Block diagonal dominant remotely operated vehicle model simulation using decentralized model predictive control

Abstract

Model predictive control on a highly coupled open-frame remotely operated vehicle system subjected to uncertain disturbances has always been a challenge. A decentralized model predictive control uses a design scaling to balance the interactions between the loops to achieve a block diagonal dominant remotely operated vehicle model is used. The numerical stability of the model predictive control algorithm improves despite the sensitivity of the control parameters on the output performance. The model predictive control gives a better two-error norm performance and more tuning options to control the velocity and position output as compared to other design scaling methods and other controllers such as proportional–integral–derivative, fuzzy logic, sliding mode, and proportional–integral plus proportional cascaded control when subjected to underwater current disturbance.

Introduction

Model predictive control (MPC) uses an explicit and a separately identifiable model for control. The concept of MPC algorithm is based on the moving horizon approach where the control action is computed to obtain the desired performance over a finite time horizon under uncertainties and constraints. ¹ The use of the constraints clearly distinguishes MPC from other control methods such as neural network–based approach,^2–8 leading to a more reliable controller and tighter control with no requirement on training data. However, the heavy online computational burden ⁹ and numerical condition of the control algorithm are the main obstacles in the application of MPC. A large computational delay and global minimum or even local minimum cannot be achieved due to time constraint in each time horizon. One approach is to obtain a closed-form optimal MPC^10,11 to circumvent the computational issue. However, it is hard to predict the system output over a long horizon since the output order is limited by the relative degree of a nonlinear model, and it can be quite unstable^12–14 for a large relative degree. With the derivatives of the control as an addition state, the robustness of the MPC under the model complexity such as highly coupled in states ¹⁵ and uncertainty when subjected to external disturbances pose an even greater challenge to the applied control community.

In the current literature, most of the applications of MPC were applied on the following systems such as ALSTOM gasifier control, ¹⁶ pH control, ¹⁷ obstacle avoidance for land robot, ¹⁸ pneumatic brake system, ¹⁹ trailing edge flapping on wind turbine blade, ²⁰ power converters, ²¹ solar energy system ²² and induction motor, ²³ and few autonomous underwater vehicles (AUVs).^24,25 The robust stability due to disturbances and the model uncertainties on the AUV were largely studied. However, the model used for the AUV control is quite decoupled in states due to its streamlined body design. The MPC on a highly coupled in states model such as a remotely operated vehicle (ROV) has not been studied. The inherent ROV’s model uncertainties due to hydrodynamic parameters, inertial nonlinearity, and coupling between the degrees of freedoms (DoFs) are worsened by the external disturbances due to underwater current. Applying the conventional MPC on this class of underwater vehicle creates a numerical stability issue on the control algorithm.

To overcome the numerical stability problem of MPC on a highly coupled ROV model, we propose the use of non-unity permutation matrices on ROV model before MPC system design to numerically condition and decentralize the model. It results in a block decentralized model for MPC. This article gives new simulation results of MPC laws and robust stability for a hover-capable ROV designed by Nanyang Technological University (NTU) in Singapore. The application of the MPC with two non-unity permutation matrices on a highly coupled system such as ROV model leads to a more decentralized predictive control system architecture design which makes the overall approach more decentralized.

The remainder of this article is organized as follows. Section “ROV model descriptions” gives the details of ROV model. Section “Block diagonal dominance model for ROV” describes the block diagonal dominance model for ROV. It is followed by the MPC design methodology in section “Decentralized MPC.” Section “Robust performance analysis” includes the robust performance analysis of ROV and comparison study of different controllers. Finally, the article ends with a conclusion.

ROV model descriptions

The ROV model used for control is the Robotics Research Centre (RRC) ROV^26,27 designed by NTU. The ROV plant model is quite coupled in the velocity states, namely, surge, sway, heave, roll, pitch, and yaw velocity in the dynamics equation. The RRC ROV (Figure 1) was designed to perform underwater pipeline inspections such as locating pipe leakages or cracks. The twin “eyeball” ROV depicted in Figure 1 has an open-frame structure. It measured 1 m long, 0.9 m wide, and 0.9 m high. It has a dry mass of 115 kg and a current operating depth of 100 m. The RRC ROV has only four thrusters to control 6 DoFs (i.e. surge, sway, heave, roll, pitch, and yaw velocity) and has a high level of cross coupling between the states. The roll and pitch motions are passive as the metacentric height is sufficient to provide adequate static stability. A brief description of the component layout of the RRC ROV is given:

• Four thrusters, each providing up to 70 N of thrust;

• Two cylindrical floats with four balancing steel weight;

• Main Pod (Pod 1) and sensors with navigational Pod (Pod 2);

• Two halogen lamps, an altimeter (for depth measurement), and a Doppler velocity log (for velocity measurement).

Figure 1.

Coordinate systems used in ROV (left) and actual ROV in pool (right).^26,27

Before the ROV modeling, the following assumptions were made:

• ROV is a rigid body and is fully submerged once in the water;

• Water is assumed to be ideal fluid that is incompressible, inviscid, and irrotational;

• ROV is slow moving during pipeline inspection;

• The earth-fixed frame of reference is inertial;

• Disturbance due to wave is neglected as it is designed for submerged operation;

• Tether dynamics attached to ROV is not modeled.

The motion of a rigid body on the body-fixed reference frame at the origin (Figure 1) is given by the equation ²⁸

$${\mathbf{M}}_{\mathbf{RB}}\stackrel{·}{\mathbf{v}}+{\mathbf{C}}_{\mathbf{RB}}(\mathbf{v})={\mathbf{\tau }}_{\mathbf{RB}}$$

(1)

where M_RB∈ℛ^6×6 is the mass inertia matrix, C_RB(v) ∈ℛ^6×6 is the Coriolis and centripetal matrix, τ_RB = [τ_RB1τ_RB2]^T∈ℛ⁶ is a vector of external forces and moments, v = [v₁ v₂]^T∈ℛ⁶ is the linear and angular velocity vector, namely: v₁ = [u v w]^T and v₂ = [p q r]^T, respectively.

The mass inertia matrix given in equation (1) can be written as follows

$${\mathbf{M}}_{\mathbf{RB}}=\left[\begin{array}{cccccc}m\hfill & 0\hfill & 0\hfill & 0\hfill & m{z}_G\hfill & -m{y}_G\hfill \\ \multicolumn{1}{c}{0}& m\hfill & 0\hfill & -m{z}_G\hfill & 0\hfill & m{x}_G\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & m\hfill & m{y}_G\hfill & -m{x}_G\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -m{z}_G\hfill & m{y}_G\hfill & I_x\hfill & -I_{\mathrm{xy}}\hfill & -I_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{m{z}_G}& 0\hfill & -m{x}_G\hfill & -I_{\mathrm{yx}}\hfill & I_y\hfill & -I_{\mathrm{yz}}\hfill \\ \multicolumn{1}{c}{-m{y}_G}& m{x}_G\hfill & 0\hfill & -I_{\mathrm{zx}}\hfill & -I_{\mathrm{zy}}\hfill & I_z\hfill \end{array}\right]$$

(2)

where x_G, y_G, and z_G are the coordinates of the center of gravity and m is the mass of the ROV. Here I_x, I_y, and I_z are the moments of inertia about the principal axes of the ROV, and I_xy = I_yx, I_xz = I_zx, and I_yz = I_zy are the products of the inertia. Here, the parameters used in equation (2) are as follows:

• $m=115.00\mathrm{kg},I_x=6.1000\mathrm{kg}{\mathrm{m}}^2,I_y=5.9800\mathrm{kg}{\mathrm{m}}^2,I_z=5.5170\mathrm{kg}{\mathrm{m}}^2$ ;

• $I_{\mathrm{xy}}=I_{\mathrm{yx}}=-0.0002\mathrm{kg}{\mathrm{m}}^2,I_{\mathrm{xz}}=I_{\mathrm{zx}}=-0.1850\mathrm{kg}{\mathrm{m}}^2,I_{\mathrm{yz}}=I_{\mathrm{zy}}=0.0006\mathrm{kg}{\mathrm{m}}^2$ .

The Coriolis and centripetal terms describing the angular motion of ROV can be expressed as follows

$${\mathbf{C}}_{\mathbf{RB}}(\mathbf{v})=\left[\begin{array}{cc}{0}_{3\times 3}\hfill & {\mathbf{C}}_{12}(v)\hfill \\ \multicolumn{1}{c}{-{\mathbf{C}}_{12}^{\mathrm{T}}(\mathbf{v})}& {\mathbf{C}}_{22}(\mathbf{v})\hfill \end{array}\right]$$

(3)

with

