Surrogate-Based Optimization Design for Air-Launched Vehicle Using Iterative Terminal Guidance

: In recent years, the penetration of low-cost air-launched vehicles for nano/micro satellites has signiﬁcantly increased worldwide. Conceptual design and overall parameters optimization of the air-launched vehicle has become an exigent task. In the present research, a modiﬁed surrogate-based sequential approximate optimization (SAO) framework with multidisciplinary simulation is proposed for overall design and parameters optimization of a solid air-launched vehicle system. In order to reduce the large computation costs of time-consuming simulation, a local density-based radial basis function is applied to build the surrogate model. In addition, an improved particle swarm algorithm with adaptive control parameters is proposed to ensure the efﬁciency and reliability of the optimization method. According to the LauncherOne air-launched vehicle, the overall optimization design problem aims to improve payload capacity with the same lift-off mass. Reasonable constraints are imposed to ensure the orbit injection accuracy and stability of the launch vehicle. The inﬂuences of the vehicle conﬁguration, optimization method, and terminal guidance are considered and compared for eight different cases. Finally, the effect on the speed of optimization convergence of employing a terminal guidance module is investigated. The payload capability of the optimized conﬁgurations increased by 27.52% and 23.35%, respectively. The ﬁnal estimated results and analysis show the signiﬁcant efﬁciency of the proposed method. These results emphasize the ability of SAO to optimize the parameters of an air-launched vehicle at a lower computation cost.


Introduction
As an emerging low-cost, quick-response, high-mobility and flexible launch system for micro/nano-satellites, air-launched vehicles are drawing worldwide attention with the rapidly increasing need for micro/nano-satellite launches. In 1990, the American corporation "Orbital Sciences" successfully launched a solid launch vehicle named "Pegasus" for the first time using a B-52, which was the world's first air-launched vehicle to make a successful orbit, and has now been successfully launched more than forty times. Currently, a Russian aerospace company is developing a two-stage air-launched vehicle named "Polyot" and, Boeing from the United States is developing an air-launched vehicle as well. In 2021, after a failed launch attempt, Virgin Orbit finally launched the "LauncherOne" liquid launch vehicle, using a Boeing 747 carrier aircraft which had a 300 kg load capacity for sun-synchronous orbit (SSO) and a 500 kg load capacity for low-earth orbit (LEO) [1]. Compared with traditional ground-launched vehicles with the same load capacity, the total mass of air-launched vehicles can be reduced by 20%~30% and the total requirement velocity can be reduced by 10%~15% [2]. The principal advantage of air-launched vehicles is that there is no need to fly through the low, dense atmosphere, the drag and dynamic pressure of which require a considerable amount of extra configuration design and mass of propellant. In addition, air-launched vehicles can utilize the initial velocity of the air launch platform (such as airplane) and have no requirement for a launching site, which directly reduces the launch cost and preparation time. On the other hand, the initial mass of the launch vehicle is limited by the payload capability of the air launch platform, which makes the mass of payload a significant indicator. Thus, it is necessary to optimize the overall parameters of the launch vehicle's configuration to obtain the maximum payload with the limited initial mass. The overall design of air-launched vehicles is a typical multidisciplinary coupling problem. The complex coupling relationship among disciplines makes it difficult for traditional individual disciplinary design optimization methods to find the optimal solution. A high-efficiency and reliable optimization method for overall design could further develop the advantages of air-launched vehicle system and potentially shorten the iteration periods of launch missions. Existing studies have of air-launched vehicles have mainly focused on their advantages compared to other options [3], system composition and analysis [4,5], propulsion system design [6,7], separation analysis [8], and multidisciplinary design and optimization studies [9][10][11]. Meanwhile, studies on the overall parameter design of air-launch vehicles using modern intelligent optimization methods are lacking.
In this study, an overall parameter optimization design for an air-launched vehicle system is based on the coupling of RBF surrogate-based SAO algorithms and multidisciplinary simulations, including aerodynamics, propulsion, trajectory, guidance, etc. The multidisciplinary simulation of air-launched vehicles could potentially decrease the huge computational costs. The key point in enhancing optimization efficiency is reducing the calling times of the original multidisciplinary model as far as possible. The advanced optimization technique is the main approach to solving problems in engineering design and application, and has been wildly studied [12][13][14][15][16][17]. Sequential approximation optimization method is known for its high efficiency, robustness [18], and generality, as it requires far fewer original model evaluations to locate the global optimum compared with evolutionary algorithms (EAs) such as genetic algorithms (GA) [19], simulated annealing (SA) [20], particle swarm optimization (PSO) [21], immune algorithms (IA) [22], and artificial bee colonies (ABC) [23]. In addition, the SAO method is specifically adaptive for overall parameter optimization based on coupling of high computational cost models including aerodynamic analysis, internal ballistic analysis, flight trajectory, guidance simulation, etc. The key step in the process of SAO is surrogate modeling, which is constructed and sequentially updated by infilling new sampling points until the results meet the terminal criterion [24]. The surrogate modeling technique is essentially a multivariate function approximation method, different varieties of which have been widely studied [25,26], including Kriging [27][28][29], Radial basis Function (RBF) [30][31][32], Support Vector Regression (SVR) [33][34][35], and Response Surface Methodology (RSM) [36][37][38]. The RBF model was originally proposed to fit the irregular topographic contours of geographical data [39], which is extensively used in SAO algorithms for its reliability in terms of accuracy and robustness [40]. As the parameters (centers, widths, and weights) of RBF models have conclusive impacts on its accuracy, several methods have been proposed to optimize the parameters of RBF for uniform samples, nonuniform samples, and sufficient samples. Recently, Han has proposed a self-organizing RBF model using an adaptive gradient multi-objective PSO (AGMPSO) method to balance the accuracy and complexity of RBF models [41]. Zhang has proposed a two-stage space division-based optimization method to optimize the width parameters of RBF models for engineering problems [42]. Li has proposed an immune algorithm system-based RBF (IAS-RBF) method to overcome the premature convergence problem [43]. Slema has proposed a model reference adaptive controller for nonlinear systems based on RBF using the error between the system response and the desired response [44].
In order to improve the optimization efficiency and performance while searching for the global optimal solution, an RBF model based on the local density of sampling points is proposed here to act as a surrogate in the time-consuming multidisciplinary model of the airlaunched vehicle. Specifically, a novel adaptive control parameter method is developed to enhance the global search capability of PSO in solving the surrogate optimization problems.
In addition, a multistage adaptive infilling strategy along with the penalty method is employed to balance the exploration and exploitation of the SAO. Compared with the conventional optimization algorithm, which calls the original model thousands of times, the main merit of SAO is that it notably reduces computation cost. The efficiency of the proposed method is evaluated in the present study, with a significant improvement in the payload capability of the air-launched vehicle system being attained.
The organization of this paper is as follows: Section 2 establishes a multidisciplinary conceptual design framework of the air-launched vehicle system. Section 3 illustrates the optimization problem, including the design variables, objective function, and constraints. Section 4 elaborates the procedures of the proposed SAO method, which uses a local density-based RBF surrogate model, and the PSO method, which uses adaptive control parameters. Section 5 analyzes the optimization results and flight performance of the studied cases. Finally, a brief conclusion of this study is provided in Section 6.

