В этом примере показано, как выбор алгоритма оценки может повлиять на результаты для нелинейной серой оценки модели поля. Мы используем данные, произведенные нелинейной системой маятника, которую схематично показывают в рисунке 1. Мы показываем в частности, как выбор решателя для дифференциальных уравнений влияет на результаты.
Рисунок 1: Схематическое представление классического маятника.
Мы запускаем наш тур моделирования путем загрузки выходных данных (данные timeseries). Данные содержат тот выход, y, который является угловым положением маятника [рад]. Угол является нулем, когда маятник указывает вниз, и это увеличивается против часовой стрелки. Существует 1 001 (симулированная) выборка точек данных, и шаг расчета составляет 0,1 секунды. Маятник затронут постоянной силой силы тяжести, но никакая другая внешняя сила не влияет на движение маятника. Чтобы исследовать эту ситуацию, мы создаем объект IDDATA:
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'pendulumdata')); z = iddata(y, [], 0.1, 'Name', 'Pendulum'); z.OutputName = 'Pendulum position'; z.OutputUnit = 'rad'; z.Tstart = 0; z.TimeUnit = 's';
Угловое положение маятника (выход) показывают в окне графика.
figure('Name', [z.Name ': output data']); plot(z);
Рисунок 2: Угловое положение маятника (выводится).
Следующий шаг должен задать структуру модели для системы маятника. Динамика его была изучена в многочисленных книгах и статьях и хорошо изучена. Если мы выбираем x1 (t) в качестве углового положения [рад] и x2 (t) как скорость вращения [rad/s] маятника, то это довольно прямо, чтобы настроить структуру пространства состояний следующего вида:
d/dt x1(t) = x2(t) d/dt x2(t) = -(g/l)*sin(x1(t)) - (b/(m*l^2))*x2(t)
y(t) = x1(t)
наличие параметров (или константы)
g - the gravity constant [m/s^2] l - the length of the rod of the pendulum [m] b - viscous friction coefficient [Nms/rad] m - the mass of the bob of the pendulum [kg]
Мы вводим эту информацию в файл MATLAB pendulum_m.m с содержимым можно следующим образом:
function [dx, y] = pendulum_m(t, x, u, g, l, b, m, varargin) %PENDULUM_M A pendulum system.
% Output equation. y = x(1); % Angular position.
% State equations.
dx = [x(2); ... % Angular position.
-(g/l)*sin(x(1))-b/(m*l^2)*x(2) ... % Angular velocity.
];
Следующий шаг должен создать типовой объект IDNLGREY для описания маятника. Мы также вводим информацию об именах и модулях входных параметров, состояний, выходных параметров и параметров. Вследствие физической действительности все параметры модели должны быть положительными, и это наложено путем установки свойства 'Minimum' всех параметров к наименьшему распознаваемому положительному значению в MATLAB®, eps (0).
FileName = 'pendulum_m'; % File describing the model structure. Order = [1 0 2]; % Model orders [ny nu nx]. Parameters = [9.81; 1; 0.2; 1]; % Initial parameters. InitialStates = [1; 0]; % Initial initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts); nlgr.OutputName = 'Pendulum position'; nlgr.OutputUnit = 'rad'; nlgr.TimeUnit = 's'; nlgr = setinit(nlgr, 'Name', {'Pendulum position' 'Pendulum velocity'}); nlgr = setinit(nlgr, 'Unit', {'rad' 'rad/s'}); nlgr = setpar(nlgr, 'Name', {'Gravity constant' 'Length' ... 'Friction coefficient' 'Mass'}); nlgr = setpar(nlgr, 'Unit', {'m/s^2' 'm' 'Nms/rad' 'kg'}); nlgr = setpar(nlgr, 'Minimum', {eps(0) eps(0) eps(0) eps(0)}); % All parameters > 0.
Прежде, чем попытаться оценить любой параметр мы симулируем систему с предполагаемыми значениями параметров. Мы делаем это для трех из доступных решателей, Эйлер вперед с фиксированной длиной шага (ode1), 23 Рунге-Кутта с адаптивной длиной шага (ode23), и 45 Рунге-Кутта с адаптивной длиной шага (ode45). Выходные параметры, полученные при использовании этих трех решателей, показывают в окне графика.
% A. Model computed with first-order Euler forward ODE solver. nlgref = nlgr; nlgref.SimulationOptions.Solver = 'ode1'; % Euler forward. nlgref.SimulationOptions.FixedStep = z.Ts*0.1; % Step size. % B. Model computed with adaptive Runge-Kutta 23 ODE solver. nlgrrk23 = nlgr; nlgrrk23.SimulationOptions.Solver = 'ode23'; % Runge-Kutta 23. % C. Model computed with adaptive Runge-Kutta 45 ODE solver. nlgrrk45 = nlgr; nlgrrk45.SimulationOptions.Solver = 'ode45'; % Runge-Kutta 45. compare(z, nlgref, nlgrrk23, nlgrrk45, 1, ... compareOptions('InitialCondition', 'model'));
Рисунок 3: Сравнение между истинным выходом и симулированными выходными параметрами трех первоначальных моделей маятника.
Как видно, результат с Эйлером передает метод (с более прекрасной сеткой, чем, что используется по умолчанию), значительно отличается от результатов, полученных с решателями Рунге-Кутта. В этом случае оказывается, что Эйлеров прямой решатель приводит к довольно хорошему результату (в терминах подгонки модели), тогда как выходные параметры, полученные с решателями Рунге-Кутта, немного ограничиваются. Однако это несколько обманывает, как будет очевидно позже, потому что начальное значение b, коэффициента вязкого трения, является вдвое более большим, чем в действительности.
Ускорение свободного падения, g, длина стержня, l, и масса боба, m, может легко быть измерено или взято из таблицы без оценки. Однако это обычно - не случай с коэффициентами трения, которые часто должны оцениваться. Мы поэтому фиксируем бывшие три параметра к g = 9.81, l = 1, и m = 1, и рассматриваем только b как свободный параметр:
nlgref = setpar(nlgref, 'Fixed', {true true false true}); nlgrrk23 = setpar(nlgrrk23, 'Fixed', {true true false true}); nlgrrk45 = setpar(nlgrrk45, 'Fixed', {true true false true});
Затем мы оцениваем b использование этих трех решателей для дифференциальных уравнений. Мы начинаем с Эйлеровой основанной на форварде структуры модели.
opt = nlgreyestOptions('Display', 'on'); tef = clock; nlgref = nlgreyest(z, nlgref, opt); % Perform parameter estimation. tef = etime(clock, tef);
Затем мы продолжаем основанное на модели на 23 решателях Рунге-Кутта (ode23).
trk23 = clock;
nlgrrk23 = nlgreyest(z, nlgrrk23, opt); % Perform parameter estimation.
trk23 = etime(clock, trk23);Мы наконец используем 45 решателей Рунге-Кутта (ode45).
trk45 = clock;
nlgrrk45 = nlgreyest(z, nlgrrk45, opt); % Perform parameter estimation.
trk45 = etime(clock, trk45);Результаты при использовании этих трех решателей получены в итоге ниже. Истинное значение b 0.1, который получен с ode45. ode23 возвращает значение вполне близко к 0,1, в то время как ode1 возвращает значение настоящее расстояние от 0,1.
disp(' Estimation time Estimated b value');Estimation time Estimated b value
fprintf(' ode1 : %3.1f %1.3f\n', tef, nlgref.Parameters(3).Value);ode1 : 6.8 0.194
fprintf(' ode23: %3.1f %1.3f\n', trk23, nlgrrk23.Parameters(3).Value);ode23: 1.7 0.093
fprintf(' ode45: %3.1f %1.3f\n', trk45, nlgrrk45.Parameters(3).Value);ode45: 2.5 0.100
Однако это не означает, что симулированные выходные параметры отличаются так очень от фактического выхода, потому что ошибки, произведенные различными решателями для дифференциальных уравнений, обычно составляются в процедуре оценки. Это - факт, который с готовностью виден путем симуляции трех предполагаемых моделей маятника.
compare(z, nlgref, nlgrrk23, nlgrrk45, 1, ... compareOptions('InitialCondition', 'model'));
Рисунок 4: Сравнение между истинным выходом и симулированными выходными параметрами трех предполагаемых моделей маятника.
Как ожидается данный этот рисунок, значения критерия Итоговой ошибки предсказания (FPE) этих моделей также скорее друг близко к другу:
fpe(nlgref, nlgrrk23, nlgrrk45);
1.0e-03 *
0.1609 0.1078 0.1056
На основе этого мы приходим к заключению, что все три модели могут получить динамику маятника, но Эйлерова основанная на форварде модель отражает базовую физику вполне плохо. Единственный способ увеличить его эффективность "физики" состоит в том, чтобы уменьшить длину шага решателя, но это также означает, что время решения значительно увеличивается. Наши эксперименты указывают, что 45 решателей Рунге-Кутта являются лучшим решателем для нежестких проблем при принятии во внимание и точности и вычислительной скорости.
45 (ode45) решателей Рунге-Кутта часто возвращают высококачественные результаты относительно быстро и поэтому выбраны в качестве решателя для дифференциальных уравнений по умолчанию в IDNLGREY. Ввод "справки idnlgrey.SimulationOptions" предоставляет еще некоторую подробную информацию о доступных решателях, и вводящий "справку nlgreyestOptions" предоставляет подробную информацию о различных опциях оценки, которые могут использоваться для моделирования IDNLGREY.
Для получения дополнительной информации об идентификации динамических систем с System Identification Toolbox™ посещают страницу информации о продукте System Identification Toolbox.