# Efficient Multivariable Generalized Predictive Control for Autonomous Underwater Vehicle in Vertical Plane

Authors: Xuliang Yao; Guangyi Yang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

This paper presents the design and simulation validation of a multivariable GPC (generalized predictive control) for AUV (autonomous underwater vehicle) in vertical plane. This control approach has been designed in the case of AUV navigating with low speed near water surface, in order to restrain wave disturbance effectively and improve pitch and heave motion stability. The proposed controller guarantees compliance with rudder manipulation, AUV output constraints, and driving energy consumption. Performance index based on pitch stabilizing performance, energy consumption, and system constraints is used to derive the control action applied for each time step. In order to deal with constrained optimization problems, a Hildreth’s QP procedure is adopted. Simulation results of AUV longitudinal control show better stabilizing performance and minimized energy consumption improved by multivariable GPC.


---

## Body

## 1. Introduction

For AUV (autonomous underwater vehicle), its motion attitude is severely affected by environmental disturbances if it is navigating near water surface. And wave disturbance is the most obvious one. Violent pitch and heave motion often discontinue normal work of AUV [1]. So, it is necessary to design an effective control pattern to the issue of AUV longitudinal motion control. For the AUV that pitch and heave motion are underactuated, external auxiliary device is needed to adjust its attitude. The bow and stern rudders are space-saving and can provide it with antipitch force and torque [2]. Meanwhile, energy consumption for longitudinal motion control is also very large. For small-scale AUV, the carried energy is limited. So, it is necessary to reduce energy consumption for driving rudders, while stabilizing performance is satisfactory.A variety of methods have been proposed to reduce AUV energy consumption. A dynamic resistance mathematical model of AUV is built in [3], which states the relationship between resistance and velocity, and they applied the terminal SMC to decrease the pitch angle of AUV and optimize the energy consumption. Energy consumption caused by overcoming the resistance and the yaw was analyzed in [4], and the self-tuning PID controller based on the MOGA method is used to optimize the performance index which has the added resistance consumption. Relevant research work based on roll additional resistance is proposed in [5]. Energy consumption used for driving fin stabilizer is added to performance in [6], GPC controller is designed, and recursive least square method is adopted for parameter estimation. An optimal path planning control of AUV based on rolling optimization is proposed in [7], in order to conduct the path of optimal energy consumption.In the vertical plane, pitch motion is dynamically coupled with surge and heave [8]. So, control system based on AUV model with single degree of freedom usually cannot attain expected pitch and heave stabilizing performance in practical application, and control systems for AUV are often designed with an implicit assumption that vehicle’s surge velocity is constant. Summing up the above, traditional SISO system controller cannot adapt to the coupled AUV dynamic model.Regarding multivariable generalized predictive control (GPC), it is one of the most popular formulations of MPC methods [9], which has been successfully implemented in many MIMO industry process control. 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 main research areas in marine control systems field. Several research works in this area were based on using GPC as the main controller for yaw [10], roll [6], depth [11], and propulsion BLDC motor [12].The implementation of GPC with constraints requires the solution of a QP (quadratic programming) problem, that is, an optimization problem discribed by a quadratic objective function and linear constraints. There are many techniques, such as the ones based on feasible direction methods [13], decreasing ellipsoid volume methods [14], pivoting methods [15], and active set methods [16]. The active set methods belong to primal methods which are commonly used in QP procedure. But the active constraints need to be identified along with the decision variable; the computational work is quite large if there are a few constraints. A dual method can systematically identify the constraints that are not active, and an algorithm called Hildreth’s QP procedure [17] leads to a very simple programming procedure for constrained optimization problems.In this paper, the main contributions are as follows. First, coupling effect of surge and heave on pitch is considered and made as a premise of issue discussion; then an augmented state space model is transformed to a CARIMA model which is suitable for GPC. In order to be more similar to the actual situation, wave disturbance is considered in vertical plane of AUV motion. Second, the energy consumption of driving rudders is calculated and added to the performance index, which expends the application range of GPC used for energy-saving and attitude control of AUV. Third, by using the novel performance index, Hildreth’s QP is adopted to deal with the constraints and obtain the optimal solution, which can be considered as the compromise between energy-saving and attitude stabilization. The rest of the paper is organized as follows. Section2 described the AUV vertical plane dynamics. GPC design is proposed in Section 3. Section 4 gives the simulation results of an AUV to verify the effectiveness of proposed approach.


