exponenta event banner

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

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

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

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