Этот пример показывает, как установить модель на установившийся в процессе оценки параметра. Установка модели к установившемуся важна во многих приложениях, таких как динамика самолета и энергосистемы. Этот пример использует модель демографической динамики.
Модели sdoPopulationInflux модели Simulink простая экология, где прирост населения организма ограничивается пропускной способностью среды:
свойственный темп роста генеральной совокупности организма.
пропускная способность среды.
Существует также приток других членов организма от соседней среды. Модель использует нормированные единицы.
Откройте модель.
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 11-Jan-2019 03:26:55
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 11-Jan-2019 03:27:20
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)