Induced Force Hovering of Spherical Robot by Under-Actuated Control of Dual Rotor

Abstract

This chapter discusses the design and modelling of a spherical flying robot. The main objective is to control its hovering and omnidirectional mobility by controlling the air mass differential pressure between two asynchronous coaxial rotors that are aligned collinearly. The spherical robot design has embedded a gyroscopic mechanism of three rings that allow the rotors’ under-actuated mobility with 3DOF. The main objective of this study is to maintain the thrust force with nearly vertical direction. The change in pressure between rotors allows to vary the rotors’ tilt and pitch. The system uses special design propellers to improve the laminar air mass flux. A nonlinear fitting model automatically calibrates the rotors’ angular speed as a function of digital values. This model is the functional form that represents the reference input to control the rotors’ speed, validated by three types controllers: P, PI, and PID. The robot’s thrust and induced forces and flight mechanics are proposed and analysed. The simulation results show the feasibility of the approach.

1. Introduction

Nowadays, unmanned aerial vehicle (UAV) robots are being deployed at an increased rate for numerous applications falling into a variety of engineering fields. There exist numerous kinds of rotary-wing-based robotic technologies in particular with active devices. Robots with rotating wings are capable to self-control over lift, propulsion, landing, hovering And take-off tasks [1–4]. Overcoming vertical flight (with minimum energy cost) is fundamental to accomplish autonomous precise tasks. One fundamental aspect in controlling and designing rotary-wing-based intelligent machines is to consider under-actuated issues to reduce the number of actuators. Under-actuated flying robots perform motion tasks more naturally, taking advantage of the inertial and gravitational forces, consequently, reducing the use of electrical energy. Biological flying birds are instances of extremely efficient under-actuated bodies. Therefore, in order to design flying machines with a reduced number of actuators, it is essential to model and understand the mechanical nature of the robot mechanics, the fluids and their physics-based relationship.

The present work has foundations on the prototype of a home-made spherical aerial robot. Some experiments can be viewed at https://www.youtube.com/watch?v=rrGH1Oh_beM. Nevertheless, the purpose of this chapter is not on showing and discussing experimental results, but on mathematically sustaining the hypothesis of the robot’s flight mechanics and control. Unlike, known spherical design approaches [5–7], rather than deploying aileron-like propellers, we proposed yaw, pitch and roll changes through under-actuation exerting an inner gyroscopic mechanism. In the present chapter, the authors are particularly interested in disclosing the physical model of a dual rotary-wing spherical robot with an under-actuated gyroscopic mechanism. The model has been divided into four major areas: the robot’s flight mechanics with direct and inverse solution, the thrusting or induced force model, the rotors control model and a proportional integral derivative (PID) based control with non-stationary reference values.

This chapter is organised as follows. In Section 2, the design and mechanical aspects of the aerial robot are presented. Section 3 presents the kinematic direct and inverse solutions of the flight mechanics. In Section 4, the acceleration components and forces involved in the robot’s aerodynamics are discussed. Section 5 presents the robot’s thrusting force model that involves two collinear induced forces. Section 6 presents the rotors’ actuator speed models that are proposed from empirical measurements, and subsequently, the analytical solution is obtained. In Section 7, the actuators’ feedback linear control is described. Finally, in Section 8, conclusions are drawn.

2. Spherical gyroscopic robot

Unlike other reported approaches [8–13], in this study, the authors have proposed an omnidirectional spherical design with two rotors vertically collinear (see Figure 1, right, and Figure 2).

Figure 1.

Left: Aerial robot real prototype. Right: Robot’s main actuated and under-actuated mechanical elements.

Figure 2.

Left: Robot’s global and local spherical coordinates. Centre: Robot’s flying space (Eqs. (1)–(3)). Right: recursive trajectory generation (Texto) and (Texto).

3. Flight mechanics model

Flight mechanics refers to the study of geometry of flight of a heavier-than-air aircraft, considering aerodynamic aspects. Expressions (1)–(3) model the three-dimension robot’s Cartesian kinematic components that describe its motion. The components, namely, x, y and z, are the space positions w.r.t. the location of the robot’s starting flight. The proposed kinematic model is constrained with initial posture as the inertial frame origin, where d is the distance between the robot’s instantaneous 3D position and its Cartesian origin. Azimuth angle ϕ₀ is w.r.t. the plane XY, and the elevation angle ϕ₁ is w.r.t. the Y-axis:

$x=d\cos \left({\varphi }_1\right)\sin \left({\varphi }_0\right)$E1

The y component (vertical) is expressed as

