Global Optimization Method to Multiple Local Optimals with the Surface Approximation Methodology and Its Application for Industry Problems

Abstract

Although generally speaking, a great number of functional evaluations may be required until convergence, it can be solved by using neural network effectively. Here, techniques to search the region of interest containing the global optimal design selected by random seeds is investigated. Also techniques for finding more accurate approximation using Holographic Neural Network (HNN) improved by using penalty function for generalized inverse matrix is investigated. Furthermore, the mapping method of extrapolation is proposed to make the technique available to general application in structural optimization. Application examples show that HNN may be expected as potential activate and feasible surface functions in response surface methodology than the polynomials in function approximations. Finally, the real design examples of a vehicle performance such as idling vibration, booming noise, vehicle component crash worthiness and combination problem between vehicle crashworthiness and restraint device performance at the head-on collision are used to show the effectiveness of the proposed method.

1. Introduction

Due to the recent developments in computer science and software technology, fluid analysis, nonlinear structural analysis and their combined analysis have become possible for large-scale models and will eventually be employed in manufacturing design. The importance of being able to make rational and correct decisions during the early stages of design is currently being emphasized. However, non-linear problems such as crashworthiness analysis require considerable time for one functional evaluation even using supercomputers. Optimization techniques which cover several design disciplines (multiple discipline optimization-MDO [1]) have been hot topics recently. Sensitivity-based optimization approaches [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], which sequentially construct local approximations, are widely employed and provide valuable contributions to linear optimization problems. Unfortunately, sensitivity approaches are not available for nonlinear optimization problems. For this purpose, response surface methodology (RSM) works effectively to achieve these design optimizations. The most important research work on RSM is to considerably reduce the extremely numerous number of function evaluations for computational cost problem and to get the robust approximation design model. That means instead of the time-consuming finite element method (FEM) analysis, robust approximation models can be constructed to evaluate responses rapidly in the design space of interest by RSM, when the response surface of the structure is infiltrated with computational noise. Early examples are such that White Jr. K.P. used RSM in their study of passenger-car crash worthiness [16] and Giunta, A.A. applied RSM to high speed civil transportation design [17]. Several key issues must be considered to develop the response surface modeling. One issue for RSM may be the selection of the activation function. Most response surfaces are expressed in terms of polynomials and among them quadratic polynomial are the most popular [18, 19, 20, 21, 22, 23, 24, 25]. Unfortunately, this is not always accepted as a proper activation function for response surfaces and may not provide a satisfactory solution, particularly for wide ranges of design space with heavy non linearity. This limits the validity of RSM to only cover a part of the design space. And so in the past, feed forward multiple layers back-propagation neural network seemed to be a potential activate function in response surface function methodology [26]. Carpenter, W. made the comparison between polynomials and multiple layers neural network [27]. There are some problems when using back-propagation neural network: (1) the learned result is influenced by the initial parameters of the network due to the nonlinear gradient descent algorithm. (2) Training of neural network costs computational time, etc. In contrast, Hagiwara et al. firstly employed the Holographic neural network (HNN) [28] for the activation function of the response surface and compared it to the multilayer feed-forward neural network and demonstrated that the former type of neural networks had a much higher approximation accuracy and the training cost of HNN was about 1/1000 the computational time compared to that of a multilayer feed-forward neural network in the case on the assumption that the designer roughly knows the location of local optimal [29, 30].

Another problem in using RSM is the difficulty in predicting the error incurred using response surface function values rather than the “true” values from the analysis code. Therefore, the extensive RSM research is focused on the selection of design trial points in design space, such as Box–Behnken Designs (BBD [20, 21, 25]) and Central Composite Design (CCD [19, 22, 23, 24]) which are two typical choices. And Central Composite Rotatable Design (CCRD [31]) and Face Central Composite Design (FCCD [18]) have also been applied to optimization studies in recent years.

By the way, now they say 3rd AI boom has come due to the born of deep learning [32]. And so recently the use of artificial neural network models to derive response surface plot has found increasing use in optimization modeling. In most of them, CCD is used such as in [33, 34, 35]. These researches mainly consist of 2 types. One is neural network itself works as RSM step [33]. The other is artificial neural network function is used as an activation function of RSM [34, 35]. Deep learning must have a lot of data. It also needs a lot of data to search the global optimization for multiple local optimal problem shown here by RSM.

Since the exact form of response surface remains unclear in general application and the accuracy of the function approximation near the global optimal design is very important, techniques to search the design space containing the global optimal design are examined. For the purpose of determining the global optimal design, techniques for searching the most-likely global optimal design (MPOD) over the entire design space and techniques for determining a more accurate approximation near the global optimal designs using the extrapolation ability of the HNN are proposed [36, 37, 38].

The MPOD is as mentioned above a sequential procedure composed of two steps. The first step is employed to determine the most likely global optimal design over the entire design space by conditioned random seeds which control the mutual distances between every two designs, to scan the wide design space rapidly. The second step is to verify and improve the accuracy of the approximation of the objective function and constraint functions at optimal design. The trade-off between the approximation accuracy and computational cost can be well balanced by the preset threshold of the mutual distance of designs. Some design examples for vehicle structure such as idling vibration, interior noise, crashworthiness and combination problem between vehicle crashworthiness and restraint device performance at the head-on collision demonstrate the validity and utility of this method. The method developed by us can omit many useless data and very unique from both of the 2 important issues for RSM such as activation function and selected design points. This approach is also useful for the execution of the deep learning. Let us get started to show this splendid method.

2. Theory of holographic neural network (HNN)

2.1 Basic theory of HNN

Here, x = {x₁, x₂, …, x_k}^T is the stimulus data, y = {y₁, y₂, …, y_m}^T is the training data. If there is learning data for n, it is presented in Eq. (1).

$$X=\left(\begin{array}{c}{X}_1^T\\ ⋮\\ X_n^T\end{array}\right)=\left(\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1k}\\ x_{21}& x_{22}& \cdots & x_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{\mathit{nk}}\end{array}\right)Y=\left(\begin{array}{c}{Y}_1^T\\ ⋮\\ Y_n^T\end{array}\right)=\left(\begin{array}{cccc}{y}_{11}& y_{12}& \cdots & y_{1m}\\ y_{21}& y_{22}& \cdots & y_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ y_{n1}& y_{n2}& \cdots & y_{\mathit{nm}}\end{array}\right)$$E1

Here, X, Y are input and output matrices. Each element of X and Y is converted to angles ${\theta }_{\mathit{ai}}\left(a=1,\cdots ,n;i=1,\cdots ,k\right)$ and ${\varphi }_{\mathit{aj}}\left(a=1,\cdots ,n;j=1,\cdots ,m\right)$.

$${\theta }_{\mathit{ai}}=f_x\left(x_{\mathit{ai}}\right){\varphi }_{\mathit{aj}}=f_y\left(y_{\mathit{aj}}\right)$$E2

Linear, Sigmoid and Arc Tangent functions are applied for f_x and f_y.

Especially, for the stimulus data, Sigmoid function is used as follows.

$${\theta }_{\mathit{ai}}=\frac{2\pi }{\left(1+e^{-c_{\mathit{ai}}\left(x_{\mathit{ai}}-{\mu }_{\mathit{ai}}\right)/{\mathrm{ϭ}}_{\mathit{ai}}}\right)}$$E3

Here, ${\mu }_{\mathit{ai}}\text{and}{ϭ}_{\mathit{ai}}$ are the mean value and standard deviation of x_ai respectively and ${\lambda }_{\mathit{ai}}$ defines a vector magnitude bound within the unit circle (0 to 1), expressing a weighting or dominance for each element of the stimulus field. If no design variable dominates, ${\lambda }_{\mathit{ai}}=1$. $c_{\mathit{ai}}$ is the coefficient to adjust the transformation slope shown in Figure 1.

Figure 1.

Transform function of phase angle: If the untrained designs are outside the trained design space, the predicted values are same values as those at trained design edges 0 or 1 for the case of $c_{\mathit{ai}}=1.0$. In contrast, the predicted function values are extended outside the trained design edges 0 or 1 for the case of $c_{\mathit{ai}}=0.5.$ In this way, the extrapolation can be done and the trained designs are not necessary on the boundary.

