## 4 Autonomous Landing Method
### Trajectory Generation and Tracking
#### 4.3.2 Landing Trajectory Generation  
A suitable trajectory must be planned to enable the FPLIR to reach the designated target location in 3D space. Two problems should be considered for the FPLIR: (1) the input vector \(\mathbf{u}\) from the system is limited, because each motor can only produce limited thrust; and (2) the FPLIR has an underdriven characteristic due to the motor's thrust along the \(z_{b}\) axis of the body. To achieve motion in 3D space, it is necessary to combine translational and rotational motion using system dynamics. Therefore, the generated trajectories must be dynamically performable for the FPLIR.  
The generation of landing maneuvers is modeled as a discrete-time multiple-shot optimization problem for \(N\) shooting points over a non-fixed time horizon \(T\), which can be expressed as follows:  
\[\underset{u_{0}\cdots u_{N-1}}{min}T\sum_{k=0}^{N}\left\|\widetilde{y}\right\| ^{2}Q_{k} \tag{19a}\] \[\text{subject to : }\mathbf{x}_{0}=\mathbf{x}_{init}\] (19b) \[T_{min}\leq T\leq T_{max}\] (19c) \[\mathbf{x}_{k+1}=f(\mathbf{x}_{k},u_{k})\forall k\in[0,N-1]\] (19d) \[z_{min}\leq p_{d}\forall k\in[0,N]\] (19e) \[0\leq\gamma\leq\gamma_{max}\forall k\in[0,N]\] (19f) \[\mathbf{u}_{min}\leq\mathbf{u}\leq\mathbf{u}_{max}\forall k\in[0,N-1] \tag{19g}\]  
The optimization problem is constructed as follows: Equation (19a) is the cost function to minimize, including final and running terms; Constraint (19c) presents the limits of the total maneuver time \(T\); Equation (19d) is the dynamics of the system; Constraint (19e) is the minimum allowable height; Constraints (19f) and (19g) are the constraints on the motor thrusts and their derivatives.  
The cost function in Equation (19a) consists of a set of errors \(\widetilde{y}\) that are dependent on the states and the inputs of the system, weighted by a diagonal matrix \(Q_{k}\) for each shooting point. Different values of \(\widetilde{y}\) are used to represent the terminal and running costs.  
\[\widetilde{\mathbf{y}}_{k}=\begin{cases}\left[\left[\gamma_{k}\frac{T}{N}+\mathbf{u} _{k}\frac{T^{2}}{2N^{2}}\right]^{T}\mathbf{\omega}_{k}^{T}\right]^{T}&k\in[0,N-1]\\ \left[\begin{array}{cc}\widetilde{\mathbf{p}}_{k}^{T}&\widetilde{\mathbf{q}}_{k}^{T }&\widetilde{\mathbf{\omega}}_{k}^{T}&\widetilde{\mathbf{\omega}}_{k}^{T}\\ \end{array}\right]^{T}&k=N\end{cases} \tag{20}\]  
The running cost minimizes the integral of the motor thrusts, as well as the angular velocities \(\mathbf{\Omega}_{k}^{T}\) of the robot. The terminal cost minimizes the position and orientation error \(\widetilde{\mathbf{p}}_{k}^{T},\widetilde{\mathbf{q}}_{k}^{T}\), as well as the final linear and angular velocity errors \(\widetilde{\mathbf{v}}_{k}^{T},\widetilde{\mathbf{\Omega}}_{k}^{T}\) at the desired landing position.  
In our study, a variable time horizon is implemented by modeling the system dynamics using a Runge-Kutta4 integration of the state space, scaling its derivative by the total time \(T\) and using an integration step of \(1/N\) seconds. Since problems where the total maneuvering time is an optimization variable suffer from poor linearization characteristics, the CasADi framework was chosen as the solver to embed a linear system solver with high numerical stability.