Mathematical optimization of a plate volume under a p-Laplace partial differential equation constraint by using standard software

A main aspect in the design of passenger cars with respect to pedestrian safety is the energy absorption capability of the engine hood. Besides that, the hood has to fulfill several other requirements. That makes it necessary to develop easy and fast to solve prediction models with little loss in accuracy for optimization purpose. Current simulation tools combined with standard optimization software are not well suited to deal with the above mentioned needs. The present paper shows the application of mathematical methods on a simplified self developed model to reduce the optimization effort. A linear and a nonlinear model are introduced and a way for solving both is pointed out. Finally it is shown, that it is possible to simplify models and get optimization results much faster by using mathematical theory. Such results can be used in support of the original problem or as an input to space mapping based optimization algorithms, such as surrogate optimization.


INTRODUCTION
Within the broad area of simulation in vehicle development there is the special field of pedestrian safety which has gained importance due to ongoing introduction of aggravated legal requirements.Therefore vehicle manufacturers are obligated to produce cars which are less harmful to pedestrians in case of collision.One of the main aspects in that area is the engine hood.The hood is one of the most important contact points of a pedestrian's head with the car during a collision with a speed up to approximately 40 km/h.That is why it is necessary to change the scenario of the head impact on the hood.There are two possibilities to take this into account.One thing is to change passive protection (geometry, stiffness, material, etc.) and the other thing is to use active protection systems (actuators that move the hood before the impact).The main focus in this paper is how to model, solve and optimize a hood-like problem in an analytical/mathematical way.This model will not include as much physics as some Finite-Element-Method (FEM) software packages but it will allow a first step to use analytical methods.Furthermore a special optimization method ("surrogate optimization") will be applicable, which is the topic of current research (see [1]).
In the following chapter 2 the physical problem is shown and a nonlinear and a linear model, represented by partial differential equations (PDEs), are introduced.The optimization problem is stated and it is shown how the gradient information for the linear problem can be derived (chapter 3).
Finally, chapter 4 depicts how a solution can be achieved for the nonlinear problem by using a fixed point iteration.Furthermore a typical optimization run and the results of a parameter study are presented.

PHYSICAL PROBLEM
For simplification a plane plate representing a vehicle hood is used (Figure 1).The boundaries of the plate are constrained with zero degrees of freedom (DOF) and each inner point has one degree of freedom, which is the movement in z-direction.Notation for Figure 1:

THE NONLINEAR AND THE LINEAR MODEL
A variational principle valid within the framework of deformation theory of plasticity (Hencky, Ilyushin and Nadai; see [2]) yields a nonlinear governing Euler-Lagrange equation for the deflection of the plane plate, shown in Figure 1. (1)

With:
Hollomon coefficient (constant material property, typically n = 0,15...0,35) plate thickness at each point of the plate Ω z-displacement of each point of the plate load at each point of the plate The solving of the nonlinear PDE in Eqn (1) is shown in chapter 4.1 by using a fixed point iteration.For easier handling and solving reason Eqn (1) was changed into a linear PDE which can be formulated as follows: (2) u = 0 on Γ This equation forms the basis to state the optimization problem in the next chapter.

OPTIMIZATION PROBLEM AND IMPLEMENTATION
We are interested to have an optimization solution for Eqn.(2) to make a first move towards the original problem mentioned in chapter 2.1.The general base for numerical optimization can be found in various literatures (e.g.[3], [4], [5] and [6]).Whereas [7] gives a brief overview on optimization in structural mechanics from an engineering point of view.In this paper we have to additionally consider the theory of optimal control of PDEs.For a general reference refer to [8], for application to [9], [10], [11], [12] and [13].
One of the main aspects for car manufacturers is the factor of weight.The weight of a component is directly proportional to the volume, which is the objective of the following optimization.