Next, the above angles are mapped on complex representation by using exponential function as follows.

$$s_{\mathit{ai}}={\lambda }_{\mathit{ai}}{e}^{\widehat{i}{\theta }_{\mathit{ai}}}{r}_{\mathit{aj}}={\gamma }_{\mathit{aj}}{e}^{\widehat{i}{\varphi }_{\mathit{aj}}}$$E4

From the above, X and Y are converted to stimulus S and response R.

$$S=\left(\begin{array}{cccc}{s}_{11}& s_{12}& \cdots & s_{1k}\\ s_{21}& s_{22}& \cdots & s_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ s_{n1}& s_{n2}& \cdots & s_{\mathit{nk}}\end{array}\right)R=\left(\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1m}\\ r_{21}& r_{22}& \cdots & r_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ r_{n1}& r_{n2}& \cdots & r_{\mathit{nm}}\end{array}\right)$$E5

H = [h₁,…, h_m], the transfer function of HNN is presented by minimizing the difference between $\text{teaching data}R\text{and}S・H\mathrm{as}\text{follows};$

$$\mathrm{H}={\left(S^{\ast }S\right)}^{-1}{S}^{\ast }R$$E6

Here, symbol * presents conjugate transpose. Output V can be predicted to be V = U・H for new input U. H is the transfer function of HNN. It is not efficient to treat Eq. (6) because it has the inverse matrix operating. When the number of activation functions is large, it will lead to considerable calculation time and may result in calculation error. The following iteration process is used to save calculation time.

$$H_0=S^{\ast }\mathrm{Y}$$E7

The enhancement of parameters will be performed by

$$∆H_{k+1}=S^{\ast }\left(Y-S{H}_k\right)$$E8

$$H_{k+1}=H_k+∆H_{k+1}$$E9

The iteration may be terminated when the error between the approximating response and the true response is orthogonal to every column vector of S in a complex domain. The iteration is terminated when the error becomes smaller than the preset threshold value, eps.

$$\mathrm{J}\left(H_{k+1}\right)={\left(\mathrm{Y}-\mathrm{S}{H}_{k+1}\right)}^{\ast }\left(\mathrm{Y}-\mathrm{S}{H}_{k+1}\right)\le \mathrm{eps}$$E10

By the way, within the HNN, limitations of mapping accuracy can be largely overcome through a pre-processing operation involving the generation of higher-order product terms. For example, for two design variables (n = 2), s = $⌈{\lambda }_1{e}^{i{\theta }_1},{\lambda }_2{e}^{i{\theta }_2}⌉$ can be increased to 12 terms (n = 12) if performed to the second order.

$$\overline{S}=⌈{\lambda }_1{e}^{i{\theta }_1},{\lambda }_2{e}^{i{\theta }_2},{\lambda }_1{e}^{-i{\theta }_1},{\lambda }_2{e}^{-i{\theta }_2},⌉$$

$${{\lambda }^2}_1{e}^{i2{\theta }_1},{{\lambda }^2}_2{e}^{i2{\theta }_2},{{\lambda }^2}_1{e}^{-i2{\theta }_1},{{\lambda }^2}_2{e}^{-i2{\theta }_2},{\lambda }_1{\lambda }_2{e}^{i\left({\theta }_1+{\theta }_2\right)},$$

$${\lambda }_1{\lambda }_2{e}^{i\left(-{\theta }_1+{\theta }_2\right)},{\lambda }_1{\lambda }_2{e}^{-i\left({\theta }_1+{\theta }_2\right)},{\lambda }_1{\lambda }_2{e}^{i\left({\theta }_1-{\theta }_2\right)}$$E11

Furthermore, the extrapolation of the approximated response function may be easily performed with trained HNN by presetting the coefficient $c_k$ in Eq. (3) to a small value to adjust the slope to be a smaller. Figure 1 shows the effect of the coefficient. In Figure 1, the experimental designs used for training the neural network are located within the range [0–1], the dotted line expresses the transform when $c_k=1.0$ and the solid line expresses the transform when $c_k$ = 0.5. We can observe from the figure that if the predicted designs are outside the trained design space, then the predicted function values outside the trained design space may be considered to be the same values as those at trained design edges 0 or 1 for the case of $c_k=1.0$. In contrast, the predicted function values outside the trained design space may keep the trend beyond the trained design edges 0 or 1 for the case of $c_k=0.5.$Figure 2 shows the effect of extrapolation. In this figure, space “I” expresses the design space inside which designs are used for approximation by training and space “D” expresses the entire feasible design space. Approximation using designs in space “I” may extrapolate the values of space “D”. This is an important issue for reducing functional evaluations and some advantages of this are: 1) design spaces with irregular shapes can be considered, 2) any number of experimental designs can be considered, and 3) searching for reasonable optimal designs over a wide space becomes possible.

Figure 2.

Explanation of extrapolation: Space “I” expresses the design space inside which designs are used for approximation by training and space “D” expresses the entire feasible design space. Approximation using designs in space “I” may extrapolate the values of space “D”.

2.2 Examination for getting the inverse matrix of transfer matrix

Although for calculation inverse matrix (${\overline{S}}^{\ast }\overline{S}$), we can use the procedure from Eqs. (7) to (10), for not so big problem, the inverse of (${\overline{S}}^{\ast }\overline{S}$) is calculated directly. If (${\overline{S}}^{\ast }\overline{S}$) has not the inverse matrix in Eq. (6), we can use Moore Penrose generalized inverse matrix ${\widehat{H}}_{\mathit{MP}}$ which is generated by operating ${‖h_k‖}^2$ (k∼1,….,m) to be the minimum. Here $h_k$ is the column vector of H.

$${\widehat{H}}_{\mathit{MP}}={\overline{S}}^{\ast }\overline{S}{\left({\overline{S}}^{\ast }\overline{S}{\overline{S}}^{\ast }\overline{S}\right)}_G^{-1}{\overline{S}}^{\ast }R$$E12

Here ${\left(\right)}_G^{-1}$ is the generalized inverse matrix.

But following is the example where it has not satisfied predicted value for Figure 3 with ${\widehat{H}}_{\mathit{MP}}.$

Figure 3.

Visualization of function used for test.

It is shown the approximation function in Figure 4 which is generated with Moore Penrose pseudo-inverse matrix by using 23 sampling points as shown ● in Figure 5. Here the base function of HNN is extended to the 3rd order where the total number of the base function is 30.

Figure 4.

Approximated function using Moore-Penrose generalized inverse: The base function of HNN is extended to the 3rd order. The surface shape is quite different from the original one because over fitting phenomena occurs.

Figure 5.

Layout of sampling point on design space.

It is recognized that the surface shape in Figure 4 is quite different from the one in Figure 3. This is because over fitting phenomena occurs which is unique to neural network. This phenomena occurs in such way that the precision of the predicted values of the points other than sampling points become worth because the values of the sampling points are forced fit to the given HNN function. In this case we cannot use HNN for the prediction.

Instead, the following penalty function was developed [39].

$$f\left(h_k\right)=\left(r_k-\right.S{h}_k){}^{\ast }\left(r_k-S{h}_k\right)+{\rho }_1{h}_k^{\ast }O{h}_k+{\rho }_2{h}_k^{\ast }{h}_k\left(k=1,\dots m\right)$$E13

Here O is the square matrix whose all elements are 1. The right side paragraph 1 means approximation error to output, the paragraph 2 means the average of regression factor h_k and the paragraph 3 means variation. The paragraphs 2 and 3 are weighted with ρ₁, ρ₂. In Eq. (13), it cannot be predicted uniquely without some constraints and so such conditions are added. Although theoretically, H is predicted uniquely without the paragraph 2, it is more stable with this paragraph. H is as follows by minimizing the penalty function (13);

$$H_P={\left(S^{\ast }·S+{\rho }_{1}O+{\rho }_{2}I\right)}^{-1}·S^{\ast }·R$$E14

