Predictive Control for a Wave-Energy Converter Array Based on an Interconnected Model

Abstract

This paper proposes a model predictive control (MPC) method based on an interconnected model to maximize the ocean wave energy extracted by a wave-energy converter (WEC) array. In the proposed method, a formally uniform interconnected model is applied to represent the dynamics of an array consisting of an arbitrary quantity of WECs, simultaneously considering the hydrodynamic interaction among all the WEC devices. First, the WEC devices and their hydrodynamic interaction are represented in an interconnected model that describes the networked dynamics of a variety of WEC arrays with distinct spatial geometry layout of the WEC devices deployed in the sea field. Second, based on the presented model, an MPC method is applied to achieve the coordinated control of the WEC array to improve its energy conversion efficiency under the constraints of buoy position and control force. Third, a hardware-in-the-loop (HIL) platform is developed to simulate the WEC array’s physical operating conditions, and the proposed method is implemented on the platform to test its performance. The test results show that the proposed MPC method using the interconnected model has a higher energy harvesting efficiency than the classic MPC method.

1. Introduction

Using ocean wave energy is a novel approach to solving the problem of resource depletion and environmental damage associated with traditional energy sources. However, the exploitation of ocean wave energy is not as well commercialized as other sustainable energy sources, such as solar and wind energy. One of the key reasons is that ocean wave energy is distributed homogeneously on the ocean surface, which makes it difficult for any wave-energy converter (WEC) to harvest a great deal of wave energy by operating a single machine over a small range. Using an array system composed of multiple point-absorber-type WECs can expand the range of wave energy harvesting options. Another advantage of WEC arrays is that they effectively reduce and eliminate power fluctuations, so that a smoother and more stable power supply for power grids and power consumers can be provided [1]. In recent years, WEC arrays have been intensively investigated due to their advantages [2].

The existence of hydrodynamic interactions among WEC devices has a significant impact on the performance of a WEC array. In [3], a wave interaction theory, including the diffraction interaction of evanescent waves, was presented to enable the calculation of the wave hydrodynamics of a multimember structure with arbitrary configurations and spacings. The effect of interactions within basic two-unit WEC arrays with different separating distances was studied in [4]. Subsequently, the influence of interactions among bodies on the overall yearly energy production of an array was assessed [5]. A CFD motion solver based on a coupled model was developed in [6] to analyze the interactions ocurring in a WEC array. Due to diffracted and radiated waves, interactions among the WECs can occur which may be constructive or destructive [7]. If such hydrodynamic interactions are effectively utilized, they have the potential to increase the energy production efficiency dramatically.

Optimizing the geometry layout allows for the effective use of constructive interaction effects. A topological method was presented in [7] to enable the design of an optimized single-row wave power plant established by a triangle array unit. An improved differential evolution algorithm was proposed in [8] to optimize the layout of a WEC array for faster speed of convergence and improved accuracy. A genetic algorithm that included a cost model was proposed in [9] to enable consideration of economic factors. For the example introduced in the paper, the distance between the devices is reduced as much as possible to diminish the maintenance cost. However, to be more efficient for wave energy harvesting in particular sea environments, it is necessary to specify an optimized spatial geometry layout of a WEC array corresponding to the ocean wave features and conditions in the area in which it is located. Therefore, a uniform modeling approach that can represent arbitrary layouts of WEC arrays is crucial for the layout optimization of WEC arrays.

To maximize the wave energy harvesting of wave energy by a WEC, efficient control of the array is necessary. A comparison between MPC and causal control was made from the perspective of performance and constraint handling in [10], with MPC being found to have significant advantages over optimal causal control. Since MPC has shown remarkable control performance when dealing with the constraints in realistic control processes, it has received increasing attention in recent years [11]. In the early literature, researchers focused mainly on the application of MPC in the control of a single WEC. A real-time non-linear model predictive controller (NMPC) was utilized to increase the WEC efficiency of a two-degrees-of-freedom WEC with highly non-linear power take-off (PTO) characteristics in [12]. An economic model predictive control (EMPC) strategy with a general economic cost function that directly reflects the economic objective of the system was proposed in [13] to maximize the energy extracted from ocean waves and to minimize the operational cost. Two MPC control strategies for a two-body point absorber were designed to extract the maximum energy from waves through unidirectional power flow in [14]. A non-quadratic piecewise discontinuous functional was applied as a cost index of NMPC to maximize the energy produced by a WEC in [15]. Various MPC strategies based on a single-unit model were successfully applied to maximize the wave energy harvested by WECs [16]. MPC methods based on WEC arrays have also attracted the attention of researchers. A model predictive control method was developed in [17], where each MPC effectively optimizes the energy absorption of its own WEC under state constraints. Nevertheless, because of the non-negligible interference generated by the movement of unit buoys, if the WEC units in the WEC array are closely spaced, implementing an MPC with an isolated unit model can result in degraded MPC control performance caused by model mismatch. Therefore, the proposed model not only takes into account the individual devices, but also the complex coupling relationships between the WEC devices, represented as hydrodynamic interactions among each pair of WEC devices. Centralized MPCs were described in [18] which led to a potential 10% improvement in converted energy. These studies investigated the impact of interference between WECs operating close to each other, which can increase wave energy harvesting efficiency. The majority of the studies reported derive sets of equations representing interference relationships directly from a particular geometry (e.g., equilateral triangles) [19]. It is challenging to directly extend the methods discussed to the modeling of WEC arrays with diverse spatial geometry layouts.

Since hydrodynamic interactions vary with local wave conditions in the ocean, it is unrealistic for a fixed model to maintain accuracy all the time. Therefore, in future commercial systems, not only does the system’s layout need to be regularly adjusted, but the model also needs to be regularly modified and adjusted to ensure the performance of MPC. This paper proposes a modeling method to rapidly adapt to arbitrary geometric layouts of a WEC array by considering a WEC array with coupled interconnection relationships. The model can be broken down into two parts: the dynamics of the WEC unit and the dynamics of the interfering effects of the coupled interconnection relationships. Generalization of the model for the WEC array can thus be addressed through the modeling method. This method reduces modeling difficulty and improves MPC survival adaptability for practical applications. We design an MPC strategy to control a WEC array to yield a globally optimal wave energy capture efficiency despite interference effects occurring between the WEC units.

There are two key contributions of this paper: First, the WEC array is depicted as an interconnected system model, i.e., a network model with dynamic edges, in which the node dynamics are the dynamics of the WEC units and the edge dynamics are the dynamics of the interfering waves formed between neighboring units. Second, underpinned by an interconnected system model, an MPC method is developed to control each WEC unit in a WEC array to maximize wave energy harvesting.

After reviewing related research in Section 1, an interconnected system model for WEC arrays is developed in Section 2. An MPC strategy for a generalized network model is presented in Section 3. In Section 4, simulations and experiments conducted are described for validating the proposed modeling and MPC methods, respectively. The results are evaluated and discussed in Section 4.2. The conclusions are presented in Section 5.

2. Model of WEC Array

2.1. Model of WEC Unit

