From faf6bb2311cab533e0fb1af37556033d87176cef Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Wed, 23 Aug 2023 12:41:15 -0700 Subject: [PATCH 1/7] Add loading section --- examples/01_Simulation.ipynb | 521 +++++++++++++++++++++++++++++++++++ 1 file changed, 521 insertions(+) diff --git a/examples/01_Simulation.ipynb b/examples/01_Simulation.ipynb index 11688e5..03ba6d6 100644 --- a/examples/01_Simulation.ipynb +++ b/examples/01_Simulation.ipynb @@ -342,6 +342,527 @@ "## Future Loading" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The previous examples feature a simple ThrownObject model, which does not have any inputs. Unlike ThrownObject, most prognostics models have some sort of [input](https://nasa.github.io/progpy/glossary.html#term-input). The input is some sort of control or loading applied to the system being modeled. In this section we will describe how to simulate a model which features an input.\n", + "\n", + "In this example we will be using the BatteryCircuit model from the models subpackage (see 3. Included Models). This is a simple battery discharge model where the battery is represented by an equivilant circuit.\n", + "\n", + "Like the past examples, we start by importing and creating the model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from prog_models.models import BatteryCircuit\n", + "m = BatteryCircuit()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can see the batteries inputs, states, and outputs (described above) by accessing these attributes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('outputs:', m.outputs)\n", + "print('inputs:', m.inputs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Consulting the [model documentation](https://nasa.github.io/progpy/api_ref/prog_models/IncludedModels.html#prog_models.models.BatteryCircuit), I see that the outputs (i.e., measurable values) of the model are temperature (`t`) and voltage (`v`). The model's input is the current (`i`) drawn from the model.\n", + "\n", + "If we try to simulate as we do above, it wouldn't work because the battery discharge is a function of the current (`i`) drawn from the battery. Simulation for a model like this requires that we define the future load. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Piecewise load\n", + "\n", + "For the first example, we define a piecewise loading profile using the `progpy.loading.Piecewise` class. This is one of the most common loading profiles. First we import the class from the loading subpackage" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from progpy.loading import Piecewise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we define a loading profile. Piecewise loader takes 3 arguments: 1. the model InputContainer, 2. times and 3. loads, each of these are explained in more detail below.\n", + "\n", + "The model input container is a class for representing the input for a model. It's a class attribute for every model, and is specific to that model. It can be found at m.InputContainer\n", + "\n", + "Times are the ending time for each load. For example, if times were [5, 7, 10], then the first load would apply until t=5, then the second load would apply for 2 seconds, following by the third load for 5 more seconds. \n", + "\n", + "Loads are a dictionary of arrays, where the keys of the dictionary are the inputs to the model (for a battery, just current `i`), and the values in the array are the value at each time in times. If the loads array is one longer than times, then the last value is the \"default load\", i.e., the load that will be applied after the last time has passed.\n", + "\n", + "For example, we might define this load profile for our battery." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "loading = Piecewise(\n", + " InputContainer=m.InputContainer,\n", + " times=[600, 900, 1800, 3000],\n", + " values={'i': [2, 1, 4, 2, 3]})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this case, the current drawn (`i`) is 2 amps until t is 600 seconds, then it is 1 for the next 300 seconds (until 900 seconds), etc. The \"default load\" is 3, meaning that after the last time has passed (3000 seconds) a current of 3 will be drawn. \n", + "\n", + "Now that we have this load profile, lets run a simulation with our model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results = m.simulate_to_threshold(loading, save_freq=100)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's take a look at the inputs to the model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = results.inputs.plot(ylabel=\"Current Draw (amps)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "See above that the load profile is piecewise, matching the profile we defined above.\n", + "\n", + "Plotting the outputs, you can see jumps in the voltage levels as the current changes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = results.outputs.plot(compact=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we simulated to threshold, loading the system using a simple piecewise load profile. This is the most common load profile and will probably work for most cases " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Moving Average" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Another common loading scheme is the moving-average load. This loading scheme assumes that the load will continue like it's seen in the past. This is useful when you dont know the exact load, but you expect it to be consistant.\n", + "\n", + "Like with Piecewise loading, the first step it to import the loading class. In this case, `progpy.loading.MovingAverage`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from progpy.loading import MovingAverage" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we create the moving average loading object, passing in the InputContainer" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "loading = MovingAverage(m.InputContainer)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The moving average load estimator requires an additional step, sending the observed load. This is done using the add_load method. Let's load it up with some observed current draws. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "measured_loads = [4, 4.5, 4.0, 4, 2.1, 1.8, 1.99, 2.0, 2.01, 1.89, 1.92, 2.01, 2.1, 2.2]\n", + " \n", + "for load in measured_loads:\n", + " loading.add_load({'i': load})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In practice the add_load method should be called whenever there's new input (i.e., load) information. The MovingAverage load estimator averages over a window of elemetns, configurable at construction using the window argument (e.g., MovingAverage(m.InputContainer, window=12))\n", + "\n", + "Now the configured load estimator can be used in simulation. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results = m.simulate_to_threshold(loading, save_freq=100)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's take a look at the resulting input current." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = results.inputs.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the loading is a constant around 2, this is because the larger loads (~4 amps) are outside of the averaging window. Here are the resulting outputs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = results.outputs.plot(compact=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the voltage and temperature curves are much cleaner. They dont have the jumps present in the piecewise loading example. This is due to the constant loading.\n", + "\n", + "In this example we simulated to threshold, loading the system using a constant load profile caluclated using the moving average load estimator. This load estimator needs to be updated with the add_load method whenever new loading data is available. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Gaussian Noise in Loading" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Typically, users have an idea of what loading will look like, but there is some uncertainty. Future load extimates are hardly ever known exactly. This is where load wrappers like the `progpy.loading.GaussianNoiseLoadWrapper` come into play. The GaussianNoiseLoadWrapper wraps around another load profile, adding a random amount of noise, sampled from a gaussian distribution, at each step. In simulation this will show some variability, but this becomes more important in prediction (see 9. Prediction).\n", + "\n", + "In this example we will repeat the Piecewise load example, this time using the GaussianNoiseLoadWrapper to represent our uncertainty in our future load estimate. \n", + "\n", + "First we will import the necessary classes and construct the Piecewise load estimation just as in the previous example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from progpy.loading import Piecewise, GaussianNoiseLoadWrapper\n", + "loading = Piecewise(\n", + " InputContainer=m.InputContainer,\n", + " times=[600, 900, 1800, 3000],\n", + " values={'i': [2, 1, 4, 2, 3]})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we will wrap this loading object in our Gaussian noise load wrapper" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "loading = GaussianNoiseLoadWrapper(loading, 0.2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this case we're adding guassian noise with a standard deviation of 0.2 to the result of the previous load estimator.\n", + "\n", + "Now let's simulate and look at the input profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results = m.simulate_to_threshold(loading, save_freq=100)\n", + "fig = results.inputs.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note the loading profile follows the piecewise shape, but with noise. If you run it again, you would get a slightly different result." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results = m.simulate_to_threshold(loading, save_freq=100)\n", + "fig = results.inputs.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here are the corresponding outputs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = results.outputs.plot(compact=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the noise in input can be seen in the resulting output plots.\n", + "\n", + "In this section we introduced the concept of NoiseWrappers and how they are used to represent uncertainty in future loading. This concept is expecially important when used with prediction (see 9. Prediction). A GuassianNoiseLoadWrapper was used with a Piecewise loading profile to demonstrate it, but NoiseWrappers can be applied to any loading object or function, including the advanced profiles introduced in the next section." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Custom load profiles" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For most applications, the standard load estimation classes can be used to represent a users uexpectation of future loading. However, there are some cases where load is some complex combination of time and state that cannot be represented by these classes. This section briefly describes a few of these cases. \n", + "\n", + "The first example is similar to the last one, in that there is gaussian noise added to some underlying load profile. In this case the magnitude of noise increases linerally with time. This is an important example, as it allows us to represent a case where loading further out in time has more uncertainty (i.e., is less well known). This is common for many prognostic usecases.\n", + "\n", + "Custom load profiles can be represented either as a function (t, x=None) -> u, where t is time, x is state, and u is input. or as a class which implements the __call__ method with the same profile as the function.\n", + "\n", + "In this case we will use the first method (i.e., the function). We will define a function that will use a defined slope (derivitive of standard deviation with time)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from numpy.random import normal\n", + "base_load = 2 # Base load (amps)\n", + "std_slope = 1e-4 # Derivative of standard deviation with time\n", + "def loading(t, x=None):\n", + " std = std_slope * t\n", + " return m.InputContainer({'i': normal(base_load, std)})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the above code is specifically for a battery, but it could be generalized to any system.\n", + "\n", + "Now let's simulate and look at the input profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results = m.simulate_to_threshold(loading, save_freq=100)\n", + "fig = results.inputs.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note how the noise in the input signal increases with time. Since this is a random process, if you were to run this again you would get a different result.\n", + "\n", + "Here is the corresponding output. Note you can see the effects of the increasingly erradic input in the voltage curve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = results.outputs.plot(compact=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the final example we will define a loading profile that considers state. In this example, we're simulating a scenario where loads are removed (i.e., turned off) when discharge event state (i.e., SOC) reaches 0.25. This emulates a \"low power mode\" often employed in battery-powered electronics.\n", + "\n", + "For simplicity the underlying load will be constant, but this same approach could be applied to more complex profiles, and noise can be added on top using a NoiseWrapper." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "normal_load = m.InputContainer({'i': 2.7})\n", + "low_power_load = m.InputContainer({'i': 1.9})\n", + "\n", + "def loading(t, x=None):\n", + " if x is not None:\n", + " # State is provided\n", + " soc = m.event_state(x)['EOD']\n", + " return normal_load if soc > 0.25 else low_power_load\n", + " return normal_load" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the above example checks if x is not None. For some models, for the first timestep, state may be None (because state hasn't yet been calculated).\n", + "\n", + "Now let's use this in simulation and take a look at the loading profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results = m.simulate_to_threshold(loading, save_freq=100)\n", + "fig = results.inputs.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, as expected, load is at normal level for most of the time, then falls to low power mode towards the end of discharge.\n", + "\n", + "Let's look at the corresponding outputs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = results.outputs.plot(compact=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice the jump in voltage at the point where the load changed. Low power mode extended the life of the battery.\n", + "\n", + "In this section we show how to make custom loading profiles. Most applications can use the standard load classes, but some may require creating complex custom load profiles using this feature." + ] + }, { "cell_type": "markdown", "metadata": {}, From 9813df3c0206dd79978014bd16078ce465aebd5e Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Wed, 23 Aug 2023 12:48:21 -0700 Subject: [PATCH 2/7] Remove old example --- examples/future_loading.py | 145 ------------------------------------- 1 file changed, 145 deletions(-) delete mode 100644 examples/future_loading.py diff --git a/examples/future_loading.py b/examples/future_loading.py deleted file mode 100644 index 78f72eb..0000000 --- a/examples/future_loading.py +++ /dev/null @@ -1,145 +0,0 @@ -# Copyright © 2021 United States Government as represented by the Administrator of the -# National Aeronautics and Space Administration. All Rights Reserved. - -""" -Example demonstrating ways to use future loading. -""" - -import matplotlib.pyplot as plt -from numpy.random import normal -from progpy.loading import Piecewise, GaussianNoiseLoadWrapper -from progpy.models import BatteryCircuit -from statistics import mean - -def run_example(): - m = BatteryCircuit() - - ## Example 1: Variable (piecewise) loading - future_loading = Piecewise( - m.InputContainer, - [600, 900, 1800, 3000, float('inf')], - {'i': [2, 1, 4, 2, 3]}) - - # Simulate to threshold - options = { - 'save_freq': 100, # Frequency at which results are saved - 'dt': 2 # Timestep - } - simulated_results = m.simulate_to_threshold(future_loading, **options) - - # Now lets plot the inputs and event_states - simulated_results.inputs.plot(ylabel = 'Variable Load Current (amps)', xlabel='time (s)') - simulated_results.event_states.plot(ylabel = 'Variable Load Event State', xlabel='time (s)') - - ## Example 2: Moving Average loading - # This is useful in cases where you are running reoccuring simulations, and are measuring the actual load on the system, - # but don't have a good way of predicting it, and you expect loading to be steady - - from progpy.loading import MovingAverage - - future_loading = MovingAverage(m.InputContainer) - - # Now lets say you have some measured loads to add - measured_loads = [10, 11.5, 12.0, 8, 2.1, 1.8, 1.99, 2.0, 2.01, 1.89, 1.92, 2.01, 2.1, 2.2] - - # We're going to feed these into the future loading eqn - for load in measured_loads: - future_loading.add_load({'i': load}) - - # Now the future_loading eqn is setup to use the moving average of whats been seen - # Simulate to threshold - simulated_results = m.simulate_to_threshold(future_loading, **options) - - # Now lets plot the inputs and event_states - simulated_results.inputs.plot(ylabel = 'Moving Average Current (amps)', xlabel='time (s)') - simulated_results.event_states.plot(ylabel = 'Moving Average Event State', xlabel='time (s)') - - # In this case, this estimate is wrong because loading will not be steady, but at least it would give you an approximation. - - # If more measurements are received, the user could estimate the moving average here and then run a new simulation. - - ## Example 3: Gaussian Distribution - # In this example we will still be doing a variable loading like the first option, but we are going to use a - # gaussian distribution for each input. - future_loading = Piecewise( - m.InputContainer, - [600, 900, 1800, 3000, float('inf')], - {'i': [2, 1, 4, 2, 3]}) - future_loading_with_noise = GaussianNoiseLoadWrapper(future_loading, 0.2) - - # Simulate to threshold - simulated_results = m.simulate_to_threshold(future_loading_with_noise, **options) - - # Now lets plot the inputs and event_states - simulated_results.inputs.plot(ylabel = 'Variable Gaussian Current (amps)', xlabel='time (s)') - simulated_results.event_states.plot(ylabel = 'Variable Gaussian Event State', xlabel='time (s)') - - # Example 4: Gaussian- increasing with time - # For this we're using moving average. This is realistic because the further out from current time you get, - # the more uncertainty there is in your prediction. - - def future_loading(t, x=None): - std = future_loading.base_std + future_loading.std_slope * (t - future_loading.t) - return {key : normal(future_loading.load[key], std) for key in future_loading.load.keys()} - future_loading.load = {key : 0 for key in m.inputs} - future_loading.base_std = 0.001 - future_loading.std_slope = 1e-4 - future_loading.t = 0 - - # Lets define another function to handle the moving average logic - window = 10 # Number of elements in window - def moving_avg(i): - for key in m.inputs: - moving_avg.loads[key].append(i[key]) - if len(moving_avg.loads[key]) > window: - del moving_avg.loads[key][0] # Remove first item - - # Update future loading eqn - future_loading.load = {key : mean(moving_avg.loads[key]) for key in m.inputs} - moving_avg.loads = {key : [] for key in m.inputs} - - # OK, we've setup the logic of the moving average. - # Now lets say you have some measured loads to add - measured_loads = [10, 11.5, 12.0, 8, 2.1, 1.8, 1.99, 2.0, 2.01, 1.89, 1.92, 2.01, 2.1, 2.2] - - # We're going to feed these into the future loading eqn - for load in measured_loads: - moving_avg({'i': load}) - - # Simulate to threshold - simulated_results = m.simulate_to_threshold(future_loading, **options) - - # Now lets plot the inputs and event_states - simulated_results.inputs.plot(ylabel = 'Moving Average Current (amps)', xlabel='time (s)') - simulated_results.event_states.plot(ylabel = 'Moving Average Event State', xlabel='time (s)') - - # In this example future_loading.t has to be updated with current time before each prediction. - - # Example 5 Function of state - # here we're pretending that input is a function of SOC. It increases as we approach SOC - - def future_loading(t, x=None): - if x is not None: - event_state = future_loading.event_state(x) - return m.InputContainer({'i': future_loading.start + (1-event_state['EOD']) * future_loading.slope}) # default - return m.InputContainer({'i': future_loading.start}) - future_loading.t = 0 - future_loading.event_state = m.event_state - future_loading.slope = 2 # difference between input with EOD = 1 and 0. - future_loading.start = 0.5 - - # Simulate to threshold - simulated_results = m.simulate_to_threshold(future_loading, **options) - - # Now lets plot the inputs and event_states - simulated_results.inputs.plot(ylabel = 'f(x) Current (amps)', xlabel='time (s)') - simulated_results.event_states.plot(ylabel = 'f(x) Event State', xlabel='time (s)') - - # In this example future_loading.t has to be updated with current time before each prediction. - - # Show plots - plt.show() - -# This allows the module to be executed directly -if __name__ == '__main__': - run_example() From 2e9ba11786919d2cb6fba9e4ccc03d7cf685afa8 Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Wed, 23 Aug 2023 13:05:38 -0700 Subject: [PATCH 3/7] add mention of seeds --- examples/01_Simulation.ipynb | 37 ++++++++++++++++++++++++++++++++---- 1 file changed, 33 insertions(+), 4 deletions(-) diff --git a/examples/01_Simulation.ipynb b/examples/01_Simulation.ipynb index 03ba6d6..8a1fedb 100644 --- a/examples/01_Simulation.ipynb +++ b/examples/01_Simulation.ipynb @@ -656,7 +656,7 @@ "metadata": {}, "outputs": [], "source": [ - "loading = GaussianNoiseLoadWrapper(loading, 0.2)" + "loading_with_noise = GaussianNoiseLoadWrapper(loading, 0.2)" ] }, { @@ -674,7 +674,7 @@ "metadata": {}, "outputs": [], "source": [ - "results = m.simulate_to_threshold(loading, save_freq=100)\n", + "results = m.simulate_to_threshold(loading_with_noise, save_freq=100)\n", "fig = results.inputs.plot()" ] }, @@ -691,7 +691,7 @@ "metadata": {}, "outputs": [], "source": [ - "results = m.simulate_to_threshold(loading, save_freq=100)\n", + "results = m.simulate_to_threshold(loading_with_noise, save_freq=100)\n", "fig = results.inputs.plot()" ] }, @@ -717,6 +717,30 @@ "source": [ "Note that the noise in input can be seen in the resulting output plots.\n", "\n", + "The seed can be set in creation of the GaussianNoiseLoadWrapper to ensure repeatable results, for example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "loading_with_noise = GaussianNoiseLoadWrapper(loading, 0.2, seed=2000)\n", + "results = m.simulate_to_threshold(loading_with_noise, save_freq=100)\n", + "fig = results.inputs.plot()\n", + "\n", + "loading_with_noise = GaussianNoiseLoadWrapper(loading, 0.2, seed=2000)\n", + "results = m.simulate_to_threshold(loading_with_noise, save_freq=100)\n", + "fig = results.inputs.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The load profiles in the two examples above are identical because they share the same random seed.\n", + "\n", "In this section we introduced the concept of NoiseWrappers and how they are used to represent uncertainty in future loading. This concept is expecially important when used with prediction (see 9. Prediction). A GuassianNoiseLoadWrapper was used with a Piecewise loading profile to demonstrate it, but NoiseWrappers can be applied to any loading object or function, including the advanced profiles introduced in the next section." ] }, @@ -900,7 +924,12 @@ "## Conclusions\n", "...\n", "\n", - "simulate_to_event(), events" + "simulate_to_event(), events\n", + "loading.plot()\n", + "guassian std dev for each event.\n", + "Guassian std dev as function of t\n", + "Switch to use generator so you can set seed\n", + "remove x=None requirement\n" ] } ], From f0e328cf540b387cda57c32b8043c55bad464b5a Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Mon, 28 Aug 2023 14:55:44 -0700 Subject: [PATCH 4/7] Resolve code review comments --- examples/01_Simulation.ipynb | 87 +++++++++++++++++++++++------------- 1 file changed, 56 insertions(+), 31 deletions(-) diff --git a/examples/01_Simulation.ipynb b/examples/01_Simulation.ipynb index 8a1fedb..95148bc 100644 --- a/examples/01_Simulation.ipynb +++ b/examples/01_Simulation.ipynb @@ -348,7 +348,7 @@ "source": [ "The previous examples feature a simple ThrownObject model, which does not have any inputs. Unlike ThrownObject, most prognostics models have some sort of [input](https://nasa.github.io/progpy/glossary.html#term-input). The input is some sort of control or loading applied to the system being modeled. In this section we will describe how to simulate a model which features an input.\n", "\n", - "In this example we will be using the BatteryCircuit model from the models subpackage (see 3. Included Models). This is a simple battery discharge model where the battery is represented by an equivilant circuit.\n", + "In this example we will be using the BatteryCircuit model from the models subpackage (see 3. Included Models). This is a simple battery discharge model where the battery is represented by an equivalent circuit.\n", "\n", "Like the past examples, we start by importing and creating the model." ] @@ -359,7 +359,7 @@ "metadata": {}, "outputs": [], "source": [ - "from prog_models.models import BatteryCircuit\n", + "from progpy.models import BatteryCircuit\n", "m = BatteryCircuit()" ] }, @@ -367,7 +367,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "You can see the batteries inputs, states, and outputs (described above) by accessing these attributes." + "You can see the battery's inputs, states, and outputs (described above) by accessing these attributes." ] }, { @@ -384,9 +384,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Consulting the [model documentation](https://nasa.github.io/progpy/api_ref/prog_models/IncludedModels.html#prog_models.models.BatteryCircuit), I see that the outputs (i.e., measurable values) of the model are temperature (`t`) and voltage (`v`). The model's input is the current (`i`) drawn from the model.\n", + "Consulting the [model documentation](https://nasa.github.io/progpy/api_ref/prog_models/IncludedModels.html#prog_models.models.BatteryCircuit), we see that the outputs (i.e., measurable values) of the model are temperature (`t`) and voltage (`v`). The model's input is the current (`i`) drawn from the battery.\n", "\n", - "If we try to simulate as we do above, it wouldn't work because the battery discharge is a function of the current (`i`) drawn from the battery. Simulation for a model like this requires that we define the future load. " + "If we try to simulate as we do above (without specifying loading), it wouldn't work because the battery discharge is a function of the current (`i`) drawn from the battery. Simulation for a model like this requires that we define the future load. " ] }, { @@ -411,11 +411,43 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next, we define a loading profile. Piecewise loader takes 3 arguments: 1. the model InputContainer, 2. times and 3. loads, each of these are explained in more detail below.\n", + "Next, we define a loading profile. Piecewise loader takes 3 arguments: 1. the model InputContainer, 2. times and 3. loads. Each of these are explained in more detail below.\n", "\n", - "The model input container is a class for representing the input for a model. It's a class attribute for every model, and is specific to that model. It can be found at m.InputContainer\n", - "\n", - "Times are the ending time for each load. For example, if times were [5, 7, 10], then the first load would apply until t=5, then the second load would apply for 2 seconds, following by the third load for 5 more seconds. \n", + "The model input container is a class for representing the input for a model. It's a class attribute for every model, and is specific to that model. It can be found at m.InputContainer. For example," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m.InputContainer" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "InputContainers are initialized with either a dictionary or an column vector, for example:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(m.InputContainer({'i': 3}))\n", + "import numpy as np\n", + "print(m.InputContainer(np.vstack((2.3, ))))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The second and third argument for the loading profile is times and loads. Together, the 'times' and 'loads' arguments specify what load is applied to the system at what times throughout simulation. The values in 'times' specify the ending time for each load. For example, if times were [5, 7, 10], then the first load would apply until t=5, then the second load would apply for 2 seconds, following by the third load for 5 more seconds. \n", "\n", "Loads are a dictionary of arrays, where the keys of the dictionary are the inputs to the model (for a battery, just current `i`), and the values in the array are the value at each time in times. If the loads array is one longer than times, then the last value is the \"default load\", i.e., the load that will be applied after the last time has passed.\n", "\n", @@ -440,7 +472,7 @@ "source": [ "In this case, the current drawn (`i`) is 2 amps until t is 600 seconds, then it is 1 for the next 300 seconds (until 900 seconds), etc. The \"default load\" is 3, meaning that after the last time has passed (3000 seconds) a current of 3 will be drawn. \n", "\n", - "Now that we have this load profile, lets run a simulation with our model" + "Now that we have this load profile, let's run a simulation with our model" ] }, { @@ -490,7 +522,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In this example we simulated to threshold, loading the system using a simple piecewise load profile. This is the most common load profile and will probably work for most cases " + "In this example we simulated to threshold, loading the system using a simple piecewise load profile. This is the most common load profile and will probably work for most cases." ] }, { @@ -504,7 +536,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Another common loading scheme is the moving-average load. This loading scheme assumes that the load will continue like it's seen in the past. This is useful when you dont know the exact load, but you expect it to be consistant.\n", + "Another common loading scheme is the moving-average load. This loading scheme assumes that the load will continue like it's seen in the past. This is useful when you don't know the exact load, but you expect it to be consistent.\n", "\n", "Like with Piecewise loading, the first step it to import the loading class. In this case, `progpy.loading.MovingAverage`" ] @@ -557,7 +589,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In practice the add_load method should be called whenever there's new input (i.e., load) information. The MovingAverage load estimator averages over a window of elemetns, configurable at construction using the window argument (e.g., MovingAverage(m.InputContainer, window=12))\n", + "In practice the add_load method should be called whenever there's new input (i.e., load) information. The MovingAverage load estimator averages over a window of elements, configurable at construction using the window argument (e.g., MovingAverage(m.InputContainer, window=12))\n", "\n", "Now the configured load estimator can be used in simulation. " ] @@ -607,9 +639,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that the voltage and temperature curves are much cleaner. They dont have the jumps present in the piecewise loading example. This is due to the constant loading.\n", + "The voltage and temperature curves are much cleaner. They don't have the jumps present in the piecewise loading example. This is due to the constant loading.\n", "\n", - "In this example we simulated to threshold, loading the system using a constant load profile caluclated using the moving average load estimator. This load estimator needs to be updated with the add_load method whenever new loading data is available. " + "In this example we simulated to threshold, loading the system using a constant load profile calculated using the moving average load estimator. This load estimator needs to be updated with the add_load method whenever new loading data is available. " ] }, { @@ -623,7 +655,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Typically, users have an idea of what loading will look like, but there is some uncertainty. Future load extimates are hardly ever known exactly. This is where load wrappers like the `progpy.loading.GaussianNoiseLoadWrapper` come into play. The GaussianNoiseLoadWrapper wraps around another load profile, adding a random amount of noise, sampled from a gaussian distribution, at each step. In simulation this will show some variability, but this becomes more important in prediction (see 9. Prediction).\n", + "Typically, users have an idea of what loading will look like, but there is some uncertainty. Future load estimates are hardly ever known exactly. This is where load wrappers like the `progpy.loading.GaussianNoiseLoadWrapper` come into play. The GaussianNoiseLoadWrapper wraps around another load profile, adding a random amount of noise, sampled from a Gaussian distribution, at each step. This will show some variability in simulation, but this becomes more important in prediction (see 9. Prediction).\n", "\n", "In this example we will repeat the Piecewise load example, this time using the GaussianNoiseLoadWrapper to represent our uncertainty in our future load estimate. \n", "\n", @@ -663,7 +695,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In this case we're adding guassian noise with a standard deviation of 0.2 to the result of the previous load estimator.\n", + "In this case we're adding Gaussian noise with a standard deviation of 0.2 to the result of the previous load estimator.\n", "\n", "Now let's simulate and look at the input profile." ] @@ -741,7 +773,7 @@ "source": [ "The load profiles in the two examples above are identical because they share the same random seed.\n", "\n", - "In this section we introduced the concept of NoiseWrappers and how they are used to represent uncertainty in future loading. This concept is expecially important when used with prediction (see 9. Prediction). A GuassianNoiseLoadWrapper was used with a Piecewise loading profile to demonstrate it, but NoiseWrappers can be applied to any loading object or function, including the advanced profiles introduced in the next section." + "In this section we introduced the concept of NoiseWrappers and how they are used to represent uncertainty in future loading. This concept is especially important when used with prediction (see 9. Prediction). A GaussianNoiseLoadWrapper was used with a Piecewise loading profile to demonstrate it, but NoiseWrappers can be applied to any loading object or function, including the advanced profiles introduced in the next section." ] }, { @@ -755,13 +787,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "For most applications, the standard load estimation classes can be used to represent a users uexpectation of future loading. However, there are some cases where load is some complex combination of time and state that cannot be represented by these classes. This section briefly describes a few of these cases. \n", + "For most applications, the standard load estimation classes can be used to represent a user's expectation of future loading. However, there are some cases where load is some complex combination of time and state that cannot be represented by these classes. This section briefly describes a few of these cases. \n", "\n", - "The first example is similar to the last one, in that there is gaussian noise added to some underlying load profile. In this case the magnitude of noise increases linerally with time. This is an important example, as it allows us to represent a case where loading further out in time has more uncertainty (i.e., is less well known). This is common for many prognostic usecases.\n", + "The first example is similar to the last one, in that there is Gaussian noise added to some underlying load profile. In this case the magnitude of noise increases linearly with time. This is an important example, as it allows us to represent a case where loading further out in time has more uncertainty (i.e., is less well known). This is common for many prognostic use-cases.\n", "\n", - "Custom load profiles can be represented either as a function (t, x=None) -> u, where t is time, x is state, and u is input. or as a class which implements the __call__ method with the same profile as the function.\n", + "Custom load profiles can be represented either as a function (t, x=None) -> u, where t is time, x is state, and u is input, or as a class which implements the __call__ method with the same profile as the function.\n", "\n", - "In this case we will use the first method (i.e., the function). We will define a function that will use a defined slope (derivitive of standard deviation with time)." + "In this case we will use the first method (i.e., the function). We will define a function that will use a defined slope (derivative of standard deviation with time)." ] }, { @@ -803,7 +835,7 @@ "source": [ "Note how the noise in the input signal increases with time. Since this is a random process, if you were to run this again you would get a different result.\n", "\n", - "Here is the corresponding output. Note you can see the effects of the increasingly erradic input in the voltage curve." + "Here is the corresponding output. Note you can see the effects of the increasingly erratic input in the voltage curve." ] }, { @@ -922,14 +954,7 @@ "metadata": {}, "source": [ "## Conclusions\n", - "...\n", - "\n", - "simulate_to_event(), events\n", - "loading.plot()\n", - "guassian std dev for each event.\n", - "Guassian std dev as function of t\n", - "Switch to use generator so you can set seed\n", - "remove x=None requirement\n" + "...\n" ] } ], From 95fd871703c3f8bdd0d3f4ce9dc01ad8b35b1e39 Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Wed, 23 Aug 2023 13:07:38 -0700 Subject: [PATCH 5/7] Add seed --- src/progpy/loading/gaussian_wrapper.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/progpy/loading/gaussian_wrapper.py b/src/progpy/loading/gaussian_wrapper.py index 25e4cdc..df7e942 100644 --- a/src/progpy/loading/gaussian_wrapper.py +++ b/src/progpy/loading/gaussian_wrapper.py @@ -11,13 +11,18 @@ class GaussianNoiseLoadWrapper(): This is a simple wrapper for future loading functions that adds gaussian noise to the inputs. It takes a future loading function and a standard deviation and returns a new future loading function that adds gaussian noise to the inputs. - Parameters + Arguments ---------- - fcn : Callable + fcn: Callable The future loading function to wrap - std : float + std: float The standard deviation of the gaussian noise to add + Keyword Arguments + ------------------- + seed: {int, SeedSequence, BitGenerator, Generator}, optional + The seed for random number generator. This can be set to make results repeatable. + Example ------- >>> from progpy.loading import GaussianNoiseLoadWrapper @@ -25,9 +30,10 @@ class GaussianNoiseLoadWrapper(): >>> future_load = GaussianNoiseLoadWrapper(future_load, STANDARD_DEV) >>> m.simulate_to_threshold(future_load) """ - def __init__(self, fcn: Callable, std: float): + def __init__(self, fcn: Callable, std: float, seed: int = None): self.fcn = fcn self.std = std + self.gen = np.random.default_rng(seed) def __call__(self, t: float, x=None): """ @@ -42,5 +48,5 @@ def __call__(self, t: float, x=None): """ input = self.fcn(t, x) for key, value in input.items(): - input[key] = np.random.normal(value, self.std) + input[key] = self.gen.normal(value, self.std) return input From 70b36908a3696d24bc08da297d0133321f887f87 Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Wed, 23 Aug 2023 13:18:05 -0700 Subject: [PATCH 6/7] Add test --- tests/test_loading.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/tests/test_loading.py b/tests/test_loading.py index 0cc4b45..f2bf688 100644 --- a/tests/test_loading.py +++ b/tests/test_loading.py @@ -2,7 +2,8 @@ import unittest -from progpy.loading import Piecewise +from progpy.loading import Piecewise, GaussianNoiseLoadWrapper + class Testloading(unittest.TestCase): def test_piecewise_construction(self): @@ -23,6 +24,27 @@ def test_piecewise_construction(self): # Should work with one more load Piecewise({}, [1, 2], {'a': [1, 2, 3]}) + + def test_gaussian_seed(self): + def loading(t, x=None): + return {'a': 10} + + # Default: two values should be different (because of randomness) + loading_with_noise = GaussianNoiseLoadWrapper(loading, 10) + load1 = loading_with_noise(10) + + loading_with_noise = GaussianNoiseLoadWrapper(loading, 10) + load2 = loading_with_noise(10) + + self.assertNotEqual(load1['a'], load2['a']) + + # Setting seed, two values should be the same now + loading_with_noise = GaussianNoiseLoadWrapper(loading, 10, seed=550) + load1 = loading_with_noise(10) + + loading_with_noise = GaussianNoiseLoadWrapper(loading, 10, seed=550) + load2 = loading_with_noise(10) + self.assertEqual(load1['a'], load2['a']) # This allows the module to be executed directly def main(): From 388d39029703c220b94a49313edba4f51915a9dc Mon Sep 17 00:00:00 2001 From: Christopher Teubert Date: Wed, 13 Sep 2023 14:41:55 -0700 Subject: [PATCH 7/7] resolve review comments --- examples/01_Simulation.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/01_Simulation.ipynb b/examples/01_Simulation.ipynb index 95148bc..88951f9 100644 --- a/examples/01_Simulation.ipynb +++ b/examples/01_Simulation.ipynb @@ -429,7 +429,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "InputContainers are initialized with either a dictionary or an column vector, for example:" + "InputContainers are initialized with either a dictionary or a column vector, for example:" ] }, { @@ -447,7 +447,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The second and third argument for the loading profile is times and loads. Together, the 'times' and 'loads' arguments specify what load is applied to the system at what times throughout simulation. The values in 'times' specify the ending time for each load. For example, if times were [5, 7, 10], then the first load would apply until t=5, then the second load would apply for 2 seconds, following by the third load for 5 more seconds. \n", + "The second and third arguments for the loading profile are times and loads. Together, the 'times' and 'loads' arguments specify what load is applied to the system at what times throughout simulation. The values in 'times' specify the ending time for each load. For example, if times were [5, 7, 10], then the first load would apply until t=5, then the second load would apply for 2 seconds, following by the third load for 3 more seconds. \n", "\n", "Loads are a dictionary of arrays, where the keys of the dictionary are the inputs to the model (for a battery, just current `i`), and the values in the array are the value at each time in times. If the loads array is one longer than times, then the last value is the \"default load\", i.e., the load that will be applied after the last time has passed.\n", "\n",