# Efficient Method of Firing Angle Calculation for Multiple Launch Rocket System Based on Polynomial Response Surface and Kriging Metamodels

Authors: Qiang Zhao; Qizhong Tang; Junli Han; Ming Yang; Zhihua Chen
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Aiming at solving the problem of firing angle calculation for the multiple launch rocket system (MLRS) under both standard and actual atmospheric conditions, an efficient method based on large sample data and metamodel is proposed. The polynomial response surface, Kriging, and the ensemble of metamodels are used to establish the functional relations between the firing angle, the maximum range angle, the maximum range, and various influencing factors under standard atmospheric conditions, and related processes are described in detail. On this basis, the initial values for the first two iterations are determined with the meteorological data being made full use of in the six degrees of freedom trajectory simulation, and then the firing angle corresponding to a specific range is automatically and iteratively calculated. The efficient method of firing angle calculation for the typical MLRS has been extensively tested with three cases. The results show that the high-order polynomial response surface, the Kriging predictors with Cubic, Gauss, and Spline correlation functions, and the ensemble of above four individual metamodels have better performances for predicting the firing angle under standard atmospheric conditions compared with those of other metamodels under identical conditions, and execution times of the above four individual metamodels with a training sample size of 9000 are all less than 0.9ms, which verifies the effectiveness and feasibility of the proposed method for calculating the firing angle under standard atmospheric conditions. Moreover, the number of iterations is effectively reduced by using the proposed iterative search approach under actual atmospheric conditions. This research can provide guidance for designing the fire control and command control system of the MLRS.


---

## Body

## 1. Introduction

As the primary suppressed weapon of the army, the MLRS plays an increasingly important role in the long-range precision strike. In the modern and future wars, it is crucial to improve the survivability of the MLRS, and reducing the launch preparation time to improve the rapid reaction capability has been becoming one of the most interesting research areas. Therefore, it is especially important to develop a rapid method for calculating the aiming parameters of the MLRS.Two primary aiming parameters including the firing angle and azimuth are determined before the fire mission [1]. The north-based azimuth can be calculated according to the given geographical positions of the launcher and target (latitude and longitude), and several existing methods available can be applied to calculate the accurate azimuth [2], such as the Bessel’s method, Rainsford’s method, Vincenty’s method, and Karney’s method. It is desirable to choose a method for calculating the azimuth according to the advantages and limitations of various methods. In addition, the whole process of calculating the azimuth is relatively simple and less time-consuming.The process of calculating the firing angle, however, is relatively complicated. There are two commonly used approaches for obtaining the firing angle: the traditional firing table approach and the iterative search approach via a trajectory program, of which the former can be realized by interpolating or approximating the tabular data under standard conditions then applying correction factors [3, 4]. Currently, the linear and quadratic interpolations are still two commonly used interpolation approaches for calculating the firing angle [5], and one of commonly used approximation approaches is the polynomial response surface (PRS) [3], particularly the polynomials of orders 1 through 4. With the exception of PRS, the feed forward neural network has also been successfully applied to predicting the firing angle of an artillery rocket [6]. However, other approximation techniques such as Kriging [7–9], radial basis function (RBF) [10], extreme learning machine (ELM) [11], and support vector machine (SVM) [12] are rarely used to predict the firing angle. The calculation result by using the firing table approach can meet the accuracy requirement within a certain range. However, when the actual atmospheric conditions are quite different from the standard atmospheric conditions, the firing angle obtained by using the firing table approach is rather inaccurate [13], which will lead to the firing mission failure of the free and simple controlled rockets with limited maneuverability.The iterative search approach is another method for calculating the firing angle by employing the trajectory program and iterative algorithm, which is applicable to the firing angle calculation under both standard and actual atmospheric conditions. The calculation efficiency depends strongly on the trajectory model, the iterative initial value, and the iterative algorithm. As one of the most well-known fire control software, the NATO Armaments Ballistic Kernel (NABK) [14, 15] searches for a firing angle corresponding to a specified range by utilizing the modified point mass or the five degrees of freedom (DOF) trajectory program [16] and the iterative algorithm. Despite the significant improvements that have been achieved, the details of the algorithms are mostly confidential and limited publications are currently available. By employing the six DOF trajectory program and an improved iterative search approach, L. J. Zhou et al. proposed a method for calculating the firing angle of the MLRS, which can reduce the number of iterations substantially [17]. D. H. Zhao et al. presented a method to determine the firing angle of the artillery based on the three-DOF trajectory program and binary search method [18]. P. Chusilp et al. came up with four iterative search algorithms and compared the calculation efficiency of four algorithms with the firing angle calculation of M107 projectile as a case study [19]. In addition, W. Charubhun et al. put forward an efficient method of firing angle calculation for the MLRS with the six DOF trajectory program and iterative binary search method by having prior knowledge of the ranges versus the firing angles [20].Without considering various types of systematic and random errors, the firing angle is determined by six influencing factors, including the latitude of launcherB0, the elevation of launcher H0, the target azimuth AT, the range between the launcher and target RG, the elevation of target HT, and the propellant temperature Ts. Moreover, the functional relation between the firing angle and six influence factors is highly nonlinear. Taking various influencing factors into consideration, conducting the research on a rapid and high-accuracy method for calculating the firing angle under standard atmospheric conditions is of great significance in the state of the art. On the one hand, the calculation result can be directly taken as the firing angle of the MLRS armed with guided rockets, which are capable of handling various types of systematic and random errors, such as aiming errors, manufacturing errors, aerodynamic errors, and wind induced errors. On the other hand, the calculation result can be used as the initial value of the iterative search approach.An efficient method for calculating the firing angle of the MLRS based on large sample data and metamodel is proposed in this paper, which can solve the problem of calculating the firing angle under both standard and actual atmospheric conditions. The PRS, Kriging, and the ensemble of metamodels (EM) [21–25] are used to establish the functional relations between the firing angle, the maximum range angle, the maximum range, and various influencing factors. Here, the metamodel is an engineering method that fits the collection of input-output pairs from runs of a simulation model, and it can be fully exploited rather than the more expensive and time-consuming simulation code [26, 27]. It turns out that the metamodel has been widely applied to the engineering design optimization problems [28–31].The outline of this paper is as follows. In the next section, the basic theories are introduced in detail including metamodels such as the PRS, Kriging and EM, the teaching-learning-based optimization algorithm [32, 33], and evaluation methods of the prediction accuracy of metamodels. Thereafter, three different kinds of problems of the firing angle calculation are defined, which includes the problem of calculating the firing angle under standard atmospheric conditions, the problem of calculating the maximum range angle under standard atmospheric conditions, and the problem of calculating the firing angle under actual atmospheric conditions. Here, the maximum range angle is the firing angle corresponding to the maximum range under specific constraint conditions. Subsequently, the results of applying the proposed method of firing angle calculation to three different kinds of problems are given. Finally, conclusions are provided in the end.


## 2. Basic Theories

### 2.1. Metamodels

#### 2.1.1. PRS

The PRS is a metamodel that conducts the function fitting via the statistical regression analysis, which has the advantages of low computational complexity and good robustness, and it is widely applied to the engineering design optimization problems. For the sake of simplicity and convenience, assuming that the response is one-dimensional, the expression of PRS can be written as(1)y^x=fxTβ=β0+∑i=1nβixi+∑i=1n∑j≥inβijxixj+···where y^x is the predicted value at the point x, xi is theith component of the vector x of lengthn, and β denotes the vector of regression parameters.


#### 2.1.2. Kriging

Kriging is a spatial interpolation method originating from the field of geostatistics [34], which is currently used to approximate deterministic data obtained by measurements or computer simulations. The mathematical formula of Kriging is given as follows:(2)y^x=fxTβ+zxwhere fxTβ is a trend function providing global approximation and z(x) is the realization of a stochastic process with mean 0 and variance σ2. The spatial covariance for design points w and x is defined as (3)covzw,zx=σ2Rθ,w,xwhere Rθ,w,x is a correlation function with the parameter vector θ, expressed as(4)Rθ,w,x=∏j=1nRjθj,wj-xj=∏j=1nRjθj,djwheredj=wj−xj and θj>0. Numerous correlation functions are available in the literature [7] including Cubic function, Exp function, Gauss function, Lin function, Spherical function, and Spline function, as shown in Table 1. Besides the most commonly used correlation functions shown in Table 1, the Matern class of correlation functions might prove more effective in practice [35], and it is possible that two cases of great interest are υ=3/2 (Matern32) and υ=5/2 (Matern52), for which (5)Rθ,w,xυ=3/2=1+3le-3l(6)Rθ,w,xυ=5/2=1+5l+5l23e-5lwhere l=∑j=1nθjdj2.Table 1
Common correlation functions.
Correlation function R j θ j , d j Cubic 1 - 3 ξ j 2 + 2 ξ j 3 ,   ξ j = min ⁡ 1 , θ j d j Exp e - θ j d j Gauss e - θ j d j 2 Lin max ⁡ 0,1 - θ j d j Spherical 1 - 1.5 ξ j + 0.5 ξ j 3 ,   ξ j = min ⁡ 1 , θ j d j Spline 1 - 15 θ j 2 d j 2 + 30 θ j 3 d j 3 0 ≤ θ j d j ≤ 0.2 1.25 1 - θ j d j 3 0.2 < θ j d j < 1 0 θ j d j ≥ 1


#### 2.1.3. EM

The EM, also known as the weighted average metamodel, has attracted much attention due to its good approximation capability [21], and it is generally expressed as (7)y^ex=∑i=1Mωiy^ixwhere y^ex is the predicted response by the EM,M denotes the number of individual metamodels, y^ix is the predicted response by theith individual metamodel, and ωi denotes the weight corresponding to theith individual metamodel. Obviously, the individual metamodel with higher prediction accuracy has larger weight and vice versa, and the sum of the weights is equal to 1; namely, ∑i=1Mωi=1.In general, the weights calculation is a key to the establishment of EM. There are currently a variety of weights calculation methods available in the literature, e.g., the heuristic calculation method, weights calculation method based on prediction variance, and combining metamodels method by minimizing the error [21–23]. The weights calculation method based on error minimization is employed in this paper, which takes the weights calculation process as an optimization problem. The design variables and optimization objective are the weights and the minimum prediction error of the EM, respectively, and the optimization problem can be formulated as(8)Find:ωiMin:Ee=Erry^eωi,y^ix,yxS.t.:∑i=1Mωi=1where Err denotes the chosen error metric for measuring the accuracy of the response predicted by the EM and yx is the actual response value at the point x. The root mean squared error is selected as the error metric Err. During the training process of the EM, there are two commonly used validation techniques, including the cross-validation technique and the validation technique with independent samples, of which the latter has lower computation cost compared with the former for large training samples. Consequently, the validation technique with independent samples is preferred, and the large testing samples are directly used as the independent validation samples in this paper. Based on the established individual metamodels, the weights of the EM are optimized by using the teaching-learning based algorithm [32, 33].