$${\mathbf{C}}_{12}(\mathbf{v})=\left[\begin{array}{ccc}m(y_{G}q+z_{G}r)\hfill & -m(x_{G}q-w)\hfill & -m(x_{G}r+v)\hfill \\ \multicolumn{1}{c}{-m(y_{G}p-w)}& m(z_{G}r+x_{G}p)\hfill & -m(y_{G}r-u)\hfill \\ \multicolumn{1}{c}{-m(z_{G}p-v)}& -m(z_{G}q-u)\hfill & m(x_{G}p+y_{G}q)\hfill \end{array}\right]$$

(4)

$${\mathbf{C}}_{22}(\mathbf{v})=\left[\begin{array}{ccc}0\hfill & -I_{\mathrm{yz}}q-I_{\mathrm{xz}}p+I_{z}r\hfill & I_{\mathrm{yz}}r+I_{\mathrm{xy}}p-I_{y}q\hfill \\ \multicolumn{1}{c}{I_{\mathrm{yz}}q+I_{\mathrm{xz}}p-I_{z}r}& 0\hfill & -I_{\mathrm{xz}}r-I_{\mathrm{xy}}q+I_{x}p\hfill \\ \multicolumn{1}{c}{-I_{\mathrm{yz}}r-I_{\mathrm{xy}}p+I_{y}q}& I_{\mathrm{xz}}r+I_{\mathrm{xy}}q-I_{x}p\hfill & 0\hfill \end{array}\right]$$

(5)

The nonlinear ROV dynamics equation and the kinematic equation can be expressed as follows

$$\mathbf{M}\stackrel{·}{\mathbf{v}}+\mathbf{C}(\mathbf{v})\mathbf{v}+\mathbf{Dv}+{\mathbf{G}}_{\mathbf{e}}(\mathbf{\eta })=\mathbf{\tau }$$

(6)

$$\stackrel{·}{\mathbf{\eta }}=\mathbf{J}({\mathbf{\eta }}_2)\mathbf{v}$$

(7)

where v = [v₁ v₂]^T = [ u v w p q r ]^T is the body-fixed velocity vector and η = [η₁η₂]^T is the earth-fixed vector, comprising the position vector η₁ = [ x y z ]^T and the orientation vector of Euler angles, η₂ = [ ϕ θ ψ ]^T. M = M_RB + M_A∈ℛ^6×6 is the sum of the rigid body inertia mass (i.e. a solid body with no deformation) and added mass (i.e. inertia added to moving ROV due to surrounding fluid being displaced) matrix, C(v) = C_RB(v) + C_A(v) ∈ℛ^6×6 is the sum of Coriolis and centripetal and the added mass forces and moments matrix, D ∈ℛ^6×6 is the hydrodynamic damping matrix, and G_e(η)∈ℛ⁶ is the gravitational and buoyancy vector. The propulsion forces and moments vector τ = Tu ∈ℛ⁶ relate the thrust output vector u = F_T $\overline{u}\in {\mathcal{R}}^4$ with the thruster configuration matrix T ∈ℛ^6×4, F_T∈ℛ^4×4 is the dynamics of each thruster and converts the input voltage command $\overline{u}\in {\mathbb{R}}^4$ into thrusts to propel the vehicle.

The dynamic model of the thruster was obtained experimentally using the test rig.^26,27 A simplified first-order thruster model was obtained as follows

$${\mathbf{F}}_{\mathbf{T}}=\frac{0.97f_T}{0.02\mathrm{s}+1}{\mathrm{I}}_4$$

(8)

where ${\mathrm{I}}_4$ is 4 × 4 identity matrix and

$$f_T=\{\begin{array}{cc}1.00\hfill & \mathrm{for\; input}=40\mathrm{V}\hfill \\ \multicolumn{1}{c}{0.82}& \mathrm{for\; input}=30\mathrm{V}\hfill \\ \multicolumn{1}{c}{0.60}& \mathrm{for\; input}=20\mathrm{V}\hfill \\ \multicolumn{1}{c}{0.33}& \mathrm{for\; input}=10\mathrm{V}\hfill \end{array}$$

As the ROV was designed to be neutrally buoyant by adding an additional float or balancing mass, the gravitational force due to the ROV weight can be made equal to the buoyancy force, $W=B=1128.15\mathrm{N}$ . By adding an extra mass on the ROV at an appropriate location, XY coordinates of the center of buoyancy coincide with the center of gravity, x_G = x_B = 0, y_G = y_B = 0. The gravitational and buoyancy vector becomes

$${\mathbf{G}}_{\mathbf{e}}(\mathbf{\eta })=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{(z_G-z_B)W\cos \theta \sin \phi }\\ \multicolumn{1}{c}{(z_G-z_B)W\sin \theta }\\ \multicolumn{1}{c}{0}\end{array}\right]$$

(9)

where $z_G-z_B=-0.03\mathrm{m}$ .

Since the ROV is symmetric about the XZ plane and close to symmetric about YZ plane, we assumed that the motions in surge, sway, pitch, and yaw are decoupled. Although it is not completely symmetric about the XY plane, the ROV is operating at a relatively low speed in which the coupling effects ²⁸ can be negligible. With this assumption, the damping matrix in equation (6) becomes

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{D}=-\mathrm{diag}\{X_u,Y_v,Z_w,K_p,M_q,N_{\mathrm{r}}\}}\\ =-\mathrm{diag}\{3.70,2.91,2.37,1.66,1.45,0.11\}\hfill \end{array}$$

(10)

For most low-speed maneuvering tasks, the off-diagonal terms in the added mass matrix are neglected ²⁸ in most applications. ${\mathbf{M}}_{\mathbf{A}}$ is simplified to a diagonal form as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{M}}_{\mathbf{A}}=-\mathrm{diag}\{X_{\stackrel{·}{u}},Y_{\stackrel{·}{v}},Z_{\stackrel{·}{w}},K_{\stackrel{·}{p}},M_{\stackrel{·}{q}},N_{\stackrel{·}{r}}\}}\\ =-\mathrm{diag}\{21.14,51.70,92.45,3.61,2.64,2.30\}\hfill \end{array}$$

(11)

Since the off-diagonal elements in ${\mathbf{M}}_{\mathbf{A}}$ are neglected, the Coriolis and centripetal added mass matrix ${\mathbf{C}}_{\mathbf{A}}(\mathbf{v})$ becomes

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{C}}_{\mathbf{A}}(\mathbf{v})=\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & -Z_{\stackrel{·}{w}}w\hfill & Y_{\stackrel{·}{v}}v\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & Z_{\stackrel{·}{w}}w\hfill & 0\hfill & -X_{\stackrel{·}{u}}u\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -Y_{\stackrel{·}{v}}v\hfill & X_{\stackrel{·}{u}}u\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -Z_{\stackrel{·}{w}}w\hfill & Y_{\stackrel{·}{v}}v\hfill & 0\hfill & -N_{\stackrel{·}{r}}r\hfill & M_{\stackrel{·}{q}}q\hfill \\ \multicolumn{1}{c}{Z_{\stackrel{·}{w}}w}& 0\hfill & -X_{\stackrel{·}{u}}u\hfill & N_{\stackrel{·}{r}}r\hfill & 0\hfill & -K_{\stackrel{·}{p}}p\hfill \\ \multicolumn{1}{c}{-Y_{\stackrel{·}{v}}v}& X_{\stackrel{·}{u}}u\hfill & 0\hfill & -M_{\stackrel{·}{q}}q\hfill & K_{\stackrel{·}{p}}p\hfill & 0\hfill \end{array}\right]}\\ =\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & -92.45w\hfill & 51.7v\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 92.45w\hfill & 0\hfill & -21.14u\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -51.70v\hfill & 21.1403u\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -92.45w\hfill & 51.70v\hfill & 0\hfill & -2.30r\hfill & 2.64q\hfill \\ \multicolumn{1}{c}{92.45w}& 0\hfill & -21.14u\hfill & 2.30r\hfill & 0\hfill & -3.61p\hfill \\ \multicolumn{1}{c}{-51.70v}& 21.14u\hfill & 0\hfill & -2.64q\hfill & 3.61p\hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

(12)

A series of experimental tests were performed to verify the hydrodynamic damping and added mass coefficients. The complete details of the numerical modeling, simulation, and experimental tests on the RRC ROV model can be found in Chin et al.^26,27

Linear ROV subsystem model

A linearized model of the ROV operating at station-keeping condition is obtained by linearizing the nonlinear ROV in equations (6) and (7) about an equilibrium point

$${\mathbf{v}}_{\mathbf{o}}(t)=[\begin{array}{cccccc}{u}_o(t)\hfill & v_o(t)\hfill & w_o(t)\hfill & p_o(t)\hfill & q_o(t)\hfill & r_o(t)\hfill \end{array}{]}^T$$

(13)

$${\mathbf{\eta }}_{\mathbf{o}}(t)=[\begin{array}{cccccc}{x}_o(t)\hfill & y_o(t)\hfill & z_o(t)\hfill & \phi (t)\hfill & {\theta }_o(t)\hfill & {\psi }_o(t)\hfill \end{array}{]}^T$$

