Huan Zhang

Zhihua Zhao
email: zhaozh@tsinghua.edu.cn

Gexue Ren

Pengxiang Hu
email: hupengxiang@tsinghua.org.cn.

Yunfei Yang

Zhongwen Pan
email: zhwpan@sina.com.

Francesco Sanfedino
email: francesco.sanfedino@isae-supaero.fr.

Ph D Candidate

Arresting-Cable System for Robust Terminal Landing of Reusable Rockets

Recent successful recovery techniques for rockets require that rockets maintain a vertical configuration with zero vertical and lateral velocities; otherwise, landings may fail. To relax this requirement, a new active-arresting system (inspired by the arresting gears used on aircraft carriers) is proposed herein to achieve a robust landing, even if the rocket deviates from the target position or has notable residual velocities and inclinations. The system consists of four deployable onboard hooks above the rocket's center of mass, an on-ground apparatus containing four arresting cables forming a square capture frame, and four buffer devices to actively catch and passively decelerate the landing rocket. To catch the rocket, the capture frame was controlled by servo motors via a simple proportional-derivative controller. After catching, the buffer devices generate decelerating forces to stop its motion. A flexible multibody model of the proposed system was built to evaluate its robust performance under various combinations of multiple uncertainties, such as noise measurement, time delay in the motor, initial conditions, and wind excitation. Using a quasi-Monte Carlo method, hundreds of deviated landing cases were generated and simulated. The results confirmed the robustness of the proposed system for achieving successful terminal landings.

D EVELOPING reusable rockets is supposed to be a promising technical path to reducing the cost of accessing and exploring space [START_REF] Reed | Performance Efficient Launch Vehicle Recovery and Reuse[END_REF][START_REF] Ragab | Launch Vehicle Recovery and Reuse[END_REF][START_REF] Donahue | Low Recurring Cost, Partially Reusable Heavy Lift Launch Vehicle[END_REF]. To achieve this goal, more than a half-century ago, scientists and engineers began exploring different decelerating and landing techniques to recover rockets or their parts [START_REF] Reed | Performance Efficient Launch Vehicle Recovery and Reuse[END_REF][START_REF] Ragab | Launch Vehicle Recovery and Reuse[END_REF][START_REF] Donahue | Low Recurring Cost, Partially Reusable Heavy Lift Launch Vehicle[END_REF][START_REF] Inatani | Concept and Preliminary Flight Testing of a Fully Reusable Rocket Vehicle[END_REF]. One of the most successful attempts is the space shuttle program, which recovered its orbiter and the two solid boosters. However, the space shuttle's hard-landing technique demands testing and maintaining the recovered parts, which increase rather than decrease launching costs [START_REF] Heinrich | Greenspace: Recovery & Reusability Scenarios for Launcher Industry[END_REF]. Nevertheless, an exciting vertical-recovery technique has recently been demonstrated by both SpaceX and Blue Origin: two space companies from the United States. By combining retropropulsion engines and landing legs to realize a soft landing, both companies have recovered and reused the first stages of their Falcon 9 [START_REF] Davis | First Stage Recovery[END_REF] and New Shepard [START_REF] Sippel | Systematic Assessment of Reusable First-Stage Return Options[END_REF] rockets, respectively. These successes not only convinced people that space can be accessed in a more economically feasible way than before, but they also encouraged scientists from other countries to develop their own techniques for reusable rockets. However, the precise and undamaged recovery of a rocket stage remains an extremely demanding task.

There are mainly two steps toward recovering a rocket stage. The first step is navigating the stage to a targeted landing place, and the next step is a soft terminal landing. In the first navigation step, significant achievements have been made by learning from various mature Mars landing technologies accumulated over the past 50 years [START_REF] Bonfiglio | Landing-Site Dispersion Analysis and Statistical Assessment for the Mars Phoenix Lander[END_REF][START_REF] Yu | Design and Optimization of Navigation and Guidance Techniques for Mars Pinpoint Landing: Review and Prospect[END_REF]. Now, based on well-developed high-precision guidance technologies, the guidance errors have been diminished from kilometers to meters with nearly minimum fuel consumption [START_REF] Wang | Optimal Rocket Landing Guidance Using Convex Optimization and Model Predictive Control[END_REF][START_REF] Simplício | Guidance of Reusable Launchers: Improving Descent and Landing Performance[END_REF][START_REF] Bojun | High-Precision Adaptive Predictive Entry Guidance for Vertical Rocket Landing[END_REF]. In the second terminallanding step, the current vertical-recovery technique requires not only advanced engines that are capable of multiple restarts and thrust modulation within a large range [START_REF] Pérez-Roca | A Survey of Automatic Control Methods for Liquid-Propellant Rocket Engines[END_REF][START_REF] Wang | Optimal Rocket Landing Guidance Using Convex Optimization and Model Predictive Control[END_REF][START_REF] Halbe | Robust Reentry Guidance of a Reusable Launch Vehicle Using Model Predictive Static Programming[END_REF] but also an accurate control system to prevent the rocket from crashing [START_REF] Blackmore | Autonomous Precision Landing of Space Rockets[END_REF][START_REF] Nonaka | Vertical Landing Aerodynamics of Reusable Rocket Vehicle[END_REF]. Moreover, the high heat flux on the base plate [START_REF] Ecker | Aerothermal Analysis of Reusable Launcher Systems During Retro-Propulsion Reentry and Landing[END_REF] and the reconfigured aerodynamic characteristics [START_REF] Nonaka | Vertical Landing Aerodynamics of Reusable Rocket Vehicle[END_REF] during the touchdown phase require thermal protection for the structure and accurate landing [START_REF] Yonezawa | Investigation of Effect of Ground Downstream of Overexpanded Dual-Bell Nozzle During Vertical Takeoff and Landing[END_REF].

In practice, managing a rocket to achieve a safe terminal landing with landing legs is a very tough task. First of all, landing a rocket is a one-shot process; and there are no second chances [START_REF] Blackmore | Autonomous Precision Landing of Space Rockets[END_REF]. Second, when a rocket touches down, it must be kept vertical with nearly zero vertical or lateral velocities. At the same time, the engines must be shut down. If the rocket is leaning, if its speed is too high, or if its engines are shut down too late, then the landing legs might fail with the potential risk of explosion. It is estimated that the vertical landing speed must be less than 2 m∕s and the inclination angle must be less than 5 deg to guarantee a successful landing [START_REF] Zhang | Adaptive Fault-Tolerant Control for a VTVL Reusable Launch Vehicle[END_REF][START_REF] Huang | Control Strategy of Launch Vehicle and Lander with Adaptive Landing Gear for Sloped Landing[END_REF]. According to a video issued by SpaceX, ‡ ‡ more than five Falcon 9 rockets have failed to land under these narrow constraints. These expensive failures have led us to think about developing an alternative landing approach that significantly expands the margin for a safe recovery, reduces the strain on the engines and control system, and achieves a higher rate of safe landings.

Currently, the field of aeronautics already offers a proven means of achieving a robust landing: the arresting gear used on aircraft carriers to recover airplanes within an incredibly short distance [START_REF] Yang | A Review of Civil Aircraft Arresting System for Runway Overruns[END_REF][START_REF] Shen | Modeling and Simulation of Arresting Gear System with Multibody Dynamic Approach[END_REF]. In a typical design, several arresting cables are strung across the deck of the carrier, and the pilot actively guides the plane to snag its tail hook on one of the cables. Then, the braking mechanism connected to the cable rapidly decelerates the airplane. This system is able to recover airplanes of various weights and velocities. Moreover, the landing planes are allowed to deviate from the orientation of the runway centerline; therefore, an airplane that is partially out of control may still be rescued by an arresting-cable system.

In this study, a novel terminal-landing technique is proposed to recover reusable rockets. The basic concept is to turn the horizontal arresting-cable systems on aircraft carriers into a vertical configuration. Since there is no pilot in the rocket, the positions of the cables are actively adjusted to catch the rocket. The proposed arresting system is mounted on towers, which can be located on ground or ocean platforms. This system offers the following three features: 1) capturing rockets under deviations in their lateral landing positions, inclination angles, lateral velocity, and angular velocity; 2) decelerating and stopping rockets with residual velocities within a prescribed maximum distance, maximum acceleration, and maximum cable stress; and 3) maintaining the recovered rocket in a vertical and stable orientation after it has come to a stop.

