Моделирование системы динамики аппарата

Этот пример показывает нелинейное моделирование серого ящика динамики аппарата. Много новых функций транспортного средства (как Электронные программы устойчивости (ESP), косвенные Системы мониторинга давления воздуха в шине (TPMS), системы мониторинга трения дорожной шины, и т.д) используют модели базовой динамики аппарата. Так называемая велосипедная модель транспортного средства является довольно простой структурой модели, которая часто используется в литературе динамики аппарата. В этом примере мы начнемся с этой структурой модели и попыткой оценить продольное и боковую жесткость шины. Фактическая работа моделирования была первоначально выполнена Эриком Нарби в его магистре наук, работают в AB Динамики NIRA, Швеция.

Моделирование динамики аппарата

Следующая фигура иллюстрирует ситуацию с моделированием транспортного средства, которая будет рассмотрена.

Рисунок 1: Схематическое представление системы динамики аппарата.

При помощи закона Ньютона движения и некоторых основных геометрических отношений, продольная скорость v_x (t), поперечная скорость v_y (t) и уровень рыскания r (t) измеренный вокруг Центра тяжести (COG) транспортного средства может быть описана следующими тремя дифференциальными уравнениями:

  d
  -- v_x(t) = v_y(t)*r(t) + 1/m*(  (F_x,FL(t)+F_x,FR(t))*cos(delta(t))
  dt                             - (F_y,FL(t)+F_y,FR(t))*sin(delta(t))
                                 + F_x,RL(t)+F_x,RR(t)
                                 - C_A*v_x(t)^2)
  d
  -- v_y(t) = -v_x(t)*r(t) + 1/m*(  (F_x,FL(t)+F_x,FR(t))*sin(delta(t))
  dt                              + (F_y,FL(t)+F_y,FR(t))*cos(delta(t))
                                  + F_y,RL(t)+F_y,RR(t))
  d
  -- r(t)   = 1/J*(  a*(  (F_x,FL(t)+F_x,FR(t))*sin(delta(t))
  dt                    + (F_y,FL(t)+F_y,FR(t))*cos(delta(t)))
                   - b*(F_y,RL(t)+F_y,RR(t)))

где индекс x используется, чтобы обозначить, что сила F действует в продольном направлении и y, что это действует в боковом направлении. Сокращения FL, FR, метка RL и RR шины: Переднее Левое, Переднее Право, Заднее Левое и Заднее Право, соответственно. Первое уравнение, описывающее продольное ускорение также, содержит термин сопротивления воздуха, который принят, чтобы быть квадратичной функцией продольной скорости транспортного средства v_x (t). Кроме того, дельта (t) (вход) является держащимся углом, J момент инерции, и a и b расстояния от центра тяжести до передних и задних осей, соответственно.

Давайте примем, что силы шины могут быть смоделированы посредством следующих линейных аппроксимаций:

  F_x,i(t) = C_x*s_i(t)
  F_y,i(t) = C_y*alpha_i(t)     for i = {FL, FR, RL, RR}

где C_x и C_y являются продольной и боковой жесткостью шины, соответственно. Здесь мы приняли, что эти параметры жесткости являются тем же самым для всех 4 шин. s_i (t) является так называемым (продольным) промахом шины i и alpha_i (t) угол промаха шины. Для управляемого транспортного средства с передними ведущими колесами (как рассмотрено здесь), промахи s_FL (t) и s_FR (t) выведены из отдельных скоростей колеса (измеренных) путем предположения, что задние колеса не показывают промаха (i.e., s_RL (t) = s_RR (t) = 0). Следовательно промахи являются входными параметрами к нашей структуре модели. Для передних колес углы промаха шины alpha_Fj (t) могут быть аппроксимированы (когда v_x (t)> 0)

alpha_Fj(t) = delta(t) - arctan((v_y(t) + a*r(t))/v_x(t))
            ~ delta(t) - (v_y(t) + a*r(t))/v_x(t)        for j = {L, R}

