# Multi-Objective Optimal Design of Dropping Shock of Series Cushioning Packaging System

Authors: Yang Xue; Wei-Sheng Song; Hong-Tao Miao; Jing Wang; Yu Li
Journal: Shock and Vibration
 (2022)

---

## Abstract

The product bears the shock load during the circulation process. The cushioning materials need to absorb the energy generated during the shock to prevent damage to the product. The cushion packaging system generally consists of products, cushioning materials, and an outside packaging carton. Therefore, it is necessary to study the factors affecting the cushion performance of the system. This paper adopts a BP neural network to develop a more precise constitutive relationship. Based on the trained constitutive network model, the dynamic equation of the system is established, and the Runge-Kutta algorithm is utilized to resolve it. The multi-objective optimization (MOP) model of the series cushioning system is established for the first time. The model considers more factors, such as the stability of the cushioning material and whether to use local cushioning. The factors affecting the performance of the cushioning system are deeply studied, and a variety of evolutionary algorithms (PESA2, SPEA2, MOPSO, MOEA/D, and NSGA-II) are utilized to solve the multi-objective optimization problem. The performances of the solution sets obtained by the evolutionary algorithms are studied based on several performance indicators. The results show that these solution sets have good convergence, and the all-around performance of NSGA-II ranks first. Finally, comparing the knee point of the non-dominated solution set calculated by NSGA-II with the parameters in the literature, it is found that the former has better defense performance and less material consumption, which can be used for future product design. Therefore, the evolutionary algorithm can favorably resolve the multi-objective optimization of the series cushioning packaging system under the dropping shock.


---

## Body

## 1. Introduction

In circulation and transportation, the product is packaged to form a package to protect the product [1]. The package usually consists of the product, cushioning material, and outer packaging box. In addition to fixing the product during transportation, the primary role of cushioning material is to reduce the damage to the product caused by external shocks [2]. Most of the outer packaging boxes used in the market are corrugated paperboard boxes [3–6] and honeycomb paperboard boxes [7] used to protect products and facilitate transportation. In order to study whether the transportation package protection design to be adopted can effectively protect the product, it is necessary to analyze the response law of packaging under external force. In logistics and transportation, most of the damage to the packages is caused by vibration [8, 9] and shock [10]. Therefore, the response analysis of the package is studied under the dropping shock load in this paper.Three relationships need to be known for analyzing the dropping shock response of a package [11]: the constitutive model of cushioning materials (i.e., the relationship between stress, strain, and strain rate); initial conditions or boundary conditions; dynamic equation of package system. Therefore, establishing a constitutive model, that is, the realization of data fitting between stress, strain, and strain rate, is necessary for package response analysis.In general, the steps for establishing a constitutive model of cushioning materials are as follows [12]. A simplified mathematical model is suggested based on experimental data. Then an available form of constitutive relationship is offered from the theoretical level. Moreover, the constitutive parameters are determined based on experiments or optimization methods. Finally, the concrete function expression of the constitutive model is established. The modeling strategy can make constructing the model clear and straightforward and ensure that the error between the constitutive equation and the experimental data is within a specific range.Many scholars currently pay more attention to data fitting based on neural networks [13–15]. BP neural network is a type of neural network that is widely used [16–18]. It uses an error backpropagation learning algorithm composed of an input layer, several hidden layers, and an output layer. It enhances the fitting accuracy of the model through error backpropagation. In addition, it does not need to know the association between input and output in advance. Thus, it shifts the research focus of the fitting problems from the nonlinear characteristics of materials to constructing the correlation between input and output data.Therefore, this article will establish the neural network model reflecting the dynamic constitutive relationship of the cushioning material based on the BP neural network.Moreover, it needs the benefit of evolutionary algorithms (MOEAs) for deeply analyzing the influence of multiple parameters in the series cushioning system on the buffer performance. MOEAs have received extensive attention from researchers since their conception [19–21]. They have demonstrated powerful problem-solving capabilities in many fields, such as function optimization [22], production scheduling [23], logistics distribution [24], and so on. We have used the population evolutionary algorithms to effectively solve the mobile robot path planning problems [25, 26].Maintaining the convergence and diversity of the population is always an essential issue in the study of population evolutionary algorithms.Multi-population co-evolutionary strategy is an effective means [27]. The sub-populations evolve independently to maintain their particular characteristics [28]. The cooperative operator is designed to enhance the interaction between sub-populations, thereby improving the efficiency of evolutionary search and the diversity of the entire population. Allowing each population corresponds to an objective to solve the fitness distribution problem in multi-objective optimization [29]. In addition, populations can update search knowledge through the external archive. An elite learning strategy and a modification strategy of the speed update equation are proposed to reform the diversity and convergence, respectively.In addition to population co-evolutionary strategies, decomposition strategies are also an effective way [30]. Direction vectors are used to decompose objective space and the convergence of each sub-objective space is guaranteed by domination decision [31]. Decomposition-based archiving methods are proposed to circumvent the shortcomings of archiving methods [32]. The entire objective space is evenly divided into multiple sub-spaces based on weight vectors. Only one non-dominated solution located in the sub-space is selected for updating the external archive in each generation to improve the diversity. Moreover, the normalized distance and Pareto dominants are used to determine the update of new solutions in the sub-space.This article considers multi-objective optimal design of dropping shock of series cushioning packaging system. The main work is shown as follows:(1)
Based on the stress, strain, and strain rate data of the cushioning materials under the dropping shock load, the BP neural network is adopted to establish the network models reflecting the dynamic constitutive relationships of the cushioning materials to enhance the accuracy of data fitting further.(2)
Establish MOP general flow of the series cushioning packaging system. The dropping shock dynamics model of the series cushioning system is constructed. Then the Runge-Kutta algorithm and the trained network are used to solve the above dynamic model. For the first time, the MOP model of the series cushioning system is proposed. The model considers more factors, such as the stability of the cushioning material and whether to use local buffering.(3)
NSGA-II is applied to solve the MOP model. The solution set is then compared with PESA2, SPEA2, MOPSO, and MOEA/D by using multiple performance indicators. The parameters in the cushioning system on the buffer performance are deeply studied, which supplies a proper selection of the amount of buffering materials. The remains of this paper are shown below. Section2 presents the related work. Section 3 establishes the dynamic model of the series cushioning system. Then Section 4 constructs the dynamic constitutive network models of the cushioning materials based on the BP algorithm. Next, Section 5 shows the MOP model of the system and the flow of NSGA-II. The experimental studies are displayed in Section 6. Subsequently, Section 7 presents the results. Eventually, Section 8 shows the conclusions.


