MPC and PSO Based Control Methodology for Path Tracking of 4 WS 4 WD Vehicles

Qifan Tan 1,* ID , Penglei Dai 2 , Zhihao Zhang 3 and Jay Katupitiya 3 1 School of Mechanical, Electronic and Control Engineering, Beijing Jiaotong University, Beijing 100044, China 2 School of Mechanical and Mechatronic Engineering, University of Technology Sydney, Sydney, NSW 2007, Australia; penglei.dai@uts.edu.au 3 School of Mechanical and Manufacturing Engineering, The University of New South Wales, Sydney, NSW 2052, Australia; zhihao.zhang1@unsw.edu.au (Z.Z.); j.katupitiya@unsw.edu.au (J.K.) * Correspondence: 14116355@bjtu.edu.cn; Tel.: +86-186-0124-6239


Introduction
With the development of Autonomous Ground Vehicles (AGVs) in the last few decades, the demand for accuracy, maneuverability and controllability in vehicle's navigation is ever increasing.For example, an AGV may be required to follow a path accurately under unstructured and uneven terrain conditions, where a significant amount of wheel slip and unpredictable disturbance forces occur at the vehicle's wheels.The 4WS4WD vehicle, with four wheels that can be steered and driven independently, is a revolutionary platform that has great potential to perform high maneuverability and flexibility in harsh environments.
The main challenge in the control of 4WS4WD control is the number of control inputs (four steering angles and four drive torques), which results in an over-actuated system, where only three outputs including its degree of freedom (DOF) in the longitudinal, lateral and angular directions of the vehicle are concerned.How to allocate all eight control inputs to achieve high path following performance has not yet been effectively solved.The control allocation is proposed to handle the control problems of over-actuated systems [1].Generally, the control allocation can be treated as an optimization problem.
For controller designing, a model of the vehicle under control is generally required to facilitate the selection of future control inputs.A dynamic model describes the states of the vehicle based on the forces applied.However, the development of a detailed vehicle dynamic model is always a challenging task due to the uncertainty of parameters and the complex disturbances from the external forces.The majority of existing control methodologies for 4WS4WD vehicles are proposed based on the linearized dynamic model, which lead to the loss of input degree of freedom [2,3].Meanwhile, some other methods use nonlinear dynamic models in the control of 4WS4WD vehicles [4,5], where control inputs are subjected to some relationship constraints to simplify the controller design.However, in practice, the four independent wheels may interact with different terrain conditions, where different slip, wheel forces and terrain disturbances are generated on the corresponding contact patches.Hence, it is desirable to make four wheels individually controlled, thereby limiting and/or overcoming different slip and disturbance on different contact patch.
There is an abundance of literature that presents kinematic modelling of ground vehicles [6][7][8], in which the vehicles are assumed to operate at low speeds to reduce the dynamic effects.Most vehicle kinematic models are developed based on the non-integrable kinematic constraints, known as non-holonomic constraints.As a result, the wheel slip has to be ignored with the assumption of zero relative velocity between wheels and terrain [9,10].To utilize the kinematic model's advantage of keeping the steering control relatively independent of velocity control [11], it is desirable to incorporate wheel slip in the vehicle kinematic model so as to facilitate the accurate vehicle control in complex terrain conditions where the no-slip assumption is not applicable.
To realize force control, a novel type of 4WS4WD vehicle with force sensors at each wheel has been designed as shown in Figure 1.This vehicle possesses the characteristics of independent steering and drive control at each wheel and force measurements.In this work, the dynamic model is partitioned into a hierarchy of three levels to facilitate the incorporation of force sensors and overcome uncertainties in the model.The force sensors at the drive unit allow the measurement of actual force data, while models of the drive unit and tire are used for simulation.Model predictive control (MPC) is selected for its ability in handling linear constraints and time-varying systems as well as its good performance in tracking problems.Particle swarm optimization (PSO) is also selected for its fast searching speed in global optimization.MPC and PSO have been successfully applied in the controller design of real-time control systems [12][13][14][15][16][17].In this work, MPC and PSO based control methods are proposed to realize the controlling aim of achieving good path following performance as well as high motion quality via vehicle steering control and independent force control at four wheels.

