Main Content

This example shows how to set a model to steady-state in the process of parameter estimation. Setting a model to steady-state is important in many applications such as power systems and aircraft dynamics. This example uses a population dynamics model.

The Simulink model `sdoPopulationInflux`

models a simple ecology where an organism population growth is limited by the carrying capacity of the environment:

is the inherent growth rate of the organism population.

is the carrying capacity of the environment.

There is also an influx of other members of the organism from a neighboring environment. The model uses normalized units.

Open the model.

```
modelName = 'sdoPopulationInflux';
open_system(modelName)
```

The file `sdoPopulationInflux_Data.mat`

contains data of the population in the environment as well as the influx of additional organisms from the neighboring environment.

load sdoPopulationInflux_Data.mat; % Timer series: Population_ts Inflow_ts hFig = figure; subplot(2,1,1); plot(Population_ts) subplot(2,1,2); plot(Inflow_ts)

The population starts in a steady state. After some time, there is an influx of organisms from the neighboring environment. Based on the measured data, we want to estimate values for the model parameters.

The parameter `R`

represents the inherent growth rate of the organism. Use 1 as the initial guess for this parameter. It is non-negative.

```
R = sdo.getParameterFromModel(modelName, 'R');
R.Value = 1;
R.Minimum = 0;
```

The parameter `K`

represents the carrying capacity of the environment. Use 2 as the initial guess for this parameter. It is known to be at least 0.1.

```
K = sdo.getParameterFromModel(modelName, 'K');
K.Value = 2;
K.Minimum = 0.1;
```

Collect these parameters into a vector.

v = [R K];

Create an Experiment object.

Exp = sdo.Experiment(modelName);

Associate `Population_ts`

with model output.

Population = Simulink.SimulationData.Signal; Population.Name = 'Population'; Population.BlockPath = [modelName '/Integrator']; Population.PortType = 'outport'; Population.PortIndex = 1; Population.Values = Population_ts;

Add `Population`

to the experiment.

Exp.OutputData = Population;

Associate `Inflow_ts`

with model input.

Inflow = Simulink.SimulationData.Signal; Inflow.Name = 'Population'; Inflow.BlockPath = [modelName '/In1']; Inflow.PortType = 'outport'; Inflow.PortIndex = 1; Inflow.Values = Inflow_ts;

Add `Inflow`

to the experiment.

Exp.InputData = Inflow;

Create a simulation scenario using the experiment, and obtain the simulated output.

```
Exp = setEstimatedValues(Exp, v); % use vector of parameters/states
Simulator = createSimulator(Exp);
Simulator = sim(Simulator);
```

Search for the model output signal in the logged simulation data.

SimLog = find(Simulator.LoggedData, ... get_param(modelName, 'SignalLoggingName') ); PopulationSim = find(SimLog, 'Population');

The model output does not match the data very well, indicating that we need to compute better estimates of the model parameters.

clf; plot(PopulationSim.Values, 'r'); hold on; plot(Population_ts, 'b'); legend('Model Simulation', 'Measured Data', 'Location', 'best');

To estimate parameters, define an objective function to compute the sum squared error between model simulation and measured data.

```
estFcn = @(v) sdoPopulationInflux_Objective(v, Simulator, Exp);
type sdoPopulationInflux_Objective.m
```

function vals = sdoPopulationInflux_Objective(v, Simulator, Exp, OpPointSetup) % Compare model output with data % % Inputs: % v - vector of parameters and/or states % Simulator - used to simulate the model % Exp - Experiment object % OpPointSetup - Object to set up computation of steady-state % operating point % Copyright 2018 The MathWorks, Inc. %Parse inputs if nargin < 4 OpPointSetup = []; end % Requirement setup req = sdo.requirements.SignalTracking; req.Type = '=='; req.Method = 'Residuals'; % Simulate the model Exp = setEstimatedValues(Exp, v); % use vector of parameters/states Simulator = createSimulator(Exp,Simulator); strOT = mat2str(Exp.OutputData(1).Values.Time); if isempty(OpPointSetup) Simulator = sim(Simulator, 'OutputOption', 'AdditionalOutputTimes', 'OutputTimes', strOT); else Simulator = sim(Simulator, 'OutputOption', 'AdditionalOutputTimes', 'OutputTimes', strOT, ... 'OperatingPointSetup', OpPointSetup); end % Compare model output with data SimLog = find(Simulator.LoggedData, ... get_param(Exp.ModelName, 'SignalLoggingName') ); OutputModel = find(SimLog, 'Population'); Model_Error = evalRequirement(req, OutputModel.Values, Exp.OutputData.Values); vals.F = Model_Error;