### 2.2. Teaching-Learning-Based Optimization Algorithm

Teaching-learning-based optimization (TLBO) algorithm is a population-based optimization method that is inspired by and based on the effect of the influence of a teacher on learners, and it is proposed by Rao et al. in 2011 [32]. In the basic TLBO algorithm, a group of learners is regarded as the current population, of which the best individual represents the teacher and the other individuals are regarded as the students. The process of basic TLBO algorithm consists of two parts, namely, ‘teacher phase’ and ‘learner phase’. The ‘teacher phase’ and ‘learner phase’ mean learning from the teacher and learning by the interaction between learners, respectively. The brief description of the basic TLBO algorithm is given as follows.


#### 2.2.1. Teacher Phase

During the teacher phase, the learners gain and update their knowledge according to the following expressions:(9)X′j,k,i=Xj,k,i+rand10,1∗Xj,kbest,i-TFMj,i(10)TF=round1+rand20,12-1where X′j,k,i and Xj,k,i denote the populations after and before learning from the teacher respectively; rand10,1 and rand20,1 are the random numbers between 0 and 1; Xj,kbest,i denotes the current best individual; TF denotes the teaching factor with the value being either 1 or 2; Mj,i is the mean individual of the population before learning from the teacher;round[X] denotes each element ofX is rounded to the nearest integer.


#### 2.2.2. Learner Phase

Having finished the teacher phase, the learners enhance their knowledge just via the interaction among themselves. A learner interacts randomly with other learners for increasing his or her knowledge. In the case of minimization problems, the updated learnerX′′j,P,i can be defined as(11)X′′j,P,i=X′j,P,i+rand30,1∗X′j,P,i-X′j,Q,ifX′j,P,i<fX′j,Q,iX′j,P,i+rand30,1∗X′j,Q,i-X′j,P,iotherwisewhere X′j,P,i and X′j,Q,i are two randomly selected learners; fX′j,P,i and fX′j,Q,i are the fitness values of the learners X′j,P,i and X′j,Q,i, respectively; rand30,1 denotes the random number between 0 and 1.


### 2.3. Evaluation Methods of the Prediction Accuracy of Metamodels

Three different metrics are applied to evaluate the prediction accuracies of various metamodels, including the maximum absolute error (MAE), maximum relative error (MRE), and root mean squared error (RMSE), expressed as, respectively,(12)MAE=maxyx-y^x(13)MRE=maxyx-y^xyx×100%(14)RMSE=1Ntest∑i=1Ntestyx-y^x2where Ntest is the number of testing samples.


## 2.1. Metamodels

### 2.1.1. PRS

The PRS is a metamodel that conducts the function fitting via the statistical regression analysis, which has the advantages of low computational complexity and good robustness, and it is widely applied to the engineering design optimization problems. For the sake of simplicity and convenience, assuming that the response is one-dimensional, the expression of PRS can be written as(1)y^x=fxTβ=β0+∑i=1nβixi+∑i=1n∑j≥inβijxixj+···where y^x is the predicted value at the point x, xi is theith component of the vector x of lengthn, and β denotes the vector of regression parameters.


### 2.1.2. Kriging

Kriging is a spatial interpolation method originating from the field of geostatistics [34], which is currently used to approximate deterministic data obtained by measurements or computer simulations. The mathematical formula of Kriging is given as follows:(2)y^x=fxTβ+zxwhere fxTβ is a trend function providing global approximation and z(x) is the realization of a stochastic process with mean 0 and variance σ2. The spatial covariance for design points w and x is defined as (3)covzw,zx=σ2Rθ,w,xwhere Rθ,w,x is a correlation function with the parameter vector θ, expressed as(4)Rθ,w,x=∏j=1nRjθj,wj-xj=∏j=1nRjθj,djwheredj=wj−xj and θj>0. Numerous correlation functions are available in the literature [7] including Cubic function, Exp function, Gauss function, Lin function, Spherical function, and Spline function, as shown in Table 1. Besides the most commonly used correlation functions shown in Table 1, the Matern class of correlation functions might prove more effective in practice [35], and it is possible that two cases of great interest are υ=3/2 (Matern32) and υ=5/2 (Matern52), for which (5)Rθ,w,xυ=3/2=1+3le-3l(6)Rθ,w,xυ=5/2=1+5l+5l23e-5lwhere l=∑j=1nθjdj2.Table 1
Common correlation functions.
Correlation function R j θ j , d j Cubic 1 - 3 ξ j 2 + 2 ξ j 3 ,   ξ j = min ⁡ 1 , θ j d j Exp e - θ j d j Gauss e - θ j d j 2 Lin max ⁡ 0,1 - θ j d j Spherical 1 - 1.5 ξ j + 0.5 ξ j 3 ,   ξ j = min ⁡ 1 , θ j d j Spline 1 - 15 θ j 2 d j 2 + 30 θ j 3 d j 3 0 ≤ θ j d j ≤ 0.2 1.25 1 - θ j d j 3 0.2 < θ j d j < 1 0 θ j d j ≥ 1


### 2.1.3. EM

The EM, also known as the weighted average metamodel, has attracted much attention due to its good approximation capability [21], and it is generally expressed as (7)y^ex=∑i=1Mωiy^ixwhere y^ex is the predicted response by the EM,M denotes the number of individual metamodels, y^ix is the predicted response by theith individual metamodel, and ωi denotes the weight corresponding to theith individual metamodel. Obviously, the individual metamodel with higher prediction accuracy has larger weight and vice versa, and the sum of the weights is equal to 1; namely, ∑i=1Mωi=1.In general, the weights calculation is a key to the establishment of EM. There are currently a variety of weights calculation methods available in the literature, e.g., the heuristic calculation method, weights calculation method based on prediction variance, and combining metamodels method by minimizing the error [21–23]. The weights calculation method based on error minimization is employed in this paper, which takes the weights calculation process as an optimization problem. The design variables and optimization objective are the weights and the minimum prediction error of the EM, respectively, and the optimization problem can be formulated as(8)Find:ωiMin:Ee=Erry^eωi,y^ix,yxS.t.:∑i=1Mωi=1where Err denotes the chosen error metric for measuring the accuracy of the response predicted by the EM and yx is the actual response value at the point x. The root mean squared error is selected as the error metric Err. During the training process of the EM, there are two commonly used validation techniques, including the cross-validation technique and the validation technique with independent samples, of which the latter has lower computation cost compared with the former for large training samples. Consequently, the validation technique with independent samples is preferred, and the large testing samples are directly used as the independent validation samples in this paper. Based on the established individual metamodels, the weights of the EM are optimized by using the teaching-learning based algorithm [32, 33].


## 2.1.1. PRS

The PRS is a metamodel that conducts the function fitting via the statistical regression analysis, which has the advantages of low computational complexity and good robustness, and it is widely applied to the engineering design optimization problems. For the sake of simplicity and convenience, assuming that the response is one-dimensional, the expression of PRS can be written as(1)y^x=fxTβ=β0+∑i=1nβixi+∑i=1n∑j≥inβijxixj+···where y^x is the predicted value at the point x, xi is theith component of the vector x of lengthn, and β denotes the vector of regression parameters.


## 2.1.2. Kriging

Kriging is a spatial interpolation method originating from the field of geostatistics [34], which is currently used to approximate deterministic data obtained by measurements or computer simulations. The mathematical formula of Kriging is given as follows:(2)y^x=fxTβ+zxwhere fxTβ is a trend function providing global approximation and z(x) is the realization of a stochastic process with mean 0 and variance σ2. The spatial covariance for design points w and x is defined as (3)covzw,zx=σ2Rθ,w,xwhere Rθ,w,x is a correlation function with the parameter vector θ, expressed as(4)Rθ,w,x=∏j=1nRjθj,wj-xj=∏j=1nRjθj,djwheredj=wj−xj and θj>0. Numerous correlation functions are available in the literature [7] including Cubic function, Exp function, Gauss function, Lin function, Spherical function, and Spline function, as shown in Table 1. Besides the most commonly used correlation functions shown in Table 1, the Matern class of correlation functions might prove more effective in practice [35], and it is possible that two cases of great interest are υ=3/2 (Matern32) and υ=5/2 (Matern52), for which (5)Rθ,w,xυ=3/2=1+3le-3l(6)Rθ,w,xυ=5/2=1+5l+5l23e-5lwhere l=∑j=1nθjdj2.Table 1
Common correlation functions.
Correlation function R j θ j , d j Cubic 1 - 3 ξ j 2 + 2 ξ j 3 ,   ξ j = min ⁡ 1 , θ j d j Exp e - θ j d j Gauss e - θ j d j 2 Lin max ⁡ 0,1 - θ j d j Spherical 1 - 1.5 ξ j + 0.5 ξ j 3 ,   ξ j = min ⁡ 1 , θ j d j Spline 1 - 15 θ j 2 d j 2 + 30 θ j 3 d j 3 0 ≤ θ j d j ≤ 0.2 1.25 1 - θ j d j 3 0.2 < θ j d j < 1 0 θ j d j ≥ 1


## 2.1.3. EM

The EM, also known as the weighted average metamodel, has attracted much attention due to its good approximation capability [21], and it is generally expressed as (7)y^ex=∑i=1Mωiy^ixwhere y^ex is the predicted response by the EM,M denotes the number of individual metamodels, y^ix is the predicted response by theith individual metamodel, and ωi denotes the weight corresponding to theith individual metamodel. Obviously, the individual metamodel with higher prediction accuracy has larger weight and vice versa, and the sum of the weights is equal to 1; namely, ∑i=1Mωi=1.In general, the weights calculation is a key to the establishment of EM. There are currently a variety of weights calculation methods available in the literature, e.g., the heuristic calculation method, weights calculation method based on prediction variance, and combining metamodels method by minimizing the error [21–23]. The weights calculation method based on error minimization is employed in this paper, which takes the weights calculation process as an optimization problem. The design variables and optimization objective are the weights and the minimum prediction error of the EM, respectively, and the optimization problem can be formulated as(8)Find:ωiMin:Ee=Erry^eωi,y^ix,yxS.t.:∑i=1Mωi=1where Err denotes the chosen error metric for measuring the accuracy of the response predicted by the EM and yx is the actual response value at the point x. The root mean squared error is selected as the error metric Err. During the training process of the EM, there are two commonly used validation techniques, including the cross-validation technique and the validation technique with independent samples, of which the latter has lower computation cost compared with the former for large training samples. Consequently, the validation technique with independent samples is preferred, and the large testing samples are directly used as the independent validation samples in this paper. Based on the established individual metamodels, the weights of the EM are optimized by using the teaching-learning based algorithm [32, 33].


