Overall Parameters Design of Air-Launched Rockets Using Surrogate Based Optimization Method

: In this paper, the design and optimization method of rocket parameters based on the surrogate model and the trajectory simulation system of the 3-DOF air-launched rockets were established. The Gaussian kernel width determination method based on the relationship between local density and width is used to ensure the efﬁciency and reliability of the optimization method, and at the same time greatly reduces the amount of calculation. An adaptive sampling point updating method was established, which includes three stages: location sampling, exploration sampling, and potential optimal sampling of the potential feasible region. The adaptive sampling is realized by the distance constraint. Based on the precision of the surrogate model, the convergence end criterion was established, which can achieve efﬁcient and reliable probabilistic global optimization. The objective function of the optimization problem was deduced to determine the maximum load mass and reasonable constraints were set to ensure that the rocket could successfully enter orbit. For solid engine rockets with the same take-off mass as Launcherone, the launch altitude and target orbit were optimized and analyzed, and veriﬁed by 3-DOF trajectory simulation. The surrogate-based optimization algorithm solved the problem of the overall parameter design optimization of the air-launched rocket and it provides support for the design of air-launched solid rockets.


Introduction
Air-launched rocket refers to relying on large aircraft (large transport aircraft, cargo aircraft, strategic bombers, aircraft, etc.) as an air-launch platform, which carries the rocket to high altitude to release and launch, after a certain attitude adjustment, the rocket engine ignites and transports the load to the preset space orbit [1]. According to the way the rocket is loaded on the aircraft, it can be divided into four kinds: hanging type, back support type, built-in type, and towed type [2]. Compared with the land launched rocket, the high-altitude launched rocket has the following advantages: (1) low dependence on ground equipment, there is no need to consider the limitation of the ground launch site and the huge infrastructure on the ground; (2) shortened launch preparation time, there is no need to consider the launch site and no need to rely on large ground facilities, so the preparation time for all launches is greatly reduced; (3) flexible launch conditions, high-altitude launch greatly reduces the impact of meteorological conditions on the launch, so it can be launched at any time on demand; (4) reduce launch consumption, the consumption of manpower and material resources for high-altitude launch is far less than that of ground launch. In particular, fuel consumption can be reduced, which can reduce launch costs [3].
In this paper, we mainly consider using large transport aircraft to launch rockets, which are launched under the condition of certain altitude and initial velocity. In 1990, America's Orbital Sciences Corporation had successfully launched its "Pegasus" rocket for the first time, this was the world's first air-launched rocket to launch a satellite. Russian aerospace company is also developing the "Flying" air-launched rocket and Boeing from the United States is also planning to develop an air-launched rocket [4]. According to relevant research data, the total mass of airborne rockets can be reduced by 20∼30% [5] compared with that of ground-launched rockets under the same load capacity. Recently, Virgin Orbit launched the Launcherone rocket with a Boeing 747, which had 300 kg of sun-synchronous orbit carrying capacity and the capacity of 500 kg of low earth orbit. Given the present research status of air-launched rockets, the research on the design and optimization of the overall parameters of air-launched rockets is still lacking. To solve this problem, this paper established the reliable air-launched rocket launch system and trajectory simulation model [6], the surrogate-based optimization method was adopted, the design and optimization of the overall parameters (engine mass, flight time, trajectory) of the solid rocket were carried out with reference to the LauncherOne rocket with a launch altitude of 10,668 m and a takeoff mass of 25.85 t. The target orbit is a 500 km Sun-synchronous orbit (SSO). The objective function is the maximum capacity of the rocket, which is analyzed by the results of algorithm optimization and three-degree-offreedom(3-DOF) trajectory simulation.
In the development of optimization algorithms in recent years, simulated annealing algorithm, genetic algorithm, ant colony algorithm, and intelligent optimization algorithm have developed more mature [7][8][9][10][11]. However, in general, this kind of optimization algorithm will appear a lot of repeated calculation of the objective function and calculation constraints, which greatly affects the optimization efficiency. So the reduction of the number of calls from the original mode is very important for the efficiency of the optimization method. For this kind of problem, response surface methods and sequence approximate optimization methods have been developed gradually. To a large extent, these methods provide a reference for optimization methods that call the original model less frequently [12,13]. More commonly in the sequence, the approximate optimization method is based on the surrogate-based optimization method, which is the common way to solve the problem of multidisciplinary optimization.
The overall design of aircraft is a typical multidisciplinary coupled problem. The complex coupling relationship among disciplines makes it difficult for traditional singledisciplinary design optimization methods to find the optimal solution. The Multidisciplinary Design Optimization (MDO) method is a methodology that fully explores and utilizes the synergistic mechanism of interaction in the system to design complex engineering systems and subsystems, it is the process of finding the optimal multidisciplinary design through optimization and compromise under each disciplinary model [14][15][16].
Cheng et al. [17] proposed a cooperative coevolution MDO algorithm based on the combination of cooperative coevolution and MDO algorithm and applied it to the design and optimization process of the overall parameters of the missile. The minimum take-off mass of the missile was taken as an indicator to optimize the parameters satisfying the aerodynamic, engine, and control requirements of the missile. Zhang et al. [18] proposed a multidisciplinary optimization method with an improved hybrid genetic algorithm. By establishing a multidisciplinary model integrating the aerodynamic, propulsion, and trajectory of the missile for analysis, and taking the range as the optimization index, the overall parameters of the missile satisfying the constraints were obtained. Although the MDO method used in the above literature considers the correlation between various disciplines in the overall design of the missile, it does not apply to the initial stage of the conceptual design of the missile due to the variety of disciplines involved and the difficulty and complexity of model establishment. To solve the above problems, Zhou et al. [19] proposed a multi-attribute comprehensive optimization design method based on variable weight TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) evaluation and particle swarm optimization algorithm fusion and took Euclidean distance as the optimization evaluation index to optimize the overall parameters of the missile. This method can effectively solve the problem of multi-disciplinary, multi-attribute, and multi-objective parameter optimization in the process of missile overall parameter design. With the innovation of approximation technology and parallel computing technology, the traditional MDO method exposes problems such as complex discipline model, large information exchange capacity, high computing cost, complex organization form, nonlinear design space, and rigid interactive environment, and so on [20]. To improve the efficiency of MDO for complex design problems, it is necessary to establish a clear discipline variable relationship network, an accurate discipline analysis model, an efficient MDO solution strategy, and a coherent discipline interface. Due to the high computational cost of some disciplines, it is necessary to replace the discipline model with advanced surrogate model technology to save optimization time and cost.
The traditional surrogate models mainly include polynomial model, Kriging model and neural network model. Polynomial model is suitable for low order nonlinear problems modeling; Kriging model is suitable for high-dimensional and low-order nonlinear problem modeling; Neural network model is suitable for high-order nonlinear problems modeling [21][22][23]. In recent years, new surrogate models suitable for actual objects have been gradually developed. By simplifying and forming low-precision calculation samples, a low-precision surrogate model can be obtained, and then a high-precision surrogate model can be obtained through high-precision calculation by updating the samples [24,25]. This can not only reduce the amount of calculation, but also make the surrogate model more accurate by using high-precision samples.