## 2. Related Work

The constitutive model of cushioning material represents the relationship between the force and deformation of the material and is generally characterized by a stress-strain curve [33]. The curve consists of three areas: linear elasticity, plateau, and densification region. Scholars have discovered that many factors affect the stress-strain relationship [34]. Under different strain rates, the stress-strain relationship of polyurethane foam has been deeply studied [35]. The relationship between stress and strain is decomposed into the product of two terms: the strain’s shape function and the strain rate’s function. Then a new constitutive model is presented to improve the accuracy of the fitting. The viscoelastic manners of polyurethane foam at different strain rates and temperatures have been studied [36]. The nonlinear manners of polyurethane foam at different conditions (strain rate, density, foaming direction) are analyzed by adapting the Frank-Brockman-Zairi elastic viscoplastic constitutive model [37].However, the above methods of constructing constitutive models require researchers to build functions to describe complex stress-strain relationships. Since the 1990s, the method of building material constitutive models based on neural networks [38] has been proposed and has received widespread attention from scholars [39–41]. The neural network constitutive model directly obtains material behavior information from the experimental data set. These pioneering works have confirmed the effectiveness of using neural network methods for constitutive model modeling. Thus this article intends to use neural network methods to establish constitutive models for common cushioning materials like corrugated paperboard and expanded polyethylene.The drop impact response of the corrugated paperboard is studied [42]. Based on the constitutive model of corrugated paperboard, a dynamic model of the drop impact process is established. Then the variational iteration method is employed to solve the kinetic equation, and the approximate solution is obtained. The approximate solution is compared with the Runge-Kutta numerical algorithm, and the result shows that the approximate solution is very close to the numerical solution. For the series cushioning packaging structure, the virtual mass method is proposed [43, 44] and used for model conversion to be suitable for the Runge-Kutta numerical method. The dynamic compression response of the trapezoidal polyethylene foam structure is studied [45]. The results obtained by the virtual mass method and Runge-Kutta algorithm are compared with the finite element method. Parameters in constitutive models are identified with the least square method [46]. The study determines the design strategy of the series cushioning packaging system during the drop impact process. However, there are still some problems in the above research: (1) It is believed that the corrugated paperboard and polyethylene have the same area, and both are smaller than the bottom area of the product. Only consider minimizing the amount of polyethylene. (2) The stability constraint of the cushioning material is not considered in the model. (3) The iteration search is performed after deterministic meshing of the parameter space (the area and thickness of the polyethylene); the obtained solution is not optimal in the single-objective problem.


## 3. Dynamic Model of Series Cushioning Packaging System

Assume that the mass and bottom area of the product arem and A, respectively. The product and n pieces of cushioning materials (the bottom cushioning material is regarded as the outer packing box) form a series cushioning system that is a package. The package falls freely from a height of H. The displacement of the product in contact with the first piece of cushioning material is set to y1. The displacement of the contact interface between the i-th (i = 1, 2, …, n − 1) piece and the i + 1th piece of cushioning material is set to yi+1. The mass, thickness, and bottom area of the i-th cushioning material are mi, hi, and Ai, respectively.Figure1 is a schematic diagram of the dropping shock of the series cushioning system. The strain εit and strain rate ε˙it of the i-th cushioning material at time t are denoted as εit=yit−yi+1t/hi and ε˙it=y˙it−y˙i+1t/hi, respectively. fiεi,ε˙i is the dynamic constitutive equation of the i-th cushioning material, which depicts the relationship between the stress, εit, and ε˙it. In order to show more simply, the strain εt and strain rate ε˙t are denoted as ε and ε˙, respectively.Figure 1
Schematic of series cushioning packaging structure.Assuming that the product is a continuous and uniform rigid body, the dynamic equation of the product is as follows:(1)my¨1+A1f1ε1,ε˙1−mg=0,where y¨1 is the acceleration at the point where the product is in contact with the first cushioning material, and A1 is the area of the first cushioning material. g is the acceleration of gravity.Force is continuous at the interface of two continuous cushioning materials. Thus, the following boundary conditions are presented:(2)A1f1ε1,ε˙1=A2f2ε2,ε˙2A2f2ε2,ε˙2=A3f3ε3,ε˙3...An−1fn−1εn−1,ε˙n−1=Anfnεn,ε˙n,where i∈1,n−1, and Ai represents the area of the i-th cushioning material. εn=yn/hn and ε˙n=y˙n/hn individually denote the strain and strain rate of the n-th cushioning material in contact with the ground in the series cushioning system.The virtual mass method is introduced in [43, 44]. Equation (2) can be transformed into an explicit differential equation by adding virtual mass m1; then the dynamic equation of the series cushioning system is as follows:(3)my¨1+A1f1ε1,ε˙1−mg=0m1y¨2−A1f1ε1,ε˙1+A2f2ε2,ε˙2=0m1y¨3−A2f2ε2,ε˙2+A3f3ε3,ε˙3=0...m1y¨n−An−1fn−1εn−1,ε˙n−1+Anfnεn,ε˙n=0.When the package falls freely, and the lowest cushioning material touches the ground, the system has the same initial velocity everywhere, with the value2gH. However, the system has not yet deformed, so the initial conditions of the system are(4)y10=y20=⋯=yn0=0y˙10=y˙20=⋯=y˙n0=2gH.When the product massm, bottom area A, drop height H, area Ai, thickness hi, and dynamic constitutive relationship fiε,ε˙ are precise, the Runge-Kutta algorithm can be used to solve (3).In addition, the area Ai and thickness hi of the cushioning materials can be changed. They are also critical parameters of MOP in this paper.


## 4. BP Neural Network for the Dynamic Constitutive Models of Cushioning Materials

