Tuning Parameters of the Fractional Order PID-LQR Controller for Semi-Active Suspension

: In order to further improve the control effect of proportion integral differential (PID) control and linear quadratic regulator (LQC) control, and improve vehicle ride comfort and enhance body stability, the 7 DOF semi-active suspension model was established, and the fractional order PI λ D µ -LQR controller was designed by combining fractional order PI λ D µ control theory and LQR control theory. The semi-active suspension model in this paper is more complex, and there are many parameters in the controller. The optimal weighting coefﬁcient of 12 vehicle smoothness evaluation indicators and parameters K p , K i , K d , λ and µ in the controller were founded by NSGA-II algorithm. After optimization, the optimized parameters were brought into the controller for random pavement simulation. Compared to the traditional passive suspension, fractional order PI λ D µ individual control and LQR separate control, the simulation results show that the effect of fractional order PI λ D µ -LQR control is very signiﬁcant. The evaluation index of vehicle smoothness has been signiﬁcantly improved, and the use of fractional order PI λ D µ -LQR control has signiﬁcantly improved the working performance of the suspension and improved the smoothness of the vehicle. At the same time, the adjusting force output of the actuator is very balanced, which inhibits the roll of the body and improves the anti-roll performance. After simulation, the excellent performance of the designed fractional PI λ D µ -LQR controller was veriﬁed, and the introduced NSGA-II algorithm played an important role in the controller parameter tuning work, which shows that the fractional order PI λ D µ -LQR controller and NSGA-II algorithm cooperate with each other to achieve good control effects.


Introduction
The research frontier for intelligent semi-active suspensions is the design of the controller, and many researchers are committed to developing a semi-active suspension controller with excellent performance to improve the working performance of the suspension and improve the smoothness of the vehicle.Ref. [1] established a 1/4 semi-active suspension model and designed a controller based on deep reinforcement learning.Based on state feedback and pre-aim feed-forward, a controller for state feedback and pre-sightforward feed-forward was designed by [2].Ref. [3] designed a semi-active suspension based on hydropneumatic mechanism and a semi-active suspension controller with a linear quadratic gaussian (LQG) optimal control algorithm.Ref. [4] developed a cloud-based parameter-adjustable adaptive semi-active suspension control method, while AMESim software was used by [5] to model the Continuous Damping Control damper adjustable damping shock absorber, analyze the dynamics of a 2-DOF semi-active suspension, and propose a genetic algorithm to improve the fuzzy rule.The control effect of the semi-active suspension controller depends largely on the parameters in the controller; for example, the control effect of the LQR controller depends on the weighting coefficient of various performance indicators, the performance of the PID controller mainly depends on K p , K i and K d parameters, and the control effect of the fuzzy controller depends on the fuzzy rules formulated, so the key to the design of the controller is the tuning of the controller parameters.For this reason, many scholars have conducted research on the tuning of parameters in the controller, and proposed different parameter tuning methods for different control modes.Ref. [6] proposed a locust optimization algorithm (GOA) and used it to adjust the parameters in the PID controller.Ref. [7] used the artificial fish swarm algorithm (AFSA) to determine the weighting coefficient in the LQR controller.Ref. [8] designed the semi-active suspension fuzzy controller and formulated the fuzzy rules of the controller, while [9] combined fuzzy control algorithm and particle swarm optimization algorithm to adjust P, I, D parameters.
The parameter adjustment in the suspension controller is highly dependent on the algorithm, and the appropriate optimization algorithm can not only improve the control effect of the controller, but also greatly improve the comfort and driving smoothness of the car.Different suspension controllers have been designed which use particle swarm algorithm, fuzzy control algorithm, deep learning algorithm, etc. to adjust the parameters in the controller, which also have good control effects [6][7][8][9].In fact, according to the "no free lunch" (NFL) theorem, for the same optimization problem, in the limited search space, different optimization algorithms can achieve different optimization effects, and no other algorithm can be better than the linear enumeration method of the search space or the pure random search algorithm [10,11].Orosz et al. [12] also mentioned that due to the "no free lunch" theorem of mathematical optimization, no one evolutionary optimization algorithm can outperform the others if all possible problems are averaged, and their article outlines the optimization techniques widely used in electric motors, summarizing the challenges and open problems of robust design optimization in applications and the development prospects of emerging technologies.
The NSGA-II algorithm belongs to a pure random search algorithm, whose principle is to imitate the "survival of the fittest" law of nature, with a strong gradient search ability.Although the NSGA-II algorithm is not universal to all optimization problems, for the semi-active suspension control system parameter optimization problem, it has a good optimization effect, which can be verified from the work of other researchers.For example, [13] constructed a rigid-flex coupling vehicle model of flexible control arm and torsion beam for vehicle dynamics simulation.The optimization goal is to control the total weight of the arm and torsion beam, ride comfort and handling stability performance indicators.Fatigue life, stiffness and modal frequency of control arms and torsion beams were used as constraints.Then, the NSGA-II was used to perform multi-objective optimization of the control arm and torsion beam to determine the lightweight scheme.Ref. [14] proposed to optimize the fuzzy logic controller (FLC) active seat suspension applied to articulated truck semi-trailer seats by using the NSGA-II algorithm, and the two objective functions optimized are seat vertical acceleration and controller control force to improve the seating comfort.Ref. [15] established a 1/4 semi-active suspension model and proposed to optimize spring rate and suspension damping by using the NSGA-II algorithm to improve vehicle comfort and ride.
The suspension model established by [13][14][15] is mainly a relatively simple 1/4 suspension model; the 1/4 suspension model can only consider the suspension performance indicators at a single axle, and it is difficult to comprehensively evaluate the performance of the whole vehicle, whereas the full vehicle model can take into account the force and movement of the vertical, pitch, roll and other degrees of freedom during the driving process of the car, optimize more indicators, evaluate the vehicle more comprehensively, and be more conducive to the study of the smoothness of the car.Secondly, because fractional order PI λ D µ control has higher design freedom, the design of fractional order PI λ D µ controller is more flexible than integer order PID controller [16,17].Therefore, according to the fractional order PI λ D µ control theory and LQR control theory, a better fractional order PI λ D µ -LQR controller is designed.Finally, the NSGA-II algorithm is improved on the basis of the traditional genetic algorithm (GA) and the first generation of non-dominated ranking genetic algorithm (NSGA), which overcomes many shortcomings of the traditional algorithm and is especially suitable for solving complex multi-objective optimization problems.The semi-active suspension system in this paper belongs to a very complex multi-objective optimization problem.Moreover, in the above literature, the two-dimensional plane representation method was used to explain the concepts of non-dominated sorting and crowding sorting of NSGA-II.For the case of many optimization goals, the concept of non-dominated sorting and crowding sorting is explained better by the three-dimensional spatial representation method, and the three-dimensional spatial solution is obtained, which is easier to understand and more clearly illustrates the problem.Therefore, this paper applies the NSGA-II algorithm to the semi-active suspension control system to solve the parameter tuning problem of the semi-active suspension fractional order PI λ D µ -LQR controller.
In summary, the 7-DOF semi-active suspension model is established, and the fractional order PI λ D µ control theory and LQR control theory are combined to design the fractional order PI λ D µ -LQR composite controller.The NSGA-II algorithm is used to find the optimal performance indicator weighting coefficient and optimal parameter K p , K i and K d , λ and µ in the fractional-order PI λ D µ -LQR controller.

