建立车辆动力学系统的模型

This example shows nonlinear grey-box modeling of vehicle dynamics. Many new vehicle features (like Electronic Stability Programs (ESP), indirect Tire Pressure Monitoring Systems (TPMS), road-tire friction monitoring systems, and so forth) rely on models of the underlying vehicle dynamics. The so-called bicycle vehicle model is a rather simple model structure that is frequently being used in the vehicle dynamics literature. In this example we will start-off with this model structure and try to estimate the longitudinal and the lateral stiffness of a tire. The actual modeling work was originally carried out by Erik Narby in his MSc work at NIRA Dynamics AB, Sweden.

Vehicle Dynamics Modeling

The following figure illustrates the vehicle modeling situation to be considered.

Figure 1: Schematic view of a vehicle dynamics system.
By the use of Newton’s law of motion and some basic geometric relationships, the longitudinal velocity v_x(t), the lateral velocity v_y(t) and the yaw rate r(t) measured around the Center Of Gravity (COG) of the vehicle can be described by the following three differential equations:

  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)))

where subscript x is used to denote that a force F acts in the longitudinal direction and y that it acts in the lateral direction. The abbreviations FL, FR, RL and RR label the tires: Front Left, Front Right, Rear Left and Rear Right, respectively. The first equation describing the longitudinal acceleration also contains an air resistance term that is assumed to be a quadratic function of the longitudinal vehicle velocity v_x(t). In addition, delta(t) (an input) is the steering angle, J a moment of inertia, and a and b the distances from the center of gravity to the front and rear axles, respectively.
Let us assume that the tire forces can be modeled through the following linear approximations:

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

where C_x and C_y are the longitudinal and lateral tire stiffness, respectively. Here we have assumed that these stiffness parameters are the same for all 4 tires. s_i(t) is the so-called (longitudinal) slip of tire i and alpha_i(t) a tire slip angle. For a front-wheel driven vehicle (as considered here), the slips s_FL(t) and s_FR(t) are derived from the individual wheel speeds (measured) by assuming that the rear wheels do not show any slip (i.e., s_RL(t) = s_RR(t) = 0). Hence the slips are inputs to our model structure. For the front wheels, the tire slip angels alpha_Fj(t) can be approximated by (when 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}
For the rear wheels, the tire slip angels alpha_Rj(t) are similarly derived and computed as
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}
With J = 1/((0.5*(a+b))^2*m) we can next set up a state-space structure describing the vehicle dynamics. Introduce the states:
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].the five measured or derived input signals
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].and the model parameters:
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].The outputs of the system are the longitudinal vehicle velocity y1(t) = x1(t), the lateral vehicle acceleration (measured by an accelerometer):
y2(t) = a_y(t) = 1/m*(  (F_x,FL(t)+F_x,FR(t))*cos(delta(t))
                    - (F_y,FL(t)+F_y,FR(t))*sin(delta(t))
                    + F_x,RL(t)+F_x,RR(t))
and the yaw rate y3(t) = r(t) (measured by a gyro).
Put together, we arrive at the following state-space model structure:
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 Vehicle Model

As a basis for our vehicle identification experiments we first need to create an IDNLGREY model file describing these vehicle equations. Here we rely on C-MEX modeling and create a vehicle_c.c model file, in which NY is set to 3. The state and output update functions of vehicle_c.c, compute_dx and compute_y, are somewhat involved and includes several standard C-defined mathematical functions, like cos(.) and sin(.) as well as pow(.) for computing the power of its argument.
The state update function compute_dx returns dx (argument 1) and uses 3 input arguments: the state vector x, the input vector u, and the six scalar parameters encoded in p (t and auxvar of the template C-MEX model file have been removed here):

  /* 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]);
  }

The output update function compute_y returns y (argument 1) and uses 3 input arguments: the state vector x, the input vector u, and five of the six parameters (the air resistance CA is not needed) encoded in 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];
}

Having a proper model structure file, the next step is to create an IDNLGREY object reflecting the modeling situation. For ease of bookkeeping, we also specify the names and units of the inputs and outputs.

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 initial states.
Ts            = 0;                                  % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ...
                'Name', 'Bicycle vehicle model', 'TimeUnit', 's', ...
                'Display', 'on');
set(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).
          'InputUnit',  {
   'ratio'; 'ratio'; 'ratio'; 'ratio'; 'rad'});
set(nlgr, 'OutputName', {
   'Long. velocity';  ...   % y(1); Longitudinal vehicle velocity
                         'Lat. accel.';   ...     % y(2); Lateral vehicle accelerat