# Multiobjective Design Optimization Framework for Multicomponent System with Complex Nonuniform Loading

Authors: Hong Zhang; Guangchen Bai; Lukai Song
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

To improve the accuracy and efficiency of multiobjective design optimization for a multicomponent system with complex nonuniform loads, an efficient surrogate model (the decomposed collaborative optimized Kriging model, DCOKM) and an accurate optimal algorithm (the dynamic multiobjective genetic algorithm, DMOGA) are presented in this study. Furthermore, by combining DCOKM and DMOGA, the corresponding multiobjective design optimization framework for the multicomponent system is developed. The multiobjective optimization design of the carrier roller system is considered as a study case to verify the developed approach with respect to multidirectional nonuniform loads. We find that the total standard deviation of three carrier rollers is reduced by 92%, where the loading distribution is more uniform after optimization. This study then compares surrogate models (response surface model, Kriging model, OKM, and DCOKM) and optimal algorithms (neighbourhood cultivation genetic algorithm, nondominated sorting genetic algorithm, archive microgenetic algorithm, and DMOGA). The comparison results demonstrate that the proposed multiobjective design optimization framework is demonstrated to hold advantages in efficiency and accuracy for multiobjective optimization.


---

## Body

## 1. Introduction

Multicomponent system is defined as the complex mechanism system comprising a plurality of rigid and flexible components, which is an indispensable part in mechanical equipment, such as excavator and loader [1–3]. For example, the carrier roller system of the excavator is a typical multicomponent system. Such multicomponent system, however, suffers from nonuniform loads caused by complex structural layout and stricter working environment, and the failure of one component will lead to the failure of the whole component system, which significantly increases the failure possibility of the multicomponent system. In addition, the multicomponent system tends to have the mutual conduction effect between components, and the design and control of one structural objective often leads to the change of another. Therefore, to reduce the load nonuniformity and improve the system performance, it is necessary to perform the multiobjective design optimization (MODO) for the complex multicomponent system. A typical MODO problem can be formulated as(1)FindxminσfixS.t.∑Fj<∑Fjglx=0au≤xu≤bu,where fi(x) denotes the i-th objective function; Fj is the j-th constraint function; [Fj] is the allowable load under j-th constraint function; glx is the equality constraint of the l-th design parameter x; xu is the u-th design parameter; and au and bu are the upper and bottom boundaries of the parameter x, respectively; i = 1, 2, …, m; j = 1, 2, …, n; u = 1, 2, …, p; l = 1, 2, …, q.To solve the MODO model shown in (1), the multiple response simulations (i.e., objective response simulation and constraint response simulation) and multiobjective optimal algorithm are unavoidable to be required. Under such circumstances, it will require thousands of calls to expensive black-box function to simulate the multiobjective and multiconstraint, which will lead to huge computing task and unaffordable computing cost. To improve the computational efficiency, differentiation algorithms, generalized speeds, Euler parameters, lower body arrays, and other dynamic numerical methods are proposed and widely applied [4–6]. Unfortunately, although these numerical methods can greatly reduce simulation time, the approximation accuracy of multiobjective and multiconstraint is unacceptable for the complex multicomponent system. Moreover, the nonlinear constrain space and multiple optimal solutions usually demand nested iterative loops, which also brings the solving difficulties of the MODO model. To improve MODO accuracy and efficiency of the complex multicomponent system, this study presents two key techniques: (1) establishment of a numerical surrogate model to calculate the multiobjective and multiconstraint; (2) development of a dynamic multiobjective algorithm to resolve the MODO model.The first technique is to establish an efficient surrogate model to approximate the multiobjective and multiconstraint. Only needing a small amount of black-box function calls, the surrogate model can be established and is promising to reduce simulation cost [7–11]. Classic surrogate models include polynomial response surface [12–14], Kriging model (KM) [15–17], artificial neural network [18–20], and support vector regression [21–23]. Among them, with integrating global nonlinear approximation ability and local precise description ability, the Kriging model possesses the potentials to approximate complex structural responses and ensure calculation accuracy [24, 25], which is suitable for complex multiobjective optimization problems of the multicomponent system. However, with regard to high-nonlinearity traits of multiobjective and multiconstraint problems, the traditional Kriging model is still insufficient to acquire satisfactory efficiency and accuracy. To address this problem, in this study, by utilizing an improved particle swarm optimization algorithm [26] to search Kriging key parameters and decomposed collaborative strategy [27–30] to simplify calculation complexity, a decomposed collaborative optimized Kriging model (DCOKM) is proposed to complete efficient approximation of multiobjective and multiconstraint.The second technique is to develop an accurate multiobjective optimal algorithm to resolve the complex MODO model. Due to the complex interaction effects between various objective functions in the multicomponent system, the global optimal solution of the MODO model is almost impossible to acquire at the same time. Therefore, an accurate optimal algorithm, which can precisely search for Pareto optimal frontier and acquire nondominated solutions, is urgently required. At present, a variety of multiobjective optimization algorithms are developed, such as multi-island genetic algorithm (MIGA) [31–33], neighbourhood cultivation genetic algorithm (NCGA) [34, 35], the second generation of nondominated sorting genetic algorithm (NSGA-II) [36–38], and archive microgenetic algorithm (AMGA) [39, 40]. Among them, owing to the virtues of nondominated sorting strategy and elitist retention strategy, NSGA-II has been widely used in MODO problems. Unfortunately, because of easy to fall into local solutions or premature convergence, the current NSGA-II algorithm is still difficult to acquire global optimal solution of the multicomponent system [41, 42]. In this case, by designing an arithmetic crossover operator and Poisson mutation operator in NSGA-II algorithm, a dynamic multiobjective genetic algorithm (DMOGA) is presented to accurately solve the MODO model and acquire nondominated solutions.The objective of this paper is to develop an efficient and accurate multiobjective design optimization framework for the multicomponent system with nonuniform loads. By fusing the proposed surrogate model (i.e., DCOKM) and optimal algorithm (i.e., DMOGA), the urgently needed multiobjective design optimization framework is established. The innovation of this paper lies in that the corresponding efficient methods are firstly proposed for the multiobjective optimization design of the multicomponent system with nonuniform loading. Regarding the multiobjective design optimization of the carrier roller system as study case, the effectiveness of the presented framework is validated. In what follows, Section2 investigates surrogate modeling with DCOKM, optimal algorithm with DMOGA, and the corresponding multiobjective design optimization framework. Section 3 performs the study case to validate the multiobjective design optimization framework. Some conclusions are summarized in Section 4.


## 2. Basic Theories and Formulations

### 2.1. Decomposed Collaborative Optimized Kriging Model, DCOKM

#### 2.1.1. Optimized Kriging Model, OKM

By integrating the Kriging model (KM) with the best unbiasedness ability and the improved particle swarm optimization (IPSO) with global search capability [43], the optimized Kriging model (OKM) is introduced as(2)yx=fTx⋅β+zx,where x = [x1, x2, …, xd] indicates the input variable; d is the dimension number of input parameters; fT(x) β signifies the global approximation trend; fT(x) = [f1(x), f2(x), …, fm(x)] is the regression basis function; β = [β1, β2, …, βm] is the regression coefficient; m is the dimension number of regression function; and z(x) is the local random deviation, which follows Gaussian distribution (0, σ2). Herein, the covariance measure of z(x) is defined as(3)Covzxi,zxi=σ2Rxi,xj,where xi and xj denote the i-th and the j-th input variables, respectively; R(·,·) is the Gauss correlation function, which is(4)Rxi,xj=exp−∑k=1nθkxik−xjk2=∏k=1nexp−θkxik−xjk2,in which n represents the dimension number of design variables; x(k) i and x(k) j are the k-th member of variables xi and xj, respectively; θk is the k-th member of correlation parameter θ.Moreover, the optimal correlation parameterθ∗ of OKM can be obtained by the following maximum likelihood estimation, i.e.,(5)θ∗=argmaxψθ=arg min12lnR+mlnσ^2,where ψ(θ) indicates the likelihood function.Then, the corresponding regression coefficientβ and variance σ2 of OKM can be estimated by(6)β^=FTR−1F−1FTR−1y,σ^2=1ny−Fβ^TR−1y−Fβ^,where F is the unit column vector and R is the correlation matrix.Clearly, the approximation accuracy of OKM mostly depends on the correlation parameterθ. Once the optimal correlation parameter θ∗ is acquired, the optimal unbiased prediction of OKM is ensured. Therefore, the high-accurate surrogate modeling problem is transformed into the searching problem of the optimal solution θ∗.Searching algorithm is a key factor to find the optimal solutionθ∗, since its searching efficiency and accuracy directly influence the approximation efficacy of the Kriging model. The PSO algorithm is an important nongradient searching algorithm, and it holds good collaborative searching ability [26]. However, the fixed inertia weight and learning factors of the ordinary PSO algorithm often fall into blind searching and trap by local optimum easily, which greatly affects the searching efficiency and accuracy. To address this issue, we propose an improved PSO algorithm by designing adaptive inertia weight and asynchronous learning factors, with which the iteration steps change. The design is to acquire better optimal solution θ∗ by balancing the global searching ability and local searching ability. The basic thought of improved PSO algorithm is summarized as follows.We assume that the particle position is assigned as the correlation parameterθ and the fitness value adopts the ψ(θ) of the optimized Kriging model. Hereinto, each particle is a potential solution for the optimal correlation parameter θ∗. In searching process, with following the direction guidance of optimal particles, the particle swarm will search for the optimal solution by updating particle individual positions, individual extreme values, and population extremum values. The renewal formulas of the improved PSO algorithm are introduced as(7)Vik+1=ωVik+c1r1Pik−Xik+c2r2Pgk−XikXik+1=Xik+Vik+1i=1,2,⋯,N,in which the adaptive inertia weight ω and asynchronous learning factors c are(8)ω=ωmax−ωmax−ωminf−fminfavg−fmin,f≤favg,ωmax,f>favg,(9)c=c1=c1s+c1n−c1s×k/K,c2=c2s+c2n−c2s×k/K,where i represents the i-th particle; k is the current iteration number; K is the largest iteration number; ω is the inertia weight; Vi is the current particle velocity; Xi is the current particle position; Pi is the current individual extremum; Pg is the current population extremum; c1 and c2 are the individual and population leaning factors, respectively; r1 and r2 are the random numbers during time domain [0,1]; ωmax is the maximum inertia weight; ωmin is the minimum inertia weight; f is the current fitness value; favg is the average fitness value; fmin is the minimum fitness value; c1s and c2s are the initial individual learning factor and initial population learning factor, respectively; and c1n and c2n are the individual learning factor and population learning factor in the largest iterations, respectively.With finite iteration times of the improved PSO algorithm, the specified searching precision is met and the optimal correlation parametersθ∗ are acquired. In light of Kriging interpolation theory, the best unbiased prediction of OKM is achieved, and the predicted response y^x˜ related to the predicted point x˜ outside sample set can be expressed by(10)y^x˜=fTx˜β^+rTx˜R−1y−Fβ^,where r=Rx˜,x1,Rx˜,x2,…,Rx˜,xnTis the correlation vector between an predicted point x˜ and training sample points (x1, x2, …, xn).With the predicted responsey^x˜, the prediction error of OKM can be evaluated by(11)φx˜=σ^21+uTx˜FTR−1F−1ux˜−rTx˜R−1rx˜s.t.ux˜=FTR−1rx˜−fx˜.Noticeably, by combining the nonlinear approximation ability of KM and the powerful global searching ability of the improved PSO algorithm, the presented OKM holds the potential to improve approximate efficiency and accuracy. The flowchart of OKM is drawn in Figure1.Figure 1
Flowchart of OKM. The establishment process of the OKM model is described.