Multidisciplinary Optimization Framework and Overall Parameters Modules
The LauncherOne air-launched vehicle is carried by a modified Boeing 747-400 carrier aircraft. Considering the flight capability of the transport aircraft, the initial launch altitude of the launch vehicle is 10.67 km, the initial launch velocity is 262 m/s, and the initial velocity inclination is 25 degrees.
In this study, according to the launch conditions and overall parameters of LauncherOne, the launch vehicle is designed in both three-stage and four-stage configurations with solid motors; its main overall parameters are shown in Table 1. The overall design and optimization studies of the air-lunched vehicle system are challenging tasks. The complex interactions among the design disciplines (aerodynamics, propulsion, structure, trajectory, etc.) make the design problem difficult to handle from the viewpoints of optimization methods and computational cost. The optimization framework and a brief illustration of the related multidisciplinary modules used in this study are shown in Figure 1.

Aerodynamic Calculation Model
Although the air-launched vehicle launches at a high altitude, which provides a benefit in terms of drag reduction (about 30% in 10 km altitude) and propellant savings, the effect of aerodynamics remains non-negligible. In this study, the main forces, including the lift force, L, and drag force, D, can be expressed as follows [45]: where C L and C D represent the coefficients of drag force and lift force, ρ represents the atmospheric density, v represents the velocity of the launch vehicle, and S M represents the pneumatic reference area. The relationship between C L and C D can be expressed as where C L 0 , C L α , and K are constants, and can be determined as follows: where C D 0 represents the zero-lift drag coefficient, the value of which changes with the Mach number, Ma. The specific value can be found in the literature [46].
The first stage of the LauncherOne launch vehicle 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 launch vehicle designed in this paper is S M 1 = 2.01 m 2 , S M 2 = 1.32 m 2 .

