# Nonlinear Reduced-Order Observer-Based Predictive Control for Diving of an Autonomous Underwater Vehicle

Authors: Xuliang Yao; Guangyi Yang; Yu Peng
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

The attitude control and depth tracking issue of autonomous underwater vehicle (AUV) are addressed in this paper. By introducing a nonsingular coordinate transformation, a novel nonlinear reduced-order observer (NROO) is presented to achieve an accurate estimation of AUV’s state variables. A discrete-time model predictive control with nonlinear model online linearization (MPC-NMOL) is applied to enhance the attitude control and depth tracking performance of AUV considering the wave disturbance near surface. In AUV longitudinal control simulation, the comparisons have been presented between NROO and full-order observer (FOO) and also between MPC-NMOL and traditional NMPC. Simulation results show the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Autonomous underwater vehicle (AUV) is an important tool for ocean exploitation. Due to its strong autonomous ability, AUV can complete the mission given by human and has become widely used in military and scientific research [1].Underactuated AUV is a kind of underwater vehicles whose actual control inputs are less than degrees of freedom (DOF) [2]. When AUV navigates near surface, diving control is one of the key research fields; a variety of methods have been proposed for AUV in vertical plane. Yan et el. [3] proposed backstepping technology and IFTSMC for surge velocity and depth control, respectively. However, the coupled interaction between the two controllers was not considered. An improved dynamic fuzzy SMC is established to enhance the depth tracking control performance in [4]. Adhami-Mirhosseini et al. [5] proposed a nonlinear dynamic controller combined with Fourier series expansion and pseudospectral ideas to achieve the depth tracking objective. In [6], an analytical SMC method is used to obtain the time optimal trajectory tracking, and the effectiveness of the proposed controller was verified via simulations. The fuzzy feedback linearization methods studied in [7] abandoned two general assumptions and used the nonlinear dynamics of depth motion directly. However, discrete-time model was not taken into consideration, and there is lack of simulation verification of the proposed algorithm. Hsu and Liu [8] investigated a steady-state depth errors problem caused by center of gravity change and considered the effects of gravity and buoyancy in AUV dynamics. By adding a switching PI controller in external-loop, the steady-state depth errors were eliminated finally. A static output feedback controller was studied in [9] to complete the diving task; however a linearized model was used.Model predictive control can calculate a sequence of future control signals in advance, in order to track a given future reference. Moreover, predictive control allows the integration of system constraints in the controller design and minimizes a cost function defined over a prediction horizon, ensuring near-optimal performance of the system. The predictive control of AUV motion is one of the promising research areas in marine control systems field. In [10], visual servoing MPC is used to bring current-influenced AUV aperiodic algorithm-generated control loop closures. A cost function minimization method is used in [11] with MPC, called least squares, whose advantage is the real-time execution ability to optimize sawtooth paths for an underwater glider. A reduced dynamical model proposed in [12] is used to control AUV in vertical and horizontal planes, respectively. The MPC is formulated as two independent optimization parts. However, by using the same initial states for both parts of linearization and neglecting the weakly coupled states, the control optimization differs from the actual situation a lot. In [13], a sampling based MPC is proposed to generate the control sequence with constraints effectively. These methods have a good performance in the field of MPC use for AUV, but they do not deal with wave disturbances. The method in [14] used MPC to predict future disturbances and counteract the wave disturbances. A linear wave theory solver is employed by estimator to approximate the fluid dynamics, but constraints were not considered. NMPC is an extension of MPC by using nonlinear systems and constraints; applications of NMPC for AUV are less frequent. In [15], NMPC is used to design the kinematic loop in a hierarchical structure for AUV DP control.Because of fewer actuators and sensors, there is an advantage for underactuated AUV to simplify the design and reduce costs. However, it brings highly coupled nonlinear models, which increases the difficulty of controller design. In general, some state variables of AUV are unmeasured due to a lack of sensors or sensor faults. So an observer-based controller becomes significant. A dual observer in [16], which combines a state and a perturbation observer, aims to solve the problem about being sensitive to external disturbance. Zhang et al. [17] considered the elevator angle constraints and proposed a controller based on MPC with artificial bee colony algorithm, and a classical linear state observer is used. In [18], a novel discrete-time proportional and integral observer is used to estimate the states, output, input, and disturbance together. Zhang et al. [19] proposed a new reduced-order observer with multiconstrained thoughts by using specific system decomposition. However, these results only use linear models.In contrast, most practical systems are nonlinear and therefore nonlinear models are required. Kinsey et al. [20] proposed a nonlinear observer based on dynamic model of AUV, which is used to estimate the vehicle’s velocity. However, in their model, the coupling terms are neglected and there was only one degree of freedom. Mahapatra et al. [21] investigated a nonlinear H∞ feedback control algorithm with a nonlinear state observer, but the error dynamics stability was not considered in observer design. Zhang et al. [22] extended an adaptive state observer to a class of nonlinear systems. However, due to the selected special Lyapunov matrix, there is a reduction of the accuracy of state estimation.At the same time, MPC need state variables to achieve desired output tracking in the minimization of cost function, so the role of observer-based MPC has a prominent position in practical implementation. However, there is a situation that some states variables of underactuated AUV can be measured. So it is not necessary to estimate all the variables. Compared with FOO, reduced-order observer (ROO) has simplified structure and good estimation performance [23, 24], which can be designed to estimate only the unmeasurable states. Recently, ROO designs used in MPC are exploited in both linear [25, 26] and nonlinear systems [27–30]. However, few research results are reported on NROO-based MPC for AUV. So the motivation of this paper is the aforementioned problems and deficiencies, and the aim is to design a NROO-based MPC with nonlinear AUV kinematics and dynamics model.In this paper, a model predictive controller with nonlinear model online-linearization based on NROO is proposed to solve the attitude control and depth tracking issues. The design process is as follows. First, the nonlinear dynamics in AUV vertical plane are described and Euler approximation is used to discrete the proposed model. Second, the nonlinear reduced-order observer is designed according to the dynamics above, to estimate the unmeasured state variables of AUV. By solving a convex optimization problem and using LMI technique, the estimation performance is enhanced and the wave disturbance influence is restricted. Third, a model predictive control based on model online-linearization is presented, in which the controller design process is independent of the NROO and it is convenient for the following simulation. Finally, simulation results show the good robustness and efficiency of the proposed methods, which can be applied to AUV attitude control and depth tracking field.The rest of the paper is organized as follows. Section2 described the AUV vertical plane dynamics. NROO design is proposed in Section 3. Section 4 introduces the MPC-NMOL strategy. Section 5 gives the simulation results of AUV to verify the effectiveness of the proposed approach.


## 2. AUV Vertical Plane Dynamics

