Optimal Design of Agricultural Mobile Robot Suspension System Based on NSGA-III and TOPSIS

Abstract

The stability of vehicles is influenced by the suspension system. At present, there are many studies on the suspension of traditional passenger vehicles, but few are related to agricultural mobile robots. There are structural differences between the suspension system of agricultural mobile robots and passenger vehicles, which requires structural simplification and modelling concerning suspension of agricultural mobile robots. This study investigates the optimal design for an agricultural mobile robot’s suspension system designed based on a double wishbone suspension structure. The dynamics of the quarter suspension system were modelled based on Lagrange’s equation. In our work, the non-dominated sorting genetic algorithm III (NSGA-III) was selected for conducting multi-objective optimization of the suspension design, combined with the Technique for Order Preference by Similarity to the Ideal Solution (TOPSIS) to choose the optimal combination of parameters in the non-dominated solution set obtained by NSGA-III. We compared the performance of NSGA-III with that of other multi-objective evolutionary algorithms (MOEAs). Compared with the second-scoring solution, the score of the optimal solution obtained by NSGA-III increased by 4.92%, indicating that NSGA-III has a significant advantage in terms of the solution quality and robustness for the optimal design of the suspension system. This was verified by simulation in Adams that our method, which utilizes multibody dynamics, NSGA-III and TOPSIS, is feasible to determine the optimal design of a suspension system for an agricultural mobile robot.

1. Introduction

The smoothness of a vehicle during driving is one of its important evaluation indicators. When a vehicle is driven in irregular road conditions, uncomfortable vibrations may be transmitted to the driver. As unpaved roads make up the majority of an agricultural mobile robot’s working environment, considerable of vibration can be transmitted to the operating equipment carried by an agricultural mobile robot, reducing its working precision and shortening its lifespan. As an important system in a vehicle that has a vibration mitigation function, a properly designed suspension system can absorb some of the vibrations and reduce the impact on the driver or working equipment, and has therefore long been the subject of research and optimization by academics. There have been many studies on conventional vehicle suspensions. The optimal design of an agricultural robot suspension system is essentially a constrained optimization problem. Lagrange multiplier methods, evolutionary algorithms and machine learning are commonly applied to solve constrained optimization problems. These methods are utilized in many areas, such as path planning [1], lesion diagnosis [2,3,4], defect detection [5,6], structural design [7] and resource allocation [8,9]. There are very few studies on the optimal design of agricultural robot suspensions. Therefore, it is helpful and instructive to carry out research on the optimal design of the suspension of agricultural robots by investigating the existing research on the suspension optimization design of passenger cars.

Ciro Moreno Ramírez et al. [10] investigated the effect of two different structures of suspension systems on vehicle dynamics, and investigated the stability of the systems using root trajectory analysis methods. Li et al. [11] used the response surface method to establish an approximate model of the MacPherson suspension systems and verified the reliability of it, and then used the NSGA-II algorithm to optimize the key hard-point coordinates of the suspension. To optimize the design of the pertinent parameters, Qian and Jin [12] examined the performance of an SUV double wishbone suspension and applied the D-optimal method based on the response surface model. Guang Li et al. [13] used the NSGA-II algorithm to optimize the suspension parameters of high-speed locomotives. The locomotive dynamics model established by SIMPACK software was used to optimize the design of six parameters by the NSGA-II algorithm, and the three evaluation indexes designed were minimum indicators. None of the above studies established an accurate non-linear mathematical model, which may lead to problems such as inaccurate solution results or reduced solution efficiency in the solution process.

Chen et al. [14] established a mathematical model of the double wishbone suspension matching an in-wheel motor with a cosine matrix method. Zhang and Li [15] established a nonlinear mathematical model of a double wishbone independent suspension using a kinematic analysis method, and then used simulation analysis in Adams software and bench testing to demonstrate that the developed mathematical model accurately expressed the dynamic properties of the suspension. The nonlinear motion equation for flexible double wishbone suspension was developed by Abdelrahman et al. [16] using the concept of imaginary displacement, but the main focus of this study was on how flexible structures respond to various road unevenness, vehicle speed, and material damping coefficients. By utilising laser scanning to create 3D models of the suspension parts, and Adams to simulate their dynamics, Prastiyo and Fiebig [17] compared the benefits of linear and progressive double wishbone suspensions and found no discernible differences between them, which showed that it is more practical to use a traditional linear double wishbone suspension. A mathematical model of the suspension system was created by the aforementioned analysis, but further optimization of the current suspension system has not been done.

A mathematical model of an independent steering-suspension guidance mechanism was developed by Chen et al. [18] based on the theory of spatial mechanics, and the structural hard-point coordinates of the suspension were optimized using sensitivity analysis. Sancibrian et al. [19] developed a model containing a large number of structural parameters, mainly the lengths of the individual links, for the structure of a double wishbone suspension. Seven functional parameters were applied to evaluate the design parameters. This study was based on the gradient descent method to find the optimal solution. Shi et al. [20] investigated the issue of optimizing the MacPherson suspension’s hard-point coordinates, established the relationship between the hard-point coordinates of the suspension and the evaluation index based on a dynamics model and support vector regression (SVR) established in the Adams/Car software, and designed a double-loop multi-objective particle swarm algorithm to optimize the hard-point coordinates of the suspension. The results showed that the improved algorithm outperformed the traditional multi-objective particle swarm algorithm and the genetic algorithm. Totu and Alexandru [21] presented a comprehensive design method based on the least squares approach for the optimal design of an innovative racing car suspension system. This study used Adams/View and Adams/Insight to build a regression model and carry out the optimization design. Zhu et al. [22] provided a sliding mode control method for the equivalent two-degrees-of-freedom model and calibrate the unknown parameters in the equivalent two-degrees-of-freedom model by parameter identification for a quarter double wishbone suspension model. In this study, the equivalent two-degrees-of-freedom suspension model was constructed by parameter identification, allowing for the analysis of the suspension response without the need to create a dynamic model. The above studies did not consider the influence of suspension performance parameters (such as damping and stiffness) on suspension performance during the optimization process.

