2.1 Kinematics model of the airship

The kinematics models of the airship⁽ Reference Chen and Duan ^{21 )} are

(1) $$\left[ {\matrix{ {\dot{x}_{g} } & {\dot{y}_{g} } & {\dot{z}_{g} } \cr } } \right]^{T} {\equals}{\bf{R}} _{{\mib{I}} }^{{\mib{B}} } \left[ {\matrix{ u & v & w \cr } } \right]^{T} $$

(2) $$\left[ {\matrix{ {\dot{\phi }} & {\dot{\theta }} & {\dot{\psi }} \cr } } \right]^{T} {\equals}{\bf{\Phi }} \left[ {\matrix{ p & q & r \cr } } \right]^{T} $$

(3) $$\openup 3pt{\bf R}_{{\mib{I}} }^{{\mib{B}} } {\equals}\left[ {\matrix{ {{\rm cos}\psi {\rm cos}\theta } & {{\rm cos}\psi {\rm sin}\theta {\rm sin}\phi {\minus}{\rm sin}\psi {\rm cos}\phi } & {{\rm cos}\psi {\rm sin}\theta {\rm cos}\phi {\plus}{\rm sin}\psi {\rm sin}\phi } \cr {{\rm sin}\psi {\rm cos}\theta } & {{\rm sin}\psi {\rm sin}\theta {\rm sin}\phi {\plus}{\rm cos}\psi {\rm cos}\phi } & {{\rm sin}\psi {\rm sin}\theta {\rm cos}\phi {\minus}{\rm cos}\psi {\rm sin}\phi } \cr {{\minus}{\rm sin}\theta } & {{\rm cos}\theta {\rm sin}\phi } & {{\rm cos}\theta {\rm cos}\phi } \cr } } \right]$$

where ϕ, θ, and ψ are the pitch, the yaw, and roll angles, respectively. $$[\matrix{ {x_{g} } & {y_{g} } & {z_{g} } \cr } ]$$ is the position in ERF. $${\mib{R}} _{{\mib{I}} }^{{\mib{B}} } $$ and Φ are the direction cosine transformation matrix and the angular transformation matrix, respectively. Using this two matrices, the velocity and the angular velocity are transferred from BRF to ERF. And both of them are detailed in (3) and (4)

(4) $${\bf \Phi }{\equals}\left[ {\matrix{ 1 & {{\rm sin}\phi {\rm tan}\theta } & {{\rm cos}\phi {\rm tan}\theta } \cr 0 & {{\rm cos}\phi } & {{\minus}{\rm sin}\phi } \cr 0 & {{\rm sin}\phi {\rm sec}\theta } & {{\rm cos}\phi {\rm sec}\theta } \cr} } \right]$$

2.2 Dynamics model of the airship

The dynamics model under wind disturbance is

(5) $${\mib{M\vec{x}}} {\equals}{\mib{F}} _{{{\mib{GB}} }} {\plus}{\mib{F}} _{{\mib{A}} } {\plus}{\mib{F}} _{{\mib{w}} } {\plus}{\mib{F}} _{{\mib{I}} } {\plus}{\mib{F}} _{{\mib{T}} } $$

where M is the inertia matrix and F _GB stands for the combined force of the buoyancy and the gravity. F _A , F _w , F _I, and F _T are the nominal aerodynamic force, the extra wind-induced aerodynamic force, the Coriolis force and the propeller force, respectively.

The mass matrix is presented as Equation (6), where I _x , I _y , I _z , I _xy , I _xz , and I _yz are the inertial moments of the airship. [x _G , y _G , z _G ] is the position of gravity center. m is the total mass of the airship and payload, and m _ii , i=1...6 are the added masses. The added mass is the inertia added to a system because when the vehicle is accelerated in a fluid, additional forces are required to increase the kinetic energy contained in the fluid. This effect appears as an apparent increase in the mass of the vehicle and is often referred to as added mass.⁽ Reference Azinheira, De Paiva and Bueno ^{22 )} The value of the added masses can be calculated by the potential flow theory.More details can be seen in Wang’s work⁽ Reference Wang ^{23 )}

(6) $${\bf{M}} {\equals}\left[ {\matrix{ {m{\plus}m_{{11}} } & 0 & 0 & 0 & {mz_{G} } & {{\minus}my_{G} } \cr 0 & {m{\plus}m_{{22}} } & 0 & {{\minus}mz_{G} } & 0 & {m_{{26}} {\plus}mx_{G} } \cr 0 & 0 & {m{\plus}m_{{33}} } & {my_{G} } & {m_{{35}} {\minus}mx_{G} } & 0 \cr 0 & {{\minus}mz_{G} } & {my_{G} } & {I_{x} {\plus}m_{{44}} } & {{\minus}I_{{xy}} } & {{\minus}I_{{yz}} } \cr {mz_{G} } & 0 & {m_{{53}} {\minus}mx_{G} } & {{\minus}I_{{xy}} } & {I_{y} {\plus}m_{{55}} } & {{\minus}I_{{xz}} } \cr {{\minus}my_{G} } & {m_{{62}} {\plus}mx_{G} } & 0 & {{\minus}I_{{xz}} } & {{\minus}I_{{yz}} } & {I_{z} {\plus}m_{{66}} } \cr } } \right]$$

(7) $$\openup 3pt{\bf F}_{{{\bf GB}}} {\equals}\left[ {\matrix{ {{\minus}(G{\minus}B){\rm sin}\theta } \cr {(G{\minus}B){\rm sin}\phi {\rm cos}\theta } \cr {(G{\minus}B){\rm cos}\phi {\rm cos}\theta } \cr {{\minus}(Gz_{G} {\minus}Bz_{B} ){\rm sin}\phi {\rm cos}\theta } \cr {{\minus}(Gz_{G} {\minus}Bz_{B} ){\rm sin}\theta {\minus}(Gx_{G} {\minus}Bx_{B} ){\rm cos}\phi {\rm sin}\theta } \cr {(Gx_{G} {\minus}Bx_{B} ){\rm sin}\phi {\rm cos}\theta } \cr } } \right]$$

The combined force of the buoyancy and the gravity is presented in Equation (7). and [x _B , y _B , z _B ] is the position of buoyancy center. G and B are the buoyancy and the gravity of the airship, respectively.

The Coriolis force is presented in Equation (8).

F _A and F _w are modeled in the next section.