Analysis of Sampling Point Updating Method
The exploration and development ability of the surrogate-based optimization method depends on the updating method of sampling points. A more reasonable updating method of sampling points is the key to ensuring the convergence speed of the optimization method under the premise of global optimization. At present, the commonly used sampling point updating methods include EI criterion [26] (Maximum expectation improvement infill criterion), potential optimal infill criterion, maximum distance infill criterion, maximum curvature criterion, and adaptive distance constrained sampling strategy, etc.
Firstly, based on the optimal Latin Hypercube (OLH), an efficient Design of Experiment (DoE) approach for constrained regions is used to form the initial sampling points from the design variable interval. Secondly, the surrogate model is established based on RPF. Finally, the multi-stage sampling point criteria and convergence criteria are used to update the sampling points.
In terms of the exploration and development ability of the optimization method, the maximum distance infill criterion, and the potential optimal infill criterion have shown superior performance, while the adaptive distance constraint strategy can balance the exploration and development ability of the surrogate-based optimization method to a certain extent. The mathematical models and performance characteristics of these three updating methods are as follows.
Optimization problems can generally be expressed in the following form: where, x is an real column vector, f (x) is the original objective function, g i (x) is inequality constraint, h j (x) is equality constraint, l and m represent the number of inequality constraints and equality constraints respectively. If they are zero, then there are no inequality constraints or equality constraints.
The sampling point of the potential optimal infilling criterion is the solution of the following optimization problem: where, f (x), g i (x) and h j (x) are the approximate models of the original objective function f (x), inequality constraint g i (x), and equality constraint h j (x). This method can accelerate the convergence of the optimization method, but it is difficult to ensure that it is the global optimal solution. The sampling point of the maximum distance criterion is the solution of the following optimization problem: where, d(x) is the minimum Euclidean distance between the new sample point and the original sample point, and N is the number of existing sample points. This method can modify the surrogate-based optimization method to a strong ability of exploration, but it is difficult to converge globally. The adaptive sample updating strategy [27,28] generally adds the optimal solution as a new sampling point to the sample by solving the specific optimization problem, and the principle is as follows: New sampling points are obtained by solving the above equation, where, δ is the minimum distance constraint among sampling points. It can be obtained by N is the number of sample points, and x i , x j are the coordinates of sampling points.

