Modified enhanced collaborative optimization method for artillery launching dynamic model with sequential process

Abstract

The artillery launch system directly influences the muzzle energy and launching accuracy, and therefore it is important to optimize the artillery launch system’s complete process to improve the launch performance. As the objective function of the artillery launch system is non-smooth with coupling parameters in sequential processes, conventional optimization methods are hard to converge for the muti-sequential process of the interaction between the projectile and the barrel. This paper develops a coupled dynamic model for artillery launching, which can predict the performance of the engraving process of the rotating band and the projectile motion in the barrel. The independent optimization problem of the artillery launch system is divided into two subspace problems, and a modified enhanced collaborative optimization (MECO) method with global search capability is proposed, in which the distance criterion and penalty design boundary method are implemented. Results show that the MECO is dedicated not only to satisfying compatibility between coupling parameters of the two sequential processes effectively but also to improving the projectile axial speed at the muzzle and launching accuracy. The MECO maintains a much stable level of convergence than the ECO when the original optimization problem is multimodal.

Introduction

Conventional artillery weapons are still playing an important role in modern warfare. In the past decades, with the development in interior ballistic research, the artillery firepower and accuracy have greatly improved. The interior ballistic system mainly consists of a propellant, a projectile, and a barrel. Artillery launch is achieved via the interaction between the above-mentioned components. The projectile-barrel interaction during launch includes an engraving process of the rotating band and the following in-bore motion of the projectile. During the launching process there exist various coupling parameters in the interior ballistic system which significantly influence the projectile-barrel interaction. Ultimately, they affect the overall launch performance. To systematically improve the launch performance of artillery, a system-integration optimization method is proposed in this paper to optimize the complete interior ballistic system.

Owing to the highly transient nature of the artillery launching process, it is difficult to accurately measure each interior ballistic parameter. Therefore, research on the interior ballistic process of artillery is mainly based on numerical simulation and supplemented by live firing experiments. First of all, as the basis of artillery launch, the propellant combustion model can be used to predict the peak pressure and projectile speed. Different research groups have proposed different models, including the lumped parameter model based on thermodynamics and Lagrangian hypothesis^1,2 and the two-phase flow model.^3,4 Nevertheless, the propellant combustion model cannot acquire the projectile attitude at the muzzle, an initial condition for calculating the exterior ballistic trajectory which can make an important difference to the launching accuracy. In recent years, with the development of three-dimensional simulation software, the combustion model has been embedded in the simulation software of the projectile-barrel contact. This greatly promotes the research of the launch dynamics. Due to different research interests, most researchers usually focus on the engraving process or the projectile motion in the chamber separately. For example, in the study of the engraving process, Keinänen et al. ⁵ set the pressure–time curve subject to the boundary conditions to control the load of the projectile. Li et al. ⁶ used three different methods (FEM, FEM-SPH, and CEL) to simulate the elastic band engraving process. Sun et al. ⁷ studied the effects of the maximum engraving resistance and the engraving pressure under different charge cases. In addition, scholars’ research interests in the engraving process also include the wear of the rifling and the influence of rifling wear degree in interior ballistic performance.^8,9 Since the engraving process is transient, lasting only a few milliseconds, the researchers who study the projectile motion in the barrel usually ignore this process. In the study of the projectile motion in the barrel, Chen ¹⁰ studied the effect of the curved barrel on projectile motion using ABAQUS®. Alexander ¹¹ established the explicit dynamic finite element model (FEM) by ABAQUS®, and calculated the displacement, angular speed, and axial acceleration of the projectile. Yu and Yang ¹² coupled the lumped parameter model with ABAQUS® to analyze the barrel dynamic response under the joint action of propellant gas pressure and projectile contact force. Therefore, the researchers only focus on a specific process of the launch dynamics without taking into full consideration the complete launching process in the chamber.

With the requirements of improving weapon system performance, in recent years, some researchers have used optimization algorithms to solve the artillery launch design problem. In optimizing the propellant combustion model, Gonzalez ¹³ first applied the augmented Lagrange multiplier method to obtain the optimal combustion model of the propellant and formulated the optimal design method. Zhang et al.^14–16 used an intelligent optimization algorithm to optimize the propellant combustion model (lumped parameter model and two-phase flow model) and charge parameters. In the optimization problem with projectile-barrel interaction, studies which have optimized the artillery launch performance only consider part of the launching process. In addition, Wang et al., ¹⁷ Xu and Yang ¹⁸ used the interval uncertainty optimization method to optimize the launching accuracy and stability of the launch system. So far, few studies have focused on the complete optimization of the different stages of the projectile-barrel interaction. It is worth noting that Wei and Wang ¹⁹ utilized Isight® to establish a multidisciplinary design optimization (MDO) model of light weapon with interior, exterior, and terminal ballistics. Despite their primary consideration of the overall performance of the full trajectory rather than of the projectile-barrel interaction, the MDO optimization method has been used to solve the sequential launch problem in the bore, which has offered an effective enlightenment for us.

While the MDO problem in artillery launching is still researched in its infancy, the MDO method has been explored for a long time. The collaborative optimization (CO) strategy, ²⁰ a typical MDO method, has been successfully applied to mathematical test problems and practical engineering problems.^21,22 However, it has also encountered several challenges, as pointed out by DeMiguel and Murry. ²³ In response to the challenges, Sobieski and Kroo ²⁴ and Zadeh et al. ²⁵ suggested using the response surface to approximate the post-optimal behavior of the subject’s sub-problems in the system’s sub-problems, given that representing a non-smooth function with a smooth face may not be an ideal approach. DeMiguel and Murray ²³ proposed the modified collaborative optimization (MCO), which uses an accurate penalty function with fixed penalty parameter values to relax troublesome constraints and adds elastic variables to maintain the smoothness of the problem. Jin et al. ²⁶ presented the design space decrease collaborative optimization and its enhanced version. Besides, Braun et al. ²⁷ suggested relaxing the system’s sub-problem equality constraints to inequalities with a relaxation tolerance, and Li et al.^28,29 promoted this approach by adaptively choosing the tolerance during the solution procedure. But these improved approaches fail to change the underlying CO problem that the system-level Jacobian is singular or that the Lagrangian multiplier is zero in the subspace problem. Roth and Kroo^30,31 reformulated the system-level and subspace problems, proposing the enhanced collaborative optimization (ECO) method. The system-level problem of ECO is an unconstrained minimization problem which is responsible only for the compatibility of shared variables between the system-level and the subspace problems, while the subspace is responsible for minimizing the global objective function. The ECO eliminates the ill-conditioning features associated with the basic CO. Based on Roth’s results,^30,31 the ECO can effectively reduce the amount of calculation by nearly an order of magnitude as compared with the CO. A distinctive feature of the ECO is that it includes linear models of the nonlocal constraints (LMNCs) and a quadratic model of the original objective function in each subspace. Nevertheless, when the LMNCs include non-local variables, the ECO needs to solve another problem, namely the constraint violation minimization (CVM). Tao et al. ³² showed that the LMNCs coincided with higher complexity of the ECO than that of the CO and ATC, and introduced the alternating direction method of multipliers (ADMM) into the ECO to simplify the LMNCs. The results of Tao’s research showed that the ECO-ADMM that simplified the LMNCs could converge with a calculation efficiency between the ECO and CO. In order to ensure the absolute convergence at system level in the ECO, the original optimization function of the subproblem was designed as a quadratic function, but it is usually a highly nonlinear multi-peak function in engineering problem. The local optimal solution may be obtained, when the gradient optimizer is still used for each subproblem. This will result in a system-level failure of the ECO to coordinate the shared variables of each subproblem to achieve consistency, and ultimately optimization failure.