Mathematical Model of Semi-Active Suspension
Taking the center of mass of the vehicle as the coordinate origin, the vehicle coordinate system is established with the longitude, lateral and vertical directions of the vehicle as the X axis, Y axis and Z axis, respectively, as shown in Figure 1.The vertical displacement of the body occurs along the Z axis, rotation around the Y axis is the pitch, rotation around the X axis is the roll, and the vertical displacement of the four unsprung masses along the Z axis are taken as the 7 DOF in the process of vehicle driving.
Electronics 2023, 12, x FOR PEER REVIEW 3 of 25 ranking genetic algorithm (NSGA), which overcomes many shortcomings of the traditional algorithm and is especially suitable for solving complex multi-objective optimization problems.The semi-active suspension system in this paper belongs to a very complex multi-objective optimization problem.Moreover, in the above literature, the two-dimensional plane representation method was used to explain the concepts of non-dominated sorting and crowding sorting of NSGA-II.For the case of many optimization goals, the concept of non-dominated sorting and crowding sorting is explained better by the threedimensional spatial representation method, and the three-dimensional spatial solution is obtained, which is easier to understand and more clearly illustrates the problem.Therefore, this paper applies the NSGA-II algorithm to the semi-active suspension control system to solve the parameter tuning problem of the semi-active suspension fractional order PI λ D µ -LQR controller.
In summary, the 7-DOF semi-active suspension model is established, and the fractional order PI λ D µ control theory and LQR control theory are combined to design the fractional order PI λ D µ -LQR composite controller.The NSGA-II algorithm is used to find the optimal performance indicator weighting coefficient and optimal parameter Kp, Ki and Kd, λ and µ in the fractional-order PI λ D µ -LQR controller.

Mathematical Model of Semi-Active Suspension
Taking the center of mass of the vehicle as the coordinate origin, the vehicle coordinate system is established with the longitude, lateral and vertical directions of the vehicle as the X axis, Y axis and Z axis, respectively, as shown in Figure 1.The vertical displacement of the body occurs along the Z axis, rotation around the Y axis is the pitch, rotation around the X axis is the roll, and the vertical displacement of the four unsprung masses along the Z axis are taken as the 7 DOF in the process of vehicle driving.In Figure 1, z s is the vertical displacement of the body center of mass; ϕ and θ represent the body pitch angle and the body roll angle, respectively; m w f 1 , m w f 2 , m wr1 , m wr2 represent the vehicle unsprung mass of the right front, left front, left rear and right rear, respectively; k w f 1 , k w f 2 , k wr1 , k wr2 represent the stiffness coefficients of right front tire, left front tire, right rear tire and left rear tire, respectively; k s f 1 , k s f 2 , k sr1 , k sr2 represent the stiffness coefficients of right front suspension, left front suspension, right rear suspension and left rear suspension coefficient, respectively; c s f 1 , c s f 2 , c sr1 , c sr2 represent the damping coefficient of right front suspension, left front suspension, right rear suspension and left rear suspension, respectively; q f 1 , q f 2 , q r1 , q r2 represent the vertical displacement of the road input to the four wheels; z w f 1 , z w f 2 , z wr1 , z wr2 represent the vertical displacement of the four wheels; z s f 1 , z s f 2 , z sr1 , z sr2 represent the vertical displacement of the body above the four wheels; F w f 1 , F w f 2 , F wr1 , F wr2 represent the right front suspension, left front suspension, right rear suspension and left rear suspension output force, respectively; a, b is the distance from the body center of mass to the front and rear wheel axle; L l , L r is the distance from the body center of mass to the center lines of the left and right wheel, respectively.The vehicle mass is expressed as m s , the pitch moment of inertia is expressed as I sy , the roll moment of inertia is expressed as I sx , the road roughness coefficient is expressed as G 0 , the experimental vehicle speed is expressed as u, and the lower stop frequency is expressed as f 0 .
By analyzing the force and moment of the vehicle's sprung and unsprung masses, the following differential equations of motion can be established [18][19][20][21]: The balance equation of the Z axis vertical force at the center of mass of the body is shown in Formula (1): ( The body roll moment balance equation is shown in Formula (3): Four unsprung mass vertical force balance equations are shown in Formula (4): When the variation range of body pitch angle and roll angle is small enough, the displacement of the suspension mass endpoints above the four wheels can be expressed as Formula (5): The state space expression (6) can be obtained by combining the above Formulas (1)-( 5): .
In Formula (6), matrix A is the system matrix, matrix B is the control matrix, matrix C is the output matrix, matrix D is the transfer matrix, matrix E is the disturbance matrix, X is the state variable, Y is the output variable, U is the control variable, and W is the perturbation variable.
A total of 18 variables including body vertical displacement, body pitch angle, body roll angle, vertical displacement of four unsprung masses, road input at four wheels, body vertical speed, body pitch angle speed, body roll angle speed and four vertical speeds of unsprung masses are selected as the state variation of the system, namely: A total of 11 variables including body vertical acceleration, body pitch angle acceleration, body roll angle acceleration, four suspension dynamic travel variables and four tire dynamic displacement variables are selected as the output variables of the system, namely: Four damping adjustment force variables, F cu f 1 , F cu f 2 , F cur1 , F cur2 , are selected as the components of the control vector U, namely: The vertical speed of the road surface at the four wheels, • q f 1 , • q f 2 , • q r1 , • q r2 , was selected as the component of the disturbance vector W, namely:

