Описание модели

Модель Simulink sdoPopulationInflux моделирует простую экологию, где прирост населения организма ограничивается пропускной способностью среды:

Существует также приток других членов организма соседней средой. Модель использует нормированные единицы.

Откройте модель.

modelName = 'sdoPopulationInflux';
open_system(modelName)

Файл sdoPopulationInflux_Data.mat содержит данные населения в среде, а также притоке дополнительных организмов соседней средой.

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)

Население запускается в устойчивом состоянии. Через какое-то время существует приток организмов соседней средой. На основе результатов измерений мы хотим оценить значения для параметров модели.

Параметр R представляет свойственный темп роста организма. Используйте 1 в качестве исходного предположения для этого параметра. Это неотрицательно.

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

Параметр K представляет пропускную способность среды. Используйте 2 в качестве исходного предположения для этого параметра. Это также неотрицательно.

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

Соберите эти параметры в вектор.

v = [R K];

Сравните результаты измерений с начальным симулированным Выходом

Создайте объект Experiment.

Exp = sdo.Experiment(modelName);

Объединенный Population_ts с выходом модели.

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

Добавьте Population к эксперименту.

Exp.OutputData = Population;

Объединенный Inflow_ts с входом модели.

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

Добавьте Inflow к эксперименту.

Exp.InputData = Inflow;

Создайте сценарий симуляции с помощью эксперимента и получите симулированный выход.

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

Ищите выходной сигнал модели в регистрируемых данных моделирования.

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

Выход модели не совпадает с данными очень хорошо, указывая, что мы должны вычислить лучшие оценки параметров модели.

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

Оцените параметры

Чтобы оценить параметры, задайте целевую функцию, чтобы вычислить квадратичную невязку суммы между симуляцией модели и результатами измерений.

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

Оцените параметры.

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

 Optimization started 30-Jan-2020 01:55:28

                                          First-order 
 Iter F-count        f(x)      Step-size  optimality
    0      5       12.485            1                                         
    1     10      1.09824        1.184        0.252
    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
    Maximum: Inf
       Free: 1
      Scale: 2
       Info: [1x1 struct]

 
1x2 param.Continuous
 

Используйте предполагаемые значения параметров в модели и получите ответ модели. Ищите выходной сигнал модели в регистрируемых данных моделирования.

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

Сравнение измеренных данных о населении с оптимизированным ответом модели показывает, что они все еще не соответствуют хорошо. Существует переходный процесс в начале ответа модели, где это заметно отличается от результатов измерений.

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

Вставьте модель, установившуюся во время оценки

Чтобы улучшить подгонку между данными модели и результатами измерений, модель должна быть установлена в установившийся, когда параметры оцениваются. Задайте спецификацию рабочей точки. Вход известен от экспериментальных данных. Поэтому (1) это не должно быть обработано как свободная переменная при вычислении установившейся рабочей точки и (2) после того, как рабочая точка найдена, ее вход не должен использоваться при симуляции модели. С другой стороны, все состояния нашли, когда вычисление рабочей точки должно использоваться при симуляции модели. Создайте sdo. OperatingPointSetup, чтобы собрать спецификацию рабочей точки, входные параметры, чтобы использовать, и состояния, чтобы использовать, таким образом, эта информация может передаваться целевой функции и использоваться при симуляции модели. Можно также предоставить четвертый аргумент sdo. OperatingPointSetup, задавая опции для вычисления рабочей точки. Например, опция 'graddescent-proj' часто используется, чтобы найти рабочую точку для систем, которые используют физическое моделирование.

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

Оцените параметры, установив модель на установившийся в процессе.

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

 Optimization started 30-Jan-2020 01:55:53

                                          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.368
    3     20  0.000692473      0.01476      0.00305
    4     25   0.00069236    5.837e-05     1.08e-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
    Maximum: Inf
       Free: 1
      Scale: 2
       Info: [1x1 struct]

 
1x2 param.Continuous
 

Используйте предполагаемые значения параметров в модели и получите ответ модели. Ищите выходной сигнал модели в регистрируемых данных моделирования.

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

Больше нет переходного процесса в начале ответа модели, и существует намного лучшее соответствие между ответом модели и результатами измерений, который также отражается более низким значением цели/функции стоимости во второй оптимизации. Все они указывают, что мы нашли хороший набор значений параметров.

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

Связанные примеры

Чтобы изучить, как поместить модели в устойчивое состояние с помощью инструмента Parameter Estimation, см. "Модель набора к Установившемуся Когда Оценка Параметров (графический интерфейс пользователя)".

Закройте модель и фигуру.

bdclose(modelName)
close(hFig)