(14)

For brevity, the linearization procedure using Taylor series will be omitted. Defining $\mathbf{x}=\mathbf{[}{\begin{array}{cc}{\mathbf{x}}_1^{\mathbf{T}}\hfill & {\mathbf{x}}_2^{\mathbf{T}}\mathbf{]}\hfill \end{array}}^{\mathbf{T}}$ , we obtain the following state-space model

$$\begin{array}{l}\left[\begin{array}{c}{\stackrel{.}{x}}_1\\ {\stackrel{.}{x}}_2\end{array}\right]=\stackrel{A(t)}{\overbrace{\left[\begin{array}{cc}-M^{-1}[C(t)+D(t)]& -M^{-1}{G}_e(t)\\ J(t)& J^{\ast }(t)\end{array}\right]}}\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]\\ +\stackrel{B}{\overbrace{\left[\begin{array}{c}{M}^{-1}T\\ 0\end{array}\right]}}u\end{array}$$

(15)

It can be written in abbreviated form as follows

$$\stackrel{·}{\mathbf{x}}=\mathbf{A}\mathbf{(}t\mathbf{)}\mathbf{x}+\mathbf{Bu}$$

(16)

During the maneuvering in a vertical plane, the steady-state linear and angular components are assumed as v_o = p_o = q_o = r_o = 0 and the equilibrium point is defined by the zero roll and pitch angles, that is ϕ_o = θ_o = 0. Hence, the time-varying matrices for the vertical plane motion $(u_o,w_o)$ with a fixed heading angle ${\psi }_o$ are simplified as shown

$$\mathbf{M}=\left[\begin{array}{cccccc}m-X_{\stackrel{·}{u}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& m-Y_{\stackrel{·}{v}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & m-Z_{\stackrel{·}{w}}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & I_{\mathrm{xx}}-K_{\stackrel{·}{p}}\hfill & -I_{\mathrm{xy}}\hfill & -I_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -I_{\mathrm{xy}}\hfill & I_{\mathrm{yy}}-M_{\stackrel{·}{q}}\hfill & -I_{\mathrm{yz}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -I_{\mathrm{xz}}\hfill & -I_{\mathrm{yz}}\hfill & I_{\mathrm{zz}}-N_{\stackrel{·}{r}}\hfill \end{array}\right]$$

(17)

$$\mathbf{D}=-\left[\begin{array}{cccccc}{X}_u\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& Y_v\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & Z_w\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & K_p\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & M_q\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & N_r\hfill \end{array}\right]$$

(18)

$$\mathbf{C}=\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & (m-Z_{\stackrel{·}{w}})w_o\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -(m-Z_{\stackrel{·}{w}})w_o\hfill & 0\hfill & (m-X_{\stackrel{·}{u}})u_o\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & -(m-X_{\stackrel{·}{u}})u_o\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& (m-Z_{\stackrel{·}{w}})w_o\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-(m-Z_{\stackrel{·}{w}})w_o}& 0\hfill & (m-X_{\stackrel{·}{u}})u_o\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -(m-X_{\stackrel{·}{u}})u_o\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(19)

$${\mathbf{G}}_{\mathbf{e}}=\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & -z_{B}B\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & -z_{B}B\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(20)

$$\mathbf{T}=\left[\begin{array}{cccc}1\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0.707\hfill & -0.707\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0.707\hfill & 0.707\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -0.293\hfill & 0.293\hfill \\ \multicolumn{1}{c}{-0.016}& -0.016\hfill & -0.012\hfill & -0.012\hfill \\ \multicolumn{1}{c}{0.310}& -0.310\hfill & 0.012\hfill & -0.012\hfill \end{array}\right]$$

(21)

If we assumed that ${\psi }_o$ is constant and ϕ_o = θ_o = 0, the kinematic transformation matrix J takes the form

$$\mathbf{J}=\left[\begin{array}{cccccc}\cos ({\psi }_o)\hfill & -\sin ({\psi }_o)\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\sin ({\psi }_o)}& \cos ({\psi }_o)\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(22)

Consequently, the linear time-invariant model for both velocity ${\mathbf{x}}_1$ and position ${\mathbf{x}}_2$ can be written in state-space form

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mathbf{x}}=\mathbf{Ax}+\mathbf{Bu}}\\ \multicolumn{1}{c}{\left[\begin{array}{c}{\stackrel{·}{\mathbf{x}}}_1\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathbf{x}}}_2}\end{array}\right]=\left[\begin{array}{cc}-{\mathbf{M}}^{-1}\mathbf{[}\mathbf{C}+\mathbf{D}\mathbf{]}\hfill & -{\mathbf{M}}^{-1}{\mathbf{G}}_{\mathbf{e}}\hfill \\ \multicolumn{1}{c}{\mathbf{J}}& 0\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{x}}_2}\end{array}\right]+\left[\begin{array}{c}{\mathbf{M}}^{-1}\mathbf{T}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]\mathbf{u}}\end{array}$$

(23)

where A and B are constant matrices.

During the station keeping, the controller is tasked to maintain the vehicle about the equilibrium position. The linearization about the equilibrium point $u_o=0.5\mathrm{m}/\mathrm{s},w_o=0.5\mathrm{m}/\mathrm{s},{\psi }_o=0.75\mathrm{rad}$ can be seen in the linear time-invariant model in equation (23). For velocity control, the first sixth-order state-space model is used. By regulating the velocity of the ROV, the position of the ROV can be indirectly controlled. The position of the ROV can be obtained via Euler’s transformation of the velocity states. By including the ROV state-space equation (with m = 6 states) and the thruster dynamic $\mathbf{u}={\mathbf{F}}_{\mathrm{T}}\overline{\mathbf{u}}$ (with q = 4 states), the ROV has n = 10 state as shown

$$\stackrel{·}{\mathbf{x}}={\mathbf{A}}_{11}\mathbf{x}+{\mathbf{B}}_1\overline{\mathbf{u}}$$

(24)

where ${\mathbf{A}}_{11}\in {\Re }^{10\times 10}$ is the system matrix, ${\mathbf{B}}_1\in {\Re }^{10\times 6}$ is the input matrix, $\mathbf{x}\in {\Re }^{10}$ is the state vector, and $\overline{\mathbf{u}}\in {\Re }^6$ is the desired input vector. The output equation can be expressed as follows

$$\mathbf{y}={\mathbf{C}}_{11}\mathbf{x}$$

(25)

where ${\mathbf{C}}_{11}$ is a 6 × 6 identity output matrix and $\mathbf{y}\in {\Re }^6$ is the velocity output vector.

Block diagonal dominance model for ROV

Before control system design, the inherent properties such as open-loop stability, right half plane (RHP) zero, and system interactions are examined. Since the RRC ROV is unlikely to exceed $u_o=0.5\mathrm{m}/\mathrm{s},w_o=0.5\mathrm{m}/\mathrm{s},{\psi }_o=1.57\mathrm{rad}$ , the open-loop RRC ROV system is evaluated at the following equilibrium: $u_o=0.5\mathrm{m}/\mathrm{s},w_o=0.5\mathrm{m}/\mathrm{s},{\psi }_o=0.75\mathrm{rad}$ . An interaction test using Rosenbrock’s row diagonal dominance with Gershgorin disks superimposed on the diagonal elements of the system frequency response was performed to show the coupling between the states of the RRC ROV. As observed in Figure 2, it indicates that the system is highly interactive over all frequencies (1–100 rad/s). It is observed from the increasing size of the disks in the diagonal elements of each row. The growing size of the disk implies that the off-diagonal terms are more dominant than its diagonal terms. To minimize the control effort required for a highly coupled system, pre-conditioning (i.e. to obtain a more diagonal dominance) on the ROV system is favorable to facilitate controller design. Note that to achieve diagonal dominance, either row or column dominance is sufficient.

Figure 2.

Direct Nyquist Array (DNA) with Gershgorin disks before scaling.

The corresponding row and column dominance measure (i.e. the sum of the off-diagonal terms in each row or column over the absolute value of the transfer function evaluated at each frequency) shows that the magnitudes are indeed not small. The two-norm values are approximately 4.9 × 10³ and 3.3 × 10³ for the row and column dominance, respectively. It further confirmed that the RRC ROV is highly coupled in motions. Also, the singular value and condition number at different frequencies plot shows that the ROV is highly coupled in motions. As shown in Figure 3, there are many peaks, in particularly at approximately 49.95 rad/s. A design scaling is therefore required to reduce the coupling and to improve the numerical condition via a more decentralized structure before MPC layout.

Figure 3.

Condition number at different frequencies.

Various design scaling methods such as Edmunds ¹² scaling, one-norm, and Perron–Frobenius (PF) scaling are used. The scaling frequency at 10 rad/s is chosen due to its lowest condition number among the scaling methods. The outputs are a non-unity diagonal pre- and post-compensator with different values. The Edmunds scaling algorithm (see Appendix 1 for the algorithm written in MATLAB) provides both scaling and input–output pairings as seen in the non-unity permutation matrices. Both the columns and rows are now reordered into the following sequence: surge, pitch, heave, yaw, sway, and roll