Pavement Simulation Model
The filtered white noise input is used as the pavement simulation model, and the time-domain expression of the vertical speed of the pavement is shown in Formula (11): where w 1 , w 2 , w 3 , w 4 is the unit Gaussian white noise with an average value of 0 and an intensity of 1, which can be directly invoked from the library when building the simulation model with Simulink.The pavement model established by Simulink is shown in Figure 2.
In Figure 2, G 0 represents the pavement unevenness coefficient of class B, u represents the speed of the vehicle, f 0 represents the frequency under the pavement, q f 1 (t) represents the vertical displacement of the road surface of the front wheel, • q f 1 represents the vertical speed of the road surface of the front wheel, and w 1 is the Gaussian white noise input to the front wheel road.
where 1 2 3 4 , , , w w w w is the unit Gaussian white noise with an average value of 0 an intensity of 1, which can be directly invoked from the library when building the simula model with Simulink.The pavement model established by Simulink is shown in Figure 2 + + ∑ f represents the frequency under the pavement, 1 ( ) resents the vertical displacement of the road surface of the front wheel, 1 f q g represent vertical speed of the road surface of the front wheel, and 1 w is the Gaussian white n input to the front wheel road.
Since the road conditions at each wheel of the vehicle do not differ much during driving process, the road model at the four wheels is the same as that at the front w 1.According to the harmonic superposition method, references [22,23], a B class 3D model with a length of 400 m and a width of 4 m is established as shown in Figure 3. Since the road conditions at each wheel of the vehicle do not differ much during the driving process, the road model at the four wheels is the same as that at the front wheel 1.According to the harmonic superposition method, references [22,23], a B class 3D road model with a length of 400 m and a width of 4 m is established as shown in Figure 3. x, y and z in Figure 3 are the length, width and rough surface thickness of th road, respectively

Fractional Order PI λ D µ Controller
Fractional order PI λ D µ control is error control.In the vehicle driving process, th hicle indicator people are most concerned about is the comfort of the body.Thus, the b vertical acceleration as the control object of fractional order PI λ D µ control, and 0 a expected value can achieve good control effect, with the actual output value of the b vertical acceleration as feedback.The difference between the feedback value and th pected value of the body's vertical acceleration is used as the error signal.The co principle of fractional order PI λ D µ control is as follows [24] (12): In Formula ( 12), ( ) e t is the error signal, ( ) is the output of the fractional o PI λ D µ controller, Kp is the scale coefficient, Ki is the integral coefficient, Kd is the differe coefficient, λ is the integral order, and µ is the differential order, where the λ and µ ca x, y and z in Figure 3 are the length, width and rough surface thickness of the 3D road, respectively.

Designing the Fractional Order PI
Fractional order PI λ D µ control is error control.In the vehicle driving process, the vehicle indicator people are most concerned about is the comfort of the body.Thus, the body vertical acceleration as the control object of fractional order PI λ D µ control, and 0 as the expected value can achieve good control effect, with the actual output value of the body vertical acceleration as feedback.The difference between the feedback value and the expected value of the body's vertical acceleration is used as the error signal.The control principle of fractional order PI λ D µ control is as follows [24] (12): In Formula ( 12), e(t) is the error signal, U PI λ D µ (t) is the output of the fractional order PI λ D µ controller, K p is the scale coefficient, K i is the integral coefficient, K d is the differential coefficient, λ is the integral order, and µ is the differential order, where the λ and µ can be any real number.
The basic structure block diagram of a fractional order PI λ D µ controller is shown in Figure 4.
In Formula ( 12), ( ) e t is the error signal, ( ) is the output of the fractional order PI λ D µ controller, Kp is the scale coefficient, Ki is the integral coefficient, Kd is the differential coefficient, λ is the integral order, and µ is the differential order, where the λ and µ can be any real number.The basic structure block diagram of a fractional order PI λ D µ controller is shown in Figure 4.The fractional order PI λ D µ control principle in Figure 4 is described as using the corresponding sensor of the vehicle to obtain the output variable y(t), and the error signal e(t) between the expected value r(t) and the output variable y(t) is used as the input of the fractional order PI λ D µ controller as the output of the controller.Fractional order PI λ D µ control extends integer integration and differentiation to fractional order when performing integration operations and differential operations, and the parameter adjustment is more flexible, which is better than the traditional integer order PID control method.
Parameter K p can adjust the rise time and adjustment time of the system, mainly affecting the response speed of the system, parameter K i ; λ can adjust the overshoot of the system, mainly affecting the adjustment accuracy of the system, parameter K d ; µ can adjust the number of oscillations of the system, mainly affecting the severity of the oscillation of the system.With reference to [25], the fractional order PI λ D µ controller is designed according to the Oustaloup algorithm.