The vertical plane dynamic and kinematic model of AUV are described in this section, and the following assumptions are considered.Assumption 1.
( 1 ) AUV navigates at a constant speed u=Uc.
( 2 ) To analyze the diving control, sway, yaw, and roll motions can be ignored, and that means y=0, v=0, and ϕ=ψ=0.
( 3 ) Center of buoyancy stays constant with no change made to external hull, but the position of center of gravity is variable.
( 4 ) The hydrodynamics drag terms whose order is higher than two are neglected.
( 5 ) Pitch angle θ has small perturbations at zero.
Therefore, the dynamic and kinematic model of AUV can be expressed in the nonlinear forms as follows:(1)mω˙-uq-xgq˙=ΣZext,Iyq˙+mzgωq-xgω˙+xguq=ΣMext,z˙=-uθ+ω,θ˙=q,where(2)ΣZext=ZHS+Zωωωω+Zω˙ω˙+Zq˙q˙+Zuquq+Zuωuω+Zδu2δ+Zwave,ΣMext=MHS+Mqqqq+Mω˙ω˙+Mq˙q˙+Muquq+Muωuω+Mδu2δ+Mwave,ZHS=W-Bcos⁡θ,MHS=zbB-zgWθ+xbB-xgW.Iy, m are the moment of inertia and the mass of AUV, respectively, (xg,yg,zg) and (xb,yb,zb) are the coordinates of AUV’s center of gravity and buoyancy, u, ω, and q are the surge, heave, and pitch angle velocity in body-fixed frame, respectively, z are depth in earth-fixed frame, δ is the elevator deflection angle, Zωωωω, Zω˙, Zq˙, Zuq, Zuω, Zδ, Mqqqq, Mω˙, Mq˙, Muq, Muω, and Mδ are the hydrodynamics coefficients of AUV, and Zwave, Mwave are the wave force and moment.The nonlinear model (1) can be simply expressed as the following system:(3)x˙ct=Acxct+gct,xt+Bcut+Dcwt,yt=Cxct,where(4)Ac=M1-1n11n1200n21n220n24100-u0100,Bc=M1-1ZδMδ00,Dc=M1-1I2×202×2,gct,xt=M1-1Zωωωω+W-B-mzgωq+Mqqqq-xgW-xbB00,M1=m11m1200m21m220000100001,C=00100001,m11=m-Zω˙,m12=-Zq˙-mxg,m21=-mxg-Mω˙,m22=Iy-Mq˙,n11=ZuωUc,n12=Zuq+mUc,n21=MuωUc,n22=Muq-mxgUc,n24=-zgW-zbB.Due to the digital controller implemented by using zero-order hold, a discrete-time vector-form nonlinear state equation is obtained by Euler approximation [31].(5)xk+1=Axk+gk,xk+Buk+Dwk,yk=Cxk,where x(k)=[x1,x2,x3,x4]T∈Rn, x1=ω, x2=q, x3=z, x4=θ, Ts is the sampling period, A=(TsAc+I), B=TsBc, D=TsDc, w(k)=[ZwaveMwave]∈Rd, u(k)∈Rm is the input δ, gk,xk=Tsgct,xt is a nonlinear-term vector, and yk∈Rp denotes the output vector. In this paper, n=4, d=2, m=1, p=2, C is a full-rank matrix, and rank(C)=2, and we assume that pair (A,C) is observable.Remark 2.
In this paper, the disturbance due to surface waves is considered as the main disturbance type and written aswk. Sea state is used to classify the sea conditions and expressed by wind speed and wave height. Wave height is affected by wind, and the wave profile is affected by water depth. For surface waves modelling, we assume the long crested and unidirectional sea state. The superposition principle is used to complete an irregular long crested wave sea state simulation. The single parameter Pierson-Moskowitz spectrum [32] is used in this model and its spectrum density is defined as follows:(6)Sω=8.1×10-3g2ω-5exp⁡-3.11Hs2ω-4,where Hs is the significant wave height (m) and ω- is wave frequency. Assuming the AUV body is small enough compared with the incoming wave and the wave speed is large enough in relation to the hull diameter D0, the transient force and moment are found by using Morison’s equation [33] and integration along the length of the hull L0 at each step is obtained below:(7)Zwave=∫L0CdρD02ωω-ω2+CmρπD204ω˙ω-ω˙2dx,Mwave=∫L0CdρD02ωω-ω2+CmρπD204ω˙ω-ω˙2xdx,where ρ is fluid density, Cd is the drag coefficient, and Cm is the added mass coefficient.In (7), (8)ωω=∑i=1Nωeiαiexp⁡-kizsin⁡ωeit+θi,ω˙ω=-∑i=1Nωei2αiexp⁡-kizcos⁡ωeit+θi,where ωei is encountering frequency, αi=2S(ω-i)dω-, N is the equal interval frequency waveband and chosen in this article to discretize the wave spectrum, ki is wave number, and θi is random phase shift (0<θi<2π) in relation to each frequency.Assumption 3.
g ( k , x ( k ) ) is a Lipschitz function [34] vector:(9)gk,x-1-gk,x-22≤Lgx-1-x-22which is assumed to satisfy the following conditions:(1)
The effect ofZHS is approximated to zero.(2)
xg≈xb=0,where Lg>0 is a Lipschitz constant and ·2 denotes the Euclidean norm. Moreover, g(0,k)=0.Remark 4.
This paper proposed MPC-NMOL design method based on NROO for a nonlinear Lipschitz model of AUV. The NROO is independent of the MPC-NMOL, which is constructed to guarantee the stability and performance of AUV diving motion. Figure1 shows the structure of the proposed system.Figure 1
Block diagram of the proposed system.


## 3. Nonlinear Reduced-Order Observer Design

### 3.1. NROO Algorithm Design

