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

Мы начинаем наш тур по моделированию с загрузки данных о выходе (timeseries данных). Данные содержат один выход, который является угловым положением маятника [рад]. Угол равен нулю, когда маятник указывает вниз, и он увеличивается против часовой стрелки. Существует 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: Сравнение истинного результата и моделируемых выходов трех начальных моделей маятника.

Как видно, результат с методом Euler forward (с более мелкой сеткой, чем используемая по умолчанию) значительно отличается от результатов, полученных с решателями Runge-Kutta. В этом случае оказывается, что передний решатель Эйлера выдает довольно хороший результат (с точки зрения модели подгонки), тогда как выходы, полученные с решателями Рунге-Кутты, немного ограничены. Однако это несколько обманчиво, как будет очевидно позже, потому что начальное значение 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. Набор текста "помогите idnlgrey. SimulationOptions "предоставляет некоторые дополнительные детали о доступных решателях, а набор" help nlgreyestOptions "предоставляет детали о различных опциях оценки, которые могут использоваться для моделирования IDNLGREY.

Дополнительная информация

Для получения дополнительной информации об идентификации динамических систем с помощью System Identification Toolbox™ посетите информационную страницу System Identification Toolbox.