Performance Evaluation for Launcher Testing Vehicle

: The paper’s purpose is to present a calculus model for a testing vehicle that can be used to validate guidance, navigation and control systems for reusable launchers in all ﬂight phases. The technical solution is based on a throttleable engine with thrust vectoring control and a reaction control system (RCS) used for roll. For calculus, we will develop a nonlinear model with six degrees of freedom, based on quaternion, extended with nonlinear equations that use pulse modulation in order to control roll. In order to synthesize the controller, we also develop a linear model similar to the launcher model. The paper analyzes two basic scenarios, ﬁrst with the ascending and the descending ﬂight phases and the second having a horizontal ﬂight interleaved between ascending and descending ﬂight phases, both scenarios being speciﬁc for reusable launchers. Based on these scenarios, the paper evaluates some performances of the proposed vehicle, namely ﬂight envelope and guidance accuracy.


Introduction
In the field of commercial launchers, nowadays, there are two directions for reducing the costs: achievement of the launcher families with modular launchers like ANGARA [1] and DELTA [2], which allow insertion on different orbits of all sorts of payloads and achievement of the reusable launchers, as it is done in the FALCON program [3,4], developed by Space X. A problem stemming from these trends is the development of an advanced navigation, control and guidance system that can solve the guidance problem for the entire launcher family or make complex maneuvers for the reusable launchers. A Launcher Testing Vehicle (LTV) can be used as an independent stage to solve these complex problems with a low development cost, with the propulsion and control system to test in flight and the launcher guidance navigation and control system (GNC).
In our vision, the Launch Testing Vehicle (LTV) is similar to the vehicle used for landing on a celestial body (Moon, Mars). Discussing technical solutions, the reference landing vehicle Apollo Lunar Module [5,6] uses a descending throttleable engine for landing and an ascent engine for docking. This approach also uses a reaction control system (RCS) having two blocks with four thrusters each for attitude control. An alternative technical solution for a Lunar Lander Demonstrator (LLD) is presented in the paper [7]. LLD has five large descending thrusters and eight small thrusters for attitude control. Many fixed thrusters allow the use of on/off modulation to control roll pitch and yaw without throttling or gimbaling for thrust vectoring control (TVC).
Among the important projects for reusable launchers, we can mention STARSHIP [8], developed by SpaceX Company(Hawthorne, USA), whose technical solution is based on the multiple Raptor Engine with gimbal capability and aerodynamic surfaces (flaps) for control. Another is the CALLISTO project [9], developed in cooperation with Centre National d'Etude Spatiales (CNES), Japan Aerospace Exploration Agency (JAXA), and Deutsches Zentrum fur Luft-und Raumfahrt (DLR), whose technical solution is based on the RV-X engine, developed by JAXA, which will have gimbaling, throttling and re-ignition capability. Moreover, the project supposes the aerodynamic surfaces (fins) and RCS with 6-8 thrusters for control.
Different from these solutions, specific for landing vehicles, in this paper, we will consider a Launch Testing Vehicle (LTV) based on thrust vectoring control (TVC) using a throttleable engine with gimbal capability for pitch and yaw, and separately, two RCS blocks with two thrusters each with Pulse Width -Pulse Frequency modulation (PWPF) [10,11] for roll control. In order to simplify the technical solution, we do not use aerodynamic surfaces, which are inefficient for low velocity.
Currently, a number of current studies [12,13] focus on determining an optimal evolution, in terms of low fuel consumption, by applying convex optimization methods on models with three degrees of freedom (3DOF) for complex evolutions.
Unlike these, the current paper proposes developing a calculation model with six degrees of freedom (6DOF), similar to that of the launchers, to underpin the vehicle's subsequent development and conduct preliminary flight tests. The model will focus on the movement's control and stability, aiming to evaluate the vehicle's performance in terms of the flight envelope and guidance accuracy. For the synthesis of the guidance system, linear models for longitudinal and lateral movement in the local frame based on the Hamilton quaternion are developed. Applying Linear Quadratic Regulator (LQR) synthesis, we will obtain an optimal regulator, and the gains used in control signals will be defined.
For performance evaluation, two simple hypothetical scenarios will be considered, one for vertical evolution, which allows the determination of the maximum altitude, and one for horizontal evolution, which allows the determination of the maximum flight distance, both cases for a fixed amount of fuel. Based on these two scenarios, for two test cases, the landing accuracy will be evaluated in the conditions of the parametric uncertainty of the model as well as the noise introduced by the sensors. Additionally, for a test case in the second scenario, the uniform and turbulent wind influence on lateral movement will be evaluated.
Taking into consideration the work objectives to implement a calculus model for LTV, which will be used to validate GNC for the launcher, the reference frames used are similar to those used for the launcher model described in detail in prior work [14,15]. This is characterized by the y-axis of the local frame oriented vertically upwards and, in the case of Euler angles, the order of rotations 3-2-1 from the local frame to the body frame.
Similar to launchers, because of their shape, the LTV has the aerodynamic focus in front of the center of mass, translating into static instability, which leads to a similar approach regarding the control problem. Different from the launcher, for LTV, some simplifying assumptions related to low velocities and distance domains will be made concerning the frames and the movement equations.
In summary, the novelties proposed by the paper, related to other similar approaches, are: A simplified technical control solution based only on propulsion. The technical solution is similar to lunar landing vehicles. Still, it is customized for the evolution of launchers by considering the aerodynamic terms in the model and the atmospheric, disturbing effects.
Formulation of the control problem based on linear coordinates in the local frame and the quaternion. Unlike our previous paper [16], we introduced the control of linear coordinates to improve guidance, especially landing accuracy, and we used quaternions to ensure the robustness of the model Evaluation of system performance in the flight envelope and landing accuracy with a highlight on the influence of uniform wind, turbulent wind, sensor noise, and uncertainty of the model parameters.
The paper's major contribution, which can be considered a contribution to the state of the art, consists of base movement analysis and linear model, both dedicated for the testing vehicle.
As will be shown during the work, for the technical solution adopted, the basic movement is close to landing on a celestial body despite the presence of the atmosphere.
Regarding the linear form of the motion equations, unlike most works [17] that use the velocity frame for the linear form of translational equations, the local frame will be used in our work, which will facilitate the control of the linear coordinates. Appendix A of the paper will present the method of obtaining the coupled linear form of the equations of motion.
In addition, for the avoidance of singularities and the system's robustness for describing the attitude, including in the linear form, the quaternion was used, which constitutes an additional element regarding the contributions of the work.
For the verification of the model but also for the design of a methodology for evaluating the vehicle's performance, two flight scenarios will be considered in both situations. The test was performed with the control based on the previously specified linear model, which was included in the 6 DOF nonlinear model.

