Combining response surface optimization and stochastic analysis for crashworthiness design – an introductory study

Objective of the design process are cost effective designs that meet certain expectations with respect to functionality and appearance. Designs are created in an iterative process where analyses of the structural behavior lead to changes in the design. The use of optimization technology makes design changes to be driven directly by analysis results. The application of optimization allows an efficient search for the right combination of design variables for a certain design. Additional use of stochastic methods in order to analyze the design from a statistical standpoint adds robustness to the design and prevents unpleasant surprises in later physical testing. This paper discusses methodology to optimize structures that undergo impact loading. Objective and constraints are transient dynamic responses. The optimization problem is solved using a sequential response surface method. An explicit finite element code is used to solve the transient dynamic problem. The optimization is not performed on results from single simulations but on statistical results from a stochastic analysis. The stochastic analysis is driven using a Monte Carlo method. Commercial software is used for the implementation of the methodology. The results from the study indicate that a combination of optimization and stochastic analysis can add safety margins to a design with respect to robustness against physical errors in the design itself and against changes in load levels and load cases. However, this initial study must be followed up by more in-depth investigations to fully understand the benefits of combined optimization-stochastic analysis.


Introduction
For the analysis of the structural behavior, numerical methods are widely used.The most acknowledged and the most flexible computational methods for structural analysis are the Finite Element Methods (FEM).Using FEM, designs of almost arbitrary complexity and physical behavior can be investigated for their structural responses.This way, large numbers of variations and design modifications can be compared effectively without physical testing.
From the results of computational simulations, conclusions can be drawn to make necessary changes to a design.The determination of such changes is usually driven by the experience of the engineer, and the engineer has normally no evidence that the design is optimal in a mathematical sense.There is also no evidence that the design shows the required robustness with respect to physical errors introduced during the manufacturing process or robustness with respect to load levels and load cases what might differ from the design loads.The objective of this study is to investigate if the combination of optimization and stochastic analysis lead to a better final design considering the above mentioned errors and uncertainties.

Design optimization
Objective of the design process are structures and structural parts that fulfill certain expectations with respect to their economy, functionality and appearance.These requirements can be formulated as an optimization problem for the studied structure.The solution to this problem is found in the design process going through multiple iteration steps.
In structural optimization it is established that complex optimization problems can be solved using sequential approximations [1,2].Beside the local ap- proximation method where the optimization problem is solved by using design sensitivities dΨ i /db j where Ψ i (b) is the vector of responses and b = {b j } is the vector of design variables.There are also methods that allow the search for an optimum in a much larger range, called global approximation.These methods are frequently used when performing optimization on highly nonlinear systems.

Global approximation
In cases where non-linearities make local approximation methods non-feasible, it is convenient to introduce higher order analytical expressions to approximate the dependence between the objective or constraint functions and the design variables.
Let Ψ i (b) be the response of interest.Then, a polynomial Ψi (b) of the degree n can be introduced such that The quantities a i0 , a ip , a ipqr are constants that are determined from the results of the system analyses by solving a linear system of equations.The number of analyses that are necessary depends on the polynomial degree n of Eq. ( 1).The function of Eq. ( 1) is called response surface.
Using Eq. ( 1), an approximation to the classical optimization problem can be established.Thus, The solution of the approximate problem yields a new set of design variables.

The sequential response surface method
A modified version of the standard response surface method, as described above, analyzes the designs as the optimization proceeds and will be called the sequential response surface method [3].The software solution allows the interfacing with any finite element solver.Considering a quadratic response surface The coefficients a i0 , a ip , a ipq are determined using a least square fit of the function Ψi on the existing designs.In the first step s = 1+m designs are analyzed and the coefficients a i0 , a ip are determined.The next m designs are analyzed to determine the coefficients a ipp .Additional designs are used to determine the remaining coefficients as it becomes necessary.After s = 1+m + (1 + m)m/2 designs, an algorithm is used to weight the designs.The optimization procedure is as follows: 1. Analyze the initial design and m perturbed designs 2. Use a least square fit to determine the polynomial coefficients for the objective function Ψo (b) and each constraint Ψi (b) The procedure is shown in Fig. 1 for a problem with one design variable.For this study the implementation of the method in Altair HyperOpt has been used [4].