The constitutive model describes the mechanical properties of materials under external loads. The dynamic constitutive model of cushioning material establishes the relationship between the material’s stress, strain, and strain rate. Generally, a simplified mathematical model is proposed based on experimental data. Then a general form of constitutive relationship is proposed from the theoretical level. Finally, the constitutive model parameters are determined based on experiments or optimization methods. So far, the concrete function expression of the constitutive model is determined.Regrettably, the above strategies require researchers to have rich experience in curve fitting. Otherwise, there will be a significant error between the fitted and experimental data. Therefore, this article builds the dynamic constitutive network models of the cushioning materials. Then the trained datafiε,ε˙ is utilized for constructing the dynamic equation of Section 3.BP neural network can be divided into network structure determination and network prediction. Figure2 is the corresponding algorithm flowchart. The steps are as follows:(1)
The inputs areε and ε˙ in the experimental data, and the output is stress. The entire data is randomly divided into training and testing data. Moreover, normalize the training and test data, respectively.(2)
The network structure is defined according to the input and output dimensions of the data. That is, it sets hidden layers and nodes in the network. Then define the network’s node transfer function, the number of iterations, learning rate and training goals, and so on.(3)
The network is trained by the training data. When the termination condition is met, the trained network is output.Figure 2
The flowchart of BP.


### 4.1. Dynamic Constitutive Network Model of Corrugated Paperboard

The package usually utilizes corrugated paperboard as the outer box and expanded polyethylene as the inner cushioning material. Therefore, in this article, the dynamic constitutive neural network models of these two cushioning materials are constructed, respectively.Dynamic stress-strain curves of corrugated paperboard and expanded polyethylene [47] are shown in Figure 3. It can be seen that the two experimental curves are nonlinear, and the accuracy of the fitted data cannot be guaranteed by using traditional curve fitting methods such as polynomial fitting.Figure 3
Dynamic stress-strain experimental data of the two cushioning materials. (a) Corrugated paperboard; (b) expanded polyethylene.
(a)(b)Accordingly, considering the influence of strain rate on the stress-strain curve, the training models of these two materials are established, respectively. The numeral of nodes in each hidden layer is 5, the hidden layer is 2, iterations are 200, the learning rate is 0.1, and the training goal is 10−5.The network prediction result of corrugated paperboard is presented in Figure4. The first two sub-graphs are the error between the expected and predicted outputs. The third sub-graph is the ratio of the error to the expected output, and the fourth sub-graph is the result of data fitting using the trained BP network. The mean square error (MSE) between the network fitting and experimental data is 5.6433e − 07. The MSE between the polynomial fitting data based on the least square method [47] and the experimental data is 0.0046. In addition, it cannot effectively reproduce the entire compression process. It can be perceived that BP has compelling modeling capability in constructing the dynamic constitutive model of corrugated paperboard.Figure 4
Prediction result of corrugated paperboard. (a) and (b) show the error between expected and predicted output; (c) shows the ratio of error to expected output; (d) shows the result of data fitting using the trained BP network.
(a)(b)(c)(d)


### 4.2. Dynamic Constitutive Network Model of Expanded Polyethylene

The prediction result of expanded polyethylene is displayed in Figure5. The first two subdiagrams are the error between the expected and predicted outputs. The third subdiagram is the ratio of the error to the expected result, and the fourth subdiagram is the outcome of data fitting using the trained BP network. MSE between the neural network fitting data and the experimental data is 5.3360e − 07. MSE between the polynomial fitting data based on the least square method [47] and the experimental data is 4.6247e − 04. It can be noted that BP has a potent modeling power in the construction of the dynamic constitutive model of expanded polyethylene.Figure 5
Prediction result of expanded polyethylene. (a) and (b) show the error between expected and predicted output; (c) shows the ratio of error to expected output; (d) shows the result of data fitting using the trained BP network.
(a)(b)(c)(d)


### 4.3. Parametric Study

To further study the influence of parameters on the fitting accuracy, based on the experimental data of corrugated paperboard, MSE between predicted and expected output is selected as an evaluation index. For the sake of fairness, each set of parameters is run 100 times to get the median of MSE.Table1 shows the median of MSE of different hidden layer nodes. The element values and lengths in the array represent the numeral of nodes in each layer and the numeral of hidden layers. When the hidden layer is single, as the number of nodes increases, the median of MSE first decreases and then increases. When the number of nodes is 13, the median is the lowest. There is a similar trend when the hidden layer is two-layer and three-layer, which is more prominent.Table 1
Comparison results of the median of MSE of different hidden layer nodes.
Hidden layerMedianHidden layerMedianHidden layerMedian[3]1.59E − 04[3, 3]1.91E − 05[3, 3, 3]7.00E − 06[4]1.02E − 04[4, 4]3.07E − 06[4, 4, 3]1.53E − 06[5]6.43E − 05[5, 5]1.68E − 06[5, 5, 5]1.56E − 06[6]2.40E − 05[6, 6]1.45E − 06[6, 6, 6]2.23E − 06[7]1.38E − 05[7, 7]1.99E − 06[7, 7, 7]2.68E − 06[8]9.38E − 06[8, 8]2.73E − 06[8, 8, 8]4.17E − 06[9]7.30E − 06[9, 9]3.69E − 06[9, 9, 9]7.00E − 06[10]6.24E − 06[10, 10]5.50E − 06[10, 10, 10]6.69E − 06[11]4.75E − 06[11, 11]8.33E − 06[11, 11, 11]9.80E − 06[12]3.90E − 06[12, 12]9.65E − 06[12, 12, 12]1.45E − 05[13]3.62E − 06[13, 13]1.63E − 05[13, 13, 13]1.38E − 05[14]3.70E − 06[14, 14]1.53E − 05[14, 14, 14]1.22E − 05[15]4.50E − 06[15, 15]1.43E − 05[15, 15, 15]1.29E − 05[20]5.05E − 06[20, 20]2.90E − 05[20, 20, 20]1.57E − 05[25]8.28E − 06[25, 25]2.78E − 05[25, 25, 25]2.08E − 05When nodes are small and the numeral of hidden layers increases, the median will decrease. However, as nodes increase, the expansion of layers will cause the median to become larger. Therefore, the numerals of nodes and layers need to be designed reasonably according to the complexity of the problem. This paper individually sets the numeral of hidden layers and nodes in each layer to 2 and 5.