Force Sensor Driving Unit
As novel contributions, the MPC methodology is applied to achieve precise path tracking of 4WS4WD vehicle.Based on the MPC theory, an offline control law is proposed to guarantee the stability of the upper controller.An sequential quadratic programming (SQP) based control allocation is developed to control the 4WS4WD vehicle in the lower controller.The inclusion of full independent force control and steering control on all four wheels enable the maximization of performance.In this work, comparison of MPC and PSO on the same vehicle model is provided, in which the proposed PSO control methodology is a further refinement of the PSO methodology presented previously in [18].The PSO control methodology in this paper simplifies the derivation and gives an algorithm in a more general form, which facilitates the comparison with other control methods.In addition, both methodologies are compared with the kinematic model based method proposed in [19].
The paper is organized as follows: Section 2 describes the 4WS4WD vehicle modelling.The MPC and PSO control methodologies are presented in Section 3. In Section 4, simulation setup and the reference path are presented.The result of the two controllers are compared.In Section 5, the discussion about two methodologies are provided based on their theoretical analysis and simulation results.Finally, conclusions are provided in Section 6.

Vehicle System Modelling
In order to develop suitable control methodology, the 4WS4WD vehicle must be modelled for controller design and simulation.The model used for this work is a fully dynamic model that consists of three components, vehicle body dynamics to estimate vehicle body movement, drive unit dynamics to model the tire force acting on the drive unit in vehicle coordinate frame and a tire dynamic model to model the force generated by the tire in a wheel coordinate frame.The separation of dynamic model allows force sensors to be incorporated easily and reduces uncertainties by obtaining actual force measurements.In the end, an offset error model is proposed for evaluating the tracking performance and is used in the controller design.

Vehicle Body Dynamic Model
As shown in Figure 2, the vehicle body dynamic model describes the motion of the vehicle by representing each drive unit as a pair of forces, with F Sl in the longitudinal and F SL in the lateral direction of the corresponding wheel.The reason for separating the dynamics of the vehicle at the drive unit boundary is to allow force sensors to be incorporated to measure the forces acting on each wheel, instead of relying on the tire model to estimate them.To facilitate the notations, the four wheels are numbered by 1, 2, 3 and 4 in circles.In particular, the force acting at each force sensor is decomposed into two components (i.e., F Sli and F SLi , i = 1, . . ., 4) normal to each other.
The equation of the vehicle body dynamics can be expressed as where the acceleration vector denoted by acc consists of longitudinal acceleration a l , radial acceleration a r and angular acceleration γ.
The mass matrix is represented by M: where M is the mass of the vehicle, J z is the vehicle body inertia, and m d is the mass of the drive unit.
Vehicle body dynamic model.
The matrix C is decided by vehicle dimension, written as The force vectors F ASi are defined as The transfer matrices R i from wheel frame to vehicle local coordinate frame The lateral components F SLi are determined by the lateral forces acting on tires, while the longitudinal components F Sli are viewed as intermediate control inputs of vehicle system, which can be achieved by controlling the torques (i.e., T i , i = 1, . . ., 4) applied on wheels.To simplify the coupled issue of left and right wheel steerings, the steering angles are constrained by Therefore, the vehicle system in this work considers six control inputs (i.e., δ f , δ r and T i , i = 1, . . ., 4) in total.

Driving Unit Dynamic Model
While the force sensors provide measurement of the current forces acting on each wheel, controller design requires prediction of future values based on control inputs.The drive unit dynamic model as shown in Figure 3 is described in this section.In each driving unit, the force sensor mounted on the wheel hub rotate with the wheel.Force measurement is in the direction of the steering angle.According to Figure 3, the dynamics of driving unit can be expressed by where a wli and a wLi denote the longitudinal and lateral accelerations of each driving unit i in the wheel coordinate system (X Wi − O Wi − Y Wi ), respectively.Using the transfer matrices R i and the accelerations a l , a r and γ given in vehicle body dynamic model, a wli and a wLi can be deduced as In Equation (7), F wi represents the force vector acting on each wheel, which is written as

