{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Welcome to the Prognostics Python Package (ProgPy) Tutorial"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"The goal of this notebook is to instruct users on how to use and extend the NASA PCoE Python Prognostics Python Package (progpy). \n",
"\n",
"First some background. The Prognostics Python Package is a python package for the modeling, simulation, state estimation, and prediction of the evolution of state for components, systems, and systems of systems, with a focus on simulating specific events. When used for prognostics, these events are typically system failures, such as a winding failure on a motor or full discharge of a battery. \n",
"\n",
"A few definitions:\n",
"* __Event__: Something that can be predicted (e.g., system failure, warning threshold). An event has either occurred or not. \n",
"* __Event State__: Progress towards event occurring. Defined as a number where an event state of 0 indicates the event has occurred and 1 indicates no progress towards the event (i.e., fully healthy operation for a failure event). For gradually occurring events (e.g., discharge) the number will progress from 1 to 0 as the event nears. In prognostics, event state is frequently called \"State of Health\" or \"SOH\"\n",
"* __Inputs__: Control applied to the system being modeled (e.g., current drawn from a battery)\n",
"* __Outputs__: Measured sensor values from a system (e.g., voltage and temperature of a battery), outputs can be estimated from system state\n",
"* __States__: Internal parameters (typically hidden states) used to represent the state of the system- can be the same as inputs/outputs but do not have to be.\n",
"* __Performance Metrics__: Performance characteristics of a system that are a function of system state, but are not directly measured. For example, a performance metric for a electric motor might be the maximum achievable torque \n",
"* __State Estimation__: The process of estimating the (possibly hidden) state of a system given sensor information on observable states\n",
"* __Prediction__: The process of predicting the evolution of a system state with time and the occurrence of events. \n",
"\n",
"The `progpy` package has the following structure\n",
"* `progpy.data_models` - package containing algorithms for data-driven models, and parent class `progpy.data_models.DataModel`\n",
"* `progpy.datasets` - package containing tools for downloading a few relevant datasets\n",
"* `progpy.loading` - package containing tools for representing different loading profiles\n",
"* `progpy.models.*` - implemented models (e.g., pump, valve, battery)\n",
"* `progpy.state_estimators.*` - State Estimation algorithms\n",
"* `progpy.predictors.*` - Prediction algorithms\n",
"* `progpy.uncertain_data.*` - Classes for representing data with uncertainty\n",
"* `progpy.utils.*` - various utility functions and classes\n",
"* `progpy.CompositeModel` - model of a system-of-systems, combining multiple interrelated models into a single model\n",
"* `progpy.EnsembleModel` - model combining multiple models of the same system into a single model\n",
"* `progpy.MixtureOfExpertsModel` - model combining multiple models of the same system into a single model where the \"best\" model is used at every step\n",
"* `progpy.LinearModel` - parent class for simple linear models\n",
"* `progpy.PrognosticsModel` - parent class for all prognostics models - defines interfaces that a model must implement, and tools for simulating said model\n",
"\n",
"In addition to the `proypy` package, the GitHub repository includes many examples illustrating how to use the package (see `examples/`), a template for implementing a new model (`prog_model_template`), a template for implementing a new state estimator (`state_estimator_template`), a template for implementing a new predictor (`predictor_template`), and this tutorial (`tutorial.ipynb`). Documentation for ProgPy can be found at <https://nasa.github.io/progpy>,\n",
"\n",
"Before you start, make sure to install progpy using the following command:\n",
"\n",
" pip install progpy\n",
"\n",
"Now let's get started with some examples"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using the included models"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"This first example is for using the included prognostics models. \n",
"\n",
"The `progpy.models` package includes implemented models, including ones for pumps, valves, batteries, and more. See <https://nasa.github.io/progpy/api_ref/progpy/IncludedModels.html> for a full description of the included models.\n",
"\n",
"First thing to do is to import the model you would like to use:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy.models import BatteryCircuit"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"This imports the BatteryCircuit model distributed with the `progpy` package. See <https://nasa.github.io/progpy/api_ref/progpy/IncludedModels.html> for details on this model.\n",
"\n",
"Next, let's create a new battery using the default settings:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"batt = BatteryCircuit()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"This creates a battery circuit model. You can also pass configuration parameters into the constructor as keyword arguments to configure the system, for example\n",
"### <center>`BatteryCircuit(qMax = 7856)`</center>\n",
"Alternatively, you can set the configuration of the system afterwards, like so:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"batt.parameters['qMax'] = 7856 \n",
"batt.parameters['process_noise'] = 0 # Note: by default there is some process noise- this turns it off. Noise will be explained later in the tutorial"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"These parameters describe the specific system (in this case battery) the model is simulating. See <https://nasa.github.io/progpy/api_ref/progpy/IncludedModels.html> for a full list of configurable model parameters. Let's use the model properties to check out some details of this model, first the model configuration:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pprint import pprint\n",
"print('Model configuration:')\n",
"pprint(batt.parameters)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"You can save or load your model configuration using pickle"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pickle\n",
"pickle.dump(batt.parameters, open('battery123.cfg', 'wb'))"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Then you can set your model configuration like below. This is useful for recording and restoring model configurations. Some users store model configuration as picked files to be restored from later."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"batt.parameters = pickle.load(open('battery123.cfg', 'rb'))"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Additionally, both pickle and JSON can be used to serialize an entire model for future use. All PrognosticsModels include functions to serialize with JSON as follows:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"save_json = batt.to_json() # Save model \n",
"serial_batt = BatteryCircuit.from_json(save_json) # Model can be called directly with serialized result"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Information is passed to and from the model using containers that function like dictionaries. The keys of the containers are specific to the model.\n",
"\n",
"Let's look at the inputs (loading) and outputs (measurements) for a BatteryCircuit model. These are the keys expected for inputs and outputs, respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('inputs: ', batt.inputs)\n",
"print('outputs: ', batt.outputs)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's look at what events we're predicting. This model only predicts one event, called EOD (End of Discharge), but that's not true for every model. See <https://nasa.github.io/progpy/models.html>"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('event(s): ', batt.events)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, let's take a look at the internal states that the model uses to represent the system:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('states: ', batt.states)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulating to a specific time"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's simulate. The following section uses the model created in the \"using the included models\" section.\n",
"\n",
"First, we define future loading. This is a function that describes how the system will be loaded as a function of time. Here we're defining a basic piecewise function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy.loading import Piecewise\n",
"# Variable (piece-wise) future loading scheme \n",
"future_loading = Piecewise(batt.InputContainer, [600, 900, 1800, 3000, float('inf')], {'i': [2, 1, 4, 2, 3]})"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that future loading can be modeled using various classes in progpy.loading or as any function f(t: float, x: StateContainer) -> InputContainer.\n",
"\n",
"At last it's time to simulate. First we're going to simulate forward 200 seconds. To do this we use the function simulate_to() to simulate to the specified time and print the results."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"time_to_simulate_to = 200\n",
"sim_config = {\n",
" 'save_freq': 20, \n",
" 'print': True # Print states - Note: is much faster without\n",
"}\n",
"(times, inputs, states, outputs, event_states) = batt.simulate_to(time_to_simulate_to, future_loading, **sim_config)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also plot the results. Here we see the temperature of the battery increase and the voltage decrease with use. This is expected. Voltage will decrease as the state of charge decreases, and temperature increases as current is drawn through the battery, until it reaches some equilibrium. Everything is very linear because the load is kept constant within the simulation window. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"inputs.plot(ylabel='Current drawn (amps)')\n",
"event_states.plot(ylabel= 'SOC')\n",
"outputs.plot(ylabel= {'v': \"voltage (V)\", 't': 'temperature (°C)'}, compact= False)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Also, note the lack of smoothness in the voltage curve. This is limited by the save_freq from sim_cfg. There is a point on the graph every 20 seconds, because that is the frequency at which we save the results. Decreasing the save frequency will result in a cleaner curve.\n",
"\n",
"The results can be further analyzed with available metrics. For example, monotonicity can be calculated for simulate_to()'s returned objects. The End of Discharge (EOD) event state (i.e., State of Charge) should be monotonic (i.e., decreasing monotonically). Note: the EOD event is specific to the battery model. Each model will simulate different events.\n",
"\n",
"The monotonicity metric indicates the degree of monotonicity where 1 is completely monotonic and 0 is perfectly non-monotonic (decreasing as much as increasing)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('monotonicity: ', event_states.monotonicity())"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Lastly, results can be stored in a container variable and be individually accessed via namedtuple syntax."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"batt_simulation = batt.simulate_to(time_to_simulate_to, future_loading, save_freq = 20)\n",
"print('times: ', batt_simulation.times) \n",
"print('\\ninputs: ', batt_simulation.inputs)\n",
"print('\\nstates: ', batt_simulation.states)\n",
"print('\\noutputs: ', batt_simulation.outputs) \n",
"print('\\nevent states: ', batt_simulation.event_states)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulating to threshold"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Instead of specifying a specific amount of time, we can also simulate until a threshold has been met using the simulate_to_threshold() method. Results can be similarly plotted and accessed via namedtuple syntax."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"options = { #configuration for this sim\n",
" 'save_freq': 100, # Frequency at which results are saved (s)\n",
" 'horizon': 5000 # Maximum time to simulate (s) - This is a cutoff. The simulation will end at this time, or when a threshold has been met, whichever is first\n",
" }\n",
"(times, inputs, states, outputs, event_states) = batt.simulate_to_threshold(future_loading, **options)\n",
"inputs.plot(ylabel='Current drawn (amps)')\n",
"event_states.plot(ylabel='SOC')\n",
"outputs.plot(ylabel= {'v': \"voltage (V)\", 't': 'temperature (°C)'}, compact= False)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"One thing to note here is that unlike the earlier plots, these plots are not smooth curves, this is because the load is piecewise, not constant. You see jumps in the plots at the times when the load changes. Also, the simulation is long enough for the temperature to reach an equilibrium. \n",
"\n",
"Default is to simulate until any threshold is met, but we can also specify which event we are simulating to (any key from model.events) for multiple event models. See `examples.sim_battery_eol` for an example of this."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Noise"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"There are three types of noise that can be added to a model in simulation, described below:\n",
"* __Process Noise__: Noise representing uncertainty in the model transition; e.g., model or model configuration uncertainty, uncertainty from simplifying assumptions. Applied during state transition\n",
"* __Measurement Noise__: Noise representing uncertainty in the measurement process; e.g., sensor sensitivity, sensor misalignments, environmental effects. Applied during estimation of outputs from states.\n",
"* __Future Loading Noise__: Noise representing uncertainty in the future loading estimates; e.g., uncertainty from incomplete knowledge of future loading. Responsibility of user to apply in supplied future loading method\n",
"\n",
"The amount of process or measurement noise is considered a property of the model and can be set using the m.parameters['process_noise'] and m.parameters['measurement_noise'], respectively.\n",
"\n",
"In this example we will use the ThrownObject model and turn off process noise. ThrownObject is a simple model of an object thrown directly into the air in a vacuum. Thrown object simulates two events: 'falling' (when the object starts falling) and 'impact' (when the object hits the ground). More details can be found later in the tutorial."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy.models import ThrownObject\n",
"\n",
"# Create an instance of the thrown object model with no process noise\n",
"m = ThrownObject(process_noise=False)\n",
"\n",
"# Define future loading\n",
"def future_load(t=None, x=None): \n",
" # The thrown object model has no inputs- you cannot load the system (i.e., effect it once it's in the air)\n",
" # So we return an empty input container\n",
" return m.InputContainer({})\n",
"\n",
"# Define configuration for simulation\n",
"config = {\n",
" 'events': 'impact', # Simulate until the thrown object has impacted the ground\n",
" 'dt': 0.005, # Time step (s)\n",
" 'save_freq': 0.5, # Frequency at which results are saved (s)\n",
"}\n",
"\n",
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, **config)\n",
"\n",
"# Print results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"print('\\nimpact time: {:.2f}s'.format(times[-1]))\n",
"# The simulation stopped at impact, so the last element of times is the impact time\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"See above the velocity decreases linearly and the position follows a clean parabola, as we would expect.\n",
"\n",
"Now with this second example we apply normal (i.e., gaussian) process noise with a standard deviation of 15 to every state. Compare the plots generated with those above- you should be able to see the effects of the noise"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = ThrownObject(process_noise = 15)\n",
"\n",
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, **config)\n",
"\n",
"# Print Results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"print('\\nimpact time: {:.2f}s'.format(times[-1]))\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also specify different amounts of noise on different states, for example here we are applying no noise to velocity but a large amount of noise to the position. Compare the plot with that above. Here you should see a smooth curve for the velocity, but a noisy curve for position."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = ThrownObject(process_noise = {'x': 50, 'v': 0})\n",
"\n",
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, **config)\n",
"\n",
"# Print Results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"print('\\nimpact time: {:.2f}s'.format(times[-1]))\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also define the shape of the noise to be uniform or triangular instead of normal. Finally, you can define your own function to apply noise for anything more complex. \n",
"\n",
"For more information see `examples.noise`"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulating - advanced"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"There are a number of advanced features that can be used in simulation. A few of these will be described below. Detail can also be found in the documentation here: https://nasa.github.io/progpy/prog_models_guide.html#simulation"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Saving results at specific points"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Earlier, we demonstrated how save_freq specifies the frequency at which the results are saved in simulation. There are occasionally circumstances when you need the results at a specific time (e.g., check points, waypoints, etc.). This is accomplished with the save_pts argument in simulation.\n",
"\n",
"For example, in the following simple case, we're reusing the example from earlier, but we're saving the results at 1 second, 5 seconds, 6 seconds, and 7 seconds. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Reset noise\n",
"m.parameters['process_noise'] = False\n",
"m.parameters['measurement_noise'] = False\n",
"\n",
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, dt=config['dt'], events='impact', save_pts=[1, 5, 6, 7])\n",
"\n",
"# Print Results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that now we only have the data at the points specified (plus initial and final). As a result, the graph looks lopsided.\n",
"\n",
"save_pts can also be used to adjust the save frequency as degradation progresses, for example by having more concentrated save pts around the end of degradation, like the example below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, dt=config['dt'], events='impact', save_pts=[1, 2, 3, 4, 5, 5.5, 6, 6.25, 6.5, 6.75, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8])\n",
"\n",
"# Print Results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the curve is much better defined as the simulation progresses. This is frequently of interest since users are mostly interested in behavior close to the event occurrence."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Keep in mind that using a fixed dt, like above, will only save the data at the closest dt increment at or after the save point. So if save_freq is less than dt, save points will be missed. Or, if save_pts are not a multiple of dt, the resulting point at which data is saved will not be exactly at the save point. Auto step sizes (see next section) can help with this."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Step Sizes"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"At the end of the save_pts example above, we noted that save_pts are often not exactly met when they are not a multiple of dt. This can be seen in the example below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, dt=0.2, events='impact', save_pts=[1, 2.5, 3, 4.5, 6, 7.5])\n",
"\n",
"# Print Results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Note above that the state was captured at 2.6, 4.6, and 7.6 instead of the requested 2.5, 4.5, and 7.5. \n",
"\n",
"This can be fixed by auto step sizes. Let's repeat this example with automatic step sizes."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, dt='auto', events='impact', save_pts=[1, 2.5, 3, 4.5, 6, 7.5])\n",
"\n",
"# Print Results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the save_pts are captured exactly, but note that the resulting behavior is very different than the results of previous simulations. This is because we used an unbounded automatic step size. The simulation took larger steps than we should, resulting in error accumulation.\n",
"\n",
"With bounded automatic step sizes we also provide the largest step it can take. This prevents the integrator from taking too large steps."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Simulate to a threshold\n",
"(times, _, states, outputs, _) = m.simulate_to_threshold(future_load, dt=('auto', 0.2), events='impact', save_pts=[1, 2.5, 3, 4.5, 6, 7.5])\n",
"\n",
"# Print Results\n",
"print('states:')\n",
"for (t,x) in zip(times, states):\n",
" print('\\t{:.2f}s: {}'.format(t, x))\n",
"\n",
"# Plot results\n",
"states.plot()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the simulation still captures every save point exactly, but it never takes a step size larger than 0.2 seconds, resulting in a better simulation."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also define more complex dynamic step sizes by providing a function (t, x)->dt. For example: `dt=lambda t, x: max(0.75 - t/100, 0.25)`"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Integration methods"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Simulation is essentially the process of integrating the model forward with time. By default, simple Euler integration is used to propagate the model forward. Advanced users can change the numerical integration method to affect the simulation quality and runtime. This is done using the `integration_method` argument in `simulate_to()`, `simulate_to_threshold()`, or the model parameters. For example:\n",
"\n",
"`m.parameters['integration_method'] = 'rk4'`\n",
"\n",
"Note: integration method can only be changed for continuous models. In this case our thrown object model is discrete, so we cannot demonstrate it on this model (see below)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m.is_continuous"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Integration methods can also be set to scipy.integrate functions, for example:\n",
"\n",
"`m.parameters['integration_method'] = scipy.integrate.Radau`\n",
"\n",
"Any keyword arguments are then saved into `m.parameters['integrator_config']` as a dictionary"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Building a new model"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"To build a model, create a seperate class to define the logic of the model. Do this by copying the file `prog_model_template.py` and replacing the functions with the logic specific to your model. \n",
"\n",
"For this example, we will model the throwing of an object directly into the air in a vacuum. This is not a particularly interesting problem, but it is simple and illustrates the basic methods of a PrognosticsModel.\n",
"\n",
"The model is illustrated below:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy import PrognosticsModel\n",
"\n",
"class ThrownObject(PrognosticsModel):\n",
" \"\"\"\n",
" Model that simulates an object thrown directly into the air (vertically) without air resistance\n",
" \"\"\"\n",
"\n",
" inputs = [] # no inputs, no way to control\n",
" states = [\n",
" 'x', # Vertical position (m) \n",
" 'v' # Velocity (m/s)\n",
" ]\n",
" outputs = [ # Anything we can measure\n",
" 'x' # Position (m)\n",
" ]\n",
" events = [ # Events that can/will occur during simulation\n",
" 'falling', # Event- object is falling\n",
" 'impact' # Event- object has impacted ground\n",
" ]\n",
"\n",
" # The Default parameters for any ThrownObject. \n",
" # Overwritten by passing parameters into constructor as kwargs or by setting model.parameters\n",
" default_parameters = {\n",
" 'thrower_height': 1.83, # Height of thrower (m)\n",
" 'throwing_speed': 40, # Velocity at which the ball is thrown (m/s)\n",
" 'g': -9.81, # Acceleration due to gravity (m/s^2)\n",
" 'process_noise': 0.0 # amount of noise in each step\n",
" } \n",
"\n",
" # First function: Initialize. This function is used to initialize the first state of the model.\n",
" # This is only needed in cases with complex initialization. \n",
" # In this case we need it because of max_x\n",
" # If not included, parameters['x0'] is returned as the initial state\n",
" # In this case we do not need input (u) or output (z) to initialize the model, \n",
" # so we set them to None, but that's not true for every model.\n",
" # u and z are Input and Output, respectively. \n",
" # Values can be accessed like a dictionary (e.g., z['x']) using the keys from inputs and outputs, respectively.\n",
" # or they can be accessed using the matrix (i.e., z.matrix)\n",
" def initialize(self, u=None, z=None):\n",
" self.max_x = 0.0\n",
" return self.StateContainer({\n",
" 'x': self.parameters['thrower_height'], # initial altitude is height of thrower\n",
" 'v': self.parameters['throwing_speed'] \n",
" })\n",
" \n",
" # Second function: state transition. \n",
" # State transition can be defined in one of two ways:\n",
" # 1) Discrete models use next_state(x, u, dt) -> x'\n",
" # 2) Continuous models (preferred) use dx(x, u) -> dx/dt\n",
" #\n",
" # In this case we choose the continuous model, so we define dx(x, u)\n",
" # This function defines the first derivative of the state with respect to time, as a function of model configuration (self.parameters), state (x) and input (u).\n",
" # Here input isn't used. But past state and configuration are.\n",
" # \n",
" # x and u are State and Input, respectively. \n",
" # Values can be accessed like a dictionary (e.g., x['x']) using the keys from states and inputs, respectively.\n",
" # or they can be accessed using the matrix (i.e., x.matrix)\n",
" def dx(self, x, u):\n",
" return self.StateContainer({\n",
" 'x': x['v'], # dx/dt = v\n",
" 'v': self.parameters['g'] # Acceleration of gravity\n",
" })\n",
" # Equivalently, the state transition could have been defined as follows:\n",
" # def next_state(self, x, u, dt):\n",
" # return self.StateContainer({\n",
" # 'x': x['x'] + x['v']*dt,\n",
" # 'v': x['v'] + self.parameters['g']*dt\n",
" # })\n",
"\n",
" # Now, we define the output equation. \n",
" # This function estimates the output (i.e., measured values) given the system state (x) and system parameters (self.parameters).\n",
" # In this example, we're saying that the state 'x' can be directly measured. \n",
" # But in most cases output will have to be calculated from state. \n",
" def output(self, x):\n",
" return self.OutputContainer({\n",
" 'x': x['x']\n",
" })\n",
"\n",
" # Next, we define the event state equation\n",
" # This is the first equation that actually describes the progress of a system towards the events.\n",
" # This function maps system state (x) and system parameters (self.parameters) to event state for each event.\n",
" # Event state is defined as a number between 0 and 1 where 1 signifies no progress towards the event, and 0 signifies the event has occurred.\n",
" # The event keys were defined above (model.events)\n",
" # Here the two event states are as follows:\n",
" # 1) falling: 1 is defined as when the system is moving at the maximum speed (i.e., throwing_speed), and 0 is when velocity is negative (i.e., falling)\n",
" # 2) impact: 1 is defined as the ratio of the current altitude (x) to the maximum altitude (max_x), and 0 is when the current altitude is 0 (i.e., impact)\n",
" def event_state(self, x): \n",
" self.max_x = max(self.max_x, x['x']) # Maximum altitude\n",
" return {\n",
" 'falling': max(x['v']/self.parameters['throwing_speed'],0), # Throwing speed is max speed\n",
" 'impact': max(x['x']/self.max_x,0) # Ratio of current altitude to maximum altitude\n",
" }\n",
"\n",
" # Finally, we define the threshold equation.\n",
" # This is the second equation that describes the progress of a system towards the events.\n",
" # Note: This function is optional. If not defined, threshold_met will be defined as when the event state is 0.\n",
" # However, this implementation is more efficient, so we included it\n",
" # This function maps system state (x) and system parameters (self.parameters) a boolean indicating if the event has been met for each event.\n",
" def threshold_met(self, x):\n",
" return {\n",
" 'falling': x['v'] < 0,\n",
" 'impact': x['x'] <= 0\n",
" }\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the model can be generated and used like any of the other provided models"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = ThrownObject()\n",
"\n",
"def future_load(t, x=None):\n",
" return m.InputContainer({}) # No loading\n",
"event = 'impact' # Simulate until impact\n",
"\n",
"(times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_load, events=[event], dt=0.005, save_freq=1)\n",
"\n",
"# Plot results\n",
"event_states.plot(ylabel= ['falling', 'impact'], compact= False)\n",
"states.plot(ylabel= {'x': \"position (m)\", 'v': 'velocity (m/s)'}, compact= False)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Note: The plots are at the resolution of save_freq (one point per second)\n",
"\n",
"See also `examples.new_model` for more information on building models"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Building a new model - advanced"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Derived Parameters \n",
"\n",
"Models can also include \"derived parameters\" (i.e., parameters that are derived from others). These can be set using the param_callbacks property. \n",
"\n",
"Let's extend the above model to include derived_parameters. Let's say that the throwing_speed was actually a function of thrower_height (i.e., a taller thrower would throw the ball faster). Here's how we would implement that"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Step 1: Define a function for the relationship between thrower_height and throwing_speed.\n",
"def update_thrown_speed(params):\n",
" return {\n",
" 'throwing_speed': params['thrower_height'] * 21.85\n",
" } # Assumes thrown_speed is linear function of height\n",
"\n",
"# Step 2: Define the param callbacks\n",
"ThrownObject.param_callbacks = {\n",
" 'thrower_height': [update_thrown_speed]\n",
" } # Tell the derived callbacks feature to call this function when thrower_height changes."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Note: Usually we would define this method within the class. For this example, we're doing it separately to improve readability. Here's the feature in action"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"obj = ThrownObject()\n",
"print(\"Default Settings:\\n\\tthrower_height: {}\\n\\tthowing_speed: {}\".format(obj.parameters['thrower_height'], obj.parameters['throwing_speed']))\n",
"\n",
"# Now let's change the thrower_height\n",
"print(\"changing height...\")\n",
"obj.parameters['thrower_height'] = 1.75 # Our thrower is 1.75 m tall\n",
"print(\"\\nUpdated Settings:\\n\\tthrower_height: {}\\n\\tthowing_speed: {}\".format(obj.parameters['thrower_height'], obj.parameters['throwing_speed']))\n",
"print(\"Notice how speed changed automatically with height\")\n",
"\n",
"# Let's delete the callback so we can use the same model in the future:\n",
"ThrownObject.param_callbacks = {}"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### State Limits\n",
"\n",
"In many cases, the values of the model states have certain physical limits. For example, temperature is limited by absolute zero. In these cases, it is useful to define those limits. In simulation, the defined limits are enforced as the state transitions to prevent the system from reaching an impossible state."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"For example, the ThrownObject model can be extended as follows:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from numpy import inf\n",
"\n",
"ThrownObject.state_limits = {\n",
" # object may not go below ground\n",
" 'x': (0, inf),\n",
"\n",
" # object may not exceed the speed of light\n",
" 'v': (-299792458, 299792458)\n",
" }"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Note: like derived parameters, these would typically be defined in class definition, not afterwards. They are defined afterwards in this case to illustrate the feature.\n",
"\n",
"State limits can be applied directly using the apply_limits function. For example:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x = {'x': -5, 'v': 3e8} # Too fast and below the ground\n",
"x = obj.apply_limits(x)\n",
"print(x)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice how the state was limited according to the model state limits and a warning was issued. The warning can be suppressed by suppressing ProgModelStateLimitWarning (`warnings.simplefilter('ignore', ProgModelStateLimitWarning)`)\n",
"\n",
"See also examples.derived_params for more information on this feature."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Parameter Estimation"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's say we don't know the configuration of the above model. Instead, we have some data. We can use that data to estimate the parameters. \n",
"\n",
"First, we define the data:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"times = [0, 1, 2, 3, 4, 5, 6, 7, 8]\n",
"inputs = [{}]*9\n",
"outputs = [\n",
" {'x': 1.83},\n",
" {'x': 36.95},\n",
" {'x': 62.36},\n",
" {'x': 77.81},\n",
" {'x': 83.45},\n",
" {'x': 79.28},\n",
" {'x': 65.3},\n",
" {'x': 41.51},\n",
" {'x': 7.91},\n",
"]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we identify which parameters will be optimized"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"keys = ['thrower_height', 'throwing_speed']"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's say we have a first guess that the thrower's height is 20m, silly I know"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m.parameters['thrower_height'] = 20"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the state of our estimation with that assumption:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('Model configuration before')\n",
"for key in keys:\n",
" print(\"-\", key, m.parameters[key])"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also calculate the error between simulated and observed data given this incorrect assumption. The calc_error functionality can calculate error using a variety of methods, including Mean Squared Error, which we use here. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(' Error: ', m.calc_error(times, inputs, outputs, dt=1e-4, method='MSE'))"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Wow, that's a large error. \n",
"\n",
"Let's run the parameter estimation to see if we can do better"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m.estimate_params([(times, inputs, outputs)], keys, dt=0.01)\n",
"\n",
"# Print result\n",
"print('\\nOptimized configuration')\n",
"for key in keys:\n",
" print(\"-\", key, m.parameters[key])\n",
"print(' Error: ', m.calc_error(times, inputs, outputs, dt=1e-4))"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Much better!\n",
"\n",
"See also examples.param_est for more details about this feature"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## UncertainData - Representing a Distribution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Uncertainty is sometimes present in data used for performing state estimations or making predictions.\n",
"\n",
"In `progpy`, data with uncertainty is represented using classes inheriting from `UncertainData`:\n",
"* `progpy.uncertain_data.MultivariateNormalDist` - Data represented by a multivariate normal distribution with mean and covariance matrix\n",
"* `progpy.uncertain_data.ScalarData` - Data without uncertainty, a single value\n",
"* `progpy.uncertain_data.UnweightedSamples` - Data represented by a set of unweighted samples. Objects of this class can be treated like a list where samples[n] returns the nth sample (Dict)\n",
"\n",
"To begin using `UncertainData`, import the type that best portrays the data. In this simple demonstration, we import the `UnweightedSamples` data type. See <https://nasa.github.io/progpy/api_ref/progpy/UncertainData.html> for full details on the available `UncertainData` types."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy.uncertain_data import UnweightedSamples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With `UnweightedSamples` imported, construct an object with samples. This object can be initialized using either a dictionary, list, or model.*Container type from prog_models (e.g., StateContainer). Let's try creating an object using a dictionary. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"samples = UnweightedSamples([{'x': 1, 'v':2}, {'x': 3, 'v':-2}])\n",
"print(samples)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given an integer value, addition and subtraction can be performed on the `UncertainData` classes to adjust the distribution by a scalar amount."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"samples = samples + 5\n",
"print(samples)\n",
"samples -= 3\n",
"print(samples)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also sample from any `UncertainData` distribution using the `sample` method. In this case it resamples from the existing samples"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(samples.sample()) # A single sample\n",
"print(samples.sample(10)) # 10 samples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see the keys present using the `.keys()` method:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(samples.keys())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and the data corresponding to a specific key can be retrieved using `.key()`"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(samples.key('x'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Various properties are available to quantify the `UncertainData` distribution"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('mean', samples.mean)\n",
"print('median', samples.median)\n",
"print('covariance', samples.cov)\n",
"print('size', samples.size)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These `UncertainData` classes are used throughout the prog_algs package to represent data with uncertainty, as described in the following sections."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## State Estimation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"State estimation is the process of estimating the system state given sensor data and a model. Typically, this is done repeatedly as new sensor data is available.\n",
"\n",
"In `progpy` a State Estimator is used to estimate the system state. \n",
"\n",
"To start, import the needed packages. Here we will import the `BatteryCircuit` model and the `UnscentedKalmanFilter` state estimator. See <https://nasa.github.io/progpy/api_ref/progpy/StateEstimator.html> for more details on the available state estimators."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy.models import BatteryCircuit\n",
"from progpy.state_estimators import UnscentedKalmanFilter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we construct and initialize the model. \n",
"\n",
"We use the resulting model and initial state to construct the state estimator. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = BatteryCircuit()\n",
"x0 = m.initialize()\n",
"\n",
"# Turn into a distribution - this represents uncertainty in the initial state\n",
"from progpy.uncertain_data import MultivariateNormalDist\n",
"from numpy import diag\n",
"INITIAL_UNCERT = 0.05 # Uncertainty in initial state (%)\n",
"# Construct covariance matrix (making sure each value is positive)\n",
"cov = diag([max(abs(INITIAL_UNCERT * value), 1e-9) for value in x0.values()])\n",
"x0 = MultivariateNormalDist(x0.keys(), x0.values(), cov)\n",
"\n",
"# Construct State estimator\n",
"est = UnscentedKalmanFilter(m, x0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use the estimator to estimate the system state."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"Prior State:\", est.x.mean)\n",
"print('\\tSOC: ', m.event_state(est.x.mean)['EOD'])\n",
"fig = est.x.plot_scatter(label='prior')\n",
"\n",
"t = 0.1\n",
"u = m.InputContainer({'i': 2})\n",
"example_measurements = m.OutputContainer({'t': 32.2, 'v': 3.915})\n",
"est.estimate(t, u, example_measurements)\n",
"\n",
"print(\"Posterior State:\", est.x.mean)\n",
"print('\\tSOC: ', m.event_state(est.x.mean)['EOD'])\n",
"est.x.plot_scatter(fig= fig, label='posterior')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As mentioned previously, this step is typically repeated when there's new data. filt.x may not be accessed every time the estimate is updated, only when it's needed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Prediction Example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Prediction is the practice of using a state estimation, a model, and estimates of future loading to predict future states and when an event will occur.\n",
"\n",
"First we will import a predictor. In this case, we will use the MonteCarlo Predictor, but see documentation <https://nasa.github.io/progpy> for a full list of predictors and their configuration parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy.predictors import MonteCarlo"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we initialize it using the model from the above example"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mc = MonteCarlo(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, let's define future loading and the first state. The first state is the output of the state estimator, and the future loading scheme is a simple piecewise function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x = est.x # The state estimate\n",
"\n",
"from progpy.loading import Piecewise\n",
"future_load = Piecewise(\n",
" m.InputContainer,\n",
" [600, 900, 1800, 3000],\n",
" {'i': [2, 1, 4, 2, 3]}\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's use the constructed mc predictor to perform a single prediction. Here we're setting dt to 0.25. Note this may take up to a minute"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mc_results = mc.predict(x, future_loading, dt=0.25, n_samples=20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The predict function returns predictions of future inputs, states, outputs, and event_states at each save point. For sample-based predictors like the monte carlo, these can be accessed like an array with the format `[sample #][time]` so that `mc_results.states[m][n]` corresponds to the state for sample `m` at time `mc_results.times[m][n]`. Alternately, use the method `snapshot` to get a single point in time. e.g., \n",
"\n",
" state = mc_results.states.snapshot(3)\n",
"\n",
"In this case the state snapshot corresponds to time `mc_results.times[3]`. The snapshot method returns type UncertainData. \n",
"\n",
"The `predict` method also returns Time of Event (ToE) as a type UncertainData, representing the predicted time of event (for each event predicted), with uncertainty.\n",
"\n",
"Next, let's use the metrics package to analyze the ToE"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"\\nEOD Predictions (s):\")\n",
"print('\\tPortion between 3005.2 and 3005.6: ', mc_results.time_of_event.percentage_in_bounds([3005.2, 3005.6]))\n",
"print('\\tAssuming ground truth 3005.25: ', mc_results.time_of_event.metrics(ground_truth = 3005.25))\n",
"from progpy.metrics import prob_success \n",
"print('\\tP(Success) if mission ends at 3005.25: ', prob_success(mc_results.time_of_event, 3005.25))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These analysis methods applied to ToE can also be applied to anything of type UncertainData (e.g., state snapshot). \n",
"\n",
"You can also visualize the results in a variety of different ways. For example, state transition"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = mc_results.states.snapshot(0).plot_scatter(label = \"t={:.0f}\".format(int(mc_results.times[0])))\n",
"for i in range(1, 4):\n",
" index = int(len(mc_results.times)/4*i)\n",
" mc_results.states.snapshot(index).plot_scatter(fig=fig, label = \"t={:.0f}\".format(mc_results.times[index]))\n",
"mc_results.states.snapshot(-1).plot_scatter(fig = fig, label = \"t={:.0f}\".format(int(mc_results.times[-1])))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Or time of event (ToE)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = mc_results.time_of_event.plot_hist()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note, for this event, there is only one event (End of Discharge). Many models have multiple events that can be predicted. For these models, ToE for each event is returned and can be analyzed.\n",
"\n",
"Alternately, a specific event (or events) can be specified for prediction. See `examples.predict_specific_event` for more details.\n",
"\n",
"Frequently the prediction step is run periodically, less often than the state estimator step"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extending - Adding a new state estimator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"New state estimators can be created by extending the state_estimator interface. As an example lets use a really dumb state estimator that adds random noise each step - and accepts the state that is closest. \n",
"\n",
"First thing we need to do is import the StateEstimator parent class"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from progpy.state_estimators import StateEstimator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we select how state will be represented. In this case there's no uncertainty- so we represent it as a scaler. Import the appropriate class"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from prog_algs.uncertain_data import ScalarData"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we construct the class, implementing the functions of the state estimator template (`state_estimator_template.py`)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import random \n",
"\n",
"class BlindlyStumbleEstimator(StateEstimator):\n",
" def __init__(self, model, x0):\n",
" self.m = model\n",
" self.state = x0\n",
"\n",
" def estimate(self, t, u, z):\n",
" # Generate new candidate state\n",
" x2 = {key : float(value) + 10*(random.random()-0.5) for (key,value) in self.state.items()}\n",
"\n",
" # Calculate outputs\n",
" z_est = self.m.output(self.state)\n",
" z_est2 = self.m.output(x2)\n",
"\n",
" # Now score them each by how close they are to the measured z\n",
" z_est_score = sum([abs(z_est[key] - z[key]) for key in self.m.outputs])\n",
" z_est2_score = sum([abs(z_est2[key] - z[key]) for key in self.m.outputs])\n",
"\n",
" # Now choose the closer one\n",
" if z_est2_score < z_est_score: \n",
" self.state = x2\n",
"\n",
" @property\n",
" def x(self):\n",
" return ScalarData(self.state)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Great, now let's try it out using the model from earlier. with an initial state of all 0s. It should slowly converge towards the correct state"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x0 = {key: 0 for key in m.states}\n",
"se = BlindlyStumbleEstimator(m, x0)\n",
"\n",
"for i in range(25):\n",
" u = m.InputContainer({'i': 0})\n",
" z = m.OutputContainer({'t': 18.95, 'v': 4.183})\n",
" se.estimate(i, u, z)\n",
" print(se.x.mean)\n",
" print(\"\\tcorrect: {'tb': 18.95, 'qb': 7856.3254, 'qcp': 0, 'qcs': 0}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extending - Adding a new Predictor"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Like the example above with StateEstimators, Predictors can be extended by subclassing the Predictor class. Copy `predictor_template.py` as a starting point."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusion\n",
"This is just the basics, there's much more to learn. Please see the documentation at <https://nasa.github.io/progpy/guide> and the examples in the `examples/` folder for more details on how to use the package, including:\n",
"* `examples.derived_params` : Example building models with derived parameters.\n",
"* `examples.state_limits`: Example building models with limits on state variables.\n",
"* `examples.param_est`: Example using the parameter estimation feature \n",
"* `examples.dynamic_step_size`: Example simulating with dynamic (i.e., changing as a function of time or state) step size.\n",
"* `examples.events`: Example extending a model to include additional events, such as warning thresholds.\n",
"* `examples.generate_surrogate`: Example generating a surrogate model\n",
"* `examples.linear_model`: Example using the new Linear Model subclass\n",
"* `examples.benchmarking`: Example benchmarking the performance of a Prognostics Model\n",
"* `examples.future_loading`: Example with complex future loading schemes\n",
"* `examples.new_model`: Example building a new model\n",
"* `examples.noise`: Example demonstrating how noise can be added in simulation\n",
"* `examples.vectorized`: Example simulating a vectorized model\n",
"* `examples.sim`, `examples.sim_battery_eol`, `examples.sim_pump`, `examples.sim_valve`, `examples.sim_powertrain`, `examples.sim_dcmotor_singlephase`: Examples using specific models from `progpy.models`\n",
"* `examples.lstm_model`, `examples.full_lstm_model`, and `examples.custom_model`: Examples using data-driven models\n",
"* `examples.basic_example` : A basic Example using prog_algs for Prognostics \n",
"* `examples.basic_example_battery` : A basic Example using prog_algs for Prognostics, using the more complex battery model\n",
"* `examples.eol_event` : An example where a model has multiple events, but the user is only interested in predicting the time when the first event occurs (whatever it is).\n",
"* `examples.measurement_eqn_example` : An example where not every output is measured or measurements are not in the same format as outputs, so a measurement equation is defined to translate between outputs and what's measured. \n",
"* `examples.new_state_estimator_example` : An example of extending StateEstimator to create a new state estimator class\n",
"* `examples.playback` : A full example performing prognostics using playback data.\n",
"* `examples.predict_specific_event` : An example where the model has multiple events, but the user is only interested in predicting a specific event (or events).\n",
"\n",
"Thank you for trying out this tutorial. Open an issue on github (https://github.com/nasa/progpy/issues) or email Chris Teubert (christopher.a.teubert@nasa.gov) with any questions or issues."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.11.0 64-bit",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.0"
},
"metadata": {
"interpreter": {
"hash": "c1e35f02e3a88578371dd5b5d88a204463a98b2d7cd5222657e170520db47be1"
}
},
"orig_nbformat": 2,
"vscode": {
"interpreter": {
"hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}