A Douyere

J-D Lan

Sun Luk

Modelling and Simulation of a Complete System of Energy Transmission

This article presents an equivalent circuit approach used for the modelling of a system consisting of an antenna and a RF/DC rectifying circuit (rectenna) at 2.45GHz. The measurements and simulations steps used for parameter extraction and model validation are presented. The major interest in this approach is the ability to study a complete system of energy transmission in a single simulator environment. This method is illustrated in the case of a microwave transmission without any modulation, that is to say the derived model can be easily integrated. The results obtained using this simple example show a relatively good agreement with the experiment and can form the basis for further developments of more complex systems.

Introduction

Wireless Power Transmission by microwaves at 2.45GHz can be used for remotely supplying power to very low power smart sensors or to transmit energy from Space Power Satellites to the Earth, [START_REF] Sasaki | Guidelines for Development of SPS 2000 Educational Model[END_REF] thus providing clean energy supply at a much wider scale. When transmitting either telecommunication signals or a cw RF beam from a transmitter to a receiver, random variations of atmospheric conditions may occur together with environmental or geometrical changes around the communicating elements. We are currently developing means to alleviate the negative effects that these random perturbations can have on the rectenna energy conversion efficiency. To design these solutions one needs to have a model of the non linear Schottky diode rectifier together with the feeding antenna. This model needs to be sufficiently precise while at the same time being simple to implement in conventional RF circuit design tools. Consequently, we have chosen to derive an equivalent circuit model of the antenna and the transmission channel following the previous work of Gerrits et al., [START_REF] Gerrits John | Modelling and simulation of a dipole antenna for UWB applications using equivalent SPICE circuit[END_REF]. The model parameters extracted are based on measurements made on the half wavelength dipole antenna used in the rectenna. The antenna+rectifier system can subsequently be simulated and optimized for performance using a harmonic balance simulator. The first and the second parts of this article present the equivalent circuit used to model the antenna, the transmission channel and the parameter extraction process. Simulation of the derived rectenna model is presented in the subsequent parts where the obtained results are compared with experiments and show a relatively good agreement. It is believed that this approach can provide a simple and yet efficient simulation framework to further study non stationary perturbation damping circuits.

Modelling of a Dipole Antenna

A half wavelength type dipole antenna with a 2.45GHz resonating frequency will be modelled. This antenna is simple to build and simulate while allowing the fabrication of large collecting arrays which are transparent to light. The latter property is of particular concern and becomes a quality when it turns to collecting large amount of power coming from Space Solar Plants at low power densities. Indeed, large rectenna arrays are needed and must not sacrifice land usage and biological life. Figure 1 shows the equivalent circuit used to model the electrical behaviour of the antenna in response to an incoming RF beam. Port 1 represents the RF beam incident power source. R C , R r , L 1 et C 1 are electrical elements whose value influences the resonating frequency and quality factor Q. Elements TLINE and C2 are included here to allow for the adjustment of the equivalent electrical length. Power dissipated in R C corresponds to the conduction losses whereas resistance R r corresponds to the radiating losses. This model will subsequently be studied in a frequency band from 1GHz to 3GHz encompassing the resonating frequency of the dipole used and corresponding to the measuring equipment at hand at the time of writing this paper.

Equivalent Circuit Parameter Estimation

A 'Gradient' method used in the APLAC software, [3] aims at minimizing a function representing the error from the simulated S 11 parameter of the dipole to the measured one in the chosen frequency band. The S 11 parameter measurement is performed using a vector network analyzer. The reference plane is determined from the calibration using 3.5mm standards. On figure 2 Because of this relatively good agreement the equivalent circuit is expected to provide a correct model of the circuit behaviour in the frequency band of interest. Moreover, with the reciprocity principle for an antenna, it should be possible to use the derived antenna model for both the reception and transmission of RF signals.

RF Transmission Model Validation

In order to validate the model of this antenna that we obtained in simulation, an experimental test setup was developed to measure the energy transmission between two identical λ/2 dipoles at various distances. Similarly, simulations were performed where an equivalent model circuit of the transmitting and receiving antenna was used, [START_REF] Gerrits John | Modelling and simulation of a dipole antenna for UWB applications using equivalent SPICE circuit[END_REF]. Path loss is accounted for by using a controlled source E 1 . The gain of the source is frequency independent but is a function of the distance x :

-0.004* 0.133

ain G x = +
. The propagation delay is accounted for by a transmission line of appropriate length. Figure 4 and 5 show the comparison between the simulated S21-parameter and the measured one.

Microwave Power Transmission System Modelling

In this part, we try to demonstrate the ability of this modelling procedure to fully describe the antennas and the transmission channel behaviour. The design and test of a rectenna (rectifier + antenna) is presented thereafter (figure 6). The purpose of this study is to experimentally estimate the DC energy collected by a constant resistive load R load =73Ω connected at the output of the rectenna, as a function of the transmitting distance and for a given transmitted power P=26dBm and to compare this data with the simulations of the same system using the equivalent circuit model approach developed in this paper. For practical reasons, this study was realized for distances ranging from 10cm to 30cm.

Figure 7 shows the measured and simulated results obtained during this study. It can be noticed that a very good agreement is obtained. The observed offset stays within the margin of error due to fabrication tolerances and measurement errors. Dependence with the distance is very well predicted by the simulator and this shows that the behaviour of the rectifier seems to be correctly taken into account by the equivalent circuit model used in this work.

Conclusion

We have presented an equivalent circuit model of a complete Wireless Power Transmission system. Models of each element of the system can be identified using experimental measurement results. The major interest in this approach is the ability to fully optimize the circuit in a single simulator environment. For example, using this equivalent circuit model, the complete system of energy transmission could be optimized for low spurious emissions and high rectenna efficiency. In a future work, impedance of the equivalent circuit representing the antenna may be optimized and the approach will require us to synthesize the antenna exhibiting the same model. This may require the modelling of more complex antenna structures.

Fig. 1 :

 1 Fig. 1: Equivalent electrical circuit of a 2.45GHz resonant frequency half wavelength dipole antenna.

Fig. 2 :

 2 Fig. 2: Simulated and measured S 11 magnitude curves compared after optimization.

Fig. 3 :

 3 Fig. 3: Smith chart representation of optimized and measured S 11 parameter.

Fig. 4 :

 4 Fig. 4: Comparison between the simulated S21parameter and the measured one.

Fig. 5 :

 5 Fig. 5: Comparison between the simulated S21- parameter and the measured one.

Fig. 6 :

 6 Fig. 6: Structure of the Rectenna. Fig. 7: DC voltage at rectenna output as a function of transmission distance.