Reference Frames, Rotation Matrix
As we mentioned earlier, the reference frames used will be similar to those used for the micro launcher model. Because LTV has a limited range, up to 10 km around the launch position, some simplifying assumptions will be added to the launcher model. So, we will consider the Earth flat, without rotation; hence we will use only the two frames described below for the equation of motion.
(a) The Start frame/ The Local frame-OX 0 Y 0 Z 0 This coordinate system denoted OX 0 Y 0 Z 0 Figure 1 has its origin on the Earth's surface at the start point of the vehicle. Axis Y 0 is oriented vertically upwards. The frame has X 0 axis orientated to the launch direction. Z 0 axis completes the right frame to the right of the launch plane. This coordinate system is considered to be an inertial reference frame.
The paper's major contribution, which can be considered a contribution to the state of the art, consists of base movement analysis and linear model, both dedicated for the testing vehicle.
As will be shown during the work, for the technical solution adopted, the basic movement is close to landing on a celestial body despite the presence of the atmosphere.
Regarding the linear form of the motion equations, unlike most works [17] that use the velocity frame for the linear form of translational equations, the local frame will be used in our work, which will facilitate the control of the linear coordinates. Appendix A of the paper will present the method of obtaining the coupled linear form of the equations of motion.
In addition, for the avoidance of singularities and the system's robustness for describing the attitude, including in the linear form, the quaternion was used, which constitutes an additional element regarding the contributions of the work.
For the verification of the model but also for the design of a methodology for evaluating the vehicle's performance, two flight scenarios will be considered in both situations. The test was performed with the control based on the previously specified linear model, which was included in the 6 DOF nonlinear model.

Reference Frames, Rotation Matrix
As we mentioned earlier, the reference frames used will be similar to those used for the micro launcher model. Because LTV has a limited range, up to 10 km around the launch position, some simplifying assumptions will be added to the launcher model. So, we will consider the Earth flat, without rotation; hence we will use only the two frames described below for the equation of motion.
(a) The Start frame/ The Local frame-This coordinate system denoted Figure 1 has its origin on the Earth's surface at the start point of the vehicle. Axis is oriented vertically upwards. The frame has axis orientated to the launch direction.
axis completes the right frame to the right of the launch plane. This coordinate system is considered to be an inertial reference frame.   (b) The Body frame-Oxyz This coordinate system originates in the center of mass of the vehicle Figure 1. The axis x is along the longitudinal symmetry axis of the body. Axis y is in a symmetry plane of the vehicle. The axis z arises, forming with the first two axes a right trihedral. Next, we use this trihedral to write dynamic rotation equations around the mass center. Also, this will be used to write thrust and aerodynamic terms.
The rotation between the start frame OX 0 Y 0 Z 0 and body frame Oxyz can be obtained using Euler's angles in order 3-2-1 or using Rodrigues parameters, or using quaternion components, as we can see from Figure 1. Next, for robustness, we will use quaternion components. First, we recall the definition of quaternion components: where l, m, n are the projection of the unit vector of the rotation axis: and σ are the principal rotation angle. The total rotation matrix from the local frame to the body frame according to the paper [10] is:

Translation Equation in Local Frame
Considering the start frame as inertial, applying the impulse theorem, we obtain: where: To the previous equation, we add the translational kinematic equations: Regarding mass variation in time, this is due to the fuel consumption: .
where I sp is the specific impulse of the rocket engine, T 0 is maximum thrust and δ T is throttling command.

Rotational Equations in Body Frame
To obtain rotational dynamic equations, we start with the angular momentum theorem. As we know, the angular momentum theorem shows that a moment applied to a body will change the angular momentum. In this case, taking into consideration that the body frame is a non-inertial frame, we can write: where: h-body angular momentum; M-applied moment; Ω-angular velocity of the body frame relative to the inertial frame (start frame/local frame). The relation can be written in matrix form: where A Ω is the antisymmetric matrix associated with angular velocity vector Ω: J represents the inertial momentum matrix. If the additional hypothesis is made that the vehicle has two symmetry planes, the inertial matrix becomes: where: L,M, N are the components of the applied moment along the body frame axis given by aerodynamics terms and thrust terms: and : denotes the rotation velocity of the body frame relative to the inertial frame expressed by projections along the body frame. Equation (9) can be written in scalar form: According to paper [10], the associated kinematic equations allow obtaining the quaternion components:

Guidance Signals
As we mentioned at the beginning of the paper, because LTV has the purpose of testing the launcher GNC, its structure will be similar to that described in [18][19][20]. The main control signals result from the quaternion error between the current and reference values. Different from the usual control system of a launcher for LTV, we will use error signals between linear coordinates.
Considering the desired attitude quaternion: q 1d q 2d q 3d q 4d T , and the current attitude quaternion q 1 q 2 q 3 q 4 T , the control error quaternion q 1e q 2e q 3e q 4e T can be obtained by using the relation indicated in papers [7,10]: With this error quaternion, as is shown in work [10], we can obtain the guidance signal for LTV evolution: where a i,j are the elements of the rotation matrix A I (3). The signal for velocity control along the body axis which will contribute to throttling command: In relation to (18), we denoted u To ; u no the basic signals corresponding to equilibrium evolution.
The control signals for linear coordinate along the local frame axis u x , u y , u z will be defined next, separately for ascending/descending evolution and for horizontal evolution.
(a) Control signals for linear coordinate during ascending/descending evolution For ascending /descending evolutions, besides the control signals previously defined, we will use the control signals for linear coordinates along the local frame axis: (b) Control signals for linear coordinate during horizontal evolution If we desire to have a horizontal evolution, the control signal for error quaternion and axial velocity is supplemented by a control signal for linear coordinates in the local frame:

Aerodynamic Terms
Due to large incidence angles (±180 • ), determination of aerodynamics terms for the vehicle is difficult. On the other hand, aerodynamical terms are important because they allow the introduction of statical instability of the vehicle due to the position of aerodynamic focus related to the center of mass. Moreover, we will analyze wind influence on vehicle flight through aerodynamic terms and define the base movement for horizontal flight.
In order to solve this problem, we can use the approach proposed in the CALLISTO project [21] by experimental activities in Wind Tunnel. This requires an advanced stage of the project when the vehicle configuration is established. In our preliminary study, we must do some simplified calculus for aerodynamics terms when the vehicle configuration is not precisely definite. Unfortunately, for slender body shapes, the theoretical results are valid only in the case of a small incidence angle [22]. In this case, we will make the following simplifying hypothesis, the vehicle has a spherical shape, and the position of the aerodynamic focus is in front of the mass center.
Based on this hypothesis, the components of the aerodynamic force in the body frame will be considered as drag force opposite to velocity components for each axis: where, due to symmetry, will have: where C d is the drag coefficient for a sphere that depends on the Reynolds number, see Figure 2, according to paper [23]. Considering aerodynamic focus placed at a distance denoted with x F in front of the mass center, the aerodynamic moments are: For LTV configuration, similar to the case of launchers x F > 0.

Thrust Vectoring Control
The guidance signals previously obtained (18) are applied to the thrust vectoring control system (TVC), which can be approximated with a first-order delay element: . δ l τ l = −δ l + u l ; . δ m τ m = −δ m + u m ; . δ n τ n = −δ n + u n ; .
which provides command terms: angular deflections or linear displacements, which are used in the thrust terms described next. Similar to the launcher, the thrust terms are given by thrust vector projection along the body frame, namely by a nozzle or entire engine tilt [24,25]. If we use two tilting angles δ n , δ m similar to the launcher, for pitch and yaw, we can write: where T 0 is the maximum thrust of the engine and δ T is the throttling command with a unitary maximum value. Considering the position of the thrust vector application at the exit of the engine nozzle, we can write the moment commands: where: in which x cm is the current position of the center of the mass from the top of the LTV, and l is the reference length (the distance between the body top and nozzle exit section).
Regarding the roll command, we suppose the existence of a Reaction Control System (RCS), which ensures the roll control through two RCS blocks with two thrusters on each one with monopropellant cold gas or another equivalent propeller, as is shown in [25].