LQR Controller
The 18 vehicle state variables in Formula (7) are taken as the feedback signals of the LQR controller.The controller uses the feedback signals as a reference to issue control instructions to the actuator, and the actuator outputs appropriate adjusting force.The control principle of LQR control is shown in Formula (13): where U LQR (t) is the output of the LQR controller, and K is the gain of the controller.Gain K can be obtained from the following formula: P is obtained from the following Riccati equation: Formulas ( 13) and ( 14) can be solved using MATLAB commands: , where: In Formula (16), Q is the weight matrix of 11 performance indicators, and R is the weight matrix of four actuator regulating forces, respectively expressed as follows: Q = diag(q 1 , q 2 , q 3 , q 4 , q 5 , q 6 , q 7 , q 8 , q 9 , q 10 , q 11 ) R = diag(q 12 , q 12 , q 12 , q 12 ) ( where: q 1 is the weighting coefficient of vertical acceleration at the body center of mass, q 2 is the weighting coefficient of pitch angle acceleration, q 3 is the weighting coefficient of roll angle acceleration, q 4 , q 5 , q 6 , q 7 are the weighting coefficients of four suspension dynamic travel variables, q 8 , q 9 , q 10 , q 11 are the weighting coefficients of four wheel dynamic displacement variables, and q 12 is the weighting coefficient of four actuator output forces.

Fractional Order PI λ D µ -LQR Controller
The fractional order PI λ D µ controller and the LQR controller are integrated into the fractional order PI λ D µ -LQR controller, and the control model is shown in Figure 5.

Fractional Order PI λ D µ -LQR Controller
The fractional order PI λ D µ controller and the LQR controller are integrated into the fractional order PI λ D µ -LQR controller, and the control model is shown in Figure 5. Figure 5 is the core of the entire control system.The fractional order PI λ D µ -LQR controller is composed of the fractional order PI λ D µ controller and the LQR controller, wherein the expected value ( ) 0 r t = of the fractional order PI λ D µ control, the controlled ob- ject is the vertical acceleration of the body, and K is the optimal feedback gain of the LQR controller.
The control principle of the fractional order PI λ D µ -LQR controller is as follows [26,27] (18): where is the comprehensive regulating force of the four actuators.

Designing the Objective Function
The process of finding an optimal solution is almost entirely dependent on the objective function, so determining the appropriate objective function is key to finding the global optimal solution.It is hoped that the body vertical acceleration, pitch angle acceleration and side inclination acceleration, as well as four suspension dynamic travel and four tire dynamic displacement indicators can be optimized.Taking the performance indicators of passive suspension as a reference, the optimization objectives as shown in Formula (19) are designed [28,29]: ( 1,2,3 10) where  Figure 5 is the core of the entire control system.The fractional order PI λ D µ -LQR controller is composed of the fractional order PI λ D µ controller and the LQR controller, wherein the expected value r(t) = 0 of the fractional order PI λ D µ control, the controlled object is the vertical acceleration of the body, and K is the optimal feedback gain of the LQR controller.
The control principle of the fractional order PI λ D µ -LQR controller is as follows [26,27] (18): where U PI λ D µ −LQR (t) is the comprehensive regulating force of the four actuators.The process of finding an optimal solution is almost entirely dependent on the objective function, so determining the appropriate objective function is key to finding the global optimal solution.It is hoped that the body vertical acceleration, pitch angle acceleration and side inclination acceleration, as well as four suspension dynamic travel and four tire dynamic displacement indicators can be optimized.Taking the performance indicators of passive suspension as a reference, the optimization objectives as shown in Formula ( 19) are designed [28,29]:

NSGA-II Algorithm to Find the Parameters of the Controller
For multi-objective optimization problems, various optimization objectives may conflict with each other and need to be weighed among various indicators, so it is necessary to select a suitable multi-objective optimization algorithm for optimization.Genetic algorithm is a group search algorithm that simulates the biological evolution process, constantly searching and comparing the best solution and eliminating the difference solution.The standard genetic algorithm will discard the last individual in the iteration process.As a classic modern genetic algorithm, the NSGA-II algorithm uses the elite retention strategy to better retain the elite individual, and solves the problem that the last generation of individuals in the standard genetic algorithm is discarded in the iteration process [30].In addition, this algorithm reduces the computational complexity of the algorithm by using a fast non-dominated sorting method.The traditional genetic algorithm is mainly used to solve single-objective optimization problems, and the NSGA-II algorithm expands the application scope of traditional genetic algorithms, especially suitable for solving multiobjective optimization problems.The advantages of the NSGA-II algorithm are mainly reflected in three aspects: (1) Elite retention strategy.Elite retention strategies expand the scope of screening when producing the next generation of individuals by mixing parent and child individuals into new groups.The elite retention strategy ensures that excellent individuals can have a greater probability of being retained, ensuring that the population is genetically better after continuous reproduction; (2) Non-dominated sorting.NSGA-II generates a set of non-dominated solution sets by ranking individuals according to their domination by other individuals, which enables NSGA-II to optimize multiple objective functions at the same time and obtain a series of optimal solutions; (3) Crowding sorting.The NSGA-II algorithm maintains population diversity by adopting the concept of crowding distance, and selects a better solution set based on this.The crowding distance is used to measure the local density of individuals in the solution space, and the distance between individuals is taken into account in the selection, so as to ensure the diversity of solution sets.
The NSGA-II algorithm was chosen to solve the multi-objective problem, and the flow chart of the algorithm is shown in Figure 6.
The NSGA-II algorithm has more operation processes than the traditional genetic algorithm (GA), as shown in part 1 (elite retention strategy) and part 2 (non-dominated sorting and crowding sorting), and more operation processes in step 1 than the first generation of nondominated ranking genetic algorithm (NSGA).The specific operation of the NSGA-II algorithm is as follows.Firstly, the initial population of a certain size is randomly generated, and the first generation population is obtained through the selection, crossover and mutation operations of the genetic algorithm after non-dominating sorting.Secondly, starting from the second generation, the parent population and the offspring population are merged, and rapid nondominant ranking is carried out.The individuals in the high-ranking non-dominant layer are retained as the parent generation as the population elite, and the individuals in the low-ranking dominant layer are directly eliminated.The crowding of the individuals in each non-dominant layer is then calculated, and the appropriate individuals are selected to form a new parent population according to the non-dominant relationship and the crowding of the individual.Finally, new offspring populations are generated through the basic operation of genetic algorithms; this process continues until the running conditions for the end of the program are met.
In Figure 6, Part 1 represents the elite retention strategy of the population, and Part 2 represents the non-dominant sorting process and the crowded sorting process of the population.
A more detailed introduction to the NSGA-II algorithm is shown in Figure 7: In Figure 6, Part 1 represents the elite retention strategy of the population, and Part 2 represents the non-dominant sorting process and the crowded sorting process of the population.
A more detailed introduction to the NSGA-II algorithm is shown in Figure 7:

Crowding distance sorting
Low non-dominance level, eliminated In Figure 6, Part 1 represents the elite retention strategy of the population, and Part 2 represents the non-dominant sorting process and the crowded sorting process of the population.
A more detailed introduction to the NSGA-II algorithm is shown in Figure 7:

Crowding distance sorting
Low non-dominance level, eliminated As shown in Figure 7, the population reproduction process P t represents the parent population, Q t represents the offspring population, and R t is the new population produced after the merger of the parent and offspring populations.The population size is set to N, and the specific steps in Figure 7 are described as follows: Step 1: Merge the parent and offspring populations into a new population R t ; nondominate the new population R t according to the dominance relationship between individuals; and divide the sorted population R t into multiple Pareto levels (here the population is divided into 6 Pareto levels: F1, F2, . . .F6).
Step 2: First put the nondominant individual (F1) with a Pareto level of F1 into a new parent set P t+1 (the darkest color set in the figure), then put the individual with a Pareto level of F2 into the new parent set P t+1 (as shown in the set with the next color depth in the figure), and so on.
Step 3: If all individuals with grades F1 and F2 are placed in the new parent set P t+1 , and the number of individuals in the set P t+1 is less than the population size N, the crowding of all individuals at level F3 is calculated and all individuals are arranged in descending order according to crowding, and then all individuals with grades greater than F3 (F4, F5, F6) are eliminated.
Step 4: Place the individuals at level F3 into the new parent set one by one in the order arranged in step 2 until the number of individuals in the parent set equals the population size N, and the remaining individuals are eliminated.
The concept of dominance can be expressed in Figure 8: Step 1: Merge the parent and offspring populations into a new population t R ; nondominate the new population t R according to the dominance relationship between individuals; and divide the sorted population t R into multiple Pareto levels (here the population is divided into 6 Pareto levels: F1, F2, … F6).
Step 2: First put the nondominant individual (F1) with a Pareto level of F1 into a new parent set 1 t P + (the darkest color set in the figure), then put the individual with a Pareto level of F2 into the new parent set 1 t P + (as shown in the set with the next color depth in the figure), and so on.
Step 3: If all individuals with grades F1 and F2 are placed in the new parent set 1 t P + , and the number of individuals in the set 1 t P + is less than the population size N, the crowding of all individuals at level F3 is calculated and all individuals are arranged in descending order according to crowding, and then all individuals with grades greater than F3 (F4, F5, F6) are eliminated.
Step 4: Place the individuals at level F3 into the new parent set one by one in the order arranged in step 2 until the number of individuals in the parent set equals the population size N, and the remaining individuals are eliminated.
The concept of dominance can be expressed in Figure 8:  ,, f f f ); the best solution is that which makes the 3 objective functions  Figure 8 is a schematic diagram of the dominant concept of the NSGA-II algorithm, in the three-dimensional spatial solution consisting of three objective functions ( f 1 , f 2 , f 3 ); the best solution is that which makes the 3 objective functions f 1 , f 2 , f 3 smallest in the Pareto frontier.As shown in Figure 6, the four points A, B, C, D are on the same plane, and the plane ABCD and the plane f 1 o f 2 are parallel.The distance from point A and point D to plane f 1 o f 3 and plane f 1 o f 2 is equal, but the distance from point A to plane f 2 o f 3 is less than the distance from point D to plane f 2 o f 3 .At this time, the non-dominant level of solution A is higher than solution D. A is the non-dominated solution, and D is the dominant solution, but here, solution A dominates solution D; thus, it can be understood that solution A is better than solution D, and in the process of population change, individual A will be retained, while individual D will be eliminated.
The crowding comparison operator is used to compare among quickly ranked peers, thus placing individuals in the quasi-Pareto realm.When the objective function obtained is s, the crowding distance can be expressed as: In Formula ( 22), i d represents the crowding distance of point i, f i−1 j represents the j objective function value of point i − 1, and f i+1 j represents the j objective function value of point i + 1.
The concept of crowding can be expressed in Figure 9.
In Formula ( 22), The concept of crowding can be expressed in Figure 9.
Figure 9 is a schematic diagram of the congestion concept of the NSGA-II algorithm.In a three-dimensional spatial solution composed of three objective functions ,, f f f , the spatial distance between point i and the adjacent two points 1 i − , 1 i + can be represented by m n t ++, by the sum of the length n , width m , and height t of the cuboid composed of these three points.In the Pareto frontier solution, the larger the spatial distance of point i , the better the solution, and point i is in the highest non-dominant position and belongs to the best solution among all solutions [31].
The optimal solution of the multi-objective optimization problem is the solution of individual i with the largest crowding degree distance in the front solution.
The main parameters of NSGA-II algorithm are shown in Table 1.  Figure 9 is a schematic diagram of the congestion concept of the NSGA-II algorithm.In a three-dimensional spatial solution composed of three objective functions f 1 , f 2 , f 3 , the spatial distance between point i and the adjacent two points i − 1, i + 1 can be represented by m + n + t, by the sum of the length n, width m, and height t of the cuboid composed of these three points.In the Pareto frontier solution, the larger the spatial distance of point i, the better the solution, and point i is in the highest non-dominant position and belongs to the best solution among all solutions [31].
The optimal solution of the multi-objective optimization problem is the solution of individual i with the largest crowding degree distance in the front solution.
The main parameters of NSGA-II algorithm are shown in Table 1.After 600 iterations, the Pareto solution is obtained, where f i (i = 1, 2, 3......10) is the i optimization objective.The solutions in three-dimensional space between each of the three optimization objectives and the projection in three planes are obtained as shown in Figure 10.From Figure 10, it can be seen that the solutions between each of the three objective functions are centrally distributed in a certain region of three-dimensional space, and Pareto frontier solutions exist in all of them, indicating that all objective functions are optimized.Among the frontier solutions in the 10-dimensional space composed of 10 objective functions, a solution with the largest crowding degree distance can be obtained as the global optimal solution.
The 12 weight variables are obtained as follows: q 1 = 80003.49,q 2 = 1000000.00,q 3 = 398530.00q 4 = 1000000.00,q 5 = 920000.00,q 6 = 760000.00q 7 = 820000.00,q 8 = 962489.99,q 9 = 880000.00q 10 = 877389.99,q 11 = 460000.00,q 12 = 1.00 The five parameters of the fractional order PI λ D µ controller are obtained as follows: To test the optimized parameters K p , K i , K d , λ and µ and whether the stability of the system can be guaranteed, the zero-pole distribution diagram of the transfer function of the closed-loop system (fractional order PI λ D µ closed-loop control part: r(t) = 0 as the expected value input, y(t) = .. z s (t) as the output, as shown in Figure 5) is obtained after the identification of the control system by using the system identification toolbox of MATLAB, as shown in Figure 11.820000.00,962489.99,880000.00877389.99,460000.001.00 The five parameters of the fractional order PI λ D µ controller are obtained as follows: 1.098 2944.466 1.336 0.0001 1.561 To test the optimized parameters Kp, Ki, Kd, λ and µ and whether the stability of the system can be guaranteed, the zero-pole distribution diagram of the transfer function of the closed-loop system (fractional order PI λ D µ closed-loop control part: ( ) 0 r t = as the ex- pected value input, ..

( ) ( )
as the output, as shown in Figure 5) is obtained after the identification of the control system by using the system identification toolbox of MATLAB, as shown in Figure 11.As can be seen from Figure 11, the zeros and poles of the closed-loop system are completely distributed in the left half of the complex plane, indicating that the optimized parameter K p , K i , K d , λ and µ can ensure system stability.