$${\mathbf{S}}_{\mathbf{pre}}^{\mathbf{N}}=\mathrm{diag}\{[\begin{array}{cccccc}19.68\hfill & 16.20\hfill & 4.324\hfill & 39.10\hfill & 213.3\hfill & 4.143\hfill \end{array}]\}$$

(26)

$${\mathbf{S}}_{\mathbf{post}}^{\mathbf{N}}=\mathrm{diag}\{[\begin{array}{cccccc}49.91\hfill & 22.52\hfill & 55.92\hfill & 2.451\hfill & 35.68\hfill & 8.666\hfill \end{array}]\}$$

(27)

$${\mathbf{S}}_{\mathbf{pre}}^{\mathbf{PF}}=\mathrm{diag}\{[\begin{array}{cccccc}0.429\hfill & 84.66\hfill & 0.409\hfill & 197.2\hfill & 24.85\hfill & 5.063\hfill \end{array}]\}$$

(28)

$${\mathbf{S}}_{\mathbf{post}}^{\mathbf{PF}}=\mathrm{diag}\{[\begin{array}{cccccc}0.001\hfill & 7.382\hfill & 0.006\hfill & 0.700\hfill & 0.004\hfill & 0.659\hfill \end{array}]\}$$

(29)

$${\mathbf{S}}_{\mathbf{post}}^{\mathbf{E}}=\left[\begin{array}{cccccc}49.91\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 35.68\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 55.92\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 8.66\hfill \\ \multicolumn{1}{c}{0}& 22.52\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 2.45\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(30)

$${\mathbf{S}}_{\mathbf{pre}}^{\mathbf{E}}=\left[\begin{array}{cccccc}19.68\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 16.20\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 4.32\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 39.10\hfill \\ \multicolumn{1}{c}{0}& 213.32\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 4.14\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(31)

The row dominance ratio of each scaling method is compared. It was found that the PF scaling (row dominance ratio, 159; column dominance ratio, 334) has the smaller dominance ratio as compared to the Edmunds scaling (row dominance ratio, 323; column dominance ratio, 497) and one-norm (row dominance ratio, 323; column dominance ratio, 497) scaling method. However, the one-norm design scaling gives the lowest condition number (or less ill-conditioned) over the frequencies as shown in Figure 4. This is also reflected in the less concentration of Gershgorin disks at the origin of the plot in Figure 5.

Figure 4.

Condition number at different frequencies.

Figure 5.

Direct Nyquist Array with Gershgorin disks after one-norm scaling.

However, the transfer function matrix (TFM) of the ROV is unlikely to be made diagonal dominant by this simple diagonal constant scaling matrix. It can be verified in Figure 6 as the gain in dB over the frequencies of interest are all greater than 6 dB.

Figure 6.

Perron–Frobenius eigenvalues of ROV system.

The Edmunds scaling using the non-unity permutation matrices is therefore chosen. It helps to numerically condition and reduces the interactions between the input-to-output interactions as it is not always possible to have a clear diagonal dominance at each row or column. The proof of the numerical stability can be seen in Edmunds. ¹² As shown in Figure 7, the ROV system can be grouped into two blocks of $3\times 3$ subsystems. Instead of looking at each element in each row (column) for row (column) diagonal dominance, a group of elements in a block can be tested for diagonal dominance. This is known as the block diagonal dominance measure of the TFM

$${\mathbf{G}}_{\mathbf{ij}}^{\mathbf{k}\times \mathbf{k}}(\mathbf{s})=\left[\begin{array}{cc}{\mathbf{G}}_{\mathbf{ii}}^{3\times 3}\hfill & {\mathbf{G}}_{\mathbf{ij}}^{3\times 3}\hfill \\ \multicolumn{1}{c}{{\mathbf{G}}_{\mathbf{ji}}^{3\times 3}}& {\mathbf{G}}_{\mathbf{jj}}^{3\times 3}\hfill \end{array}\right]$$

(32)

where superscript $3\times 3$ and subscript $i,j$ refer to the size of the partitioned block and elements in the block, respectively. The matrix $\mathbf{G}\mathbf{(}\mathbf{s}\mathbf{)}$ is considered to be block diagonal dominant if the gain of the diagonal blocks dominates the gain of the off-diagonal blocks. To test the block diagonal dominance of $\mathbf{G}\mathbf{(}\mathbf{s}\mathbf{)}$ , the gains of the sub-blocks of ${\mathbf{G}}_{\mathbf{ii}}$ are determined as follows

$$\begin{array}{c}\underset{¯}{\sigma }({\mathbf{G}}_{\mathbf{ii}})=\underset{x\ne 0}{min}\frac{\Vert {\mathbf{G}}_{\mathbf{ii}}\mathbf{x}\Vert }{\Vert \mathbf{x}\Vert }\hfill \\ =\mathrm{smallest\; singular\; value\; of}{\mathbf{G}}_{\mathbf{ii}}\hfill \end{array}$$

(33)

Figure 7.

Decoupled step response.

To verify that the RRC ROV system is now block diagonal dominant, Figure 8 shows the block dominance measure of the system before and after the Edmunds scaling using $\underset{¯}{\sigma }(\mathbf{G})$ . As shown, the ROV model becomes more block diagonal dominance. The block dominance measure of TFM shows a decrease from 1.6 × 10⁻² to 3.5 × 10⁻⁵. It implies that the TFM of the ROV is now decoupled into two blocks of each matrix of size $3\times 3$ .

Figure 8.

Block diagonal dominance measure before (left) and after (right) Edmunds scaling.

After achieving a better numerically conditioned system, another important part of a multivariable design is the input–output pairing. Applying the relative gain array (RGA) to the decentralized ROV system at the steady-state condition as $\mathbf{\Lambda }(0)=\mathbf{G}\mathbf{*}(\mathbf{G}{)}^{-\mathbf{T}}$ . Here, the operation * denotes element-by-element multiplication (Hadamard or Schur product). Since the RGA is applied at the steady state, or zero frequency, it is important to examine the RGA at other frequencies. An RGA number by Bristol is used and has the following form

$$\mathrm{RGA}\mathrm{number}=\Vert \mathbf{\Lambda }(\mathbf{G})-\mathbf{I}{\Vert }_{\mathrm{sum}}$$

(34)

By evaluating the RGA number across the frequencies of interest after the row permutation (by Edmunds scaling), there is a slight decrease in the RGA number after ordering across the frequencies as seen in Figure 9. It implies that the selected I/O pairs have an effect on the diagonal dominance of the system. A lower RGA number corresponds to lower condition number of the ROV system. It is a “good” feature of the control system design.

Figure 9.

RGA number of ROV.

The proposed Edmunds scaling algorithm on the TFM sets all the row and column sums to unity using simple iterative steps. The algorithm equalizes and maximizes the geometric mean of the row dominance and geometric mean of the column dominance for given inputs and outputs. It improves the numerical stability in MPC system design where the solution of the algebraic Riccati equation is needed. Also, the design scaling makes the individual loops more independent by balancing the interactions between the loops. As a result, the scaled state-space model of ROV becomes

$$\stackrel{·}{\mathbf{x}}={\mathbf{A}}_{11}\mathbf{x}+{\mathbf{B}}_1{\mathbf{S}}_{\mathbf{pre}}\overline{\mathbf{u}}$$

(35)

$$\mathbf{y}={\mathbf{S}}_{\mathbf{post}}{\mathbf{C}}_{11}\mathbf{x}$$

(36)

where ${\mathbf{S}}_{\mathbf{pre}}\in {\Re }^{6\times 6}$ is the pre-compensator and ${\mathbf{S}}_{\mathbf{post}}\in {\Re }^{6\times 6}$ is the post-compensator using Edmunds design scaling method in equations (30) and (31), respectively.

Decentralized MPC

In the research literature, MPC is formulated in state space. The MPC is obtained by minimizing the cost function by application of the quadratic programming algorithms ¹⁴ and the relevant references cited in this text. The model of the ROV is described by the linear difference equations (24) and (25). As the input $\overline{\mathbf{u}}$ cannot affect the output $\mathbf{y}$ at the same time, D = 0 is assumed. The following auxiliary variables are defined to include the state-space model in an augmented state-space model

$$\mathbf{z}(t)=\stackrel{·}{\mathbf{x}}(t)$$

(37)

$$\mathbf{y}(t)={\mathbf{C}}_{11}\mathbf{x}(t)$$

(38)

With a new state vector ${\mathbf{x}}_{\mathbf{a}}(t)=[\mathbf{z}(t{)}^{\mathrm{T}}\mathbf{y}(t{)}^{\mathrm{T}}{]}^{\mathrm{T}}$ , the augmented state-space model is defined as follows

$$\left[\begin{array}{c}\stackrel{·}{\mathbf{z}}(t)\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\mathbf{y}}(t)}\end{array}\right]=\stackrel{\mathbf{A}}{\overbrace{\left[\begin{array}{cc}{\mathbf{A}}_{11}\hfill & 0_{n\times m}\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_{\mathbf{post}}{\mathbf{C}}_{11}}& 0_{m\times m}\hfill \end{array}\right]}}\left[\begin{array}{c}\mathbf{z}(t)\hfill \\ \multicolumn{1}{c}{\mathbf{y}(t)}\end{array}\right]+\stackrel{\mathbf{B}}{\overbrace{\left[\begin{array}{c}{\mathbf{B}}_1{\mathbf{S}}_{\mathbf{pre}}\hfill \\ \multicolumn{1}{c}{0_{m\times m}}\end{array}\right]}}\stackrel{·}{\overline{\mathbf{u}}}(t)$$