Mass Estimate Module
The mass estimate module is a basic module of the multidisciplinary optimization framework. It considers a series of equations for the air-launched vehicle, including the mass of payload, structure, propellant, apparatus, etc. The gross mass of the air-launched vehicle can be expressed as follows [45]: where m 0 is the gross mass of the launch vehicle, m pd is the mass of the payload, m a f t is the mass of the afterskirt and wings of the launch vehicle, m f is the mass of the fairing, m c is the mass of apparatus cabin, including payload adapter, orbital maneuvering system, guidance set, etc., n is the number of the stages, m pi is the propellant mass of the i th stage motor, m si is the structure mass of the i th stage, and m ai is the apparatus mass of the i th stage of the launch vehicle, including the mass of the inter-stage structure, raceway, cable assembly, etc.
The mass ratio of the i th stage motor, α i , is defined as The relationship between the gross mass of the i th sub-stage launch vehicle, m 0i , and the gross mass of the i + 1 th sub-stage launch vehicle, m 0(i+1) , is

Propulsion Module
Simplifying the thrust of the solid rocket motor (SRM) in each stage to average the thrust is a simple and reliable method of thrust characteristic calculation that can be used effectively during the overall design optimization phase. The average thrust can be determined by the working time of the motor, the specific impulse, and the propellant mass of the SRM: where P i represents the average thrust of the i th stage SRM, t i represents the working time of the i th stage SRM, and I sp represents the specific impulse of the i th stage SRM. Therein, the average thrust of the first stage SRM is set as the 10 km altitude average thrust, and the average thrust of the second and subsequent stage SRMs is set as the vacuum average thrust.

Trajectory Module
A three-degrees-of-freedom (3-DOF) trajectory simulation using the direct orbital injection method is implemented and tested in the trajectory module. In this module, the state parameters, including velocity, flight angle, position, forces, and mass, are computed by integrating the equations of motion in the launch coordinate [45]. The basic 3-DOF trajectory equations of the launch vehicle are as follows: where α is the angle of attack, β is the angle of sideslip, θ is the flight path angle, σ is the heading angle, υ is the bank angle, p is the motor thrust, X is the drag force, Y is the lift force, Z is the yaw force, and x, y, z are the positional coordinates in the launch coordinate system. The flight process of the air-launched vehicle system is shown in Figure 2 (a three-stage configuration is shown in this example), which can include several phases: take-off, separation and ignition, powered flight and unpowered glide phases in each stage, terminal guidance, and payload injection.

Terminal Guidance Module
As the altitude of the objective orbit is above 80 km, an autonomous terminal guidance method for direct orbit injection outside the atmosphere is proposed in order to meet the for accuracy and reliability requirements of the launch mission and as well as to accelerate the convergence of the objective function in the optimization process. The guidance process of the air-launched vehicle system is shown below. Step 1-Launch Initialization: The velocity (V * A f ), geocentric distance (r * A f ), flight path angle (γ * ), and orbit inclination (i * ) of the objective orbit injection point should be calculated as the terminal constraints. In addition, the initial parameters include the average vacuum thrust (T vac ), second consumption (dm), geographic latitude of the launch point (B 0 ), launch azimuth (A 0 ), initial estimate time of the terminal guidance phase (t f ), and initial estimate valve of the costate variable ( → λ 0 ). As for the orbital injection problem, the main constraints are the injection point constraints, which include the position vector and velocity vector of the injection point. According to the optimal control theory, the optimal control problem is an unsolvable movable end-point problem, and the position vector and the velocity vector of the injection point are taken as the constraints. As there is no requirement for the right ascension of the ascending node of the objective orbit, a four-constraint problem is proposed here; the constraints are shown below: where → r A f is the dimensionless terminal position vector of the launch vehicle relative to the geocenter, → V A f is the dimensionless terminal velocity vector of the launch vehicle in the launch inertial coordinate system, r * A f is the geocentric distance of the injection point, V * A f is the velocity of the injection point in the launch inertial coordinate system, i * is the orbital inclination angle of the objective orbit, γ * is the flight path angle of the injection point, and → I N A is the unit vector of the north pole, which can be expressed in the launch inertial coordinate system as Step 2-Iteration of State Parameters: In each period of the terminal guidance phase, the current velocity vector, → V A0 , and geocentric position vector, → r A0 , of the launch vehicle are updated and nondimensionalized in the launch inertial coordinate system. The satisfied costate variables, → λ 0 , and nondimensional initial time, τ f , can be obtained by iteration using the Newton down-hill method: the initial value of the velocity vector, → V A , in the launch inertial system, and (12) where → I b is the unit vector of the body axis and g r0 = µ/r 2 0 , where µ is the gravitational constant. Following iteration, the optimal body axis, I * b (t), in the launch inertial system can be obtained and resolved in order to obtain the command attitude angles.
Step 3-Termination Criterion: The guidance process is terminated when the remaining flight time to the predetermined injection point obtained by iteration in the previous step is less than or equal to 5 s. The launch vehicle maintains the current attitude angle for the remaining flight time and then shuts down. On the contrary, if the remaining flight time is more than 5 s, the guidance process continues, then proceeds to Step 4.
Step 4-Next Guidance Period: The costate variables, → λ 0 , and dimensionless time, τ f , obtained during iteration can be regarded as the initial values of the next guidance period; then, Steps 2-3 are repeated for the next guidance period when it reaches the next guidance point.

