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

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

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

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

Функция 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 11-Jan-2019 03:25:58

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

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