I Modelling of MultiRotor Uav
The rigid body dynamics of the multirotor are given by
(1) 
where is the mass of the multirotor, is the position in the earth fixed frame,
is the threedimensional translational force vector generated by the multirotor,
is the rotation matrix from the body frame to earth fixed frame, is an attitude of the multirotor in the earth fixed frame, is the thrust force vector in the body frame, is the magnitude of the total thrust, and is a gravity vector. The parameter is the moment of inertia (MOI) of the multirotor, is an angular velocity vector defined in the body frame, and is an attitude control torque vector. For attitude dynamics, simplified dynamics of(2) 
is commonly used, taking into account the small operation range of roll and pitch angle of multirotor and negligible Coriolis term [2, 5, 11].
Ii Translational Force/acceleration Control
In order to control the translational force/acceleration of the multirotor, we need to convert the target acceleration into the target attitude and the target thrust . Throughout this paper, notation denotes the desired value of the variable . Also, we assume that the yaw of always remains zero through a wellbehaved independent controller to simplify the discussion. Now, we define as a set of states that needs to be controlled for generating the desired translational acceleration of the multirotor.
Once we choose as a set of state variables to control the translational force/acceleration of multirotor, our next task should be finding a way to convert the desired acceleration to . To figure out how to convert the signal, let us first investigate the relationship between and .
Iia Relationship between and
In Equation (1), we have discussed the dynamics of the translational motion of multirotor. Going into detail, the corresponding translational dynamics are expressed as
(3) 
where is the yaw rotation matrix. Now, let us define a vector of state variables named the pseudoacceleration vector as
(4) 
Applying Equation (4) to (3), we obtain the following relationship between and :
(5) 
IiB Calculation of from considering dynamics
From Equation (5), we begin a discussion on how to calculate based on . First, Equation (5) yields the following expression on :
(6) 
Equation (6) represents the required states to generate such translational acceleration. From this, one might try to find the input to the controller to create the desired acceleration by replacing and with and , respectively, as follows.
(7)  
(8)  
(9) 
However, this method can severely degrade control performance when multirotor is larger than a certain size as we discuss below.
Fig. 2 shows the internal structure between and . In this figure, we can see that and are realized to and through attitude controller, rotor dynamics, and attitude dynamics. In contrast, only passes through the rotor dynamics to become . Here, we treat , where , since rotor dynamics are mostly negligible. Assuming that the attitude controller is properly designed, we can model the relationship between and as the following equation:
(10) 
Here, are timevarying nonnegative delay factors. Applying Equation (10) into (5), we have
(11) 
In Equation (11), the desired attitude and total thrust are realized asynchronously due to and . Applying Equation (9) to Equation (11), the result is follows.
(12) 
In the equation of Equation (12), the parenthesized part can continuously change if and are too large to be ignored. This indicates that 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 multirotor increases, such as large multirotor or multirotor with large cargo. When the Zdirectional control performance degrades, a highlevel controller (e.g., position controller) or the operator may need to constantly modify the value to correct the poor Zdirectional control performance. As a result, this degrades the X and Y direction control performance because the values in parentheses of the and 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.
IiB1 Solution candidate 1
The first candidate is to timesynchronize the attitude and total thrust output by adding an artificial time delay to in Equation (9) as
(13) 
Here, is a delay element deliberately applied to . Applying Equation (13) to Equation (11), the equation of motion is changed from Equation (12) to
(14) 
Through Equations (7), (8) and (10), and can be described as
(15) 
Let us assume that and have the same value of since most multirotors have nearly the same roll and pitch behavior due to the symmetrical mechanical structure. Then, Equation (14) with Equations (10) and (15) becomes as
(16) 
Now, we can solve the problem in Equation (12) by setting equal to . However, this method is not easily applicable in a realworld situation because it is difficult to determine the value of that changes continuously during the flight. Therefore, the control method through Equation (13) cannot be a practical method.
IiB2 Solution candidate 2
Alternatively, we can find a reasonable solution that is applicable in the real world by selectively delaying and in Equation (13) by and , but keeping at zero. As we can see from Equation (10), the values of and delayed by and seconds are and . Applying this idea to Equation (13), we can obtain as
(17) 
where the values and can be measured from the builtin inertial measurement unit (IMU) sensor. Then, by setting to zero, we can determine the input/output relationship of the translational accelerations dynamics of the multirotor as
(18) 
where we assume . 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 threedimensional translational acceleration control method applicable to actual multirotor control. In order to compare the performance of multirotor 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 multirotor has small value of 0.1. On the other hand, when the MOI of the multirotor 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 Zdirectional 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 Zdirectional control performance on the system is evident when controlling the position of the multirotor. The bottom set of figures is the situation where the highlevel 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 Zdirectional control performance remains constant regardless of the MOI of the platform, stabilizing the X and Ydirectional control performance faster than Case 1.
This phenomenon can be understood in other ways by considering the role of the denominator term of the equation in Equation (6), which is to compensate for the reduction of the vertical thrust component in the sense of inertial coordinates when the multirotor is tilted. When is calculated based on the desired attitude as Equation (9), the situation is similar to compensating for the future event after 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 multirotor 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 [18][10] have already been proposed, this section concerns only translational disturbances applied to the system for straightforward discussion and analysis.
Iiia An overview of the disturbancemerged overall system
Fig. 4 shows the overall configuration of the system. First, the position controller generates the target acceleration input . This signal is then transformed into the target force input through the following forceacceleration relationship:
(19) 
Then, signal passes through block to transform the signal into the (refer Equation (4)). The signal then passes through the block, which converts the target acceleration to , the input to the multirotor controller, based on Equations (7), (8) and (17). Once passes through the dynamics described in Fig. 2 and outputs , it passes through block to produce (refer Equation (4) and (5)). Right after is generated, the external disturbance force immediately compromises the thrust and results in and .
IiiB Disturbance observer
In Fig. 4, the translational force disturbance is combined with to become . However, canceling is only possible by adding an appropriate disturbance cancellation term to the signal. Therefore, it is preferable to assume that there is an equivalent input disturbance that has the same effect on the system as [21]. Then is replaced by , making . As we can see in Fig. 4, the signal is merged into , which is the translational acceleration control input with disturbance cancellation signal. Now, let us construct the DOB based on the above settings.
IiiB1 estimation algorithm
For the estimation of , we first estimate the sum of and by
(20) 
We can easily achieve the signal from Equation (19) where is measured by the IMU sensor. The transfer function is the nominal model of , and is the representation of the estimation of signal throughout this paper. Once we estimate , we then obtain by
(21) 
The signal passes through the block, which is basically a low pass filter, to overcome both the causality violation issue due to the improperness of and the potential instability issue caused by the nonminimum phase characteristic of . The filter is used to match the phase with signal. In the end, we generate a disturbancecompensating control input by
(22) 
This makes become
(23) 
The most important factor in the estimation process is the proper design of and . Of these, 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 and then explain the structure of the filter.
IiiB2 Nominal model
The internal structure of is described as in Fig. 5, all of which are simple conversion blocks except for the block. The block is the relationship between and depicted in Fig. 2. The is constructed from two parts: attitude and thrust dynamics. We denote these as and respectively.
As we see from Fig. 2, 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
(24) 
where represent , , axis, respectively. For attitude control, PD control in the following form is used.
(25) 
The parameters , represent control gains in each attitude component. Then, the overall transfer function between desired and current attitude becomes
(26) 
In the case of , the only dynamics involved is rotor dynamics, which we decided to neglect. Thus, it can be expressed as
(27) 
Now, we can construct the transfer matrix for using Equations (26) and (27) as
(28) 
Equation (28) is a detailed representation of the relationship between and , which was introduced in Equation (10). On the other hand, , which defines the nominal relationship between and (or and ), was introduced in Equation (18) (refer Equation (4) for the relationship between and ). Here, we can see that both Equations (10) and (18) have the same input/output characteristics with time delay of for the first and second channels and no time delay for the third channel. Therefore, we can conclude that in Equation (28) is also the transfer function between and as well as between and , which is
(29) 
IiiB3 filter design
In filter design, we choose to make , which is now identical to , a proper function with relative degree of 1. Since is composed of three channels in , and directions, we need to design three separate filters. As shown in Equation (26), () and () among the three transfer functions of are systems with a relative degree of 1. The thrust transfer function () has a relative degree of 0, as can be seen from Equation (27). Therefore, the filters for making with a relative degree of 1 are designed as
(30) 
(31) 
(32) 
where and are filters corresponding to the horizontal () and vertical () models respectively. The symbol is the time constant and is the damping ratio of the filter. The filter is designed to have a gain of 1 when [8]. The filter is set to , to easily achieve the purpose of phase matching.
Iv Stability Analysis
The design of 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 and the actual model , due to various uncertainties and applied assumptions.
Although the smallgain theorem (SGT) [13]
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
analysis [7, 20, 9], to reflect the combined effects of uncertainties.Iva Modeling of considering uncertainties
The multirotor’s actual transfer function between and in Fig. 4 is
(33) 
Here, , and represent the input/output translational force relationship in the , , and directions, respectively. This research considers a small but nonzero DCgain error, parametric error and phase shift error between and . Then each can be expressed as the following equation:
(34) 
where represent , , axis. The symbols represent the uncertain variable gain and time delay parameters, respectively. The nominal transfer function can be replaced by based on Equation (29). The portion containing only the parametric uncertainty is denoted by , and the time delay uncertainty is denoted by .
In Equation (34), each contains three uncertain variables, which are , and . In the case of , we define as
(35) 
where is the error value of . In the case of , determining the actual value of is difficult compared to other physical quantities. We also define in the same manner as for the convenience of analysis as
(36) 
where are the nominal and error values of . Because the term containing is of an irrational form that is not suitable for analysis, we use an analytic approximation of the uncertain timedelay to a rational function with unmodeled dynamic uncertainty [9]. First, we change the representation of the model to a multiplicative uncertainty form that combines parametric uncertainties and unmodeled timedelay uncertainty as follows:
(37) 
A complex unstructured uncertainty corresponds to unknown time delay , and is the maximum uncertainty that can be caused by . Here, we can obtain using Equation (34) as
(38) 
The maximum value of for each can be found using Euler’s formula as
(39) 
where
(40) 
As a result of analyzing a large amount of actual experimental data, we confirmed that the time delay between () and does not exceed 0.1 second in all three channels. We put 20 percent margin so that . Fig. 6 is multiple Bode magnitude plots of generated by varying from to . From Fig. 6, we can extract
(41) 
for all , which is the upper boundary of sets represented by the red solid line.
IvB determination through analysis
IvB1 robust stability analysis
In [6], the structured singular value is defined as
(44) 
where is a complex structured blockdiagonal unmodeled uncertainty block which gathers all model uncertainties [26]. Following the common notation, the symbol represents a set of all stable transfer matrices with the same structure (full, blockdiagonal, or scalar blocks) and nature (real or complex) as . The is the maximum singular value of uncertainty block . The matrices and are defined by collapsing the simplified overall system to upper LFT uncertainty description as
(45) 
where is the known part of the system, is a reference input and is an output of the overall system. In the theory of the analysis, it is wellknown that the system is robustly stable if satisfies the following conditions
(46) 
The analysis is performed separately for each channel of , , thanks to the structure of the platform described by Equation (28), but since and 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 (45) by using MATLAB’s Robust Control Toolbox, where and in our case. As a reminder, subscript refers to each channel of , , and . Also, structured uncertainty is constructed as
(47) 
which includes unmodeled MOI uncertainty, time and gain uncertainty in our system.
IvB2 Results of analysis
Name  Value  Name  Value 