STATING THE OPTIMIZATION PROBLEM
(3) lower an upper bound fpr the plate thickness upper bound for the z-displacement In this subsection we formulate the optimization problem and derive first order necessary optimality conditions.For more details refer to [1].
The Lagrange multiplier for the state constraints is a measure, therefore we penalize the state restriction on u.This turns Eqn.(3) into the following optimization problem with an augmented objective functional.To deduce the gradient for the objective functional of Eqn.(4) it is necessary to express the Lagrange function as follows.For better readability the arguments (x, y) are omitted.
(5) which leads us to the directional derivative (6) for δu: Ω → R and' • 'the scalar product of two vectors.
By the optimality condition we are able to derive the adjoint equation ( 8) Furthermore we need another directional derivative (9) for This leads to the variational inequality (10) Now we can state the optimality system which consists of Eqns.( 2), ( 8) and (10).

Optimality system
If (t * , u * ) is a solution for Eqn. ( 4), then there is a p * that satisfies the following relations: state equation: variational inequality: Thus we have the gradient of the extended objective functional (4): Therefore we have to solve the state equation (10) first to calculate the gradient of u * .Afterwards the adjoint equation has to be solved (under respect of the actual solution u * ) to gain the solution for p * .With both solutions it is possible to obtain the wanted solution for the gradient in Eqn.(14).

IMPLEMENTATION DETAILS
The optimality system in Eqn. ( 10) is implemented and solved in Matlab combined with COMSOL.Matlab is used as a general development environment whereas COMSOL serves as PDE solving tool.
A trust region algorithm "dogleg" was chosen for optimization purpose.This algorithm is a part of the Matlab optimization toolbox.Trust region algorithms are well suited to deal with nonlinear medium scale problems.An overview of these algorithms can be found in [14] and [6].
The required algorithm is called by the Matlab function fmincon.Standard settings of this function did not provide appropriate results.It was possible to obtain better results by using modified options (GradObj: on, DerivativeCheck: off, LargeScale: off, MaxIter: off, Display: iter, Diagnostics: on).
Typical values for the weight factors were chosen with k = 1 and α = 100.000.COMSOL as well as the Matlab function fmincon offer the possibility to numerically calculate the gradient information.A self programmed central difference scheme did not make much difference to the result but helped understanding the behavior of the gradients of u * and p * .

SOLVING THE NONLINEAR PROBLEM BY FIXED POINT ITERATION
The solution of Eqn. ( 2) can be calculated by COMSOL, but solving the original Eqn.(1) did not work the same way.
A first approach is the Banach fixed point theorem which leads to a fixed point iteration (see both in [15]) to use the solution of Eqn.(2) on solving Eqn.(1).The necessary iteration rule is stated as follows: (15) with u k + 1 = 0 at Γ.
Notice the slight change (additional ε-term) which was made from Eqn. (2) to Eqn. (15).This is necessary to avoid singularities during the PDE solving for gradients equal to zero.The implementation of the iteration rule leads to the following results.Using the parameters n = 0, 22, f = 5N (const.)and a randomized thickness distribution (Figure 2) leads to the convergence of the fixed point iteration according to Table 1.
Ω, κ = 0 ( ) Mathematical optimization of a plate volume under a p-Laplace partial differential equation constraint by using standard software Figure 2 Randomized thickness distribution.Table 1 shows, that there is quite some relative change in the maximal displacement (u k max ) until iteration step 17.From iteration step 18 to iteration step 40 the absolute value of the relative change in the maximal displacement is smaller than 1e-6.That shows, that accuracy is not much gaining during 22 iterations.Hence, 20 iteration steps seem to be accurate enough to obtain a solution for Eqn.(1)