## 2. AUV Vertical Plane Dynamics

The vertical plane model of AUV is linearized for predefined operating conditions to extract the linear model. First, the matrix equations of motion can be given as follows:(1)m-Zw˙-mxG+Zq˙00-mxG+Mw˙Iyy-Mq˙0000100001w˙q˙z˙θ˙-ZwmU+Zq00Mw-mxGU+Mq0-mgh100-U0100wqzθ=ZδbMδb00ZδsMδs00δbδs+ZwaveMwave00,which can be simply expressed as discretized state space form(2)xmt+1=Amxmt+Bmut+Dmdt,yt=Cmxmt.An augmented state space model can be used if the model input isΔu(t)=u(t)-u(t-1), and Δxm(t)=xm(t)-xm(t-1). ε(t)=d(t)-d(t-1) is a zero-mean white noise sequence.Notice that, from (2), the following difference equation can be given as (3)xmt=Amxmt-1+Bmut-1+Dmdt-1, and then subtracting (3) from (2) leads to(4)Δxmt+1=AmΔxmt+BmΔut+Dmεt.In order to relate the outputy(t) to the state variable Δxm, we deduce that (5)Δyt+1=CmΔxmt+1=CmAmΔxmt+CmBmΔut+CmDmεt,where Δy(t+1)=y(t+1)-y(t).Choosing a new state variable vectorx(t)=Δxm(t)Ty(t)TT, we have the following:(6)Δxmt+1yt+1=AmomTCmAmIq×qΔxmtyt+BmCmBmΔut+DmCmDmεt,yt=omIq×qΔxmtyt,where Iq×q is the identity matrix with dimensions q×q, which is the number of outputs, and om is a q×n1 zero matrix. In (6), Am, Bm, Cm, and Dm have dimensions n1×n1, n1×m, q×n1, and n1×2, respectively.For notational simplicity, we denote (6) by(7)xt+1=Adxt+BdΔut+Pεt,yt=Cxt,where Ad, Bd, P, and C are matrices corresponding to the forms given in (7). In the following, the dimensionality of the augmented state space equation is taken to be n=n1+q.Remark 1.
The disturbance due to surface waves is the main disturbance type which is considered in this paper and written asd(t). Wave height is affected by wind, and the wave profile is affected by water depth. Sea state is used to classify the sea conditions and expressed by wind speed and wave height. In surface waves modeling, we assumed the long crested and unidirectional sea state, which means all the waves travel in a primary direction. In order to complete an irregular long crested wave sea state simulation, the superposition principle is used here. The single parameter Pierson-Moskowitz spectrum [18] is used in this model and its spectrum density is defined by the following equation:(8)Sω=8.1×10-3g3ω5exp⁡-3.11Hs2ω4,where Hs is the significant wave height (m), ω is wave frequency. Assuming the AUV body is small compared with the incoming wave and the wave speed is large enough in relation to the body diameter D0. The transient force and moment are found by the use of Morison’s equation [19] and integration along the length of the body L0 at each step as follows:(9)Zwave=∫L0CdρD02ww-w2+CmρπD024w˙w-w˙2dx,Mwave=∫L0CdρD02ww-w2+CmρπD024w˙w-w˙2xdx,where ρ is fluid density, Cd is the drag coefficient, and Cm is the added mass coefficient.
In (9), (10)ww=∑i=1Nωeiαiexp⁡-kizsin⁡ωeit+θi,w˙w=-∑i=1Nωei2αiexp⁡-kizcos⁡ωeit+θi, where N is the equal interval frequency waveband and chosen in this thesis to discretize the wave spectrum, ωei is encountering frequency, αi=2S(ωi)dω, ki is wave number, and θi is random phase shift (0<θi<2π) in relation to each frequency.


## 3. Multivariable GPC Controller Design

### 3.1. Multivariable GPC