The rest of this paper is organized as follows: Sec. II introduces the overall design of the terminal-landing system; Sec. III establishes a flexible multibody dynamic model of the proposed system; Sec. IV provides simulation results and statistical results to verify the accuracy, feasibility, and applicability of the proposed system; and Sec. V concludes with the pros and cons of the proposed terminal-landing device.

II. Overall Design of the Terminal-Landing System

Figure 1a shows the overall structure of the proposed original arresting-cable system for the recovery of a reusable rocket stage. The system has four supporting truss towers arranged in a square, four guide rails connecting these towers, and four sliders moving along the rails. Along the rails are four identical arresting-cable subsystems. As schematically shown in Fig. 1b, each subsystem has a single cable connecting two opposing sliders after passing through several fixed passive pulleys, a fixed active pulley, and several moving pulleys. If the active fixed pulley is rotated, one slider can be moved along the guide rail and dragged by one cable of a subsystem. The moving pulleys are connected to a buffer device that consists of a counterweight and a friction brake. The recoverable rocket itself is equipped with four hooks that are symmetrically distributed above its center of mass.

It is assumed the rocket stage is navigated to a target position, within the supporting truss frame. With the process of the stage decent, the arresting cables are actively moved to catch the hooks. Once these hooks are caught by the cables, the downward motion of the rocket will drag the buffer devices, which in turn generates decelerating forces that slow the rocket. Because the hooks are located above the rocket's center of mass, the resulting pendulum motion ensures that the rocket is stable. In summary, the four arresting-cable subsystems, together with the hooks, are designed to track, capture, decelerate, and maintain the stability of the rocket during its landing phase.

In this design, the entire landing process can be separated into two sequential phases: 1) the active-tracking/capturing phase, and 2) the passive-decelerating phase. During the first phase, as shown in Fig. 2a, the counterweights are locked by friction, and the servo motors collectively reel in or reel out the cables to move the capture frame enclosed by the pre-tensioned arresting cables (as shown in purple); this is done to catch the rocket. Catching occurs when the rocket falls to a designated height and ends when the hooks engage with the arresting cables. In the second phase, as shown in Fig. 2b, all motors are powered down and the hooks drag the arresting cables downward, which in turn moves the counterweights upward and results in frictional forces from the brakes. During this process, the kinematic energy of the rocket is partially transferred into gravitational energy of the counterweights and partially dissipated by the friction of the brakes. Designing this landing device would require two central tasks: creating a control algorithm for the motors in the active-tracking/capturing phase, and designing the counterweights and brakes of the buffer devices for the passive-deceleration phase.

A. Designing the Active-Tracking/Capturing Phase

Figure 3a shows the active-control system, which contains specific hardware (i.e., laser-measuring instruments, filters, a calculation unit, a controller, and four servo motors) and a proportional-differential (PD) algorithm for controlling the motors to actively track and capture the landing rocket. This control system aims to ensure that 1) the center of the capture frame coincides with the target landing position of the rocket, 2) the bottom of the rocket goes through the hollow capture frame, and 3) the hooks on the rocket catch the arresting cables.

In designing the PD algorithm, the center position of the capture frame is needed. However, obtaining that information is not a straightforward task due to the fact that the frame is surrounded by flexible cables. Fortunately, as shown in Fig. 3a, if the cables are sufficiently pre-tensioned, then the position of the capture frame could be approximated by positions of the driven sliders. Therefore, the algorithm alternately controls the center positions of the sliders, which were selected as control variables in the form of d x s 1 ; x s 2 ; y s 3 ; y s 4 ⊤ . The effect of cable flexibility on the rocket capturing task is evaluated in Sec. IV. Next, the target landing position of the rocket is determined. During the rocket stages' decent, as shown in Fig. 3a, if the entire rocket was above the top plane oxy, the target position of the rocket r t was set to be a projection of the rocket's bottom on plane oxy. If the rocket crosses this plane, the target position was set to be the intersection point between the central axis of the rocket and plane oxy. Here, the position vector of the head and the bottom of the rocket are measured by a laser-measuring instrument in the arresting coordinate system oxyz. The measured position vectors of the head and bottom of the rocket are r m h and r m b , respectively. The real positions are r 0 h and r 0 b , which differ from the measured value by the measurement noises of δr h and δr b :

r m h r 0 h δr h ; r m b r 0 b δr b (1)
To eliminate measurement noises, filters based on a moving average method were adopted; and the measured positions were filtered as

r h t 0.3r m h t -Δt 0.4r m h t -2Δt 0.3r m h t -3Δt (2a) r b t 0.3r m b t -Δt 0.4r m b t -2Δt 0.3r m b t -3Δt (2b)
where r h x h ; y h ; z h ⊤ and r b x b ; y b ; z b ⊤ , and Δt is the sampling time period. The measured noises were assumed to be threedimensional isotropic and ergodic noises, and they have zero mean values and nonzero root mean squares (RMSs) σ r . These filtered positions, r h and r b , were used to calculate the target position of the rocket r t x t ; y t ; z t ⊤ , depending on whether the rocket was above or across the top plane oxy:

r t 2 6 4 x b y b 0 3 7 5; if z b ≥ 0; r t 2 6 6 6 6 4 x b - x b -x h z b -z h z b y b - y b -y h z b -z h z b 0 3 7 7 7 7 5 ; if z b < 0 (3)
Then, the control target variable was set as d target x t ; x t ; y t ; y t ⊤ , for which each component corresponds to the target displacement of the sliders from one to four.

As shown in Fig. 3b, e d targetd was used to denote the error between the target and the controlled positions. The controlled reelout velocities of active pulleys or servo motors v control P were designed by a conventional PD-control algorithm:

v control P a p e _ e ( 4
)
where a p is the proportional gain with a positive constant value. However, this simple PD-control law has a startup problem. At the initial tracking time of t 0 s, the capture frame was placed at the center of the truss frame, and the target positions were distant from the center. Therefore, the PD-control law required a large starting value in terms of the pulley reelout velocity. Meanwhile, the pulley is at rest before tracking so that there is a significant velocity discontinuity at the beginning. To avoid this problem, a saturation function of time κt was multiplied to the original PD-control law v control P to generate a smooth control command for the reelout velocities of the active pulley or servo motors, i.e.,

2t∕t c t < t c ∕2 1 t ≥ t c ∕2 (6)
and t c is the time duration of the active-tracking/capturing phase. Continuously inputting the control command v smooth P into the servo motors, the motors' time-delay feature makes the real output velocity v P t and displacement d P t differ from the inputted target, assuming they satisfy the following discrete time-delay relationship [START_REF] Jirsa | Will a Large Complex System with Time Delays be Stable?[END_REF]:

v P t -d P t Z t-τ 0 v smooth P t dt ( 7
)
where τ is the time-delay parameter. Additionally, the pre-tensioned cables were wound tightly on the pulleys and are reeled in or out by the motor; therefore, the reelout displacement of the active pulleys was the same as the displacement of the sliders, namely, d P d.

Each of the four servo motors has an angular-displacement sensor mounted on its shaft, enabling it to obtain the positions of the four sliders as well as the reelout displacement of the active pulleys. Thus, the active-control system can be regarded as a closed control loop, as shown in Fig. 3b. The transfer function Gs from the error vector e to the real reelout velocities of the active pulley v P can be obtained by doing the Laplace transform to Eqs. (4-7), which can be expressed as follows:

Gs ≜ V P s Es a p se -τs Ks s 1 ( 8
)
where s is the Laplace variable; and V P s, Es, and Ks are the Laplace transforms of v p t, et, and κt. Based on this equation, the effectiveness of the aforementioned designed control algorithm on tracking the landing rocket was evaluated. The evaluation was separated into two steps: τ 0 and τ ≠ 0. First, assuming τ 0, then the theoretical solution of the activecontrol system for t < t c ∕2 is as follows:

d d target 1 -e -a p ∕t c t 2 (9)
Hence, the relative error is e d target e -a p ∕t c t 2 [START_REF] Wang | Optimal Rocket Landing Guidance Using Convex Optimization and Model Predictive Control[END_REF] This response asymptotically converges to zero as long as a p is positive. The larger the proportional gain of a p , the quicker the converged result will be obtained. To obtain a fast convergence, such as letting the relative error be less than e -5 at t t c ∕2, Eq. ( 10) says that the gain a p can be selected as

a p ≥ 20∕t c (11)
Second, a case of τ 0.1 s ≠ 0 was evaluated. The other parameters of the case are the RMS of the measurement noises of σ r 0.1 m, the duration of the active-tracking/capturing phase of t c 10 s, a p 2, and e 17. Figure 4 compares the performance of the control system with and without the saturation function, using slider 1 and pulley 1 as examples. It could be conjectured that the saturation function reduces the maximum reelout velocity of the motor at the initial time.

B. Designing the Passive-Decelerating Phase

Just before one of hooks touches an arresting cable, the engines on the rocket need to be shut down and the servo motors powered off. Then, the passive-decelerating phase begins. As illustrated in Fig. 5, the falling rocket drags all four buffer devices, and its mechanical energy is transferred to the counterweights and dissipated as heat from braking friction. Therefore, the central task of designing the decelerating phase is to select the mass of the counterweight m w and the friction force F f based on the following analysis of the dynamic process.

For a rough design estimate, a reference condition was adopted to analyze the dynamic process of the passive-deceleration phase. It was assumed that the deformation of the cables was negligible, the rocket was landing vertically at the center of top plane oxy, and all of the buffer devices worked identically. In this case, the entire system had only one degree of freedom, and can therefore be characterized by a generalized coordinate, the angle θ between the arresting cable, and the horizontal line, as shown in Fig. 5. Kinematically, if the rocket moves downward from top plane oxy by distance z l 0 tan θ, where l 0 is half the length of the supporting truss. The counterweights are correspondingly raised up by the following equation:

w αl 0 sec θ -1 α z 2 l 2 0 q -l 0 ( 12
)
where α is the constant transfer rate of the moving pulley group. Energywise, the variation of the system's mechanical energy equals the work done by the frictional forces. This yields the equation of motion for the decelerating process, which is expressed as

1 2 m r v 2 0 m r gz 1 2 m r _ z 2 1 2 4m w _ w 2 4m w g F f w (13)
Here, m r is the mass of the rocket, v 0 is the rocket's initial vertical velocity when z 0, g is the standard gravitational acceleration, and a dot over a symbol represents a differentiation with respect to time t.

a
To make the analysis clearer and more concise, all quantities were nondimensionalized and denoted by a bar over the symbols. Here, the reference values of length were chosen as l 0 , of time as l 0 ∕g p , and of mass as m r . Thus, quantities of velocity were normalized by gl 0 p , of acceleration by g, and of force by m r g. Then, the equation of motion [Eq. [START_REF] Pérez-Roca | A Survey of Automatic Control Methods for Liquid-Propellant Rocket Engines[END_REF]] could be reformulated into another dimensionless form as

v 2 0 2 z z 02 2 m t α 2 w 02 2 F t α w ( 14
)
where the prime symbol indicates a differentiation with respect to dimensionless time t ≜ t∕ l 0 ∕g p , and

m t ≜ 4α 2 m w ; F t ≜ 4α m w F f ( 15
)
were introduced to simplify the expression. In addition,

w αsec θ -1; z tan θ (16)
To obtain an explicit expression for the velocity of the rocket, the kinematic expressions of Eq. ( 16) were differentiated with respect to the dimensionless time t. Also, u 0 , z 0 , and θ 0 are related as

w 0 α sin θ z 0 ; θ 0 z 0 cos 2 θ ( 17
)
Substituting Eqs. ( 16) and ( 17) into Eq. ( 14) directly links the velocity of the rocket with the generalized coordinate θ:

z 02 2 1 m t sin 2 θ " v 2 0 2 tan θ - F t sec θ -1 # (18)
Equation ( 18) implies that the dimensionless velocity is completely determined by two design parameters: namely, m t and F t . Their selection fulfills the following three requirements:

1) The rocket is stopped within an allowed dimensionless distance z m .

2) The maximum acceleration of the landing rocket is less than an accepted dimensionless value a m .

3) The maximum tension in the cables is smaller than a dimensionless tension

T m to leave a sufficient margin for the cable strength. The proposed arresting system was designed for a rocket mass of mass m r 25 ton and a residual vertical velocity of v 0 20 m∕s. Geometrically, the chosen half-length of the truss frame was l 0 25 m, and the rocket was fully stopped when its displacement reached the maximum of z m 10 m. For this system, the corresponding dimensionless parameters are

v 0 v 0 ∕ gl 0 p 20∕ 9.81 × 25 p
1.28 and θ m arctan10∕25 0.38 rad. Substituting the preceding parameters and z 0 0 into Eq. ( 18) yields a requirement on

F t : F t v 2 0 ∕2 tan θ m sec θ m -1 1.28 2 ∕2 tan 0.38 sec 0.38-1 15.86 (19)
To obtain the acceleration value of the rocket, differentiating Eq. ( 18) with respect to the dimensionless time t and making use of the relationships of Eqs. ( 17) and ( 19) yield the following:

z 0 0 1 - F t sin θ 1 m t sin 2 θ - m t sin2θcos 2 θ 1 m t sin 2 θ 2 v 2 0 2 tan θ - F t sec θ -1 (20)
Since both v 0 and F t are given, the maximum acceleration of the rocket only depends on the parameter m t . Figure 6 gives the j z 0 0 j maxm t curve as a red solid line; it is not a monotonic function. Setting a m 4 and then demanding j z 0 0 j max ≤ a m require m t ≥ 1.55 [START_REF] Huang | Control Strategy of Launch Vehicle and Lander with Adaptive Landing Gear for Sloped Landing[END_REF] To obtain the tension in the arresting cables, Newton's second law was applied to the rocket's vertical direction, yielding

m r z m r g -8T sin θ (22)
This equation can be nondimensionalized as

T 1 8 sin θ 1 - z 0 0 (23)
As with the dimensionless acceleration of the rocket, the maximum cable tension only depends on the parameter m t . Figure 6 shows the T maxm t curve as a blue dashed line.

Fig. 6 The relationship between j z 0 0 j max , T max , and m t . Here, the arresting cable is the same type as that is used in aircraft carriers. Its radius was 20 mm, and the allowed stress was 800 MPa such that T m 800 × π × 20 2 ∕25 × 10 3 × 9.81 4.1. Therefore, requiring

T max ≤ T m requires m t ≤ 10.2 (24)
In fulfilling the requirements of Eqs. ( 21) and ( 24), m t was set to two since this design is concerned more with cable tension than the acceleration of the rocket. Then, the mass of the counterweight m w and the frictional force F f were determined by the transfer rate α according to Eq. [START_REF] Halbe | Robust Reentry Guidance of a Reusable Launch Vehicle Using Model Predictive Static Programming[END_REF].

The design and structural parameters of the system are listed in Table 1. The center of mass of the landing rocket was generally below the geometric center, and a 12.5 ton landing rocket mounted by an engine of 12.5 ton at its bottom was selected as the object of this study. For simplicity, the mass of the rocket's thin-walled structure was uniformly distributed along the axial direction.

III. Dynamical Simulation Model

Although the design process described in the previous section assumed the rocket to be descending vertically with zero lateral velocity and along the center of the supporting truss, these assumptions are not needed for the proposed system. In particular, when the rocket deviates from these working conditions (i.e., it is off center from the truss frame, it is tilted, or it is moving sideways), then its dynamics are much too complex to be computed analytically. Therefore, a dynamical simulation model was built to evaluate the landing performance of the proposed arresting-cable system under different landing conditions.

A. Governing Equations of the System

The model was based on a multibody approach, as shown in Fig. 7. The truss structure, guide rails, pulleys, and sliders are all simplified as massless rigid bodies because they are not important to the dynamic process. For each buffer device, the counterweight was considered as a point mass, and the brake was replaced by a stick-slip friction force.