## 4.1. Dynamic Constitutive Network Model of Corrugated Paperboard

The package usually utilizes corrugated paperboard as the outer box and expanded polyethylene as the inner cushioning material. Therefore, in this article, the dynamic constitutive neural network models of these two cushioning materials are constructed, respectively.Dynamic stress-strain curves of corrugated paperboard and expanded polyethylene [47] are shown in Figure 3. It can be seen that the two experimental curves are nonlinear, and the accuracy of the fitted data cannot be guaranteed by using traditional curve fitting methods such as polynomial fitting.Figure 3
Dynamic stress-strain experimental data of the two cushioning materials. (a) Corrugated paperboard; (b) expanded polyethylene.
(a)(b)Accordingly, considering the influence of strain rate on the stress-strain curve, the training models of these two materials are established, respectively. The numeral of nodes in each hidden layer is 5, the hidden layer is 2, iterations are 200, the learning rate is 0.1, and the training goal is 10−5.The network prediction result of corrugated paperboard is presented in Figure4. The first two sub-graphs are the error between the expected and predicted outputs. The third sub-graph is the ratio of the error to the expected output, and the fourth sub-graph is the result of data fitting using the trained BP network. The mean square error (MSE) between the network fitting and experimental data is 5.6433e − 07. The MSE between the polynomial fitting data based on the least square method [47] and the experimental data is 0.0046. In addition, it cannot effectively reproduce the entire compression process. It can be perceived that BP has compelling modeling capability in constructing the dynamic constitutive model of corrugated paperboard.Figure 4
Prediction result of corrugated paperboard. (a) and (b) show the error between expected and predicted output; (c) shows the ratio of error to expected output; (d) shows the result of data fitting using the trained BP network.
(a)(b)(c)(d)


## 4.2. Dynamic Constitutive Network Model of Expanded Polyethylene

The prediction result of expanded polyethylene is displayed in Figure5. The first two subdiagrams are the error between the expected and predicted outputs. The third subdiagram is the ratio of the error to the expected result, and the fourth subdiagram is the outcome of data fitting using the trained BP network. MSE between the neural network fitting data and the experimental data is 5.3360e − 07. MSE between the polynomial fitting data based on the least square method [47] and the experimental data is 4.6247e − 04. It can be noted that BP has a potent modeling power in the construction of the dynamic constitutive model of expanded polyethylene.Figure 5
Prediction result of expanded polyethylene. (a) and (b) show the error between expected and predicted output; (c) shows the ratio of error to expected output; (d) shows the result of data fitting using the trained BP network.
(a)(b)(c)(d)


## 4.3. Parametric Study

To further study the influence of parameters on the fitting accuracy, based on the experimental data of corrugated paperboard, MSE between predicted and expected output is selected as an evaluation index. For the sake of fairness, each set of parameters is run 100 times to get the median of MSE.Table1 shows the median of MSE of different hidden layer nodes. The element values and lengths in the array represent the numeral of nodes in each layer and the numeral of hidden layers. When the hidden layer is single, as the number of nodes increases, the median of MSE first decreases and then increases. When the number of nodes is 13, the median is the lowest. There is a similar trend when the hidden layer is two-layer and three-layer, which is more prominent.Table 1
Comparison results of the median of MSE of different hidden layer nodes.
Hidden layerMedianHidden layerMedianHidden layerMedian[3]1.59E − 04[3, 3]1.91E − 05[3, 3, 3]7.00E − 06[4]1.02E − 04[4, 4]3.07E − 06[4, 4, 3]1.53E − 06[5]6.43E − 05[5, 5]1.68E − 06[5, 5, 5]1.56E − 06[6]2.40E − 05[6, 6]1.45E − 06[6, 6, 6]2.23E − 06[7]1.38E − 05[7, 7]1.99E − 06[7, 7, 7]2.68E − 06[8]9.38E − 06[8, 8]2.73E − 06[8, 8, 8]4.17E − 06[9]7.30E − 06[9, 9]3.69E − 06[9, 9, 9]7.00E − 06[10]6.24E − 06[10, 10]5.50E − 06[10, 10, 10]6.69E − 06[11]4.75E − 06[11, 11]8.33E − 06[11, 11, 11]9.80E − 06[12]3.90E − 06[12, 12]9.65E − 06[12, 12, 12]1.45E − 05[13]3.62E − 06[13, 13]1.63E − 05[13, 13, 13]1.38E − 05[14]3.70E − 06[14, 14]1.53E − 05[14, 14, 14]1.22E − 05[15]4.50E − 06[15, 15]1.43E − 05[15, 15, 15]1.29E − 05[20]5.05E − 06[20, 20]2.90E − 05[20, 20, 20]1.57E − 05[25]8.28E − 06[25, 25]2.78E − 05[25, 25, 25]2.08E − 05When nodes are small and the numeral of hidden layers increases, the median will decrease. However, as nodes increase, the expansion of layers will cause the median to become larger. Therefore, the numerals of nodes and layers need to be designed reasonably according to the complexity of the problem. This paper individually sets the numeral of hidden layers and nodes in each layer to 2 and 5.


## 5. Multi-Objective Optimization Model of Series Cushioning Packaging System