It is known that the GPC algorithm is based on using a CARIMA model(11)Az-1yt=Bz-1ut-1+ξz-1ηtΔ,where the unmeasurable disturbance is given by a white noise colored by η(t).We consider the most usual case whenξ(z-1)=In×n, because the coloring polynomials are very difficult to estimate with sufficient accuracy in practice, especially in the multivariable case.Multiply (11) by Δ(12)A-z-1yt=Bz-1Δut-1+ηt.With zero initial conditions, it can be seen that both (7) and (11) are equivalent if (13)CzI-Ad-1Bd=A-z-1-1Bz-1z-1,CzI-Ad-1Pdt=A-z-1-1ηt.By makingd(t)=η(t), we can get the matrix A-(z-1) and B(z-1).The vector output prediction equation can be calculated and expressed in condensed form which predicts the future dynamic behavior of the AUV longitudinal motion along the horizonNp:(14)y=Gu+f,where u=Δδ=[Δδ(t+1)T⋯Δδ(t+Nu)T]T, and other procedures are detailed in [9].Notice that matrixA(z-1) is diagonal; the problem is reduced to the recursion of n scalar Diophantine equations, which become much simpler to program and require less computation. The computation of G and f is also considerably simplified.Conventional GPC performance index can be written as follows:(15)J=y-wTQy-w+uTRu,where w is a future reference sequence and R and Q are positive definite weighting matrices.Remark 2.
CARIMA (controlled auto regressive integrated moving average) model [20] has been long used in GPC. The main reason we use the I-O model (CARIMA model) rather than state space model (7) is that we can directly use the input and output signals of the plant to avoid the use of a state observer, and the complexity of the controller’s front stage design is reduced. However, system identification technology is usually used in GPC to obtain the required I-O model (GPC can be considered as an adaptive control when the system identification is carried out online). Apparently, it will increase the difficulty of controller design by using system identification. State space model of AUV is widely used in many classic literatures [21, 22] and can be modeled from AUV dynamics equations. Therefore, it is wise and simple to obtain an I-O model by equivalent transformation of state space model.


### 3.2. Improved Performance Index

The energy consumption used for driving bow and stern rudders can be expressed as(16)E0=λ0·Δδ·Mr,where λ0 is selected according to energy consumption used for pitch stabilizing and depth tracking and Δδ is variation of rudder angle. Mr is the driving moment of rudder.Then the Lagrange dynamics equation is chosen to analyze flat wing-type rudder with thin foil condition as shown in Figure1. Rudder length is 2a, the distance from midpoint to shaft is c, and cross-section minor axis length is b.Figure 1
Flat wing-type rudder with thin foil condition.Lagrange functionL is defined as the difference between flat wing kinetic energy T and potential energy P. By setting the flap wing rotational angle δb,s(t), without considering the effects of the potential energy (under the thin foil condition, the quality of flat wing can be ignored, and the flat wing flaps within a small range), L can be expressed as(17)L=T=12Jδ˙b,s2t.The second Lagrange equation of flat wing can be described as(18)ddt∂L∂δ˙b,st-∂L∂δb,st=Mb,st.From (17) we can get (19)ddt∂L∂δ˙b,st=Jδ¨b,st,∂L∂δb,st=0.Substitute (19) into (18); then the driving moment of rudder can be obtained as follows:(20)Mb,st=Jδ¨b,st=Jω˙t,where J=πρ[a2c2+(a2-b2)2/8] is the additional moment inertia and ω(t) is the rotation angular velocity, without loss of generality, and for simplification of the following discussion, we set the control horizon Nu=2, and then(21)Mr=Mb,st+1TMb,st+2TT.Equation (21) and Δδ are substituted into (16); it follows that(22)E0=λ0Δδb,st+1TΔδb,st+2TMb,st+1Mb,st+2=λ0Δδb,st+1TMb,st+1+Δδb,st+2TMb,st+2,where Δδb,s(t+1)=δb,s(t+1)-δb,s(t)=ω(t+1+d)·Ts, Ts is sampling period and d is the delay between the angular velocity and driving rudder signal. Considering internal complicated mechanical structure of rudders, d is set as 1, and (23) is obtained:(23)ωt+2=Δδb,st+1Ts.Denote(24)Mb,st+1=Jω˙t+1=Jωt+2-ωt+1Ts.Substitute (23) into (24), and (25) is obtained:(25)Mb,st+1=JΔδb,st+1-Δδb,stTs2.Substituting (25) into (22), we obtain (26)E0=λnΔδTΔδ-F0Δδ,where λn=λ0J/Ts2, F0=c0c1, c0=Δδ(t) is the last step variation of rudder angle which is known, and c1=Δδ(t+1) is unknown.For a brief expression, (15) can be rewritten with (26) as follows: (27)J=12uTEu+uTFn+f0,s.t.Mu≤γ, where(28)E=2GTQG+R+λnI,Fn=F-F0,f0=f-wTf-w,F=2GTQf-w,u=Δδ.The constraints of the rudders deflection and pitch angle, corresponding to the inputs and outputs of the GPC, are characterized byM and γ and described next.


