Robust Translational Force Control of Multi-Rotor UAV for Precise Acceleration Tracking

TL;DR

This paper presents a robust translational force control method for multi-rotor UAVs that accurately tracks desired accelerations by considering multi-rotor dynamics and using disturbance observers for disturbance rejection.

Contribution

The paper introduces a novel force control approach that incorporates multi-rotor dynamics into the kinematic inversion process and employs DOB for disturbance robustness.

Findings

Enhanced acceleration tracking accuracy demonstrated in simulations.

Robustness against external disturbances validated through experiments.

Stable control performance confirmed via mu-analysis.

Abstract

In this paper, we introduce a translational force control method with disturbance observer (DOB)-based force disturbance cancellation for precise three-dimensional acceleration control of a multi-rotor UAV. The acceleration control of the multi-rotor requires conversion of the desired acceleration signal to the desired roll, pitch, and total thrust. But because the attitude dynamics and the thrust dynamics are different, simple kinematic signal conversion without consideration of those difference can cause serious performance degradation in acceleration tracking. Unlike most existing translational force control techniques that are based on such simple inversion, our new method allows controlling the acceleration of the multi-rotor more precisely by considering the dynamics of the multi-rotor during the kinematic inversion. By combining the DOB with the translational force system that includes the improved conversion technique, we achieve robustness with respect to the external force disturbances that hinders the accurate acceleration control. mu-analysis is performed to ensure the robust stability of the overall closed-loop system, considering the combined effect of various possible model uncertainties. Both simulation and experiment are conducted to validate the proposed technique, which confirms the satisfactory performance to track the desired acceleration of the multi-rotor.

Tables

Table 1. TABLE I: PHYSICAL QUANTITIES AND CONTROLLER GAINS.

Name	Value	Name	Value
$P_{\varphi ,\theta }$	3	Mass	3.24 $\mathrm{Kg}$
$D_{\varphi ,\theta }$	1	${\overline{J}}_{1,2}$	0.82 $\mathrm{Kg}\cdot m^2$
$\ddot{𝐗}$ Limit	$\mathrm{pm3}m/s^2$	${\overline{J}}_3$	1.49 $\mathrm{Kg}\cdot m^2$
$\max \left|{\delta }_j\right|$	$0.12$	$\max \left|J_{\Delta ,j}\right|$	$0.3$
$\max \left|K_{\Delta ,j}\right|$	$0.1$	$\zeta$	0.707

Equations