## 2.2. Teaching-Learning-Based Optimization Algorithm

Teaching-learning-based optimization (TLBO) algorithm is a population-based optimization method that is inspired by and based on the effect of the influence of a teacher on learners, and it is proposed by Rao et al. in 2011 [32]. In the basic TLBO algorithm, a group of learners is regarded as the current population, of which the best individual represents the teacher and the other individuals are regarded as the students. The process of basic TLBO algorithm consists of two parts, namely, ‘teacher phase’ and ‘learner phase’. The ‘teacher phase’ and ‘learner phase’ mean learning from the teacher and learning by the interaction between learners, respectively. The brief description of the basic TLBO algorithm is given as follows.


### 2.2.1. Teacher Phase

During the teacher phase, the learners gain and update their knowledge according to the following expressions:(9)X′j,k,i=Xj,k,i+rand10,1∗Xj,kbest,i-TFMj,i(10)TF=round1+rand20,12-1where X′j,k,i and Xj,k,i denote the populations after and before learning from the teacher respectively; rand10,1 and rand20,1 are the random numbers between 0 and 1; Xj,kbest,i denotes the current best individual; TF denotes the teaching factor with the value being either 1 or 2; Mj,i is the mean individual of the population before learning from the teacher;round[X] denotes each element ofX is rounded to the nearest integer.


### 2.2.2. Learner Phase

Having finished the teacher phase, the learners enhance their knowledge just via the interaction among themselves. A learner interacts randomly with other learners for increasing his or her knowledge. In the case of minimization problems, the updated learnerX′′j,P,i can be defined as(11)X′′j,P,i=X′j,P,i+rand30,1∗X′j,P,i-X′j,Q,ifX′j,P,i<fX′j,Q,iX′j,P,i+rand30,1∗X′j,Q,i-X′j,P,iotherwisewhere X′j,P,i and X′j,Q,i are two randomly selected learners; fX′j,P,i and fX′j,Q,i are the fitness values of the learners X′j,P,i and X′j,Q,i, respectively; rand30,1 denotes the random number between 0 and 1.


## 2.2.1. Teacher Phase

During the teacher phase, the learners gain and update their knowledge according to the following expressions:(9)X′j,k,i=Xj,k,i+rand10,1∗Xj,kbest,i-TFMj,i(10)TF=round1+rand20,12-1where X′j,k,i and Xj,k,i denote the populations after and before learning from the teacher respectively; rand10,1 and rand20,1 are the random numbers between 0 and 1; Xj,kbest,i denotes the current best individual; TF denotes the teaching factor with the value being either 1 or 2; Mj,i is the mean individual of the population before learning from the teacher;round[X] denotes each element ofX is rounded to the nearest integer.


## 2.2.2. Learner Phase

Having finished the teacher phase, the learners enhance their knowledge just via the interaction among themselves. A learner interacts randomly with other learners for increasing his or her knowledge. In the case of minimization problems, the updated learnerX′′j,P,i can be defined as(11)X′′j,P,i=X′j,P,i+rand30,1∗X′j,P,i-X′j,Q,ifX′j,P,i<fX′j,Q,iX′j,P,i+rand30,1∗X′j,Q,i-X′j,P,iotherwisewhere X′j,P,i and X′j,Q,i are two randomly selected learners; fX′j,P,i and fX′j,Q,i are the fitness values of the learners X′j,P,i and X′j,Q,i, respectively; rand30,1 denotes the random number between 0 and 1.


## 2.3. Evaluation Methods of the Prediction Accuracy of Metamodels

Three different metrics are applied to evaluate the prediction accuracies of various metamodels, including the maximum absolute error (MAE), maximum relative error (MRE), and root mean squared error (RMSE), expressed as, respectively,(12)MAE=maxyx-y^x(13)MRE=maxyx-y^xyx×100%(14)RMSE=1Ntest∑i=1Ntestyx-y^x2where Ntest is the number of testing samples.


## 3. An Efficient Method for Calculating the Firing Angle of the MLRS Based on Metamodels

### 3.1. Firing Angle Calculation under Standard Atmospheric Conditions

Under standard atmospheric conditions, the actual process of firing angle calculation for the MLRS based on the metamodel is shown in Figure1, and the following two specific steps are involved.Figure 1
Calculation procedure of the firing angle under standard atmospheric conditions.Step 1 (build the trajectory model).
The accuracy of the firing angle calculation is directly determined by the trajectory model. We apply the six DOF trajectory model, and the effects of the oblateness and rotation of the Earth are taken into account. Moreover, the fudge factors are used. Consequently, the calculation results by utilizing the six DOF trajectory model agree well with flight test results.Step 2.
Establish the functional relation between the firing angle andB0, H0, AT, RG, HT, Ts via metamodels, and the whole process mainly includes obtaining a dataset through Latin hypercube sampling (LHS) [36] and computer simulations, selecting the metamodel type, building the metamodel, and validating the metamodel. Here, the LHS is especially suitable for very large designs.
Once the geographical positions of launcher and target are given, the detailed process of the firing data calculation under standard atmospheric conditions consists of the following steps,Step 1.
CalculateAT and RG according to the given latitude, longitude, and elevation of the launcher and target.Step 2.
Calculate the firing angle by using the established functional relation between the firing angle andB0, H0, AT, RG, HT, Ts.
With regard to the firing angles corresponding to various design points under standard atmospheric conditions, they are obtained via the iterative search approach and six DOF trajectory program, and the following six specific steps are involved.Step 1.
Calculate the firing angleθ1 corresponding to the given RG and H0 by using the Lagrange method [37] and take θ1 as the initial value of the iteration.Step 2.
Calculate the rangeR1 corresponding to the θ1 via the six DOF trajectory program, as well as the difference between R1 and RG by d1=RG-R1.Step 3.
Estimate the step lengthΔθ1 according to d1 and calculate the firing angle of the second iteration by θ2=θ1+Δθ1.Step 4.
Calculate the rangeR2 corresponding to the θ2 via the six DOF trajectory program.Step 5.
Determine the firing angle with higher accuracy according to the following iteration formula:(15)θk=θk-1+θk-1-θk-2Rk-1-Rk-2RG-Rk-1where θk, θk-1, and θk-2 are the firing angles at thekth, (k-1)th and (k-2)th iterations, respectively; Rk-1 and Rk-2 are the ranges corresponding to firing angles θk-1 and θk-2, respectively. Based on formula (15), the iterative calculation of the firing angle is conducted with the six-DOF trajectory program, until the following convergence criterion is met:(16)Rk-RG≤4RG10000where Rk is the range corresponding to θk. Thereafter, RG is replaced by Rk in the sample.Step 6.
Iterate from Step 1 to Step 5, until the calculation processes corresponding to all design points have been finished.There is a special case thatRG is greater than the maximum range corresponding to the given B0, H0, AT, HT, and Ts. In this case, it only needs to iterate the firing angle to the maximum range angle by the iterative method in Section 3.2, and then replace RG by the maximum range.


### 3.2. Maximum Range Angle Calculation under Standard Atmospheric Conditions

Without considering various types of systematic and random errors, the maximum range angle and maximum range under standard atmospheric conditions are determined byB0, H0, AT, HT and Ts. Prior to establishing the functional relations for the maximum range angle and maximum range, the training samples and testing samples are obtained, which consists of the following two specific steps.Step 1.
Calculate the maximum range angle with variable step sizes, which consists of four rounds of iterations with step sizes of 1.000, 0.100, 0.010, and 0.001, respectively. In theith round, let θi1, ui, and θi2+ be the initial value of the iteration, the step size, and the assumed firing angle of the second iteration, respectively, where θi2+=θi1+ui, and the ranges corresponding to θi2+ and θi1 are Xi2+ and Xi1, respectively. Then thekth firing angle θik can be expressed as(17)θik=θi1+kuiXi2+>Xi1θi1-kuiotherwisewherek denotes the number of iterations. For each round, iterate formula (17), until the maximum range angle is reached. For the subsequent round, θi1 is initialized by the final output of the (i-1)th round. After four rounds of iterative calculations, output the corresponding maximum range.Step 2.
Iterate Step 1, until the calculation processes corresponding to all design points have been finished.After obtaining the training samples and testing samples, the functional relations between the maximum range angle, the maximum range, and various influencing factors are established by using different metamodels.


### 3.3. Firing Angle Calculation under Actual Atmospheric Conditions