Reaction Control System in Roll
In order to form the roll command, we start from the roll guidance signal u l . This signal is applied to the actuator, which in this case is RCS.
As we can see from Figure 3, there are two blocks of the reactive elements located in the yoz plane of the LTV body, which allows obtaining the roll command with the reactive elements (δ 1 , δ 2 , δ 3 ; δ 4 ).
in which is the current position of the center of the mass from the top of the LTV, and is the reference length (the distance between the body top and nozzle exit section).
Regarding the roll command, we suppose the existence of a Reaction Control System (RCS), which ensures the roll control through two RCS blocks with two thrusters on each one with monopropellant cold gas or another equivalent propeller, as is shown in [25].

Reaction Control System in Roll
In order to form the roll command, we start from the roll guidance signal . This signal is applied to the actuator, which in this case is RCS.
As we can see from Figure 3, there are two blocks of the reactive elements located in the plane of the LTV body, which allows obtaining the roll command with the reactive elements ( , , ; ). To formulate the problem, we will consider the as oriented vectors along nozzles ( Figure 3) with the module in binary form, having a value "1" if the reactive element is working or "0" if the reactive element is not working.
As shown in [10,11,26], the system uses Pulse Width-Pulse Frequency modulation (PWPF) for roll control (Figure 4). PWPF comprises a linear integrator that allows a system supplementary adjustment and a nonlinear element (N), Type Schmidt Trigger [11,25]. To follow the input, a feedback reaction loop is used. The nonlinear element which switches between reactive elements can be equivalent to the relay type with hysteresis and dead zone, as shown in Figure 5. To formulate the problem, we will consider the δ i as oriented vectors along nozzles ( Figure 3) with the module in binary form, having a value "1" if the reactive element is working or "0" if the reactive element is not working.
As shown in [10,11,26], the system uses Pulse Width-Pulse Frequency modulation (PWPF) for roll control (Figure 4). PWPF comprises a linear integrator that allows a system supplementary adjustment and a nonlinear element (N), Type Schmidt Trigger [11,25]. To follow the input, a feedback reaction loop is used.
in which is the current position of the center of the mass from the top of the LTV, and is the reference length (the distance between the body top and nozzle exit section).
Regarding the roll command, we suppose the existence of a Reaction Control System (RCS), which ensures the roll control through two RCS blocks with two thrusters on each one with monopropellant cold gas or another equivalent propeller, as is shown in [25].

Reaction Control System in Roll
In order to form the roll command, we start from the roll guidance signal . This signal is applied to the actuator, which in this case is RCS.
As we can see from Figure 3, there are two blocks of the reactive elements located in the plane of the LTV body, which allows obtaining the roll command with the reactive elements ( , , ; ). To formulate the problem, we will consider the as oriented vectors along nozzles ( Figure 3) with the module in binary form, having a value "1" if the reactive element is working or "0" if the reactive element is not working.
As shown in [10,11,26], the system uses Pulse Width-Pulse Frequency modulation (PWPF) for roll control (Figure 4). PWPF comprises a linear integrator that allows a system supplementary adjustment and a nonlinear element (N), Type Schmidt Trigger [11,25]. To follow the input, a feedback reaction loop is used. The nonlinear element which switches between reactive elements can be equivalent to the relay type with hysteresis and dead zone, as shown in Figure 5. The nonlinear element which switches between reactive elements can be equivalent to the relay type with hysteresis and dead zone, as shown in Figure 5. As shown in Figure 5, the dimension of the dead zone is 2 , hysteresis width is , and the output in saturation is ±1.
For passing from the continuous command to the binary values of the elements we first form the intermediate signals: The signal is integrated, obtaining the signal : The signal , is passed through the nonlinear element shown in Figure 5, obtaining signal Using the reactive elements, we obtain the roll equivalent command: We can observe that, due to the binary structure of the elements , the equivalent roll command can have values "±1" or "0". Finally, using the Fourier transform [27], we can obtain a linear differential equation: similar to relation (25). Given the previous reactive force of the RCS element and body diameter , we can obtain the roll command moment given by:

Base Movement, Linear Form of the Equation of Motion, Stability and Command Matrices
To define the base movement, we start from the dynamic equations for translation (4) and rotation (9). Due to the initial hypothesis for the equation of motion, the Earth is considered a flat surface that does not rotate; the gravity acceleration does not depend on the orientation of the local frame, which can be arbitrarily chosen. In this case, we can suppose, without reducing the problem generality, that the base movement takes place in As shown in Figure 5, the dimension of the dead zone is 2a, hysteresis width is b, and the output in saturation is ±1.
For passing from the continuous command u l to the binary values of the elements δ i we first form the intermediate signals: The signal u is integrated, obtaining the signal x: .
The signal x, is passed through the nonlinear element shown in Figure 5, obtaining signal y y = N(x), where we denoted with N(x) the nonlinear element applied to x signal in the open loop: Using the reactive elements, we obtain the roll equivalent command: We can observe that, due to the binary structure of the elements δ i , the equivalent roll command δ l can have values "±1" or "0".
Finally, using the Fourier transform [27], we can obtain a linear differential equation: . δ l τ l ∼ = −δ l + u l (34) similar to relation (25). Given the previous reactive force of the RCS element R and body diameter d, we can obtain the roll command moment given by:

Base Movement, Linear Form of the Equation of Motion, Stability and Command Matrices
To define the base movement, we start from the dynamic equations for translation (4) and rotation (9). Due to the initial hypothesis for the equation of motion, the Earth is considered a flat surface that does not rotate; the gravity acceleration does not depend on the orientation of the local frame, which can be arbitrarily chosen. In this case, we can suppose, without reducing the problem generality, that the base movement takes place in a vertical plane X 0 OY 0 of a local frame, the lateral motion being null. Simultaneously, we suppose a configuration with two symmetry planes adopted for the inertia matrix in dynamic rotation equations, and we can choose the roll angle null. That means in the base movement; we have only q 3 as an independent parameter, the other quaternions components being: In this case, the rotation matrix (5) becomes: From the translation equation, we will consider the relations in the vertical plane: where we denoted G = mg, and from rotation equations, we hold the pitch equation: If we cancel the velocity derivatives, it results: and from the pitch equation, we cancel the angular acceleration, obtaining: The nonlinear system (40), (41) can be solved with Newton's iterative algorithm. If we impose as a base movement the horizontal flight at a specified altitude and velocity, from Equations (40) and (41), we obtain the component q 3 , and thrust components X T , Y T , the thrust moment in pitch N T being linked to the component Y T of the thrust given by Equation (27).
From thrust components, we can obtain the throttling command δ T and the angular pitch thrust command δ n , which means the equilibrium commands: Taking into consideration that the actuators in Equation (25) have unitary gain, the guidance signals equilibrium commands u T0 and u n0 from the relation (18) can be obtained with the relations (42).
Having base movement defined, considering as commands: throttling command δ T the angular deflections δ m , δ n , and equivalent command δ l the stability matrix and the command matrix can be obtained. For cross-checking, the stability and command matrices were determined both numerically and by analytical relations, obtaining identical results. For clarity, next, we present their analytical expressions.
Because in the base movement, the longitudinal motion is decupled from lateral motion, we will separately analyze these two motions.
Starting from the coupled linear development presented in Appendix A, for longitudinal motion, having the states Vx, Vy, r, q 3 , x, y and the commands δ T , δ n , the linear form of the equations are: where K = S 2 ρC d and the drag dependence of the velocity and the altitude is given by: Denoting: we obtain for the longitudinal motion the stability matrix in Table 1 and the command matrix in Table 2.
Because in longitudinal motion, the controlled states differ from the ascending/descending evolutions to the horizontal evolution, we will define a separate stability matrix and command matrix for each type of longitudinal evolution.
For longitudinal motion in ascending/descending evolution, the states are Vx, Vy, r, q 3 , x, and the commands are δ T , δ n . Because in this particular case, the basic movement has supplementary conditions: the stability matrix has the form from Table 3, Table 3. Stability matrix for longitudinal motion in ascending/descending evolution A. and the command matrix has the form from Table 4. Table 4. Command matrix for longitudinal motion in ascending/descending evolution B.
where: For longitudinal motion in horizontal evolution, the states are Vx, Vy, r, q 3 , y and the commands are δ T , δ n . Because in this particular case, the basic movement has the supplementary condition: The stability matrix has the form from Table 5, and the command matrix has the form from Table 6. Table 6. Command matrix for longitudinal motion in horizontal evolution B. where: Starting from the coupled linear development presented in Appendix A, for lateral motion, the states are V z , p, q, q 1 , q 2 , z and the commands are δ l , δ m the linear form of the equations are: . ∆q 1 = q 4 2 ∆p − q 3 2 ∆q; . ∆q 2 = q 3 2 ∆p + q 4 2 ∆q; . ∆z 0 = ∆V z . Denoting: we obtain for the lateral motion the stability matrix shown in Table 7   Table 7. Stability matrix for lateral motion A.
and the command matrix shown in Table 8.
For lateral motion in ascending/descending evolution, considering condition (45), the matrix elements become: For lateral motion in horizontal evolution, considering condition (47), the matrix elements become: (52)

Linear Form of the Guidance Signals, the Regulator's Design
The linearization of guidance signals will be done related to the system's states: Velocity in local frame: V x ; V y ; V z ; Coordinate in local frame: x 0 ; y o ; z 0 ; Angular velocity in body frame: p; q; r; Attitude (quaternion): q 1 ; q 2 ; q 3 . First, we linearized the error quaternion given by relation (17). Considering for the basic movement, the desired quaternion is equal to the basic quaternion, the relation (17) becomes: Separating the deviation of the desired quaternion as an input signal and considering: we obtain: where the input signals are: Next, we linearized the axial velocity control (19). Considering the basic movement, the linear form of the relation (19) becomes: where ∆ f u is an input signal containing desired axial velocities. On the other hand, for basic movement, the elements of the rotation matrix (3) become: This means that the relationship (57) becomes: Having the quaternion-related functions linearized, we will look for the linear form of the guidance signal for linear coordinates in different evolutions.
For ascending/descending evolution, the linear coordinate signal in local frame (20) becomes: where ∆ f z ; ∆ f y is the input signal containing the desired linear coordinates. Then, for q 3 = √ 2/2 corresponding base movement for θ = 90 • , we add axial velocity control (59): For horizontal evolution, the linear coordinate in the local frame (21) became: Then, for V y = 0 we add axial velocity control (58), (59): If we neglected the time delay introduced by actuator Equations (25) and (34), we obtain a direct link between the linear form of the guidance signals and commands: In this case, starting from the scalar form of the guidance signals (18), we can write the linear form of the command components directly, separately, for the two motions.
The command for longitudinal motion in ascending/descending evolution becomes: Denoting: k we obtain the regulator matrix indicated in Table 9. Table 9. The regulator matrix for longitudinal motion in ascending/descending evolution K.
The command for longitudinal motion in horizontal evolution becomes: Denoting: we obtain the regulator matrix shown in Table 10. The commands for lateral motion in all evolutions are: Denoting: the regulator matrix for lateral motion is shown in Table 11. Table 11. The regulator matrix for lateral motion K.
To obtain the matrix terms, as shown in work [28], we can apply an LQR synthesis to an optimal regulator. The problem implies solving an algebraic Riccati equation.
Finally, the control coefficients introduced in control signals relations (18)-(21) can be obtained from the regulator's elements: For longitudinal motion in ascending/descending evolution, k 3 , k 8 , k 9 , are obtained directly, while the others are given by: For longitudinal motion in horizontal evolution, k 3 is obtained directly, while the others are given by: For lateral motion, k 1 , k 5 , k 8 , k 9 are obtained directly, while the others are given by:

Model Parameters
To build a performance evaluation methodology, including the quality of the guidance system, it is necessary to define a realistic LTV model. Since the data for such test vehicles is uncertain, we will refer to the information for classical launcher stages. For this purpose, we consider the Lox-Kerosene pair as fuel for the rocket engine, which is a common solution found in a series of classic (Soyuz) but also current launchers (Atlas V, Antares, Falcon 1, Falcon 9) If we adopt an initial take-off mass, the ratio between the structural mass and the mass at the start is defined by the structural coefficient [29]: where m s is the structural mass and m p is the propellant mass. As shown in [29], this coefficient measures the vehicle designer's skill and evolved over time for different launchers. Thus, if for the stages of the Soyuz launcher, the structural coefficient is around 9% for Falcon 1, it becomes 6%, and for Falcon 9, it becomes 3%, which shows the technological evolution of the launchers. In order to have a realistic model for LTV, we will consider a structural coefficient of 30%, which is much oversized and will allow the introduction of additional control and testing elements, specific to a test vehicle, into the structural mass.
Based on the previous statements, we will consider a hypothetical LTV with an initial mass of 100 kg, from which propellant is 70 kg, equipped with a liquid rocket engine with a specific impulse I sp = 230[s], specific for the regular liquid engine (Kerosene + LOX), with maximum thrust T 0 = 1200 [N] and reactive force of one RCS element R = 5 [N].
The vehicle has the shape and the dimension presented in Figure 6. From a dimensional point of view, we can check if the volume of the imposed fuel is compatible with the dimensions indicated in Figure 6.
For this, we start from the density of LOX of 1141 kg/m 3 and Kerosene of 825 kg/m 3 . Considering the LOX/Kerosene mixture ratio = 2.72, we obtain an average propellant density of 1056 kg/m 3 , which for a mass of 70 kg leads to a hypothetical spherical tank with a radius of 0.25 m, a size compatible with the geometric dimensions of the vehicle indicated in Figure 6.
For the considered configuration, the reference length is l = 1 m, the reference surface S = 0.283 m 2 and the position of aerodynamic focus is x F = 0.4 m.
Based on the previous statements, we will consider a hypothetical LTV with an initial mass of 100 kg, from which propellant is 70 kg, equipped with a liquid rocket engine with a specific impulse = 230[ ], specific for the regular liquid engine (Kerosene + LOX), with maximum thrust = 1200 [ ] and reactive force of one RCS element = 5 [ ]. The vehicle has the shape and the dimension presented in Figure 6.  Figure 6.
For this, we start from the density of LOX of 1141 kg/m 3 and Kerosene of 825 kg/m 3 . Considering the LOX/Kerosene mixture ratio = 2.72, we obtain an average propellant density of 1056 kg/m 3 , which for a mass of 70 kg leads to a hypothetical spherical tank with a radius of 0.25 m, a size compatible with the geometric dimensions of the vehicle indicated in Figure 6 For the considered configuration, the reference length is = 1 , the reference surface = 0.283 and the position of aerodynamic focus is = 0.4 .

