Оценка параметров модели с использованием нескольких экспериментов (код)

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

Этот пример показывает, как оценить параметры модели из нескольких наборов экспериментальных данных. Вы оцениваете параметры системы масса-пружина-демпфер.

Откройте модель и получите экспериментальные данные

Этот пример использует 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')

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