Notebook

2. Parameter Estimation¶

Parameter estimation is used to tune the parameters of a general model so that its behavior matches that of a specific system. For example, the parameters of a battery model can be tuned to configure the model to more accurately describe the behavior of a specific battery.

Generally, parameter estimation is done by tuning the parameters of the model so that the simulation (see 01 Simulation) best matches the behavior observed in some available data. This is done using a mixture of data, knowledge (e.g., from system specs), and intuition. For large, complex models, it can be VERY difficult and computationally expensive.

In ProgPy, parameter estimation is done using the progpy.PrognosticsModel.estimate_params() method. This method takes input and output data from one or more runs, and uses scipy.optimize.minimize function to estimate the parameters of the model. For more information, refer to the documentation here.

A few definitions:

keys (list[str]): Parameter keys to optimize
times (list[float]): Array of times for each run
inputs (list[InputContainer]): Array of input containers where inputs[x] corresponds to times[x]
outputs (list[OutputContainer]): Array of output containers where outputs[x] corresponds to times[x]
method (str, optional): Optimization method. See scipy.optimize.minimize
tol (int, optional): Tolerance for termination. Depending on the provided minimization method, specifying tolerance sets solver-specific options to tol
error_method (str, optional): Method to use in calculating error. See calc_error for options
bounds (tuple or dict, optional): Bounds for optimization in format ((lower1, upper1), (lower2, upper2), ...) or {key1: (lower1, upper1), key2: (lower2, upper2), ...}
options (dict, optional): Options passed to optimizer. See scipy.optimize.minimize for options

Simple Example¶

In this example, we estimate the model parameters from data. In general, the data will usually be collected from the physical system or from a different model (model surrogacy). In this case, we will use example data.

In [ ]:

times = [0, 1, 2, 3, 4, 5, 6, 7]
inputs = [{}] * 8
outputs = [
    {"x": 1.83},
    {"x": 36.5091999066245},
    {"x": 60.05364349596605},
    {"x": 73.23733081022635},
    {"x": 76.47528104941956},
    {"x": 69.9146810161441},
    {"x": 53.74272753819968},
    {"x": 28.39355725512131},
]

First, we will import a model from the ProgPy Package. For this example, we will be using the simple ThrownObject model.

In [ ]:

from progpy.models import ThrownObject

We can now build a model with a best guess for the parameters. We will guess that our thrower is 20 meters tall, has a throwing speed of 3.1 $m/s$, and that acceleration due to gravity is 15 $m/s^2$. However, given our times, inputs, and outputs, we can clearly tell this is not true! Let's see if parameter estimation can fix this.

In [ ]:

m = ThrownObject(thrower_height=20, throwing_speed=3.1, g=15)

Next we will define specific parameters that we want to estimate. We can pass the desired parameters to our keys keyword argument.

In [ ]:

keys = ["thrower_height", "throwing_speed", "g"]

To really see what estimate_params() is doing, we will print out the state before executing the estimation.

In [ ]:

# Printing state before
print("Model configuration before")
for key in keys:
    print("-", key, m[key])
print(" Error: ", m.calc_error(times, inputs, outputs, dt=0.1))

Notice that the error is quite high. This indicates that the parameters are not accurate.

Now, we will run estimate_params() with the data to correct these parameters.

In [ ]:

m.estimate_params(times=times, inputs=inputs, outputs=outputs, keys=keys, dt=0.1)

Let's see what the new parameters are after estimation.

In [ ]:

print("\nOptimized configuration")
for key in keys:
    print("-", key, m[key])
print(" Error: ", m.calc_error(times, inputs, outputs, dt=0.1))

Sure enough, parameter estimation determined that the thrower's height wasn't 20m. Instead, it was closer to 1.8m, a much more reasonable height. Parameter estimation also correctly estimated g as ~-9.81 $m/s^2$ and throwing speed at around 40 $m/s$, the values used to generate our example data.

Using Tol¶

An additional feature of the estimate_params() function is the tolerance feature, or tol. The exact function that the tol argument uses is specific to the method used. For example, the tol argument for the Nelder-Mead method is the change in the lowest error and its corresponding parameter values between iterations. The difference between iterations for both of these must be below tol for parameter estimation to converge.