Gabriele Perozzi

A toolbox for quadrotors: from aerodynamic science to control theory

In this paper the quadrotor unmanned aerial vehicle (UAV) is studied. Aerodynamic forces and moments are explained together with aerodynamic coefficients and a complete quadrotor model is illustrated. After that, feedback control, trajectory generation and proportional integrative derivative (PID) control algorithms are applied. Notions of controllability, observability and area of attraction are also introduced. It is always very hard to conciliate aerodynamic science with control theory because of different hypothesis. Thus, the main focus of this paper is to group them to create a basic guide that lets the reader to understand how quadrotors work and how to build a first linear control effectively, covering aerodynamics and control background together. Important remarks are done to the simplifications of aerodynamic coefficients and flight dynamic equations useful for control theory. Two quadrotor Simulink models are provided together with the article. One is built using the highly nonlinear aerodynamic coefficients, the other with some simplifying hypothesis. Both of them use only blocks provided with basic Simulink environment to ensure a better compatibility with most of Simulink versions, avoiding integrated matlab functions and any other auxiliary toolbox.

Introduction

The presented work is based on quadrotors having the configuration as in Fig. 1. Quadrotors are underactuated systems with 6 degrees of freedom and 4 controls. The paper outline is as follows. In Section 2 the inertial and body frames are introduced, the quadrotor flight dynamics and the aerodynamic coefficients are derived. In Section 3 some admissible simplifications are allowed to make the system easier to study, the control-oriented model is illustrated, notions of controllability, observability and area of attraction are introduced, then linear feedback state control, PID control and trajectory generation are explained for a quadrotor respecting its physical constraints. The Section 4 is dedicated to the explanation of the Simulink quadrotor models, provided together with this paper. Appendices A, B, C conclude the paper.

2 Quadrotor model

Body and inertial reference systems

A quadrotor is a vehicle that can fly in the airspace. Hence to be properly studied, two different frames must be considered. The body frame of the quadrotor, and an inertial frame which can be represented by the earth. Forces acting on the quadrotor are computed in body frame, and than transposed in inertial frame (see Fig. 2) thanks to the rotation matrix R. R is the final rotational matrix defined as

R = R z (ψ) R y (θ) R x (φ) =   c ψ c θ -s ψ c φ + c ψ s θ s φ s φ s ψ + c ψ s θ c φ s ψ c θ c ψ c φ + s ψ s θ s φ -c ψ s φ + s ψ s θ c φ -s θ c θ s φ c θ c φ   ,
where c ψ = cos(ψ), s ψ = sin(ψ) and similarly. Then the passage from the earth inertial frame (R i ) to the body frame (R b ) is

[X T ] R b = [X T ] Ri R.
(1)

Flight dynamics

For the Parrot AR Drone 2.0 or similar shape, the rotors gyroscopic effects and inertial counter torques can be neglected since they are very small. Then the translation dynamics of the drone in the body frame yield

m u + m × u = F aero + mR T g, ( 2
)
where m is the mass of the UAV, u = [u v w] T is its linear velocity expressed in body frame, = [p q r] T is its angular velocity in body frame, F aero = [F Xaero F Y aero F Zaero ] T is the vector of the external aerodynamic forces in body frame, g = [0 0 g] T is the gravity acceleration in inertial frame, R is the rotational matrix. The rotational dynamics of the drone with respect to inertial earth frame are

I ˙ = -× I + L aero , (3)
where I is the inertia matrix defined as: the UAV, and τ aero = [L aero M aero N aero ] T is the external aerodynamic moments in the body frame. The relation between angular velocities and Euler angles φ = p + tan θ(q sin φ + r cos φ), θ = q cos φ -r sin φ,

