В этом примере показано, как установить модель в стационарное состояние в процессе оценки параметров. Перевод модели в стационарное состояние важен во многих областях применения, таких как энергетические системы и динамика самолета. В этом примере используется модель динамики заполнения.
Модель 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 в качестве начального предположения для этого параметра. Известно, что она составляет по меньшей мере 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. Настройка операционной точки для сбора спецификации операционной точки, входных данных для использования и состояний для использования, чтобы эту информацию можно было передать целевой функции и использовать при моделировании модели. Можно также предоставить четвертый аргумент для sdo. Настройка операционной точки, задающая параметры для вычисления рабочей точки. Например, опция «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');
Сведения о том, как перевести модели в устойчивое состояние с помощью приложения «Оценка параметров», см. в разделе Установка модели в устойчивое состояние при оценке параметров (GUI).
Закройте модель и фигуру.
bdclose(modelName) close(hFig)