Для задних колес углы промаха шины alpha_Rj (t) так же выведены и вычислены как

alpha_Rj(t) = - arctan((v_y(t) - b*r(t))/v_x(t))
            ~ - (v_y(t) - b*r(t))/v_x(t)                 for j = {L, R}

С J = 1 / ((0.5* (a+b)) ^2*m), мы можем затем настроить структуру пространства состояний, описывающую динамику аппарата. Введите состояния:

  x1(t) = v_x(t)      Longitudinal velocity [m/s].
  x2(t) = v_y(t)      Lateral velocity [m/s].
  x3(t) = r(t)        Yaw rate [rad/s].

эти пять измеренных или выведенных входных сигналов

  u1(t) = s_FL(t)     Slip of Front Left tire [ratio].
  u2(t) = s_FR(t)     Slip of Front Right tire [ratio].
  u3(t) = s_RL(t)     Slip of Rear Left tire [ratio].
  u4(t) = s_RR(t)     Slip of Rear Right tire [ratio].
  u5(t) = delta(t)    Steering angle [rad].

и параметры модели:

  m                   Mass of the vehicle [kg].
  a                   Distance from front axle to COG [m].
  b                   Distance from rear axle to COG [m].
  Cx                  Longitudinal tire stiffness [N].
  Cy                  Lateral tire stiffness [N/rad].
  CA                  Air resistance coefficient [1/m].

Выходные параметры системы являются продольной скоростью транспортного средства y1 (t) = x1 (t), боковое ускорение транспортного средства (измеренный акселерометром):

y2(t) = a_y(t) = 1/m*(  (F_x,FL(t) + F_x,FR(t))*sin(delta(t))
                      + (F_y,FL(t) + F_y,FR(t))*cos(delta(t))
                      + F_y,RL(t) + F_y,RR(t))

и уровень рыскания y3 (t) = r (t) (измеренный гироскопом).