Design Variables
According to the multidisciplinary optimization framework of air-launched vehicles, the design variables and their boundaries are shown in Table 2.

Objective Function
In order to determine the design direction and emphasis, it is necessary to obtain a reasonable objective function. As a key index in the launch vehicle design process, payload capacity partly represents the performance of the launch vehicle. Thus, it has generally been considered as the design objective in the optimization process. In this study, the objective of the optimization method for overall parameters design is to seek the maximum payload capacity with the same gross mass as the air-launched vehicle. The objective function can be expressed as

Constraints
The constraints of the optimization problem consist of inequality constraints and equality constraints. Considering the safety and stability of the launch vehicle and the launch requirements during flight, the vertical overload (n y ), dynamic pressure (q sep ) at the separation of the first stage motor, and Mach number (Ma) at the end time of the control method are constrained. Thus, the inequality constraints, g i (x), are shown below: where n ymax represents the upper limit of normal overload, q max sep represents the maximum dynamic pressure at the time of separation, and Ma min represents the minimum Mach number. Considering that the angle of attack should remain zero at transonic velocity, in this study, n ymax = 0.05, q max sep = 10 kPa, and Ma min = 1.1. As the objective orbit is a 500-km sun-synchronous orbit, the equality constraints, Ceq i (x), can be expressed as where h, v, e, and i are the altitude, velocity, eccentricity, and inclination of the orbital injection point, respectively, and h obj , v obj , e obj , and i obj are the altitude, velocity, eccentricity, and inclination of the objective orbit, respectively. In this study, h obj = 500 km, v obj = 7612.57 m/s, e obj = 0, and i obj = 96.67 • .
As the equality constraints are time-consuming and extremely difficult to satisfy in the optimization process, it is necessary, efficient, and simple to transform the equality constraints to inequality constraints. As the orbital inclination, i, is mainly affected by the launch azimuth, A 0 , it is easy to meet Ceq 4 (x) in the optimization process by controlling the launch azimuth A 0 in 190 degrees, while Ceq 1 (x), Ceq 2 (x), and Ceq 3 (x) are the constraints for height, velocity, and eccentricity, which are complex and sensitive to other design variables. In order to simplify this problem, an energy-based method is proposed during the optimization process. As the objective function is to search for the maximum payload capacity with the same gross mass of the launch vehicle, a larger payload capacity would lead to the launch vehicle being able to reach a lower orbital altitude. Thus, a constraint on the lower limit value of the perigee, h p , of the injection point can fully represent the effect of Ceq 1 (x), Ceq 2 (x), and Ceq 3 (x). Finally, the global optimal solution is supposed to be a circular orbit of 500 km; thus, the constraint conditions of the air-launched vehicle can be obtained as follows:

Surrogate-Based SAO Method
In this study, the overall parameters of air-launched vehicles are optimized for a practical launching mission based on the coupling of numerical optimization algorithms and multidisciplinary simulations. As the computation cost of multidisciplinary simulations can be potentially very high, the key point in accomplishing successful overall parameters optimization is to increase the efficiency of the optimization methods while reducing the calling times of multidisciplinary simulations as minimal as possible during the optimization procedure. SAOs are known for their lower computational costs, generality, robustness, and accuracy [47]. SAO algorithms require much lower times in the evaluation of original models to locate the global optimum when compared to evolutionary algorithms, such as genetic algorithms, simulated annealing, particle swarm optimization algorithms, immune algorithms, artificial bee colony algorithms, etc. Objectively, SAO algorithms are a particularly good fit for overall parameters optimization based on coupling with multidisciplinary simulations.
The surrogate model is the basis of the SAO method, the essence of which is multivariate function approximation aiming to satisfy certain accuracy conditions for optimal design and analysis, instead of high precision numerical simulation. In the process of SAO, a reasonable surrogate model can effectively reduce the calling times of high-precision simulation models and improve optimization design efficiency. Among the existing surrogate models, the RBF surrogate model is widely used in SAO thanks to its simple principle and ease of use. However, in practical applications a good surrogate model with good performance needs to be obtained by optimizing the shape parameters of the basis function, and the complexity of the shape parameter optimization increases sharply with the increase in the number of sample points; therefore, the reasonable and efficient determination of the shape parameters of the basis function represents a difficult research area in the field of RBF surrogate models.
In this section, a shape parameter determination method is proposed for the radial basis function based on the local density of sampling points, which unifies the determination of multiple shape parameters into one parameter and solves the problem of the complexity of shape parameter optimization increasing sharply with the sample size. In addition, a modified PSO using adaptive inertial weight is proposed to search the optimal shape parameter and improve the efficiency of constructing the surrogate model.