Tire Model
The actual generation of wheel forces are from the contact between tires and the ground.The wheel forces depend on the surface friction, load, and the intrinsic properties of the tire.To analyze the wheel forces F wi in Equation ( 9), the tire model is built as shown in Figure 4.
Considering the generation of wheel forces and external disturbances, F wi can be expressed as where F li and F Li (i = 1, . . ., 4) are longitudinal and lateral forces caused by wheel slip, respectively.F gli and F gLi denote the terrain disturbances corresponding to F li and F Li .
According to [20], the longitudinal force can be considered proportional to the slip ratio in small range, which is expressed by where F Ni is the weight force on wheel i.In the simulation, the load transfer due to the longitudinal and lateral accelerations and roll angle are taken into consideration.F Ni can be calculated by considering longitudinal and lateral load transfer according to [21].The longitudinal slip stiffness k li is determined by the tire type and terrain condition.ς i is the longitudinal slip ratio of wheel i presented in [22]: where R w and ω i denote the tire radius and angular velocity of wheel i as shown in Figure 4b.v wi represents the actual velocity of the wheel i, which can be calculated by the vehicle geometry as well as velocities V l , V r and Ω shown in Figure 2. Note that the model is proposed only considering the vehicle is moving forward.Thus, v wi and ω i are set to be nonnegative.As is discussed in Equation ( 12), the slip ratio is zero when the vehicle is resting.Finally, as per Figure 4b, the wheel dynamic equation indicating the relationship between T i and F li can be written as where J w represents the wheel inertia, and F ri is the rolling resistance calculated by where K r0 = 0.015 and K r1 = 7 × 10 −6 s 2 /m 2 are the parameters for common car tires [23].The steering motion of wheels result in the lateral slip velocities v si , which causes the actual direction of wheel velocity v wi to differ from the wheel center plane by the slip angle α i .According to the tire lateral characteristic curve [24], a linear model for F Li is given in where k Li is the lateral slip stiffness of the wheel i.Note that, due to the sign conventions in a negative slip angle causes a positive lateral force and vice versa.
As shown in Figure 4a, the slip angle α i can be calculated by steering angle δ i and side slip angle β i , where the expression of β i is available in our previous work presented in [25].

Offset Model
An offset error model is proposed for the evaluation of the performance of the controllers.The reference position (RP) of the vehicle on the reference path is defined as the normal projection of the heading of the vehicle at the centre of gravity on to the reference path, as shown in Figure 5.The tracking error in vehicle coordinate frame consists of lateral offset error l os and heading error θ os , expressed as The position errors are always in the lateral direction of the vehicle.The longitudinal position error is set to 0. The heading error can be found by where θ re f is the reference heading direction of the RP, tangential to the reference path.
The velocity error states vel is the difference between vehicle velocity vel and reference velocity states vel re f in the vehicle coordinate frame.vel can be expressed as where V lerr , V rerr and Ω err are the errors in longitudinal, lateral, and angular velocity, respectively.Reference velocity states are defined as where reference angular velocity Ω re f is the rate of change of the heading of RP on the reference path.
The reference longitudinal velocity V lre f and lateral velocity V rre f in vehicle coordinate frames can be derived from the reference path coordinate frame by the equations: where V l p is the reference longitudinal velocity, taken as tangential to the reference path.V rp is the reference radial velocity normal to the reference path at RP. Finally, the acceleration error states acc are written as: where acc is the difference between vehicle acceleration acc and reference acceleration acc re f .The reference acceleration vector acc re f is defined as where a lre f and a rre f are, respectively, the reference longitudinal and lateral acceleration at point RP, measured in the direction and normal to the vehicle heading, and γ re f is the reference angular acceleration of the point RP.
The reference acceleration in vehicle coordinate frame can be derived from the path coordinate frame by a lre f = a l p cos(θ os ) − a rp sin(θ os ), where a l p and a rp are the reference acceleration tangential and normal to the reference path.

Control Methodology
Two methods of finding the optimal control inputs for steering angles (i.e., δ f and δ r ) and drive forces F Sli , (i = 1, ..., 4) are developed.The first method used model predictive control (MPC) with a sequential quadratic programming (SQP) solver.The second method is particle swarm optimization (PSO).The objective of the controller is to accurately follow a predetermined path through multiple terrain types.