RESULTS FOR A TYPICAL OPTIMIZATION RUN
Consider the optimization problem in Eqn. ( 4) for which the equations in the optimality system of Eqn. ( 14) have to be solved.The available parameters are set as follows: t = 0,2 mm, t -=10 mm, u -= 0,045 mm, α =10 6 , κ = 1, n = 0, 22, f = 5N (const.)and a starting thickness of 1,75mm.With these parameters the optimization process generates the iteration steps in Figure 3.There you can see the change in the thickness for discrete FE nodes.Blue denotes small and red large thickness values.For the optimization run it was necessary to solve the PDE from Eqn. (2) 316 times.For comparison, consider, that an optimization without adjoint equation (Eqn.( 14)) and analytic gradient information takes more evalutations.In this example, the optimization with a numerically differentiated gradient needs 5461 evaluations of Eqn.(2).The advantage of using mathematical methods to derive an analytical gradient results in an 17 times lower PDE evaluation effort for an optimization result.

PARAMETER STUDIES
The influence of different parameters on the optimization solution and optimization convergence was examined by a parameter study.Parameters are described in the following part, whereas values and/or settings for the parameters are shown in Table 2.

•
Norm of the gradient: Because of the discrete scalar product it is necessary to define a norm.One possibility is to use an Euclidean norm, another one to use a h-2-Norm (weighted Euclidean norm) which can be stated as follows: With υ an arbitrary vector of length N and T the area of the entire plate.Additionally N depicts the degree of freedom of the finite element discretization.

•
The plate thickness starting value in each FE-node (t start ) • Hollomon coefficient n: A material parameter used in the material law of Ludwik-Hollomon [16].It describes the relation between effective stress and effective strain.
With σ effective stress, κ and n material constants and ϕ the effective strain.

•
Weighting factor κ: Weighting of the volume in the objective function • Weighting factor α: Weighting of the penalty term in the objective function : , • Gradient norm: Both norms differ solely by a scaling factor, which explains why the optimization results do not qualitatively differ.Accurate results are only derived with the h-2-norm, but the directional vectors differ in terms of length only.An implementation of a stepsize rule (e.g.Armijo stepsize rule [17]) would speed up the convergence of the optimization.

•
The plate thickness starting value: Values, far beyond the optimum value led to a slower convergence of the optimization but resulted in the same optimal thickness parameters.The convergence rate could also be increased by the above mentioned stepsize rule.There was no other observed influence on the optimization result.

•
Hollomon coefficient n: Variation of this parameter showed only small effects on the result of the state equation.This is can be explained by Eqns.( 11) and ( 17).

•
Weighting factor α: Because of the small value of the restricted maximal displacement (0,045mm) this factor was important for the weighting of the penalty term in Eqn. ( 4).It has to set the penalty term to a reasonable level, which allows the objective function to consider the maximal displacement and to don't exceed it.A value of 1e6 showed good results.

•
Weighting factor κ : This factor, combined with α , accounts for the order of magnitude of the objective function.So it has a similar role as α.It (κ ) was set to 1 to achieve the effective plate volume.Hence there is just a left for setting the order of magnitude correctly.

CONCLUSION
The physically motivated problem was to optimize the thickness distribution of a vehicle hood.It was mathematically described and stated as a linear and a nonlinear model.Using Lagrange functions it was possible to derive a solution for the gradient of the stated optimization problem.This caused a significant decrease of function evaluations.By using a fixed point iteration it was shown, that the solution of the nonlinear PDE (p-Laplace equation) can be computed.Parameter studies showed the influence of different parameters on the optimization with respect to the optimization result and convergence.
Ongoing and future research is to incorporate the models from this paper into the space mapping technique and surrogate models (see [18], [19] and [20]).This optimization method is promising to be used with sophisticated models and time-consuming simulations.First results will be published in [1].Moreover, some efforts will be done to further enhance the solving of the PDEs (Eqns.( 1) and ( 2)) by a multigrid method according to [21] to improve calculation time.
-...lower and upper bound for the plate thickness u -...upper bound for the z-displacement κ, a...positive weight functions for the volume and the penalty term

Table 1
Convergence of the Fixed Point Iteration

Table 2
Values and settings for the parameter study