For example, if in the nth iteration of the optimizer above the best error was 2e-5 and the cooresponding values were thrower_height=1.8, throwing_speed=40, and g=-9.8 and at the n+1th iteration the best error was 1e-5 and the cooresponding values were thrower_height=1.85, throwing_speed=39.5, and g=-9.81, then the difference in error would be 1e-5 and the difference in parameter values would be

$$\sqrt{(1.85 - 1.8)^2 + (40 - 39.5)^2 + (9 - 9.81)^2} = 0.5025932749$$

In this case, error would meet a tol of 1e-4, but the parameters would not, so optimization would continue. For more information, see the scipy.optimize.minimize documentation.

In our previous example, note that our total error was roughly 6e-10 after the estimate_params() call, using the default tol of 1e-4. Now, let us see what happens to the parameters when we pass a tolerance of 1e-6.

In [ ]:

m = ThrownObject(thrower_height=20, throwing_speed=3.1, g=15)
m.estimate_params(
    times=times, inputs=inputs, outputs=outputs, keys=keys, dt=0.1, tol=1e-6
)
print("\nOptimized configuration")
for key in keys:
    print("-", key, m[key])
print(" Error: ", m.calc_error(times, inputs, outputs, dt=0.1))

As expected, reducing the tolerance leads to a decrease in the overall error, resulting in more accurate parameters.

Note that if we were to set a high tolerance, such as 10, our error would consequently be very high! Also note that the tol value is for scipy minimize. It is different but strongly correlated to the result of calc_error. For more information on how the tol feature works, please refer to scipy's minimize() documentation.

You can also adjust the metric that is used to estimate parameters by setting the error_method to a different calc_error() method (see example below). Default is Mean Squared Error (MSE). See calc_error() method for list of options.

In [ ]:

m["thrower_height"] = 3.1
m["throwing_speed"] = 29

# Using MAE, or Mean Absolute Error instead of the default Mean Squared Error.
m.estimate_params(
    times=times,
    inputs=inputs,
    outputs=outputs,
    keys=keys,
    dt=0.1,
    tol=1e-9,
    error_method="MAX_E",
)
print("\nOptimized configuration")
for key in keys:
    print("-", key, m[key])
print(" Error: ", m.calc_error(times, inputs, outputs, dt=0.1, method="MAX_E"))

Note that MAX_E is frequently better at capturing tail behavior in many prognostic models.

Handling Noise with Multiple Runs¶

In the previous two examples, we demonstrated how to use estimate_params() using a clearly defined ThrownObject model. However we assumed that there would be no noise in the data used to estimate parameters. This is almost never the case in real life.

In this example, we'll show how to use estimate_params() with noisy data. First, let's repeat the previous example, this time generating data from a noisy model.

In [ ]:

m = ThrownObject(process_noise=1)
results = m.simulate_to_threshold(save_freq=0.5, dt=("auto", 0.1))

# Resetting parameters to their incorrectly set values.
m["thrower_height"] = 20
m["throwing_speed"] = 3.1
m["g"] = 15
keys = ["thrower_height", "throwing_speed", "g"]

In [ ]:

m.estimate_params(
    times=results.times, inputs=results.inputs, outputs=results.outputs, keys=keys
)
print("\nOptimized configuration")
for key in keys:
    print("-", key, m[key])
print(" Error: ", m.calc_error(results.times, results.inputs, results.outputs))

Here, the error from calc_error is low. To have an accurate estimation of the error, we should actually be manually measuring the Absolute Mean Error rather than using calc_error().

The reason being is simple. calc_error() is calculating the error between the simulated and observed data. However, the observed and simulated data in this case are being generated from a model that has noise. In other words, we are comparing the difference of noise to noise, which can lead to inconsistent results.

Let's create a helper function to calculate the Absolute Mean Error between our original and estimated parameters.

In [ ]:

# Creating a new model with the original parameters to compare to the model with noise.
true_Values = ThrownObject()


# Function to determine the Absolute Mean Error (AME) of the model parameters.
def AME(m, keys):
    error = 0
    for key in keys:
        error += abs(m[key] - true_Values[key])
    return error