#### 2.1.2. Mathematical Modeling of DCOKM

Multiobjective design optimization of the multicomponent system refers to multiobjective (i.e., the force standard deviation in the first componentσF1, in the third component σF3, and in all rollers σMF, etc.) and multiconstraint (i.e., the sum of mean force in all rollers SMF and the sum of maximum force in all rollers SFmax, etc.), which is difficult to perform the nonlinear approximation with high efficiency and acceptable accuracy. To address the multiobjective and multiconstraint (MOMC) problem, we develop a decomposed collaborative optimized Kriging model (DCOKM) with respect to the decomposed collaborative strategy and optimized Kriging model, to further address the strong nonlinearity issues and enhance approximation efficacy. For the MOMC approximation with the presented DCOKM, an entire complex MOMC system with all input variables and output responses is decomposed into many single-objective single-constraint subsystems, in which each submodel has few input variables and output responses. Considering Latin hypercube sampling technique and OKM, the decomposed OKM for many single-objective single-constraints is established. Subsequently, considering collaborative sampling technique and decomposed OKM, the DCOKM is mathematically constructed. The mathematical modeling process of DCOKM is summarized.Assuming that the multiobjective and multiconstraint problem involvesr objectives and l constraints, when x(t) indicates the t-th input sample, then the corresponding t-th output response Y(t) is(12)Yt=fxt,xt=x1t,x2t,…,xntT,t=1,2,…,p,p+1,…,l,where t = 1, 2, …, p indicates output objective, while t = p + 1, p + 2, …, l represents the constraint response.With a given predicted pointx˜t for (10), the corresponding output response Y˜t is(13)Y˜t=fTx˜tβ^t+rTx˜tRt−1yt−Fβ^t.Equation (13) is the decomposed OKM under p-th objective response and l-th constraint response, in which β^t, R(t), and rTx˜t indicate the optimal regression coefficient, the correlation model, and the correlation vector in decomposed OKM, respectively. From the above analysis, the OKM (10) of multiobjective and multiconstraint models is divided into a series of OKM submodels under single-objective and single-constraint such as (12) and (13), and the approximation complexity is significantly reduced. As the efficient and accurate DCOKM is developed, the MOMC problem is promising to be addressed for multiobjective design optimization.


### 2.2. Dynamic Multiobjective Genetic Algorithm, DMOGA

The resolving algorithm is a critical factor for solving the MODO model and obtaining the optimal solutions. NSGA-II is a vital nondominated fast sorting algorithm in multiobjective optimal design, which can effectively control the population distribution and reduce the resolving complexity. However, due to randomness and uncertain fluctuation effects in simulated binary crossover and polynomial mutation operator, its convergence precision is still unacceptable for the complex MODO model. To tackle with this problem, by designing the arithmetic crossover operator and Poisson mutation operator, we propose a dynamic multiobjective genetic algorithm (DMOGA) based on the traditional NSGA-II. The objective of the DMOGA is to reduce the uncertain fluctuation effects of operators and gain ideal Pareto optimal frontier. The design thought of the arithmetic crossover operator and Poisson mutation operator is presented as follows.To improve the population diversity and avoid over propagation of excellent solutions, an arithmetic crossover operator is presented by combining the information of nondominated sorting levels. The arithmetic crossover operator can keep the high-ranking individuals in parent generation and increase the low-ranking individuals in offspring population. Therefore, the quality of offspring population and the diversity of whole population shall be improved. For two given parent vectorsP1=v11,v21,⋯,vm1 and P2=v12,v22,⋯,vm2, two offspring Cz = (z1, z2, …, zm) and Cw=w1,w2,…,wm can be obtained by the designed arithmetic crossover operator:(14)zi=αivi1+1−αivi2,wi=αivi2+1−αivi1,where vi1and vi2 indicate the i-th parent individual of the first and second generations; zi and wi are the child individual; i = 1, 2, …, m; and αi is the random number in [0,1], which is expanded as(15)αi=BrAr+Br,Ar≠Br,AdAd+Bd,Ar=Br,where Ar and Br indicate the nondominated sorting level of individual A and individual B, respectively and Ad and Bd are the crowding distance of individual A and individual B, respectively.To change the disadvantage of low convergence speed caused by polynomial mutation parameters, a new Poisson distribution mutation operator is proposed, which can overstep the local optimum and thereby is conducive to search for global solutions. Assuming that the individual before mutation is expressed asx = (x1, x2, …, xn) and the individual after compilation is expressed as x′=x1′,x2′,⋯,xn′, then the Poisson mutation operator can be designed as(16)xi′=xi+xiu−xilδk,where xi′ means the i-th child individual after mutation; xi is the i-th parent individual; xui and xli are the upper and bottom boundaries of gene value at the mutation point i, respectively; i = 1, 2, …, m; and δk is the Poisson random number, which is introduced as(17)δk=λke−λk!,where λ denotes the distribution parameter and k is the variable value.To obtain the Pareto optimal frontier and acquire nondominated solutions, the DMOGA algorithm is proposed by designing the arithmetic crossover operator and Poisson mutation operator. The essential process of DMOGA algorithm is summarized as follows:(1)
Generate initialN individuals of population Pn, and set iteration times n = 0(2)
Evaluate fitness values of all individuals, and rank them with Pareto dominance and crowding distances(3)
Select individuals fromPn by the binary tournament method(4)
Obtain the childCn by the arithmetic crossover operator and Poisson mutation operator(5)
Generate populationQn by merging parent and child individuals, evaluate fitness value of all individuals, and perform fast nondominated sorting of Qn(6)
Calculate crowding degree and crowding distance, and select individuals to form the new species groupPn+1(7)
Terminate the algorithm when the accuracy requirement is met; otherwise, the algorithm will back to (2)(8)
Acquire Pareto optimal surface and global optimal solutionBy absorbing the global search capability of the arithmetic crossover operator and fast convergence ability of the Poisson mutation operator, the proposed DMOGA can accurately solve the multiobjective design optimization model and quickly obtain the Pareto solution set. To achieve the tradeoff effect between multiple objectives and optimal solution, a Utopia-Pareto directing adaptive (UPDA) search scheme is adopted by capturing plenty of characteristics and utilizing the ordering information of Pareto solutions [35, 44], which is conducive to find the best optimal solution for multiobjective design optimization of the multicomponent system.


### 2.3. Multiobjective Design Optimization Framework with DCOKM and DMOGA

For the multicomponent system optimal design problem, we develop a multiobjective design optimization framework based on a surrogate model (i.e., DCOKM) to approximate multiple responses and an optimal algorithm (i.e., DMOGA) to acquire optimal solution set. Herein, the improved PSO algorithm and Kriging model are first absorbed into decomposed collaborative strategy, to enhance the approximation efficacy of the surrogate model; then, the arithmetic crossover operator and Poisson mutation operator are designed, to acquire the efficient and accurate convergence of the optimal algorithm. Therefore, by combining the DCOKM and DMOGA, the calculation accuracy and efficiency of the multiobjective design optimization model is promising to be greatly improved. The multiobjective design optimization framework is shown in Figure2.Figure 2
Multiobjective design optimization framework.