### 3.3. Consideration of Constraints

In practical application, considering the physical limitations of the driving device, and in order to ensure the correct operation of the AUV, the constraints of the rudders deflection and system output used to defineM and γ for each time step t can be defined as follows: (29)δmin≤δ≤δmax,Δδmin≤Δδ≤Δδmax,ymin≤y≤ymax,and then from (29) we can get(30)M1M2M3Δδ≤N1N2N3, where (31)M1=-C2C2,M2=-II,M3=-GG,N1=-δmin+C1δt-1δmax-C1δt-1,N2=-ΔδminΔδmax,N3=-Ymin+fYmax-f,C1=I2×2I2×2,C2=I2×2I2×20I2×2.For compactness of expression, we describe (30) by(32)MΔδ≤γ,which is equivalent to the constraints in (27).


### 3.4. Optimization with Equality Constraints

To minimize the quadratic function subject to (32), the performance index here becomes a QP (quadratic programming) problem. Let us consider the expression which contains the Lagrange multipliers, that is, a QP problem subject to equality constraints Mu=γ; we construct (33)J=12uTEu+uTFn+f0+λTMu-γ.The minimization ofJ is to take the first partial derivatives with respect to u and λ and make the results equal to zero. We obtain (34)∂J∂u=Eu+MTλ+∂uTFn∂u=0,(35)∂J∂λ=Mu-γ=0.Denote(36)uT=u11u12u21u22,Fn=F-F0=F1F2-c0c1=FcF2-c1=f11f12f21-u11f22-u12,uTFn=u11f11+u12f12+u21f21+u22f22-u21u11-u22u12.Then, in (34),(37)∂uTFn∂u=f11-u21f12-u22f21-u11f22-u12=FcF2-0010000110000100u=Fu-Su.The minimization ofJ can be made by finding the optimal u and λ via (34) and (35) where (38)λ=-ME-S-1MT-1γ+ME-S-1Fu,(39)u=-E-S-1MTλ+Fu=η-E-S-1MTλ,where η=-E-S-1Fu is the global optimal solution.


### 3.5. Optimization with Inequality Constraints

In (32), the inequality constraints may comprise active and inactive constraints. If Miu=γi, an inequality constraint Miu≤γi is said to be active, and if Miu<γi, it is inactive, where Mi and γi together form the ith inequality constraint. The Kuhn-Tucker conditions [13] are widely used to define the active and inactive constraints in terms of λ. If the active set were known, the original problem could become an equality constrains problem in (33).The conventional active set method in [9] belongs to the primal methods, in which the solutions are based on u (called decision variables). If the MIMO system has too many constraints, the calculation is quite large and it is not a straightforward work.To identify the constraints which are inactive systematically, a dual method can be used. 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 (39) in (33), the dual problem is written as (40)minλ≥012λTHλ+λTK+P1,where the matrices H and K are given by (41)P1=12FuTE-S-12Fn-Fu+λTME-S-1Fn-Fu-f0,H=ME-S-1MT,K=γ+ME-S-1Fu.The set ofλ which minimize the dual performance index(42)J=12λTHλ+λTK+P1,subject to λ≥0, are denoted as λ∗. Here Hildreth’s QP procedure [23] was used to solve the dual problem. The method can be written as (43)λit+1=max⁡0,βit+1,with(44)βit+1=-1hiiki+∑j=1i-1hijλjt+1+∑j=i+1Nkhijλjt,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 two sets of λ in one iterative cycle, λj(t+1) and λj(t). So, we set λj(0)=0, at t=1, and the iterative procedure converged to λ∗ as a result; by (39) we have (45)u=η-E-S-1MTλ∗,where λ∗=[β1(t+1),β2(t+1),…,βn(t+1)]T.According to the receding horizon control in GPC, the first two elements (bow and stern rudders deflection) inu are taken to form the Δu(t).Remark 3.
Because Hildreth’s QP is an element-by-element search based algorithm, there is no need to calculate any matrix inversion. However, if the active constraints are linearly dependent or the number of the active constraints is more than the number ofu, then λ will not converge to λ∗ and the iteration will terminate at the largest value of the iterative counter. But the algorithm will continue because there is no matrix inversion calculation. In this case, the algorithm will get a near-optimal solution with constraints if the constraints are violated. This is a main reason for using Hildreth’s QP in this paper, because of the ability to automatically recover from a deterioration constrained process.