Using our new AME function, we see that the error isn't as great as we thought.

In [ ]:

AME(m, keys)

Note that the error changes every time due to the randomness of noise.

In [ ]:

for count in range(10):
    m = ThrownObject(process_noise=1)
    results = m.simulate_to_threshold(save_freq=0.5, dt=("auto", 0.1))

    # Resetting parameters to their originally incorrectly set values.
    m["thrower_height"] = 20
    m["throwing_speed"] = 3.1
    m["g"] = 15

    m.estimate_params(
        times=results.times,
        inputs=results.inputs,
        outputs=results.outputs,
        keys=keys,
        dt=0.1,
    )
    error = AME(m, ["thrower_height", "throwing_speed", "g"])
    print(f"Estimate Call Number {count} - AME Error {error}")

This issue with noise can be overcome with more data. Let's repeat the example above, this time using data from multiple runs. First, let's generate the data.

In [ ]:

times, inputs, outputs = [], [], []
m = ThrownObject(process_noise=1)
for count in range(20):
    results = m.simulate_to_threshold(save_freq=0.5, dt=("auto", 0.1))
    times.append(results.times)
    inputs.append(results.inputs)
    outputs.append(results.outputs)

Next, let's reset the parameters to our incorrect values.

In [ ]:

m["thrower_height"] = 20
m["throwing_speed"] = 3.1
m["g"] = 15

Finally, we will call estimate_params() with all the collected data.

In [ ]:

m.estimate_params(times=times, inputs=inputs, outputs=outputs, keys=keys, dt=0.1)
print("\nOptimized configuration")
for key in keys:
    print("-", key, m[key])
error = AME(m, ["thrower_height", "throwing_speed", "g"])
print("AME Error: ", error)

Notice that by using data from multiple runs, we are able to produce a lower AME Error than before. This is because we are able to simulate the noise multiple times, which in turn, allows our estimate_params() to produce a more accurate result since it is given more data to work with.

Simplified Battery¶

The previous examples all used a simple model, the ThrownObject. For large, complex models, it can be VERY difficult and computationally expensive.

In this example, we will estimate the parameters for the simplified battery model. This model is more complex than the ThrownObject model but is still a relatively simple model. This example demonstrates some approaches useful for estimating parameters in complex models, like estimating parameter subsets on data selected to highlight specific features.

Let's prepare some data for parameter estimation. We will be using the datasets subpackage in progpy for this.

Data Prep¶

In [ ]:

from progpy.datasets import nasa_battery

(desc, data) = nasa_battery.load_data(1)

The dataset includes 4 different kinds of runs: trickle, step, reference, random walk. We're going to split the dataset into one example for each of the different types for use later.

Let's take a look at the trickle discharge run first.

In [ ]:

trickle_dataset = data[0]
print(trickle_dataset)
trickle_dataset.plot(
    y=["current", "voltage", "temperature"], subplots=True, xlabel="Time (sec)"
)

Now let's do the same for a reference discharge run (5).

In [ ]:

reference_dataset = data[5]
reference_dataset.plot(
    y=["current", "voltage", "temperature"], subplots=True, xlabel="Time (sec)"
)

Now we will do it for the step runs. Note that this is actually multiple runs that we need to combine. relativeTime resets for each "run". So if we're going to use multiple runs together, we need to stitch these times together.

In [ ]:

data[7]["absoluteTime"] = data[7]["relativeTime"]
for i in range(8, 32):
    data[i]["absoluteTime"] = (
        data[i]["relativeTime"] + data[i - 1]["absoluteTime"].iloc[-1]
    )

Next, we should combine the data into a single dataset and investigate the results.

In [ ]:

import pandas as pd

step_dataset = pd.concat(data[7:32], ignore_index=True)
print(step_dataset)
step_dataset.plot(
    y=["current", "voltage", "temperature"], subplots=True, xlabel="Time (sec)"
)

Finally, let's investigate the random walk discharge. Like the step discharge, we need to stitch together the times and concatenate the data.

In [ ]:

data[35]["absoluteTime"] = data[35]["relativeTime"]
for i in range(36, 50):
    data[i]["absoluteTime"] = (
        data[i]["relativeTime"] + data[i - 1]["absoluteTime"].iloc[-1]
    )

