An improved multi-objective artificial physics optimization algorithm based on R2 indicator and its application

Abstract

Artificial Physics Optimization (APO) algorithms in solving constrained multi-objective problems often encounter challenges such as uneven population distribution and imbalances between global and local search capabilities. To address these issues, we propose the Constrained Rank Multi-Objective Artificial Physical Optimization based on the R2 indicator (R2-ICRMOAPO) algorithm. This algorithm integrates non-dominated sorting with the R2 indicator and updates the external storage set using the contribution value derived from the R2 indicator formula, ensuring both set distribution and convergence. It also dynamically adjusts inertia weights and gravitational factors to enhance its global and local search capabilities. To evaluate the performance of the R2-ICRMOAPO algorithm, we compared it with four other multi-objective optimization algorithms using standard test functions. The results indicate that it demonstrates superior distribution and optimization performance. Furthermore, we applied it to optimize the parameters of a hydro-pneumatic suspension system. The experimental results show that it can reduce the root mean square value of car body vertical acceleration by approximately 21.4% and the root mean square value of dynamic tire load by about 19.6%. This reduction effectively enhances vehicle smoothness within a reasonable range. Consequently, these results confirm the feasibility of the R2-ICRMOAPO algorithm to solve practical problems.

Keywords

Artificial physical optimization multi-objective optimization R2 indicator hydro-pneumatic suspension

Introduction

In recent years, many problems in engineering fields are finally transformed to solve the optimal value of constrained multi-objective optimization problems (CMOPs). With the emergence of a multi-objective optimization problem (MOP), scholars began to propose traditional optimization methods. They typically adopt certain techniques to convert the MOP into a single-objective optimization problem (SOP) and uses SOP methods to solve MOP. Although existing SOP methods are well-established and can assist in solving MOPs, they necessitate a comprehensive understanding of the relevant background knowledge. Moreover, each optimization process yields only one Pareto-optimal solution, which necessitates multiple iterations to fully explore the solution space. Swarm Intelligence (SI) algorithms take inspiration from natural group behaviors, simulating collaboration and information exchange among individuals to identify the best solution to a problem. Consequently, they exhibit significant potential and advantages in addressing complex real-world challenges. A variety of innovative optimization approaches have emerged, including Ant Colony Algorithm (ACO),¹ Particle Swarm Optimization (PSO),² Firefly Algorithm,³ Grey Wolves Algorithm,⁴ Differential Evolution Algorithm,⁵ among others. ACO and PSO are regarded as the most classical SI algorithms, known for their robustness and powerful search capabilities. However, ACO can be computationally demanding and time-consuming to solve, while PSO is highly sensitive to parameter settings and susceptible to getting trapped in local optima.⁶ Artificial Physics Optimization (APO) is an emerging SI algorithm proposed in 2010. It can prevent the convergence to local optima that PSO algorithms often encounter during global optimization. APO has been effectively applied to solving SOPs, demonstrating strong global search capabilities, stability, rapid convergence, and good robustness. The APO algorithm is analogous to the PSO algorithm, which has been widely adopted and extensively applied in solving MOPs. Consequently, APO can be effectively applied to solve MOPs.⁷ Wang et al.⁷ analyzed the MOP from the perspective of virtual forces and employed fast sorting to rank non-dominated solution sets. Sun et al.⁸ introduced concepts from unconstrained multi-objective optimization into the constrained domain, proposing an improved APO algorithm that incorporates an ordinal value and congestion degree based on the constraint violation degree judgment criterion. However, the comparison process between new individuals and archived solution sets often leads to the loss of high-quality solutions and fails to identify those that are more widely and uniformly distributed. Consequently, further improvements to the algorithm are necessary.

Currently, researchers have proposed various strategies to design improved algorithms for MOPs, which can be categorized into three types: Pareto-dominated optimization algorithms, decomposition-based optimization algorithms, and indicator-based optimization algorithms. Pareto-dominated optimization algorithms may face selection difficulties at the initial iteration as the number of non-dominated solutions increases with the number of iterations, such as Non-dominated Sorting Genetic Algorithm-II.⁹ Decomposition-based optimization algorithms obtain a single objective function for single-objective optimization by aggregating the sub-problem functions together using weight aggregation, Chebyshev and cross-boundary methods, such as Multi-objective Evolutionary Algorithm Based on Decomposition.¹⁰ The core of indicator-based optimization algorithms is to design or select the appropriate performance evaluation indicator to improve algorithm performance, such as an IGD∼+ Performance Indicator Based Particle Swarm Optimizer for Multi-objective Optimization.¹¹During the process of improving the design of algorithms for MOPs, it is common to start with the use of some metrics capable of evaluating algorithm performance. Common indicators include Hypervolume (HV)¹² and R2.¹³ Although HV metrics are widely used in the current research phase, the time consumed in the solution selection grows exponentially with the number of objectives. In contrast, the R2 metric is less computationally intensive and is characterized by weak monotonicity. Therefore, the R2 metric can serve as a selection mechanism for updating the algorithm. Zitzler et al.¹⁴ defined the R2 metric in 2008 and gave the formula for the standard R2 metric. Wei et al.¹⁵ proposed a hybrid multi-objective PSO algorithm based on R2 metric by introducing R2 metric, S-shaped function mapping method, and simulated binary crossover operator to improve the algorithm. Liu et al.¹⁶ used a selection strategy guided by the R2 metric and reference vectors to propose a two-stage evolutionary algorithm based on the R2 metric. In this paper, a constrained rank multi-objective artificial physical optimization based on the R2 indicator (R2-ICRMOAPO) algorithm is proposed to address the shortcomings of the ICRMOAPO algorithm. This algorithm leverages the R2 indicator for a comprehensive evaluation of the solution set’s characteristics. It integrates non-dominated sorting with the R2 metric to update and maintain the external storage set. This process is guided by the contribution value of the R2 indicator, ensuring the retention of a uniformly distributed set of elite solutions. The algorithm balances global and local search capabilities by dynamically adjusting inertia weights and gravity factors based on the iteration count. It has demonstrated superior convergence and distribution performance compared to advanced multi-objective algorithms on standard test functions, proving its effectiveness in dealing with complex CMOPs. Moreover, its practicality is confirmed through successful application to a real-world engineering problem.