$y=d\sin \left({\varphi }_1\right)$E2

And, the z component is expressed as

$z=d\cos \left({\varphi }_1\right)\cos \left({\varphi }_0\right)$E3

For further purpose, the inverse solution is obtained by an algebraic arrangement of derivatives. The first-order derivatives of Eqs. (1)–(3) are obtained and shown in Eqs.(4)–(6),

$\dot{x}=\dot{d}\cos \left({\varphi }_1\right)\sin \left({\varphi }_0\right)+d\cos \left({\varphi }_1\right)\cos \left({\varphi }_0\right){\dot{\varphi }}_0-d\sin \left({\varphi }_1\right)\sin \left({\varphi }_0\right){\dot{\varphi }}_1$E4

The vertical component is expressed as

$\dot{y}=\dot{d}\sin \left({\varphi }_1\right)+d\cos \left({\varphi }_1\right)\dot{\phi }$E5

And, the z component is expressed as

$\dot{z}=\dot{d}\cos \left({\varphi }_1\right)\cos \left({\varphi }_0\right)+d\cos \left({\varphi }_1\right)\cos \left({\varphi }_0\right){\dot{\varphi }}_0$E6

Figure 2 centre depicts the robot’s flying space, which is spherical with the Cartesian origin at robot’s starting flying task.

Expressing in the matrix form, the first-order derivative dp/dt w.r.t. time is

$\dot{p}=\left(\begin{array}{c}\dot{x}\\ \dot{y}\\ \dot{z}\end{array}\right)=\left(\begin{array}{ccc}d{C}_1{C}_0& -d{S}_1{S}_0& C_1{S}_0\\ 0& d{C}_1& S_1\\ -d{C}_1{C}_0& -d{S}_1{C}_0& C_1{C}_0\end{array}\right)\stackrel{\cdot }{}\left(\begin{array}{c}{\dot{\varphi }}_0\\ {\dot{\varphi }}_1\\ \dot{d}\end{array}\right)$E7

By simplifying, the linear equation of direct kinematic is denoted by the Jacobian matrix J and the first-order vector of independent variables,

$\dot{p}=J\dot{\Phi }$E8

In order to obtain a recursive functional form equivalent to previous state variables, the derivatives are expressed in the following manner:

$\frac{d}{dt}p=J\frac{d}{dt}\Phi$E9

Time differentials are eliminated, and the integration operators complete the remaining differentials dp and dϕ:

${\displaystyle {\int }_{{\dot{p}}_1}^{{\dot{p}}_2}dp=J}{\displaystyle {\int }_{{\dot{\varphi }}_1}^{{\dot{\varphi }}_2}d\Phi }$E10

Thus, to solve for the robot’s Cartesian position a recursive form is obtained by algebraically reordering. Next, robot’s position p_t+1 is obtained by successive approximations of ϕ_t until ϕ_f:

$p_{t+1}=p_t+J\left({\varphi }_f-{\varphi }_t\right)$E11

In addition, the kinematic inverse solution requires the inverse-squared Jacobian matrix, assuming that it is an invertible and non-singular matrix, with non-zero determinant:

$J^{-1}=\frac{1}{d^2{C}_1\left({\left[C_1{C}_0\right]}^2+\left[S_0{S}_1\right]\right)}\left(\begin{array}{ccc}d{C}_0& 0& -d{S}_0\\ d{C}_1{S}_0{S}_1& d{C}_1^2& -d{C}_1{C}_0{S}_1\\ {\left(d{C}_1\right)}^2{S}_0& -d{C}_1{C}_0{S}_1& {\left(d{C}_1\right)}^2{C}_0\end{array}\right)$E12

Therefore, the first-order inverse kinematic is obtained by an algebraic approach:

$\dot{\Phi }=J^{-1}\dot{p}$E13

As described earlier, we complete the differentials dp and dϕ:

${\displaystyle {\int }_{{\dot{\Phi }}_1}^{{\dot{\Phi }}_2}d\Phi =J^{-1}}{\displaystyle {\int }_{{\dot{p}}_1}^{{\dot{p}}_2}dp}$E14

By integrating both sides of the equality, respectively, a recursive inverse solution in terms of the rotor’s angular speed is obtained,

${\dot{\Phi }}_{t+1}={\dot{\Phi }}_t+J^{-1}\left({\dot{p}}_f-{\dot{p}}_t\right)$E15

where p_t is the actual robot 3D position and p_f represents the final desired position in space. To achieve such location, the robot recursively approximates the next desired rotor’s controlled velocity. The next figure depicts a simulation result where the aerial robot successively approached the final desired position, staring from the Cartesian origin.