Соедините, мы прибываем в следующую структуру модели в пространстве состояний:

  d
  -- x1(t) = x2(t)*x3(t) + 1/m*(  Cx*(u1(t)+u2(t))*cos(u5(t))
  dt                            - 2*Cy*(u5(t)-(x2(t)+a*x3(t))/x1(t))*sin(u5(t))
                                + Cx*(u3(t)+u4(t))
                                - CA*x1(t)^2)
  d
  -- x2(t) = -x1(t)*x3(t) + 1/m*(  Cx*(u1(t)+u2(t))*sin(u5(t))
  dt                             + 2*Cy*(u5(t)-(x2(t)+a*x3(t))/x1(t))*cos(u5(t))
                                 + 2*Cy*(b*x3(t)-x2(t))/x1(t))
  d
  -- x3(t) = 1/((0.5*(a+b))^2)*m)*(  a*(  Cx*(u1(t)+u2(t)*sin(u5(t))
  dt                                    + 2*Cy*(u5(t) - (x2(t)+a*x3(t))/x1(t))*cos(u5(t)))
                                   - 2*b*Cy*(b*x3(t)-x2(t))/x1(t))
     y1(t) = x1(t)
     y2(t) = 1/m*(  Cx*(u1(t)+u2(t))*sin(u5(t))
                  + 2*Cy*(u5(t)-(x2(t)+a*x3(t))/x1(t))*cos(u5(t))
                  + 2*Cy*(b*x3(t)-x2(t))/x1(t))
     y3(t) = x3(t)

Модель транспортного средства IDNLGREY

Когда базис для нашей идентификации транспортного средства экспериментирует, мы сначала должны создать файл модели IDNLGREY, описывающий эти уравнения транспортного средства. Здесь мы используем моделирование C-MEX и создаем vehicle_c.c файл модели, в котором Нью-Йорк установлен в 3. Состояние и выходные функции обновления vehicle_c.c, compute_dx и compute_y, несколько включены, и включает несколько стандартных математических функций C-defined, как потому что(.) и sin (.), а также голова (.) для вычисления степени его аргумента.

Функция обновления состояния compute_dx возвращает дуплекс (аргумент 1) и использует 3 входных параметра: вектор состояния x, входной вектор u и шесть скалярных параметров, закодированных в p (t и auxvar файла модели шаблона C-MEX были удалены сюда):

  /* State equations. */
  void compute_dx(double *dx, double *x, double *u, double **p)
  {
      /* Retrieve model parameters. */
      double *m, *a, *b, *Cx, *Cy, *CA;
      m  = p[0];   /* Vehicle mass.                    */
      a  = p[1];   /* Distance from front axle to COG. */
      b  = p[2];   /* Distance from rear axle to COG.  */
      Cx = p[3];   /* Longitudinal tire stiffness.     */
      Cy = p[4];   /* Lateral tire stiffness.          */
      CA = p[5];   /* Air resistance coefficient.      */
      /* x[0]: Longitudinal vehicle velocity. */
      /* x[1]: Lateral vehicle velocity. */
      /* x[2]: Yaw rate. */
      dx[0] = x[1]*x[2]+1/m[0]*(Cx[0]*(u[0]+u[1])*cos(u[4])
              -2*Cy[0]*(u[4]-(x[1]+a[0]*x[2])/x[0])*sin(u[4])
              +Cx[0]*(u[2]+u[3])-CA[0]*pow(x[0],2));
      dx[1] = -x[0]*x[2]+1/m[0]*(Cx[0]*(u[0]+u[1])*sin(u[4])
              +2*Cy[0]*(u[4]-(x[1]+a[0]*x[2])/x[0])*cos(u[4])
              +2*Cy[0]*(b[0]*x[2]-x[1])/x[0]);
      dx[2] = 1/(pow(((a[0]+b[0])/2),2)*m[0])
              *(a[0]*(Cx[0]*(u[0]+u[1])*sin(u[4])
              +2*Cy[0]*(u[4]-(x[1]+a[0]*x[2])/x[0])*cos(u[4]))
              -2*b[0]*Cy[0]*(b[0]*x[2]-x[1])/x[0]);
  }

Выходная функция обновления compute_y возвращает y (аргумент 1) и использует 3 входных параметра: вектор состояния x, входной вектор u и пять из этих шести параметров (CA сопротивления воздуха не нужен), закодированный в p:

  /* Output equations. */
  void compute_y(double *y, double *x, double *u, double **p)
  {
      /* Retrieve model parameters. */
      double *m  = p[0];   /* Vehicle mass.                    */
      double *a  = p[1];   /* Distance from front axle to COG. */
      double *b  = p[2];   /* Distance from rear axle to COG.  */
      double *Cx = p[3];   /* Longitudinal tire stiffness.     */
      double *Cy = p[4];   /* Lateral tire stiffness.          */
      /* y[0]: Longitudinal vehicle velocity. */
      /* y[1]: Lateral vehicle acceleration. */
      /* y[2]: Yaw rate. */
      y[0] = x[0];
      y[1] = 1/m[0]*(Cx[0]*(u[0]+u[1])*sin(u[4])
             +2*Cy[0]*(u[4]-(x[1]+a[0]*x[2])/x[0])*cos(u[4])
             +2*Cy[0]*(b[0]*x[2]-x[1])/x[0]);
      y[2] = x[2];
  }

Имея соответствующий файл структуры модели, следующий шаг должен создать объект IDNLGREY, отражающий ситуацию с моделированием. Для простоты бухгалтерии мы также задаем имена и модули вводов и выводов.

FileName      = 'vehicle_c';                        % File describing the model structure.
Order         = [3 5 3];                            % Model orders [ny nx nu].
Parameters    = [1700; 1.5; 1.5; 1.5e5; 4e4; 0.5];  % Initial parameters.
InitialStates = [1; 0; 0];                          % Initial value of initial states.
Ts            = 0;                                  % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ...
                'Name', 'Bicycle vehicle model', 'TimeUnit', 's');
        nlgr.InputName =  {'Slip on front left tire';               ...   % u(1).
                         'Slip on front right tire';              ...   % u(2).
                         'Slip on rear left tire';                ...   % u(3).
                         'Slip on rear right tire';               ...   % u(4).
                         'Steering angle'};                       ...   % u(5).
          nlgr.InputUnit =  {'ratio'; 'ratio'; 'ratio'; 'ratio'; 'rad'};

          nlgr.OutputName = {'Long. velocity';  ...   % y(1); Longitudinal vehicle velocity
                         'Lat. accel.';   ...     % y(2); Lateral vehicle acceleration
                         'Yaw rate'};                             ...   % y(3).
          nlgr.OutputUnit = {'m/s'; 'm/s^2'; 'rad/s'};