(a) Mechanical data
For mechanical data, we will consider two cases. The first will be the LTV at the start, with full fuel, and the second at the end of evolution, without fuel. Because the fuel consumption defined through Equation (7) depends on flight conditions, the intermediate values of the other mechanical data will be obtained through interpolation as a function of LTV mass. The mechanical characteristics of LTV are shown in Table 12.

(b) Time constants and controller gains
The time constants for actuators are: Model parameters defined will be used for the development of subsequent applications. For the evaluation of the guidance precision, an uncertainty of these parameters will be considered, which will lead to a dispersion of the evaluated trajectories, including the impact point. For mechanical data, we will consider two cases. The first will be the LTV at the start, with full fuel, and the second at the end of evolution, without fuel. Because the fuel consumption defined through Equation (7) depends on flight conditions, the intermediate values of the other mechanical data will be obtained through interpolation as a function of LTV mass. The mechanical characteristics of LTV are shown in Table 12.

(b) Time constants and controller gains
The time constants for actuators are: τ T = 0.5 s; τ n = 0.1 s; τ l = 0.1 s; τ m = 0.1 s and the gains used in control signals are: k 1 = 3.9; k 2 = 2.5; k 3 = 2.0; k 4 = 5.5; k 5 = 2.2; k 6 = 6.9; k 7 = 0.32; k 8 = 0.7; k 9 = 1. Model parameters defined will be used for the development of subsequent applications. For the evaluation of the guidance precision, an uncertainty of these parameters will be considered, which will lead to a dispersion of the evaluated trajectories, including the impact point.

Flight Scenarios
Although in some papers [12,13], the optimal trajectories for reaching a desired final position are analyzed, in the initial phase of LTV testing, it is necessary to define some simple trajectories to verify the vehicle's controllability as its field of use by defining the flight envelope.
For this purpose, we will further define two flight scenarios, one vertical ascent/descent type, and the second having a horizontal evolution interspersed between the vertical evolutions of ascent/descent.
In the first scenario, LTV starts at an altitude y 0 = 0.1 [m], with an initial velocity V y0 = 0.1 ms −1 and initial pitch angle: θ 0 = 90 • . Next, the ascending phase follows, which lasts until LTV achieves the desired altitude y d , evolving at an imposed pitch angle θ d1 = 90 • and imposed ascensional velocity V yd1 = V. After ascending phase, the descending phase follows with imposed pitch angle θ d2 = 90 • , evolution being vertical with descending constant velocity V yd2 = −V, which lasts until it reaches the breaking altitude y f . After y f braking occurs, and the speed decreases until V y4 = 1 ms −1 ensuring a smooth landing for the LTV.
In the second scenario, LTV starts at an altitude y 0 = 0.1 [m], with an initial velocity V y0 = 0.1 ms −1 and initial pitch angle: θ 0 = 90 • . Next, the ascending phase follows, which lasts until LTV achieves the desired altitude y d evolving at an imposed pitch angle θ d1 = 90 • and imposed ascensional velocity V yd1 = V. After the ascending phase, the horizontal phase follows at an imposed altitude y d , with V yd2 = 0, until LTV achieves the desired abscissa x d , with imposed velocity V xd = V. After the horizontal phase, the descending phase follows with imposed pitch angle θ d2 = 90 • , evolution being vertical with descending constant velocity V yd3 = −V which lasts until it reaches the breaking altitude y f . After y f braking occurs, and the speed decreases until V y4 = 1 ms −1 ensuring a smooth landing for the LTV.
To synthesize the flight evaluation, we chose to have the same module of velocity V in the entire flight evolution, except for the breaking phase when the value is 1 ms −1 . The value of flight velocity V will be particularized for the test case.

(a) First test case
In order to exemplify the first scenario, we will consider an evolution with V = 10 ms −1 , with the desired altitude y d = 400 [m]. Figure 7 presents the vertical trajectory obtained by LTV for this first test case. The trajectory contains two flight phases: ascending phase and descending phase.    Figure 8 shows the equilibrium solutions of Equations (40) and (41) (q3_e, dt_e, dnt_e) compared to the solutions obtained through differential equations integration for the 6DOF model in the first scenario (q3, dt, dnt). Figure 9 shows eigenvalues of the stability matrix A for longitudinal motion in ascending/descending evolution corresponding to Table 3. We can observe a pair of complex eigenvalues with real part positive due to static instability (position aerodynamics focus in front of the mass center). The variation of the eigenvalues corresponds with the variation of the flight parameters in the first scenario.    Figure 9 shows eigenvalues of the stability matrix for longitudinal motion in ascending/descending evolution corresponding to Table 3. We can observe a pair of complex eigenvalues with real part positive due to static instability (position aerodynamics focus in front of the mass center). The variation of the eigenvalues corresponds with the variation of the flight parameters in the first scenario.  Figure 10 shows the eigenvalues of the regulated matrix A-BK for longitudinal motion in ascending/descending evolution. Considering the variation of the eigenvalues for the open loop system in Figure 9, one can observe that the real part of all eigenvalues for the closed-loop system is negative, resulting in stability for the first scenario evolution. Figure 9. Eigenvalues of the stability matrix A for longitudinal motion in ascending/descending evolution. Figure 10 shows the eigenvalues of the regulated matrix A − BK for longitudinal motion in ascending/descending evolution. Considering the variation of the eigenvalues for the open loop system in Figure 9, one can observe that the real part of all eigenvalues for the closed-loop system is negative, resulting in stability for the first scenario evolution. Figure 11 shows the eigenvalues of the stability matrix A for the uncontrolled lateral motion corresponding to Table 7. We can observe a pair of complex eigenvalues with real part positive due to static instability (position aerodynamics focus in front of the mass center). The variation of the eigenvalues corresponds to the variation of the flight parameters in the first scenario. Due to the symmetry of the configuration and evolution, the shown eigenvalues of the stability matrix for the lateral motion correspond to the eigenvalues of the stability matrix A for longitudinal motion presented in Figure 9. The difference consists in the number of eigenvalues. For lateral motion, we have six values; for longitudinal, we only have five values. Figure 9. Eigenvalues of the stability matrix A for longitudinal motion in ascending/descending evolution. Figure 10 shows the eigenvalues of the regulated matrix A-BK for longitudinal motion in ascending/descending evolution. Considering the variation of the eigenvalues for the open loop system in Figure 9, one can observe that the real part of all eigenvalues for the closed-loop system is negative, resulting in stability for the first scenario evolution. Figure 10. Eigenvalues of the matrix A-BK for longitudinal motion in ascending/descending evolution. Figure 11 shows the eigenvalues of the stability matrix for the uncontrolled lateral motion corresponding to Table 7. We can observe a pair of complex eigenvalues with real part positive due to static instability (position aerodynamics focus in front of the mass center). The variation of the eigenvalues corresponds to the variation of the flight parameters in the first scenario. Due to the symmetry of the configuration and evolution, the shown eigenvalues of the stability matrix for the lateral motion correspond to the eigenvalues of the stability matrix for longitudinal motion presented in Figure 9. The   Figure 12 shows eigenvalues of the regulated matrix A − BK for lateral motion in all evolutions. One can observe that the real part of all eigenvalues has real negative parts, implying the stability of this evolution. Due to the symmetry of the configuration and evolution, the shown eigenvalues of the regulated matrix for lateral motion correspond to the eigenvalues of the regulated matrix A − BK for longitudinal motion presented in Figure 10. The difference consists in the number of the eigenvalues. For lateral motion, there are six values, while for longitudinal only five.