ψ = q sin φ + r cos φ cos θ , (4)
out the singularity at θ = π 2 , which is never reached during flight (it is a reasonable assumption if we don't want to reach aggressive maneuvers). Hence the full model is given by the equations ( 2), (3), (4). The aerodynamic forces F Xaero , F Y aero , F Zaero , moments L aero , M aero , N aero , and related coefficients are derived below using a combination of momentum and blade element theory in helicopters. Aerodynamic forces and moments for each rotor, where subscript j indicates the j th rotor, are derived as

F X j = -ρAR 2 u j -u w (u j -u w ) 2 + (v j -v w ) 2 C H j ω 2 j , F Y j = -ρAR 2 v j -v w (u j -u w ) 2 + (v j -v w ) 2 C H j ω 2 j , F Z j = -ρAR 2 C T j ω 2 j , L j = -sign ω j ρAR 3 u j -u w (u j -u w ) 2 + (v j -v w ) 2 C Rmj ω 2 j , M j = -sign ω j ρAR 3 v j -v w (u j -u w ) 2 + (v j -v w ) 2 C Rmj ω 2 j , N j = -sign ω j ρAR 3 C Q j ω 2 j ,
where ρ is the air density, A is the rotor area, R is the rotor radius, u w , v w , w w are the wind velocities with respect to the earth in body frame, C H j is the hub force coefficient, C T j is the rotor thrust coefficient, ω j is the rotor angular velocity, C Q j is the rotor drag moment coefficient, C Rmj is the rotor rolling moment coefficient. Total aerodynamic forces are

F Xaero = 4 j=1 F X j , F Y aero = 4 j=1 F Y j , F Zaero = 4 j=1 F Z j . (5)
Total aerodynamic moments are

L aero = 4 j=1 (L j + F Z j ls j -hF Y j ), M aero = 4 j=1 (M j -F Z j lc j + hF X j ), N aero = 4 j=1 (N j + F Y j lc j -F X j ls j ), ( 6
)
where h is the distance between rotors plane and the center of gravity of the UAV, l is the arm length, and with

c j = cos π 2 (j -1) + , s j = sin π 2 (j -1) + ,
where in our UAV configuration we have = π 4 . Thus, for vectors c j and s j we have cosines and sinus of the angles

[ π 4 , 3 4 π, 5 4 π, 7 4 π].
Linear rotors speed in body frame are computed as a function of the state in body frame

  u j v j w j   =   p q r   ×   lc j ls j h   +   u v w   .

Aerodynamic coefficients

The blade flapping is modeled as a fictitious hinge spring as in Fig. 3. With the hypothesis that β is small, the equation is

I x β + ω 2 I x + M blade x g eR 2 β = M a -k β β,
where, ω is the rotor angular velocity, I x is the inertia of the blade, x g R is the position of the CoG of the blade, M blade is the mass of the blade, and M a is the aerodynamic momentum of each element of the blade. Then the blade flapping equation can be rewritten as

d 2 β dψ 2 + 1 8 γ(1 + 4 3 µ sin ψ) dβ dψ + (1 + δ) + 1 8 γ( 4 3 µ cos ψ + µ 2 sin 2ψ) β = - 1 8 γ θ 0 (1 + 8 3 µ sin ψ + 2µ 2 sin 2 ψ) -θ tw ( 4 5 + 2µ sin ψ + 4 3 µ 2 sin 2 ψ) - 4 3 λ -2µλ sin ψ ,
where ψ is the azimuth angle of the blade. If the transition phase is stable and the blade movement is periodic in ψ the equation becomes

β = a 0 + a 1 sin ψ + b 1 cos ψ, (7)
with

a 0 = - γ 8(1 + δ) θ 0 (1 + µ 2 ) -θ tw 4 5 + 2 3 µ 2 - 4 3 λ , δa 1 = - γ 4 µ 4 3 θ 0 -θ tw -λ + γ 8 b 1 1 - µ 2 2 , δb 1 = - γ 6 µa 0 - γ 8 a 1 1 + µ 2 2 , δ = M blade x g eR 2 I x , γ = ρacR 4 I x ,
where c is the average length of the blade chord.

The advance ratio indicates the component of the UAV velocity parallel to the rotor disk with respect to the blade tip speed:

µ j = (u j -u w ) 2 + (v j -v w ) 2 R|ω j | .
The inflow ratio is the ratio between the component of UAV velocity perpendicular to the rotor disk with respect to the blade tip speed:

λ j = σa (1 + 3 2 µ 2 j ) θ0 6 -(1 + µ 2 j ) θtw 8 - λj 4 2 µ 2 j + λ 2 j + w w -w j + qlc j -pls j R|ω j | ,
where σ is the rotor solidity ratio, a is the lift curve slope of the blade section, θ 0 is the angle of attack of the root profile, θ tw is the twist angle of the blades.

The rolling moment of a propeller exists in forward flight when the advancing blade is producing more lift than the retreating one and it is the integration over the entire rotor of the lift of each section acting at a given radius. Its coefficient is:

C Rmj = σa µ j 8 (λ j - 4 3 θ 0 + θ tw ) + b 1 16 (1 - µ 2 j 2 ) .
The thrust is the resultant of the vertical forces acting on all the blade elements. Its coefficient is:

C T j = σa (1 + 3 2 µ 2 j ) θ 0 6 -(1 + µ 2 j ) θ tw 8 - λ j 4 .
The hub forces is the resultant of the horizontal forces acting on all the blade elements. Its coefficient is:

C H j = C H P j + C H ij , C H P j σ = µ j 4 (C D0 + C Di θ 2 0 ) + C Di µ j 24 (3θ 2 tw -8θ 0 θ tw ) + θ 0 24 (3µ 2 j b 1 -12λ j µ j -4b 1 ) - θ tw 16 (µ 2 j b 1 -4λ j µ j -2b 1 ) + C Di µ 2 j 8 a 0 a 1 + µ j 16 (a 2 1 -b 2 1 ) + 1 4 λ j b 1 , C H ij σa = θ 0 4 (λ j µ j + 2 3 b 1 ) - θ tw 8 (λ j µ j + b 1 ) + µ j 8 (a 2 0 + b 2 1 ) - 3 8 λ j b 1 + 1 12 a 0 a 1 ,
where C D0 is the drag coefficient of the blade section, C Di is the induced drag coefficient of the blade section.

The drag moment about the rotor shaft is caused by the aerodynamic forces acting on the blade elements, the horizontal forces acting on the rotor are multiplied by the moment arm and integrated over the rotor. Its coefficient is:

C Q j = C Q P j + C Q ij , C Q P j σ = 1 8 (C D0 + C Di θ 2 0 )(1 + µ 2 j ) -C Di θ 0 θ tw ( 1 5 + µ 2 j 6 ) + C Di θ 2 tw ( 1 12 + µ 2 j 16 ) -C Di λ j × ( θ 0 3 - θ tw 4 ) + C Di µ 2 j 8 (a 2 0 + a 2 1 4 + 3b 2 1 4 ) + 1 16 (a 2 1 + b 2 1 ) + λ 2 j 4 + µ j 6 a 0 a 1 - λ j µ j 4 b 1 C Q ij σa = λ j ( θ 0 6 - θ tw 8 - λ j 4 ) - µ 2 j 8 (a 2 0 + a 2 1 4 + 3b 2 1 4 ) - 1 16 (a 2 1 + b 2 1 ) - µ j 6 a 0 a 1 + λ j µ j 4 b 1 ,

Rotors dynamics

The rotors are driven by DC-motors, which binds electrical and mechanical quantities. Rotors can be represented as in Fig. 4, and described by the dynamic system

L di dt = u -Ri -K e ω m J dωm dt = τ m -τ d ,
where u is the motor input, R is the motor inertial resistance, L is the motor inductance, K e is the motor constant corresponding to the back EMF constant, τ m is the motor torque, ω m is the motor angular speed, τ d is the motor load, J is the motor moment of inertia, i is the current. Considering a small motor with a very low inductance

Figure 4: Rotor model the given system can be simplified and approximated, then linearized around an operating point ω 0 , and further identified as:

G(s) = 1 bs + 1 ,
where the time constant b = 0.1 can be found considering Parrot Drones. Time constant slightly varies based on the considered motor and thus on the considered UAV mass. Precisely, in the Parrot Drones, the time constant is less than 0.1 for the increasing rotor angular velocities, and it is almost equal to 0.1 for decreasing rotor angular velocities.

3 Trajectory and Control

Model simplification

The presented quadrotor model is highly nonlinear and thus it is not compatible to be studied in control theory. A trade-off should be found between the aerodynamic science (which allows highly nonlinearities to better represent the real quadrotor behaviour) and the control theory (which allows the system to be as much simple as possible).

Identification results thanks to indoor experiments on Parrot drone at low/moderate velocity, validates the previous UAV model and allows some simplifications to be accepted:

• θ tw = 0;

• C Di = 0;

• λ j = λ stat -4 σa K z wj -ww R|ωj | , K z , λ stat ≥ 0; • C T j = C T stat + K z wj -ww
R|ωj | , C T stat ≥ 0, such models of λ j and C T j are rather precise in the climbing phase, but less accurate in descent phase, since the model tends to slightly overestimate the propulsion in the descent phase;

• C H j = K D µ j , K D ≥ 0;
the UAV drag is modeled as ρAR 2 C H j ω 2 j , corresponding to the rotors at low UAV speed, otherwise at higher speed we need to add the body drag effect since it depends on square of velocity. However the constant K D is identified with tests considering the whole UAV body and rotors;

• Blade flapping is neglected.

Control-oriented model

The system can be rewritten in the state-space form

Ẋ = f (X, U, d) (8)
where f is expressed in (2), (3), (4) with the simplifying hypothesys. The state vector X is chosen as

X = [x y z ẋ ẏ ż φ θ ψ p q r] T , the control input is U = [U z U θ U φ U ψ ] T ,
where its relation with the rotor velocities is defined by an invertible matrix

    U z U θ U φ U ψ     =     K f K f K f K f K f lc j K f lc j K f lc j K f lc j -K f ls j -K f ls j -K f ls j -K f ls j K m -K m K m -K m         ω 2 1 ω 2 2 ω 2 3 ω 2 4     , (9)
with ω min ≤ ω j ≤ ω max , where

K f = ρAR 2 C T stat and K m = ρAR 3 σC D0 8 + λ stat σa( θ0 6 -λstat 4
) . In this work, the control inputs are selected to be proportional to the terms with ω 2 j . Thus, expanding ( 2) and ( 3), the other terms dependent linearly on ω j and wind velocities are considered as disturbances. Since we do not know in advance the wind perturbations, then we cannot use these terms in controls. Such a decomposition of thrust (which is proportional to ω j and ω 2 j ) and selection of disturbances are almost exact in the hover flight, where we have (p, q, r) ≈ ( φ, θ, ψ). Linear rotors velocities in body frame are computed as a function of the state in body frame

  u j v j w j   =   p q r   ×   lc j ls j h   +   u v w   .
We will call signals, whose influence is compensated by robust abilities of the control, as "disturbances", which include external inputs (the wind velocities), state components and controls. We have the disturbance d

= [d x , d y , d z , d φ , d θ , d ψ ]
, where d x , d y , d z are expressed in body frame, explicited as

d x = 4 j=1 -ρARK D (u j -u w ) |ω j |, d y = 4 j=1 -ρARK D (v j -v w ) |ω j |, d z = 4 j=1 -ρARK z (w j -w w ) |ω j |, d φ = 4 j=1 ω j ρAR 2 (u j -u w ) σa 2 ( θ 0 3 - λ stat 4 ) + |ω j |ρAR hK D (v j -v w ) -lK z (w j -w w ) s j + sign (ω j ) 1 2 ρARK z (u j -u w ) (w j -w w ) , d θ = 4 j=1 ω j ρAR 2 (v j -v w ) σa 2 ( θ 0 3 - λ stat 4 ) + |ω j |ρAR -hK D (u j -u w ) + lK z (w j -w w ) c j + sign (ω j ) 1 2 ρARK z (v j -v w ) (w j -w w ) , d ψ = 4 j=1 ω j ρAR 2 K z (w j -w w ) ( 2θ 0 3 -2λ stat ) -|ω j |ρlARK D (v j -v w ) c j -(u j -u w ) s j -sign (ω j ) ρAR σC D0 8 (u j -u w ) 2 + (v j -v w ) 2 - 4 σa K 2 z (w j -w w ) 2 .
In earth frame, using (1), the linear dynamics disturbances become

d xe =c ψ c θ d x + (c ψ s θ s φ -s ψ c φ )d y + (s φ s ψ + c ψ s θ c φ )d z d ye =s ψ c θ d x + (s ψ s θ s φ + c ψ c φ )d y + (s ψ s θ c φ -c ψ s φ )d z d ze = -s θ d x + c θ s φ d y + c θ c φ d z
Under the hypothesis of quasi-hover the system is explicited in compact form as

                                                 ẋ1 = x 4 ẋ4 = (s φ s ψ + c ψ s θ c φ ) 1 m (U z -d xe ) ẋ2 = x 5 ẋ5 = (s ψ s θ c φ -c ψ s φ ) 1 m (U z -d ye ) ẋ3 = x 6 ẋ6 = g -(c φ c θ ) 1 m (U z + d ze ) ẋ7 = x 10 ẋ10 = x 11 x 12 Ixx-Izz Ixx + 1 Ixx (U φ + d φ ) ẋ8 = x 11 ẋ11 = x 10 x 12 Izz-Ixx Iyy + 1 Iyy (U θ + d θ ) ẋ9 = x 12 ẋ12 = x 10 x 11 Ixx-Iyy Izz + 1 Izz (U ψ + d ψ ) (10)

State controllability

The state controllability property of a generic system describes the ability of an external input, such as the control signal, to move the state of the system from an initial point to a given final point in a finite time.

Let's consider the LTI system:

Ẋ = A X + B U Y = C X + D U
Then the system is state controllable if the controllability matrix

R = [B AB A 2 B . . . A n-1 B]
has rank n, where n is the dimension of the matrix A. In matlab, the function rank(ctrb(A,B)), gives the rank of the controllability matrix, which should be equal to the dimension of the matrix A, in our case n = 12. This concept is useful to be sure that the system can be fully controllable with the given control input. It means for example that, to control properly the quadrotor system the full control vector must be considered, to control the attitude sub-system only 3 elements of the control vector can be considered (U φ , U θ , U ψ ), and so on.

Observability

A system is observable if the behavior of the entire system can be determined from the output. It means that the values of the state can be reconstructed using the output values generally in the sensors. Let's consider again the same system for the controllability, then if the row rank of the following observability matrix

O = [C CA CA 2 . . . CA n-1 ] T
is equal to n, then the system is observable. This is a key concept in estimation algorithms. For example, if we want to estimate the wind velocity using the quadrotor as wind sensor, the quadrotor must be observable in earth frame using cameras, since the wind estimation algorithms use information related to the drone velocity with respect to the wind in earth frame. In matlab, the function rank(obsv(A,C)), gives the rank of the observability matrix, which should be equal to the dimension of the matrix A, in our case n = 12.

Area of Attraction

Another key concept in control algorithms, is the Area of Attraction (AoA). The AoA is the area of the state-space where the state of the system, from any initial point, can reach a final point, called equilibrium point. Any point initially placed on the unstable equilibrium point, inside this area, will reach the equillibrium point in finite-time (stability notion) or in infinite time (asymptotical stability notion). The AoA concept is strictly correlated with the control algorithms, because different controls have different AoAs. For example the linear control, such as PID, have small AoA around the chosen equilibrium point. Many equilibrium points can be chosen in the state-space and each of them has a different small local AoA. While the nonlinear controllers, has a larger domain of stability.

For complex controls and with high dimension systems, AoA can be found only empirically in simulations because of the high computational complexity. However, the AoA of linear controls, such as feedback control in Section 3.7 for the system (8), can be found analytically following the steps (see code in Appendix A):

1. Given the wind d w = (u w , v w , w w ) T , find the equilibrium point (X 0 , U 0 ) such that f (X 0 , U 0 , d w ) = 0.

2. Find the matrices A, B.

3. Given the matrices A, B, find the gain K such that A -B K is Hurwitz.

4. Find P = P T > 0 such that (A -BK) T P + P (A -BK) < 0.

5. Find δx such that dV = 2δx T P f (X 0 + δx, U 0 -Kδx, d w ) ≤ 0.

6. δx is the area in 12 dimensions of the state-space, inside this area the system is stable. We finally found the Area of Attraction.

7. Repeat the entire process for another chosen wind d w .

Physical constraints and trajectory generation

Since UAVs are physical objects, then they are subjected to physical constraints, which correspond to their maximal linear and angular velocities and accelerations. Consequently, a smooth position trajectory is desired. In literature many articles are presented to generate desired trajectories admissible to the system dynamics. One method is to use a third order filter to obtain realistic velocity and acceleration trajectories

ξ f ξ wp = 1 (1 + τ s) 3 ,
where ξ wp is the way-point step reference, ξ f is the filtered signal, τ is the arbitrary weight to create the desired smooth trajectory. The reference signal is third order filtered as in Fig. 5, with G 1 = 1/(τ pos τ vel τ acc ), G 2 = τ pos τ vel + τ acc τ vel + τ pos τ acc , G 3 = τ pos + τ vel + τ acc where τ pos , τ vel , τ acc are the three weighted parameters for position, velocity and acceleration respectively. This method allows to consider the limitations, inside the integrals. Results are illustrated in Figures 6, 7, 8. Different gains allow to obtain different trajectories. respectively for linear velocity and acceleration, where T max is 90% of the maximal thrust achievable by the UAV to ensure that we still have some available thrust to proper control the drone in case of safe maneuvers; T min is the minimal thrust. The reference angles φ des , θ des are saturated in an admissible range of φ min , θ min minimal values and φ max , θ max maximal values, to avoid the singularities. Physical constraints for angles and torque controls are imposed as

∆ min =   φ min θ min ψ min   , ∆ max =   φ max θ max ψ max   ; Ω min =   φmin θmin ψmin   , Ω max =   φmax θmax ψmax   ; U φ,θ,ψmin =   -lρAR 2 C T stat sin( )(ω 2 max -ω 2 min ) -lρAR 2 C T stat cos( )(ω 2 max -ω 2 min ) -2ρAR 3 C Q (ω 2 max -ω 2 min )   , U φ,θ,ψmax =   lρAR 2 C T stat sin( )(ω 2 max -ω 2 min ) lρAR 2 C T stat cos( )(ω 2 max -ω 2 min ) 2ρAR 3 C Q (ω 2 max -ω 2 min )   .

Feedback state control

The first simplest control is the stabilization of the drone in different operating points (see Appendix B). To this end, let's consider the system: Supposing the trim point is (X 0 , U 0 , d 0 ) (important remark is that the angles depend by the linear quadrotor velocity and the wind, thus the effect of the wind should be considered while computing the trim point), let denote by δX, δU and δd the deviations from that point of X, U and d, respectively, then:

Ẋ = f (X, U, d),
X = X 0 + δX, U = U 0 + δU, d = d 0 + δd.
Expanding the model equations in a truncated Taylor series of first order the following is obtained

δX = f (X 0 , U 0 , d 0 ) + ∂f ∂X (X 0 , U 0 , d 0 )δX + ∂f ∂U (X 0 , U 0 , d 0 )δU + ∂f ∂d (X 0 , U 0 , d 0 )δd Denote constant matrices: A = ∂f ∂X (X 0 , U 0 , d 0 )δX, B = ∂f ∂U (X 0 , U 0 , d 0 )δU, D = ∂f ∂d (X 0 , U 0 , d 0 )δd, G = F (X 0 , U 0 , d 0 ),
then the linearized equations of the system around the point (X 0 , U 0 , d 0 ) can be rewritten as follows

δ Ẋ = AδX + BδU + Dδd + G.
To stabilize the quadrotor around the chosen trim condition, a simple feedback state control can be applied

U = K 0 X + K 1 ,
where K 1 stabilizes the UAV against the gravity force. Thus, if the system is chosen to be stabilized around stationary case (but any desired operating point can be chosen), the control gain can be computed as

K 1 = -g b1 0 0 0 T ,
where b 1 = ρAR 2 C T stat is the term proportional to the control U z = -mg in ż dynamics for the stationary flight.

To find K 0 the closed-loop system must be analyzed

Ẋ = (A + B K 0 )X + BK 1 + G Y = X
It means to find K 0 such that the eigenvalues are at the desired location in the roots graph. Hence, this linear system is stable if and only if the eigenvalues of A + B K 0 are either real and negative or else complex with non-positive real part, and it is asymptotically stable if the complex eigenvalues of A + B K 0 have strictly negative real part.

Since the system is fully controllable the gain can be tuned using the Matlab function place.

PID control

An enhanced version of the previous feedback state control is the PID. While the feedback state allows only the proportional term, the PID can provide better stability using the integrative and the derivative terms too. The proposed PID is valid for linearized system around hover flights, so under the following hypothesis:

U z = -mg, cos φ = 1, cos θ = 1, sin φ = φ, sin θ = θ, φ = p, θ = q, ψ = r,
p des = q des = r des = 0.

Since the x, y cannot be directly controlled, then auxiliary controls U x , U y are added, which will be used to find the controlled angles. Position controller is structured as

U x =K P x (x des -x) + K Ix (x des -x) + K Dx ( ẋdes -ẋ) + ẍdes , U y =K P y (y des -y) + K Iy (y des -y) + K Dy ( ẏdes -ẏ) + ÿdes , U z =g -K P z (z des -z) + K Iz (z des -z) + K Dz ( żdes -ż) + zdes ,
where K P i , K Ii , K Di for i = x, y, z are the tuning parameters for proportional, integrative and derivative terms. Attitude controller is structured as

U φ =I xx K P φ (φ c -φ) + K Iφ (φ c -φ) -K Dφ p , U θ =I yy K P θ (θ c -θ) + K Iθ (θ c -θ) -K Dθ q , U ψ =I zz K P ψ (ψ des -ψ) + K Iψ (ψ des -ψ) -K Dψ r ,
where K P i , K Ii , K Di for i = φ, θ, ψ are the tuning parameters. The controlled angles φ c , θ c are derived from U x , U y and using the desired value of ψ des

φ c = - 1 g (U x sin ψ des -U y cos ψ des ) , θ c = - 1 g (U x cos ψ des + U y sin ψ des ) .
The generic scheme of hierarchical control algorithm presented above is given in Fig. 9. The attitude control loop usually works 10 times faster than the position control loop. • A nonlinear model which includes all the nonlinear aerodynamic coefficients (called complex);

• A nonlinear model which includes the identification in Section 3.1 (called simple).

The best strategy to work with these models is to use at first the simple one, which is more suitable to be studied for control and estimation techniques, then validate the algorithms with the complex model, checking the differences due to the neglected nonlinear terms. Their peculiarities are that:

• Only basic simulink blocks are used, avoid the necessity of interpreted or external matlab functions, allowing better compatibility with most of simulink versions. In particular the used blocks are: Constant, Mux, Demux, Gain, Integrator, In, Out, Product, Terminator, Switch, Compare To Constant, Relational Operator, Logical Operator, Abs, Add, Divide, Math Function, Product of Elements, Sqrt, Subtract, Sum of Elements, Unary Minus, Unit Delay (it is only used in the complex model), Trigonometric Function, Subsystem, While Iterator Subsystem (it is only used in the complex model), From, Goto, Stop Simulation, To Workspace, Cross vector.

• Auxiliary matlab/simulink toolbox, such as the purchasable Aerospace toolbox provided by Simulink, is not needed except of the basic simulink package;

• They are based on a detailed model ensuring a good quality simulator close to the reality;

• They are fully configurable changing the parameter values, allowing to study any kind of quadrotor having the shape similar to the commercial Parrot AR 2.0 Drone;

• The complex model can recognize the vortex ring state and stop the simulation accordingly. When the variable "test", in matlab workspace, is equal to 1.

Input and Output of the model are defined as (see Figures 12,13):

• Input -Rotors velocity: One vector of angular velocities expressed in rad/s, where the 4 rotors are indicated with their signs according to the velocity direction, ex: +, -, + -;

-Wind: One vector of the wind velocity expressed in earth frame of the 3 components in x, y, z axes in m/s.

• Output -State: One vector of 12 components of the state of the system. In order: 3 linear velocities in earth frame expressed in m/s, 3 angular velocities expressed in rad/s, 3 angles expressed in rad, 3 position expressed in m.

Aerodynamic coefficients, forces, moments and flight dynamics are illustrated in the schemes in Figures 14,15.

The following parameter values can be changed to simulate quadrotors with different shapes and in different environment conditions using the higly nonlinear coefficients (see Fig. 16):

• Environment:

-Gravity acceleration;

-Air density.

• Quadrotor body:

-Mass of the drone;

-Moment of Inertia on x axes;

-Moment of Inertia on y axes;

-Moment of Inertia on x-z axes;

-Angle (see Fig. 1);

-Arm length;

-Distance between the center of gravity of the drone and the rotors plane.

• Rotors:

-Angle of attack of the blade at the root profile;

-Twist of the blade;

-Rotor radius;

-Lift curve slope of the blade section;

-Blades number;

main chord length of the blades;

-Drag coefficient of the blade section;

-Induced drag coefficient of the blade section.

• Blade flapping equation (see ( 7)):

-Position of the spring (e);

-Center of gravity of the blade;

-Mass of the blade;

-Inertia of the blade.

• Auxiliary parameter to solve the λ equation:

-Maximum number of allowed iterations to solve the nonlinear λ equation, after that the system is in unstable position and automatically blocks the simulation: It means we are in vortex ring state.

• UAV initial conditions:

-Initial linear velocity in earth frame;

-Initial position in earth frame;

-Initial angular velocity;

-Initial angles.

The following parameter values can be changed to simulate quadrotors with different shapes and in different environment conditions using the simplifying hypothesis based on identification in Section 3.1 (see Fig. 17):

• Environment:

-Gravity acceleration;

-Air density.

• Quadrotor body:

-Mass of the drone;

-Moment of Inertia on x axes;

-Moment of Inertia on y axes;

-Angle (see Fig. 1);

-Arm length;

-Distance between the center of gravity of the drone and the rotors plane.

• Rotors:

-Angle of attack of the blade at the root profile;

-Rotor radius;

-Lift curve slope of the blade section;

-Blades number;

-Main chord length of the blades;

-Rotors drag coefficient (K D );

-Rotors thrust coefficient (K z );

-Rotors thrust coefficient in stationary phase (C T stat ).

• UAV initial conditions:

-Initial linear velocity in earth frame;

-Initial position in earth frame;

-Initial angular velocity;

-Initial angles.

An auxiliary rotor model is also provided with a configurable mask, which allows to transform the control input to the rotors velocity taking into account the rotor transfer function. Simulink blocks used are: Constant, Sqrt, Product, Divide, Abs, Gain, Discrete-Time Integrator, In, Out, Sum. Input and Output of the model are defined as (see Fig. 18):

• Input:

control: vector of 4 components (U z , U φ , U θ , U ψ )

• Output:

omega: vector of 4 components. Rotors angular velocity in rad/s.

The following parameter values can be changed to simulate rotors with different shapes and in different environment conditions using the mixer matrix (9) (see Fig. 19):

• Rotor transfer function:

-Denominator;

-Numerator;

-Rotor angular velocity in hover phase.

• Mixer Matrix:

-Air density;

-Rotor radius;

-Rotor thrust coefficient in hover flight (C T stat );

-Drag coefficient of the blade section;

-Angle of attack at the root profile of the blade section;

-Rotors position ( , see Fig. 1);

-Blades number;

-Main chord length of the blade section;

-Drone arm length.

A Area of Attraction in matlab %Area of Attraction syms omega1_2 omega2_2 omega3_2 omega4_2 U v w p q r phi theta psi t1 t2 t3 t4 t5 t6 ∆x P = sdpvar(9,9); matrix=-double(vpa(subs(de2bi(0:511),0,-1)))'; 7))+x(6) * cos(x(7))); f8=x(5) * cos(x(7))-x(6) * sin(x( 7)); f9=(x(5) * sin(x( 7))+x(6) * cos(x(7)))/cos(x(8)); FU=[f1;f2;f3;f4;f5;f6;f7;f8;f9];

∆ =
%-------------------------------------------------------------------------------------- %wind effect
t0=[0 0 363.666^2 363.666^2 363.666^2 363.666^2]; t = fsolve(@(tt) fcn_change(tt,FU,t1,t2,t3,t4,t5,t6),t0); u=[omega1_2;omega2_2;omega3_2;omega4_2];

%System equations--------------------------------------------------------------------------------------
x=[U,v,w,p,q,r,phi,theta,psi]; %NonLinear System F(x0+∆x,u0-K∆x,d) %System equations-------------------------------------------------------------------------------------- %--------------------------------------------------------------------------------------%linear system matrices AM = jacobian ([f1, f2, f3, f4, f5, f6, f7, f8, f9], [x(1), x(2), x(3), x(4), x(5), x(6), x(7), ... x(8), x(9)]); BM = jacobian ([f1, f2, f3, f4, f5, f6, f7, f8, f9], [u(1), u(2), u(3), u(4)]); %variables to evaluate jacobian U0=x0(1);v0=x0(2);w0=x0(3);p0=x0(4);q0=x0(5);r0=x0(6);phi0=t(1);theta0=t(2);psi0=x0(9); ,t1,t2,t3,t4,t5,t6) F=double (vpa(subs(F1,[t1,t2,t3,t4,t5,t6],[tt(1),tt(2),tt(3),tt(4),tt(5),tt(6)]))); end