To consider the elastic vibration and small deformation of the rocket, the rocket was modeled by lumped rigid bodies connected by masslessbeam elements [START_REF] Hu | Multibody Dynamics of Flexible Liquid Rockets with Depleting Propellant[END_REF]. The rigid hooks were hinged to the rocket and were also connected to the rocket by springs and dampers. The arresting cables were meshed with variable-length cable elements based on an arbitrary Lagrangian-Euler (ALE) description, and the contacts between the hooks and cables were accounted for with Hertz contact forces [START_REF] Liu | A Model Order Reduction Method for the Simulation of Gear Contacts Based on Arbitrary Lagrangian Eulerian Formulation[END_REF][START_REF] Liu | An Arbitrary Lagrangian-Eulerian Formulation of a Geometrically Exact Timoshenko Beam Running Through a Tube[END_REF]. The control of the active-tracking/capturing phase was integrated into the dynamic equation of the system in the form of state space. The details are presented in the following subsections. In summary, the governing equations of the arresting-cable system and the landing rocket form a set of differential algebraic equations (DAEs) [START_REF] Hairer | Solving Ordinary Differential Equations. II: Stiff and Differential-Algebraic Problems[END_REF]:

8 > > > > > > < > > > > > > : Mq q Φ ⊤ q λ -Qq; _ q; t 0 Φq; _ q; v P ; t 0 e eq; t _ x Ax Be v P Cx De (25)
Here, q is the generalized coordinate vector of the system; Mq is the mass matrix of the system; Qq; _ q; t is the generalized force vector of the forces applied to the system; Φq; _ q; v P ; t is the vector of the kinematic and geometric constraints; Φ q is the gradient of Φ with respect to q; λ is the corresponding Lagrange multiplier; fA; B; C; Dg are the state-space matrices of the controller; x is the state variables; and e and v P are the control input and output of the controller, respectively. The whole system model was built using our laboratory's in-house multibody code, and the obtained DAEs were numerically integrated through an index-3 backward differentiation formula.

B. Modeling the Arresting Cables

The main difficulty in modeling the proposed recovery system is overcoming the problem of efficiency in modeling flexible arresting cables with a large displacement, as well as the moving points of contacts between the hooks and cables. The traditional finite element modeling method, based on the Lagrange formulation, meshes the cable as a fixed grid of constant-length elements and calculates the contact force between these cable elements and the hooks.

For contact detections, this method demands that the size of the cable elements that are potentially contacting with the hooks is much smaller than the hooks' radius. Therefore, very fine elements are necessary to mesh the cables if a large portion of them might make contact. This fact significantly increases the calculation scale and time.

To substantially improve the calculation efficiency, a new modeling method was established to handle the cables and their contact with the hooks. As shown in Fig. 7a, the basic idea was to dynamically adjust the cable meshes to guarantee that only a few elements contact the hooks and that these cable elements be finely meshed while the others are roughly meshed. Thus, the number of elements is notably reduced and the contact detections are focused on a limited range, thus accelerating the calculation. Technically speaking, this meshing was implemented by combining a kind of variable-length cable element and a length-adjusting algorithm.

Variable-Length Cable Elements

As shown in Fig. 8a, to vary the length of an element, the generalized coordinates of a cable node q 1 contain not only the global nodal position r 1 but also a material coordinate s 1 or

q 1 r ⊤ 1 ; s 1 ⊤ (26)
Taking the simplest two-node straight element as an example, its generalized coordinates are

q e r ⊤ 1 ; s 1 ; r ⊤ 2 s 2 ⊤ (27)
where subscripts 1 and 2 represent the node numbers. The element length L e s 2s 1 is variable, which is in sharp contrast to the traditional fixed length cable element. Since the material coordinate of the node is variable, this cable element is actually based on an ALE description [START_REF] Hong | Dynamic Modeling of Mass-Flowing Linear Medium with Large Amplitude Displacement and Rotation[END_REF][START_REF] Chen | Equivalence of Lagrange's Equations for Non-Material Volume and the Principle of Virtual Work in ALE Formulation[END_REF] and is called an ALE element; and the corresponding nodes are called ALE nodes.

The internal point position of the element was linearly interpolated as follows: rs; t N r sq r t;

q r r ⊤ 1 ; r ⊤ 2 ⊤ ; N r 1 -ξ∕2I 3 ; 1 ξ∕2I 3 ( 28
)
where N r is the shape function matrix, I 3 is the identity matrix of order 3, and ξ 2ss 1s 2 ∕s 2s 1 is a dimensionless material coordinate.

For an infinitely small cable segment with length ds around the material coordinate s ∈ s 1 ; s 2 , the elastic force F e , external force f, and inertial force -ρA r are balanced, where ρA is the mass density of the cable per arc length. Therefore, the following integration of the virtual work is zero:

δW Z s 2 s 1 δr ⊤ f F e -ρA r ds 0 ( 29
)
Taking the second time derivatives of Eq. ( 28) and substituting the obtained result into Eq. ( 29) yield the governing equation of the cable element:

M e q e Q s Q E Q f ( 30
)
where

M e s 2 -s 1 2 Z 1 -1 ρAN ⊤ Ndξ Q s - s 2 -s 1 2 Z 1 -1 ρAN ⊤ r s dξ; Q E - s 2 -s 1 2 Z 1 -1 ∂ϵ ∂q e ⊤ EAkϵ β_ ϵ dξ; Q f s 2 -s 1 2 Z 1 -1 N ⊤ f dξ; N N r ∂N r ∂s 1 q r ∂N r ∂s 2 q r ; ϵ 1 2 ∂r ⊤ ∂s ∂r ∂s -1 ; k 1 ϵ ≥ 0; 0 ϵ < 0 ; r s 2 ∂N r ∂s 1 _ s 1 ∂N r ∂s 2 _ s 2 _ q r ∂ 2 N r ∂s 2 1 _ s 2 1 2 ∂ 2 N r ∂s 1 ∂s 2 _ s 1 _ s 2 ∂ 2 N r ∂s 2 2 _ s 2 2 ! q r (31)
where M e is the generalized mass matrix of the element, Q s is the inertial force vector caused by the mass flow, Q f is the generalized external force vector, Q E includes the generalized elastic and damping forces, ϵ is the axial Green strain, k represents the cable can resist only tension and not compression, and β is the stiffness proportional to the Rayleigh damping factor. The cables in the four subsystems are adopted using the same modeling method and the same number of elements. Each arresting cable is modeled by n 1 ALE elements, and the transfer cable is modeled by n 2 ALE elements. There are 4n 1 4n 2 ALE elements in the proposed system.

Mesh-Adjusting Algorithm

As shown in Fig. 7, all the cables in the terminal-landing system were meshed with the aforementioned variable-length cable elements, and the mesh was dynamically adjusted to ensure that only a few predefined elements have the potential to be in contact with the rocket hooks.

As for the mesh, let us take a single arresting cable as an example (indicated in red in Figs. 7 and8b). This cable was separated into several segments by key points along its length, including its two ends, the points coinciding with the pulleys (ignoring the radius), and the intersection points with the adjacent orthogonal arresting cables. For each key point, a cable node was assigned; these are called boundary nodes, and they are marked as solid black dots. Any cable segment surrounded by two boundary nodes was divided into a b) Mesh-adjusting constraint a) A typical two-node ALE cable element number of elements by uniformly distributing internal nodes, marked as hollow black dots.

The positions of the mesh nodes were adjusted and determined by introducing the following three kinds of constraint equations into the model. First, if a boundary node is at the end of a cable, located on a slider, then its position is the same as that of the slider and its material coordinate is a constant since there is no material flow. The equation is expressed as follows:

Φ 1 " r node -r slider s node -s 0 node # 0 ( 32
)
where s 0 node is a constant for the material coordinate of the node at the initial time.

Second, if a boundary node coincides with a pulley or an intersection point, then its position r node equals the corresponding pulley's or point's position r p ; and the constraints relate only to position [START_REF] Hong | A Modeling of Sliding Joint on One-Dimensional Flexible Medium[END_REF][START_REF] Fu | From Multiscale Modeling to Design of Synchronization Mechanisms in Mesh Antennas[END_REF] or are expressed as

Φ 2 r node -r p 0 (33)
Third, consider a cable segment bounded by two boundary nodes, b 1 and b 2 , whose material coordinates are s b 1 and s b 2 . The segment was meshed into N elements by N -1 internal nodes, and the material coordinates of these internal nodes were distributed to equally split the length of the segment. Then, the material coordinate s i of the ith node, numbered from the boundary node b 1 , was determined to be