## 2.1. Decomposed Collaborative Optimized Kriging Model, DCOKM

### 2.1.1. Optimized Kriging Model, OKM

By integrating the Kriging model (KM) with the best unbiasedness ability and the improved particle swarm optimization (IPSO) with global search capability [43], the optimized Kriging model (OKM) is introduced as(2)yx=fTx⋅β+zx,where x = [x1, x2, …, xd] indicates the input variable; d is the dimension number of input parameters; fT(x) β signifies the global approximation trend; fT(x) = [f1(x), f2(x), …, fm(x)] is the regression basis function; β = [β1, β2, …, βm] is the regression coefficient; m is the dimension number of regression function; and z(x) is the local random deviation, which follows Gaussian distribution (0, σ2). Herein, the covariance measure of z(x) is defined as(3)Covzxi,zxi=σ2Rxi,xj,where xi and xj denote the i-th and the j-th input variables, respectively; R(·,·) is the Gauss correlation function, which is(4)Rxi,xj=exp−∑k=1nθkxik−xjk2=∏k=1nexp−θkxik−xjk2,in which n represents the dimension number of design variables; x(k) i and x(k) j are the k-th member of variables xi and xj, respectively; θk is the k-th member of correlation parameter θ.Moreover, the optimal correlation parameterθ∗ of OKM can be obtained by the following maximum likelihood estimation, i.e.,(5)θ∗=argmaxψθ=arg min12lnR+mlnσ^2,where ψ(θ) indicates the likelihood function.Then, the corresponding regression coefficientβ and variance σ2 of OKM can be estimated by(6)β^=FTR−1F−1FTR−1y,σ^2=1ny−Fβ^TR−1y−Fβ^,where F is the unit column vector and R is the correlation matrix.Clearly, the approximation accuracy of OKM mostly depends on the correlation parameterθ. Once the optimal correlation parameter θ∗ is acquired, the optimal unbiased prediction of OKM is ensured. Therefore, the high-accurate surrogate modeling problem is transformed into the searching problem of the optimal solution θ∗.Searching algorithm is a key factor to find the optimal solutionθ∗, since its searching efficiency and accuracy directly influence the approximation efficacy of the Kriging model. The PSO algorithm is an important nongradient searching algorithm, and it holds good collaborative searching ability [26]. However, the fixed inertia weight and learning factors of the ordinary PSO algorithm often fall into blind searching and trap by local optimum easily, which greatly affects the searching efficiency and accuracy. To address this issue, we propose an improved PSO algorithm by designing adaptive inertia weight and asynchronous learning factors, with which the iteration steps change. The design is to acquire better optimal solution θ∗ by balancing the global searching ability and local searching ability. The basic thought of improved PSO algorithm is summarized as follows.We assume that the particle position is assigned as the correlation parameterθ and the fitness value adopts the ψ(θ) of the optimized Kriging model. Hereinto, each particle is a potential solution for the optimal correlation parameter θ∗. In searching process, with following the direction guidance of optimal particles, the particle swarm will search for the optimal solution by updating particle individual positions, individual extreme values, and population extremum values. The renewal formulas of the improved PSO algorithm are introduced as(7)Vik+1=ωVik+c1r1Pik−Xik+c2r2Pgk−XikXik+1=Xik+Vik+1i=1,2,⋯,N,in which the adaptive inertia weight ω and asynchronous learning factors c are(8)ω=ωmax−ωmax−ωminf−fminfavg−fmin,f≤favg,ωmax,f>favg,(9)c=c1=c1s+c1n−c1s×k/K,c2=c2s+c2n−c2s×k/K,where i represents the i-th particle; k is the current iteration number; K is the largest iteration number; ω is the inertia weight; Vi is the current particle velocity; Xi is the current particle position; Pi is the current individual extremum; Pg is the current population extremum; c1 and c2 are the individual and population leaning factors, respectively; r1 and r2 are the random numbers during time domain [0,1]; ωmax is the maximum inertia weight; ωmin is the minimum inertia weight; f is the current fitness value; favg is the average fitness value; fmin is the minimum fitness value; c1s and c2s are the initial individual learning factor and initial population learning factor, respectively; and c1n and c2n are the individual learning factor and population learning factor in the largest iterations, respectively.With finite iteration times of the improved PSO algorithm, the specified searching precision is met and the optimal correlation parametersθ∗ are acquired. In light of Kriging interpolation theory, the best unbiased prediction of OKM is achieved, and the predicted response y^x˜ related to the predicted point x˜ outside sample set can be expressed by(10)y^x˜=fTx˜β^+rTx˜R−1y−Fβ^,where r=Rx˜,x1,Rx˜,x2,…,Rx˜,xnTis the correlation vector between an predicted point x˜ and training sample points (x1, x2, …, xn).With the predicted responsey^x˜, the prediction error of OKM can be evaluated by(11)φx˜=σ^21+uTx˜FTR−1F−1ux˜−rTx˜R−1rx˜s.t.ux˜=FTR−1rx˜−fx˜.Noticeably, by combining the nonlinear approximation ability of KM and the powerful global searching ability of the improved PSO algorithm, the presented OKM holds the potential to improve approximate efficiency and accuracy. The flowchart of OKM is drawn in Figure1.Figure 1
Flowchart of OKM. The establishment process of the OKM model is described.


### 2.1.2. Mathematical Modeling of DCOKM

Multiobjective design optimization of the multicomponent system refers to multiobjective (i.e., the force standard deviation in the first componentσF1, in the third component σF3, and in all rollers σMF, etc.) and multiconstraint (i.e., the sum of mean force in all rollers SMF and the sum of maximum force in all rollers SFmax, etc.), which is difficult to perform the nonlinear approximation with high efficiency and acceptable accuracy. To address the multiobjective and multiconstraint (MOMC) problem, we develop a decomposed collaborative optimized Kriging model (DCOKM) with respect to the decomposed collaborative strategy and optimized Kriging model, to further address the strong nonlinearity issues and enhance approximation efficacy. For the MOMC approximation with the presented DCOKM, an entire complex MOMC system with all input variables and output responses is decomposed into many single-objective single-constraint subsystems, in which each submodel has few input variables and output responses. Considering Latin hypercube sampling technique and OKM, the decomposed OKM for many single-objective single-constraints is established. Subsequently, considering collaborative sampling technique and decomposed OKM, the DCOKM is mathematically constructed. The mathematical modeling process of DCOKM is summarized.Assuming that the multiobjective and multiconstraint problem involvesr objectives and l constraints, when x(t) indicates the t-th input sample, then the corresponding t-th output response Y(t) is(12)Yt=fxt,xt=x1t,x2t,…,xntT,t=1,2,…,p,p+1,…,l,where t = 1, 2, …, p indicates output objective, while t = p + 1, p + 2, …, l represents the constraint response.With a given predicted pointx˜t for (10), the corresponding output response Y˜t is(13)Y˜t=fTx˜tβ^t+rTx˜tRt−1yt−Fβ^t.Equation (13) is the decomposed OKM under p-th objective response and l-th constraint response, in which β^t, R(t), and rTx˜t indicate the optimal regression coefficient, the correlation model, and the correlation vector in decomposed OKM, respectively. From the above analysis, the OKM (10) of multiobjective and multiconstraint models is divided into a series of OKM submodels under single-objective and single-constraint such as (12) and (13), and the approximation complexity is significantly reduced. As the efficient and accurate DCOKM is developed, the MOMC problem is promising to be addressed for multiobjective design optimization.


## 2.1.1. Optimized Kriging Model, OKM