There are three main factors considered in the design of the series cushioning packaging system:(1)
The shock protection design is reliable. The product’s maximum acceleration is less than the fragility, and it should be less than the given allowable value for safety.(2)
The amounts of cushioning materials are the most economical.(3)
Each layer of cushioning materials must meet the stability constraints, i.e.,Ai>1.33hi2.Therefore, the series cushioning packaging system design is a typical MOP, and the mathematical model is shown in(5)minV1,V2,…,Vn,maxy¨1V1=A1h1,V2=A2h2,…,Vn=Anhnsubject tomaxy¨1≤Gmns0≤Ai≤A,A<An<A00<hi<h0Stability Constraint,where Vi is the amount of the i-th cushioning material. y¨1 is the acceleration of the product. A is the bottom area of the product. An is the bottom area of the corrugated paperboard. A0 is the maximum allowable bottom area of the package. h0 is the maximum allowable thickness of cushioning material. hn is the corrugated paperboard’s thickness, which has specific requirements in industrial production and cannot be changed at pleasure. Gm is the product’s fragility, and ns >1 is the safety factor. Except for corrugated paperboard, the area Ai of other cushioning materials should be less than or equal to the bottom size A of the product, which means partial buffer or full buffer. The stability constraint is presented as follows:(6)if full buffer,Ai>1.33hi2if partial buffer,14Ai>1.33hi2An>1.33hn2,Due to the limited space, we chose NSGA-II for a brief introduction. The pseudocode is as follows:Step 1: randomly generate individuals (or chromosomes) in the feasible parameter space, and the individuals are encoded in real numbers. The set of chromosomes forms the initial population (Pop).Step 2: compute objectives of the individual in population (Pop), that is, the amount of expanded polyethylene, corrugated paperboard, and the maximum response acceleration of the product under this set of parameters.Step 3: generate an empty set (NDPop) to stock non-dominated individuals discovered.Step 4: while the iteration condition has not been reached, do(1)
Use the selection operator to randomly choose individuals from the population (Pop) to form a new set (Pops).(2)
Randomly pick two individuals from the set (Pops) in turn, and use the crossover operator to generate two new individuals. Repeat this process until the individuals in the set (Pops) have performed the selection operator. Store the new individuals into a new set (Popc).(3)
Use the mutation operator on individuals in the population (Pop) to generate new individuals and store the new individuals into a new set (Popm).(4)
Merge the two sets (Popc) and (Popm) into the population (Pop). Perform non-dominated sorting and crowding measurement on the population (Pop). Retain individuals with high nondominant ranking and low crowding in the population (Pop), and keep the population (Pop) volume unchanged.(5)
Non-dominated individuals of population (Pop) are backed up to set (NDPop). Perform non-dominated sorting of set (NDPop). Retain non-dominated individuals of set (NDPop).(6)
Increase the number of population iterations by 1.In addition, it should be noted that, in the entire algorithm execution process, each newly generated individual must make the stability judgment. The individuals who do not meet the stability constraints are corrected to enhance the computational efficiency of the algorithm.So far, the resolution of MOP has been completed. The steps of the whole algorithm are shown as follows:Step 1: construct dynamic constitutive network models of the materials based on experimental dataStep 2: the dynamic equation of the series cushioning system is established. Then, the Runge–Kutta algorithm is adopted to solve the dynamic equation. In the end, the law of the product acceleration changing with time can be obtained.Step 3: establish the MOP model of series cushioning system and utilize MOEA to solve. Then, non-dominated solutions can be received.


## 6. Experimental Studies

### 6.1. Experimental Settings

The parameters of the multi-objective optimization problem are defined as follows: the product mass is set to 6 kg, the bottom area is 0.01 m2, the fragility value Gm is 85 g, the safety factor ns is 1.2, and the maximum allowable packaging bottom area is 0.04 m2. The thickness of corrugated paperboard is generally a fixed value in actual production, which is set to 6.7 mm in this paper. The product, the expanded polyethylene, and the corrugated paperboard form the series cushioning packaging system and then fall from a height of 0.5 m.This paper compares the results computed by various MOEAs. These algorithms include a non-dominated sorting genetic algorithm (NSGA-II [48]), modified Pareto envelope-based selection algorithm (PESA2 [49]), strength Pareto evolutionary algorithm (SPEA2 [50]), multi-objective particle swarm optimization (MOPSO [51]), and multi-objective evolutionary algorithm based on decomposition (MOEA/D [52]). These algorithms are chosen because they perform pleasingly on most MOP problems. Therefore, it makes comparisons more believable.Table2 shows the parameters of algorithms. For fairness, the capacity of non-dominated solutions is set to 100.Table 2
Parameters of algorithms.
AlgorithmsParametersPESA2MaxIt = 100, nPop = 80, nArchive = 100, nGrid = 7,pc = 0.8, pm = 0.1SPEA2MaxIt = 100, nPop = 80, nArchive = 100,K = 13, pc = 0.8, pm = 0.1MOPSOMaxIt = 100, nPop = 80, nArchive = 100,c1 = 1, c2 = 2, ω = 0.5 ⟶ 0.2MOEA/DMaxIt = 100, nPop = 80, nArchive = 100, nNeigh = 4,pc = 0.8, pm = 0.1NSGA-IIMaxIt = 100, nPop = 80, nArchive = 100,pc = 0.8, pm = 0.1


### 6.2. Performance Metric

Performance evaluation is required to evaluate the approximate Pareto solution sets obtained by MOP algorithms. The performance indicators can be divided into three categories:(1)
Only evaluate the approaching extent of Pareto approximate solution set to Pareto front (PF), that is, the convergence (e.g., C-metric [53] and GD [54]).(2)
Only evaluate the distribution of the approximate solution set on PF, that is, diversity (e.g., ∆p [55] and CPF [56]). Diversity assessment is further divided into uniformity (e.g., spacing [57]) and spread (e.g., PD [58]) of the solution set distribution.(3)
Comprehensive indicators measure convergence and diversity (e.g., IGD [59] and HVR [60]).C-metric compares convergence of two solution sets and presents below:(7)CA,B=u∈B|∃v∈A:vweakly dominatesuB.The value ranges from 0 to 1. C(A, B) = 1 implies that some individuals of set A dominate each individual in set B. That is, the convergence of B is worse than A. Conversely, C(A, B) = 0 implies that any solutions of set A do not dominate all solutions in set B. That is, the convergence of B is upper than A. The larger the C(A, B), the upper the convergence of A. Note that the sum of C(A, B) and C(B, A) does not have to be 1.Generational distance (GD) calculates the average minimum distance from each solution of set S to setP∗. Set P∗ contains reference points, which are uniformly spread on PF. S is a set of solutions calculated by the algorithm.(8)GDS,P∗=∑i=1Sdi2S,where |S| is the numeral of solutions in set S, and di denotes Euclidean distance from the i-th solution of set S to the closest solution on set P∗. The less the GD is, the higher the convergence of S is, and the closer it stands to the whole PF.∆p is the average Hausdorff distance from set S to PF. The smaller ∆p, the higher the diversity of S.Coverage over Pareto Front (CPF) replaces each pointx of set S with the point closest to x on reference solution set and forms a new solution set S∗. CPF is the set S∗ coverage on reference point set. The larger CPF, the higher the diversity of S.Spacing evaluates diversity by computing the standard deviation of minimum Euclidean distance between each point of set S and others. The less the spacing, the better the uniformity S. The calculation is as follows:(9)SpacingS=1S−1∑i=1Sd¯−di2,where di denotes the minimum distance, and d¯ represents the average of all di.Pure distance (PD) evaluates diversity by calculating dissimilarity, which is described as follows:(10)PDS=maxx∈SPDS\x+miny∈S/xfx−fyp,where fx−fyp represents the Lp-norm-based distance. PD mainly considers the spread of S. The higher the PD, the better the diversity.Inverted Generational Distance (IGD) calculates the average of minimum distances betweenP∗ and S. The smaller the IGD, the higher the algorithm’s all-around performance.(11)IGDS,P∗=∑x∈P∗distx,SP∗.Hypervolume Ratio (HVR) refers to the ratio of the hypervolume of S to PF. The higher the HVR, the more the power of the comprehensive performance of solution set S.For fairness, each algorithm runs 30 times independently and then calculates performance indicators’ mean and standard deviation. The best in each row is emphasized with bold font. In addition, the Wilcoxon rank-sum test with a significance level ofα = 0.05 is utilized.