Φ 3i s i -s b 1 s b 2 -s b 1 i∕N; 0 (34)
where i 1; 2; : : : ; N -1. The number N is set to be large for the cable segments that contact the hooks and small for the other segments [START_REF] Zhang | An Efficient Multibody Dynamic Model of Arresting Cable Systems Based on ALE Formulation[END_REF].

C. Flexible Model of a Landing Rocket

Since the rocket's propellant was nearly exhausted in the terminallanding phase, its structures were modeled as a flexible multibody model in which a cluster of lumped rigid bodies was connected by massless-beam elements [START_REF] Hu | Multibody Dynamics of Flexible Liquid Rockets with Depleting Propellant[END_REF] and the hooks and engines (as well as the counterweights and sliders) were modeled by the rigid bodies and constraints.

Rigid Bodies

As shown in Fig. 9, the generalized coordinates of a rigid body i in quaternion representation are

q i r ⊤ i ; θ ⊤ i ⊤ ( 35
)
where r i are the position coordinates of the center of mass of the rigid body i; and θ i θ 0 i ; θ 1 i ; θ 2 i ; θ 3 i ⊤ are the orientation quaternion parameters of the ith rigid body, which satisfies θ ⊤ i θ i 1. Taking the first time derivatives of Eq. ( 35), calculating the kinetic energy of the ith rigid body T i 1∕2 _ q ⊤ i M i _ q i , and substituting the obtained result into the second kind of Lagrangian equation in analytical mechanics yield the following governing equation of the ith rigid body:

M i q i Q q Q f (36)
Here, Q q -_ M i _ q i ∂T i ∕∂q i ⊤ is the generalized inertial force vector, and Q f H i f is the generalized external force vector. In addition,

M i " m i I 3×3 0 0 4G ⊤ i J i G i # ; H i " I 3×3 0 3×4 0 4×3 2G T i # (37)
G i 2 6 6 4 -θ 1 i θ 0 i θ 3 i -θ 2 i -θ 2 i -θ 3 i θ 0 i θ 1 i -θ 3 i θ 2 i -θ 1 i θ 0 i 3 7 7 5 (38)
where J i is the inertial tensor of the ith rigid body. The external force vector f was determined by a combination of gravity, the massless-beam force, aerodynamics, the spring-damper force, the contact force, and the friction of different rigid bodies in this system as follows:

f 8 > > < > > : f grav f beam f aero f sprdmp rocket f grav f contact f aero f sprdmp hook f grav f friction counterweight (39)
This will be discussed in detail in the next subsection. The rocket was clustered by n b rigid bodies and connected by n b -1 masslessbeam elements. The four hooks, the four sliders, and the four counterweights were all modeled as rigid bodies and were constrained to only one degree of freedom, as described in the following. Thus, there are n b 12 rigid bodies in this system.

Constraints of Rigid Bodies

In the arresting-cable system, the hooks on the rocket were designed to only rotate relative the hinged pins, whereas the counterweights and sliders can only slide along the slideway and guide rails. These constraints were modeled by a set of algebraic constraint equations.

For a revolute constraint between the rocket and a hook, revolution was only allowed around the z axis of the two objects, as shown in the upper dotted box of Fig. 7a; and the constraint equations are given as follows:

Φ 4 2 6 6 4 r 1 -r 2 z ⊤ 1 x 2 z ⊤ 1 y 2 3 7 7 5 0 ( 40
)
where r 1 x 1 ; y 1 ; z 1 and r 2 x 2 ; y 2 ; z 2 are center positions of the relatively two revolute bodies, and x 1 , y 1 , z 1 and x 2 , y 2 , z 2 are orientation axes of the two bodies.

For a translational constraint between a counterweight and a slideway, movement was only allowed along the z axis of two objects, as shown in the bottom dotted box of Fig. 7a, and the constraint equations are given as follows:

x 1 -x 2 y 1 -y 2 z ⊤ 1 x 2 z ⊤ 1 y 2 x ⊤ 1 y 2 3 7 7 7 7 7 7 7 7 5 0 (41)
where subscripts 1 and 2 represent the two sliding bodies.

D. Loads

The forces exerted on the system will be presented as generalized forces in the system's governing equations in Eqs. [START_REF] Chen | Equivalence of Lagrange's Equations for Non-Material Volume and the Principle of Virtual Work in ALE Formulation[END_REF] and [START_REF] Dick | Higher Order Quasi-Monte Carlo Integration for Bayesian PDE Inversion[END_REF]. There are six types of forces: 1) the gravitational force modeled as a constant force; 2) the massless-beam force between two rigid bodies in the rocket modeled in the form of a force matrix [START_REF] Hu | Multibody Dynamics of Flexible Liquid Rockets with Depleting Propellant[END_REF]; 3) the spring damper between the hooks and the rocket modeled using a linear function of relative displacement and velocity between the hooks and the rocket; 4) the stick-slip friction on the counterweight modeled as a unified form of a sticking state and a sliding state based on the Coulomb friction law [START_REF] Gu | A Fully Coupled Elastohydrodynamic Model for Static Performance Analysis of Gas Foil Bearings[END_REF]; 5) the contacts on cables and hooks; and 6) the aerodynamic force on the cables and the rocket along with the cable or beam elements. More details on the latter two types are described in the following two sections:

Modeling Contacts on Cables and Hooks

The geometry of the rocket and hooks is shown in Fig. 1c. Using Hertz contact theory, the contact force between the hook and the arresting-cable element is

f contact f n n f τ τ ( 42
)
where n is the unit normal vector of the contact surface, τ is the corresponding unit tangential vector, f n is the normal collision force, and f τ is the tangential friction force. Collision forces are exerted on the nodes of the cable element. More details regarding these forces can be found in the Appendix.

Modeling Aerodynamic Forces

Both the descending rocket and the cables of the arresting-cable system are affected by wind excitation, which consists of steady wind and gust. Then, steady wind speed can be represented by the mean wind speed v ws . The gust speed v wg can be described by a normal stationary random process with zero mean and the nonzero RMS σ wg , which is obtained by the mean wind velocity v ws multiplied by a fluctuation intensity λ wg (i.e., σ wg λ wg v ws ). The wind direction is represented by the clockwise angle θ w with respect to the north. In summary, the wind field is determined by three parameters: the mean wind speed v ws , the fluctuation intensity λ wg , and the wind angle θ w .

As shown in Fig. 10, the velocity of an arbitrary point M of the cable or rocket relative to the wind, in the local coordinate system o m x m y m z m , can be expressed as

v A m a v M v w ( 43
)
where v M is the velocity of an arbitrary point M in the arresting coordinate system oxyz, A m a is the transformation matrix from the arresting coordinate system oxyz to the local coordinate system o m x m y m z m , and v w is the wind speed consisting of steady wind speed v ws and gust speed v wg (i.e., v w v ws v wg ).

The aerodynamic force on an arbitrary point M of the rocket-beam element or cable element in the arresting coordinate system oxyz can be obtained as follows:

f aero x; t A a m 2 6 6 4 -C D x; tS ref q a x; t C N x; tS ref q a x; tα a x; t -C N x; tS ref q a x; tβ a x; t 3 7 7 5 ( 44
)
where C D x; t and C N x; t are the drag and lift coefficients obtained from the wind-tunnel test, and S ref is the reference area. The angle of attack α a x; t, sideslip angle β a x; t, and dynamic pressure q a x; t of point M in the local coordinate system o m x m y m z m can be expressed as

α a arctan v y v x ; β a arcsin v z k vk ; q a 1 2 ρ a k vk 2 (45)
where ρ a is the atmospheric density at the current flight altitude. The drag forces are applied on the cables in the same way they are applied on the rocket.

E. Time-Domain State-Space Model of the Active-Tracking/Capturing

Algorithm