(8) $$\openup 3pt{\bf F}_{{\bf I}} {\equals}\left[ {\matrix{ {{\minus}(m{\plus}m_{{33}} )wq{\plus}(m{\plus}m_{{22}} )vr{\minus}mz_{G} pr{\plus}mx_{G} (r^{2} {\plus}q^{2} ){\minus}my_{G} pq} \cr {{\minus}(m{\plus}m_{{11}} )ur{\plus}(m{\plus}m_{{33}} )wp{\minus}mz_{G} qr{\plus}my_{G} (r^{2} {\plus}p^{2} ){\minus}mx_{G} pq} \cr {{\minus}(m{\plus}m_{{22}} )vp{\plus}(m{\plus}m_{{11}} )uq{\minus}my_{G} qr{\plus}mz_{G} (p^{2} {\plus}q^{2} ){\minus}mx_{G} pr} \cr {(m_{{55}} {\minus}m_{{66}} {\minus}I_{z} {\plus}I_{y} )qr{\minus}my_{G} (pv{\minus}qu){\plus}mz_{G} (ru{\minus}pw)} \cr {(m_{{66}} {\minus}m_{{44}} {\minus}I_{x} {\plus}I_{z} )pr{\minus}mz_{G} (qw{\minus}rv){\plus}mx_{G} (pv{\minus}qu)} \cr {(m_{{44}} {\minus}m_{{55}} {\minus}I_{y} {\plus}I_{z} )qp{\minus}mx_{G} (ru{\minus}pw){\plus}my_{G} (qw{\minus}rv)} \cr } } \right]$$

2.3 Aerodynamic force model

If the airship works in an ideal circumstance and there exist no wind field during its operation, the aerodynamic force is defined as a nominal aerodynamic force. According to Liu,⁽ Reference Liu ^{24 )} it can be presented as

(9) $$\openup 3pt{\mib{F}} _{{\mib{A}} } {\equals}\left[ {\matrix{ {{1 \over 2}\rho S_{{ref}} V^{2} C_{x} cos\varphi } \cr {{1 \over 2}\rho S_{{ref}} V^{2} C_{x} sin\varphi } \cr {{1 \over 2}\rho S_{{ref}} V^{2} C_{z} } \cr {{\minus}{1 \over 2}\rho S_{{ref}} L_{{ref}} V^{2} C_{{my}} sin\varphi {\minus}\xi _{p} p} \cr {{1 \over 2}\rho S_{{ref}} L_{{ref}} V^{2} C_{{my}} cos\varphi {\minus}\xi _{q} q} \cr {{1 \over 2}\rho S_{{ref}} L_{{ref}} V^{2} C_{{mz}} {\minus}\xi _{r} r} \cr } } \right]$$

where C _x , C _z , C _my, and C _mz are the aerodynamic coefficients. It can be observed that the aerodynamic coefficient in the x-axis and the y-axis is always the same, and also the rotation aerodynamic coefficient in the x-axis and the y-axis are the same, this is because the shape of the airship is symmetrical with respect to the z-axis. φ is the angle between the x-axis and the horizontal component of the inertial airship speed, so $${\rm tan}\varphi {\equals}{u \over v}$$ . ρ is the air density. S _ref and L _ref are the reference area and the length, respectively. The presentations are S _ref =Vol ^2/3 and L _ref =Vol ^1/3, where Vol is the volume of the airship. And $$V{\equals}\sqrt {u^{2} {\plus}v^{2} {\plus}w^{2} } $$ is the inertial airship speed in BRF, ξ _p , ξ _q , ξ _r are the rotation damping coefficients. The damping coefficients are often proposed by the experience of the designer. In this paper, these are obtained by experiments.

The aerodynamic coefficients are often obtained by wind-tunnel experiment, the symmetrical structure of the airship can help to reduce the cost of the experiment because the number of the coefficients is reduced. The wind-tunnel experiment is often conducted under different angle of attacks. In this paper, following Chen’s work,⁽ Reference Chen, Zhang and Duan ^{25 )} the wind-tunnel experiment is conducted every 3° in [−180°,180°], and using linear interpolation, the coefficients of other angles are obtained. The results show that C _mz ≈10⁻⁵ is close to zero and can be ignored, and C _x , C _z , C _my are shown in Fig. 3.

Figure 3 (Colour online) Aerodynamic coefficients under different angle of attacks.

If the wind velocity is given as $$[u_{w} \quad v_{w} \quad w_{w} ]^{T} $$ , the real aerodynamic force can then be reconsidered. Define the wind-induced aerodynamic force as the difference between the real aerodynamic force and the nominal one, thus the wind-induced aerodynamic force can be presented as

(10) $$\openup 3pt{\mib{F}} _{{\mib{w}} } {\equals}\left[ {\matrix{ {{1 \over 2}\rho S_{{ref}} \left( {V_{a}^{2} {\rm cos}\varphi _{w} {\minus}V^{2} {\rm cos}\varphi } \right)C_{x} } \cr {{1 \over 2}\rho S_{{ref}} \left( {V_{a}^{2} {\rm sin}\varphi _{w} {\minus}V^{2} {\rm sin}\varphi } \right)C_{x} } \cr {{1 \over 2}\rho S_{{ref}} \left( {V_{a}^{2} {\minus}V^{2} } \right)C_{z} } \cr {{\minus}{1 \over 2}\rho S_{{ref}} L_{{ref}} \left( {V_{a}^{2} {\rm sin}\varphi _{w} {\minus}V^{2} {\rm sin}\varphi } \right)C_{{my}} } \cr {{1 \over 2}\rho S_{{ref}} L_{{ref}} \left( {V_{a}^{2} {\rm cos}\varphi _{w} {\minus}V^{2} {\rm cos}\varphi } \right)C_{{my}} } \cr {{1 \over 2}\rho S_{{ref}} L_{{ref}} \left( {V_{a}^{2} {\minus}V^{2} } \right)C_{{mz}} } \cr } } \right]$$

where $$V_{a} {\equals}\sqrt {(u{\minus}u_{w} )^{2} {\plus}(v{\minus}v_{w} )^{2} {\plus}(w{\minus}w_{w} )^{2} } $$ is the airspeed and tanφ _w =(u−u _w )/(v−v _w ). φ _w is the angle between the x-axis and the horizontal component of the airspeed.

2.4 Propeller force model

The propellers of the airship can be rotated in the vertical plane with a deflection angle μ _i , and the propulsion forces are denoted by f _i , i=1, 2, 3, 4. The detailed definition of deflection angle μ _i can be seen in Fig. 4, the thrust can rotate in the range that μ _i ∈[−180^°,180^°]. The deflection angle is 0 when the thrust is in the negative direction of the z-axis, the deflection angle is positive when the thrust rotates clockwise, and the deflection angle is negative when the thrust rotation is counterclockwise. These are the direct control variables given to the controller. In BRF, the thrust force is presented as