Consider thatC is full-row rank, and there always exists a matrix C⊥∈R(n-p)×n such that C⊥C∈Rn×n is nonsingular (C⊥ can be given as an orthogonal basis of the null-space of C). Then under the nonsingular coordinate transformation z(k)=Tx(k), in which T=C⊥C is used, the dynamics (5) can be expressed into the form(10)z1k+1z2k+1=A11A12A21A22z1kz2k+C⊥Cgk,T-1zk+B1B2uk+D1D2wk,yk=0p×n-pIpz1kz2k,where z1(k)∈Rn-p and z2(k)∈Rp are new state vectors, Ip is identity matrix whose dimension is p×p, and (11)TAT-1=A11A12A21A22,TB=B1B2,TD=D1D2,CT-1=0p×n-pIp,zk=Txk=z1kz2k,and then (12)z2k=yk,z1k+1=A11z1k+A12yk+C⊥gk,T-1zk+B1uk+D1wk,yk+1=A21z1k+A22yk+Cgk,T-1zk+B2uk+D2wk.By using the hypothesis inputμ(k) and output ρ(k): (13)μk=A12yk+B1uk,ρk=yk+1-A22yk-B2ukwe get (14)z1k+1=A11z1k+C⊥gk,T-1zk+μk+D1wk,ρk=A21z1k+Cgk,T-1zk+D2wk.From the dynamics (14), NROO is constructed as (15)z^1k+1=A11z^1k+C⊥gk,T-1z^k+μk-Gρ^k-ρk,(16)ρ^k=A21z^1k+Cgk,T-1z^k,where z^1(k)∈Rn-p is the observer state of the reduced-order system (14), ρ^(k)∈Rp is the output of reduced-order observer, and G∈R(n-p)×p is the observer gain matrix.Denote(17)ek=z^1k-z1k,Gk=gk,T-1z^k-gk,T-1zkand then, from (14)–(16), the state estimation error dynamics are described by the following:(18)ek+1=A11-GA21ek+C⊥-GCGk+GD2-D1wk.Equation (18) is rewritten as the following form:(19)ek+1=A-1-GA-2ek+GC2-C1Gk+GD2-D1wk, where A-1=A11, A-2=A21, C1=-C⊥, and C2=-C.Theorem5 gives a constrained NROO design, whose state estimation performance is specified via a H∞ performance index.Theorem 5.
Let a prescribedH∞ performance level γ1>0, and if there exist a symmetric positive definite matrix P-∈R(n-p)×(n-p), a matrix Y-∈R(n-p)×p, and a ε1>0 such that the following condition holds:(20)-P-P-A-1-Y-A-2Y-C2-P-C1Y-D2-P-D10∗-P-+ε1Lg2In-p00In-p∗∗-ε1λminTTTIn00∗∗∗-γ1Id0∗∗∗∗-γ1In-p≤0,where Y-=P-G and ∗ denotes the symmetric elements in a matrix, then the error dynamics (18) satisfy H∞ performance index e(k)2≤γ1w(k)2, and the NROO gain matrix can be obtained by G=P--1Y-.Proof.
Choose the following Lyapunov function:(21)Vk=eTkP-ek.
According to error dynamics (18) and (21), the difference ΔV(k) is (22)ΔVk=Vk+1-Vk=eTk+1P-ek+1-eTkP-ek=eTkA-1-GA-2TP-A-1-GA-2ek+GTkGC2-C1TP-GC2-C1Gk+wTkGD2-D1TP-GD2-D1wk+2eTkA-1-GA-2TP-GC2-C1Gk+2eTkA-1-GA-2TP-GD2-D1wk+2GTkGC2-C1TP-GD2-D1wk-eTkP-ek.Since y(k) is measurable, z^2(k) can be substituted by y(k), and the observer state x^(k) can be written as(23)x^k=T-1z^1kyk.
One gets(24)x^k-xk=T-1z^1kyk-T-1z1kyk=T-1ek0.
SinceG(k)=g(k,T-1z^(k))-g(k,T-1z(k)) and g(k,x(k)) satisfies the Lipschitz condition (25)Gk2≤Lgx^k-xk2for a positive scalar ε1 we have (26)ε1GTkGk≤ε1Lg2x^k-xkTx^k-xk=ε1Lg2eTk0T-1TT-1ek0.
The above inequality is multiplied byλmin(TTT) on both sides, and we have (27)ε1λminTTTGTkGk≤ε1Lg2λminTTTeTk0T-1TT-1ek0≤ε1Lg2eTk0TT-1TTTT-1ek0=ε1Lg2eTk0ek0=ε1Lg2eTkekwhich is (28)ε1Lg2eTkek-ε1λminTTTGTkGk≥0,where λmin(·) is the smallest eigenvalue matrix.
We define(29)J1=∑k=0K1γ1eTkek-γ1wTkwk.
Under zero initial conditions, one gets(30)J1≤∑k=0K1γ1eTkek-γ1wTkwk+ΔVk.
Substituting (22) and (28) into (30), it follows that(31)J1≤∑k=0K1γ1eTkek-γ1wTkwk+ε1Lg2eTkek-ε1λminTTTGTkGk+ΔVk=∑k=0KeTkGTkwTkΩekGkwk,where (32)Ω=c11A-1-GA-2TP-GC2-C1A-1-GA-2TP-GD2-D1∗c22GC2-C1TP-GD2-D1∗∗c33,c11=A-1-GA-2TP-A-1-GA-2-P-+ε1Lg2+1γ1In-p,c22=GC2-C1TP-GC2-C1-ε1λminTTTIn,c33=GD2-D1TP-GD2-D1-γ1Id.
By using the Schur complement lemma,Ω≤0 is equivalent to (33)c11-1γ1In-pA-1-GA-2TP-GC2-C1A-1-GA-2TP-GD2-D1In-p∗c22GC2-C1TP-GD2-D10∗∗c330∗∗∗-γ1In-p≤0.
Furthermore, (33) can be rewritten as(34)A-1-GA-2TP-GC2-C1TP-GD2-D1TP-0P--1P-A-1-GA-2P-GC2-C1P-GD2-D10+-P-+ε1Lg200In-p∗-ε1λminTTTIn00∗∗-γ1Id0∗∗∗-γ1In-p≤0.
Let us use Schur complement lemma again; finally, we obtain (20), which guarantees e(k)2≤γ1w(k)2. So if (20) holds, then (19) is exponentially stable with a H∞ performance index e(k)2≤γ1w(k)2.Remark 6.
We use the measurable outputy(k)=z2(k) and the inequality λmin(P)aTa≤aTPa for any matrix P>0 and real vector a, to make the Lipschitz condition transform into ε1λmin(TTT)GT(k)G(k)≤ε1Lg2eT(k)e(k), which guarantees the diagonal element (3,3) is nonzero and negative and helps to form constraint (20).Remark 7.
In Theorem5, the LMI constraint (20) estimates the states by suppressing the influence of the disturbance term w(k) on e(k). By solving the convex optimization problem, the level γ1 can be obtained, in order to minimize e(k). Apparently, the NROO order n-p is less than the FOO [35] order n, and the transfer function from w(k) to e(k) is simplified. Therefore, NROO has the advantages of simple construction and better performance, which is compared to FOO.


### 3.2. NROO Existence Condition

On one hand, for the existence of NROO, the necessity of conditionrank(CD)=rank(D) has been proved in [36]. It is called observability with unknown input, which means the necessity of condition can guarantee the observability of systems with unknown disturbances w(k).On the other hand, for the existence of FOO, the necessity of condition requires the pair(A,C) to be observable [35]; that means (35)rankλIn-p-A11-A12-A21λIp-A220p×n-pIp=rankλIn-p-A11A21+p=n.It is not difficult to prove that the formula above is equivalent to the necessity of condition in NROO, because the nonsingular coordinate transformation which will be introduced in next section cannot transform the observability of the system, which means the NROO has the same range of application as FOO in this paper.


### 3.3. State Estimation

The NROO gain matrixG has been derived in the previous section. However, hypothesis output ρ(k) includes future output y(k+1) which is not available in practical implementation. In this section, a novel expression will be introduced for state estimation.Substituting (16) and (23) into (15), one gets (36)z^1k+1=A11-GA21z^1k+C⊥-GCgk,T-1z^1kyk+In-pμk+Gρk=A-1-GA-2z^1k+GC2-C1gk,T-1z^1kyk+In-pμk+Gρk.Then, substitutingμ(k) and ρ(k) into (36) yields(37)z^1k+1=A-1-GA-2z^1k+GC2-C1gk,T-1z^1kyk+In-pA12-GA22yk+Gyk+1+In-pB1-GB2uk.We denoteχ(k)=z^1(k)-Gy(k) in order to eliminate y(k+1) as follows:(38)χk+1=A-1-GA-2χk+Gyk+GC2-C1gk,T-1z^1kyk+In-pA12-GA22yk+In-pB1-GB2uk=A-1-GA-2χk+GC2-C1gk,T-1z^1kyk+A-1-GA-2+In-pA12-GA22yk+In-pB1-GB2uk.Finally, the state estimator is (39)x^k=T-1χk+Gykyk.