The purpose of this paper is to develop a modified enhanced collaborative optimization (MECO) method to implement the distance criterion and penalty design boundary strategy in order to solve the optimization problem with a two-sequential artillery launching process.

The remaining part of the paper is organized as follows. Section 2 presents the dynamic model of artillery launching coupled with the propellant combustion model. The coupled model can simulate the engraving process and the projectile motion process. Section 3 presents the MECO method and calculation procedure, derives the artillery launch system’s optimization formulation in the form of MECO. Section 4 reports the results. Finally, Section 5 makes further discussion and draws a brief conclusion.

Artillery launching mechanism and model

In this section, we first review and develop the artillery launching mechanism, the propellant combustion model, and a coupled dynamic model of artillery launch. Next, we exhibit the structural parameters and introduce the response indexes and design variables of the artillery launch system. Finally, we describe the optimization problems of artillery launch system to be solved in this paper.

Artillery launching mechanism

The artillery launch system consists of a propellant, a projectile, and a barrel, as illustrated in Figure 1. The combustion of the propellant is confined to the area enclosed by the wall of the barrel and the base of the projectile. During the launching process, the barrel’s bottom is assumed to be fixed and the projectile is pushed forward by the propellant gas pressure. The outer diameter of the rotating band is larger than the inner diameter of the rifling so that the combustion chamber can be sealed. Therefore, the rifling completes the engraving process before guiding the rotation of the projectile, which ensures gyro stability for the projectile during its flight through the air.

Figure 1.

Artillery launching mechanism and parameters.

As shown in Figure 1, there are many coupling parameters of the propellant, projectile, and barrel, which affect not only each other geometrically but also the engraving process of the rotating band and the projectile motion in the barrel. For example, the coupling parameters of the propellant include the propellant thickness $e_1$ and the propellant chamber ${\omega }_0$ . The coupling parameters in the projectile (contain the rotation band) and the barrel include the rifling depth $t_{sh}$ , the bourrelet-to-rifling distance $l_{bd}$ , the mass eccentricity of the projectile $l_c$ , rotating band width $W$ , rotating band location $l_{RB}$ , forcing cone angle $\alpha$ , and groove width $b$ .

The rifling depth $t_{sh}$ is defined in equation (1), where $d_{RB}$ is the outer diameter of the rotating band and $d_G$ is the inner diameter of the rifling

$t_{sh}=\frac{d_{RB}-d_G}{2}$

(1)

And the bourrelet-to-rifling distance $l_{bd}$ is defined in equation (2)

$l_{bd}=\frac{d_G-d_B}{2}$

(2)

where $d_B$ is the diameter of bourrelet.

Therefore, a coupled FEM model, including the propellant combustion model and the dynamic model of the projectile-barrel mechanical interaction, should be developed to study the performance of artillery launching. The two models exchange parameters with each other in their respective computing processes.

Propellant combustion model

In our study, a lumped parameter model is adapted to implement the combustion model as shown in equations (3) and (4). ³³ The combustion model can describe the engraving process and the projectile motion