In [ ]:

random_walk_dataset = pd.concat(data[35:50], ignore_index=True)
print(random_walk_dataset)
random_walk_dataset.plot(
    y=["current", "voltage", "temperature"], subplots=True, xlabel="Time (sec)"
)

Now the data is ready for this tutorial, let's dive into it.

Set Up Model¶

In [ ]:

from progpy.models import SimplifiedBattery

m = SimplifiedBattery()

Parameter Estimation¶

Let's take a look at the parameter space.

In [ ]:

m.parameters

We will now test how well it fits the random walk dataset. First, let's prepare the data and future load equation.

In [ ]:

times_rw = random_walk_dataset["absoluteTime"]
inputs_rw = [
    elem[1]["voltage"] * elem[1]["current"] for elem in random_walk_dataset.iterrows()
]
outputs_rw = [{"v": elem[1]["voltage"]} for elem in random_walk_dataset.iterrows()]

import numpy as np


def future_load_rw(t, x=None):
    power = np.interp(t, times_rw, inputs_rw)
    return {"P": power}

We can now evaluate how well the battery matches the data.

In [ ]:

result = m.simulate_to(
    random_walk_dataset["absoluteTime"].iloc[-1], future_load_rw, dt=1, save_freq=100
)
from matplotlib import pyplot as plt

plt.figure()
plt.plot(times_rw, [z for z in random_walk_dataset["voltage"]])
plt.plot(result.times, [z["v"] for z in result.outputs])
plt.xlabel("Time (sec)")
plt.ylabel("Voltage (volts)")
fig = result.event_states.plot()

This is a terrible fit. Clearly, the battery model isn't properly configured for this specific battery. Reading through the paper, we see that the default parameters are for a larger battery pouch present in a UAV, much larger than the 18650 battery that produced our dataset.

To correct this, we need to estimate the model parameters.

There are 7 parameters to set (assuming initial SOC is always 1). We can start with setting a few parameters we know. We know that $v_L$ is about 4.2 (from the battery specs). We also expect that the battery internal resistance is the same as that in the electrochemistry model (which also uses an 18650). Finally, we know that the capacity of this battery is significantly smaller than the default values for the larger pouch battery.

In [ ]:

m["v_L"] = 4.2  # We know this
from progpy.models import BatteryElectroChemEOD

m["R_int"] = BatteryElectroChemEOD.default_parameters["Ro"]
m["E_crit"] /= 4  # Battery capacity is much smaller

Now, let's take a look at the model fit again and see where that got us.

In [ ]:

result_guess = m.simulate_to(
    random_walk_dataset["absoluteTime"].iloc[-1], future_load_rw, dt=1, save_freq=5
)
plt.plot(times_rw, [z for z in random_walk_dataset["voltage"]])
plt.plot(result_guess.times, [z["v"] for z in result_guess.outputs])
plt.xlabel("Time (sec)")
plt.ylabel("Voltage (volts)")

Much better, but not there yet. Next, we need to use the parameter estimation feature to estimate the parameters further. Let's prepare some data. We'll use the trickle, reference, and step datasets for this. These are close enough temporally that we can expect aging effects to be minimal.

NOTE: It is important to use a different dataset to estimate parameters as to test

In [ ]:

times_trickle = trickle_dataset["relativeTime"]
inputs_trickle = [
    {"P": elem[1]["voltage"] * elem[1]["current"]}
    for elem in trickle_dataset.iterrows()
]
outputs_trickle = [{"v": elem[1]["voltage"]} for elem in trickle_dataset.iterrows()]

times_ref = reference_dataset["relativeTime"]
inputs_ref = [
    {"P": elem[1]["voltage"] * elem[1]["current"]}
    for elem in reference_dataset.iterrows()
]
outputs_ref = [{"v": elem[1]["voltage"]} for elem in reference_dataset.iterrows()]

times_step = step_dataset["relativeTime"]
inputs_step = [
    {"P": elem[1]["voltage"] * elem[1]["current"]} for elem in step_dataset.iterrows()
]
outputs_step = [{"v": elem[1]["voltage"]} for elem in step_dataset.iterrows()]