Имена и модули (начальных) состояний и параметров модели заданы через SETINIT. Мы также используем эту команду, чтобы указать, что первое начальное состояние (продольная скорость) должно быть строго положительно для модели быть допустимым и указать, что все параметры модели должны быть строго положительными. Эти ограничения будут впоследствии соблюдаться при выполнении оценки параметра модели и/или начального состояния.

nlgr = setinit(nlgr, 'Name', {'Longitudinal vehicle velocity'        ... % x(1).
                       'Lateral vehicle velocity'             ... % x(2).
                       'Yaw rate'});                          ... % x(3).
nlgr = setinit(nlgr, 'Unit', {'m/s'; 'm/s'; 'rad/s'});
nlgr.InitialStates(1).Minimum = eps(0);   % Longitudinal velocity > 0 for the model to be valid.
nlgr = setpar(nlgr, 'Name', {'Vehicle mass';                         ... % m.
                      'Distance from front axle to COG';      ... % a
                      'Distance from rear axle to COG';       ... % b.
                      'Longitudinal tire stiffness';          ... % Cx.
                      'Lateral tire stiffness';               ... % Cy.
                      'Air resistance coefficient'});         ... % CA.
nlgr = setpar(nlgr, 'Unit', {'kg'; 'm'; 'm'; 'N'; 'N/rad'; '1/m'});
nlgr = setpar(nlgr, 'Minimum', num2cell(eps(0)*ones(6, 1)));   % All parameters > 0!

Четыре из шести параметров этой структуры модели могут с готовностью быть получены через таблицу данных рассматриваемого транспортного средства:

  m  = 1700 kg
  a  = 1.5 m
  b  = 1.5 m
  CA = 0.5 or 0.7 1/m (see below)

Следовательно мы не оценим эти параметры:

nlgr.Parameters(1).Fixed = true;
nlgr.Parameters(2).Fixed = true;
nlgr.Parameters(3).Fixed = true;
nlgr.Parameters(6).Fixed = true;

С этим текстовые сводные данные вводимой структуры модели IDNLGREY получены через PRESENT можно следующим образом.

present(nlgr);
                                                                                          