Here, I is identity matrix. The approximation function by using this H_P is depicted in Figure 6 where 23 sampling data are used as above, ${\rho }_1,{\rho }_2$, are 1.0. It can be said it does not occur the over fitting phenomena in Figure 6. The approximation functions by both methods are estimated by mean square error.

$$E=\sqrt{\frac{1}{N}{\sum }_{k=1}^N{\left(y_k-{\widehat{y}}_k\right)}^2}$$E15

Figure 6.

Approximated function using penalty coefficient method: The over fitting phenomena does not occur.

Here N is number of estimation points, $y_k$ is the real function value and ${\widehat{\mathrm{y}}}_{\mathrm{k}}$ is the value of approximation function. Here it is called “approximation error” which is estimated by using 23 sampling points and called “prediction error” which is estimated by 2601 points, generated by dividing design region into each of the two axes 50 equally separated. From Table 1, it can be estimated in the case of Moore Penrose, the approximation error is very small but the prediction error gets worth. On the other hand, in the case of proposed method, the approximation error is not smaller than in the case of Moore Penrose but the prediction error is moderately suppressed. In the case of proposed one, it is important how much ${\rho }_1,{\rho }_2\text{to}\mathrm{be}\mathrm{set}.$ If they are very small, it approaches to the approximation function by Moore Penrose because the part of right side paragraph 1 of Eq. (14) is larger. On the other hand, if they are too large, the part of right side paragraph 1 of Eq. (14) is too small to be used effectively the information of sampling data. And so when the number of sampling data is large and the sampling data are selected equally over the design region, ${\rho }_1,{\rho }_2$ are small so that the information of sampling data can be used effectively. On the other hand, if it is difficult to have enough sampling data, it is better ${\rho }_1,{\rho }_2$ are made reasonably large so that design region is roughly approximated by small information in order not to have over fitting phenomena.

	${\widehat{{H}}}_{{MP}}$	${\widehat{{H}}}_{{P}}$
Error in learning	0.0447	0.4993
Error in prediction	1.8681	0.9373

Table 1.

Comparison between two methods by square-mean-error.

In the case of Moore Penrose, the approximation error is very small but the prediction error gets worth. In the case of Penalty Coefficient Method, the approximation error is not smaller than in the case of Moore Penrose but the prediction error is moderately suppressed.

3. Optimization method with holographic neural network

3.1 Procedure for MPOD

A two-step sequential optimization is proposed to determine the optimal design as shown in Figure 7. The first step is to scan the entire design space to determine where the global optimum design exists, most probably by designs selected using conditioned random seeds. The second step is to increase the accuracy of the most-likely global optimum design by imputing additional information where required. The details regarding these two steps are described below. The optimum problem is defined as.

$$Maxf\left(x\right)$$E16

$$\text{Subject to}{g}_j\left(x\right)\le 0,\left(j=1,2,3,\dots .,\mathrm{m}\right)$$E17

Figure 7.

Flowchart of MPOD: The first step is to scan the entire design space to determine where the global optimum design exists, most probably by designs selected using conditioned random seeds. The second step is to increase the accuracy of the most-likely global optimum design by imputing additional information where required.

Here m is the number of constraints.

Step 1: To search the most-likely global optimum design over the entire design space. Selecting initial trial designs x_k (k = 1, 2, …., N) where N is the initial number of design points. Although one should understand that this is not the only choice, it is generally recommended N = (n+1)(n+2)/2, n is the number of design variables. Actually the determination of the N value depends on the scale of the problem such as the number of design variables and the computation time of each functional evaluation. For instance, as far as the formulation of CCD (Center Composite Design), points N=2n+1 is recommended for problems with a large number of designs. It should be emphasized that the designs produced by conditioned random seeds are employed to search the entire design space. The uniformly distributed conditioned random seeds are used to provide a good valance of search in the design space. Design conditions are required to satisfy Eq. (18).

$$d_m=min\left\{\frac{‖x_i-x_j‖}{‖x_{\mathit{CG}}‖}\right\}\left(1\le i,j\le N,i\ne j,m=1,2,\dots \dots ,p\right)$$E18

Where $x_{\mathit{CG}}$ is the normalization parameter which can be selected to be the geometrical center of design $x_k$ = $\left\{x_{1,k},x_{2,k},\dots .,x_{n,k}\right\}$, (k=1,……,N), $d_m\ge d_{\mathit{min}}$ (m=1,2,….,p), $d_{\mathit{min}}$ is the threshold of distance, p is the total number of design combinations p = ${{}_{N}C}_2$. If no new design is obtained after 5000 random trials within the threshold distance, the new threshold distance must be decreased by $d_{\mathit{min}}=\mathrm{\beta }{d}_{\mathit{min}}$ (0 $<\beta <$ 1). Designs satisfying Eq. (18) and their objective function f($x_i)\left(i=1,\dots \dots ,N\right)$. and constraint $g_j\left(x_i\right)$ (i=1,……,N, j=1,……,m) values are trained using neural networks to approximate the response surface $\overline{f}$ (x) and $\overline{g}{}_j\left(x\right).$ A global optimization design $x_{\mathit{op}}\left(0\right)$ with an approximation function can be solved by optimization codes. Then, more S designs are created with the same conditions $d_m\ge d_{\mathit{min}}$ (m=1,2,….,p),which are used to improve the accuracy of the approximation. Here, S is generally set to n. But it is also not the only choice. The random S designs are selected as the Gaussian normal distribution with the mean value of $x_{\mathit{op}}\left(0\right)$ and the deviation must be carefully chosen. The approximation accuracy can be judged by the change in the movement of approximation optimal design location. Therefore, the deviation should change with the movement distance σ=αΔd. Here, Δd is the distance of designs calculated by Eq. (19) and α is the coefficient which controls the deviation value and affects the approximation accuracy and the number of functional evaluations, as discussed below. The kth sequential improved approximation functions by the neural network is to be used to obtain new optimal design $x_{\mathit{op}}\left(k\right)$. The convergent conformation for the kth sequential approximation is taken to terminate the random seeds operation

$$\mathit{\Delta d}=\frac{‖x_{\mathit{op}}\left(\mathrm{k}\right)-x_{\mathit{op}}\left(\mathrm{k}-1\right)‖}{‖x_{\mathit{CG}}‖}\le {\mathrm{\varepsilon }}_1$$E19

Thus far, the above sequential approximation is used to determine the approximation global optimal location.

Step 2: To verify and improve the accuracy of optimal design and its performance. The accuracy check is based on