### 6.3. Experimental Results

Table3 presents the comparison results. ∆p values of all comparison algorithms are the same as the IGD values, which means that all the evolutionary algorithms involved in the comparison have good convergence performance [56]. Moreover, the GD of NSGA-II is the smallest, but it is similar to other algorithms, confirming the above point of view. Both ∆p and CPF of NSGA-II rank first, implying NSGA-II has the best diversity performance. IGD and HVR of NSGA-II are both the first, indicating NSGA-II has the best comprehensive performance. NSGA-II also ranks first in all indicators overall.Table 3
Comparison results of performance metric.
MetricPESA2SPEA2MOPSOMOEA/DNSGA-IIGDMean3.40E − 033.40E − 033.30E − 033.20E − 033.10E − 03Std2.99E − 042.25E − 042.46E − 043.16E − 042.40E − 04Rank5−4−3−2≈1∆pMean3.15E − 011.06 E + 001.10 E + 001.09 E+001.97E − 01Std1.36E − 012.74E − 014.27E − 015.17E − 013.69E − 02Rank2−3−5−4−1CPFMean4.71E − 012.04E − 014.45E − 011.53E − 015.33E − 01Std6.07E − 025.04E − 025.36E − 023.81E − 025.03E − 02Rank2−4−3−5−1SpacingMean2.81E − 011.31E − 011.40E − 012.89E − 011.94E − 01Std1.37E − 011.92E − 021.90E − 021.45E − 012.41E − 02Rank4−1+2+5−3PDMean5.85 E + 034.45 E + 036.40 E + 034.24 E + 035.69 E + 03Std6.47 E + 025.30 E + 026.98 E + 026.94 E + 025.67 E + 02Rank2≈4−1+5−3IGDMean3.15E − 011.06 E + 001.10 E + 001.09 E + 001.97E − 01Std1.36E − 012.74E − 014.27E − 015.17E − 013.69E − 02Rank2−3−5−4−1HVRMean9.73E − 018.89E − 018.89E − 019.21E − 019.95E − 01Std2.18E − 021.24E − 021.87E − 022.73E − 021.00E − 03Rank2−4−5−3−1Total rank1923242811Table4 shows the C-metric comparison results. The comparisons indicate the convergence of NSGA-II is the same as SPEA2 and MOPSO. The convergence of NSGA-II is a little poor than PESA2 and somewhat better than MOEA/D. This conclusion conflicts with the ordering of GD values. It also shows a certain degree of deviation between the C-metric and the GD measurement results.Table 4
Comparison results of C-metric.
MetricPESA2SPEA2MOPSOMOEA/DC (NSGA-II, −)Mean6.57E − 029.03E − 021.41E − 013.65E − 01Std3.35E − 028.79E − 021.14E − 011.65E − 01C (−, NSGA-II)Mean2.02E − 019.51E − 021.06E − 014.47E − 02Std7.46E − 026.52E − 025.99E − 024.36E − 02Figure6 illustrates non-dominated solutions after running each algorithm 30 times. Solutions received by NSGA-II and PESA2 have better distribution, with MOPSO and MOEA/D coming in second. The solution set distribution calculated by SPEA2 is the worst.Figure 6
Final non-dominated solutions of series cushioning system. (a) PESA2. (b) SPEA2. (c) MOPSO. (d) MOEA/D. (e) NSGA-II.
(a)(b)(c)(d)(e)


### 6.4. Impacts of Parameter Settings

Research the influence of population size, algebra, crossover rate, and mutation rate in NSGA-II on algorithm accuracy. The likelihood of optimality (Lopt) [26, 61] is utilized to metric non-dominated solutions under diverse parameters. Based on the same parameter settings, the algorithm runs n times independently, and the number of times the optimal solution is found in each dimension is m1, m2, ..., md. The Loptk of parameters k is m1+m2+⋯+md/d×n. In this article, n is set to 30.Figure7 shows the Lopt comparison results of different parameters. It can be seen that when the population algebra exceeds 60, the population size exceeds 60, the crossover rate exceeds 0.4, and the mutation rate exceeds 0.1, and Lopt is already close to 0.9, which implies that the algorithm has excellent stability. Therefore, in this paper, the population is 80, the generations are 100, and the selection, crossover, and mutation are 1, 0.8, and 0.1, respectively.Figure 7
TheLopt comparison results of different parameters: (a) population generation; (b) population size; (c) crossover rate; (d) mutation rate.
(a)(b)(c)(d)


## 6.1. Experimental Settings