### 3.6. Steps of Multivariable GPC with Constraints

(a)
Set values ofNP, Nu, and specify Q, R. Calculate matrix G and vector f; then E and F can be calculated. (b)
Sample current depthz(t) and pitch angle θ(t) and update constraints matrix M and γ by using δ(t-1). (c)
Check if the global optimal solutionη satisfies the constraints. If so, make λ∗ equal to zero vector and go to (e). If not, go to (d). (d)
Calculate matricesH and K; then the dual variable λ∗ can be calculated from (35). (e)
ReturnΔu(t) from optimal solution u. (f)
Go to step (b).


## 3.1. Multivariable GPC

It is known that the GPC algorithm is based on using a CARIMA model(11)Az-1yt=Bz-1ut-1+ξz-1ηtΔ,where the unmeasurable disturbance is given by a white noise colored by η(t).We consider the most usual case whenξ(z-1)=In×n, because the coloring polynomials are very difficult to estimate with sufficient accuracy in practice, especially in the multivariable case.Multiply (11) by Δ(12)A-z-1yt=Bz-1Δut-1+ηt.With zero initial conditions, it can be seen that both (7) and (11) are equivalent if (13)CzI-Ad-1Bd=A-z-1-1Bz-1z-1,CzI-Ad-1Pdt=A-z-1-1ηt.By makingd(t)=η(t), we can get the matrix A-(z-1) and B(z-1).The vector output prediction equation can be calculated and expressed in condensed form which predicts the future dynamic behavior of the AUV longitudinal motion along the horizonNp:(14)y=Gu+f,where u=Δδ=[Δδ(t+1)T⋯Δδ(t+Nu)T]T, and other procedures are detailed in [9].Notice that matrixA(z-1) is diagonal; the problem is reduced to the recursion of n scalar Diophantine equations, which become much simpler to program and require less computation. The computation of G and f is also considerably simplified.Conventional GPC performance index can be written as follows:(15)J=y-wTQy-w+uTRu,where w is a future reference sequence and R and Q are positive definite weighting matrices.Remark 2.
CARIMA (controlled auto regressive integrated moving average) model [20] has been long used in GPC. The main reason we use the I-O model (CARIMA model) rather than state space model (7) is that we can directly use the input and output signals of the plant to avoid the use of a state observer, and the complexity of the controller’s front stage design is reduced. However, system identification technology is usually used in GPC to obtain the required I-O model (GPC can be considered as an adaptive control when the system identification is carried out online). Apparently, it will increase the difficulty of controller design by using system identification. State space model of AUV is widely used in many classic literatures [21, 22] and can be modeled from AUV dynamics equations. Therefore, it is wise and simple to obtain an I-O model by equivalent transformation of state space model.


## 3.2. Improved Performance Index