In the multibody dynamic model of the proposed system, the real position vector of the head and the bottom of the rocket, r 0 h and r 0 b , can be obtained by the position coordinates of the rigid bodies at the two ends of the flexible rocket model. The measured errors are modeled by the white noise with zero mean values and nonzero RMS σ r . After filtering by Eq. ( 2) with the sampling time period and conversion by Eq. (3), the target vector of the active-control can be obtained as a function of the generalized coordinates of the system q, i.e., d target d target q; t. Similarly, the control vector can be obtained by the position coordinates of the rigid bodies corresponding to the sliders, i.e., d dq; t. Therefore, the error vector was e eq; t.

The real control variables were the reelout velocities of the active pulleys as the actuator, i.e., v P . The active-control system was introduced into the multibody model in the form of state space by converting the transfer function [Eq. ( 8)] to the state space:

_ x Ax Be v P Cx De ( 46
)
where x is the state variables, and fA; B; C; Dg are state-space matrices.

When an active pulley reels in or out, it changes the length of the wound cable. This effect was modeled by constraining the material coordinate velocity of the corresponding ALE node:

Φ 6i _ s i -v p i 0 ( 47
)
where i 1; 2; 3; 4 represents the ith active pulley; and v p i and _ s i denote the reelout velocity and the material coordinate velocity of the active pulley i and the corresponding ALE node i, respectively. With this constraint equation and the preceding state-space equations, the major effect of the active-control system is formulated.

F. Validation of the Dynamic Modeling Method

To validate the aforementioned proposed modeling method, a numerical model was constructed with the parameters listed in Table 1, and an ideal working case in which the rocket lands vertically Fig. 10 The angle of attack α a and sideslip angle β a at point P.

along the center of the truss frame was simulated. In this case, the analytical solution, detailed in Sec. II.B, has already been obtained in Eqs. ( 20) and ( 23), assuming that the cables did not stretch. However, the cables have a certain amount of flexibility, which should affect the dynamics of deceleration.

Figure 11a shows several time snapshots of the landing process; and Figs. 11b-11e compare the analytical results (drawn in black solid lines) with the results of the numerical model (drawn in black dashed lines). Overall, the two results agree closely, suggesting the correctness of the numerical model. Also, it is not surprising to find that the simulation results oscillate around the analytical results due to the flexibility of the cables.

IV. Robust Landing Performance of the Arresting-Cable System

The proposed arresting-cable system was designed to both catch and decelerate rockets, even if they deviate from the ideal landing condition. In other words, a rocket is allowed to have residual linear and angular velocities, a landing position that deviates from the center of the truss frame, and an attitude that may depart from the vertical.

To quantify the deviations of a descending rocket, a local coordinate system OXYZ was attached to it. As shown in Fig. 1c, the XY plane passes through the four top points of the deployed hooks and Z is the rotation axis. Referring to the truss-frame coordinate system oxyz, the x and y positions as well as the three velocities of point O, together with the attitude and angular velocities of OXYZ, were adopted to describe the rocket deviations. For attitude, rocket coordinate system OXYZ was obtained by rotating the reference coordinate system oxyz with respect to a unit vector n n 1 ; n 2 ; n 3 ⊤ by an angle θ ∈ -π; π. Defining the rotation vector θ ≜ θn θ x ; θ y ; θ z ⊤ , its first two components (θ x and θ y ) represent the inclination angles of the rocket, for which θ z represents the rolling angle. The angular velocity is ω ω x ; ω y ; ω z ⊤ . Since the vertical landing velocity margin of the rocket was given as 20 m∕s, there are 10 uncertain parameters to describe the deviations of the rocket states, as listed in the first five rows of Table 2.

Except for the deviation of the rocket states, noises from the measurement, and the time delay of the motor system, the wind excitations will also affect the terminal-landing performance of the proposed arresting-cable system. To quantify them, five additional uncertain parameters were introduced, as stated in the previous sections: measurement noise by the standard deviation σ r , the time delay by τ in Sec. II.A, the steady wind by the stationary velocity v ws , the wind angle θ w , and the lateral gust fluctuating intensity λ wg . Together with the 10 state deviation variables for the rocket, there were a total of 15 uncertain parameters in the simulated dynamic system, which are indicated by the dark red boxes in Fig. 12. The robust evaluations for the proposed arresting-cable system will be sequentially conducted for the two working phases. For the active-tracking/capturing phase, the initial state deviation ranges of the rocket for simulation are listed in Table 3, and simulations were used to validate the system's capacity to capture the landing rocket and expand the range of landing state deviation. For the passivedecelerating phase, the initial state deviation ranges of the rocket were the captured state ranges obtained from the first active-tracking phase. The process of the entire two-phase simulation is shown in Fig. 12; and the transition between the active-tracking/capturing phase and the passive-decelerating phase was triggered by the contact force and the heights of the hooks.

To comprehensively evaluate the performance of the proposed system, three scenarios corresponding to three maximal steady wind speeds were studied, as listed in Table 4. For each scenario, a random quasi Monte Carlo method [START_REF] Hellekalek | Open Type Quasi-Monte Carlo Integration Based on Halton Sequences in Weighted Sobolev Spaces[END_REF][START_REF] Dick | Higher Order Quasi-Monte Carlo Integration for Bayesian PDE Inversion[END_REF][START_REF] Camberos | On Direct Simulation Quasi-Monte Carlo Methods[END_REF] was adopted to generate 100 cases to assess whether the designed system can capture and decelerate the rockets with the uncertain parameters listed in Table 2. To achieve the fast convergence of multidimensional random variables, the Halton sequence [START_REF] Chi | On the Optimal Halton Sequence[END_REF][START_REF] Hinrichs | Tractability Properties of the Weighted Star Discrepancy of the Halton Sequence[END_REF], provided by the haltonset function in MATLAB 2017a, was used to generate low-discrepancy and quasirandom sequences in the multidimensional virtual space by deterministic combinations of uncertain parameters.

A. Dynamics of a Deviated Landing Case

To provide an overview of the system's dynamic behaviors, a random landing case with all 15 nonzero deviations was simulated and analyzed in the sequential active-tracking/capturing and passivedecelerating phases. The detailed parameters of this case are x 13 m, y 13 m, v x v y 0.5 m, θ x θ y 5 deg, θ z 25 deg, ω x ω y ω z 0.1 deg∕s, τ 0.1 s, σ r 0.1 m, θ w 90 deg, and λ wg 10% in scenario I (i.e., v ws 5 m∕s, h 200 m, t c 10 s, and a P 2).

As shown in Fig. 13a, assume a rocket descends freely from a point over the truss frame at t 1 . At the same time, the capture frame, driven by the active-control system, begins to track the target point of the rocket. Then, at t 2 , the center of the capture frame, indicated by the positions of the sliders, coincides with the target point of the rocket as shown in Figs. 13b and14a, which satisfies the design of Eq. [START_REF] Simplício | Guidance of Reusable Launchers: Improving Descent and Landing Performance[END_REF]. Until t 3 , the hooks engage the arresting cables with the nonzero contact force as shown in Fig. 14c, and the passive-decelerating phase takes over the active-tracking/capturing phase. Then, the rocket moves downward, and its kinetic energy is partially transferred into the gravitational energy of the counterweights and partially dissipated by the friction of the brakes. At t 4 , the vertical velocities of the rocket and counterweights are zero. After this moment, the counterweights are stationary and the rocket moves toward the center until it is stationary and stable at t 5 .

In addition, Figs. 13b and14a show that, in the passive-decelerating phase, if a rocket was captured with a deviated state, the asymmetrical forces generated from the four arresting cables could gradually correct the position of the rocket to the center of arresting plane oxy. Similarly, Fig. 14b shows the roll angle also approaches a stable value. These results imply that the system has an automatic-correction capability similar to the arresting gears for aircraft carriers [START_REF] Zhang | An Efficient Multibody Dynamic Model of Arresting Cable Systems Based on ALE Formulation[END_REF]. More details regarding the dynamic process of this deviated case can be found in Supplemental Video S1.

A successful terminal landing requires the success of both the first active-tracking/capturing phase and the second passive-decelerating phase. There are two criteria to judge the success of the first active phase. The first is that the locus error kek between the center of the capture frame and the target point of the rocket (as shown in Fig. 13b) approaches zero, indicating a good tracking performance. The second is that the contact forces between the four hooks and the arresting cables are nonzero (as shown in Fig. 14c), indicating the frame captured the rocket.

