Установите модель в установившееся состояние при оценке параметров (кода)

Открыть скрипт

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

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

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

$\frac{d y}{dt} = R (1 - \frac{y}{K}) (y + 10)$

- неотъемлемая скорость роста населения организма.

- грузоподъемность окружения.

Также происходит приток других представителей организма из соседнего окружения. В модели используются нормированные модули.

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

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 в качестве начального предположения для этого параметра. Известно, что она составляет по меньшей мере 0,1.

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

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

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 27-Jan-2021 14:02:22

                                          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

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

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 27-Jan-2021 14:02:32

                                          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

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

Документация

Установите модель в установившееся состояние при оценке параметров (кода)

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

Сравнение измеренных данных с исходным моделируемым выходом

Оценка параметров

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

Похожие примеры

Simulink Design Optimization

Поддержка