Model Description

The Simulink model sdoPopulationInflux models a simple ecology where an organism population growth is limited by 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];

Compare Measured Data to Initial Simulated Output

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');

Estimate Parameters

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 30-Jul-2020 18:27:29

                                          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');

Put Model in Steady-State During Estimation

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 30-Jul-2020 18:27:39

                                          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');

Related Examples

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)