(11) $$\openup 3pt{\mib{F}} _{{\mib{T}} } {\equals}\left[ {\matrix{ {f_{2} {\rm sin}\mu _{2} {\minus}f_{4} {\rm sin}\mu _{4} } \cr {{\minus}f_{1} {\rm sin}\mu _{1} {\plus}f_{3} {\rm sin}\mu _{3} } \cr {{\minus}f_{1} {\rm cos}\mu _{1} {\minus}f_{2} {\rm cos}\mu _{2} {\minus}f_{3} {\rm cos}\mu _{3} {\minus}f_{4} {\rm cos}\mu _{4} } \cr {\left( {f_{4} {\rm cos}\mu _{4} {\minus}f_{2} {\rm cos}\mu _{2} } \right)R_{p} } \cr {\left( {f_{1} {\rm cos}\mu _{1} {\minus}f_{3} {\rm cos}\mu _{3} } \right)R_{p} } \cr {\left( {f_{1} {\rm sin}\mu _{1} {\plus}f_{2} {\rm sin}\mu _{2} {\plus}f_{3} {\rm sin}\mu _{3} {\plus}f_{4} {\rm sin}\mu _{4} } \right)R_{p} } \cr } } \right]$$

Figure 4 Definition of the deflection angle of the thrust.

Note in Equation (11), nonlinear couplings between f _i and μ _i are demonstrated, and control allocation must be conducted to obtain f _i and μ _i after F _T is obtained, the control allocation is often completed by choosing a control allocation matrix. The chosen of the control allocation matrix is very important in dynamic inverse control allocation to avoid a singular value. Here, the indirect control vector u _in is introduced as

(12) $${\mib{u}} _{{{\mib{in}} }} {\equals}[f_{1} {\rm sin}\mu _{1} ,\,f_{2} {\rm sin}\mu _{2} ,\,f_{3} {\rm sin}\mu _{3} ,\,f_{4} {\rm sin}\mu _{4} ,\,f_{1} {\rm cos}\mu _{1} ,\,f_{2} {\rm cos}\mu _{2} ,\,f_{3} {\rm cos}\mu _{3} ,\,f_{4} {\rm cos}\mu _{4} ]^{T} $$

This will result in a constant control allocation matrix. Thus, the inverse matrix is also a constant matrix, that is, a solution is always available in the control allocation process. Therefore, the propeller force is then decomposed as

(13) $${\mib{F}} _{{\mib{T}} } {\equals}{\mib{Bu}} _{{{\mib{in}} }} $$

the resulting control allocation matrix B is presented as

$${\mib{B}} {\equals}\left[ {\matrix{ 0 & 1 & 0 & {{\minus}1} & 0 & 0 & 0 & 0 \cr {{\minus}1} & 0 & 1 & 0 & 0 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & {{\minus}1} & {{\minus}1} & {{\minus}1} & {{\minus}1} \cr 0 & 0 & 0 & 0 & 0 & {{\minus}R_{p} } & 0 & {R_{p} } \cr 0 & 0 & 0 & 0 & {R_{p} } & 0 & {{\minus}R_{p} } & 0 \cr {{\minus}R_{p} } & {{\minus}R_{p} } & {{\minus}R_{p} } & {{\minus}R_{p} } & 0 & 0 & 0 & 0 \cr } } \right]$$

The direct control variables can then be solved by

