Asymptotic tracking position control with active oscillation damping of a multibody Mars vehicle using two artificial augmentation approaches

The Valles Marineris Explorer Cooperative Swarm navigation, Mission and Control research project aims to explore the Valles Marineris canyon system on Mars with, among others, multibody rotary-wing unmanned aerial vehicles (UAVs) comprising of a hexrotor system and a helium-filled balloon being attached to it by means of a rope. In this paper, we develop a high-fidelity closed-loop control system in MATLAB® and Simulink™ to present the application of an adequate flight controller guaranteeing (1) asymptotic tracking position control of the multibody flight system, (2) suppression of the balloon’s swinging motion in forward flight case, and (3) stabilization of the rope angle around its equilibrium for steady-state conditions. Applying feedback linearization for the outer loop and analytical backstepping for the inner loop of a nonlinear cascaded control design model of the hexrotor system, we propose an extension of the baseline flight controller by two artificial augmentation approaches to cope with the balloon dynamics. Basically, by utilizing oscillation damping feedbacks of the uncertain plant which are applied as additional commands to either the inner or the outer loop’s reference model. Simulation results are presented for an eight-shaped flight maneuver at the bottom of Valles Marineris proving that the augmentation units increase the flight controller capabilities to suppress modeling errors artificially—without changing the baseline control laws. The augmentation units actively damp the balloon motion in the forward flight case for non-steady-state conditions to counteract the rope oscillations and finally stabilize the rope angle around its equilibrium, so that the Mars vehicle is able to reach a steady-state in position when its extraterrestrial mission profile is successfully completed.

Cross product of two vectors • Dot product of two vectors (⋅) −1 Inverse of a matrix (⋅) + Moore-Penrose inverse of a matrix (⋅) T Transpose of a matrix (⋅) Creation of a diagonal matrix d(⋅)

Introduction
Deeper questions that arise in countless scientific reports [1,2] or [3] deal with the origin and evolution of life: How did life begin and where did it originate if not on Earth? These are questions that perplexed scientists for a long time and, eventually on account of that, triggered a wide search for extinct organisms on other planets in our solar system. The most promising terrestrial planet to solve this gap in understandings of humanity is Mars. Valles Marineris Explorer Cooperative Swarm navigation, Mission and Control, in short called VaMEx-CoSMiC or VaMEx, is a German research project which was started in 2012 to focus on the exploration of the red planet or rather the Valles Marineris. Until today, the Valles Marineris is stated to be the largest canyon system in our solar system providing geological peak values of ~3000 km length, ~600 km width, and ~11 km depth in total. It is located in the southern hemisphere of the planet parallel to its equatorial axis and expands roughly from 0° to 20° S and from 50° to 90° W excluding the Tharsis region in the West. Due to the deep cliffs and fractured rock layers, geologists are hoping to reveal new insights into the history of Mars and find evidence for microbial life or even petrified microorganism being hidden into the Martian crust. The key purpose of VaMEx is therefore to develop, compare, and evaluate technologies for an unmanned autonomous robotic swarm mission exploring the Valles Marineris. Within the project, the Institute of Flight System Dynamics (FSD, TU Munich) is, inter alia, responsible for the design of rotarywing platforms being used as swarm participants and their embedment into a high-fidelity simulation environment. This incorporates a model-based plant design of the flying platform and the derivation of a suitable flight controller both forming the closed-loop control system [4][5][6].
Regarding the thin atmosphere on Mars, the development of a fully operational multicopter concept is quite a difficult task. Since all propellers' require high-density values to produce sufficient thrust, only the deepest points are investigated for an extraterrestrial mission profile leading to a geographic location of 14.035° S and 58.5° W at the bottom of Valles Marineris with an average height h avg of −4907 m below the Martian geoid. At the FSD, conceptual studies for this given design point were conducted by [7] to determine a first suitable solution for a Martian rotary-wing unmanned aerial vehicle (UAV) being depicted in Fig. 1. A MATLAB ® -based scalable, parametric design tool [8] was invented afterwards to reevaluate different vehicle configurations. As an ultimate version, a 6-rotor UAV with a balloon being attached to it by means of a flexible rope was chosen as most appropriate solution. The balloon provides additional buoyancy for the vehicle leading to an enhancement of the actuators' battery lifetime and thus to a total flight time improvement.

