exponenta event banner

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

Этот пример показывает нелинейное моделирование серого ящика динамики аппарата. Много новых функций транспортного средства (как Электронные программы устойчивости (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) выведены из отдельных скоростей колеса (измеренных) путем предположения, что задние колеса не показывают промаха (т.е. 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.                                                   

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

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