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

,


Introduction
In circulation and transportation, the product is packaged to form a package to protect the product [1]. e 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][4][5][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]. erefore, the response analysis of the package is studied under the dropping shock load in this paper. ree 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. erefore, 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. en 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. e 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][14][15]. BP neural network is a type of neural network that is widely used [16][17][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. us, it shifts the research focus of the fitting problems from the nonlinear characteristics of materials to constructing the correlation between input and output data. erefore, 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][20][21]. ey 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]. e sub-populations evolve independently to maintain their particular characteristics [28]. e 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 multiobjective 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]. e 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.
is article considers multi-objective optimal design of dropping shock of series cushioning packaging system. e 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. e dropping shock dynamics model of the series cushioning system is constructed.
en 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. e 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. e solution set is then compared with PESA2, SPEA2, MOPSO, and MOEA/D by using multiple performance indicators. e parameters in the cushioning system on the buffer performance are deeply studied, which supplies a proper selection of the amount of buffering materials. e remains of this paper are shown below. Section 2 presents the related work. Section 3 establishes the dynamic model of the series cushioning system. en 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. e experimental studies are displayed in Section 6. Subsequently, Section 7 presents the results. Eventually, Section 8 shows the conclusions.

Related Work
e 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]. e 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]. e relationship between stress and strain is decomposed into the product of two terms: the strain's shape function and the strain rate's function.
en a new constitutive model is presented to improve the accuracy of the fitting. e viscoelastic manners of polyurethane foam at different strain rates and temperatures have been studied [36]. e 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][40][41].
e neural network constitutive model directly obtains material behavior information from the experimental data set. ese pioneering works have confirmed the effectiveness of using neural network methods for constitutive model modeling.
us this article intends to use neural network methods to establish constitutive models for common cushioning materials like corrugated paperboard and expanded polyethylene.
e 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. en the variational iteration method is employed to solve the kinetic equation, and the approximate solution is obtained. e 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. e dynamic compression response of the trapezoidal polyethylene foam structure is studied [45]. e 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]. e 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) e stability constraint of the cushioning material is not considered in the model. (3) e 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.

Dynamic Model of Series Cushioning Packaging System
Assume that the mass and bottom area of the product are m and A, respectively. e 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. e package falls freely from a height of H. e displacement of the product in contact with the first piece of cushioning material is set to y 1 . e 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 y i+1 . e mass, thickness, and bottom area of the i-th cushioning material are m i , h i , and A i , respectively. Figure 1 is a schematic diagram of the dropping shock of the series cushioning system. e strain ε i (t) and strain rate _ ε i (t) of the i-th cushioning material at time t are denoted as is the dynamic constitutive equation of the i-th cushioning material, which depicts the relationship between the stress, ε i (t), and _ ε i (t). In order to show more simply, the strain ε(t) and strain rate _ ε(t) are denoted as ε and _ ε, respectively. Assuming that the product is a continuous and uniform rigid body, the dynamic equation of the product is as follows: where y 1 is the acceleration at the point where the product is in contact with the first cushioning material, and A 1 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.
us, the following boundary conditions are presented: where i ∈ (1, n − 1), and A i represents the area of the i-th cushioning material. ε n � y n /h n and _ ε n � _ y n /h n individually denote the strain and strain rate of the n-th cushioning material in contact with the ground in the series cushioning system. e 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: When the package falls freely, and the lowest cushioning material touches the ground, the system has the same initial velocity everywhere, with the value ���� 2gH. However, the system has not yet deformed, so the initial conditions of the system are When the product mass m, bottom area A, drop height H, area A i , thickness h i , and dynamic constitutive relationship f i (ε, _ ε) are precise, the Runge-Kutta algorithm can be used to solve (3).

Shock and Vibration
In addition, the area A i and thickness h i of the cushioning materials can be changed. ey are also critical parameters of MOP in this paper.

BP Neural Network for the Dynamic Constitutive Models of Cushioning Materials
e constitutive model describes the mechanical properties of materials under external loads. e 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. en 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. erefore, this article builds the dynamic constitutive network models of the cushioning materials. en the trained data f i (ε, _ ε) is utilized for constructing the dynamic equation of Section 3.
BP neural network can be divided into network structure determination and network prediction. Figure 2 is the corresponding algorithm flowchart. e steps are as follows: (1) e inputs are ε and _ ε in the experimental data, and the output is stress. e entire data is randomly divided into training and testing data. Moreover, normalize the training and test data, respectively.
(2) e network structure is defined according to the input and output dimensions of the data. at is, it sets hidden layers and nodes in the network. en define the network's node transfer function, the number of iterations, learning rate and training goals, and so on. (3) e network is trained by the training data. When the termination condition is met, the trained network is output.

Dynamic Constitutive Network Model of Corrugated
Paperboard. e package usually utilizes corrugated paperboard as the outer box and expanded polyethylene as the inner cushioning material. erefore, 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.
Accordingly, considering the influence of strain rate on the stress-strain curve, the training models of these two materials are established, respectively. e 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 . e network prediction result of corrugated paperboard is presented in Figure 4. e first two sub-graphs are the error between the expected and predicted outputs. e 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. e mean square error (MSE) between the network fitting and experimental data is 5.6433e − 07. e 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.

Dynamic Constitutive Network Model of Expanded
Polyethylene.
e prediction result of expanded polyethylene is displayed in Figure 5. e first two subdiagrams are the error between the expected and predicted outputs. e 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.

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. Table 1 shows the median of MSE of different hidden layer nodes. e 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. ere is a similar trend when the hidden layer is two-layer and three-layer, which is more prominent.
When 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. erefore, the numerals of nodes and layers need to be designed reasonably according to the complexity of the problem. is paper individually sets the numeral of hidden layers and nodes in each layer to 2 and 5.