$$\left\{\begin{array}[]{lr}m\ddot{\mathbf{X}}=\mathbf{F}+m\mathbf{g}=R(\mathbf{q})\mathbf{T}_{t}+m\mathbf{g}\\ J\dot{\mathbf{\Omega}}=\mathbf{T}_{r}-\mathbf{\Omega}\times J\mathbf{\Omega}\;,\end{array}\right.$$

$$J\ddot{\mathbf{q}}={\mathbf{T}}_r$$

$$m\ddot{\mathbf{X}}=-R(\psi )\left[\begin{array}{c}\cos \varphi \sin \theta \\ -\sin \varphi \\ \cos \varphi \cos \theta \end{array}\right]T_t+m\mathbf{g},$$

$$\tilde{\ddot{\mathbf{X}}}=\begin{bmatrix}\tilde{\ddot{x}}\\ \tilde{\ddot{y}}\\ \tilde{\ddot{z}}\end{bmatrix}=R^{-1}(\psi)\left(\ddot{\mathbf{X}}-\mathbf{g}\right)=R^{-1}(\psi)\big{(}{{1}\over{m}}\mathbf{F}\big{)}.$$

$$m\widetilde{\ddot{\mathbf{X}}}=-h(\varphi ,\theta )T_t=-\left[\begin{array}{c}\cos \varphi \sin \theta \\ -\sin \varphi \\ \cos \varphi \cos \theta \end{array}\right]T_t.$$

$$\mathbf{r}=\left[\begin{array}{c}\theta \\ \varphi \\ T_t\end{array}\right]=\left[\begin{array}{c}\arctan \left(\frac{\tilde{\ddot{x}}}{\tilde{\ddot{z}}}\right)\\ \arctan \left(-\frac{\tilde{\ddot{y}}\cos \theta }{\tilde{\ddot{z}}}\right)\\ -\frac{m\tilde{\ddot{z}}}{\cos \varphi \cos \theta }\end{array}\right].$$

$${\theta }_d$$

$${\varphi }_d$$

$$T_{t,d}$$

$$\mathbf{r}(t)=\left[\begin{array}{c}\theta (t)\\ \varphi (t)\\ T_t(t)\end{array}\right]\approx \left[\begin{array}{c}{\theta }_d\left(t-{\gamma }_{\theta }\right)\\ {\varphi }_d\left(t-{\gamma }_{\varphi }\right)\\ T_{t,d}(t)\end{array}\right].$$

$$m\tilde{\ddot{\mathbf{X}}}=h\Big{(}\phi_{d}\big{(}t-\gamma_{\phi}\big{)},\theta_{d}\big{(}t-\gamma_{\theta}\big{)}\Big{)}T_{t,d}(t).\\$$

$$\left[\begin{array}{c}\tilde{\ddot{x}}(t)\\ \tilde{\ddot{y}}(t)\\ \tilde{\ddot{z}}(t)\end{array}\right]=\left[\begin{array}{c}\left(\frac{\cos {\varphi }_d(t-{\gamma }_{\varphi })}{\cos {\varphi }_d(t)\cos {\theta }_d(t)}{\tilde{\ddot{z}}}_d(t)\right)\sin {\theta }_d(t-{\gamma }_{\theta })\\ -\left(\frac{1}{\cos {\varphi }_d(t)\cos {\theta }_d(t)}{\tilde{\ddot{z}}}_d(t)\right)\sin {\varphi }_d(t-{\gamma }_{\varphi })\\ \left(\frac{\cos {\varphi }_d(t-{\gamma }_{\varphi })\cos {\theta }_d(t-{\gamma }_{\theta })}{\cos {\varphi }_d(t)\cos {\theta }_d(t)}\right){\tilde{\ddot{z}}}_d(t)\end{array}\right]$$

$$T_{t,d}=-\frac{m{\tilde{\ddot{z}}}_d(t-{\gamma }_v)}{\cos {\varphi }_d(t-{\gamma }_{\varphi })\cos {\theta }_d(t-{\gamma }_{\theta })}.$$

$$\left[\begin{array}{c}\tilde{\ddot{x}}(t)\\ \tilde{\ddot{y}}(t)\\ \tilde{\ddot{z}}(t)\end{array}\right]=\left[\begin{array}{c}\tan {\theta }_d(t-{\gamma }_{\theta })\\ -\frac{\tan {\varphi }_d(t-{\gamma }_{\varphi })}{\cos {\theta }_d(t-{\gamma }_{\theta })}\\ 1\end{array}\right]{\tilde{\ddot{z}}}_d(t-{\gamma }_v).$$

$$\left\{\begin{array}[]{lr}\theta(t)=\arctan{\left({{\tilde{\ddot{x}}_{d}(t-\gamma_{\theta})}\over{\tilde{\ddot{z}}_{d}(t-\gamma_{\theta})}}\right)}\\ \phi(t)=\arctan{\left(-{{\tilde{\ddot{y}}_{d}(t-\gamma_{\phi})\cos{\theta_{d}(t-\gamma_{\phi})}}\over{\tilde{\ddot{z}}_{d}(t-\gamma_{\phi})}}\right)}.\end{array}\right.$$

$$\begin{bmatrix}\tilde{\ddot{x}}(t)\\ \tilde{\ddot{y}}(t)\\ \tilde{\ddot{z}}(t)\end{bmatrix}=\begin{bmatrix}\tilde{\ddot{x}}_{d}(t-\gamma_{h})\big{(}{{{\tilde{\ddot{z}}_{d}(t-\gamma_{v})}\over{\ddot{z}}_{d}(t-\gamma_{h})}}\big{)}\\ \tilde{\ddot{y}}_{d}(t-\gamma_{h})\big{(}{{{\tilde{\ddot{z}}_{d}(t-\gamma_{v})}\over{\ddot{z}}_{d}(t-\gamma_{h})}}\big{)}\\ \tilde{\ddot{z}}_{d}(t-\gamma_{v})\end{bmatrix}.$$

$$T_{t,d}=-\frac{m{\tilde{\ddot{z}}}_d(t)}{\cos \varphi (t)\cos \theta (t)},$$

$$\widetilde{\ddot{\mathbf{X}}}=\left[\begin{array}{c}{\tilde{\ddot{x}}}_d(t-{\gamma }_h)\left(\frac{{\tilde{\ddot{z}}}_d(t)}{{\tilde{\ddot{z}}}_d(t-{\gamma }_h)}\right)\\ {\tilde{\ddot{y}}}_d(t-{\gamma }_h)\left(\frac{{\tilde{\ddot{z}}}_d(t)}{{\tilde{\ddot{z}}}_d(t-{\gamma }_h)}\right)\\ {\tilde{\ddot{z}}}_d(t)\end{array}\right]\approx \left[\begin{array}{c}{\tilde{\ddot{x}}}_d(t-{\gamma }_h)\\ {\tilde{\ddot{y}}}_d(t-{\gamma }_h)\\ {\tilde{\ddot{z}}}_d(t)\end{array}\right],$$

$$\mathbf{F}=m(\ddot{\mathbf{X}}-\mathbf{g}).$$

$$\hat{\kappa }(s)={\widehat{\tilde{\mathbf{F}}}}_d(s)+{\hat{\mathbf{d}}}_{EID}(s)=P_n^{-1}(s)\tilde{\mathbf{F}}(s).$$

$${\hat{\mathbf{d}}}_{EID}=Q_1(s)\hat{\kappa }(s)-Q_2(s){\tilde{\mathbf{F}}}_d(s).$$

$${\tilde{\mathbf{F}}}_d={\mathbf{F}}_d-{\hat{\mathbf{d}}}_{EID}.$$

$$\kappa ={\tilde{\mathbf{F}}}_d+{\mathbf{d}}_{EID}={\mathbf{F}}_d-{\hat{\mathbf{d}}}_{EID}+{\mathbf{d}}_{EID}\approx {\mathbf{F}}_d.$$

$$\frac{{\mathbf{q}}_i(s)}{{\mathbf{T}}_{r,i}(s)}=\frac{1}{J_i{s}^2},$$

$$\frac{{\mathbf{T}}_{r,i}(s)}{{\mathbf{q}}_{i,d}(s)-{\mathbf{q}}_i(s)}=P_i+D_{i}s$$

$${\Lambda }_{n,a,i}(s)=\frac{{\mathbf{q}}_i(s)}{{\mathbf{q}}_{i,d}(s)}=\frac{D_{i}s+P_i}{J_i{s}^2+D_{i}s+P_i}.$$

$${\Lambda }_{n,t}(s)=1.$$

$${\Lambda }_n(s)=\left[\begin{array}{ccc}{\Lambda }_{n,a,2}(s)& 0& 0\\ 0& {\Lambda }_{n,a,1}(s)& 0\\ 0& 0& {\Lambda }_{n,t}(s)\end{array}\right].$$

$$P_n(s)={\Lambda }_n(s).$$

$${Q_{1}}(s)=diag\big{(}{Q_{1,h}}(s),{Q_{1,h}}(s),{Q_{1,v}}(s)\big{)},$$

$$Q_{1,h}(s)=\frac{1}{({\tau }_{1}s{)}^2+\zeta ({\tau }_{1}s)+1},$$

Full text

Note to Practitioners:

Precise acceleration tracking is a fundamental requirement of multi-rotor unmanned aerial vehicles (UAVs) for widening their application area beyond basic autonomous flight. For such an objective, we need an accurate three-dimensional force control and a robust rejection method of translational motion disturbance.

First, for accurate force control, the target force command must be converted to the appropriate target attitude and thrust value, because the multi-rotor generates the three-dimensional translational forces by the combination of the current attitude and the total propeller thrust [1]. Once the target attitude and the total thrust command are determined, each value passes through attitude and thrust dynamics that are quite different from each other : the process of achieving actual attitude involves feedback attitude control[2], torque generation by the combination of motor’s thrust, followed by the rotation of the fuselage that has larger moment of inertia than the propellers. Due to such difference, simple kinematic conversion of the force signal without consideration of the actual attitude and thrust dynamics can cause unsynchronized realization of the attitude and the total thrust, which degrades acceleration tracking performance given that the acceleration of the multirotor is determined by the combination of attitude and thrust. To the best of our knowledge, however, many studies have not investigated this issue. In [1] and [3], the target thrust signal was computed without considering attitudinal dynamics while treating Z-directional translational dynamics as a separate channel to other horizontal dynamics. All three axes of translational dynamics have been simultaneously considered in [4] during the conversion process, but they also did not reflect the different characteristics of attitude and thrust dynamics. Those differences become noticeable in multi-rotors that have large moment of inertia, due to significant time delay between input and output attitude.

Second, for a satisfactory level of translational disturbance rejection, we need a controller that estimates and offsets the effect of the disturbance[5, 6]. As a way to achieve this goal, we can consider constructing the Disturbance Observer (DOB)-based robust control algorithm [7]. However, although several studies applied the Disturbance Observer (DOB) robust control technique to their controllers [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], most of them [8, 9, 10, 11, 12, 13, 14, 15] aimed to maintain the nominal attitude control performance against torque disturbance. Therefore, this approach has a limitation in overcoming translational movement disturbances. Only a few studies exist on applying DOB to overcome the translational force disturbances including [16], [17], and [18]. In [16], however, only the estimation method of the disturbance is introduced and no specific control method for overcoming the disturbance using the estimated disturbance is proposed. In [17], inverse kinematics rather than inverse dynamics is used in the process of generating disturbance compensation signal. This approach can cause severe degradation in disturbance estimation performance as the dynamics is not negligible. In [18], which is the preliminary research of this paper, the structure of DOB to cope with translational force disturbance is proposed. However, accurate translational acceleration control is not achieved because of the error in converting the target translational force command to the target attitude and total thrust. Also, the nominal model used in DOB is based on inaccurate desired acceleration-to-desired states conversion technique.

In this paper, we present a new accurate three-dimensional translational acceleration tracking control that overcomes the limitations of the previous studies. The contributions of the proposed acceleration control technique are as follows. First, we introduce a new conversion method that reflects the difference between attitude dynamics and thrust dynamics when computing the target attitude and total thrust command from the translational force command (i.e., ‘Converter #2’ block of Fig. 1). Second, we model the translational force system (i.e., the shaded part of Fig. 1) that includes the new command conversion method, and design a DOB-based robust controller (i.e., ‘Disturbance Estimator’ block of Fig. 1) that overcomes translational force disturbance based on our new model. In the DOB controller design process, we perform $\mu$-analysis to systematically take into account the complex effects of various uncertainty. By presenting simulation and experimental results, we demonstrate the target acceleration tracking performance of the proposed conversion technique and the ability to overcome the external translational force disturbance of the designed DOB controller.

This paper is organized as follows. In Section II, we discuss the mathematical model of the multi-rotor used in the controller design. Section III deals with the force control of multi-rotor, and Section IV describes how DOB is applied to the force control. Section V provides the stability analysis to determine the range of DOB parameters that guarantee the stability of the designed system even in the presence of various uncertain elements. In Section VI, we demonstrate the empirical validity through simulations and actual experiments.

I MODELLING OF MULTI-ROTOR UAV

The rigid body dynamics of the multi-rotor are given by

[TABLE]

where $m$ is the mass of the multi-rotor, $\mathbf{X}=[x\ y\ z]^{T}\in\mathbb{R}^{3\times 1}$ is the position in the earth fixed frame, $\mathbf{F}=[F_{x}\ F_{y}\ F_{z}]^{T}\in\mathbb{R}^{3\times 1}$ is the three-dimensional translational force vector generated by the multi-rotor, $R(\mathbf{q})$ is the rotation matrix from the body frame to earth fixed frame, $\mathbf{q}=[\phi\ \theta\ \psi]^{T}\in\mathbb{R}^{3\times 1}$ is an attitude of the multi-rotor in the earth fixed frame, $\mathbf{T}_{t}=[0\ 0\ -T_{t}]^{T}\in\mathbb{R}^{3\times 1}$ is the thrust force vector in the body frame, $T_{t}\in\mathbb{R}$ is the magnitude of the total thrust, and $\mathbf{g}=[0\ 0\ g]^{T}\in\mathbb{R}^{3\times 1}$ is a gravity vector. The parameter $J\in\mathbb{R}^{3\times 3}$ is the moment of inertia (MOI) of the multi-rotor, $\mathbf{\Omega}=[p\ q\ r]^{T}\in\mathbb{R}^{3\times 1}$ is an angular velocity vector defined in the body frame, and $\mathbf{T}_{r}=[\tau_{r}\ \tau_{p}\ \tau_{y}]^{T}\in\mathbb{R}^{3\times 1}$ is an attitude control torque vector. For attitude dynamics, simplified dynamics of

[TABLE]

is commonly used, taking into account the small operation range of roll and pitch angle of multi-rotor and negligible Coriolis term [1, 17, 19].

II TRANSLATIONAL FORCE/ACCELERATION CONTROL

In order to control the translational force/acceleration of the multi-rotor, we need to convert the target acceleration $\ddot{\mathbf{X}}_{d}$ into the target attitude $\mathbf{q}_{d}$ and the target thrust $T_{t,d}$. Throughout this paper, notation $(*)_{d}$ denotes the desired value of the variable $*$. Also, we assume that the yaw $\psi$ of $\mathbf{q}$ always remains zero through a well-behaved independent controller to simplify the discussion. Now, we define $\mathbf{r}=[\theta\ \phi\ T_{t}]^{T}\in\mathbb{R}^{3\times 1}$ as a set of states that needs to be controlled for generating the desired translational acceleration of the multi-rotor.

Once we choose $\mathbf{r}=[\theta\ \phi\ T_{t}]^{T}$ as a set of state variables to control the translational force/acceleration of multi-rotor, our next task should be finding a way to convert the desired acceleration $\ddot{\mathbf{X}}_{d}$ to $\mathbf{r}_{d}$. To figure out how to convert the signal, let us first investigate the relationship between $\mathbf{r}$ and $\ddot{\mathbf{X}}$.

II-A Relationship between $\mathbf{r}$ and $\tilde{\ddot{\mathbf{X}}}$

In Equation (1), we have discussed the dynamics of the translational motion of multi-rotor. Going into detail, the corresponding translational dynamics are expressed as

[TABLE]

where $R(\psi)\in\mathbb{R}^{3\times 3}$ is the yaw rotation matrix. Now, let us define a vector of state variables $\tilde{\ddot{\mathbf{X}}}$ named the pseudo-acceleration vector as

[TABLE]

Applying Equation (4) to (3), we obtain the following relationship between $\mathbf{r}$ and $\tilde{\ddot{\mathbf{X}}}$:

[TABLE]

II-B Calculation of $\mathbf{r}_{d}$ from $\tilde{\ddot{\mathbf{X}}}_{d}$ considering dynamics

From Equation (5), we begin a discussion on how to calculate $\mathbf{r}_{d}$ based on $\tilde{\ddot{\mathbf{X}}}_{d}$. First, Equation (5) yields the following expression on $\mathbf{r}$:

[TABLE]

Equation (6) represents the required states $\mathbf{r}$ to generate such translational acceleration. From this, one might try to find the input to the controller to create the desired acceleration by replacing $\tilde{\ddot{\mathbf{X}}}$ and $\mathbf{r}$ with $\tilde{\ddot{\mathbf{X}}}_{d}$ and $\mathbf{r}_{d}$, respectively, as follows.

[TABLE]

However, this method can severely degrade control performance when multi-rotor is larger than a certain size as we discuss below.

Fig. 2 shows the internal structure between $\mathbf{r}_{d}$ and $\mathbf{r}$. In this figure, we can see that $\phi_{d}$ and $\theta_{d}$ are realized to $\phi$ and $\theta$ through attitude controller, rotor dynamics, and attitude dynamics. In contrast, $T_{t,d}$ only passes through the rotor dynamics to become $T_{t}$. Here, we treat $\mathbf{u}_{d}=\mathbf{u}$, where $\mathbf{u}=[\mathbf{T}_{r}\ T_{t}]^{T}\in\mathbb{R}^{4\times 1}$, since rotor dynamics are mostly negligible. Assuming that the attitude controller is properly designed, we can model the relationship between $\mathbf{r}_{d}$ and $\mathbf{r}$ as the following equation:

[TABLE]

Here, $\gamma_{*}\in[0,\infty)$ are time-varying non-negative delay factors. Applying Equation (10) into (5), we have

[TABLE]

In Equation (11), the desired attitude and total thrust are realized asynchronously due to $\gamma_{\phi}$ and $\gamma_{\theta}$. Applying Equation (9) to Equation (11), the result is follows.

[TABLE]

In the $\tilde{\ddot{z}}(t)$ equation of Equation (12), the parenthesized part can continuously change if $\gamma_{\phi}$ and $\gamma_{\theta}$ are too large to be ignored. This indicates that $z$-directional control performance can be significantly reduced if the delay between the desired and actual attitude signals becomes large, for example in situations when the MOI of the multi-rotor increases, such as large multi-rotor or multi-rotor with large cargo. When the Z-directional control performance degrades, a high-level controller (e.g., position controller) or the operator may need to constantly modify the $\tilde{\ddot{z}}_{d}$ value to correct the poor Z-directional control performance. As a result, this degrades the X and Y direction control performance because the values in parentheses of the $\tilde{\ddot{x}}(t)$ and $\tilde{\ddot{y}}(t)$ equations in (12) also constantly change. The decline in control performance due to this control scheme will be shown in Fig. 3.

To address this issue, we next consider two candidate solutions.

II-B1 Solution candidate 1

The first candidate is to time-synchronize the attitude and total thrust output by adding an artificial time delay to $T_{t,d}$ in Equation (9) as

[TABLE]

Here, $\gamma_{v}$ is a delay element deliberately applied to $\tilde{\ddot{z}}_{d}$. Applying Equation (13) to Equation (11), the equation of motion is changed from Equation (12) to

[TABLE]

Through Equations (7), (8) and (10), $\phi(t)$ and $\theta(t)$ can be described as

[TABLE]

Let us assume that $\gamma_{\phi}$ and $\gamma_{\theta}$ have the same value of $\gamma_{h}$ since most multi-rotors have nearly the same roll and pitch behavior due to the symmetrical mechanical structure. Then, Equation (14) with Equations (10) and (15) becomes as

[TABLE]

Now, we can solve the problem in Equation (12) by setting $\gamma_{v}$ equal to $\gamma_{h}$. However, this method is not easily applicable in a real-world situation because it is difficult to determine the value of $\gamma_{h}$ that changes continuously during the flight. Therefore, the control method through Equation (13) cannot be a practical method.

II-B2 Solution candidate 2

Alternatively, we can find a reasonable solution that is applicable in the real world by selectively delaying $\phi_{d}(t)$ and $\theta_{d}(t)$ in Equation (13) by $\gamma_{\phi}$ and $\gamma_{\theta}$, but keeping $\gamma_{v}$ at zero. As we can see from Equation (10), the values of $\phi_{d}(t)$ and $\theta_{d}(t)$ delayed by $\gamma_{\phi}$ and $\gamma_{\theta}$ seconds are $\phi(t)$ and $\theta(t)$. Applying this idea to Equation (13), we can obtain $T_{t,d}$ as

[TABLE]

where the values $\phi(t)$ and $\theta(t)$ can be measured from the built-in inertial measurement unit (IMU) sensor. Then, by setting $\gamma_{v}$ to zero, we can determine the input/output relationship of the translational accelerations dynamics of the multi-rotor as

[TABLE]

where we assume ${{\tilde{\ddot{z}}_{d}(t)}\over{\tilde{\ddot{z}}_{d}(t-\gamma_{h})}}\approx 1$. This assumption is valid in most cases, except in situations where the change in target vertical acceleration is abnormally large and rapid.

Through the control techniques of solution candidate 2 (Equations (7), (8) and (17)), we obtained a three-dimensional translational acceleration control method applicable to actual multi-rotor control. In order to compare the performance of multi-rotor control using Equations (9) and (17), a brief simulation is conducted as shown in Fig. 3.

The simulation shows the comparison of the target acceleration tracking performance of Case 1 with Equations (7), (8), (9) and Case 2 with Equations (7), (8), (17). The upper set of figures show the acceleration tracking performance of Cases 1 and 2 with arbitrary acceleration command. Here, we can see that there are no differences in performance between Cases 1 and 2 when MOI of the multi-rotor has small value of 0.1. On the other hand, when the MOI of the multi-rotor increases, both Cases 1 and 2 show delayed responses in the X and Y direction acceleration tracking as expected. However, we can observe that the Z-directional performance of the Case 2 remains the same regardless of the magnitude of the MOI, unlike Case 1 where the performance degradation is observed. The effect of the decline in Z-directional control performance on the system is evident when controlling the position of the multi-rotor. The bottom set of figures is the situation where the high-level position controller generates the desired acceleration command to track the predefined trajectory. In Case 1, we can observe a decrease in acceleration tracking performance in both the X and Y directions as well as the Z direction as the MOI increases. On the other hand, in Case 2, the Z-directional control performance remains constant regardless of the MOI of the platform, stabilizing the X and Y-directional control performance faster than Case 1.

This phenomenon can be understood in other ways by considering the role of the denominator term of the $T_{t}$ equation in Equation (6), which is to compensate for the reduction of the vertical thrust component in the sense of inertial coordinates when the multi-rotor is tilted. When $T_{t,d}$ is calculated based on the desired attitude as Equation (9), the situation is similar to compensating for the future event after $\gamma_{h}$ seconds. Instead, it is intuitive to use the current attitude as in Equation (17) to correct the vertical thrust reduction. From the flight results using Equation (17) in Fig. 3, we can confirm that the control performance in all directions is satisfactory.

III DISTURBANCE OBSERVER

External disturbances applied to multi-rotor act not only in the form of translational disturbances but also in the form of rotational torques. However, given that a number of solutions for overcoming the rotational torque disturbances [8]$\sim$[15] have already been proposed, this section concerns only translational disturbances applied to the system for straightforward discussion and analysis.

III-A An overview of the disturbance-merged overall system

Fig. 4 shows the overall configuration of the system. First, the position controller ${C(t)}$ generates the target acceleration input $\ddot{\mathbf{X}}_{d}$. This signal is then transformed into the target force input $\mathbf{F}_{d}$ through the following force-acceleration relationship:

[TABLE]

Then, $\mathbf{F}_{d}$ signal passes through ${{1}\over{m}}R^{-1}(\psi)$ block to transform the signal into the $\tilde{\ddot{\mathbf{X}}}_{d}$ (refer Equation (4)). The signal $\tilde{\ddot{\mathbf{X}}}_{d}$ then passes through the $\mathbf{B}\scriptscriptstyle{[\tilde{\ddot{\mathbf{X}}}_{d}\rightarrow\mathbf{r}_{d}]}$ block, which converts the target acceleration $\tilde{\ddot{\mathbf{X}}}_{d}$ to $\mathbf{r}_{d}$, the input to the multi-rotor controller, based on Equations (7), (8) and (17). Once $\mathbf{r}_{d}$ passes through the dynamics described in Fig. 2 and outputs $\mathbf{r}$, it passes through $\mathbf{B}\scriptscriptstyle{[\mathbf{r}\rightarrow\mathbf{F}]}$ block to produce $\mathbf{F}$ (refer Equation (4) and (5)). Right after $\mathbf{F}$ is generated, the external disturbance force $\mathbf{d}_{actual}$ immediately compromises the thrust and results in $\tilde{\mathbf{F}}$ and $\ddot{\mathbf{X}}$.

III-B Disturbance observer

In Fig. 4, the translational force disturbance $\mathbf{d}_{actual}$ is combined with $\mathbf{F}$ to become $\tilde{\mathbf{F}}$. However, canceling $\mathbf{d}_{actual}$ is only possible by adding an appropriate disturbance cancellation term to the $\mathbf{F}_{d}$ signal. Therefore, it is preferable to assume that there is an equivalent input disturbance $\mathbf{d}_{EID}$ that has the same effect on the system as $\mathbf{d}_{actual}$ [20]. Then $\mathbf{d}_{actual}$ is replaced by $\mathbf{d}_{EID}$, making $\mathbf{F}=\tilde{\mathbf{F}}$. As we can see in Fig. 4, the $\mathbf{d}_{EID}$ signal is merged into $\tilde{\mathbf{F}}_{d}$, which is the translational acceleration control input with disturbance cancellation signal. Now, let us construct the DOB based on the above settings.

III-B1 $\mathbf{d}_{EID}$ estimation algorithm

For the estimation of $\mathbf{d}_{EID}$, we first estimate $\kappa$ the sum of $\tilde{\mathbf{F}}_{d}$ and $\mathbf{d}_{EID}$ by

[TABLE]

We can easily achieve the $\tilde{\mathbf{F}}$ signal from Equation (19) where $\ddot{\mathbf{X}}$ is measured by the IMU sensor. The transfer function ${P_{n}(s)}$ is the nominal model of ${P(s)}$, and $\hat{(*)}$ is the representation of the estimation of $(*)$ signal throughout this paper. Once we estimate $\hat{\kappa}$, we then obtain $\hat{\mathbf{d}}_{EID}$ by

[TABLE]

The signal $\hat{\kappa}(s)$ passes through the ${Q_{1}}$ block, which is basically a low pass filter, to overcome both the causality violation issue due to the improperness of ${P_{n}^{-1}}(s)$ and the potential instability issue caused by the non-minimum phase characteristic of ${P_{n}}(s)$. The filter ${Q_{2}}(s)$ is used to match the phase with ${Q_{1}}(s)\hat{\kappa}(s)$ signal. In the end, we generate a disturbance-compensating control input $\tilde{\mathbf{F}}_{d}$ by

[TABLE]

This makes $\kappa$ become

[TABLE]

The most important factor in the $\mathbf{d}_{EID}$ estimation process is the proper design of ${P_{n}}$ and ${Q}$. Of these, ${Q}$ is deeply related to the stability of the system and will be discussed in more detail in the next section. In the remainder of this section, we first discuss the design of the nominal model ${P_{n}}$ and then explain the structure of the $Q$-filter.

III-B2 Nominal model ${P_{n}}$

The internal structure of ${P_{n}^{-1}}(t)$ is described as in Fig. 5, all of which are simple conversion blocks except for the ${{\it{\Lambda}}}^{-1}_{n}{\it(t)}$ block. The block ${\it{\Lambda}}(t)$ is the relationship between $\mathbf{r}_{d}$ and $\mathbf{r}$ depicted in Fig. 2. The $\it{\Lambda_{n}(s)}$ is constructed from two parts: attitude and thrust dynamics. We denote these as $\it{\Lambda_{n,a}}(s)$ and $\it{\Lambda_{n,t}}(s)$ respectively.

As we see from Fig. 2, $\it{\Lambda}_{n,a}(s)$ is constructed with attitude controller, rotor dynamics and attitudinal dynamics. Since rotor dynamics can be ignored, we only need to find the transfer function of the attitudinal dynamics and attitude controller. For attitude dynamics, let us refer to Equation (2) and express it as

[TABLE]

where $i=1,2,3$ represent $\phi$, $\theta$, $\psi$ axis, respectively. For attitude control, PD control in the following form is used.

[TABLE]

The parameters $P_{i}$, $D_{i}$ represent control gains in each attitude component. Then, the overall transfer function ${\it{\Lambda}}_{n,a,i}$ between desired and current attitude becomes

[TABLE]

In the case of ${\it{\Lambda}}_{n,t}$, the only dynamics involved is rotor dynamics, which we decided to neglect. Thus, it can be expressed as

[TABLE]

Now, we can construct the transfer matrix for ${\it{\Lambda}}_{n}=diag({\it{\Lambda}}_{n,1},{\it{\Lambda}}_{n,2},{\it{\Lambda}}_{n,3})$ using Equations (26) and (27) as

[TABLE]

Equation (28) is a detailed representation of the relationship between $\mathbf{r}_{d}$ and $\mathbf{r}$, which was introduced in Equation (10). On the other hand, $P_{n}$, which defines the nominal relationship between $\kappa$ and $\mathbf{F}$ (or $\tilde{\ddot{\mathbf{X}}}_{d}$ and $\tilde{\ddot{\mathbf{X}}}$), was introduced in Equation (18) (refer Equation (4) for the relationship between $\tilde{\ddot{\mathbf{X}}}$ and $\mathbf{F}$). Here, we can see that both Equations (10) and (18) have the same input/output characteristics with time delay of $\gamma_{h}$ for the first and second channels and no time delay for the third channel. Therefore, we can conclude that ${\it{\Lambda}}_{n}(s)$ in Equation (28) is also the transfer function between $\kappa$ and $\mathbf{F}$ as well as between $\mathbf{r}_{d}$ and $\mathbf{r}$, which is

[TABLE]

III-B3 $Q$-filter design

In $Q$-filter design, we choose to make ${Q_{1}}(s){\it{\Lambda}}_{n}^{-1}(s)$, which is now identical to ${Q_{1}}(s){P}_{n}^{-1}(s)$, a proper function with relative degree of 1. Since ${P}_{n}(s)$ is composed of three channels in $X$, $Y$ and $Z$ directions, we need to design three separate $Q$-filters. As shown in Equation (26), ${\it{\Lambda}}_{n,1}(s)$($={\it{\Lambda}}_{n,a,2}(s)$) and ${\it{\Lambda}}_{n,2}(s)$($={\it{\Lambda}}_{n,a,1}(s)$) among the three transfer functions of ${\it{\Lambda}}_{n}(s)$ are systems with a relative degree of 1. The thrust transfer function ${\it{\Lambda}}_{n,3}(s)$($={\it{\Lambda}}_{n,t}(s)$) has a relative degree of 0, as can be seen from Equation (27). Therefore, the $Q$-filters for making ${Q_{1}}(s){\it{\Lambda}}_{n}^{-1}(s)$ with a relative degree of 1 are designed as

[TABLE]

[TABLE]

[TABLE]

where ${Q_{1,h}}$ and ${Q_{1,v}}$ are $Q$-filters corresponding to the horizontal ($\it{\Lambda}_{n,a}$) and vertical ($\it{\Lambda}_{t}$) models respectively. The symbol $\tau$ is the time constant and $\zeta$ is the damping ratio of the filter. The filter ${Q_{1}}$ is designed to have a gain of 1 when $s=0$ [21]. The filter ${Q_{2}}$ is set to ${Q_{2}}={Q_{1}}$, to easily achieve the purpose of phase matching.

IV STABILITY ANALYSIS

The design of $Q$-filter in the DOB structure should be based on rigorous stability analysis to ensure the overall stability. In particular, we note that there is always a difference between the nominal model $P_{n}(s)$ and the actual model $P(s)$, due to various uncertainties and applied assumptions.

Although the small-gain theorem (SGT) [18] can still be a tool for stability analysis, the SGT analysis based on the largest singular value among uncertainties is likely to yield overly conservative results especially if multiple uncertain elements are involved. Instead, we use structured singular value analysis, or $\mu$-analysis [22, 23, 24], to reflect the combined effects of uncertainties.

IV-A Modeling of ${P}(s)$ considering uncertainties

The multi-rotor’s actual transfer function ${P}(s)$ between $\kappa$ and $\mathbf{F}$ in Fig. 4 is

[TABLE]

Here, $P_{1}$, $P_{2}$ and $P_{3}$ represent the input/output translational force relationship in the $X$, $Y$, and $Z$ directions, respectively. This research considers a small but nonzero DC-gain error, parametric error and phase shift error between ${P}_{n}(s)$ and ${P}(s)$. Then each ${P}_{j}(s)$ can be expressed as the following equation:

[TABLE]

where $j=1,2,3$ represent $X$, $Y$, $Z$ axis. The symbols $K_{j},\delta_{j}\in\mathbb{R}$ represent the uncertain variable gain and time delay parameters, respectively. The nominal transfer function ${P}_{n,j}$ can be replaced by $\it{\Lambda}_{n,j}$ based on Equation (29). The portion containing only the parametric uncertainty is denoted by ${P}_{p,j}(s)=K_{j}{\it{\Lambda}}_{n,j}(s)$, and the time delay uncertainty is denoted by $\Gamma_{j}(s)=e^{-\delta_{j}s}$.

In Equation (34), each ${P}_{j}(s)$ contains three uncertain variables, which are $K_{j}$, $J_{j}$ and $\delta_{j}$. In the case of $K_{j}$, we define $K_{j}$ as

[TABLE]

where $K_{\Delta,j}\in\mathbb{R}$ is the error value of $K_{j}$. In the case of $J_{j}$, determining the actual value of $J_{j}$ is difficult compared to other physical quantities. We also define $J_{j}$ in the same manner as $K_{j}$ for the convenience of analysis as

[TABLE]

where $\bar{J}_{j},J_{\Delta,j}\in\mathbb{R}$ are the nominal and error values of $J_{j}$. Because the term $\Gamma_{j}(s)$ containing $\delta_{j}$ is of an irrational form that is not suitable for analysis, we use an analytic approximation of the uncertain time-delay $\Gamma_{j}(s)$ to a rational function with unmodeled dynamic uncertainty [24]. First, we change the representation of the $P_{j}(s)$ model to a multiplicative uncertainty form that combines parametric uncertainties and unmodeled time-delay uncertainty as follows:

[TABLE]

A complex unstructured uncertainty $\Delta_{\delta,j}\in\mathbb{C}$ corresponds to unknown time delay $\delta_{j}$, and $W_{\delta,j}(s)$ is the maximum uncertainty that can be caused by $\Gamma_{j}(s)$. Here, we can obtain $W_{\delta,j}(s)$ using Equation (34) as

[TABLE]

The maximum value of $|\Gamma_{j}(j\omega)-1|$ for each $\omega$ can be found using Euler’s formula as

[TABLE]

where

[TABLE]

As a result of analyzing a large amount of actual experimental data, we confirmed that the time delay between ${P}_{n}(s)$($={\it{\Lambda}}_{n}(s)$) and ${P}(s)$ does not exceed 0.1 second in all three channels. We put 20 percent margin so that $|\delta_{j}|\leq 0.12$. Fig. 6 is multiple Bode magnitude plots of $|\Gamma(s)-1|$ generated by varying $\delta$ from $-0.12$ to $+0.12$. From Fig. 6, we can extract

[TABLE]

for all $j$, which is the upper boundary of $|\Gamma_{j}(s)-1|$ sets represented by the red solid line.

The uncertainties of $K_{j}$ and $J_{j}$ can also be modeled in the same manner as in Equation (37) as

[TABLE]

where $||{\Delta}_{K,j}||_{\infty},||{\Delta}_{J,j}||_{\infty}\leq 1$. The transfer function ${\it{\Lambda}}_{n,n,j}$ is basically the same as ${\it{\Lambda}}_{n,j}$, except that $J$ in Equation (26) is replaced to the nominal MOI value $\bar{J}$. The transfer functions $W_{K,j}$ and $W_{J,j}$ are

[TABLE]

IV-B $\tau$-determination through $\mu$-analysis

IV-B1 $\mu$-robust stability analysis

In [25], the structured singular value $\mu$ is defined as

[TABLE]

where $\Delta$ is a complex structured block-diagonal unmodeled uncertainty block which gathers all model uncertainties [26]. Following the common notation, the symbol $\mathbf{\Delta}$ represents a set of all stable transfer matrices with the same structure (full, block-diagonal, or scalar blocks) and nature (real or complex) as $\Delta$. The \bar{\sigma}\big{(}\Delta\big{)} is the maximum singular value of uncertainty block $\Delta$. The matrices $M$ and $\Delta$ are defined by collapsing the simplified overall system to upper LFT uncertainty description as

[TABLE]

where $M$ is the known part of the system, $r$ is a reference input and $y$ is an output of the overall system. In the theory of the $\mu$-analysis, it is well-known that the system is robustly stable if $\mu$ satisfies the following conditions

[TABLE]

[22][25].

The $\mu$-analysis is performed separately for each channel of $X$, $Y$, $Z$ thanks to the structure of the platform described by Equation (28), but since $X$ and $Y$ channels are composed of the same structure, they share the identical analysis result. As we can see from Fig. 7, the system is collapsed in the form of Equation (48) by using MATLAB’s Robust Control Toolbox${}^{\textrm{TM}}$, where $r_{j}=[\mathbf{F}_{d,j}\ \mathbf{d}_{EID,j}]^{T}\in\mathbb{R}^{2\times 1}$ and $y_{j}=\tilde{\mathbf{F}}_{j}\in\mathbb{R}$ in our case. As a reminder, subscript $j$ refers to each channel of $X$, $Y$, and $Z$. Also, structured uncertainty ${\Delta}_{j}(s)\in\mathbb{C}^{3\times 3}$ is constructed as

[TABLE]

which includes unmodeled MOI uncertainty, time and gain uncertainty in our system.

IV-B2 Results of analysis

Table I shows the multi-rotor’s physical quantities and controller gains used both in the simulation and the experiment. The gains $P_{\phi,\theta}$ and $D_{\phi,\theta}$ are predefined values set during the primary gain-tuning process to obtain the ability to control the attitude of the platform. The translational acceleration limit is set to prevent flight failure due to excessive acceleration control inputs and is set at $\pm 3\$\mathrm{m}\mathrm{/}\mathrm{s}^{2} to have a roll and pitch limit of approximately $\pm 0.3\ $\mathrm{r}\mathrm{a}\mathrm{d} in level flight condition. As previously mentioned, the unmodeled time delay $\delta_{j}$ is set to 0.1, and the gain error is assumed to be a maximum error of 10 percent. For MOIs that are difficult to estimate, we assumed a wider 30 percent uncertainty. The damping ratio $\zeta$ of the second order filter is set to 0.707, which is the critical damping ratio, to balance the overshoot and late response. Fig. 8 shows the results of $\mu$-analysis. From the analysis, we can see that the system is stable when $\tau_{1}>0.12$ and $\tau_{2}>0.09$.

Fig. 9 shows the results of the SGT-based stability analysis, performed in the same manner as [18]. The analysis is based on the following model:

[TABLE]

where all uncertainties due to $\delta_{j}$, $K_{j}$ and $J_{j}$ are lumped using the functions $W_{l,h}(=W_{l,1},W_{l,2})$ and $W_{l,v}(=W_{l,3})$, whose magnitude increases over frequency as shown in blue curves of Fig. 9. The stability condition of the SGT-based analysis in this case is

[TABLE]

[18, 22, 27]. In the SGT-based analysis, the bode plots of the $Q$-filter with $\tau_{1}=0.12$ and $\tau_{2}=0.09$ indicate that system with those $\tau$ values could be unstable. However, through the $\mu$-analysis, those $\tau$ values are still in the stable region. From this, we can confirm that the $\mu$-analysis provides more rigorous $\tau$ boundary values than SGT-based analysis.

V SIMULATION AND EXPERIMENTAL RESULT

This section reports simulation and experimental results to validate the performance of our three-dimensional force controller and the disturbance cancellation performance of the DOB technique. The comparison of the acceleration tracking performance of the force control methods according to the MOI variation is already shown in the simulation of Fig. 3. Therefore, in this section, we provide

experimental result to demonstrate the performance of the proposed force control technique for the actual plant, and 2. 2.

simulation and experimental results to demonstrate the capability of the DOB in overcoming the translational force disturbance.

Based on the results from the previous section, the cutoff frequencies of the $Q$-filter are set to $\tau_{1}=0.15$ and $\tau_{2}=0.12$ in both simulation and actual experiment with additional margins to ensure additional stability.

V-A Validation of acceleration tracking performance

In the experiment, arbitrary desired acceleration commands for $X$ and $Y$ directions are given by the operator-controlled radio controller. Fig. 10 shows the multi-rotor accurately following the target acceleration. From this result, we can confirm that our three-dimensional translational acceleration control technique functions effectively even in the actual flight.

V-B Validation of DOB performance

V-B1 Simulation result

In the simulation, the multi-rotor follows a circular trajectory with radius of 3 $\mathrm{m}$ and height of 5 $\mathrm{m}$. Meanwhile, the multi-rotor is exposed to periodic disturbances with accelerations up to 5.5 $\mathrm{m}\mathrm{/}\mathrm{s}^{2}$ in each axis. Fig. 11 compares the multi-rotor’s position tracking performance before and after applying DOB. On the left graphs of Fig. 11, the target trajectory tracking results are not smooth due to the unexpected disturbances, whereas the trajectory deviation is drastically reduced in the right graphs where the DOB algorithm is applied.

V-B2 Experimental Result

In the experiment, the multi-rotor is commanded to hover at a specific point in three-dimensional space but connected to the translational force measurement sensor via the tether to measure the applied disturbance force. As we can see in Fig. 13, the operator aligns the force sensor in the $X$-axis and pulls and releases the force sensor periodically to apply a disturbance to the multi-rotor.

Fig. 12 is a comparison of hovering performance before (left) and after (right) applying the DOB algorithm. The graphs in the left column are the case when the DOB is not applied, which has a larger $X$ directional position shift than other axes. Unlike the DOB-off case, the DOB-on case shows a significant reduction in position error.

Two graphs at the forth row shows the acceleration tracking results. When DOB is not applied, an acceleration signal is generated by the position error, but we can see that the target acceleration cannot be followed due to the disturbance. Meanwhile, we can see that the acceleration of the platform (yellow solid line) well tracks the target acceleration (blue dash-single dotted line). This is because the well-behaved DOB algorithm generated control input including the disturbance compensation signal (orange dash-single dotted line) and applied to the platform. The effect of the DOB can be confirmed by significantly reduced position error. Four graphs at the bottom of the figure show the difference between the signal $\vec{\mathcal{F}}_{d}$ and $\tilde{\vec{\mathcal{F}}}_{d}$ (fifth row), and the comparison between $d_{actual,x}$ measured by force sensor and $\hat{d}_{EID,x}$ estimated by DOB algorithm (sixth row). When DOB is not applied, $\hat{d}_{EID}$ estimation process is working internally but the signal is not merged into $\tilde{\vec{\mathcal{F}}}_{d}$ signal, making $\vec{\mathcal{F}}_{d}$ and $\tilde{\vec{\mathcal{F}}}_{d}$ have the same value. On the other hand, we can see the difference between the $\vec{\mathcal{F}}_{d}$ and the $\tilde{\vec{\mathcal{F}}}_{d}$ signal when DOB is applied, because the $\hat{d}_{EID}$ signal is merged into the $\tilde{\vec{\mathcal{F}}}_{d}$ signal. Two graphs in the last row show the comparison between the measured disturbance and the estimated disturbance, and we can confirm that the estimates are fairly accurate in both cases.

An extra flight experiment is conducted under wind disturbance to validate the DOB performance in a more realistic environment. As we can see in Fig. 14, the target location of the multi-rotor is set on the centerline of a strong wind generator that generates wind speed of 7 $\mathrm{m}\mathrm{/}\mathrm{s}$. The performance of DOB is visualized by comparing the position difference between DOB-on and DOB-off situations. the multi-rotor has a position error of about 1 $\mathrm{m}$ in the DOB-off case and about 0.3 $\mathrm{m}$ in the DOB-on case. Through the experiment, we can confirm that the proposed DOB algorithm works effectively even against a wind disturbance.

VI CONCLUSION

In this paper, we introduced 1) a new method of converting the target acceleration command to the desired attitude and total thrust, and 2) a DOB method for overcoming the disturbance that obstructs the translational motion, to more accurately control the translational acceleration of the multi-rotor UAV. In the control input conversion process, we reflect the different dynamic characteristics of attitude and thrust, so that more precise control is possible than the existing methods. Then, by using the DOB-based robust control algorithm based on the nominal translational force system, the magnitude of the disturbance force applied to the fuselage is estimated and compensated. For the robust stability guarantee, the $Q$-filter of the DOB is designed based on the $\mu$-stability analysis. The validity of the proposed method is confirmed through simulation and actual experiments.

The proposed technique is useful in various applications such as aerial parcel delivery service or drone-based industrial operations where precise acceleration control is required. For example, in a multi-rotor-based parcel delivery service, the proposed DOB algorithm can maintain the nominal flight performance by considering the additional force due to the weight of the cargo attached to the multi-rotor as a disturbance to be estimated. Also, the proposed algorithm is suitable for situations that require precise trajectory tracking performance even in windy conditions such as maritime operations or human-rescue missions. For industrial applications involving collaborative flight of multiple multi-rotors, the proposed algorithm can be used to estimate and stabilize internal forces caused in between physically-coupled multi-rotors.