Skip to content

Commit

Permalink
Merge tutorials
Browse files Browse the repository at this point in the history
  • Loading branch information
@teubert
teubert committed Oct 20, 2023
1 parent ce5659e commit a2e6d05
Showing 1 changed file with 369 additions and 0 deletions.
369 changes: 369 additions & 0 deletions tutorial.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1374,6 +1374,368 @@
"## 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",
Expand All @@ -1395,6 +1757,13 @@
"* `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."
]
Expand Down

0 comments on commit a2e6d05

Please sign in to comment.