By integrating the Kriging model (KM) with the best unbiasedness ability and the improved particle swarm optimization (IPSO) with global search capability [43], the optimized Kriging model (OKM) is introduced as(2)yx=fTx⋅β+zx,where x = [x1, x2, …, xd] indicates the input variable; d is the dimension number of input parameters; fT(x) β signifies the global approximation trend; fT(x) = [f1(x), f2(x), …, fm(x)] is the regression basis function; β = [β1, β2, …, βm] is the regression coefficient; m is the dimension number of regression function; and z(x) is the local random deviation, which follows Gaussian distribution (0, σ2). Herein, the covariance measure of z(x) is defined as(3)Covzxi,zxi=σ2Rxi,xj,where xi and xj denote the i-th and the j-th input variables, respectively; R(·,·) is the Gauss correlation function, which is(4)Rxi,xj=exp−∑k=1nθkxik−xjk2=∏k=1nexp−θkxik−xjk2,in which n represents the dimension number of design variables; x(k) i and x(k) j are the k-th member of variables xi and xj, respectively; θk is the k-th member of correlation parameter θ.Moreover, the optimal correlation parameterθ∗ of OKM can be obtained by the following maximum likelihood estimation, i.e.,(5)θ∗=argmaxψθ=arg min12lnR+mlnσ^2,where ψ(θ) indicates the likelihood function.Then, the corresponding regression coefficientβ and variance σ2 of OKM can be estimated by(6)β^=FTR−1F−1FTR−1y,σ^2=1ny−Fβ^TR−1y−Fβ^,where F is the unit column vector and R is the correlation matrix.Clearly, the approximation accuracy of OKM mostly depends on the correlation parameterθ. Once the optimal correlation parameter θ∗ is acquired, the optimal unbiased prediction of OKM is ensured. Therefore, the high-accurate surrogate modeling problem is transformed into the searching problem of the optimal solution θ∗.Searching algorithm is a key factor to find the optimal solutionθ∗, since its searching efficiency and accuracy directly influence the approximation efficacy of the Kriging model. The PSO algorithm is an important nongradient searching algorithm, and it holds good collaborative searching ability [26]. However, the fixed inertia weight and learning factors of the ordinary PSO algorithm often fall into blind searching and trap by local optimum easily, which greatly affects the searching efficiency and accuracy. To address this issue, we propose an improved PSO algorithm by designing adaptive inertia weight and asynchronous learning factors, with which the iteration steps change. The design is to acquire better optimal solution θ∗ by balancing the global searching ability and local searching ability. The basic thought of improved PSO algorithm is summarized as follows.We assume that the particle position is assigned as the correlation parameterθ and the fitness value adopts the ψ(θ) of the optimized Kriging model. Hereinto, each particle is a potential solution for the optimal correlation parameter θ∗. In searching process, with following the direction guidance of optimal particles, the particle swarm will search for the optimal solution by updating particle individual positions, individual extreme values, and population extremum values. The renewal formulas of the improved PSO algorithm are introduced as(7)Vik+1=ωVik+c1r1Pik−Xik+c2r2Pgk−XikXik+1=Xik+Vik+1i=1,2,⋯,N,in which the adaptive inertia weight ω and asynchronous learning factors c are(8)ω=ωmax−ωmax−ωminf−fminfavg−fmin,f≤favg,ωmax,f>favg,(9)c=c1=c1s+c1n−c1s×k/K,c2=c2s+c2n−c2s×k/K,where i represents the i-th particle; k is the current iteration number; K is the largest iteration number; ω is the inertia weight; Vi is the current particle velocity; Xi is the current particle position; Pi is the current individual extremum; Pg is the current population extremum; c1 and c2 are the individual and population leaning factors, respectively; r1 and r2 are the random numbers during time domain [0,1]; ωmax is the maximum inertia weight; ωmin is the minimum inertia weight; f is the current fitness value; favg is the average fitness value; fmin is the minimum fitness value; c1s and c2s are the initial individual learning factor and initial population learning factor, respectively; and c1n and c2n are the individual learning factor and population learning factor in the largest iterations, respectively.With finite iteration times of the improved PSO algorithm, the specified searching precision is met and the optimal correlation parametersθ∗ are acquired. In light of Kriging interpolation theory, the best unbiased prediction of OKM is achieved, and the predicted response y^x˜ related to the predicted point x˜ outside sample set can be expressed by(10)y^x˜=fTx˜β^+rTx˜R−1y−Fβ^,where r=Rx˜,x1,Rx˜,x2,…,Rx˜,xnTis the correlation vector between an predicted point x˜ and training sample points (x1, x2, …, xn).With the predicted responsey^x˜, the prediction error of OKM can be evaluated by(11)φx˜=σ^21+uTx˜FTR−1F−1ux˜−rTx˜R−1rx˜s.t.ux˜=FTR−1rx˜−fx˜.Noticeably, by combining the nonlinear approximation ability of KM and the powerful global searching ability of the improved PSO algorithm, the presented OKM holds the potential to improve approximate efficiency and accuracy. The flowchart of OKM is drawn in Figure1.Figure 1
Flowchart of OKM. The establishment process of the OKM model is described.


## 2.1.2. Mathematical Modeling of DCOKM

Multiobjective design optimization of the multicomponent system refers to multiobjective (i.e., the force standard deviation in the first componentσF1, in the third component σF3, and in all rollers σMF, etc.) and multiconstraint (i.e., the sum of mean force in all rollers SMF and the sum of maximum force in all rollers SFmax, etc.), which is difficult to perform the nonlinear approximation with high efficiency and acceptable accuracy. To address the multiobjective and multiconstraint (MOMC) problem, we develop a decomposed collaborative optimized Kriging model (DCOKM) with respect to the decomposed collaborative strategy and optimized Kriging model, to further address the strong nonlinearity issues and enhance approximation efficacy. For the MOMC approximation with the presented DCOKM, an entire complex MOMC system with all input variables and output responses is decomposed into many single-objective single-constraint subsystems, in which each submodel has few input variables and output responses. Considering Latin hypercube sampling technique and OKM, the decomposed OKM for many single-objective single-constraints is established. Subsequently, considering collaborative sampling technique and decomposed OKM, the DCOKM is mathematically constructed. The mathematical modeling process of DCOKM is summarized.Assuming that the multiobjective and multiconstraint problem involvesr objectives and l constraints, when x(t) indicates the t-th input sample, then the corresponding t-th output response Y(t) is(12)Yt=fxt,xt=x1t,x2t,…,xntT,t=1,2,…,p,p+1,…,l,where t = 1, 2, …, p indicates output objective, while t = p + 1, p + 2, …, l represents the constraint response.With a given predicted pointx˜t for (10), the corresponding output response Y˜t is(13)Y˜t=fTx˜tβ^t+rTx˜tRt−1yt−Fβ^t.Equation (13) is the decomposed OKM under p-th objective response and l-th constraint response, in which β^t, R(t), and rTx˜t indicate the optimal regression coefficient, the correlation model, and the correlation vector in decomposed OKM, respectively. From the above analysis, the OKM (10) of multiobjective and multiconstraint models is divided into a series of OKM submodels under single-objective and single-constraint such as (12) and (13), and the approximation complexity is significantly reduced. As the efficient and accurate DCOKM is developed, the MOMC problem is promising to be addressed for multiobjective design optimization.


## 2.2. Dynamic Multiobjective Genetic Algorithm, DMOGA

The resolving algorithm is a critical factor for solving the MODO model and obtaining the optimal solutions. NSGA-II is a vital nondominated fast sorting algorithm in multiobjective optimal design, which can effectively control the population distribution and reduce the resolving complexity. However, due to randomness and uncertain fluctuation effects in simulated binary crossover and polynomial mutation operator, its convergence precision is still unacceptable for the complex MODO model. To tackle with this problem, by designing the arithmetic crossover operator and Poisson mutation operator, we propose a dynamic multiobjective genetic algorithm (DMOGA) based on the traditional NSGA-II. The objective of the DMOGA is to reduce the uncertain fluctuation effects of operators and gain ideal Pareto optimal frontier. The design thought of the arithmetic crossover operator and Poisson mutation operator is presented as follows.To improve the population diversity and avoid over propagation of excellent solutions, an arithmetic crossover operator is presented by combining the information of nondominated sorting levels. The arithmetic crossover operator can keep the high-ranking individuals in parent generation and increase the low-ranking individuals in offspring population. Therefore, the quality of offspring population and the diversity of whole population shall be improved. For two given parent vectorsP1=v11,v21,⋯,vm1 and P2=v12,v22,⋯,vm2, two offspring Cz = (z1, z2, …, zm) and Cw=w1,w2,…,wm can be obtained by the designed arithmetic crossover operator:(14)zi=αivi1+1−αivi2,wi=αivi2+1−αivi1,where vi1and vi2 indicate the i-th parent individual of the first and second generations; zi and wi are the child individual; i = 1, 2, …, m; and αi is the random number in [0,1], which is expanded as(15)αi=BrAr+Br,Ar≠Br,AdAd+Bd,Ar=Br,where Ar and Br indicate the nondominated sorting level of individual A and individual B, respectively and Ad and Bd are the crowding distance of individual A and individual B, respectively.To change the disadvantage of low convergence speed caused by polynomial mutation parameters, a new Poisson distribution mutation operator is proposed, which can overstep the local optimum and thereby is conducive to search for global solutions. Assuming that the individual before mutation is expressed asx = (x1, x2, …, xn) and the individual after compilation is expressed as x′=x1′,x2′,⋯,xn′, then the Poisson mutation operator can be designed as(16)xi′=xi+xiu−xilδk,where xi′ means the i-th child individual after mutation; xi is the i-th parent individual; xui and xli are the upper and bottom boundaries of gene value at the mutation point i, respectively; i = 1, 2, …, m; and δk is the Poisson random number, which is introduced as(17)δk=λke−λk!,where λ denotes the distribution parameter and k is the variable value.To obtain the Pareto optimal frontier and acquire nondominated solutions, the DMOGA algorithm is proposed by designing the arithmetic crossover operator and Poisson mutation operator. The essential process of DMOGA algorithm is summarized as follows:(1)
Generate initialN individuals of population Pn, and set iteration times n = 0(2)
Evaluate fitness values of all individuals, and rank them with Pareto dominance and crowding distances(3)
Select individuals fromPn by the binary tournament method(4)
Obtain the childCn by the arithmetic crossover operator and Poisson mutation operator(5)
Generate populationQn by merging parent and child individuals, evaluate fitness value of all individuals, and perform fast nondominated sorting of Qn(6)
Calculate crowding degree and crowding distance, and select individuals to form the new species groupPn+1(7)
Terminate the algorithm when the accuracy requirement is met; otherwise, the algorithm will back to (2)(8)
Acquire Pareto optimal surface and global optimal solutionBy absorbing the global search capability of the arithmetic crossover operator and fast convergence ability of the Poisson mutation operator, the proposed DMOGA can accurately solve the multiobjective design optimization model and quickly obtain the Pareto solution set. To achieve the tradeoff effect between multiple objectives and optimal solution, a Utopia-Pareto directing adaptive (UPDA) search scheme is adopted by capturing plenty of characteristics and utilizing the ordering information of Pareto solutions [35, 44], which is conducive to find the best optimal solution for multiobjective design optimization of the multicomponent system.


