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

Открыть сценарий

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

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

В этом примере используется sdoMassSpringDamper модель. Модель включает в себя два интегратора для моделирования скорости и положения массы в системе масса-пружина-демпфер.

open_system('sdoMassSpringDamper');

Загрузите данные эксперимента.

load sdoMassSpringDamper_ExperimentData

Переменные texp1, yexp1, texp2, и yexp2 загружаются в рабочую область. yexp1 и yexp2 описывать массовое положение для времени texp1 и texp2 соответственно.

Определение экспериментов по оценке

Создайте двухэлементный массив объектов эксперимента для хранения измеренных данных для двух экспериментов.

Создайте объект эксперимента для первого эксперимента.

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];

Создайте двухэлементный массив экспериментов. Поскольку два эксперимента идентичны, за исключением ожидаемых выходных данных, скопируйте первый эксперимент дважды.

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);

sdoMassSpringDamper_Objective функция:

Имеет один входной аргумент, определяющий массу, постоянную пружины и значение демпфера, а также исходное положение массы.

Имеет один входной аргумент, указывающий объект эксперимента, содержащий измеренные данные.

Возвращает вектор ошибок между моделируемыми и экспериментальными выходами.

sdoMassSpringDamper_Objective для функции требуется два входа, но sdo.optimize требуется функция с одним входным аргументом. Чтобы обойти это, estFcn является анонимной функцией с одним входным аргументом, v, но это вызывает sdoMassSpringDamper_Objective используя два входных аргумента, v и Exp.

Дополнительные сведения об анонимных функциях см. в разделе Анонимные функции.

sdo.optimize команда минимизирует возвращаемый аргумент анонимной функции estFcn, то есть остаточные ошибки, возвращенные sdoMassSpringDamper_Objective. Дополнительные сведения о том, как написать функцию цели/ограничения для использования с 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')

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