Stochastic analysis
Our world is characterized by non-deterministic processes.Uncertainties and non-deterministic behavior are essential in all physical processes that are known as the core applications of explicit finite element analysis, i.e. chocks, explosions and impacts.Car crash safety and occupant safety testing as parts of the impact-type problems are highly non-deterministic, e.g.repeated collision tests of the same vehicle type will always lead to different results.For this study the explicit finite element code LS-Dyna has been used [5].Simulation models are, however, deterministic.Two simulations of the same input deck will yield identical results provided they are run on the same computer under the same conditions.Thus, the computer environment fail to mirror the physical reality in this sense.In order to increase the robustness and reliability of a design with respect to physical deviations from the ideal simulation model, and to assure its functionality for loads and load cases which deviate from the design load, we must perform a large number of simulations of the design with perturbations from the ideal model.The results from such study must be evaluated from a statistical standpoint in order to make judgements about the reliability and robustness of the design.
The general idea of the Monte Carlo method is as follows.Let the vector of the output responses y i , depend on a vector of random variables x i , whose marginal probabilistic descriptions are known or possible to estimate.By using well tested sampling algorithms for the random variables, a large population of sample sets which contain a sample of each random variable can be generated.The sample sets are chosen such that the histograms for each random variable approaches its marginal probabilistic distribution.An equally large population of output variable sets y i can be obtained by using the analysis code in a deterministic way to generate a set of output responses for each set of samples of the random variables x i .Two issues are of vital importance when defining a Monte-Carlo simulation.
1.A marginal probabilistic distribution must be defined for each random variable in the random variable vector x i .2. A sampling algorithm must be defined for which a limited population of sample sets will lead to histograms for each random variable that approach the marginal probabilistic distributions.
The PROMENVIR computer environment provides a generic environment for performing Monte-Carlo analyses using any analysis code [6].It provides several choices for defining marginal probabilistic distributions such as Three sampling algorithms for choosing the samples from the distributions are available: -Latin Hyper-cube sampling -Descriptive sampling -Antithetic variates All above sampling algorithms are so-called reduced sampling algorithms that optimize the content of the population for a certain size without sacrifying the quality of the statistical description.

Objective
The objective of this study has been to perform an initial investigation if optimization performed on the result of a stochastic analysis will yield different results than a pure optimization of a single deterministic ideal model.The two cases have been compared with respect to: 1. Robustness of final design to changes in load case and/or load level; 2. Differences in design variable values and objective function of final design, i.e. the optimization result for the deterministic model.

Approach
A software solution for an optimization method using response surfaces usually provides a general environment for performing response surface optimization.The analysis input file with the substituted design variables is passed to a shell script which is responsible for all steps going from the analysis input file to an output file which contains objective and constraints (Fig. 3).When the objective and the constraints have been calculated from the analysis code, the execution is passed back to the main optimizer again.
Performing optimization on the output from a stochastic analysis adds another layer to the simulation environment.The solver call of the response surface method no longer calls the finite element solution directly directly, but starts a Monte-Carlo session.The Monte-Carlo driver executes a project that contains the complete description of the stochastic analysis to be performed on the input file of the finite element solver for the current values of the design variables (Fig. 4).
In this case, the objective and the constraints are given in terms of statistical values instead of direct simulation results, e.g.mean value of a chosen response, maximum or minimum values, standard deviations or confidence intervals.The Monte Carlo analysis provides the statistical information necessary to drive the optimization process towards a converged final design.

Applications from crash and occupant analysis 4.3.1. Example 1: Validation of Hybrid III neck model
The Hybrid III dummy which is used for assessing passenger safety during frontal vehicle collision tests must regularly undergo standardized tests.The dummy neck must for instance be tested in a so-called neck pendulum flexion test that will certify the stiffness and damping of the neck during forward bending (Fig. 5).
The test procedure is as follows: The neck and head of the dummy are mounted at the lower end of a pendulum arm.The pendulum is rotated to an initial position and released to swing from left to right in Fig. 5(A).In its vertical position, the pendulum strikes a honeycomb block which retardates the pendulum.This retardation causes the neck and head to bend forward and then rebound, Fig. 5(B).The test conditions for the neck flexion test are specified in Table 1.Validations of such simulation model have been published earlier [7,8].
The objective of Example 1 is to investigate if stochastic analysis can make the simulation model of the neck more robust considering that the test speci-