## 2.3. Multiobjective Design Optimization Framework with DCOKM and DMOGA

For the multicomponent system optimal design problem, we develop a multiobjective design optimization framework based on a surrogate model (i.e., DCOKM) to approximate multiple responses and an optimal algorithm (i.e., DMOGA) to acquire optimal solution set. Herein, the improved PSO algorithm and Kriging model are first absorbed into decomposed collaborative strategy, to enhance the approximation efficacy of the surrogate model; then, the arithmetic crossover operator and Poisson mutation operator are designed, to acquire the efficient and accurate convergence of the optimal algorithm. Therefore, by combining the DCOKM and DMOGA, the calculation accuracy and efficiency of the multiobjective design optimization model is promising to be greatly improved. The multiobjective design optimization framework is shown in Figure2.Figure 2
Multiobjective design optimization framework.


## 3. Case Study

In this section, a carrier roller system from the track driving device of the excavator is selected as an engineering case to verify the effectiveness of the presented multiobjective design optimization framework. It should be noted that all computations are performed on an Inter(R) Core (TM) Desktop Computer (i7-9700K CPU 3.6 GHz and 16 GB RAM).


### 3.1. Multiobjective Optimal Modeling

A typical carrier roller system mainly includes three components, and each component endures multidirectional nonuniform loading. The change of structural size for one component will cause the change of the contact force between the whole track and the roller, and then influence the load of other components, which will lead to the load change among the multiple component systems. The load transfer path of the carrier roller system is shown in Figure3. During the working process, each carrier roller bears the load influence of three directions at the same time, and the load fluctuation in each direction is indicated in Figure 4. To confirm the important influencing factors of complex nonuniform loading, the load distribution of the carrier roller system is quantified in Figures 5–8. It reveals that the resultant force of three components shows V-shape distribution. Moreover, to further find the reasons of the great dispersion on resultant force, the load distribution of standard deviation in each direction is investigated and illustrated in Figures 9–11. Meanwhile, the stiffness analysis of the carrier roller system is performed and the results are depicted in Figure 12. We observe that with the increase of stiffness, the force standard deviation in the first component σF1, in the third component σF3, and in all rollers σMF (i.e., the sum of maximum force) appears significant fluctuations, while the sum of mean force in all rollers SMF and the sum of maximum force in all rollers SFmax show little variations.Figure 3
Load transfer path of the multicomponent system.Figure 4
Load analysis of the first carrier roller: (a) forward forceX, (b) radial force Y, (c) axial force Z, and (d) three directional forces.
(a)(b)(c)(d)Figure 5
Load distribution of forward forceX of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 6
Load distribution of radial forceY of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 7
Load distribution of axial forceZ of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 8
Load distribution of resultant force FM of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 9
Load distribution of standard deviationX of the carrier roller system: (a) simulation history. and (b) boxplot curve.
(a)(b)Figure 10
Load distribution of standard deviationY of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 11
Load distribution of standard deviationZ of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 12
Load fluctuation with different stiffness value. (a) Force standard deviations. (b) Sum of mean force and maximum force.
(a)(b)Therefore, by minimizing the responses ofσF1, σF3, and σMF, the nonuniform appearance would be reduced and the load distribution will become more uniform. Therefore, to minimize the nonuniform loading of the multicomponent system, by regarding the force standard deviations (i.e., σF1, σF3, and σMF) as optimal objectives, force mean values (i.e., SMF and SFmax) and other constraints as constraint functions, and the structural sizes (i.e., height of first component x1, height of third component x3, and displacement condition x2) as design variables, the multiobjective design optimization model is established as illustrated in the following equation:(18)Findx1,x2,x3,minσF1x=σF1x1,x2,x3,minσF3x=σF3x1,x2,x3,minσMFx=σMFx1,x2,x3,S.t.SMF<SMF,SFmax<SFmax,407.5≤x1≤527.5,10≤x2≤70,407.5≤x3≤527.5,where [SMF] and [SFmax] indicate the allowable sum of mean force and allowable sum of maximum force in all rollers, respectively. In this study, [SMF] and [SFmax] are set as 59013 N and 131842 N, respectively.


### 3.2. Surrogate Modeling by DCOKM

#### 3.2.1. Sample Generations

To obtain the objective responses and constraint function values of the multiobjective design optimization model, the kinematic and dynamic equations are solved as follows.Based on Lagrange multiplier modeling technique [45, 46], the dynamic equation of the carrier roller system in the generalized coordinate system is established as(19)ddt∂T∂q˙−∂T∂qT+fqTρ+gq˙Tμ−Q=0,where T denotes the system kinetic energy; q is the generalized coordinates, which is related to the variable x1, x2, and x3; q˙ is the generalized velocities, which is related to the variable x1, x2, and x3; Q is the generalized forces; ρ is the Lagrange multiplier array under holonomic constraints; and μ is the Lagrange multiplier array under nonholonomic constraints. Herein, the holonomic and nonholonomic constraints in (18) are expressed as(20)fq,t=0,gq,q˙,t=0.Subsequently, by solving the constraint equation with the gear prediction-correction algorithm, the velocity and acceleration at timetn can be obtained as(21)∂ϕ∂qq˙=−∂ϕ∂t,(22)∂ϕ∂qq¨=−∂2ϕ∂t2+∑k=1n∑l=1n∂2ϕ∂qk∂qlq˙kq˙l+∂∂t∂ϕ∂qq˙+∂∂q∂ϕ∂tq˙,where qk and ql are the generalized coordinates of the k-th and l-th iterations, respectively, and q˙k and q˙l are the generalized velocities of the k-th and l-th iterations, respectively.Consequently, by solving the Lagrange multiplier equation, the constraint and response reaction force in (18) are acquired:(23)∂ϕ∂qTλ=−ddt∂T∂q˙T+∂T∂qT+Q.Obviously, it is time-consuming to solving the nonlinear dynamic equations for thousands of times. Therefore, to build the surrogate model of the above dynamic equations, a handful of the input variables (x1, x2, and x3) are extracted by Latin hypercube sampling, and the corresponding output responses (σF1, σF3, σMF, SMF, and SFmax) are obtained by running the above dynamic equations. The extracted samples would be considered as training samples for surrogate modeling, which are shown in Table 1.Table 1
Extracted training samples.
Samplex1x2x3σF1σF3σMFSMFSFmax1390394465478.366088.8662096.7534708.4870669.7324331064626752.666670.624389.7849454.3293768.563409484064213.394146.5394882.8434017.6663289.524418824535718.215777.5631403.2743320.4485526.515393864596297.806242.5552301.8541561.7789275.856449914485431.015322.8615860.5947852.889216.5674241404518671.587605.7782878.4851373.97113861.8084731084887110.486871.81312331.3657903.74112841.709485764706382.036362.14611143.5952054.4593507.65104571284094517.625344.15310078.7449212.41112440.5111451234324282.405153.784734.3037107.0376721.8612423554436500.605537.3991825.1939707.0295494.29134811234003658.484148.8013708.9551039.6199463.7514458984033675.753830.379571.4544760.4185159.2315441604856235.715748.826150.3446056.87111631164791364508958.538488.7612654.4858537.9011964317472313952550.093079.388237.6736571.7265595.7518462694646185.576892.617400.7747623.44105196.8019429904345715.606410.483553.8043720.3985185.6320426474685295.955673.302924.8240827.1577965.02214141073922459.883156.218467.2639682.1178960.6822413304524183.204876.16520.43036114.3663734.8623456563974436.304444.577942.6538563.3583004.92244281034876741.817158.786310.2151060.0795940.3025477244654936.896018.777260.1041893.8772833.50263991294729084.016987.263386.8249349.60117109.3027475264404903.206094.306697.0539920.4567536.7928430974104272.045065.086294.4942083.3287140.5429431573983964.194139.246453.21236342.5070549.1930480634045216.056164.2210156.8742814.12120803.00


#### 3.2.2. Fitting DCOKM

By importing the generated samples into the presented DCOKM method, we perform the surrogate modeling for multiobjective and multiconstraint. Considering the improved PSO algorithm, the optimal correlation parameterθ∗ is searched and substituted into OKM, and the DCOKM can be established. The optimal correlation parameter θ∗ of DCOKM is shown in Table 2. When the population size N = 100, number of iteration K = 80, maximum inertia weight ωmax = 0.9, minimum inertia weight ωmin = 0.4, initial individual learning factor c1s = 2.5, individual learning factor in the largest iterations c1n = 0.5, initial population learning factor c2s = 0.5, and population learning factor in the largest iterations c2n = 2.5, the training process of DCOKM is shown in Figure 13, which shows that the DCOKM modeling holds fast convergence speed and accuracy. Moreover, to verify the approximation accuracy of the established surrogate model, the absolute errors between the predicted objective/constraint and the actual objective/constraint are calculated as shown in Figures 14 and 15. It reveals that the absolute errors of the built DCOKM are nearly equal to zero, which is sufficient to meet the accuracy requirements of surrogate modeling.Table 2
Optimal correlation parameter of DCOKM.
ResponseRelated parametersθ∗θ1θ2θ3σF138.8348840.3881942.71815σF33.0864333.2010693.394138σMֿF48.5207350.4612853.37277SMF49.9766351.9756954.97429SFmax49.9893451.9889154.98827Figure 13
Training process of DCOKM.Figure 14
Absolute errors of three objectives: (a)σF1., (b) σF3., and (c) σMֿF.
(a)(b)(c)Figure 15
Absolute errors of two constraints: (a)SMF and (b) SFmax.
(a)(b)