(14) $$\eqalignno{ & f_{i} {\equals}\sqrt {\left( {f_{i} {\rm sin}\mu _{i} } \right)^{2} {\plus}\left( {f_{i} {\rm cos}\mu _{i} } \right)^{2} } \cr & \mu _{i} {\equals}{\rm arctan}{{f_{i} {\rm sin}\mu _{i} } \over {f_{i} {\rm cos}\mu _{i} }} $$

The dynamics equation can be transformed into a MIMO affine system. It is noted that in station-keeping control, the useful output is the position of the airship in ERF. Recall the dynamics model (5), the resulting affine system is

(15) $$\eqalignno{ & {\dot{\mib{x}}} {\equals}{\mib{f}} \left( {\mib{x}} \right){\plus}{\mib{g}} _{{\rm 1}} \left( {\mib{x}} \right){\mib{u}} _{{{\mib{in}} }} {\plus}{\mib{d}} \cr & {\mib{y}} {\equals}[x_{g} \quad y_{g} \quad z_{g} ]^{T} $$

where d = M ⁻¹ F _w is the wind disturbance, g ₁ ( x )= M ⁻¹ B , f ( x )= M ⁻¹( F _GB + F _A + F _I ) and y is the system output.

3.1 Explicit nonlinear MPC

For a nonlinear MIMO affine system such as the airship, it is well known that after continuously differentiating the output for a specific number of times, the control input appears in the expressions.

The relative degree is a vector, and for station-keeping control of the airship, it can be denoted as ρ _r =[ρ _r1, ρ _r2, ρ _r3]. Since the indirect control vector appears after some differentiations, the relative degree is first specified in the controller design. Before that, the following assumption is made to ensure the existence of these differentiations.

Under Assumption 1, the relative degree and $${\mib{y}} ^{{[\rho _{r} ]}} $$ are obtained by the following differentiations.⁽ Reference Slegers and Costello ^{27 )} The first derivative of y is the kinematics Equation (1).

The second derivative of y is

(16) $$\left[ {\matrix{ {\ddot{x}_{g} } \cr {\ddot{y}_{g} } \cr {\ddot{z}_{g} } \cr } } \right]{\equals}{\bf R}_{{\mib{I}} }^{{\mib{B}} } \left[ {\matrix{ {\dot{u}{\plus}qw{\minus}rv} \cr {\dot{v}{\plus}ru{\minus}pw} \cr {\dot{w}{\plus}pv{\minus}qu} \cr } } \right]$$

Then u _in appears in the expression of $$\dot{u},\,\dot{v},\,\dot{w}$$ in Equation (5). Thus, the relative degree is ρ _r1=ρ _r2=ρ _r3=2.

Substituting Equation (15) into Equation (16) gives

(17) $$\left[ {\matrix{ {\ddot{x}_{g} } \cr {\ddot{y}_{g} } \cr {\ddot{z}_{g} } \cr } } \right]{\equals}{\bf R}_{{\mib{I}} }^{{\mib{B}} } {\mib{C}} \left( {{\mib{f}} \left( {\mib{x}} \right){\plus}{\mib{g}} _{{\rm 1}} \left( {\mib{x}} \right){\mib{u}} _{{{\mib{in}} }} {\plus}{\mib{d}} } \right){\plus}{\mib{F}} _{{{\mib{qw}} }} $$

where $${\mib{F}} _{{{\mib{qw}} }} {\equals}{\bf R}_{{\mib{I}} }^{{\mib{B}} } \left[ {qw{\minus}rv,\,ru{\minus}pw,\,pv{\minus}qu} \right]^{T} $$ and $${\mib{C}} {\equals}\left[ {\matrix{ {{\mib{I}} _{{3{\times}3}} } & {{\mib{O}} _{{3{\times}3}} } \cr } } \right]$$ .

Separating u _in from Equation (17) gives

(18) $$\left[ {\matrix{ {\ddot{x}_{g} } & {\ddot{y}_{g} } & {\ddot{z}_{g} } \cr } } \right]^{T} {\equals}{\mib{L}} ({\mib{x}} ){\plus}{\bf R}_{{\mib{I}} }^{{\mib{B}} } {\mib{Cd}} {\plus}{\mib{A}} \left( {\mib{x}} \right){\mib{u}} _{{{\mib{in}} }} $$

where $${\mib{L}} \left( {\mib{x}} \right){\equals}{\bf R}_{{\mib{I}} }^{{\mib{B}} } {\rm (}{\mib{Cf}} {\rm (}{\mib{x}} {\rm )}{\plus}{\mib{F}} _{{{\mib{qw}} }} {\rm )}$$ and $${\mib{A}} \left( x \right){\equals}{\bf R}_{I}^{B} {\mib{Cg}} _{1} ({\mib{x}} )$$ .

To achieve station-keeping of the airship, we need to design a controller such that the output y tracks the reference position w _d . In MPC strategy, this can be achieved by minimising the following cost function:

(19) $${\mib{J}} {\equals}{\int}_0^{N_{p} } {\bf e}_{{\bf y}} \left( {t{\plus}\tau } \right)^{T} {\bf e}_{{\bf y}} \left( {t{\plus}\tau } \right)d\tau $$

where $${\bf e}_{{\bf y}} (t{\plus}\tau ){\equals}{\,\,\mib{{\vskip -3pt \hat}{\hskip -2pt y}}} (t{\plus}\tau ){\minus}{\,\,\mib{{\vskip -3pt \hat}{\hskip -3.5pt w}}} _{{\mib{d}} } (t{\plus}\tau )$$ , N _p is the predictive horizon and 0≤τ≤ N _p . $${\,\,\mib{{\vskip -3pt \hat}{\hskip -2pt y}}} (t{\plus}\tau )$$ , and $${\,\,\mib{{\vskip -3pt \hat}{\hskip -3.5pt w}}} _{{\mib{d}} } (t{\plus}\tau )$$ are the approximations of the output and reference position at time t+τ respectively.

In general, the cost function is optimised by MPC to obtain u _in at each sampling time. In this paper, $${\,\,\mib{{\vskip -3pt \hat}{\hskip -2pt y}}} (t{\plus}\tau )$$ and $${\,\,\mib{{\vskip -3pt \hat}{\hskip -4pt w}}} _{{\mib{d}} } (t{\plus}\tau )$$ are approximated by Taylor series expansion up to ρ _r th order

The approximation using Taylor serial expansion of $${\,\,\mib{{\vskip -3pt \hat}{\hskip -2pt y}}} (t{\plus}\tau )$$ can be presented as

(20) $${\,\,\mib{{\vskip -4pt \hat}{\hskip -2pt y}}} (t{\plus}\tau ){\equals}\left[ {\matrix{ {x_{g} (t{\plus}\tau )} \cr {y_{g} (t{\plus}\tau )} \cr {z_{g} (t{\plus}\tau )} \cr } } \right]{\equals}\left[ {\matrix{ {x_{g} (t)} \cr {y_{g} (t)} \cr {z_{g} (t)} \cr } } \right]{\plus}\tau \left[ {\matrix{ {\dot{x}_{g} (t)} \cr {\dot{y}_{g} (t)} \cr {\dot{z}_{g} (t)} \cr } } \right]{\plus}{{\tau ^{2} } \over 2}\left[ {\matrix{ {\ddot{x}_{g} (t)} \cr {\ddot{y}_{g} (t)} \cr {\ddot{z}_{g} (t)} \cr } } \right]$$

To separate the unknown time parameter τ from the approximation, Equation (20) is rearranged as

(21) $${\,\,\mib{{\vskip -4pt \hat}{\hskip -2pt y}}} (t{\plus}\tau ){\equals}{\bf \Gamma } \left( \tau \right){\mib{Y}} \left( t \right)$$

where Γ(τ)=[ T _f , Ts ], $${\mib{T}} _{{\bf s}} {\equals}diag\left[ {{{\tau ^{2} } \over 2},{{\tau ^{2} } \over 2},{{\tau ^{2} } \over 2}} \right]$$ , T _f =blkdiag ([1,τ], [1,τ], [1,τ]), it is the unkonwn matrix with the unknown τ. $${\mib{Y}} \left( t \right){\equals}[{\mib{\bar{Y}}} _{1} \left( t \right),{\mib{\bar{Y}}} _{2} \left( t \right),{\mib{\bar{Y}}} _{3} \left( t \right),\ddot{x}_{g} ,\ddot{z}_{g} ,\ddot{z}_{g} ]^{T} ,$$ $${\mib{\bar{Y}}} _{1} \left( t \right){\equals}[x_{g} ,\dot{x}_{g} ],{\mib{\bar{Y}}} _{2} \left( t \right){\equals}[y_{g} ,\dot{y}_{g} ]\!\,, {\mib{\bar{Y}}} _{3} \left( t \right){\equals}[z_{g} ,\dot{z}_{g} ].$$

Similar to $${\mib{\hat{y}}} \left( {t{\plus}\tau } \right),\,{\mib{\hat{w}}} _{{\mib{d}} } \left( {t{\plus}\tau } \right)$$ can be approximated as

(22) $${\mib{\hat{w}}} _{d} \left( {t{\plus}\tau } \right){\equals}{\bf \Gamma }\left( \tau \right){\mib{W}} _{d} \left( t \right)$$

where $${\mib{W}} _{d} \left( t \right){\equals}[{\mib{\bar{W}}} _{1} \left( t \right),{\mib{\bar{W}}} _{2} \left( t \right),{\mib{\bar{W}}} _{3} \left( t \right),\ddot{x}_{d} ,\ddot{y}_{d} ,\ddot{z}_{d} ]^{T} ,\,{\mib{\bar{W}}} _{1} \left( t \right){\equals}[x_{d} ,\dot{x}_{d} ],\,{\bf \bar{W}}_{3} \left( t \right){\equals}[z_{d} ,\dot{z}_{d} ],\,{\mib{\bar{W}}} _{2} \left( t \right){\equals}[y_{d} ,\dot{y}_{d} ],$$

Substituting Equations (21) and (22) into the Equation (19) gives

(23) $${\mib{J}} {\equals}\left[ {{\mib{Y}} \left( t \right){\minus}{\mib{W}} _{{\mib{d}} } \left( t \right)} \right]^{T} \left[ {\matrix{ {{\bf \Lambda }_{1} } & {{\bf \Lambda }_{2} } \cr {{\bf \Lambda }_{2}^{T} } & {{\bf \Lambda }_{3} } \cr } } \right]\left[ {{\mib{Y}} \left( t \right){\minus}{\mib{W}} _{{\mib{d}} } \left( t \right)} \right]$$

where $${\bf \Lambda }_{1} {\equals}{\int}_0^{N_{p} } {\mib{T}} _{{\mib{f}} } ^{T} {\mib{T}} _{{\mib{f}} } d\tau ,\,{\bf \Lambda }_{2} {\equals}{\int}_0^{N_{p} } {\mib{T}} _{{\mib{f}} } ^{T} {\bf T}_{{\mib{s}} } \left| {d\tau } \right.,\,{\bf \Lambda }_{3} {\equals}{\int}_0^{N_{p} } {\mib{T}} _{{\mib{s}} } ^{T} {\mib{T}} _{{\mib{s}} } d\tau .$$

Equation (23) can be minimised with respect to u _in (t + τ). The necessary condition for the optimality is

(24) $${{\partial {\mib{J}} } \over {\partial {\mib{u}} _{{{\mib{in}} }} }}{\equals}{\rm 2}\left[ {{\mib{Y}} \left( {\mib{t}} \right){\minus}{\mib{W}} _{{\mib{d}} } \left( {\mib{t}} \right)} \right]^{T} \left[ {\matrix{ {{\mib{\Lambda }} _{{\mib{1}} } } & {{\mib{\Lambda }} _{{\mib{2}} } } \cr {{\mib{\Lambda }} _{{\mib{2}} }^{{\mib{T}} } } & {{\mib{\Lambda }} _{{\mib{3}} } } \cr } } \right]{{\partial {\mib{Y}} } \over {\partial {\mib{u}} _{{{\mib{in}} }} }}{\equals}0$$

Recall Y (t), W _d and define $${\mib{M}} _{{\rho _{r} }} {\equals}[{\mib{\bar{Y}}} _{1} {\minus}{\mib{\bar{W}}} _{1} \quad {\mib{\bar{Y}}} _{2} {\minus}\bar{W}_{2} \quad {\mib{\bar{Y}}} _{3} {\minus}{\mib{\bar{W}}} _{3} ]^{T} $$ , then Equation (24) can be

(25) $$2[{\mib{M}} _{{\rho _{r} }}^{T} \quad (\ddot{x}_{g} {\minus}\ddot{x}_{d} ,\ddot{y}_{g} {\minus}\ddot{y}_{d} ,\ddot{z}_{g} {\minus}\ddot{z}_{d} )]\,\left[ {\matrix{ {{\bf \Lambda }_{1} } & {{\bf \Lambda }_{2} } \cr {{\bf \Lambda }_{2}^{T} } & {{\bf \Lambda }_{3} } \cr } } \right]\,[{\rm 0}^{T} ,{\rm 0}^{T} ,{\mib{A}} ({\mib{x}} )^{T} ]^{T} {\equals}0$$

It gives that

(26) $${\mib{M}} _{{\rho _{r} }}^{T} {\bf \Lambda }_{2} {\plus}(\ddot{x}_{g} {\minus}\ddot{x}_{d} ,\ddot{y}_{g} {\minus}\ddot{y}_{d} ,\ddot{z}_{g} {\minus}\ddot{z}_{d} ){\bf \Lambda }_{3} {\equals}0$$

Recall Equation (18), and note Λ ₃ is symmetrical, then the optimal control can be obtained by

(27) $${\mib{u}} _{{{\mib{in}} }} {\bf (}{\mib{t}} {\plus}\tau {\bf )}^{{\bf {\asterisk}}} {\equals}{\minus}{\mib{A}} \left( {\mib{x}} \right)^{{\plus}} \left[ {{\mib{KM}} _{{\rho _{r} }} {\plus}{\mib{M}} _{1} {\plus}{\mib{R}} _{I}^{B} {\mib{Cd}} } \right]$$

where A ( x ) ⁺ is the Moore–Penrose pseudoinverse of $${\mib{A}} \left( {\mib{x}} \right),{\mib{M}} _{1} {\equals} {\mib{L}} ({\bf x}){\minus}[\ddot{x}_{d} \quad \ddot{y}_{d} \quad \ddot{z}_{d} ]^{T} ,\,{\mib{K}} {\equals}{\bf \Lambda }_{3}^{{{\minus}1}} {\bf \Lambda }_{2}^{T} .$$

After solving the nonlinear equation (27), we can obtain the optimal control u _in ^*(t+τ). In this paper, only the current control in the control profile is implemented. The explicit solution is u _in ^*= u _in (t+τ), for τ=0.

In addition, it can be noted that in Equation (27), the disturbance d is still unknown in the control law, so an observer is designed in the next section for disturbance rejection. After that, d is replaced by its estimation value $${\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}} $$ , then the influence of wind field is compensated. Therefore, the station-keeping control law is