A point-absorber-type WEC is primarily composed of a buoy and a PTO unit, as shown in Figure 1a. The WEC array is represented as an interconnected system model in which the topology structure describes the interfering relationship between all WEC units, as shown in Figure 1b. The symbol ∂ in Figure 1b is the angle of the wave propagation direction, and ${\ell }_1,\cdots ,{\ell }_6$ denote the directed interfering relationship between each pair of the three units, Unit 1, Unit 2 and Unit 3. Buoys are used to harvest the kinetic energy of ocean waves, and PTO converts the kinetic energy of the buoy into electric energy. Meanwhile, PTO produces electric power and a damping force, which affects the motion of the buoy simultaneously excited by ocean waves.

The dynamics of a point-absorber WEC are expressed as

$$m_f\ddot{z}\left(t\right)=F_w\left(t\right)+F_m\left(t\right),$$

(1)

where $m_f$ is the mass of the buoy, $z\left(t\right)$ is the vertical displacement of the buoy, $F_w\left(t\right)$ is the hydrodynamic force, and $F_m\left(t\right)$ is the mechanical force.

The variables $F_w\left(t\right)$ and $F_m\left(t\right)$ satisfy

$$\left\{\begin{array}{c}{F}_w\left(t\right)=F_h\left(t\right)+F_r\left(t\right)+F_e\left(t\right),\hfill \\ F_m\left(t\right)=F_o\left(t\right)+F_u\left(t\right),\hfill \end{array}\right.$$

(2)

where $F_h\left(t\right)$, $F_r\left(t\right)$, $F_e\left(t\right)$, $F_o\left(t\right)$ and $F_u\left(t\right)$ are the buoyancy, radiation force, excitation force, friction force of the mechanical device, and reaction force, respectively. The buoyancy $F_h\left(t\right)$ is

$$F_h\left(t\right)=S_h\left[\mathrm{w}\left(t\right)-z\left(t\right)\right],$$

(3)

where $\mathrm{w}\left(t\right)$ and $S_h$ are the vertical height of the incident wave and constant buoyancy coefficient, respectively. The radiation force $F_r\left(t\right)$ is

$$F_r\left(t\right)=F_D\left(t\right)-m_{\infty }\ddot{z}\left(t\right),$$

(4)

where $m_{\infty }$ denotes the additional mass related to the quantity of the water being moved and is used as the coefficient of the linear term of the radiation force in the case of infinite frequency. Inspired by [20], the convolution term $F_D\left(t\right)$ is approximated as $D·\left(\dot{\mathrm{w}}\left(t\right)-\dot{z}\left(t\right)\right)$, where D is a constant damping coefficient to roughly approximate the Fourier transform $\widehat{D}\left(jw\right)$ of $F_D\left(t\right)$. $\widehat{D}\left(jw\right)$ can be regarded as a frequency-dependent damping coefficient.

The excitation force $F_e\left(t\right)$ is a non-causal force that is approximated by a fifth-order state-space representation of the relationship between the wave height and excitation force [21]. Since the influence of $F_o\left(t\right)$ is relatively slight, we assume $F_o\left(t\right)=0$ [22]. The proposed control algorithm calculates the reaction force $F_u\left(t\right)$.

Let $m=m_f+m_{\infty }$ and $\xi \left(t\right)={\left[F_h\left(t\right),\dot{z}\left(t\right)\right]}^T$. Then, the WEC unit model is achieved as

$$\begin{array}{cc}\hfill & \dot{\xi }\left(t\right)=A\xi \left(t\right)+B_w\dot{\mathrm{w}}\left(t\right)+B\left[F_u\left(t\right)+F_e\left(t\right)\right],\hfill \\ \hfill & y\left(t\right)=C\xi \left(t\right),\hfill \\ \hfill & y_z\left(t\right)=C_z\xi \left(t\right),\hfill \end{array}$$

(5)

where

$$\begin{array}{c}A=\left[\begin{array}{cc}0& -S_h\\ -\frac{1}{m}& -\frac{D}{m}\end{array}\right],B_w=\left[\begin{array}{c}{S}_h\\ \frac{D}{m}\end{array}\right],B=\left[\begin{array}{c}0\\ \frac{1}{m}\end{array}\right],C=\left[\begin{array}{cc}0\hfill & 1\hfill \end{array}\right],\mathrm{and}{C}_z=\left[\begin{array}{cc}1\hfill & 0\hfill \end{array}\right].\end{array}$$

Remark 1.

According to (5), it holds that $y\left(t\right)=\dot{z}\left(t\right)$, and $y_z\left(t\right)=F_h\left(t\right)$.

2.2. WEC Array

2.2.1. WEC Units in Array

When the WEC units are close together, their interactions are enhanced, thus affecting the movement of the buoys more significantly. Radiation waves are the main source of interactions and therefore cannot be ignored. Considering the interaction effect, the i-th ($i=1,\dots ,N$) WEC unit receives the cumulative effect of the radiation wave velocity ${\dot{\mathrm{w}}}_{ij}$ generated by the other j-th ($j=1,\dots ,N-1,j\ne i$) WEC units, in addition to the effect of the incident wave velocity ${\dot{\mathrm{w}}}_i\left(t\right)$ [23]. The model of the i-th WEC unit is formulated as follows,

$$\begin{array}{cc}\hfill & {\dot{\xi }}_i\left(t\right)=A_i{\xi }_i\left(t\right)+B_{wi}\left[{\dot{\mathrm{w}}}_i\left(t\right)+\sum _{j\ne i}^{N-1}{\dot{\mathrm{w}}}_{ij}\left(t\right)\right]+B_i\left(F_{ui}\left(t\right)+F_{ei}\left(t\right)\right),\hfill \\ \hfill & y_i\left(t\right)=C_i{\xi }_i\left(t\right),\hfill \\ \hfill & y_{zi}\left(t\right)=C_{zi}{\xi }_i\left(t\right),\hfill \end{array}$$

(6)

where ${\xi }_i\left(t\right)={[F_{hi}\left(t\right),{\dot{z}}_i\left(t\right)]}^T$ is the state of the i-th WEC unit, $F_{hi}\left(t\right)$ is the buoyancy of the i-th unit, ${\dot{z}}_i\left(t\right)$ is the velocity of the i-th buoy. In this model, the input to the WEC unit is the total of the velocity of the incident wave and the radiated waves generated by the other WEC units.

Moreover, the array composed of N WEC units is considered as a continuous interconnected system modeled as follows:

$$\begin{array}{cc}\hfill & \dot{\tilde{\xi }}\left(t\right)={\mathbb{A}}_c\tilde{\xi }\left(t\right)+{\mathbb{B}}_{wc}\dot{\tilde{\mathrm{w}}}\left(t\right)+{\mathbb{B}}_{wc}{\dot{\tilde{\mathrm{w}}}}_{\mathrm{ij}}\left(t\right)+{\mathbb{B}}_c\left[{\tilde{F}}_u\left(t\right)+{\tilde{F}}_e\left(t\right)\right],\hfill \\ \hfill & \tilde{y}\left(t\right)={\mathbb{C}}_c\tilde{\xi }\left(t\right),\hfill \\ \hfill & {\tilde{y}}_z\left(t\right)={\mathbb{C}}_{zc}\tilde{\xi }\left(t\right),\hfill \end{array}$$

(7)

where

$$\begin{array}{c}\tilde{\xi }\left(t\right)=\left[\begin{array}{c}{\xi }_1\\ {\xi }_2\\ ⋮\\ {\xi }_N\end{array}\right],\dot{\tilde{\mathrm{w}}}\left(t\right)=\left[\begin{array}{c}{\dot{\mathrm{w}}}_1\\ {\dot{\mathrm{w}}}_2\\ ⋮\\ {\dot{\mathrm{w}}}_N\end{array}\right],{\dot{\tilde{\mathrm{w}}}}_{\mathrm{ij}}=\left[\begin{array}{c}\sum _{j\ne 1}^{N-1}{\dot{\mathrm{w}}}_{1j}\\ \sum _{j\ne 2}^{N-1}{\dot{\mathrm{w}}}_{2j}\\ ⋮\\ \sum _{j\ne N}^{N-1}{\dot{\mathrm{w}}}_{Nj}\end{array}\right],{\tilde{F}}_u\left(t\right)=\left[\begin{array}{c}{F}_{u1}\\ F_{u2}\\ ⋮\\ F_{uN}\end{array}\right],\\ {\tilde{F}}_e\left(t\right)=\left[\begin{array}{c}{F}_{e1}\\ F_{e2}\\ ⋮\\ F_{eN}\end{array}\right],\tilde{y}\left(t\right)=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_N\end{array}\right],{\tilde{y}}_z\left(t\right)=\left[\begin{array}{c}{y}_{z1}\\ y_{z2}\\ ⋮\\ y_{zN}\end{array}\right],\\ {\mathbb{A}}_c={\mathbb{I}}_N\otimes A,{\mathbb{B}}_{wc}={\mathbb{I}}_N\otimes B_w,{\mathbb{B}}_c={\mathbb{I}}_N\otimes B,{\mathbb{C}}_c={\mathbb{I}}_N\otimes C,{\mathbb{C}}_{zc}={\mathbb{I}}_N\otimes C_z,\end{array}$$