The energy consumption used for driving bow and stern rudders can be expressed as(16)E0=λ0·Δδ·Mr,where λ0 is selected according to energy consumption used for pitch stabilizing and depth tracking and Δδ is variation of rudder angle. Mr is the driving moment of rudder.Then the Lagrange dynamics equation is chosen to analyze flat wing-type rudder with thin foil condition as shown in Figure1. Rudder length is 2a, the distance from midpoint to shaft is c, and cross-section minor axis length is b.Figure 1
Flat wing-type rudder with thin foil condition.Lagrange functionL is defined as the difference between flat wing kinetic energy T and potential energy P. By setting the flap wing rotational angle δb,s(t), without considering the effects of the potential energy (under the thin foil condition, the quality of flat wing can be ignored, and the flat wing flaps within a small range), L can be expressed as(17)L=T=12Jδ˙b,s2t.The second Lagrange equation of flat wing can be described as(18)ddt∂L∂δ˙b,st-∂L∂δb,st=Mb,st.From (17) we can get (19)ddt∂L∂δ˙b,st=Jδ¨b,st,∂L∂δb,st=0.Substitute (19) into (18); then the driving moment of rudder can be obtained as follows:(20)Mb,st=Jδ¨b,st=Jω˙t,where J=πρ[a2c2+(a2-b2)2/8] is the additional moment inertia and ω(t) is the rotation angular velocity, without loss of generality, and for simplification of the following discussion, we set the control horizon Nu=2, and then(21)Mr=Mb,st+1TMb,st+2TT.Equation (21) and Δδ are substituted into (16); it follows that(22)E0=λ0Δδb,st+1TΔδb,st+2TMb,st+1Mb,st+2=λ0Δδb,st+1TMb,st+1+Δδb,st+2TMb,st+2,where Δδb,s(t+1)=δb,s(t+1)-δb,s(t)=ω(t+1+d)·Ts, Ts is sampling period and d is the delay between the angular velocity and driving rudder signal. Considering internal complicated mechanical structure of rudders, d is set as 1, and (23) is obtained:(23)ωt+2=Δδb,st+1Ts.Denote(24)Mb,st+1=Jω˙t+1=Jωt+2-ωt+1Ts.Substitute (23) into (24), and (25) is obtained:(25)Mb,st+1=JΔδb,st+1-Δδb,stTs2.Substituting (25) into (22), we obtain (26)E0=λnΔδTΔδ-F0Δδ,where λn=λ0J/Ts2, F0=c0c1, c0=Δδ(t) is the last step variation of rudder angle which is known, and c1=Δδ(t+1) is unknown.For a brief expression, (15) can be rewritten with (26) as follows: (27)J=12uTEu+uTFn+f0,s.t.Mu≤γ, where(28)E=2GTQG+R+λnI,Fn=F-F0,f0=f-wTf-w,F=2GTQf-w,u=Δδ.The constraints of the rudders deflection and pitch angle, corresponding to the inputs and outputs of the GPC, are characterized byM and γ and described next.


## 3.3. Consideration of Constraints

In practical application, considering the physical limitations of the driving device, and in order to ensure the correct operation of the AUV, the constraints of the rudders deflection and system output used to defineM and γ for each time step t can be defined as follows: (29)δmin≤δ≤δmax,Δδmin≤Δδ≤Δδmax,ymin≤y≤ymax,and then from (29) we can get(30)M1M2M3Δδ≤N1N2N3, where (31)M1=-C2C2,M2=-II,M3=-GG,N1=-δmin+C1δt-1δmax-C1δt-1,N2=-ΔδminΔδmax,N3=-Ymin+fYmax-f,C1=I2×2I2×2,C2=I2×2I2×20I2×2.For compactness of expression, we describe (30) by(32)MΔδ≤γ,which is equivalent to the constraints in (27).


## 3.4. Optimization with Equality Constraints

To minimize the quadratic function subject to (32), the performance index here becomes a QP (quadratic programming) problem. Let us consider the expression which contains the Lagrange multipliers, that is, a QP problem subject to equality constraints Mu=γ; we construct (33)J=12uTEu+uTFn+f0+λTMu-γ.The minimization ofJ is to take the first partial derivatives with respect to u and λ and make the results equal to zero. We obtain (34)∂J∂u=Eu+MTλ+∂uTFn∂u=0,(35)∂J∂λ=Mu-γ=0.Denote(36)uT=u11u12u21u22,Fn=F-F0=F1F2-c0c1=FcF2-c1=f11f12f21-u11f22-u12,uTFn=u11f11+u12f12+u21f21+u22f22-u21u11-u22u12.Then, in (34),(37)∂uTFn∂u=f11-u21f12-u22f21-u11f22-u12=FcF2-0010000110000100u=Fu-Su.The minimization ofJ can be made by finding the optimal u and λ via (34) and (35) where (38)λ=-ME-S-1MT-1γ+ME-S-1Fu,(39)u=-E-S-1MTλ+Fu=η-E-S-1MTλ,where η=-E-S-1Fu is the global optimal solution.


## 3.5. Optimization with Inequality Constraints