Multi-Stage Adaptive Sampling Point Updating Method Model
The multi-stage adaptive sampling points updating method proposed in this paper is based on the surrogate model and the sequence approximate optimization method, which can greatly improve the optimization efficiency of the optimization method under the premise of ensuring the probabilistic global optimum. The main steps are as follows: (1) The optimization problem of the multi-stage adaptive sampling points updating method is established. The known optimization problems are as follows: where, N is the sample number; p d (x) is the distance constraint function; δ constraints the distance between the new sample points and the existing sample points, and the result is as follows: where, dim is the dimension of the optimization problem; 1/ dim √ N is the average distance of the existing sample points; n is the number of steps updating at the sampling points; x a and x b are coordinates of the updating sampling points. According to the above equation, the distance constraint δ will gradually decrease with the increase in samples and can be kept within a reasonable value range. In this way, the self-adaptive adjustment of the distance constraint between the new sampling points and the existing sample points can be realized.
When transforming constrained optimization problems into unconstrained optimization problems, the multiplier method or penalty function method is generally considered [29], as shown below: where, α and β are penalty factors, ϕ, ψ and φ are continuous functions satisfying the following conditions: where, y represents equality constraint or inequality constraint or distance constraint. y ≥ 0 and y = 0 mean that the constraints are satisfied, y < 0 and y = 0 mean that these constraints are not satisfied. The optimization method proposed in this paper is to add new sampling points to the original sample, which are the optimal solutions of the unconstrained optimization problem corresponding to Equation (7). The unconstrained optimization problem also includes the penalty function term α ∑ i ϕ[ g i (x)]+∑ j ψ h j (x) of constraints g i (x), and h j (x), and also the penalty function term βφ[p d (x)] of adaptive distance constraint p d (x) corresponding to the original optimization problem. The ability to explore and develop the surrogate-based optimization can be changed by adjusting the penalty factor α and β.
(2) Establish the accuracy evaluation model of the surrogate model of the objective function and the constraint value relative to the original model Suppose x N+1 is the N + 1 th sampling point obtained by updating the sampling points. The approximate models based on the first N sample points are respectively f (x), g i (x) and h j (x), then the accuracy evaluation model of the approximate model of the objective function and the constraint value relative to the original model is as follows: In the Formula (9), when the independent variable x is greater than or equal to zero, η, λ and µ are non-negative functions, and η(0) = 0, λ(0) = 0, µ(0) = 0.
(3) First sampling stage-Location sampling of potential feasible region The value of penalty factor α and β are set to make 10 3 to 10 5 times, when sampling is conducted in the potential feasible region.
The main steps for updating the N th sampling point in the potentially feasible region are as follows: (a) firstly, an approximate surrogate model based on the objective function and constraint values of the existing sample points is constructed; (b) according to Equation (7), the optimal solution is obtained, namely, the coordinates of the new sampling point x 1n are obtained; (c) calculate the objective function and constraint value according to the coordinates of the new sampling points, the surrogate model and the original model. The constraint value can be expressed as ε(x 1n ). The sampling updating process is repeated until the new sampling point satisfies any of the following criteria, and then the localization sampling phase of the potentially feasible region is completed.
In the above criteria, ε(x 1n ) is calculated according to Equation (9), ε 1 , λ 1 , and I 1 are preset parameters in the location sampling stage and related to the criteria at the end. Criterion I means that ε(x 1n ) value of the nearest λ 1 sample points are less than or equal to ε 1 ; Criterion II means that the number of sampling points in the potential feasible region is more than or equal to I 1 .
In the sampling stage of the potential feasible region, the sampling points first need to satisfy the adaptive distance constraints p d (x) ≥ 0. On the premise of satisfying the distance constraints, the sampling points will tend to the approximate model satisfying the inequality constraints g i (x) ≤ 0 and the equality constraints h j (x) = 0, which indicates that the sampling points will be distributed in the sparse region and tend to the potential feasible region. The sampling stage of the potential feasible region is to explore the whole solution domain in the sampling process, and gradually tends to the potential feasible region. In this way, the approximate accuracy of the entire surrogate model in the design space can be gradually improved, and the positioning of the potential feasible region can be well realized. The end criterions of sampling stage based on the objective function and the constraint value of the approximate model accuracy have the characteristics of strong maneuverability and can use a reasonable number of sampling points to end the sampling stage of the potential feasible region, which greatly improves the optimization efficiency of the surrogate-based optimization method.
(4) Second sampling stage-potential feasible region exploration sampling The second sampling stage is the exploratory sampling under the potential feasible region. In this process, the value of penalty factor α is set to constant value. The value of have the same order of magnitude through adjusting β.
The judgment basis for the end of sampling at the N th sampling point in the exploration sampling stage of potential feasible region is similar to that in the first stage. The main sampling and updating process is as follows: (a) firstly, an approximate surrogate model of objective function and constraint values is constructed based on the existing sample points; (b) according to Equation (7), the optimal solution is obtained, which means the coordinates of the new sampling points x 2n are obtained; (c) calculate the objective function and constraint value according to the coordinates of the new sampling points, the approximate model and the original model. The constraint value can be expressed as: ε(x 2n ). The sampling updating process is repeated until the new sampling point satisfies any of the following criteria, and the exploratory sampling phase of the potential feasible region is ended.
where, ε(x n ) is calculated according to Equation (9), ε 2 , λ 2 , I 2 are preset parameters in the exploratory sampling stage and related to the criteria at the end. The sampling points in the exploration and sampling stage of the potential feasible region are mainly distributed in the potential feasible region, and the new sampling points will keep a certain distance from the original sample points. The main sampling mode is to explore and sample the whole feasible region, which can greatly improve the approximate accuracy of the approximate surrogate model in the feasible region.
(5) Third sampling stage-potentially optimal sampling In the potentially optimal sampling stage, the penalty factors in Equation (7) are adjusted as α remains unchanged, β = 0.
The sampling process in this stage is the same as the main steps in the previous two stages. When the sampling points are updated for the N th time, the following steps are mainly carried out: (a) firstly, an approximate surrogate model of the objective function and constraint values is constructed based on the existing sample points. (b) According to Equation (7), the optimal solution is obtained, which means the coordinates of the new sampling points x 3n are obtained. (c) Calculate the objective function and constraint values according to the coordinates of the new sampling points, the approximate model, and the original model, and the constraint value can be expressed as ε(x 3n ). The process of sampling updating is repeated until the new sampling point satisfies any of the following criteria, and sampling stage of potential optimal sampling can be completed. In this algorithm, the end of convergence in the potential optimal sampling stage is the end of convergence of the whole algorithm. The criteria are as follows: where, ε x , ε 3 and λ 3 are the preset parameters required by the convergence criterion. The meaning of criterion I is that the current feasible solution is x N 3 , which satisfies the inequality and equality constraints. Criterion II means the Euclidean distance between the current feasible solution x N 3 and the recent λ 3 sampling points is less than or equal to ε x . Criterion III means the absolute value of error between the current sample points surrogate approximation model f (x N 3 ) and its original function f (x N 3 ) is less than or equal to ε 3 . The potential optimal sampling points are in the complete development mode, and the update of sampling points will converge quickly to the optimal point, which ensures the optimization convergence efficiency to a great extent.

