Этот пример показывает, как выбор алгоритма оценки может повлиять на результаты для нелинейной оценки модели серого прямоугольника. Мы используем данные, полученные с помощью нелинейной маятниковой системы, которая схематически показана на рис. 1. В частности, мы покажем, как выбор решателя дифференциальных уравнений влияет на результаты.
Рис. 1: Схематический вид классического маятника.
Мы начинаем наше моделирующее путешествие с загрузки выходных данных (данных временных рядов). Данные содержат один выход, y, который является угловым положением маятника [рад]. Угол равен нулю, когда маятник направлен вниз и увеличивается против часовой стрелки. Имеется 1001 (смоделированная) выборка точек данных, и время выборки составляет 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) в качестве углового положения [rad] и 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]
Мы вводим эту информацию в pendulum_m.m файла MATLAB со следующим содержимым:
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.
Перед попыткой оценить какой-либо параметр мы моделируем систему с предполагаемыми значениями параметра. Мы делаем это для трех доступных решателей, Euler вперед с фиксированной длиной шага (ode1), Runge-Kutta 23 с адаптивной длиной шага (ode23) и Runge-Kutta 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 с помощью трех дифференциальных решателей уравнений. Мы начинаем со структуры модели на основе Euler.
opt = nlgreyestOptions('Display', 'on'); tef = clock; nlgref = nlgreyest(z, nlgref, opt); % Perform parameter estimation. tef = etime(clock, tef);
Затем мы продолжаем модель на основе решателя Runge-Kutta 23 (ode23).
trk23 = clock;
nlgrrk23 = nlgreyest(z, nlgrrk23, opt); % Perform parameter estimation.
trk23 = etime(clock, trk23);Наконец, мы используем решатель Runge-Kutta 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
На основании этого мы делаем вывод, что все три модели способны улавливать динамику маятника, но модель, основанная на Эйлере вперёд, отражает лежащую в основе физику довольно плохо. Единственный способ увеличить его «физические» характеристики - уменьшить длину шага решателя, но это также означает, что время решения значительно увеличивается. Наши эксперименты показывают, что решатель Runge-Kutta 45 является лучшим решателем для нежестких задач при учете как точности, так и вычислительной скорости.
Решатель Runge-Kutta 45 (ode45) часто возвращает результаты высокого качества относительно быстро и поэтому выбирается в качестве решателя дифференциальных уравнений по умолчанию в IDNLGREY. Ввод "help idnlgrey. В разделе «Параметры» содержатся дополнительные сведения о доступных решателях, а при вводе «help nlgreyestOptions» - сведения о различных вариантах оценки, которые можно использовать для моделирования IDNLGREY.
Дополнительные сведения об идентификации динамических систем с помощью Toolbox™ System Identification см. на информационной странице System Identification Toolbox.