Surrogate-based multi-objective optimization of firing accuracy and firing stability for a towed artillery

2.1. Multi-body rigid-flexible dynamic modeling

Fig. 1 shows a diagram of the towed artillery. When it is firing the gunpowder combustion produces a large number of high pressure gas, which promotes the projectile a forward movement along the gun bore with a great acceleration. At the same time the recoiling parts, which including the breechblock, recoil mechanism, barrel and muzzle brake, make a backward recoiling and then forward counter recoiling motion under the joint action of the forces of the gun bore and recoil resistance.

Fig. 1Diagram of firepower part of a towed artillery: 1 – spade, 2 – undercarriage, 3 – breechblock, 4 – trunnion, 5 – recoil mechanism, 6 – cradle, 7 – barrel, 8 – muzzle brake, 9 – front bushing, 10 –balancer, 11 – top carriage, 12 – running wheel, 13 – base

Basic hypothesis: regardless of the coupling effect between projectile and barrel in the process of firing. And the artillery was in static equilibrium position before firing.

The flexible parts, such as the barrel, cradle and top carriage, were created based on finite element modal neutral file. Their dynamic response was calculated by using the superposition method. The rest parts, which did not affect artillery dynamic response much, were established as rigid parts. Then proper joint relationship and boundary condition was set up between all the parts and the ground. The launching load and forces of recoil mechanism were written in FORTRAN program as a dynamic link library for ADAMS. The topological structure of this dynamic model is presented as Fig. 2.

Fig. 2Topological structure of the dynamic model

2.2. Model verification

In order to verify the artillery launching dynamics model, a series of dynamic verification tests were performed. Usually, there are two aspects need to verify. The first one, testing the recoil displacement and recoil resistance to inspect the rationality of the modeling process and artillery recoil movement. The second aspect, measuring the values of firing dynamic parameters like the angular displacement of the muzzle to evaluate the tactical level of the artillery.

As showed in Fig. 3, a high-speed photography test system was put in place to measure the recoil movement of the towed artillery when firing. It is mainly composed of a high-speed camera, optical lens, a fixed bracket, a computer and an image tracking processing software, and other parts. The Phantom V710 high-speed camera has a wide shooting frequency range, from 6,000 fps to 1,400,000 fps, its minimum exposure duration is 1 us. Optical lens is adopting the Nikkor AF-S 400 mm f/2.8D ED lens. Image tracking software using ProAnalyst software developed by company Xcitex. This test system can catch the white-on-black test marking (showed in Fig. 3) precisely in the whole firing process. Thus, we can get kinematic characteristics of the launching artillery system.

As showed in Fig. 4, a testing device for muzzle dynamic response was set. It is mainly composed of a SDI-ARG-720 angular rate gyroscope and a DEWETRON 1201 data acquisition system.

Fig. 3High-speed photography test system

Fig. 4The muzzle responses test