Compared with various types of systematic and random errors, the meteorological data has a greater influence on the accuracy of the firing angle of the MLRS armed with free and simple controlled rockets. Consequently, it is essential to make full use of the meteorological data available, including the actual atmospheric pressure, virtual temperature, wind velocity and wind direction.There are a number of iterative formulas available to calculate the accurate firing angle under actual atmospheric conditions [19]. The formulas have been extensively tested for the 155mm M107 projectile as a case study, by comparison with results from other formulas. It has been proven that formula (15) is preferred to other iterative formulas, and it is just the iterative formula (15) that we use in this paper. Once the iterative formula is determined, the initial values for the first two iterations are of great importance to the minimization of the number of iterations, of which the firing angle θ1 for the first iteration is obtained by using the established functional relation between the firing angle and various influencing factors.LetR1 be the range corresponding to θ1 under actual atmospheric conditions, and θ1∗ is the firing angle corresponding to R1 under standard atmospheric conditions. Let θ2 be the firing angle corresponding to the range R2 under actual atmospheric conditions, and R2∗ is the range corresponding to θ2 under standard atmospheric conditions. Two basic assumptions are given as follows.Assumption 1.
The difference betweenR1 and RG is equal to the difference between RG and R2∗; namely, R1-RG=RG-R2∗.Assumption 2.
The difference betweenθ1 and θ1∗ is equal to the difference between θ2 and θ1; namely, θ1-θ1∗=θ2-θ1.However, there is an unavoidable case thatRG is near the minimum or maximum ranges, which may result in R1 and R2∗ being less than the minimum range RMin or greater than the maximum range RMax under standard atmospheric conditions. In this case, if metamodels are still used to predict θ2, it may lead to the sharp decrease of the prediction accuracy. Aiming at solving this problem, an efficient method for calculating θ2 according to the ranges R1 and R2∗ is proposed, as shown below.If there isR1≤RMin or R1≥RMax such that RMin≤R2∗≤RMax, then follow Assumption 1 to calculate R2∗ by R2∗=2RG-R1 as well as the corresponding firing angle θ2∗ under standard atmospheric conditions, and subsequently θ2 is obtained by θ=θ2∗. Similarly, if there exists R2∗≤RMin or R2∗≥RMax such that RMin≤R1≤RMax, then follow Assumption 2 to calculate the θ2 directly by θ2=2θ1-θ1∗. Regarding other cases where RMin≤R2∗≤RMax and RMin≤R1≤RMax, θ2 is then calculated by θ2=θ2∗+2θ1-θ1∗/2 based on Assumptions 1 and 2. To sum up, we get the following expression for θ2,(18)θ2=θ2∗R1<RMinorR1>RMax2θ1-θ1∗R2∗<RMinorR2∗>RMaxθ2∗+2θ1-θ1∗2otherwiseThe detailed steps for calculating the firing angle under actual atmospheric conditions are summarized as follows.Step 1.
Based on the established functional relation between the firing angle and various influencing factors, calculate the firing angleθ1 corresponding to the given B0, H0, AT, RG, HT, and Ts.Step 2.
By utilizing the six DOF trajectory program, calculate the rangeR1 corresponding to θ1 under actual atmospheric conditions.Step 3.
Based on the established functional relation between the maximum range and various influencing factors, calculate the maximum rangeRMax corresponding to the given B0, H0, AT, HT, and Ts.Step 4.
Check whetherR1 satisfies RMin≤R1≤RMax, and if yes, calculate the firing angle θ1∗ corresponding to R1 based on the established functional relation between the firing angle and various influencing factors.Step 5.
Calculate the distance difference betweenR1 and RG by ΔR1=R1-RG, and then calculate the range R2∗ by R2∗=RG-ΔR1.Step 6.
Check whetherR2∗ satisfies RMin≤R2∗≤RMax, and if yes, calculate the firing angle θ2∗ corresponding to the R2∗ based on the established functional relation between the firing angle and various influencing factors.Step 7.
Calculateθ2 by the following equation: (19)θ2=θ2∗R1<RMinorR1>RMax2θ1-θ1∗R2∗<RMinorR2∗>RMaxθ2∗+2θ1-θ1∗2otherwiseand then calculate the range R2 corresponding to θ2 under actual atmospheric conditions by utilizing the six DOF trajectory program.Step 8.
Having obtained the iterative initial valuesθ1, θ2 and their corresponding ranges R1, R2, based on the six DOF trajectory program, iterate formula (15) until the convergence criterion (16) is met.


## 3.1. Firing Angle Calculation under Standard Atmospheric Conditions

Under standard atmospheric conditions, the actual process of firing angle calculation for the MLRS based on the metamodel is shown in Figure1, and the following two specific steps are involved.Figure 1
Calculation procedure of the firing angle under standard atmospheric conditions.Step 1 (build the trajectory model).
The accuracy of the firing angle calculation is directly determined by the trajectory model. We apply the six DOF trajectory model, and the effects of the oblateness and rotation of the Earth are taken into account. Moreover, the fudge factors are used. Consequently, the calculation results by utilizing the six DOF trajectory model agree well with flight test results.Step 2.
Establish the functional relation between the firing angle andB0, H0, AT, RG, HT, Ts via metamodels, and the whole process mainly includes obtaining a dataset through Latin hypercube sampling (LHS) [36] and computer simulations, selecting the metamodel type, building the metamodel, and validating the metamodel. Here, the LHS is especially suitable for very large designs.
Once the geographical positions of launcher and target are given, the detailed process of the firing data calculation under standard atmospheric conditions consists of the following steps,Step 1.
CalculateAT and RG according to the given latitude, longitude, and elevation of the launcher and target.Step 2.
Calculate the firing angle by using the established functional relation between the firing angle andB0, H0, AT, RG, HT, Ts.
With regard to the firing angles corresponding to various design points under standard atmospheric conditions, they are obtained via the iterative search approach and six DOF trajectory program, and the following six specific steps are involved.Step 1.
Calculate the firing angleθ1 corresponding to the given RG and H0 by using the Lagrange method [37] and take θ1 as the initial value of the iteration.Step 2.
Calculate the rangeR1 corresponding to the θ1 via the six DOF trajectory program, as well as the difference between R1 and RG by d1=RG-R1.Step 3.
Estimate the step lengthΔθ1 according to d1 and calculate the firing angle of the second iteration by θ2=θ1+Δθ1.Step 4.
Calculate the rangeR2 corresponding to the θ2 via the six DOF trajectory program.Step 5.
Determine the firing angle with higher accuracy according to the following iteration formula:(15)θk=θk-1+θk-1-θk-2Rk-1-Rk-2RG-Rk-1where θk, θk-1, and θk-2 are the firing angles at thekth, (k-1)th and (k-2)th iterations, respectively; Rk-1 and Rk-2 are the ranges corresponding to firing angles θk-1 and θk-2, respectively. Based on formula (15), the iterative calculation of the firing angle is conducted with the six-DOF trajectory program, until the following convergence criterion is met:(16)Rk-RG≤4RG10000where Rk is the range corresponding to θk. Thereafter, RG is replaced by Rk in the sample.Step 6.
Iterate from Step 1 to Step 5, until the calculation processes corresponding to all design points have been finished.There is a special case thatRG is greater than the maximum range corresponding to the given B0, H0, AT, HT, and Ts. In this case, it only needs to iterate the firing angle to the maximum range angle by the iterative method in Section 3.2, and then replace RG by the maximum range.


## 3.2. Maximum Range Angle Calculation under Standard Atmospheric Conditions

Without considering various types of systematic and random errors, the maximum range angle and maximum range under standard atmospheric conditions are determined byB0, H0, AT, HT and Ts. Prior to establishing the functional relations for the maximum range angle and maximum range, the training samples and testing samples are obtained, which consists of the following two specific steps.Step 1.
Calculate the maximum range angle with variable step sizes, which consists of four rounds of iterations with step sizes of 1.000, 0.100, 0.010, and 0.001, respectively. In theith round, let θi1, ui, and θi2+ be the initial value of the iteration, the step size, and the assumed firing angle of the second iteration, respectively, where θi2+=θi1+ui, and the ranges corresponding to θi2+ and θi1 are Xi2+ and Xi1, respectively. Then thekth firing angle θik can be expressed as(17)θik=θi1+kuiXi2+>Xi1θi1-kuiotherwisewherek denotes the number of iterations. For each round, iterate formula (17), until the maximum range angle is reached. For the subsequent round, θi1 is initialized by the final output of the (i-1)th round. After four rounds of iterative calculations, output the corresponding maximum range.Step 2.
Iterate Step 1, until the calculation processes corresponding to all design points have been finished.After obtaining the training samples and testing samples, the functional relations between the maximum range angle, the maximum range, and various influencing factors are established by using different metamodels.


## 3.3. Firing Angle Calculation under Actual Atmospheric Conditions

Compared with various types of systematic and random errors, the meteorological data has a greater influence on the accuracy of the firing angle of the MLRS armed with free and simple controlled rockets. Consequently, it is essential to make full use of the meteorological data available, including the actual atmospheric pressure, virtual temperature, wind velocity and wind direction.There are a number of iterative formulas available to calculate the accurate firing angle under actual atmospheric conditions [19]. The formulas have been extensively tested for the 155mm M107 projectile as a case study, by comparison with results from other formulas. It has been proven that formula (15) is preferred to other iterative formulas, and it is just the iterative formula (15) that we use in this paper. Once the iterative formula is determined, the initial values for the first two iterations are of great importance to the minimization of the number of iterations, of which the firing angle θ1 for the first iteration is obtained by using the established functional relation between the firing angle and various influencing factors.LetR1 be the range corresponding to θ1 under actual atmospheric conditions, and θ1∗ is the firing angle corresponding to R1 under standard atmospheric conditions. Let θ2 be the firing angle corresponding to the range R2 under actual atmospheric conditions, and R2∗ is the range corresponding to θ2 under standard atmospheric conditions. Two basic assumptions are given as follows.Assumption 1.
The difference betweenR1 and RG is equal to the difference between RG and R2∗; namely, R1-RG=RG-R2∗.Assumption 2.
The difference betweenθ1 and θ1∗ is equal to the difference between θ2 and θ1; namely, θ1-θ1∗=θ2-θ1.However, there is an unavoidable case thatRG is near the minimum or maximum ranges, which may result in R1 and R2∗ being less than the minimum range RMin or greater than the maximum range RMax under standard atmospheric conditions. In this case, if metamodels are still used to predict θ2, it may lead to the sharp decrease of the prediction accuracy. Aiming at solving this problem, an efficient method for calculating θ2 according to the ranges R1 and R2∗ is proposed, as shown below.If there isR1≤RMin or R1≥RMax such that RMin≤R2∗≤RMax, then follow Assumption 1 to calculate R2∗ by R2∗=2RG-R1 as well as the corresponding firing angle θ2∗ under standard atmospheric conditions, and subsequently θ2 is obtained by θ=θ2∗. Similarly, if there exists R2∗≤RMin or R2∗≥RMax such that RMin≤R1≤RMax, then follow Assumption 2 to calculate the θ2 directly by θ2=2θ1-θ1∗. Regarding other cases where RMin≤R2∗≤RMax and RMin≤R1≤RMax, θ2 is then calculated by θ2=θ2∗+2θ1-θ1∗/2 based on Assumptions 1 and 2. To sum up, we get the following expression for θ2,(18)θ2=θ2∗R1<RMinorR1>RMax2θ1-θ1∗R2∗<RMinorR2∗>RMaxθ2∗+2θ1-θ1∗2otherwiseThe detailed steps for calculating the firing angle under actual atmospheric conditions are summarized as follows.Step 1.
Based on the established functional relation between the firing angle and various influencing factors, calculate the firing angleθ1 corresponding to the given B0, H0, AT, RG, HT, and Ts.Step 2.
By utilizing the six DOF trajectory program, calculate the rangeR1 corresponding to θ1 under actual atmospheric conditions.Step 3.
Based on the established functional relation between the maximum range and various influencing factors, calculate the maximum rangeRMax corresponding to the given B0, H0, AT, HT, and Ts.Step 4.
Check whetherR1 satisfies RMin≤R1≤RMax, and if yes, calculate the firing angle θ1∗ corresponding to R1 based on the established functional relation between the firing angle and various influencing factors.Step 5.
Calculate the distance difference betweenR1 and RG by ΔR1=R1-RG, and then calculate the range R2∗ by R2∗=RG-ΔR1.Step 6.
Check whetherR2∗ satisfies RMin≤R2∗≤RMax, and if yes, calculate the firing angle θ2∗ corresponding to the R2∗ based on the established functional relation between the firing angle and various influencing factors.Step 7.
Calculateθ2 by the following equation: (19)θ2=θ2∗R1<RMinorR1>RMax2θ1-θ1∗R2∗<RMinorR2∗>RMaxθ2∗+2θ1-θ1∗2otherwiseand then calculate the range R2 corresponding to θ2 under actual atmospheric conditions by utilizing the six DOF trajectory program.Step 8.
Having obtained the iterative initial valuesθ1, θ2 and their corresponding ranges R1, R2, based on the six DOF trajectory program, iterate formula (15) until the convergence criterion (16) is met.