Local Density-Based RBF Shape Parameter Determination Method
Given a set of training samples [x i , y i ], i = 1, 2, · · · , n, the output prediction of the RBF surrogate model at any design variable can be expressed aŝ where x is the design variable, ϕ i (r) is the basis function of the i th sample, and ω i is the coefficient of the basis function. In this study, the Gaussian basis function is applied as the basis function, which can be expressed as where c i is the shape parameter of the i th sample. As for any non-zero shape parameter c i , after substituting the training samples into Equation (17), the basis function coefficients, ω i , can be obtained by introducing the interpolation condition,f (x i ) = y i where Φ ij = ϕ j (x i ) . According to the basic theory of RBF interpolation, for a training sample set with different sampling points, the corresponding basis function coefficients can be obtained by solving linear equations for any shape parameter c in Equation (19), after which the surrogate model can be constructed. In practice, different shape parameters have important impacts on the prediction accuracy at non-sampling points, which is shown in Figure 3. The influence of the shape parameters is mainly reflected in the following two aspects: (1) The global prediction accuracy of the RBF surrogate model first increases as the shape parameter increases from near zero, then the surrogate model oscillates and decreases in prediction accuracy as the shape parameter increases.
(2) There is no uniform criterion for when the appropriate shape parameter appears, and it depends on the property of the original output function.
Therefore, a reasonable determination method for the shape parameters of the basis function can effectively improve the global prediction accuracy of the surrogate model, which in turn reduces the number of training samples and improves the efficiency of optimization design. The basic idea of the surrogate model is to predict the output law of the whole design space through the output of training samples in the design space.
Thus, the output response of any sampling point is jointly influenced by the response of all sample points, each sample point has the dominant influence on itself, and the influence gradually decays on the region further away from the point. Furthermore, in the same form of the basis function, the shape parameter determines the rate of the decay of influence to each sample point in the surrounding space. A larger shape parameter means that the influence of the point can be propagated further in the entire output, while a smaller shape parameter determines a smaller influence range. Figure 4 shows the influence of the distance between sampling points and different shape parameters on the shape of the basis function. According to the distribution of sample points in the space, the influence range of each sample point should be appropriately small for the region with dense distribution of sample points to prevent overfitting. For a region with sparse distribution of sample points, the input parameters are provided arbitrarily given as there is no more information around the point to predict the output, and the output of the point can only be predicted by increasing the influence of the surrounding sample points, which means appropriately increasing the value of the shape parameter. Based on the above analysis, a determination method for the shape parameter is proposed based on the local density of sample points.
In this method, the local density, ρ(x), can be expressed as where c i represents the influence extent of the i th sample point on the contribution of local density. Too large or too small a value of c i will cause the local density to lose the local characteristics. Without loss of generality, in the unit hypercube design space the value of c i can be expressed as where n is the number of sample points and d is the number or dimension of design variables. The local density of the one-dimensional problem in Equation (20) can be obtained according to Equation (21). Figure 5 shows an illustrative example of density in one dimension. The example shows that the function has larger values in the region with a dense distribution of sampling points and smaller values in the region with a sparse distribution of sampling points. This method can reasonably quantify the local density of discrete sample points. As the shape parameter can be considered as the influence extent, the influence extent should be inversely proportional to the local density of each sample point, which can be expressed as follows: where V i represents the influence volume of the i th sample point and ρ i represents the local density at the i th sample point. The reference value, c i , of the shape parameter can be obtained by normalizing the total volume of influence, which can be expressed as Because the values of the optimal shape parameters vary with the model characteristics, the reference value of shape parameters should be scaled based on the same relative size of each shape parameter; the values of each shape parameter can be obtained as follows: where λ is the scale coefficient, which determines the effective range of the surrogate model. Thus, as the shape parameters are simplified into one parameter, λ, the RBF surrogate model can be efficiently constructed.