$$\frac{‖f(x_{\mathit{OP}}\left(\mathrm{k}\right)-\overline{f}(x_{\mathit{OP}}\left(\mathrm{k}\right)‖}{‖f\left(x_{\mathit{OP}}\left(\mathrm{k}\right)\right)‖}\le {\mathrm{\varepsilon }}_2$$E20

where $x_{\mathit{OP}}\left(k\right)$ is the optimal design in Step 1, ${\varepsilon }_2$ is the threshold of the termination and $f\left(x_{\mathit{OP}}\left(\mathrm{k}\right)\right)\text{and}\overline{f}\left(x_{\mathit{OP}}\left(\mathrm{k}\right)\right)$are the corresponding exact and approximation function values. If Eq. (20) is not satisfied, n learning points are added by normal distribution random number of which center is $x_{\mathit{OP}}\left(\mathrm{k}\right)$ until Eq. (20) is satisfied. Here n is number of design values and variance ${\sigma }_i$ is decided in Eq. (21).

$${\sigma }_i={\alpha }_i\left(x_{\mathit{OP}}\left(k\right)-x_{\mathit{OP}}\left(k-1\right)\right)\left(i=\mathit{1},\mathit{2},\dots \dots \dots ,n\right)$$E21

3.2 Necessary designs used for optimal design by MPOD in step 1

The convergence and efficiency of the proposed method are based on the deviation of the Gaussian random seeds. The probability of determining the global optimal design is increased when a larger deviation is used which can search a wider feasible design space resulting in an increase of the functional evaluations. In contrast, if the deviation is small, the number of functional evaluations will decrease but the possibility of searching the smaller cluster will reduce. As illustrated in Figure 8, the scatter of design points is adjusted by the deviation value. This method is sufficiently flexible to satisfy any user’s requirements through the adjustment of coefficient α. The necessary designs used for solving the global optimal can definitely be estimated. For sake of convergence and without a loss of generality, a one-dimensional example is employed to express the estimator as shown in Figure 9. Assume the design space $\mathcal{R}$ is d₀ and the size of the entire feasible design space is D.

$$d_0=\frac{2\mathit{min}‖x_e-x_c‖}{‖x_c‖}$$E22

$$f\left(x_i\right)>\mathrm{f}\left(x_j\right);x_i\in \mathcal{R};x_j\notin \mathcal{R}$$E23

We define the space $\mathcal{R}$ as an optimal cluster, which is mathematically defined by Eq. (22).which indicates that the function values within the optimal cluster $\mathcal{R}$ as an optimal cluster are longer than those outside the optimal cluster. x_e is the edge point of the optimal cluster. x_c is the geometrical center of optimal cluster and d₀ expresses the size of the optimal cluster. If the size of the optimal cluster is known, it is estimated that the maximum number of designs for determining the optimum cluster with probability unity is decided by Eq. (24).

Figure 8.

Random data distribution according to standard deviation: The scatter of design points is adjusted by the deviation value.

Figure 9.

Explanation of necessary designs to find global optimal: The design space $\mathcal{R}$ is d₀ and the size of the entire feasible design space is D.

$$P=int\left(\frac{2\mathit{max}‖x_E-x_{\mathit{CG}}‖/‖x_{\mathit{CG}}‖}{d_0}\right)$$E24

Where $x_E$ expresses the group points on the edge of design space, $x_{\mathit{CG}}$ is the geometrical center of the design space and int. expresses the rounded integer. In this case, if ${\mathrm{\epsilon }}_1$ = $d_0$ in Step 1, it can be moved to Step 2 with less than P designs.

4. Numerical examples

4.1 Comparison of approximation by holographic neural network and polynomials

The comparison of the HNN with polynomials for the function approximation is performed by the following multiple peaks function, which is shown in Figure 9.

$$f\left(\mathrm{x}\right)=10x^5\left(x^2-2\right){\mathrm{e}}^{-x^2},-4.0\le \mathrm{x}\le 3.0$$E25

For the above numerical example, the number of training designs is 15 and they are evenly spaced over the design space. The approximation accuracy of the HNN and polynomials is shown in Figure 10. In this figure, the horizontal axis expresses the orders of polynomials (a) or orders of the HNN (b), the vertical axis expresses the standard variance calculated with the 200 designs evenly spaced in the design space [−4.0,3.0]. We can see from the results that the HNN is much more accurate than the polynomials and maintains good accuracy after the ninth order of expansion. The standard deviation of the HNN at the ninth order is 0.3, while those of polynomials at the 14th order is higher than 1.0. In other words, a larger number of coefficients must be selected for the polynomials and thus, more experimental designs must be used compared to that of the HNN. Therefore, the HNN is a more suitable activation response function for this problem than polynomials.

Figure 10.

Performance comparison of polynomial and holographic neural networks: The horizontal axis is the orders of polynomials (left) or number of expansions for the HNN (right), the vertical axis is the standard variance calculated with the 200 designs evenly spaced in the design space [−4.0,3.0].

4.2 Numerical example with holographic neural network for extrapolation application

To demonstrate the efficiency of this method, a numerical example of multiple local maximums with constraints is used to search the global optimal design in the feasible design space. The optimization problem is defined as

$$maxf\left(x_1,x_2\right)$$

$$f\left(x_1,x_2\right)=3{\left(1-x_1\right)}^2{\mathrm{e}}^{-{x_1}^2-{\left(x_2+1\right)}^2}-10\left(\frac{x_1}{5}-{x_1}^3-{x_2}^5\right){\mathrm{e}}^{-\left({x_1}^2+{x_2}^2\right)}–$$

$$\frac{1}{3}{\mathrm{e}}^{-{\left(1+x_1\right)}^2-x_2^2}$$

$$\mathit{\text{Submit to}}\left\{-2,-2\right\}\le \left\{x_1,x_2\right\}\le \left\{2,2\right\}$$

$$g_1\left(\mathrm{x}\right)=\frac{15}{80}{x}_1^2-\frac{5}{8}{x}_1+x_2-2\le 0$$

$$g_2\left(\mathrm{x}\right)=\frac{1}{2}{x}_1-x_2-1\le 0$$

Three local optimal designs, marked ①, ②, ③, exist within the feasible design space depicted in Figure 11. It is a contour plot, and the feasible design space is enclosed within the shadowed area. The following two cases are performed to investigate the performances of convergent accuracy and the number of functional evaluations.

1. Case of large deviation (σ=5Δd): a sequential approximation is performed using the proposed two steps as shown in Figure 7. Value of β is set to 0.6. The method of training the HNN is: to start with expansion terms of 2, and training accuracy of 0.02,which if not satisfied after 1000 epochs, is gradually increased to a higher-order expansion by two orders. The thresholds of Eqs. (19) and (20) are ${\varepsilon }_1$ =0.15, ${\varepsilon }_2$ =0.01. Step 1 is terminated after nine iterations (corresponding to 22 functional evaluations) of random design searches. Then, Eq. (20) is satisfied after two iterations and the total number of functional evaluations until convergence is 24. The approximation optimal design is $x_{\mathit{OP}}\left(\mathrm{k}\right)$ = (−0.013, 1.583), $\overline{f}\left(x_{\mathit{OP}}\left(\mathrm{k}\right)\right)$ =8.106. The history of each iteration is illustrated in Figure 12, in which the horizontal axis stands for the distance (ds) between the original position (0,0) and the approximation optimal design.

2. Case of small deviation (σ=2Δd): The training procedure and stop thresholds are similar to the case of large deviation. Step 1 is terminated after six iterations (corresponding to 16 functional evaluations) of random design searches. Then, Eq. (20) is satisfied after two iterations and the total number of functional evaluations until convergence is 18. The approximation optimal design is =8.105. The history of each iteration is illustrated in Figure 13, in which the horizontal and vertical axes have the same meaning as those of Figure 12. Further inspection of the global approximation between the above two cases is performed using 1600 (40×40) lattice design points covering the entire design space. The deviations of cases (1) and (2) are 1.12 and 1.32, respectively. That is, a larger deviation of Gauss random designs gives slower convergence but higher accuracy of approximation over the entire design space. On the other hand, a smaller deviation of Gauss random designs gives rapid convergence but lower accuracy of approximation over the entire design space. Let us confirm the effect of Eqs. (22) and (23). If the threshold termination in Step 1 is set to ε₁=0.15, the number of designs is 12 and 22 in case α=2 and α=5, respectively. If the threshold termination in Step 1 is set sufficiently small, say, ε₁=1.0×10⁻³, then the history of convergence is depicted in Figure 14 and the number of functional evaluations in Step 1 becomes 36. However, on the basis of the result shown in Figure 14 we know that there is no further improvement in determining the global optimal design after 36 functional evaluations, except to increase the probability. On the other hand, d₀=0.9 can be determined for this problem, and Step 1 is terminated when ε₁=0.9 at eight designs. In practice, however, this is not known; therefore, if one attempts to search for the global optimal design including small clusters, one is expected to select a stricter termination threshold in Step 1 and a larger deviation of Gauss random designs.

Figure 11.

Contour plot of objective function and constraints: It is a contour plot, and the feasible design space is enclosed within the shadowed area. Three local optimal designs, marked ①, ②, ③ exist within the feasible design space.

Figure 12.

Convergence of large standard deviation: The horizontal axis stands for the distance (ds) between the original position (0,0) and the approximation optimal design.

Figure 13.

Convergence of small standard deviation: The horizontal axis stands for the distance (ds) between the original position (0,0) and the approximation optimal design.

Figure 14.

Change of distance in large standard deviation: The number of functional evaluations in step 1 becomes 36.

5. Application to vehicle problems

5.1 Vehicle idling vibration problem

5.1.1 Modeling for optimal problem

NVH such as noise, vibration and harshness is very important for merchantability of vehicle. Here idling vibration is treated as one of NVH. Many studies have been done related to idling vibration where two types of modeling exist. One modeling is that body is separated from powertrain and engine mount [40]. The other one consists of powertrain, engine mount and body together where body is modeled very simply such as solid or as only 1st bending and twist modes in most of studies [41, 42]. Here powertrain and engine mount are modeled together with rather precise body model.

Here powertrain and engine mount are modeled as shown inFigure 15where weight of powertrain is supported by 3 main mounts. Each mount is modeled with 3 springs in x, y, z directions and damping. The weight of each mount is so small compared with body and powertrain that its weight is ignored.

Figure 15.

Mounts model: Weight of powertrain is supported by 3 main mounts. Each mount is modeled with 3 springs in x, y, z directions and damping.

Powertrain is so stiff compared with body that it is modeled as 6 rigid body modes. Body is presented with 44 modes under 50 Hz which are calculated by eigenvalue analysis for the FEM model in Figure 16. Total model is generated by component mode synthesis [43] by using spring parameters, damping for engine mount and by vectors, eigenvalues and eigen damping for body and powertrain. This total model has 50 numbers of freedom which is much lower than that of physical coordinates. Only engine mount is changed in repeated optimal calculations.

Figure 16.

Vehicle FEM model: Body is presented with 44 modes under 50 Hz.

5.1.2 Setting up optimization problems

5.1.2.1 Position arrangement of engine mount

It is shown 34 candidates for position arrangement of engine mount in Figure 17. They are arranged so that powertrain is enclosed. 3 points are selected among these 34 points. Among the 3 selected points $p_{t1}$, $p_{t2}$, $p_{t3}$, the next relationship is held:

$$1\le p_{t1}<p_{t2}<p_{t3}\le 34$$E26

Figure 17.

Candidates for location of mounts: 34 candidates for position arrangement of engine mount are arranged so that powertrain is enclosed. 3 points are selected among these 34 points.

When these 3 mounts are placed closely, support for powertrain becomes unstable.

To deal with this, the following restrictions are added:

$$p_{t2}-p_{t1}\ge 5$$

$$p_{t3}-p_{t2}\ge 5$$E27

$$p_{t1}+34-p_{t3}\ge 5$$

Spring and damping parameters are continuous variables. Following constraints are given to guarantee other characteristics than idling vibration.

• Front and rear stiffness for dynamic vibration and acceleration shock when the engine starts:

$$\sum _{i=1}^3\mathrm{kxi}=35.0\left[\mathrm{kgf}/\mathrm{mm}\right]$$E28

• Left and right stiffness for handling and stability

$$\sum _{i=1}^3\mathrm{kyi}=45.0\left[\mathrm{kgf}/\mathrm{mm}\right]$$E29

• Up and down stiffness for engine shake

$$\sum _{i=1}^3\mathrm{kzi}=45.0\left[\mathrm{kgf}/\mathrm{mm}\right]$$E30

Here k_di (d=x, y, z, i=1,2,3) are stiffness parameters of number i, direction d.

All parameters have upper and lower limits as follows:

$$7.0\le {\mathrm{k}}_{\mathrm{xi}}\le 21.0$$

$$9.0\le {\mathrm{k}}_{\mathrm{yi}}\le 27.0$$E31

$$9.0\le {\mathrm{k}}_{\mathrm{zi}}\le 27.0$$

$$0.02\le {\mathrm{c}}_{\mathrm{di}}\le 0.1\left(\mathrm{d}=\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{i}=\mathrm{1,2,3}\right)$$

Damping parameter is also constrained as follows:

$$\sum _{i=1}^3{c}_{\mathit{di}}=0.1\left(\mathrm{d}=\mathrm{x},\mathrm{y},\mathrm{z}\right)$$E32

Total number of design values is 21 such as mount location (3×3), stiffness parameters of 3 directions (3×3) and damping parameter (3).

From Eqs. (27) and (30), follows are induced:

$${\mathrm{K}}_{\text{d3}}=35.0-\left({\mathrm{k}}_{\text{d1}}+{\mathrm{k}}_{\text{d2}}\right)$$

$$14.0\le {\mathrm{k}}_{\text{d1}}+{\mathrm{k}}_{\text{d2}}\le 28.0\left(\mathrm{d}=\mathrm{x},\mathrm{y},\mathrm{z}\right)$$E33

This means that 21 design values and 6 linear equation constraints induce 15 design values (as shown in Table 2) and 6 boundary values (as shown in Table 3).

	Mount1	Mount2	Mount3
Location	pt1	pt2	pt3
Stiffness (x)	kx1	kx2	kx3
Stiffness (y)	ky1	ky2	ky3
Stiffness (z)	kz1	kz2	kz3
Damper (x)	Cx1	Cx2	Cx3
Damper (y)	Cy1	Cy2	Cy3
Damper (z)	Cz1	Cz2	Cz3

Table 2.

21 design variable candidates. Since there are 6 boundary conditions as shown in Table 3, 15 of them are design variables.

Lower		Upper
14.0	kx1 + kx2	28.0
18.0	Ky1 + ky2	36.0
18.0	Kz1+ kz2	36.0
0.04	Cd1 + Cd2(d=x,y,z)	0.08

Table 3.

Boundary conditions.

5.1.2.2 Settings of objective function

Idling vibration is estimated at 2 observer points where the occupant seat is installed on the floor. These points are intersection points between left and right side sills and cross member where the seat is installed. The frequency vertical responses of these 2 points are shown in Figure 18. Idling vibration is the phenomena during engine rpm 600–900 and so it is estimated by frequency response from 20 Hz to 30 Hz. As a result, the objective function is the sum of 2 integral functions with frequency in this range. Trapezoidal rule is used for integration with integral interval 0.2 Hz.

Figure 18.

Acceleration response for frequency at observation point: Idling vibration is the phenomena during engine rpm 600–900 and so it is estimated by frequency response from 20 Hz to 30 Hz.

5.1.3 Optimization results

First, with N initial samplings, response values are calculated. N means at least necessary number for the approximation with secondary polynomials;

N= (n+1)(n+2)/2 (n is the number of design variables, here, n=15. So N=136). Each sampling point is selected among feasible region which is given from Eqs. (25), (26), (30) and Tables 2 and 3. Both of ${\varepsilon }_1$, ${\varepsilon }_2$ are set up 0.05 and ${\alpha }_i$ is set up 1.5. Both of GA and SQP (Sequential Quadratic Programming) are used to explore the optimal value on the approximated response surfaces.

Three types of calculation are performed as far as initial point for optimal analysis by SQP. In 1st one, it is selected the point where it is the optimal point using GA. In 2nd one, last optimal point is selected as the initial point. In the last one, the best solution point among sample points is selected as the initial point. In each step, these 3 calculations are performed. However, because the solution by SQP is continuous value, at last stage discrete value is searched around the solution point. It is the optimal point which has the least value among three types calculation in each step.

5.1.3.1 Comparison with GA

To certify the effectiveness of MPOD, it is compared with GA. Here general GA is selected although there are many types of GA. Real-coded Genetic Algorithms is adopted with elite preservation strategy where population size is 50 and number of generations is 100. It is shown in Table 4 the comparison between GA and MPOD.

		MPOD	GA
A	Response of optimum design	0.244	0.240
B	Number of sampling points	189	5000
C	Total calculation time [sec]	1571.8	1256.8
D	Calculation time of response for all sampling points [sec]	43.7	1150.0
E	Rate of D for C [%]	2.78	91.50

Table 4.

Comparison between result of MPOD and one of GA.

As far as the problem where it takes much time to calculate objective function, it is very effective to use MPOD.

In term A, it means the precision of optimal value is almost same. In term B, it is compared the number of calculation of objective function. It is fewer in MPOD than in GA. But in term C, total time is shorter in GA than in MPOD. The calculation time of objective function for one selected point is 0.23 sec. In term D, all calculation times for object functions are filled out. In term E, in MPOD, the calculation time of objective function accounts for about 13% in the total calculation time although in GA it accounts about 90%. In other words, the total calculation time in GA is more influenced by the calculation time of object function than in MPOD. On the other hand, the most time-spent process in MPOD is to search the optimal solution on the approximation function updated every time. In GA, this process is finished within few seconds although in SQP, it takes sometimes few minutes with some initial values. And so in MPOD, it influences very much to the calculation time how the optimal value can be searched on the approximation function. As far as the problem where it takes much time to calculate objective function, it is very effective to use MPOD.

5.2 Booming noise problem for simplified vehicle model

5.2.1 Modeling for optimal problem

Here booming noise for wagon-type automobile is treated. This noise is from 40 to 60 hertz when a vehicle gets over a small protrusion.

This is the coupled acoustic-structural vibration problem [44, 45, 46, 47, 48].

It is shown the simplified model simulated vehicle system in Figure 19 which has average dimensions such as width 1.5 m, length 3 m, height 1.2 m. It is supposed that the nodal shape box is body and inside the box is interior space. The FEM model is generated such that the body has 1980 rectangular shell elements with length of one side 0.1 m, thickness 1 mm and inside the box has 5400 cube acoustic elements with length of one side 0.1 m. To treat the coupled forward and rear direction acoustic-upper panel structural phenomena, the 5 faces except upper panel are high stiffness by 100,000 times Young modulus of these 5 faces such that the 1st eigen-frequency of acoustic system and second panel bending one are near. The mode attenuation ratio is 3.5% both for the structural system and acoustic system.

Figure 19.

Simplified vehicle FE model: To treat the coupled forward and rear direction acoustic-upper panel structural phenomena, the 5 faces except upper panel are high stiffness by 100,000 times young modulus of these 5 faces such that the 1st eigen-frequency of acoustic system and second panel bending one are near.

5.2.2 Setting up optimization problems

The problem is set as follows:

The objective function: the average of 8 frequency integral values of each sound pressure response absolute value at observation points shown in Figure 19. The eight coordinate positions are 1 m from the front or the back panel in Z direction, 0.3 m from the right or the left panel in Y direction and 0.3 m from the upper or the bottom panel in X direction symmetrically. Integral values of each sound pressure response absolute value are calculated by trapezoid official over from 40 to 50 hertz divided into 200 pieces:

Design values: $X_1$,$X_2$.

X₁: the position coordinate of the roof bow on the roof panel

0.5 $\le X_1\le 2.5$, 0.1 m increments, 21 point discrete values

$X_2$: panel thickness 0.6 $\le X_2\le 1.8$, 0.3 mm increments, 5 point discrete values. Input force: Sinusoid forced displacement in forward and rear direction with amplitude 0.2 mm which is equivalent to the force which is generated when the vehicle gets over a small protrusion.

5.2.3 Optimization results

Because the number of design values is 2, 6 initial learning points are selected through uniform random number by keeping longer than $d_{\mathit{min}}$. Here $d_{\mathit{min}}$ is 0.33, ${\varepsilon }_{1.}$ in Eq. (19) is 0.05, ${\varepsilon }_2$ in Eq. (20) is 0.01 and ${\alpha }_1$ = ${\alpha }_2$ in Eq. (21) are 2.0. The distribution of learning points is shown in Figure 20. In step 1, 6 points marked 1st and 3 points marked 7th–9th are set over entire design space. On the other hand, 3 points marked 10th–12th in step 2 are concentrated on one part of the space. In step 1, although after 3 repetition, Eq. (19) is satisfied with $X_1$ =2.5 and $X_2$ =0.6, Eq. (20) is not satisfied. And so in step 2 it is repeated. The first optimal point in step 2 is $X_1$ =1.3 and $X_2$ =0.6 apart from $X_1$ =2.5 and $X_2$ =0.6. And finally, Eq. (20) is satisfied with $X_1$ =1.5 and $X_2$ =0.9. This is because there are 2 design points which have almost same minimum values, one design point is $X_1$ =1.5, the other is $X_1$ =2.5. And so the optimal point moves between 2 minimum points on the way to be improved accuracy on response surface. As a result, it converges to the optimal solution $X_1$ =1.5, $X_2$ =0.9 and Y=95.59 with 12 learning points. In such way, optimal value is searched by repeating step 1where global solution is searched and by repeating step 2 where approximation precision is improved on the response surface. The approximation response surface is shown in Figure 21. Although referring to Figure 19, originally, the response surface must be symmetric with respect to $X_1$ =1.5, it is asymmetric because the learning points have biased distribution. This shows we can abbreviate unnecessary learning points for getting minimum value by taking many learning points in neighborhood of optimal point.

Figure 20.

Learning points distribution in MPOD: In step 1, 6 points marked 1st and 3 points marked 7th–9th are set over entire design space. On the other hand, 3 points marked 10th–12th in step 2 are concentrated on one part of the space.

Figure 21.

Approximate response surface by MPOD: The response surface must be symmetric with respect to $X_1$ =1.5, it is asymmetric because the learning points have biased distribution. This shows we can abbreviate unnecessary learning points for getting minimum value by taking many learning points in neighborhood of optimal point.

5.2.4 Comparison between MPOD and experimental planning method using orthogonal table (EPO)

As shown in Figures 22, 25 points at 5 levels are selected by quartering upper and lower design limits. Figure 23 shows the response surface approximated by 4th order direct polynomial. It is high in the symmetric level of response surface with respect to $X_1$ =1.5 because the design points are uniformly selected. It has 3 valleys in the neighborhood of $X_1$ =0.5, 1.5, 2.5 and 2 mountains in the neighborhood of $X_1$ =1.0, 2.0. Gradient of response surface for $X_2$ is extremely smaller than the one for $X_1.$ The optimal solution by this response surface is $X_1$ =1.5, $X_2$ =0.6 and Y=94.36 which are slightly different from $X_1$ =1.5, $X_2$ =0.9 and Y=95.59 by the MPOD. The reason for this difference is that although in EPO the trial points must be placed uniformly, in MPOD they are placed adaptably such as the distance between any 2 points is longer than d_min The minimum point by EPO exits in the region within d_min =0.33 where it is out of this problem.

Figure 22.

Designing point distribution in EPO: 25 points at 5 levels are selected by quartering upper and lower design limits.

Figure 23.

Approximate response surface by EPO: This response surface is approximated by 4^th order direct polynomial. It is high in the symmetric level of response surface with respect to $X_1$ =1.5 because the design points are uniformly selected.

5.3 Booming noise problem for real vehicle model

5.3.1 Modeling for optimal problem

Following the previous section, it is applied the MPOD also to the booming coupled noise and structural problem for a real wagon-type vehicle where the dynamic damper should be set. It is shown the overview of vehicle FEM model in Figure 24 where the structural system is modeled with shell elements, the acoustic system is modeled with solid elements and the dynamic damper is modeled with lumped mass and one degree of freedom scalar spring.

Figure 24.

(a) Structural FE model for trimmed body. (b) Acoustic FE model of cabin. Overview of vehicle FEM model: The structural system is modeled with shell elements, the acoustic system is modeled with solid elements and the dynamic damper is modeled with lumped mass and one degree of freedom scalar spring.

5.3.2 Setting up optimization problems

Let us set up the objective function Y and the design variables X1, X2 as follows:

Y: Minimize the maximum SPL (Sound Pressure Level) during 30 and 60 Hertz at the observation point ● which is set ear height position of the center of the front seat under the condition such that forward and rear unit input is loaded at the connected point of the arm to the front suspension for giving the same-phase vibration on the left and right. This load is generated when the vehicle gets over a small protrusion.

X1: the position of the damper which is selected among 22 candidate from ➀ to ㉒ in Figure 25.

Figure 25.

Candidates of dynamic damper location: The position of the damper which is selected among 22 candidate from ➀ to ㉒.

X2: Spring constant which is selected among 5 candidates: 7.0,9.0,11.0,13.0,15.0 kgf/mm. The weight and the structural damping coefficient are constant such as is 1.5 kg and 0.1 each other.

5.3.3 Optimization results

As far as the convergence judgment for the optimal, ε₁ is 0.05, ε₂ is 0.01 and also as far as the control coefficient, α₁ and α₂ are 2.0. From the approximation response surface with the initial 6 points, the 1st optimal is X1=16, X2=11.0, Y=0.659. The 2nd optimal is X1=15, X2=7.0, Y=0.567 by adding 7th, 8th 2 points. The step 1 is terminated because the 3rd optimal is X1=15, X2=7.0 that is same as last time by adding 9th, 10th 2 points. It is finished because the optimal solution X1=14, X2=7.0, Y=0.556 are gained by adding the 11th point and this satisfies with Eq. (20). It is shown the learning distribution points in Figure 26 and the approximated response surface using 11 learning points in Figure 27. This means the minimal solution over the total design region is gained with very few learning points compared with total 110 design points. As shown in Figure 28. The objective function Y is reduced 6 dB around 45Herz by the optimal analysis. As shown in Figure 29, the coupled acoustic structural mode at 44.3 Hz gives the entire car body elastic deformation induced by the deformation of the upper panel and rear floor. Moreover the acoustic system gives 1st mode of which node is a little rear from the center of passenger cabin and of which anti-node is front to back end of passenger cabin. It is shown the difference of the panel vibration acceleration along the center line of vehicle body in Figure 30 with the optimally optimum arranged damper and without damper. From Figure 30, the effect of optimal arrangement of damper can be interpreted to be the vibration-absorbing effect for the global coupled acoustic-structural mode at 44.3 Hz rather than the local vibration reduction around the front roof rail.

Figure 26.

Learning points distribution: The minimal solution over the total design region is gained with very few learning points compared with total 110 design points.

Figure 27.

Response surface by 11 learning points.

Figure 28.

Estimation of optimized dynamic damper effect: The objective function Y is reduced 6 dB around 45Herz by the optimal analysis.

Figure 29.

Fluid–structure coupled eigen mode (44.3 Hz): The coupled acoustic structural mode at 44.3 Hz gives the entirecar body elastic deformation induced by the deformation of the upper panel and rear floor.

Figure 30.

Acceleration of vehicle body panel (44.3 Hz): Optimal arrangement of damper can give the vibration-absorbing effect for the global coupled acoustic-structural mode at 44.3 Hz rather than the local vibration reduction around the front roof rail.

To certify the effect of the optimal damper arrangement, bench test using smooth drum is performed with the optimal arrangement damper and without damper. It is shown the noise level–frequency relationship in Figure 31. Maximum approximately 5 dB reduction of SPL is realized among 45∼50 Hz. Although it needs many man-hours to find out the optimal arrangement of dynamic damper by experiment trial and error, it has been certified the method for booming noise reduction can be realized efficiently by MPOD.

Figure 31.

Validation of optimized dynamic damper effect: Maximum approximately 5 dB reduction of SPL is realized among 45∼50 Hz.

5.4 Vehicle crash problem

5.4.1 Setting up optimization problems

The vehicle side member (component parallel to the central axis of vehicle) plays a big role in energy absorption while the frontal collision and energy absorption of the vehicle is determined by its size, shape and welding [49, 50, 51, 52, 53]. In this example, the uniform hollow square-cross section, a perfectly straight side member with uniform thickness is investigated and reinforcement of the component is considered as a way to increase the energy absorption. The size and material property of the member are ρ = 7.85E-3 g/m${\mathrm{m}}^3$, E = 210GPa, $E_t$ = 2.5GPa, ${\sigma }_y$ = 220 MPa, v = 0.3, and the form of reinforcement is depicted in Figure 32.

Figure 32.

Analysis model of crashworthiness optimization: The member is perfectly straight with uniform hollow square-cross section and uniform thickness.

All degrees of freedom at the bottom end are rigidly fixed, and at the top end, a rigid mass of 500 kg, and velocity of 54 km/h are used as a load to simulate a crush. The load–displacement behavior of the member while crushing is calculated by FEM solver LS-DYNA3D [54], in which 2700 shell elements and 2754 nodes are used. The design variables include the thickness of the base plate (t₁), the upper plate (t₂) and the location of reinforcement (z). The total weight of the component is constant (w=780 g), therefore only two design variables are independent. The objective function is the energy absorption of the component, 10 ms after crushing. The mathematical definition of the problem is as follows;

$$Maxf\left(t_1,t_2,z\right)$$

$$\text{Subject to}:0.5\le t_1,t_2\le 1.0,1.5\le z\le 250$$

$$t_2=\left(1.5+\left(1-t_1\right)\right)/2.0$$

5.4.2 Optimization results

The training procedure and end thresholds are identical to the case of large deviation, that is, α=5, β=0.6 and ε₁=0.15, ε₂=0.01. Step 1 is terminated after four iterations (corresponding to 12 functional evaluations) of random seed searches. Then, Eq. (20) is satisfied after two iterations and the total number of functional evaluations until convergence, is 14. The approximation optimal is (1.21, 1.39,160), f (x_op(k))=8484.11 kJ, which is 11.6% higher than that of the original design (1.00,1.50,0.0). The contour plot of the approximation function is depicted in Figure 33. In the figure, the open circle ○ indicates the optimal design using the proposed MPOD method and the plus sign+ expresses the designs used for approximation. The comparison between the approximation function and FEM values with t₁=1.21 and t₂=1.39 are fixed while the changing reinforcement location z is plotted in Figure 34. It can be seen from the figure, and also mentioned in Refs. [49, 50, 51, 52, 53], that the FEM value varies, because of the heavy nonlinearity of the crashworthiness problems [49, 50, 51, 52, 53]. A smooth and robust approximation function is obtained by the HNN approximation and an optimal design is realized successfully. The MPOD approach based on the HNN has a robust property against calculation noise. However, the deviation of the coefficient is dependent on the design distribution. This topic will be discussed further in future works.

Figure 33.

Approximation function of the present method: The open circle ◯ indicates the optimal design using the proposed MPOD method and the plus sign+ expresses the designs used for approximation.

Figure 34.

Comparison of FEM and HNN when design variable z is changed: The FEM value varies, because of the heavy nonlinearity of the crashworthiness problems. A smooth and robust approximation function is obtained by the HNN approximation and an optimal design is realized successfully.

5.5 Combination problem between vehicle crashworthiness and restraint device performance at the head-on collision

5.5.1 Setting up optimization problems

Here the problem is the head-on collision where the vehicle collides perpendicular to a concrete wall. This phenomena is that the vehicle finishes the moving after the crush energy is absorbed completely by the vehicle structure [55]. During this phenomena, the occupant moves length L+I in the axial direction where L is the amount of crush and I is the amount of movement of the occupant in the cabin. The longer L+I is the smaller the average deceleration of the occupant is. And so it is necessary to raise L+I. On the other hand, if the amount of crush of the vehicle front-end is too large, the possible amount of occupant moving becomes smaller. Also if the movement amount of the occupant to the vehicle front side, it arises the possibilities of the crush between the occupant and the steering wheel, etc. In such way, the most important issue is to get the optimal valance between the vehicle crashworthiness and occupant restraints performances. So let us start the feasibility study concerned with this. Ever before, representative 4 patterns have been considered as the best vehicle performance for the occupant. As shown in Figure 35, these 4 patterns of vehicle crush performance have been referred to the standard numbers which are generated when the small passenger cars collide perpendicular to a concrete wall with 35mph based on the criterion of NHTSA (National Highway Traffic Safety Administration). Here it is selected among these 4 vehicle crush performances as the design of vehicle. That means the design vehicle value becomes integer for the vehicle performance. As far as the occupant restraint, it is supposed that there is not the contact between the occupant and the steering wheel for simplicity so that only seat belt is considered. For the design values of seat belt, the magnification $x_1$ of standard value 2.5GPa of seat belt Young modulus and the magnification $x_2$ of standard value 3000 N of working load of the load limiter are selected. Here the load limiter is the device which relieves the impact on the chest by sending seat belt with the load remained constant when the load applied to the seat belt winder reaches to the working load as shown in Figure 36. Here the send-off amount of load limiter is 100 mm. Next it is shown in Figure 37 the analysis model of dummy system which consists of dummy, seat belt and sled. The dummy is Hybrid III whose size is the average of adult male. The number of total nodal points is 3700. The number of shell elements is 3203 for dummy and 116 for sled. Seat belt is basically modeled with 58 beams. But the parts of the seat belt for the chest and the waist are modeled with 192 shell elements so that the contact between dummy and seat belt can be considered precisely. The dummy and the sled are solid and seat belt is elastic. Contact by penalty method is defined for the friction between dummy and seat. The coefficient of friction is 0.2. The model for seat is comparatively coarse because it is only for the in-plane force to get friction as shown in this figure. The objective function is Eq. (35) by using HIC in Eq. (34) and chest G which are generated by head-on collision with crush velocity 35mph referred to the criterion of NHTSA.

$$\mathrm{HIC}={\left[\frac{1}{t_{2-t_1}}{\int }_{t_1}^{t_2}\alpha \mathit{dt}\right]}^{2.5}\left(t_2-t_1\right)$$E34

$$Min\mathrm{f}\left(\mathrm{n},x_1,x_2\right)=\sqrt{{\left(\frac{\mathit{HIC}-200}{200}\right)}^2+{\left(\frac{\mathit{\text{ChestG}}-20}{10}\right)}^2}$$E35

$$\text{Subject to}1\le \mathrm{n}\le 4$$

$$0.7\le x_1\le 1.3$$

$$0.5\le x_2\le 1.5$$E36

$$H_{\mathit{\text{disp}}}\left(\mathrm{n},x_1,x_2\right)\le 420\mathrm{mm}$$

Figure 35.

Type of deceleration curve of vehicle: These 4 patterns of vehicle crush performance have been referred to the standard numbers which are generated when the small passenger cars collide perpendicular to a concrete wall with 35mph based on the criterion of NHTSA.

Figure 36.

Safety belt pull out vs. load of load limiter: The load limiter is the device which relieves the impact on the chest by sending seat belt with the load remained constant when the load applied to the seat belt winder reaches to the working load.

Figure 37.

Hybrid III dummy and sled FEM model: The analysis model of dummy system which consists of dummy, seat belt and sled.

Here, $\alpha$ in Eq. (34) is deceleration of head. The unit of deceleration for HIC and chest G is used g(1 g=9.8 m/ $s^2$). Chest G is the maximum deceleration which continuously generates longer than 3 ms.

The objective function in Eq. (35) is rated on a 5-point stage. Here $H_{\mathit{\text{disp}}}$ in Eq. (36) means the moving amount of head to the front and it is shorter than 420 mm so that the occupant does not crush to the steering.

5.5.2 Optimization results

6 points are selected for the initial learning points. ${\varepsilon }_1=0.02$ in Eq. (19), ${\varepsilon }_2=0.05$ in Eq. (20), ${\alpha }_1$ = ${\alpha }_2={\alpha }_3$ =1.0 in Eq. (21). Firstly step 1 is terminated with 5 iterations which use total 22 points. In step 2, it is terminated with 1 iteration which uses 2 points. So 24 points are used. The optimal design point is n=3, $x_1=0.70,x_2=1.20.$ The value of the objective function is 1.34 (HIC=355.8、Chest G=30.9) of which estimation is 5☆ in Figure 38. The other local optimal design point is n=3, $x_1=1.30,x_2=1.15.$ of which objective function is 3.12 (HIC=620.2, Chest G=43.1). This corresponds to 4☆ in Figure 38. It is shown the contour of the objective function with n=3 in Figure 39. In this figure, A is the given optimal point, B is the local optimal point and the shadow part means the one where it cannot be designed. Although all trial points are given within the region, A is on the boundary. This means the optimal point A is never found without the MPOD’s intrinsic extrapolation ability. Such result is never found only with the conventional interpolation unless the trial points are included also on the boundary. In such style, it becomes very efficiency. It is shown the motion behavior of the dummy model at optimal design in Figure 40. The HIC becomes the maximum when the jaw contacts the chest around 100 ms as shown in this figure. But the maximum HIC is rather small because the contact between the jaw and the chest becomes loose around 100 ms by the operation of the load limiter. Generally speaking, lots of studies for response surface optimal problem have been done for only one peak problem. On the other hand, MPOD gives us plural optimal combination simultaneously without the knowledge such as how many and where the maximum points exist. Now present vehicle also has air back in addition to seat belt as restraint device. And so it has more local maximum points in the optimal combination problem among vehicle and restraint device than this problem mentioned here. That means the best combination has not been realized between vehicle performance and restraint device one. Only the MPOD can be applied to this multi-peak optimal problem.

Figure 38.

Safety criteria using NHTSA: The optimal design point is n=3, $x_1=0.70,x_2=1.20.$ The value of the objective function is 1.34 (HIC=355.8、chest G=30.9) estimated as 5☆. The other local optimal design point is n=3, $x_1=1.30,x_2=1.15.$ The objective function is 3.12 (HIC=620.2, chest G=43.1) estimated as 4☆.

Figure 39.

Contour of approximation response surface at n=3: A is the given optimal point. Although all trial points are given within the region, a is on the boundary. This means the optimal point a is never found without the MPOD’s intrinsic extrapolation ability.

Figure 40.

Motion behavior of the dummy model at optimal design: The HIC becomes the maximum when the jaw contacts the chest around 100 ms. but the maximum HIC is rather small because the contact between the jaw and the chest becomes loose around 100 ms by the operation of the load limiter.

6. Conclusion

Here, it is presented a novel method to solve the global optimal design within a feasible design space using response surface methodology based on HNN activation. The sequential approximation procedure and extrapolation of the HNN are proposed for irregular design space and arbitrary designs. The detailed summaries of the issues raised in this paper are as follows.

1. RSM has become more and more important especially for non-linear problems. There are 2 important issues in RSM such as the selection of activation function where quadratic polynomials are most used and the selection of the design trial points where CCD is most utilized. Also artificial neural network function is applied where the multi-layered one is generally applied. In contrast we have developed MPOD where both of the selections of design trial points and the activation function are different from others.

2. In MPOD the number of design trial points can be reduced much more than others. And it is very efficient when the calculation time for one design trial point is long in real big and non-linear problems. Generally the more layers there are the higher the precision becomes in multi-layer neural network. In contrast it is shown that higher orders of the HNN can be realized by the products of fundamental terms to improve the approximation capability.

3. In higher order case with HNN, it becomes often necessary to use pseudo inverse matrix where we have developed a superior method using penalty function to Moore–Penrose generalized inverse matrix.

4. It is estimated that the HNN is more feasible to approximate the multiple peaks function than the polynomials.

5. The extrapolation can be performed by adjusting the transformation coefficient which expresses the slope of the transformation function of the HNN.

6. It is illustrated that the proposed sequential approximation procedure involves two steps. Step 1 aims to conduct a rapid scan of the entire design space with fewer designs, Step 2 focuses on increasing and confirming the approximation accuracy near the optimal design determined in Step 1. The proposed sequential optimal design method can be applied to design problems with regular and irregular boundaries of design space by the benefit of extrapolation.

7. Rapid approximation of the entire design space can be performed by presetting the threshold d_min to be larger, however, it can also feasibly be used to search for relatively small clusters by adjusting the deviation of random seeds and presetting the threshold d_min to be smaller in Step 1. However, the deviation of random seeds is assigned smaller values in Step 1 to increase the probability of solving the design within the small cluster. In other words, a trade-off between the accuracy over the entire design space and computation expense can be comprehensively achieved.

8. Comparing with GA, the proposed method can suppress the number of samples and has the same level of accuracy as GA. That means MPOD is much superior to GA in the case of long calculation time problem for one sample such as non-linear problem.

9. Applications to the vehicle problems such as idling vibration, booming noise, crashworthiness design and crash occupant restraint problem with multiple peaks reveals that the proposed method is a feasible and practical approach.

10. Some points still need to be discussed, such as the robustness of the optimal design and the confident regions of estimated coefficients of holographic neural network, which will be addressed in the future.

Acknowledgments

This research has been done by many co-authors. Especially this author thanks to Dr. Shi, Q., Dr. Kamada, Y., Dr. Tokura, S., Mr. Takasima, F. and Mr. Fukushima, H.