(39)

$$\mathbf{y}(t)=\stackrel{\mathbf{C}}{\overbrace{\left[\begin{array}{cc}{0}_{m\times n}\hfill & {\mathbf{I}}_{m\times m}\hfill \end{array}\right]}}\left[\begin{array}{c}\mathbf{z}(t)\hfill \\ \multicolumn{1}{c}{\mathbf{y}(t)}\end{array}\right]$$

(40)

where ${\mathbf{I}}_{m\times m}$ is the identity matrix m × m; $0_{m\times m}$ and $0_{m\times n}$ are zero matrices of dimension m × m and m × n, respectively.

The information of the state vector at a time $t_i$ is written as $\mathbf{x}(t_i)$ . The state feedback control with gain matrix K is valid for $0\le \tau \le T_p$ where $\tau$ is the time variable within the window from $t_i$ to prediction horizon $T_p$ . The derivative of the control within this window is therefore written as follows

$$\stackrel{·}{\overline{\mathbf{u}}}(\tau )=-\mathbf{Kx}(\tau )$$

(41)

The predicted future state at time $\tau$

$$\mathbf{x}(t_i+\tau |t_i)=e^{\mathbf{A}\tau }\mathbf{x}(t_i)+\mathbf{\Phi }(\tau {)}^{\mathrm{T}}\mathbf{\eta }$$

(42)

where

$$\mathbf{\Phi }(\tau {)}^{\mathrm{T}}=\underset{0}{\overset{\tau }{\int }}{e}^{\mathbf{A}(\tau -\gamma )}[\begin{array}{cccc}{\mathbf{B}}_1{\mathbf{L}}_1{(\gamma )}^{\mathrm{T}}\hfill & {\mathbf{B}}_2{\mathbf{L}}_2{(\gamma )}^{\mathrm{T}}\hfill & \dots \hfill & {\mathbf{B}}_{\mathbf{m}}{\mathbf{L}}_{\mathbf{m}}{(\gamma )}^{\mathrm{T}}\hfill \end{array}]d\gamma ;{\mathbf{\eta }}^{\mathrm{T}}=\left[\begin{array}{cccc}{\mathbf{\eta }}_1^{\mathrm{T}}\hfill & {\mathbf{\eta }}_2^{\mathrm{T}}\hfill & \dots \hfill & {\mathbf{\eta }}_{\mathrm{m}}^{\mathrm{T}}\hfill \end{array}\right]$$

(43)

$${\mathbf{L}}_{\mathrm{i}}(\tau {)}^{\mathrm{T}}=[\begin{array}{cccc}{l}_1^i(\tau )\hfill & l_2^i(\tau )\hfill & \dots \hfill & l_{N_p}^i(\tau )\hfill \end{array}]=e^{{\mathbf{A}}_{\mathbf{p}}\tau }\mathbf{L}(0)$$

(44)

$${\mathbf{A}}_{\mathbf{p}}=\left[\begin{array}{cccc}-{\rho }_A\hfill & 0\hfill & \dots \hfill & 0\hfill \\ \multicolumn{1}{c}{-2{\rho }_A}& -{\rho }_A\hfill & \dots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \dots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{-2{\rho }_A}& \cdots \hfill & \cdots \hfill & -2{\rho }_A-{\rho }_A\hfill \end{array}\right]$$

(45)

$$\mathbf{L}(0)=\sqrt{2{\rho }_A}{\left[\begin{array}{cccc}1\hfill & 1\hfill & \dots \hfill & 1\hfill \end{array}\right]}^{\mathrm{T}}$$

(46)

where ${\rho }_A$ is the time scaling factor and $N_p$ is the Laguerre parameter for each input.

The predicted output at a time $\tau$ can be computed as follows

$$\mathbf{y}(t_i+\tau |t_i)=\mathbf{C}[e^{\mathbf{A}\tau }\mathbf{x}(t_i)+\mathbf{\Phi }(\tau {)}^{\mathrm{T}}\mathbf{\eta }]$$

(47)

where

$$\mathbf{\Omega }=\underset{0}{\overset{T_p}{\int }}[\mathbf{\Phi }(t_i)\mathbf{Q}\mathbf{\Phi }(t_i{)}^{\mathrm{T}}+\mathbf{R}]d\tau ;\mathbf{Q},\mathbf{R}\ge 0$$

(48)

$$\mathbf{\Psi }=\underset{0}{\overset{T_p}{\int }}\mathbf{\Phi }(t_i)\mathbf{Q}{e}^{\mathbf{A}\tau }d\tau$$

(49)

The optimal $\mathbf{\eta }$ becomes

$$\mathbf{\eta }=-{\mathbf{\Omega }}^{-1}\mathbf{\Psi }\mathbf{x}(t_i)$$

(50)

that minimizes $\mathbf{J}$ as

$${\mathbf{J}}_{min}=\mathbf{x}(t_i{)}^{\mathrm{T}}[\underset{0}{\overset{T_p}{\int }}{e}^{{\mathbf{A}}^{\mathrm{T}}\tau }\mathbf{Q}{e}^{\mathbf{A}\tau }d\tau -{\mathbf{\Psi }}^{\mathrm{T}}{\mathbf{\Omega }}^{-1}\mathbf{\Psi }]\mathbf{x}(t_i)$$

(51)

The optimal control signal is written as follows

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\overline{\mathbf{u}}}(t)=\left[\begin{array}{cccc}{L}_1{(0)}^{\mathrm{T}}\hfill & 0_2\hfill & \dots \hfill & 0_m\hfill \\ \multicolumn{1}{c}{0_1}& L_2{(0)}^{\mathrm{T}}\hfill & \dots \hfill & 0_m\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \dots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0_1}& 0_2\hfill & \cdots \hfill & L_m{(0)}^{\mathrm{T}}\hfill \end{array}\right]\left[\begin{array}{c}{\eta }_1\hfill \\ \multicolumn{1}{c}{{\eta }_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\eta }_m}\end{array}\right]}\\ =-\mathbf{L}(\tau {)}^{\mathrm{T}}{\mathbf{\Omega }}^{-1}\mathbf{\Psi }\mathbf{x}(t)\hfill \\ =-{\mathbf{K}}_{\mathbf{mpc}}\mathbf{x}(t)\hfill \end{array}$$

(52)

The first term is the derivative of the input state $\mathbf{x}(t)$ , while the error signal $\mathbf{y}(t)-\mathbf{r}(t)$ is the second term of the state feedback as shown

$$\stackrel{·}{\overline{\mathbf{u}}}(t)=-\stackrel{{\mathbf{K}}_{\mathbf{mpc}}}{\overbrace{[\begin{array}{cc}{\mathbf{K}}_{\mathbf{x}}\hfill & {\mathbf{K}}_{\mathbf{y}}\hfill \end{array}]}}\left[\begin{array}{c}\stackrel{·}{\mathbf{x}}(t)\hfill \\ \multicolumn{1}{c}{\mathbf{y}(t)-\mathbf{r}(t)}\end{array}\right]$$

(53)

where $\mathbf{r}(t)$ is the set-point signal.

The control signal computation is carried out by simply approximating the derivative $\stackrel{·}{\overline{\mathbf{u}}}(t)$ with

$$\stackrel{·}{\overline{\mathbf{u}}}(t_i)\approx \frac{\overline{\mathbf{u}}(t_i)-\overline{\mathbf{u}}(t_i-\Delta t)}{\Delta t}\Rightarrow \overline{\mathbf{u}}(t_i)=\overline{\mathbf{u}}(t_i-\Delta t)+\stackrel{·}{\overline{\mathbf{u}}}(t_i)\Delta t$$

(54)

The actual control signal $\overline{\mathbf{u}}(t)$ will update itself repeatedly from its past value