Offset Model Linearization
In order to be applied in the MPC controller design, the offset model is linearized by a feedback linearization method.By defining the state vector x = [x 1 , x 2 , x 3 , x 4 , x 5 ], as x 1 = V lerr , x 2 = l os , x 3 = los , x 4 = θ os and x 5 = θos , the offset model can be written in the following form: In this model, a r , a l and γ are considered as the inputs while other constants including reference accelerations and velocities can be obtained from the vehicle states.The nonlinearity of the model comes from trigonometric terms of x 4 .For accurate path tracking problems, the yaw error is assumed to vary smoothly within [−10°10°].Then, the model can be linearized by feedback linearization.Define the input vector u = [a r a l γ] T ; then, the new input vector at time t k is obtained as and the offset model is expressed as where It can be seen that the model is time-varying but can be treated as a linear model at each sampling step.The output vector is denoted by y c and the equation is written as where

Model Predictive Control
To develop the MPC controller, the offset model as given in Equations ( 27) and ( 28) needs to be discretized and its discrete-time state vector form can be written as are defined, thereby an augmented state-space model can be expressed as where O is a zero matrix, and I represents the identity matrix.Defining the new state vector T , the augmented model can be written in the following matrix form: where Theorem 1.Given a discrete time system following the form of Equation (31), the asymptotic stabilization of the closed-loop system can be realized by substituting the first item of ∆U * as the control input ∆u(k) , when ∆U * is the optimal solution of the following optimization problem: where Y denotes the predicted output sequence, ∆U denotes the future input sequence, and N p is the prediction horizon.R s is the sequence of the control target vector.y N p (∆U) represents the error between the final predicted output and target.
Proof of Theorem 1.To prove the stability, the Lyapunov function V(x k ) is defined equal to the value of the objective function J k subjected to its optimal solution, which can be expressed as where, ∆u k , ..., ∆u k+N p −1 are obtained by the optimal solution of future inputs, and y k , ..., y k+N p represent the corresponding error sequence between future outputs and target.r w is the nonnegative gain matrix.
According to the definition in Equation (33), the Lyapunov function at the next sample k + 1 is written as To facilitate the comparison between two neighboring Lyapunov function values, a intermediate function V is defined, which is formed by evaluating V(x k+1 ) with a defined inputs sequence, which is obtained by shifting the optimal inputs sequence of V(x ki ) one step forward, and setting its last input ∆u k+N p as zero.It is obvious that the objective function value of non-optimal inputs sequence has to be no less than V(x ki+1 ), which can be expressed as Since V shares the same future inputs sequence and the predictive outputs sequence with V(x k ) for the sample time k + 1, ..., k + N p − 1, it can be easily derived that the difference between these two functions is As is given in Equation (32), the optimization problem is subjected to the constraint y k+N p = 0.Then, it can be obtained that Then, the monotonicity of the Lyapunov function can be obtained by which can prove the asymptotic stability of the system.
The model predictive control algorithm is realized by receding optimization.In order to apply the MPC efficiently, we assume that the predicted outputs sequence is in a finite prediction horizon N p and the inputs sequence is in a control horizon N c , which is less than N p .The sequences mentioned above can be expressed in the matrix form: where y(k + n|k) denotes the predicted outputs at time k + n based on the states at time k.Based on Theorem 1 and the assumption above, the corollary about finite-time unconstrained MPC can be obtained as follows.
Corollary 1.Given the system without input and output constraints, the prediction and control horizon are N p and N c , respectively.Then, the following feedback control law ∆u(k) = Kx(k) can asymptotically stabilize the closed-loop system, where and, Proof of Corollary 1.According to Equations ( 31) and (40), the predicted output sequence Y can be expressed by Meanwhile, the final item of Y can be expressed as By substituting Equations ( 43) and (44) into Theorem 1, the optimization problem turns into the following form: arg min where R is a weighting matrix.Note that in the application of offset model based path tracking, the target vector R s should be zero all the time.It can be seen that J meets an equality constrained quadratic programming.Then, the objective function is expanded by Lagrange expression and simplified by omitting the constant term, where ξ is the Lagrange multiplier.
According to the Lagrange multiplier method, the optimal control input vector ∆U * can be found by solving the following equation system.The solution is obtained by taking the first partial derivatives of J with respect to the vectors ∆U and λ, and then equating these derivatives to zero: Then, its optimal solution can be obtained: where According to the receding horizon control principle, the first increment of ∆U * is applied as the control inputs.Then, the control law in Equation (41) can be obtained.
It can be seen that Corollary 1 gives an offline solution of the MPC algorithm, which significantly improves its computing efficiency.Based on Corollary 1, the desired accelerations A = [ a r , a l , γ] T can be obtained by integrating the optimal solution ∆u(k) * .

