最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

車両運動システムのモデル化

この例は、車両運動の非線形グレーボックス モデリングを示しています。車両の機能 (横滑り防止装置 (ESP)、間接的タイヤ プレッシャー モニタリング システム (TPMS)、ロード タイヤ摩擦モニタリング システムなど) は、基礎として車両運動のモデルを利用しています。いわゆる自転車車両モデルはシンプルなモデル構造で、車両運動の資料ではよく使用されています。この例では、このモデル構造から始めて、タイヤの前後方向および横方向の剛性を推定してみます。実際のモデル化作業は、Erik Narby が NIRA Dynamics AB (スウェーデン) で修士課程の研究として行ったものです。

車両運動のモデル化

以下の図は、検討する車両モデル化状況を示しています。

図 1: 車両運動システムの概略図

ニュートンの運動の法則と一部の基本的な幾何学的関係を使用して、車両の重心 (COG) で測定した前後方向の速度 v_x(t)、横方向の速度 v_y(t)、ヨーレート r(t) を、次の 3 つの微分方程式で記述できます。

  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 はそれぞれ次のタイヤを表します。前左 (Front Left)、前右 (Front Right)、後左 (Rear Left)、後右 (Rear Right)。前後方向の加速度を表現する最初の方程式にも空気抵抗の項があり、これは前後方向の車両速度 v_x(t) の 2 次関数と仮定されます。また、delta(t) (入力) はステアリング角度、J は慣性モーメント、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].

5 つの測定または導出された入力信号

  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 モデル ファイルを作成します。ここで、NY は 3 に設定されています。vehicle_c.c の状態および出力更新関数の compute_dx と compute_y は多少複雑で、いくつかの標準 C 定義数学関数 cos(.) および sin(.) や、引数のべき乗を計算する pow(.) が含まれています。

状態更新関数 compute_dx は dx (引数 1) を返し、次の 3 つの入力引数を使用します。状態ベクトル x、入力ベクトル u、および p でエンコードされた 6 つのスカラー パラメーター (テンプレート C MEX モデル ファイルの t と auxvar はここでは削除されています)。

  /* 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、p でエンコードされた 6 つのパラメーターのうちの 5 つ (空気抵抗 CA は不要)。

  /* 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!

このモデル構造の 6 つのパラメーターのうちの 4 つは、該当する車両のデータ シートからすぐに取得できます。

  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.                                            

入出力データ

この時点で、使用できる入出力データを読み込みます。このファイルには、次の 3 つの異なる実験からのデータが含まれています。

  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. シミュレートした高剛性タイヤ データを使用したシステム同定

最初の車両同定実験では、高剛性タイヤをシミュレーションしたデータを検討します。モデル構造 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';

入出力は、2 つのプロット ウィンドウで表示されます。

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: 高剛性タイヤの初期車両モデルの実際の出力とシミュレーションの出力の比較

モデルの一致を高めるため、2 つのタイヤ剛性パラメーター Cx と Cy を次に推定し、推定されたモデルの新しいシミュレーションを実行します。

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. シミュレートした低剛性タイヤ データを使用したシステム同定

2 番目の実験で、最初の実験のモデル化を繰り返しますが、今度はシミュレートした低剛性タイヤ データを使用します。

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';

入出力は、2 つのプロット ウィンドウで表示されます。

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: 低剛性タイヤの初期車両モデルの実際の出力とシミュレーションの出力の比較

次の 2 つの剛性パラメーターを推定します。

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. 測定した Volvo V70 データを使用したシステム同定

最後の実験では、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';

入出力は、2 つのプロット ウィンドウで表示されます。見てわかるように、測定されたデータにはノイズが含まれています。

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 車両モデルの測定された出力とシミュレーションの出力の比較

次に、タイヤ剛性パラメーター Cx と Cy を推定します。ここでは、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

推定された Volov 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.                                                   

結果のまとめ

タイヤ剛性パラメーターの推定は実際には、複雑な問題です。最初に、上記で説明したモデル構造の近似が有効な運転範囲はかなり狭く、高加速度、ブレーキなどの間のデータは使用できません。剛性は環境状況によっても異なります。周辺温度、タイヤの温度、道路表面の条件などで、これらは使用したモデル構造には反映されていません。2 番目に、剛性パラメーターの推定は、運転方法に大きく依存します。3 番目の同定実験のようにほとんど直線に進む場合、剛性パラメーター (特に横方向) の推定は困難になります。つまり、パラメーターの不確かさが高まります。