Literature overview
In the literature, several conceptual studies besides VaMEx exist to determine the optimal design of a flying platform on Mars as in [9]. In addition to rotary-wing UAV concepts, the Martian Autonomous Rotary-Wing Vehicle (MARV) [10], the Mars Helicopter Scout [11], or a meso-scaled 4-rotors UAV concept called Mesicopter [12] were published among many other concepts.
Flight control approaches for those Martian UAVs as well as their dynamic modeling into a simulation environment are Fig. 1 VaMEx project emblem [6] rare to find in literature-not to mention multibody rotarywing concepts like the proposed Mars vehicle of Sect. 1. Only limited data have been published [11,12] which does not vary widely from terrestrial procedures. Based on these findings, we conclude that terrestrial methods for modeling UAV dynamics, incorporated control strategies, and their simulative setup, e.g., within MATLAB ® and Simulink™, can likewise be used to build a virtual flight control system for Martian rotorcraft. As state-of-the-art simulation framework, the FSD Simulation Environment (FSDSE), which incorporates a high-fidelity six degree of freedom rigid-body model, is studied and used to develop a Martian simulation environment which meets the requirements of the VaMEx-CoSMiC project. Successful applications of the FSDSE are, e.g., given in [13,14] and a comprehensive model description is given in [15]. A summary of commonly used approaches to design, model, and control terrestrial hybrid UAVs was conducted by [16]. Our main distinctions within the stated methods are reflected by three key issues: First, we address the embedment of a realistic atmospheric model which is built on high-resolution atmospheric measurement data being smoothed and prepared by a meteorological database of Mars to best fit 14.035° S and 58.5° W, the vehicle's geographic take-off location at Valles Marineris. Due to the low Reynolds numbers, propellers are not able to rotate as fast as on Earth which was also the most challenging factor presented in [10].
Second, we propose a procedure to develop a nine degrees of freedom multibody UAV model which relies on the classical Newtonian formulations in which a six degrees of freedom rigid-body model for the multirotor system and a three degrees of freedom point-mass model for the balloon are afflicted with one kinematic constraint each. The suspension line, or rather rope, is enforced with one kinematic constraint at each coupling point being obtained by a massless springdamper oscillator model. This model significantly aids in the design of multibody dynamics which are, in general, derived with analytical mechanics leading to complex first-or second-order differential equations. For instance, the governing motions for comparable dynamic models as, e.g., a multibody Parafoil-UAV are obtained using either Lagrange's equations, shown in [17], or the Hamiltonian procedure, published in [18]. In [19], Lagrange's equations are likewise used to obtain the dynamic model system formulations of a multibody quadrotor UAV with cable-suspended payload.
Third, and main contribution of this paper, we present a flight controller which is capable to handle multibody flight dynamics guaranteeing (i) Asymptotic tracking position control of the multirotor system, while damping the suspension line oscillations simultaneously.
(ii) Suppression of the balloon's swinging motion in forward flight case. (iii) Stabilization of rope angle around its equilibrium for steady-state conditions.
Among a greater diversity of nonlinear flight controllers, Backstepping [20] and Nonlinear Dynamic Inversion (NDI) [21] are two of the most widely used control design methodologies for agile multirotor systems. Examples of successful applications are, e.g., given in [14,[22][23][24][25][26]. Essentially, Backstepping is a model-based approach to obtain, based on control-Lyapunov functions [20], asymptotical closedloop stability of a feedback control system which provides a so-called strict feedback form [20], whereas NDI is a state transformation which ensures, without any approximations, linear input-output dynamics of a nonlinear system. It is therefore often called Feedback Linearization [21]. A comparison between NDI and analytical backstepping was conducted by [27] showing that both methods ensure the same satisfying controller performance for agile multirotor systems. Thus, in this paper, we propose a nonlinear flight controller which is a cascaded version of a backsteppingbased attitude controller nestling inside a superior position controller, which is based on NDI, to meet objective (i).
The idea of both backstepping and NDI methodologies is, in general, to make a dynamic system or rather its state trajectory follow a desired reference trajectory by determining an appropriate control law. In case of uncertain dynamic systems, however, modeling errors are inevitable within the control design model which leads to a mismatch between the state and reference trajectory. A common approach to cope with modeling errors is to augment the baseline flight controller architecture for addressing secondary controller capabilities as stated in (ii) and (iii).
The first, and perhaps most intuitive, approach is to augment the baseline flight controller by input shapers [28] which are, e.g., used in multibody quadrotor UAVs with suspended load [29] or in helicopter slung load systems [30] to damp the swinging motion. Input shaping is a powerful feedforward technique where the natural frequency and damping ratios are estimated from the linearized well-known plant dynamics to determine an oscillation damping command being applied to the reference model. Since the actual reference model remains open-loop, robustness of this technique is often achieved by adding a delayed feedback of the uncertain plant to the closed-loop control system which is related to the theory of delayed feedback control [31]. This feedback is then added to the feedforward reference trajectory which results in an active vibration-damping and also delayed reference signal. Similar to the input shaping technique, control design parameters have to be determined from a linearized state-space model, where the actual time-delay needs to be modeled. Examples of successful applications are, e.g., given in [32] for container cranes or in [31,33,34] for helicopter slung load systems to actively damp the swing motion. In [24], a comparison between the usage of input shaping, delayed feedback control, and a combination of both techniques showed that the input shaping is indeed advantageous in suppressing unintended swing motions for a multibody hexrotor UAV with cable-suspended load, but also, that a delayed feedback controller pose a powerful standalone augmentation approach to address the controller capabilities (ii) and (iii).
In this paper, we summarize the nonlinear dynamic effects of the plant, being represented by oscillating motion of the suspension line, as uncertainties and suppress them by augmenting either the inner loop or outer loop's reference model with oscillation damping feedbacks which are applied as additional commands to solve objective (ii) and (iii). Although this method is in the spirit of delayed feedback controllers [31], the augmentation signals are not subjected to time-delays. Additionally, we directly augment the input signals of the reference model which generates a closed-loop including both the uncertain plant and the reference model. The proposed artificial augmentation units are therefore assigned to the classical theory of closed-loop reference models [35,36] which have their origin in the framework of adaptive control [37,38]. In literature, augmentation approaches for multibody UAVs which explicitly use the theory of closed-loop reference models to satisfy the controller's capabilities (ii) and (iii) could not have been identified.

Outline
The remainder of this paper is organized as follows: The plant design of the multibody Mars vehicle is presented in Sect. 2. Section 3 reviews the nonlinear control design model and derives the baseline flight controller. Two oscillation damping feedback augmentation approaches are presented in Sect. 4, both extending the baseline flight controller separately. Simulation results are stated in Sect. 5 to validate the control design including the artificial augmentation units. Section 6 summarizes the stated findings and provides a short conclusion.

Plant design
The Mars vehicle is modularly implemented into MATLAB ® and the Simulink™ toolbox as multibody UAV forming the plant of the closed-loop control system. It consists of two nonlinear UAV models, a six degrees of freedom (6DoF) rigid-body model for the multicopter and a 3DoF pointmass model for the balloon. Both models are additionally afflicted with kinematic constraints to not only generate a linked connection, but also to support a modular design process. The general structure is depicted in Fig. 2. The highfidelity simulation model consists of four main submodels, the environment submodel, the equations of motion (EoM), the multibody system formalism, and the airframe submodel. Its baseline architecture is inspired by the FSDSE [15].
Using the state-space model representation, the multirotor system is input affine in the system's commanded input vector, containing the square of the actuators' commanded rotational speeds c,i ∈ ℝ . Its system dynamics can be described by as well as the autonomous balloon system dynamics are described by To uniquely describe the current situation of both dynamic systems for an arbitrary point in time, the multicopter state vector x ∈ ℝ 12 includes the position r R ∈ ℝ 3 , the kinematic velocity v R ∈ ℝ 3 , the angular rates IB ∈ ℝ 3 , and the Euler Angles ∈ ℝ 3 as attitude representation. The balloon state vector x Bal ∈ ℝ 6 only includes the position r Q ∈ ℝ 3 and the kinematic velocity v Q ∈ ℝ 3 . Summarizing both state vectors to x P = x x Bal T , the total plant state vector, it is possible to write the plant system dynamics to so that the system output equation is given by x Bal = f x, x Bal . To determine the current plant state vector x P (t + ∆t) only by means of the last known plant state vector x P (t), Simulink™ provides a broad range of higher order ordinary differential equation (Ode) solvers like Ode4 which is based on the fourth-order Runge-Kutta method [39]. In contrast to first-order Ode solvers, where Euler's integration [39] serves as a base, the Ode4 solver determines an averaging value for the plant state vector's derivative ̇x P (t) at exactly four weighted sampling points. This not only ensures a higher accuracy to solve, respectively, propagate the UAVs' nonlinear state equations numerically, but also upholds the fact that ̇x P (t) is usually not constant during a time increment of ∆t. To build a time-history simulation [39], Ode4 is used as numerical integration method with a fixed sampling time ∆t of 0.001 s to step the plant state vector forward.

Environment submodel
To ensure that the multibody UAV can be simulated under realistic environmental conditions of planet Mars, a Martian Standard Atmosphere (MSA) is implemented inside the environment submodel which is related to the multicopter reference point R.
This atmospheric model is built on high-resolution atmospheric measurement data being smoothed and prepared by the Martian climate database (MCD) V5.2 [40]. For altitude definitions, a Martian reference ellipsoid with a geocentric radius r M of 3394.6 km is established which is based on the concept of the World Geodetic System 84 [41]. The obtained ellipsoid best fits the geographic location of 14.035° S and 58.5° W at the Valles Marineris, while its semi-major axis a equals 3396.2 km and the square of its first eccentricity e 2 equals 0.0117. This reference surface also defines the normal gravity potential g 0 ∶= 3.717 m∕s 2 being required for the 1D quadratic gravitation model [42,43].
In addition to seasonal and diurnal changes, the MSA is related to a Martian solar longitude of 359 • and a local true solar time of 12 ∶ 00 in Mars year 24. Regarding the vertical extension of the Valles Marineris, the crucial part of the MSA forms the Martian troposphere which can be considered as real polytropic up to a height of ~ 6.4 km above the ground. Meteorological studies yield that the decrease of atmospheric temperature with altitude can be represented by two major dry adiabatic lapse rates [44] (DALRs) [40].
very precisely, so that the Martian troposphere is divided into two layers, the lower Martian troposphere (index l) and the upper Martian troposphere (index up). Due to both constant DALRs, the barometric formula [42] can be solved analytically which, in terms of the ideal gas law [44], leads to all formulas for atmospheric changes relative to the geopotential height [42] where h R is the geodetic height of the multirotor system above the reference ellipsoid. A summary is listed in Tables 1 and 2 where T, p , and denoting the atmospheric temperature, pressure, and density as well as represents the dynamic viscosity on Mars. A proof, that the atmospheric model maps the MSA V5.2 measurement data with high accuracy is shown in Fig. 3. The only rough distinctions can be identified at ground level of Valles Marineris, due to the temperature ground effect as a result of cosmic radiation on Mars, and for heights beyond the upper Martian troposphere. These model uncertainties of the MSA are excluded within the flight envelope of the Mars vehicle, since simulation results are only evaluated from an average height h avg of −4907 m below the Martian geoid. Therefore, they can be neglected for further investigations.

Equations of motion
For describing the state equations of the multirotor system, we use a body-fixed (B) frame and a North-East-Down (I) frame as reference frames being depicted in Fig. 4. The multicopter reference point R denotes the origin of the B-frame and is assumed to be congruent with its center of gravity G. This yields decoupled differential equations for the translation and rotation dynamics. The I-frame is located on the surface of Valles Marineris and has a nonrelocatable placement at the multirotor system's point of departure. Its x-y-plane is parallel to the local surface whereby the x I -axes points to the Martian north pole. By assuming that planet Mars is non-rotating and flat, without any elliptical shape, the I-frame is considered as inertial and can be used to apply Newton's second law of motion [39]. The state propagation equations for position, translation, rotation, and attitude can thus be formulated by whereby I ∈ ℝ 3x3 is the identity matrix and H( ) ∈ ℝ 3x3 is given by In Eq. (8), m = 2.317 kg denotes the total mass of the multirotor system including 0.2 kg of payload, F R ∈ ℝ 3 and M R ∈ ℝ 3 the total amount of all external forces and moments related to R, and J R represents the multicopter's mass moment of inertia tensor which is approximately given by diag 0.131 kg m 2 0.131 kg m 2 0.261 kg m 2 ∈ ℝ 3x3 . To describe the attitude of the multirotor system, represents the pitch angle, and Φ ∈ [− ; ] symbolizes the roll angle. Together, they constitute the rotation matrix which maps a vector from the B-frame into the I-frame. To avoid the inherent singularity of Eq. (9) for Θ = ± π∕2 , the simulation model can be switched-based on the configuration values (see Fig. 2)-to quaternions [15] being used as 4-dimensional representation to describe the orientation of the multirotor system. In addition to the environment submodel of Sect. 2.1, the multicopter's differential equation for its geodetic height is given by where ̇h R denotes the time derivative of the geodetic height with respect to the I-frame being required for altitude definitions above the reference ellipsoid.
To describe the sate propagation equations of the balloon point-mass model, we choose Q as reference point and use a second body-fixed coordinate frame, the B Bal -frame, to formulate its dynamics. The state propagation equations for position and translation are thus given by where m Bal = 0.06 kg denotes the balloon mass and F Q ∈ ℝ 3 the total amount of all external forces related to Q. In relation to Sect. 2.1, the balloon's differential equation to propagate its geodetic height h Q can be formulated to where ̇h Q denotes the time derivative of h Q with respect to the I-frame being required for altitude definitions above the reference ellipsoid.

Multibody system formalism
Complex first-or second-order differential equations do usually occur while working with multibody dynamics. These equations do not only require a lot of computing power, it is also much more difficult to maintain an analytical formulation using, e.g., the theory of virtual work [45] or the Lagrange equations [45], respectively, solving them. In this paper, we use the classical Newtonian formulations to impose kinematic constraints on both dynamic systems ̇x and ̇x Bal as massless spring-damper oscillator (SDO) model to approximate a realistic rope connection. This SDO model triggers F R rp ∈ ℝ 3 and F Q rp ∈ ℝ 3 , the rope reaction forces, and adds them-as part of the airframe submodel-to the total amount of applied forces F R and F Q acting on both dynamic systems. For describing the kinematics of the SDO model, the multicopter rope (R * ) frame as well as the balloon rope (R * Bal ) frame are introduced as reference frames. Both R * -and R * Bal -frame are moving freely with the UAVs' reference point while serving the following properties (see Fig. 4):

Bal
-axis and the z R * ∕z R * Bal -axis forming a righthand coordinate system.
To describe the orientation of the resulting rope force relative to the B-and B Bal -frame, R * and R * , respectively, , are defined as rope angles. Together, they constitute the rotation matrix are given through the geometric vector chains and (13) whereby I ∈ ℝ 3x3 is the identity matrix. The multicopter rope angles can thus be calculated to as well as the balloon rope angles to

Design of spring force
To design the spring force for the SDO model, it is assumed that a spring produces a restoring force proportional to its deflection l . In case of the multicopter, l can be calculated through ‖ � r RQ � B ‖ as well as l is determined by ‖ � r QR � Bal ‖ for the balloon UAV. Using Eqs. (13) and (14), this leads to a spring force of [45] for the multirotor system and likewise for the balloon, where l 0 = 5 m implies the virtual spring's norm length. Since the rope should act in a naturally flexible manner, the spring design may not be too stiff. The virtual spring stiffness k sp is therefore chosen to equal 0.5 kg∕s 2 which will, in the end, confer a smooth oscillating behavior to the SDO model. Regarding Eqs. (17) and (18), , or rahter, Bal are both logical operators ensuring that only positive spring forces are transmitted by the rope. Otherwise, Fig. 4 System architecture of multibody Mars vehicle with SDO model and I-frame pushing forces would act on both UAVs for l being smaller than l 0 which is not desired. Thus can be defined for the multicopter as well as for the balloon UAV.

Design of damper force
To design the damping force for the SDO model, it is assumed that a damper always produces a force acting in opposite direction of its movement [45]. The damping force depends therefore on the UAVs' change of position over time in rope direction which are given by for the multirotor system and for the balloon. Otherwise, kinematic velocity elements of multicopter or balloon in non-rope direction would wrongly lead to a damping force which is not desired. Based on the R * -and R * Bal -frame properties and coordinate definitions, the damping force of the SDO model can be formulated to for the multicopter and likewise for the balloon, where , or rather, Bal are already predefined by Eqs. (19) and (20) ensuring that damping forces only occur, while the rope is oscillating due to a positive spring force. Furthermore, d dp implies the virtual damping constant which is chosen to equal 0.75 kg∕s . This choice of parameter is synergetic with the spring design of Sect. 2.3.1, so that a smooth oscillating behavior is damped in a short time period.

Rope force of SDO model
To design the rope force for the SDO model we simply use Eqs. (17) and (23), or rather, (18) and (24) and add them as part of a vector summation. In addition to the multicopter, it directly follows: for its rope force being notated in the B-frame, as well as for the balloon's rope force being notated in the B Bal -frame.

Airframe submodel
The airframe is the most extensive submodel in the plant. Its outputs are strongly interrelated to the EoM submodel of Sect. 2.2, since the state equations of multicopter and balloon can only be propagated with the total sum of all external forces and moments acting on the UAVs. The airframe's task is to provide and prepare them. As the Martian density is most crucial to produce buoyancy for the balloon UAV, or rather, has a linear dependency to the thrust equation of the propellers, the airframe submodel is mainly depending on the MSA being embedded in the environment submodel (see Sect. 2.1).
In case of the multirotor system, all external forces and moments are given by and being subdivided into their physical types of origin which are propulsion (P), aerodynamics (A), gravity (g), ground contact (C), and rope (rp). To formulate the multicopter's gravity force, we use the 1D quadratic gravitation model [42] Since R is assumed to be congruent with the multicopter's center of gravity, the gravitational force as well as the rope force, given by Eq. (25), induce no moments in Eq. (28) and therefore only contribute to the balance of forces in Eq. (27). To design a simplified behavior for the multirotor system's ground contact, four contact points r C i ∈ ℝ 3 , each of them given by the vector summation of r R and r RC i , are defined as landing gear exerting the total contact force as well as the total contact moment on the vehicle if, and only if, the vehicle is in direct contact with the valley surface of Valles Marineris. The particular contact force of each contact point F C i ∈ ℝ 3 is composed of frictional forces [45] in the x B -y B plane and an impact force in z B -direction both counteracting the vehicle's movement.
Since we assume that the landing gear carries the vehicle mass equally, i.e., each landing gear leg carries one quarter of the vehicle's gravitational force, given by Eq. (29), the contact force of each contact point can be formulated as [45] where fr represents the frictional coefficient which equals 100. Applying Euler's classical treatment of vector analysis [39], the kinematic velocity of each contact point is given by The propulsive forces and moments are produced by six propellers. Each of them provides two blades, a calculated density of 240 kg∕m 3 , a blade radius r bl of 0.2055 m , and a mass inertia J zz = 8.9 ⋅ 10 −4 kg∕m 2 around its rotational axis [8]. To avoid an inherent minus sign within the propulsion unit modeling, six propeller-fixed (P i ) frames are introduced, where the z P i -axis is perpendicular to the x B -y B plane and points in direction of positive thrust. The corresponding rotation matrix from the body-fixed B-frame into the particular P i -frame is therefore simply given by 3 and, while using the P i -frame as notation frame, the rotational speeds vector of each propeller can be written as P i Rot = 0 0 act,i ⋅ i T ∈ ℝ 3 depending on act,i , the actuators' predefined direction of rotation with respect to the z P i -axis. For rotor one, three, and five, act,i equals 1 , whereas for rotor two, four and six, act,i equals −1 . Assuming that the position of each propeller r T i ∈ ℝ 3 is given by the vector summation of r R and r RT i , the total propulsive force in Eq. (27) can be formulated to as well as the corresponding total propulsive moment in Eq. (28) can be determined by In Eqs. (34) and (35), the thrust-generating force is the aerodynamic force of the propeller which is indicated by an index A. We derive this force from the Blade Element Momentum Theory [46] while assuming the propeller can be approximated by a circular disc providing an elliptical circulation distribution and therefore a constant induced velocity distribution in radial direction. Based on [46], the aerodynamic force of each propeller can be written to where the thrust coefficient C T equals 0.074 [8]. In Eq. (35) Here, C n denotes the torque coefficient and is given by 0.0164 [8]. Based on Euler's classical treatment of vector analysis [39], the gyroscopic moment can be written to where J T i = diag 0 0 J zz ∈ ℝ 3x3 is the propeller's mass moment of inertia tensor and ̇ P i Rot the propeller's angular acceleration which is related to the actuator model.
A common approach to reproduce the actuator's dynamic behavior is to filter the propeller's angular rate by a PT1 element instead of modeling the actual electric drive. This PT1 filter, based on [47], is given as where t PT1 is the filter's time constant which equals 0.05 s and P i Rot c ∈ ℝ 3 is the actuator's commanded rotational speed vector being related to the vector elements of the system's commanded input u c , given by Eq. (1), through where 2 c,i ∈ u c . As the multirotor system has no lift-generating surfaces, except the thrust-producing propellers, the aerodynamic force in Eq. (27) is solely approximated by the body drag. Based on the static MSA of Sect. 2.1, the aerodynamic velocity of the multicopter is equal to its kinematic velocity. Since the aerodynamic reference point is assumed to be congruent with the multirotor system's center of gravity, or rather R, the aerodynamic force is, based on [46], written to where the dimensionless drag coefficients C Dx and C Dy are assumed to equal −0.06 , respectively C Dz to equal −0.13 . In Eq. (41), the reference area S r is given by 1 m 2 . Similar to the rope and gravitational force, the aerodynamic force induces no moments in Eq. (28).
For the balloon UAV the total sum of all external forces simplifies, in contrast to the multirotor system, to where the rope force is provided through Eq. (26). Since the balloon has no thrust-producing propulsion unit, the aerodynamic force F Q A ∈ ℝ 3 is the solely force acting on the UAV. In Eq. (42), one part of this force is generated by the balloon's buoyancy and the other part by its body drag. For the sake of simplicity, we assume that the balloon can be approximated by a spherical shape and that unheated helium ( R He = 2077 J∕kg K) is used as carrier gas which is less dense than the surrounding atmosphere in the Valles Marineris. Since the modeling of the balloon includes an envelope which generates a gravitational force due to its mass, the balloon's net buoyant force can be formulated to Here, evl = 0.003128 kg∕m 2 is the envelope's areal density which corresponds to a polyethylene hull with an estimated film thickness of 0.0034 mm [48], and S Bal is the balloon surface which is calculated to 15.09 m 2 based on its spherical shape approximation. As per [48], the balloon's buoyant force can be written as (40) where V Bal denotes the total volume of the balloon being predefined by 5.51 m 3 . Based on the static MSA of Sect. 2.1, the aerodynamic velocity of the balloon is equal to its kinematic velocity. Thus, the body drag force in Eq. (42) is, based on [46], calculated to where the dimensionless drag coefficient C D equals −0.2 and the reference area S r,Bal is predefined by 3.77 m 2 .

Control design
The goal of this section is to obtain a simplified mathematical description of the plant, so that the derived model can be used as a basis for the flight control design of the proposed Mars vehicle. In comparison to the high-fidelity simulation model in Sect. 2, this procedure only takes aspects and properties of physical significance into account. Thus, it is assumed that both multirotor system and balloon UAV can be decoupled in terms of flight characteristics so that the balloon dynamics and the SDO model are neglectable within the controller model.
To solve objective (i), the flight controller is designed using a nonlinear cascaded control design model (CDM) which is inspired by [26] and [49].
For describing the main dynamic aspects of the multirotor system, another body-fixed (b) frame is introduced which provides similar properties as the B-frame of Sect. 2.2 and is illustrated in Fig. 5. Referring to Eqs. (36) and (37), it is assumed that all force and moments of the propulsion unit are directly Note that, since only main physical effects are considered, gyroscopic moments are neglected within the CDM. Equations (46) and (47) illustrate that all aerodynamic forces and moments are explicitly related to the commanded input vector of the plant, given by Eq. (1), which contains the square of the actuators' commanded rotational speeds. Since the CDM's output vector v d ∈ ℝ 4 contains the desired virtual control effort and this effort is represented through the spatial distribution of all desired moments M R d and through the total desired force T R d , both acting on the multirotor system's reference point R, v d is directly correlated to the plant's commanded input vector by Note that the mapping between the desired virtual control and the set of the actuators' squared commanded rotational speeds is described by the control matrix B ∈ ℝ 4×6 . Based on Eqs. (46) and (47), this matrix is composed of B M ∈ ℝ 3×6 , which describes the linear mapping for the multicopter's total propulsive moment M R b ∶= B M u c being defined as as well as B T ∈ ℝ 1×6 , which is describing the mapping for the multicopter's total propulsive thrust T R ∶= B T u c being defined as In Eq. (49), l denotes the multicopter arm length which is predefined by 0.4315 m . By assuming a unique mapping similar to quadcopter configurations where v d and the (46) number of actuators have the same dimension ( B −1 is quadratic and does exist), the commanded input vector of the plant can be calculated by introducing the Moore-Penrose pseudoinverse B + ∈ ℝ 6×4 , so that [50] To ensure that the desired virtual control is mapped into the physically attainable set of the actuators' squared commanded rotational speeds, Eq. (51) is further used to introduce the actuators' rate constraints by saturating the vector elements of u c to 2 c,i,min = 0 rad 2 ∕s 2 as well as 2 c,i,max = (6000 ⋅ ∕30) 2 rad 2 ∕s 2 . Euler's first and second law [45] serves as a basis for the controller model to derive the translational and rotational dynamics of the multirotor system which are both related to the I-frame used as an inertial reference frame. The translational dynamics can thus be written to where m = 2.317 kg denotes the multirotor system mass, mg R ∈ ℝ 3 its gravitational force including g 0 , t R ∈ ℝ 3 its total propulsive force, and ̇v R II its absolute acceleration which is the time derivative of the multicopter's kinematic velocity with respect to the I-frame. Applying Euler's classical treatment of vector analysis [39] as well as using the b-frame as notation frame, the rotational dynamics can be written to where J R ∈ ℝ 3×3 represents the multirotor system's mass moment of inertia tensor, which is constantly chosen as diag 0.131 kg m 2 0.131 kg m 2 0.261 kg m 2 . To obtain an attitude representation, Euler angles are introduced-similar to the plant design-where Ψ ∈ [− ; ] represents the azimuth angle, Θ ∈ [− ; ] symbolizes the pitch angle, and Φ ∈ [− ; ] is the roll angle. Together, they constitute the rotation matrix R Ib = R z (−Ψ)R y (−Θ)R x (−Φ) ∈ SO 3 which maps a vector from the b-frame into the I-frame. The first two CDM states are defined by the multicopter's kinematic velocity v R as well as Ib which is the vector of kinematic angular body rates. Their time-depending dynamic behavior is given by Eqs. (52) and (53). To control the position of the multirotor system, r R ∈ ℝ 3 is defined as the CDM's position state whereby its kinematic relation to the multicopter velocity is provided through ̇r R I = v R I . Instead of choosing ∈ ℝ 3 , the vector of Euler angles, as the CDM's attitude state, we chose the multicopter's total propulsive force t R ∈ ℝ 3 as attitude state which can be calculated, while using the b-frame as notation frame, to (51) This thrust vector is composed of the total propulsive thrust T R , given by Eq. (50), and z b = 0 0 1 T , which defines the unit vector of the body z b -axis. To obtain the time derivative of Eq. (54), we use the I-frame as reference frame while applying Euler's classical treatment of vector analysis [39] which yields For a more compact notation of Eq. (55), the thrust vector dynamics can be reformulated into a matrix-vector form of notation [49] which includes T ∈ ℝ 3×3 , the thrustmagnitude matrix, as well as T ∈ ℝ 3×3 , the reduced thrust-magnitude matrix: Note that this approach is based on a reduced attitude parametrization [51] where the multicopter's yaw rate is excluded from the thrust vector dynamics in Eq. (55) or (56). Hence, a rotation around the z b -axis is not correlated to a change of the CDM's attitude state based on Eq. (55) or (56), so that the multirotor system's yaw control can be designed independently.

Baseline controller
Based on the CDM of Sect. 3, a flight controller for the Mars vehicle can be derived which is primarily inspired by [26] and [49]. The flight controller is a cascaded version of a backstepping-based attitude controller nestling inside a superior position controller which is based on Feedback Linearization [52]. Both nonlinear controllers are denoted as baseline controller, since they serve as a base for the artificial augmentation strategies presented in Sect. 4. Figure 6 summarizes the overall controller architecture.
Based on its cascaded structure, the flight controller is divided into an inner attitude loop being enclosed by an outer position loop. Since the dynamics of the inner loop are much faster than the dynamics of the outer loop, we conclude that both loops can be designed independently being justified by the theory of time-scale separation [53]. Inside the position loop, r R and v R are both controlled by (54) stabilizing second-order error dynamics. Inside the attitude loop, t R is controlled using a two-step backstepping control law as well as Ib z is controlled using a one-step backstepping control law. The position as well as the yaw rate controller are fed by r R c and Ib z,c , the pilot's position and yaw rate command, whereas the attitude controller is fed by a commanded thrust vector t R c as feedforward. Since the overall controller testbed is treated as a modelin-the-loop (MIL) simulation, the controller states are initialized through the plant state vector x P . Feasible reference trajectories, indicated by the subscript r, are not only provided by linear first-and second-order reference models (PT1 and PT2 filters), but also by a second-order nonlinear reference model [49], to define the controllers' tracking objectives.

Baseline position controller
To derive the feedback law for the baseline position controller, position error dynamics are stabilized by the outer loop's pseudo-control t R c around its zero equilibrium. The control objective has been reached when r R tracks a desired reference trajectory r R r ∈ ℝ 3 , so that the position tracking error e p ∈ ℝ 3 and its time derivatives converge to zero: We choose the I-frame as notation frame and insert Eq. (52) into (59) yielding being notated in the I-frame. Next, we assume that a total commanded thrust t R c does establish which denotes a desired behavior of ë p . All matrices K v , K x , K i are Hurwitz and of the set ℝ 3×3 . If the total thrust t R is selected according to the pseudo-control law the desired position error dynamics in Eq. (61) and the real position error dynamics in Eq. (60) become equal so that the CDM's position state r R and the translation state v R approach their reference trajectories exponentially fast.

Baseline attitude controller
The control objective of the baseline attitude controller is achieved when t R tracks a desired reference trajectory t R r ∈ ℝ 3 , so that the attitude tracking error e t ∈ ℝ 3 and its time derivative converge to zero: We choose the I-frame as notation frame and insert Eq. (56) into (64) by using the rotation matrix R Ib : Here, ( * ) defines the vector of pseudo-control variables for the thrust vector tracking (TVT). To derive the control law of this first backstepping step, a positive definite Lyapunov function (LF) candidate V 1 [53] is chosen in quadratic form as whereby the time derivative of V 1 ∈ ℝ is calculated to To obtain an exponential decay behavior of ̇e t , Eq. (67) must ensure negative semi-definiteness. Therefore, let T t e t + 1 2 e T tė t = e T tė t .̇e t,d = −K t e t denote a desired decay behavior of the real TVT error dynamics. Then, using Eq. (67) once more a desired time derivative of V 1 can likewise be written to where K t ∈ ℝ 3×3 is a Hurwitz matrix acting as feedback gain of e t . We assume that V 1,d can be generated by a desired pseudo-control input u 1 = u 1x u 1y u 1z T ∈ ℝ 3 , so that for ( * ) = u 1 , the desired time derivative V 1,d and the time derivative of the actual LF V 1 become equal. A comparison of Eqs. (67) and (68), while using the I-frame as notation frame, yields the pseudo-control law for the TVT after inserting Eq. (65) into Eq. (67). Note that, since the thrustmagnitude matrix is quadratic and provides full column and row rank, its inverse T −1 ∈ ℝ 3×3 can be calculated to In Eq. (70), a problem of singularity arises for the inverse thrust-magnitude matrix if the total propulsive thrust T R is equal to zero. To avoid that T −1 becomes singular we saturate T R to a minimum permissible propulsive thrust of 0.1 N . Concerning the general backstepping methodology [53], the error cross-coupling term lacks in Eq. (69) which illustrates the main distinction in deriving the baseline controllers' pseudo-control laws compared to [26]. The attitude loop is stabilized independent of e p being justified by the theory of time-scale separation [53]. The baseline position controller shall lose some of its performance, since the inner loop and therefore also the plant's commanded input vector u c receive less information of the outer loop.
Proceeding with the general backstepping methodology, the comparison between the desired pseudo-control law for the TVT, given by (69), and the current vector of pseudocontrol variables ( * ) (see Eq. (65)) leads to sub-targets of the primary control objective. Since both angular body rates Ib x and Ib y are related to the rotational dynamics, given by Eq. (53), they cannot be manipulated directly in a manner that Ib xy =

Ib x
Ib y T ∈ ℝ 2 matches the desired pseudo-control u 1xy = u 1x u 1y T ∈ ℝ 2 at any time. Hence, to ensure angular rates tracking (ART), the ART error e xy = e x e y T ∈ ℝ 2 and its time derivative with respect to the b-frame must converge to zero: (71) e xy = u 1xy − Ib xy , To satisfy TVT and ART at the same time, the ART error must be included within the TVT error dynamics, so that Eq. (71) is reformulated to Ib xy = u 1xy − e xy and inserted into Eq. (65) which yields Since the actual pseudo-control variable Ṫ R becomes an input for the inner loop while propagating it to a firstorder integrator chain (see Fig. 6), and additionally, since we assume that every demanded, or rather desired, input for the inner loop can be generated in such a short period of time that it is always available, Ṫ R matches its desired pseudo-control u 1z at any time, so that no additional subordinated control error arises. For the sake of clarity, we rewrite Eq. (73) to To derive the control law of this second backstepping step, the valid LF candidate V 1 , given by Eq. (66), is extended to whereby the time derivative of V 2 ∈ ℝ can be determined to We use the I-frame to notate the TVT error and its dynamics as well as the b-frame to notate the ART error and its dynamics. (73) (78) retransformed, so that the ART error is explicitly included in the vector space ℝ 2 . For this purpose, a weighted matrix W ∈ ℝ 3×2 is introduced which satisfies that e x e y 0 T b can be replaced by W e xy b . Since V 2 , respectively its time derivative, is a scalar quantity, we transpose the error crosscoupling term e t T I R Ib TW e xy b and rewrite Eq. (78) to Note that Eq. (79) includes the angular body accelerations ̇I b x and ̇I b y which are directly affected by the rotational dynamics, given by Eq. (53). Thus, ̇ Ib xy defines the vector of pseudo-control variables for the ART. To render Eq. (79) negative semi-definite, which proofs stability of both tracking errors e t and e xy around their zero equilibria 0 ∈ ℝ 3 and 0 ∈ ℝ 2 , let ̇e xy,d = −K xy e xy denote a desired decay behavior of the real ART error dynamics. Then, using Eq. (77) once more, a desired time derivative of V 2 can likewise be written to where K xy ∈ ℝ 2×2 is a Hurwitz matrix acting as feedback gain of e xy . We assume that V 2,d can be generated by a desired pseudo-control input u 2xy = u 2x u 2y T ∈ ℝ 2 , so that for ̇ Ib xy = u 2xy , the desired time derivative V 2,d and the time derivative of the actual LF V 2 become equal which yields the pseudo-control law for the ART: Note that the pseudo-control law for the ART includes ̇u 1x and ̇u 1y , the time derivatives of the TVT's desired pseudocontrols. To obtain the dynamics of u 1xy with respect to the b-frame, we use Eq. (69) and apply Euler's classical treatment of vector analysis [39] which yields In Eq. (82), the time derivative of the inverse thrust-magnitude matrix with respect to the b-frame can be calculated to [26] (79) and the time derivative of the rotation matrix R bI with respect to the b-frame is given by the following kinematics [26]: being related to the Strapdown equation [39]. To finally constitute the pseudo-control law for the overall ART in the vector space ℝ 3 we use Eq. (81) and, additionally, assume that u 2z ∈ ℝ , which represents the pseudo-control law of the one-step backstepping yaw rate controller, is already given from Sect.

leading to
Since all vector elements of u 2 ∈ ℝ 3 can be interpreted as a desired angular body acceleration around the x b -,y b -, and z b -axis of the CDM, the relation between the desired moments M R d and u 2 is given by inserting Eq. (85) into the CDM's rotational dynamics (see Eq. (53)): Since the CDM's output vector v d , given by Eq. (48), also represents, conversely, the input to stabilize the inner loop of the CDM, Eq. (86) emphasizes that the angular body accelerations match their purposed desired values, being comprised in u 2 , if v d is mapped into the attainable set of the actuators' squared commanded rotational speeds, established through Eq. (51). This mapping also involves the total desired force T R d , or rather T R , which can be determined by propagating the desired pseudo-control u 1z to a first-order integrator chain: With regard to Eq. (70), the initial condition of Eq. (87) is predefined by T R t 0 = 0.1 N . The final result for the baseline flight controller's desired virtual control vector can then be summarized to

Yaw rate controller
Using the thrust vector t R as reduced attitude state allows an independent yaw control design for the CDM, since the angular body rate Ib z ∈ ℝ is entirely decoupled from the attitude kinematics, given through Eqs. (55) and (56). The control objective of the yaw rate controller is achieved when Ib z tracks a desired reference trajectory Ib z,r ∈ ℝ, so that the yaw motion tracking error e z ∈ ℝ and its time derivative converge to zero: To derive the control law for the yaw rate tacking, we choose again a positive definite LF candidate V 3 [53] in quadratic form as whereby the time derivative of V 3 ∈ ℝ is calculated to To enforce an exponential decay behavior of ̇e z , Eq. (92) must ensure negative definiteness. Therefore, let ̇e z,d = −k ⋅ e z denote a desired decay behavior of the real yaw rate error dynamics. Then, we use Eq. (92) once more and write a desired time derivative of V 3 as where the feedback gain k z ∈ ℝ is positive definite and constant. We assume that V 3,d is generated by a desired pseudocontrol input u 3 ∈ ℝ such that for ̇I b z = u 3 , Eqs. (92) and (93) become equal which yields the pseudo-control law for the yaw rate tracking after inserting Eq. (90) into Eq. (92). We conclude from Eq. (94) that yaw rate tracking is guaranteed under the assumption that the angular body acceleration is adjusted in a manner that ̇I b z equals u 3 at any time, so that ̇e z declines to zero exponentially fast. For the sake of consistency, we write u 2z = u 3 .

Reference models
Within this section, reference models are presented which are used to generate feasible trajectories for the baseline controllers that shall ultimately be tracked by the states of the plant. Feasible means that the reference trajectories shall be sufficient smooth, so that they can physically be achieved by the Mars vehicle. Otherwise, the baseline controllers' tracking objectives would be doomed to fail in the long term [54]. Usually, linear reference models are used to obtain the r -trajectories since they provide the simplest approach to plan a desired trajectory by only using a n-dimensional filter. To obtain the desired trajectory for the yaw rate controller, the pilot's yaw rate command Ib z,c is filtered by a PT1 element. This PT1 filter represents a linear first-order reference model and, based on [47], is given by whereby the time constant t z is constantly chosen as 1.25 s . This emphasizes that the attainable vector space for a desired moment around the body z b -axis is reduced to a minimum (while the yaw rate controller is active) being beneficial with regard to the low-density conditions in Valles Marineris (see Sect. 2.1). In case of the baseline position controller, the linear second-order reference model [47] 1

Stability analysis of closed-loop system
To provide a valid baseline flight controller, closed-loop stability of the CDM must be guaranteed which is shown in this section. To track the reference position r R r as well as the reference velocity v R r , both given by Eq. (96), we define the total error of the CDM's position loop e pos ∈ ℝ 6 using Eqs. (57) and (58) to As we are aiming to stabilize the dynamics of Eq. (109) around its zero equilibrium, we use the I-frame as notation and reference frame while applying Euler's classical treatment of vector analysis [39] which yields where ë p is given by Eq. (60). To obtain a matrix-vector notation, the position loop error dynamics from Eq. (110) are reformulated into where the matrices M ∈ ℝ 6x6 , A ∈ ℝ 6x6 , and B ∈ ℝ 6x3 are given by including I ∈ ℝ 3×3 as identity matrix and 0 ∈ ℝ 3×3 as zero matrix. Since it is assumed that every demanded input for the position loop is, based on the theory of time-scale separation [53], always available and disposable, t R equals t R c (106) For the desired position error dynamics ë p,d , given by Eq. (61), we only consider the feedback matrices K x and K v in the closed-loop system. By multiplying Eq. (113) with M −1 , it follows: with denoting the error matrix of the position loop. All inherent nonlinearities of the CDM's outer loop are canceled out by the pseudo-control law (62) which yields a linear statespace representation for ̇e pos . To achieve asymptotic stability, meaning that the linear error dynamics (114) show an exponential decay behavior, so that e pos converges at its zero equilibrium 0 ∈ ℝ 6 , E must be a stability matrix, and thus, Eq. (115) is only allowed to exhibit poles in the left-half complex plane. To examine E , we use the so-called Lyapunov equation [53] which manifests that a positive definite LF V pos ∈ ℝ must be found satisfying V pos e pos > 0 and V pos (0) = 0 for all e pos ≠ 0 . We obtain a valid LF candidate in quadratic form as whereby the time derivative of V pos is calculated to Equations (116) and (117) both include the matrix P ∈ ℝ 6×6 which we can chose arbitrarily. The only primary constraint is that P must be positive definite as well as symmetric, so that P = P T is always satisfied [53]. After inserting the linear error dynamics, given by Eq. (114), into Eq. (117), we find so that the matrix Q ∈ ℝ 6×6 is given by the Lyapunov equation [53]: While choosing, e.g., P as identity matrix I ∈ ℝ 6×6 , we conclude from Eq. (119) that Q becomes positive definite and therefore that E is Hurwitz [53]. Thus, Eq. (114) is a stable linear system while choosing the feedback gain matrices K x and K v positive definite which guarantees position and velocity tracking at any given time. Since the dynamics of the CDM's inner loop are considerable faster than the dynamics of the outer loop, we conclude that the closed-loop system is stable under the pseudo-control inputs u 1 and u 2 , given by Eqs. (69) and (85), provided that the feedback gain matrices K t , K xy , and k z are positive definite and constant. Note that the desired pseudo-control law for the TVT always exists, since we saturate T R to a minimum permissible propulsive thrust of 0.1 N to avoid that the inverse thrust-magnitude matrix T −1 , given through Eq. (70), becomes singular [53].

Artificial augmentation units
In this section, artificial augmentation units are developed to extend the baseline flight controller's architecture, shown in Fig. 6, to increase its capabilities in suppressing modeling errors artificially. This procedure is advantageous to achieve objective (ii) and (iii), since the units work independent from the baseline controller. Thus, the control laws derived in Sect. 3.1 are still valid.
Due to the uncontrollability of the balloon UAV, unintended rope oscillations are causing steady-state errors [47] in the multirotor system's position. Additionally, when the rope touches the rotor plane, or rather the propeller tips, the flight mission is endangered and may fail. The idea is therefore to damp the balloon motion artificially within the horizontal plane by generating a closed-loop reference model (CRM) [35] either for the inner or the outer loop of the CDM. Based on the developed methodology, we call our first approach the thrust vector augmentation (TVA), whereas our second approach is called the position command augmentation (PCA).
The general idea of both TVA and PCA is to embed a feedback of the uncertain plant into the reference model so that a closed-loop results. If this feedback is defined as an error, it tries to pull the reference model towards the plant so that both meet "half-way" and the error is, at least, reduced [35]. We use the balloon's relative position to the multicopter r RQ xy = x RQ y RQ T ∈ ℝ 2 , given by Eq. (13), and the balloon's relative velocity v RQ xy = u RQ v RQ T ∈ ℝ 2 , given by the vector subtraction of v R xy = u R v R T ∈ ℝ 2 from v Q xy = u Q v Q T ∈ ℝ 2 , as error signals and embed them into the baseline controllers' open-loop reference models (ORM). These error signals are then scaled by the Hurwitz matrices D Bal and K Bal , the so-called Luenberger gains [35], which are both of the set ℝ 2×2 .

Thrust vector augmentation
We implement the TVA approach by augmenting the commanded thrust vector t R c which represents the pseudo-control input of the CDM's outer loop. The augmented thrust vector is then defined as whereby t R c is given by Eq. (62) and the weighted matrix W ∈ ℝ 3×2 is given through Eq. (79). Due to the TVA, Eq. (98) is redefined, so that the unit vector z b,c indicates the pointing direction of t R c without violating that ‖z b,c ‖ = 1 at every time instant. Hence, while the Mars vehicle operates on a mission, a balance between a desired attitude command and a reduction of the balloon error signals is obtained. Since the error signals are only defined within the horizontal plane and, by scaling K Bal and D Bal to a lower level, the inner loop's control objective remains and, finally, concludes in the stabilization of the rope angle R * around 90 • for steadystate conditions (see Fig. 4). The main advantage of the TVA is based on its robustness against external disturbances. While the TVA is activated, the baseline control target-to stabilize the tracking error e p around its equilibrium 0 ∈ ℝ 3 -cannot be bypassed. Only a significant growth of e p is obtained for non-steady-state conditions since a mismatch between the desired pseudo-control input t R c and the actual augmented pseudo-control input t R c of the CDM's outer loop is produced on purpose to stabilize the relative position of the balloon in the forward flight case.

Position command augmentation
To realize the PCA approach and generate a CRM inside the position loop, we augment the initial position command r R c , so that the new pilot command is defined as whereby the weighted matrix W ∈ ℝ 3×2 is given through Eq. (79). In case of the PCA, the position tracking error in Eq. (57) is redefined to ẽ p =r R r − r R , since the r-trajectories of the second-order reference model, given by Eq. (96), are fed by the augmented position command r R c as feedforward. Hence, while the Mars vehicle operates on a mission, a balance between a desired position command and a reduction of the balloon error signals is obtained. Since the error signals are only defined within the horizontal plane and, by scaling K Bal and D Bal to a lower level, the outer loop's control objective remains and, finally, concludes in the stabilization of the rope angle R * around 90 • for steady-state conditions (see Fig. 4). Although a disproportionate growth , of ẽ p is prevented, by pulling the reference model towards the plant, the overall controller target, to reach a steadystate in the vehicle's position, so that e p converges to zero, is only satisfied in the absence of stationary external disturbances. In case of, e.g., static wind, the error feedback ‖r RQ xy ‖ is always positive and low-frequent. Thus, reference model and plant do constantly meet "half-way", so that the overall control target is misled by the deviation between r R r and r R r . To overcome this problem, the linear error term K Bal r RQ xy is passed through a washout filter [39] satisfying that the outer loop's control objective remains even under steady wind by only using the transient rate of the feedback. We embed the washout filter as second-order system and only use the first derivative for Eq. (121) which will solely produce non-zero augmentation signals when the linear error term is highfrequent and not steady. Figure 7 shows the schematic structure of the PCA approach including the washout filter where 0 = 0.6 1 s , = 1 and k s = 1 define its system parameters.

Simulation results
This section summarizes the simulation results of testing the CDM of Sect. 3 on the high-fidelity plant design model of Sect. 2. Since the overall simulation testbed is treated as MIL, the controller states are initialized through the plant state vector x P so that, inter alia, the CDM's b-frame and the plant's B-frame are congruent from scratch. To demonstrate the controller performance, the closedloop system response for the Mars vehicle's position state r R is evaluated on a digital-8-flight-maneuver at the bottom of Valles Marineris. The flight path is obtained by a reference trajectory being in shape of a long curved "8" representing a suitable extraterrestrial mission profile to, e.g., scan the Martian surface from above and detect prominent spots where swarm participants should perform scientific experiments.
The mission has a total flight time of 180 s . The actual maneuver starts eastwards and is performed within 120 s . This provides a time slot of 60 s to not only reach a steadystate in the vehicle position, but also to eliminate steady-state errors [47]. The take-off location is initialized according to the environment submodel of Sect. 2.1 with a geodetic altitude of -4906 m blow the Martian geoid. Since R * describes the radial distance between the rope and the x B -y B plane of the multirotor system, R * ,ne = 20 • is defined as never exceed rope angle for a quantitative mission assessment. To make the handling more intuitive and, to provide proper smooth trajectories, the pilot's position command arises from an integrated joystick-velocity command v R cmd . For the sake of comparability, v R cmd is constantly given by ±3 ±3 ±1 T , so that the pilot's position command is determined as follows: Equation (122) also ensures that, for every time instant, the multirotor system does not accelerate upwards faster than the balloon. For the baseline controllers, the feedback gain matrices K v = diag 2 2 1 , K x = diag 0.5 0.5 0.5 ,  tracks the reference position r R r (red) even under influence of the balloon UAV and rope dynamics very accurately which implies a satisfying robustness for all three controller configurations. In terms of the TVA, the oscillating offset between r R and r Q is reduced, so that an additional offset between r R and r R r is caused on purpose, due to the mismatch between t R c and t R c , to reduce the rope oscillations. A significant reduction of this additional offset is achieved in terms of the PCA by smoothly shifting the aimed position tracking objective from the ORM to the CRM for non-steady-state conditions. Due to the fact that, while the PCA is activated, the augmented reference position r R r does distinguish from the nominal reference position r R r , the Mars vehicle's position tracking for the flight controller with activated PCA shows the nominal trajectory for the reference position r R r , being depicted as black dashed line, and additionally, shows the actual augmented trajectory for the reference position r R r , which is tracked, in red. For the sake of clarity, the position tracking for the flight controller with activated TVA also shows the nominal trajectory for the reference position r R r (122)  Figure 9 depicts the vector elements of the position tracking error e x,p , e y,p , and e z,p ( ẽ x,p , ẽ y,p and ẽ z,p for the PCA) in x I -, y I -, and z I -direction of the I-frame over time. In case of the TVA, the position errors are mainly capped between ∓0.5 m for e x,p and even ∓2 m for e y,p , whereas for the baseline controller and for the baseline controller with PCA e x,p and e y,p , respectively, ẽ x,p and ẽ y,p , are settled below ∓0.2 m and ∓0.6 m . Due to the artificial augmentation units, the Mars vehicle is able to not only reach a steady-state in position, but also to eliminate all steady-state errors, so that ‖e p ‖ , or rather ‖ẽ p ‖ , converge to zero at the end of the simulation period. The baseline controller itself is not able to eliminate the steady-state errors which is caused by the modeling errors inside the CDM (neglection of balloon and rope dynamics).
For the sake of completeness, Fig. 10 illustrates the attitude tracking (TVT) in x I -, y I -, and z I -direction over time for the baseline controller during the digital-8-maneuver. The thrust vector elements t R x , t R y , and t R z track the reference thrust vector elements t R x,r , t R y,r , and t R z,r very accurately which implies a good controller performance similar to the position tracking simulation results presented before. This is expected since t R x,c , t R y,c , and t R z,c , the vector elements of the desired pseudo-control input for the CDM's outer loop, are only used as virtual control to stabilize the position loop exponentially. The initial control objective to track the ORM remains without stabilizing the balloon dynamics in In Fig. 11, the attitude tracking (TVT) in x I -, y I -, and z I -direction over time for the baseline controller with PCA is presented during the digital-8-maneuver. The controller performance of tracking the reference thrust vector elements t R x,r , t R y,r , and t R z,r by the vehicle's estimated thrust vector elements t R x , t R y , and t R z are comparably good as the simulation results of the TVT for solely using the baseline controller, without PCA, as shown in Fig. 10. The only difference can be recognized in the end of the simulation period. Here, all disturbing rope oscillations, illustrated by a swirling of t R x , t R y , and t R z , since the rope pulls the rotor plane back and forth, are canceled out entirely being established through t R x,c , t R y,c , and t R z,c . Based on the outer loop's CRM, the desired pseudo-control t R c propagates the augmented trajectories r R r and ṽ R r for the position and velocity reference to the nonlinear reference model of the CDM's inner loop which concludes in a "PCA-dominated" trajectory design for the reference thrust vector t R r . In Fig. 12, the attitude tracking (TVT) in x I -, y I -, and z I -direction of the I-frame over time for the baseline controller with TVA is presented during the digital-8-maneuver. The elements of the estimated thrust vector t R x , t R y , and t R are tracking the reference thrust vector elements t R x,r , t R y,r , and t R z,r with high accuracy which implies a good controller performance. Similar to the simulation results of the TVT while using the baseline controller with PCA, shown in Fig. 11, the rope oscillations do completely vanish in the end of the flight maneuver. However, for the TVT, the augmented thrust vector command t R c does distinguish from the nominal thrust vector command t R c . That is why, the attitude tracking in Fig. 12 shows the nominal trajectories for the thrust vector command t R x,c , t R y,c , and t R z,c , and, additionally, shows the actual augmented trajectories for the thrust vector command t R x,c , t R y,c , and t R z,c to emphasize that the TVT objective is smoothly shifted from the ORM to the CRM for non-steady-state conditions. The elements t R x,c and t R y,c do, besides their initial task to control r R xy and v R xy by stabilizing second-order position error dynamics, also stabilize the balloon dynamics around their zero equilibria, so that the rope oscillations are canceled out entirely at the end of the simulation period. In addition to that, the TVT in z I -direction of the I-frame over time shows that t R z,c and t R z,c are equal, since the augmentation unit solely damps the balloon motion artificially within the horizontal plane. The latter two statements lead to the fact that the position error e x,p as well as e y,p in x I -and y I -direction of the I-frame must increase on average which can be confirmed by Fig. 9.
To provide a quantitative mission assessment, Fig. 13 shows the simulation results for the rope angle variation over time during the digital-8-maneuver for all three controller configurations. For the baseline controller, the rope angle R * is capped between 90 • and 30 • . For the baseline controller with PCA, R * is capped between 90 • and 40 • which is further reduced to a range between 90 • and 50 • in terms of the baseline controller with TVA. The flight mission can be stated as successful and practicable for all of the three controller configurations since R * ,ne is not undershot. In terms of a smallest possible rope angle range, the baseline controller with TVA shows the best performance, or rather, is most successful in damping the balloon motion artificially within the horizontal plane. Nevertheless, the rope oscillations can not only be reduced significantly, but also canceled out entirely at the end of the flight maneuver when the baseline controller is extended artificially by either the PCA or the TVA approach. At this time, the balloon is perfectly arranged above the multirotor system which cannot be achieved by solely using the baseline controller without the artificial augmentation units.

Summary, conclusion, and outlook
In this paper, we presented a nonlinear flight controller for a multibody Mars vehicle guaranteeing asymptotic tracking position control with active oscillation damping. Relating thereto, as main contribution, we have studied two artificial augmentation approaches for the baseline controller: • Thrust vector augmentation (TVA), • Position command augmentation (PCA), which both increase the controller capabilities to suppress modeling errors artificially without changing the baseline control laws. Basically, by utilizing oscillation damping feedbacks of the uncertain plant which are applied as additional commands to either the inner or the outer loop of the controller's reference model. We pose that both TVA and PCA approach are advantageous to cope with plant uncertainties and that they can be added to an existing controller architecture (linear/nonlinear).
By utilizing the classical Newtonian formulations, a nine degrees of freedom multibody UAV model was derived including a six degrees of freedom rigid-body model for a hexrotor system and a three degrees of freedom point-mass model for a balloon UAV to setup a first stage high-fidelity simulation model. Kinematic constraints were imposed on both UAV models as massless spring-damper oscillator model to approximate a realistic rope connection. To simulate the extraterrestrial flight mission under realistic conditions, an atmospheric model was built on high-resolution measurement data being smoothed and prepared by a meteorological database of Mars to best fit 14.035° S and 58.5° W, the vehicle's geographic takeoff location at Valles Marineris. A nonlinear control design model of the hexrotor system was created and considered as baseline controller providing a cascaded structure. Feedback linearization was applied for the outer loop and analytical backstepping for the inner loop to derive all feedback laws. To cope with uncertain dynamic effects of the plant, being represented by oscillating motion of the suspension line, two artificial augmentation units were developed (TVA and PCA approach) both extending the baseline flight controller's architecture separately. Simulation results were presented for an eight-shaped flight maneuver at the bottom of Valles Marineris proving that both augmentation units actively damp the balloon motion in the forward flight case for non-steady-state conditions to counteract the rope oscillations and finally stabilize the rope angle around its equilibrium while keeping the system within its limits. The Mars vehicle was only able to reach a steady state in position including either the TVA or PCA approach in the closed-loop control system underlining their importance for the extraterrestrial mission simulation. However, further analysis of both proposed approaches is still required, since stability or rather semi-global stability (concerning closed-loop reference models for uncertain plants whose states are accessible) was not in the focus of this paper.
Future work should therefore hone the artificial augmentation units in a second stage high-fidelity simulation model including sensor models and state estimators for both hexrotor system and balloon UAV. Additionally, proper trajectory planning should be examined to ensure that (a) smoother rope angle trajectories are provided by adding additional feedforward terms to the baseline controller's architecture and (b) the closed-loop control system stays within predefined safety bounds. To investigate the performance of the controller when subjected to exogenous disturbances, different wind scenarios shall be applied to the controlled plant to reproduce mission-critical environmental conditions such as dust storms on Mars, and especially, in the Valles Marineris.
Acknowledgements The results shown in this paper have been mostly developed from July 2017 until Dec. 2018 during the mid-term and master thesis of the first author at the Institute of Flight System Dynamics, Technical University Munich. The idea to only use the transient rate of the relative position error signal for the closed-loop reference model within the control design model's outer loop has been developed from the first author during his research at the Institute of Air Transportation Systems, Technical University Hamburg, where he continued his study in control theory.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.