Statistical results
Objective and constr.1.
The design variables are scale factors on the initial load curves or stiffnesses.All design variables were set to b i = 1.0 for the initial design.A normal distribution was selected to describe the random variables.Limits for the distributions were set at the test limits according to Table 1.
Mean values from a Monte-Carlo analysis for the maximum head/neck moment and head angle were chosen as output responses for the optimization.The objective function minimized in the combined analysis is defined by Eq. ( 4).
The function σ(X) is the standard deviation of X.In Eq. ( 4) M is the maximum head/neck moment and A is the maximum head angle.The value n is the total number of samples in the Monte-Carlo analysis and the counter is the sample number.The superscripts S and P stand for simulation results and ideal physical results, respectively.The quantitiesare weighting factors.The optimization was constrained by upper and lower bounds on maximum moment and angle.These values were allowed to vary up to 5% from their respective ideal test results.
A population of n = 50 samples was used for the Monte-Carlo analysis to get average values for all responses.A direct minimization of the function X without Monte-Carlo analysis was also performed in order to compare the final designs with and without Monte-Carlo analysis.The final designs from the optimizations with and without Monte-Carlo analysis were finally subjected to a 150 sample Monte-Carlo analysis.The results from these analyses were evaluated and compared.
Table 2 shows the results from the investigation.The two final designs from optimization with and without Monte-Carlo simulation are compared with the physical test requirements.The following designs are compared: Design 1: Direct optimization of the function X, followed by a 150-sample Monte-Carlo simulation on the final design.
Design 2: Optimization of function Ψ MC 0 , followed by a 150-sample Monte-Carlo simulation on the final design.
The results from the Monte-Carlo optimization were generally worse compared to the direct optimization.Two statistical measures were, however, reduced: the standard deviations of the objective function and of the maximum head angle.Thus, an objective function based on the mean value of the function X in Eq. ( 5) from the samples adjusted for the standard deviation was not sufficient.In order to drive the optimization process to generally improved statistical results compared to a direct optimization of X, the statistical measures to minimize need to have a stronger direct influence on the objective function than what was used in this example.