Model Simulation
In order to highlight the control effect of fractional order PI λ D µ -LQR composite control, the simulation results are compared with those of passive suspension, the fractional order PI λ D µ control and the LQR control.The simulation time is 10 s, the simulation step length range is [0.01, 0.001], the fixed speed is 20 m/s, the road grade is branch road, and the value range of the road unevenness coefficient under this road condition is 5 × 10 −7 ~3 × 10 −5 , with a mean of 5 × 10 −6 .The time domain curve of each performance indicator is obtained, and the power spectral density curve of each performance indicator can be obtained by Fourier transform of the output time domain curve.The range of horizontal coordinates is 0-5 s. Figure 12 shows the time domain curve of each performance indicator: order PI D control and the LQR control.The simulation time is 10 s, the simulation step length range is [0.01, 0.001], the fixed speed is 20 m/s, the road grade is branch road, and the value range of the road unevenness coefficient under this road condition is 5 × 10 −7 ~ 3 × 10 −5 , with a mean of 5 × 10 −6 .The time domain curve of each performance indicator is obtained, and the power spectral density curve of each performance indicator can be obtained by Fourier transform of the output time domain curve.The range of horizontal coordinates is 0-5 s. Figure 12  As can be seen from Figure 12, the effect of fractious-order PI λ D µ -LQR composite control is more obvious than that obtained by passive suspension, fractious-order PI λ D µ control alone and LQR control alone.The peak value of the output performance indicator curve of the fractional order PI λ D µ -LQR control significantly decreases, but the roll angle acceleration slightly worsens.
Figure 13 shows the power spectral density curves corresponding to the time-domain As can be seen from Figure 12, the effect of fractious-order PI λ D µ -LQR composite control is more obvious than that obtained by passive suspension, fractious-order PI λ D µ control alone and LQR control alone.The peak value of the output performance indicator curve of the fractional order PI λ D µ -LQR control significantly decreases, but the roll angle acceleration slightly worsens.
Figure 13 shows the power spectral density curves corresponding to the time-domain curves of each performance indicator.Considering that the excitation of the branch road generally does not exceed 10 Hz, the frequency range of the horizontal coordinate is 0-15 Hz.As can be seen from Figure 13, compared with passive suspension, the fractional order PI λ D µ control and the LQR control, after adopting the fractional order PI λ D µ -LQR composite control, the power spectral density of body vertical acceleration, pitch angle acceleration, suspension dynamic travel and tire dynamic displacement significantly decreases in the frequency range of 0-5 Hz, while the change is not obvious in the frequency range of 5-15 Hz.This indicates that after active control, the actuator can well suppress the lowfrequency impact from the road surface, but the suppression effect on the high-frequency impact is not obvious.

