Оценка значений параметров модели (код)

Этот пример показывает, как использовать экспериментальные данные для оценки значений параметров модели.

Модель самолета

Модель Simulink, sdoAircraftEstimation, моделирует систему продольного управления полетом самолета.

open_system('sdoAircraftEstimation')

Задача оценки

Вы используете измеренные данные для оценки параметров и состояний модели самолета.

Измеренные выходные данные:

  • Пилотная сила G, выход Pilot G-force calculation блок

  • Угол атаки, четвертый выход Aircraft Dynamics Model блок

Параметры:

  • Постоянная времени привода, Ta, используемый в Actuator Model блок

  • Вертикальная скорость, Zd, используемый в Aircraft Dynamics Model блок

  • Коэффициент усиления тангажа, Md, используемый в Aircraft Dynamics Model блок

Состояние:

  • Начальное состояние модели привода первого порядка, sdoAircraftEstimation/Actuator Model

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

Получите измеренные данные.

[time,iodata] = sdoAircraftEstimation_Experiment;

The sdoAircraftEstimation_Experiment функция возвращает измеренные выходные данные, iodata, и соответствующий временной вектор. Первый столбец iodata - сила G пилота, а второй столбец - угол атаки.

Чтобы увидеть код для этой функции, введите edit sdoAircraftEstimation_Experiment.

Создайте объект эксперимента, чтобы сохранить измеренные входные/выходные данные.

Exp = sdo.Experiment('sdoAircraftEstimation');

Создайте объект, чтобы сохранить измеренный выход G-Force пилота.

PilotG = Simulink.SimulationData.Signal;
PilotG.Name      = 'PilotG';
PilotG.BlockPath = 'sdoAircraftEstimation/Pilot G-force calculation';
PilotG.PortType  = 'outport';
PilotG.PortIndex = 1;
PilotG.Values    = timeseries(iodata(:,2),time);

Создайте объект, чтобы сохранить измеренный угол атаки (альфа) на выходе.

AoA = Simulink.SimulationData.Signal;
AoA.Name = 'AngleOfAttack';
AoA.BlockPath = 'sdoAircraftEstimation/Aircraft Dynamics Model';
AoA.PortType  = 'outport';
AoA.PortIndex = 4;
AoA.Values    = timeseries(iodata(:,1),time);

Добавьте измеренную G-силу пилота и данные по углу атаки к эксперименту в качестве ожидаемых выходных данных.

Exp.OutputData = [...
    PilotG; ...
    AoA];

Добавьте начальное состояние для Actuator Model блок к эксперименту. Установите его Free поле к true так, чтобы оно было оценено.

Exp.InitialStates = sdo.getStateFromModel('sdoAircraftEstimation','Actuator Model');
Exp.InitialStates.Minimum = 0;
Exp.InitialStates.Free    = true;

Сравнение измеренного выхода и начального моделируемого выхода

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

Simulator    = createSimulator(Exp);
Simulator    = sim(Simulator);

Поиск G-силы пилота и угла сигналов атаки в записанных данных моделирования.

SimLog       = find(Simulator.LoggedData,get_param('sdoAircraftEstimation','SignalLoggingName'));
PilotGSignal = find(SimLog,'PilotG');
AoASignal    = find(SimLog,'AngleOfAttack');

Постройте график измеренных и моделируемых данных.

Как ожидалось, реакция модели не совпадает с экспериментальными выходными данными.

plot(time, iodata, ...
    AoASignal.Values.Time,AoASignal.Values.Data,'--', ...
    PilotGSignal.Values.Time,PilotGSignal.Values.Data,'-.');
title('Simulated and Measured Responses Before Estimation')
legend('Measured angle of attack',  'Measured pilot g force', ...
     'Simulated angle of attack', 'Simulated pilot g force');

Задайте параметры для оценки

Выберите параметры модели, которые описывают систему приведения в действие управления полетом. Задайте границы для предполагаемых значений параметров на основе нашего понимания системы приведения в действие.

p = sdo.getParameterFromModel('sdoAircraftEstimation',{'Ta','Md','Zd'});
p(1).Minimum = 0.01;  %Ta
p(1).Maximum = 1;
p(2).Minimum = -10;   %Md
p(2).Maximum = 0;
p(3).Minimum = -100;  %Zd
p(3).Maximum = 0;

Получите начальное значение состояния привода, которое должно быть оценено из эксперимента.

s = getValuesToEstimate(Exp);

Сгруппируйте параметры модели и начальные состояния, которые будут оценены вместе.

v = [p;s]
 
v(1,1) =
 
       Name: 'Ta'
      Value: 0.5000
    Minimum: 0.0100
    Maximum: 1
       Free: 1
      Scale: 0.5000
       Info: [1x1 struct]

 
v(2,1) =
 
       Name: 'Md'
      Value: -1
    Minimum: -10
    Maximum: 0
       Free: 1
      Scale: 1
       Info: [1x1 struct]

 
v(3,1) =
 
       Name: 'Zd'
      Value: -80
    Minimum: -100
    Maximum: 0
       Free: 1
      Scale: 128
       Info: [1x1 struct]

 
v(4,1) =
 
       Name: 'sdoAircraftEstimation/Actuator...'
      Value: 0
    Minimum: 0
    Maximum: Inf
       Free: 1
      Scale: 1
    dxValue: 0
     dxFree: 1
       Info: [1x1 struct]

 
4x1 param.Continuous
 

Определите целевую функцию оценки

Создайте целевую функцию оценки, чтобы оценить, насколько близко выход симуляции, сгенерированный с использованием оцененных значений параметров, совпадает с измеренными данными.