and the subscript c indicates the term of continuous system. $\mathbb{I}$ is an identity matrix.

Remark 2.

In this paper, a system connecting several independent WEC units with a particular network structure is called an interconnected system, where the WEC units are coupled with each other after being connected by the network. The primary distinction between the network of an interconnected system and a typical network is that the former processes dynamical edges, which can obtain system inputs based on node outputs and transmit outputs to the inputs of other nodes [24].

2.2.2. Interactions among WEC Units

According to Equation (6), ${\dot{\mathrm{w}}}_{ij}\left(t\right)$ indicates a variation in the vertical displacement of WEC unit i excited by the interference wave generated by WEC unit j. ${\dot{\mathrm{w}}}_{ij}\left(t\right)$ is associated to the locations of WEC unit i and j, and the heave motion ${\dot{z}}_j\left(t\right)$ of WEC unit j. Using experimental data, and ignoring the non-linearity relationship between ${\dot{\mathrm{w}}}_{ij}\left(t\right)$ and ${\dot{z}}_j\left(t\right)$, this transfer function ${\mathrm{H}}_{ij}\left(s\right)$ of a linear dynamic system with ${\dot{z}}_j\left(t\right)$ and ${\dot{\mathrm{w}}}_{ij}\left(t\right)$ as input and output, respectively, can be derived by standard linear system identification techniques [23]. In addition to being an edge with dynamic properties in the network topology of the interconnected system, this linear dynamical model can be expressed simultaneously as a frequency characteristics model ${\mathrm{H}}_{ij}\left(\mathrm{j}\omega \right)$ (simplified as ${\mathrm{H}}_{ij}$) as follows:

$${\dot{w}}_{ij}\left(\omega \right)={\mathrm{H}}_{ij}{\dot{z}}_j\left(\omega \right),$$

(8)

$${\mathrm{H}}_{ij}=\mathrm{j}\omega \sqrt{\frac{{\pi \left|\omega \right|}^4}{8g^3}}{a}^2\frac{1}{\sqrt{d_{ij}}}{e}^{-\mathrm{j}\left\{\frac{{\omega }^2}{g}{d}_{ij}+\frac{3\pi }{4}\right\}},$$

(9)

where ${\dot{z}}_j\left(\omega \right)$ and ${\dot{w}}_{ij}\left(\omega \right)$ are defined as the Fourier transform of ${\dot{z}}_j\left(t\right)$ and ${\dot{\mathrm{w}}}_{ij}\left(t\right)$; g, a and $d_{ij}$ are the gravity constant, the diameter of the buoy, and the distance between the i-th WEC unit and the j-th WEC unit, respectively; $\omega$ is the frequency and j is the imaginary number. For $\omega >0$, the frequency characteristics model Equation (8) is equivalently formulated as [25]: 

$${\mathrm{H}}_{ij}=\mathrm{A}\left(\omega \right)e^{-\mathrm{j}\varphi \left(\omega \right)}$$

(10)

where $\mathrm{A}\left(\omega \right)={\omega }^3{a}^2\sqrt{\frac{\pi }{8g^3{d}_{ij}}}\mathrm{and}\varphi \left(\omega \right)=\frac{d_{ij}{\omega }^2}{g}+\frac{\pi }{4}$. The frequency characteristics model Equation (10) as a linear system is rewritten as a state-space formula as follows:

$$\begin{array}{cc}\hfill {\dot{\zeta }}_{\ell }\left(t\right)& ={\mathcal{E}}_{\ell }{\zeta }_{\ell }\left(t\right)+{\mathcal{F}}_{\ell }{u}_{\ell }\left(t\right),\hfill \\ \hfill v_{\ell }\left(t\right)& ={\mathcal{H}}_{\ell }{\zeta }_{\ell }\left(t\right),\hfill \end{array}$$

(11)

where ${\mathcal{E}}_{\ell }\in {\mathbb{R}}^{\left(5\times 5\right)}$, ${\mathcal{F}}_{\ell }\in {\mathbb{R}}^{\left(5\times 1\right)}$, ${\mathcal{H}}_{\ell }\in {\mathbb{R}}^{\left(1\times 5\right)}$, ${\zeta }_{\ell }\left(t\right)\in {\mathbb{R}}^{\left(5\times 1\right)}$ are the system matrix, input matrix, output matrix and state, respectively; $u_{\ell }\left(t\right)$ and $v_{\ell }\left(t\right)$ are the input and output of the ℓ-th interaction state-space model, respectively. Equation (11) denotes the ℓ-th ($\ell =1,\dots ,2M,M=\frac{N(N-1)}{2}$) directed interaction effect between any two WEC units. This interaction effect is the weight of the corresponding edge of the interconnected system, as shown in Figure 1b.

In the WEC array, the integrated interaction model is expressed as

$$\begin{array}{c}{\dot{\tilde{\zeta }}}_{\ell }\left(t\right)={\mathbb{E}}_c{\tilde{\zeta }}_{\ell }\left(t\right)+{\mathbb{F}}_c{\tilde{u}}_{\ell }\left(t\right),\hfill \\ {\tilde{v}}_{\ell }\left(t\right)={\mathbb{H}}_c{\tilde{\zeta }}_{\ell }\left(t\right),\hfill \end{array}$$