## 4. Case Studies

The case studies considered here seek to reveal the effectiveness of the firing angle calculation method for the MLRS under both standard and actual atmospheric conditions. Table2 provides the influencing factors and their variation ranges, where the mil is a unit of angular measurement used in artillery equal to 1/6400.00 of a complete revolution.Table 2
Influencing factors and their variation ranges.
Factors Units Lower bound Upper bound B 0 ° 0.00 70.00 H 0 m 0.00 5000.00 A T mil 0.00 6400.00 R G km 80.00 300.00 H T m 0.00 5000.00 T s °C -45.00 55.00


### 4.1. Case 1:Firing Angle Calculation under Standard Atmospheric Conditions

The six DOF trajectory simulation is based on the International Standard Atmosphere (ISO 2533:1975). The design of experiments is carried out for six influencing factors by using the LHS, and the total numbers of training samples and testing samples are 9000 and 5000, respectively. The firing angles corresponding to various design points are obtained by employing the iterative search approach and six DOF trajectory program described in Section3.1.Table3 shows the prediction accuracies of the PRS with different orders. As for the orders of PRS ranging from 1 to 9, the prediction accuracy can be improved significantly with the increase of the order. It is easily seen that the prediction power is relatively low for the commonly used orders of PRS ranging from 1 to 4. Regarding the prediction of the firing angle, the MAE, MRE, and RMSE of the 4th-order PRS are 7.055mil, 0.701%, and 1.424mil, respectively, which indicates that the low-order PRS lacks the capability of approximating the highly complicated functional relation. Compared with the low-order PRS, the high-order PRS has higher prediction accuracy of the firing angle, and the MAE, MRE, and RMSE of the 9th PRS are 1.047mil, 0.101%, and 0.103mil, respectively. Furthermore, in order to further evaluate the prediction accuracies of the PRS with different orders, the range error induced by the prediction error of the firing angle is analyzed simultaneously, which has the same changing trend as that of the firing angle, and the MAE, MRE, and RMSE of the 9th PRS are 291.49m, 0.185%, and 51.08m, respectively. Thus, high-order PRS is preferable to the establishment of a highly nonlinear functional relation between the firing angle and B0, H0, AT, RG, HT, Ts.Table 3
Effect of the order of PRS on prediction accuracy of firing angle.
Order Firing angle Range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) 1 58.774 5.562 13.584 18751.03 12.476 7165.87 2 27.953 3.535 8.638 12062.80 9.439 4817.94 3 11.159 1.383 2.721 4410.95 3.394 1463.24 4 7.055 0.701 1.424 2893.34 1.362 752.29 5 4.789 0.463 0.673 1436.73 0.695 336.77 6 2.976 0.288 0.372 898.42 0.512 180.78 7 2.122 0.204 0.215 589.93 0.232 102.04 8 1.509 0.146 0.137 422.01 0.244 64.78 9 1.047 0.101 0.103 291.49 0.185 51.08 10 2.052 0.248 0.166 1308.71 0.501 94.99In terms of the prediction of firing angle, the correlation function significantly affects the prediction accuracy of Kriging with the trend function of 4th PRS, as illustrated in Table4. Obviously, the Kriging predictors with different correlation functions all have higher prediction accuracies by comparison with the prediction results of the low-order PRS under identical conditions. Nevertheless, the correlation function still plays a significant role in Kriging. The Gauss, Spline and Cubic correlation functions of great interest demonstrate a promising capability to predict the firing angle, which have a parabolic behavior near the origin. In particular, the Gauss and Spline correlation functions are quite similar in terms of the prediction accuracy. However, if Linear, Exp, and Spherical functions with a linear behavior near the origin are selected as the correlation functions, the prediction performances are relatively poor.Table 4
Effect of the correlation function of Kriging on prediction accuracy of firing angle.
Correlation function Firing angle Range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) Cubic 0.809 0.080 0.099 286.63 0.108 51.45 Exp 2.242 0.223 0.227 664.52 0.250 85.75 Gauss 0.865 0.083 0.084 243.33 0.094 39.71 Linear 2.426 0.233 0.233 674.98 0.262 93.47 Matern32 1.667 0.161 0.150 465.00 0.180 60.81 Matern52 1.192 0.115 0.121 332.74 0.152 53.66 Spherical 2.330 0.224 0.221 649.45 0.251 83.96 Spline 0.738 0.074 0.087 262.95 0.103 44.97The high-order PRS (PRSho, order=9) and Kriging predictors with Cubic (KRcub), Exp (KRexp), Gauss (KRgau), Linear (KRlin), Matern32 (KRmat32), Matern52 (KRmat52), Spherical (KRsph), and Spline (KRspl) correlation functions are used to construct the EM. The weights of the EM are calculated by using the TLBO algorithm and initialize the optimization parameters of the population sizeNP=30 and the number of generationsNG=600. The optimized weights corresponding to the PRSho, KRcub, KRexp, KRgau, KRlin, KRmat32, KRmat52, KRsph, and KRspl are 0.139072, 0.146959, 0.000076, 0.363127, 0.000073, 0.000255, 0.000681, 0.000082, and 0.349675, respectively, and the sum of the weights corresponding to the PRSho, KRcub, KRgau, and KRspl is 0.998833, which indicates that the above four metamodels play a key role in the construction of the EM. With regard to the prediction of firing angle, the MAE, MRE, and RMSE of the EM are 0.698 mil, 0.068%, and 0.071mil, respectively. In terms of the range error induced by the prediction error of firing angle, the MAE, MRE, and RMSE of the EM are 195.08m, 0.076%, and 35.92m, respectively. Obviously, the EM outperforms other individual metamodels. In light of the above analysis, we can conclude that the PRSho, KRcub, KRgau, KRspl, and EM have higher prediction accuracies.The distributions of absolute errors of the firing angle and range of five metamodels are shown in Figures2 and 3. According to Figure 2, for each of the five metamodels, the testing samples with the absolute error of firing angle being less than 0.15mil and 0.30mil account for over 91.16% and 98.24%, respectively. As illustrated in Figure 3, the proportions of testing samples with the absolute error of range being less than 50m and 100m makes up over 70.00% and 95.48%, respectively. The above results further verify the effectiveness of the firing angle calculation method under standard atmospheric conditions.Figure 2
Distributions of the absolute error of firing angle of five metamodels.
(a)
Overall distribution (b)
Partial amplification of P1 portionFigure 3
Distributions of the absolute error of range of five metamodels.
(a)
Overall distribution (b)
Partial amplification of P2 portionThe execution time of a program depends strongly on the programming environment and the computing platform used. Based on the established metamodels PRSho, KRcub, KRgau, KRspl, and EM for the prediction of firing angle, five programs are written in Microsoft Visual C++ 6.0. Each program operates on a PC with an Intel(R) Core(TM) i5-6300HQ CPU @ 2.30GHz, 4GB RAM, Windows 7 (64bit). The execution times of various metamodels for a calculation of the firing angle are given in Figure4, where the execution time of each metamodel is the average time of the execution times of 5000 testing cases. It is quite obvious that the execution times of four individual metamodels are all less than 0.9ms, and the PRSho and KRspl are the fastest and slowest metamodels with execution times being 0.03676ms and 0.82525ms, respectively. Regarding the execution time of the EM, it is 2.23317ms, which is still far less than the required time of the fire control system that is measured in seconds. Therefore, based on the above analysis of the rapidity of the firing angle prediction, the feasibility of the firing angle calculation method under standard atmospheric conditions is verified.Figure 4
Execution times of various metamodels for a prediction of firing angle.Figures5 and 6 clearly illustrate how the prediction performances of metamodels depend on the training sample size. The training samples with different sizes are generated independently based on the LHS. The PRSho predictors corresponding to the training sample size ranging from 3000 to 21000 with an increment of 3000 have the orders of 7, 8, 9, 9, 10, 10, and 10, respectively. The EMs of the firing angle corresponding to different training samples are constructed in the same way as the EM with 9000 training samples, and thus the optimized weights corresponding to different training samples are different. Regarding the prediction error of the firing angle and its inducing range error, not surprisingly, the RMSE of each metamodel shows a downward trend with the training sample size. In particular, if the number of the training samples is less than 9000, the RMSE of each metamodel decreases significantly with the increase of training sample size, which indicates that the prediction performance can be improved remarkably with an increasing training sample size, and the PRSho has the most significant decrease. However, while the training sample size is more than 9000, the RMSE of each metamodel decreases relatively slowly, which indicates that 9000 training samples are enough to achieve reasonably good prediction performance. The MAE and MRE of each metamodel have the similar changing trends, and the KRgau is then taken as an example. Regarding the prediction of firing angle, the MAE and MRE of the KRgau are 1.494mil and 0.150%, respectively, at a training sample size of 3000, but they decrease to 0.865mil and 0.083% with 9000 training samples and further decrease to 0.451mil and 0.043% with 21000 training samples.Figure 5
Variation of the prediction accuracy of firing angle with the training sample size.Figure 6
Variation of the range error due to the prediction error of firing angle with the training sample size.


### 4.2. Case 2: Calculations of Maximum Range Angle and Maximum Range under Standard Atmospheric Conditions