In (32), the inequality constraints may comprise active and inactive constraints. If Miu=γi, an inequality constraint Miu≤γi is said to be active, and if Miu<γi, it is inactive, where Mi and γi together form the ith inequality constraint. The Kuhn-Tucker conditions [13] are widely used to define the active and inactive constraints in terms of λ. If the active set were known, the original problem could become an equality constrains problem in (33).The conventional active set method in [9] belongs to the primal methods, in which the solutions are based on u (called decision variables). If the MIMO system has too many constraints, the calculation is quite large and it is not a straightforward work.To identify the constraints which are inactive systematically, a dual method can be used. 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 (39) in (33), the dual problem is written as (40)minλ≥012λTHλ+λTK+P1,where the matrices H and K are given by (41)P1=12FuTE-S-12Fn-Fu+λTME-S-1Fn-Fu-f0,H=ME-S-1MT,K=γ+ME-S-1Fu.The set ofλ which minimize the dual performance index(42)J=12λTHλ+λTK+P1,subject to λ≥0, are denoted as λ∗. Here Hildreth’s QP procedure [23] was used to solve the dual problem. The method can be written as (43)λit+1=max⁡0,βit+1,with(44)βit+1=-1hiiki+∑j=1i-1hijλjt+1+∑j=i+1Nkhijλjt,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 two sets of λ in one iterative cycle, λj(t+1) and λj(t). So, we set λj(0)=0, at t=1, and the iterative procedure converged to λ∗ as a result; by (39) we have (45)u=η-E-S-1MTλ∗,where λ∗=[β1(t+1),β2(t+1),…,βn(t+1)]T.According to the receding horizon control in GPC, the first two elements (bow and stern rudders deflection) inu are taken to form the Δu(t).Remark 3.
Because Hildreth’s QP is an element-by-element search based algorithm, there is no need to calculate any matrix inversion. However, if the active constraints are linearly dependent or the number of the active constraints is more than the number ofu, then λ will not converge to λ∗ and the iteration will terminate at the largest value of the iterative counter. But the algorithm will continue because there is no matrix inversion calculation. In this case, the algorithm will get a near-optimal solution with constraints if the constraints are violated. This is a main reason for using Hildreth’s QP in this paper, because of the ability to automatically recover from a deterioration constrained process.


## 3.6. Steps of Multivariable GPC with Constraints

(a)
Set values ofNP, Nu, and specify Q, R. Calculate matrix G and vector f; then E and F can be calculated. (b)
Sample current depthz(t) and pitch angle θ(t) and update constraints matrix M and γ by using δ(t-1). (c)
Check if the global optimal solutionη satisfies the constraints. If so, make λ∗ equal to zero vector and go to (e). If not, go to (d). (d)
Calculate matricesH and K; then the dual variable λ∗ can be calculated from (35). (e)
ReturnΔu(t) from optimal solution u. (f)
Go to step (b).


## 4. Simulation Results