Im
From Figure 13, we can observe that vertical velocity Vy increases during ascending phase to 10 m/s and becomes negative (−10 m/s) during the descending phase and finally has the value of −1 m/s during the breaking phase. Moreover, we can observe the vertical coordinate y that increases from 0.1 m to 400 m and then decreases to zero.
From Figure 14, we observe that the thrust throttling command dT decreases during ascending and descending phases, following the corresponding decrease in LTV mass. The command reaches peak values during transition phases. Figure 12 shows eigenvalues of the regulated matrix − for lateral motion in all evolutions. One can observe that the real part of all eigenvalues has real negative parts, implying the stability of this evolution. Due to the symmetry of the configuration and evolution, the shown eigenvalues of the regulated matrix for lateral motion correspond to the eigenvalues of the regulated matrix − for longitudinal motion presented in   From Figure 13, we can observe that vertical velocity increases during ascending phase to 10 m/s and becomes negative (−10 m/s) during the descending phase and finally has the value of −1 m/s during the breaking phase. Moreover, we can observe the vertical coordinate that increases from 0.1 m to 400 m and then decreases to zero. From Figure 14, we observe that the thrust throttling command dT decreases during ascending and descending phases, following the corresponding decrease in LTV mass. The command reaches peak values during transition phases. In terms of lateral and rolling movement, the values of states and commands are null during all flight phases.

Vy[m/s] y [m]
The lateral and roll channels have no state or command variations (their derivatives remain null) during the first flight scenario.

(b) Second test case
In order to exemplify the second scenario, we will consider an evolution with V = 7 [m/s], at altitude y d = 300 m with the imposed distance of the horizontal flight x d = 1000 m. Figure 15 presents the vertical trajectory obtained by LTV. We can observe three flight phases: ascending, horizontal, and descending. We can also observe that the flight attitude for LTV in all three phases is with the pitch angle close to 90 degrees (q 3 = 0.707) corresponding to the basic movement shown in Figure 16. From Figure 14, we observe that the thrust throttling command dT decreases during ascending and descending phases, following the corresponding decrease in LTV mass. The command reaches peak values during transition phases.  In order to exemplify the second scenario, we will consider an evolution with = 7 [m/s], at altitude = 300 m with the imposed distance of the horizontal flight = 1000 m. Figure 15 presents the vertical trajectory obtained by LTV. We can observe three flight phases: ascending, horizontal, and descending. We can also observe that the flight attitude for LTV in all three phases is with the pitch angle close to 90 degrees ( = 0.707) corresponding to the basic movement shown in Figure 16.   In order to exemplify the second scenario, we will consider an evolution with = 7 [m/s], at altitude = 300 m with the imposed distance of the horizontal flight = 1000 m. Figure 15 presents the vertical trajectory obtained by LTV. We can observe three flight phases: ascending, horizontal, and descending. We can also observe that the flight attitude for LTV in all three phases is with the pitch angle close to 90 degrees ( = 0.707) corresponding to the basic movement shown in Figure 16.    scenario (q3, dt, dnt). Figure 17 shows eigenvalues of the stability matrix for longitudinal motion in hor-  Figure 16 shows the equilibrium solutions of Equations (40) and (41) (q3_e, dt_e, dnt_e) compared to the solution obtained through the integration of differential equations for the 6 DOF model in the second scenario (q3, dt, dnt). Figure 17 shows eigenvalues of the stability matrix A for longitudinal motion in horizontal evolution corresponding to Table 5. We can observe a pair of complex eigenvalues with real part positive due to static instability (position of the aerodynamics focus in front of the mass center). The variation of the eigenvalues corresponds to the variation of the flight parameters in horizontal evolution in the second scenario.   Figure 17, one can observe that the real part of all eigenvalues for the closed-loop system is negative, resulting in stability for the entire evolution for the second scenario.   Figure 17, one can observe that the real part of all eigenvalues for the closed-loop system is negative, resulting in stability for the entire evolution for the second scenario.   Figure 17, one can observe that the real part of all eigenvalues for the closed-loop system is negative, resulting in stability for the entire evolution for the second scenario.  Figure 19 shows vertical velocity and vertical coordinate (altitude). We can observe that vertical velocity Vy increases during ascending phase until 7 m/s and is null in horizontal evolution, then becomes negative V y = −7 m/s during the descending phase and finally has the value V y = −1 m/s during the breaking phase. Moreover, we can observe the vertical coordinate y that increases from 0.  From Figure 20, we can observe that horizontal velocity Vx, which is zero during the ascending phase, has a value of 7 [m/s] in horizontal evolution and then becomes zero during the descending and breaking phases. Moreover, we can observe the abscissa x that increases from 0 m until 1000 m.  From Figure 20, we can observe that horizontal velocity Vx, which is zero during the ascending phase, has a value of 7 [m/s] in horizontal evolution and then becomes zero during the descending and breaking phases. Moreover, we can observe the abscissa x that increases from 0 m until 1000 m.  From Figure 20, we can observe that horizontal velocity Vx, which is zero during the ascending phase, has a value of 7 [m/s] in horizontal evolution and then becomes zero during the descending and breaking phases. Moreover, we can observe the abscissa x that increases from 0 m until 1000 m.   Figure 21 shows the velocity components in the local frame. We can observe that vertical velocity V y increases during the ascending phase, null in the horizontal evolution and becomes negative during the descending phase. The horizontal velocity Vx is null during the ascending phase, becomes constant during the horizontal phase, and null during the descending phase.  Even though a quaternion approach is used in the 6DOF model, we prefer to present the results using Euler angles for a better understanding. Figure 22 shows the desired pitch angle (ted) and the achieved pitch angle (te). We can observe that in ascending phase, the pitch angle has the value θ d1 = 90°, during horizontal evolution, follows the value of the equilibrium pitch angle, which becomes smaller while the mass decreases, and in descending phase, it takes the value θ d2 = 90°. Figure 23 shows a detail of Figure 22.   Even though a quaternion approach is used in the 6DOF model, we prefer to present the results using Euler angles for a better understanding. Figure 22 shows the desired pitch angle (ted) and the achieved pitch angle (te). We can observe that in ascending phase, the pitch angle has the value θ d1 = 90 • , during horizontal evolution, follows the value of the equilibrium pitch angle, which becomes smaller while the mass decreases, and in descending phase, it takes the value θ d2 = 90 • . Figure 23 shows a detail of Figure 22.  Even though a quaternion approach is used in the 6DOF model, we prefer to present the results using Euler angles for a better understanding. Figure 22 shows the desired pitch angle (ted) and the achieved pitch angle (te). We can observe that in ascending phase, the pitch angle has the value θ d1 = 90°, during horizontal evolution, follows the value of the equilibrium pitch angle, which becomes smaller while the mass decreases, and in descending phase, it takes the value θ d2 = 90°. Figure 23 shows a detail of Figure 22.    Figure 24 shows the longitudinal commands. We observe that thrust throttling command dT has maximum value during the ascending phase, after which it decreases continuously during horizontal and descending evolutions. Pitch deflection command dn is null during the ascending phase, becomes negative during horizontal evolution, and is null during descending evolution.  Figure 25 shows the incidence angle alfa and the pitch angle te. We can observe that incidence reaches the value of the pitch angle after the ascending phase, then maintains a value equal to the value of the pitch angle during horizontal evolution, and during the descending phase, it gets close to 180 degrees due to vertical descending evolution. The incidence angles and related pitch angle values prove that vehicle maintains a vertical attitude with the tip up in all flight phases.  Figure 24 shows the longitudinal commands. We observe that thrust throttling command dT has maximum value during the ascending phase, after which it decreases continuously during horizontal and descending evolutions. Pitch deflection command dn is null during the ascending phase, becomes negative during horizontal evolution, and is null during descending evolution.  Figure 24 shows the longitudinal commands. We observe that thrust throttling command dT has maximum value during the ascending phase, after which it decreases continuously during horizontal and descending evolutions. Pitch deflection command dn is null during the ascending phase, becomes negative during horizontal evolution, and is null during descending evolution.  Figure 25 shows the incidence angle alfa and the pitch angle te. We can observe that incidence reaches the value of the pitch angle after the ascending phase, then maintains a value equal to the value of the pitch angle during horizontal evolution, and during the descending phase, it gets close to 180 degrees due to vertical descending evolution. The incidence angles and related pitch angle values prove that vehicle maintains a vertical attitude with the tip up in all flight phases.  Figure 25 shows the incidence angle alfa and the pitch angle te. We can observe that incidence reaches the value of the pitch angle after the ascending phase, then maintains a value equal to the value of the pitch angle during horizontal evolution, and during the descending phase, it gets close to 180 degrees due to vertical descending evolution. The incidence angles and related pitch angle values prove that vehicle maintains a vertical attitude with the tip up in all flight phases. In terms of lateral and roll motion, the values of states and commands are null during all the flight phases.
During the second flight scenario, the lateral and roll channels have no state or command variations (their derivatives remain null).