The research on calculation methods of the maximum range angle and maximum range under standard atmospheric conditions is carried out by using the PRS, Kriging, and EM. The LHS is selected as the method for the design of experiments, and the variation ranges of influencing factorsB0, H0, AT, HT, and Ts are shown in Table 2. Then the training samples and testing samples are obtained with the sizes of 500 and 1000, respectively. Thereafter, the functional relations between maximum range angle, maximum range, and B0, H0, AT, HT, Ts are established.The effect of the order of PRS on prediction accuracy of maximum range angle is illustrated in Table5. Obviously, for the order within the range of 1 to 5, the prediction performance can be improved substantially with an increasing order, and the 5th PRS has the best prediction performance. In terms of the prediction of maximum range angle, the MAE, MRE, and RMSE of 5th PRS are 0.949mil, 0.090%, and 0.122mil, respectively. Meanwhile, the maximum range error induced by the prediction error of maximum range angle is analyzed to further evaluate the prediction performance, and the MAE, MRE, and RMSE of 5th PRS are 202.86m, 0.089%, and 29.81m, respectively. Compared with prediction results of 5th PRS, the prediction accuracy of 6th PRS decreases significantly. While the order of PRS is more than 6, the available training samples are not enough, resulting in an inaccurate estimation of regression parameters. Therefore, based on the analysis above, we draw the conclusion that the appropriate order is of significant importance to the improvement of the prediction accuracy of PRS.Table 5
Effect of the order of PRS on prediction accuracy of maximum range angle.
Order Maximum range angle Maximum range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) 1 11.441 1.124 4.231 3328.52 1.362 1041.31 2 8.509 0.873 3.077 2501.27 1.059 767.22 3 3.690 0.350 0.692 864.47 0.382 171.82 4 2.605 0.259 0.499 697.64 0.300 124.37 5 0.949 0.090 0.122 202.86 0.089 29.81 6 1.787 0.168 0.174 479.80 0.195 43.12Considering the actual need of predicting the maximum range corresponding to the givenB0, H0, AT, HT, and Ts under standard atmospheric conditions, the functional relation between the maximum range and the above five influencing factors is established directly by using the PRS, and the prediction results are illustrated in Table 6. It is easily seen that the prediction accuracy with the order has the same changing trend as that of PRS for predicting the maximum range angle, and the 5th PRS has the highest prediction accuracy with the MAE, MRE, and RMSE being 746.63m, 0.314%, and 77.39m, respectively.Table 6
Effect of the order of PRS on prediction accuracy of maximum range.
Order MAE (m) MRE (%) RMSE (m) 1 5929.71 2.249 1995.12 2 4084.98 1.456 1285.65 3 2431.36 0.988 369.05 4 1102.10 0.392 220.33 5 746.63 0.314 77.39 6 1028.23 0.410 110.31Regarding the prediction of maximum range angle, Table7 gives the influence of the correlation function on prediction accuracy of the Kriging with the trend function of 2th PRS. It is easily noticed that the Kriging predictors with different correlation functions all demonstrate good prediction performance. In terms of various correlation functions, the Spline ranks first in prediction performance, which illustrates that the Spline correlation function is the most preferable for the prediction of maximum range angle, and it is followed Cubic, Gauss, Matern52, and Matern32 in turn. Furthermore, each of Kriging predictors with above five correlation functions has better prediction performance by comparison with that of 5th PRS under identical conditions. However, the Kriging predictors with the Exp, Linear, and Spherical correlation functions demonstrate similar prediction capability to the 5th PRS.Table 7
Effect of the correlation function of Kriging on prediction accuracy of maximum range angle.
Correlation function Maximum range angle Maximum range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) Cubic 0.318 0.032 0.062 93.58 0.037 15.53 Exp 0.720 0.074 0.121 213.22 0.088 30.26 Gauss 0.385 0.036 0.073 95.81 0.039 18.29 Linear 0.753 0.079 0.122 223.23 0.090 30.40 Matern32 0.613 0.064 0.099 182.36 0.075 24.29 Matern52 0.427 0.040 0.076 103.39 0.043 18.87 Spherical 0.686 0.072 0.121 207.54 0.084 29.95 Spline 0.309 0.029 0.060 58.02 0.025 14.71Table8 shows how the performance of Kriging for predicting the maximum range depends on the correlation function. Obviously, the correlation function has an important influence on the prediction performance, and the prediction performance with the correlation function has the same changing trend as that of Kriging for predicting the maximum range angle. Compared with the 5th PRS, the Kriging predictors with eight correlation functions all have better prediction performances, and especially the prediction performances of Kriging predictors with Cubic and Spline correlation functions are much better than those of Kriging predictors with other correlation functions, which illustrates that the Kriging predictors with Spline and Cubic correlation functions are quite capable of predicting the maximum range.Table 8
Effect of correlation function of Kriging on prediction accuracy of maximum range.
Correlation function MAE (m) MRE (%) RMSE (m) Cubic 61.45 0.025 11.53 Exp 363.70 0.124 59.98 Gauss 135.14 0.057 31.80 Linear 372.36 0.127 60.67 Matern32 286.50 0.104 46.79 Matern52 146.48 0.061 33.42 Spherical 351.18 0.147 59.27 Spline 50.56 0.021 10.56The PRSho (order=5) and Kriging predictors with eight correlation functions are used to construct the EM of maximum range angle, and the corresponding weights are obtained by using the TLBO algorithm with the initial parameters ofNP=22 andNG=500. The optimized weights corresponding to the PRSho, KRcub, KRexp, KRgau, KRlin, KRmat32, KRmat52, KRsph, and KRspl are 0.004214, 0.288490, 0.004953, 0.199488, 0.003182, 0.017866, 0.173068, 0.005553, and 0.303186, respectively. Thus, the KRspl ranks first in the proportion, and it is followed by KRcub, KRgau, and KRmat52 in turn. Moreover, the sum of the weights corresponding to the KRcub, KRgau, KRmat52, and KRspl is 0.964232, which indicates that the above four metamodels play a major role in the construction of the EM. Regarding the prediction of maximum range angle, the MAE, MRE, and RMSE of the EM are 0.259mil, 0.026%, and 0.045mil, respectively, and the maximum range error induced by the prediction error of the maximum range angle is analyzed simultaneously with the MAE, MRE, and RMSE being 50.07m, 0.021%, and 10.93m, respectively. Thus, the EM is clearly superior to other individual metamodels. Based on the above analysis, we can conclude that the KRcub, KRgau, KRmat52, KRspl, and EM have higher prediction accuracies.Similarly, the EM of the maximum range is constructed in the above way, and the optimized weights are obtained, among which the optimized weights corresponding to the KRcub and KRspl are 0.292025 and 0.696638, respectively. With the exception of the KRcub and KRspl, the sum of the optimized weights corresponding to the PRSho and Kriging predictors with other correlation functions is 0.011337, which is mainly due to the fact that the KRcub and KRspl have obviously higher prediction accuracies than those of PRSho and Kriging predictors with other correlation functions, and consequently the prediction accuracy of EM mainly depends on prediction accuracies of the KRcub and KRspl. The MAE, MRE, and RMSE of the EM are 46.84m, 0.019%, and 10.47m, respectively. Obviously, the prediction accuracy of EM is improved slightly compared with that of the KRspl.


### 4.3. Case 3: Firing Angle Calculation under Actual Atmospheric Conditions

Based on the monthly meteorological data below 30km provided by a dozen observation stations during the period from 1990 to 1995, the case study on the firing angle calculation under actual atmospheric conditions is conducted to validate the method described in Section3.3. The meteorological data are adequate in spatial coverage, and it includes various meteorological variables such as temperature, atmospheric pressure, relative humidity, wind velocity, wind direction, virtual temperature, and solar radiation. Under each of varieties of actual atmospheric conditions, the latitude, elevation, and ground surface temperature of the position of the meteorological observation station are selected as the latitude of launcher, the elevation of launcher, and the propellant temperature, respectively, and the LHS is adopted to generate 100 design points for influencing factors AT, RG, and HT with variation ranges being shown in Table 2. If the given RG is greater than the maximum range under actual atmospheric conditions, it only needs to calculate the number of iterations corresponding to the maximum range. Based on the established functional relations in Sections 4.1 and 4.2 with the sample sizes of 9000 and 500, respectively, of which the functional relation between maximum range and B0, H0, AT, HT, Ts is established by using the KRspl, and the distributions of the number of iterations are finally obtained, as illustrated in Figures 7 and 8.Figure 7
Distributions of the number of iterations merely considering the atmospheric pressure and virtual temperature.
(a)
Overall distribution (b)
Partial amplification of P3 portionFigure 8
Distributions of the number of iterations considering the atmospheric pressure, virtual temperature, wind velocity, and wind direction.
(a)
Overall distribution (b)
Partial amplification of P4 portionFigure7 describes the distributions of the number of iterations of five metamodels merely considering the meteorological variables of atmospheric pressure and virtual temperature. Regarding the firing angle calculation corresponding to each of five metamodels, the maximum number of iterations is 4, and the proportions corresponding to the number of iterations 1, 2, 3, and 4 are over 1.49%, 78.19%, 18.05%, and 0.08%, respectively. Thus, the number of iterations 2 and 3 rank first and second in the proportion, and they are followed by the number of iterations 1 and 4 in turn. However, as for the commonly used 4th PRS, the proportions corresponding to the number of iterations 1, 2, 3, and 4 are 1.13%, 47.48%, 51.16%, and 0.23%, respectively. Obviously, the calculation efficiency of each of five metamodels is much higher than that of 4th PRS, and the order of the calculation efficiency from high to low in turn is EM, KRgau, KRspl, PRSho, and KRcub.Figure8 shows the distributions of the number of iterations of five metamodels considering the atmospheric pressure, virtual temperature, wind velocity, and wind direction. It is easily seen that the maximum number of iterations is 5 for each of five metamodels and the proportions corresponding to the number of iterations less than 4 are over 99.88%. Furthermore, the number of iterations 2 and 3 account for the largest proportion (over 94.81%), whereas the proportions corresponding to the number of iterations 1 and 5 are quite small, representing less than 0.96% and 0.12% of all testing samples, respectively. By contrast, with regard to the 4th PRS, the proportions corresponding to the number of iterations 1, 2, 3, 4, and 5 are 0.78%, 32.45%, 60.29%, 6.29%, and 0.19%, respectively. Thus, the efficiency of the firing angle calculation of each of five metamodels is improved significantly by comparison with that of 4th PRS, and the order of the calculation efficiency from high to low in turn is still EM, KRgau, KRspl, PRSho, and KRcub.On the other hand, by comparing the distributions of the number of iterations shown in Figure7, we can see that the maximum number of iterations increases for each of five metamodels if the wind velocity and wind direction are taken into consideration, and the proportions corresponding to the numbers of iterations 1 and 2 are markedly reduced, of which the latter has a reduction of over 30.01%. However, the proportions corresponding to the numbers of iterations 3 and 4 are significantly increased, of which the former has an increment of over 26.41%.


## 4.1. Case 1:Firing Angle Calculation under Standard Atmospheric Conditions