10 groups firing test were conducted, but 2 of them were get rid since incorrect collecting data. By comparison between the simulation results and the average of 8 groups testing results in Table 1. The less relative error proves the multi-body rigid-flexible coupling dynamic model can simulate the fire process well. Where, $S_r$ is the recoil stroke length, $F_{R\_\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maxmum of the recoil resistance, and ${\theta }_z$ is the elevation jump angle at muzzle center.

3.1. Optimization objectives and design variables

Considering the artillery tactical and technical requirements especially for firing accuracy and firing stability, some of the firing dynamic parameters, whose dynamic responses can be used to represent the artillery firing tactical and technical properties, were selected as optimization objectives. The elevation and lateral jump angles at muzzle center, the horizontal displacement of the spade and the vertical displacement of base are included. For all the above parameters, we considered their values at the moment when the projectile came out of the muzzle.

In the multi-body rigid-flexible model of the towed artillery, direct structural sensitivity analysis for flexible parties is not easy because it is difficult to parametric modeling for them which are based on the modal neutral file from FEA. So, we used the related multi-rigid body model to do structural sensitivity analysis, then picked five design variables according to the analysis results and special design requirements. They are the vertical and lateral spectral centroid shift$e_y$, $e_z$ of recoiling parts, the front cradle bushing’s axial direction offset $\mathrm{\Delta }{l}_x$, the trunnion center’s vertical offset ${\mathrm{\Delta }h}_y$ and the gap $c_b$ between cradle bushing and barrel. Where, $e_y$, $e_z$ and $h_y$ are three of the most important general parameters for artillery designs, $l_x$ and $c_b$ also play the key roles in firing accuracy and firing stability.

The initial values and value ranges of design variables are presented as Table 2. The initial values come from the original design of this towed artillery, and we determine value ranges of design variables by comprehensive considering expert system and design experience in the past years.

3.2. Establishment of sample library

The design of experiment is the sampling plan in design variable space. The key question is how we assess the goodness of such designs, considering the number of samples is severely limited by the computational expense of each sample.

In order to reduce the size of the training database while keeping the sample representative, optimal Latin Hypercube Sampling (LHS) is used. LHS is a stratified sampling approach with the restriction that each of the input variables has all portions of its distribution represented by input values, and optimal LHS have been proposed to overcome the potential lack of uniformity of LHS. Based on optimal LHS algorithm, according to the initial values and value ranges of design variables in Table 2, a MATLAB program was implemented to calculate 130 training samples and 15 validate samples for design variables. Partial training samples and validate samples are listed below in Table 3.

Based on the above dynamic model in Section 2, according to the values of training samples above, we modified multi-body coupling model respectively. Every flexible body part was modified by altering its finite element mesh model to regenerate a new modal neutral file according to the sample parameters, and importing it into the dynamics model to replace the original component. Then, we can run dynamics simulation to get the relative sample output values. All in this way, the training sample, test sample input and output values are obtained, the sample library was created.

3.3. RBF-BP series combine artificial neural network

In the surrogate model construction stage, we use systematic parameter variation and a selection process that takes into consideration both accuracy and architectural simplicity to find and train the most robust surrogate design for each of the considered target variables.

A genetic optimization RBF-BP series combine neural network was built to create the approximate model. It is a series combine neural network with a RBF sub-network and a BP sub-network double hidden layers. According to the related structures of artificial neural network, the output and input relationship equations can be derived by the artificial neural network theory knowledge [15] as below.

The transfer function of the first RBF hidden layer nodes is Gaussian function, it can be written as:

where $X$ is a $l$ dimensional input vector, $c_i$is the center of the $k$th radial basis function, and it have the same dimension with $X$, ${\sigma }_i$ is the width of the $k$th radial basis function of hidden layer neurons, ${‖X-c_i‖}^2$ represents the Euclidean norm between $X$ and $c_i$, with the increase of it, $Ф\left(X\right)$ will gradually decay till 0. The number of hidden layer neurons in RBF sub-network set as $M$, then the output of the $i$th neurons in the RBF sub-network follows that:

where $b_{1i}$ presents threshold value of the $i$th neurons in RBF sub-network, $w_{1i}$($i=$1, 2, …, $M$) presents the weight from hidden layer to output layer of the RBF sub-network. The input of the $k$th neurons in BP sub-network is the sum of the threshold value and the related weighted $w_{2j}$ of the output in all the neurons in RBF sub-network, that is:

where $b_{2j}$ presents threshold value of the $j$th neurons in BP sub-network. The transfer function of the BP sub-network is Tansig-type function, it can be written as:

So, the output of the $k$th neurons in BP sub-network:

Then the output of the $n$th neurons in RBF-BP series combine network is the sum of the threshold value and the related weighted $w_{3n}$ of the output in all the neurons in BP network, that is:

where $b_{3n}$ presents threshold value of the $n$th neurons of output layer in RBF-BP series combine network. For the given sample input $x_i$($x_i\in R^l$) and expected output $d_i$($d_i\in R^m$), the output error is defined as $\epsilon =d_i-{\tilde{y}}_i$.

In this manuscript, the RBF-BP series combine neural network structure is five inputs, four outputs and two hidden layers, the number of hidden layer neurons $r$ has the approximate relationship with the number l of neurons in the input layer, which is $r=2l+1$, that is 11. But in order to get a good precision, trial method was used to define the number of the two hidden layer neurons from 9 to 20, and found that when the first and second hidden layer neurons was 15, 15 respectively, a total of 154 weights and 34 thresholds can get better approximating precisions.

The data were initially pre-processed to remove outliers. Additional inputs were formed by means of aggregating operators. The input and output variables were rescaled in the [0.1, 0.9] range to facilitate ANN training. After normalization of the input sample $X$, first training through sub-RBF neural network, the training results as input to the training of BP sub-network, and finally get the training results. And the network has the ability of error back to study, when the training results cannot meet the precision demand, reverse changes of the neural network weights and threshold until the training results meet the accuracy requirements, finally finished the training.

3.4. Validation of surrogate model

When surrogate model constructed, precision evaluations must pass to ensure the validity of the model. The evaluation includes two aspects: the reappearance of sample points and prediction ability of non-samples space points. Here we taken the most commonly validation method, the coefficient of determination $R^2$, to test this surrogate model. The expression presented below:

where $y_i$, ${\tilde{y}}_i$and $\stackrel{-}{y}$ are dynamic simulation value (true value) of the $i$th test sample, predicted approximate value by RBF-BP series combine network of the $i$th test sample and the mean value of all the dynamic simulation values of the test samples. $n$ is the number of the test samples.

The more of the $R^2$ close to 1, means the more precision of surrogate model. An $R^2$ of 1 means the dependent variables can be predicted without error from the independent variables.

In the process of practical validation on surrogate model, it is usually randomly generated the additional test points in the design for model precision reason. Thus, we use the additional 15 test samples and its dynamic simulation response mentioned in Section 3.3.2 to evaluate the constructed RBF-BP series combine ANN surrogate model. The inspection results show in Table 4. According to that, it can be regarded the established approximate model have a good generalization ability and high prediction accuracy.

4.1. Multi-objective optimization

In MO optimization problems, several conflicting objective functions have to be minimized concurrently. Specifically, in our towed artillery optimization, the elevation jump angle ${\theta }_z$ and lateral jump angle ${\theta }_y$ at muzzle center, the horizontal displacement ${\tau }_x$ of the spade and the vertical displacement ${\lambda }_y$of the base are selected as the optimization objectives. The first two are responsible for firing accuracy and the rest two are for firing stability. We considered their values in the moment of the projectile came out of the muzzle. They can be illustrated below:

where, $X$ are the design variables, their initial value and value range are mentioned in Table 2 previously. $c_1$, $c_2$, $c_3$and $c_4$ are the experience weight coefficients in the evaluation of firing accuracy and firing stability by certain evaluation criterions, here we selected as 0.6, 0.4, 0.3 and 0.7 respectively. ${\theta }_{z0}$, ${\theta }_{y0}$, ${\tau }_{x0}$ and ${\lambda }_{y0}$ are the initial objective value in the original design. And $a_i$, $b_i$ are the lower limit and upper limit of the design variables $X$, which were also presented in Table 2.

One of the most appropriate strategies to determine a good approximate solution of a NP-hard multi-objective problem is a metaheuristic algorithm based on population which evolves along the solution space to find a set of non-dominated solutions. In this paper, we propose a methodology based on the non-dominated sorting genetic algorithm-II (NSGA-II) first introduced by Deb K. [16]. The main objective of algorithm NSGA-II is to find a set of solutions ordered by fronts under the concept of Pareto dominance.

In our work, the proposed algorithm has been implemented in MATLAB, the parameters used for the proposed algorithm are: the rate of recombination is 95 %, the number of iterations is 300, and the size $NN$ of the initial population is 80. Once we have an efficient ANN and a well-set NSGA-II algorithm, we can perform the optimization. Fig. 5 shows the Pareto fronts found by the proposed methodology, it offers a range of possibilities depending on the objectives for the proposed problem. In this case, the best solution depends on the specific interests of the decision maker. He would pick solutions close to side A, if pay more attention to firing accuracy, or close to side D for firing stability.

The first part (A to B) and the second part (B to C) are almost linear with a persistent fall of $F_{obj2}$ for lower values of $F_{obj1}$, but the slope of the Pareto front is greater for the first part than the second. Indeed, in the first part, the most appropriate axial direction offset of the front cradle bushing’s ($\mathrm{\Delta }{l}_x$) will be the optimal solution. But in the second part, it just has little influence on both optimization objectives. And lateral spectral centroid shift of recoiling parts${(e}_z$) and the trunnion center’s vertical offset ${\mathrm{\Delta }h}_y$ influence a lot. The third part (C to D) corresponds to a relatively gentle fall of $F_{obj2}$ for higher values of $F_{obj1}$.

Table 5 indicates that the change of optimization objectives really due to the five design variables, but $F_{obj1}$ and $F_{obj2}$ are conflicting, they have different change trends. Specifically, just in the numerical view, with the decrease of$e_y$, ${\mathrm{\Delta }h}_y$ and $c_b$, and the increase of $e_z$ and ${\mathrm{\Delta }l}_x$, the firing accuracy ($F_{obj1}$) will become poorer but the firing stability ($F_{obj2}$) will be better.

When got these values of the design variables, a optimization towed arillery dynamic model can be obtained by modifying the relative parameters. For the changes of the $e_y$ and $e_y$, usually making a rearrangement of the affix in the breech and a slight change of the mass of the barrel and muzzel brake can reach the goals. For ${\mathrm{\Delta }l}_x$, ${\mathrm{\Delta }h}_y$ and $c_b$, we need to modify the structure dimensions of the cradle and its related adjacent parts, and all these dimension modifications is easy to do. And in flexible multi-body dynamic modeling level, just as description in Section 3.2, we can modify the finite elements models to change the modal neutral files of the related components and parts. In this way a optimization towed arillery dynamic model can be obtained.

Fig. 5Optimal Pareto front

4.2. Selecting a solution of the Pareto front

We got a range of possible optimization results above, and the decision maker can choose a best one for their specific interests. However, there are different ways to select the optimum solution Pareto front.

Max-min criterion is used to select a solution from the Pareto front, in order to show solutions that represent each set of data. This method finds a solution that is equidistant from the ends of each objective. The criterion is defined as follows:

where ${FO}_{1,\mathrm{m}\mathrm{a}\mathrm{x}}$, ${FO}_{1,\mathrm{m}\mathrm{i}\mathrm{n}}$ are the maximum and minimum values of the objective function 1, ${FO}_{2,\mathrm{m}\mathrm{a}\mathrm{x}}$, ${FO}_{2,\mathrm{m}\mathrm{i}\mathrm{n}}$ are the maximum and minimum values of the objective function 2, and ${FO}_{1,i}$, ${FO}_{2,i}$ are the values of the objective function 1 and function 2 for solution $i$ respectively.

This criterion selects from the Pareto front a specific solution (enclosed in box in Fig. 5), whose characteristic and their relative objectives are shown in Table 5. And the comparisons of optimization objectives between the selected solution and the original design are also showed in Table 7 and Fig. 6.

It should be pointed out that our surrogate model only can calculate the parameter values at the moment when the projectile came out of the muzzle, so in the Table 7, we compare the parameter values at that moment in different conditions. While in the Fig. 6, we compare the different parameter values at whole period from firing second to the projectile came out of the muzzle moment, so the selected solution value in surrogate model showed in Fig. 6 just an instantaneous value which are marked as C. And the curve A and curve B indicate the original solution values in ADAMS model and the selected solution values in ADAMS model, respectively.

Fig. 6Comparisons of optimization objectives between the selected solution and the original design: A – optimization objectives in original design of ADAMS model, B – optimization objectives in selected solution of ADAMS model, C – optimization objectives in selected solution of surrogate model

a)${\theta }_y$