4. Aerodynamic robot’s model

The aerodynamic robot’s model refers to the application of the Newton’s second law of motion in three dimensions to infer the thrusting force T and other involved forces that produce the Cartesian accelerations:

$\begin{array}{c}\ddot{x}=\ddot{d}{c}_1{s}_0-\dot{d}{s}_1{s}_0{\dot{\varphi }}_1+\dot{d}{c}_1{c}_0{\dot{\varphi }}_0+\dot{d}{c}_1{c}_0{\varphi }_1-d{s}_1{c}_0{\dot{\varphi }}_0{\dot{\varphi }}_1-d{c}_1{c}_0{\dot{\varphi }}_0^2\\ -d{c}_1{c}_0{\ddot{\varphi }}_0-\dot{d}{s}_1{s}_0{\dot{\varphi }}_1-d{c}_1{s}_0{\dot{\varphi }}_1^2-d{s}_1{c}_0{\dot{\varphi }}_1{\dot{\varphi }}_0-d{s}_1{s}_0{\dot{\varphi }}_1\end{array}$E16

$\ddot{y}=\ddot{d}{s}_1+d{c}_1{\dot{\varphi }}_1+\dot{d}{c}_1{\varphi }_1-d{s}_1{\dot{\varphi }}_1^2+d{c}_1{\ddot{\varphi }}_1$E17

$\begin{array}{c}\ddot{z}=\ddot{d}{c}_1{c}_0-\dot{d}{s}_1{c}_0{\dot{\varphi }}_1-\dot{d}{c}_1{s}_0{\dot{\varphi }}_0-d{s}_1{c}_0{\varphi }_1-d{c}_1{c}_0{\dot{\varphi }}_1^2+d{s}_1{s}_0{\dot{\varphi }}_0{\dot{\varphi }}_1\\ -d{s}_1{c}_0{\ddot{\varphi }}_1-\dot{d}{c}_1{s}_0{\dot{\varphi }}_0+\dot{d}{s}_1{s}_0{\dot{\varphi }}_0{\dot{\varphi }}_1-d{c}_1{c}_0{\dot{\varphi }}_0^2-d{c}_1{s}_0{\ddot{\varphi }}_0\end{array}$E18

In the matrix form, the following equation represents the direct kinematic solution for the Cartesian accelerations:

$\begin{array}{c}\ddot{p}=\left(\begin{array}{ccc}d{c}_1{c}_0& -d{s}_1{s}_0& c_1{s}_0\\ 0& d{c}_1& s_1\\ -d{c}_1{c}_0& -d{s}_1{c}_0& c_1{s}_0\end{array}\right)\left(\begin{array}{c}{\dot{\varphi }}_0\\ {\dot{\varphi }}_1\\ \dot{d}\end{array}\right)+\left(\begin{array}{ccccc}-d{c}_1{c}_0& -d{c}_1{c}_0& -2d{c}_1{c}_0& 2d{c}_1{c}_0& -2d{s}_1{s}_0\\ 0& -d{s}_1& 0& 0& 2c_1\\ -d{c}_1{c}_0& -d{c}_1{c}_0& -2d{s}_1{s}_0& 2d{c}_1{s}_0& -2d{s}_1{c}_0\end{array}\right)\left(\begin{array}{c}{\dot{\varphi }}_0^2\\ {\dot{\varphi }}_1^2\\ {\dot{\varphi }}_0{\dot{\varphi }}_1\\ \dot{d}{\dot{\varphi }}_0\\ \dot{d}{\dot{\varphi }}_1\end{array}\right)\end{array}$E19

or

$\ddot{p}=J_1\cdot \dot{\Phi }+J_2\cdot \ddot{\Phi }$

From the previous expression, the vector acceleration is substituted into the Newton’s second law of motion,

$f=m\cdot \ddot{p}$E20

where the Cartesian force components defined in robot’s local inertial frame are expressed as follows:

$f_x=T\cos \left({\theta }_1\right)\sin \left({\theta }_0\right)$E21

The force component along the robot’s local y component is expressed as

$f_y=T\sin \left({\theta }_1\right)$E22

And, the force component along the robot’s local Z component is expressed as

$f_z=T\cos \left({\theta }_1\right)\cos \left({\theta }_0\right)$E23

By simplifying the Cartesian force components in the matrix form

$f=T\left(\begin{array}{c}{c}_1{s}_0\\ s_1\\ c_1{c}_0\end{array}\right)$E24