Issa and Samn [23] used the Harris Hawk Optimization (HHO) algorithm to optimize the damping and stiffness coefficients of passive suspension, which they simplified to a two-degrees-of-freedom model. The result showed that the optimized passive suspension performance was improved. Huang et al. [24] simplified the double wishbone suspension to a single-wheel two-degrees-of-freedom model and employed evolutionary algorithms to find its parameters. Based on the suspension equivalent two-degrees-of-freedom model and the seat-passenger equivalent eight-degrees-of-freedom model, Papaioannou and Koulocheris [25] improved the solution speed and quality of evolutionary algorithms by dividing the optimization targets into primary optimization objectives and auxiliary ones. Drehmer et al. [26] modeled the motion of the whole vehicle by introducing the influence on the driver’s seat to establish an 8-degrees-of-freedom motion model. Particle swarm optimization algorithms and sequential quadratic programming algorithms were used to optimize the damping and stiffness of the four independent suspensions of the complete vehicle under different road conditions. Gobbi et al. [27] used a 2-degrees-of-freedom linear model to analytically describe the dynamic behaviour of a vehicle travelling on a randomly contoured road, and respectively optimized the damping and stiffness of the passive suspension and the damping, stiffness and controller gain of the active suspension based on multi-objective planning theory and robust design theory. Kwon et al. [28] developed a mathematical model of a hydro-pneumatic suspension system and a whole vehicle model for the design of hydro-pneumatic suspension in heavy vehicles. This study developed an agent model to reduce the computational effort in the optimal design process, and used it for multi-objective optimization solutions. Zheng et al. [29] designed an active, tuned inertial damper (TID) suspension systems based on a combination of active actuator and inertializer, and proposed a parameter optimization method based on an analytical solution. Stability algebraic analysis was carried out using Hurwitz’s criterion and it was verified that the parameters obtained with this method could guarantee the stability of the suspension. Yang et al. [30] used a decomposition-based multidisciplinary optimization approach for the optimal design of passive suspension systems based on an equivalent two-degrees-of-freedom model. In terms of the selection of design parameters, this study focused on the refinement of the parameters affecting damping and stiffness in the suspension system, such as the number of coils of the spring, the spring diameter and the piston diameter, among other structural parameters. Li et al. [31] proposed a dimensionless hybrid index based on safety probability to evaluate suspension performance for the optimal design of energy-harvesting suspension systems, transformed the multi-objective optimization problem into a single-objective optimization problem, and then used a genetic algorithm to solve the multi-parameter optimization problem. Truong and Dao [32] proposed a hybrid HNSGA-III & MOPSO algorithm based on the MOPSO algorithm and NSGA-III for the optimal design of the stiffness and damping of a powertrain suspension model. HNSGA-III & MOPSO algorithms showed better efficiency and solution quality than MOPSO and NSGA-III in this study. Grotti et al. [33] proposed a multi-objective archive-based Quantum Particle Swarm Optimizer (MOQPSO) algorithm for the optimization of suspension systems. The results were compared with NSGA-II algorithm and COGA-II algorithm. It was shown that the MOQPSO algorithm proposed in this study could obtain better non-dominated set solutions than the NSGA-II and COGA-II algorithms. Bingul and Yildiz [34] carried out a multi-objective optimized design of a non-linear suspension system for an electric vehicle based on the NSGA-II algorithm. Three suspension systems, including passive suspension, active suspension based on PD control, and active suspension based on FL control, were optimized. This study designed seven minimum optimization objectives to evaluate the degree of influence of the active suspension on the driver under road vibration conditions. Prasad et al. [35] optimized and improved the stiffness and PID control parameters of a semi-active suspension system in a quarter suspension model. They proposed a fast convergence optimization algorithm for saving computational costs. Compared with NSGA-II algorithm, the geometry-inspired genetic algorithm proposed in this study converges quickly under small population conditions. The above studies did not consider the influence of suspension structural parameters on suspension performance during the optimization process.

Gialleonardo et al. [36] used the NSGA-II algorithm for optimal design of control strategy parameters for active suspensions. They designed six minimum targets to evaluate the parameters of three different active control strategies and to compare the performance of them. Xu et al. [37] used a multi-objective particle swarm algorithm to optimize the design of key parameters of a straw back-throwing device. This study showed that the multi-objective evolutionary algorithm is effective and feasible for the optimal design of mechanical structures. Jiang et al. [38] used the Kriging model and NSGA-II algorithm for the multi-objective lightweight design of two parts for a control arm and torsion beam, which are widely used in passenger car suspensions. Dang et al. [39] used NSGA-III to optimize the parameters of the rotational speed and opening length of a fertiliser spreader with the objectives of accuracy, uniformity, adjustment time, and crushing rate. They chose the single-objective evolutionary algorithm GA and the multi-objective evolutionary algorithm MOEA-D-DE for comparison. The evaluation metrics showed that NSGA-III has significant advantages in solving this kind of problem. Li et al. [40] proposed a multi-objective optimal control method for an active suspension system aimed at solving the negative vibration problem generated by the in-wheel motor of an electric vehicle. An integrated model considering the electromechanical coupling between the electromagnetic excitation of the motor and the transient dynamics of the vehicle was established and developed. The Pareto solution set for the optimal parameters of the active suspension system was solved using the particle swarm optimization algorithm. Chen et al. [41] conducted a study on the multi-objective optimization problem of high-speed train suspension systems using the NSGA-II algorithm to optimize four suspension parameters with three minimum objectives. The research above did not offer a clear strategy on how to select the optimal solution from the Pareto solution set. Some studies [13,31,32,33,34,36,37,41] lacked comparisons with other multi-objective optimization algorithms.

To address the above problems, our study proposes an optimal design method for the suspension of the agricultural mobile robot based on NSGA-III and TOPSIS. First, a dynamics model of the suspension system of the agricultural mobile robot was established based on the first-class Lagrange equation, which includes structural and performance parameters of the suspension system. Four minimum objectives were established for evaluating the performance of the suspension system. Second, the dynamics model was multi-objective optimized using eight multi-objective optimization algorithms including NSGA-III. The optimal solution in the Pareto solution set obtained by the multi-objective optimization algorithm was selected using the TOPSIS method. Finally, the optimal solutions set obtained by the multiple evolutionary algorithms was scored and ranked using the TOPSIS method. The results show that the NSGA-III algorithm obtained the optimal parameter combinations in this study. The top five parameter combinations were simulated by Adams software to verify the feasibility and effectiveness of the method.

This work aimed to find the optimum design of an agricultural robot passive suspension system using the NSGA-III algorithm and TOPSIS method. The main contributions and innovations of this study are as follows:

• A suspension system dynamics model based on the first type of Lagrangian equations containing structural and performance parameters was developed to facilitate the optimization process for the overall optimal design of both types of parameters.
• NSGA-III was used to compute the multi-objective optimization problem, resulting in a more significant diversity and convergence of feasible solutions. A comparison with various multi-objective evolutionary algorithms was also made, showing the advantages of NSGA-III in solving this problem.
• The optimal solution in the Pareto solution set was selected by scoring it using the TOPSIS method. In the comparison to multiple algorithms, the optimal solution set composed was again scored and ranked using the TOPSIS method, and the scoring was more reasonable.

2. Materials and Methods

2.1. Suspension Dynamics Model

2.1.1. Mechanical Model Simplification