Sequential Quadratic Programming Based Control Allocation
Considering the vehicle body with the mass M and inertia J in this work, the command force vector τ is defined as where τ l , τ r and τ γ are the longitudinal and lateral forces and the moment about a vertical axis of the vehicle, respectively.According to the vehicle body dynamic model given in Equation ( 1), the actual actuating forces come from the longitudinal forces F Sl and lateral forces F SL on each wheel.Let the command forces τ produced jointly by the wheel forces and steering angles be expressed as where the i-th column of B u (δ) and B w (δ) can be written as (L xi , L yi ) represents the location of each wheel in a coordinate system with its origin at the centre of gravity and positive x-axis forward.
On a 4WS4WD vehicle, the drive forces F Sl and the steering angles δ are the direct inputs while the lateral forces F SL obtained from sensors are considered as a measured disturbance.In this work, the steering angles δ are composed of δ f and δ r , which represent the front and rear steering inputs, respectively.Thus, the control problem is reduced to obtaining the feasible solution of Equation (50).In order to facilitate the computation, a slack vector s is defined by which denotes the error between the commanded and actual generalized forces.The slack variable s guarantees that there always exists a feasible solution in the following optimization [26].
In order to solve this control allocation problem, the objective function is defined with respect to F Sl , δ and s, where B s , B F Sl and B δ are the search bounds of each variable.
In this function, the first term minimizes the magnitudes of the drive forces; the second term is used to ensure the steering angle to search around its previous value.By penalizing the slack variable s in the third term, the actual generalized force vector coincides as much as possible with the commanded forces τ.The matrix Q f ∈ I 4×4 , Q δ ∈ I 2×2 and Q s ∈ I 3×3 are used to tune the objective.The search bounds of all variables (i.e.F Sl , δ and s) are specified by the constraints.
Based on the objective function presented above, the control allocation is converted to a nonlinear constrained optimization problem.Using the sequential quadratic programming (SQP), the optimal solution can be computed efficiently and reliably by standard numerical software.

PSO-Based Method
For the PSO-based method, an objective function including vehicle error states is proposed using the sliding surfaces.In Sliding Mode Control (SMC), the time-varying sliding surface is normally defined by the scalar equation s(x; t) = 0, in which s(x; t) is expressed by [27], where λ is a positive constant and x is the error state vector.
According to the idea of SMC, the problem of maintaining x = 0 is transformed into keeping s = 0.In this work, the scalar quantities composed of vehicle error states are introduced into the definition of objective function.Instead of designing the switching control law in SMC, the optimization is used to maintain the scalar quantities on the sliding surface s = 0.The vehicle error state vectors including pos, vel and acc are following the definitions in the offset model.Therefore, to track the trajectories of the RP, the vector s = [s l , s r , s a ] needs to be defined correspondingly.Note that vector s follows a different definition than that in SQP.It is defined to facilitate the notations in the local boundary part.
For pos presented in Equation ( 17), its first component, which represents the position error in the longitudinal direction, is always zero.Therefore, according to Equation (53), the longitudinal scalar quantity s l can be obtained as where a lerr and V lerr are the longitudinal components of acc and vel, which are given in Equations ( 22) and ( 19), respectively.
To maintain the vehicle on the RP geometrically, another first order scalar quantity s r is designed as which aims to minimize the offset errors l os and V rerr .The angular acceleration error γ is included to consider the effects of forces and yaw movements of vehicle, and thus a second order scalar quantity is chosen as where θ os , Ω err and γ err are given in Equations ( 17), ( 19) and ( 22), respectively.
Then, the problem of following the trajectories of RP is transformed into maintaining s at 0. Using the linear scalarization, an objective function is defined as where C l , C r and C a are the weighting coefficients strictly positive and constrained by The variables of objective function in Equation (57) consist of the forces F Sli , steering angles δ f and δ r , which can be written as a variable vector, First invented by Kennedy and Eberhart (1995), PSO has been successfully applied to solve problems featuring nonlinearity, non-differentiable, and multiple optima.PSO is found to be capable of generating high quality solutions with more stable and faster convergence characteristics, and shorter calculation times than other stochastic methods [17].For standard PSO at time t, the updating velocity v i (t) and position x i (t) of the i-th particle are presented in the following equations: where v i (t) and x i (t) are vectors in multi-dimensional space.p b and g b denote the local optimal position and the global optimal position, respectively.The particle inertia weight is represented by ζ.
The particle cognitive acceleration and social acceleration are denoted by φ 1 and φ 2 , which are defined as positive constants.η 1 and η 2 are two stochastic parameters within [0 1].The search space of PSO in this work is defined as a six-dimensional space corresponding to the dimension of v obj .Therefore, the particle position vector x i (i) in Equation (60) represents a possible solution of the objective function in Equation (57).