In this paper, simulations are presented to demonstrate the effectiveness of multivariable GPC. The methods are used in a given depth control of DOLPHIN MARK II AUV which is developed by the University of British Columbia. The parameters about DOLPHIN can be found from [21]; here, the physical parameters which will be used in the simulation are shown in Tables 1 and 2.Table 1
Physical parameters of DOLPHIN AUV.
Description Values Units m 4390 kg Length 8.534 m Beam 1 m I x x 1315 Kg⋅m2 I y y 5900 Kg⋅m2 I z z 5057 Kg⋅m2 Location of CG ( 0,0 , 0 ) m Location of BG ( 0 . 086 , 0 , - 0 . 137 ) mTable 2
Parameters of rudders.
Description Type Chord Semispan Bow rudders NACA0025 0.35 m 0.968 m Stern rudders NACA0025 0.35 m 1.016 mCase 1.
For simulation, the parameters of GPC are chosen asNp=30, Nu=2, Q=INp×Np, and R=0.5I4×4. The initial values of state variables are all zero. The waves disturbance exposed on DOLPHIN 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 U=1.832(m/s), the desired depth is 3.5(m), and pitch angle is 0∘. Input and output constraints are(46)δbmax=-δbmin=30∘,δsmax=-δsmin=25∘,Δδb,smax=-Δδb,smin=5∘s,-5∘≤θ≤5∘.The purpose of the designed GPC with constraints is to achieve a given depth and pitch angle within the input and output constraints. Figures2 and 4 show the output and control signal of AUV with Hildreth’s QP procedure. It is seen that the depth converges to the desired value after 130 sampling instants, and pitch angle converges to zero after 190 sampling instants. And we note that the optimal control signal has large amplitude close to minimum limit but the output does not overshoot. From Figures 3 and 5, assume that the control signal amplitudes are limited directly without Hildreth’s QP procedure; when the control signals exceed their limits, the closed-loop performance deteriorates; as a result, the depth has a significant overshoot and the control signals become oscillatory after coming out of the saturation. So the simulation results above demonstrate that multivariable GPC with Hildreth’s QP procedure has the ability to deal with the control saturation and presents satisfactory performance of AUV in vertical plane.Figure 2
Output with Hildreth’s QP procedure.Figure 3
Output without Hildreth’s QP procedure.Figure 4
Control signal with Hildreth’s QP procedure.Figure 5
Control signal without Hildreth’s QP procedure.Figure6 shows the wave force and moment. The wave force and moment result from the calculation of (9), 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 of P-M spectrum t to superimpose the irregular waves.Figure 6
Wave force and moment.All these results demonstrate that AUV could achieve the desired depth and pitch angle under the wave disturbance. In addition, the control signals in this method are smooth and without control signal saturation.Case 2.
Energy consumption (EC) is expressed according to the following equation:(47)EC=2∑t=1NsimMb,stTΔδt=2JT2ΔδtTΔδt-c0.In order to demonstrate the energy-saving effectiveness of AUV in vertical plane, EC based on (15) is compared with EC based on (27) in the same conditions as shown in Figures 7 and 8. Figure 9 shows an EC comparison between GPC with and without Hildreth’s QP. Value results of simulation are shown in Table 3. The percentages ηe in parentheses show the differences between ECi and EC1. And the formula of percentage calculation is shown as follows:(48)ηe=ECi-EC1EC1×100%,where i=2,3, and the meaning of EC1, EC2, and EC3 is EC with E0 and Hildreth’s QP, without E0, without Hildreth’s QP, respectively.Table 3
Comparison of AUV diving EC.
Sampling instant EC1 (J) EC2 (J) EC3 (J) (ηe) (ηe) 30 2.7926 3.3425 (19.69%) 3.8365 (37.38%) 68 3.3593 4.7944 (29.93%) 5.8165 (73.15%) 100 14.6156 30.3467 (107.63%) 42.6749 (191.98%) 200 94.3112 165.1662 (75.13%) 216.9435 (129.71%) 500 372.6770 624.8213 (67.66%) 792.4392 (112.63%)Figure 7
EC comparison with and withoutE0.Figure 8
EC comparison of 0–100 sampling instants.Figure 9
EC comparison with and without Hildreth’s QP.It is seen from Figure8 that EC is constant when the AUV is diving because the control signals reach the minimum limit. It takes less EC to dive based on performance index with E0. Figure 7 shows the overall EC which includes the diving and fixed depth stage. From Table 3 we can see, at 68 sampling instants, EC increase as AUV begins to flap rudders to keep attitude constant. As shown in Figure 9 and Table 3, EC3 is particularly higher than EC1, because of the control saturation, different frequency, and amplitude of rudder angle shown in Figures 4 and 5, which lead to more frequent and significant stern rudder flapping.Moreover, the actual simulation time is 141 ms under the condition whereTs=0.1s and total sample instant is 500 (using Intel i5-4460 3.2 GHz PC). We note that the algorithm represented in this paper is fast enough within one sample instant (282 μs).The simulation results demonstrate the energy-saving effectiveness of proposed method that means multivariable GPC is an efficient and meaningful method for the controller design of AUV in vertical plane and for energy-saving.


## 5. Conclusions

In this paper, a multivariable GPC controller for AUV in vertical plane is presented, which controls the depth and pitch angle. This design uses the linearized AUV model and its augmented form to get the equivalent CARIMA model used in GPC. The proposed controller incorporates an energy consumption term in the performance index providing the system with significant energy-saving effectiveness. The design of the controller also takes into account the practical rudders deflection constraints and output constraints. By using a Hildreth’s QP procedure, the constraints can be simple handled. Making use of the proposed multivariable GPC controller, the AUV can navigate in vertical plane with different depth references and pitch angle. It is robust against rough wave disturbance near surface. The simulations carried out provide the validation of the proposed multivariable GPC controller, presenting fast dynamical response and energy consumption reduction. The utilization of carried energy is raised by means of driving rudders energy optimization.



---

*Source: 4650380-2016-11-23.xml*