The study was carried out concerning the suspension system of the agricultural mobile robot. The quarter model of the suspension system is shown in Figure 1, where the slider is fixed to the ABC bar, indicating that the robot’s frame is able to move freely in the vertical direction. The BD and CF bars indicate the wishbone in the double wishbone suspension, the AE bar indicates the spring damped shock absorber, the DEF bar indicates the steering system, and G indicates the wheels.

For analysis of the quarter suspension system model shown in Figure 1a, the BD bar and CF bar motion forms and motion parameters are the same, so a simplified model of the quarter suspension system can be obtained, as shown in Figure 1b.

2.1.2. Dynamic Model

In Figure 1b, H is the BC midpoint. For the motion analysis of the simplified model shown in Figure 1b, the generalised coordinates can be taken as $q=x,x_E,y_E,y_H,y_A$ (m). The active member is the frame AH bar, the wheel bracket E is constrained by the ground, the centre of mass of each member is $s_i(i=2,4,5)$, the mass of each member is $m_i(i=2,4,5)$ (kg), the velocity of motion at the centre of mass of each member is $v_i(i=2,4,5)$ (ms⁻¹), the angular velocity at the centre of mass of each member is ${\omega }_i(i=2,4)$ (rads⁻¹), the spring deformation is x, the wheel bracket E is subjected to the support reaction force F(N) from the ground, and its motion analysis is shown in Figure 2.

The dynamics of the suspension system can be modelled on the basis of the first type of Lagrange equation as follows:

$$\frac{d}{dt}(\frac{\partial L}{\partial {\dot{q}}_i})-\frac{\partial L}{\partial q_i}+{\displaystyle \sum _{a=1}^4{\lambda }_a\frac{\partial C_a}{\partial q_i}}=Q_i$$

(1)

where L is the Lagrangian function of the system, $L=E_k-E_p$.

$E_k$ is the kinetic energy of the system.

$E_p$ is the potential energy of the system.

$q_i$ are generalized coordinates.

$Q_i$ for broad forces.

The structural constraints are

$$C_1:{(x_E-x_H)}^2+{(y_E-y_H)}^2=l_2{}^2$$

(2)

$$C_2:y_A-y_H=l_{ah}$$

(3)

$$C3:{(x_E-x_A)}^2+{(y_E-y_A)}^2={(l-x)}^2$$

(4)

$$C_4:y_E=Y_E$$

(5)

The Lagrangian function of the system is

$$\begin{array}{cc}\hfill & L=-\frac{g}{2}\left(m_2(y_E+y_H)+m_4(y_A+y_E)+m_5(y_A+y_H)\right)-\frac{1}{2}k{x}^2\hfill \\ \hfill & +\frac{m_4{(l-x)}^2{\left(\frac{{\dot{y}}_A-{\dot{y}}_E}{x_E}-\frac{{\dot{x}}_E(y_A-y_E)}{{x_E}^2}\right)}^2}{6{\left(\frac{{(y_A-y_E)}^2}{{x_E}^2}+1\right)}^2}+\frac{{l_2}^2{m}_2{\left(\frac{{\dot{y}}_E-{\dot{y}}_H}{x_E}-\frac{{\dot{x}}_E(y_E-y_H)}{{x_E}^2}\right)}^2}{6{\left(\frac{{(y_E-y_H)}^2}{{x_E}^2}+1\right)}^2}\hfill \\ \hfill & +\frac{1}{2}{m}_2\left(\frac{1}{4}{{\dot{x}}_E}^2+\frac{1}{4}{\left({\dot{y}}_E+{\dot{y}}_H\right)}^2\right)+\frac{1}{2}{m}_4\left(\frac{1}{4}{{\dot{x}}_E}^2+\frac{1}{4}{\left({\dot{y}}_A+{\dot{y}}_E\right)}^2\right)+\frac{1}{8}{m}_5{\left({\dot{y}}_A+{\dot{y}}_H\right)}^2\hfill \end{array}$$

(6)

Generalized forces $Q_i$

$$\left.\begin{array}{l}{Q}_1=-c\dot{x}\\ Q_2=0\\ Q_3=(m_2+m_4+m_5+m_b)g\\ Q_4=0\\ Q_5=0\end{array}\right\}$$

(7)

According to Equations (1)–(7), we can obtain

$$M_{11}{\dot{x}}_E{\dot{y}}_A+M_{12}{\dot{y}}_E{\dot{y}}_A+M_{13}{\dot{x}}_E{\dot{y}}_E+M_{14}{\dot{x}}_E{}^2+M_{15}{\dot{y}}_A{}^2+M_{16}{\dot{y}}_E{}^2+kx+2{\lambda }_3(l-x)=-c\dot{x}$$

(8)

$$M_{21}+M_{22}\dot{x}+M_{23}{\dot{x}}_E+M_{24}\left({\dot{y}}_A-{\dot{y}}_E\right)+\frac{(m_2+m_4)}{4}{\ddot{x}}_E+2({\lambda }_1+{\lambda }_3)x_E+M_{25}\left({\dot{y}}_E-{\dot{y}}_H\right)=0$$

(9)

$$\left.\begin{array}{l}{M}_{31}+M_{32}{\dot{x}}_E{}^2+M_{33}{\dot{x}}_E{}^2+M_{34}{\dot{y}}_A{}^2+M_{35}{\dot{x}}_E{\dot{y}}_E+M_{36}{\dot{y}}_E{}^2+M_{37}{\dot{x}}_E\dot{x}+M_{38}{\dot{y}}_A\dot{x}\\ +M_{39}{\dot{y}}_E\dot{x}+M_{310}{\dot{x}}_E{\dot{y}}_A+M_{311}{\dot{y}}_E{\dot{y}}_A+M_{312}{\dot{x}}_E{\dot{y}}_H+M_{313}{\dot{y}}_E{\dot{y}}_H+M_{314}{\ddot{y}}_A\\ +M_{315}{\ddot{y}}_E+M_{316}{\dot{y}}_H{}^2+M_{317}{\ddot{x}}_E+M_{318}{\ddot{y}}_H=g(m_2+m_4+m_5+m_b)\end{array}\right\}$$

(10)

$$M_{41}+M_{42}{\dot{x}}_E{}^2+M_{43}({\dot{y}}_E{}^2+{\dot{y}}_H{}^2)+M_{44}{\ddot{y}}_H+M_{45}{\ddot{x}}_E+M_{46}{\ddot{y}}_E+\frac{m_5}{4}{\ddot{y}}_A+M_{47}{\dot{x}}_E{\dot{y}}_H+M_{48}{\dot{y}}_E{\dot{y}}_H+M_{49}{\dot{x}}_E{\dot{y}}_E=0$$