For the second passive phase, there are three criteria to judge its success. First, the maximum arresting distance is less than 10 m; namely, jz max j < 10 m, as shown in Fig. 14d. Second, the maximum 3 Captured states in Table 8 Time delay

τ max 0.1 s -- -- -- Measurement noise σ max r 0.1 m -- -- --
arresting acceleration is less than 5 g; namely, ja max j < 5 g where g 9.80 m∕s 2 , as shown in Fig. 14e. Third, the maximum arrestingcable stress is less than 800 MPa; namely, jσ max j < 800 MPa, as shown in Fig. 14f. These criteria were used in the subsequent performance evaluations to judge whether the terminal landings were successful. Although the arresting-cable system is expected to be able to retrieve rockets that deviate from the ideal landing condition, this capability has its limits. To evaluate this capability, 300 simulations (100 cases for each scenario in Table 4) under the 15 uncertain parameters in Table 2 were carried out. The ranges of the initial rocket states in the simulations were the ranges of the terminalguidance states in Table 3. The process from the rocket reaching the terminal-guidance point to the rocket being stopped and stabilized was simulated.

As shown in Figs. 15a-17a, the black points represent the position deviation; and the red, green, and blue arrows represent the velocity, angle, and angular velocity deviation vectors, respectively. The upper and lower yellow areas indicate the terminal-guidance point range and the available capture area in plane oxyz, respectively. The lines between the two areas indicate the trajectories of the rockets of 100 shooting cases for each scenario.

The robustness of the active-tracking/capturing phase was evaluated by the capability of expanding the states of the rocket from the terminal-guidance states to the captured state. This capability can be described by the extension from the upper area to the lower area in Figs. 15a-17a and the extension of the deviation range from Table 3 to Table 5. In addition, the lateral velocity deviation and angular velocity deviation of the rocket increased with the increasing wind speed. The tracking errors kek (indicated in Fig. 13b capture frame coincides with the target point of the rocket at design moment t c ∕2 (i.e., 5, 3, and 2.5 s for scenarios I, II, and III, respectively). The tracking errors are small enough that all the hooks engaged with the arresting cables. According to the success criteria of the active-tracking/capturing phase, all 300 cases successfully captured the rockets.

After the rockets were captured, the landing process entered into the passive-decelerating phase, and the dynamic behaviors of the rockets and arresting cables are shown in Figs. 181920, corresponding to the three scenarios. Then, the three parameters (i.e., the maximum buffering distance jz max j, the maximum acceleration of rocket a max , and the maximum stress of the cables σ max ) in Table 6 verify that all 3 × 100 cases were successfully stopped because they met the three success criteria. The landing success of the proposed system includes the success of the active-tracking/capturing phase and the passivedecelerating phase. The success rates of statistical simulation cases are shown in Table 7 for the three scenarios and validate the robustness of the proposed system under various combinations of multiple uncertainties.

C. Robust Performance of the Passive-Decelerating Phase

To verify the robust capability of the proposed system at the passive-decelerating phase, decelerating processes were simulated for various different initial conditions when the rocket was arrested by the capture frame. The ranges of these initial conditions covered the final states of the rockets when the first active-tracking phase ended, as shown in Table 4. The expanded ranges, shown in Table 8, were selected to initiate 100 quasi Monte Carlo simulation cases for each of the three scenarios so that a total of 300 cases were simulated.

The initial conditions of these 100 different cases are shown in Fig. 21a.

For these 300 examples, the statistical results of the three criteria parameters (i.e., the maximum buffering distance jz max j, the maximum acceleration of the rocket a max , and the maximum stress of the cables σ max ) are shown in Table 9. The simulated results show that these three parameters all meet the success criteria of the passivedecelerating phase; therefore, the rocket can be successfully landed. The arresting-cable system can achieve robust deceleration with different landing states of rockets.

In addition, the distributions of the three criteria in the arresting plane oxyz are shown in Figs. 21b-21d, respectively. For the maximum buffering distances kz max k, its maximum value appears in the case of central capturing and, the farther away from the center of the arresting area, the smaller its value was. Thus, the height of the truss frame can be designed to refer to the maximum buffering distances in the central landing case. The maximum accelerations of the rocket a max and maximum stress of the cables σ max were not highly correlated to the captured position deviation. The effects of the wind upon the buffering performance of the system were negligible. The results also supported the robustness of the proposed system.

V. Conclusions

In this study, a robust terminal-landing system was proposed to safely land rockets under various deviations from ideal landing conditions, including uncertainties in the system parameters, wind excitation, and initial states. The system consists of two parts: one on board and one on ground. The onboard system consists of four hooks on the rocket designed to catch the on-ground arresting cables. The on-ground system consists of an active-arresting-cable system comprising four movable arresting cables forming a capture frame and four buffer devices; it is responsible for actively catching the hooks and then decelerating the rocket.

A flexible multibody model was built to evaluate the interaction dynamics between the landing rocket and the arresting-cable system. Based on this model, quasi Monte Carlo simulations confirmed that a rocket can be safely landed even if its landing state deviates from the ideal landing state. The proposed system does not need additional complex onboard equipment other than retractable hooks. The advantages of the proposed system include a potentially higher Table 6 The decelerating capability of the arresting-cable system under the rocket initial states in success rate for vertical landings and less propellant usage for landings, since the system can manage a higher terminal vertical-landing velocity. Another potential advantage is the acceptance of greater wind excitation.

Appendix: Formulation of the Contact Force Between the Hook and the Arresting Cable

The expression for a normal collision force is as follows:

f n kδ e c _ δ
where δ is the penetration depth of the collision point; _ δ is the change rate of this penetration depth; k is the stiffness coefficient; c is the damping coefficient of the collision; and e is the index of the nonlinear collision force, which is related to the shape of the contact surface.

The tangential frictional force is calculated based on a modified Coulomb's frictional force, and it is given by f τ μvf n Here, the friction coefficient μ is expressed as STEPt; t 0 ; h 0 ; t 1 ; h 1 8 > > < > > :

h 0 t ≤ t 0 h 0 aη 2 3 -2η t 0 < t < t 1 h 1 t ≥ t 1
with a h 1h 0 , η tt 0 ∕t 1t 0 .

Table 9 The decelerating capability of the arresting-cable system under the rocket's initial states in

0 =

 0 the arresting cable, m 2 a p = proportional gain c h = damping coefficient of the spring damper, N ⋅ s∕m d h = width of a hook, m E = modulus of the arresting cable, GPa F f = sliding friction, N H = height of the supporting truss frame, m h h = size of a hook, m J engine k = moment of inertia of the rocket's engine relative to OXYZ (k xx, yy), kg ⋅ m 2 J engine zz = polar moment of inertia of the rocket's engine relative to OXYZ, kg ⋅ m 2 k h = stiffness coefficient of the spring damper, N∕m L c = length from the center of mass to the top of the rocket stage, m L h = length from the origin of OXYZ to the top of the rocket stage, m L r = length of the rocket stage, m l h = length of a hook, m l 0 = half-width of arresting area, m l 1 = half-width of capture area, m m engine = mass of rocket's engine, ton m h = mass of a hook, kg m w = mass of a counterweight, ton R p = radius of the active pulley, m R r = radius of the rocket stage, m v p i = reelout velocity of the active pulley i, m∕s v wg = stationary wind speed, m∕s α = transfer rate of the pulley group, dimensionless β = the stiffness proportional to the Rayleigh damping factor θ = angle between the arresting cable and the horizontal line in the simplified theoretical model, rad θ h = initial deployment angle between a hook and the rocket stage, rad θ w = wind angle, rad λ wg = ratio of root mean square of gust speed with respect to stationary wind speed, % ρ = density of the arresting cable, kg∕m 3 σ r = root mean square of measurement noise, m τ = time delay of servo motor, s

Fig. 1

 1 Fig. 1 Sketch of the proposed arresting-cable system.

Fig. 2

 2 Fig.2Sketches of the two phases of the proposed arresting-cable system.

Fig. 3

 3 Fig. 3 Design of the active-tracking/capturing control system.

  κt

) Displacement response of slider 1 b) Reelout velocity of active pulley 1 Fig. 4

 114 Fig.[START_REF] Inatani | Concept and Preliminary Flight Testing of a Fully Reusable Rocket Vehicle[END_REF] Controlling results of the slider's displacement and motor's reelout velocity with and without the saturation function κt.