b)${\theta }_z$

c)${\tau }_x$

d)${\lambda }_y$

4.3. Discussions

According to the artillery system dynamics [17] and artillery design theory [18], $e_y$, $e_z$ and $h_y$ are three of the most important general parameters for artillery designs, $l_x$ and $c_b$ also play the key roles in firing accuracy and firing stability. To be specific, the vertical and lateral spectral centroid shift $e_y$, $e_z$ of recoiling parts should be as smaller as possible, and trunnion center’s height $h_y$ plays a very important and complex role in firing accuracy and firing stability since the cradle can be rotate slightly around it when firing, its position has key relation to the movement of the artillery because it’s the position of the turning center of the recoiling parts. Similarly, for the gap $c_b$ between cradle bushing and barrel, it also has complex and important effect on firing accuracy and firing stability, it is a complicated contact/collision problem under dynamic impact loading and a thermal phenomenon, many artillery researchers worked on this but it’s hard to decide it in a different situation. For the distance $l_x$ between front and back bushings, since the acting position of the relative sliding between the recoiling parts and cradle are front and back bushings, so theoretically the larger distance $l_x$, the smaller vibration of the muzzle, and that means the better firing accuracy, however, when the distance $l_x$ is too long in some way, it will influence the overturning moment of the whole artillery, and then the displacement of the spade and the base, that means worse firing stability. So a proper distance $l_x$ is needed to balance the firing accuracy and firing stability. In our optimization design, the front cradle bushing’s axial direction offset $\Delta l_x$ is changed to change the distance lx.

Compare with the parameter values in Table 2 and Table 6, it can be found that in our selected solution the vertical and lateral spectral centroid shift $e_y$, $e_z$ of recoiling parts became smaller, the distance $l_x$ and the gap $c_b$ became larger slightly. Combine with artillery system dynamics and artillery design theory, it is not difficult to explain and verify our optimization results. Actually, it can be seen from Fig. 6 that at the end of the curves B (present the selected solution) have smaller values compare to curves A (present the original design), and curves B also get smaller maximum value at the whole interior period, that means our selected solution has a better tactical and technical properties.

What’s more, the curves B have the similar growth pattern with curves A, they are all similar to the classical situations like reference [3] and [13], and it offers another perspective on the validity of our ADAMS simulation model. And the optimization objectives (curves A) in selected solution of ADAMS model are very close to whose at the surrogate model (point C), this indicates our RBF-BP series combine neural network surrogate model has strong generalization ability and high prediction accuracy, it can reflect the artillery launch system well. All these indicate our optimization method effective and it can verify our optimization results well.

To conclude on these comparisons, the selected solution considered all the four objectives at the same time, they all meet an optimal results. In this way, we can get a relatively balanced solution. Moreover, these also claimed our methodology’s feasible and effective.