## 3.1. NROO Algorithm Design

Consider thatC is full-row rank, and there always exists a matrix C⊥∈R(n-p)×n such that C⊥C∈Rn×n is nonsingular (C⊥ can be given as an orthogonal basis of the null-space of C). Then under the nonsingular coordinate transformation z(k)=Tx(k), in which T=C⊥C is used, the dynamics (5) can be expressed into the form(10)z1k+1z2k+1=A11A12A21A22z1kz2k+C⊥Cgk,T-1zk+B1B2uk+D1D2wk,yk=0p×n-pIpz1kz2k,where z1(k)∈Rn-p and z2(k)∈Rp are new state vectors, Ip is identity matrix whose dimension is p×p, and (11)TAT-1=A11A12A21A22,TB=B1B2,TD=D1D2,CT-1=0p×n-pIp,zk=Txk=z1kz2k,and then (12)z2k=yk,z1k+1=A11z1k+A12yk+C⊥gk,T-1zk+B1uk+D1wk,yk+1=A21z1k+A22yk+Cgk,T-1zk+B2uk+D2wk.By using the hypothesis inputμ(k) and output ρ(k): (13)μk=A12yk+B1uk,ρk=yk+1-A22yk-B2ukwe get (14)z1k+1=A11z1k+C⊥gk,T-1zk+μk+D1wk,ρk=A21z1k+Cgk,T-1zk+D2wk.From the dynamics (14), NROO is constructed as (15)z^1k+1=A11z^1k+C⊥gk,T-1z^k+μk-Gρ^k-ρk,(16)ρ^k=A21z^1k+Cgk,T-1z^k,where z^1(k)∈Rn-p is the observer state of the reduced-order system (14), ρ^(k)∈Rp is the output of reduced-order observer, and G∈R(n-p)×p is the observer gain matrix.Denote(17)ek=z^1k-z1k,Gk=gk,T-1z^k-gk,T-1zkand then, from (14)–(16), the state estimation error dynamics are described by the following:(18)ek+1=A11-GA21ek+C⊥-GCGk+GD2-D1wk.Equation (18) is rewritten as the following form:(19)ek+1=A-1-GA-2ek+GC2-C1Gk+GD2-D1wk, where A-1=A11, A-2=A21, C1=-C⊥, and C2=-C.Theorem5 gives a constrained NROO design, whose state estimation performance is specified via a H∞ performance index.Theorem 5.
Let a prescribedH∞ performance level γ1>0, and if there exist a symmetric positive definite matrix P-∈R(n-p)×(n-p), a matrix Y-∈R(n-p)×p, and a ε1>0 such that the following condition holds:(20)-P-P-A-1-Y-A-2Y-C2-P-C1Y-D2-P-D10∗-P-+ε1Lg2In-p00In-p∗∗-ε1λminTTTIn00∗∗∗-γ1Id0∗∗∗∗-γ1In-p≤0,where Y-=P-G and ∗ denotes the symmetric elements in a matrix, then the error dynamics (18) satisfy H∞ performance index e(k)2≤γ1w(k)2, and the NROO gain matrix can be obtained by G=P--1Y-.Proof.
Choose the following Lyapunov function:(21)Vk=eTkP-ek.
According to error dynamics (18) and (21), the difference ΔV(k) is (22)ΔVk=Vk+1-Vk=eTk+1P-ek+1-eTkP-ek=eTkA-1-GA-2TP-A-1-GA-2ek+GTkGC2-C1TP-GC2-C1Gk+wTkGD2-D1TP-GD2-D1wk+2eTkA-1-GA-2TP-GC2-C1Gk+2eTkA-1-GA-2TP-GD2-D1wk+2GTkGC2-C1TP-GD2-D1wk-eTkP-ek.Since y(k) is measurable, z^2(k) can be substituted by y(k), and the observer state x^(k) can be written as(23)x^k=T-1z^1kyk.
One gets(24)x^k-xk=T-1z^1kyk-T-1z1kyk=T-1ek0.
SinceG(k)=g(k,T-1z^(k))-g(k,T-1z(k)) and g(k,x(k)) satisfies the Lipschitz condition (25)Gk2≤Lgx^k-xk2for a positive scalar ε1 we have (26)ε1GTkGk≤ε1Lg2x^k-xkTx^k-xk=ε1Lg2eTk0T-1TT-1ek0.
The above inequality is multiplied byλmin(TTT) on both sides, and we have (27)ε1λminTTTGTkGk≤ε1Lg2λminTTTeTk0T-1TT-1ek0≤ε1Lg2eTk0TT-1TTTT-1ek0=ε1Lg2eTk0ek0=ε1Lg2eTkekwhich is (28)ε1Lg2eTkek-ε1λminTTTGTkGk≥0,where λmin(·) is the smallest eigenvalue matrix.
We define(29)J1=∑k=0K1γ1eTkek-γ1wTkwk.
Under zero initial conditions, one gets(30)J1≤∑k=0K1γ1eTkek-γ1wTkwk+ΔVk.
Substituting (22) and (28) into (30), it follows that(31)J1≤∑k=0K1γ1eTkek-γ1wTkwk+ε1Lg2eTkek-ε1λminTTTGTkGk+ΔVk=∑k=0KeTkGTkwTkΩekGkwk,where (32)Ω=c11A-1-GA-2TP-GC2-C1A-1-GA-2TP-GD2-D1∗c22GC2-C1TP-GD2-D1∗∗c33,c11=A-1-GA-2TP-A-1-GA-2-P-+ε1Lg2+1γ1In-p,c22=GC2-C1TP-GC2-C1-ε1λminTTTIn,c33=GD2-D1TP-GD2-D1-γ1Id.
By using the Schur complement lemma,Ω≤0 is equivalent to (33)c11-1γ1In-pA-1-GA-2TP-GC2-C1A-1-GA-2TP-GD2-D1In-p∗c22GC2-C1TP-GD2-D10∗∗c330∗∗∗-γ1In-p≤0.
Furthermore, (33) can be rewritten as(34)A-1-GA-2TP-GC2-C1TP-GD2-D1TP-0P--1P-A-1-GA-2P-GC2-C1P-GD2-D10+-P-+ε1Lg200In-p∗-ε1λminTTTIn00∗∗-γ1Id0∗∗∗-γ1In-p≤0.
Let us use Schur complement lemma again; finally, we obtain (20), which guarantees e(k)2≤γ1w(k)2. So if (20) holds, then (19) is exponentially stable with a H∞ performance index e(k)2≤γ1w(k)2.Remark 6.
We use the measurable outputy(k)=z2(k) and the inequality λmin(P)aTa≤aTPa for any matrix P>0 and real vector a, to make the Lipschitz condition transform into ε1λmin(TTT)GT(k)G(k)≤ε1Lg2eT(k)e(k), which guarantees the diagonal element (3,3) is nonzero and negative and helps to form constraint (20).Remark 7.
In Theorem5, the LMI constraint (20) estimates the states by suppressing the influence of the disturbance term w(k) on e(k). By solving the convex optimization problem, the level γ1 can be obtained, in order to minimize e(k). Apparently, the NROO order n-p is less than the FOO [35] order n, and the transfer function from w(k) to e(k) is simplified. Therefore, NROO has the advantages of simple construction and better performance, which is compared to FOO.