nlgr =                                                                                    
Continuous-time nonlinear grey-box model defined by 'vehicle_c' (MEX-file):               
                                                                                          
   dx/dt = F(t, u(t), x(t), p1, ..., p6)                                                  
    y(t) = H(t, u(t), x(t), p1, ..., p6) + e(t)                                           
                                                                                          
 with 5 input(s), 3 state(s), 3 output(s), and 2 free parameter(s) (out of 6).            
                                                                                          
 Inputs:                                                                                  
    u(1)  Slip on front left tire(t) [ratio]                                              
    u(2)  Slip on front right tire(t) [ratio]                                             
    u(3)  Slip on rear left tire(t) [ratio]                                               
    u(4)  Slip on rear right tire(t) [ratio]                                              
    u(5)  Steering angle(t) [rad]                                                         
 States:                                         Initial value                            
    x(1)  Longitudinal vehicle velocity(t) [m/s]   xinit@exp1   1   (fixed) in ]0, Inf]   
    x(2)  Lateral vehicle velocity(t) [m/s]        xinit@exp1   0   (fixed) in [-Inf, Inf]
    x(3)  Yaw rate(t) [rad/s]                      xinit@exp1   0   (fixed) in [-Inf, Inf]
 Outputs:                                                                                 
    y(1)  Long. velocity(t) [m/s]                                                         
    y(2)  Lat. accel.(t) [m/s^2]                                                          
    y(3)  Yaw rate(t) [rad/s]                                                             
 Parameters:                                 Value                                        
    p1   Vehicle mass [kg]                         1700      (fixed) in ]0, Inf]          
    p2   Distance from front axle to COG [m]        1.5      (fixed) in ]0, Inf]          
    p3   Distance from rear axle to COG [m]         1.5      (fixed) in ]0, Inf]          
    p4   Longitudinal tire stiffness [N]         150000      (estimated) in ]0, Inf]      
    p5   Lateral tire stiffness [N/rad]           40000      (estimated) in ]0, Inf]      
    p6   Air resistance coefficient [1/m]           0.5      (fixed) in ]0, Inf]          
                                                                                          
Name: Bicycle vehicle model                                                               
                                                                                          
Status:                                                                                   
Created by direct construction or transformation. Not estimated.                          
More information in model's "Report" property.                                            

Данные ввода - вывода

На данном этапе мы загружаем доступные данные ввода - вывода. Этот файл содержит данные из трех различных экспериментов:

  A. Simulated data with high stiffness tires [y1 u1].
  B. Simulated data with low stiffness tires  [y2 u2].
  C. Measured data from a Volvo V70           [y3 u3].

Во всех случаях, шаг расчета Ts = 0,1 секунды.

load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'vehicledata'));

A. System Identification Используя симулированные высокие данные о жесткости шины

В нашем первом транспортном средстве экспериментирует идентификация, мы считаем симулированными, высоко утомляют данные о жесткости. Копия структуры модели nlgr и объекта IDDATA z1 отражение этой конкретной ситуации с моделированием сначала создается. Эти 5 входных сигналов хранятся в u1 и этих 3 выходных сигналах в y1. Входные параметры промаха (сгенерированный от сигналов скорости колеса) для передних колес были выбраны, чтобы быть синусоидальными с постоянным смещением; уровень рыскания был также синусоидальным, но с различной амплитудой и частотой. В действительности это - несколько искусственная ситуация, потому что каждый редко волнует транспортное средство так в боковом направлении.

nlgr1 = nlgr;
nlgr1.Name = 'Bicycle vehicle model with high tire stiffness';
z1 = iddata(y1, u1, 0.1, 'Name', 'Simulated high tire stiffness vehicle data');
z1.InputName = nlgr1.InputName;
z1.InputUnit = nlgr1.InputUnit;
z1.OutputName = nlgr1.OutputName;
z1.OutputUnit = nlgr1.OutputUnit;
z1.Tstart = 0;
z1.TimeUnit = 's';

Вводы и выводы показывают на двух рисунках графика.

h_gcf = gcf;
set(h_gcf,'DefaultLegendLocation','southeast');
h_gcf.Position = [100 100 795 634];
for i = 1:z1.Nu
   subplot(z1.Nu, 1, i);
   plot(z1.SamplingInstants, z1.InputData(:,i));
   title(['Input #' num2str(i) ': ' z1.InputName{i}]);
   xlabel('');
   axis tight;
end
xlabel([z1.Domain ' (' z1.TimeUnit ')']);

Рисунок 2: Входные параметры к системе транспортного средства с высокой жесткостью шины.

for i = 1:z1.Ny
   subplot(z1.Ny, 1, i);
   plot(z1.SamplingInstants, z1.OutputData(:,i));
   title(['Output #' num2str(i) ': ' z1.OutputName{i}]);
   xlabel('');
   axis tight;
end
xlabel([z1.Domain ' (' z1.TimeUnit ')']);