Используйте анонимную функцию с одним входным параметром, который вызывает sdoAircraftEstimation_Objective функция. Передаем анонимную функцию в sdo.optimize, который оценивает функцию при каждой итерации оптимизации.

estFcn = @(v) sdoAircraftEstimation_Objective(v,Simulator,Exp);

The sdoAircraftEstimation_Objective функция:

  • Имеет один входной параметр, который задает значения параметров привода и начальное состояние привода.

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

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

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

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

The sdo.optimize команда минимизирует возвращаемый аргумент анонимной функции estFcn, то есть остаточные ошибки, возвращенные sdoAircraftEstimation_Objective. Для получения дополнительной информации о том, как написать функцию objective/constraint для использования со sdo.optimize команда, введите help sdoExampleCostFunction в командной строке MATLAB.

Чтобы изучить целевую функцию оценки более подробно, введите edit sdoAircraftEstimation_Objective в командной строке MATLAB.

type sdoAircraftEstimation_Objective
function vals = sdoAircraftEstimation_Objective(v,Simulator,Exp) 
%SDOAIRCRAFTESTIMATION_OBJECTIVE
%
%    The sdoAircraftEstimation_Objective function is used to compare model
%    outputs against experimental data.
%
%    vals = sdoAircraftEstimation_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,
%    sdoAircraftEstimation_cmddemo
%
 
% 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.
%
Simulator = createSimulator(Exp,Simulator);
Simulator = sim(Simulator);

SimLog       = find(Simulator.LoggedData,get_param('sdoAircraftEstimation','SignalLoggingName'));
PilotGSignal = find(SimLog,'PilotG');
AoASignal    = find(SimLog,'AngleOfAttack');

PilotGError = evalRequirement(r,PilotGSignal.Values,Exp.OutputData(1).Values);
AoAError    = evalRequirement(r,AoASignal.Values,Exp.OutputData(2).Values);

%%
% Return the residual errors to the optimization solver.
%
vals.F = [PilotGError(:); AoAError(:)];
end

Оцените параметры

Используйте sdo.optimize функция для оценки значений параметров привода и начального состояния.

Задайте опции оптимизации. Функция оценки sdoAircraftEstimation_Objective возвращает невязки ошибок между моделируемыми и экспериментальными данными и не включает никаких ограничений, что делает эту задачу идеальной для решателя 'lsqnonlin'.

opt = sdo.OptimizeOptions;
opt.Method = 'lsqnonlin';

Оцените параметры.

vOpt = sdo.optimize(estFcn,v,opt)
 Optimization started 27-Jan-2021 14:00:21

                                          First-order 
 Iter F-count        f(x)      Step-size  optimality
    0      8      27955.2            1                                         
    1     17      10121.6       0.4744     5.68e+04
    2     26       3127.5       0.3845     1.24e+04
    3     35      872.577       0.4293     2.81e+03
    4     44      238.548       0.5156          618
    5     53      71.5731       0.4939          147
    6     62      16.9301       0.4222         44.3
    7     71      1.82188       0.3019         11.5
    8     80    0.0426814       0.1356         1.35
    9     89   0.00113417      0.02555        0.243
   10     98  0.000247412     0.008296      0.00764
Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.
 
vOpt(1,1) =
 
       Name: 'Ta'
      Value: 0.0500
    Minimum: 0.0100
    Maximum: 1
       Free: 1
      Scale: 0.5000
       Info: [1x1 struct]

 
vOpt(2,1) =
 
       Name: 'Md'
      Value: -6.8849
    Minimum: -10
    Maximum: 0
       Free: 1
      Scale: 1
       Info: [1x1 struct]

 
vOpt(3,1) =
 
       Name: 'Zd'
      Value: -63.9989
    Minimum: -100
    Maximum: 0
       Free: 1
      Scale: 128
       Info: [1x1 struct]

 
vOpt(4,1) =
 
       Name: 'sdoAircraftEstimation/Actuator...'
      Value: 1.3254e-04
    Minimum: 0
    Maximum: Inf
       Free: 1
      Scale: 1
    dxValue: 0
     dxFree: 1
       Info: [1x1 struct]

 
4x1 param.Continuous
 

Сравнение измеренного выхода и итогового моделируемого выхода

Обновите эксперименты с оцененными значениями параметров.

Exp = setEstimatedValues(Exp,vOpt);

Симулируйте модель с помощью обновленного эксперимента и сравните моделируемый выход с экспериментальными данными.

Реакция модели с использованием оцененных значений параметров тесно совпадает с выходными данными эксперимента.

Simulator    = createSimulator(Exp,Simulator);
Simulator    = sim(Simulator);
SimLog       = find(Simulator.LoggedData,get_param('sdoAircraftEstimation','SignalLoggingName'));
PilotGSignal = find(SimLog,'PilotG');
AoASignal    = find(SimLog,'AngleOfAttack');

plot(time, iodata, ...
    AoASignal.Values.Time,AoASignal.Values.Data,'-.', ...
    PilotGSignal.Values.Time,PilotGSignal.Values.Data,'--')
title('Simulated and Measured Responses After Estimation')
legend('Measured angle of attack',  'Measured pilot g force', ...
    'Simulated angle of attack', 'Simulated pilot g force');

Обновите значения параметров модели

Обновите модель с оцененными значениями параметров привода. Не обновляйте начальное состояние привода модели (четвертый элемент vOpt), так как он зависит от эксперимента.

sdo.setValueInModel('sdoAircraftEstimation',vOpt(1:3));

Похожие примеры

Чтобы узнать, как оценить параметры модели с помощью приложения Parameter Estimator, смотрите Estimate Model Parameter Values (GUI).

Закройте модель.

bdclose('sdoAircraftEstimation')