## 3.2. NROO Existence Condition

On one hand, for the existence of NROO, the necessity of conditionrank(CD)=rank(D) has been proved in [36]. It is called observability with unknown input, which means the necessity of condition can guarantee the observability of systems with unknown disturbances w(k).On the other hand, for the existence of FOO, the necessity of condition requires the pair(A,C) to be observable [35]; that means (35)rankλIn-p-A11-A12-A21λIp-A220p×n-pIp=rankλIn-p-A11A21+p=n.It is not difficult to prove that the formula above is equivalent to the necessity of condition in NROO, because the nonsingular coordinate transformation which will be introduced in next section cannot transform the observability of the system, which means the NROO has the same range of application as FOO in this paper.


## 3.3. State Estimation

The NROO gain matrixG has been derived in the previous section. However, hypothesis output ρ(k) includes future output y(k+1) which is not available in practical implementation. In this section, a novel expression will be introduced for state estimation.Substituting (16) and (23) into (15), one gets (36)z^1k+1=A11-GA21z^1k+C⊥-GCgk,T-1z^1kyk+In-pμk+Gρk=A-1-GA-2z^1k+GC2-C1gk,T-1z^1kyk+In-pμk+Gρk.Then, substitutingμ(k) and ρ(k) into (36) yields(37)z^1k+1=A-1-GA-2z^1k+GC2-C1gk,T-1z^1kyk+In-pA12-GA22yk+Gyk+1+In-pB1-GB2uk.We denoteχ(k)=z^1(k)-Gy(k) in order to eliminate y(k+1) as follows:(38)χk+1=A-1-GA-2χk+Gyk+GC2-C1gk,T-1z^1kyk+In-pA12-GA22yk+In-pB1-GB2uk=A-1-GA-2χk+GC2-C1gk,T-1z^1kyk+A-1-GA-2+In-pA12-GA22yk+In-pB1-GB2uk.Finally, the state estimator is (39)x^k=T-1χk+Gykyk.


## 4. Model Predictive Control with Nonlinear Model Online Linearization

### 4.1. Nonlinear Predictive Model Online Linearization

In this section, MPC-NMOL was proposed to deal with the nonlinear state equation (5). Due to the time-consuming and computational complex problem, nonlinear optimization is converted into a quadratic optimization, by considering an online linearized method. Using gain scheduling technique [12], the online linearized dynamics at current sampling instance are given as (40)xk+1=A+g-k,xxk+Buk+Dwk,yk=Cxk,where the nonlinear-term g-(k,x) is linearized from g(k,x(k)), and one gets (41)g-x,k=∂gk,xk∂xkxk=x^kwhere the current operating state x(k) is defined by the estimation state x^(k) at kth sampling instance.We denote the general expression as follows:(42)gk,xk=g11x1x1+g12x2x2+g13x1x2g21x1x1+g22x2x2+g23x1x200,where the nonlinear-term coefficient in (42) can be calculated by M1, Ts, and gc(t,x(t)) mentioned before.Substitute (42) into (41), and we can get (43)g-k,x=g-1g-200g-3g-40000000000,where (44)g-1=g11x^1+g11x^1sgn⁡x^1+g13x^2,g-2=g12x^2+g12x^2sgn⁡x^2+g13x^1,g-3=g21x^1+g21x^1sgn⁡x^1+g23x^2,g-4=g22x^2+g22x^2sgn⁡x^2+g23x^1.We assume that(45)wk=wk+1∣k=⋯=wk+Np∣k,where Np is the predictive horizon. Then, considering that w(k) is the difference between the estimation state x^(k) and x(k) (calculated by x(k-1), u(k-1), and g-(k-1,x)), this means(46)wk=x^k-A+g-k-1,xxk-1+Buk-1.Here the wave disturbance model can be considered as the state disturbance and we can use it to estimate the external disturbance in vertical plane of AUV motion.Denote(47)A-=A+g-k-1,xand then the iterative predictive states over Np at step k is (48)x~k+1∣k=A-x^k+Bu~k∣k+wk,x~k+2∣k=A-x^k+1+Bu~k+1∣k+wk,⋮x~k+Np∣k=A-x^k+Np-1+Bu~k+Np-1∣k+wk,where u~k+Np∣k is the predictive input at step k.When we implement MPC algorithm, an incremental predictive model is always required, and therefore the inputu~k+i∣k can be replaced by Δu~k+i∣k which means(49)u~k+i∣k=u~k+i-1+Δu~k+i∣kand here we assume that u~k+i∣k will change just at every step (i<Nu) and remain constant after step k+i (Nu≤i<Np-1).With the assumption that the predictive value ofw(k) at sample instance k is zero, the vector output prediction equation can be calculated and expressed in condensed form which predicts the future dynamic behavior of the AUV longitudinal motion over the horizon Np:(50)y=Ψx^k+Γuk-1+Θu,where (51)y=y~k+1∣k⋮y~k+Np∣k,u=Δu~k∣k⋮Δu~k+Nu-1∣k,Ψ=CA-CA-2⋮CA-Np,Γ=CBCA-B+B⋮∑i=0Np-1CA-iB,Θ=CB⋯0CA-B+B⋯0⋮⋱⋮∑i=0Np-1CA-iB⋯∑i=0Np-NuCA-iB.


### 4.2. Consideration of Constraints

In order to guarantee the correct operation of the AUV, the constraints of the input elevator deflectionδ∈RNu over control horizon Nu and system output y∈R2Np over predictive horizon Np are used to define M and γ for each time step k considering the physical limitations of the driving device in practical implementation, and one obtains (52)δmin≤δ≤δmax,Δδmin≤Δδ≤Δδmax,ymin≤y≤ymax.For simplification of the following discussion, and without loss of generality, we set the control horizonNu=2, and then from (52) we can get(53)M1M2M3Δδ≤N1N2N3,where (54)M1=-Λ2Λ2,N1=-δmin+Λ1δk-1δmax-Λ1δk-1,M2=-I2×2I2×2,N2=-ΔδminΔδmax,M3=-ΘΘ,N3=-ymin+Ψx^k+Γδk-1ymax-Ψx^k-Γδk-1,Δδ=ΔδkΔδk+1,Λ1=11,Λ2=1011.We describe (53) by (55)MΔδ≤γ, which is equivalent to the constraints in the next section, and the constraints of the elevator deflection angle δ and pitch angle θ, corresponding to the inputs and outputs of the MPC, are characterized by M and γ.


### 4.3. Optimization with Constraints