(11)

$$\left.\begin{array}{l}{M}_{51}+M_{52}{\dot{x}}_E{}^2+M_{53}{\dot{y}}_A{}^2+M_{54}{\dot{x}}_E{\dot{y}}_E+M_{55}{\dot{y}}_E{}^2+M_{56}{\dot{x}}_E\dot{x}+M_{57}{\dot{y}}_A\dot{x}+M_{58}{\dot{y}}_E\dot{x}\\ +M_{59}{\dot{x}}_E{\dot{y}}_A+M_{510}{\dot{y}}_E{\dot{y}}_A+M_{511}{\dot{x}}_E{\dot{y}}_A+M_{512}{\ddot{x}}_E+M_{513}{\ddot{y}}_A+M_{514}{\ddot{y}}_E+\frac{m_5}{4}{\ddot{y}}_H=0\end{array}\right\}$$

(12)

where

$$M_{13}=\frac{2m_4(l-x)(y_A-y_E)}{3C_2{}^2{x}_E{}^3},M_{14}=\frac{m_4(l-x){(y_A-y_E)}^2}{3C_2{}^2{x}_E{}^4}$$

$$M_{15}=\frac{m_4(l-x)}{3C_2{}^2{x}_E{}^2},M_{16}=\frac{m_4(l-x)}{3C_2{}^2{x}_E{}^2}$$

$$\begin{array}{cc}\hfill M_{21}=& -\frac{2{C_1}^2{m}_4{(l-x)}^2{(y_A-y_E)}^2}{3{C_2}^3{x_E}^3}+\frac{2C_1{C}_4{m}_4{(l-x)}^2(y_A-y_E)}{3{C_2}^3{x_E}^2}\hfill \\ \hfill & -\frac{C_1{C}_3{m}_4{(l-x)}^2}{3{C_2}^2}-\frac{C_{10}{l_2}^2{m}_2(y_E-y_H)}{3{C_6}^2{x_E}^2}\hfill \\ \hfill & -\frac{C_9{m}_4{(l-x)}^2(y_A-y_E)}{3{C_2}^2{x_E}^2}-\frac{2{C_5}^2{l_2}^2{m}_2{(y_E-y_H)}^2}{3{C_6}^3{x_E}^3}\hfill \\ \hfill & -\frac{C_5{C}_7{l_2}^2{m}_2}{3{C_6}^2}+\frac{2C_5{C}_8{l_2}^2{m}_2(y_E-y_H)}{3{C_6}^3{x_E}^2}\hfill \end{array}$$

$$\begin{gathered} M_{22}=\frac{2C_1{m}_4(l-x)(y_A-y_E)}{3C_2{}^2{x}_E{}^2} \\ {M}_{23}=\frac{2C_5{l}_2{}^2{m}_2(y_E-y_H)}{3C_6{}^2{x}_E{}^3}+\frac{2C_1{m}_4{(l-x)}^2(y_A-y_E)}{3C_2{}^2{x}_E{}^3} \end{gathered}$$

$$M_{24}=-\frac{C_1{m}_4{(l-x)}^2}{3C_2{}^2{x}_E{}^2},M_{25}=-\frac{C_5{l}_2{}^2{m}_2}{3C_6{}^2{x}_E{}^2}$$

$$M_{31}=\frac{g}{2}\left(m_2+m_4\right)+2{\lambda }_1(y_E-y_H)+2{\lambda }_3(y_E-y_A)+{\lambda }_4$$

$$M_{32}=(\frac{{(y_A-y_E)}^2}{C_2{x}_E{}^2}-1)\frac{2m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^2{x}_E{}^4}$$

$$M_{33}=(1-\frac{{(y_E-y_H)}^2}{C_6{x}_E{}^2})\frac{2l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^2{x}_E{}^4}$$

$$M_{34}=\frac{2m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^3{x}_E{}^4}$$

$$M_{35}=(\frac{2{(y_A-y_E)}^2}{C_2{x}_E{}^2}-1)\frac{2m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^3}+(\frac{2{(y_E-y_H)}^2}{C_6{x}_E{}^2}-1)\frac{2l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^3}$$

$$M_{36}=\frac{2m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^3{x}_E{}^4}-\frac{2l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^3{x}_E{}^4}$$

$$M_{37}=-\frac{2m_4(l-x)(y_A-y_E)}{3C_2{}^2{x}_E{}^3},M_{38}=\frac{2m_4(l-x)}{3C_2{}^2{x}_E{}^2},M_{39}=-\frac{2m_4(l-x)}{3C_2{}^2{x}_E{}^2}$$

$$M_{310}=(1-\frac{2{(y_A-y_E)}^2}{C_2{x}_E{}^2})\frac{2m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^3},M_{311}=-\frac{4m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^3{x}_E{}^4}$$

$$M_{312}=(1-\frac{2{(y_E-y_H)}^2}{C_6{x}_E{}^2})\frac{2l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^3},M_{313}=\frac{4l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^3{x}_E{}^4}$$

$$M_{314}=\frac{m_4}{4}-\frac{m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^2},M_{315}=\frac{m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^2}+\frac{l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^2}+\frac{m_2}{4}+\frac{m_4}{4}$$

$$\begin{gathered} M_{316}=-\frac{2l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^3{x}_E{}^4} \\ {M}_{317}=\frac{m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^2{x}_E{}^3}-\frac{l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^2{x}_E{}^3} \end{gathered}$$

$$M_{318}=\frac{m_2}{4}-\frac{l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^2}$$

$$\begin{gathered} M_{41}=\frac{g}{2}\left(m_2+m_5\right)-{\lambda }_2-2{\lambda }_1(y_E-y_H) \\ {M}_{42}=(\frac{{(y_E-y_H)}^2}{C_6{x}_E{}^2}-1)\frac{2l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^2{x}_E{}^4} \end{gathered}$$

$$\begin{gathered} M_{43}=\frac{2l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^3{x}_E{}^4},M_{44}=\frac{l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^2}+\frac{m_2}{4}+\frac{m_5}{4} \\ {M}_{45}=\frac{l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^2{x}_E{}^3},M_{46}=\frac{m_2}{4}-\frac{l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^2} \end{gathered}$$

$$M_{47}=(\frac{2{(y_E-y_H)}^2}{C_6{x}_E{}^2}-1)\frac{2l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^3},M_{48}=-\frac{4l_2{}^2{m}_2(y_E-y_H)}{3C_6{}^3{x}_E{}^4}$$

$$M_{49}=(1-\frac{2{(y_E-y_H)}^2}{C_6{x}_E{}^2})\frac{2l_2{}^2{m}_2}{3C_6{}^2{x}_E{}^3}$$