The six DOF trajectory simulation is based on the International Standard Atmosphere (ISO 2533:1975). The design of experiments is carried out for six influencing factors by using the LHS, and the total numbers of training samples and testing samples are 9000 and 5000, respectively. The firing angles corresponding to various design points are obtained by employing the iterative search approach and six DOF trajectory program described in Section3.1.Table3 shows the prediction accuracies of the PRS with different orders. As for the orders of PRS ranging from 1 to 9, the prediction accuracy can be improved significantly with the increase of the order. It is easily seen that the prediction power is relatively low for the commonly used orders of PRS ranging from 1 to 4. Regarding the prediction of the firing angle, the MAE, MRE, and RMSE of the 4th-order PRS are 7.055mil, 0.701%, and 1.424mil, respectively, which indicates that the low-order PRS lacks the capability of approximating the highly complicated functional relation. Compared with the low-order PRS, the high-order PRS has higher prediction accuracy of the firing angle, and the MAE, MRE, and RMSE of the 9th PRS are 1.047mil, 0.101%, and 0.103mil, respectively. Furthermore, in order to further evaluate the prediction accuracies of the PRS with different orders, the range error induced by the prediction error of the firing angle is analyzed simultaneously, which has the same changing trend as that of the firing angle, and the MAE, MRE, and RMSE of the 9th PRS are 291.49m, 0.185%, and 51.08m, respectively. Thus, high-order PRS is preferable to the establishment of a highly nonlinear functional relation between the firing angle and B0, H0, AT, RG, HT, Ts.Table 3
Effect of the order of PRS on prediction accuracy of firing angle.
Order Firing angle Range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) 1 58.774 5.562 13.584 18751.03 12.476 7165.87 2 27.953 3.535 8.638 12062.80 9.439 4817.94 3 11.159 1.383 2.721 4410.95 3.394 1463.24 4 7.055 0.701 1.424 2893.34 1.362 752.29 5 4.789 0.463 0.673 1436.73 0.695 336.77 6 2.976 0.288 0.372 898.42 0.512 180.78 7 2.122 0.204 0.215 589.93 0.232 102.04 8 1.509 0.146 0.137 422.01 0.244 64.78 9 1.047 0.101 0.103 291.49 0.185 51.08 10 2.052 0.248 0.166 1308.71 0.501 94.99In terms of the prediction of firing angle, the correlation function significantly affects the prediction accuracy of Kriging with the trend function of 4th PRS, as illustrated in Table4. Obviously, the Kriging predictors with different correlation functions all have higher prediction accuracies by comparison with the prediction results of the low-order PRS under identical conditions. Nevertheless, the correlation function still plays a significant role in Kriging. The Gauss, Spline and Cubic correlation functions of great interest demonstrate a promising capability to predict the firing angle, which have a parabolic behavior near the origin. In particular, the Gauss and Spline correlation functions are quite similar in terms of the prediction accuracy. However, if Linear, Exp, and Spherical functions with a linear behavior near the origin are selected as the correlation functions, the prediction performances are relatively poor.Table 4
Effect of the correlation function of Kriging on prediction accuracy of firing angle.
Correlation function Firing angle Range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) Cubic 0.809 0.080 0.099 286.63 0.108 51.45 Exp 2.242 0.223 0.227 664.52 0.250 85.75 Gauss 0.865 0.083 0.084 243.33 0.094 39.71 Linear 2.426 0.233 0.233 674.98 0.262 93.47 Matern32 1.667 0.161 0.150 465.00 0.180 60.81 Matern52 1.192 0.115 0.121 332.74 0.152 53.66 Spherical 2.330 0.224 0.221 649.45 0.251 83.96 Spline 0.738 0.074 0.087 262.95 0.103 44.97The high-order PRS (PRSho, order=9) and Kriging predictors with Cubic (KRcub), Exp (KRexp), Gauss (KRgau), Linear (KRlin), Matern32 (KRmat32), Matern52 (KRmat52), Spherical (KRsph), and Spline (KRspl) correlation functions are used to construct the EM. The weights of the EM are calculated by using the TLBO algorithm and initialize the optimization parameters of the population sizeNP=30 and the number of generationsNG=600. The optimized weights corresponding to the PRSho, KRcub, KRexp, KRgau, KRlin, KRmat32, KRmat52, KRsph, and KRspl are 0.139072, 0.146959, 0.000076, 0.363127, 0.000073, 0.000255, 0.000681, 0.000082, and 0.349675, respectively, and the sum of the weights corresponding to the PRSho, KRcub, KRgau, and KRspl is 0.998833, which indicates that the above four metamodels play a key role in the construction of the EM. With regard to the prediction of firing angle, the MAE, MRE, and RMSE of the EM are 0.698 mil, 0.068%, and 0.071mil, respectively. In terms of the range error induced by the prediction error of firing angle, the MAE, MRE, and RMSE of the EM are 195.08m, 0.076%, and 35.92m, respectively. Obviously, the EM outperforms other individual metamodels. In light of the above analysis, we can conclude that the PRSho, KRcub, KRgau, KRspl, and EM have higher prediction accuracies.The distributions of absolute errors of the firing angle and range of five metamodels are shown in Figures2 and 3. According to Figure 2, for each of the five metamodels, the testing samples with the absolute error of firing angle being less than 0.15mil and 0.30mil account for over 91.16% and 98.24%, respectively. As illustrated in Figure 3, the proportions of testing samples with the absolute error of range being less than 50m and 100m makes up over 70.00% and 95.48%, respectively. The above results further verify the effectiveness of the firing angle calculation method under standard atmospheric conditions.Figure 2
Distributions of the absolute error of firing angle of five metamodels.
(a)
Overall distribution (b)
Partial amplification of P1 portionFigure 3
Distributions of the absolute error of range of five metamodels.
(a)
Overall distribution (b)
Partial amplification of P2 portionThe execution time of a program depends strongly on the programming environment and the computing platform used. Based on the established metamodels PRSho, KRcub, KRgau, KRspl, and EM for the prediction of firing angle, five programs are written in Microsoft Visual C++ 6.0. Each program operates on a PC with an Intel(R) Core(TM) i5-6300HQ CPU @ 2.30GHz, 4GB RAM, Windows 7 (64bit). The execution times of various metamodels for a calculation of the firing angle are given in Figure4, where the execution time of each metamodel is the average time of the execution times of 5000 testing cases. It is quite obvious that the execution times of four individual metamodels are all less than 0.9ms, and the PRSho and KRspl are the fastest and slowest metamodels with execution times being 0.03676ms and 0.82525ms, respectively. Regarding the execution time of the EM, it is 2.23317ms, which is still far less than the required time of the fire control system that is measured in seconds. Therefore, based on the above analysis of the rapidity of the firing angle prediction, the feasibility of the firing angle calculation method under standard atmospheric conditions is verified.Figure 4
Execution times of various metamodels for a prediction of firing angle.Figures5 and 6 clearly illustrate how the prediction performances of metamodels depend on the training sample size. The training samples with different sizes are generated independently based on the LHS. The PRSho predictors corresponding to the training sample size ranging from 3000 to 21000 with an increment of 3000 have the orders of 7, 8, 9, 9, 10, 10, and 10, respectively. The EMs of the firing angle corresponding to different training samples are constructed in the same way as the EM with 9000 training samples, and thus the optimized weights corresponding to different training samples are different. Regarding the prediction error of the firing angle and its inducing range error, not surprisingly, the RMSE of each metamodel shows a downward trend with the training sample size. In particular, if the number of the training samples is less than 9000, the RMSE of each metamodel decreases significantly with the increase of training sample size, which indicates that the prediction performance can be improved remarkably with an increasing training sample size, and the PRSho has the most significant decrease. However, while the training sample size is more than 9000, the RMSE of each metamodel decreases relatively slowly, which indicates that 9000 training samples are enough to achieve reasonably good prediction performance. The MAE and MRE of each metamodel have the similar changing trends, and the KRgau is then taken as an example. Regarding the prediction of firing angle, the MAE and MRE of the KRgau are 1.494mil and 0.150%, respectively, at a training sample size of 3000, but they decrease to 0.865mil and 0.083% with 9000 training samples and further decrease to 0.451mil and 0.043% with 21000 training samples.Figure 5
Variation of the prediction accuracy of firing angle with the training sample size.Figure 6
Variation of the range error due to the prediction error of firing angle with the training sample size.


## 4.2. Case 2: Calculations of Maximum Range Angle and Maximum Range under Standard Atmospheric Conditions