The parameters of the multi-objective optimization problem are defined as follows: the product mass is set to 6 kg, the bottom area is 0.01 m2, the fragility value Gm is 85 g, the safety factor ns is 1.2, and the maximum allowable packaging bottom area is 0.04 m2. The thickness of corrugated paperboard is generally a fixed value in actual production, which is set to 6.7 mm in this paper. The product, the expanded polyethylene, and the corrugated paperboard form the series cushioning packaging system and then fall from a height of 0.5 m.This paper compares the results computed by various MOEAs. These algorithms include a non-dominated sorting genetic algorithm (NSGA-II [48]), modified Pareto envelope-based selection algorithm (PESA2 [49]), strength Pareto evolutionary algorithm (SPEA2 [50]), multi-objective particle swarm optimization (MOPSO [51]), and multi-objective evolutionary algorithm based on decomposition (MOEA/D [52]). These algorithms are chosen because they perform pleasingly on most MOP problems. Therefore, it makes comparisons more believable.Table2 shows the parameters of algorithms. For fairness, the capacity of non-dominated solutions is set to 100.Table 2
Parameters of algorithms.
AlgorithmsParametersPESA2MaxIt = 100, nPop = 80, nArchive = 100, nGrid = 7,pc = 0.8, pm = 0.1SPEA2MaxIt = 100, nPop = 80, nArchive = 100,K = 13, pc = 0.8, pm = 0.1MOPSOMaxIt = 100, nPop = 80, nArchive = 100,c1 = 1, c2 = 2, ω = 0.5 ⟶ 0.2MOEA/DMaxIt = 100, nPop = 80, nArchive = 100, nNeigh = 4,pc = 0.8, pm = 0.1NSGA-IIMaxIt = 100, nPop = 80, nArchive = 100,pc = 0.8, pm = 0.1


## 6.2. Performance Metric

Performance evaluation is required to evaluate the approximate Pareto solution sets obtained by MOP algorithms. The performance indicators can be divided into three categories:(1)
Only evaluate the approaching extent of Pareto approximate solution set to Pareto front (PF), that is, the convergence (e.g., C-metric [53] and GD [54]).(2)
Only evaluate the distribution of the approximate solution set on PF, that is, diversity (e.g., ∆p [55] and CPF [56]). Diversity assessment is further divided into uniformity (e.g., spacing [57]) and spread (e.g., PD [58]) of the solution set distribution.(3)
Comprehensive indicators measure convergence and diversity (e.g., IGD [59] and HVR [60]).C-metric compares convergence of two solution sets and presents below:(7)CA,B=u∈B|∃v∈A:vweakly dominatesuB.The value ranges from 0 to 1. C(A, B) = 1 implies that some individuals of set A dominate each individual in set B. That is, the convergence of B is worse than A. Conversely, C(A, B) = 0 implies that any solutions of set A do not dominate all solutions in set B. That is, the convergence of B is upper than A. The larger the C(A, B), the upper the convergence of A. Note that the sum of C(A, B) and C(B, A) does not have to be 1.Generational distance (GD) calculates the average minimum distance from each solution of set S to setP∗. Set P∗ contains reference points, which are uniformly spread on PF. S is a set of solutions calculated by the algorithm.(8)GDS,P∗=∑i=1Sdi2S,where |S| is the numeral of solutions in set S, and di denotes Euclidean distance from the i-th solution of set S to the closest solution on set P∗. The less the GD is, the higher the convergence of S is, and the closer it stands to the whole PF.∆p is the average Hausdorff distance from set S to PF. The smaller ∆p, the higher the diversity of S.Coverage over Pareto Front (CPF) replaces each pointx of set S with the point closest to x on reference solution set and forms a new solution set S∗. CPF is the set S∗ coverage on reference point set. The larger CPF, the higher the diversity of S.Spacing evaluates diversity by computing the standard deviation of minimum Euclidean distance between each point of set S and others. The less the spacing, the better the uniformity S. The calculation is as follows:(9)SpacingS=1S−1∑i=1Sd¯−di2,where di denotes the minimum distance, and d¯ represents the average of all di.Pure distance (PD) evaluates diversity by calculating dissimilarity, which is described as follows:(10)PDS=maxx∈SPDS\x+miny∈S/xfx−fyp,where fx−fyp represents the Lp-norm-based distance. PD mainly considers the spread of S. The higher the PD, the better the diversity.Inverted Generational Distance (IGD) calculates the average of minimum distances betweenP∗ and S. The smaller the IGD, the higher the algorithm’s all-around performance.(11)IGDS,P∗=∑x∈P∗distx,SP∗.Hypervolume Ratio (HVR) refers to the ratio of the hypervolume of S to PF. The higher the HVR, the more the power of the comprehensive performance of solution set S.For fairness, each algorithm runs 30 times independently and then calculates performance indicators’ mean and standard deviation. The best in each row is emphasized with bold font. In addition, the Wilcoxon rank-sum test with a significance level ofα = 0.05 is utilized.


## 6.3. Experimental Results