$\{\begin{array}{c}\psi =\chi Z\left(1+\lambda Z+\mu Z^2\right)\\ \frac{dZ}{dt}=\frac{u_1{p}^n}{e_1}\\ \phi m_p\frac{d{v}_p}{dt}=S_{a}p\\ S_{a}p\left(l_{\psi }+l\right)=f\omega \psi +f_{ig}{\omega }_{ig}-\frac{1}{2}\theta \phi m_p{v}_p^2\end{array}$

(3)

where

$l_{\psi }=l_0\left[1-\frac{\mathit{\Delta }}{{\rho }_p}\left(1-\psi \right)-\alpha \mathit{\Delta }\psi \right]$

(4)

$\psi$ is the percentage of propellant burned; $\chi$ , $\lambda$ and $\mu$ are the shape parameters of the propellant; $Z$ is the relative burned thickness, $u_1$ is the burning rate, and $n$ is the burning rate index; $\phi$ is the secondary power factor, ${\phi }_1$ is the resistance coefficient; $p$ is the powder gas pressure, $m_p$ is the projectile mass, $l$ is the projectile travel, and $v_p$ is the projectile speed in the axial direction of the artillery barrel; $S_a$ is the artillery barrel cross-sectional area, $\theta$ is the heat ratio, $l_0\left(l_{\psi }\right)$ is the chamber (free) volume-to-bore area ratio; $f$ is the propellant force, $\omega$ is the propellant mass, ${\rho }_p$ is the propellant density, $f_{ig}$ is the igniting powder force, ${\omega }_{ig}$ is the igniting powder mass, $e_1$ is the propellant thickness, $\alpha$ is the powder gas volume, and $\mathit{\Delta }$ is the loading density.

Through the above equations, the computing program for the propellant combustion model is developed based on the Runge–Kutta method.

The coupled dynamic model of artillery launch

In this paper, the coupled dynamic model of artillery launch is implemented in ABAQUS®. The projectile-barrel interaction can be figured out by the control function of explicit dynamics, as shown in equation (5)

$\mathbf{M}\mathbf{Ü}+\mathbf{C}\dot{\mathbf{U}}+\mathbf{K}\mathbf{U}=\mathbf{F}$

(5)

where $\mathbf{U}$ , $\dot{\mathbf{U}}$ , and $\ddot{\mathbf{U}}$ are the matrices of displacement, velocity, and acceleration, respectively, in the FEM; $\mathbf{M}$ , $\mathbf{C}$ , and $\mathbf{K}$ are the matrices of mass, damping, and stiffness, respectively; $\mathbf{F}$ is the load vectors.

In the artillery launch FEM, the projectile base pressure drives the projectile to move in the bore. The function relationship between the projectile base pressure $p_d$ and the average pressure $p$ can be calculated by the propellant combustion model, as expressed in equation (6)

$p_d=\frac{3\phi m}{\omega +3\phi m}p$

(6)

The coupling between the FEM and the propellant combustion model is realized through the VUAMP subroutine in ABAQUS®. First, the VUAMP coupled with the propellant combustion model calculates the instantaneous projectile base pressure $p_d\left(t\right)$ and records the instantaneous relative burned thickness $Z\left(t\right)$ . Next, the FEM of artillery launch exerts the pressure $p_d\left(t\right)$ on the projectile, calculates the position $l_0\left(t+\mathit{\Delta }t\right)$ and axial speed $v_x\left(t+\mathit{\Delta }t\right)$ of the projectile after the time increment $\mathit{\Delta }t$ , and sends the position and speed to the VUAMP. This step loops until the end of the engraving process or the projectile exits the muzzle. One distinction between the engraving process and the projectile movement in the barrel is that the engraving process needs a failure model for the rotating band in the FEM, as shown in Figure 2.

Figure 2.

The coupled dynamic model of the artillery launch.

Response indexes and design variables of the coupled dynamic model

Through the coupled dynamic model of artillery launch introduced in the previous section, the response indexes can be calculated from the engraving process of the rotating band and the movement process of the projectile in the barrel, respectively. These two processes are named “sub a” and “sub b,” respectively.

The engraving process of the rotating band calculated one response parameter: the engraving resistance $f_R$ , and the movement process of the projectile in the barrel calculated six response parameters: the maximum chamber pressure $p_m$ , the projectile axial speed at the muzzle $v_{xm}$ , the projectile’s lateral angular displacement at the muzzle ${\mathit{\Omega }}_{ym}$ , the vertical angular displacement at the muzzle ${\mathit{\Omega }}_{zm}$ , the projectile’s lateral angular speed at the muzzle ${\mathit{\Phi }}_{ym}$ , and the vertical angular speed at the muzzle ${\mathit{\Phi }}_{zm}$ .

Many scholars have studied the parameters or sensitivity that affects artillery firing performance in previous literatures.^3,8,16,34 In this paper, 11 high-sensitivity parameters are selected as the design variables. According to the description of the coupled dynamic model, five design variables are grouped as shared variables $\mathbf{z}=\left[b t_{sh} W e_1 \omega \right]$ , where $b$ is the groove width, $t_{sh}$ is the rifling depth, $W$ is the rotating band width, $e_1$ is the propellant thickness, and $\omega$ is the propellant mass. The rest design variables are distributed to the two subspaces as local variables, such as the subspace of rotating band engraving has two local variables ${\mathit{x}}_{La}=\left[\alpha d_{RB}\right]$ , and the subspace of projectile motion has four local variables ${\mathit{x}}_{Lb}=\left[l_c l_{RB} l_{bd} W_0\right]$ , where $\alpha$ is the forcing cone angle, $d_{RB}$ is the outer diameter of the rotating band, $l_c$ is the mass eccentricity of the projectile, $l_{RB}$ is the rotating band location, $l_{bd}$ is the bourrelet-to-rifling distance, and $W_0$ is the propellant volume.

Problem description of the artillery launch system

Usually, artillery launch aims to maximize the muzzle energy and firing accuracy in safety.

Global objective function

To be distinguished from the objective function in the MDO problem in Section 3, the “global objective function” defined by Roth and Kroo ³⁰ is used to represent the overarching design goal.

There are many indicators to evaluate artillery launch. In this paper, both the projectile axial speed normalization coefficient ${\eta }_1$ and the projectile disturbance coefficient ${\eta }_2$ at the muzzle are considered as the global objective functions. Obviously, this is a multi-objective optimization (MOO) problem. The linear weighted sum method ³⁵ is used to convert these two evaluation indicators into a single global objective function which can be calculated by equation (7)

$\begin{array}{c}{F}_0={\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2\\ ={\lambda }_1\frac{{\varphi }_{v_{xm}}}{v_{xm}}+{\lambda }_2\left(\frac{abs\left({\mathrm{\Omega }}_{ym}\right)}{{\varphi }_{{\mathit{\Omega }}_{ym}}}+\frac{abs\left({\mathrm{\Omega }}_{zm}\right)}{{\varphi }_{{\mathit{\Omega }}_{zm}}}+\frac{abs\left({\mathrm{\Phi }}_{ym}\right)}{{\varphi }_{{\mathit{\Phi }}_{ym}}}+\frac{abs\left({\mathrm{\Phi }}_{zm}\right)}{{\varphi }_{{\mathit{\Phi }}_{zm}}}\right)\end{array}$

(7)

where $v_{xm}$ is the projectile axial speed at the muzzle, which is a set of response indexes, ${\varphi }_{\ast }$ denotes the regularization factor, ${\lambda }_1$ and ${\lambda }_2$ are the weight coefficient, the weight coefficients are ${\lambda }_1=0.6$ and ${\lambda }_2=0.1$ , which were calculated by the method of analytic hierarchy process (AHP), ³⁶ ${\mathit{\Omega }}_{ym}$ and ${\mathit{\Omega }}_{zm}$ are the angular displacements of the projectile moving in lateral and vertical at the muzzle, respectively; ${\mathit{\Phi }}_{ym}$ and ${\mathit{\Phi }}_{zm}$ are the angular velocities of the projectile moving in lateral and vertically at the muzzle, respectively; and the regularization factors set as below: ${\varphi }_{v_{xm}}=950$ , ${\varphi }_{\mathrm{\Omega }ym}=4.75\times {10}^{-5}$ , ${\varphi }_{\mathrm{\Omega }zm}=1.95\times {10}^{-3}$ , ${\varphi }_{\mathrm{\Phi }ym}=\mathbf{ }6\mathbf{.}5$ , ${\varphi }_{\mathrm{\Phi }zm}=\mathbf{ }4\mathbf{.}8$ ; $abs(\ast )$ denotes the absolute value function.

Constraints

To ensure launching safety and service life, three constraints should be considered for the artillery launch problem.

(1) Engraving resistance

From the perspective of launching safety, it is necessary to limit the engraving resistance $f_R$ ∈ [900, 1600] kN ⁷ during the engraving process. The constraint function is expressed in equation (8)

$g_{fR}=350-\left|f_R\left(b,t_{sh},W,e_1,\omega ,\alpha ,d_{RB}\right)-1250\right|\ge 0$

(8)

(2) Maximum chamber pressure

The maximum chamber pressure $p_m$ , which occurs during the movement of the projectile in the bore, should be limited to no more than 320 MPa. ¹⁷ The constraint function is expressed in equation (9)

$g_{pm}=320-p_m\left(b,t_{sh},W,e_1,\omega ,l_c,l_{RB},l_{bd},W_0\right)\ge 0$

(9)

(3) Ultimate life of artillery

The ultimate life of artillery is defined as the time it takes for the projectile’s axial speed at the muzzle to drop by less than 10% of the initial design value when the rifling is worn. ³⁷ The design variables $l_c$ and $t_{sh}$ are considered subject to rifling wear in this paper. According to Ma, ³⁸ the function relationship between the amount of rifling wear $\mathit{\Delta }{t}_{sh}$ and the number of shots $N$ from the experimental data is defined in equation (10), and the rifling depth $t_{sh}^{\prime }$ upon reaching the ultimate life can be calculated by equation (11)

$\mathit{\Delta }{t}_{sh}=\frac{1}{1.322}\mathit{ln}\frac{N}{116.74}$

(10)

$t_{sh}^{\prime }=t_{sh}-\mathit{\Delta }{t}_{sh}=t_{sh}-\frac{1}{1.322}\mathit{ln}\frac{N}{116.74}$

(11)

From equations (1) and (11) the bourrelet-to-rifling distance at the limit life $l_{bd}^{\prime }$ is defined by equation (12)

$l_{bd}^{\prime }=\frac{d_{RB}-d_B}{2}+t_{sh}-\frac{1}{1.322}\mathit{ln}\frac{N}{116.74}$

(12)

The constraint function $\mathit{\Delta }{v}_{xm}$ is the difference between the initial axial speed of the projectile at the muzzle $v_{xm}$ and the projectile axial speed at the muzzle when reach the limit life $v_{lxm}$ and it can be calculated by equation (13)

$g_{\mathrm{\Delta }vm}={\lambda }_p{p}_m\left(b,t_{sh},W,e_1,\omega ,\mathbf{ }{l}_c,l_{RB},l_{bd},W_0\right)-\mathit{\Delta }{p}_m\left(b,t_{sh}^{\prime },W,e_1,\omega ,\mathbf{ }{l}_c,l_{RB},l_{bd}^{\prime },W_0\right)\ge 0$

(13)

where $t_{sh}^{\prime }$ is the rifling depth at the limit life; $l_{bd}^{\prime }$ is the bourrelet-to-rifling distance at the limit life; $v_{lxm}$ is the projectile axial speed at the muzzle at the limit life, which can be calculated by the coupled dynamic model of artillery launch in chapter 2 with the parameters $t_{sh}^{\prime }$ and $l_{bd}^{\prime }$ ; and ${\lambda }_p$ is the life coefficient.

The US military has verified the life of the M1 artillery with the results showing that the ultimate life of the artillery is 15,000 rounds. ³⁷ Considering the boundary of the surrogate model, the amount of wear and the number of shots need to be redetermined. Therefore, the projectile’s axial speed at the muzzle upon reaching the ultimate life should not exceed by the life coefficient ${\lambda }_p$ = 1.11% by equation (11) when the number of shots N = 200.

Methodology

In this section, the ECO method is reviewed. Next, the MECO method is proposed for the problem of the non-smooth global objective function. Finally, the optimization problem with a two-sequential artillery launching process is derived in the form of MECO.

Overview of the ECO method

The ECO solves the MDO problem with a bi-level structure, that is, a system-level problem and $n$ subspace problems. The system-level of ECO is an unconstrained minimized problem and does not have the global objective function $F_0$ . The system-level is only responsible for the compatibility of shared design variables between the system level $\mathbf{z}$ and the subspace ${\mathit{x}}_i$ . The system-level problem is expressed in equation (14)

$\begin{array}{rr}\hfill \mathrm{F}\mathrm{i}\mathrm{n}\mathrm{d}& \mathbf{z}\hfill \\ \hfill \underset{\mathbf{z}}{\mathbf{min}}& J_{sys}={\displaystyle \sum }_{i=1}^n\Vert \mathbf{z}-{\mathit{x}}_i{\Vert }_2^2\hfill \\ \hfill \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathbf{\hspace{0.5em}}\mathrm{t}\mathrm{o}& \mathrm{N}\mathrm{o}\mathbf{\hspace{0.5em}}\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\hfill \end{array}$

(14)

In contrast with the system-level problem, the subspace problems optimize the global objective functions. The objective function of the subspace consists of three parts: the global objective function, a quadratic measure of compatibility ${\lambda }_C(\Vert \mathbf{z}-{\mathit{x}}_i{\Vert }_2^2)$ , and a set of slack variables ${\lambda }_F{\displaystyle \sum \mathit{s}}$ . The constraint functions of subspace include the local constraints and linear models of the nonlocal constraints (LMNCs). The $i^{th}$ subspace problem is formulated in equation (15)

$\begin{array}{rr}\hfill \mathrm{F}\mathrm{i}\mathrm{n}\mathrm{d}& \mathbf{x}=\left[{\mathit{x}}_i\mathbf{ }{\mathit{x}}_{Li}\mathbf{ }{\mathit{s}}_j\right]\hfill \\ \hfill \underset{\mathrm{x}}{\mathbf{min}}& J_{subi}={\tilde{F}}_0\left({\mathit{x}}_i,{\mathit{x}}_{Li}\right)+{\lambda }_{Ci}(\parallel \widehat{\mathbf{z}}-{\mathit{x}}_{\mathit{i}}{\parallel }_{\mathit{2}}^{\mathit{2}})+{\lambda }_{Fi}{{\displaystyle \sum }}_{j=1,j\ne i}^n{{\displaystyle \sum }}_{k=1}^{n_s}{s}_{jk}\hfill \\ \hfill \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathbf{\hspace{0.5em}}\mathrm{t}\mathrm{o}& {\mathbf{g}}_i\left({\mathit{x}}_i,{\mathit{x}}_{Li}\right)\ge \mathbf{0}\hfill \\ \hfill & {\tilde{\mathbf{g}}}_j\left({\mathit{x}}_i\right)+{\mathit{s}}_j\ge \mathbf{0},\hspace{0.17em}j=\left\{1,\cdots ,n\right\},\hspace{0.17em}j\ne i\hfill \\ \hfill & {\mathbf{s}}_j\ge \mathbf{0},\hspace{0.17em}j=\left\{1,\cdots ,n\right\},\hspace{0.17em}j\ne i\hfill \end{array}$

(15)

where ${\tilde{F}}_0({\mathit{x}}_i,{\mathit{x}}_{Li})$ is a quadratic model of the global objective function, ${\mathit{x}}_i$ is the shared design variables in subspace $i$ , ${\mathit{x}}_{Li}$ is the local design variables in subspace $i$ , $\widehat{\mathbf{z}}$ is a copy of shared variables provided by the system-level, ${\mathit{g}}_i({\mathit{x}}_i,{\mathit{x}}_{Li})$ is the local constraint function in subspace $i$ , ${\tilde{\mathit{g}}}_j\left({\mathit{x}}_i\right)$ is the linear models of nonlocal constraint (LMNCs), as formulated in equation (16), also means that the constraints used in subspace $j$ , ${\mathit{s}}_j$ is the vector of slack variables for ${\tilde{\mathit{g}}}_j\left({\mathit{x}}_i\right)$ , ${\lambda }_{Ci}$ and ${\lambda }_{Fi}$ are the penalty parameters for the consistency and the slack variables feasible, respectively, $n_s$ is the number of LMNCs

${\tilde{g}}_j\left(x_i\right)={g_j|}_{\widehat{\mathbf{z}},{\mathit{x}}_L^{\ast }}+{\displaystyle \sum }_{i=1}^{n_s}\left[{\frac{\partial g_j}{\partial x_i}|}_{\widehat{\mathbf{z}},{\mathit{x}}_L^{\ast }}\left(x_i-\widehat{z}\right)\right]$

(16)

where $\mathbf{ }\partial g_j/\partial x_i$ is the constraint model coefficient. The constraint model coefficients are obtained through the central difference method and pre-stored in the program in this study.

Since the parameters ${\tilde{\mathit{x}}}_{Li}$ are used in the constraint function ${\tilde{\mathit{g}}}_j\left({\mathit{x}}_i\right)$ used in each subspace contains, the parameters ${\tilde{\mathit{x}}}_{Li}$ are the copy of local design variables provided by the constraint violation problem (CVM) and cannot be directly called from the subspace $i$ . The CVM is expressed in equation (17)

$\begin{array}{rr}\hfill \mathrm{F}\mathrm{i}\mathrm{n}\mathrm{d}& \mathrm{x}=\left[x_{Li}\mathbf{ }{s}_j\right]\hfill \\ \hfill \underset{\mathrm{x}}{\mathbf{min}}& J_{CVM\mathbf{ }i}={\displaystyle \sum }_{j=1}^{n_s}{s}_j\hfill \\ \hfill \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathbf{\hspace{0.5em}}\mathrm{t}\mathrm{o}& {\mathbf{g}}_i\left(x_{Li},s_j\right)\ge \mathbf{0}\hfill \\ \hfill & {\mathit{s}}_i\ge \mathbf{0}\hfill \end{array}$

(17)

where $n_s$ is the numbers of the constraint ${\mathit{g}}_i(x_{Li},s_j)$ in the subspace $i$ .

MECO method

One of the distinctive features in ECO is that a quadratic model of the global objective function ${\tilde{F}}_0(\cdot )$ is used in the subspace for absolute convergence. Tao et al. ³² reported that designers had the option to replace the quadratic model ${\tilde{F}}_0(\cdot )$ with the true model $F_0(\cdot )$ when the true global objective function $F_0(\cdot )$ can be evaluated directly by subspaces. But the global objective function $F_0(\cdot )$ is usually non-smooth, especially in the artillery launch problem to be solved in this paper. If the gradient optimizer remains used in each subspace, the extremely different local solutions will be obtained.

The intelligent optimization algorithms have been widely used in the global optimization for non-smooth functions. ³⁹ Therefore, the intelligent optimizer is used in the subspaces to solve the non-smooth global objective problem and hence the artillery launching problem. In fact, the simple combination of the intelligent optimizer and the ECO may not converge in the system level. This is because the subspace uses the target values provided at the system level as the starting point for optimization when the gradient optimizer is applied in the ECO while the intelligent optimizer uses a group of random populations.

The ECO-based MECO proposed in this study has two distinctive characteristics: (1) A subspace loop is added by the distance criterion method to choose the optimal solution and (2) the penalty design variable boundary is used to improve the convergence.

The distance criterion method for the subspace loop

Based on the ECO strategy, a loop is added for each subspace and a distance criterion is adopted to choose the optimal solution to maintain the stability of the subspace in the main loop. The distance criterion is used to estimate the optimal solution’s Euclidean distance between the subspaces and the system level, and to compare it with the previous distance to decide whether the subspace needs to be reoptimized. The criterion is formulated in equation (18)

$D=\{\begin{array}{cc}1,\hfill & k=1\hfill \\ \parallel \mathbf{z}-{\mathit{x}}_k{\parallel }_2,\hfill & k>1\hfill \end{array}$

(18)

where ${\mathit{x}}_k$ is the optimal solution in the subspace at the $k^{\mathbf{th}}$ loop, and $D$ is the Euclidean distance.

The penalty design variable boundary method

Since the only system-level responsibility is to handle variable consistency among all the subspaces, the system-level guidance direction tends to deviate from the previous cycle when the optimal solutions in each subspace are extremely different. Keeping the system-level guidance stable is the key to ensuring convergence.

We set the penalty design variable boundary of the subspaces to maintain the stability of orientation after the system-level direction is altered. The criterion $\mathit{\Delta }d$ introduced in this paper is defined as the ratio of the system-level objective function value to that in the previous loop. When $\mathit{\Delta }d$ > $\epsilon$ , the system-level direction is altered by the optimal solutions of a subspace and the design variable boundary to the system-level solution should be penalized, as shown in equation (19)

$\{\begin{array}{c}{x}_k^B=x_k^B, \mathit{\Delta }{d}_j\le \epsilon \\ x_k^B=x_k^z, \mathit{\Delta }{d}_j>\epsilon \end{array}$

(19)

where $\mathit{\Delta }{d}_j=d_{kj}/d_{k(j-1)}$ , $j$ is the system-level iteration numbers, and $j\ge 2$ , $d_{kj}=\Vert z_k-x_k{\Vert }_2$ is the distance of the kth design variable, $B$ denotes the boundary of design variables, $z$ is the system-level value, and $\epsilon$ is the threshold, taken as 1.5 in this study.

Calculation procedure of the MECO

Compared with the ECO, it is obvious that the subspace and system-level loop are modified in the MECO. The flowchart of the MECO is shown in Figure 3. The main steps of the MECO method are outlined as follows:

Figure 3.

Flowchart of MECO.

Step 1: Initialize/update the parameters or design variables and their boundaries.

Step 2: Solve the constraint violation model (CVM) for the non-local variables of the subspace to get the constraint model coefficients (CMCs) of each subspace by the gradient optimizer.

Step 3: Obtain the solutions of the subspace by the intelligent optimizer, such as the generational particle swarm optimization (GPSO). Calculate the space distance $D_i$ by equation (18). Set the initial distance $D_0=D$ for the first loop; otherwise, continue the subspace loop until $D=D_0$ or the end condition is satisfied.

Step 4: Update the design variables of each subspace after all subspaces are completed.

Step 5: Obtain the system-level feasible solution by the gradient optimizer and calculate the distance $\mathit{\Delta }{d}_j$ of each design variable by equation (19) and figure out the new design variable boundaries by penalty design variable boundary method of equation (19) for $\mathit{\Delta }{d}_j$ < $\epsilon$ .

Step 6: Check whether the system level satisfies the convergence tolerance or reaches the maximum iterations. If it does, end the iteration and obtain the optimal solution and global objectives; else substitute $T=T+1$ , update the system-level shared variables and all boundaries, and return to Step 1.

MECO-based optimization model of artillery launch

The optimization problem with a two-sequential artillery launch is divided into two subspace problems. The form of MECO formulation consists of five parts: the CVM of sub a, the CVM of sub b, the optimization problem of sub a, the optimization problem of sub b, and the system-level optimization problem.

The CVM a for parameters ${\mathit{x}}_{La}=[\alpha \hspace{0.17em}{d}_{RB}]$ is expressed in equation (20)

$\begin{array}{rr}\hfill \mathrm{Find}& \mathbf{x}=[\alpha \hspace{0.17em}{d}_{RB}\hspace{0.17em}{s}_1]\hfill \\ \hfill \underset{\mathbf{x}}{\min }& J_{CVM\hspace{0.17em}a}=s_1\hfill \\ \hfill \mathrm{subject}\hspace{0.17em}\mathrm{to}& g_{a1}=350-|f_R(\widehat{b},{\widehat{t}}_{sh},\widehat{W},{\widehat{e}}_1,\widehat{\omega },\alpha ,d_{RB})-1250|+s_1⩾0\hfill \\ & s_1⩾0\hfill \end{array}$

(20)

where $J_{CVMa}$ is the CVM objective function in the subspace $a$ , and $\widehat{\ast }$ means it is a parameter provided by the system-level.

And the CVM b for parameters ${\mathit{x}}_{\mathit{Lb}}=\left[l_c{l}_{RB}{l}_{bd}{W}_0\right]$ is shown in equation (21)

$\begin{array}{rr}\hfill \mathrm{Find}& \mathbf{x}=\left[l_c\hspace{0.17em}{l}_{RB}\hspace{0.17em}{l}_{bd}\hspace{0.17em}W\hspace{0.17em}{s}_1\hspace{0.17em}{s}_2\right]\hfill \\ \hfill \underset{\mathbf{x}}{\min }& J_{CVM\hspace{0.17em}\hspace{0.17em}b}=s_1+s_2\hfill \\ \hfill \mathrm{subject}\hspace{0.17em}\mathrm{to}& g_{b1}=320-p_m\left(\widehat{b},{\widehat{t}}_{\mathit{sh}},\widehat{W},{\widehat{e}}_1,\widehat{\omega },l_c,l_{RB},l_{bd},W_0\right)+s_1⩾0\hfill \\ & g_{b2}={\lambda }_p{p}_m\left(\widehat{b},{\widehat{t}}_{\mathit{sh}},\widehat{W},{\widehat{e}}_1,\widehat{\omega },l_c,l_{RB},l_{bd},W_0\right)\hfill \\ & -{\mathit{\Delta }}_{pm}\left(\widehat{b},t_{sh}^{\prime },\widehat{W},{\widehat{e}}_1,\widehat{\omega },l_c,l_{RB},l_{bd}^{\prime },W_0\right)+s_2⩾0\hfill \\ & s_1,s_2⩾0\hfill \end{array}$

(21)

where $J_{CVMa}$ is the CVM objective function in the subspace $b$ , $\widehat{\ast }$ means it’s a parameter provided by the system-level.

The optimization problem for the subspace of rotating band engraving process (sub a) is formulated in equation (22)

$\begin{array}{rr}\hfill \mathrm{Find}& \mathbf{x}=\left[{\mathit{x}}_a\hspace{0.17em}{\mathit{x}}_{La}\hspace{0.17em}{\mathit{s}}_a\hspace{0.17em}\right]=\left[b\hspace{0.17em}{t}_{sh}\hspace{0.17em}W\hspace{0.17em}{e}_1\hspace{0.17em}\omega \hspace{0.17em}\alpha \hspace{0.17em}{d}_{RB}\hspace{0.17em}{s}_{a1}\hspace{0.17em}{s}_{a2}\right]\hfill \\ \hfill \underset{\mathbf{X}}{\min }& J_{sub\hspace{0.17em}a}=F_0+{\lambda }_c\left({\left(b-\widehat{b}\right)}^2+{\left(t_{sh}-{\widehat{t}}_{sh}\right)}^2+{\left(W-\widehat{W}\right)}^2+{\left(e_1-{\widehat{e}}_1\right)}^2+{\left(\omega -\widehat{\omega }\right)}^2\right)\hfill \\ & +{\lambda }_F\left(s_{a1}+s_{a2}\right)\hfill \\ \hfill \mathrm{subject}\hspace{0.17em}\mathrm{to}& g_{a1}=350-|f_R\left(b,t_{sh},W,e_1,\omega ,\alpha ,d_{RB}\right)-1250|⩾0\hfill \\ & {\tilde{g}}_{b1}=320-p_m\left(b,t_{sh},W,e_1,\omega ,{\widehat{l}}_c,{\widehat{l}}_{RB},{\widehat{l}}_{bd},{\widehat{W}}_0\right)+\frac{\partial g_{b1}}{\partial b}\left(b-\widehat{b}\right)+\frac{\partial g_{b1}}{\partial t_{sh}}\left(t_{sh}-{\widehat{t}}_{sh}\right)\hfill \\ & +\frac{\partial g_{b1}}{\partial \omega }\left(\omega -\widehat{\omega }\right)+\frac{\partial g_{b1}}{\partial e_1}\left(e_1-{\widehat{e}}_1\right)+\frac{\partial g_{b1}}{\partial \omega }\left(\omega -\widehat{\omega }\right)+s_{a1}⩾0\hfill \\ & {\tilde{g}}_{b2}={\lambda }_p{p}_m\left(b,t_{sh},W,e_1,\omega ,{\widehat{l}}_c,{\widehat{l}}_{RB},{\widehat{l}}_{bd},{\widehat{W}}_0\right)-\mathit{\Delta }{p}_m\left(b,t_{\mathit{sh}}^{\prime },W,e_1,\omega ,{\widehat{l}}_c,{\widehat{l}}_{RB},l_{bd}^{\prime },{\widehat{W}}_0\right)\hspace{0.17em}+\frac{\partial g_{b2}}{\partial b}(b-\widehat{b})+\frac{\partial g_{b2}}{\partial t_{sh}}(t_{sh}-{\widehat{t}}_{sh})+\frac{\partial g_{b2}}{\partial W}(W-\widehat{W})+\frac{\partial g_{b2}}{\partial e_1}(e_1-{\widehat{e}}_1)\hfill \\ & +\frac{\partial g_{b2}}{\partial \omega }\left(\omega -\widehat{\omega }\right)+s_{a2}⩾0\hfill \\ & s_{a1}⩾0,s_{a2}⩾0\hfill \end{array}$

(22)

where $J_{sub\hspace{0.17em}a}$ is the subspace objective function in the subspace $a$ , $\widehat{\ast }$ means it’s a parameter, provided by the system-level or CVM b, ${\lambda }_C$ is the penalty parameters for the consistency and ${\lambda }_F$ is the penalty parameters for slack variables feasible.

And the optimization problem for the subspace of projectile motion (sub b) is formulated in equation (23)

$\begin{array}{rr}\hfill \hspace{0.17em}\mathrm{Find}& \mathbf{x}=\left[{\mathit{x}}_b\hspace{0.17em}{\mathit{x}}_{Lb}\hspace{0.17em}{\mathit{s}}_b\right]\hspace{0.17em}=\left[b\hspace{0.17em}{t}_{sh}\hspace{0.17em}W\hspace{0.17em}{e}_1\hspace{0.17em}\omega \hspace{0.17em}{l}_c\hspace{0.17em}{l}_{RB}\hspace{0.17em}{l}_{bd}\hspace{0.17em}{W}_0{s}_{b1}\right]\hfill \\ \hfill \hspace{0.17em}\underset{\mathbf{X}}{\min }& J_{sub\hspace{0.17em}b}=F_0+{\lambda }_c\left({\left(b-\widehat{b}\right)}^2+{\left(t_{sh}-{\widehat{t}}_{sh}\right)}^2+{\left(W-\widehat{W}\right)}^2+{\left(e_1-{\widehat{e}}_1\right)}^2+{\left(\omega -\widehat{\omega }\right)}^2\right)\hfill \\ & +{\lambda }_F{s}_{b1}\hfill \\ \hfill \mathrm{subject}\hspace{0.17em}\mathrm{to}& g_{b1}=320-p_m\left(b,t_{sh},W,e_1,\omega ,l_c,\hspace{0.17em}{l}_{RB},\hspace{0.17em}{l}_{bd},\hspace{0.17em}{W}_0\right)⩾0\hfill \\ & g_{b2}={\lambda }_p{p}_m\left(b,t_{sh},W,e_1,\omega ,l_c\hspace{0.17em},l_{RB}\hspace{0.17em},l_{bd},\hspace{0.17em}{W}_0\right)\hfill \\ & -\mathit{\Delta }{p}_m\left(b,t{\prime }_{sh},W,e_1,\omega ,l_c,\hspace{0.17em}{l}_{RB},\hspace{0.17em}{l}_{bd}^{\prime }\hspace{0.17em}{W}_0\right)⩾0\hfill \\ & {\tilde{g}}_{a1}=350-|f_R\left(b,t_{sh},W,e_1,\omega ,\widehat{\alpha },{\widehat{d}}_{RB}\right)-1250|+\frac{\partial g_{a1}}{\partial b}\left(b-\widehat{b}\right)\hspace{0.17em}+\frac{\partial g_{a1}}{\partial t_{sh}}\left(t_{sh}-{\widehat{t}}_{sh}\right)+\frac{\partial g_{a1}}{\partial W}\left(W-\widehat{W}\right)+\frac{\partial g_{a1}}{\partial e_1}\left(e_1-{\widehat{e}}_1\right)\hspace{0.17em}+\frac{\partial g_{a1}}{\partial \omega }\left(\omega -\widehat{\omega }\right)+s_{b1}⩾0\hfill \\ & s_{b1}⩾0\hfill \end{array}$

(23)

where $J_{sub\hspace{0.17em}\hspace{0.17em}b}$ is the subspace objective function in the subspace $b$ , $\widehat{\ast }$ means it is a parameter provided by the system-level or CVM a, ${\lambda }_C$ is the penalty parameters for the consistency and ${\lambda }_F$ is the penalty parameters for slack variables feasible.

System-level problem is expressed in equation (24)

$\begin{array}{rr}\hfill \mathrm{Find}\hspace{0.17em}\mathbf{z}& \mathbf{z}=\left[b\hspace{0.17em}{t}_{\mathit{sh}}\hspace{0.17em}W\hspace{0.17em}{e}_1\hspace{0.17em}\omega \right]\hfill \\ \hfill \underset{\mathbf{X}}{\min }& J_{sys}={\left(b-{\widehat{b}}_a\right)}^2+{\left(t_{\mathit{sh}}-{\widehat{t}}_{\mathit{s}\mathit{ha}}\right)}^2+{\left(W-{\widehat{W}}_a\right)}^2+{\left(e_1-{\widehat{e}}_{1a}\right)}^2+{\left(\omega -{\widehat{\omega }}_a\right)}^2\hfill \\ & \hspace{0.17em}+{\left(b-{\widehat{b}}_b\right)}^2+{\left(t_{\mathit{sh}}-{\widehat{t}}_{s\mathit{hb}}\right)}^2+{\left(W-{\widehat{W}}_b\right)}^2+{\left(e_1-{\widehat{e}}_{1b}\right)}^2+{\left(\omega -{\widehat{\omega }}_b\right)}^2\hfill \\ \hfill \mathrm{subject}\hspace{0.17em}\mathrm{to}& \mathrm{No}\hspace{0.17em}\mathrm{constraints}\hfill \end{array}$

(24)

where $J_{sys}$ is the system-level objective function, $\widehat{\ast }a$ and $\widehat{\ast }b$ are parameters provided by the subspace a and subspace b, respectively.

Results

In this paper, the complete optimization for the two processes of the projectile-barrel interaction is realized by implementing the MECO method with global search capability.

First, the design variables and boundaries of the subspace of the engraving process and the subspace of the projectile movement process in the barrel are listed in Tables 1 and 2, respectively.

Table 1.

Optimal solutions of shared design variables and boundaries.

Variable	Groove width $b$ [mm]	Rifling depth $t_{sh}$ [mm]	Rotating band width $W$ [mm]	Propellant thickness $e_1$ [mm]	Propellant mass $\omega$ [kg]
Initial	6.33	2.33	55.52	2.40	17.40
Lower	6.10	2.33	54.00	2.10	17.40
Upper	6.33	2.80	62.00	2.40	18.00
Optimal	6.33	2.33	61.98	2.39	17.58

Table 2.

Optimal solutions of local variables for the subspaces and boundaries.

Variable	Subspace of rotating band engraving		Subspace of projectile motion
	Forcing cone angle $\alpha$ [deg]	Outer diameter of the rotating band $d_{RB}$ [mm]	Mass eccentricity of the projectile $l_c$ [mm]	Rotating band location $l_{RB}$ [mm]	Bourrelet-to-rifling distance $l_{bd}$ [mm]	Propellant chamber $W_0$ [mm3]
Initial	2.25	81.00	0.25	182.00	0.25	25.02
Lower	1.78	79.80	0.00	172.00	0.10	25.00
Upper	5.71	81.00	0.30	192.00	0.60	27.00
Optimal	2.40	80.78	0.14	174.00	0.12	26.26

The Latin hypercube sampling (LHS) method is applied for DOEs. Each indicator acquires 35 training sets data (15 test sets data) from the model of the engraving process, and 70 training sets data (30 test sets data) from the projectile movement process in the barrel.

The surrogate models of BP neural network are obtained by implementing the TensorFlow® training of the training sets data, and the accuracy of the surrogate models is validated by determination coefficient $R^2$ of the test sets. The accuracy of the surrogate models is listed in Table 3. It is shown that the values of determination coefficient $R^2$ are all greater than 0.90.

Table 3.

Accuracy of the surrogate models.

Response	Engraving resistance $f_R$	Projectile speed at the muzzle $v_{xm}$	Maximum chamber pressure $p_{xm}$	Projectile lateral angular displacement at the muzzle ${\mathrm{\Omega }}_{ym}$	Projectile vertical angular displacements at the muzzle ${\mathrm{\Omega }}_{zm}$	Projectile lateral angular speed at the muzzle ${\mathrm{\Phi }}_{\mathit{ym}}$	Projectile vertically angular speed at the muzzle ${\mathrm{\Phi }}_{\mathit{zm}}$
$R^2$	0.9587	0.9737	0.9795	0.9328	0.9269	0.9347	0.9216

To realize the complete optimization for the subspace engraving process and subspace projectile movement process, the MECO strategy is implemented in Python 3. The sequential least squares programming (SLSQP) and the generational particle swarm optimization (GPSO) from pygmo2 are used as the gradient optimizer and the global optimizer, respectively. pygmo2 is a scientific library providing a large number of optimization problems and algorithms. ⁴⁰

We set ${\lambda }_F=10$ and ${\lambda }_C=0.1$ ,^30–32 the system-level consistency tolerance ${\epsilon }_0=1\times {10}^{-6}$ , and the threshold for the penalty design variable boundary $\epsilon =1.5$ , and computation is performed by the ECO method, the method of distance criterion with ECO, and the MECO method, respectively. The system-level convergence chart is shown in Figure 4.

Figure 4.

Convergence chart of system-level objective function values.

Under the implementation of the MECO method, the convergence curve of shared design variables is shown in Figure 5, which also displays the penalization process of the boundary of some design variables and the dynamic variation range of design variables in the subspace.

Figure 5.

Convergence curves and boundaries of shared design variables.

The values of local constraints with the iteration from the MECO method are listed in Figure 6.

Figure 6.

Values of the three system-level constraint functions with MECO.

The indicators and global objective function by the MECO method are shown in Figure 7.

Figure 7.

Convergence curves of two indicators and the system-level global objective function with MECO.

Finally, the optimal design variables (including the shared variables and local variables of each subspace) are brought into the coupled dynamic model of artillery launch introduced in Section 2. The original and optimal speeds are compared in Table 4.

Table 4.

Projectile axial speeds at the muzzle $v_{xm}$ before and after optimization.

Initial value [m/s]	Optimal value [m/s]	Improved [m/s]	Validated in FEM [m/s]	Error [%]
982.0	994.0	12.0	992.0	0.2

Discussion and Conclusion

Discussion

This study has presented an MECO method that optimizes the artillery’s maximum projectile axial speed and minimum projectile attitude disturbance at the muzzle while considering the two different stages of the projectile-barrel interaction.

Figure 4 shows the system-level convergence curves, which are implemented by different methods. The following insights are provided: First, the value of the system-level computing through the ECO strategy exhibits an oscillation characteristic, making it difficult to converge. Second, the solution converges when the distance criterion is adopted, but remains unstable in the first 8 iterations. Finally, the system-level convergence is improved when the boundary penalty of the subspace is further implemented, through only 13 iterations of MECO. It can be found that the boundaries of the design variables are narrowed by the MECO method, as shown in Figure 5. Compared with the other two methods, the MECO demonstrates higher efficiency and stability.

One interesting phenomenon observed is that the constraints of the maximum chamber pressure are broken through and the global objective function is also abnormal at the 6th iteration. For example, $p_{m(i=6)}=365$ MPa >320 MPa in Figure 6, and $F_{0(i=6)}=1.3878>$ $F_{0(i=5)}=0.6950$ in Figure 7. This phenomenon can be understood in the following way. Although the subspaces are subject to the constraints, the abnormal solutions that are incompatible with other subspaces are still produced since the objective function is non-smooth. This is contrary to the purpose of LMNCs of increasing each subspace’s “awareness” of their influence on other subspace. Furthermore, the system-level problem is an unconstrained optimization problem, and the system-level solution fluctuates drastically upon receiving the subspace’s abnormal solutions, as shown in Figure 5. The violation phenomenon happens when the system-level solution is projected into the constraint function and the objective function. At this time, the system-level design variables are deprived of the dynamic “move limit” function and the “guide direction” function, which ultimately affects the convergence.

It can be seen from Figure 5 that the system-level and each subspace’s shared variables all converge under the guidance of the dynamic “move limit” function. In this study, the distance criterion has been used in the subspace loop on the one hand, and the penalty design variables boundary method has been implemented on the other hand. The boundaries of the groove width $b$ , the rifling depth $t_{sh}$ , and the propellant thickness $e_1$ are narrowed as these variables change drastically at the 6th iteration ( $\mathit{\Delta }{d}_j$ > 1.5), as can be observed in Figure 5. At the 7th iteration the boundaries of these design variables are penalized by the MECO method. The system-level objective function then decreases rapidly and finally converges. The results show that the MECO method proposed in this study can solve the non-smooth optimization problem of the artillery launching process.

It can observe from Figure 7 that the curves of two indicators and the global objective function change drastically at the 6th iteration, due to the default of the maximum chamber pressure at the 6th iteration, as shown in Figure 6. The value of the global objective function has been improved from 0.7241 to 0.7140. The projectile axial speed $v_{xm}$ at the muzzle and the projectile disturbance coefficient ${\eta }_2$ concerned are not so convergent as the system-level variables, but these indicators have improved the result. They fluctuate slightly around the optimal value ( $v_{xm}$ = 994.0 m/s and ${\eta }_2$ = 0.9810) and coordinate the shared design variables of the propellant, the projectile, and the barrel used in the two-sequential process. As Roth and Kroo^30,31 reported, the subspace with the ECO strategy is objectively focused not only on satisfying compatibility but also on reducing the global objective function.

As shown in Figure 6 and Figure 7, the projectile axial speed $v_{xm}$ at the muzzle and the launching accuracy have been improved under the three constraints of the engraving resistance $f_R$ , the maximum chamber pressure $p_m$ , and the launching life. For example, $v_{xm}$ has increased to 994.0 m/s by 12 m/s, and ${\eta }_2$ has decreased from 1.0089 to 0.9810. At the same time, the projectile axial speed error between the optimal solution and the verified value in the FEM after optimization is less than 0.2%. The MECO method proves accurate and effective in solving the multi-stage complete optimization problem of artillery launch.

Conclusion

In this study, the dynamic model for artillery launch coupled with the propellant combustion has been developed to predict the performance of the engraving process of the rotating band and the movement process of the projectile in the barrel. The modified method has been developed with distance criterion and penalty design boundary in the MECO. To solve the optimization problem of the dynamic model for artillery launch, the model has been formulated in the form of MECO. The optimal results show that the system-level objective function presents oscillation characteristic that can hardly converge when the ECO strategy is implemented. The number of iterations of the MECO strategy is 13, less than that of the MECO strategy, showing higher efficiency and stability. The projectile axial speed at the muzzle has increased by 1.22% and the projectile disturbance coefficient has decreased by 2.77% under the constraints of the engraving resistance, the maximum chamber pressure, and the launching life. The projectile axial speed at the muzzle and the launching accuracy have been improved while the two-stage optimization problem has been solved effectively by the MECO strategy. Unlike most of the studies on artillery launch found in the literatures, the MECO is focused not only on satisfying the compatibility between the coupled parameters in the two sequential stages but also on reducing the global objective function, as in this paper from 0.7241 to 0.7140. Future work will involve the following aspects: (1) to generalize the deterministic optimization to uncertain problems based on the MECO strategy, (2) to generalize the optimization of artillery launch in bore to the full trajectory, and (3) to consider the dynamic weight problem of multi-objective optimization, and to couple it with MDO.