(28) $${\mib{u}} _{{{\mib{in}} }} {\equals}{\minus}{\mib{A}} \left( {\mib{x}} \right)^{{\plus}} \left[ {{\mib{KM}} _{{\rho _{r} }} {\plus}{\mib{M}} _{1} {\plus}{\mib{R}} _{I}^{B} {\mib{C }{\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}} \right]$$

The direct control variables can then be calculated by Equation (14).

3.2 Nonlinear DOB

As disscused in the previous section, it is clear that $${\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}} $$ is indispensable in controller design. In this section, a nonlinear DOB is proposed to estimate the disturbance $${\mib{d}} $$ .

The nonlinear DOB for MIMO affine system (15) is designed as

(29) $$\eqalignno{ & {\dot{\mib{z}}} {\equals}{\minus}{\mib{l}} \left( {\mib{x}} \right){\mib{z}} {\minus}{\mib{l}} \left( {\mib{x}} \right)\left( {{\mib{p}} \left( {\mib{x}} \right){\plus}{\mib{f}} \left( {\mib{x}} \right){\plus}{\mib{g}} _{{\rm 1}} \left( {\mib{x}} \right){\mib{u}} _{{{\mib{in}} }} } \right) \cr & {\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}} {\equals}{\mib{z}} {\plus}{\mib{p}} \left( {\mib{x}} \right) $$