Example 2:
Optimization and stochastic analysis of vehicle front structure Two front rails have been connected to a cross member to form a simplified vehicle front structure (Fig. 6).Both rails are clamped at their right ends and the left end.The structure is subjected to a straight impact from a rigid planar wall with the mass of 1500 kg and a velocity of 3 m/s.The rails have an initial shell thickness of 2.0 mm.The cross member has the initial shell thickness of 0.5 mm.All parts of the front structure are modeled using linear elastic-plastic material behavior.
The structure was subjected to a combined optimization and Monte-Carlo simulation.The objective of Example 2 was to investigate if a more robust design could be found by using stochastic analysis when considering changes in impact angle of the rigid wall.The following design variables and random variables were used for the combined analysis: Design variables b 1 : Shell thickness of the first part of the bend for both rails, 0.2 mm < b 1 < 3.0 mm b 2 : Shell thickness of the second part of the bend for both rails, 0.2 mm < b 2 < 3.0 mm b 3 : Shell thickness of the cross member, 0.2 mm < b 3 < 3.0 mm Random variable x 1 : The angle of the rigid wall relative to the front structure, 0 deg < x 1 < 18.9 deg A half normal distribution was selected for the random variable x 1 with bounds at 0.0 deg and 18.9 deg, respectively (Fig. 7).
The objective function to be maximized in the combined analysis is defined as The functionis σ(X) the standard deviation of X. E int is the total internal energy absorbed by the front structure at the time 40 ms.A population of n = 20 samples was used in the MC analysis to calculate the objective function.A direct optimization of the function X was also performed in order to compare the final designs with and without the use of stochastic analysis.The final designs from optimization with and without MC analysis were finally subjected to a 40 sample MC analysis where a more narrow corridor was used for the random variable.The range for the random variable The thickness of the cross member has no influence since a straight impact will lead to zero net forces in the lateral direction.Thus, comparisons were performed using three different designs.
Design 1: Final design from direct optimization of function, X.A conclusion was made that the direct optimization of the function X was not successful for this optimization setup.A comparison with the optimal design from the linear case shows only small differences for with zero impact angle.However, when the angle is increased stepwise from zero to 14.0 degrees, large differences appear as depicted in Fig. 8.
The Monte-Carlo optimization yields a much more stable design, which is less sensitive to the impact angle.The energy absorption remains constant for angles from 3.58 to 14.0 deg and is generally larger than the energy absorption for Design 3. Design 3 shows an unstable behavior where energy absorption capability changes rapidly when the impact angle changes.The reason for this behavior can be seen in Fig. 9.The design from the Monte-Carlo optimization results in a localized buckling behavior (Fig. 9 rect optimization of function X, on the other hand, collapses completely in the longitudinal/vertical direction (Fig. 9(A)).

Conclusions
Stochastic analysis using the Monte Carlo method has been combined with a response surface optimization in order to investigate if an increased robustness of the final optimized design can be achieved compared to a final design from a direct optimization.As a measure of the robustness, the sensitivity of the final design to changes in load level and load case was used.Two examples have been given where optimization and stochastic analysis were combined: one validation example and one design optimization example.
The study shows the potential of the presented method but also the difficulties.Even more care must be taken when defining objective functions in order to drive the optimization in the desired direction, and the number of samples included in the Monte-Carlo simulation must be carefully selected in order to reach convergence for the optimization.The results from the examples show these difficulties; in one case the standard deviations of important responses were reduced, and in the other a better general solution and a more robust buckling behavior was found.In both examples a completely different setup of design variables was found when using Monte-Carlo optimization.Thus, no clear trends in the results could be found between the two examples.Nevertheless, the results from the study indicate that stochastic analysis can be a helpful tool in order to increase the robustness of designs.
This study should be looked upon as an initial study of how optimization can be combined with stochastic analysis.Further investigations are necessary in order to fully understand the benefits from the usage of this tool.Nevertheless, the authors believe that this methodology, as presented here or in a modified form, will be an important tool in the future design process.
x 1 was changed to 0 deg < x 1 < 7.1 deg.Since no constraints were used for the optimization, the solution to the linear version using rigid-plastic material behavior of this problem is known.The maximum absorbed energy at any given time is reached with the maximum stiffness of the model.Thus, an optimization with linear behavior would give the design b 1 = b 2 = 3.0 mm, b 3 = 0.5 mm = initial value.

Fig. 7 .
Fig. 7. Marginal probabilistic distribution for random variable x 1 .Design 2: Final design from the optimization of function Ψ MC 0 .Design 3: Optimal design from the linear case, d 1 = d 2 = 3.0 mm, d 3 = 0.5 mm.The models were compared with respect to final objective function for the deterministic model and to the results of the final 40 sample Monte-Carlo simulation of designs 1 and 2. Comparisons were also made for deterministic simulations of the final designs 1, 2, and 3 where the hit angle from the wall was varied from zero to 14.0 deg.The optimizations performed with and without Monte-Carlo simulation show completely different results.The Monte-Carlo optimization shows a better behavior from all aspects considered in this analysis.A conclusion was made that the direct optimization of the function X was not successful for this optimization setup.A comparison with the optimal design from the linear case shows only small differences for with zero impact angle.However, when the angle is increased stepwise from zero to 14.0 degrees, large differences appear as depicted in Fig.8.The Monte-Carlo optimization yields a much more stable design, which is less sensitive to the impact angle.The energy absorption remains constant for angles from 3.58 to 14.0 deg and is generally larger than the energy absorption for Design 3. Design 3 shows an unstable behavior where energy absorption capability changes rapidly when the impact angle changes.The reason for this behavior can be seen in Fig.9.The design from the Monte-Carlo optimization results in a localized buckling behavior (Fig.9(B)).Design 3 experiences global buckling in the lateral direction and the energy absorption capability decreases as the front structure deforms (Fig.9(C)).The design from the di-

Table 1
Specifications for the neck pendulum test in flexion, the results are examples

Table 2
Comparison of results for Example 1 (Using a 150-sample Monte-Carlo simulation)