The research on calculation methods of the maximum range angle and maximum range under standard atmospheric conditions is carried out by using the PRS, Kriging, and EM. The LHS is selected as the method for the design of experiments, and the variation ranges of influencing factorsB0, H0, AT, HT, and Ts are shown in Table 2. Then the training samples and testing samples are obtained with the sizes of 500 and 1000, respectively. Thereafter, the functional relations between maximum range angle, maximum range, and B0, H0, AT, HT, Ts are established.The effect of the order of PRS on prediction accuracy of maximum range angle is illustrated in Table5. Obviously, for the order within the range of 1 to 5, the prediction performance can be improved substantially with an increasing order, and the 5th PRS has the best prediction performance. In terms of the prediction of maximum range angle, the MAE, MRE, and RMSE of 5th PRS are 0.949mil, 0.090%, and 0.122mil, respectively. Meanwhile, the maximum range error induced by the prediction error of maximum range angle is analyzed to further evaluate the prediction performance, and the MAE, MRE, and RMSE of 5th PRS are 202.86m, 0.089%, and 29.81m, respectively. Compared with prediction results of 5th PRS, the prediction accuracy of 6th PRS decreases significantly. While the order of PRS is more than 6, the available training samples are not enough, resulting in an inaccurate estimation of regression parameters. Therefore, based on the analysis above, we draw the conclusion that the appropriate order is of significant importance to the improvement of the prediction accuracy of PRS.Table 5
Effect of the order of PRS on prediction accuracy of maximum range angle.
Order Maximum range angle Maximum range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) 1 11.441 1.124 4.231 3328.52 1.362 1041.31 2 8.509 0.873 3.077 2501.27 1.059 767.22 3 3.690 0.350 0.692 864.47 0.382 171.82 4 2.605 0.259 0.499 697.64 0.300 124.37 5 0.949 0.090 0.122 202.86 0.089 29.81 6 1.787 0.168 0.174 479.80 0.195 43.12Considering the actual need of predicting the maximum range corresponding to the givenB0, H0, AT, HT, and Ts under standard atmospheric conditions, the functional relation between the maximum range and the above five influencing factors is established directly by using the PRS, and the prediction results are illustrated in Table 6. It is easily seen that the prediction accuracy with the order has the same changing trend as that of PRS for predicting the maximum range angle, and the 5th PRS has the highest prediction accuracy with the MAE, MRE, and RMSE being 746.63m, 0.314%, and 77.39m, respectively.Table 6
Effect of the order of PRS on prediction accuracy of maximum range.
Order MAE (m) MRE (%) RMSE (m) 1 5929.71 2.249 1995.12 2 4084.98 1.456 1285.65 3 2431.36 0.988 369.05 4 1102.10 0.392 220.33 5 746.63 0.314 77.39 6 1028.23 0.410 110.31Regarding the prediction of maximum range angle, Table7 gives the influence of the correlation function on prediction accuracy of the Kriging with the trend function of 2th PRS. It is easily noticed that the Kriging predictors with different correlation functions all demonstrate good prediction performance. In terms of various correlation functions, the Spline ranks first in prediction performance, which illustrates that the Spline correlation function is the most preferable for the prediction of maximum range angle, and it is followed Cubic, Gauss, Matern52, and Matern32 in turn. Furthermore, each of Kriging predictors with above five correlation functions has better prediction performance by comparison with that of 5th PRS under identical conditions. However, the Kriging predictors with the Exp, Linear, and Spherical correlation functions demonstrate similar prediction capability to the 5th PRS.Table 7
Effect of the correlation function of Kriging on prediction accuracy of maximum range angle.
Correlation function Maximum range angle Maximum range MAE (mil) MRE (%) RMSE (mil) MAE (m) MRE (%) RMSE (m) Cubic 0.318 0.032 0.062 93.58 0.037 15.53 Exp 0.720 0.074 0.121 213.22 0.088 30.26 Gauss 0.385 0.036 0.073 95.81 0.039 18.29 Linear 0.753 0.079 0.122 223.23 0.090 30.40 Matern32 0.613 0.064 0.099 182.36 0.075 24.29 Matern52 0.427 0.040 0.076 103.39 0.043 18.87 Spherical 0.686 0.072 0.121 207.54 0.084 29.95 Spline 0.309 0.029 0.060 58.02 0.025 14.71Table8 shows how the performance of Kriging for predicting the maximum range depends on the correlation function. Obviously, the correlation function has an important influence on the prediction performance, and the prediction performance with the correlation function has the same changing trend as that of Kriging for predicting the maximum range angle. Compared with the 5th PRS, the Kriging predictors with eight correlation functions all have better prediction performances, and especially the prediction performances of Kriging predictors with Cubic and Spline correlation functions are much better than those of Kriging predictors with other correlation functions, which illustrates that the Kriging predictors with Spline and Cubic correlation functions are quite capable of predicting the maximum range.Table 8
Effect of correlation function of Kriging on prediction accuracy of maximum range.
Correlation function MAE (m) MRE (%) RMSE (m) Cubic 61.45 0.025 11.53 Exp 363.70 0.124 59.98 Gauss 135.14 0.057 31.80 Linear 372.36 0.127 60.67 Matern32 286.50 0.104 46.79 Matern52 146.48 0.061 33.42 Spherical 351.18 0.147 59.27 Spline 50.56 0.021 10.56The PRSho (order=5) and Kriging predictors with eight correlation functions are used to construct the EM of maximum range angle, and the corresponding weights are obtained by using the TLBO algorithm with the initial parameters ofNP=22 andNG=500. The optimized weights corresponding to the PRSho, KRcub, KRexp, KRgau, KRlin, KRmat32, KRmat52, KRsph, and KRspl are 0.004214, 0.288490, 0.004953, 0.199488, 0.003182, 0.017866, 0.173068, 0.005553, and 0.303186, respectively. Thus, the KRspl ranks first in the proportion, and it is followed by KRcub, KRgau, and KRmat52 in turn. Moreover, the sum of the weights corresponding to the KRcub, KRgau, KRmat52, and KRspl is 0.964232, which indicates that the above four metamodels play a major role in the construction of the EM. Regarding the prediction of maximum range angle, the MAE, MRE, and RMSE of the EM are 0.259mil, 0.026%, and 0.045mil, respectively, and the maximum range error induced by the prediction error of the maximum range angle is analyzed simultaneously with the MAE, MRE, and RMSE being 50.07m, 0.021%, and 10.93m, respectively. Thus, the EM is clearly superior to other individual metamodels. Based on the above analysis, we can conclude that the KRcub, KRgau, KRmat52, KRspl, and EM have higher prediction accuracies.Similarly, the EM of the maximum range is constructed in the above way, and the optimized weights are obtained, among which the optimized weights corresponding to the KRcub and KRspl are 0.292025 and 0.696638, respectively. With the exception of the KRcub and KRspl, the sum of the optimized weights corresponding to the PRSho and Kriging predictors with other correlation functions is 0.011337, which is mainly due to the fact that the KRcub and KRspl have obviously higher prediction accuracies than those of PRSho and Kriging predictors with other correlation functions, and consequently the prediction accuracy of EM mainly depends on prediction accuracies of the KRcub and KRspl. The MAE, MRE, and RMSE of the EM are 46.84m, 0.019%, and 10.47m, respectively. Obviously, the prediction accuracy of EM is improved slightly compared with that of the KRspl.


## 4.3. Case 3: Firing Angle Calculation under Actual Atmospheric Conditions

Based on the monthly meteorological data below 30km provided by a dozen observation stations during the period from 1990 to 1995, the case study on the firing angle calculation under actual atmospheric conditions is conducted to validate the method described in Section3.3. The meteorological data are adequate in spatial coverage, and it includes various meteorological variables such as temperature, atmospheric pressure, relative humidity, wind velocity, wind direction, virtual temperature, and solar radiation. Under each of varieties of actual atmospheric conditions, the latitude, elevation, and ground surface temperature of the position of the meteorological observation station are selected as the latitude of launcher, the elevation of launcher, and the propellant temperature, respectively, and the LHS is adopted to generate 100 design points for influencing factors AT, RG, and HT with variation ranges being shown in Table 2. If the given RG is greater than the maximum range under actual atmospheric conditions, it only needs to calculate the number of iterations corresponding to the maximum range. Based on the established functional relations in Sections 4.1 and 4.2 with the sample sizes of 9000 and 500, respectively, of which the functional relation between maximum range and B0, H0, AT, HT, Ts is established by using the KRspl, and the distributions of the number of iterations are finally obtained, as illustrated in Figures 7 and 8.Figure 7
Distributions of the number of iterations merely considering the atmospheric pressure and virtual temperature.
(a)
Overall distribution (b)
Partial amplification of P3 portionFigure 8
Distributions of the number of iterations considering the atmospheric pressure, virtual temperature, wind velocity, and wind direction.
(a)
Overall distribution (b)
Partial amplification of P4 portionFigure7 describes the distributions of the number of iterations of five metamodels merely considering the meteorological variables of atmospheric pressure and virtual temperature. Regarding the firing angle calculation corresponding to each of five metamodels, the maximum number of iterations is 4, and the proportions corresponding to the number of iterations 1, 2, 3, and 4 are over 1.49%, 78.19%, 18.05%, and 0.08%, respectively. Thus, the number of iterations 2 and 3 rank first and second in the proportion, and they are followed by the number of iterations 1 and 4 in turn. However, as for the commonly used 4th PRS, the proportions corresponding to the number of iterations 1, 2, 3, and 4 are 1.13%, 47.48%, 51.16%, and 0.23%, respectively. Obviously, the calculation efficiency of each of five metamodels is much higher than that of 4th PRS, and the order of the calculation efficiency from high to low in turn is EM, KRgau, KRspl, PRSho, and KRcub.Figure8 shows the distributions of the number of iterations of five metamodels considering the atmospheric pressure, virtual temperature, wind velocity, and wind direction. It is easily seen that the maximum number of iterations is 5 for each of five metamodels and the proportions corresponding to the number of iterations less than 4 are over 99.88%. Furthermore, the number of iterations 2 and 3 account for the largest proportion (over 94.81%), whereas the proportions corresponding to the number of iterations 1 and 5 are quite small, representing less than 0.96% and 0.12% of all testing samples, respectively. By contrast, with regard to the 4th PRS, the proportions corresponding to the number of iterations 1, 2, 3, 4, and 5 are 0.78%, 32.45%, 60.29%, 6.29%, and 0.19%, respectively. Thus, the efficiency of the firing angle calculation of each of five metamodels is improved significantly by comparison with that of 4th PRS, and the order of the calculation efficiency from high to low in turn is still EM, KRgau, KRspl, PRSho, and KRcub.On the other hand, by comparing the distributions of the number of iterations shown in Figure7, we can see that the maximum number of iterations increases for each of five metamodels if the wind velocity and wind direction are taken into consideration, and the proportions corresponding to the numbers of iterations 1 and 2 are markedly reduced, of which the latter has a reduction of over 30.01%. However, the proportions corresponding to the numbers of iterations 3 and 4 are significantly increased, of which the former has an increment of over 26.41%.


## 5. Conclusion

Based on the large sample data, the functional relation between the firing angle and six influencing factors is established by using the PRS, Kriging and EM, which effectively deals with the problem of the rapid calculation of firing angle for the MLRS under standard atmospheric conditions. On this basis, together with the established functional relation between the maximum range and five influencing factors, by utilizing the proposed iterative search approach and six DOF trajectory program, the problem of firing angle calculation under actual atmospheric conditions has been effectively solved as well. The main conclusions of this study are as follows.( 1 ) Regarding the firing angle calculation under standard atmospheric conditions, the PRSho, KRcub, KRgau, KRspl, and EM have higher prediction accuracies by comparison with those of other metamodels, and the order of the prediction accuracy at a training sample size of 9000 from high to low in turn is EM, KRgau, KRspl, KRcub, and PRSho. In addition, for the training sample size less than 9000, the prediction accuracy of each of the above five metamodels improves significantly with the increase of training sample size, whereas it improves slightly when the training sample size is over 9000.( 2 ) The KRcub, KRgau, KRmat52, KRspl, and EM have better performances for predicting the maximum range angle under standard atmospheric conditions compared with those of other metamodels at a training sample size of 500, and the KRcub and KRspl are especially suitable for predicting the maximum range.( 3 ) The methods for calculating the firing angle and maximum range angle based on metamodels feature high accuracy and fast speed, which can provide a new solution to the calculation of the firing data for the MLRS under standard atmospheric conditions.( 4 ) The proposed method of firing angle calculation under actual atmospheric conditions can effectively reduce the number of iterations, and thus the rapid reaction capability of the MLRS is improved significantly.With the rapid development of technologies such as artificial intelligence and machine learning, the research and application of the firing angle calculation for the MLRS will also gain significant developments in the near future.



---

*Source: 7689860-2019-04-04.xml*