$$M_{51}=\frac{g}{2}\left(m_4+m_5\right)+{\lambda }_2-2{\lambda }_3(y_E-y_A)$$

$$\begin{gathered} M_{52}=(1-\frac{{(y_A-y_E)}^2}{C_2{x}_E{}^2})\frac{2m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^2{x}_E{}^4} \\ {M}_{53}=-\frac{2m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^3{x}_E{}^4} \end{gathered}$$

$$M_{54}=(1-\frac{2{(y_A-y_E)}^2}{C_2{x}_E{}^2})\frac{2m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^3},M_{55}=-\frac{2m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^3{x}_E{}^4}$$

$$M_{56}=\frac{2m_4(l-x)(y_A-y_E)}{3C_2{}^2{x}_E{}^3},M_{57}=-\frac{2m_4(l-x)}{3C_2{}^2{x}_E{}^2},M_{58}=\frac{2m_4(l-x)}{3C_2{}^2{x}_E{}^2}$$

$$M_{59}=\frac{4m_4{(l-x)}^2{(y_A-y_E)}^2}{3C_2{}^3{x}_E{}^5},M_{510}=\frac{4m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^3{x}_E{}^4}$$

$$M_{511}=-\frac{2m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^3},M_{512}=-\frac{m_4{(l-x)}^2(y_A-y_E)}{3C_2{}^2{x}_E{}^3}$$

$$M_{513}=\frac{m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^2}+\frac{m_4}{4}+\frac{m_5}{4},M_{514}=\frac{m_4}{4}-\frac{m_4{(l-x)}^2}{3C_2{}^2{x}_E{}^2}$$

$$C_1=\frac{{\dot{y}}_A-{\dot{y}}_E}{x_E}-\frac{\left(y_A-y_E\right){\dot{x}}_E}{x_E{}^2}{C}_2=1+\frac{{\left(y_A-y_E\right)}^2}{x_E{}^2}$$

$$C_3=\frac{2\left(y_A-y_E\right){\dot{x}}_E}{x_E{}^3}-\frac{{\dot{y}}_A-{\dot{y}}_E}{x_E{}^2}$$

$$C_4=\frac{2\left(y_A-y_E\right)\left({\dot{y}}_A-{\dot{y}}_E\right)}{x_E{}^2}-\frac{2{\left(y_A-y_E\right)}^2{\dot{x}}_E}{x_E{}^3}$$

$$C_5=\frac{{\dot{y}}_E-{\dot{y}}_H}{x_E}-\frac{\left(y_E-y_H\right){\dot{x}}_E}{x_E{}^2}$$

$$C_6=1+\frac{{\left(y_E-y_H\right)}^2}{x_E{}^2}$$

$$C_7=\frac{2\left(y_E-y_H\right){\dot{x}}_E}{x_E{}^3}-\frac{{\dot{y}}_E-{\dot{y}}_H}{x_E{}^2}$$

$$C_8=\frac{2\left(y_E-y_H\right)\left({\dot{y}}_E-{\dot{y}}_H\right)}{x_E{}^2}-\frac{2{\left(y_E-y_H\right)}^2{\dot{x}}_E}{x_E{}^3}$$

$$C_9=\frac{2\left(y_A-y_E\right){\dot{x}}_E{}^2}{x_E{}^3}-\frac{2{\dot{x}}_E\left({\dot{y}}_A-{\dot{y}}_E\right)}{x_E{}^2}-\frac{\left(y_A-y_E\right){\ddot{x}}_E}{x_E{}^2}+\frac{{\ddot{y}}_A-{\ddot{y}}_E}{x_E}$$

$$C_{10}=\frac{2\left(y_E-y_H\right){\dot{x}}_E{}^2}{x_E{}^3}-\frac{2{\dot{x}}_E\left({\dot{y}}_E-{\dot{y}}_H\right)}{x_E{}^2}-\frac{\left(y_E-y_H\right){\ddot{x}}_E}{x_E{}^2}+\frac{{\ddot{y}}_E-{\ddot{y}}_H}{x_E}$$

$$m_2={\rho }_a{l}_2$$

$$m_4={\rho }_{d}l$$

$$m_5={\rho }_f{l}_{ah}$$

$$l_{ah}=\sqrt{l^2-l_2{}^2}$$

The relevant symbols in the quarter suspension dynamics model are defined as shown in

Table 1. Parameter setting of the quarter suspension model.

Parameter	Symbol	Value	Unit
Body mass	$m_b$	50.000	$kg$
Length of damper	$l$	- *	$m$
Length of horizontal arm	$l_2$	- *	$m$
Damping coefficient	$c$	- *	$N\cdot s\cdot m^{-1}$
Stiffness coefficient	$k$	- *	$N\cdot m^{-1}$
Mass coefficient of damper	${\rho }_d$	3.566	$kg\cdot m^{-1}$
Mass coefficient of horizontal arm	${\rho }_a$	1.614	$kg\cdot m^{-1}$
Mass coefficient of frame	${\rho }_f$	2.400	$kg\cdot m^{-1}$

* These values are given by the optimization procedure.

.

Suspension system parameters include damper length ($l$), wishbone length ($l_2$), damping factor (c) and stiffness factor (k). The selection of the suspension system parameters affects the objectives of the mobile robot’s own weight, stability, the degree of fluctuation and the deformation of the shock absorber. These four objectives are described as follows.

1.

The total mass of the quarter suspension system is used as a measure of the degree of lightness, as shown in Equation (13).

$$J_1=Min(M_{sum})$$

(13)

where $M_{sum}={\displaystyle \sum m_i}$.

2.

The height of the centre of mass at steady state is used as a measure of suspension stability, as shown in Equation (14).

$$J_2=Min(Y_{cm})$$

(14)

where $Y_{cm}=\frac{{\displaystyle \sum m_i{y}_i}}{{\displaystyle \sum m_i}}$.

3.

The standard deviation of the suspension centre-of-mass height time response curve is used as a measure of the degree of system fluctuation, as shown in Equation (15).

$$J_3=Min(\sigma )$$

(15)

where $\sigma =\sqrt{\frac{{\displaystyle \sum _{t=1}^T{[{(Y_{cm})}_t-{\overline{Y}}_{cm}]}^2}}{T-1}}$.

4.

The ratio of the deformation of the shock absorber x to the wishbone $l_2$ is chosen as the deformation factor to measure the deformation of the shock absorber as shown in Equation (16).

$$J_4=Min(DF)$$

(16)

where $DF=\frac{x}{l_2}$.

The focus of the optimal design problem for the suspension system of an agricultural mobile robot is the selection of the optimal combination of parameters for the quarter suspension model damper length ($l$), wishbone length ($l_2$), damping factor (c) and stiffness factor (k) to achieve the optimal objective function, i.e., the optimal suspension own weight, stability, degree of system fluctuation and damper deformation.