(12)

where ${\mathbb{E}}_c={\mathbb{I}}_{2M}\otimes {\mathcal{E}}_{\ell }$, ${\mathbb{F}}_c={\mathbb{I}}_{2M}\otimes {\mathcal{F}}_{\ell }$, ${\mathbb{H}}_c={\mathbb{I}}_{2M}\otimes {\mathcal{H}}_{\ell }$, and

$$\begin{array}{c}{\tilde{\zeta }}_{\ell }\left(t\right)=\left[\begin{array}{c}{\zeta }_1\\ {\zeta }_2\\ ⋮\\ {\zeta }_{2M}\end{array}\right],{\tilde{u}}_{\ell }\left(t\right)=\left[\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_{2M}\end{array}\right],{\tilde{v}}_{\ell }\left(t\right)=\left[\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ v_{2M}\end{array}\right].\end{array}$$

2.2.3. WEC Array System

Since a WEC array is an interconnected system, it can be represented as a graph $\mathcal{G}$. The incidence matrix $\mathcal{D}$ of the graph $\mathcal{G}$ is applied to denote the interconnected system with N nodes and $2M$ directed edges. The element located at the i-th row and ℓ-th column of $\mathcal{D}$ is defined as

$${\mathcal{D}}_{i\ell }=\left\{\begin{array}{cc}+1,\hfill & \mathrm{node}i\mathrm{is}\mathrm{the}\mathrm{head}\mathrm{of}\mathrm{edge}\ell ,\hfill \\ -1,\hfill & \mathrm{node}i\mathrm{is}\mathrm{the}\mathrm{tail}\mathrm{of}\mathrm{edge}\ell ,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(13)

The incidence matrix $\mathcal{D}$ describes the connection relationship between the edges and nodes in graph $\mathcal{G}$, and is used to determine the edge on which the output of the node is located and the node at which the input of the edge arrives. From the perspective of the node inputs and outputs, the incidence matrix satisfies

$$\mathcal{D}={\mathbb{D}}_{\mathrm{O}}-{\mathbb{D}}_{\mathrm{I}}^{\mathrm{T}},$$

(14)

where ${\mathbb{D}}_{\mathrm{O}}$ and ${\mathbb{D}}_{\mathrm{I}}$ are the input and output matrices of the nodes, respectively, as

$$\begin{array}{c}{\mathbb{D}}_{\mathrm{O}i\ell }=\left\{\begin{array}{cc}1,\hfill & \mathrm{edge}\ell \mathrm{is}\mathrm{the}\mathrm{input}\mathrm{of}\mathrm{node}i,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.,{\mathbb{D}}_{\mathrm{I}\ell i}=\left\{\begin{array}{cc}1,\hfill & \mathrm{edge}\ell \mathrm{is}\mathrm{the}\mathrm{output}\mathrm{of}\mathrm{node}i,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

Furthermore, considering Equations (6), (11) and (13), ${\sum }_{j\ne i}^{N-1}{\dot{\mathrm{w}}}_{ij}\left(t\right)$ in Equation (6) and $u_{\ell }\left(t\right)$ in Equation (11) are denoted as

$$\left\{\begin{array}{cc}\hfill \sum _{j\ne i}^{N-1}{\dot{\mathrm{w}}}_{ij}\left(t\right)& =\sum _{\ell =1}^M{\mathbb{D}}_{\mathrm{O}i\ell }{v}_{\ell }\left(t\right),i=1,2,\dots ,N,\hfill \\ \hfill u_{\ell }\left(t\right)& =\sum _{i=1}^N{\mathbb{D}}_{\mathrm{I}\ell i}{y}_i\left(t\right),\ell =1,2,\dots ,2M.\hfill \end{array}\right.$$

(15)

Moreover, Equation (15) can be written in matrix form as follows

$$\left\{\begin{array}{c}{\dot{\tilde{\mathrm{w}}}}_{\mathrm{ij}}\left(t\right)={\mathbb{D}}_{\mathrm{O}}{\tilde{v}}_{\ell }\left(t\right),\hfill \\ {\tilde{u}}_{\ell }\left(t\right)={\mathbb{D}}_{\mathrm{I}}\tilde{y}\left(t\right).\hfill \end{array}\right.$$

(16)

For example, for the three-unit WEC array interconnected system shown in Figure 1b, the incidence matrix $\mathcal{D}$ is

$$\mathcal{D}=\left[\begin{array}{cccccc}1& 0& -1& -1& 0& 1\\ -1& 1& 0& 1& -1& 0\\ 0& -1& 1& 0& 1& -1\end{array}\right].$$

(17)

Considering Equation (14), the incidence matrix $\mathcal{D}$ is decomposed into ${\mathbb{D}}_{\mathrm{O}}$ and ${\mathbb{D}}_{\mathrm{I}}$. According to Equations (15) and (16), we can obtain

$${\dot{\tilde{\mathrm{w}}}}_{\mathrm{ij}}\left(t\right)=\left[\begin{array}{c}\sum _{j\ne 1}^{N-1}{\dot{\mathrm{w}}}_{1j}\left(t\right)\\ \sum _{j\ne 2}^{N-1}{\dot{\mathrm{w}}}_{2j}\left(t\right)\\ \sum _{j\ne 3}^{N-1}{\dot{\mathrm{w}}}_{3j}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{v}_1+v_6\\ v_2+v_4\\ v_3+v_5\end{array}\right],{\tilde{u}}_{\ell }\left(t\right)=\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right]=\left[\begin{array}{c}{y}_2\\ y_3\\ y_1\\ y_1\\ y_2\\ y_3\end{array}\right].$$

(18)

By substituting Equation (16) into Equations (7) and (12), the model of a three-unit WEC array as an interconnected system is obtained. More generally, a formally uniform state-space model is achieved to express a WEC array system with N units as

$$\begin{array}{c}\left[\begin{array}{c}\dot{\tilde{\xi }}\left(t\right)\\ {\dot{\tilde{\zeta }}}_{\ell }\left(t\right)\end{array}\right]=\left[\begin{array}{cc}{\mathbb{A}}_c& {\mathbb{B}}_{wc}{\mathbb{D}}_{\mathrm{O}}{\mathbb{H}}_c\\ {\mathbb{F}}_c{\mathbb{D}}_{\mathrm{I}}{\mathbb{C}}_c& {\mathbb{E}}_c\end{array}\right]\left[\begin{array}{c}\tilde{\xi }\left(t\right)\\ {\tilde{\zeta }}_{\ell }\left(t\right)\end{array}\right]+\left[\begin{array}{c}{\mathbb{B}}_{wc}\\ \mathbf{0}\end{array}\right]\dot{\tilde{\mathrm{w}}}\left(t\right)+\left[\begin{array}{c}{\mathbb{B}}_c\\ \mathbf{0}\end{array}\right]\left({\tilde{F}}_u\left(t\right)+{\tilde{F}}_e\left(t\right)\right),\\ \tilde{y}\left(t\right)=\left[{\mathbb{C}}_c,\mathbf{0}\right]\left[\begin{array}{c}\tilde{\xi }\left(t\right)\\ {\tilde{\zeta }}_{\ell }\left(t\right)\end{array}\right],\\ {\tilde{y}}_z\left(t\right)=\left[{\mathbb{C}}_{zc},\mathbf{0}\right]\left[\begin{array}{c}\tilde{\xi }\left(t\right)\\ {\tilde{\zeta }}_{\ell }\left(t\right)\end{array}\right],\end{array}$$

(19)

where $\dot{\tilde{w}}\left(t\right)=\dot{\tilde{\mathrm{w}}}\left(t\right)+{\dot{\tilde{\mathrm{w}}}}_{\mathrm{ij}}\left(t\right)$.

3. MPC Algorithm Design

3.1. Control Problem Statement

Equation (19) is discretized by a zero-order hold method as

$$\begin{array}{c}\left[\begin{array}{c}\dot{\tilde{\xi }}(k+1)\\ {\tilde{\zeta }}_{\ell }(k+1)\end{array}\right]=\left[\begin{array}{cc}\mathbb{A}& {\mathbb{B}}_w{\mathbb{D}}_{\mathrm{O}}\mathbb{H}\\ \mathbb{F}{\mathbb{D}}_{\mathrm{I}}\mathbb{C}& \mathbb{E}\end{array}\right]\left[\begin{array}{c}\tilde{\xi }\left(k\right)\\ {\tilde{\zeta }}_{\ell }\left(k\right)\end{array}\right]+\left[\begin{array}{c}{\mathbb{B}}_w\\ \mathbf{0}\end{array}\right]\dot{\tilde{\mathrm{w}}}\left(k\right)+\left[\begin{array}{c}\mathbb{B}\\ \mathbf{0}\end{array}\right]\left({\tilde{F}}_u\left(k\right)+{\tilde{F}}_e\left(k\right)\right),\\ \tilde{y}\left(k\right)=[\mathbb{C},\mathbf{0}]\left[\begin{array}{c}\tilde{\xi }\left(k\right)\\ {\tilde{\zeta }}_{\ell }\left(k\right)\end{array}\right],\\ {\tilde{y}}_z\left(k\right)=\left[{\mathbb{C}}_z,\mathbf{0}\right]\left[\begin{array}{c}\tilde{\xi }\left(k\right)\\ {\tilde{\zeta }}_{\ell }\left(k\right)\end{array}\right],\end{array}$$

(20)

where

$$\begin{array}{cc}\hfill & \mathbb{A}=e^{{\mathbb{A}}_{c}T},\mathbb{B}={\int }_0^T{e}^{{\mathbb{A}}_c\tau }d\tau {\mathbb{B}}_c,{\mathbb{B}}_w={\int }_0^T{e}^{{\mathbb{A}}_c\tau }d\tau {\mathbb{B}}_{wc},\mathbb{C}={\mathbb{C}}_c,{\mathbb{C}}_z={\mathbb{C}}_{zc},\mathbb{E}=e^{{\mathbb{E}}_{c}T},\hfill \\ \hfill & \mathbb{F}={\int }_0^T{e}^{{\mathbb{E}}_c\tau }d\tau {\mathbb{F}}_c,\mathbb{H}={\mathbb{H}}_c,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \dot{\tilde{\mathrm{w}}}\left(k\right)=\left[\begin{array}{c}{\dot{\mathrm{w}}}_1\left(k\right)\\ {\dot{\mathrm{w}}}_2\left(k\right)\\ ⋮\\ {\dot{\mathrm{w}}}_N\left(k\right)\end{array}\right],\widetilde{F_u}\left(k\right)=\left[\begin{array}{c}{F}_{u1}\left(k\right)\\ F_{u2}\left(k\right)\\ ⋮\\ F_{uN}\left(k\right)\end{array}\right],\widetilde{F_e}\left(k\right)=\left[\begin{array}{c}{F}_{e1}\left(k\right)\\ F_{e2}\left(k\right)\\ ⋮\\ F_{eN}\left(k\right)\end{array}\right],\tilde{y}\left(k\right)=\left[\begin{array}{c}{y}_1\left(k\right)\\ y_2\left(k\right)\\ ⋮\\ y_N\left(k\right)\end{array}\right],{\tilde{y}}_z\left(k\right)=\left[\begin{array}{c}{y}_{z1}\left(k\right)\\ y_{z2}\left(k\right)\\ ⋮\\ y_{zN}\left(k\right)\end{array}\right].\hfill \end{array}$$

In Equation (20), $\tilde{\xi }\left(k\right)\in {\mathbb{R}}^{\left(\alpha \times N\right)\times 1}$ is the states set of all the WEC units at instant k, where $\alpha$ is the model order of each WEC unit. In addition, $\widetilde{{\zeta }_{\ell }}\left(k\right)\in {\mathbb{R}}^{\left(\beta \times 2M\right)\times 1}$ is the states set of all interactions at time instant k, where $\beta$ is the model order of each interaction. The wave excitation force and wave velocity are assumed to be predictable in this paper.

The MPC method controls the WEC array to maximize the harvesting of wave energy and convert it into electrical energy while satisfying the constraints. Since the buoy of the point-absorbing WEC is directly connected to PTO, the inner product of the buoy velocity and the corresponding reaction force is equal to the total power of the wave energy captured by the WEC array. As the velocity of the buoy is opposite to the positive direction of the reaction force, it holds that

$$P\left(k\right)=-{\tilde{F}}_u^{\mathrm{T}}\left(k\right)\tilde{y}\left(k\right),$$

(21)

and

$$\underset{{\tilde{F}}_{uk}^{k+N_p}}{max}\sum _{T=0}^{N_p}P\left(k+T\right)=\underset{{\tilde{F}}_{uk}^{k+N_p}}{min}\sum _{T=0}^{N_p}{\widetilde{F_u}}^{\mathrm{T}}\left(k+T\right)\tilde{y}\left(k+T\right),$$

(22)

where $N_p$ is the predicted time domain, and T is the time step. To ensure the convexity of the optimization problem, Equation (22) is modified as the cost function

$$J=\sum _{T=0}^{N_p}\left[{\tilde{F}}_u^{\mathrm{T}}(k+T)\tilde{y}(k+T)+r{\tilde{F}}_u^{\mathrm{T}}(k+T){\tilde{F}}_u(k+T)+q{\tilde{y}}_z^{\mathrm{T}}(k+T){\tilde{y}}_z(k+T)\right],$$

(23)

where r and q are two artificially determined weighting factors. The first part of the cost function Equation (23) is the wave energy harvesting power in Equation (22). The second part of Equation (23) represents a soft constraint on the control force, primarily to guarantee that the solution of the MPC is to solve a convex optimization problem and the control force is within the PTO’s operating range. The third part represents the penalty term of the relative displacement of the buoy and the wave surface, ensuring the PTO operates within an effective operating range. To guarantee the safe operation of the WEC array and reduce the loss and maintenance cost, it is necessary to restrict the amplitude of the reaction force. Moreover, the buoy needs to be constrained to keep the proportional displacement between the midpoint of the buoy and the wave surface within a safe range. These constraints are formulated as follows,

$$\begin{array}{cc}\hfill & \left|F_u\left(k+T\right)\right|\le F_{u,max},\mathrm{for}T=0,1,\dots ,N_p,\hfill \\ \hfill & \left|y_z\left(k+T\right)\right|\le y_{z,max},\mathrm{for}T=0,1,\dots ,N_p,\hfill \end{array}$$

(24)

where $F_{u,max}$ and $y_{z,max}$ are the maximum amplitude of the control force and buoyancy, respectively.

The MPC aims to solve the optimal control input sequence that minimizes the cost function Equation (23) in the prediction time domain $N_p$, under the premise that the constraint Equation (24) is satisfied.

Remark 3.

The weighting factors r and q in Equation (23) are specified as two specific values independently to ensure that the MPC is a convex optimization problem.

3.2. Prediction Model

The cost function Equation (23) can be expressed in matrix form as

$$J={\left({\tilde{U}}_k^{k+N_p}\right)}^{\mathrm{T}}{\tilde{Y}}_k^{k+N_p}+{\left({\tilde{U}}_k^{k+N_p}\right)}^{\mathrm{T}}R{\tilde{U}}_k^{k+N_p}+{\left({\tilde{Y}}_{zk}^{k+N_p}\right)}^{\mathrm{T}}Q{\tilde{Y}}_{zk}^{k+N_p},$$

(25)

where $R=r\times I_{N\times \left(N_P+1\right)},Q=q\times I_{N\times \left(N_P+1\right)}$,

$${\tilde{U}}_k^{k+N_p}=\left[\begin{array}{c}{\tilde{F}}_u\left(k\right)\\ {\tilde{F}}_u(k+1)\\ ⋮\\ {\tilde{F}}_u\left(k+N_p\right)\end{array}\right],{\tilde{Y}}_k^{k+N_p}=\left[\begin{array}{c}\tilde{y}\left(k\right)\\ \tilde{y}(k+1)\\ ⋮\\ \tilde{y}\left(k+N_p\right)\end{array}\right],{\tilde{Y}}_{zk}^{k+N_p}=\left[\begin{array}{c}{\tilde{y}}_z\left(k\right)\\ {\tilde{y}}_z(k+1)\\ ⋮\\ {\tilde{y}}_z\left(k+N_p\right)\end{array}\right],$$

and I is a unit diagonal matrix.

In the cost function, it is necessary to know the buoy velocity and buoyancy of each unit in the predicted time domain. By the discrete model of the interconnected system, we derive a prediction model from k to $k+N_p$ as follows:

$$\begin{array}{c}{\tilde{Y}}_k^{k+N_p}=\Lambda X\left(k\right)+{\Phi }_u\left({\tilde{U}}_k^{k+N_p}+{{\tilde{F}}_{ek}}^{k+N_p}\right)+{\Phi }_w{\dot{\tilde{W}}}_k^{k+N_p},\hfill \\ {\tilde{Y}}_{zk}^{k+N_p}={\Lambda }_{z}X\left(k\right)+{\Phi }_{uz}\left({\tilde{U}}_k^{k+N_p}+{{\tilde{F}}_{ek}}^{k+N_p}\right)+{\Phi }_{wz}{\dot{\tilde{W}}}_k^{k+N_p},\hfill \end{array}$$

(26)

where

$$X\left(k\right)=\left[\begin{array}{c}\tilde{\xi }\left(k\right)\\ {\tilde{\zeta }}_{\ell }\left(k\right)\end{array}\right],{\tilde{F}}_{ek}^{k+N_p}=\left[\begin{array}{c}{\tilde{F}}_e\left(k\right)\\ {\tilde{F}}_e(k+1)\\ ⋮\\ {\tilde{F}}_e\left(k+N_p\right)\end{array}\right],{\dot{\tilde{W}}}_k^{k+N_p}=\left[\begin{array}{c}\dot{\tilde{w}}\left(k\right)\\ \dot{\tilde{w}}(k+1)\\ ⋮\\ \dot{\tilde{w}}\left(k+N_p\right)\end{array}\right],\Lambda =\left[\begin{array}{c}\mathcal{C}\\ \mathcal{C}\mathcal{A}\\ \mathcal{C}{\mathcal{A}}^2\\ ⋮\\ \mathcal{C}{\mathcal{A}}^{N_p}\end{array}\right],$$

$${\Phi }_u=\left[\begin{array}{ccccc}\mathbf{0}& & & & \\ \mathcal{C}{\mathcal{B}}_u& \mathbf{0}& & & \\ \mathcal{C}\mathcal{A}{\mathcal{B}}_u& \mathcal{C}{\mathcal{B}}_u& \mathbf{0}& & \\ ⋮& ⋮& \ddots & \ddots & \\ \mathcal{C}{\mathcal{A}}^{N_p-1}{\mathcal{B}}_u& \mathcal{C}{\mathcal{A}}^{N_p-2}{\mathcal{B}}_u& \cdots & \mathcal{C}{\mathcal{B}}_u& \mathbf{0}\end{array}\right],$$

$${\Phi }_w=\left[\begin{array}{ccccc}\mathbf{0}& & & & \\ \mathcal{C}{\mathcal{B}}_w& \mathbf{0}& & & \\ \mathcal{C}\mathcal{A}{\mathcal{B}}_w& \mathcal{C}{\mathcal{B}}_w& \mathbf{0}& & \\ ⋮& ⋮& \ddots & \ddots & \\ \mathcal{C}{\mathcal{A}}^{Np-1}{\mathcal{B}}_w& \mathcal{C}{\mathcal{A}}^{N_p-2}{\mathcal{B}}_w& \cdots & \mathcal{C}{\mathcal{B}}_w& \mathbf{0}\end{array}\right],$$

$$\Lambda \in {\mathbb{R}}^{\left(N\times \left(N_p+1\right)\right)\times (\alpha N+2\beta M)},{\Phi }_u\in {\mathbb{R}}^{\left(N\times \left(N_p+1\right)\right)\times \left(N\times \left(N_p+1\right)\right)},$$

$${\Phi }_w\in {\mathbb{R}}^{\left(N\times \left(N_p+1\right)\right)\times \left(N\times \left(N_p+1\right)\right)},$$

with

$$\mathcal{A}=\left[\begin{array}{cc}\mathbb{A}& {\mathbb{B}}_w{\mathbb{D}}_{i\ell }\mathbb{H}\\ \mathbb{F}{\mathbb{D}}_{\ell i}\mathbb{C}& \mathbb{E}\end{array}\right],{\mathcal{B}}_u=\left[\begin{array}{c}\mathbb{B}\\ \mathbf{0}\end{array}\right],{\mathcal{B}}_w=\left[\begin{array}{c}{\mathbb{B}}_w\\ \mathbf{0}\end{array}\right],\mathcal{C}=[\mathbb{C},\mathbf{0}],{\mathcal{C}}_z=[{\mathbb{C}}_z,\mathbf{0}].$$

Similarly, ${\Lambda }_z$, ${\Phi }_{uz}$ and ${\Phi }_{wz}$ need to replace $\mathcal{C}$ in the corresponding matrix element by ${\mathcal{C}}_z$. Using prediction models Equation (24) to represent their constraints in the predicted time domain achieves

$$\begin{array}{ll}\hfill & \left|{\tilde{U}}_k^{k+N_p}\right|\le U_{max},\\ \hfill & \left|{\Lambda }_{z}x\left(k\right)+{\Phi }_{uz}\left({\tilde{U}}_k^{k+N_p}+{{\tilde{F}}_{ek}}^{k+N_p}\right)+{\Phi }_{wz}{\dot{\tilde{W}}}_k^{k+N_p}\right|\le Y_{z,max},\end{array}$$

(27)

with

$$U_{max}=\left[\begin{array}{c}{F}_{u,max}\\ F_{u,max}\\ ⋮\\ F_{u,max}\end{array}\right]\in \mathbb{R}\left(N\times N_p\right)\times 1,Y_{z,max}=\left[\begin{array}{c}{y}_{z,max}\\ y_{z,max}\\ ⋮\\ y_{z,max}\end{array}\right]\in \mathbb{R}\left(N\times N_p\right)\times 1.$$

3.3. Quadratic Programming Problem

We need to solve the optimal reaction force sequence so that the cost function J can obtain the minimum value as

$${\tilde{U}}_k^{*k+N_p}=arg\underset{{\tilde{U}}_k^{k+N_p}}{\min }J.$$

(28)

Equation (25) is reorganized into a general form to utilize the standard QP solver.

Substituting the prediction model Equation (26) into the cost function Equation (25), and ignoring terms unrelated to ${\tilde{U}}_k^{k+N_p}$, leads to

$$J=\frac{1}{2}{\left({\tilde{U}}_k^{k+N_p}\right)}^{\mathrm{T}}\mathcal{H}{\tilde{U}}_k^{k+N_p}+{\mathcal{F}}^{\mathrm{T}}{\tilde{U}}_k^{k+N_p},$$

(29)

where

$$\begin{array}{c}\mathcal{H}={\Phi }_u+{\Phi }_u^{\mathrm{T}}+2R+2{\Phi }_{uz}^{\mathrm{T}}Q{\Phi }_{uz},\\ \mathcal{F}=\left(\Lambda +2{\Phi }_{uz}^{\mathrm{T}}Q{\Lambda }_z\right)X\left(k\right)+\left({\Phi }_w+2{\Phi }_{uz}^{\mathrm{T}}Q{\Phi }_{wz}\right){\dot{\tilde{W}}}_k^{k+N_p}+\left({\Phi }_u+2{\Phi }_{uz}^{\mathrm{T}}Q{\Phi }_{uz}\right){\tilde{F}}_{ek}^{k+N_p}.\end{array}$$

(30)

The constraints Equation (27) are reorganized into

$$\mathcal{S}{\tilde{U}}_k^{k+N_p}\le b_u,$$

(31)

where $\mathcal{S}=\left[\begin{array}{c}I\\ -I\\ {\Phi }_{uz}\\ -{\Phi }_{uz}\end{array}\right]$, and $b_u=\left[\begin{array}{c}{U}_{max}\\ U_{max}\\ Y_{z,max}-{\Lambda }_{xz}x\left(k\right)-{\Phi }_{wz}{\dot{\tilde{W}}}_k^{k+N_p}-{\Phi }_{uz}{\tilde{F}}_{ek}^{k+N_p}\\ Y_{z,max}+{\Lambda }_{xz}x\left(k\right)+{\Phi }_{wz}{\dot{\tilde{W}}}_k^{k+N_p}+{\Phi }_{uz}{\tilde{F}}_{ek}^{k+N_p}\end{array}\right].$

Finally, the general form of the QP problem is summarized as

$$\begin{array}{c}{\tilde{U}}^{*}=arg\underset{U_0^k}{min}\frac{1}{2}{\left({\tilde{U}}_k^{k+N_p}\right)}^{\mathrm{T}}\mathcal{H}{\tilde{U}}_k^{k+N_p}+{\mathcal{F}}^{\mathrm{T}}{\tilde{U}}_k^{k+N_p},\\ \mathrm{S}.\mathrm{t}.\mathcal{S}{\tilde{U}}_k^{k+N_p}\le b_u.\end{array}$$

(32)

In each control horizon, the first term ${\tilde{U}}_k^{k+1}$ in the sequence of optimal control forces ${\tilde{U}}_k^{*k+N_p}$ is adopted as the optimal control force at time k.

4. Experimental Results and Discussion

4.1. Simulation of Wave Energy Harvesting

We designed a semi-physical hardware-in-the-loop (HIL) experiment to simulate the wave energy capture more realistically than the digital simulation, and to reduce the high cost and complexity of the real ocean or large pool environment. The semi-physical experimental platform consisted of a front-end part and a back-end part as shown in Figure 2a. The front-end part comprised a desktop computer running MATLAB${}^{\circledR }$ and an RT-LAB simulator (a rapid control prototyping (RCP) platform from OPAL-RT Technologies Canada), as shown in the upper part of Figure 2a. The front-end part was mainly dedicated to simulating realistic sea conditions to implement a fully digital simulation of the WEC array to validate the established interconnected system model and the MPC control strategy mapped to the back-end part for the physical experiment bench operation.

The back-end part was a physical experiment bench as shown in the lower part of Figure 2a,b. The rotating motor was controlled by a commercial drive (Copley Controls ADP-090-36) and was used to imitate the rippling motion of the WEC buoys in the array under the effect of wave force. As the PTO device of the WEC, the linear generator (see

Table 1. Specifications for wave energy linear generators.

Electrical Parameters	Value	Mechanical Parameters	Value
Rated power (W)	50	Mover stroke (cm)	70
Rated current (A)	20	Mover speed (m/s)	$3.5$ m
Rated voltage (V)	20	Mover length (cm)	18
Conversion efficiency	$85\%$	Mover weight (kg)	$1.65$

for the specific parameters) was dragged by a belt to convert the motion of the simulated buoy into an AC power output. The AC power generated by three PTO devices was converted to DC power using an AC-DC converter and rectifier (modules 6 and 7 in Figure 2b) developed in our laboratory and fed to the DC bus and electronic loads.

The purpose of constructing the above experimental system was to perform the following procedures:

(i)

Verifying the correctness of the interconnected model and the control performance of the proposed MPC algorithm;

(ii)

Comparing the difference in the motion characteristics affected by the interaction and the wave energy harvesting efficiency between the interconnected model and the isolated model of the WEC array, and exploring how to affect the wave energy harvesting efficiency by change in the incident wave propagation direction in the WEC interconnected system.

4.2. Simulation Parameter Setting

We built a three-unit WEC array in the front end with reference to the PB3 device described in [26,27] developed by Ocean Power Technologies Inc. All the parameters are shown in

Table 2. Simulation parameters.

Description	Notation	Value
Diameter	a	$2.65$ m
Height	h	$3.5$ m
Damping	D	$2\times {10}^5$ Nm/s
Buoy mass	m	8300 kg
Added mass	$m_{\infty }$	6100 kg
Stiffness	$S_h$	$2.21\times {10}^5$ N/m
Density of sea water	$\rho$	1020 kg/m${}^3$
Gravity	g	$9.81$ m/s${}^2$
Maximum heave motion	$y_{z,max}$	$1.75$ m
Maximum input magnitude	$F_{u,max}$	$1\times {10}^5$ N

. The three WEC units were arranged in an equilateral triangle with a distance of 25 m from each other.

The JONSWAP function generated the incident wave [28]; the wave significant height was 2 m, and the peak period was 3 s. The wave velocity and the corresponding wave excitation force $F_e$ are shown in Figure 3. To investigate the effect that the angle ∂ of wave propagation direction waves had on the WEC array, three sets of experiments with different wave propagation directions $\partial ={30}^{\circ },{45}^{\circ },{60}^{\circ }$ were set up.

4.3. Results and Discussion

The prediction horizon and the control horizon were both $10\mathrm{s}$, the sampling time was $0.02\mathrm{s}$, and the simulation time was $50\mathrm{s}$. The weight value r was set to $2.5\times {10}^{-6}$ to ensure that the Hermitian matrix was a positive definite matrix, and Q was also set to $2.5\times {10}^{-6}$. When the radius of the buoy was $5\times a\mathrm{m}$ and the distance $d_{ij}$ was $25\mathrm{m}$, the matrices of the state-space representation of the edge model were

$$\begin{array}{c}{\mathcal{E}}_{\ell }=\left[\begin{array}{ccccc}-1.3192& -2.6404& -1.8147& -1.3222& -0.2090\\ 1& & & & \\ & 1& & & \\ & & 1& & \\ & & & 1& \end{array}\right],\hfill \\ {\mathcal{F}}_{\ell }={\left[\begin{array}{ccccc}1& 0& 0& 0& 0\end{array}\right]}^{\mathrm{T}},\mathrm{and}{\mathcal{H}}_{\ell }=\left[\begin{array}{ccccc}0.1470& -0.1931& -0.1043& -0.0100& -0.0053\end{array}\right].\hfill \end{array}$$

Figure 4 shows the velocity of each unit’s radiation wave in the three layout cases. In particular, when $\partial ={30}^{\circ }$, Units 2 and 3 were affected by the same interference. This was because of the symmetrical structure of the equilateral triangle, and the propagation direction of the incident wave overlapping the symmetric axis, so that the waveform of the incident wave received by Units 2 and 3 was synchronized, and the identical edge model caused the same interference. Figure 5 shows that, under the control of the MPC algorithm, the control force of each WEC unit was limited to a range of $\pm 1\times {10}^5$ N, which satisfies the input constraint in Table 1.

In Figure 6, the waveform of each buoy velocity is almost the same as that of the corresponding incident wave velocity, and resonance is almost realized between the buoy motion and the wave motion. Due to the existence of other terms in the cost function, the execution delay, and the limitation of the control force and displacement, the resonance between the buoy and the incident wave is not fully realized. However, the wave energy capture efficiency is optimized by the control of the MPC algorithm.

To show the impact of interactions in the interconnected system, three units of isolated systems were set up as baseline cases. The isolated system had no interaction among the three WEC units. Thus, the radiation waves caused by the oscillating motion of units did not affect the neighborhood units. In the mathematical model Equation (20), the induction matrix $\mathbb{D}$ is an all-0 matrix. Because the interaction exists in the WEC array, the buoy’s motion velocity is different between the interconnected and isolated systems without interaction. Furthermore, to adapt to the existence of the interaction, the MPC algorithm solves the optimal reaction force to make the units coordinate with each other to achieve the global optimal efficiency.

Figure 7 shows the variation in the total power generated by the corresponding three WEC units between the interconnected and the isolated system when the wave incidence angle $\partial ={45}^{\circ }$. Each curve in Figure 7 illustrates the difference in wave energy captured by a unit, taking into account the impact of interfering waves between units compared with those not taking into account interfering waves. The energy error curves for the three units in Figure 7 highlight the various effects of interactions within the array on the energy harvesting of each unit. Figure 8 shows the waveform of the radiation wave and the incident wave when the incident wave angle $\partial ={45}^{\circ }$. Note that in Figure 7, Unit 2 of the interconnected system harvests more wave energy than Unit 2 of the isolated system. In Figure 8, the waveforms of the incident wave received by Unit 2 are nearly the same as those of the radiation wave generated by the other units. The wave energy harvested by Units 1 and 3 of the interconnected system is not as high as the corresponding unit in the isolated system, and the waveforms of the incident wave and the radiation wave are significantly asynchronous.

When the waveform of the radiation wave received by the WEC unit is close to that of the incident wave, the unit can generate more energy. When their phase error is close to zero, so that the radiation wave synchronizes with the incident wave, just as when the buoy motion synchronizes with the wave motion, the efficiency of the wave energy captured by the WEC is optimized. At this time, the radiation wave also enables maximum optimization of the wave energy captured by the corresponding unit of the interconnected system.

To show the influence of the interaction in different cases, we define the index $\varrho$ as

$$\varrho =\frac{P_{\mathrm{array},i}}{{\displaystyle \sum _{j=1}^N}{P}_{\mathrm{isolate},i}^j},$$

(33)

where $P_{\mathrm{array},i}$ and $P_{\mathrm{isolate},i}^j$ are the power of the WEC array as an interconnected system and the power of each isolated unit at $t=i$, respectively.

In the three cases, the interaction of the interconnected system had a positive or negative effect at different instances within 50 s. The interaction in the first and second cases had more positive than negative effects, and $\varrho >1$ was satisfied, respectively, $62.72\%$ and $52.08\%$ of the total time. In the third case, the negative effect was more than the positive effect, and $\varrho >1$ was satisfied only $29.08\%$ of the total time.

The results in Figure 9 show that in the first and second cases, the interconnected system controlled by the MPC algorithm had higher wave energy harvesting efficiency than the isolated system. In the third case, the interconnected system harvested less wave energy than the isolated system. In addition, when the incident waves had different propagation directions, the increase in the wave energy harvesting efficiency of the WEC interconnected system changed accordingly.

5. Conclusions

In this paper, we proposed an interconnected system model of a WEC array which graphically integrates a WEC unit model and an interference model. The proposed model with wave interference dynamics describes the WEC array system more accurately. Each part of the interconnected system model has a well-defined structure that can better reflect uniformity when the array configuration changes. The interconnected model accurately describes the wave energy harvesting process under the influence of wave interference. The MPC algorithm based on the interconnected system model is applied to enable the WEC array to achieve optimal global performance of wave energy capture control, subject to several constraints. Through experimental tests, it was found that the phase error between the radiation and incident waves significantly affected the wave energy capture efficiency, and the wave energy capture efficiency was effectively improved as the error decreased. The experimental results show that the MPC with the interconnected system model had a higher wave energy capture efficiency compared with the isolated system model. However, the performance of the MPC relies on accurate and sufficient information on wave excitation forecasts. Future research will concentrate on more realistic hydrodynamics of the buoys and their coupling in the ocean environment, as well as inaccurate or insufficient wave excitation forecasts. Strengthening the robustness and adaptability of the MPC controller by combining intelligent control techniques is another promising direction.