### 3.3. Optimal Solution by DMOGA

By employing the built DCOKM to approximate the multiobjective and multiconstraint, the DMOGA is executed to solve the established MODO model shown in (18). Herein, the parent population size, crossover probability, mutation probability, and iteration number of DMOGA algorithm are set as 100, 0.9, 0.1, and 1000, respectively. After 1000 iterative times, the Pareto optimal frontier is obtained as shown in Figure 16. According to the Pareto curve in Figure 16, the optimization objectives σMF and σF1, and σMF and σF3 both show negative correlation and σF1 and σF3 show positive correlation. To achieve the reasonable balance between three objectives, we sort σF1 and σF3 as ascending order and σMF as descending order. Then, the middle point in the Pareto front curve is selected as the optimal solution of multiobjective optimization [35, 44]. In this case, the design variables x1, x2, and x3 corresponding to the optimal solution are chosen as 412, 10, and 459, respectively.Figure 16
Pareto optimal frontier of MODO.To validate the optimal results, we compare the optimal parameters before and after optimization in Table3, which shows that all objective responses are reduced after optimization, particularly the σMF decreased by 92%. To check the loading uniformity after optimization, we further compare the loading distribution before and after optimization in Figures 17 and 18. The loading distribution results indicate that the system loads are significantly reduced and the loading after optimization becomes more uniform than that of before optimization. Moreover, through the multiobjective design optimization of the carrier roller system, the average failure ratio of track driving device is reduced from 31% to 8%.Table 3
Comparison of results before and after optimization.
Optimization indexBefore optimizationAfter optimizationReduce the proportion (%)Remarkx1467412——x24010——x3467459——σF18735290467Reduce track run out and impact loadσF38156438246Reduce track run out and impact loadσMF1098792492Better load uniformitySMF590133392455Reduce wear failure rateSFmax1318425747743Reduce fracture failure rateFigure 17
Simulation history of three carrier rollers before and after optimization: (a) before optimization and (b) after optimization.
(a)(b)Figure 18
Resultant force of three carrier rollers before and after optimization: (a) before optimization and (b) after optimization.
(a)(b)


### 3.4. Validations of Methods

Because of the particularity of the study case (practical engineering problem), there is no standard solution for reference. To compare the effectiveness of the proposed method, some advanced surrogate models and optimization algorithms are selected to perform the multiobjective design of the carrier roller system and then compared with the presented model (DCOKM) and algorithm (DMOGA). Herein, to support the feasibility of the proposed surrogate model (i.e., DCOKM), the multiobjective and multiconstraint approximation of the carrier roller system MODO is also studied based on the modeling methods of the direct Monte Carlo simulation (MCS), response surface (RS) model, KM, OKM, and the presented DCOKM; Furtherly, to validate the effectiveness of the developed optimal algorithm (i.e., DMOGA), the multiobjective design optimization of the carrier roller system is also performed based on NCGA algorithm, AMGA algorithm, NSGA-II algorithm, and DMOGA algorithm. Note that to ensure the rationality of algorithm comparison, all of optimal algorithms are using DCOKM to approximate the objectives and constraints in optimal algorithm comparisons. The approximation efficacy of different models is listed in Table4, and the computing efficiency and accuracy of different algorithms are listed in Table 5.Table 4
Approximation efficacy of different models.
MethodsFitting surrogate modelsTesting surrogate modelsPrecisionFitting numberFitting time (s)Testing numberTesting time (s)MCS————1RS2452.16 × 10720015.6230.9451KM1801.55 × 10720013.5320.9792OKM1301.12 × 10720011.3670.9875DCOKM907.78 × 1062009.24760.9952Table 5
Computing efficiency and accuracy of different algorithms.
MethodsBefore optimization (×104 N)After optimization (×104 N)ReductionTime (s)Y1Y2Y3Y1Y2Y3Y1Y2Y3NCGA0.87350.81561.09870.33320.46200.10230.610.430.9018.23AMGA0.87350.81561.09870.31460.46570.12750.640.430.8816.52NSGA-II0.87350.81561.09870.30650.46400.11120.650.430.8915.72DMOGA0.87350.81561.09870.29040.43820.09240.670.460.9213.31Note. Y1 indicates the σF1, Y2 indicates the σF3, and Y3 indicates the σMF.As illustrated in Table4, the fitting number and fitting time of DCOKM are less than RS, KM, OKM, and the testing time of KM, OKM, and DCOKM are far less than the MCS method. Therefore, the DCOKM holds the highest approximation efficiency. Moreover, the DCOKM is more accurate than RS, KM, and OKM and is almost consistent with the direct MCS. The superior performance of DCOKM is induced by (i) the proposed improved PSO training algorithm can rapidly obtain the optimal parameters of the optimized Kriging model, which improves the nonlinear fitting ability and generalization ability of DCOKM; (ii) the decomposed collaborative sampling strategy in DCOKM brings the distributed parallel computation, which is conducive to reduce nonlinear features and computational task. By integrating the strengths of the optimized Kriging model and decomposed collaborative strategy, DCOKM can efficiently and accurately approximate the multiresponse and multiconstraint. Therefore, the DCOKM is a high-accuracy and high-efficiency surrogate model in multiobjective design optimization.As revealed in Table5, the optimization time of the DMOGA algorithm is shorter than NSGA-II, AMGA, and NCGA, and the optimization efficacy of DMOGA is better than NSGA-II, AMGA, and NCGA. The computing advantages of DMOGA is induced by (i) the proposed arithmetic crossover operator can combine nondominated sorting level information, which improves the quality of the offspring population and ensure the diversity of the whole population; (ii) the Poisson mutation operator can generate random numbers in a wider range and overstep the local optimum efficiently, which ensures a faster convergence speed to search for global solutions. Therefore, the DMOGA algorithm is a high-efficiency algorithm for resolving the MODOTherefore, with the integration of the surrogate model (i.e., DCOKM) and optimal algorithm (i.e., DMOGA), the proposed multiobjective design optimization framework is validated to improve optimal efficiency while maintaining acceptable optimal accuracy and thereby is an effective way for the multiobjective optimization design of the multicomponent system.


## 3.1. Multiobjective Optimal Modeling

A typical carrier roller system mainly includes three components, and each component endures multidirectional nonuniform loading. The change of structural size for one component will cause the change of the contact force between the whole track and the roller, and then influence the load of other components, which will lead to the load change among the multiple component systems. The load transfer path of the carrier roller system is shown in Figure3. During the working process, each carrier roller bears the load influence of three directions at the same time, and the load fluctuation in each direction is indicated in Figure 4. To confirm the important influencing factors of complex nonuniform loading, the load distribution of the carrier roller system is quantified in Figures 5–8. It reveals that the resultant force of three components shows V-shape distribution. Moreover, to further find the reasons of the great dispersion on resultant force, the load distribution of standard deviation in each direction is investigated and illustrated in Figures 9–11. Meanwhile, the stiffness analysis of the carrier roller system is performed and the results are depicted in Figure 12. We observe that with the increase of stiffness, the force standard deviation in the first component σF1, in the third component σF3, and in all rollers σMF (i.e., the sum of maximum force) appears significant fluctuations, while the sum of mean force in all rollers SMF and the sum of maximum force in all rollers SFmax show little variations.Figure 3
Load transfer path of the multicomponent system.Figure 4
Load analysis of the first carrier roller: (a) forward forceX, (b) radial force Y, (c) axial force Z, and (d) three directional forces.
(a)(b)(c)(d)Figure 5
Load distribution of forward forceX of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 6
Load distribution of radial forceY of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 7
Load distribution of axial forceZ of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 8
Load distribution of resultant force FM of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 9
Load distribution of standard deviationX of the carrier roller system: (a) simulation history. and (b) boxplot curve.
(a)(b)Figure 10
Load distribution of standard deviationY of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 11
Load distribution of standard deviationZ of the carrier roller system: (a) simulation history and (b) boxplot curve.
(a)(b)Figure 12
Load fluctuation with different stiffness value. (a) Force standard deviations. (b) Sum of mean force and maximum force.
(a)(b)Therefore, by minimizing the responses ofσF1, σF3, and σMF, the nonuniform appearance would be reduced and the load distribution will become more uniform. Therefore, to minimize the nonuniform loading of the multicomponent system, by regarding the force standard deviations (i.e., σF1, σF3, and σMF) as optimal objectives, force mean values (i.e., SMF and SFmax) and other constraints as constraint functions, and the structural sizes (i.e., height of first component x1, height of third component x3, and displacement condition x2) as design variables, the multiobjective design optimization model is established as illustrated in the following equation:(18)Findx1,x2,x3,minσF1x=σF1x1,x2,x3,minσF3x=σF3x1,x2,x3,minσMFx=σMFx1,x2,x3,S.t.SMF<SMF,SFmax<SFmax,407.5≤x1≤527.5,10≤x2≤70,407.5≤x3≤527.5,where [SMF] and [SFmax] indicate the allowable sum of mean force and allowable sum of maximum force in all rollers, respectively. In this study, [SMF] and [SFmax] are set as 59013 N and 131842 N, respectively.