Boundary Definition
Given that both methods are reduced to the optimization problems, the definition of the search space determines the quality of the solution.

Global Boundaries
In this work, the global boundaries can be assigned based on the properties of each actuator, which is defined as where δ min and δ max are the minimum and maximum steering angles of each wheel.F d min and F d max are the minimum and maximum drive forces provided by the driving unit, which can be calculated by where T max is the maximum torque that can be achieved by the driving unit of the vehicle.

Local Boundaries
To realize real-time optimization, the computing time obtaining optimal values for variables [F Sl , δ] needs to be constrained within the sample time of controller T s .According to the properties of each actuator, the maximum variations of [F Sl , δ] within T s can be obtained.In this work, to improve the computing efficiency, the local boundaries are also determined by the states of s.
According to a r in Equation (1), when s l < 0, the forces F Sli need to be increased, thereby dragging s l in Equation (54) towards the surface s l = 0. Similarly, when s l > 0, the forces F Sli need to be decreased.Therefore, the local boundaries of F Sli can be written as where ∆F d max is the maximum change of F Sli that can be achieved within T s .At the time t + T s , F Sli (t + T s ) are the possible optimal solutions.For the steering angles δ f and δ r in v obj , they affect the vehicle lateral and angular motions in a coupled way.To solve this problem, an allocating method is specified that the vehicle lateral error determines searching direction of δ f , while δ r is related to vehicle angular error.When s r < 0, δ f needs to be increased to drive s r back to the surface s r = 0, and vice versa.Thus, the local boundary of δ f is summarized as, where ∆δ max is the maximum change of steering angle that can be achieved within T s .At the time t + T s , the possible optimal solution of front steering angle is δ f (t + T s ).
Similarly, based on s r , the local boundary of δ r is defined as For the PSO in particular, its particle velocity boundaries define the range of speed that particles can achieve to search for the optimal solution.To improve search performance in PSO, the absolute maximum particle velocity is normally set as a certain percentage of particle position range [28].According to Equations ( 63)-(65), the particle velocity boundaries can be obtained as where σ F , σ δ f and σ δ r are the particle velocity coefficients for F Sli , δ f and δ r , respectively.

Simulation Setup
The working process of the vehicle control system is illustrated in Figure 6.The dashed box at the top shows the operation of the vehicle dynamic model.Based on the actual inputs (i.e., T i and δ i ) as well as the vehicle velocity states vel, the tire model can generate the wheel forces, which is then used in driving unit model.The driving unit model considers the wheel inertia and gives the actual forces acting on the vehicle body.The dynamic states of the vehicle such as positions and velocities are simulated by running the proposed vehicle body dynamic model in Matlab (R2017a, MathWorks, Natick, MA, United States) using a Runge-Kutta method based solver with a time step of 8.33 µs.In the simulation, both external disturbances and state measurement noises are involved in validating the robustness of the proposed control methodology.The measured states are compared with the reference profiles to generate the error states.Using the method proposed in Section 3.3, the search boundaries of the drive forces and steering angles are obtained.Meanwhile, the error states are transfered to the offset model.Virtual inputs calculated by the MPC algorithm are delivered to the control allocation modual.Then, the actual control inputs including drive forces and steering angles are generated and substituted into the dynamic model for the next iteration.The dashed box at the bottom represents the main controller.PSO controller can be substituted in place of the MPC controller.
The PSO and MPC controllers are applied to drive the simulated vehicle to track the reference path shown in Figure 7a.In the simulation, the varying terrain conditions are considered to verify the validity and robustness of each controller.As presented in Figure 7b, the reference path is divided into ten sections, which have different tire stiffnesses k li and k Li acting on each wheel.To simulate the terrain disturbances, the disturbances F gi applied to all four wheels were always in the vehicle longitudinal direction.As shown in Figure 8, the disturbances are modeled as step signals and their magnitudes are generated randomly to show the uncertainties.
As a comparison, a kinematic model based control method is applied to drive the same vehicle.According to relative kinematic path tracking research [8,19], the side slip is the main disturbance that leads to the unpredicted tracking errors.Using the dynamic model based observer, the side slip of the vehicle can be predicted with an error of 10% to 30%.In this simulation, a random side slip velocity that is less than 10% of its longitudinal velocity is added in the kinematic model.

