В этом примере показано, как установить модель на установившийся в процессе оценки параметра. Установка модели к установившемуся важна во многих приложениях, таких как динамика самолета и энергосистемы. Этот пример использует модель демографической динамики.
Модель 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 26-Jul-2019 20:15:35
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 26-Jul-2019 20:16:01
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)