Improvement of Performance Indicators
In order to more clearly reflect the advantages of the fractional order PI λ D µ -LQR composite control, the time-domain curves of each performance indicator are squared and the average root values are shown in Table 2.As can be seen from Table 2, after the fractional order PI λ D µ -LQR composite control is adopted, compared with the passive suspension, the other performance indicators are As can be seen from Figure 13, compared with passive suspension, the fractional order PI λ D µ control and the LQR control, after adopting the fractional order PI λ D µ -LQR composite control, the power spectral density of body vertical acceleration, pitch angle acceleration, suspension dynamic travel and tire dynamic displacement significantly decreases in the frequency range of 0-5 Hz, while the change is not obvious in the frequency range of 5-15 Hz.This indicates that after active control, the actuator can well suppress the low-frequency impact from the road surface, but the suppression effect on the highfrequency impact is not obvious.

Improvement of Performance Indicators
In order to more clearly reflect the advantages of the fractional order PI λ D µ -LQR composite control, the time-domain curves of each performance indicator are squared and the average root values are shown in Table 2.As can be seen from Table 2, after the fractional order PI λ D µ -LQR composite control is adopted, compared with the passive suspension, the other performance indicators are greatly optimized except for the slight deterioration of the roll angle acceleration.In order to display the improved condition of the optimized performance indicators more directly, a bar chart of the improvement of the indicators is shown in Figure 14: By analyzing the data, it can be found that the control effect of the fractional order PI λ D µ -LQR composite control is remarkable.All performance indicators are optimized to the greatest extent on the basis of the fractional order PI λ D µ control and the LQR control.