2.2. Suspension Structure Optimization Based on Non-Dominated Sorting Genetic Algorithm III (NSGA-III) and Technique for Order Preference by Similarity to Ideal Solution (TOPSIS)

In single-objective optimization problems, since only one objective function needs to be considered optimal, the optimal solution can be found by some common mathematical methods. In multi-objective optimization problems, however, there may be certain constraints between the various objective functions. When one of the objectives reaches optimality, the performance of the other objectives may be unacceptable. Therefore, solving the multi-objective problem for the optimal solution is very difficult. Further, in the solution of multi-objective optimization problems, the non-dominated solutions set, the Pareto set, is usually used to represent the acceptable better solution. The solutions of the Pareto set are non-dominated by any other solution, and are considered as equally optimal solutions. With the TOPSIS method, it is possible to select the optimal solution needed from the non-dominated solution set.

2.2.1. Non-Dominated Sorting Genetic Algorithm III (NSGA-III)

The non-dominated genetic sorting algorithm is a multi-objective evolutionary algorithm based on a genetic algorithm that simulates the evolutionary behaviour of natural populations in terms of “survival of the fittest” and the reproduction of offspring.

In natural selection, individuals in a population that are well adapted to nature are more likely to survive and thus reproduce their offspring, passing on the genes for good traits from generation to generation. The process of reproduction produces offspring with different traits from their parents due to crossover and variation, ensuring a diversity of traits in the population. Through multiple generations of evolution, the individuals that eventually survive in the population are those that are well adapted to nature.

Compared to NSGA, NSGA-II reduces the computational complexity of nondominated sorting, introduces and elitism strategy which could help prevent the loss of good solutions, and doesnot require specifying the sharing parameter. With the above improvements, the iterative convergence speed of NSGA-II is improved, and the computational complexity is reduced from $O(M{N}^3)$ to $O(M{N}^2)$. NSGA-II, when faced with high-dimensional multi-objective optimization problems with more than three objectives, suffers from the shortcomings of convergence. Compared to NSGA-II, NSGA-III, based on the reference point selection mechanism, is effective in reducing computational complexity and improving convergence for high-dimensional multi-objective optimization problems with a large number of objectives. NSGA-III improves on the population update selection mechanism of NSGA-II by providing and adaptively updating a number of well-distributed reference points to help maintain diversity among population members. With these improvements, the computational complexity of NSGA-III is the greater of $O(M{N}^2)$ or $O(N^2{\log }^{M-2}N)$.

2.2.2. Technique for Order Preference by Similarity to Ideal Solution (TOPSIS)

The TOPSIS method is a ranking method that approximates an ideal solution and makes full use of the information from raw data to accurately reflect the gaps between individual evaluation objects. It is a method of ranking a limited number of evaluation objects according to their proximity to an idealised target, and is an evaluation of the relative merits of the available objects. Therefore, TOPSIS is a common and effective method for multi-objective decision analysis.

In this study, since there are a finite number of evaluation objects in the non-dominated solution set, and each object has four known evaluation indicators, the TOPSIS method of scoring is as follows.

Assume that there are n solutions in the non-dominated solution set, i.e., there are n evaluation objects and four evaluation indicators, consisting of a standardised matrix.

$$Z=\left[\begin{array}{cccc}{z}_{11}& z_{12}& z_{13}& z_{14}\\ z_{21}& z_{22}& z_{23}& z_{24}\\ ⋮& ⋮& ⋮& ⋮\\ z_{n1}& z_{n2}& z_{n3}& z_{n4}\end{array}\right]$$

(17)

Define the minimum value

$Z^{-}=(Z_1{}^{-},Z_2{}^{-},Z_3{}^{-},Z_4{}^{-})=(\min \{z_{11},z_{21},\cdots ,z_{n1}\},\min \{z_{12},z_{22},\cdots ,z_{n2}\},\min \{z_{13},z_{23},\cdots ,z_{n3}\},\min \{z_{14},z_{24},\cdots ,z_{n4}\})$

Define the distance of the $i(i=1,2,\cdots ,n)$ evaluation object from the minimum value $D_i{}^{-}=\sqrt{{\displaystyle \sum _{j=1}^4{(Z_j{}^{-}-z_{ij})}^2}}$

The unnormalised score of the $i(i=1,2,\cdots ,n)$ evaluation object can be calculated $S_i=\frac{1}{D_i{}^{-}}$

Normalised scores for the $i(i=1,2,\cdots ,n)$ evaluation objects ${\overline{S}}_i=\frac{S_i}{{\displaystyle \sum _{i=1}^n{S}_i}}$.

The individual with the highest score based on the normalised score is selected as the optimal solution.

2.2.3. Parameters Optimization of the Suspension System Based on NSGA-III and TOPSIS

This section describes the use of NSGA-III to optimize the structural parameters of a robotic suspension system and to select the optimal solution in a non-dominated solutions set using the TOPSIS method. The objective of the optimized robotic suspension system problem is to find optimal values for the system parameters (damper length, cross arm length, stiffness and damping coefficient), aiming to minimise the four evaluation metrics developed (total suspension mass, suspension centre of mass height, suspension fluctuation coefficient and suspension compression coefficient) to obtain the highest score. Algorithm 1 describes the process of optimizing the optimum parameter values for the quarter passive robot suspension system.

Algorithm 1: Parameters Optimization of the Suspension System Based on NSGA-III and TOPSIS

1:  Input the ranges of ($l$), ($l_2$), (k) and (c)

2:  Initialize (N) search agents of NSGA-III with four values of $l^{(1:N)}$, $l_2{}^{(1:N)}$, $k^{(1:N)}$ and $c^{(1:N)}$

3:  Calculate the indicators ($M_{sum}{}^{(1:N)}$, $Y_{CM}{}^{(1:N)}$, ${\sigma }^{(1:N)}$, $D{F}^{(1:N)}$) of search agents using the multibody dynamic modal.

4:  Determine the nondominated solution set ($l^{F_1}$, $l_2{}^{F_1}$, $k^{F_1}$ and $c^{F_1}$) that have the lowest indicators ($M_{sum}{}^{(1:N)}$, $Y_{CM}{}^{(1:N)}$, ${\sigma }^{(1:N)}$, $D{F}^{(1:N)}$).

5:  Update the search agents values ($l^{(1:N)}$, $l_2{}^{(1:N)}$, $k^{(1:N)}$ and $c^{(1:N)}$) based on NSGA-III genetic operations.

6:  Repeat from step 3 for T iterations.

7:  Output the nondominated solution set ($l^{F_1}$, $l_2{}^{F_1}$, $k^{F_1}$ and $c^{F_1}$).