There are six parts to this paper. Section 2 provides background knowledge, including the concepts of CMOP and APO algorithms. Section 3 introduces the improved algorithm: R2-ICRMOAPO. Section 4 details the experimental setup and presents numerical results on various benchmark problems. Section 5 discusses the application of R2-ICRMOAPO in hydro-pneumatic suspension. The final section draws conclusions and suggests directions for future studies.

Background

The objective of this paper is to address the CMOP, and the solution we employ is an improved APO algorithm. We begin by introducing the CMOP and the APO algorithm, followed by the presentation of our enhanced algorithm.

Constrained multi-objective problems

Most engineering design problems require that multiple indicators simultaneously achieve optimal values within a specified range, and such problems are termed CMOPs.

The mathematical model of CMOPs can be expressed as follows:

\begin{matrix} min F (x) = [f_{1} (x), f_{2} (x), \dots, f_{s} (x)] \\ s . t . p_{i} (x) < 0, i = 1, 2, . . ., m \\ q_{j} (x) = 0, j = 1, 2, . . ., n \\ x = (x_{1}, x_{2}, \dots, x_{k}) \\ F (x) \in R^{s}, x \in R^{k} \end{matrix}

(1)

Where $x$ is the decision variable, $F (x)$ is the objective function, and $p (x)$ and $q (x)$ are inequality and equation constraints, respectively, which together define the boundaries of the feasible region.

For CMOPs, the infeasibility of each solution is:

ϕ (x_{i}) = \sum_{j = 1}^{k} {[min {0, g_{i} (x_{i})}]}^{2} + \sum_{p = 1}^{l} {[h_{p} (x_{i})]}^{2}

(2)

An infeasible domain consists of all points that do not satisfy the constraints. Conversely, when the infeasibility measure for a solution $x_{i}$ is 0, $x_{i}$ is considered a feasible solution.

Artificial physical optimization algorithm

APO algorithm is an intelligent optimization algorithm based on population search, its essence is the simulation of Newton’s second law. Each individual is treated as a particle in physics. By calculating the individual’s mass, the rules of the interaction force between individuals, and acceleration, the algorithm updates each individual’s position to solve optimization problems. It is characterized by a simple principle, species diversity, and high search efficiency. The algorithm can be described as follows:

Suppose the population size is $M$ . The mass function of $x_{i}$ can be expressed as:

m_{i} = e^{\frac{f (x_{best}) - f (x_{i})}{f (x_{worst}) - f (x_{best})}}

(3)

Where $f (x)$ is the value of the objective function for the decision variable.

The virtual force formula is derived from the force of gravity. Consequently, the virtual force that individual $i$ is subjected to in the $k - th$ dimension is given by:

