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

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 onpitch 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 andminimized energy consumption improved by multivariable GPC.


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 2 Mathematical Problems in Engineering 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 energysaving and attitude stabilization.The rest of the paper is organized as follows.Section 2 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.

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: which can be simply expressed as discretized state space form An augmented state space model can be used if the model input is Δ() = ()−(−1), and Δ  () =   ()−  (−1).() = () − ( − 1) is a zero-mean white noise sequence.
Choosing a new state variable vector () = [Δ  ()  ()  ]  , we have the following: where  × is the identity matrix with dimensions ×, which is the number of outputs, and   is a  ×  1 zero matrix.In (6) where   ,   , , and  are matrices corresponding to the forms given in (7).In the following, the dimensionality of the augmented state space equation is taken to be  =  1 + .
Remark 1.The disturbance due to surface waves is the main disturbance type which is considered in this paper and written as ().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: where   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  0 .The transient force and moment are found by the use of Morison's equation [19] and integration along the length of the body  0 at each step as follows: where  is fluid density,   is the drag coefficient, and   is the added mass coefficient.In (9), where  is the equal interval frequency waveband and chosen in this thesis to discretize the wave spectrum,   is encountering frequency,   = √2(  ),   is wave number, and   is random phase shift (0 <   < 2) in relation to each frequency.

Multivariable GPC.
It is known that the GPC algorithm is based on using a CARIMA model where the unmeasurable disturbance is given by a white noise colored by ().
We consider the most usual case when ( −1 ) =  × , because the coloring polynomials are very difficult to estimate with sufficient accuracy in practice, especially in the multivariable case.
Multiply (11) by Δ With zero initial conditions, it can be seen that both ( 7) and ( 11) are equivalent if By making () = (), we can get the matrix A( −1 ) and B( −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 horizon   : where u = Δ = [Δ( + 1)  ⋅ ⋅ ⋅ Δ( +   )  ]  , and other procedures are detailed in [9].Notice that matrix A( −1 ) is diagonal; the problem is reduced to the recursion of  scalar Diophantine equations, which become much simpler to program and require less computation.The computation of  and f is also considerably simplified.
Conventional GPC performance index can be written as follows: where w is a future reference sequence and  and  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.

Improved Performance Index.
The energy consumption used for driving bow and stern rudders can be expressed as where  0 is selected according to energy consumption used for pitch stabilizing and depth tracking and Δ is variation of rudder angle.  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 Figure 1.Rudder length is 2, the distance from midpoint to shaft is , and cross-section minor axis length is .
Lagrange function  is defined as the difference between flat wing kinetic energy  and potential energy .By setting the flap wing rotational angle  , (), 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),  can be expressed as The second Lagrange equation of flat wing can be described as From ( 17) we can get Substitute ( 19) into (18); then the driving moment of rudder can be obtained as follows: where  = [ 2  2 + ( 2 −  2 ) 2 /8] is the additional moment inertia and () is the rotation angular velocity, without loss of generality, and for simplification of the following discussion, we set the control horizon   = 2, and then Equation ( 21) and Δ are substituted into (16); it follows that , ( + 2) ] =  0 [Δ , ( + 1)   , ( + 1) where Δ , ( + 1) =  , ( + 1) −  , () = ( + 1 + )⋅   ,   is sampling period and  is the delay between the angular velocity and driving rudder signal.Considering internal complicated mechanical structure of rudders,  is set as 1, and ( 23) is obtained: Denote Substitute ( 23) into (24), and ( 25) is obtained: Substituting (25) into (22), we obtain where is the last step variation of rudder angle which is known, and  1 = Δ( + 1) is unknown.
For a brief expression, (15) can be rewritten with (26) as follows: where The constraints of the rudders deflection and pitch angle, corresponding to the inputs and outputs of the GPC, are characterized by  and  and described next.

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 define  and  for each time step  can be defined as follows: and then from (29) we can get where For compactness of expression, we describe (30) by which is equivalent to the constraints in (27).

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 u = ; we construct The minimization of  is to take the first partial derivatives with respect to u and  and make the results equal to zero.We obtain Denote The minimization of  can be made by finding the optimal u and  via (34) and ( 35) where where  = −( − ) −1   is the global optimal solution.

Optimization with Inequality
Constraints.In (32), the inequality constraints may comprise active and inactive constraints.If   u =   , an inequality constraint   u ≤   is said to be active, and if   u <   , it is inactive, where   and   together form the th 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  (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 where the matrices  and  are given by The set of  which minimize the dual performance index 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 with where ℎ  is the th element in ,   is the th element in , and   is the number of rows of .In this method there are two sets of  in one iterative cycle,   ( + 1) and   ().So, we set   (0) = 0, at  = 1, and the iterative procedure converged to  * as a result; by (39) we have where  * = [ 1 ( + 1),  2 ( + 1), . . .,   ( + 1)]  .According to the receding horizon control in GPC, the first two elements (bow and stern rudders deflection) in u are taken to form the Δ().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 of , 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.(b) Sample current depth () and pitch angle () and update constraints matrix  and  by using ( − 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 matrices  and ; then the dual variable  * can be calculated from (35).
(f) Go to step (b).

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. surge speed is  = 1.832 (m/s), the desired depth is 3.5 (m), and pitch angle is 0 ∘ .Input and output constraints are The purpose of the designed GPC with constraints is to achieve a given depth and pitch angle within the input and output constraints.Figures 2 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 6 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.
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.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   in parentheses show the differences between EC and EC1.And the formula of percentage calculation is shown as follows: where  = 2, 3, and the meaning of EC1, EC2, and EC3 is EC with  0 and Hildreth's QP, without  0 , without Hildreth's QP, respectively.It is seen from Figure 8 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  0 .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 where   = 0.1 s 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.

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.

Figure 1 :
Figure 1: Flat wing-type rudder with thin foil condition.

Table 1 :
Physical parameters of DOLPHIN AUV.
(a) Set values of   ,   , and specify , .Calculate matrix  and vector ; then  and  can be calculated.

Table 3 :
Comparison of AUV diving EC.