$$\begin{array}{c}\multicolumn{1}{c}{\overline{\mathbf{u}}(t_i)={\overline{\mathbf{u}}}_{\mathbf{actual}}(t_i-\Delta t)+\stackrel{·}{\overline{\mathbf{u}}}(t_i)\Delta t}\\ ={\overline{\mathbf{u}}}_{\mathbf{actual}}(t_i-\Delta t)\hfill \\ +\left[\begin{array}{cccc}{L}_1{(0)}^{\mathrm{T}}\hfill & 0_2\hfill & \dots \hfill & 0_m\hfill \\ \multicolumn{1}{c}{0_1}& L_2{(0)}^{\mathrm{T}}\hfill & \dots \hfill & 0_m\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \dots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0_1}& 0_2\hfill & \cdots \hfill & L_m{(0)}^{\mathrm{T}}\hfill \end{array}\right]\left[\begin{array}{c}{\eta }_1\hfill \\ \multicolumn{1}{c}{{\eta }_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\eta }_m}\end{array}\right]\Delta t\hfill \end{array}$$

(55)

The optimal solution with the exponentially weighted cost function is identical to the original optimal solution without exponential weighting when the prediction horizon is sufficiently large. However, without modification on the pair of weight matrices $\mathbf{Q}$ and $\mathbf{R}$ , the solution does not guarantee exponential decay of the original variable, $\mathbf{x}(t_i+\tau |t_i)$ within the horizon. The weight matrix $\mathbf{Q}$ becomes ${\mathbf{Q}}_{\mathbf{\alpha }}=\mathbf{Q}+2(\alpha +\beta )\mathbf{P};\alpha ,\beta >0,\mathbf{R}>0$ . $\mathbf{P}$ is the solution of the steady-state Riccati equation remains unchanged. The optimal control $\stackrel{·}{\overline{\mathbf{u}}}(\tau )$ is thus obtained. The proposed approach is numerically sound as it allows the use of a sufficiently large prediction horizon using ${\mathbf{A}}_{\mathbf{\alpha }}=\mathbf{A}-\alpha \mathbf{I}$ to guarantee asymptotical stability.

The control amplitude constraints for the state-space model at a time t can be written in the form where ${\overline{\mathbf{u}}}_{\mathbf{min}}$ and ${\overline{\mathbf{u}}}_{\mathbf{max}}$ are vectors containing the required lower and upper limits of the control amplitude, respectively

$${\overline{\mathbf{u}}}_{\mathbf{min}}\le \overline{\mathbf{u}}(t)\le {\overline{\mathbf{u}}}_{\mathbf{max}}$$

(56)

The constraints on $\stackrel{·}{\overline{\mathbf{u}}}(t)$ take the form

$${\stackrel{·}{\overline{\mathbf{u}}}}_{\mathbf{min}}\le \stackrel{·}{\overline{\mathbf{u}}}(t)\le {\stackrel{·}{\overline{\mathbf{u}}}}_{\mathbf{max}}$$

(57)

where ${\stackrel{·}{\overline{\mathbf{u}}}}_{\mathbf{min}}$ and ${\stackrel{·}{\overline{\mathbf{u}}}}_{\mathbf{max}}$ correspond to the minimum and maximum changes in the thruster rate, respectively. The constraints are applied to the system inputs on both the absolute limits and the speed of change. The rate of the thrust input is chosen to be lower than the maximum rate (this is an estimate as the thrust is proportional to square of its speed). If the controller requires a significant change of thrust beyond the thruster specification, then the thruster dynamics will perform like a low-pass filter to improve the controller performance. The rate of change of set points for each thruster is set to ±0.01 N/s. The maximum (minimum) absolute thrust value is limited to 70 N (−70 N), corresponding to a thruster speed of approximately 1000 r/min.

The implementation of MPC law also requires direct measurement of all elements in $[\stackrel{·}{\mathbf{x}}(t)\mathbf{y}(t)-\mathbf{r}(t){]}^{\mathrm{T}}$ at each time instant. The state vector component $\stackrel{·}{\mathbf{x}}(t)$ can be estimated. However, there will be noise in $\mathbf{x}(t)$ causing the computed difference $\stackrel{·}{\mathbf{x}}(t)$ to amplify the noise. Thus, an observer is used to estimate $\stackrel{·}{\mathbf{x}}(t)$ as

$$\widehat{\stackrel{·}{\mathbf{x}}}(t)=\mathbf{A}\hat{\mathbf{x}}(t)+\mathbf{B}\stackrel{·}{\overline{\mathbf{u}}}(t)+{\mathbf{K}}_{\mathbf{ob}}[\mathbf{y}(t)-\mathbf{C}\hat{\mathbf{x}}(t)]$$

(58)

where $\hat{\mathbf{x}}(\mathbf{t}\mathbf{)}$ is the estimate of $\mathbf{x}(t)$ and ${\mathbf{K}}_{\mathbf{ob}}$ is the observer gain matrix. ${\mathbf{K}}_{\mathbf{ob}}$ is designed using discrete-linear quadratic regulator theory with weighting matrices ${\mathbf{Q}}_{\mathbf{ob}}$ and ${\mathbf{R}}_{\mathbf{ob}}$ . ${\mathbf{K}}_{\mathbf{ob}}$ is selected such that the closed-loop observer error system matrix is stable. The effects of measurement noise in $\mathbf{y}(t)$ can be reduced by choosing a large value of ${\mathbf{R}}_{\mathbf{ob}}$ . The initial state vector of the observer is set to zero for output vector $\mathbf{y}(t)=0$ .

In summary, the steps to design the decentralized MPC with Edmunds scaling are as follows:

1. Determine the non-unity permutation matrices ( ${\mathbf{S}}_{\mathbf{post}}^{\mathbf{E}}$ and ${\mathbf{S}}_{\mathbf{pre}}^{\mathbf{E}}$ ) of Edmunds scaling in equations (30) and (31) to achieve least interaction and block diagonal dominance on the TFM (32);

2. Form the augmented matrix for MPC using Edmunds scaling in equations (39) and (40);

3. Set parameters for continuous MPC-time scaling parameter $({\rho }_A)$ , number of Laguerre parameters $(N_p)$ , prediction horizon $(T_p)$ , and plane shifting parameters on matrix $\mathbf{A}(\alpha ,\beta >0)$ ;

4. Solve Riccati equation for $\mathbf{P}$ using ${\mathbf{Q}}_{\mathbf{\alpha }}=\mathbf{Q}+2(\alpha +\beta )\mathbf{P}>0,\mathbf{R}>0$ and ${\mathbf{A}}_{\alpha }=\mathbf{A}-\mathbf{\alpha }\mathbf{I}$ ;

5. Define constraint $\mathbf{L}(\tau )$ using equation (44) and input constraints on thrusters: ${\overline{\mathbf{u}}}_{\mathbf{min}}\le \overline{\mathbf{u}}(t)\le {\overline{\mathbf{u}}}_{\mathbf{max}}$ ;

6. Solve MPC for $\mathbf{\Omega }$ and $\mathbf{\Psi }$ using equations (48) and (49), respectively;

7. Compute observer gain matrix, ${\mathbf{K}}_{\mathbf{ob}}$ with weighting matrices ${\mathbf{Q}}_{\mathbf{ob}}$ and ${\mathbf{R}}_{\mathbf{ob}}$ ;

8. Compute the derivative state-feedback gain for $\stackrel{·}{\overline{\mathbf{u}}}(t)=-{\mathbf{K}}_{\mathbf{mpc}}\mathbf{x}(t)$ where ${\mathbf{K}}_{\mathbf{mpc}}=\mathbf{L}(\tau {)}^{\mathrm{T}}{\mathbf{\Omega }}^{-1}\mathbf{\Psi }$ using equation (52);

9. Compute the state feedback gain $\overline{\mathbf{u}}(t)$ ;

10. Plot $\hat{\mathbf{y}}(t)$ and control input $\overline{\mathbf{u}}(t)$ .

Robust performance analysis

As the ROV is subjected to uncertainties, its robustness performance needs to be examined. The ROV subjected to changing underwater current effect will be simulated. It is followed by comparing the results to various scaling methods. For completeness, the proposed MPC will be compared with different types of controllers under a fixed and varying controller gains.

Robust analysis on ROV subjected to underwater current

The robustness analysis of the decentralized MPC on the ROV is performed. The uncertain perturbations are formed in a block diagonal matrix of stable perturbations, $\Delta$ where each element is real. The common $\mathbf{N}\mathbf{\Delta }$ structure is used to analyze the ROV system (Figure 10), $\mathbf{P}$ with MPC controller, ${\mathbf{K}}_{\mathbf{mpc}}$ and uncertainty. The lower linear fractional transformation (LTF) is denoted by

$$\mathbf{N}={\mathbf{P}}_{11}+{\mathbf{P}}_{12}{\mathbf{K}}_{\mathbf{mpc}}(\mathbf{I}-{\mathbf{P}}_{22}{\mathbf{K}}_{\mathbf{mpc}}{)}^{-1}{\mathbf{P}}_{21}$$

(59)

Figure 10.