Table3 presents the comparison results. ∆p values of all comparison algorithms are the same as the IGD values, which means that all the evolutionary algorithms involved in the comparison have good convergence performance [56]. Moreover, the GD of NSGA-II is the smallest, but it is similar to other algorithms, confirming the above point of view. Both ∆p and CPF of NSGA-II rank first, implying NSGA-II has the best diversity performance. IGD and HVR of NSGA-II are both the first, indicating NSGA-II has the best comprehensive performance. NSGA-II also ranks first in all indicators overall.Table 3
Comparison results of performance metric.
MetricPESA2SPEA2MOPSOMOEA/DNSGA-IIGDMean3.40E − 033.40E − 033.30E − 033.20E − 033.10E − 03Std2.99E − 042.25E − 042.46E − 043.16E − 042.40E − 04Rank5−4−3−2≈1∆pMean3.15E − 011.06 E + 001.10 E + 001.09 E+001.97E − 01Std1.36E − 012.74E − 014.27E − 015.17E − 013.69E − 02Rank2−3−5−4−1CPFMean4.71E − 012.04E − 014.45E − 011.53E − 015.33E − 01Std6.07E − 025.04E − 025.36E − 023.81E − 025.03E − 02Rank2−4−3−5−1SpacingMean2.81E − 011.31E − 011.40E − 012.89E − 011.94E − 01Std1.37E − 011.92E − 021.90E − 021.45E − 012.41E − 02Rank4−1+2+5−3PDMean5.85 E + 034.45 E + 036.40 E + 034.24 E + 035.69 E + 03Std6.47 E + 025.30 E + 026.98 E + 026.94 E + 025.67 E + 02Rank2≈4−1+5−3IGDMean3.15E − 011.06 E + 001.10 E + 001.09 E + 001.97E − 01Std1.36E − 012.74E − 014.27E − 015.17E − 013.69E − 02Rank2−3−5−4−1HVRMean9.73E − 018.89E − 018.89E − 019.21E − 019.95E − 01Std2.18E − 021.24E − 021.87E − 022.73E − 021.00E − 03Rank2−4−5−3−1Total rank1923242811Table4 shows the C-metric comparison results. The comparisons indicate the convergence of NSGA-II is the same as SPEA2 and MOPSO. The convergence of NSGA-II is a little poor than PESA2 and somewhat better than MOEA/D. This conclusion conflicts with the ordering of GD values. It also shows a certain degree of deviation between the C-metric and the GD measurement results.Table 4
Comparison results of C-metric.
MetricPESA2SPEA2MOPSOMOEA/DC (NSGA-II, −)Mean6.57E − 029.03E − 021.41E − 013.65E − 01Std3.35E − 028.79E − 021.14E − 011.65E − 01C (−, NSGA-II)Mean2.02E − 019.51E − 021.06E − 014.47E − 02Std7.46E − 026.52E − 025.99E − 024.36E − 02Figure6 illustrates non-dominated solutions after running each algorithm 30 times. Solutions received by NSGA-II and PESA2 have better distribution, with MOPSO and MOEA/D coming in second. The solution set distribution calculated by SPEA2 is the worst.Figure 6
Final non-dominated solutions of series cushioning system. (a) PESA2. (b) SPEA2. (c) MOPSO. (d) MOEA/D. (e) NSGA-II.
(a)(b)(c)(d)(e)


## 6.4. Impacts of Parameter Settings

Research the influence of population size, algebra, crossover rate, and mutation rate in NSGA-II on algorithm accuracy. The likelihood of optimality (Lopt) [26, 61] is utilized to metric non-dominated solutions under diverse parameters. Based on the same parameter settings, the algorithm runs n times independently, and the number of times the optimal solution is found in each dimension is m1, m2, ..., md. The Loptk of parameters k is m1+m2+⋯+md/d×n. In this article, n is set to 30.Figure7 shows the Lopt comparison results of different parameters. It can be seen that when the population algebra exceeds 60, the population size exceeds 60, the crossover rate exceeds 0.4, and the mutation rate exceeds 0.1, and Lopt is already close to 0.9, which implies that the algorithm has excellent stability. Therefore, in this paper, the population is 80, the generations are 100, and the selection, crossover, and mutation are 1, 0.8, and 0.1, respectively.Figure 7
TheLopt comparison results of different parameters: (a) population generation; (b) population size; (c) crossover rate; (d) mutation rate.
(a)(b)(c)(d)


## 7. Results

From the optimal solutions found by NSGA-II presented in Figure6, the following conclusions can be observed:(1)
From a vertical perspective, when the amount of expanded polyethylene is fixed, the increase in the amount of corrugated paperboard will reduce the maximum response acceleration of the product to a certain extent. That is, the cushioning effect of the outer packing box cannot be ignored.(2)
From a horizontal view, when the amount of the outer packaging box is fixed, increasing the amount of expanded polyethylene will reduce the maximum response acceleration of the product to a certain degree, that is, improve the protection performance of the series cushioning system.(3)
From the trend in the yellow area, it can be seen that, under the same packaging protection and safety, reducing the amount of any cushioning material will inevitably lead to the increase of other cushioning materials.(4)
When the drop height is changed to 0.7 m and other parameters remain unchanged, the solution set computed with NSGA-II is displayed in Figure8. The non-dominated solutions set moves towards the right in the solution space when the product quality increases. Because the thickness of corrugated paperboard is fixed, when the fall altitude increases, increasing the area of the corrugated paperboard is not enough to meet the safety requirements. The cushioning performance can only be improved by increasing the amount of expanded polyethylene. It also reflects the limited cushioning performance of corrugated paperboard compared with expanded polyethylene materials.Figure 8
The non-dominated solutions set of the series cushioning system when the drop height is 0.7 m.The results obtained are compared with experimental results [46] to verify the method’s effectiveness proposed in this paper. The settings of the relevant parameters are the same as those in the literature. The product mass is 6 kg, the fragility value Gm is 100 g, the safety factor ns is 1, the lower surface area of the product is 0.011 m2, and the drop height is 0.6 m. The results obtained by the iteration search method [46] are shown as follows: The bottom area of the expanded polyethylene is 0.007 m2, and the thickness is 0.038 m. The parameters of the knee point of the non-dominated solution set are as follows. The bottom area of the expanded polyethylene is 0.0101 m2, and the thickness is 0.0208 m. The amount of cushioning material has been saved by 21.02%. The reason is that the iteration search method is to divide the parameter space (area and thickness of expanded polyethylene) deterministically and then perform the iteration search. The obtained solution is optimal on the current grid-scale, not the optimal global solution. In comparison, the solution set of the evolutionary algorithm is approached towards PF. In other words, the solution set calculated by NSGA-II has outstanding overall performance.


## 8. Conclusions

This article builds a multi-objective optimization model of the series cushioning packaging system. BP is adopted to establish materials' dynamic constitutive network models to achieve higher-precision data prediction. Various MOEAs are employed to resolve the MOP model, and multiple performance indicators are used to conduct in-depth research on the performance. The comparisons indicate optimal solution set obtained by NSGA-II has good comprehensive performance. In addition, the parameters corresponding to the knee point can effectively save the amount of expanded polyethylene while satisfying the protection safety. Therefore, the method offered in this article is viable and has specific guiding significance for MOP protection design of transport packaging.Due to the limited conditions, the experimental data are all from the literature. The parameters corresponding to the knee point have only been simulated to verify the feasibility, and no experimental verification has been taken. The following action is to build a test platform for verification.For the experimental data in this article, the BP neural network has achieved an excellent fitting effect. However, further verification is required for more complex data in the future.The next step is to establish a more general dynamic constitutive network model that reflects more factors (temperature and density).



---

*Source: 9491102-2022-03-25.xml*