We can print the keys and the error beforehand for reference. The error here is what is used to estimate parameters.

In [ ]:

inputs_reformatted_rw = [
    {"P": elem[1]["voltage"] * elem[1]["current"]}
    for elem in random_walk_dataset.iterrows()
]
all_keys = ["v_L", "R_int", "lambda", "gamma", "mu", "beta", "E_crit"]
print("Model configuration")
for key in all_keys:
    print("-", key, m[key])
error_guess = m.calc_error(
    times=times_rw.to_list(), inputs=inputs_reformatted_rw, outputs=outputs_rw
)
print("Error: ", error_guess)

Next, let's set the bounds on each of the parameters.

For $v_L$ and $R_{int}$, we're defining some small bounds because we have an idea of what they might be. For the others we are saying it's between 0.1 and 10x the default battery. We also are adding a constraint that E_crit must be smaller than the default, since we know it's a smaller battery.

In [ ]:

bounds = {
    "v_L": (3.75, 4.5),
    "R_int": (
        BatteryElectroChemEOD.default_parameters["Ro"] * 0.5,
        BatteryElectroChemEOD.default_parameters["Ro"] * 2.5,
    ),
    "lambda": (0.046 / 10, 0.046 * 10),
    "gamma": (3.355 / 10, 3.355 * 10),
    "mu": (2.759 / 10, 2.759 * 10),
    "beta": (8.482 / 10, 8.482 * 10),
    "E_crit": (202426.858 / 10, 202426.858),  # (smaller than default)
}

Finally, we'll estimate the parameters. See the Parameter Estimation section in the ProgPy documentation for more details.

We can throw all of the data into estimate parameters, but that will take a long time to run and is prone to errors (e.g., getting stuck in local minima). For this example, we will split characterization into parts.

First, we try to capture the base voltage ($v_L$). If we look at the equation above, $v_L$ is the only term that is not a function of either SOC or power. So, for this estimation we use the trickle dataset, where power draw is the lowest, and we only use the first section where SOC can be assumed to be about 1.

In [ ]:

keys = ["v_L"]
m.estimate_params(
    times=trickle_dataset["relativeTime"].iloc[:10].to_list(),
    inputs=inputs_trickle[:10],
    outputs=outputs_trickle[:10],
    keys=keys,
    dt=1,
    bounds=bounds,
)

We will now run the simulation and plot the result.

In [ ]:

print("Model configuration")
for key in all_keys:
    print("-", key, m[key])
error_fit1 = m.calc_error(
    times=times_rw.to_list(), inputs=inputs_reformatted_rw, outputs=outputs_rw
)
print(f"Error: {error_guess}->{error_fit1} ({error_fit1 - error_guess})")

result_fit1 = m.simulate_to(
    random_walk_dataset["absoluteTime"].iloc[-1], future_load_rw, dt=1, save_freq=5
)
plt.plot(times_rw, [z for z in random_walk_dataset["voltage"]], label="ground truth")
plt.plot(result_guess.times, [z["v"] for z in result_guess.outputs], label="guess")
plt.plot(result_fit1.times, [z["v"] for z in result_fit1.outputs], label="fit1")
plt.legend()
plt.xlabel("Time (sec)")
plt.ylabel("Voltage (volts)")

plt.figure()
plt.plot([0, 1], [error_guess, error_fit1])
plt.xlabel("Parameter Estimation Run")
plt.ylabel("Error")
plt.ylim((0, 0.25))

A tiny bit closer, but not significant. Our initial guess (from the packaging) must have been pretty good.

The next step is to estimate the effect of current on the battery. The Parameter $R_{int}$ (internal resistance) effects this. To estimate $R_{int}$ we will use 2 runs where power is not minimal (ref and step runs). Again, we will use only the first couple steps so EOL can be assumed to be 1.

In [ ]:

keys = ["R_int"]
m.estimate_params(
    times=[times_ref.iloc[:5].to_list(), times_step.iloc[:5].to_list()],
    inputs=[inputs_ref[:5], inputs_step[:5]],
    outputs=[outputs_ref[:5], outputs_step[:5]],
    keys=keys,
    dt=1,
    bounds=bounds,
)

Let's look at what that got us.