Control system structure for uncertain ROV system.

The uncertain closed-loop transfer function is related to $\mathbf{N}$ and $\mathbf{\Delta }$ by an upper LTF

$$\mathbf{F}={\mathbf{N}}_{22}+{\mathbf{N}}_{21}\mathbf{\Delta }\mathbf{(}\mathbf{I}-{\mathbf{N}}_{11}\mathbf{\Delta }{)}^{-1}{\mathbf{N}}_{12}$$

(60)

As shown in Figure 10, the ROV closed-loop system has multiplicative input uncertainty, ${\mathbf{\Delta }}_{\mathbf{I}}$ . Here, ${\mathbf{W}}_{\mathbf{p}}$ is a performance weight and ${\mathbf{W}}_{\mathbf{I}}$ is a normalized weight for uncertainty. The generalized ROV plant is written as follows

$$\mathbf{P}=\left[\begin{array}{ccc}0\hfill & 0\hfill & {\mathbf{W}}_{\mathbf{I}}\hfill \\ \multicolumn{1}{c}{{\mathbf{W}}_{\mathbf{p}}\mathbf{G}}& {\mathbf{W}}_{\mathbf{p}}\hfill & {\mathbf{W}}_{\mathbf{p}}\mathbf{G}\hfill \\ \multicolumn{1}{c}{-\mathbf{G}}& -\mathbf{I}\hfill & -\mathbf{G}\hfill \end{array}\right]$$

(61)

The partition of $\mathbf{P}$ becomes

$$\begin{array}{l}{P}_{11}=\left[\begin{array}{cc}0& 0\\ W_{p}G& W_p\end{array}\right];P_{12}=\left[\begin{array}{c}{W}_I\\ W_{p}G\end{array}\right];\\ P_{21}=\left[\begin{array}{cc}-G& -I\end{array}\right];P_{22}=-G\end{array}$$

(62)

Equation (59) becomes

$$\mathbf{N}=\left[\begin{array}{cc}-{\mathbf{W}}_{\mathbf{I}}\mathbf{KG}\mathbf{(}\mathbf{I}+\mathbf{KG}{)}^{-1}\hfill & -{\mathbf{W}}_{\mathbf{I}}\mathbf{K}\mathbf{(}\mathbf{I}+\mathbf{GK}{)}^{-1}\hfill \\ \multicolumn{1}{c}{{\mathbf{W}}_{\mathbf{p}}\mathbf{G}\mathbf{(}\mathbf{I}+\mathbf{KG}{)}^{-1}}& {\mathbf{W}}_{\mathbf{p}}{(\mathbf{I}+\mathbf{GK}\mathbf{)}}^{-1}\hfill \end{array}\right]$$

(63)

The weighting function is used to shape the sensitivity of the closed-loop system, $\mathbf{N}$ to the desired level. A typical type of performance function used is a low-pass filter, high-pass filter, or constant weight. The tuning of the weighting function, irrespective of the kind used, is performed iteratively:

• ${\mathbf{W}}_1$ . From the frequency response of the disturbance (Figure 11), the disturbance is dominant up to 0.0001 rad/s. To reduce the disturbance up to this frequency, ${\mathbf{W}}_1$ is selected to be a low-pass filter at the frequency range

  $${\mathbf{W}}_1=\frac{1\times {10}^{-5}s+1\times {10}^{-6}}{s+1\times {10}^{-4}}{\mathbf{I}}_6$$

  (64)

• ${\mathbf{W}}_{\mathbf{p}}$ . A high-pass filter was chosen to shape the input sensitivity of the closed-loop system. A similar first-order high-pass filter was used in each channel with a corner frequency of around 0.001 rad/s, to limit the input magnitudes at high frequencies and thereby limit the closed-loop bandwidth as seen in Figure 11

  $${\mathbf{W}}_{\mathbf{p}}=\frac{1\times {10}^{-2}s+2.5\times {10}^{-7}}{1\times {10}^{5}s+25}I$$

  (65)

Figure 11.

Weighting function for output performance and input.

For robust stability, the following conditions for the ROV in the uncertainty set need to satisfy for robust performance:

• N is internally stable (all eigenvalues of N is negative as seen in Figures 12 and 13);

• $\Vert \mathbf{N}\Vert =2\times {10}^{-7}<1$ ; moreover, $\mathbf{N}$ is internally stable;

• F is stable $\forall \Delta =\mathrm{diag}\{[0.4850.4850.780.750.770.74]\}$ , $\Vert \Delta {\Vert }_{\infty }=0.78\le 1$ , and $\mathbf{N}$ are internally stable;

• $\Vert \mathbf{F}\Vert =2\times {10}^{-7}<1$ , $\forall \Delta =\mathrm{diag}\{[0.4850.4850.780.750.770.74]\}$ , $\Vert \Delta {\Vert }_{\infty }=0.78\le 1$ , and $\mathbf{N}$ are internally stable.

Figure 12.

Sensitivity plot of $\mathbf{N}$ .

Figure 13.

Pole and zeros of closed-loop ROV system.

As shown in Figure 13, all eigenvalues (including $\mathbf{N}$ ) lie on the left-hand plane. The closed-loop ROV model with random input disturbance $\Vert \Delta {\Vert }_{\infty }=0.78\le 1$ is stable. In addition, the singular value of the sensitivity function $\mathbf{S}=(1+{\mathbf{GK}\mathbf{)}}^{-1}$ is small that indicates the closed-loop system with a disturbance signal that has less influence on the system output. The system has a good stability margin. As seen in Figure 12, the input sensitivity $\mathbf{KS}=\mathbf{K}\mathbf{(}1+\mathbf{GK}{)}^{-1}$ plot, the effect of the input disturbances is quite negligible on the output response. The robust stability due to parametric uncertainty is below the margin of 0 dB. It indicates the ROV system is robust against parametric uncertainty. As observed in $\mathbf{T}=\mathbf{KG}\mathbf{(}1+\mathbf{KG}{)}^{-1}$ plot, the gain over the frequencies is less than the 0 dB $(=1)$ margin. It indicates that the system output will be quite small for any significant inputs.

The current velocity is included as a relative velocity, ${\mathbf{v}}_{\mathbf{r}}=\mathbf{v}-{\mathbf{v}}_{\mathbf{c}}$ where ${\mathbf{v}}_{\mathbf{c}}=[\begin{array}{cccccc}{u}_c\hfill & v_c\hfill & w_c\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}{]}^{\mathrm{T}}$ is the current velocity vector acting on the ROV’s body-fixed frame. The Earth-fixed current velocity vector is given by ${\mathbf{v}}_{\mathbf{c}}^{\mathbf{E}}=[\begin{array}{cccccc}{u}_c^E\hfill & v_c^E\hfill & w_c^E\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}{]}^{\mathrm{T}}$ , while the body-fixed velocity can be determined by the following Euler’s transformation

$$u_c^E=u_c\cos {\alpha }_c\cos {\beta }_c$$

(66)

$$v_c^E=v_{\mathbf{c}}\sin {\beta }_c$$

(67)

$$w_c^E=w_{\mathbf{c}}\sin {\alpha }_c\cos {\beta }_c$$

(68)

where α_c is the angle of attack and β_c is the sideslip angle.

During the simulation, both angles of attack and sideslip angle are set to 45°. The current speed was incremented by 1 m/s until 10 m/s (for testing the robustness of the system). However, these settings can be adjusted to reflect the actual current experienced by the ROV. As shown in Figure 14, the H-infinity norm of the linear velocity increases with the sea current. The increase in the output linear velocities begins gradually with less impact on the angular velocities (as observed from the smaller magnitude in the H-infinity norm). The ROV system is quite robustly stable under the underwater current disturbances for underwater current not more than 3 m/s. As observed in Figure 14, the hinf-norm increases exponentially after 3 m/s.

Figure 14.

H-infinity norm of output velocity (subjected to underwater current).

Robust analysis on ROV using various scaling methods

The time response simulation of the decentralized MPC on the ROV was performed under the MATLAB/Simulink™ environment. The parameters used in the scaling methods and MPC are varied to study the robust performance of each method. From the simulation results, the parameters employed in the simulation were determined through multiple trials till the system becomes unstable. From the previous simulations, the control horizon $T_p$ has the least significant effect on ROV performance. When a higher value of time scaling factor ${\rho }_A$ is used, all DoFs exhibit a lower steady-state error. Although higher Laguerre parameter $N_p$ allows the ROV to perform better, the computational time increases linearly as control horizon $T_p$ and ${\rho }_A$ increase. There is no significant improvement in the output response by increasing $T_p$ , $\alpha$ , and $\beta$ . However, the stability of the response can be affected by increasing the difference between the $\beta$ and $\alpha$ . The response in $\mathbf{\Delta }\overline{\mathbf{u}}$ and $\mathbf{\Delta }\stackrel{·}{\overline{\mathbf{u}}}$ becomes quite oscillatory although it shows some convergence to its steady-state value. For a high $\beta$ value (i.e. by shifting the eigenvalues to further RHP), the steady-state error increases and the velocity responses become quite unstable.

