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.
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 launcher , the elevation of launcher , the target azimuth , the range between the launcher and target , the elevation of target , and the propellant temperature . 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 aswhere is the predicted value at the point , is the ith component of the vector of length n, 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:where is a trend function providing global approximation and is the realization of a stochastic process with mean 0 and variance . The spatial covariance for design points and is defined as where is a correlation function with the parameter vector , expressed aswhere dj=wj−xj and . 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 (Matern32) and (Matern52), for which where .
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 where is the predicted response by the EM, M denotes the number of individual metamodels, is the predicted response by the ith individual metamodel, and denotes the weight corresponding to the ith 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, .
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 aswhere denotes the chosen error metric for measuring the accuracy of the response predicted by the EM and is the actual response value at the point . The root mean squared error is selected as the error metric . 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: where and denote the populations after and before learning from the teacher respectively; and are the random numbers between 0 and 1; denotes the current best individual; denotes the teaching factor with the value being either 1 or 2; 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 learner can be defined aswhere and are two randomly selected learners; and are the fitness values of the learners and , respectively; 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,where 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 Figure 1, and the following two specific steps are involved.

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 and , , , , , 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. Calculate and 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 and , , , , , .
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 corresponding to the given and by using the Lagrange method [37] and take as the initial value of the iteration.
Step 2. Calculate the range corresponding to the via the six DOF trajectory program, as well as the difference between and by .
Step 3. Estimate the step length according to and calculate the firing angle of the second iteration by .
Step 4. Calculate the range corresponding to the via the six DOF trajectory program.
Step 5. Determine the firing angle with higher accuracy according to the following iteration formula:where , , and are the firing angles at the kth, (k-1)th and (k-2)th iterations, respectively; and are the ranges corresponding to firing angles and , 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:where is the range corresponding to . Thereafter, is replaced by 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 that is greater than the maximum range corresponding to the given , , , , and . 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 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 by , , , and . 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 the ith round, let , , and be the initial value of the iteration, the step size, and the assumed firing angle of the second iteration, respectively, where , and the ranges corresponding to and are and , respectively. Then the kth firing angle can be expressed aswhere k denotes the number of iterations. For each round, iterate formula (17), until the maximum range angle is reached. For the subsequent round, 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 for the first iteration is obtained by using the established functional relation between the firing angle and various influencing factors.
Let be the range corresponding to under actual atmospheric conditions, and is the firing angle corresponding to under standard atmospheric conditions. Let be the firing angle corresponding to the range under actual atmospheric conditions, and is the range corresponding to under standard atmospheric conditions. Two basic assumptions are given as follows.
Assumption 1. The difference between and is equal to the difference between and ; namely, .
Assumption 2. The difference between and is equal to the difference between and ; namely, .
However, there is an unavoidable case that is near the minimum or maximum ranges, which may result in and being less than the minimum range or greater than the maximum range under standard atmospheric conditions. In this case, if metamodels are still used to predict , it may lead to the sharp decrease of the prediction accuracy. Aiming at solving this problem, an efficient method for calculating according to the ranges and is proposed, as shown below.
If there is or such that , then follow Assumption 1 to calculate by as well as the corresponding firing angle under standard atmospheric conditions, and subsequently is obtained by . Similarly, if there exists or such that , then follow Assumption 2 to calculate the directly by . Regarding other cases where and , is then calculated by based on Assumptions 1 and 2. To sum up, we get the following expression for ,
The 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 corresponding to the given , , , , , and .
Step 2. By utilizing the six DOF trajectory program, calculate the range corresponding to under actual atmospheric conditions.
Step 3. Based on the established functional relation between the maximum range and various influencing factors, calculate the maximum range corresponding to the given , , , , and .
Step 4. Check whether satisfies , and if yes, calculate the firing angle corresponding to based on the established functional relation between the firing angle and various influencing factors.
Step 5. Calculate the distance difference between and by , and then calculate the range by .
Step 6. Check whether satisfies , and if yes, calculate the firing angle corresponding to the based on the established functional relation between the firing angle and various influencing factors.
Step 7. Calculate by the following equation: and then calculate the range corresponding to under actual atmospheric conditions by utilizing the six DOF trajectory program.
Step 8. Having obtained the iterative initial values , and their corresponding ranges , , 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. Table 2 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.
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 Section 3.1.
Table 3 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 , , , , , .
In 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 Table 4. 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.
The 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 size NP=30 and the number of generations NG=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 Figures 2 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.

(a) Overall distribution

(b) Partial amplification of P1 portion

(a) Overall distribution

(b) Partial amplification of P2 portion
The 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 Figure 4, 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.

Figures 5 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.


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 factors , , , , and 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 , , , , are established.
The effect of the order of PRS on prediction accuracy of maximum range angle is illustrated in Table 5. 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.
Considering the actual need of predicting the maximum range corresponding to the given , , , , and 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.
Regarding the prediction of maximum range angle, Table 7 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 8 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.
The 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 of NP=22 and NG=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 Section 3.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 , , and with variation ranges being shown in Table 2. If the given 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 , , , , is established by using the KRspl, and the distributions of the number of iterations are finally obtained, as illustrated in Figures 7 and 8.

(a) Overall distribution

(b) Partial amplification of P3 portion

(a) Overall distribution

(b) Partial amplification of P4 portion
Figure 7 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.
Figure 8 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 Figure 7, 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.
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.
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.
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.
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.
Data Availability
We regret that access to data is restricted for the time being. The data used to support the findings of this study have not been made available because the protection of technical privacy and confidentiality.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by Postgraduate Research & Practice Innovation Program of Jiangsu Province (Grant no. KYCX17_0392).