Fig. 5

 5 Fig. 5 Schematic of the quarter-model of the arresting-cable system.

a)Fig. 7

 7 Fig.7Simulation model of the active arresting-cable system.

Fig. 8

 8 Fig.8Simulation model of the ALE cable element.

Fig. 9

 9 Fig.9Simulation model of the rocket.

Fig. 11

 11 Fig. 11 Performance of the decelerating process: a) snapshots of the simulated solutions, and b-e) comparisons of the numerical and analytical solutions.

Fig. 12

 12 Fig. 12 Topology diagram of the dynamics and control loop of the system.

a)

  Fig. 13 Snapshots and trajectories of a deviated landing case.

  Fig. 15 State deviations, and tracking errors in scenario I.

  Fig. 17 State deviations, trajectories, and tracking errors in scenario III.

a)Fig. 18

 18 Fig.[START_REF] Ecker | Aerothermal Analysis of Reusable Launcher Systems During Retro-Propulsion Reentry and Landing[END_REF] Dynamic behaviors of the rocket and arresting cables after the rocket was captured in scenario I.

Fig. 20

 20 Fig. 20 Dynamic behaviors of the rocket and arresting cables after the rocket was captured in scenario III.

  Fig.[START_REF] Huang | Control Strategy of Launch Vehicle and Lander with Adaptive Landing Gear for Sloped Landing[END_REF] Performance of the decelerating phase in the 3 × 100 landing cases of the quasi Monte Carlo simulations.

  μ d jvj > v d signvSTEPjvj; v d ; μ d ; v s ; μ s v s ≤ jvj ≤ v d STEPv; -v s ; -μ s ; v s ; μ s -v s < v < v swhere v is the relative tangential velocity of the two contacting bodies; v s and v d are the static and dynamic friction transition velocities, respectively; and μ s and μ d are the corresponding static and dynamic friction coefficients. The signv function extracts the sign of the velocity v. The STEP function is a second-order derivation continuous function defined as follows:

Table 1

 1 Physical parameters of the proposed terminal-landing system

	Item		Value
		Rocket stage
	R r		2.00 m
	L r		40.00 m
	L c		30.00 m
	m engine		12.50 tons
	J engine zz		2.50 × 10 4 kg ⋅ m 2
	J engine xx	J engine yy	4.16 × 10 3 kg ⋅ m 2
	EA		2.63 × 10 9 N
	ρA		501.40 kg∕m
	EJ		5.24 × 10 9 N∕m
			Hook
	m h		45.00 kg
	l h		3.50 m
	h h		1.00 m
	θ h		π∕3 rad
	d h		0.15 m
	k h		1.50 × 10 5 N∕m
	c h		4.00 × 10 5 N ⋅ s∕m
	L h		20.00 m
		Arresting-cable system
	H		≥ 30.00 m
	2l 0		50.00 m
	2l 1		7.00 m
	EA		2.64 × 10 8 N
	ρA		9.80 kg∕m
	R p		0.10 m
	m w		0.92 ton
	F

f 3.13 × 10 5 N α 3

Table 2

 2 Deviation ranges of the uncertain parameters of the simulated system

	Numbers	Parameters	Range [min max] Distribution	Description
	1-2	r x ; r y , m	r min r max	Halton	Lateral position deviation
	3-4	v x ; v y , m∕s	v min v max	Halton	Lateral velocity deviation
	5-6	θ x ; θ y , deg	θ min i	θ max i	Halton	Inclination angle deviation
	7	θ z , deg	θ min r	θ max r	Halton	Roll angle deviation
	8-10	ω x ; ω y ; ω z , deg ∕s	ω min ω max	Halton	Angular velocity deviation
	11	θ w , deg	0 θ max w	Halton	Deviation of wind angle
	12	v ws , m∕s	0 v max ws	Halton	Stationary wind speed deviation
	13	λ wg , %	λ min wg	λ max wg	Halton	Fluctuation intensity deviation of gust
	14	τ, s	0 τ max	Halton	Time delay
	15	σ r , m	0 σ max r	Halton	RMS of measurement noise

Table 3

 3 The deviation ranges of the rocket's terminal-guidance states

	Numbers	Parameters	Min Max
	1-2	r x ; r y , m	-13.0 13.0
	3-4	v x ; v y , m∕s	-0.1	0.1
	5-6	θ x ; θ y , deg	-5.0	5.0
	7	θ z , deg	-25.0 25.0
	8-9	ω x ; ω y , deg ∕s	-0.1	0.1
	10	ω z , deg ∕s	-0.1	0.1

Table 4

 4 Description of three scenarios for robustness evaluation

		Three scenarios for active-tracking/ capturing phase	Three scenarios for passive-decelerating phase
	Parameter	I	II	III	I	II	III
	v max ws , m∕s	5	10	20	5	10	20
	h, m	200	120	100	--	--	--
	t c , s	10	6	5	--	--	--
	a P	2	5	6	--	--	--
	Wind excitation		θ max w 360 deg, λ min wg 10%, λ max wg 20%
	Rocket initial states	Terminal-guidance states in Table

Table 5

 5 Deviation ranges of the rockets' captured states after the active-tracking/capturing phase

			Scenario I	Scenario II	Scenario III
	Numbers	Parameters	Min	Max	Min	Max	Min	Max
	1-2	r x ; r y , m	-16.84 16.94 -15.90 16.27 -15.70 16.97
	3-4	v x ; v y , m∕s	-0.97 0.74 -1.43 1.09 -1.86 1.59
	5-6	θ x ; θ y , deg -15.72 15.90 -12.89 12.70 -14.44 15.26
	7	θ z , deg	-26.06 23.42 -25.59 23.14 -25.41 23.67
	8-9	ω x ; ω y , deg ∕s -2.37 3.02 -4.07 4.77 -5.82 5.98
	10	ω z , deg ∕s	-0.22 0.19 -0.34 0.63 -0.49 0.46

Table 5

 5 MPa 505.11 617.98 563.52 524.64 604.12 563.72 534.74 615.33 564.60

				Scenario I			Scenario II			Scenario III
	Numbers	Parameters	Min	Max	Mean	Min	Max	Mean	Min	Max	Mean
	1	jz max j, m	7.96	9.63	9.16	7.96	9.63	9.16	8.25	9.59	9.12
	2	a max , 9.8 m∕s 2	3.58	4.71	4.05	3.58	4.71	4.04	3.51	4.88	4.21
	3	σ max ,

Table 7

 7 Statistics of the quasi Monte Carlo cases

	Numbers	Scenario I Scenario II Scenario III
	Simulation cases	100	100	100
	Successfully captured cases	100	100	100
	Successfully stopped cases	100	100	100
	Success rate, %	100	100	100

Table 8

 8 The deviation ranges of the rocket's captured states for the passive-decelerating phase

	Numbers	Parameters	Min	Max
	1-2	r x ; r y , m	-17.0	17.0
	1-2	v x ; v y , m∕s	-2.0	2.0
	3-4	θ x ; θ y , deg	-16.0	16.0
	5-6	θ z , deg	-26.0	26.0
	7	ω

x ; ω y , deg ∕s -6.

Table 8

 8 MPa 544.79 642.41 584.83 544.79 642.87 584.50 544.79 644.48 584.81

				Scenario I			Scenario II			Scenario III
	Numbers	Parameters	Min	Max	Mean	Min	Max	Mean	Min	Max	Mean
	1	jz max j, m	7.99	9.78	8.80	7.99	9.79	8.83	8.00	9.77	8.82
	2	a max , 9.8 m∕s 2	2.96	4.33	3.56	2.98	4.71	3.56	2.81	4.71	3.54
	3	σ max ,

Acknowledgments

This research was supported by the National Natural Science Foundation of China (grant nos. 11872221 and 11302114) and the Major State Basic Research Development Program (grant no. 2012CB821203). All the authors acknowledge Yongpeng Gu and Jianqiao Guo from Tsinghua University, as well as Daniel Alazard (AIAA Associate Fellow) from the University of Toulouse for several useful discussions. Many thanks are extended to the Editors and Reviewers of this paper.