Particle Swarm Optimization Method Using Adaptive Control Parameters
As the RBF surrogate model of the air-launched vehicle system is sequentially established and updated, the computation cost is gradually reduced. The optimization problem with a surrogate model can be expressed as where f (x), g i (x), h j (x) are the surrogate models of the objective function, inequality constraints, and equality constraint, respectively, d(x) = min( (x − x i ) T (x − x i )) (i = 1, 2, · · · , n) is the minimum Euclidean distance between the new sampling point and former sampling points, and δ is the minimum distance between the existing sampling points, which gradually decreases during the SAO iteration process.
In order to solve this optimization problem, a modified PSO method is proposed in this section. Standard PSO methods cannot adaptively balance the performance of global search (exploration) and local search (development). It is easy to fall into local optima and low search accuracy. This paper improves the standard PSO from the perspective of control parameters.
The inertia weight, ω, and learning factor, c 1 and c 2 , are independent control parameters in the PSO method. Therein, the inertia weight, ω, is a key parameter to balance the global and local search capabilities of the particles. Larger inertia weight can make particles have better exploration capabilities, while smaller inertia weight can make particles have better development capabilities. In general, the particle swarm needs better exploration ability at the early stage and better development ability at the later stage in order to improve search efficiency. Thus, an adaptive inertia weight strategy should be applied in the optimization process. Existing strategies need to set the maximum number of optimization steps in advance, and the inertia weight is related to the maximum number. In this study, an adaptive inertia weight strategy is proposed based on the distribution of the current optimal solution, p best , found by a particle P i .
In this method, the design space is first normalized as follows: where X U i and X L i are the upper and lower bounds of the i th dimension of the original design space and n is the dimension of the design space.
The distribution of p best can be characterized by its standard deviation, σ pbest . At the beginning of the search, p best is approximately considered as a uniform distribution in the design space, of which the standard deviation, σ pbest , is 0.2887. Along with optimization process, σ pbest continues to decrease, and approaches 0 when the optimization process is convergent. Thus, the inertia weight can be expressed as where i is the current iteration number and ω max and ω min are the initial value and final value of the inertia weight, ω.
In addition, the learning factor is an important factor in balancing the exploration and development capabilities of PSO. In the early stage of optimization, the learning factor should have greater self-learning ability and weaker swarm learning ability to strengthen the global search. In the later stage of optimization, the learning factor should have greater swarm learning ability and weaker self-learning ability to accelerate the convergence to the global optimal solution. In [48], a learning factor strategy is proposed based on sine and cosine functions, although it relies on information about the maximum iteration number: where Iter _max is the maximum iteration number, i is the current iteration number, and ρ = 2, δ = 0.5. In this study, the modified learning factor strategy is based on the information of the standard deviation, σ pbest , which is shown below: where i is the current iteration number and ρ = 2, δ = 0.5.

Modified Sequence Approximate Optimization Method
In this study, the modified surrogate-based SAO method mainly consists of a design of experiment (DoE) stage, an approximation stage, and a termination criteria and infilling stage. The procedure of the proposed surrogate-based SAO method is shown in Figure 6.

DoE Stage
In this stage, the design space of the optimization problem is mapped to the unit hypercube space and the optimal Latin hypercube design (OLHD) method is applied to generate the initial sample points in the unit hypercube. The original multidisciplinary model is employed to calculate the objective function and constraint value of the initial sample points. According to the OLHD method, the number of initial sampling points is supposed to be the multiple of the dimension of the design variables. In addition, as the computation cost of the original multidisciplinary model is substantial, the calling times of the multidisciplinary model should be limited. Thus, the number of initial sampling points is set to 28. The initial sample set, including the initial samples and corresponding responses, are then obtained and applied in order to establish the surrogate models in the Approximation Stage.

Approximation Stage
In this stage, the proposed local-density-based method is applied based on the initial sample set to determine the influence area of each sample point and calculate the shape parameters of RBF. Then, the RBF surrogate models of the objective function and constraints are constructed and sequentially updated when new sample points are infilled using the infilling strategy in the Infilling Stage.

Termination Criteria
Termination of the proposed SAO is determined by the following criteria: (1) If the relative distance between the optimal solutions of two successive iterations is less than 1‰, then evaluate criterion (2). Otherwise, the SAO proceeds to the Infilling Stage.
(2) If the relative error between the two successive iterations' optimal objective functions with the constraints imposed using a penalty method is less than 0.0001, then evaluate criterion (3). Otherwise, the SAO proceeds to the Infilling Stage.
(3) If the relative error between the objective functions of the surrogate model and the true model is less than 0.0001, convergence is reached and the proposed SAO algorithm is terminated. Otherwise, the SAO proceeds to the Infilling Stage.