And, by substituting the functional form of the vector force f into the Newton’s second law,

$m\cdot \ddot{p}=T\left(\begin{array}{c}{c}_1{s}_0\\ s_1\\ c_1{c}_0\end{array}\right)$E25

Thus, by reordering the previous equation, we substitute the vector constraints w^T =(C₁S₀, S₁, C₁C₀). The acceleration vector d²p/dt² is a function of the next position p_t+1 and the rotors variables:

$Tw=m\ddot{p}\left(p_{t+1},{\Phi }_f,{\dot{\Phi }}_t\right)$E26

By dropping off the induced robot’s force T,

$T=w^{-1}m\ddot{p}\left(p_{t+1},{\Phi }_f,{\dot{\Phi }}_t\right)$E27

So far, in this expression, the total thrusting induced force T represents the robot’s global flying force. Thus, T is an arithmetic result produced by the sum of the top rotor’s induced force T₁ and the below rotor’s induced force T₂ according to the following governing constraints:

1. (a)

   For T₁ = T₂, the inflow air mass is same throughout both rotors, hence T = T₁ + T₂.

2. (b)

   For T₁ > T₂, speed and air mass below rotor 1 are greater than rotor 2 inflow, T = T₁ + α₂ T₂.

3. (c)

   For T₁<T₂, opposed to constraint (b), then T₂ = α₂T₁ and T = T₁(1+α₂).

Here, the numerical factors α_1,2 are gains denoting rotors’ speed-rate differences. For either constraint (b) or (c), the gyroscopic mechanism angles’ tilt and pitch are affected, consequently changing the robot’s azimuth and elevation angles.

5. Induced force model

According to the depictions of Figure 3, the rotors are continuously pushing the air down. As per Newton’s third law, an equal and opposite reaction force, denoted as rotor thrust, is acting on the rotor due to air. The induced force model refers to the thrusting force exerted to accelerate the robot. And at a constant velocity the quasi-static hovering is achieved [1–4].

Figure 3.

Rotors’ flow conditions in the slip streams.

The momentum conservation is obtained by relating the induced force T₂ to the rate of momentum change. It is the mass rate and the far-field wake-induced velocity v_w below rotor 2, where dm/dt=ρAv₂ and the rotor disk area A = πR². Thus, the moment conservation

$T=\dot{m}{v}_w$E28

The energy conservation per unit time

$Tv=\frac{1}{2}\dot{m}{v}_w^2$E29

To obtain a relationship between v and v_w, let us substitute T and dm/dt,

$\dot{m}{v}_{w}v=\frac{1}{2}\dot{m}{v}_w^2$E30

Algebraically simplifying,

$v=\frac{1}{2}{v}_w;v_w=2v$E31

Hence, substituting v_w into T,

$T=\dot{m}{v}_w=2\dot{m}v$E32

Then,

$T=2\rho A{v}^2$E33

Dropping off v, the following expression is obtained;

$v=\sqrt{\frac{T}{2\rho A}}$E34

The propulsive power P_w is the thrusting force T capable to move the robot at a given velocity (distance over time):

$P_w=Tv=T\sqrt{\frac{T}{2\rho A}}$E35

The induced power per unit thrust for a hovering rotor can be written as

$\frac{P_w}{T}=v=\sqrt{\frac{T}{2\rho A}}$E36

The above expression indicates that, for a low inflow velocity, the efficiency is higher. This is possible if the rotor has a low disk loading (T/A). Note that the parameter determining the induced power is essentially T/(ϱA). Therefore, the effective disk loading increases with an increase in altitude and temperature.

From previous analysis, let us now precisely define the thrusting force for rotors 1 and 2, according to Figure 3 (left side). For rotor 1, the air mass flow is

${\dot{m}}_1=\rho A_1{v}_1$E37

Hence, considering only rotor 1, the induced force is

$T_1={\dot{m}}_1{v}_2$E38

The energy conservation principle for rotor 1 is expressed as

$T_1{v}_1=\frac{1}{2}{\dot{m}}_1{v}_2^2$E39

And finding a relationship between v₁ and v₂ in accordance with Figure 3 (left side),

$v_2=2v_1$E40

Thus,

$T_1=2{\rho }_1{A}_1{v}_1^2$E41

The induced air velocity induced by rotor 1 is modelled by

$v_1=\sqrt{\frac{T_1}{2{\rho }_1{A}_1}}$E42

Similarly, modelling both the induced force T₂ and velocity v₂ for rotor 2, the following analysis is developed. The airflow rate,