where z is the state vector of DOB, l ( x ) is the observer gain matrix, and p ( x ) is the nonlinear function to be designed. The observer gain matrix can be determined by

(30) $${\mib{l}} \left( {\mib{x}} \right){\equals}{{\partial {\mib{p}} \left( {\mib{x}} \right)} \over {\partial {\mib{x}} }}$$

(31) $${\dot{\mib{e}}} {\equals} \,\, {{ \vskip -7pt \dot} \hskip 2pt{{\mib {{{\vskip -5pt \hat}}{\hskip -3pt d}}}}} &{\minus}{\dot{\mib{d}}}$$

Substituting Equation (29) into Equation (31) gives

(32) $$\eqalignno{ {\dot{\mib{e}}} {\equals} {\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}} &{\equals}{\dot{\mib{z}}} {\plus}{{\partial {\mib{p}} \left( {\mib{x}} \right)} \over {\partial {\mib{x}} }}{\dot{\mib{x}}} \cr &{\equals}{\minus}{\mib{l}} \left( {\mib{x}} \right){\rm (}{\mib{z}} {\minus}{\dot{\mib{x}}} {\rm )}{\minus}{\mib{l}} \left( {\mib{x}} \right)\left( {{\mib{p}} \left( {\mib{x}} \right){\rm {\plus}}{\mib{f}} \left( {\mib{x}} \right){\rm {\plus}}{\mib{g}} _{{\rm 1}} \left( {\mib{x}} \right){\mib{u}} _{{{\mib{in}} }} } \right) $$

Substituting Equations (15) and (29) into Equation (32) gives

(33) $$\eqalignno{ \dot{e}&{\equals} {\minus}{\mib{l}} \left( {\mib{x}} \right){\rm (}{\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}} {\minus}{\mib{p}} \left( {\mib{x}} \right){\minus}{\mib{f}} \left( {\mib{x}} \right){\minus}{\mib{g}} _{{\mib{1}} } \left( {\mib{x}} \right){\mib{u}} _{{{\mib{in}} }} {\minus}{\mib{d)}} {\minus}{\mib{l}} \left( {\mib{x}} \right)\left( {{\mib{p}} \left( {\mib{x}} \right){\mib{{\plus}f}} \left( {\mib{x}} \right){\mib{{\plus}g}} _{{\mib{1}} } \left( {\mib{x}} \right){\mib{u}} _{{{\mib{in}} }} } \right) \cr &{\rm {\equals}}{\minus}{\mib{l}} \left( {\mib{x}} \right){\mib{e}} $$

This demonstrates that $${\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}} $$ approaches d exponentially if l ( x ) is Hurwitz.

3.3 Actuator saturation suppression based on TD

The amplitude of the propulsion force is severely limited in the realistic scenario, while the analytical control law (28) can handle no saturation. Therefore, to preserve the analytical control law and decrease the online computation burden, the TD employing time-varying design parameter is proposed.

TD is first proposed by Han,⁽ Reference Han and Wang ^{28 )} and it uses a numerical method based on integrity to track the differentiation of the signal. In handling the saturation, TD always calculate a smooth govern path within the saturation. The main advantage of TD is the control of rapidness and strong robustness to the disturbance. The most common used TD⁽ Reference Han ^{20 )} can be presented as

(34) $$\left\{ {\matrix{ {\dot{x}_{1} {\equals}x_{2} } \cr {\dot{x}_{2} {\equals}fhan\left( {x_{1} ,x_{2} ,R_{{TD}} } \right)} \cr } } \right.$$

where x ₁ and x ₂ are the states of the TD, they represent the position and velocity when respecting to the station-keeping control. R _TD is the design parameter that represents the acceleration of the transition process. It is the only parameter that determines the transition time.⁽ Reference Han and Yuan ^{29 )} fhan denotes the optimal control synthesis function. This function can avoid the chattering in steady state caused by the non-smooth structure of TD. It can be presented as