8:  Calculate the scores of the nondominated solution set ($l^{F_1}$, $l_2{}^{F_1}$, $k^{F_1}$ and $c^{F_1}$) by the indicators ($M_{sum}{}^{F_1}$, $Y_{CM}{}^{F_1}$, ${\sigma }^{F_1}$, $D{F}^{F_1}$) based on TOPSIS.

9:  Select the highest score solution as the best global solution ($l^{gbest}$, $l_2{}^{gbest}$, $k^{gbest}$ and $c^{gbest}$).

3. Results and Discussion

To investigate the advantages of multi-objective evolutionary algorithms in multi-objective optimization problems, we compared the performance of NSGA-III with some multi-objective evolutionary algorithms. At the same time, we built a quarter suspension model in Adams software to perform a dynamics analysis of the parameter optimization results and compare the time response of the parameter combinations in the optimal solution set for the suspension dynamics analysis.

3.1. Comparison of Multi-Objective Evolutionary Algorithms’ (MOEA) Results

The non-dominated ranking genetic algorithm II (NSGA-II), the covariance matrix adaptive evolution strategy (CMA-ES), the third generation generalised difference algorithm (GDE3), the decomposition-based multi-objective evolution algorithm (MOEAD), the multi-objective particle swarm optimization algorithm (OMOPSO), the speed-constrained multi-objective particle swarm optimization algorithm (SMOPSO) and the strength Pareto Evolutionary Algorithm (SPEA2) as representatives of multi-objective evolutionary algorithms were compared with NSGA-III. The parameter settings for the experimental tests are shown in

Table 2. Parameter setting of the optimization test.

Parameter	Lower Bound	Upper Bound	Unit
$l$	0.150	0.450	$m$
$l_2$	0.035	0.300	$m$
$c$	1.000	1000.000	$N\cdot s\cdot m^{-1}$
$k$	0.100	30,000.000	$N\cdot m^{-1}$

.

The differential equation model based on Equations (8)–(12) was optimally solved within the parameter intervals shown in Table 2 using the Wolfram Language NDSolve function and the multi-objective evolutionary algorithm described above. NDSolve function is a general-purpose numerical solver for differential equations. It can solve many ordinary differential equations (ODEs) and some partial differential equations (PDEs). The set of Pareto solutions obtained from each evolutionary algorithm solution was selected using the TOPSIS method to find the optimal solution. The optimal solutions obtained by the eight evolutionary algorithms were integrated and scored again using the TOPSIS method.

The experimental results of the optimization of the multi-objective dynamics model using NSGA-III on a quarter suspension model compared to other multi-objective evolutionary algorithms are shown in

Table 3. Optimum score of the NSGA-III in comparison with other evolutionary algorithms.

Algorithm	Best	Worst	Average	St. Dev
CMA-ES	0.02879	0.01780	0.02394	0.00424
GDE3	0.03361	0.01547	0.02198	0.00648
MOEAD	0.03835	0.01957	0.02782	0.00655
OMOPSO	0.02954	0.01375	0.02094	0.00537
SMPSO	0.03858	0.01252	0.02681	0.00896
SPEA2	0.03278	0.01483	0.02250	0.00753
NSGA-II	0.02954	0.01387	0.02141	0.00554
NSGA-III	0.04048	0.03017	0.03460	0.00381

,

Table 4. Based on TOPSIS method, the optimal solution sets integrated by eight evolutionary algorithms were comprehensively scored, and the top five solutions are listed.

Rank	Algorithm	Score	l (m)	l2 (m)	k (Nm−1)	c (Nsm−1)
1	NSGA-III	0.04048	0.21790	0.08181	19,201.10593	385.76418
2	SMPSO	0.03858	0.21371	0.08497	16,474.23182	188.98510
3	MOEAD	0.03835	0.22376	0.07908	21,444.41076	300.06606
4	NSGA-III	0.03749	0.21944	0.09227	14,775.16239	237.24490
5	GDE3	0.03361	0.23349	0.06862	24,567.76130	404.21663

and

Table 5. The top five scoring solutions correspond to the four design objective values based on Equations (13)–(16), respectively.

Rank	Algorithm	Msum (kg)	Ycm (m)	σ (×10−5 m)	DF
1	NSGA-III	1.39373	0.33824	4.33187	0.31139
2	SMPSO	1.36980	0.33203	6.27938	0.34897
3	MOEAD	1.42793	0.34342	3.33967	0.29398
4	NSGA-III	1.40929	0.33234	6.18176	0.32603
5	GDE3	1.47897	0.35391	2.31888	0.28621

. Each algorithm was run five times and standard deviations analysed to determine the stability of the algorithm. The experimental results show that the NSGA-III achieved optimal evaluation metrics for best, worst and average results with the lowest standard deviation. Meanwhile, among the 40 optimal solution sets obtained by the eight multi-objective evolutionary algorithms, the optimal solution obtained by the NSGA-III achieved the highest score in the TOPSIS composite score. This indicates that the NSGA-III outperforms other comparative algorithms in terms of solution quality and robustness in solving multi-objective suspension optimization design problems. A comparison of the computational time of NSGA-III with other multi-objective evolutionary algorithms is presented in

Table 6. Computational time of the NSGA-III in comparison with other evolutionary algorithms.

Algorithm	Shortest (s)	Longest (s)	Average (s)	St. Dev
CMA-ES	37.86	43.88	40.76	2.04
GDE3	37.58	54.15	47.44	5.54
MOEAD	174.59	216.50	193.76	14.05
OMOPSO	40.62	52.04	44.05	4.35
SMPSO	45.54	62.82	51.32	6.49
SPEA2	41.24	54.02	46.28	4.22
NSGA-II	42.95	50.79	46.85	2.51
NSGA-III	173.21	195.65	183.31	8.78

. All optimizations were conducted on a laptop with Intel i7-6700HQ (2.60 GHz) and 16 GB RAM. Both NSGA-III and MOEAD took longer than the others in terms of computation time. However, there was no real-time requirement for the task of suspension design. Therefore, the average computation time of 183.31 s for NSGA-III is acceptable, which does not affect the excellent performance of NSGA-III in accomplishing the design problem targeted in our study.

3.2. Simulation Experiments

Further comparative validation of the above optimization results was carried out. A parametric quarter suspension model was built in Adams software and the dynamics of the optimized parameter combinations were analysed to compare results.

The Adams model used in the dynamic simulation is shown in Figure 3, and its relevant parameters are shown in

Table 7. Main parameters of the Adams parametric model for dynamics simulation.

