Выходные данные

Мы запускаем наш тур моделирования путем загрузки выходных данных (данные 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 :  4.8               0.194

fprintf('   ode23:  %3.1f               %1.3f\n', trk23, nlgrrk23.Parameters(3).Value);

   ode23:  1.9               0.093

fprintf('   ode45:  %3.1f               %1.3f\n', trk45, nlgrrk45.Parameters(3).Value);

   ode45:  3.0               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.