equations-------------------------------------------------------------------------------------

Figure 1 :

 1 Figure 1: Quadrotor configuration and reference systems.

Figure 2 :

 2 Figure 2: Quadrotor configuration and reference systems.

I

  xx I xy I xz I yx I yy I yz I zx I zy I zz   ,

Figure 3 :

 3 Figure 3: Fictious spring modeling the blade flapping

Figure 5 :Figure

 5 Figure 5: Third order filter scheme for trajectory

Figure 8 :

 8 Figure 8: Quadrotor filtered acceleration

Figure

  Figure 9: Hierachical control

  %change only the first three components if we are in different operating point wind=[1;2;3]; %change wind value to check the AoA of a wind range: 3 cycles for uu=[t3; t4; t5; t6]; xx=[x0(1); x0(2); x0(3); x0(4); x0(5); x0(6); t1; t2; x0(9)]; %System equations----------------------------------------------------------------------------------x(3)-wind(3))./(R. * abs(sqrt(u))); %inflow ratio lambda=4. * (-CT./(sigma * a)+theta_0/6); %coeff rolling moment CRm=sigma * a. * ((mu./8). * (lambda-(4/3) * theta_0)); %coeff torque CQi=sigma. * a. * lambda. * (theta_0/6-lambda./4); CQp=sigma./8. * (CD_0+CD_i * theta_0^2). * (1+mu. * mu)-CD_i. * (lambda. * theta_0/3-lambda. * lambda./4(epsilon); sin(pi/2+epsilon); sin(pi+epsilon); sin(3 * pi/2+epsilon)]; coseno=[cos(epsilon); cos(pi/2+epsilon); cos(pi+epsilon); cos(3 * pi/2+epsilon)-1;1;-1]. * rho. * A. * R^3. * CQ. * u; else Fxj=-rho. * A. * R^2. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fyj=-rho. * A. * R^2. * ((x(2)-wind(2))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fzj=-rho. * A. * R^2. * u. * CT; Lj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CRm. * u; Mj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(2)-wind(2))/sqrt((x(2)-wind(2))^2+(x(2)-wind(2))^2)). * CRm. * u; Nj=[1;-1;1;-1]. * rho. * A. * R^3. * CQ. * u; end %aeroynamics F=[sum(Fxj);sum(Fyj);sum(Fzj)]; J=[sum(Lj)-h * sum(Fyj)+l * sum(Fzj. * seno);... sum(Mj)+h * sum(Fxj)-l * sum(Fzj. * coseno);... sum(Nj)+l * sum(Fyj. * coseno)-l * sum(Fxj. * seno)Ixx-x(5) * x(6) * (Izz-Iyy)/Ixx; f5=J(2)/Iyy-x(4) * x(6) * (Ixx-Izz)/Iyy; f6=J(3)/Izz-x(4) * x(5) * (Iyy-Ixx)/Izz; f7=x(4)+tan(x(8)) * (x(5) * sin(x(

  a. * lambda. * (theta_0/6-lambda./4); CQp=sigma./8. * (CD_0+CD_i * theta_0^2). * (1+mu. * mu)-CD_i. * (lambda. * theta_0/3-lambda. * lambda./4); CQ=CQi+CQp; %sin and cos seno=[sin(epsilon); sin(pi/2+epsilon); sin(pi+epsilon); sin(3 * pi/2+epsilon)]; coseno=[cos(epsilon); cos(pi/2+epsilon); cos(pi+epsilon); cos(3 * pi/2+epsilon)]; %aerodynamic forces and momentum %avoid singularity Fxj=-rho. * A. * R^2. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fyj=-rho. * A. * R^2. * ((x(2)-wind(2))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fzj=-rho. * A. * R^2. * u. * CT; Lj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CRm. * u; Mj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(2)-wind(2))/sqrt((x(2)-wind(2))^2+(x(2)-wind(2))^2)). * CRm. * u; Nj=[1;-1;1;-1]. * rho. * A. * R^3. * CQ. * u; %aeroynamics F=[sum(Fxj);sum(Fyj);sum(Fzj)]; J=[sum(Lj)-h * sum(Fyj)+l * sum(Fzj. * seno);... sum(Mj)+h * sum(Fxj)-l * sum(Fzj. * coseno);... sum(Nj)+l * sum(Fyj. * coseno)-l * sum(Fxj. * seno)]; %dynamic equations f1=F(1)/m-g * sin(x(8))+(x(3) * x(5)-x(6) * x(2)); f2=F(2)/m+cos(x(8)) * sin(x(7)) * g-(x(6) * x(1)-x(3) * x(4)); f3=F(3)/m+cos(x(8)) * cos(x(7)) * g-(x(3) * x(2)-x(5) * x(2)); f4=J(1)/Ixx-x(5) * x(6) * (Izz-Iyy)/Ixx; f5=J(2)/Iyy-x(4) * x(6) * (Ixx-Izz)/Iyy; f6=J(3)/Izz-x(4) * x(5) * (Iyy-Ixx)/Izz; f7=x(4)+tan(x(8)) * (x(5) * sin(x(7))+x(6) * cos(x(7))); f8=x(5) * cos(x(7))-x(6) * sin(x(7)); f9=(x(5) * sin(x(7))+x(6) * cos(x(7)))/cos(x(8));%--------------------------------------------------------------------------------------%linear system matrices AM = jacobian([f1, f2, f3, f4, f5, f6, f7, f8, f9], [x(1), x(2), x(3), x(4), x(5), x(6), x(7), ... x(8), x(9)]); BM = jacobian([f1, f2, f3, f4, f5, f6, f7, f8, f9], [u(1), u(2), u(3), u(41, -2, -3, -4, -5, -6, -7, -8, -9]; K=place(AM,BM,eigenvalues); %find P and optimize Linear system epsilon = 1e-3; G = [P ≥ epsilon * eye(9) , (AM-BM * K)' * P+P * (AM-BM * K)≤ -epsilon * eye(9)]; optimize(G,[],sdpsettings('solver','sdpt3')); %find ∆ and optimize NonLinear system for ii=1:512 ∆x=matrix(:,ii). * ∆;

  %aerodynamic forces and momentum %avoid singularity Fxj=-rho. * A. * R^2. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fyj=-rho. * A. * R^2. * ((x(2)-wind(2))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fzj=-rho. * A. * R^2. * u. * CT; Lj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CRm. * u; Mj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(2)-wind(2))/sqrt((x(2)-wind(2))^2+(x(2)-wind(2))^2)). * CRm. * u; Nj=[1;-1;1;-1]. * rho. * A. * R^3. * CQ. * u; %aeroynamics F=[sum(Fxj);sum(Fyj);sum(Fzj)]; J=[sum(Lj)-h * sum(Fyj)+l * sum(Fzj. * seno);... sum(Mj)+h * sum(Fxj)-l * sum(Fzj. * coseno);... sum(Nj)+l * sum(Fyj. * coseno)-l * sum(Fxj. * seno)]; %dynamic equations f1=F(1)/m-g * sin(x(8))+(x(3) * x(5)-x(6) * x(2)); f2=F(2)/m+cos(x(8)) * sin(x(7)) * g-(x(6) * x(1)-x(3) * x(4)); f3=F(3)/m+cos(x(8)) * cos(x(7)) * g-(x(3) * x(2)-x(5) * x(2)); f4=J(1)/Ixx-x(5) * x(6) * (Izz-Iyy)/Ixx; f5=J(2)/Iyy-x(4) * x(6) * (Ixx-Izz)/Iyy; f6=J(3)/Izz-x(4) * x(5) * (Iyy-Ixx)/Izz; f7=x(4)+tan(x(8)) * (x(5) * sin(x(7))+x(6) * cos(x(7))); f8=x(5) * cos(x(7))-x(6) * sin(x(7)); f9=(x(5) * sin(x(7))+x(6) * cos(x(7)))/cos(x(8)); F=[f1;f2;f3;f4;f5;f6;f7;f8;f9];%--------------------------------------------------------------------------------------L = [∆x≥epsilon * ones(9,1), 2. * ∆x' * P * F≤-epsilon * ones(9,1)]; optimize(L,[],sdpsettings('solver','sdpt3')); %∆ value for each iteration: first element of each column ∆_value(ii)=abs(∆x(1)); end %radius of the AoA: minimum of a vector which has 512 elements AoA_radius=min(∆_value); B Feedback control in simulink/matlab function [K_u_x] = fun_wind_effect(vector) syms omega1_2 omega2_2 omega3_2 omega4_2 U0 v0 w0 p0 q0 r0 phi0 theta0 psi0 t1 t2 t3 t4 t5 t6

Figure 10 :

 10 Figure 10: Overall simulink scheme of the feedback control

Fyj=0

  . * rho. * A. * R^3. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CRm. * u; Mj=0. * [1;1;1;1]; Nj=[1;-1;1;-1]. * rho. * A. * R^3. * CQ. * u; else Fxj=-rho. * A. * R^2. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fyj=-rho. * A. * R^2. * ((x(2)-wind(2))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fzj=-rho. * A. * R^2. * u. * CT; Lj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CRm. * u; Mj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(2)-wind(2))/sqrt((x(2)-wind(2))^2+(x(2)-wind(2))^2)). * CRm. * u; Nj=[1;-1;1;-1]. * rho. * A. * R^3. * CQ. (Fxj);sum(Fyj);sum(Fzj)]; J=[sum(Lj)-h * sum(Fyj)+l * sum(Fzj. * seno);... sum(Mj)+h * sum(Fxj)-l * sum(Fzj. * coseno);... sum(Nj)+l * sum(Fyj. * coseno)-l * sum(Fxj. * seno)(@(tt) fcn_change(tt,FU,t1,t2,t3,t4,t5,t6),t0); %System equations--------------------------------------------------------------------------------------u=[omega1_2;omega2_2;omega3_2;omega4_2]; x=[U0,v0,w0,p0,q0,r0,phi0,theta0,psi0]; mu=sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)./(R. * abs(sqrt(u))); a. * lambda. * (theta_0/6-lambda./4); CQp=sigma./8. * (CD_0+CD_i * theta_0^2). * (1+mu. * mu)-CD_i. * (lambda. * theta_0/3-lambda. * lambda./4(epsilon); sin(pi/2+epsilon); sin(pi+epsilon); sin(3 * pi/2+epsilon)]; coseno=[cos(epsilon); cos(pi/2+epsilon); cos(pi+epsilon); cos(3 * pi/2+epsilon)]; %aerodynamic forces and momentum %avoid singularity Fxj=-rho. * A. * R^2. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fyj=-rho. * A. * R^2. * ((x(2)-wind(2))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CH. * u; Fzj=-rho. * A. * R^2. * u. * CT; Lj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(1)-wind(1))/sqrt((x(1)-wind(1))^2+(x(2)-wind(2))^2)). * CRm. * u; Mj=[1;-1;1;-1]. * rho. * A. * R^3. * ((x(2)-wind(2))/sqrt((x(2)-wind(2))^2+(x(2)-wind(2))^2)). * CRm. * u; Nj=[1;-1;1;-1]. * rho. * A. * R^3. * CQ. * u; %aeroynamics F=[sum(Fxj);sum(Fyj);sum(Fzj)]; J=[sum(Lj)-h * sum(Fyj)+l * sum(Fzj. * seno);... sum(Mj)+h * sum(Fxj)-l * sum(Fzj. * coseno);... sum(Nj)+l * sum(Fyj. * coseno)-l * sum(Fxj. * seno)sin(x(7))+x(6) * cos(x(7)))/cos(x(8));

Figure 15 :

 15 Figure 15: "Flight dynamics" simulink submodel

Figure 17 :

 17 Figure 17: Fully configurable mask of the implemented simple simulink model. (Values are an example).

Figure 18 :

 18 Figure 18: From control to rotors.

Figure 19 :

 19 Figure 19: Fully configurable mask of the model from control to rotors angular velocity. (Values are an example).

C Quadrotor models in simulink