In [ ]:

print("Model configuration")
for key in all_keys:
    print("-", key, m[key])
error_fit2 = m.calc_error(
    times=times_rw.to_list(), inputs=inputs_reformatted_rw, outputs=outputs_rw
)
print(f"Error: {error_fit1}->{error_fit2} ({error_fit2 - error_fit1})")

result_fit2 = m.simulate_to(
    random_walk_dataset["absoluteTime"].iloc[-1], future_load_rw, dt=1, save_freq=5
)
plt.plot(times_rw, [z for z in random_walk_dataset["voltage"]], label="ground truth")
plt.plot(result_guess.times, [z["v"] for z in result_guess.outputs], label="guess")
plt.plot(result_fit1.times, [z["v"] for z in result_fit1.outputs], label="fit1")
plt.plot(result_fit2.times, [z["v"] for z in result_fit2.outputs], label="fit2")
plt.legend()
plt.xlabel("Time (sec)")
plt.ylabel("Voltage (volts)")

plt.figure()
plt.plot([0, 1, 2], [error_guess, error_fit1, error_fit2])
plt.xlabel("Parameter Estimation Run")
plt.ylabel("Error")
plt.ylim((0, 0.25))

Much better, but not there yet! Finally we need to estimate the effects of SOC on battery performance. This involves all of the rest of the parameters. For this we will use all the rest of the parameters. We will not be using the entire reference curve to capture a full discharge.

Note that we're using the error_method MAX_E, instead of the default MSE. This results in parameters that better estimate the end of the discharge curve and is recommended when estimating parameters that are combined with the event state.

In [ ]:

keys = ["lambda", "gamma", "mu", "beta", "E_crit"]
m.estimate_params(
    times=times_ref.to_list(),
    inputs=inputs_ref,
    outputs=outputs_ref,
    keys=keys,
    dt=1,
    bounds=bounds,
    error_method="MAX_E",
)

We will now run the simulation and plot the result.

In [ ]:

print("Model configuration")
for key in all_keys:
    print("-", key, m[key])
error_fit3 = m.calc_error(
    times=times_rw.to_list(), inputs=inputs_reformatted_rw, outputs=outputs_rw
)
print(f"Error: {error_fit2}->{error_fit3} ({error_fit3 - error_fit2})")

result_fit3 = m.simulate_to(
    random_walk_dataset["absoluteTime"].iloc[-1], future_load_rw, dt=1, save_freq=5
)
plt.plot(times_rw, [z for z in random_walk_dataset["voltage"]], label="ground truth")
plt.plot(result_guess.times, [z["v"] for z in result_guess.outputs], label="guess")
plt.plot(result_fit1.times, [z["v"] for z in result_fit1.outputs], label="fit1")
plt.plot(result_fit2.times, [z["v"] for z in result_fit2.outputs], label="fit2")
plt.plot(result_fit3.times, [z["v"] for z in result_fit3.outputs], label="fit3")
plt.legend()
plt.xlabel("Time (sec)")
plt.ylabel("Voltage (volts)")

plt.figure()
plt.plot([0, 1, 2, 3], [error_guess, error_fit1, error_fit2, error_fit3])
plt.xlabel("Parameter Estimation Run")
plt.ylabel("Error")
plt.ylim((0, 0.25))

This is even better. Now we have an "ok" estimate, ~150 mV (for the sake of a demo). The estimate could be refined further by setting a lower tolerance (tol parameter), or repeating the 4 parameter estimation steps, as shown above.

Parameter estimation is also limited by the model itself. This is a simplified battery model, meaning there were some simplifying assumptions made. It will likely not be able to capture the behavior of a model as well as a higher fidelity model (e.g., BatteryElectroChemEOD).

Conclusion¶

This chapter introduced the concept of parameter estimation, through which the parameters of a physics-based model are estimated. This is done using a mixture of data, knowledge (e.g., from system specs), and intuition. For large, complex models, it can be VERY difficult and computationally expensive. Fortunately, in this case we have a relatively simple model.

In ProgPy a models estimate_params method is used to estimate the parameters. See Parameter Estimation Docs for more details.

In the next notebook, we will be exploring (see 03 Existing Models).