## 3.2. Surrogate Modeling by DCOKM

### 3.2.1. Sample Generations

To obtain the objective responses and constraint function values of the multiobjective design optimization model, the kinematic and dynamic equations are solved as follows.Based on Lagrange multiplier modeling technique [45, 46], the dynamic equation of the carrier roller system in the generalized coordinate system is established as(19)ddt∂T∂q˙−∂T∂qT+fqTρ+gq˙Tμ−Q=0,where T denotes the system kinetic energy; q is the generalized coordinates, which is related to the variable x1, x2, and x3; q˙ is the generalized velocities, which is related to the variable x1, x2, and x3; Q is the generalized forces; ρ is the Lagrange multiplier array under holonomic constraints; and μ is the Lagrange multiplier array under nonholonomic constraints. Herein, the holonomic and nonholonomic constraints in (18) are expressed as(20)fq,t=0,gq,q˙,t=0.Subsequently, by solving the constraint equation with the gear prediction-correction algorithm, the velocity and acceleration at timetn can be obtained as(21)∂ϕ∂qq˙=−∂ϕ∂t,(22)∂ϕ∂qq¨=−∂2ϕ∂t2+∑k=1n∑l=1n∂2ϕ∂qk∂qlq˙kq˙l+∂∂t∂ϕ∂qq˙+∂∂q∂ϕ∂tq˙,where qk and ql are the generalized coordinates of the k-th and l-th iterations, respectively, and q˙k and q˙l are the generalized velocities of the k-th and l-th iterations, respectively.Consequently, by solving the Lagrange multiplier equation, the constraint and response reaction force in (18) are acquired:(23)∂ϕ∂qTλ=−ddt∂T∂q˙T+∂T∂qT+Q.Obviously, it is time-consuming to solving the nonlinear dynamic equations for thousands of times. Therefore, to build the surrogate model of the above dynamic equations, a handful of the input variables (x1, x2, and x3) are extracted by Latin hypercube sampling, and the corresponding output responses (σF1, σF3, σMF, SMF, and SFmax) are obtained by running the above dynamic equations. The extracted samples would be considered as training samples for surrogate modeling, which are shown in Table 1.Table 1
Extracted training samples.
Samplex1x2x3σF1σF3σMFSMFSFmax1390394465478.366088.8662096.7534708.4870669.7324331064626752.666670.624389.7849454.3293768.563409484064213.394146.5394882.8434017.6663289.524418824535718.215777.5631403.2743320.4485526.515393864596297.806242.5552301.8541561.7789275.856449914485431.015322.8615860.5947852.889216.5674241404518671.587605.7782878.4851373.97113861.8084731084887110.486871.81312331.3657903.74112841.709485764706382.036362.14611143.5952054.4593507.65104571284094517.625344.15310078.7449212.41112440.5111451234324282.405153.784734.3037107.0376721.8612423554436500.605537.3991825.1939707.0295494.29134811234003658.484148.8013708.9551039.6199463.7514458984033675.753830.379571.4544760.4185159.2315441604856235.715748.826150.3446056.87111631164791364508958.538488.7612654.4858537.9011964317472313952550.093079.388237.6736571.7265595.7518462694646185.576892.617400.7747623.44105196.8019429904345715.606410.483553.8043720.3985185.6320426474685295.955673.302924.8240827.1577965.02214141073922459.883156.218467.2639682.1178960.6822413304524183.204876.16520.43036114.3663734.8623456563974436.304444.577942.6538563.3583004.92244281034876741.817158.786310.2151060.0795940.3025477244654936.896018.777260.1041893.8772833.50263991294729084.016987.263386.8249349.60117109.3027475264404903.206094.306697.0539920.4567536.7928430974104272.045065.086294.4942083.3287140.5429431573983964.194139.246453.21236342.5070549.1930480634045216.056164.2210156.8742814.12120803.00


### 3.2.2. Fitting DCOKM

By importing the generated samples into the presented DCOKM method, we perform the surrogate modeling for multiobjective and multiconstraint. Considering the improved PSO algorithm, the optimal correlation parameterθ∗ is searched and substituted into OKM, and the DCOKM can be established. The optimal correlation parameter θ∗ of DCOKM is shown in Table 2. When the population size N = 100, number of iteration K = 80, maximum inertia weight ωmax = 0.9, minimum inertia weight ωmin = 0.4, initial individual learning factor c1s = 2.5, individual learning factor in the largest iterations c1n = 0.5, initial population learning factor c2s = 0.5, and population learning factor in the largest iterations c2n = 2.5, the training process of DCOKM is shown in Figure 13, which shows that the DCOKM modeling holds fast convergence speed and accuracy. Moreover, to verify the approximation accuracy of the established surrogate model, the absolute errors between the predicted objective/constraint and the actual objective/constraint are calculated as shown in Figures 14 and 15. It reveals that the absolute errors of the built DCOKM are nearly equal to zero, which is sufficient to meet the accuracy requirements of surrogate modeling.Table 2
Optimal correlation parameter of DCOKM.
ResponseRelated parametersθ∗θ1θ2θ3σF138.8348840.3881942.71815σF33.0864333.2010693.394138σMֿF48.5207350.4612853.37277SMF49.9766351.9756954.97429SFmax49.9893451.9889154.98827Figure 13
Training process of DCOKM.Figure 14
Absolute errors of three objectives: (a)σF1., (b) σF3., and (c) σMֿF.
(a)(b)(c)Figure 15
Absolute errors of two constraints: (a)SMF and (b) SFmax.
(a)(b)


## 3.2.1. Sample Generations

To obtain the objective responses and constraint function values of the multiobjective design optimization model, the kinematic and dynamic equations are solved as follows.Based on Lagrange multiplier modeling technique [45, 46], the dynamic equation of the carrier roller system in the generalized coordinate system is established as(19)ddt∂T∂q˙−∂T∂qT+fqTρ+gq˙Tμ−Q=0,where T denotes the system kinetic energy; q is the generalized coordinates, which is related to the variable x1, x2, and x3; q˙ is the generalized velocities, which is related to the variable x1, x2, and x3; Q is the generalized forces; ρ is the Lagrange multiplier array under holonomic constraints; and μ is the Lagrange multiplier array under nonholonomic constraints. Herein, the holonomic and nonholonomic constraints in (18) are expressed as(20)fq,t=0,gq,q˙,t=0.Subsequently, by solving the constraint equation with the gear prediction-correction algorithm, the velocity and acceleration at timetn can be obtained as(21)∂ϕ∂qq˙=−∂ϕ∂t,(22)∂ϕ∂qq¨=−∂2ϕ∂t2+∑k=1n∑l=1n∂2ϕ∂qk∂qlq˙kq˙l+∂∂t∂ϕ∂qq˙+∂∂q∂ϕ∂tq˙,where qk and ql are the generalized coordinates of the k-th and l-th iterations, respectively, and q˙k and q˙l are the generalized velocities of the k-th and l-th iterations, respectively.Consequently, by solving the Lagrange multiplier equation, the constraint and response reaction force in (18) are acquired:(23)∂ϕ∂qTλ=−ddt∂T∂q˙T+∂T∂qT+Q.Obviously, it is time-consuming to solving the nonlinear dynamic equations for thousands of times. Therefore, to build the surrogate model of the above dynamic equations, a handful of the input variables (x1, x2, and x3) are extracted by Latin hypercube sampling, and the corresponding output responses (σF1, σF3, σMF, SMF, and SFmax) are obtained by running the above dynamic equations. The extracted samples would be considered as training samples for surrogate modeling, which are shown in Table 1.Table 1
Extracted training samples.
Samplex1x2x3σF1σF3σMFSMFSFmax1390394465478.366088.8662096.7534708.4870669.7324331064626752.666670.624389.7849454.3293768.563409484064213.394146.5394882.8434017.6663289.524418824535718.215777.5631403.2743320.4485526.515393864596297.806242.5552301.8541561.7789275.856449914485431.015322.8615860.5947852.889216.5674241404518671.587605.7782878.4851373.97113861.8084731084887110.486871.81312331.3657903.74112841.709485764706382.036362.14611143.5952054.4593507.65104571284094517.625344.15310078.7449212.41112440.5111451234324282.405153.784734.3037107.0376721.8612423554436500.605537.3991825.1939707.0295494.29134811234003658.484148.8013708.9551039.6199463.7514458984033675.753830.379571.4544760.4185159.2315441604856235.715748.826150.3446056.87111631164791364508958.538488.7612654.4858537.9011964317472313952550.093079.388237.6736571.7265595.7518462694646185.576892.617400.7747623.44105196.8019429904345715.606410.483553.8043720.3985185.6320426474685295.955673.302924.8240827.1577965.02214141073922459.883156.218467.2639682.1178960.6822413304524183.204876.16520.43036114.3663734.8623456563974436.304444.577942.6538563.3583004.92244281034876741.817158.786310.2151060.0795940.3025477244654936.896018.777260.1041893.8772833.50263991294729084.016987.263386.8249349.60117109.3027475264404903.206094.306697.0539920.4567536.7928430974104272.045065.086294.4942083.3287140.5429431573983964.194139.246453.21236342.5070549.1930480634045216.056164.2210156.8742814.12120803.00