\begin{matrix} F_{ij, k} = \\ {\begin{matrix} G m_{i} m_{j} (x_{j, k} - x_{i, k}), f (x_{i}) > f (x_{j}) \\ - G m_{i} m_{j} (x_{j, k} - x_{i, k}), f (x_{i}) \leq f (x_{j}) \end{matrix} \forall i \neq j, i \neq best \end{matrix}

(4)

Where $x_{j, k}$ represent the position of the individual $j$ in the $k - th$ dimension, and $G$ denote the gravitational factor.

The formula for calculating the resultant force acting on $x_{i}$ in the $k - th$ dimension is given by:

F_{i, k} = \sum_{\binom{j = 1}{i \neq j}}^{n} F_{ij, k}, \forall i \neq best

(5)

The velocity and position update of individual $x_{j}$ is given by the following equations:

v_{j, k} (t + 1) = w v_{j, k} (t) + \frac{λ F_{j, k}}{m_{j}}, v_{i} \in [v_{min}, v_{max}], λ \in [0, 1]

(6)

x_{j, k} (t + 1) = x_{j, k} (t) + v_{j, k} (t + 1), x_{i} \in [x_{min}, x_{max}]

(7)

Where $v_{j, k}$ represent the velocity of individual $j$ in the $k - th$ dimension, $w$ be the inertia weight, and $λ$ be a random variable that follows a certain distribution.

The flow of the APO algorithm is as follows:

Algorithm 1 APO
Input: population size $M$ and archive size Q.
Output: Q
1: initialize parameters 2: randomized generation of individual initial velocities and positions
3: for $i = 1$ to $M$ 4: calculate the objective function
5: end for
6: identify non-dominant individuals to be stored in the archive
7: while iter< $Ma x_{iterations}$ do
8: calculate the individual forces $F_{i, k}$ according to (5)
9: calculate the velocity and position of individuals according to (6) and (7)
10: update the archive
11: end while
12: returnQ

Improved algorithm: R2-ICRMOAPO algorithm

With the widespread application of CMOPs, corresponding optimization algorithms have emerged. The selection of non-dominated solutions is crucial to the algorithm’s performance on MOPs. This paper introduces the R2-ICRMOAPO algorithm in response to the deficiencies of APO algorithms, which include uneven population distribution and an imbalance between global and local search capabilities when solving MOPs.

The R2 indicator updates the external storage set

The R2 indicator is based on a utility function, which can be used both to estimate the relative mass of individuals and to distinguish the strengths and weaknesses of candidate solution schemes. In this context, the weighted Chebyshev function is used as the utility function for the R2 indicator. With the chosen utility function and a specific reference point, the quality of a set of individuals relative to the reference point can be evaluated, and the candidate solution with greater utility is selected. The R2 indicator not only comprehensively evaluates algorithms but also obtains a uniformly distributed solution set through rapid calculation. Therefore, it has been applied to solve multi-objective problems (MOPs).

The utility function expression is as follows:

\begin{matrix} R 2 (A, W, q^{*}) = \frac{1}{| W |} \sum_{μ \in W} min_{x \in A} {max_{i = 1, 2, \dots, m} μ_{i} | f_{i} (x) - q^{*} |}, \\ μ = (μ_{1}, μ_{2}, . . ., μ_{m}) \end{matrix}

(8)

Where $A$ is a group of solutions, $W$ is a set of weight vectors uniformly distributed in the target space, and $q^{*}$ is a reference point.

When different vectors of weights determine the functions within the utility function, The individual’s contribution to the R2 indicator (CR2) is:

CR 2 (x, A, W, q^{*}) = R 2 (A, W, q^{*}) - R 2 (A \ {x}, W, q^{*})

(9)

In the population, individuals can be fully sorted based on their CR2 values. A larger CR2 indicates that the solution set is better distributed, a higher quality of the solution, and a higher probability that it will be retained in the next generation. The size of non-dominated solutions increases with the size of iterations of the algorithm and may exceed the preset solution set size. If the external storage space is more than a predefined size, the non-dominant solution of the external storage set is evaluated by using (9), and individuals with smaller CR2 values are deleted.

The time complexity of R2

Assume that the size of the candidate set is M and the number of targets is n.

The R2 time complexity is the R2 metric for calculating each individual in the candidate set, that is $O (| M | | V | (\log | M | + n))$ .

Where |V| is the number of Chebyshev’s weights.¹⁷

Key elements of the algorithm

Assume that $P (t)$ is the $t - th$ generation population.

The ordinal value of an individual $i$ at generation $t$ is defined as:

r_{i} (t) = 1 + n_{i} (t)

(10)

Where $n_{i} (t)$ denotes the total number of dominant individuals $i$ in $P (t)$ .

When the ordinal values are equal, it is necessary to introduce a neighborhood radius. To prevent the neighborhood radius from being too small in the algorithm, it can be defined as follows:

r = \frac{{(\sum_{i = 1}^{n} max {(| E_{i} |)}^{2})}^{\frac{1}{2}}}{archivesize - 1}

(11)

Where $| E_{i} |$ is the Euclidean distance between individuals on dimension $i$ at initialization.

To take into account the dominant relationship between individuals and the radius of the neighborhood in which the individuals are located the degree of crowding. The mass function is⁸:

m_{i} (t) = e^{\frac{- ran k_{i} (t)}{M}}

(12)

The magnitude of the inertia weights $w$ is related to the current position affected by the individual’s velocity. A smaller value of $w$ indicates that the algorithm is more capable of localized search. The selection of $w$ will affect the performance and efficiency of the algorithm, which is crucial for achieving a balance between population development and exploration. Generally, large inertia weights are needed at the beginning of the iteration to maintain the global search ability and ensure convergence to the global optimum more quickly and accurately. At the end of the iteration, a small inertia weight is chosen to ensure accurate convergence to the optimal solution. Therefore, in the R2-ICRMOAPO algorithm, the inertia weight is decreased linearly, which allows the inertia weight to decrease dynamically as the iteration progresses:

w (iter) = w_{i} - \frac{(w_{i} - w_{f}) (iter - 1)}{Ma x_{iterations}}

(13)

Where $w_{i}$ and $w_{f}$ indicate the initial and final values of the inertia weights $w$ , and $iter$ indicates the size of iterations. In this algorithm, $w_{i} = 0.9$ and $w_{f} = 0.4$ .

The gravity coefficient significantly influences the performance of the algorithm. To some extent, it determines the number of individuals who make convergent motions. In the standard APO algorithm, the gravitational factor is typically set to a value of ten. The R2-ICRMOAPO algorithm improves this by making the gravitational factor change dynamically with the evolution of the population. Equation (14) is the formula for calculating the gravitational factor in the R2-ICRMOAPO algorithm.

G (iter) = G_{i} - \frac{(G_{i} - G_{f}) (iter - 1)}{Ma x_{iteration}}

(14)

Where $G_{i}$ and $G_{f}$ are the initial and final values of the gravitational factors $G$ , respectively. Given the large initial value of the inertia weight $w$ in this paper, it is necessary to set a correspondingly large initial value for the gravity factor $G$ to balance the global and local search capabilities of the algorithm. Thus, $G_{i} = 70$ and the final value of $G$ from iteration to the end is $G_{f} = 1$ .

Algorithm flow

The selection of non-dominated individuals significantly influences the convergence and distribution of the algorithm. As iterations proceed, the number of non-dominated solutions may increase and potentially exceed the pre-set size of the solution set. Consequently, updating and maintaining the external archive is crucial in the solution process. This paper introduces an approach that combines non-dominated ranking with the R2 indicator. Initially, the Pareto dominance relationship is used to identify a set of non-dominated solution sets, and the R2 indicator contribution is calculated for each. Subsequently, the algorithm selects the candidate solution with better utility by deleting individuals with smaller R2 index contributions. During the update process, the algorithm balances its global and local search capabilities by linearly decreasing the inertia weights and gravity coefficients. The update mechanism of the external archive, referred to as CR2, ensures both the renewal of non-dominated solutions and the preservation of their distribution. Guided by these concepts, the R2-ICRMOAPO algorithm is constructed, offering an innovative method for addressing Multi-Objective Problems (MOPs). Figure 1 illustrates the main design concept of the R2-ICRMOAPO algorithm.

Figure 1.

The main design idea for R2-ICRMOAPO.

The detailed steps for the R2-ICRMOAPO algorithm are described below:

Algorithm 2 R2-ICRMOAPO
Input: population size M, archive size Q, reference point $q^{*}$
Output: Q
1: randomized generation of individual initial velocities and positions 2: initialize $q^{*}$
3: calculate r by (11) 4: identify non-dominant individuals to be stored in the archive
6: calculate the individual mass $m_{i}$ by (12)
7: calculate the individual virtual force $F_{ij, k}$ by (4)
9: update the velocity $v_{i}$ by (6)
10: update the position $x_{i}$ by (7)
11: update archive according to individual CR2
12: while iter< $Ma x_{iteratives}$ do 13: calculate r and sort the individuals 14: update $v_{i}$ and $x_{i}$ 15: end 16: return Q

Experimental study and analysis

Test functions and comparison algorithm

All operational tests are performed at Matlab-R2020b. In this paper, we evaluate the performance of the R2-ICRMOAPO by comparing it with CMOEAD,¹⁸ NSGA-II,⁹ CMOAPO,⁷ and C-TAEA.¹⁹ Experiments were conducted on the CF²⁰ test functions and DC-DTLZ²¹ test functions. By comparing the decomposition-based optimization algorithm CMOEAD and the non-dominated ordering-based optimization algorithm NSGA-II, we demonstrate the feasibility of the indicator-based approach of R2-ICRMOAPO. Furthermore, by comparing the constrained multi-objective-based optimization algorithms CMOAPO and C-TAEA, we analyze the role of improved strategies in R2-ICRMOAPO. The DC-DTLZ functions, classical test functions, and CF benchmark functions, an improved set of test functions, are widely used for evaluating algorithm performance. Experimentation with these functions provides a comprehensive assessment of an algorithm’s capabilities.

For consistency, the population size is set to 100 for all algorithms, with a maximum archive size of 100. Table 1 shows the parameter settings for the R2-ICRMOAPO algorithm and the comparison algorithms.

Table 1.

Algorithm parameter setting.

Algorithm name	Parameter
CMOEAD	$p_{c} = 1$ , $p_{m} = 1 / n$ , $η_{c} = 20$ , $η_{m} = 20$
NSGA-II	$p_{c} = 0.9$ , $p_{m} = 1 / n$ , $η_{c} = 20$ , $η_{m} = 20$
CTAEA	$p_{c} = 1$ , $p_{m} = 1 / n$
CMOAPO	$G = 10$
R2-ICRMOAPO	$G_{i} = 70$ , $G_{f} = 1$

Performance indicator

In the experiment, HV and Inverted Generational Distance (IGD) are used as performance evaluation metrics for algorithms. The algorithm is then comprehensively evaluated from the perspectives of convergence and distribution, illustrating its effectiveness.

HV measures the quality of the algorithm to identify a non-dominated solution set within the solution space. A larger HV value indicates better convergence and distribution of the solution set. HV can be expressed as:

HV = δ (⋃_{i = 1}^{| S |} v_{i})

(15)

Where $δ$ denotes the Lebesgue measure, $s$ denotes the non-dominated solutions, and $v_{i}$ represents the hypervolume, which is composed of reference points with the ith solution in the collection of solutions.

The IGD metric allows for a comprehensive evaluation of the solution set’s performance by measuring the average of the smallest distances from the true Pareto optimal solutions to the approximate solution set. The IGD can be expressed as:

IGD = \frac{(\sum_{i = 1}^{n} \bar{d_{i}})}{n}

(16)

Where $\bar{d_{i}}$ is the Euclidean distance between each individual in the target space and the nearest individual in the true Pareto front, and $n$ denotes the number of individuals in the target space.

Experimental results

To graphically illustrate the optimization ability of the R2-ICRMOAPO algorithm for solving multi-objective problems, Figure 2 shows the Pareto boundary graphs obtained by the R2-ICRMOAPO algorithm and the comparison algorithm for CF1, CF2, and CF6 test problems. Figure 3 plots the Pareto frontier graphs obtained by five algorithms on the CD1-DTLZ1, DC2-DTLZ2, and DC3-DTLZ3 test problems. In Figure 2, the CMOEAD and CMOAPO algorithms demonstrated uneven distribution on the CF1 and CF6 test problems and failed to converge to the Pareto Front (PF). For the CF6 test problem, no algorithm fully converged to the PF upon completion. Nonetheless, the R2-ICRMOAPO algorithm exhibited superior distribution compared to the comparison algorithm. As can be seen from Figure 3, on the DC1-DTLZ1 test problem, the C-TAEA algorithm showed an uneven population distribution with significant gaps, while the NAGAII algorithm displayed a more uniform distribution than the others. The R2-ICRMOAPO algorithm demonstrated enhanced distribution capabilities in the DC2-DTLZ1 and DC3-DTLZ3 test problems. A comprehensive analysis reveals that the R2-ICRMOAPO algorithm optimizes the test problems with a frontier that approaches the true PF more closely and offers a more uniform distribution than the other algorithms.

Figure 2.

The compared Pareto frontiers of the CF1, CF2, and CF6 test functions.

Figure 3.

The compared Pareto frontiers of the DC1-DTLZ1, DC2-DTLZ2, and DC3-DTLZ3 test functions.

Without loss of generality, 30 independent experiments were conducted on 13 test problems, recording the IGD and HV values for each run. Tables 2 and 3 detail the HV and IGD comparison results between the R2-ICRMOAPO algorithm and other comparison algorithms on the CF and DC-DTLZ family of test functions, with the best results for each test function highlighted in bold. Moreover, the Wilcoxon rank sum test has been applied to determine if there is a significant difference between the R2-ICRMOAPO algorithm and the other algorithms. The symbols +, −, and ≈ indicate that the comparative algorithms are better, worse, or equivalent to the R2-ICRMOAPO algorithm in terms of performance, respectively. The final statistics are presented in the last row of the table.

Table 2.

HV indicator test results for the test function.

	CMOEAD	NSGA-II	C-TAEA	CMOAPO	R2-ICRMOAPO
CF1	4.2290e-1(1.74e-2)-	4.6000e-1(9.57e-3)-	4.8576e-1(9.07e-3)-	3.9895e-1(1.13e-2)-	5.1612e-1(4.12e-3)
CF2	5.4359e-1(3.38e-2)-	5.5276e-1(3.51e-2)-	5.7826e-1(2.47e-2)-	5.3555e-1(2.56e-2)-	5.8452e-1(2.16e-2)
CF3	6.3113e-2(2.90e-2)-	7.6051e-2(3.53e-2)-	6.7550e-2(3.32e-2)-	3.1841e-2(3.66e-2)-	8.4917e-2(1.38e-2)
CF4	3.2828e-2(4.89e-2)-	3.4033e-1(3.57e-2)-	3.7302e-1(4.23e-2)≈	2.9734e-1(4.17e-2)-	3.7699e-1(3.48e-2)
CF5	1.5796e-1(5.81e-2)-	2.0747e-1(7.09e-2)≈	1.9502e-1(6.98e-2)-	1.9911e-1(6.57e-2)-	2.0480e-1(5.63e-2)
CF6	5.4580e-1(4.79e-2)-	5.8870e-1(3.03e-2)-	6.0115e-1(2.76e-2)+	5.4338e-1(3.02e-2)-	5.9765e-1(2.17e-2)
CF7	3.0061e-1(7.83e-2)-	3.7028e-1(7.56e-2)-	3.6286e-1(9.54e-2)≈	3.3451e-1(8.31e-2)-	3.9140e-1(7.53e-2)
DC1-DTLZ1	2.8527e-1(2.51e-1) ≈	3.7260e-1(2.43e-1)+	1.0167e-2(6.33e-2)-	1.7838e-1(7.95e-1)-	2.6038e-1(2.76e-1)
DC1-DTLZ3	2.3559e-3(1.19e-2)-	1.2893e-3(5.37e-3)-	1.4637e-3(4.87e-3)-	2.2035e-2(6.18e-2)-	4.6292e-2(4.09e-3)
DC2-DTLZ1	3.6754e-1(6.24e-2)-	4.8702e-1(5.32e-1)-	3.1482e-1(9.90e-2)-	6.3555e-1(9.89e-2)-	7.9291e-1(4.11e-2)
DC2-DTLZ3	1.0796e-1(2.07e-2)-	6.7578e-2(3.17e-3)-	3.4591e-2(4.68e-3)-	2.5720e-2(4.76e-3)-	1.4341e-1(1.95e-3)
DC3-DTLZ1	2.7385e-1(1.40e-2)+	1.5671e-2(8.07e-2)-	1.0805e-1(4.99e-3)-	6.7180e-2(3.79e-3)-	1.1737e-1(1.72e-1)
DC3-DTLZ3	3.8247e-2(3.15e-4)-	2.2731e-2(2.57e-3)-	9.5074e-3(2.38e-3)-	1.9324e-2(7.13e-2)-	3.9835e-2(1.66e-3)
+/−/≈	1/11/1	1/11/1	1/10/2	0/12/0

Table 3.

IGD indicator test results for test functions.

	CMOEAD	NSGA-II	C-TAEA	CMOAPO	R2-ICRMOAPO
CF1	7.1329e-1(5.41e-3)-	6.9409e-2(7.39e-3)-	6.7077e-2(8.81e-3)-	9.2755e-2(9.31e-3)-	4.1288e-2(3.50e-3)
CF2	1.0110e-1(2.80e-2)-	8.4006e-2(2.26e-2)-	6.8745e-2(2.01e-2)-	9.3894e-2(2.66e-2)-	6.7028e-2(1.31e-2)
CF3	5.3041e-1(1.22e-1)-	3.4817e-1(8.94e-2)≈	4.0735e-1(1.20e-1)-	5.3117e-1(1.19e-1)-	3.4754e-1(8.21e-2)
CF4	1.8667e-1(7.58e-2)-	1.2408e-1(3.42e-2)+	1.2949e-1(3.60e-2)≈	1.5402e-1(3.64e-2)-	1.3210e-1(3.52e-2)
CF5	4.0264e-1(9.94e-2)≈	3.7692e-1(1.09e-1)-	3.9753e-1(9.77e-2)≈	3.8748e-1(9.92e-2)≈	3.5117e-1(6.49e-2)
CF6	2.1306e-1(4.26e-2)-	1.0684e-1(3.98e-2)-	1.3217e-1(4.86e-2)-	1.2033e-1(3.96e-2)-	1.0632e-1(2.55e-2)
CF7	1.9278e+0(2.74e-1)-	3.5723e-1(1.21e-1)-	3.8489e-1(1.41e-1)-	3.9145e-1(1.35e-1)-	3.1439e-1(1.13e-1)
DC1-DTLZ1	2.2387e-1(2.23e-2)≈	1.4556e-1(1.80e-2)+	8.5587e-1(6.20e-1)-	1.6028e+0(8.40e-1)-	2.4065e-1(2.92e-1)
DC1-DTLZ3	1.4199e+1(5.50e-1)-	9.6957e-1(4.31e-1)≈	3.7401e+0(6.15e-1)-	4.3728e+0(5.09e-1)-	9.1825e-1(4.22e-1)
DC2-DTLZ1	8.9051e-2(9.44-2)-	2.3138e-1(5.51e-2)-	3.7783e-2(1.54e-2)+	2.0851e-1(3.72e-2)≈	1.5572e-1(8.96e-2)
DC2-DTLZ3	9.3187e-1(4.41e-1)-	9.6028e-1(4.82e-2)-	6.3144e-1(3.97e-1)-	9.0247e-1(4.11e-1)-	6.2733e-1(3.65e-2)
DC3-DTLZ1	7.1497e-1(1.41e-1)-	8.8718e-1(3.17e-1)-	8.0256e-1(3.22e-1)-	6.7477e-1(5.68e-1)-	1.4101e-2(2.47e-1)
DC3-DTLZ3	7.6057e-1(5.62e-1)-	4.0305e-1(5.32e-1)-	6.0729e-1(6.18e-1)-	1.9086e+0(7.00e-1)-	2.0251e-1(1.88e-2)
+/−/≈	0/11/2	2/19/2	0/13/0	1/11/1

From Table 2, it is evident that the R2-ICRMOAPO algorithm demonstrates superior performance compared to other algorithms in the CF1-5, CF7, DC1-DTLZ3, DC2-DTLZ1, DC2-DTLZ3, and DC3-DTLZ3. This suggests that R2-ICRMOAPO has a better overall performance with respect to the HV metric.

Regarding Wilcoxon’s rank-sum test results, the R2-ICRMOAPO algorithm significantly outperformed its competitors in more than half of the test functions. The statistical results clearly indicate that, for the IGD indicator, R2-ICRMOAPO achieves better overall performance.

To specifically analyze the convergence of the R2-ICRMOAPO algorithm, the upward trend curves of the HV metrics for each algorithm under various test problems are plotted in Figure 4. The horizontal and vertical coordinates indicate the maximum assessment number and HV index value for each run, respectively. As shown in Figure 4, the R2-ICRMOAPO algorithm demonstrates an iterative tendency to outperform the comparison algorithms on most test problems. Specifically, in CF1, CF3, CF5, DC2-DTLZ1, and DC3-DTLZ3, the R2-ICRMOAPO algorithm not only rises most rapidly but also reaches a stable state when the number of evaluations is approximately 2 × 10⁴. Furthermore, it achieves the largest HV value at the end of the evaluations. In contrast, the CMOAPO algorithm achieves the minimum HV value when the evaluation times are maximized, thereby highlighting the effectiveness of the improvement measures adopted by the R2-ICRMOAPO algorithm.

Figure 4.

HV indicator convergence curve.

Although no optimization algorithm can solve all problems optimally, the R2-ICRMOAPO algorithm demonstrates high competitiveness in addressing MOPs. The experimental results indicate that R2-ICRMOAPO often finds the best solution among the tested problems, suggesting its effectiveness in solving MOPs.

Multi-objective optimization of hydro-pneumatic suspension parameters

As a piece of indispensable mechanical equipment in modern engineering projects, wheeled excavators play an irreplaceable role in both civil municipal and military construction. Increasing the machine’s speed offers significant practical benefits for military and civilian applications. However, higher travel speeds can amplify road surface input excitation, leading to excessive vibration. Currently, wheeled excavators rely solely on cab and seat suspension for comfort, with no suspension between axles for vibration damping. This design is fine at low speeds, but at speeds above 50 km/h, it generates excessive road vibration so a suspension system must be introduced to dampen it.²²

Hydro-pneumatic suspension, utilizing gas as the elastic medium and liquid as the force transmission medium, can effectively dampen vibrations, and stabilize the excavator at high speeds. It is necessary to enhance the suspension performance evaluation indices based on the excavator’s smoothness and driving comfort under varying road conditions. However, these performance evaluation indicators are not only constrained by each other but improving one or two of them is likely to cause deterioration in the other.^23,24 Therefore, it is necessary to seek the best balance between the contradictions by optimizing the design and adjusting its parameters appropriately. This optimization process of finding an optimal solution among multiple mutually constrained objective functions is called a MOP. In this work, a mathematical model is developed using the vertical body acceleration and tire dynamic load as objective functions. The R2-ICRMOAPO algorithm is used to optimize the model in the MATLAB environment. The optimization of the main parameters of the hydro-pneumatic suspension is realized interactively with Simulink.

Establishment of multi-objective optimization model

The structure of the single-chamber hydro-pneumatic suspension system is shown in Figure 5 below.

Figure 5.

The Single-chamber hydro-pneumatic suspension structure.

During the compression stroke, a portion of oil in the piston chamber enters and fills the annular chamber through a damping orifice and check valve, while the remaining oil simultaneously compresses the gas. Conversely, during the recovery stroke, the check valve closes. The gas flowing out of the annular chamber through the damping orifice is insufficient to fill the space created by the cylinder’s elongation, causing the gas to expand. The compression and expansion of the nitrogen chamber then act as a gas spring, and the flow of hydraulic oil in the orifice is equivalent to a damping element, which will rapidly attenuate the vibration energy transmitted to the body and transform it into heat energy.

The design goal of hydro-pneumatic suspension is to minimize the car-body vertical acceleration and dynamic tire load as much as possible within a limited working space. In this paper, half of the Front axle model of a high-speed wheel excavator is abstractly simplified, and the Two degree-of-freedom suspension vibration system of a high-speed wheel excavator is obtained as shown in Figure 6.

Figure 6.

Two-degree-of-freedom suspension vibration system structure.

Newton’s law of motion is adopted to derive the motion equation of the vibration system, as shown below:

M_{w} {z ″}_{1} = K_{t} (x_{0} - z_{1}) - K_{s} (z_{1} - z_{2}) - C_{s} (z_{1}' - z_{2}')

(17)

M_{b} {z ″}_{1} = K_{s} (z_{1} - z_{2}) + C_{s} (z_{1}' - z_{2}')

(18)

Where $K_{t}$ represents the tire stiffness, and $C_{s}$ denotes the damping coefficient.

In light of the practical challenges encountered in production design, the external structural parameters of the hydro-pneumatic suspension cylinders remain unchanged. Instead, the initial gas volume and pressure, the damping hole diameter, and the check valve diameter are selected as design variables.^25,26 That is $x = [P_{0}, V_{0}, d_{k}, d_{f}]$ .

According to the actual situation and experience, the value range of design variables is as follows (Table 4):

Table 4.

Design variable value range.

Variable	$P_{0}$ (L)	$V_{0}$ (MPa)	$d_{k}$ (mm)	$d_{f}$ (mm)
Initial value	5.76	0.50	3.00	6.00
Upper limit	6.34	0.55	5.00	8.00
Lower limit	5.18	0.45	2.00	4.00

Car-body vertical acceleration, dynamic tire load, and suspension dynamic forgiveness are the three primary indicators for assessing vehicle driving smoothness. However, effectively unifying these indicators for road excitation is challenging due to the inherent constraints between them. Therefore, this paper proposes the establishment of a suspension parameter optimization model. In this model, car-body vertical acceleration and dynamic tire load are considered as the objective functions, while suspension dynamic forgiveness serves as the constraint condition. This approach aims to find the best compromise among the competing objectives, leading to an optimized design that enhances overall vehicle smoothness.

The objective function expression of suspension multi-objective optimization is given by:

min (a_{w}, σ_{f})

(19)

Where $a_{w}$ represents the vertical body acceleration, and $σ_{f}$ is the dynamic tire load.

When a wheel excavator operates at high speeds, it is crucial to maintain sufficient tire-ground contact to prevent the risk of the machine becoming airborne. To mitigate this risk, the dynamic tire load is typically restricted to one-third of the total suspension mass. Studies indicate that under these conditions, the likelihood of the wheel dynamic load surpassing the safety threshold is only 0.3%, and the chance of wheels losing contact with the ground is a mere 0.15%.

σ_{f} \leq \frac{1}{3} G

(20)

The dynamic deflection of the suspension system is restricted to prevent the suspension cylinder from hitting the limit device. This constraint is necessary because when the suspension dynamic stroke is excessive, the likelihood of the suspension striking the limit block increases. Such a situation could occur if the amplitude of suspension motion becomes too large.

σ_{n} \leq \frac{1}{3} [f_{d}]

(21)

Where $[f_{d}]$ is the allowable value of the maximum dynamic deflection of the suspension system, with a selected limit of 120 mm.

Therefore, the calculation model for the optimization algorithm of the hydro-pneumatic suspension is as follows:

\begin{matrix} min f_{1} = a_{w} \\ min f_{2} = σ_{f} \\ s . t . σ_{f} \leq \frac{1}{3} G \\ σ_{n} \leq \frac{1}{3} [f_{d}] \\ 5.18 \leq P_{0} \leq 6.34 \\ 0.45 \leq V_{0} \leq 0.55 \\ 2 \leq d_{k} \leq 5 \\ 4 \leq d_{f} \leq 8 \end{matrix}

(22)

Using equation (22), we establish the corresponding Simulink model and invoke the R2-ICRMOAPO algorithm optimization program written in the MATLAB environment to obtain the optimization results.

Optimization results

The R2-ICRMOAPO algorithm is employed to estimate the optimal values of the vehicle suspension system parameters, effectively minimizing the cost toward its optimal value. The vehicle model parameters utilized in the study are outlined in Table 5.

Table 5.

Vehicle suspension system parameters.

Name	Unit	Numerical
Sprung mass	kg	2950
Unsprung mass	kg	650
Initial gas volume	L	0.5
Initial gas pressure	MPa	5.76
Piston rod outer diameter	mm	100
Cylinder inner diameter	mm	80
Check valve equivalent diameter	mm	6
Damping hole apertures	mm	3
Tire stiffness	N/m	1,000,000

The simulation model of a dual-degree-of-freedom suspension system was established in Simulink, as depicted in Figure 7, based on the established multi-objective model.

Figure 7.

Simulink model.

According to the model, the excitation generated by the wheeled excavator traveling at a speed of 15 m/s on a Class D road surface is considered as the road disturbance. The corresponding R2-ICRMOAPO algorithm optimization program has been developed, and the optimization results obtained after running the program are shown in Table 6.

Table 6.

Parameter optimization results.

Variable	$P_{0}$ /L	$V_{0}$ /MPa	$d_{k}$ /mm	$d_{f}$ /mm
Initial value	5.76	0.50	3.00	6.00
Optimization value	5.49	0.54	2.63	5.86

The optimized parameters of the hydro-pneumatic suspension were incorporated into the suspension system model, and the ride comfort was evaluated. Comparisons of the ride comfort before and after optimization, including the vertical acceleration, dynamic deflection, and dynamic load on the bodywork of the vehicle, were made under identical operating conditions (Class D road surface, with a driving speed of 15 m/s). The results are presented in Figures 8 to 10.

Figure 8.

Optimize front and rear car-body vertical acceleration curves: (a) before optimization and (b) after optimization.

Figure 9.

Optimize front and rear dynamic tire load curves: (a) before optimization and (b) after optimization.

Figure 10.

Optimize front and rear suspension dynamic forgiveness curves: (a) before optimization and (b) after optimization.

Figures 8 to 10 illustrates that the vertical acceleration of the car-body, dynamic tire load, and suspension dynamic forgiveness of the hydro-pneumatic suspension are improved after the algorithm is optimized. Table 7 depicts the root mean square comparison of the metrics before and after optimizing the hydro-pneumatic suspension.

Table 7.

Comparison of root mean square values of each indicator before and after optimization.

Indicators	Car-body vertical acceleration (m/s²)	Dynamic tire load (N)	Suspension dynamic forgiveness (mm)
Before optimization	1.351	4839.671	7.868
After optimization	1.062	3890.274	7.492
Performance enhancements	21.4%	19.6%	4.72%

Table 7 shows the improvements achieved after optimizing the wheeled excavator’s suspension system. The vertical acceleration of the car body was significantly reduced, 1.351 to 1.062 m/s², representing a 21.4% decrease. Similarly, the dynamic wheel load experienced a 19.6% reduction, dropping from 4839.671 to 3890.3274 N. The suspension dynamic forgiveness also saw a 4.72% reduction, decreasing from 7.868 to 7.492 mm. The optimization not only significantly reduces the vertical acceleration and tire dynamic load but also ensures that the suspension’s dynamic camber remains within a reasonable range. The results suggest that the R2-ICRMOAPO algorithm effectively enhances the smoothness of wheeled excavators, offering a theoretical foundation for further suspension system optimization.

Conclusions and future studies

The APO algorithm has received widespread attention as an emerging SI optimization algorithm. The R2-ICRMOAPO algorithm proposed in this paper obtains well-distributed solutions. Based on APO’s high search efficiency, the algorithm takes advantage of the dominant relationship between individuals and the crowding degree of individuals, selects individuals based on the contribution value of R2 indicators. This approach enhances population distribution by rationally ordering individuals from the search space for retention and updating in the external storage set. Furthermore, the inertia weights and gravity factors are dynamically adjusted to optimize the global and local search capability of the R2-ICRMOAPO algorithm. By analyzing the test results in CF and DC-DTLZ test functions and Wilcoxon rank sums, it is shown that the algorithm is highly competitive.

We utilized the R2-ICRMOAPO algorithm to further optimize the parameter design of the hydro-pneumatic suspension system, formulating a constrained multi-objective optimization model. In this model, our primary objectives were to minimize the vertical acceleration of the car body and the dynamic loads on the tires. The main structural parameters of the hydro-pneumatic suspension system were considered as variables. The results show that the parameters optimized by the R2-ICRMOAPO algorithm reduce the root mean square value of car-body vertical acceleration by 21.4%, the root mean square value of dynamic tire load by 19.6%, and the smoothness of the optimized performance is greatly improved. The optimization results may offer a theoretical foundation that can be used for improving the hydro-pneumatic suspension.

In the future, we are prepared to further extend the algorithm to the high-dimensional domain and apply it to a broader range of engineering fields. This is because the R2 index introduced in this paper is particularly suitable for solving high-dimensional MOP problems due to its weak Pareto dominance property and computational simplicity. Meanwhile, numerous problems in engineering involve high dimensions.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by National Natural Science Foundation of China (Grant No. 52302475), Basic Research Programs of Shanxi Province (Grant No. 202203021221150).

ORCID iD

Yihong Chang

Data availability statement

Data sharing not applicable to this article as no datasets.

References

Soheili

Zoka

Abachizadeh

Tuned mass dampers for the drift reduction of structures with soil effects using ant colony optimization. Adv Struct Eng 2021; 24(4): 771–783.

Han

Chen

Ling

, et al. Multi-objective particle swarm optimization with adaptive strategies for feature selection. Swarm Evol Comput 2021; 62: 100847.

Wei

, et al. A survey on firefly algorithms. Neurocomputing 2022; 500: 662–678.

Barman

Dewang

Mewada

Facial recognition using grey wolf optimization. Mater Today Proc 2022; 58: 273–285.

Molina-Pérez

Mezura-Montes

Portilla-Flores

, et al. A differential evolution algorithm for solving mixed-integer nonlinear programming problems. Swarm Evol Comput 2024; 84: 101427.

Lin

Chen

, et al. Survey of novel swarm intelligence optimization algorithms. Comput Eng Appl 2018; 54(12): 1–9.

Wang

Zeng

Tan

. An artificial physics optimization algorithm for multi-objective problems based on virtual force sorting proceedings. In: Swarm, evolutionary, and memetic computing: first international conference on swarm, evolutionary, and memetic computing, SEMCCO 2010, Chennai, India, December 16-18, 2010. Proceedings 1, 2010, pp.615–622. Berlin, Heidelberg: Springer.

Sun

, et al. Artificial physics multi-objective algorithm based on sequence value and crowding degree. Syst Eng Electron 2014; 36(12): 2442–2448.

Deb

Pratap

Agarwal

, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans Evol Comput 2002; 6(2): 182–197.

10.

Tan

Zheng

, et al. A decomposition method based on random objective division for MOEA/D in many-objective optimization. IEEE Access 2020; 8: 103550–103564.

11.

Dong

Liu

, et al. An IGD+ Performance Indicator Based Particle Swarm Optimizer For Multi-objective Optimization. In: 2021 33rd Chinese Control and Decision Conference (CCDC), pp. 3633–3638. IEEE, 2021.

12.

Dhiman

Soni

Pandey

, et al. A novel hybrid hypervolume indicator and reference vector adaptation strategies based evolutionary algorithm for many-objective optimization. Eng Comput 2021; 37(4): 3017–3035.

13.

Jiang

, et al. Multi-objective particle swarm optimization with R2 indicator and adaptive method. Complex Intell Syst 2021; 7(5): 2697–2710.

14.

Zitzler

Knowles

Thiele

Quality assessment of pareto set approximations. Multiobject Optim 2008; 5252(4): 373–404.

15.

Wei

Fan

, et al. A hybrid multiobjective particle swarm optimization algorithm based on R2 indicator. IEEE Access 2018; 6: 14710–14721.

16.

Liu

, et al. An R2 indicator and weight vector-based evolutionary algorithm for multi-objective optimization. Soft Comput 2020; 24(7): 5079–5100.

17.

Cheng

Liu

, et al. A two-stage R2 indicator based evolutionary algorithm for many-objective optimization. Appl Soft Comput 2018; 67: 245–260.

18.

Cai

Wang

A multiobjective optimization-based evolutionary algorithm for constrained optimization. IEEE Trans Evol Comput 2006; 10(6): 658–675.

19.

Zou

Shi

Guo

, et al. Solving multi-stage weapon target assignment problems by C-TAEA. In: 2020 39th Chinese control conference (CCC), Shenyang, China, 2020, pp.1593–1598. DOI: 10.23919/CCC50068.2020.9188559

20.

Zhu

Zhang

Lin

A constrained multiobjective evolutionary algorithm with detect-and-escape strategy. IEEE Trans Evol Comput 2020; 24(5): 938–947.

21.

Shan

XY.

Two-population evolution algorithm for coordinated constrained multi-objective optimization. University of Electronic Science and Technology of China, 2021. DOI: 10.27005/d.cnki.gdzku.2021.005504.

22.

Sun

Liu

Zhang

, et al. Analysis on the research status of oil and gas suspension system control technology of heavy-duty engineering vehicles. Trans China Soc Eng Mach 2021; 19(06): 482–486.

23.

Yang

, et al. Innovative Design Method of hydro-pneumatic suspension for large high-clearance sprayer based on improved NSGA-II algorithm. Agriculture 2023; 13(5): 1071.

24.

Tian

Xia

XX . Vibration optimization and control of seven-axle hybrid oil-gas suspension vehicle. Mech Strength 2022; 44(01): 19–28.

25.

Tao

Fan

You

Investigation on non-linear characteristics of hydropneumatic suspension modeling and simulation for engineering vehicles. Acta Simulata Syst. Sin 2005; 1(2): 55.

26.

Tian

Dou

XX.

Modelling and characteristics of hybrid coupling hydro-pneumatic suspension for a seven-axle vehicle. Proc IMechE, Part D: J Automobile Engineering 2021; 235(7): 1892–1901.