Flight Envelopes
We will evaluate the vehicle's performance based on the previously presented scenarios. According to [19], the performance of the guided vehicle means the flight envelope and the guidance precision. Next, we will evaluate the flight envelope using the previous described scenarios by determining maximum altitude and distance. Then, considering the uncertainty of the model parameters and sensor noise, we will evaluate the trajectory dispersions.

(a) Maximum altitude
Considering the first scenario, the vertical evolution, and using different ascending/descending velocities, we can obtain the maximum altitude of the vehicle as a velocity function, presented in Figure 26.
From Figure 26, we can observe that maximum altitude increases with velocity until it reaches a peak after that decreases. The maximum altitude obtained by the vehicle was   Figure 25. Diagram of the pitch angle (te) and the incidence angle (alfa).
In terms of lateral and roll motion, the values of states and commands are null during all the flight phases.
During the second flight scenario, the lateral and roll channels have no state or command variations (their derivatives remain null).

Flight Envelopes
We will evaluate the vehicle's performance based on the previously presented scenarios. According to [19], the performance of the guided vehicle means the flight envelope and the guidance precision. Next, we will evaluate the flight envelope using the previous described scenarios by determining maximum altitude and distance. Then, considering the uncertainty of the model parameters and sensor noise, we will evaluate the trajectory dispersions.

(a) Maximum altitude
Considering the first scenario, the vertical evolution, and using different ascending/descending velocities, we can obtain the maximum altitude of the vehicle as a velocity function, presented in Figure 26.

(b) Maximum distance
Considering the second scenario, with horizontal evolution, and using different velocities and different altitudes for the horizontal flight, we can obtain the maximum distance of the vehicle as a velocity and altitude function is presented in Figure 27. From Figure 26, we can observe that maximum altitude increases with velocity until it reaches a peak after that decreases. The maximum altitude obtained by the vehicle was Hmax = 2.4 km, corresponding to V = 21 [m/s]. Moreover, from Figure 26, we can observe three restrictions that limit the flight envelope: fuel consumption; exceeding maximum acceleration during maneuver a max = 30 m/s 2 ; exceeding final velocity (V f inal max = 2 [m/s]).

(b) Maximum distance
Considering the second scenario, with horizontal evolution, and using different velocities and different altitudes for the horizontal flight, we can obtain the maximum distance of the vehicle as a velocity and altitude function is presented in Figure 27.

(b) Maximum distance
Considering the second scenario, with horizontal evolution, and using different velocities and different altitudes for the horizontal flight, we can obtain the maximum distance of the vehicle as a velocity and altitude function is presented in Figure 27.

Wind Influence
For wind influence, we consider the second scenario, uniform wind and turbulence, as described in [17]. First, we consider lateral wind, with the velocity Wz = 1 m/s, the influence of lateral parameters being presented in the next figure ( Figure 28).
Further, we consider a layer of turbulence between the altitude of 200-250 m. This layer of turbulence influences the motion in the ascending and descending phases, as seen in Figure 29.
The results prove that the vehicle is stable when considering the wind influence, lateral deviation in z coordinate, and Vz velocity in both cases being insignificant. exceedance of maximum acceleration ( = 30 [m/s ]), the same restriction as in the case of maximum altitude.

Wind Influence
For wind influence, we consider the second scenario, uniform wind and turbulence, as described in [17]. First, we consider lateral wind, with the velocity Wz = 1 m/s, the influence of lateral parameters being presented in the next figure ( Figure 28). Further, we consider a layer of turbulence between the altitude of 200-250 m. This layer of turbulence influences the motion in the ascending and descending phases, as seen in Figure 29. For wind influence, we consider the second scenario, uniform wind and turbulence, as described in [17]. First, we consider lateral wind, with the velocity Wz = 1 m/s, the influence of lateral parameters being presented in the next figure (Figure 28). Further, we consider a layer of turbulence between the altitude of 200-250 m. This layer of turbulence influences the motion in the ascending and descending phases, as seen in Figure 29.

Vz[m/s] z [m]
Medium intensity turbulence Figure 29. Influence of turbulence on the lateral coordinate (z) and lateral velocity (Vz).