Conventional MPC performance index can be written as(56)J=y-yrefTQy-yref+uTRu,where yref is a future reference vector and R and Q are positive definite weighted matrices.To simplify the expression, (56) can be rewritten as follows: (57)J=12uTEu+uTF+f0,s.t.Mu≤γ,where(58)E=2ΘTQΘ+R,F=2ΘTQΨx^k+Γuk-1-yref,f0=Ψx^k+Γuk-1-yrefTΨx^k+Γuk-1-yref.To minimize the quadratic function subject to (55), a QP (Quadratic Programming) problem has come out. Let us consider the expression which contains the Lagrange multipliers, that is, a QP problem subject to equality constraints Mu=γ, by the formula below: (59)J=12uTEu+uTF+f0+λTMu-γ.The minimization ofJ is to take the first partial derivatives with respect to u and λ, and we make them equal to zero and obtain the formula below:(60)∂J∂u=Eu+MTλ+F=0,(61)∂J∂λ=Mu-γ=0.The minimization ofJ can be made by finding the optimal u and λ via (60) and (61) where (62)λ=-ME-1MT-1γ+ME-1F,(63)u=-E-1MTλ+F=η-E-1MTλ, where η=-E-1F is the global optimal solution.The inequality constraints may comprise active constraints and inactive constraints in (55). We use both Mi and γi to form the ith inequality constraint. If Miu=γi, an inequality constraint Miu≤γi can be considered as active, and if Miu<γi, it is inactive. Here we use the Kuhn-Tucker conditions [37] to define the active and inactive constraints in terms of λ. If the active set were known, the original problem could become equality constrains problem in (59).In the conventional active set method [38], which belongs to the primal methods, the solutions are based on u (called decision variables). If the MIMO system has too many constraints, the calculations are complex and it is not a straightforward work.A dual method can be used to identify the constraints which are inactive systematically. The inactive constraints can be eliminated in the solution, andλ are called dual variables here. For constrained minimization problem, this method is a very simple programming procedure. The dual problem is derived from original primal problem as follows. Substituting (63) in (59), the dual problem is written as (64)minλ≥012λTHλ+λTK+12γTE-1γ,where the matrices H and K are given by (65)H=ME-1MT,K=γ+ME-1F.Subject toλ≥0, we minimize the dual performance index(66)J=12λTHλ+λTK+12γTE-1γ,where the set of λ are denoted as λ∗. By using Hildreth’s QP procedure [39] the dual problem is solved and the method can be written as (67)λik+1=max⁡0,βik+1with(68)βik+1=-1hiiki+∑j=1i-1hijλjk+1+∑j=i+1Nkhijλjk,where hij is the ijth element in H, ki is the ith element in K, and Nk is the number of rows of K. In this method, there are λj(k+1) and λj(k) in one iterative cycle. And we set λj(0)=0, at k=1, and the iterative procedure will converge to λ∗ as a result. Substitute λ∗ into (63) and we have (69)u=η-E-1MTλ∗,where λ∗=[β1(k+1),β2(k+1),…,βNk(k+1)]T.According to the receding horizon control in MPC, the first elements (elevator deflectionΔδ(k)) in u are taken to construct Δu(t).Remark 8.
Because Hildreth’s QP is a search-based point-by-point algorithm, there is no matrix inversion calculation. However, if the number of the active constraints is more than the number ofu’s or the active constraints are linearly dependent, then λ will not converge to λ∗ and the iteration will terminate at the largest value of the iterative counter. But the algorithm will not end because there is no matrix inversion calculation. In this case, finally, the algorithm will end in a near-optimal solution with the violation of constraints. This is the reason why we use Hildreth’s QP here, for its good ability to automatically recover from a deterioration constrained process.


### 4.4. Steps of NROO-Based MPC-NMOL with Constraints

(a)
Set values ofNP, Nu, and specify Q, R. (b)
Get the estimation of current statex^(k), sample current depth z(t), and pitch angle θ(t). (c)
Calculate matrixA- by online-linearization at current operating point which is defined by x^(k), to get matrixes Ψ, Γ, and Θ. Update constraints matrixes M and γ by using δ(t-1), and then E, F, and f0 can be calculated. (d)
Check if the global optimal solutionη satisfies the constraints. If so, make λ∗equal to zero vector and go to (f). If not, go to (e). (e)
Calculate matricesH and K, and then the dual variable λ∗ can be calculated from (61). (f)
GetΔu(t) from optimal solution u. (g)
Go to step (b).


## 4.1. Nonlinear Predictive Model Online Linearization

In this section, MPC-NMOL was proposed to deal with the nonlinear state equation (5). Due to the time-consuming and computational complex problem, nonlinear optimization is converted into a quadratic optimization, by considering an online linearized method. Using gain scheduling technique [12], the online linearized dynamics at current sampling instance are given as (40)xk+1=A+g-k,xxk+Buk+Dwk,yk=Cxk,where the nonlinear-term g-(k,x) is linearized from g(k,x(k)), and one gets (41)g-x,k=∂gk,xk∂xkxk=x^kwhere the current operating state x(k) is defined by the estimation state x^(k) at kth sampling instance.We denote the general expression as follows:(42)gk,xk=g11x1x1+g12x2x2+g13x1x2g21x1x1+g22x2x2+g23x1x200,where the nonlinear-term coefficient in (42) can be calculated by M1, Ts, and gc(t,x(t)) mentioned before.Substitute (42) into (41), and we can get (43)g-k,x=g-1g-200g-3g-40000000000,where (44)g-1=g11x^1+g11x^1sgn⁡x^1+g13x^2,g-2=g12x^2+g12x^2sgn⁡x^2+g13x^1,g-3=g21x^1+g21x^1sgn⁡x^1+g23x^2,g-4=g22x^2+g22x^2sgn⁡x^2+g23x^1.We assume that(45)wk=wk+1∣k=⋯=wk+Np∣k,where Np is the predictive horizon. Then, considering that w(k) is the difference between the estimation state x^(k) and x(k) (calculated by x(k-1), u(k-1), and g-(k-1,x)), this means(46)wk=x^k-A+g-k-1,xxk-1+Buk-1.Here the wave disturbance model can be considered as the state disturbance and we can use it to estimate the external disturbance in vertical plane of AUV motion.Denote(47)A-=A+g-k-1,xand then the iterative predictive states over Np at step k is (48)x~k+1∣k=A-x^k+Bu~k∣k+wk,x~k+2∣k=A-x^k+1+Bu~k+1∣k+wk,⋮x~k+Np∣k=A-x^k+Np-1+Bu~k+Np-1∣k+wk,where u~k+Np∣k is the predictive input at step k.When we implement MPC algorithm, an incremental predictive model is always required, and therefore the inputu~k+i∣k can be replaced by Δu~k+i∣k which means(49)u~k+i∣k=u~k+i-1+Δu~k+i∣kand here we assume that u~k+i∣k will change just at every step (i<Nu) and remain constant after step k+i (Nu≤i<Np-1).With the assumption that the predictive value ofw(k) at sample instance k is zero, the vector output prediction equation can be calculated and expressed in condensed form which predicts the future dynamic behavior of the AUV longitudinal motion over the horizon Np:(50)y=Ψx^k+Γuk-1+Θu,where (51)y=y~k+1∣k⋮y~k+Np∣k,u=Δu~k∣k⋮Δu~k+Nu-1∣k,Ψ=CA-CA-2⋮CA-Np,Γ=CBCA-B+B⋮∑i=0Np-1CA-iB,Θ=CB⋯0CA-B+B⋯0⋮⋱⋮∑i=0Np-1CA-iB⋯∑i=0Np-NuCA-iB.


