Research Article | Open Access
Multifidelity Multidisciplinary Design Optimization of Integral Solid Propellant Ramjet Supersonic Cruise Vehicles
Integral solid propellant ramjet (ISPR) supersonic cruise vehicles share the characteristic that they are highly integrated configurations. The traditional design of vehicles cannot achieve a balance between computational expense and accuracy. A multifidelity multidisciplinary design optimization (MDO) platform has been developed in this study. The focus of the platform is on ISPR supersonic cruise vehicles. Firstly, codes of discipline with different levels of fidelity (LoF) were established, such as geometry, aerodynamics, radar cross-section calculations, propulsion, mass, and trajectory discipline codes. Secondly, two MDO frameworks were constructed through discipline codes. A low LoF MDO framework is suitable for conceptual design, and a medium LoF MDO framework is suitable for preliminary design. Finally, taking the optimization problem with the minimum overall detection probability of flight trajectory as an example, the low LoF framework first explores the entire design space to achieve the mission requirements, and then, the medium LoF MDO framework accepts the low LoF framework optimization parameters. Hence, the optimization target is reached with more detailed parameters and higher fidelity. Additionally, an example for a solid propellant missile with minimum total mass is tested by the platform. The study results show that the multifidelity MDO framework not only exploits interactions between the disciplines but also improves the accuracy of optimization results and reduces the iteration time.
With the development of flight vehicle design, the complex and nonlinear system is no longer the prestige of single disciplines, and the relationship between disciplines is increasingly complicated. Designers must search for optimal design plans, starting from the overall concept, by considering the interaction between disciplines. Multidisciplinary design optimization (MDO) is a methodology to solve this complex and nonlinear problem.
MDO is a methodology for design of systems where the interaction between several disciplines or subsystems must be considered and where the designer is free to significantly affect the system performance in more than one discipline or component . Flight vehicle MDO has developed into a mature system. One of the first applications of MDO was aircraft wing design, and now, MDO has been increasingly implemented in the field; as a result, many MDO frameworks have been applied to different types of vehicles. Baker et al. developed the IHAT system focusing high-speed airbreathing flight vehicles [2, 3]. Alston et al. developed the HFMDO system for the next-generation supersonic aircraft .
To a certain extent, the fidelity of the MDO framework depends on the fidelity of disciplinary codes. If the disciplinary codes are not well developed, the MDO system may be useless or worse . Sobieszczanski-Sobieski and Haftka and Roshanian et al. have highlighted that the two main challenges of MDO are computational expense and organizational complexity [1, 6] and that optimization solution of engineering problems mainly suffer from the curse of dimensionality . Therefore, multifidelity MDO was developed by researchers. Geiselhart et al. developed integration of multifidelity MDO codes for a supersonic aircraft, but those codes only were discussed in conceptual level design . Roshanian et al. developed a multifidelity MDO system for small solid propellant launch vehicles by conceptual design and preliminary design . Feng et al. proposed a multiobjective optimization-based effective global optimization for addressing the balance between local and global explorations . Bataleblu et al. introduced computational intelligence to optimize a robust trajectory .
Viana et al. think that the use of metamodeling techniques in MDO has evolved in the 29 years, and multifidelity approximations of metamodels are emphasized . By processing a large quantity of low level of fidelity (LoF) data and a small amount of high LoF data, by exploiting the advantages of the high LoF method with high accuracy and the low LoF method with low computational cost, a multifidelity method can be used to address the computational cost problem of MDO. Some MDO frameworks have adopted multifidelity methods to improve computational accuracy [6, 8]. Multifidelity optimization strategies are usually established within these frameworks. Low LoF strategies are used to calculate part of the optimization results and then refine the optimization results with medium or high LoF strategies.
Radar cross-section (RCS), as an important factor in the design and mission planning of a vehicle, has been applied in the trajectory design of unmanned combat aerial vehicles, where the problem of unmanned combat aerial vehicle path planning was formulated as an optimal control problem by Kabamba et al. . Low observability motion planning formulations and the pseudospectral multiphase optimal control are proposed . As an important technical indicator for supersonic cruise vehicles, the RCS plays an important role in the design of the vehicle shape and trajectory. Due to its importance, this study incorporates RCS into the MDO frameworks.
Compared with solid propellant launch vehicles, the relationships between disciplines are more complicated for integral solid propellant ramjet (ISPR) supersonic cruise vehicles. Firstly, the propulsion codes are equivalent to those containing solid propellant motor codes and solid ramjet codes, and secondly, the ISPR working conditions are closely related to the flight conditions, such as height, Mach number, angle of attack, air capture of the inlet, and throat radius. Meanwhile, the air capture of the inlet is closely related to the inlet size, which also significantly affects the aerodynamic characteristics. In the traditional preliminary design stage, each iteration in every discipline takes several days or even weeks. MDO helps to reduce the complexity of designing the process and shortens the development cycle of these vehicles.
Figure 1 shows the baseline supersonic cruise vehicle configuration used in this study. It also shows the vehicle geometry parameterization. The vehicle is divided into a payload cabin, an equipment cabin, and ISPR. The inlets and rudders are arranged outside the vehicle. The vehicle has a plane symmetry shape and twin inlet 90° configuration. The internal packaging layout and configuration are determined by external geometry. Changing the length and diameter can significantly change the engine charge and aerodynamic characteristics.
Based on the multifidelity strategy, we began with the conceptual and preliminary designs of the vehicle then established two MDO frameworks. At the system design level, different fidelity strategies were adopted. Based on the characteristics of the conceptual and preliminary designs, different fidelity MDO frameworks satisfying these two different design stages were established. Different from the multifidelity MDO system designed by Roshanian et al. , the MDO frameworks combine Variable Complex Modeling (VCM) with multilevel of fidelity MDO. VCM is a method which can replace low LoF codes by medium LoF codes. Different from VCM, not only low LoF disciplinary codes but also high LoF disciplinary codes are contained in the low LoF MDO framework in this system. There is the same condition in the medium LoF MDO framework. At the discipline level, the fidelity of codes was increased by using different multifidelity methods according to the characteristics of the disciplines. Geometry, aerodynamics, propulsion, mass, and trajectory disciplinary codes were used. RCS calculations were also incorporated into the MDO frameworks. The focus of this study was the ISPR supersonic cruise vehicle. An MDO system that satisfies the conceptual and preliminary designs of the vehicle was established, and an example was optimized. The geometry, propulsion, and trajectory parameters were adjusted to solve the problem of the minimum detection probability of the flight path.
The paper is organized as follows. Section 2 introduces the basic idea of multifidelity strategies and methods. Section 3 presents disciplinary codes depending on different fidelities. The calculation time and accuracy of those codes also are shown in this section. Two examples are presented in Section 4. One of them is an application for the supersonic cruise vehicle using ISPR. Another is a simple example for a small solid propellant missile. Finally, the conclusions are reported in Section 5.
2. Multifidelity Calculation
2.1. Concept of Multifidelity Strategy
Fidelity refers to the validity and reliability of the entire MDO process and optimization results. It depends on the fidelity of each discipline and the complexity of the optimization process; the more complex the optimization process and the higher the discipline fidelity, the higher fidelity of the optimization.
Multifidelity, also known as variable fidelity and variable complexity, provides a way to minimize high LoF models at reduced computational cost . Alexandrov and Lewis give an overview of first-order approximation/model management optimization and a study of variable-fidelity physics models for airfoil . Gano et al. used low-fidelity models and a scaling function to approximate the high-fidelity model; a Kriging-based scaling function was introduced to better approximate the high-fidelity response . A multifidelity technique was employed for a gust performance evaluation problem by Berci et al., where a metamodel of the high-fidelity model was built based on a tuned low-fidelity aeroelastic model . Balabanov and Venter introduced a new approach to multifidelity optimization in reference, where correlation between the results of the low LoF and high LoF analyses is not required . An approach for combining conceptual design and preliminary design for wing optimization was introduced by Hutchison et al.; multifidelity design strategies were used to combine conceptual and preliminary designs . Similarly, different fidelity models are used to establish conceptual and preliminary design models in this work.
The MDO process depends on the design stage of the vehicle. Conceptual design of the vehicle refers to an initial demonstration based on technical indicators, forming a basic feasible concept and scheme and determining the overall parameters of the vehicle. Preliminary design is based on the concept design results, which are used as the initial iterative value, and the technical indicators and overall parameters in order to demonstrate feasibility and form a more detailed design. This allows more parameters to be determined, not only the overall parameters but also the parameters for all of the disciplines. In these two stages, the application of MDO can shorten the prophase design cycle and reduce the design cost. The MDO process in conceptual design is simpler and therefore less reliable.
There are two cases of discipline fidelity calculations. The first case is in the aerodynamic and RCS calculation codes. The high LoF method uses the same input as the low LoF method, but it requires too much calculation time and computing resources. As such, a small amount of high LoF data and a large amount of low LoF data are combined to form a medium LoF calculation method in order to increase the fidelity of the codes. The other case is in mass and propulsion codes, where the required inputs of the high LoF method are more detailed and accurate. However, the computational cost does not increase significantly in this case. The conceptual design stage cannot provide such detailed data, so the low LoF method is used for conceptual design after satisfying the calculation requirements of the high LoF method. This method performs an iterative calculation in the preliminary design phase.
A multifidelity strategy was adopted for MDO in this study. Multifidelity strategy refers to use of different fidelity discipline codes to establish different fidelity MDO frameworks for different application scenarios. Different from other multifidelity MDO frameworks, the low LoF MDO framework does not mean that all disciplinary codes are low LoF, while medium LoF MDO framework does not mean that all codes are medium LoF. Similar to the cask effect, the fidelity of the framework is determined by the lowest fidelity disciplinary codes. In the conceptual design phase, a low LoF calculation strategy is adopted. For aerodynamic and RCS calculations, a low LoF calculation method is adopted; the low LoF calculation method is still used for other discipline codes such as propulsion. Through iterative optimization calculations, the MDO results at the conceptual design stage can be formed. Once the conceptual design is complete, the preliminary design begins. At this time, a medium LoF strategy is adopted. The multifidelity calculation method is adopted for all disciplines, and the initial value for the iterative optimization is the optimum value of the conceptual design.
The architectures of MDO are known as the all-at-once (AAO) problem, which is a general formulation according to Martins and Lambe . Figure 2 shows the relationship between two multifidelity strategies. Figure 2(a) shows the low LoF MDO strategy in the conceptual design phase, and Figure 2(b) represents the medium LoF MDO strategy in the preliminary design phase. The low LoF MDO strategy contains the codes: low LoF geometry, low LoF aerodynamics, low LoF radar cross-section calculations, low LoF propulsion, low LoF mass, and trajectory. The medium LoF MDO strategy contains the codes: high LoF geometry, medium LoF aerodynamics, medium LoF radar cross-section calculations, high LoF propulsion, high LoF mass, and trajectory. The orange full line box presents the high LoF method, the green narrow dotted line presents the medium LoF method, and the blue wide dotted line presents the low LoF method.
2.2. Multifidelity Methods
The multifidelity calculation method forms medium LoF codes by combining low LoF discipline codes with high LoF discipline codes.
For two different fidelity methods solving the same problem, the high LoF method can provide high confidence calculation results but often requires a longer calculation time and more computing resources; although the low LoF method is less accurate, its computational efficiency is much higher than that of the high LoF method. Using the characteristics of both methods, the multifidelity calculation method can make up for the poor accuracy of the low LoF method and can significantly reduce the time by the high LoF method.
A commonly used multifidelity method is the polynomial fitting method by Doyle et al.  and Le et al.  applied in building disciplinary codes, which uses polynomial interpolation to fit the high LoF data with the low LoF data. However, the error has been shown to be large. The Co-Kriging method is an improved method based on the Kriging method; it is presented by Myers  and can be used to construct a multifidelity surrogate model with low LoF data and high LoF data. Its accuracy has been established and validated in aerodynamic calculations by Forrester et al.  and Huang et al. . Furthermore, an extended Co-Kriging method with higher accuracy has been developed by Xiao et al. . Co-Kriging approximation models using an automated Euler and Navier-Stokes-based method had been developed by Chung and Alonso . Kim et al. made a comparison study on the accuracy of the technique of nonconvex functions and got the conclusion that Kriging can create a more accurate metamodel than support vector regression . Yuan and Bai made a comparison of the Kriging method and neural network and got a conclusion that Kriging approximation is likely to be more accurate . A variety of multifidelity calculation methods were compared in this study, and the Co-Kriging method, which has the highest accuracy, was used as the multifidelity calculation method for aerodynamics and other disciplines.
The specific operation steps for the multifidelity method are as follows:
Step 1: select the state to be calculated according to the specific task of the discipline
Step 2: use the low LoF method for this discipline to calculate the data under all states to generate a low LoF database
Step 3: select the characteristic calculation state from all the calculation states using the uniform test or Latin hypercubic experiment design method. Use the selected characteristic calculation state and select the characteristic database based on the low LoF database, where the characteristic database can represent natural features of the low LoF database
Step 4: use the high LoF method to calculate the data under the characteristic calculation state and then generate a high LoF database
Step 5: multifidelity data can be generated with higher accuracy
The multifidelity data by the specific operation steps also are shown in Figure 3.
In this section, we describe in detail the three types of multifidelity calculation methods: polynomial interpolation fitting, delta-Kriging, and Co-Kriging.
2.2.1. Polynomial Interpolation Fitting
The polynomial interpolation fitting method is used in many MDO frameworks [5, 6]. This method is easy to implement and has a certain degree of accuracy. It is based on the difference value between the high LoF data and the characteristic data. It carries on the polynomial to merge the interpolation value. The interpolation results and the sum of the low LoF data represent medium LoF data. Polynomial interpolation fitting is done using the equation where we used the cubic spline interpolation method.
The greater the amount of high LoF data that is available, the higher the accuracy we obtained for the medium LoF data. Therefore, a sufficient amount of high LoF data is needed for this method, which means that a longer calculation time is required. For a computational module with time constraints, using this method means that a certain amount of computational accuracy is lost.
2.2.2. Delta-Kriging Method
This method is based on the traditional surrogate model method, to model the difference between the high and low LoF codes. This difference is used to modify the low LoF code errors, in order to improve accuracy. The difference surrogate model method uses the equations where is the surrogate model method using . We used the Kriging method to build the surrogate model.
As with the polynomial interpolation fitting method, the delta-Kriging method deals with the difference between the high and low LoF data. The latter uses the surrogate model to predict the difference; it is sensitive to the selection of the characteristic flight state, but it has higher precision than the latter; the former uses polynomial interpolation to fit the difference, which is less sensitive to flight status selection but has low precision.
2.2.3. Co-Kriging Method
The Co-Kriging method is an extension of the Kriging method, using the data of two independent high and low LoF data sets, and can make full use of the existing sample data to excavate the high LoF code information, while reducing the quantity of sample data and the computational cost. The Co-Kriging method is based on the following equation:
Compared with the traditional Kriging method, this method has two different correlation functions; because two separate sets of sample data are used, the relevant parameters , , , , and can be calculated [24, 30].
The fundamental point which makes the Co-Kriging method different from the traditional Kriging method is . In the traditional Kriging method, it can be regarded as 1, but in the Co-Kriging method, it must be obtained through global optimization (using a genetic algorithm, for example) with other related parameters and is not usually 1.
2.2.4. Two-Variable Demonstration
For certain shape of the vehicle, the variables of aerodynamic data are two or three flight conditions, which are , , and . The same situation happens with the RCS calculation of and .We used an example with two variables to test the feasibility of the three methods and evaluate their errors. Let us assume that the high LoF model could be represented by the function , the low LoF model by , and the sampling method was a Latin hypercubic sampling method. A total of 100 low LoF sampling points were extracted. Figure 4 shows the high and low LoF models.
Figures 5–7 show the results of the three methods with a different number of high LoF sampling points. They are a result of the change in the confidence when the numbers of high confidence points are 4, 9, or 16. The stars represent the data at the low LoF sampling points, and the squares represent the high LoF data sampling points. The closer the model is to the high LoF model, the further the model is from the stars, the higher the fidelity. When there are 4 high LoF data points, the Co-Kriging method fitting effect is significantly higher than the other two methods. The delta-Kriging and polynomial methods produced similar results, but the polynomial fitting result is close to the low confidence data. When the number of high LoF points was increased to 9, the fidelity of the Co-Kriging method was higher than it was with 4 data points, and the polynomial method is slightly better than the delta-Kriging method. When the number of high confidence sampling points was increased to 16, none of the methods improve their precision.
The delta-Kriging method and the polynomial method are similar in principle. Both of them deal with the difference between the high and low LoF, so it is reasonable that results were similar.
The mean relative square error (MRSE) and the maximum relative error (MRE) are used as criteria to test the error prediction of the multifidelity methods. Table 1 shows the MRSE of the three methods, which were established with different numbers of high LoF sampling points and the high LoF model. Table 2 shows the MRE of the three methods. These results further indicate that the Co-Kriging method was the best of the three methods. The delta-Kriging method and the polynomial method have similar results. Therefore, we used the Co-Kriging method when constructing a multifidelity MDO framework:
3. Math Discipline Codes
The relationship between the discipline codes is shown in Figure 8. The coupling relationship between codes is complex, and they affect each other. The main function of the geometry codes is to establish a parameterized vehicle shape and internal structure and to provide the geometric data for the aerodynamics, RCS calculation, propulsion, and mass codes. The main function of the aerodynamic discipline codes is in the defined geometric shape parameters used to calculate aerodynamic coefficients, and RCS discipline codes are used to calculate RCS data for trajectory disciplines in the defined geometric shape parameters. The propulsion codes calculate parameters such as thrust, specific impulse, and engine quality of ISPR, then provide these parameters to the mass and trajectory disciplines. The mass discipline codes provide the vehicle mass characteristics and equipment parameters based on the geometric shape and mass parameters of the ISPR for the trajectory discipline. Finally, the trajectory discipline facilitates the calculation of a specific flight trajectory and parameters using the aerodynamic, RCS, mass, and propulsion data; it also provides overload data to the mass codes. Table 3 shows the interaction parameters between all of the discipline codes in Figure 8.
3.1. Aerodynamic Codes
Common aerodynamic calculation methods include the engineering estimation method and the Computational Fluid Dynamics (CFD) method. As a low LoF method, the advantage of engineering estimation is the fast calculation speed which means that dozens of computing states can be completed in a short time; the disadvantage is that the calculation results are not accurate enough. The component combination method and surface element method are common engineering estimation methods. CFD can be divided into two parts: the viscous flow Navier-Stokes equation and the inviscid Euler equation . Viscous flow involves the flow of transport phenomena, such as friction, heat conduction, and mass diffusion. These transport phenomena are dissipative, and they always increase the entropy of the fluid. It is obvious that the definition of inviscid flow ignores dissipation, viscous transport, mass diffusion, and the flow of heat through conduction. The assumption of aerodynamic calculation codes is based on the inviscid flow conditions. As a high LoF method, the advantage of CFD is that the accuracy of the calculation results is very high. The disadvantage is that the calculation process is time-consuming. The calculation often takes hours or even days. The calculation of the viscosity is often more time-consuming than the inviscid calculation.
3.1.1. Low LoF Aerodynamic Codes
The low LoF codes do not require highly reliable methods but do require high computational efficiency. To find the supersonic cruising vehicle’s speed characteristics, the low LoF method uses a theoretical calculation. The theoretical calculation method establishes a simplified mathematical model and a control equation based on the study of the basis of air motion law and assumptions. The body is solved using the rotating body theory. When the angle of attack is 0-5°, irrotational stream theory is used to find a solution. When the angle of attack is 5-25°, cross-flow theory is used to find a solution. The supersonic aerodynamic characteristics of the wing are based on cone-shaped fluid theory. The wing-body assembly and tail-body assembly are solved according to the component combination method. Therefore, the aerodynamic characteristics of the whole vehicle are those of the wing-body assembly and the tail-body combination section [32, 33].
3.1.2. High LoF Aerodynamic Codes
The high LoF method needs to be able to accurately describe the aerodynamic status of the aircraft. Therefore, it requires a high-order aerodynamic calculation method. The three-dimensional inviscid Euler equation was used in this study. The finite volume method was used based on the integral form of Euler’s equation. Firstly, the computing area was approximately discretized into a finite number of nonoverlapping grids. Secondly, a series of nonoverlapping control units around each grid point were selected with only one node in each control unit. Thirdly, the flow rate to be determined was set on the grid node, and each control unit was integrated using the law of the conservation of the flow rate to derive a set of discrete formats. Finally, the solution was determined, and the numerical solution of the flow was obtained [31, 34].
In this method, the Runge-Kutta method was used for time discretization, and the spatial discretization is in an upwind format and total variation diminishing (TVD) scheme difference scheme.
For certain shape, we chose the flight conditions at and . Mesh independence is described. Different mesh numbers such as 0.5 million, 0.75 million, and 1 million are tested. The results are shown in Table 4. The results indicate that mesh number of 0.75 million can meet the calculation accuracy requirements.
For a simple aircraft shape and a moderate mesh size, the calculation time of this method is short, making it suitable for MDO calculation.
3.1.3. Medium LoF Aerodynamic Codes
For a certain step in the optimization iteration, the MDO system automatically generates an aircraft shape, which is required to calculate the aerodynamic table based on this shape and flight conditions. During each round of optimization, for the certain shape, the aerodynamic table is a function of flight conditions. The process for medium LoF aerodynamics is as follows:
Step 1. Getting the calculation status. The calculation status required for trajectory calculation was obtained, generally including the flight Mach number, angle of attack, and flight altitude. These statuses are the basis for the formation of aerodynamic tables
Step 2. Generating the low LoF aerodynamic table. For all calculation statuses, a low LoF aerodynamic table can be formed in a very short time by using the theoretical calculation method. This table is identical to the table used in the low LoF MDO framework
Step 3. Selecting the characteristic status. The characteristic status are sampled in all calculation status by the Latin hypercube sampling method. In general, if the number of factors in calculation status is 2, the sampling point is no less than 9, and if the number of factors is 3, the sampling point is no less than 16 
Step 4. Generating the high LoF aerodynamic table. For all characteristic status, a high LoF aerodynamic table can be formed using the high LoF calculation method. This table is identical to the table used in the high LoF MDO framework
Step 5. Forming the final medium LoF aerodynamic table or model. The premise of generating a medium LoF aerodynamic table is a low LoF table containing a large amount of data and a high LoF table containing a small amount of data. The high LoF data is used to correct the low LoF data. Using the Co-Kriging method, the two kinds of tables can be used as input to the method to quickly generate a medium LoF model. A medium LoF aerodynamic table can be generated by bringing in all calculation statuses. The medium LoF model also can be an input of the trajectory calculation
RCS Calculation Codes
The fundamental solution to high-frequency scattering RCS analysis is to solve Maxwell’s equations. Solving the equations directly provides a precise solution for the spatial electromagnetic field and the RCS. The main ways of doing this are the moment method, finite element method, time domain finite element method, fast multisubstage method, and multilayer fast multisubstage method. Although all of these numerical methods are highly reliable, they require a substantial amount of calculation.
The RCS approximate calculation method uses approximation and asymptotic analysis based on the strict numerical analysis of Maxwell’s equations and mathematical operation skills. Mature methods for doing this include geometric optics, physical optics, the geometrical theory of diffraction, physical theory of diffraction, incremental length diffraction coefficients, and method of equivalent current . Research shows that these approximate calculation methods can meet the requirements of the preliminary design phase and complement each other. The results are more accurate when multiple methods are used in combination .
3.1.4. Low LoF RCS Calculation Codes
The low LoF RCS method uses the physical optics method [38, 39]. The electromagnetic induction field of the target surface is approximated, and the scattering field is obtained by integration. To facilitate the solution, Maxwell’s equations are transformed into the Stratton-Chu integral equations : where ψ is Green’s function, given by
The following basic assumptions are made for the Stratton-Chu equation: (a) electromagnetic waves do not diffract on the object surface; they only reflect; (b) the total scattered field of the target shadow is zero; (c) the distance from the scatterer to the radar receiver is much greater than the scatterer size; and (d) when the radius of curvature of the target surface is much larger than the wavelength, the tangent plane is used for calculation instead of the original surface.
Based on the assumptions of the physical optics method without considering the mutual interference of the scattering fields of each element, for a specified incident field direction, the total scattered electric field can be obtained by successively calculating the scattered electric field of each surface element.
The steps of the low LoF RCS calculation are shown in Figure 9. Step 1: reading the geometry model of the vehicle. Step 2: setting up calculation conditions: range of view angle and incident power. Step 3: calculation of the surface feature parameter, element physical optical current, and surface element scattered electric field. Step 4: loop computing in all surface elements, summing all scattered electric field. Step 5: loop computing all direction of the view angle, calculating RCS data. Step 6: RCS data output, where the view angle is an angle of the sea level and the line from radar to vehicle.
3.1.5. High LoF RCS Calculation Codes
A high LoF method is needed to describe the RCS characteristics more precisely than the low LoF method. The time domain finite element method within the numerical analysis method is used to solve Maxwell’s equations . The basic idea is as follows: Maxwell’s equations are differentially discretized; then, they are converted from differential equations to difference equations using the center difference method; starting from the wave source and directly describing the electromagnetic wave propagation process at the moment , the finite-difference calculation of each spatial grid allows waves to propagate in a simulated space. Data is sampled continuously within a certain volume and period in a continuous electromagnetic field, so the computational process completely simulates the electromagnetic wave propagation and interaction with the target object. Eventually, an accurate solution of the spatial electromagnetic field is obtained; the RCS can then be found using the near-far field conversion.
3.1.6. Medium LoF RCS Calculation Codes
Once the vehicle coating material has been determined, its RCS data depends on its shape. The following basic assumptions are made for the medium LoF RCS calculation: in the MDO process, the topological structure of the vehicle’s shape will not change, and the shape parameters will change within a certain range, so the RCS data of different aircraft shapes will not change greatly. The following medium LoF method is used:
Step 1. Getting the calculation status. The calculation status required for trajectory calculation was obtained, generally including all view angles. For a certain shape, these statuses are the basis for the formation of RCS tables
Step 2. Generating the low LoF RCS table. For all calculation status, a low LoF RCS table can be formed in a very short time by using the physical optics method. This table is identical to the table used in the low LoF MDO framework
Step 3. Selecting the characteristic status. The characteristic status is sampled in all calculation statuses by the Latin hypercube sampling method
Step 4. Generating the high LoF RCS table. For all characteristic statuses, a high LoF RCS table can be formed using the time domain finite element method. This table is identical to the table used in the high LoF MDO framework
Step 5. Forming the final medium LoF RCS table or model. Using the Co-Kriging method, the two kinds of tables can be used as input to the method to quickly generate a medium LoF RCS model
Step 6. Gauss filtering. The RCS data model can be processed by Gauss filtering to improve the computational speed in trajectory calculation
3.2. Propulsion Codes
The propulsion codes are used for the internal trajectory calculation model of the ISPR. The ISPR can be considered in two parts: the solid propellant motor and the solid fuel ramjet, and their calculation models are established accordingly. When a solid propellant motor model is used to replace a booster model in the ISPR, a solid fuel ramjet model is used instead of the gas generator model. The low LoF code is used in the conceptual design phase, which cannot build a detailed model; it simplifies the ISPR codes because fewer parameters are available for the vehicle and ISPR. A high LoF code is used in the preliminary design stage. Compared with the previous stage, more parameters are available for the ISPR. This is sufficient to establish a more detailed model of the power system and can describe the working state of the ISPR more comprehensively. It can also improve the accuracy of the trajectory codes.
The following basic assumptions are made for propulsion codes: (a) the internal gas flow in the engine is one-dimensional flow; (b) there is no energy dissipation; (c) there is no ablation of the thermal insulation layer; and (d) the flow rate of the ISPR can be continuously adjusted.
3.2.1. Low LoF Propulsion Codes
The low LoF model for a solid propellant motor assumes that the engine charge a specific impulse, and the flow rate is constant. Therefore, the engine can provide a fixed thrust. The thrust model is
The low LoF model of the solid fuel ramjet assumes that the specific impulse remains unchanged, but it can adjust the flow rate. Hence, thrust is a variable, and the specific value is where
3.2.2. High LoF Propulsion Codes
In the high LoF model of the booster motor, the theoretical specific impulse is calculated using the theoretical extrapolation function. The influence of the nozzle expansion ratio and the pressure of the combustor are considered when the interpolation is carried out. On the basis of the thermodynamic calculation, the evaluation of the interpolated function  is
The relative error between the result of this extrapolation formula and the result of the thermal calculation does not exceed 0.1%. At the same time, the actual specific impulse is also calculated by the Solid Performance Program (SPP) empirical formula . The flow rate is determined by its charge and burning time.
The high LoF codes for solid fuel ramjet are calculated using the nondesign point performance. The total pressure recovery coefficient and inlet flow coefficient are calculated by the following equations :
The size of the gas generator and the pressure of combustor are dependent on the shape of the vehicle and mass of fuel-rich propellant. After determining the type of propellant, thermodynamic calculations can be carried out. The composition, temperature, and other thermodynamic parameters of combustion products of combustor can be obtained. According to those thermodynamic parameters, parameters of afterburning combustor can be obtained, which are , , and . The total pressure recovery coefficient of the nozzle can be obtained by experience in preliminary design. A side air inlet internal trajectory calculation program is used . The thrust of the high LoF model is
3.3. Mass Codes
At different stages of vehicle design, different mass codes can be used for calculation purposes. The low LoF model has a poor calculation accuracy but requires fewer parameters and is generally used in the conceptual design stage. Conversely, the high LoF model used in the preliminary design stage has a high calculation accuracy, but it requires more complicated parameters. However, the computational challenge of them are almost equal and simple.
The following basic assumptions are made for mass codes: (a) the number of internal parts is known and no longer increases; (b) the mass, mass center, and shape of all internal parts are known; and (c) since the aircraft is in a supersonic state, the mass loss caused by thermal ablation is ignored.
3.3.1. Low LoF Mass Codes
The low LoF method simplifies the vehicle shell into a thin shell of even thickness, so the mass of the vehicle can be determined. The mass discipline codes can be simplified and divided into four parts:
The mass center of the vehicle is where and are calculated from the low LoF propulsion codes. and are dependent on equipment in the vehicle, which are fixed. and are calculated by the thin shell of the simplified vehicle shell.
3.3.2. High LoF Mass Codes
The high LoF codes model the mass in detail. Firstly, the vehicle is divided into different structural sections according to their function. Then, the models of each section are established by structural and geometrical parameters. Finally, according to the materials and density of each section, we can get the cabin mass and position of the center of mass. The vehicle mass and mass center can be quickly calculated using the information above: where is the number of sections.
3.4. Trajectory Codes
The supersonic cruise aircraft is a plane-symmetrical aircraft, so it uses Bank-To-Turn (BTT) for control. During the flight, the slide angle is 0, and the trajectory is controlled through the angle of attack and roll angle. The simulation uses a trajectory with three degrees of freedom [46–48], and the state variables include velocity, path angle, deflection angle, vehicle position, and mass.
The following basic assumptions are made for propulsion codes: (a) the earth’s rotation and curvature are neglected; (b) the rotational motion caused by the moment of inertia is neglected; (c) the slide angle is 0; and (d) the aircraft is in transient equilibrium all the time: where
In order to demonstrate the influence of the RCS on the trajectory, we need to add the low detectability constraint, that is, to increase the radar set point on the route of the trajectory. Constraints can include the instantaneous detection probability of the radar effective area or the overall detection probability of the entire flight. The control variables of the aircraft during flight are the angle of attack and roll angle. If the vehicle is required to reach a specified target point, the end position is added as a constraint condition. The model of the instantaneous detection probability of the radar effective area is
The model of the overall detection probability is
The complicated constraints on trajectory codes mean that that the derivation process is cumbersome, and the convergence domain is small when using Pontryagin’s extremum principle. The pseudospectral method  is used to transform the problem from a continuous system into a nonlinear programming problem. In order to increase the convergence rate of the pseudospectral method, polynomial fitting is carried out on the aerodynamic data, and the RCS data is Gauss filtered before polynomial fitting. By using the pseudospectral method, the optimal control problem of the trajectory is solved during each iteration of MDO.
3.5. Analysis of Discipline Codes
For the vehicle in this study, the computational time and accuracy of each fidelity disciplinary codes are compared. The computational time is shown in Table 5. The computational time of low LoF aerodynamics and RCS codes is lower than those with high LoF. And an average of 19 times of the time is reduced by using medium LoF codes, which can solve the problem of computational expense. Due to the fact that there are some nonexistent codes, the computational time of them is not counted in Table 5.
The accuracy of each codes are shown in Table 6 in which data used for comparison are as follows: in the status of and ; in the status of and ; of ISPR in , , and m; and .
The calculation time of the two MDO frameworks in Figure 2 is tested. A high LoF MDO framework consisting entirely of high LoF disciplinary codes is also tested. All of the test data is listed in Table 7. The conclusion can be drawn from the test that the medium LoF framework can reduce calculation time and improve calculation accuracy in acceptable limits. However, the low LoF framework of the conceptual design has a low fidelity.
4. Example Problem
We used the different fidelity codes of different disciplines to build the two MDO frameworks in Figure 2, meeting the conceptual design and preliminary design requirements. Through the low LoF MDO framework, the conceptual design of a supersonic cruise vehicle was completed.
4.1. Optimization Problem and Constraints
The optimization object is an ISPR supersonic cruise vehicle. The launch platform is a supersonic aircraft, and the mission target position of the flight trajectory is known. A noncooperative radar is located along the longitudinal flight trajectory. The optimization problem is to reduce the overall detection probability for the entire flight by modifying the shape, ISPR, and flight trajectory parameters. There are 34 input parameters, of which 14 are related to the shape, 11 to the ISPR, and 9 to the trajectory. All input parameters and their initial values are shown in Table 8.
The sensitivity analysis of all 34 parameters for the overall detection probability is calculated, and the results are shown in Figure 10. In this figure, all the parameters are arranged according to the sensitivity. The most influential parameters are gas generator total impulse, cone length, column length, diameter, mass of gas generator, total time, and mass of booster. In the MDO process, special attention needs to be paid to the influence of the changes of these parameters on the overall detection probability during the flight of the aircraft.
4.2. Solutions and Conclusion
Genetic algorithms were used as optimization algorithms for this example. The conceptual design took 65 hours. After 243 generations of optimization calculations, the overall detection probability was reduced by 66.9%, and the effect was considerable. The preliminary design optimization took a total of 221 hours. After 108 generations of optimization calculations, the overall detection probability of the mission trajectory was reduced by 5.76%. A comparison of all output parameters is shown in Table 9. Figure 11 shows the lateral trajectory before and after optimization, including the initial trajectory, the optimal trajectory of conceptual design, and the optimal trajectory of preliminary design.
For simplifying the calculation, we set the booster mass and booster special impulse to fixed values. From Table 9, we can conclude that (a) the cone length is increased, the column length is reduced, the total length is reduced relative to its initial value, and the result is that the total mass and the average RCS are reduced; (b) while the total length and the gas generator length are decreasing, the length of powder and the diameter are increased, and finally, the gas generator is increased; and (c) from the results, we can see that after optimization, the optimal trajectory can avoid the radar further the initial trajectory; the overall detection probability is significantly reduced.
In this case, increasing the cone length, reducing the column length, and reducing the total length can improve the aerodynamic performance and reduce the RCS. Since the range need to be increased to stay away from the radar, the mass and length of gas generator are increased. A combination of factors leads to the total mass reduction. Finally, the overall detection probability is significantly reduced.
After the optimal design values have been achieved at the conceptual design stage, the models with higher LoF are adopted. In the entire design exploration space, the preliminary design stage can achieve optimal results with only a small number of computational iterations. Although the results of preliminary design are not much different from that of conceptual design, the orientation of vehicle design and optimization is provided by these results.
4.3. Another Example
There is an additional example for the multifidelity MDO framework using a small solid propellant missile. The design target is the minimum total mass with the constraint of 350 km for range. Aerodynamic, propellant, mass, and trajectory codes are mentioned before. A comparison of all output parameters is shown in Table 10. In this case, total mass of the missile is reduced significantly. Due to the decrease of cone length, column length, and total length, total mass is reduced. To meet the predetermined requirement of range, and both are increased, and then, is increased. This simple example shows that the multifidelity MDO platform is not only suitable for the ISPR supersonic cruise vehicle but also suitable for the solid propellant missile.
The MDO method has a substantial influence on the design of the flight vehicle, and it is convenient when designing supersonic cruise vehicles.
We introduced multifidelity MDO for a supersonic cruise vehicle considering the RCS. Geometry, aerodynamics, RCS calculations, propulsion, mass, and trajectory discipline codes were used. Initially, two MDO frameworks were used for this complex system with coupling disciplines to find the optimal solution. The conceptual design MDO framework adopted a low LoF strategy, which provided the initial value for the preliminary design MDO framework. The preliminary design adopted a medium LoF strategy, which provided sufficient design parameters to improve the accuracy of the MDO. Secondly, the aerodynamic, RCS calculation, propulsion, and mass disciplines adopted a multifidelity model. In comparison to previous models which used empirical methods, the accuracy of model calculations improved. Finally, the designed MDO frameworks were used to explore the design space, provide weight targets for vehicle optimization, reduce the design workload, and shorten the design cycle.
By adopting multifidelity design at both the discipline and the system levels, it is possible to customize the LoF for each design stage in different disciplines. The optimization cases demonstrated that this multifidelity MDO framework performed well in minimizing the total mass or the overall detection probability of flight trajectory. This multifidelity MDO framework is a general approach for solving flight vehicles’ design problems.
|:||Inlet capture area|
|:||Parameters of radar|
|:||Lateral force coefficient|
|:||The low LoF database or method|
|:||The high LoF database or method|
|:||The characteristic database or method|
|:||A medium LoF database, also known as a calibration database|
|:||Predicted value of the surrogate model|
|:||Spline interpolation method, the gravitational acceleration|
|:||Actual specific impulse|
|:||Theoretical specific impulse|
|:||Theoretical air volume|
|:||Mass of the vehicle|
|:||Unit normal vector of the unit cell|
|:||Pressure of environment|
|:||Pressure of the combustor|
|:||Distance from the surface element to the observation point, distance from the radar to the flight vehicle|
|:||Initial time of trajectory simulation|
|:||End time of trajectory simulation|
|:||Flight status that generates low LoF data|
|:||Characteristic flight status (that generates the characteristic and the high LoF data)|
|:||Medium LoF surrogate model|
|:||Surrogate model constructed by the low LoF data|
|:||Surrogate model constructed by the difference between high and characteristic data|
|:||Angle of attack, residual gas coefficient|
|:||Angle of slide slip|
|:||Shock angle of the inlet|
|:||Nozzle expansion ratio|
|:||Relevant parameters of|
|:||Relevant parameters of|
|:||Turning angle of inlet|
|:||Wavelength, velocity coefficient|
|:||Amplification factor of the Co-Kriging method|
|:||Radar cross-section (RCS)|
|:||Total pressure recovery coefficient|
|:||Heating coefficient of the chamber|
|:||Inlet flow coefficient|
|:||Combined quantity of the chamber|
|:||Gas quality increase coefficient.|
|:||Throat of inlet|
|1:||Inlet section of inlet|
|5:||Outlet section of the nozzle|
|:||Integral solid propellant ramjet (ISPR)|
|:||Payload and equipment|
|:||Wings or rudders.|
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
The work was partially sponsored by the National Natural Science Foundation of China (No. 11772353). The authors would like to acknowledge the support of Zhiwei Feng and the rest of the team, without whose contributions this work would not be possible.
- J. Sobieszczanski-Sobieski and R. T. Haftka, “Multidisciplinary aerospace design optimization: survey of recent developments,” Structural Optimization, vol. 14, no. 1, pp. 1–23, 1997.
- M. Baker, K. Alston, and M. Munson, “Integrated Hypersonic Aeromechanics Tool (IHAT),” in 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, pp. 4–6, Atlanta, Georgia, September 2002.
- M. L. Baker, M. J. Munson, G. W. Hoppus, and K. Y. Alston, “The Integrated Hypersonic Aeromechanics Tool (IHAT), Build 4,” in 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Albany, New York, August 2004.
- P. K. Alston, S. Doyle, T. Winter, H. Kim, and S. Ragon, “High fidelity multidisciplinary optimization (HFMDO),” in 13th AIAA/ISSMO Multidisciplinary Analysis Optimization Conference, Fort Worth, Texas, September 2010.
- I. Kroo and I. Kroo, “Multidisciplinary optimization applications in preliminary design - status and directions,” in 38th Structures, Structural Dynamics, and Materials Conference, Kissimmee,FL,USA, April 1997.
- J. Roshanian, J. Jodei, M. Mirshams, R. Ebrahimi, and M. Mirzaee, “Multi-level of fidelity multi-disciplinary design optimization of small, solid-propellant launch vehicles,” Transactions of the Japan Society for Aeronautical and Space Sciences, vol. 53, no. 180, pp. 73–83, 2010.
- J. Roshanian, A. A. Bataleblu, and M. Ebrahimi, “A novel evolution control strategy for surrogate-assisted design optimization,” Structural and Multidisciplinary Optimization, vol. 58, no. 3, pp. 1255–1273, 2018.
- K. Geiselhart, L. Ozoroski, J. Fenbert, E. Shields, and W. Li, “Integration of multifidelity multidisciplinary computer codes for design and analysis of supersonic aircraft,” in 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, p. 465, Orlando, Florida, January 2011.
- Z. Feng, Q. Zhang, Q. Zhang, Q. Tang, T. Yang, and Y. Ma, “A multiobjective optimization based framework to balance the global exploration and local exploitation in expensive optimization,” Journal of Global Optimization, vol. 61, no. 4, pp. 677–694, 2015.
- A. A. Bataleblu and J. Roshanian, “Robust trajectory optimization of space launch vehicle using computational intelligence,” in 2015 IEEE Congress on Evolutionary Computation (CEC), pp. 3418–3425, Sendai, Japan, May 2015.
- F. A. C. Viana, T. W. Simpson, V. Balabanov, and V. Toropov, “Special section on multidisciplinary design optimization: metamodeling in multidisciplinary design optimization: how far have we really come?” AIAA Journal, vol. 52, no. 4, pp. 670–690, 2014.
- P. T. Kabamba, S. M. Meerkov, and F. H. Zeitz, “Optimal path planning for unmanned combat aerial vehicles to defeat radar tracking,” Journal of Guidance, Control, and Dynamics, vol. 29, no. 2, pp. 279–288, 2006.
- H. Liu, S. Chen, Y. Zhang, and J. Chen, “Modelling radar tracking features and low observability motion planning for UCAV,” in 2012 4th International Conference on Intelligent Human-Machine Systems and Cybernetics, vol. 2, pp. 162–166, Nanchang, China, August 2012.
- T. Robinson, K. Willcox, M. Eldred, and R. Haimes, “Multifidelity optimization for variable-complexity design,” in 11th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Portsmouth, Virginia, September 2006.
- N. M. Alexandrov and R. M. Lewis, “An overview of first-order model management for engineering optimization,” Optimization and Engineering, vol. 2, no. 4, pp. 413–430, 2001.
- S. E. Gano, J. E. Renaud, and B. Sanders, “Hybrid variable fidelity optimization by using a kriging-based scaling function,” AIAA Journal, vol. 43, no. 11, pp. 2422–2433, 2005.
- M. Berci, V. Toropov, R. Hewson, and P. Gaskell, “Metamodelling based on high and low fidelity model interaction for UAV gust performance optimization,” in 50th Aiaa/Asme/Asce/Ahs/Asc Structures, Structural Dynamics, and Materials Conference, Palm Springs, California, May 2009.
- V. Balabanov and G. Venter, “Multi-fidelity optimization with high-fidelity analysis and low-fidelity gradients,” in 10th Aiaa/issmo Multidisciplinary Analysis and Optimization Conference, Albany, New York, August 2004.
- M. G. Hutchison, E. R. Unger, W. H. Mason, B. Grossman, and R. T. Haftka, “Variable-complexity aerodynamic optimization of a high-speed civil transport wing,” Journal of Aircraft, vol. 31, no. 1, pp. 110–116, 1994.
- J. R. R. A. Martins and A. B. Lambe, “Multidisciplinary design optimization: a survey of architectures,” AIAA Journal, vol. 51, no. 9, pp. 2049–2075, 2013.
- S. Doyle, K. Alston, and T. Winter, “Developing the aerodynamics module for the integrated multidisciplinary optimization object system,” in 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, p. 8, Orlando, Florida, January 2011.
- A.-T. Le, K. Gray, and M. L. Baker, “Building the aerodynamics module for the integrated hypersonic aeromechanics tool,” in 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, January 2004.
- D. E. Myers, “Matrix formulation of co-kriging,” Journal of the International Association for Mathematical Geology, vol. 14, no. 3, pp. 249–257, 1982.
- A. I. J. Forrester, A. Sóbester, and A. J. Keane, “Multi-fidelity optimization via surrogate modelling,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 463, no. 2088, pp. 3251–3269, 2007.
- L. Huang, Z. Gao, and D. Zhang, “Research on multi-fidelity aerodynamic optimization methods,” Chinese Journal of Aeronautics, vol. 26, no. 2, pp. 279–286, 2013.
- M. Xiao, G. Zhang, P. Breitkopf, P. Villon, and W. Zhang, “Extended Co-Kriging interpolation method based on multi-fidelity data,” Applied Mathematics and Computation, vol. 323, pp. 120–131, 2018.
- H. S. Chung and J. Alonso, “Design of a low-boom supersonic business jet using cokriging approximation models,” in 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, pp. 2002–5598, Atlanta, Georgia, September 2002.
- B. S. Kim, Y. B. Lee, and D. H. Choi, “Comparison study on the accuracy of metamodeling technique for non-convex functions,” Journal of Mechanical Science and Technology, vol. 23, no. 4, pp. 1175–1181, 2009.
- R. Yuan and G. Bai, “Comparison of neural network and kriging method for creating simulation-optimization metamodels,” in 2009 Eighth IEEE International Conference on Dependable, Autonomic and Secure Computing, pp. 815–821, Chengdu, China, December 2009.
- D. Marcotte, “Cokriging with Matlab,” Computers & Geosciences, vol. 17, no. 9, pp. 1265–1280, 1991.
- J. D. Anderson and J. Wendt, Computational Fluid Dynamics, Springer, 1995.
- H. Yan, J. Chen, and Y. Feng, Aerodynamic Characteristics Analysis and Engineering Calculation of Aircraft, Northwestern Polytechnical University Press, 1990.
- Q. Shi, “The engineering prediction method for aerodynamic characteristics of the symmetrical tactical missile configuration with inlet,” Journal of Experiments in Fluid Mechanics, vol. 18, no. 2, pp. 98–101, 2004.
- D. Zhang, A Course in Computational Fluid Dynamics, Higher Education Press, 2010.
- L. Huang, Z. Gao, and D. Zhang, “Aerodynamic optimizaton based on multi-fidelity surrogate,” ACTA Aerodynamics Sinica, vol. 31, no. 6, pp. 783–788, 2013.
- G. Yuan, Computation of Electromagnetic Characteristics of Complex Targets, Beijing Jiaotong University, 2007.
- J. Fan, Study on the Algorithm of Wideband Scattering by Electrical Large Complex Target, Xidian University, 2015.
- Y. An, D. Wang, and R. Chen, “Improved multilevel physical optics algorithm for fast computation of monostatic radar cross section,” IET Microwaves, Antennas & Propagation, vol. 8, no. 2, pp. 93–98, 2014.
- S. Blume and G. Kahl, “The physical optics radar cross section of an elliptic cone,” IEEE Transactions on Antennas and Propagation, vol. 35, no. 4, pp. 457–460, 1987.
- Y. Wang, Calculation on Radar Cross Section of Complex Targets, Xidian University, 2008.
- Q. Wang, Aircraft Three Dimensional Reconstruction and Stealth Characterstic, Najing University of Aeronautics and Astronautics, 2009.
- J. Yang and R. Chen, “Intergal optimization design for tactical rocket and it’s srm by utilizing the synthesis of internal ballistics and trajectory,” Journal of Astronautics, vol. 20, no. 1, pp. 21–27, 1999.
- G. Fang, “Solid rocket motor overall optimization design,” Journal of Solid Rocket Technology, vol. 1990, no. 2, pp. 1–8, 1990.
- J. Hu, W. Zhang, Z. Xia, and L. Hunag, Ramjet Propulsion Technology, NUDT Press, Changsha, 2013.
- F. Bao, X. Huang, and Z. Zhang, Integral Solid Propellant Ramjet Rocket Motor, China Aetronautic Publishing House, 2006.
- K. Chen, L. Liu, and Y. Meng, Launch Vehicle Flight Dynamics and Guidance, National Defense Industry Press, 2014.
- X. Qian, R. Lin, and Y. Zhao, Missile Flight Mechanics, Beijing Institute of Technology Press, 2008.
- B. Chen, Y. Liu, H. Shen, H. Lei, and Y. Lu, “Performance limitations in trajectory tracking control for air-breathing hypersonic vehicles,” Chinese Journal of Aeronautics, vol. 32, no. 1, pp. 167–175, 2019.
- M. A. Patterson and A. V. Rao, “GPOPS-II: a MATLAB software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming,” ACM Transactions on Mathematical Software, vol. 41, no. 1, pp. 1–37, 2014.
Copyright © 2019 Xiaojian Sun et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.