Flight Parameters Dispersion
To evaluate dispersion of flight parameters we consider the dispersion with normal distribution of some uncertain model parameter, defined in Section 2.5 and Section 5.1, as a percentage of the nominal value, as follows: Dispersions of aerodynamic drag coefficient and position of the aerodynamic focus related reference length: Dispersion of mass and inertial moments: For LTV, the following values for standard deviations were considered: The uncertainty of the model parameters described previously are considered input data and do not change over time during the flight.
Supplementary an additive noise affecting sensors measurement having uniform distribution was considered: sensors for angular attitude ±1/60 deg; sensors for linear position ±0.3 m.
The additive noise of the sensor will be variable in time but maintain uniform distribution.
To obtain the trajectory dispersion, we ran 500 tests for each flight case presented previously, the values for uncertain sizes and noise for the sensor being obtained by random number generators [30]. Figure 30 shows a trajectory beam in a vertical plane corresponding to the first case.
The uncertainty of the model parameters described previously are considered input data and do not change over time during the flight.
Supplementary an additive noise affecting sensors measurement having uniform distribution was considered: sensors for angular attitude ± 1 60 ⁄ ; sensors for linear position ±0.3 m.
The additive noise of the sensor will be variable in time but maintain uniform distribution.
To obtain the trajectory dispersion, we ran 500 tests for each flight case presented previously, the values for uncertain sizes and noise for the sensor being obtained by random number generators [30]. Figure 30 shows a trajectory beam in a vertical plane corresponding to the first case.   Synthesizing the results for landing positions in Table 13, we can observe that the final coordinate average does not exceed the desired values (X = 0, Z = 0) with more than a few centimeters, and standard deviations are also small values. Regarding final velocity, this is very close to the desired values (V = 1 m/s), and the standard deviation has a small value (12 cm/s).  Figure 32 shows a trajectory beam in a vertical plane corresponding to the second case. Synthesizing the results for landing positions in Table 13, we can observe that the final coordinate average does not exceed the desired values (X = 0, Z = 0) with more than a few centimeters, and standard deviations are also small values. Regarding final velocity, this is very close to the desired values (V = 1 m/s), and the standard deviation has a small value (12 cm/s).       Synthesizing the results for landing positions (Tables 13 and 14), we can observe that the average of the final coordinate does not exceed the desired values (X = 500 m, Z = 0 m) with more than a few centimeters, and standard deviations are also small values. Regarding final velocity, this is very close to the desired value (V = 1 m/s), and the standard deviation has a very small value (0.8 cm/s).

Conclusions
In order to build the equation of motion for the testing vehicle, we use two frames. The first one is the Local Frame/Start Frame, which allows us to write the translational equations; the second is the Body Frame, which allows us to write rotational dynamic equations. In order to correlate the model and the result of the testing vehicle with the Synthesizing the results for landing positions (Tables 13 and 14), we can observe that the average of the final coordinate does not exceed the desired values (X = 500 m, Z = 0 m) with more than a few centimeters, and standard deviations are also small values. Regarding final velocity, this is very close to the desired value (V = 1 m/s), and the standard deviation has a very small value (0.8 cm/s).

Conclusions
In order to build the equation of motion for the testing vehicle, we use two frames. The first one is the Local Frame/Start Frame, which allows us to write the translational equations; the second is the Body Frame, which allows us to write rotational dynamic equations. In order to correlate the model and the result of the testing vehicle with the model and results of the launcher, we use the Start Frame with the y-axis up and quaternion Hamilton. Further, we write the relations that describe the guidance command, which allows for building a 6DOF guided model. To complete the definition, aerodynamics, and thrust terms were introduced. Starting from the 6DOF nonlinear model, we obtained the basic movement, the linear form of the motion equation, and the stability and command matrices. From the analysis of the basic movement, it was found that for the technical solution adopted, the flight attitude is always vertical. Figure 15 shows that it leads to horizontal evolution with a high incidence angle, close to 90 degrees. Using the model, two flight scenarios were evaluated. The first one contains only ascending and descending evolution, and the second contains three phases: ascending, horizontal, and descending. For each scenario, a flight case was defined and analyzed, and the flight envelope of the testing vehicle was defined. Supplementary, the influence of the uncertainty of the model parameters and sensor noise on trajectory dispersion was evaluated. Moreover, the influence of uniform wind and turbulence on lateral trajectory deviation was analyzed. The model developed must be improved by using experimental measurements of the dynamic regime's aerodynamic terms and thrust characteristics. After completing ground measurements, the 6DOF model can be used for flight experiments design. In order to achieve the control system, it is necessary to use the inertial measurement unit (IMU) combined with GPS measurements to obtain information on the position, as well as the attitude and angular rate of the vehicle. The developed model for LTV, which is similar to the launcher model, can be used to validate the reusable launcher GNC in ascending/descending and horizontal flight phases.
In summary, the novelty element of the paper, and at the same time, the contributions in the field are:

Appendix A. Coupled form of the Linear Equations of Motion
The objective of the appendix is to obtain the coupled linear form of the equations of motion. The basic movement considered is a translation in a vertical plane. This appendix uses the same nomenclature as the main paper.

A.1. Translation Movement
Starting from the nonlinear equations for translation in local frame: The dynamic equations: The kinematic equations: The linear form of the translation equations becomes: The dynamic equations: (A3) The kinematic equations: For the kinematic equations, the previous relation becomes the final form.
Regarding the dynamic equations, aerodynamic terms can be written as: where K = S 2 ρC d depends on altitude y 0 and velocity V. The linear form of the aerodynamic term becomes: The form of the aerodynamic term is: where: and: In order to evaluate the influence of altitude and velocity, based on the hypothesis of the spherical shape of the vehicle, for the considered velocity and altitude regime, through regression, the following approximation functions were obtained: where: Next, we resume the development of the translational aerodynamic term (A7), which becomes: Because in the basic movement V z = 0, the translational aerodynamic term becomes: Returning to the dynamic equations of translation (A3), the linear form of the propulsive term is: For the first term of the propulsive force development, we start from the nonlinear expression: which can be expressed in the linear form: Because in the basic movement, the yaw command is null δ m = 0, we will obtain: For the second term of relation (A13), relying on the inverse rotation matrix: q 2 4 + q 2 1 − q 2 2 − q 2 3 2(q 1 q 2 − q 3 q 4 ) 2(q 3 q 1 + q 2 q 4 ) 2(q 1 q 2 + q 3 q 4 ) q 2 4 + q 2 2 − q 2 3 − q 2 1 2(q 2 q 3 − q 4 q 1 ) 2(q 3 q 1 − q 2 q 4 ) 2(q 2 q 3 + q 4 q 1 ) q 2 4 + q 2 3 − q 2 1 − q 2 2   , the derivatives with respect to the quaternion components are: In the basic movement, we have: q 1 = 0; q 2 = 0; δ m = 0, from where: where: P 31 = q 3 q 4 cosδ n − q 2 4 sinδ n , P 32 = −q 2 4 cosδ n − q 3 q 4 sinδ n , P 13 = −2q 3 q 4 cosδ n + 1 − 2q 2 3 sinδ n ,P 23 = 1 − 2q 2 3 cosδ n + 2q 3 q 4 sinδ n .
Cumulating the two terms, the linear form of the propulsive term from the translation Equation (A13) is obtained: (A22)

A.2. Rotational Movement
In this section, we will obtain the linear form of the rotational equations in the body frame: The nonlinear rotational equations are: The dynamic equations: The kinematic equations: If we consider as a basic movement the translation in the vertical plane, we have: q 1 = 0; q 2 = 0; q 4 = ± 1 − q 2 3 , with zero angular velocity: p = 0, q = 0, r = 0, the linear form of the rotational equations become: The dynamic equations: For the kinematic equations, the previous relation becomes the final form.
Regarding the dynamic equations, linear forms of the applied moments must be expressed.
For the propulsive moments of (A25), we resume the linear form of the translation terms (A16): and we consider that the command propulsive moments are: Regarding the roll moment command, it is considered as: where δ l is the equivalent roll deflection.