3  Mass  3.24  
1  0.82  
Limit  1.49  
0.707 
Table I shows the multirotor’s physical quantities and controller gains used both in the simulation and the experiment. The gains and are predefined values set during the primary gaintuning 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 to have a roll and pitch limit of approximately in level flight condition. As previously mentioned, the unmodeled time delay 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 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 analysis. From the analysis, we can see that the system is stable when and .
Fig. 9 shows the results of the SGTbased stability analysis, performed in the same manner as [13]. The analysis is based on the following model:
(48) 
where all uncertainties due to , and are lumped using the functions and , whose magnitude increases over frequency as shown in blue curves of Fig. 9. The stability condition of the SGTbased analysis in this case is
(49) 
[13, 7, 22]. In the SGTbased analysis, the bode plots of the filter with and indicate that system with those values could be unstable. However, through the analysis, those values are still in the stable region. From this, we can confirm that the analysis provides more rigorous boundary values than SGTbased analysis.
V Simulation and Experimental Result
This section reports simulation and experimental results to validate the performance of our threedimensional 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

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 filter are set to and in both simulation and actual experiment with additional margins to ensure additional stability.
Va Validation of acceleration tracking performance
In the experiment, arbitrary desired acceleration commands for and directions are given by the operatorcontrolled radio controller. Fig. 10 shows the multirotor accurately following the target acceleration. From this result, we can confirm that our threedimensional translational acceleration control technique functions effectively even in the actual flight.
VB Validation of DOB performance
VB1 Simulation result
In the simulation, the multirotor follows a circular trajectory with radius of 3 and height of 5 . Meanwhile, the multirotor is exposed to periodic disturbances with accelerations up to 5.5 in each axis. Fig. 11 compares the multirotor’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.
VB2 Experimental Result
In the experiment, the multirotor is commanded to hover at a specific point in threedimensional 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 axis and pulls and releases the force sensor periodically to apply a disturbance to the multirotor.
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 directional position shift than other axes. Unlike the DOBoff case, the DOBon 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 dashsingle dotted line). This is because the wellbehaved DOB algorithm generated control input including the disturbance compensation signal (orange dashsingle 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 and (fifth row), and the comparison between measured by force sensor and estimated by DOB algorithm (sixth row). When DOB is not applied, estimation process is working internally but the signal is not merged into signal, making and have the same value. On the other hand, we can see the difference between the and the signal when DOB is applied, because the signal is merged into the 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 multirotor is set on the centerline of a strong wind generator that generates wind speed of 7 . The performance of DOB is visualized by comparing the position difference between DOBon and DOBoff situations. the multirotor has a position error of about 1 in the DOBoff case and about 0.3 in the DOBon 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 multirotor 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 DOBbased 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 filter of the DOB is designed based on the 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 dronebased industrial operations where precise acceleration control is required. For example, in a multirotorbased 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 multirotor 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 humanrescue missions. For industrial applications involving collaborative flight of multiple multirotors, the proposed algorithm can be used to estimate and stabilize internal forces caused in between physicallycoupled multirotors.
References
 [1] (2015) Time domain disturbance observer based control of a quadrotor unmanned aerial vehicle. In XXV International Conference on Information, Communication and Automation Technologies (ICAT), pp. 1–6. Cited by: Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [2] (2008) Quadrotor dynamics and control rev 0.1. Cited by: §I, Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [3] (2014) Robust trajectory tracking for underactuated vtol aerial vehicles: extended for adaptive disturbance compensation. IFAC Proceedings Volumes 47 (3), pp. 3184–3189. Cited by: Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [4] (2016) Comparison of various quaternionbased control methods applied to quadrotor with disturbance observer and position estimator. Robotics and Autonomous Systems 79, pp. 87–98. Cited by: Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [5] (2014) Highperformance trajectory tracking control of a quadrotor with disturbance observer. Sensors and Actuators A: Physical 211, pp. 67–77. Cited by: §I, Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [6] (1982) Analysis of feedback systems with structured uncertainties. In IEE Proceedings DControl Theory and Applications, Vol. 129, pp. 242–250. Cited by: §IVB1.
 [7] (2002) Advanced techniques for clearance of flight control laws. Vol. 283, Springer Science & Business Media. Cited by: §IVB1, §IVB2, §IV.
 [8] (1999) Disturbance observer and feedforward design for a highspeed directdrive positioning table. IEEE Transactions on control systems Technology 7 (5), pp. 513–526. Cited by: §IIIB3.
 [9] (2018) A robust control approach for hydraulic excavators using synthesis. International Journal of Control, Automation and Systems 16 (4), pp. 1615–1628. Cited by: §IVA, §IV.
 [10] (2018) Robust control of an equipmentadded multirotor using disturbance observer. IEEE Transactions on Control Systems Technology 26 (4). Cited by: §III, Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [11] (2017) Trajectory tracking control of multirotors from modelling to experiments: a survey. International Journal of Control, Automation and Systems 15 (1), pp. 281–292. Cited by: §I.
 [12] (2014) Nonlinear disturbance observer based robust attitude tracking controller for quadrotor uavs. International Journal of Control, Automation and Systems 12 (6), pp. 1266–1275. Cited by: Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [13] (2016) Robust acceleration control of a hexarotor uav with a disturbance observer. In IEEE 55th Conference on Decision and Control (CDC), pp. 4166–4171. Cited by: §IVB2, §IV, Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [14] (2016) An lmi approach to adaptive robust tracker design for uncertain nonlinear systems with timedelays and input nonlinearities. Nonlinear Dynamics 85 (3), pp. 1965–1978. Cited by: Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [15] (2018) Composite nonlinear feedback integral sliding mode tracker design for uncertain switched systems with input saturation. Communications in Nonlinear Science and Numerical Simulation 65, pp. 173–184. Cited by: Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [16] (2018) Adaptive sliding mode control for finitetime stability of quadrotor uavs with parametric uncertainties. ISA transactions 72, pp. 1–14. Cited by: Robust Translational Force Control of MultiRotor UAV for Precise Acceleration Tracking*.
 [17] (1996) Motion control for advanced mechatronics. IEEE/ASME transactions on mechatronics 1 (1), pp. 56–67. Cite
Comments
There are no comments yet.