Parameters of Simulation
The first two sections of Table 1 list the parameters of vehicles used in simulation.The common parameters in the objective functions are listed in the third section of Table 1, in which the vehicle mass and inertia values used in the objective function are different from actual vehicle values to simulate parametric uncertainties.This helps to verify the robustness of the vehicle control method in dealing with variations in the system.In order to run in real time, both optimization methods are limited to 15 ms of search time per iteration, less than the system sample time T s = 20 ms.In optimization, the optimal solution is searched within the boundaries until time reaches TI or the error gradient of 0.001 is achieved by the objective function.The last two sections present the particular parameters for each method.

-
I i×i represents the i-th-order identity matrix.MPC: model predictive control; SQP: sequential quadratic programming; PSO: particle swarm optimization.

Simulation Results
Figure 9 shows path tracking results by different controllers.The MPC-SQP controller has superior performance with offset error less than 3.2 cm, and heading error less than 2 • .The PSO controller has maximum offset error of 6.1 cm and heading error of 3.2 • .Both dynamic model based control methods perform better than the kinematic model based one in constraining the offset errors.The drive torques and steering angles applied on four wheels are presented in Figure 10.At each cornering, the torques applied at the two outside wheels (i.e., T 1 and T 4 ) firstly increase and then decrease, which are contrary to that of the two inside wheels (i.e., T 2 and T 3 ).The difference in torques is used to compensate the insufficient angular accelerations of vehicle in corners.The steering angles of MPC-SQP and PSO controller follow the same varying trend in which steering curves of MPC-SQP show smoother variation.Thus, the MPC-SQP controller can provide more stable steerings compared with the PSO controller.In this work, another aim of control is to maintain high vehicle motion quality including its velocity and acceleration performances.Figure 11a presents the longitudinal velocity curves achieved by MPC-SQP and PSO controllers.Starting at 0.5 m/s, the two curves reach 3 m/s at 6.5 s and 4.8 s , respectively.Thus, the trajectory and error curves of the MPC controller responses are slower than the PSO controller in Figure 9.Both controllers are capable of maintaining the longitudinal velocity around 3 m/s in the rest process.In Figure 11b, the longitudinal acceleration curve of PSO increases faster than MPC-SQP and has a peak value of 1 m/s 2 .The curve of MPC-SQP controller reaches 0.6 m/s 2 at a maximum and has smoother variation during the whole process.
From the offset model, lateral velocity is undesirable as it causes tracking error.As shown in Figure 11c, both controllers maintain the lateral velocities fluctuating around zero.The MPC-SQP controller has a maximum lateral velocity of 0.08 m/s, whereas PSO has a maximum value of 0.17 m/s.The lateral velocities are well constrained, which provides higher accuracy of path tracking.The lateral acceleration curves follow the same varying trend in Figure 11d.It can be seen the PSO controller causes a lateral acceleration oscillation of 0.6 m/s 2 relative to the MPC-SQP.The smoother change of lateral acceleration achieved by MPC-SQP controller reduces the jerk effect, decreasing the deviation from the reference path.The angular velocities are demonstrated in Figure 11e, in which the MPC-SQP has higher accuracy and a smoother manner, compared with that of the PSO.In Figure 11f, the angular acceleration curves are presented.MPC-SQP and PSO curves vary with oscillations of 12 °/s 2 and 20 °/s 2 , respectively.The undesirable oscillations are eliminated in the process of the integration, which can be demonstrated in Figure 11e.From all the acceleration plots given in Figure 11, the MPC-SQP controller performs better in constraining oscillations of accelerations, which provides better stability in vehicle motion control.
Considering that the optimization based method may lead to an expensive computation, it is essential to analyze the computing efficiency for each controller.As shown in Figure 12, the computing efficiency of each controller is compared.Figure 12a shows the computing time used by the controllers in each sampling time T s .It can be seen that both controllers achieve finishing computing within T s , which validates the feasibility of proposed controllers.In the box plot in Figure 12b, both controllers have an average computing time of 8 ms.The PSO controller gives a stable variation range from 4 ms to 10 ms, as the particle velocity vector is decided by the limitation of the sampling time, which guarantees relatively high quality solutions within a short time.The computing time of MPC-SQP controller substantially increases and reaches 15 ms to 20 ms during each cornering because the control allocation may encounter a complex optimization problem when the lateral forces start to vary.It indicates that the PSO controller performs better in computing efficiency.

