From 06e7da05a00de01f2c65541eb611454c877e25af Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Mon, 30 Oct 2023 09:28:38 -0700 Subject: [PATCH] Add back example #4 --- examples/04_New Models.ipynb | 1661 ++++++++++++++++++++++++++++++++++ 1 file changed, 1661 insertions(+) create mode 100644 examples/04_New Models.ipynb diff --git a/examples/04_New Models.ipynb b/examples/04_New Models.ipynb new file mode 100644 index 0000000..0cd7256 --- /dev/null +++ b/examples/04_New Models.ipynb @@ -0,0 +1,1661 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 4. Defining new Physics-Based Prognostic Models" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All of the past sections describe how to use an existing model. In this section we will describe how to create a new model. This section specifically describes creating a new physics-based model. For training and creating data-driven models see 5. Data-driven Models." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Linear Models" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The easiest model to build is a linear model. Linear models are defined as a linear time series, which can be defined by the following equations:\n", + "\n", + "\n", + "\n", + "**The State Equation**:\n", + "$$\n", + "\\frac{dx}{dt} = Ax + Bu + E\n", + "$$\n", + "\n", + "**The Output Equation**:\n", + "$$\n", + "z = Cx + D\n", + "$$\n", + "\n", + "**The Event State Equation**:\n", + "$$\n", + "es = Fx + G\n", + "$$\n", + "\n", + "$x$ is `state`, $u$ is `input`, $z$ is `output`, and $es$ is `event state`\n", + "\n", + "Linear Models are defined by creating a new model class that inherits from progpy's LinearModel class and defines the following properties:\n", + "* $A$: 2-D np.array[float], dimensions: n_states x n_states. The state transition matrix. It dictates how the current state affects the change in state dx/dt.\n", + "* $B$: 2-D np.array[float], optional (zeros by default), dimensions: n_states x n_inputs. The input matrix. It dictates how the input affects the change in state dx/dt.\n", + "* $C$: 2-D np.array[float], dimensions: n_outputs x n_states. The output matrix. It determines how the state variables contribute to the output.\n", + "* $D$: 1-D np.array[float], optional (zeros by default), dimensions: n_outputs x 1. A constant term that can represent any biases or offsets in the output.\n", + "* $E$: 1-D np.array[float], optional (zeros by default), dimensions: n_states x 1. A constant term, representing any external effects that are not captured by the state and input.\n", + "* $F$: 2-D np.array[float], dimensions: n_es x n_states. The event state matrix, dictating how state variables contribute to the event state.\n", + "* $G$: 1-D np.array[float], optional (zeros by default), dimensions: n_es x 1. A constant term that can represent any biases or offsets in the event state.\n", + "* __inputs__: list[str] - `input` keys\n", + "* __states__: list[str] - `state` keys\n", + "* __outputs__: list[str] - `output` keys\n", + "* __events__: list[str] - `event` keys" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will now utilize our LinearModel to model the classical physics problem throwing an object into the air. This is a common example model, the non-linear version of which (`progpy.examples.ThrownObject`) has been used frequently throughout the examples. This version of ThrownObject will behave nearly identically to the non-linear ThrownObject, except it will not have the non-linear effects of air resistance.\n", + "\n", + "We can create a subclass of LinearModel which will be used to simulate an object thrown, which we will call the ThrownObject Class.\n", + "\n", + "First, some definitions for our Model:\n", + "\n", + "**Events**: (2)\n", + "* `falling: The object is falling`\n", + "* `impact: The object has hit the ground`\n", + "\n", + "**Inputs/Loading**: (0)\n", + "* `None`\n", + "\n", + "**States**: (2)\n", + "* `x: Position in space (m)`\n", + "* `v: Velocity in space (m/s)`\n", + "\n", + "**Outputs/Measurements**: (1)\n", + "* `x: Position in space (m)`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, for our keyword arguments:\n", + "\n", + "* __thrower_height : Optional, float__\n", + " * Height of the thrower (m). Default is 1.83 m\n", + "* __throwing_speed : Optional, float__\n", + " * Speed at which the ball is thrown (m/s). Default is 40 m/s" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With our definitions, we can now create the ThrownObject Model.\n", + "\n", + "First, we need to import the necessary packages." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from progpy import LinearModel" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we'll define some features of a ThrownObject LinearModel. Recall that all LinearModels follow a set of core equations and require some specific properties (see above). In the next step, we'll define our inputs, states, outputs, and events, along with the $A$, $C$, $E$, and $F$ values.\n", + "\n", + "First, let's consider state transition. For an object thrown into the air without air resistance, velocity would decrease linearly by __-9.81__ \n", + "$\\dfrac{m}{s^2}$ due to the effect of gravity, as described below:\n", + "\n", + " $$\\frac{dv}{dt} = -9.81$$\n", + "\n", + " Position change is defined by velocity (v), as described below:\n", + " \n", + " $$\\frac{dx}{dt} = v$$\n", + "\n", + " Note: For the above equation x is position not state. Combining these equations with the model $\\frac{dx}{dt}$ equation defined above yields the A and E matrix defined below. Note that there is no B defined because this model does not have any inputs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject(LinearModel):\n", + " events = ['impact']\n", + " inputs = [] \n", + " states = ['x', 'v']\n", + " outputs = ['x']\n", + " \n", + " A = np.array([[0, 1], [0, 0]])\n", + " C = np.array([[1, 0]])\n", + " E = np.array([[0], [-9.81]])\n", + " F = None" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that we defined our `A`, `C`, `E`, and `F` values to fit the dimensions that were stated at the beginning of the notebook! Since the parameter `F` is not optional, we have to explicitly set the value as __None__.\n", + "\n", + "Next, we'll define some default parameters for our ThrownObject model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject(ThrownObject): # Continue the ThrownObject class\n", + " default_parameters = {\n", + " 'thrower_height': 1.83,\n", + " 'throwing_speed': 40,\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the following cells, we'll define some class functions necessary to perform prognostics on the model.\n", + "\n", + "The `initialize()` function sets the initial system state. Since we have defined the `x` and `v` values for our ThrownObject model to represent position and velocity in space, our initial values would be the thrower_height and throwing_speed parameters, respectively." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject(ThrownObject):\n", + " def initialize(self, u=None, z=None):\n", + " return self.StateContainer({\n", + " 'x': self.parameters['thrower_height'],\n", + " 'v': self.parameters['throwing_speed']\n", + " })" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For our `threshold_met()`, we define the function to return True for event 'falling' when our thrown object model has a velocity value of less than 0 (object is 'falling') and for event 'impact' when our thrown object has a distance from of the ground of less than or equal to 0 (object is on the ground, or has made 'impact').\n", + "\n", + "`threshold_met()` returns a _dict_ of values, if each entry of the _dict_ is __True__, then our threshold has been met!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject(ThrownObject):\n", + " def threshold_met(self, x):\n", + " return {\n", + " 'falling': x['v'] < 0,\n", + " 'impact': x['x'] <= 0\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, for our `event_state()`, we will calculate the measurement of progress towards the events. We normalize our values such that they are in the range of 0 to 1, where 0 means the event has occurred." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject(ThrownObject):\n", + " def event_state(self, x): \n", + " x_max = x['x'] + np.square(x['v'])/(9.81*2)\n", + " return {\n", + " 'falling': np.maximum(x['v']/self.parameters['throwing_speed'],0),\n", + " 'impact': np.maximum(x['x']/x_max,0) if x['v'] < 0 else 1\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With these functions created, we can now use the `simulate_to_threshold()` function to simulate the movement of the thrown object in air. For more information, see 1. Simulation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = ThrownObject()\n", + "save = m.simulate_to_threshold(print=True, save_freq=1, threshold_keys='impact', dt=0.1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Note__: Because our model takes in no inputs, we have no need to define a future loading function. However, for most models, there would be inputs, and thus a need for a future loading function. For more information on future loading functions and when to use them, please refer to the future loading section in 1. Simulation.\n", + "\n", + "Let's take a look at the outputs of this model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = save.outputs.plot(title='generated model')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that that plot resembles a parabola, which represents the position of the ball through space as time progresses!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For more information on Linear Models, see the [Linear Model](https://nasa.github.io/progpy/api_ref/prog_models/LinearModel.html) Documentation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## New State Transition Models" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the previous section, we defined a new prognostic model using the LinearModel class. This can be a powerful tool for defining models that can be described as a linear time series. Physics-based state transition models that cannot be described linearly are constructed by subclassing [progpy.PrognosticsModel](https://nasa.github.io/progpy/api_ref/prog_models/PrognosticModel.html#prog_models.PrognosticsModel). To demonstrate this, we'll create a new model class that inherits from this class. Once constructed in this way, the analysis and simulation tools for PrognosticsModels will work on the new model.\n", + "\n", + "For this example, we'll create a simple state-transition model of an object thrown upward into the air without air resistance. Note that this is the same dynamic system as the linear model example above, but formulated in a different way. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we'll import the necessary packages to create a general prognostics model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from progpy import PrognosticsModel" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we'll define our model class. PrognosticsModels require defining [inputs](https://nasa.github.io/progpy/glossary.html#term-input), [states](https://nasa.github.io/progpy/glossary.html#term-state), [outputs](https://nasa.github.io/progpy/glossary.html#term-output), and [event](https://nasa.github.io/progpy/glossary.html#term-event) keys. As in the above example, the states include position (`x`) and velocity(`v`) of the object, the output is position (`x`), and the events are `falling` and `impact`. \n", + "\n", + "Note that we define this class as `ThrownObject_ST` to distinguish it as a state-transition model compared to the previous linear model class. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(PrognosticsModel):\n", + " \"\"\"\n", + " Model that simulates an object thrown into the air without air resistance\n", + " \"\"\"\n", + "\n", + " inputs = [] # no inputs, no way to control\n", + " states = [\n", + " 'x', # Position (m) \n", + " 'v' # Velocity (m/s)\n", + " ]\n", + " outputs = [ # Anything we can measure\n", + " 'x' # Position (m)\n", + " ]\n", + " events = [\n", + " 'falling', # Event- object is falling\n", + " 'impact' # Event- object has impacted ground\n", + " ]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we'll add some default parameter definitions. These values can be overwritten by passing parameters into the constructor. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(ThrownObject_ST):\n", + "\n", + " default_parameters = {\n", + " 'thrower_height': 1.83, # default height \n", + " 'throwing_speed': 40, # default speed\n", + " 'g': -9.81, # Acceleration due to gravity (m/s^2)\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All prognostics models require some specific class functions. We'll define those next. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we'll need to add the functionality to set the initial state of the system. There are two ways to provide the logic to initialize model state. \n", + "\n", + "1. Provide the initial state in `parameters['x0']`, or \n", + "2. Provide an `initialize` function \n", + "\n", + "The first method here is preferred. If `parameters['x0']` are defined, there is no need to explicitly define an initialize method, and these parameter values will be used as the initial state. \n", + "\n", + "However, there are some cases where the initial state is a function of the input (`u`) or output (`z`) (e.g. a use-case where the input is also a state). In this case, an explicitly defined `initialize` method is required. \n", + "\n", + "Here, we'll set our initial state by defining an `initialize` function. In the code below, note that the function can take arguments for both input `u` and output `z`, though these are optional. \n", + "\n", + "Note that for this example, defining initialize in this way is not necessary. We could have simply defined `parameters['x0']`. However, we choose to use this method for ease when using the `derived_params` feature, discussed in the next section. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(ThrownObject_ST):\n", + "\n", + " def initialize(self, u=None, z=None):\n", + " return self.StateContainer({\n", + " 'x': self.parameters['thrower_height'], # Thrown, so initial altitude is height of thrower\n", + " 'v': self.parameters['throwing_speed'] # Velocity at which the ball is thrown - this guy is a professional baseball pitcher\n", + " })" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, the PrognosticsModel class requires that we define how the state transitions throughout time. For continuous models, this is defined with the method `dx`, which calculates the first derivative of the state at a specific time. For discrete systems, this is defined with the method `next_state`, using the state transition equation for the system. When possible, it is recommended to use the continuous (`dx`) form, as some algorithms will only work on continuous models.\n", + "\n", + "Here, we use the equations for the derivatives of our system (i.e., the continuous form)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(ThrownObject_ST):\n", + "\n", + " def dx(self, x, u):\n", + " return self.StateContainer({\n", + " 'x': x['v'], \n", + " 'v': self.parameters['g']}) # Acceleration of gravity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we'll define the `output` method, which will calculate the output (i.e., measurable values) given the current state. Here, our output is position (`x`). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(ThrownObject_ST):\n", + " \n", + " def output(self, x):\n", + " return self.OutputContainer({'x': x['x']})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The next method we define is [`event_state`](https://nasa.github.io/progpy/glossary.html#term-event-state). As before, \n", + "`event_state` calculates the progress towards the events. Normalized to be between 0 and 1, 1 means there is no progress towards the event and 0 means the event has occurred. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(ThrownObject_ST):\n", + " \n", + " def event_state(self, x): \n", + " # Use speed and position to estimate maximum height\n", + " x_max = x['x'] + np.square(x['v'])/(-self.parameters['g']*2)\n", + " # 1 until falling begins\n", + " x_max = np.where(x['v'] > 0, x['x'], x_max)\n", + " return {\n", + " 'falling': max(x['v']/self.parameters['throwing_speed'],0), # Throwing speed is max speed\n", + " 'impact': max(x['x']/x_max,0) # 1 until falling begins, then it's fraction of height\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At this point, we have defined all necessary information for the PrognosticsModel to be complete. There are other methods that can additionally be defined (see the [PrognosticsModel](https://nasa.github.io/progpy/api_ref/prog_models/PrognosticModel.html) documentation for more information) to provide further configuration for new models. We'll highlight some of these in the following sections. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As an example of one of these, we additionally define a `threshold_met` equation. Note that this is optional. Leaving `threshold_met` empty will use the event state to define thresholds (threshold = event state == 0). However, this implementation is more efficient, so we include it. \n", + "\n", + "Here, we define `threshold_met` in the same way as our linear model example. `threshold_met` will return a _dict_ of values, one for each event. Threshold is met when all dictionary entries are __True__. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(ThrownObject_ST):\n", + "\n", + " def threshold_met(self, x):\n", + " return {\n", + " 'falling': x['v'] < 0, # Falling occurs when velocity becomes negative\n", + " 'impact': x['x'] <= 0 # Impact occurs when the object hits the ground, i.e. position is <= 0\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With that, we have created a new ThrownObject state-transition model. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's can test our model through simulation. First, we'll create an instance of the model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m_st = ThrownObject_ST()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll start by simulating to impact. We'll specify the `threshold_keys` to specifically indicate we are interested in impact. For more information on simulation, see 1. Simulation. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Simulate to impact\n", + "event = 'impact'\n", + "simulated_results = m_st.simulate_to_threshold(threshold_keys=[event], dt=0.005, save_freq=1, print = True)\n", + "\n", + "# Print result: \n", + "print('The object hit the ground in {} seconds'.format(round(simulated_results.times[-1],2)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To summarize this section, we have illustrated how to construct new physics-based models by subclassing from [progpy.PrognosticsModel](https://nasa.github.io/progpy/api_ref/prog_models/PrognosticModel.html#prog_models.PrognosticsModel). Some elements (e.g. inputs, states, outputs, events keys; methods for initialization, dx/next_state, output, and event_state) are required. Models can be additionally configured with additional methods and parameters.\n", + "\n", + "Note that in this example, we defined each part one piece at a time, recursively subclassing the partially defined class. This was done to illustrate the parts of the model. In reality, all methods and properties would be defined together in a single class definition. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Derived Parameters" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the previous section, we constructed a new model from scratch by subclassing from [progpy.PrognosticsModel](https://nasa.github.io/progpy/api_ref/prog_models/PrognosticModel.html#prog_models.PrognosticsModel) and specifying all of the necessary model components. An additional optional feature of PrognosticsModels is derived parameters, illustrated below. \n", + "\n", + "A derived parameter is a parameter (see parameter section in 1. Simulation) that is a function of another parameter. For example, in the case of a thrown object, one could assume that throwing speed is a function of thrower height, with taller throwing height resulting in faster throwing speeds. In the electrochemistry battery model (see 3. Included Models), there are parameters for the maximum and minimum charge at the surface and bulk, and these are dependent on the capacity of the battery (i.e. another parameter, qMax). When such derived parameters exist, they must be updated whenever the parameters they depend on are updated. In PrognosticsModels, this is achieved with the `derived_params` feature. \n", + "\n", + "This feature can also be used to cache combinations of parameters that are used frequently in state transition or other model methods. Creating lumped parameters using `derived_params` causes them to be calculated once when configuring, instead of each time step in simulation or prediction. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this example, we will use the `ThrownObject_ST` model created in the previous section. We will extend this model to include a derived parameter, namely `throwing_speed` will be dependent on `thrower_height`.\n", + "\n", + "To implement this, we must first define a function for the relationship between the two parameters. We'll assume that `throwing_speed` is a linear function of `thrower_height`. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def update_thrown_speed(params):\n", + " return {\n", + " 'throwing_speed': params['thrower_height'] * 21.85\n", + " } \n", + " # Note: one or more parameters can be changed in these functions, whatever parameters are changed are returned in the dictionary" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we'll define the parameter callbacks, so that `throwing_speed` is updated appropriately any time that `thrower_height` changes. The following effectively tells the derived callbacks feature to call the `update_thrown_speed` function whenever the `thrower_height` changes. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_ST(ThrownObject_ST):\n", + "\n", + " param_callbacks = {\n", + " 'thrower_height': [update_thrown_speed]\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also have more than one function be called when a single parameter is changed. You would do this by adding the additional callbacks to the list (e.g., 'thrower_height': [update_thrown_speed, other_fcn])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have now added the capability for `throwing_speed` to be a derived parameter. Let's try it out. First, we'll create an instance of our class and print out the default parameters. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "obj = ThrownObject_ST()\n", + "print(\"Default Settings:\\n\\tthrower_height: {}\\n\\tthrowing_speed: {}\".format(obj.parameters['thrower_height'], obj.parameters['throwing_speed']))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, let's change the thrower height. If our derived parameters work correctly, the thrower speed should change accordingly. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "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']))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As we can see, when the thrower height was changed, the throwing speed was re-calculated too. \n", + "\n", + "In this example, we illustrated how to use the `derived_params` feature, which allows a parameter to be a function of another parameter. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Direct Models" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the previous sections, we illustrated how to create and use state-transition models, or models that use state transition differential equations to propagate the state forward. In this example, we'll explore another type of model implemented within ProgPy - Direct Models. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Direct models estimate the time of event directly from the system state and future load, rather than through state transitions. This approach is particularly useful for physics-based models where the differential equations of state transitions can be explicitly solved, or for data-driven models that map sensor data directly to the time of an event. When applicable, using a direct model approach provides a more efficient way to estimate the time of an event, especially for events that occur later in the simulation. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To illustrate this concept, we will extend the state-transition model, `ThrownObject_ST`, defined above, to create a new model class, `DirectThrownObject`. The dynamics of a thrown object lend easily to a direct model, since we can solve the differential equations explicitly to estimate the time at which the events occur. \n", + "\n", + "Recall that our physical system is described by the following differential equations: \n", + "\\begin{align*}\n", + "\\frac{dx}{dt} &= v \\\\ \\\\\n", + "\\frac{dv}{dt} &= -g \n", + "\\end{align*}\n", + "\n", + "which can be solved explicity, given initial position $x_0$ and initial velocity $v_0$, to get:\n", + "\\begin{align*}\n", + "x(t) &= -\\frac{1}{2} gt^2 + v_0 t + x_0 \\\\ \\\\ \n", + "v(t) &= -gt + v_0\n", + "\\end{align*}\n", + "\n", + "Setting these equations to 0 and solving for time, we get the time at which the object hits the ground and begins falling, respectively. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To construct our direct model, we'll extend the `ThrownObject_ST` model to additionally include the method [time_to_event](https://nasa.github.io/progpy/api_ref/prog_models/PrognosticModel.html#prog_models.PrognosticsModel.time_of_event). This method will calculate the time at which each event occurs (i.e., time when the event threshold is met), based on the equations above. `time_of_event` must be implemented by any direct model. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class DirectThrownObject(ThrownObject_ST):\n", + " \n", + " def time_of_event(self, x, *args, **kwargs):\n", + " # calculate time when object hits ground given x['x'] and x['v']\n", + " # 0 = x0 + v0*t - 0.5*g*t^2\n", + " g = self.parameters['g']\n", + " t_impact = -(x['v'] + np.sqrt(x['v']*x['v'] - 2*g*x['x']))/g\n", + "\n", + " # 0 = v0 - g*t\n", + " t_falling = -x['v']/g\n", + " \n", + " return {'falling': t_falling, 'impact': t_impact}\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With this, our direct model is created. Note that adding `*args` and `**kwargs` is optional. Having these arguments makes the function interchangeable with other models which may have arguments or keyword arguments. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's test out this capability. To do so, we'll use the `time` package to compare the direct model to our original timeseries model. \n", + "\n", + "Let's start by creating an instance of our timeseries model, calculating the time of event, and timing this computation. Note that for a state transition model, `time_of_event` still returns the time at which `threshold_met` returns true for each event, but this is calculated by simulating to threshold." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time \n", + "\n", + "m_timeseries = ThrownObject_ST()\n", + "x = m_timeseries.initialize()\n", + "print(m_timeseries.__class__.__name__, \"(Direct Model)\" if m_timeseries.is_direct else \"(Timeseries Model)\")\n", + "tic = time.perf_counter()\n", + "print('Time of event: ', m_timeseries.time_of_event(x, dt = 0.05))\n", + "toc = time.perf_counter()\n", + "print(f'execution: {(toc-tic)*1000:0.4f} milliseconds')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's do the same using our direct model implementation. In this case, when `time_to_event` is called, the event time will be estimated directly from the state, instead of through simulation to threshold. \n", + "\n", + "Note that a limitation of a direct model is that you cannot get intermediate states (i.e., save_pts or save_freq) since the time of event is calculated directly. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m_direct = DirectThrownObject()\n", + "x = m_direct.initialize() # Using Initial state\n", + "# Now instead of simulating to threshold, we can estimate it directly from the state, like so\n", + "print('\\n', m_direct.__class__.__name__, \"(Direct Model)\" if m_direct.is_direct else \"(Timeseries Model)\")\n", + "tic = time.perf_counter()\n", + "print('Time of event: ', m_direct.time_of_event(x))\n", + "toc = time.perf_counter()\n", + "print(f'execution: {(toc-tic)*1000:0.4f} milliseconds')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that execution is significantly faster for the direct model. Furthermore, the result is actually more accurate, since it's not limited by the timestep (see dt section in 1. Simulation). These observations will be even more pronounced for events that occur later in the simulation. \n", + "\n", + "It's important to note that this is a very simple example, as there are no inputs. For models with inputs, future loading must be provided to `time_of_event` (see the Future Loading section in 1. Simulation). In these cases, most direct models will encode or discretize the future loading profile to use it in a direct estimation of time of event." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the example provided, we have illustrated how to use a direct model. Direct models are a powerful tool for estimating the time of an event directly from the system state. By avoiding the process of state transitions, direct models can provide more efficient event time estimates. Additionally, the direct model approach is not limited to physics-based models. It can also be applied to data-driven models that can map sensor data directly to the time of an event. \n", + "\n", + "In conclusion, direct models offer an efficient and versatile approach for prognostics modeling, enabling faster and more direct estimations of event times. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Matrix Data Access Feature" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the above models, we have used dictionaries to represent the states. For example, in the implementation of `ThrownObject_ST` above, see how `dx` is defined with a StateContainer dictionary. While all models can be constructed using dictionaries in this way, some dynamical systems allow for the state of the system to be represented with a matrix. For such use-cases, ProgPy has an advanced *matrix data access feature* that provides a more efficient way to define these models.\n", + "\n", + "In ProgPy's implementation, the provided model.StateContainer, InputContainer, and OutputContainers can be treated as dictionaries but use an underlying matrix. This is important for some applications like surrogate and machine-learned models where the state is represented by a tensor. ProgPy's *matrix data access feature* allows the matrices to be used directly. Simulation functions propagate the state using the matrix form, preventing the inefficiency of having to convert to and from dictionaries. Additionally, this implementation is faster than recreating the StateContainer each time, especially when updating inplace." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example, we'll illustrate how to use the matrix data access feature. We'll continue with our ThrownObject system, and create a model to simulate this using matrix notation (instead of dictionary notation as in the standard model, seen above in `ThrownObject_ST`). The implementation of the model is comparable to a standard model, except that it uses matrix operations within each function, as seen below. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, the necessary imports." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from progpy import PrognosticsModel" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To use the matrix data access feature, we'll subclass from our state-transition model defined above, `ThrownObject_ST`. Our new model will therefore inherit the default parameters and methods for initialization, output, threshold met, and event state. \n", + "\n", + "To use the matrix data access feature, we'll use matrices to define how the state transitions. Since we are working with a discrete version of the system now, we'll define the `next_state` method, and this will override the `dx` method in the parent class. \n", + "\n", + "In the following, we will use the matrix version for each variable, accessed with `.matrix`. We implement this within `next_state`, but this feature can also be used in other functions. Here, both `x.matrix` and `u.matrix` are column vectors, and `u.matrix` is in the same order as model.inputs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class ThrownObject_MM(ThrownObject_ST):\n", + "\n", + " def next_state(self, x, u, dt):\n", + "\n", + " A = np.array([[0, 1], [0, 0]]) # State transition matrix\n", + " B = np.array([[0], [self.parameters['g']]]) # Acceleration due to gravity\n", + " x.matrix += (np.matmul(A, x.matrix) + B) * dt\n", + "\n", + " return x" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our model is now specified. Let's try simulating with it." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we'll create an instance of the model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m_matrix = ThrownObject_MM()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, let's simulate to threshold. We'll also time the simulation so we can compare with the non-matrix state-transition model below. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time \n", + "\n", + "tic_matrix = time.perf_counter()\n", + "# Simulate to threshold \n", + "m_matrix.simulate_to_threshold(\n", + " print = True, \n", + " threshold_keys = 'impact', \n", + " dt = 0.1, \n", + " save_freq = 1)\n", + "toc_matrix = time.perf_counter()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our matrix notation was successful in simulating the thrown object's behavior throughout time. \n", + "\n", + "Finally, let's simulate the non-matrix version to compare computation time. \n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tic_st = time.perf_counter()\n", + "m_st.simulate_to_threshold(\n", + " print = True, \n", + " threshold_keys = 'impact', \n", + " dt = 0.1, \n", + " save_freq = 1)\n", + "toc_st = time.perf_counter()\n", + "\n", + "print(f'Matrix execution: {(toc_matrix-tic_matrix)*1000:0.4f} milliseconds')\n", + "print(f'Non-matrix execution: {(toc_st-tic_st)*1000:0.4f} milliseconds')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As we can see, for this system, using the matrix data access feature is computationally faster than a standard state-transition matrix that uses dictionaries. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As illustrated here, the matrix data access feature is an advanced capability that represents the state of a system using matrices. This can provide efficiency for use-cases where the state is easily represented by a tensor and operations are defined by matrices." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## State Limits" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In real-world physical systems, there are often constraints on what values the states can take. For example, in the case of a thrown object, if we define our reference frame with the ground at a position of $x=0$, then the position of the object should only be greater than or equal to 0, and should never take on negative values. In ProgPy, we can enforce constraints on the range of each state for a state-transition model using the [state limits](https://nasa.github.io/progpy/prog_models_guide.html#state-limits) attribute. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To illustrate the use of `state_limits`, we'll use our thrown object model `ThrownObject_ST`, created in an above section. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m_limits = ThrownObject_ST()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before adding state limits, let's take a look at the standard model without state limits. We'll consider the event of `impact`, and simulate the object to threshold." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "event = 'impact'\n", + "simulated_results = m_limits.simulate_to_threshold(threshold_keys=[event], dt=0.005, save_freq=1)\n", + "\n", + "print('Example: No State Limits')\n", + "for i, state in enumerate(simulated_results.states):\n", + " print(f'State {i}: {state}')\n", + "print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that at the end of the simulation, the object's position (`x`) is negative. This doesn't make sense physically, since the object cannot fall below ground level (at $x=0$)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To avoid this, and keep the state in a realistic range, we can change the `state_limits` attribute of the model. The `state_limits` attribute is a dictionary that contains the state limits for each state. The keys of the dictionary are the state names, and the values are tuples that contain the lower and upper limits of the state. \n", + "\n", + "In our Thrown Object model, our states are position, which can range from 0 to infinity, and velocity, which we'll limit to not exceed the speed of light." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import inf\n", + "from math import inf\n", + "\n", + "m_limits.state_limits = {\n", + " # object position may not go below ground height\n", + " 'x': (0, inf),\n", + "\n", + " # object velocity may not exceed the speed of light\n", + " 'v': (-299792458, 299792458)\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we've specified the ranges for our state values, let's try simulating again. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "event = 'impact'\n", + "simulated_results = m_limits.simulate_to_threshold(threshold_keys=[event], dt=0.005, save_freq=1)\n", + "\n", + "print('Example: With State Limits')\n", + "for i, state in enumerate(simulated_results.states):\n", + " print(f'State {i}: {state}')\n", + "print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that now the position (`x`) becomes 0 but never reaches a negative value. This is because we have defined a state limit for the `x` state that prevents it from going below 0. Also note that a warning is provided to notify the user that a state value was limited. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's try a more complicated example. This time, we'll try setting the initial position value to be a number outside of its bounds. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x0 = m_limits.initialize(u = {}, z = {})\n", + "x0['x'] = -1 # Initial position value set to an unrealistic value of -1\n", + "\n", + "simulated_results = m_limits.simulate_to_threshold(threshold_keys=[event], dt=0.005, save_freq=1, x = x0)\n", + "\n", + "# Print states\n", + "print('Example 2: With -1 as initial x value')\n", + "for i, state in enumerate(simulated_results.states):\n", + " print('State ', i, ': ', state)\n", + "print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that the simulation stops after just two iterations. In this case, the initial position value is outside the state limit. On the first iteration, the position value is therefore adjusted to be within the appropriate range of 0 to $\\infty$. Since we are simulating to impact, which is defined as when position is 0, the threshold is immediately satisfied and the simulation stops. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, note that limits can also be applied manually using the `apply_limits` function. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = {'x': -5, 'v': 3e8} # Too fast and below the ground\n", + "print('\\t Pre-limit: {}'.format(x))\n", + "\n", + "x = m_limits.apply_limits(x)\n", + "print('\\t Post-limit: {}'.format(x))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In conclusion, setting appropriate [state limits](https://nasa.github.io/progpy/prog_models_guide.html#state-limits) is crucial in creating realistic and accurate state-transition models. It ensures that the model's behavior stays within the constraints of the physical system. The limits should be set based on the physical or practical constraints of the system being modeled. \n", + "\n", + "As a final note, state limits are especially important for state estimation (to be discussed in the State Estimation section), as it will force the state estimator to only consider states that are possible or feasible. State estimation will be described in more detail in section 08. State Estimation. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Custom Events" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the examples above, we have focused on the simple event of a thrown object hitting the ground or reaching `impact`. In this section, we highlight additional uses of ProgPy's generalizable concept of `events`. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The term [events](https://nasa.github.io/progpy/prog_models_guide.html#events) is used to describe something to be predicted. Generally in the PHM community, these are referred to as End of Life (EOL). However, they can be much more. \n", + "\n", + "In ProgPy, events can be anything that needs to be predicted. Systems will often have multiple failure modes, and each of these modes can be represented by a separate event. Additionally, events can also be used to predict other events of interest other than failure, such as special system states or warning thresholds. Thus, `events` in ProgPy can represent End of Life (EOL), End of Mission (EOM), warning thresholds, or any Event of Interest (EOI). " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are a few components of the model that must be specified in order to define events:\n", + "\n", + "1. The `events` property defines the expected events \n", + "\n", + "2. The `threshold_met` method defines the conditions under which an event occurs \n", + "\n", + "3. The `event_state` method returns an estimate of progress towards the threshold \n", + "\n", + "Note that because of the interconnected relationship between `threshold_met` and `event_state`, it is only required to define one of these. However, it is generally beneficial to specify both. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To illustrate this concept, we will use the `BatteryElectroChemEOD` model (see section 03. Included Models). In the standard implementation of this model, the defined event is `EOD` or End of Discharge. This occurs when the voltage drops below a pre-defined threshold value. The State-of-Charge (SOC) of the battery is the event state for the EOD event. Recall that event states (and therefore SOC) vary between 0 and 1, where 1 is healthy and 0 signifies the event has occurred. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Suppose we have the requirement that our battery must not fall below 5% State-of-Charge. This would correspond to an `EOD` event state of 0.05. Additionally, let's add events for two warning thresholds, a $\\text{\\textcolor{yellow}{yellow}}$ threshold at 15% SOC and a $\\text{\\textcolor{red}{red}}$ threshold at 10% SOC. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To define the model, we'll start with the necessary imports." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "from progpy.loading import Piecewise\n", + "from progpy.models import BatteryElectroChemEOD" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, let's define our threshold values. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "YELLOW_THRESH = 0.15 # 15% SOC\n", + "RED_THRESH = 0.1 # 10% SOC\n", + "THRESHOLD = 0.05 # 5% SOC" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we'll create our model by subclassing from the `BatteryElectroChemEOD` model. First, we'll re-define `events` to include three new events for our two warnings and new threshold value, as well as the event `EOD` from the parent class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class BattNewEvent(BatteryElectroChemEOD):\n", + " events = BatteryElectroChemEOD.events + ['EOD_warn_yellow', 'EOD_warn_red', 'EOD_requirement_threshold']\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we'll override the `event_state` method to additionally include calculations for progress towards each of our new events. We'll add yellow, red, and failure states by scaling the EOD state. We scale so that the threshold SOC is 0 at their associated events, while SOC of 1 is still 1. For example, for yellow, we want `EOD_warn_yellow` to be 1 when SOC is 1, and 0 when SOC is 0.15 or lower. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class BattNewEvent(BattNewEvent):\n", + " \n", + " def event_state(self, state):\n", + " # Get event state from parent\n", + " event_state = super().event_state(state)\n", + "\n", + " # Add yellow, red, and failure states by scaling EOD state\n", + " event_state['EOD_warn_yellow'] = (event_state['EOD']-YELLOW_THRESH)/(1-YELLOW_THRESH) \n", + " event_state['EOD_warn_red'] = (event_state['EOD']-RED_THRESH)/(1-RED_THRESH)\n", + " event_state['EOD_requirement_threshold'] = (event_state['EOD']-THRESHOLD)/(1-THRESHOLD)\n", + "\n", + " # Return\n", + " return event_state" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we'll override the `threshold_met` method to define when each event occurs. Based on the scaling in `event_state` each event is reached when the corresponding `event_state` value is less than or equal to 0. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class BattNewEvent(BattNewEvent):\n", + " def threshold_met(self, x):\n", + " # Get threshold met from parent\n", + " t_met = super().threshold_met(x)\n", + "\n", + " # Add yell and red states from event_state\n", + " event_state = self.event_state(x)\n", + " t_met['EOD_warn_yellow'] = event_state['EOD_warn_yellow'] <= 0\n", + " t_met['EOD_warn_red'] = event_state['EOD_warn_red'] <= 0\n", + " t_met['EOD_requirement_threshold'] = event_state['EOD_requirement_threshold'] <= 0\n", + "\n", + " return t_met" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With this, we have defined the three key model components for defining new events. \n", + "\n", + "Let's test out the model. First, create an instance of it. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = BattNewEvent()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Recall that the battery model takes input of current. We will use a piecewise loading scheme (see 01. Simulation)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Variable (piecewise) future loading scheme\n", + "future_loading = Piecewise(\n", + " m.InputContainer,\n", + " [600, 900, 1800, 3000],\n", + " {'i': [2, 1, 4, 2, 3]})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can simulate to threshold and plot the results. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "simulated_results = m.simulate_to_threshold(future_loading, threshold_keys=['EOD'], print = True)\n", + "\n", + "simulated_results.event_states.plot()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, we can see the SOC plotted for the different events throughout time. The yellow warning (15% SOC) reaches threshold first, followed by the red warning (10% SOC), new EOD threshold (5% SOC), and finally the original EOD value. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this section, we have illustrated how to define custom [events](https://nasa.github.io/progpy/prog_models_guide.html#events) for prognostics models. Events can be used to define anything that a user is interested in predicting, including common values like Remaining Useful Life (RUL) and End of Discharge (EOD), as well as other values like special intermediate states or warning thresholds. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Serialization " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "ProgPy includes a feature to serialize models, which we highlight in this section. \n", + "\n", + "Model serialization has a variety of purposes. For example, serialization allows us to save a specific model or model configuration to a file to be loaded later, or can aid us in sending a model to another machine over a network connection. Some users maintain a directory or repository of configured models representing specific systems in their stock.\n", + "\n", + "In this section, we'll show how to serialize and deserialize model objects using `pickle` and `JSON` methods. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we'll import the necessary modules." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import pickle\n", + "import numpy as np\n", + "from progpy.models import BatteryElectroChemEOD\n", + "from progpy.loading import Piecewise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this example, we'll use the BatteryElectroChemEOD model. We'll start by creating a model object. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "batt = BatteryElectroChemEOD()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we'll serialize the model in two different ways using 1) `pickle` and 2) `JSON`. Then, we'll plot the results from simulating the deserialized models to show equivalence of the methods. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To save using the `pickle` package, we'll serialize the model using the `dump` method. Once saved, we can then deserialize using the `load` method. In practice, deserializing will likely occur in a different file or in a later use-case, but here we deserialize to show equivalence of the saved model. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pickle.dump(batt, open('save_pkl.pkl', 'wb')) # Serialize model\n", + "load_pkl = pickle.load(open('save_pkl.pkl', 'rb')) # Deserialize model " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we'll serialize using the `to_json` method. We deserialize by calling the model directly with the serialized result using the `from_json` method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "save_json = batt.to_json() # Serialize model\n", + "json_1 = BatteryElectroChemEOD.from_json(save_json) # Deserialize model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the serialized result can also be saved to a text file and uploaded for later use. We demonstrate this below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "txtFile = open(\"save_json.txt\", \"w\")\n", + "txtFile.write(save_json)\n", + "txtFile.close()\n", + "\n", + "with open('save_json.txt') as infile: \n", + " load_json = infile.read()\n", + "\n", + "json_2 = BatteryElectroChemEOD.from_json(load_json)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have now serialized and deserialized the model using `pickle` and `JSON` methods. Let's compare the resulting models. To do so, we'll use ProgPy's [simulation](https://nasa.github.io/progpy/prog_models_guide.html#simulation) to simulate the model to threshold and compare the results. \n", + "\n", + "First, we'll need to define our [future loading profile](https://nasa.github.io/progpy/prog_models_guide.html#future-loading) using the PiecewiseLoad class. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Variable (piecewise) future loading scheme\n", + "future_loading = Piecewise(\n", + " batt.InputContainer,\n", + " [600, 1000, 1500, 3000],\n", + " {'i': [3, 2, 1.5, 4]})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, let's simulate each model to threshold using the `simulate_to_threshold` method. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Original model \n", + "results_orig = batt.simulate_to_threshold(future_loading, save_freq = 1)\n", + "# Pickled version \n", + "results_pkl = load_pkl.simulate_to_threshold(future_loading, save_freq = 1)\n", + "# JSON versions\n", + "results_json_1 = json_1.simulate_to_threshold(future_loading, save_freq = 1)\n", + "results_json_2 = json_2.simulate_to_threshold(future_loading, save_freq = 1)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's plot the results for comparison." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "voltage_orig = [results_orig.outputs[iter]['v'] for iter in range(len(results_orig.times))]\n", + "voltage_pkl = [results_pkl.outputs[iter]['v'] for iter in range(len(results_pkl.times))]\n", + "voltage_json_1 = [results_json_1.outputs[iter]['v'] for iter in range(len(results_json_1.times))]\n", + "voltage_json_2 = [results_json_2.outputs[iter]['v'] for iter in range(len(results_json_2.times))]\n", + "\n", + "plt.plot(results_orig.times,voltage_orig,'-b',label='Original surrogate') \n", + "plt.plot(results_pkl.times,voltage_pkl,'--r',label='Pickled serialized surrogate') \n", + "plt.plot(results_json_1.times,voltage_json_1,'-.g',label='First JSON serialized surrogate') \n", + "plt.plot(results_json_2.times, voltage_json_2, '--y', label='Second JSON serialized surrogate')\n", + "plt.legend()\n", + "plt.xlabel('Time (sec)')\n", + "plt.ylabel('Voltage (volts)')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All of the voltage curves overlap, showing that the different serialization methods produce the same results. \n", + "\n", + "Additionally, we can compare the output arrays directly, to ensure equivalence. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "# Check if the arrays are the same\n", + "are_arrays_same = np.array_equal(voltage_orig, voltage_pkl) and \\\n", + " np.array_equal(voltage_orig, voltage_json_1) and \\\n", + " np.array_equal(voltage_orig, voltage_json_2)\n", + "\n", + "print(f\"The simulated results from the original and serialized models are {'identical. This means that our serialization works!' if are_arrays_same else 'not identical. This means that our serialization does not work.'}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To conclude, we have shown how to serialize models in ProgPy using both `pickle` and `JSON` methods. Understanding how to serialize and deserialize models can be a powerful tool for prognostics developers. It enables the saving of models to a disk and the re-loading of these models back into memory at a later time. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conclusions" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In these examples, we have described how to create new physics-based models. We have illustrated how to construct a generic physics-based model, as well as highlighted some specific types of models including linear models and direct models. We highlighted the matrix data access feature for using matrix operations more efficiently. Additionally, we discussed a few important components of any prognostics model including derived parameters, state limits, and events. \n", + "\n", + "With these tools, users are well-equipped to build their own prognostics models for their specific physics-based use-cases. In the next example, we'll discuss how to create data-driven models." + ] + } + ], + "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.9.6" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}