## 3.2.2. Fitting DCOKM

By importing the generated samples into the presented DCOKM method, we perform the surrogate modeling for multiobjective and multiconstraint. Considering the improved PSO algorithm, the optimal correlation parameterθ∗ is searched and substituted into OKM, and the DCOKM can be established. The optimal correlation parameter θ∗ of DCOKM is shown in Table 2. When the population size N = 100, number of iteration K = 80, maximum inertia weight ωmax = 0.9, minimum inertia weight ωmin = 0.4, initial individual learning factor c1s = 2.5, individual learning factor in the largest iterations c1n = 0.5, initial population learning factor c2s = 0.5, and population learning factor in the largest iterations c2n = 2.5, the training process of DCOKM is shown in Figure 13, which shows that the DCOKM modeling holds fast convergence speed and accuracy. Moreover, to verify the approximation accuracy of the established surrogate model, the absolute errors between the predicted objective/constraint and the actual objective/constraint are calculated as shown in Figures 14 and 15. It reveals that the absolute errors of the built DCOKM are nearly equal to zero, which is sufficient to meet the accuracy requirements of surrogate modeling.Table 2
Optimal correlation parameter of DCOKM.
ResponseRelated parametersθ∗θ1θ2θ3σF138.8348840.3881942.71815σF33.0864333.2010693.394138σMֿF48.5207350.4612853.37277SMF49.9766351.9756954.97429SFmax49.9893451.9889154.98827Figure 13
Training process of DCOKM.Figure 14
Absolute errors of three objectives: (a)σF1., (b) σF3., and (c) σMֿF.
(a)(b)(c)Figure 15
Absolute errors of two constraints: (a)SMF and (b) SFmax.
(a)(b)


## 3.3. Optimal Solution by DMOGA

By employing the built DCOKM to approximate the multiobjective and multiconstraint, the DMOGA is executed to solve the established MODO model shown in (18). Herein, the parent population size, crossover probability, mutation probability, and iteration number of DMOGA algorithm are set as 100, 0.9, 0.1, and 1000, respectively. After 1000 iterative times, the Pareto optimal frontier is obtained as shown in Figure 16. According to the Pareto curve in Figure 16, the optimization objectives σMF and σF1, and σMF and σF3 both show negative correlation and σF1 and σF3 show positive correlation. To achieve the reasonable balance between three objectives, we sort σF1 and σF3 as ascending order and σMF as descending order. Then, the middle point in the Pareto front curve is selected as the optimal solution of multiobjective optimization [35, 44]. In this case, the design variables x1, x2, and x3 corresponding to the optimal solution are chosen as 412, 10, and 459, respectively.Figure 16
Pareto optimal frontier of MODO.To validate the optimal results, we compare the optimal parameters before and after optimization in Table3, which shows that all objective responses are reduced after optimization, particularly the σMF decreased by 92%. To check the loading uniformity after optimization, we further compare the loading distribution before and after optimization in Figures 17 and 18. The loading distribution results indicate that the system loads are significantly reduced and the loading after optimization becomes more uniform than that of before optimization. Moreover, through the multiobjective design optimization of the carrier roller system, the average failure ratio of track driving device is reduced from 31% to 8%.Table 3
Comparison of results before and after optimization.
Optimization indexBefore optimizationAfter optimizationReduce the proportion (%)Remarkx1467412——x24010——x3467459——σF18735290467Reduce track run out and impact loadσF38156438246Reduce track run out and impact loadσMF1098792492Better load uniformitySMF590133392455Reduce wear failure rateSFmax1318425747743Reduce fracture failure rateFigure 17
Simulation history of three carrier rollers before and after optimization: (a) before optimization and (b) after optimization.
(a)(b)Figure 18
Resultant force of three carrier rollers before and after optimization: (a) before optimization and (b) after optimization.
(a)(b)


## 3.4. Validations of Methods

Because of the particularity of the study case (practical engineering problem), there is no standard solution for reference. To compare the effectiveness of the proposed method, some advanced surrogate models and optimization algorithms are selected to perform the multiobjective design of the carrier roller system and then compared with the presented model (DCOKM) and algorithm (DMOGA). Herein, to support the feasibility of the proposed surrogate model (i.e., DCOKM), the multiobjective and multiconstraint approximation of the carrier roller system MODO is also studied based on the modeling methods of the direct Monte Carlo simulation (MCS), response surface (RS) model, KM, OKM, and the presented DCOKM; Furtherly, to validate the effectiveness of the developed optimal algorithm (i.e., DMOGA), the multiobjective design optimization of the carrier roller system is also performed based on NCGA algorithm, AMGA algorithm, NSGA-II algorithm, and DMOGA algorithm. Note that to ensure the rationality of algorithm comparison, all of optimal algorithms are using DCOKM to approximate the objectives and constraints in optimal algorithm comparisons. The approximation efficacy of different models is listed in Table4, and the computing efficiency and accuracy of different algorithms are listed in Table 5.Table 4
Approximation efficacy of different models.
MethodsFitting surrogate modelsTesting surrogate modelsPrecisionFitting numberFitting time (s)Testing numberTesting time (s)MCS————1RS2452.16 × 10720015.6230.9451KM1801.55 × 10720013.5320.9792OKM1301.12 × 10720011.3670.9875DCOKM907.78 × 1062009.24760.9952Table 5
Computing efficiency and accuracy of different algorithms.
MethodsBefore optimization (×104 N)After optimization (×104 N)ReductionTime (s)Y1Y2Y3Y1Y2Y3Y1Y2Y3NCGA0.87350.81561.09870.33320.46200.10230.610.430.9018.23AMGA0.87350.81561.09870.31460.46570.12750.640.430.8816.52NSGA-II0.87350.81561.09870.30650.46400.11120.650.430.8915.72DMOGA0.87350.81561.09870.29040.43820.09240.670.460.9213.31Note. Y1 indicates the σF1, Y2 indicates the σF3, and Y3 indicates the σMF.As illustrated in Table4, the fitting number and fitting time of DCOKM are less than RS, KM, OKM, and the testing time of KM, OKM, and DCOKM are far less than the MCS method. Therefore, the DCOKM holds the highest approximation efficiency. Moreover, the DCOKM is more accurate than RS, KM, and OKM and is almost consistent with the direct MCS. The superior performance of DCOKM is induced by (i) the proposed improved PSO training algorithm can rapidly obtain the optimal parameters of the optimized Kriging model, which improves the nonlinear fitting ability and generalization ability of DCOKM; (ii) the decomposed collaborative sampling strategy in DCOKM brings the distributed parallel computation, which is conducive to reduce nonlinear features and computational task. By integrating the strengths of the optimized Kriging model and decomposed collaborative strategy, DCOKM can efficiently and accurately approximate the multiresponse and multiconstraint. Therefore, the DCOKM is a high-accuracy and high-efficiency surrogate model in multiobjective design optimization.As revealed in Table5, the optimization time of the DMOGA algorithm is shorter than NSGA-II, AMGA, and NCGA, and the optimization efficacy of DMOGA is better than NSGA-II, AMGA, and NCGA. The computing advantages of DMOGA is induced by (i) the proposed arithmetic crossover operator can combine nondominated sorting level information, which improves the quality of the offspring population and ensure the diversity of the whole population; (ii) the Poisson mutation operator can generate random numbers in a wider range and overstep the local optimum efficiently, which ensures a faster convergence speed to search for global solutions. Therefore, the DMOGA algorithm is a high-efficiency algorithm for resolving the MODOTherefore, with the integration of the surrogate model (i.e., DCOKM) and optimal algorithm (i.e., DMOGA), the proposed multiobjective design optimization framework is validated to improve optimal efficiency while maintaining acceptable optimal accuracy and thereby is an effective way for the multiobjective optimization design of the multicomponent system.


## 4. Conclusions

The purpose of this study is to develop an efficiently and accurately multiobjective optimization design framework for the multicomponent system. Firstly, to efficiently approximate the multiobjective and multiconstraint, a surrogate model named as the decomposed collaborative optimized Kriging model (DCOKM) is proposed, by absorbing the strengths of the improved PSO algorithm and Kriging model into decomposed collaborative strategy. Then, to accurately solve the multiobjective design optimization model, an optimal algorithm named as the dynamic multiobjective genetic algorithm (DMOGA) is developed, by designing the arithmetic crossover operator and Poisson mutation operator. Regarding the carrier roller system MODO as a study case, the validity and feasibility of multiobjective optimization design framework is proved. Some conclusions are summarized as follows:(1)
From the multiobjective design optimization of the carrier roller system, we discover that the total standard deviation of three carrier rollers is reduced by 92 %, and loading distribution is more uniform after optimization(2)
DCOKM has higher accuracy and efficiency than the ordinary Kriging model and optimized Kriging model and is an effective surrogate model for the multiobjective design optimization(3)
The optimal solutions of DMOGA are better than the NSGA-II algorithm due to the superior properties of the arithmetic crossover operator and Poisson mutation operator(4)
The comparison of methods reveals that the proposed multiobjective optimization design framework (i.e., DCOKM and DMOGA) possesses high-efficiency and high-accuracy in multiobjective optimal design of the multicomponent system



---

*Source: 7695419-2020-08-19.xml*