Multi-Objective Optimization Model of Series
Cushioning Packaging System ere are three main factors considered in the design of the series cushioning packaging system: (1) e shock protection design is reliable. e product's maximum acceleration is less than the fragility, and it should be less than the given allowable value for safety. (2) e amounts of cushioning materials are the most economical. (3) Each layer of cushioning materials must meet the stability constraints, i.e., erefore, the series cushioning packaging system design is a typical MOP, and the mathematical model is shown in where V i 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. A n is the bottom area of the corrugated paperboard. A 0 is the maximum allowable bottom area of the package. h 0 is the maximum allowable thickness of cushioning material. h n is the corrugated paperboard's thickness, which has specific requirements in industrial production and cannot be changed at pleasure. G m is the product's fragility, and n s >1 is the safety factor. Except for corrugated paperboard, the area A i 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. e stability constraint is presented as follows: if full buffer, Due to the limited space, we chose NSGA-II for a brief introduction. e pseudocode is as follows: Step 1: randomly generate individuals (or chromosomes) in the feasible parameter space, and the     Shock and Vibration individuals are encoded in real numbers. e 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 nondominated individuals discovered.
Step 4: while the iteration condition has not been reached, do   In addition, it should be noted that, in the entire algorithm execution process, each newly generated individual must make the stability judgment. e 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. e steps of the whole algorithm are shown as follows: Step 1: construct dynamic constitutive network models of the materials based on experimental data Step 2: the dynamic equation of the series cushioning system is established. en, 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. en, non-dominated solutions can be received.

Experimental Settings.
e parameters of the multiobjective optimization problem are defined as follows: the product mass is set to 6 kg, the bottom area is 0.01 m 2 , the fragility value G m is 85 g, the safety factor n s is 1.2, and the maximum allowable packaging bottom area is 0.04 m 2 . e thickness of corrugated paperboard is generally a fixed value in actual production, which is set to 6.7 mm in this paper. e product, the expanded polyethylene, and the corrugated paperboard form the series cushioning packaging system and then fall from a height of 0.5 m. is paper compares the results computed by various MOEAs. ese 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 multiobjective evolutionary algorithm based on decomposition (MOEA/D [52]). ese algorithms are chosen because they perform pleasingly on most MOP problems. erefore, it makes comparisons more believable. Table 2 shows the parameters of algorithms. For fairness, the capacity of non-dominated solutions is set to 100.

Performance Metric.
Performance evaluation is required to evaluate the approximate Pareto solution sets obtained by MOP algorithms. e 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]).
Generational distance (GD) calculates the average minimum distance from each solution of set S to set P * . Set P * contains reference points, which are uniformly spread on PF. S is a set of solutions calculated by the algorithm.
where |S| is the numeral of solutions in set S, and d i denotes Euclidean distance from the i-th solution of set S to the closest solution on set P * . e 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. e smaller ∆p, the higher the diversity of S. Coverage over Pareto Front (CPF) replaces each point x 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. e 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. e less the spacing, the better the uniformity S. e calculation is as follows: where d i denotes the minimum distance, and d represents the average of all d i . Pure distance (PD) evaluates diversity by calculating dissimilarity, which is described as follows: where ‖f(x) − f(y)‖ p represents the L p -norm-based distance. PD mainly considers the spread of S. e higher the PD, the better the diversity. Inverted Generational Distance (IGD) calculates the average of minimum distances between P * and S. e smaller the IGD, the higher the algorithm's all-around performance.
Hypervolume Ratio (HVR) refers to the ratio of the hypervolume of S to PF. e 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. e 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. Table 3 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 4 shows the C-metric comparison results. e comparisons indicate the convergence of NSGA-II is the same as SPEA2 and MOPSO. e convergence of NSGA-II is a little poor than PESA2 and somewhat better than MOEA/ D. is 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. Figure 6 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. e solution set distribution calculated by SPEA2 is the worst.

Impacts of Parameter Settings.
Research the influence of population size, algebra, crossover rate, and mutation rate in NSGA-II on algorithm accuracy. e 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 m 1 , m 2 , ..., m d . e Lopt k of parameters k is (m 1 + m 2 + · · · + m d )/d × n. In this article, n is set to 30. Figure 7 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. erefore, in this paper, the population is

Results
From the optimal solutions found by NSGA-II presented in Figure 6, 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. at 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 Figure 8. e 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. e 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. e results obtained are compared with experimental results [46] to verify the method's effectiveness proposed in this paper. e settings of the relevant parameters are the same as those in the literature. e product mass is 6 kg, the fragility value G m is 100 g, the safety factor n s is 1, the lower surface area of the product is 0.011 m 2 , and the drop height is 0.6 m. e results obtained by the iteration search method [46] are shown as follows: e bottom area of the expanded polyethylene is 0.007 m 2 , and the thickness is 0.038 m. e parameters of the knee point of the non-dominated solution set are as follows. e bottom area of the expanded polyethylene is 0.0101 m 2 , and the thickness is 0.0208 m. e amount of cushioning material has been saved by 21.02%. e 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. e obtained solution is optimal on the current gridscale, 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.

Conclusions
is 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. e 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. erefore, 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. e parameters corresponding to the knee point have only been simulated to verify the feasibility, and no experimental verification has been taken. e 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. e next step is to establish a more general dynamic constitutive network model that reflects more factors (temperature and density).

Data Availability
Data are contained within the article.

Conflicts of Interest
e authors declare no conflicts of interest.