## 4.2. Consideration of Constraints

In order to guarantee the correct operation of the AUV, the constraints of the input elevator deflectionδ∈RNu over control horizon Nu and system output y∈R2Np over predictive horizon Np are used to define M and γ for each time step k considering the physical limitations of the driving device in practical implementation, and one obtains (52)δmin≤δ≤δmax,Δδmin≤Δδ≤Δδmax,ymin≤y≤ymax.For simplification of the following discussion, and without loss of generality, we set the control horizonNu=2, and then from (52) we can get(53)M1M2M3Δδ≤N1N2N3,where (54)M1=-Λ2Λ2,N1=-δmin+Λ1δk-1δmax-Λ1δk-1,M2=-I2×2I2×2,N2=-ΔδminΔδmax,M3=-ΘΘ,N3=-ymin+Ψx^k+Γδk-1ymax-Ψx^k-Γδk-1,Δδ=ΔδkΔδk+1,Λ1=11,Λ2=1011.We describe (53) by (55)MΔδ≤γ, which is equivalent to the constraints in the next section, and the constraints of the elevator deflection angle δ and pitch angle θ, corresponding to the inputs and outputs of the MPC, are characterized by M and γ.


## 4.3. Optimization with Constraints

Conventional MPC performance index can be written as(56)J=y-yrefTQy-yref+uTRu,where yref is a future reference vector and R and Q are positive definite weighted matrices.To simplify the expression, (56) can be rewritten as follows: (57)J=12uTEu+uTF+f0,s.t.Mu≤γ,where(58)E=2ΘTQΘ+R,F=2ΘTQΨx^k+Γuk-1-yref,f0=Ψx^k+Γuk-1-yrefTΨx^k+Γuk-1-yref.To minimize the quadratic function subject to (55), a QP (Quadratic Programming) problem has come out. Let us consider the expression which contains the Lagrange multipliers, that is, a QP problem subject to equality constraints Mu=γ, by the formula below: (59)J=12uTEu+uTF+f0+λTMu-γ.The minimization ofJ is to take the first partial derivatives with respect to u and λ, and we make them equal to zero and obtain the formula below:(60)∂J∂u=Eu+MTλ+F=0,(61)∂J∂λ=Mu-γ=0.The minimization ofJ can be made by finding the optimal u and λ via (60) and (61) where (62)λ=-ME-1MT-1γ+ME-1F,(63)u=-E-1MTλ+F=η-E-1MTλ, where η=-E-1F is the global optimal solution.The inequality constraints may comprise active constraints and inactive constraints in (55). We use both Mi and γi to form the ith inequality constraint. If Miu=γi, an inequality constraint Miu≤γi can be considered as active, and if Miu<γi, it is inactive. Here we use the Kuhn-Tucker conditions [37] to define the active and inactive constraints in terms of λ. If the active set were known, the original problem could become equality constrains problem in (59).In the conventional active set method [38], which belongs to the primal methods, the solutions are based on u (called decision variables). If the MIMO system has too many constraints, the calculations are complex and it is not a straightforward work.A dual method can be used to identify the constraints which are inactive systematically. The inactive constraints can be eliminated in the solution, andλ are called dual variables here. For constrained minimization problem, this method is a very simple programming procedure. The dual problem is derived from original primal problem as follows. Substituting (63) in (59), the dual problem is written as (64)minλ≥012λTHλ+λTK+12γTE-1γ,where the matrices H and K are given by (65)H=ME-1MT,K=γ+ME-1F.Subject toλ≥0, we minimize the dual performance index(66)J=12λTHλ+λTK+12γTE-1γ,where the set of λ are denoted as λ∗. By using Hildreth’s QP procedure [39] the dual problem is solved and the method can be written as (67)λik+1=max⁡0,βik+1with(68)βik+1=-1hiiki+∑j=1i-1hijλjk+1+∑j=i+1Nkhijλjk,where hij is the ijth element in H, ki is the ith element in K, and Nk is the number of rows of K. In this method, there are λj(k+1) and λj(k) in one iterative cycle. And we set λj(0)=0, at k=1, and the iterative procedure will converge to λ∗ as a result. Substitute λ∗ into (63) and we have (69)u=η-E-1MTλ∗,where λ∗=[β1(k+1),β2(k+1),…,βNk(k+1)]T.According to the receding horizon control in MPC, the first elements (elevator deflectionΔδ(k)) in u are taken to construct Δu(t).Remark 8.
Because Hildreth’s QP is a search-based point-by-point algorithm, there is no matrix inversion calculation. However, if the number of the active constraints is more than the number ofu’s or the active constraints are linearly dependent, then λ will not converge to λ∗ and the iteration will terminate at the largest value of the iterative counter. But the algorithm will not end because there is no matrix inversion calculation. In this case, finally, the algorithm will end in a near-optimal solution with the violation of constraints. This is the reason why we use Hildreth’s QP here, for its good ability to automatically recover from a deterioration constrained process.


## 4.4. Steps of NROO-Based MPC-NMOL with Constraints

(a)
Set values ofNP, Nu, and specify Q, R. (b)
Get the estimation of current statex^(k), sample current depth z(t), and pitch angle θ(t). (c)
Calculate matrixA- by online-linearization at current operating point which is defined by x^(k), to get matrixes Ψ, Γ, and Θ. Update constraints matrixes M and γ by using δ(t-1), and then E, F, and f0 can be calculated. (d)
Check if the global optimal solutionη satisfies the constraints. If so, make λ∗equal to zero vector and go to (f). If not, go to (e). (e)
Calculate matricesH and K, and then the dual variable λ∗ can be calculated from (61). (f)
GetΔu(t) from optimal solution u. (g)
Go to step (b).


## 5. Simulation Results

In this paper, simulations are presented to demonstrate the effectiveness of NROO-based MPC-NMOL. The method is used in a given depth control of REMUS AUV which is developed by MIT (Massachusetts Institute of Technology). The values of nonlinear model parameters are shown as follows:(70)A=0.99210.006300.00020.04530.99610-0.00710.0101-0.015100.0101,B=-0.0065-0.0400,C=00100001,D=0.0002000.00120.0002000.0012,gk,xk=M-1-0.02x1x1-0.0066x2x20.0046x1x1-0.2272x2x200.Here the physical parameters of REMUS AUV, which can be found from [40], are shown in Table 1.Table 1
Physical Parameters of REMUS AUV.
Description Values Units m 30.48 kg Length 1.33 m Beam 0.26 m Ixx 0.177 Kg⋅m2 Iyy 3.45 Kg⋅m2 Izz 3.45 Kg⋅m2 Location of CG (0,0,0.0196) m Location of BG (−0.611,0,0) m


### 5.1. Nonlinear Reduced-Order Observer Design