Influence of Actuator Output Force on Roll Angle Acceleration
By analyzing Formula (3), it can be found that while the left and right suspension By analyzing the data, it can be found that the control effect of the fractional order PI λ D µ -LQR composite control is remarkable.All performance indicators are optimized to the greatest extent on the basis of the fractional order PI λ D µ control and the LQR control.

Influence of Actuator Output Force on Roll Angle Acceleration
By analyzing Formula (3), it can be found that while the left and right suspension stiffness are the same, the roll angle acceleration indicator is only related to the dynamic travel of the four suspensions and the regulating force of the four actuators.Combined with Table 1, it is noted that the dynamic travel of the left and right suspension of the passive suspension is exactly the same, and the left and right roll torque is balanced, so the square mean root value of the roll angle acceleration of the passive suspension is very small.After active control is adopted, when the dynamic travel of left and right suspension is not much different, the acceleration indicator of side inclination is mainly affected by the adjustment force of four actuators.
When the weighting coefficients of the four actuators in the matrix R of Formula ( 17) are the same, the regulating force output by the four actuators is the balanced force; when the weighting coefficients of the four actuators are different, the regulating force output is the unbalanced force.The adjusting force output by the actuator in two cases is shown in Figures 15 and 16  In order to visually compare the influence of the adjusting force output by the four actuators on the roll angle acceleration, the time-domain curve of the roll angle acceleration is obtained as shown in Figure 17:  In order to visually compare the influence of the adjusting force output by the four actuators on the roll angle acceleration, the time-domain curve of the roll angle acceleration is obtained as shown in Figure 17: In order to visually compare the influence of the adjusting force output by the four actuators on the roll angle acceleration, the time-domain curve of the roll angle acceleration is obtained as shown in Figure 17: In order to visually compare the influence of the adjusting force output by the actuators on the roll angle acceleration, the time-domain curve of the roll angle acce tion is obtained as shown in Figure 17: Obviously, the balanced adjustment force output by four actuators can more e tively restrain the roll of the car body and improve the stability of the car body.

Conclusions
This paper introduces theoretical research on the analysis, modeling and simula of semi-active suspension control systems of vehicles.Firstly, through the force ana of the semi-active suspension dynamic model of the whole vehicle, the differential e tion of motion of the whole vehicle in the direction of seven degrees of freedom is der Obviously, the balanced adjustment force output by four actuators can more effectively restrain the roll of the car body and improve the stability of the car body.

Conclusions
This paper introduces theoretical research on the analysis, modeling and simulation of semi-active suspension control systems of vehicles.Firstly, through the force analysis of the semi-active suspension dynamic model of the whole vehicle, the differential equation of motion of the whole vehicle in the direction of seven degrees of freedom is derived, and it is expressed in the form of state space.Secondly, according to the fractional PI λ D µ control theory and LQR control theory, the fractional PI λ D µ -LQR controller was designed, and the simulation model of the fractional PI λ D µ -LQR control system was built by MATLAB/simulink.A non-dominant ranking genetic algorithm (NSGA-II) is then introduced, and the concepts of dominance, crowding and elite retention strategy of the NSGA-II algorithm are explained in detail.Because of the advantages embodied by the NSGA-II algorithm, it was used to optimize the weighting coefficient of 12 performance indicators and parameters K p , K i , K d , λ and µ in the fractional order PI λ D µ -LQR controller.Finally, the parameters after tuning are brought into the controller for simulation, and at the same time, in order to reflect the advantages of the proposed fractional PI λ D µ -LQR control method, compared with passive suspension, the fractional order PI λ D µ control and the LQR control, the comparison results show that the fractional PI λ D µ -LQR compound control is better than the fractional order PI λ D µ control and LQR control, and the performance indicators are improved to the greatest extent.Through simulation, it is found that the balanced adjustment force output by four actuators can effectively restrain the body roll.The fractional order PI λ D µ control and the LQR control are combined to give full play to the advantages of the two control methods and improve the control effect.