Рисунок 3: Выходные параметры от системы транспортного средства с высокой жесткостью шины.

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

clf
compare(z1, nlgr1, [], compareOptions('InitialCondition', 'model'));

Рисунок 4: Сравнение между истинными выходными параметрами и симулированными выходными параметрами первоначальной модели транспортного средства с высокой жесткостью шины.

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

nlgr1 = nlgreyest(z1, nlgr1);

Сравнение между истиной и симулированными выходными параметрами (с предполагаемой моделью) показывают в окне графика.

compare(z1, nlgr1, [], compareOptions('InitialCondition', 'model'));

Рисунок 5: Сравнение между истинными выходными параметрами и симулированными выходными параметрами предполагаемой модели транспортного средства с высокой жесткостью шины.

Эффективность симуляции предполагаемой модели довольно хороша. Предполагаемые параметры жесткости также близко к тем используемым в Simulink®, чтобы сгенерировать истинные выходные данные:

disp('                        True      Estimated');
fprintf('Longitudinal stiffness: %6.0f    %6.0f\n', 2e5, nlgr1.Parameters(4).Value);
fprintf('Lateral stiffness     : %6.0f    %6.0f\n', 5e4, nlgr1.Parameters(5).Value);
                        True      Estimated
Longitudinal stiffness: 200000    198517
Lateral stiffness     :  50000     53752

B. System Identification Используя симулированные низкие данные о жесткости шины

Во втором эксперименте мы повторяем моделирование из первого эксперимента, но теперь с симулированными низкими данными о жесткости шины.

nlgr2 = nlgr;
nlgr2.Name = 'Bicycle vehicle model with low tire stiffness';
z2 = iddata(y2, u2, 0.1, 'Name', 'Simulated low tire stiffness vehicle data');
z2.InputName = nlgr2.InputName;
z2.InputUnit = nlgr2.InputUnit;
z2.OutputName = nlgr2.OutputName;
z2.OutputUnit = nlgr2.OutputUnit;
z2.Tstart = 0;
z2.TimeUnit =  's';

Вводы и выводы показывают на двух рисунках графика.

clf
for i = 1:z2.Nu
   subplot(z2.Nu, 1, i);
   plot(z2.SamplingInstants, z2.InputData(:,i));
   title(['Input #' num2str(i) ': ' z2.InputName{i}]);
   xlabel('');
   axis tight;
end
xlabel([z2.Domain ' (' z2.TimeUnit ')']);

Рисунок 6: Входные параметры к системе транспортного средства с низкой жесткостью шины.

clf
for i = 1:z2.Ny
   subplot(z2.Ny, 1, i);
   plot(z2.SamplingInstants, z2.OutputData(:,i));
   title(['Output #' num2str(i) ': ' z2.OutputName{i}]);
   xlabel('');
   axis tight;
end
xlabel([z2.Domain ' (' z2.TimeUnit ')']);

Рисунок 7: Выходные параметры от системы транспортного средства с низкой жесткостью шины.

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

clf
compare(z2, nlgr2, [], compareOptions('InitialCondition', 'model'));

Рисунок 8: Сравнение между истинными выходными параметрами и симулированными выходными параметрами первоначальной модели транспортного средства с низкой жесткостью шины.

Два параметра жесткости затем оцениваются.

nlgr2 = nlgreyest(z2, nlgr2);

Сравнение между истиной и симулированными выходными параметрами (с предполагаемой моделью) показывают в окне графика.

compare(z2, nlgr2, [], compareOptions('InitialCondition', 'model'));

Рисунок 9: Сравнение между истинными выходными параметрами и симулированными выходными параметрами предполагаемой модели транспортного средства с низкой жесткостью шины.

Эффективность симуляции предполагаемой модели снова действительно хороша. Даже с той же начальной точкой параметра, как использовался в высоком случае жесткости шины, предполагаемые параметры жесткости также здесь близко к тем используемым в Simulink, чтобы сгенерировать истинные выходные данные:

disp('                        True      Estimated');
fprintf('Longitudinal stiffness: %6.0f    %6.0f\n', 1e5, nlgr2.Parameters(4).Value);
fprintf('Lateral stiffness     : %6.0f    %6.0f\n', 2.5e4, nlgr2.Parameters(5).Value);
                        True      Estimated
Longitudinal stiffness: 100000     99573
Lateral stiffness     :  25000     26117

C. System Identification Используя измеренную Volvo V70 Data

В итоговом эксперименте мы считаем данные собранными в Volvo V70. Как выше, мы делаем копию типового объекта модели транспортного средства nlgr и создаем новый объект IDDATA, содержащий результаты измерений. Здесь мы также увеличили коэффициент сопротивления воздуха с 0,50 до 0,70, чтобы лучше отразить ситуацию Volvo V70.

nlgr3 = nlgr;
nlgr3.Name = 'Volvo V70 vehicle model';
nlgr3.Parameters(6).Value = 0.70;   % Use another initial CA for the Volvo data.
z3 = iddata(y3, u3, 0.1, 'Name', 'Volvo V70 data');
z3.InputName = nlgr3.InputName;
z3.InputUnit = nlgr3.InputUnit;
z3.OutputName = nlgr3.OutputName;
z3.OutputUnit = nlgr3.OutputUnit;
z3.Tstart = 0;
z3.TimeUnit = 's';

Вводы и выводы показывают на двух рисунках графика. Как видно, результаты измерений являются довольно шумными.

clf
for i = 1:z3.Nu
   subplot(z3.Nu, 1, i);
   plot(z3.SamplingInstants, z3.InputData(:,i));
   title(['Input #' num2str(i) ': ' z3.InputName{i}]);
   xlabel('');
   axis tight;
end
xlabel([z3.Domain ' (' z3.TimeUnit ')']);

Рисунок 10: Измеренные входные параметры от транспортного средства Volvo V70.

clf
for i = 1:z3.Ny
   subplot(z3.Ny, 1, i);
   plot(z3.SamplingInstants, z3.OutputData(:,i));
   title(['Output #' num2str(i) ': ' z3.OutputName{i}]);
   xlabel('');
   axis tight;
end
xlabel([z3.Domain ' (' z3.TimeUnit ')']);

Рисунок 11: Измеренные выходные параметры от транспортного средства Volvo V70.

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

nlgr3 = setinit(nlgr3, 'Value', {18.7; 0; 0});   % Initial value of initial states.
clf
compare(z3, nlgr3);

Рисунок 12: Сравнение между измеренными выходными параметрами и симулированными выходными параметрами первоначальной модели транспортного средства Volvo V70.

Параметры жесткости шины Ккс и Сай затем оцениваются, в этом случае с помощью метода поиска Levenberg-Marquardt, после чего новая симуляция с предполагаемой моделью, выполняются. Кроме того, мы здесь оцениваем начальное значение продольной скорости, тогда как начальные значения поперечной скорости и уровня рыскания сохранены фиксированными.

nlgr3 = setinit(nlgr3, 'Fixed', {false; true; true});
nlgr3 = nlgreyest(z3, nlgr3, nlgreyestOptions('SearchMethod', 'lm'));

Сравнение между истиной и симулированными выходными параметрами (с предполагаемой моделью) показывают в окне графика.

compare(z3, nlgr3);

Рисунок 13: Сравнение между измеренными выходными параметрами и симулированными выходными параметрами первой предполагаемой модели транспортного средства Volvo V70.

Предполагаемые параметры жесткости итоговой модели Volvo V70 разумны, все же это здесь неизвестно, каковы их действительные значения.

disp('                        Estimated');
fprintf('Longitudinal stiffness: %6.0f\n', nlgr3.Parameters(4).Value);
fprintf('Lateral stiffness     : %6.0f\n', nlgr3.Parameters(5).Value);
                        Estimated