From the result thatrank(CD)=rank(D)=2 and (A,C) are observable, we can easily verify the existence of NROO.By using Schmidt orthogonalization, we can get the matrixT from the combination of C1 and its standard orthogonal basis C2. Obviously, the nonsingular matrix T has only one form which is (71)T=1000010000100001.We chooseLg=0.001 as the value of the Lipschitz constant, and then the matrices parameter values in (19) can be calculated as follows:(72)A-1=0.99210.00630.04530.9961,A-2=0.01000.01,C1=-10000-100,C2=00-10000-1,D1=0.0002000.0012,D2=0.0002000.0012.With the help of MATLAB LMI toolbox, condition (20) is solved to obtain γ1=2.3561×103; at the same time, other results are found as follows:(73)P-=144.6253-54.7977-54.797733.5797,Y-=42.6923-1.27540.495741.8309,ε1=4.3871×103.One obtains (74)G-=P--1Y-=0.78801.21351.30073.2259.Figures2 and 3 show the estimation of x1(k) and x2(k) by using NROO and FOO. Although both observers can make the state estimation error converge asymptotically, compared to FOO, NROO has a better performance of state estimation.Figure 2
Comparison ofx1 and x^1 (a) using the NROO and (b) using the FOO.
(a) (b)Figure 3
Comparison ofx2 and x^2 (a) using the NROO and (b) using the FOO.
(a) (b)


### 5.2. Implementation of MPC-NMOL

The parameters of MPC are chosen asNp=30, Nu=2, Q=INp×Np, and R=10. The initial values of state variables are all zero. The waves disturbance exposed on REMUS is assumed at a level 3 sea state, Hs=0.88(m), β=45∘, Cd=0.65, Cm=1.95, and N=271. The surge speed is Uc=1.51(m/s), the desired depth is 5(m), and pitch angle is 0∘. Input and output constraints are(75)δmax=-δmin=30∘,Δδmax=-Δδmin=5∘/s,-90∘≤θ≤90∘.Case 1.
First, we assume that there is no disturbance in simulation process, and NMPC is used to compare with the proposed method. Both of the two methods (MPC-NMOL and NMPC) have all state variables measurable. Figure4 compares MPC-NMOL and NMPC simulation results of depth output. Figure 5 compares the simulation results of pitch output. Figure 6 compares the elevator deflection angle input.Figure 4
Comparison of depth with no disturbance.Figure 5
Comparison of pitch angle with no disturbance.Figure 6
Comparison of elevator deflection angle with no disturbance.Case 2.
Next, it is assumed that wave disturbance affects the state process, and the other condition is the same as Case1. Figures 7–9 show the comparison of depths, pitch angles, and elevator deflection angles with wave disturbance.Figure 7
Comparison of depth with disturbance.Figure 8
Comparison of pitch angle with disturbance.Figure 9
Comparison of elevator deflection angle with disturbance.Figure10 shows the wave force and moment, which can be seen as state process disturbance and assumed to be zero-mean white noise sequence. The wave force and moment are calculated from (7), so they have the same form but different amplitude. Furthermore, when we simulate the wave force and moment, we choose multiple influential frequencies which are near the given main frequency ofP-M spectrum to superimpose the irregular waves.Figure 10
Wave force and moment.All these results in Cases1 and 2 demonstrate that AUV could achieve the desired depth and pitch angle under the wave disturbance. In addition, the input signals in MPC-NMOL are smooth and without control signal saturation.


## 5.1. Nonlinear Reduced-Order Observer Design

From the result thatrank(CD)=rank(D)=2 and (A,C) are observable, we can easily verify the existence of NROO.By using Schmidt orthogonalization, we can get the matrixT from the combination of C1 and its standard orthogonal basis C2. Obviously, the nonsingular matrix T has only one form which is (71)T=1000010000100001.We chooseLg=0.001 as the value of the Lipschitz constant, and then the matrices parameter values in (19) can be calculated as follows:(72)A-1=0.99210.00630.04530.9961,A-2=0.01000.01,C1=-10000-100,C2=00-10000-1,D1=0.0002000.0012,D2=0.0002000.0012.With the help of MATLAB LMI toolbox, condition (20) is solved to obtain γ1=2.3561×103; at the same time, other results are found as follows:(73)P-=144.6253-54.7977-54.797733.5797,Y-=42.6923-1.27540.495741.8309,ε1=4.3871×103.One obtains (74)G-=P--1Y-=0.78801.21351.30073.2259.Figures2 and 3 show the estimation of x1(k) and x2(k) by using NROO and FOO. Although both observers can make the state estimation error converge asymptotically, compared to FOO, NROO has a better performance of state estimation.Figure 2
Comparison ofx1 and x^1 (a) using the NROO and (b) using the FOO.
(a) (b)Figure 3
Comparison ofx2 and x^2 (a) using the NROO and (b) using the FOO.
(a) (b)


## 5.2. Implementation of MPC-NMOL

The parameters of MPC are chosen asNp=30, Nu=2, Q=INp×Np, and R=10. The initial values of state variables are all zero. The waves disturbance exposed on REMUS is assumed at a level 3 sea state, Hs=0.88(m), β=45∘, Cd=0.65, Cm=1.95, and N=271. The surge speed is Uc=1.51(m/s), the desired depth is 5(m), and pitch angle is 0∘. Input and output constraints are(75)δmax=-δmin=30∘,Δδmax=-Δδmin=5∘/s,-90∘≤θ≤90∘.Case 1.
First, we assume that there is no disturbance in simulation process, and NMPC is used to compare with the proposed method. Both of the two methods (MPC-NMOL and NMPC) have all state variables measurable. Figure4 compares MPC-NMOL and NMPC simulation results of depth output. Figure 5 compares the simulation results of pitch output. Figure 6 compares the elevator deflection angle input.Figure 4
Comparison of depth with no disturbance.Figure 5
Comparison of pitch angle with no disturbance.Figure 6
Comparison of elevator deflection angle with no disturbance.Case 2.
Next, it is assumed that wave disturbance affects the state process, and the other condition is the same as Case1. Figures 7–9 show the comparison of depths, pitch angles, and elevator deflection angles with wave disturbance.Figure 7
Comparison of depth with disturbance.Figure 8
Comparison of pitch angle with disturbance.Figure 9
Comparison of elevator deflection angle with disturbance.Figure10 shows the wave force and moment, which can be seen as state process disturbance and assumed to be zero-mean white noise sequence. The wave force and moment are calculated from (7), so they have the same form but different amplitude. Furthermore, when we simulate the wave force and moment, we choose multiple influential frequencies which are near the given main frequency ofP-M spectrum to superimpose the irregular waves.Figure 10
Wave force and moment.All these results in Cases1 and 2 demonstrate that AUV could achieve the desired depth and pitch angle under the wave disturbance. In addition, the input signals in MPC-NMOL are smooth and without control signal saturation.


## 6. Conclusions

In this paper, a NROO-based model predictive controller with nonlinear model online linearization for AUV in vertical plane is presented, which controls the depth and pitch angle. This design uses the NROO to estimate the states used in MPC. The design process of the controller also takes into account the practical elevator deflection constraints and output constraints. By using a Hildreth’s QP procedure, the constraints can be simply handled. Making use of the proposed MPC methods, the AUV can navigate in vertical plane with desired depth and pitch angle. It is robust against rough wave disturbance near surface. The simulations carried out provide the validation of the proposed methods, presenting fast dynamical response and strong robustness to external disturbances. Accurate control and state estimation can also be achieved.



---

*Source: 4394571-2017-01-03.xml*