Discussion
Based on the theoretical analysis and simulation results provided above, the discussion between MPC based method and PSO method are given as follows:

D1
Comparing the simulations, it is obvious that dynamic-based methods proposed in this paper perform better than kinematic model based methods.In the kinematic model based methods, it is difficult to measure the side slip directly.The majority of research works are trying to design observers to predict its side slip.However, this kind of method can only get an approximate estimate.Thus, it is not feasible to use a kinematic model based controller to completely eliminate the error due to side slip.Both MPC-SQP and PSO methods are proposed based on dynamic models and controlled by drive forces.In the dynamic model, the lateral forces can be obtained easily using force sensors.This kind of method gives a practical way to avoid the side slip estimation and achieve precise tracking along curved paths.

D2
Both MPC-SQP and PSO methods are proposed partly or completely based on optimizations.
As is compared in the simulation, the PSO controller achieves obtaining the solution with a stable computing time.The algorithm can optimize the quality of solutions and the calculation times at the same time, which guarantees the feasibility and capability of the controller.The MPC-SQP controller may encounter the calculation timeout in the simulation results.This is because the SQP solver spends more time when calculating a complex Hessian matrix.Thus, it is necessary to set a sufficient sampling time when applying the MPC-SQP method.

D3
The MPC-SQP method gives a stability proof of the control system while the PSO method has difficulty with the mathematical proof due to the limitation of intelligent optimization.
The stability analysis of MPC-SQP method provides a possibility to analyze the stable range of its control parameters and margin of errors.

D4
Both control algorithms are reduced to constrained optimization problems with six input variables.PSO is a global algorithm that performs better with searching for a global optimal solution while SQP is a reliable solver but may be held in local optimal solutions.Thus, for each method, it is important to define a proper search range.In this paper, the boundaries are calculated based on the vector s, which improves the efficiency and stability of the optimization process.

Conclusions
This paper presented two control methodologies applicable to four wheel steering and four wheel drive vehicle systems to track paths accurately.The system model is a nonlinear coupled dynamic model that is over-actuated.The first methodology determines the drive force inputs and steer inputs using an MPC based method coupled with a control allocation based on SQP.The second methodology is proposed as a fully optimization based method that is developed by refining a previously proposed PSO based method.The performance of the MPC-SQP method has been compared with the PSO based control method and a kinematic model based control method.In the simulation, the path tracking results have proved the superior performance of the MPC-SQP based controller.In the cases when computing times are considered, the PSO based controller offers more efficiency and better stability in computing compared with MPC-SQP based controller.

Figure 1 .
Figure 1.The 4WS4WD (four wheel steering and four wheel drive) vehicle.

Figure 4 .
Figure 4. Tire model.(a) Top view of the tire; (b) lateral view of the tire.

Figure 6 .
Figure 6.Flowchart of the vehicle control system.

Figure 12 .
Figure 12.Computing efficiency comparison.(a) computing time; (b) box plot of computing time.
where x d (k), y d (k) and u d (k) are the state vector, output vector and input vector, respectively.The coefficient matrices A d , B d and C d are updated after discretization.To eliminate undesirable oscillations, embedded integrator vectors ∆x d

Table 1 .
Parameters used in the simulation.