%Define options for the optimization. % opts = sdo.OptimizeOptions; opts.Method = 'lsqnonlin';

Estimate the parameters.

vOpt = sdo.optimize(estFcn, v, opts); disp(vOpt)

Optimization started 03-Oct-2020 18:49:24 First-order Iter F-count f(x) Step-size optimality 0 5 12.485 1 1 10 1.09824 1.184 0.244 2 15 0.9873 1.088 0.0259 3 20 0.952948 1.217 0.00624 4 25 0.946892 0.9151 0.00197 5 30 0.946484 0.3541 0.00153 Local minimum possible. lsqnonlin stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. (1,1) = Name: 'R' Value: 5.5942 Minimum: 0 Maximum: Inf Free: 1 Scale: 1 Info: [1x1 struct] (1,2) = Name: 'K' Value: 3.2061 Minimum: 0.1000 Maximum: Inf Free: 1 Scale: 2 Info: [1x1 struct] 1x2 param.Continuous

Use the estimated parameter values in the model, and obtain the model response. Search for the model output signal in the logged simulation data.

Exp = setEstimatedValues(Exp, vOpt); Simulator = createSimulator(Exp,Simulator); Simulator = sim(Simulator); SimLog = find(Simulator.LoggedData, ... get_param(modelName, 'SignalLoggingName') ); PopulationSim = find(SimLog, 'Population');

Comparing the measured population data with the optimized model response shows that they still do not match well. There is a transient at the beginning of the model response, where it is markedly different from the measured data.

clf; plot(PopulationSim.Values, 'r'); hold on; plot(Population_ts, 'b'); legend('Model Simulation', 'Measured Data', 'Location', 'best');

To improve the fit between the model and measured data, the model needs to be set to steady-state when parameters are estimated. Define an operating point specification. The input is known from experimental data. Therefore, (1) it should not be treated as a free variable when computing the steady-state operating point, and (2) after the operating point is found, its input should not be used when simulating the model. On the other hand, all the states found when computing the operating point should be used when simulating the model. Create an sdo.OperatingPointSetup to collect the operating point specification, inputs to use, and states to use, so this information can be passed to the objective function and used when simulating the model. You can also provide a fourth argument to sdo.OperatingPointSetup, specifying options for computing the operating point. For example, the option 'graddescent-proj' is often used to find the operating point for systems that use physical modeling.

opSpec = operspec(modelName); opSpec.Inputs(1).Known = true; inputsToUse = []; statesToUse = 1:numel(opSpec.States); OpPointSetup = sdo.OperatingPointSetup(opSpec, inputsToUse, statesToUse);

Estimate the parameters, setting the model to steady-state in the process.

estFcn = @(v) sdoPopulationInflux_Objective(v, Simulator, Exp, OpPointSetup); vOpt = sdo.optimize(estFcn, v, opts); disp(vOpt)

Optimization started 03-Oct-2020 18:49:38 First-order Iter F-count f(x) Step-size optimality 0 5 11.1517 1 1 10 0.025674 0.5732 0.045 2 15 0.00239293 0.3451 0.357 3 20 0.000692478 0.0148 0.00301 4 25 0.00069236 6.539e-05 1.16e-07 Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance. (1,1) = Name: 'R' Value: 0.5953 Minimum: 0 Maximum: Inf Free: 1 Scale: 1 Info: [1x1 struct] (1,2) = Name: 'K' Value: 3.0988 Minimum: 0.1000 Maximum: Inf Free: 1 Scale: 2 Info: [1x1 struct] 1x2 param.Continuous

Use the estimated parameter values in the model, and obtain the model response. Search for the model output signal in the logged simulation data.

Exp = setEstimatedValues(Exp, vOpt); Simulator = createSimulator(Exp,Simulator); Simulator = sim(Simulator, 'OperatingPointSetup', OpPointSetup); SimLog = find(Simulator.LoggedData, ... get_param(modelName, 'SignalLoggingName') ); PopulationSim = find(SimLog, 'Population');

There is no more transient at the beginning of the model response, and there is a much better match between the model response and measured data, which is also reflected by the lower objective/cost function value in the second optimization. All these indicate that we have found a good set of parameter values.

clf; plot(PopulationSim.Values, 'r'); hold on; plot(Population_ts, 'b'); legend('Model Simulation', 'Measured Data', 'Location', 'best');

To learn how to put models in a steady state using the **Parameter Estimator** app, see Set Model to Steady-State When Estimating Parameters (GUI).

Close the model and figure.

bdclose(modelName) close(hFig)