Longitudinal stiffness: 108873
Lateral stiffness     :  29964

Дополнительная информация о предполагаемой модели транспортного средства Volvo V70 получена через PRESENT. Здесь интересно отметить, что неопределенность, связанная с предполагаемой боковой жесткостью шины, довольно высока (и значительно выше, чем для продольной жесткости шины). Эта неопределенность происходит частично, от которого поперечное ускорение варьируется так мало во время тест-драйва.

present(nlgr3);
                                                                                                 
nlgr3 =                                                                                          
Continuous-time nonlinear grey-box model defined by 'vehicle_c' (MEX-file):                      
                                                                                                 
   dx/dt = F(t, u(t), x(t), p1, ..., p6)                                                         
    y(t) = H(t, u(t), x(t), p1, ..., p6) + e(t)                                                  
                                                                                                 
 with 5 input(s), 3 state(s), 3 output(s), and 2 free parameter(s) (out of 6).                   
                                                                                                 
 Inputs:                                                                                         
    u(1)  Slip on front left tire(t) [ratio]                                                     
    u(2)  Slip on front right tire(t) [ratio]                                                    
    u(3)  Slip on rear left tire(t) [ratio]                                                      
    u(4)  Slip on rear right tire(t) [ratio]                                                     
    u(5)  Steering angle(t) [rad]                                                                
 States:                                         Initial value                                   
    x(1)  Longitudinal vehicle velocity(t) [m/s]   xinit@exp1   17.6049   (estimated) in ]0, Inf]
    x(2)  Lateral vehicle velocity(t) [m/s]        xinit@exp1         0   (fixed) in [-Inf, Inf] 
    x(3)  Yaw rate(t) [rad/s]                      xinit@exp1         0   (fixed) in [-Inf, Inf] 
 Outputs:                                                                                        
    y(1)  Long. velocity(t) [m/s]                                                                
    y(2)  Lat. accel.(t) [m/s^2]                                                                 
    y(3)  Yaw rate(t) [rad/s]                                                                    
 Parameters:                                    ValueStandard Deviation                          
    p1   Vehicle mass [kg]                          1700            0   (fixed) in ]0, Inf]      
    p2   Distance from front axle to COG [m]         1.5            0   (fixed) in ]0, Inf]      
    p3   Distance from rear axle to COG [m]          1.5            0   (fixed) in ]0, Inf]      
    p4   Longitudinal tire stiffness [N]          108873      26.8501   (estimated) in ]0, Inf]  
    p5   Lateral tire stiffness [N/rad]          29963.5      217.877   (estimated) in ]0, Inf]  
    p6   Air resistance coefficient [1/m]            0.7            0   (fixed) in ]0, Inf]      
                                                                                                 
Name: Volvo V70 vehicle model                                                                    
                                                                                                 
Status:                                                                                          
Termination condition: Maximum number of iterations reached..                                    
Number of iterations: 20, Number of function evaluations: 41                                     
                                                                                                 
Estimated using Solver: ode45; Search: lm on time domain data "Volvo V70 data".                  
Fit to estimation data: [-374.2;29.74;34.46]%                                                    
FPE: 2.362e-07, MSE: 0.3106                                                                      
More information in model's "Report" property.                                                   

Заключительные замечания

Оценка параметров жесткости шины является на практике довольно сложной проблемой. Во-первых, приближения, введенные в структуре модели выше, допустимы для довольно узкой области операции и данных во время высоких ускорений, торможение, и т.д., не может использоваться. Жесткость также меняется в зависимости от состояния окружающей среды, e.g., окружающая температура, температура в шинах и условиях дорожного покрытия, которые не составляются в используемой структуре модели. Во-вторых, оценка параметров жесткости полагается в большой степени на ведущий стиль. При больше всего движении прямо вперед как в третьем идентификационном эксперименте, становится трудно оценить параметры жесткости (особенно боковой) или поместить иначе, неопределенность параметра становится довольно высокой.