${\dot{m}}_2={\rho }_2{A}_2{v}_3$E43

And the induced force T₂ considers the inflow air mass and the far-field wake-induced velocity v_w,

$T_2={\dot{m}}_2{v}_w$E44

The energy conservation for the second rotor is expressed as

$T_2{v}_3=\frac{1}{2}{\dot{m}}_2{v}_w^2$E45

The relationship between far-field wake-induced air velocity v_w and v₃ is given by v_w=2v₃, and

$T_2=2{\rho }_2{A}_2{v}_3^2$E46

Therefore, the second rotor’s air induced velocity is

$v_3=\sqrt{\frac{T_2}{2{\rho }_2{A}_2}}$E47

From the three previous postulates of Figure 3, let us deduce the conditions when both rotors, although asynchronous, simultaneously induce the airflow equally, when v₁=v₃ (Figure 3a). For this case, the total robot’s thrusting force T=T₁+T₂,

$T=2{\rho }_1{A}_1{v}_1^2+2{\rho }_2{A}_2{v}_3^2$E48

For case 1, let us assume A₁=A₂ and v₂=v₃ through the second rotor’s disc area. And,

$v_w=2v_3;\text{likewise},v_2=2v_1;\text{hence}{v}_w=2\left(2v_1\right)=4v_1$E49

Rewriting total T as a function of v₁, thus we have

${\rho }_2=\frac{m_c}{v_c}$E50

The air mass variation is denoted as the mass derivative w.r.t. time,

$\dot{m}={\rho }_{2}A{V}_2$E51

Now, for the case v₁>v₃,

$T=T_1+{\alpha }_2{T}_2$E52

Let us substitute our T₁ and T₂ models previously analysed,

$T=2{\rho }_1{A}_1{v}_1^2+2{\rho }_2{A}_2{v}_3^2$E53

In addition, A₁=A₂, v₂=2v₁ and v_w=2v₃, but v₂>v₃ hence v₃=α₂v₂. Thus, v_w=2(α₂ 2(2v₁)), and we obtain the relationship between the far-field wake-induced velocity and v₁:

$v_w=6{\alpha }_2{v}_1$E54

By substituting v₁ into expression T, developing and factorising algebraically,

$T=2A{v}_1^2\left({\rho }_1+4{\alpha }_2{\rho }_2\right)$E55

Similarly, for the case when v₁<v₃, T=T₁(1+α₂),

$T=2{\rho }_1{A}_1{v}_1^2+2{\rho }_2{A}_2{v}_3^2$E56

for this case, A₁=A₂, v₂=2v₁ and v_w=2v₃. However, just above and below rotor 2, v₂<v₃, and therefore v₂=α₂v₃, then v₃=2v₂/α₂.

$T=2{\rho }_{1}A{v}_1^2+2{\rho }_{2}A{\left(2\frac{v_1}{{\alpha }_2}\right)}^2$E57

Factorizing and algebraically arranging,

$T=2A{v}_1^2\left({\rho }_1+\frac{8{\rho }_2}{{\alpha }_2^2}\right)$E58

Solving the parameter α₂ for each of the three cases,