Infilling Stage
In order to reasonably infill new sample points in the design space, a multistage adaptive infilling strategy is employed to balance the exploitation and exploration capacity of SAO [49]. In this stage, the minimum distance, δ, between the sampling points in the sample set is initially calculated; then, based on the surrogate models constructed in the approximation stage, an adaptive infilling with penalty method is implemented by solving the optimization problem using the proposed PSO algorithm.
The optimization problem (25) can be converted to an unconstrained optimization problem with penalties as follows: where | f | max , |g i | max and h j max are the roughly maximum absolute values of the objective function, inequality, and equality constraints, respectively, and α i , β j , and τ are the penalty coefficients. The orders of magnitude for α i and β j are 10 4 , τ is set to different values in order to achieve sampling goals for different stages, and the global minimum point, x s f , of Formula (31) is solved using the proposed PSO algorithm with adaptive control parameters.
Then, x s f is added to the sample set in order to update the surrogate models in the next iteration. The termination conditions of the sampling stages are determined based on the approximation accuracy of the response surfaces. The details of the multistage infilling strategy are elaborated as follows: (1) Potentially feasible region locating stage (Stage 1) In this stage, τ is two orders greater than α i and β j , and reaches the magnitude order as 10 6 . The global minimum point, x s f , lies in the sparse regions, which clearly improves the global accuracy of the response surfaces. On the other hand, x s f inclines towards the feasible regions, which is beneficial for discovery of potentially feasible regions separate from each other. This sampling stage terminates automatically when one of the following two conditions is met: (a) the number of sampling iteration steps in Stage 1 reaches a presupposed parameter I 1 ; or (b) the relative error between the response surfaces and the original models satisfies the following inequation for the last successive λ 1 sampling iterations: where ε is the maximum relative error of the response surfaces of the objective function and constraints at the latest sampling point of each iteration and ε 1 is the maximum permissible error of Stage 1.
(2) Exploring in potentially feasible regions stage (Stage 2) In this stage, τ is set at the same magnitude order as α i and β j . The global minimum point, x s f , obtained from Formula (31) lies in the potentially feasible regions, and maintains distance with the existing sampling points, which clearly improves the response surfaces accuracy in the feasible regions. Sampling stage 2 terminates similarly to stage 1, and the parameters I 1 , λ 1 , and ε 1 are replaced with I 2 , λ 2 , and ε 2 , respectively. In this stage, τ is set as 0, and the optimum point of the response surfaces is sampled. Termination of Stage 3, as well as the whole SAO algorithm, occurs according to the termination criteria.
Finally, the optimal solution, together with its objective function and the constraints of the original model, are added to the sample set in order to update the surrogate models in the next iteration.

Numerical Examples
In order to provide a comparison with other optimization methods in the literature, five classic benchmark functions with continuous design variables are utilized here, namely, the well-known surrogate-based global optimization algorithms SOCE [50], EGO [51], HAM [52], and KMS [50], which have been widely cited. Table 3 shows the five benchmark functions. Table 3. Benchmark functions.

Benchmark Function Number of Design Variables
Design Space Global Optimum The optimization results of the benchmark functions are shown in Table 4, which states the best values obtained by all SBO algorithms; "Best, Median, Worst" represents the best, median and worst results obtained in ten tests within 300 function evaluations (NFE), while "MNFE" represents the mean NFE of the ten results. In Table 4, all five surrogate-based optimization algorithms perform well in lowdimension functions. However, the KMS, EGO, and HAM algorithms have problems with convergence in different cases and cannot solve high-dimension problems well. The proposed SAO and SOCE algorithms show similar performance in both accuracy and function evaluations, while the proposed SAO can achieve better results with fewer evaluations. Consequently, the proposed SAO method shows high robustness on the test functions and better performance in searching for the true global optimum, and proves competitive with the other surrogate-based optimization algorithms.