Parameter	Symbol	Value	Unit
Mass of body	$m_b{}^{Sim}$	50.000	$kg$
Length of damper	$l^{Sim}$	- *	$m$
Length of horizontal arm	$l_2{}^{Sim}$	- *	$m$
Damping coefficient	$c^{Sim}$	- *	$N\cdot s/m$
Stiffness coefficient	$k^{Sim}$	- *	$N/m$
Density of damper	${\rho }_d{}^{Sim}$	11,350.931	$kg/m^3$
Density of horizontal arm	${\rho }_a{}^{Sim}$	29,494.779	$kg/m^3$
Density of frame	${\rho }_f{}^{Sim}$	5328.001	$kg/m^3$

* These values relate to the optimization results.

. In the Adams model, we used RigidBody: Link to represent frame and horizontal arm in the suspension system; RigidBody: Cylinder to represent damper, and RigidBody: Box to represent the body, under-spring mass (e.g., wheel) and the ground. Artificial masses of all entities were used to ensure that the model parameters in the Adams simulation were the same as the PDE model parameters, involving the parameterized frame, horizontal arm and damper corresponding RigidBody: Link density calculated by linear density shown in Table 1.

Table 8. Parameters of the globally optimal solution (obtained by NSGA-III) in the partial differential equation model compared with those in the Adams model.

Parameter	PDE Value (kg)	Adams Value (kg)
Mass of body	50.00000	50.00000
Mass of damper	0.77703	0.77703
Mass of horizontal arm	0.13204	0.13204
Mass of frame	0.55860	0.55860

compares the parameters of the globally optimal solution (obtained from NSGA-III) in the PDE model with those in the Adams model. Since the parameters of the Adams model were calculated based on the parameters of the PDE model, the quality parameters of the two models were exactly the same. In Adams View, add fixed, revolute and translational joints to the model at the corresponding position.

A step signal was applied at the wheel-ground contact position, the amplitude of which was the sum of the suspension gravity and the design load of the individual wheels. Three parameters: suspension centre of mass acceleration, suspension centre of mass height and damper spring deflection coefficient, were evaluated for comparison over a period of 1.5 s. The parameter response curves over time are shown in Figure 4, Figure 5 and Figure 6.

Analysis of Figure 4 shows that the optimal combination of suspension parameters derived from NSGA-III can quickly reduce the body acceleration while keeping the acceleration within a small fluctuation range, which has a positive impact on maintaining smoothness in the work of agricultural mobile robots. By analysing Figure 5, it can be seen that the NSGA-III optimal combination of suspension parameters allows the suspension to reach a steady state in the shortest possible stabilisation time, while keeping the suspension’s centre of mass variation range to a minimum, which is beneficial to the mobile robot’s ability to effectively maintain and quickly recover a steady state when encountering bumps. Analysis of Figure 6 shows that the optimal combination of suspension parameters obtained by NSGA-III enables the spring deflection coefficient of the suspension to be kept at a low value. The deflection coefficient was chosen as a measure of spring deflection rather than a direct measure of deflection because an increase in the length of the cross-arm causes an increase in spring deflection when other parameters are held constant. The use of a deflection factor avoids this problem affecting the optimization process. The results of the combined dynamics simulations show that the NSGA-III combined with the TOPSIS method is indeed the best solution in the Pareto solution set of the multi-objective evolutionary algorithm for the multi-objective suspension optimization problem.

Combining Table 4 and Table 5 and Figure 4, Figure 5 and Figure 6, it is obvious that no solution is optimal for all evaluation criteria. For example, when we wish to reduce the deformation factor (DF) of the damper, we can do so by increasing the length of the horizontal arm (l₂). However, at the same time, increasing the horizontal length (l₂) leads to an increase in the total mass of the suspension (M_sum). The conflict between the objectives makes it difficult to select an optimal solution from the non-dominated solution set, which is why we introduced the TOPSIS method for optimal solution selection. The TOPSIS method measured the advantages and disadvantages of different non-dominated solutions in a scientific and rational way and scored them. Table 4 shows that one of the solutions obtained by NSGA-III after analysis using the TOPSIS method received the highest score. The four objective values corresponding to the top five scoring solutions are shown in Table 5. We identified that both of the four objectives of the highest scoring NSGA-III solution were not the smallest, i.e., optimal. The smallest total mass (M_sum) and the lowest centre-of-mass height (Y_cm) appeared in the solution obtained by SMPSO and the smallest degree of centre-of-mass fluctuation(σ) and the smallest degree of damper deformation (DF) appeared in the solution obtained by GDE3. Meanwhile, the time response curves for the NSGA-III derived solution in Figure 4, Figure 5 and Figure 6, respectively, were not optimal. In the centre-of-mass acceleration time response curves and the damper deformation factor time response curves, the GDE3 solution exhibited better attenuation than NSGA-III, while both the peak acceleration and deformation factor were smaller than the NSGA-III solution. In the centre-of-mass height time response curve, another solution derived by NSGA-III (ranked 4th) was able to obtain the lowest centre-of-mass height, which meant that it has the best stability.

The above results do not seem to indicate that the solution found by NSGA-III is optimal, but when the four objectives were considered collectively, the solution found by NSGA-III showed superiority over the other solutions. We calculated the difference between each objective and their minimum of the sample set according to the TOPSIS method. This meant that any solution will only receive a higher score if all four objectives are as close as possible to their minimum of the sample set. The SMPSO derived solutions had the largest degree of centre-of-mass fluctuation (σ) and damper deformation factor (DF), the GDE3 derived solutions had the second largest total mass (M_sum) and the largest height of mass (Y_cm); therefore, their scores were weakened. Under the comprehensive scoring of the TOPSIS method, the solution obtained by NSGA-III had the highest score, i.e., the optimal solution.

4. Conclusions

We proposed a multi-objective optimization design method for the suspension of agricultural robots that can balance different performance indicators and obtain the optimal parameter combination of the suspension. We utilized the Lagrangian equation to establish the partial differential equation(PDE) model of the agricultural robot double wishbone suspension, including structural parameters and performance parameters, and utilized the multi-objective evolutionary algorithm NSGA-III and TOPSIS method to carry out the agricultural robot double wishbone optimum design of the suspension. We established four evaluation indicators to evaluate the performance of the suspension, including the total mass of the suspension system, the barycenter’s height of the suspension in a stable state, the fluctuation degree of the suspension under the step response, and the deformation degree of the damper. We selected eight typical multi-objective evolutionary algorithms to solve the multi-objective suspension optimization design problem, took advantage of the TOPSIS method to score the non-dominated solution set, and selected the solution with the highest score as the optimal solution. The results show that the optimal solution obtained by NSGA-III and TOPSIS method has advantages in a comprehensive performance. Therefore, we come to the conclusion that NSGA-III combined with the TOPSIS method can effectively obtain a high-quality design of an agricultural mobile robot’s suspension system.