The Proof of Global Optimal
Since the penalty function method has been used to transform the constrained optimization problem into an unconstrained optimization problem, the unconstrained optimization problem is considered in the proof of global convergence.
Define function: f (X), where X is an independent variable. Then the objective function can be defined as: min f (X) (13) Construct an auxiliary objective function as follows: where, when X is in the feasible area, P(X) = 0. When X is not in the feasible area, P(X) > 0. Then there is the following theorem: Let M k > 0 increase and tend to infinity, then the unconstrained problem can be described as Define the global optimal solution of the above optimization problem be X * k (k = 1, 2, . . .). If lim k→∞ X * k = X * , then X * is the global optimal solution of the unconstrained problem.
Proof. Because X * is the global optimal solution of unconstrained problem minF(X, Replace M k with M k+1 at the right end of Equation (15) to obtain: Adding Equations (15) and (16) and a transfer merge can get Combining Equations (15) and (19) gives: Then take any point X ∈ S (Attention: Equations (17) and (20) illustrate F(X * k , M k ) and f (X * k ) are incremental. In Equation (21), it is shown that they are upper bounded, so the limit lim k→∞ F(X * k , M k ) and the limit lim k→∞ f (X * k ) exist. Further notes that: So the limit lim k→∞ M k P(X * k ) exist, thus lim k→∞ P(X * k ) = 0. Then due to the continuance of P(X), we can get This can illustrate X * ∈ S. Pay more attention to both ends of Equation (21), ∀X ∈ S (Due to the continuance of f (X)) So X * is the global optimal solution of the unconstrained problem.
Note 1: Although the global convergence of this method has been proved in this section, considering the influence of interval selection of optimization variables, number of iteration steps, and accuracy of the approximate model, the optimization method based on surrogate model can only ensure the global optimization of probability.

1
This method realizes the adaptive adjustment of the distance between the new sampling points and the existing sample points, which can ensure that the distance constraint value is always within a reasonable range in the process of updating the sampling points. 2 This method makes each stage of sampling have a unified form and adjusts the exploration and development ability of the surrogate-based optimization method by changing the penalty factor, which can realize the sampling purpose of different stages. 3 Based on the approximate model of the objective function and the constraint value, the convergence end criterion of each sampling stage is established, which is conducive to terminating the sampling update after a reasonable number of samples and achieving the purpose of improving the optimization efficiency. 4 Decision science, convenient operation, clearly defined objective and process controllability are the salient characteristics of this method, which ensures the reliability and global optimality of the surrogate-based optimization method, with good convergence efficiency.

Improve Surrogate-Based Optimization Method
The flow chart of the surrogate-based optimization method proposed in this paper is as follows: The blue lines in Figure 1 mean that adaptive sampling will be carried out in each stage. The difference is that the values of penalty factors α and β in Formula (7) are different in different stages of adaptive sampling, which can make the exploration and development capabilities of different stages different.
The improved surrogate-based optimization method is mainly divided into the initialization stage, the approximation stage, the convergence judgment stage, and the adaptive sampling stage: (1) Initialization stage The design space is mapped to the unit hypercube space, and the optimized Latin hypercube method is used to determine the initial sample. The sample number is determined by the following equation: where, dim is the dimension of the design space, and in this paper dim= 9. The original model is used to calculate the objective function and constraint value of the sample, and the initial sample set can be obtained: where, x i represents the coordinate of the sampling points, y i represents the output response of the original model. Generally, it is a vector composed of the objective function and all constraints.
(2) The approximate stage The improved radial basis function is used to establish the approximate model, and the particle swarm optimization algorithm is used to solve the optimal solution of the approximate model. (3) The convergence criterion The convergence criterion is consistent with the above in the three stages of location sampling, exploration sampling, and potential optimal sampling in the potential feasible region.
(4) The updating stage of sampling points The adaptive infill strategy proposed [30] in this paper is as follows: The information of new sampling points can be obtained by solving the optimization problem of the above equation.

Mathematical Modeling and Verification of Overall Parameters
The LauncherOne air-launched rocket referred to in this paper adopts the hanging loading mode, and its carrier is a modified Boeing 747-400 large passenger aircraft. Considering the flight capability of the carrier aircraft, the initial launch altitude of the rocket is 10.67 km, the initial launch velocity is 262 m/s, and the initial launch velocity inclination is 25 degrees.
In this paper, according to the launching conditions and overall parameter of LauncherOne, the rocket (original proposal) is designed with a three-stage solid engine and a four-stage solid engine, and its main overall parameters are shown in the Table 1.  Table 1 mainly shows the relevant overall parameters of the original scheme, which provides support for the overall parameters optimal design below.
Rocket trajectory simulation is an important module to judge whether the optimization results meet the requirements. For the trajectory simulation of the air-launched rocket, the dynamic equation is the basis of trajectory simulation. As the mass of the rocket will change with the consumption of propellant and the falling off of the engine during flight, it is necessary to establish a mass and thrust model. As the falling off of the engine changes the aerodynamic reference area of the rocket and the change in air density, which will affect the aerodynamic force of the rocket, it is necessary to establish its aerodynamic model. Different flight timing affects the integration results of the rocket trajectory dynamic equation, so the flight timing model is also established. Various models are shown in the following sections:

3-DOF Trajectory Simulation Model
In the launch coordinate, the 3-DOF trajectory equation of the rocket [31] are as follows: where, α, β are angle of attack and angle of sideslip; θ, σ, υ are flight path angle, heading angle and bank angle; F is motor thrust; X, Y, Z are drag, lift force, and yawing force; x, y, zare the coordinate in the launching coordinate system; V is the speed of the rocket; mis the mass of the rocket; g is the gravitational acceleration; dt is the derivative of time.

Mass Calculation Model
The mass model of the rocket mainly includes load mass, equipment mass, structure mass, and propellant mass. The calculation of the take-off weight of the rocket is the basis of the trajectory simulation and the optimization of the overall parameter design. The derived calculation equation [32] of the takeoff mass of the rocket can be expressed as: In the above formula, m 0 is the takeoff mass of the rocket; m f is the mass of the rocket fairing; m b is the quality of instruments and equipment on the rocket; m pi is the mass of stage i motor propellant; m eni is the structural mass of stage i engine; m si is the mass of stage i motor interstate section (connecting mechanism, missile wing, tail fin, etc.).
The relation between the take-off mass of the i sub-stage rocket m 0(i) and the take-off mass of the i + 1 sub-stage rocket m 0(i+1) is The mass coefficient of the i stage rocket motor structure is denoted by α i (refers to the ratio of structural mass and charge quantity, α i = α en = m en /m p ), and it can be concluded that m pi + m eni = (1 + α i )m pi (31)

Thrust Calculation Model
As the overall scheme in this paper belongs to the demonstration stage of the rocket carrying capacity, the thrust model of the motor can be simplified to the average thrust. Set the average thrust of the i stage-level rocket motor as F i , the working time of thrust as t i , and specific impulse is defined as I sp , then the following relation can be obtained: That is, the thrust multiplied by the working time is equal to the total impulse. In the design scheme of this paper, considering the height of the motor ignition moment of each level, the average thrust of the primary motor is set as the ground average thrust, and the average thrust of the secondary and higher engines is the vacuum average thrust.

Aerodynamic Calculation Model
Like other objects, when the rocket moves relative to the atmosphere, the atmosphere will form forces on the surface of the rocket, among which the main forces are lift force and drag force, which are expressed as follows [31]: In the above formula, C D and C L are drag and lift coefficients respectively; ρ is atmospheric density; S M is the pneumatic reference area. The relationship between lift coefficient and drag coefficient is expressed as follows: In the above formula, C L 0 , C L α and K are constants, and their values are as follows: C D 0 is the zero-lift drag coefficient, and its value changes with the change of Mach number. The specific value can be referred to in the literature [33].
The first stage of the LauncherOne rocket is known to have a diameter of 1.6 m and the second stage is 1.3 m. Therefore, the reference area of the solid rocket designed in this paper is: S M 1 = 2.01 m 2 , S M 2 = 1.32 m 2 .

Flight Program Design
LauncherOne is an air-launched two-stage rocket. Its flight process is mainly divided into five stages, namely, the first powered flight stage, the first unpowered glide stage, the second powered flight stage, the second unpowered glide stage, and the second re-ignition (final boost) stage. The flight procedure is as follows: (1) The take-off stage of the carrier aircraft. After the final push stage was shut down, the load was sent to the preset orbit. According to the basic timing sequence and overall parameters of LauncherOne, the simulation of 500 km solar synchronous orbit entry was realized by adjusting the relevant parameters, and the carrying capacity was at least 300 kg. The simulation results are as shown in Figure 2. It can be seen from the above Figure 2 that the LauncherOne rocket designed can send an object with a total mass of 450 kg into a solar synchronous orbit of 500 km, which proved that the 3-DOF trajectory simulation module in this paper had certain reliability and provided support for the design optimization of the next step of solid rocket Based on the flight time sequence of LauncherOne, this paper designed the flight time sequence of an air-launched rocket based on a solid rocket engine, which adopted the way of direct entry into orbit. It was mainly divided into five stages, such as aircraft and rocket separation stage, first-stage flight, second-stage flight, third-stage flight, and final boost stage. The details are as follows [31]: (1) Aircraft and rocket separation stage Referring to the parameters information of Virgin air-launched rocket and carrier aircraft when separating, the launch parameters of the designed air-launched rocket in the 3-DOF trajectory simulation were consistent with the LauncherOne, that was, the launch altitude was 10.67 km, the speed was 262 m/s, and the launch inclination was 25 degrees.
(2) First stage flight After separation, the rocket flew at an altitude of nearly 11 km, fired at an inclination of 25 degrees, and then ignited the first stage motor and performed a programmed angleof-attack pull-up flight after transonic velocity, and maintained zero-angle-of-attack flight after reaching a certain altitude and velocity inclination. The flight procedure is as follows: In the above formula, t 1 is the time to start the pull-up of the angle of attack. Since the initial velocity of the rocket after separation had reached 0.9 Ma, it can be set as a few seconds. t 2 is The moment to pull up the end for the angle of attack; t 3 is the end time of the first gliding stage; t 4 is the end time of separation between the first and second stages. Write the gliding time of first stage as: t g1 , then there is t g1 + t 2 = t 3 . α m is the maximum absolute value of attack angle; t m is the time corresponding to the maximum angle of attack.

(3) Secondary and tertiary flight stage
The rocket is already flying in the rarefied atmosphere during the second stage, then the aerodynamic for pitch Angle program choice is negligible. In general, the simplified form of the secondary flight procedure is adopted in engineering applications. In this case, the linear expression of the pitch angle program in terms of time is where ϕ 0 is the pitching angle program of the starting point of the stage, which is equal to the pitching program angle of the termination point of the previous stage, .
ϕ is the change rate of the pitching program angle of this stage, and the initial value of the optimal pitching program angle and other trajectory parameters can be provided for the next flight stage by designing . ϕ and ∆t of the gliding stage according to ϕ 0 . In the optimization theory, the optimal solution obtained by the numerical method is very close to the linear relationship. Equation (27) applies not only to the second flight stage in the vacuum section but also to all flight stages above the second flight stage in the vacuum section.
The flight procedure [28] at this stage is designed as follows: In the above equation, ϕ(t 4 ) is the program pitch angle at the end of the first stage flight. The separation between the first stage and the second stage is zero angle of attack, ϕ(t 4 ) is mainly determined by the flight path angle θ. The optimal initial pitching angle of the second stage program can be obtained by designing the first stage gliding time t g1 . t 5 is the shutdown time of the second-stage engine; . ϕ 2 and . ϕ 2g are respectively the change rate of pitch angle in the second-stage power flight stage and the change rate of pitch angle in the second-stage glide stage. Set the second-stage gliding time as t g2 , thent g2 + t 5 = t 6 , t 6 is the end moment of the second-stage gliding; t 7 is the end time of fairing separation.

(4) Final boost stage
The program pitch angle of the final boost stage is designed as follows: In the above formula, t 8 is the shutdown time of the final stage engine, and . ϕ 3 is the change rate of the pitch angle of the final boost stage. The final stage motor shutdown time should make the rocket reach the conditions of orbit.
Take the three-stage rocket as an example, the flight process is as shown in Figure 3.  The program pitch angle of the final boost stage is designed as follows: In the above formula, 8 t is the shutdown time of the final stage engine, and 3 ϕ  is the change rate of the pitch angle of the final boost stage. The final stage motor shutdown time should make the rocket reach the conditions of orbit. Take the three-stage rocket as an example, the flight process is as shown in Figure 3.

Design Variables
Based on the rocket flight procedure described above, the design variables selected for rocket optimization are shown below [34]:

Design Variables
Based on the rocket flight procedure described above, the design variables selected for rocket optimization are shown below [34]: where, F i and ∆t i are respectively thrust and working time of each motor; A 0 is the launch azimuth; t g12 , t g23 and t g34 are respectively the sliding time of the working clearance between the first and second motors, the second and third motors, the third and fourth motors; ∆m e is the variable of the load mass, that is, the separation mass at the time of take-off; ∆t p and α p are respectively the duration of the attack angle and the value of attack angle. The range of optimization design variables are shown in the following Table 2.

Objective Function
The objective of the optimization algorithm for trajectory optimization is to seek the maximum payload mass whereas the takeoff mass remains constant. Therefore, this paper selects the maximum payload mass as the optimization objective function [28]. In the optimization process, the change of load mass at ground take-off time is taken as the objective function, that is, the minimum separation mass is taken as the optimization objective.

Constraint Conditions
In this paper, the optimization problem needs to meet inequality constraints and equality constraints. To ensure the safety of the rocket body and the launching requirements during the flight, the constraints of normal overload n y , dynamic pressure during the separation of the first and second engines q sep and Mach number Ma before the attack angle pulling up are constrained. The inequality constraints are as follows: where n ymax is the upper limit of normal overload, q max sep is the maximum dynamic pressure at the time of separation, Ma min is the minimum Mach number. Considering that the angle of attack should keep zero at transonic velocity, in this paper, n ymax = 0.05, q max sep = 10kpa and Ma min = 1.1.
According to a sun-synchronous circular orbit with a target orbit of 500km, the following equality constraints can be obtained: where, h, v, e, i are respectively orbital altitude, orbital velocity, orbital eccentricity, and orbital inclination; the target orbit parameter are h obj = 500km, v obj = 7612.57m/s, e obj = 0, i obj = 96.67 • .
Since the orbital inclination is mainly affected by the launch azimuth A 0 , it is easy to satisfy in the optimization process and take A 0 = 190deg. The constraints for height, velocity, and orbital eccentricity are complex, which are sensitive to more variables, and difficult to obtain the global optimal solution under the three constraints. To simplify the optimization problem, Ceq 1 (x), Ceq 2 (x), and Ceq 3 (x) are transformed to some other inequality constraints. Consider calculating the perigee height h p at the shutdown point of the final boost stage, the perigee height h p ≥ 500km is taken as the constraint. Since the objective function is the maximum load mass, the global optimal solution must be a circular orbit of 500 km. Thus, the constraint conditions of the air-launched rocket in this paper can be obtained as follows:

Optimization Results and Analysis
The trajectory of the air-launched rocket was optimized by using the surrogate-based optimization method. Firstly, the overall parameters of the solid rocket were designed and adjusted to meet the constraint indexes (defined as the original scheme), and the 3-DOF trajectory model was applied for simulation analysis. The overall parameters were optimized to maximize the carrying capacity, and 100 initial samples were selected. The results of the optimization for the three-stage air-launched and four-stage air-launched rockets are shown in the following Table 3. Its trajectory simulation is shown in the following Figure 4. The optimized convergence curves are as follows: The results show that the objective function of the three-stage rocket converges at 1750 iterations, and the optimal feasible solution is −79.35 kg at the 1730th time, (the separation mass at launch is −79.35 kg) which is equivalent to the payload mass of the rocket can be increased by 79.35 kg. Compared with the original scheme, the carrying capacity is increased by 26.45%, and all the indicators met the constraints. The objective function of the four-stage rocket converges at 1020 iterations, and the optimal feasible solution is −141.48 kg at 1012th iteration, which means the load capacity is improved by 141.48 kg. The load capacity is improved by 47.16% compared with the original scheme and all the indexes also meet the constraints.
From the above comparative results of air-launched three-stage and four-stage solid rockets, it is clear that the four-stage rocket has stronger launch capability, but generally has a longer flight time and a more complex structure. To further compare the air-launched rocket for the improvement of the launch capacity, the four-stage rocket is selected for the optimization of the maximum launch capacity for ground vertical launch, and the results are as shown in Table 4 and Figure 5.  The results show that the objective function of the three-stage rocket converg 1750 iterations, and the optimal feasible solution is -79.35 kg at the 1730th time, (the aration mass at launch is -79.35 kg) which is equivalent to the payload mass of the ro can be increased by 79.35 kg. Compared with the original scheme, the carrying capaci increased by 26.45%, and all the indicators met the constraints. The objective functio the four-stage rocket converges at 1020 iterations, and the optimal feasible solution 141.48 kg at 1012th iteration, which means the load capacity is improved by 141.48 kg. load capacity is improved by 47.16% compared with the original scheme and all th dexes also meet the constraints.
From the above comparative results of air-launched three-stage and four-stage s rockets, it is clear that the four-stage rocket has stronger launch capability, but gene has a longer flight time and a more complex structure. To further compare the launched rocket for the improvement of the launch capacity, the four-stage rocket i lected for the optimization of the maximum launch capacity for ground vertical lau  The optimized convergence curves and the height, velocity, local flight path angle, and thrust curves are as follows:  Figure 5 show that the objective function converges after 1279 iterations, and the optimal solution is 13.71 kg after the 1028 iteration, which is equivalent to a 13.71 kg reduction in carrying capacity. Compared with the original scheme, the carrying capacity is 4.57% lower. According to the four-stage rocket model proposed in this paper, the carrying capacity of the air-launched rocket is about 51.49% higher than that of the ground-launched rocket. It further explained the superiority of an air-launched rocket.

Conclusions
The surrogate-based optimization algorithm proposed in this paper has the characteristics of high global optimization efficiency and high precision. The precision of the surrogate model is significantly improved by the adaptive sampling point updating strategy and the kernel width estimation of the radial basis function, which balances the sampling explo-ration and development ability of the optimization method. Through the establishment of the 3-DOF trajectory simulation model, mass calculation model, thrust calculation model, aerodynamic calculation model, and flight program design to determine the overall parameter design simulation model, and through the replication of LauncherOne rocket to verify the correctness of the design model, lay a foundation for the next design optimization.
Based on the proposed method, the overall parameter design and optimization problem of the air-launched rocket is illustrated, and the load capacity is optimized with the amount of change in payload mass (separation mass) as the objective function and analyzed by 3-DOF trajectory simulations. Based on the overall parameters and basic flight time sequence of the LauncherOne air-launched rocket, the overall parameters of the air-launched solid rocket in this paper were designed, and the original scheme was formed through 3-DOF trajectory simulation. Then the optimization calculation was carried out based on the original scheme by designing the optimization variables, determining the optimization objective function, and defining the constraint conditions.
Compared with the traditional optimization algorithm, which calls the original model thousands of times, this method greatly improves the optimization efficiency and optimization precision. From the optimization results, the optimal feasible solution of the air-launched rocket increases the capacity respectively by 26.45% and 47.16%, compared with the original solution, which demonstrates the reliability of the optimization method in this paper. For the four-stage rocket, compared with the ground launch, the carrying capacity of the optimal solution is improved by 51.49%. This further illustrates the advantages of air-launched rockets.
In summary, the surrogate-based model optimization algorithm proposed in this paper solves the main difficulties of multidisciplinary optimization algorithms. The method was applied to the conceptual design and parameter optimization of the air-launched rocket system. The results show that the method is effective and reliable for the optimal design of the launch vehicle, which can improve the launch capacity and reduce the amount of calculation under the same take-off mass. In addition, the calculation results of the fourstage rocket are compared. The results show that the air-launched rocket can significantly improve the carrying capacity and indirectly reduce the economic cost of the launch mission. Therefore, combined with the financial cost and reliability model, this research work can find the ideal overall parameters for the air-launched rocket system, and provide support for the development of low financial cost, fast response, and high-reliability satellite launching device. In general, it has a certain reference significance for the development of the overall parameter optimization of air-launched rockets.
Author Contributions: S.C. is responsible for model construction and simulation verification, J.L. is responsible for learning and developing optimization theory, S.Z. and X.B. are responsible for grasping the structure and content of the article and providing technical support, and D.S. is responsible for the final typesetting and revision of the article. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All data, models, and code generated or used during the study appear in the submitted article.