Results and Discussion
In this study, eight design cases are set up in order to explore the effect on speeding up the optimization convergence rate by employing the terminal guidance module in a conceptual design process of air-launched vehicles, exploring the efficiency between genetic algorithm (PSO) and surrogate-based methods (SAO). The specific illustration of the design cases is shown in Table 5. Based on the multidisciplinary optimization design framework, the optimal design results are shown in Table 6. The optimal objective function (payload capacity) and certain constraints (orbital perigee and apogee altitude) of each case are shown in Figures 7 and 8.  For the objective function, the results (Case 1|2 and Case 5|6) show that employing the terminal guidance module in the multidisciplinary framework can significantly improve the global search ability of classic evolutionary algorithm PSO; for the surrogate-based SAO method, the results (Case 3|4 and Case 7|8) show less effect from the terminal guidance module, as the construction of surrogate model mainly relies on the input and response of the original model without considering the intrinsic connection of the multidisciplinary framework, which means the surrogate-based SAO method can reduce the requirement of the terminal guidance module in the conceptual design process of air-launched vehicles. In addition, the results (Case 1|3, Case 2|4, Case 5|7 and Case 6|8) show the better global search ability of SAO than PSO, whether in a three-stage or four-stage configuration.
As for the constraints, the results (Case 1|3, Case 2|4, Case 5|7 and Case 6|8) show that the surrogate-based SAO has much better performance and stability in dealing with the constraints compared to PSO; the eccentricity is controlled on the order of magnitude 10 −4 , while the PSO shows unstable performance in dealing with the constraints.
The optimization convergence procedures of the design cases are shown in Figure 9; the control parameters of the PSO, including inertia weight, learning factor, and number of populations, are set as ω = 0.2, c = 0.8, N pop = 20. The results (Case 1|2,5|6) show that employing the terminal guidance module in the conceptual design process can improve both the convergence speed and the strength of the global search capability of the objective function. The oscillation in the convergence procedure in cases 3, 4, 7, and 8 show the effect of the multistage infilling strategy applied in the SAO method, as it tried several times to jump out of the local optimal solutions. In addition, the run time of PSO is over ten times longer than that of SAO, which demonstrates the great performance of the proposed SAO.

Performance Analysis
The flight simulation and performance of the air-launched vehicles using the optimal solutions of the overall parameters in two configurations are shown in Figure 10. The flight characteristic parameters include the altitude, Mach number, velocity, dynamic pressure, thrust, and gross mass of the launch vehicle. The title and units of each curve are shown in the legend of each figure with the same color. Figure 10a,c shows the flight altitude, thrust, Mach number, and velocity of the optimized air-launched vehicles. Both of the altitude curves show that the launch vehicles successfully reach the target injection point at 500 km, which means that the optimization results meet the constraint of the altitude. The trust curves show the time sequence and average trust of each stage. Figure 10b,d shows the local slope angle, weight, and dynamic pressure of the optimized air-launched vehicles. The curves of the local slope angle show the flight attitude of the launch vehicles after separating from the carrier aircraft. Furthermore, the optimization results meet the conditions of all the constraints with an error rate under 0.05%.

Conclusions
This paper proposes a conceptual design and overall parameters optimization framework for air-launched vehicle systems using a modified surrogate-based SAO method. A multidisciplinary model including aerodynamic, propulsion, mass, trajectory, and terminal guidance was established in order to simulate and analyze the flight performance of the airlaunched vehicle system. In order to improve the optimization efficiency and performance in searching for the global optimal solution, the calling times of the time-consuming multidisciplinary model were reduced by a proposed local density-based RBF surrogate model. The surrogate-based optimization problem can be solved by the proposed PSO method using adaptive control parameters. In addition, a multistage adaptive infilling strategy with a penalty method was employed to balance the exploration and exploitation of SAO. According to the results of eight design cases, employing the terminal guidance module in the optimization process can improve both the convergence speed and the global search strength and capability of PSO. The proposed SAO method only relies on the input and response of the original model, without requiring the internal information, which means that it can reduce the requirements of the terminal guidance system in the conceptual design process of air-launched vehicles. For the overall parameters optimization problem concerning air-launched vehicles, the computational costs for locating optimum results are reduced by an order of magnitude when the proposed SAO is used. The optimal solution increased the payload capacity of the three-stage air-launched vehicle by 82.55 kg (27.52%) and the payload capacity of the four-stage air-launched vehicle by 70.06 kg (23.35%) with the same lift-off mass of the original scheme (300 kg for 500 km SSO). Thus, combined with a surrogate optimization technique and reliable terminal guidance module, the proposed SAO method and multidisciplinary optimization framework are suggested for solving complicated aerospace engineering problems for spacecraft chasing high payload capacity, rapid response, and high reliability.
Author Contributions: J.L. was responsible for multidisciplinary modeling, verification, and developing optimization methods; D.W. was responsible for providing resources and software technical support; W.Z. was responsible for the structure and content of the article and for providing technical support. All authors have read and agreed to the published version of the manuscript.