Figure 1 .
Figure 1.Seven degrees of freedom dynamic model of semi-active suspension.In Figure1, sz is the vertical displacement of the body center of mass; ϕ andθ represent the body pitch angle and the body roll angle, respectively; 1

Figure 1 .
Figure 1.Seven degrees of freedom dynamic model of semi-active suspension.

Figure 2 .
Figure 2. Pavement model of front wheel 1.In Figure 2, 0 G represents the pavement unevenness coefficient of class B, u re sents the speed of the vehicle, 0

Figure 4 .
Figure 4. Structure block diagram of fractional order PI λ D µ controller.Figure 4. Structure block diagram of fractional order PI λ D µ controller.

Figure 4 .
Figure 4. Structure block diagram of fractional order PI λ D µ controller.Figure 4. Structure block diagram of fractional order PI λ D µ controller.

Figure 5 .
Figure 5. Simulation model of semi-active suspension system.
the square mean root value of 10 performance indicators output by semi-active suspension, and is the square mean root value of 10 performance indicators output by passive suspension corresponding to the output of semiactive suspension.The 10 performance indicators are vertical acceleration of the body, pitch angle acceleration, four suspension dynamic travel variables and four tire dynamic

Figure 5 .
Figure 5. Simulation model of semi-active suspension system.

Figure 8
Figure 8 is a schematic diagram of the dominant concept of the NSGA-II algorithm, in the three-dimensional spatial solution consisting of three objective functions ( 1 2 3

13 f of and plane 12 f
smallest in the Pareto frontier.As shown in Figure 6, the four A, B, C, D are on the same plane, and the plane ABCD and the plane 12 f of are parallel.The distance from point A and point D to plane of is equal, but the distance from point A to plane 23 f of is less
the j objective function value of point

1 f
is used to represent the vertical acceleration indicator of the body; 2 f is used to represent the pitch angle acceleration indicator; 3 4 5 , , f f f and 6 f are used to represent the four suspension dynamic travel indicators; 7 8 9 , , f f f and 10 f are used to represent the four tire dynamic displacement indicators; (a), (b), (c), (d) represent the Pareto solution of each of the three optimization indicators and the projection of these solutions in the spatial Cartesian coordinate system; and the projection represents the local optimal solution between the corresponding two indicators.The local optimal solution between each of the three performance indicators in Figure 10 is circled by the red dashed line as the Pareto leading solution, and the global optimal solution is obtained in the local optimal solution.

Figure 11 .
Figure 11.Zero-pole distribution diagram of closed-loop system.Figure 11.Zero-pole distribution diagram of closed-loop system.

Figure 11 .
Figure 11.Zero-pole distribution diagram of closed-loop system.Figure 11.Zero-pole distribution diagram of closed-loop system.

Figure 12 .
Figure 12.(a) Time domain curve of vehicle body vertical acceleration; (b) time domain curve of pitch angle acceleration; (c) time domain curve of roll angle acceleration; (d) time domain curve of left front suspension dynamic travel; (e) time domain curve of left rear suspension dynamic travel; (f) time domain curve of right front tire dynamic displacement; (g) time domain curve of left rear tire dynamic displacement.

Figure 12 .
Figure 12.(a) Time domain curve of vehicle body vertical acceleration; (b) time domain curve of pitch angle acceleration; (c) time domain curve of roll angle acceleration; (d) time domain curve of left front suspension dynamic travel; (e) time domain curve of left rear suspension dynamic travel; (f) time domain curve of right front tire dynamic displacement; (g) time domain curve of left rear tire dynamic displacement.

Figure 13 .
Figure 13.(a) Power spectral density curve of the vertical acceleration of the vehicle body; (b) power spectral density curve of pitch angle acceleration; (c) power spectral density curve of roll angle acceleration; (d) power spectral density curve of left front suspension dynamic travel; (e) power spectral density curve of left rear suspension dynamic travel; (f) power spectral density curve of right front tire dynamic displacement; (g) power spectral density curve of the left rear tire dynamic displacement.

Figure 14 .
Figure 14.(a) Improvement of performance indicators of the fractional order PI λ D µ -LQR control compared with passive suspension; (b) performance improvement of the fractional order PI λ D µ -LQR control compared with fractional order PI λ D µ control; (c) fractional order PI λ D µ -LQR control compared with the LQR control.

Figure 14 .
Figure 14.(a) Improvement of performance indicators of the fractional order PI λ D µ -LQR control compared with passive suspension; (b) performance improvement of the fractional order PI λ D µcompared fractional order PI λ D µ control; (c) fractional order PI λ D µ -LQR control compared with the LQR control.

Figure 17 .
Figure 17.Influence of adjusting force of the actuator on roll angle acceleration.

Figure 17 .
Figure 17.Influence of adjusting force of the actuator on roll angle acceleration.
suspension dynamic travel variables and four tire dynamic displacement variables.In order to achieve better optimization results, a penalty function is given as in Equation (20):At that time, if RMS semi−active−i > RMS passive−i (i = 1, 2, 3 • • • 10), then: • • • 10) is the square mean root value of 10 performance indicators output by semi-active suspension, and is the square mean root value of 10 performance indicators output by passive suspension corresponding to the output of semi-active suspension.The 10 performance indicators are vertical acceleration of the body, pitch angle acceleration, four

Table 1 .
Main parameters of NSGA-II algorithm.

Table 1 .
Main parameters of NSGA-II algorithm.

Table 2 .
Average root value of each performance indicator.

Table 2 .
Average root value of each performance indicator.