${\alpha }_2^2=\{\begin{array}{l}4,\\ 1,\\ \sqrt{2},\end{array}\begin{array}{l}{v}_1=v_3\\ v_1>v_3\\ v1<v_3\end{array}|$E59

Therefore, we synthesise the total thrusting force T for all cases by

$T=2A{v}_1^2\left({\rho }_1+{\alpha }_2^2{\rho }_2\right)$E60

Therefore, from Eq. (60), now we have a functional form for the thrusting force T. Thus, to reach a controlled rotors’ velocity, we must establish a relationship between the induced velocity v₁ and the rotor angular velocity dϕ/dt using the tip speed of the rotor blade as reference. The rotor inflow is represented in non-dimensional form as

$\lambda =\frac{v_3}{\dot{\varphi }R}=\sqrt{\frac{C_T}{2}}$E61

where C_T is the thrust coefficient modelled by

$C_T=\frac{T}{{\rho }_{2}A{\left({\dot{\varphi }}_{t}R\right)}^2}$E62

Therefore, the following equality allow us to deduce C_T with more

$\dot{\varphi }=\frac{v_1}{R\sqrt{\frac{C_T}{2}}}$E63

Thus, substituting (Texto) into v₁,

$v_1={\left(\dot{\varphi }R\sqrt{\frac{C_T}{2}}\right)}^2$E64

And subsequently, v₁² is substituted into T,

$T=A\left({\dot{\varphi }}^2{R}^2{C}_T\right)\left({\rho }_1+{\alpha }_2^2{\rho }_2\right)$E65

Therefore, the induced force T arising from the fluid mechanics equation (65) is equated with the induced force of the flight mechanics equation (27), the following expression is obtained:

$2A\left({\dot{\varphi }}^2{R}^2\frac{C_T}{2}\right)\left({\rho }_1+{\alpha }_2^2{\rho }_2\right)=w^{-1}m\ddot{p}$E66

Our objective is to find an analytical solution for the rotor speed required to reach the induced velocity v₁ and the induced force T, which is governed by the flight mechanics law. The induced air mass velocity v₁ can be expressed in terms of the rotor speed that is controlled to obtain the desired angular velocity,

$\dot{\varphi }=\sqrt{\frac{m{w}^{-1}\ddot{p}}{A{R}^2{C}_T\left({\rho }_1+{\alpha }_2^2{\rho }_2\right)}}$E67

This model represents the independent variable to control the motors’ Speed of Eq. (73) or Eq. (74):

Thus, it follows a set of empirical temperature measurements where some experimental hovering experiments were carried out. The plots in Figure 4 depict how the air pressure is affected as the temperature varies over time.

Figure 4.

Left: Thrusting force w.r.t. induced velocity. Right: Air mass density as temperature varies.

6. Actuators’ speed model

The actuators’ self-calibration speed model is discussed in this section. The real rotary velocity in a range from minimal to maximal values approached a logarithmic model. It considered the empirical set of angular speed measurements w.r.t. digital controls. The inherent physical variations, such as temperature, air pressure, density and air dust particles, affected the actuators’ performance. Since the angular speed value capable to hover the spherical robot’s body is disturbed, a self-calibration is required to maintain position control as accurate as possible. From experiments, the empirical models that obtained (Figure 5) φ vs. d are fitted according to the next model. The parameters A and β are unknown and must fit the speed measurements φ, w.r.t to digital word d,

$\dot{\varphi }=A+\beta \ln \left(d\right)$E68

We temporally substitute d’=ln(d), and thus the rotary speed

$\dot{\varphi }=A=\beta d\text{'}$E69

To estimate the unknown parameter β, a linear mean-squared method is applied,

$\beta =\frac{n\left({{\displaystyle \sum }}^{\text{}}d\text{'}\dot{\varphi }-\left({{\displaystyle \sum }}^{\text{}}d\text{'}\right)\left({{\displaystyle \sum }}^{\text{}}\dot{\varphi }\right)\right)}{n\left({{\displaystyle \sum }}^{\text{}}d{\text{'}}^2\right)-{\left({{\displaystyle \sum }}^{\text{}}d\text{'}\right)}^2}$E70

Subsequently, the previous parameter solution is used in the next expression A,

$A=\frac{{\displaystyle \sum \dot{\varphi }}}{n}-\beta \frac{{{\displaystyle \sum }}^{\text{}}d\text{'}}{n}$E71

By substituting the parameters numerical values, the rotary speed model is obtained as follows:

$\dot{\varphi }=10838.57637\ln \left(d\right)-46776.63059$E72

In addition, in order to obtain the inverse solution, we algebraically drop off the variable d from Eq. (72),

$d=e^{\frac{\dot{\varphi }-A}{\beta }}$E73

And numerical parametric values from Eq. (72) are substituted into (73) to obtain

$d=e^{\frac{\dot{\varphi }-46776.6306}{10838.5764}}$E74

In order to compare how our theoretical model fits the empirical model, both inverse and direct operation control modes are depicted in Figure 5.

Figure 5.

Actuator’s raw measurements. Left: Direct solution. Right: Inverse solution.

7. Rotor’s speed control

From the previous section, the actuators’ speed model is now used to formulate a feedback linear control. Let us assume that a rotor control variable (i.e., angle, velocity and acceleration) should ideally be equal to the real sensed control variable, as expressed by the next equality (75). Nevertheless, in a realistic scenario real and ideal control variables are different due to a number of factors, such as frictions, inertial and gravity forces, motor electromagnetic performance and so on. Thus, both variables are approached by a multiplicative gain or factor alpha, which approximates both numerical values according to the relation

$\phi =k\widehat{\dot{\phi }}$E75

Assuming an arbitrary actual derivative order, the equation is equivalently expressed as

$\frac{d}{dt}\phi =k\frac{d}{dt}\widehat{\phi }$E76

The time differentials are eliminated and the remaining differentials are obtained by solving the following definite integrals,

${\displaystyle {\int }_{{\phi }_1}^{{\phi }_2}d\phi =k}{\displaystyle {\int }_{{\widehat{\phi }}_1}^{{\widehat{\phi }}_2}d\widehat{\phi }}$E77

By solving the definite integrals,

${\phi }_2-{\phi }_1=k\left({\widehat{\phi }}_2-{\widehat{\phi }}_1\right)$E78

In this case, ${\widehat{\phi }}_2$ is the expected or the reference value to be reached ideally, ${\widehat{\phi }}_2={\phi }^{\text{ref}}$. Thus, by adjusting the times sub-interval labels, and algebraically reordering, the next recursive numerical successive approximation equation is expressed as follows:

${\phi }_{t+1}={\phi }_t+k\left({\phi }^{\text{ref}}-{\phi }_t\right)={\phi }_t+k{e}_t$E79

From the previous expression, the error ke_t is spanned into the past (angular displacement), the actual (rotary speed) and the future (angular accelerations) errors in order to cover the whole error history. And the general constant gain k is proportional to k_p, k_I and k_d. Thus,

${\phi }_{t+1}={\phi }_t+\left(k_p{e}_p+k_I{e}_I+k_d{e}_d\right)$E80

Therefore, the next feedback proportional error e_p (rad/s) with proportional gain k_p (dimensionless) is obtained with the observation ϕ_t measured online. And the reference model (Texto) is established in terms of the instantaneous control word δ_t:

$k_p{e}_p=k_p\left({\dot{\varphi }}^{\text{ref}}-{\dot{\varphi }}_t\right)=k_p\left(a\ln \left(\delta \right)-b-{\dot{\varphi }}_t\right)$E81

Figure 6.

Left: Measurement and reference proportional models. Right: Proportional error.

For illustrative purpose, an accelerative rotor’s task to exert robot’s propulsion was performed. Figure 6 (left) depicts both the reference model ϕ^ref and the observation ϕ(t).

In addition, Figure 6 (right) shows the proportional error behaviour ϕ^ref-ϕ(t) without k_p. Furthermore, the feedback integral error e_I (rad/s) with integral gain k_I (dimensionless) is expressed by the time integration of the difference of (d²ϕ/dt² – d²ϕ(t)/dt²).

$k_I{e}_I=k_I\left({{\displaystyle \int }}^{\text{}}\left({\ddot{\varphi }}^{\text{ref}}-{\ddot{\varphi }}_t\right)dt\right)$E82

The observation model d²ϕ/dt² was obtained online by the numerical derivatives of the optical encoder according to the following relationship:

$\ddot{\varphi }\left(t\right)=\frac{{\dot{\varphi }}_{t_2}-{\dot{\varphi }}_{t_1}}{t_2-t_1}$E83

In addition, since the observation model inherently poses perturbations, an analytical reference model d²ϕ^ref/dt² was obtained using a nonlinear regressive fitting process for parameters identification,

$\left(\begin{array}{c}{\displaystyle \sum {\ddot{\varphi }}_i}\\ {\displaystyle \sum \ddot{\varphi }{t}_i}\\ {\displaystyle \sum {\ddot{\varphi }}_i{t}_i^2}\end{array}\right)=\left(\begin{array}{ccc}n& {{\displaystyle \sum }}^{\text{}}{t}_i& {{\displaystyle \sum }}^{\text{}}{t}_i^2\\ {{\displaystyle \sum }}^{\text{}}{t}_i& {{\displaystyle \sum }}^{\text{}}{t}_i^2& {{\displaystyle \sum }}^{\text{}}{t}_i^3\\ {{\displaystyle \sum }}^{\text{}}{t}_i^2& {{\displaystyle \sum }}^{\text{}}{t}_i^3& {{\displaystyle \sum }}^{\text{}}{t}_i^4\end{array}\right)\stackrel{\cdot }{}\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)$E84

where the previous expression is similarly expressed as

$y=A\stackrel{\cdot }{}x$E85

And by solving vector x,

$x=A^{-1}\stackrel{\cdot }{}y$E86

Hence, the reference model is a theoretical nonlinear function of time,

${\ddot{\varphi }}^{\text{ref}}\left(t\right)=34.7950-1.2832t+0.0145t^2$E87

Therefore, for the sake of the integral control u_I (rad/s), the angular acceleration reference model (Texto) is substituted next in its general form:

$u_I=k_I\left({{\displaystyle \int }}^{\text{}}\left(a_0-a_{1t}+a_2{t}^2\right)-\ddot{\varphi }\left(t\right)dt\right)$E88

Finally, in order to keep data homogeneity (numerical data subtraction), time integration is obtained by the trapezoid rule for numerical integration,

${\dot{\varphi }}_k=\frac{t_f-t_0}{2n}\left[\left({\ddot{\varphi }}_0^{\text{ref}}-{\ddot{\varphi }}_0\left(t\right)\right)+2{\displaystyle \sum \left({\ddot{\varphi }}_t^{\text{ref}}-{\ddot{\varphi }}_t\right)}+\left({\ddot{\varphi }}_{t_f}^{\text{ref}}\left(t\right)-{\ddot{\varphi }}_{t_f}\left(t\right)\right)\right]$E89

Figure 7.

Left: Measurement and reference integral models. Right: Integral error.

Figure 7 (left) depicts both reference and empirical models integrated with time. Figure 7 (right) shows the integral error behaviour.

In addition, the derivative control u_d=k_de_d (rad/s) with feedback derivative error e_d, and with derivative gain k_d (dimensionless), improves the closed-loop stability as follows:

$k_d{e}_d=k_d\left(\frac{d}{dt}{\varphi }^{\text{ref}}-\frac{d}{dt}{\varphi }_t\right)$E90

In order to obtain the time derivative observation model, the rotor’s angle evolution ϕ(t) (rad) is observed online using an optical encoder during the time slot where velocity and acceleration are also measured. As the measurements are read with noise, the analytical reference model is fitted as a nonlinear polynomial of the following form:

${\varphi }^{\text{ref}}\left(t\right)=-3709.6273+2116.9484t+27.9289t^2$E91

where the derivative error general form is expressed as

$k_d{e}_d=k_d\frac{d}{dt}\left(-b_0+b_{1t}+b_2{t}^2-\varphi \left(t\right)\right)$E92

Figure 8 (left) depicts both reference ϕ^ref(t) and observation ϕ(t) models together. Although both curves are apparently fitted, the vertical scale is provided in thousands of radians. Figure 8 (right) shows the derivative error scale.

Figure 8.

Left: Measurement and reference derivative models. Right: Derivative error.

Generally, the PID controller is expressed by the next expression

$u_t=k_p\left(\left(a\ln \left(d\right)-b\right)-{\dot{\varphi }}_i\right)+k_I\left(\underset{t}{{\displaystyle \int }}\left(a_0-a_{1t}+a_{2t}^2-{\ddot{\varphi }}_t\right)dt\right)+k_d\left(\frac{d}{dt}\left(-b_0+b_{1t}+b_{2t}^2-\varphi \left(t\right)\right)\right)$E93

Therefore, the controlled rotor’s velocity that is recursively calculated by dφ_t+1/dt= dφ_t/dt + u_t is applied, and we established the following controller choices: proportional (P), proportional integral (PI) and proportional-integral-derivative. Figure 9 (left) depicts the rotors’ angular speed without control and with three types of controllers. The constant parameters were adjusted accordingly to obtain such results. We can see that after 25 s the responses P and PI gradually converge w.r.t. the raw rotor’s speed (Figure 9, left).

Figure 9.

Controllers P, PI and PID. Left: Rotor’s angular speed. Right: Induced Cartesian forces.

In addition, with the controlled rotor’s speed output, the induced force is iteratively calculated by

$T\left({\dot{\phi }}_{t+1}\right)=A{\dot{\phi }}_{t+1}^2{R}^2{C}_T\left({\rho }_1+{\alpha }_2^2{\rho }_2\right)$E94

Thus, Figure 9 (right) depicts the induced component forces that are produced using three types of controllers.

8. Conclusions

This work briefly introduced the design of an aerial spherical robot with under-actuated gyroscopic mechanism. Although the purpose of this study was not to describe the robot flying and physical capabilities, the main objective was to demonstrate the induced force model deploying dual collinear rotary wings with no steering actuators. This study describes the following four major areas: the robot flight mechanics, the model of the induced thrusting forces, the self-calibration actuators’ Speed model validated with three types of controllers (P, PI, PID) to drive the rotors’ motor speed. The proposed aerodynamic mechanism poses neither ailerons nor propellers for steering control. Although the platform is a home-made laboratory prototype with special arrangements, the main focus of this chapter is to model the hypothesis of controlling the robot’s directions by varying rotors’ asynchronous speed. To achieve this, the under-actuated gyroscopic mechanism provides the ability to control its inner yaw, pitch and roll angles. Until this stage, the robot development is under an early control capability. However, this study presents mathematical solutions and simulation results to demonstrate the proposed aerodynamic hypothesis.