As observed, the output responses are constant with higher $\alpha$ (i.e. by shifting the eigenvalues further into the left-hand plane). There is no effect of using $\beta$ to improve the numerical condition of the Riccati equation. Hence, a lower $\beta$ was adopted such that $\beta -\alpha$ is kept small. However, with a lower $\alpha (<\beta )$ , the output responses improve with less oscillatory in motion and steady-state error. As a result of multiple simulations, the chosen parameters are as follows:

• One-norm scaling: ${\rho }_A=10,N_p=6,T_p=4,{\overline{\mathrm{u}}}_{\max }=70,\Delta {\overline{\mathrm{u}}}_{\max }=0.01,\beta =0.2,\alpha =0.005$ ;

• Edmunds scaling: ${\rho }_A=10,N_p=6,T_p=4,{\overline{\mathrm{u}}}_{\max }=70,\Delta {\overline{\mathrm{u}}}_{\max }=0.01,\beta =0.01,\alpha =0.005$ ;

• PF scaling: ${\rho }_A=10,N_p=6,T_p=4,{\overline{\mathrm{u}}}_{\max }=70,\Delta {\overline{\mathrm{u}}}_{\max }=0.01,\beta =0.01,\alpha =0.005$ ;

• Without scaling: ${\rho }_A=10,N_p=6,T_p=4,{\overline{\mathrm{u}}}_{\max }=70,\Delta {\overline{\mathrm{u}}}_{\max }=0.01,\beta =0.01,\alpha =0.005$ .

Figures 15–17 show the comparison results of using different scaling methods when the ROV is moving at $u=v=1\mathrm{m}/\mathrm{s}$ while trying to maintain a zero velocity for the remaining DoF. For Edmunds scaling, the surge and sway velocity converge quickly with minimal overshoot as compared to other scaling methods. However, there is steady-state error in the surge direction. As compared to the one-norm and PF scaling design method, the steady-state error is much higher. Due to the constraints on ${\overline{\mathrm{u}}}_{\max }$ and $\Delta {\overline{\mathrm{u}}}_{\max }$ , Edmunds scaling gives the least control effort. It can be recognized by the decoupled control structure of the ROV as a result of the non-unity permutation matrices used. When compared to the one-norm and PF scaling method, the rate of the control effort required is much larger. The velocity response of the PF scaling method is observed to be smaller than the rest. It is because the amount of control effort required is much lower. By examining the MPC design without any scaling, the velocity responses are unable to achieve the desired steady-state values.

Figure 15.

Velocity time response of various design scaling methods.

Figure 16.

Control input rate of various design scaling methods.

Figure 17.

Control input of various design scaling methods.

All scaling methods exhibit quite a high sway velocity (but within the constraint) as compared to surge direction. However, as the vehicle moves diagonally to the set point, there is a higher sway position response as compared to the surge direction. Fortunately, the remaining DoFs are not excited as the ROV system matrix is now less coupled. It shows that the ROV system has a faster response to the desired set point target under the input constraints. If the ROV encountered a significant disturbance in the inputs, the maximum thrust could help the ROV to maintain its closed-loop stability and reject the disturbance quickly.

The simulated results provide evidence that the decentralized model predictive controller using Edmunds design scaling is capable to control both surge and sway velocities to the desired set point. The roll, pitch, and yaw rate are quite small in the design scaling methods.

Comparisons with different controllers under fixed control parameters

MPC has not been applied to a highly coupled open-frame ROV system subjected to uncertain disturbance. The conventional control schemes such as proportional–integral–derivative (PID), ²⁹ fuzzy logic, ³⁰ sliding-mode control (SMC), ³¹ and proportional–integral plus proportional (P-PI) cascaded control ²⁶ are compared with the proposed decentralized MPC. Note that this is not the complete list of controllers used in the current literature. However, it serves to compare with the proposed MPC. For consistency, the same input commands, $u=v=1\mathrm{m}/\mathrm{s}$ on ROV model are used. The sea current disturbance was set to 2 m/s for testing the robustness of all the controllers. As observed in Figure 18, the PID controller tries to regulate the desired velocity, but the response is quite oscillatory with high pitch rate.

Figure 18.

Proposed decentralized MPC comparisons with other control schemes.

On the other hand, the cascaded control scheme with PI-P gives less oscillatory behavior. However, the overshoot is quite significant as compared to the conventional PID single-loop control structure. Fuzzy logic control exhibits less oscillatory motion but exhibits steady-state error of approximately 50% with oscillatory in the angular velocities and a prominent response of 1 m/s in heave motion. SMC has demonstrated quite a good steady-state value. It can attenuate the disturbance with a small steady-state error. Finally, the proposed decentralized MPC exhibits good overall response and is capable of achieving a good steady-state value under the disturbances as compared to the controllers. As observed, the roll and yaw motion are quite stable for most control schemes. As shown in Table 1, the proposed MPC controller with a decentralized scheme using Edmunds scaling performs slightly better (as seen in the two-norm and infinity norm of the velocity response) than other controllers.

Table 1.

Two and infinity norm comparison between the controller’s velocity response.

Control schemes	Parameters	Two-norm	Infinity-norm
PID control	$k_p=31,k_i=11,k_d=0.5$	234.7	39.10
MPC-Edmunds scaling	${\rho }_A=10,N_p=6,T_p=4,{\overline{\mathrm{u}}}_{\max }=70$ $\beta =0.01$ $\alpha =1$	9.900	3.400
Fuzzy-Sugeno control		16.10	6.900
SMC	$k_s=0.01,k_v=20,k_d=0.05$ BL = 0.05, lambda = 0.01	10.60	6.200
Cascaded control	$k_p=31$ $k_p=31,k_i=6$	45.90	371.5

PID: proportional–integral–derivative; MPC: model predictive control; SMC: sliding-mode control; BL: boundary layer.

Comparisons with different controllers under varying control parameters

The effects of different controller gains on the ROV’s output velocity responses are examined. We reviewed the change in only one controller gain by a fixed increment at any one time. The two-norm errors of the velocity outputs are then computed. The decentralized MPC using the Edmunds scaling is first examined. Initially, ${\rho }_A=10,N_p=6,T_p=4$ , $\beta =0.01$ , and $\alpha =1$ are used. As shown in Figure 19, the parameters such as ${\rho }_A,N_p$ , $T_p$ , $\beta$ , and $\alpha$ are plotted against the two-norm errors for each velocity. It can be seen that the decentralized MPC is quite sensitive to changes in the control parameters such as $N_p,\beta$ , and $\alpha$ relatively less impact on $T_p$ and ${\rho }_A$ . However, the computational time increases when a larger Laguerre parameter $N_p$ is used. As compared to PID controller in Figure 20, the control parameters have less influence on the output velocities. In SMC shown in Figure 21, the increased control parameters $k_d$ and boundary layer (BL) thickness have some impacts on the velocity responses. It has a reasonable two-norm error as compared to the smallest error norm in the proposed decentralized MPC. As compared to the cascaded P-PI control scheme in Figure 22, the impacts on the controller gains can be seen by the large proportional gain used in both the inner and outer loops. The two-norm errors are the highest as compared to the above-mentioned control schemes. In summary, each controller has its advantages and disadvantages. The proposed MPC can be quite sensitive to control parameters changes, but at the same time, it gives more tuning options for better performance as reflected by the lowest two-norm error for certain set of control parameters.

Figure 19.

Two-norm error of proposed decentralized MPC comparisons at different control parameters.

Figure 20.

Two-norm error of PID controller at different control parameters.

Figure 21.

Two-norm error of SMC at different control parameters.

Figure 22.

Two-norm error of P-PI cascaded control at various control parameters.

Conclusion

The robust decentralized MPC of a ROV is proposed. The proposed Edmunds design scaling algorithm on the TFM equalizes and maximizes the geometric mean of the row dominance for given inputs and outputs. The design scaling makes the individual control loops more independent, hence able to achieve a decentralized ROV model for MPC. The simulation results show that the decentralized MPC using Edmunds scaling exhibits a better performance as compared to other design scaling methods.

The proposed decentralized MPC behaves quite well in the velocity response as compared with PID, fuzzy logic, SMC, and P-PI cascaded control scheme. The simulated results have shown some evidences that the decentralized model predictive controller using Edmunds design scaling is capable of controlling the ROV. From the results, the decentralized MPC algorithm can be implemented with a less numerical error. Although the proposed decentralized MPC with Edmunds scaling is quite sensitive to some control parameters, it offers more options to tune the ROV’s output performance to achieve a lower two-norm error.

Future works will include using a time-varying nonlinear ROV model subjected to parametric and multiple uncertain external disturbances. Comparisons with more controllers such as adaptive neural network control will be performed on ROV system.