Этот пример показывает, как оценить параметры модели из нескольких наборов экспериментальных данных. Вы оцениваете параметры системы масса-пружина-демпфер.
Этот пример использует sdoMassSpringDamper модель. Модель включает два интегратора, чтобы смоделировать скорость и положение массы в системе масса-пружина-демпфер.
open_system('sdoMassSpringDamper');
Загрузите данные эксперимента.
load sdoMassSpringDamper_ExperimentData
Переменные texp1, yexp1, texp2, и yexp2 загружаются в рабочую область. yexp1 и yexp2 опишите положение массы для времени texp1 и texp2 соответственно.
Создайте 2-элементный массив объектов эксперимента, чтобы сохранить измеренные данные для двух экспериментов.
Создайте объект эксперимента для первого эксперимента.
Exp = sdo.Experiment('sdoMassSpringDamper');
Создайте объект для хранения выходов измерения положения массы.
MeasuredPos = Simulink.SimulationData.Signal; MeasuredPos.Values = timeseries(yexp1,texp1); MeasuredPos.BlockPath = 'sdoMassSpringDamper/Position'; MeasuredPos.PortType = 'outport'; MeasuredPos.PortIndex = 1; MeasuredPos.Name = 'Position';
Добавьте измеренные данные о положении массы к эксперименту в качестве ожидаемых выходных данных.
Exp.OutputData = MeasuredPos;
Создайте объект, чтобы задать начальное состояние для Velocity блок. Начальная скорость массы составляет 0 м/с.
sVel = sdo.getStateFromModel('sdoMassSpringDamper','Velocity'); sVel.Value = 0; sVel.Free = false;
sVel.Free установлено в false поскольку начальная скорость известна и не нуждается в оценке.
Создайте объект, чтобы задать начальное состояние для Position блок. Задайте предположение начального положения массы. Установите Free поле объекта начального положения, чтобы true так, чтобы оно было оценено.
sPos = sdo.getStateFromModel('sdoMassSpringDamper','Position'); sPos.Free = true; sPos.Value = -0.1;
Добавьте начальные состояния к эксперименту.
Exp.InitialStates = [sVel;sPos];
Создайте 2-элементный массив экспериментов. Поскольку два эксперимента идентичны, за исключением ожидаемых выходных данных, скопируйте первый эксперимент дважды.
Exp = [Exp; Exp];
Измените ожидаемые выходные данные второго объекта эксперимента в Exp.
Exp(2).OutputData.Values = timeseries(yexp2,texp2);
Создайте сценарий симуляции с помощью первого эксперимента и получите моделируемый выход.
Simulator = createSimulator(Exp(1)); Simulator = sim(Simulator);
Поиск сигнала положения в записанных данных моделирования.
SimLog = find(Simulator.LoggedData,get_param('sdoMassSpringDamper','SignalLoggingName')); Position = find(SimLog,'Position');
Получите имитированный сигнал положения для второго эксперимента.
Simulator = createSimulator(Exp(2),Simulator); Simulator = sim(Simulator); SimLog = find(Simulator.LoggedData,get_param('sdoMassSpringDamper','SignalLoggingName')); Position(2) = find(SimLog,'Position');
Постройте график измеренных и моделируемых данных.
Реакция модели не совпадает с экспериментальными выходными данными.
subplot(211) plot(... Position(1).Values.Time,Position(1).Values.Data, ... Exp(1).OutputData.Values.Time, Exp(1).OutputData.Values.Data,'--') title('Experiment 1: Simulated and Measured Responses Before Estimation') ylabel('Position') legend('Simulated Position','Measured Position','Location','SouthEast') subplot(212) plot(... Position(2).Values.Time,Position(2).Values.Data, ... Exp(2).OutputData.Values.Time, Exp(2).OutputData.Values.Data,'--') title('Experiment 2: Simulated and Measured Responses Before Estimation') xlabel('Time (seconds)') ylabel('Position') legend('Simulated Position','Measured Position','Location','SouthEast')
Выберите массу m, коэффициент упругости k, и коэффициент демпфирования b параметры из модели. Задайте, что оценочные значения для этих параметров должны быть положительными.
p = sdo.getParameterFromModel('sdoMassSpringDamper', {'b', 'k', 'm'}); p(1).Minimum = 0; p(2).Minimum = 0; p(3).Minimum = 0;
Получите начальные значения состояния положения, которые будут оценены в эксперименте.
s = getValuesToEstimate(Exp);
s содержит два начальных объекта состояния, оба для Position блок. Каждый объект соответствует эксперименту в Exp.
Сгруппируйте параметры модели и начальные состояния, которые будут оценены вместе.
v = [p;s]
v(1,1) =
Name: 'b'
Value: 100
Minimum: 0
Maximum: Inf
Free: 1
Scale: 128
Info: [1x1 struct]
v(2,1) =
Name: 'k'
Value: 500
Minimum: 0
Maximum: Inf
Free: 1
Scale: 512
Info: [1x1 struct]
v(3,1) =
Name: 'm'
Value: 8
Minimum: 0
Maximum: Inf
Free: 1
Scale: 8
Info: [1x1 struct]
v(4,1) =
Name: 'sdoMassSpringDamper/Position'
Value: -0.1000
Minimum: -Inf
Maximum: Inf
Free: 1
Scale: 0.1250
dxValue: 0
dxFree: 1
Info: [1x1 struct]
v(5,1) =
Name: 'sdoMassSpringDamper/Position'
Value: -0.1000
Minimum: -Inf
Maximum: Inf
Free: 1
Scale: 0.1250
dxValue: 0
dxFree: 1
Info: [1x1 struct]
5x1 param.Continuous
Создайте целевую функцию оценки, чтобы оценить, насколько близко выход симуляции, сгенерированный с использованием оцененных значений параметров, совпадает с измеренными данными.
Используйте анонимную функцию с одним входным параметром, который вызывает sdoMassSpringDamper_Objective функция. Передаем анонимную функцию в sdo.optimize, который оценивает функцию при каждой итерации оптимизации.
estFcn = @(v) sdoMassSpringDamper_Objective(v,Simulator,Exp);
The sdoMassSpringDamper_Objective функция:
Имеет один входной параметр, который задает массу, коэффициент упругости и значения демпфера, а также начальное положение массы.
Имеет один входной параметр, который задает объект эксперимента, содержащий измеренные данные.
Возвращает вектор ошибок между моделируемым и экспериментальным выходами.
The sdoMassSpringDamper_Objective функция требует двух входов, но sdo.optimize требует функцию с одним входным параметром. Чтобы обойти это, estFcn является анонимной функцией с одним входным параметром, v, но он звонит sdoMassSpringDamper_Objective используя два входных параметров, v и Exp.
Для получения дополнительной информации об анонимных функциях см. «Анонимные функции».
The sdo.optimize команда минимизирует возвращаемый аргумент анонимной функции estFcn, то есть остаточные ошибки, возвращенные sdoMassSpringDamper_Objective. Для получения дополнительной информации о том, как написать функцию objective/constraint для использования со sdo.optimize команда, введите help sdoExampleCostFunction в командной строке MATLAB.
Чтобы изучить целевую функцию оценки более подробно, введите edit sdoMassSpringDamper_Objective в командной строке MATLAB.
type sdoMassSpringDamper_Objective
function vals = sdoMassSpringDamper_Objective(v,Simulator,Exp)
%SDOMASSSPRINGDAMPER_OBJECTIVE
%
% The sdoMassSpringDamper_Objective function is used to compare model
% outputs against experimental data.
%
% vals = sdoMassSpringDamper_Objective(v,Exp)
%
% The |v| input argument is a vector of estimated model parameter values
% and initial states.
%
% The |Simulator| input argument is a simulation object used
% simulate the model with the estimated parameter values.
%
% The |Exp| input argument contains the estimation experiment data.
%
% The |vals| return argument contains information about how well the
% model simulation results match the experimental data and is used by
% the |sdo.optimize| function to estimate the model parameters.
%
% see also sdo.optimize, sdoExampleCostFunction
%
% Copyright 2012-2015 The MathWorks, Inc.
%%
% Define a signal tracking requirement to compute how well the model output
% matches the experiment data. Configure the tracking requirement so that
% it returns the tracking error residuals (rather than the
% sum-squared-error) and does not normalize the errors.
%
r = sdo.requirements.SignalTracking;
r.Type = '==';
r.Method = 'Residuals';
r.Normalize = 'off';
%%
% Update the experiments with the estimated parameter values.
%
Exp = setEstimatedValues(Exp,v);
%%
% Simulate the model and compare model outputs with measured experiment
% data.
%
Error = [];
for ct=1:numel(Exp)
Simulator = createSimulator(Exp(ct),Simulator);
Simulator = sim(Simulator);
SimLog = find(Simulator.LoggedData,get_param('sdoMassSpringDamper','SignalLoggingName'));
Position = find(SimLog,'Position');
PositionError = evalRequirement(r,Position.Values,Exp(ct).OutputData.Values);
Error = [Error; PositionError(:)];
end
%%
% Return the residual errors to the optimization solver.
%
vals.F = Error(:);
end
Используйте sdo.optimize функция для оценки значений параметров привода и начального состояния.
Задайте опции оптимизации. Функция оценки sdoMassSpringDamper_Objective возвращает невязки ошибок между моделируемыми и экспериментальными данными и не включает никаких ограничений, что делает эту задачу идеальной для решателя 'lsqnonlin'.
opt = sdo.OptimizeOptions;
opt.Method = 'lsqnonlin';
Оцените параметры. Заметьте, что начальное положение массы оценивается дважды, один раз для каждого эксперимента.
vOpt = sdo.optimize(estFcn,v,opt)
Optimization started 27-Jan-2021 14:01:44
First-order
Iter F-count f(x) Step-size optimality
0 11 0.777696 1
1 22 0.00413099 3.696 0.00648
2 33 0.00118327 0.3194 0.00243
3 44 0.0011106 0.06718 5.09e-05
Local minimum found.
Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
vOpt(1,1) =
Name: 'b'
Value: 58.1959
Minimum: 0
Maximum: Inf
Free: 1
Scale: 128
Info: [1x1 struct]
vOpt(2,1) =
Name: 'k'
Value: 399.9452
Minimum: 0
Maximum: Inf
Free: 1
Scale: 512
Info: [1x1 struct]
vOpt(3,1) =
Name: 'm'
Value: 9.7225
Minimum: 0
Maximum: Inf
Free: 1
Scale: 8
Info: [1x1 struct]
vOpt(4,1) =
Name: 'sdoMassSpringDamper/Position'
Value: 0.2995
Minimum: -Inf
Maximum: Inf
Free: 1
Scale: 0.1250
dxValue: 0
dxFree: 1
Info: [1x1 struct]
vOpt(5,1) =
Name: 'sdoMassSpringDamper/Position'
Value: 0.0994
Minimum: -Inf
Maximum: Inf
Free: 1
Scale: 0.1250
dxValue: 0
dxFree: 1
Info: [1x1 struct]
5x1 param.Continuous
Обновите эксперименты с оцененными значениями параметров.
Exp = setEstimatedValues(Exp,vOpt);
Получите моделируемый выход для первого эксперимента.
Simulator = createSimulator(Exp(1),Simulator); Simulator = sim(Simulator); SimLog = find(Simulator.LoggedData,get_param('sdoMassSpringDamper','SignalLoggingName')); Position(1) = find(SimLog,'Position');
Получите моделируемый выход для второго эксперимента.
Simulator = createSimulator(Exp(2),Simulator); Simulator = sim(Simulator); SimLog = find(Simulator.LoggedData,get_param('sdoMassSpringDamper','SignalLoggingName')); Position(2) = find(SimLog,'Position');
Постройте график измеренных и моделируемых данных.
Реакция модели, использующая оцененные значения параметров, хорошо соответствует выходным данным для экспериментов.
subplot(211) plot(... Position(1).Values.Time,Position(1).Values.Data, ... Exp(1).OutputData.Values.Time, Exp(1).OutputData.Values.Data,'--') title('Experiment 1: Simulated and Measured Responses After Estimation') ylabel('Position') legend('Simulated Position','Measured Position','Location','NorthEast') subplot(212) plot(... Position(2).Values.Time,Position(2).Values.Data, ... Exp(2).OutputData.Values.Time, Exp(2).OutputData.Values.Data,'--') title('Experiment 2: Simulated and Measured Responses After Estimation') xlabel('Time (seconds)') ylabel('Position') legend('Simulated Position','Measured Position','Location','SouthEast')
Обновите модель m, k, и b значений параметров. Не обновляйте начальное значение положения модели, так как это зависит от эксперимента.
sdo.setValueInModel('sdoMassSpringDamper',vOpt(1:3));
Закройте модель
bdclose('sdoMassSpringDamper')