(35) $$fhan\left( {x_{1} ,x_{2} ,R_{{TD}} } \right){\equals}\left\{ {\matrix{ {R_{{TD}} a\quad \quad \,\mid\,a\,\mid\leq d_{t} } \cr {{\minus}R_{{TD}} {\rm sign}(a)\quad \quad \,\mid\,a\,\mid\gt\,d_{t} } \cr } } \right.$$

where d _t =R _TD h ₀, h ₀ is the sampling time interval, and a is presented as

(36) $$a{\equals}\left\{ {\matrix{ {x_{2} {\plus}c\,/\,h_{0} \quad \quad \,\mid\,c\,\mid\leq d_{{t0}} } \cr {x_{2} {\plus}sign(c)(a_{0} {\minus}d_{t} )\,/\,2\quad \quad \,\mid\,c\,\mid\,\gt\,d_{{t0}} } \cr } } \right.$$

where d _t0=d _t h ₀, c=x ₁ + x ₂ h ₀ and $$a_{0} {\equals}\sqrt {d_{t}^{2} {\plus}8R_{{TD}} \,\mid\,c\,\mid\,^{2} } $$ .

To better understanding the basic theory and tracking performance of TD in path planning, the forward flight of the airship is taken as an example. The flight acceleration, transition time and reference position in the x-axis are defined as R _TD , T _0, and x _d , respectively. The acceleration has the shape of the square wave, and velocity having the triangle wave. To illustrate it in detail, in the first half of the process, the airship moves with a constant acceleration R _TD , reaching the position $${{x_{d} } \over 2}$$ and velocity $${{T_{0} } \over 2}R_{{TD}} $$ . While in the latter half, the airship moves with the acceleration −R _TD , reaching the position x _d and velocity 0.

Figure 5 shows the position, velocity, and acceleration under different reference positions 10 m, 20 m, and 30 m, respectively, with R _TD =0.5 m/s ². It can be seen clearly that once the acceleration is determined, the controller will always achieve the same performance in handling the saturation, no matter how the reference positions change. Different reference positions will only lead to different transition times.

Figure 5 (Colour online) Path planning of TD.

In this paper, the main motivation of saturation handling strategy of TD is to choose a proper R _TD so that the direct control variables will always be within its boundary. We propose a time-varying R _TD choosing method based on the dynamics model of the airship via the relationship between the maximum flight acceleration and maximum propulsion force. Before that, the following assumption is made to simplify and decouple the dynamics model.

Define * _max as the maximum value of *. Recall the dynamics model (5) and replace d by $${\mib{\hat{d}}} $$ . The maximum acceleration of the airship can be expressed by

(37) $${\mib{M}} {\dot{\mib{x}}}_{{max}} {\equals}{\mib{F}} _{{{\mib{T}} max}} {\plus}{\mib{F}} _{{\mib{A}} } {\plus}{\mib{F}} _{{{\mib{GB}} }} {\plus}{\mib{F}} _{{\bf I}} {\plus}{\mib{M{\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}}} $$

Under Assumption 4, only the first two rows of Equation (37) are considered. $$\dot{p},\,\dot{q}\,{\rm and}\,\dot{r}$$ are not taken into consideration because of their small values. Define d _x and d _y as the disturbances in x- and y-axes, $$\hat{d}_{x} \,{\rm and}\,\hat{d}_{y} $$ as their estimations, respectively. Then the time-varying accelerations in x- and y-axes are

(38) $$\eqalignno{ & \left( {m_{0} {\plus}m_{{11}} } \right)\left( {\dot{u}_{{{\rm max}}} {\minus}\hat{d}x} \right){\equals}f_{{{\rm max}}} \left( {{\rm sin}\mu _{2} {\minus}{\rm sin}\mu _{4} } \right){\plus}{\mib{F}} _{{\mib{A}} } \left( 1 \right){\plus}{\mib{F}} _{{\bf I}} \left( 1 \right){\plus}{\mib{F}} _{{{\mib{GB}} }} \left( 1 \right) \cr & \left( {m_{0} {\plus}m_{{22}} } \right)\left( {\dot{v}_{{{\rm max}}} {\minus}\hat{d}y} \right){\equals}f_{{{\rm max}}} \left( {{\rm sin}\mu _{1} {\minus}{\rm sin}\mu _{3} } \right){\plus}{\mib{F}} _{{\mib{A}} } \left( 2 \right){\plus}{\mib{F}} _{{\mib{I}} } \left( 2 \right){\plus}{\mib{F}} _{{{\mib{GB}} }} \left( 2 \right) $$

Thus the maximum accelerations in x and y axes are

(39) $$\eqalignno{ & \dot{u}_{{max}} {\equals}{{f_{{max}} \left( {sin\mu _{2} {\minus}sin\mu _{4} } \right){\plus}{\mib{F}} _{{\mib{A}} } \left( 1 \right){\plus}{\mib{F}} _{{\mib{I}} } (1){\plus}{\mib{F}} _{{{\mib{GB}} }} \left( 1 \right)} \over {m0{\plus}m_{{11}} }}{\plus}\hat{d}_{x} \cr & \dot{v}_{{max}} {\equals}{{f_{{max}} \left( {sin\mu _{1} {\minus}sin\mu _{3} } \right){\plus}{\mib{F}} _{{\mib{A}} } \left( 2 \right){\plus}{\mib{F}} _{{\mib{I}} } {\plus}{\mib{F}} _{{{\mib{GB}} }} \left( {\bf 2} \right)\left( 2 \right)} \over {m_{0} {\plus}m_{{22}} }}{\plus}\hat{d}_{y} $$

Then $$\dot{u}_{{{\rm max}}} $$ and $$\dot{v}_{{{\rm max}}} $$ are used as the time-varying design parameters in Equation (34).

4.1 Controllerstructure

Based on the previous section, the station-keeping controller consists of three modules: explicit NMPC, nonlinear DOB, and TD. Its structure can be seen in Fig. 6.

Figure 6 Structure of the station-keeping controller.

In Fig. 6, w _d is the reference position of the station-keeping, and $${\mib{w'}} _{{\mib{d}} } {\equals}\left[ {\matrix{ {x'_{{d}} } & {y'_{{d}} } & {z'_{{d}} } \cr } } \right]^{T} $$ is the command position planned by TD, and the other parameters remain the same as previous.

4.2 Stability analysis

From the structure of the controller, it can be seen that TD has no effect on the DOB-based NMPC controller. It only plans a command path in the station-keeping, so the stability relies on DOB and explicit NMPC. Now we focus on the stability of the DOB-based NMPC controller, and the stability analysis is given as follows.

The position tracking error of station-keeping control is defined as

(40) $${\mib{z}} _{{\mib{p}} } ^{0} {\equals}\left[ {\matrix{ {x_{g} {\minus}x{\hbox{'}} _{d} } & {y_{g} {\minus}y \hbox{'} _{d} } & {z_{g} {\minus}z \hbox{'}_{d} } \cr } } \right]^{T} $$

Its first derivative is

(41) $${{\dot{\mib z}}} _{p}^{0} {\equals}{\mib{z}} _{{\mib{p}} } ^{1} {\equals}\left[ {\matrix{ {\dot{x}_{g} {\minus}\dot{x}\hbox{'}_{d} } & {\dot{y}_{g} {\minus}\dot{y} \hbox{'} _{d} } & {\dot{z}_{g} {\minus}\dot{z} \hbox{'} _{d} } \cr } } \right]^{T} $$

And the second derivative is

(42) $${{\dot{\mib z}}} _{{\mib{p}} }^{{\rm 1}} {\equals}\left[ {\matrix{ {\ddot{x}_{g} {\minus}\ddot{x} \hbox{'} _{d} } & {\ddot{y}_{g} {\minus}\ddot{y} \hbox{'} _{d} } & {\ddot{z}_{g} {\minus}\ddot{z} \hbox{'} _{d} } \cr } } \right]^{T} $$

Define

(43) $${\mib{z}} _{{\mib{p}} }^{{\rm 2}} {\equals}\left[ {\matrix{ {\ddot{x}_{g} } & {\ddot{y}_{g} } & {\ddot{z}_{g} } \cr } } \right]^{T} $$

Recall $$M_{{\rho _{r} }} $$ , M ₁. Substituting Equations (18) and (28) into Equation (42) gives

(44) $$\eqalignno{ {{\dot{\mib z}}} _{p}^{1} & {\equals} {\mib{M}} _{1} {\plus}{\mib{R}} _{I}^{B} {\mib{Cd}} {\plus}{\mib{A}} \left( {\mib{x}} \right)\left[ {{\minus}{\mib{A}} \left( {\mib{x}} \right)^{{\plus}} \left( {{\mib{KM}} _{{\rho _{r} }} {\plus}{\mib{M}} _{1} {\plus}{\mib{R}} _{{\mib{I}} }^{{\mib{B}} } {\mib{C{\,\,\mib{{\vskip -5pt \hat}{\hskip -3pt d}}}}} } \right)} \right] \cr &{\equals}{\minus}{\mib{KM}} _{{{\mib{\rho }} _{{\mib{r}} } }} {\mib{{\plus}R}} _{{\mib{I}} }^{{\mib{B}} } {\mib{Ce}} $$

where K =blkdiag ([k ₁₁, k ₁₂], [k ₂₁, k ₂₂], [k ₃₁, k ₃₂])

Recall e and $${\mib{z}} _{{\mib{p}} } ^{{\rm 2}} {{.\,\dot{\mib z}}} _{{\mib{p}} } ^{1} $$ can be expressed as

(45) $${{\dot{\mib z}}} _{{\mib{p}} }^{{\mib{1}} } {\equals}{\mib{z}} _{{\mib{p}} } ^{2} {\plus}{\mib{R}} _{{\mib{I}} }^{{\mib{B}} } {\mib{Ce}} $$

Recall the presentation of K and $${\mib{M}} _{{\rho _{r} }} $$ , and rearrange them. Recall Equations (41), (42), and (45). The tracking error of the control system is

(46) $$\left[ {\matrix{ {{\dot {\mib z}} _{{\mib{p}} }^{{\rm 0}} } \cr {{\dot {\mib z}} _{{\mib{p}} }^{{\rm 1}} } \cr } } \right]{\equals}\left[ {\matrix{ {{\rm 0}_{{3{\times}3}} } & {{\mib{I}} _{{3{\times}3}} } \cr {{\mib{K}} _{1} } & {{\mib{K}} _{2} } \cr } } \right]\left[ {\matrix{ {{\mib{z}} _{{\mib{p}} }^{{\rm 0}} } \cr {{\mib{z}} _{{\mib{p}} }^{{\rm 1}} } \cr } } \right]{\plus}\left[ {\matrix{ {{\rm 0}_{{3{\times}3}} } \cr {\epsilon} \cr } } \right]$$

where K _i =diag (k _1i , k _2i , k _3i ), i=1, 2,3, and $${\epsilon}{\equals}{\bf{R}} _{{\mib{I}} }^{{\mib{B}} } {\mib{Ce}} $$ .

The tracking error of the station-keeping controller can then be expressed as

(47) $${{\dot{\mib z}}} _{{\mib{p}} } {\equals}{\mib{f}} _{{\mib{A}} } {\mib{z}} _{{\mib{p}} } {\plus}{\epsilon}$$

From Lemma 1, we know that ϵ will converge to zeros if l ( x ) is Hurwitz. Therefore, the tracking error can be reduced to a linear system $${\dot{\mib{z}}} _{{\mib{p}} } {\equals}{\mib{f}} _{{\mib{A}} } {\mib{z}} _{{\mib{p}} } $$ , and its global asymptotic stability can be guaranteed by properly choosing K such that $${{\mib{f}}_{{\mib{A}}$$ is Hurwitz. In this paper, it can be achieved by adjusting the predictive horizon such that $${{\mib{f}}_{{\mib{A}}$$ is Hurwitz.

5.1 Parameters analysis of nonlinear DOB and explicit NMPC

In this section, the parameters l and N _p are investigated. Some comparisons on the performance are carried out by applying different l and N _p , and comments are given. Then the parameters with better performance are chosen for the next section to illustrate the effectiveness of the proposed station-keeping controller.

The DOB is first examined. The predictive horizon is given as N _p =1 s. The estimation performances with different ls (1, 10, and 80) are given in Fig. 7. From lemma 1, we know that $$\hat{d}$$ approaches d exponentially if l ( x ) is Hurwitz. It means that the lager gain in DOB leads to faster converge rate. It is consistent with the simulation result in Fig. 7 (see when l=1 and l=10). When l increases, overshoot can be observed (see when l=80), but the converge rate does not change too much. Therefore, the tradeoff between overshoot and converge rate has to be made. Thus, l=10 is choosen.

Figure 7 (Colour online) Disturbance observer performance under different l.

To verify remark 1, a random wind and a cosine wind field are also implemented in the simulation, and the parameters are l=10 and N _p =1 s. The velocity of the random wind is given as [10 + 0.5ξ,10 + 0.5ξ,0] m/s, where ξ is a random value in [−1,1]. The cosine wind is given as [10 + sin(0.2t),10 + cos(0.2t),0] m/s. Figs. 8 and 9 show the results of the estimations. It can be seen that the estimations are both converged to the wind disturbance. The remark 1 is verified.

Figure 8 (Colour online) Disturbance observer performance of cosine wind.

The explicit NMPC is then examined, and l=10 are fixed. The position responses under different N _p s(1, 2, and 3 s) are given in Fig. 10. Figure 10 shows that the larger N _p is, the larger the transition time of the system is. And larger predictive horizon also leads to larger overshoot( see in the z-axis). Thus, N _p =1 s is chosen (Fig. 10).

Figure 9 (Colour online) Positions under different predictive horizons.

Figure 10 (Colour online) Positions under different predictive horizons.

5.2 The controller performance

In this section, taking the wind disturbance and thrust saturation into consideration, we implement our method to the airship. According to the previous section, l=10 and N _p =1 s are fixed in the following simulations.

The nonlinear DOB is illustrated with its ability for disturbance rejection. The states of the airship with or without nonlinear DOB are given in Fig. 11. The result shows that the controller with DOB performs much better. Specifically, the DOB can restrain the states’ variation degree, especially rolling angular velocity (its amplitude decreases significantly from 7.8 to 0.8 radians/s). It is able to drive the airship to the reference position (otherwise, the position of the airship cannot converge to the reference position). Since the disturbance is mainly caused by wind field, we can conclude that wind field (disturbance) has huge effects on the station-keeping performance. In addition, DOB has the ability of disturbance rejection, thus, the station-keeping performance has been significantly improved. Those interested in disturbance estimation performance can refer to Fig. 7 when l=10.

Figure 11 (Colour online) Velocities, angular velocities, and positions with or without nonlinear DOB.

Finally, the ability to handle the thrust saturation of TD is studied, and the propulsion forces and deflection angles are given in Fig. 12. The propulsion forces and deflection angles without using TD are given in Fig. 13 as a comparison. The two figures demonstrate that when the controller equipped with TD, the maximum value of propulsion forces decreases from 3685N (in propeller 4 at 5.8 s) to 502N (in propeller 3 at time 13.2 s). The propulsion force slightly exceeds the given maximum value. From the comparison, it can be concluded that the TD can handle thrust saturation efficiently.

Figure 12 Propulsion forces and deflection angles with TD.

Figure 13 Propulsion forces and deflection angles without TD.