Model-Predictive Control Strategy for an Array of Wave-Energy Converters

To facilitate the commercialization of wave energy in an array or farm environment, effective control strategies for improving energy extraction efficiency of the system are important. In this paper, we develop and apply model-predictive control (MPC) to a heaving point-absorber array, where the optimization problem is cast into a convex quadratic programming (QP) formulation, which can be efficiently solved by a standard QP solver. We introduced a term for penalizing large slew rates in the cost function to ensure the convexity of this function. Constraints on both range of the states and the input capacity can be accommodated. The convex formulation reduces the computational hurdles imposed on conventional nonlinear MPC. For illustration of the control principles, a point-absorber approximation is adopted to simplify the representation of the hydrodynamic coefficients among the array by exploiting the small devices to wavelength assumption. The energy-capturing capabilities of a two-cylinder array in regular and irregular waves are investigated. The performance of the MPC for this two-WEC array is compared to that for a single WEC, and the behavior of the individual devices in head or beam wave configuration is explained. Also shown is the reactive power required by the power takeoff system to achieve the performance.


Introduction
In the broad subject of wave-energy extraction, various control strategies have been proposed and implemented on single or stand-alone wave-energy converters, including resistive control, approximate complex-conjugate control, and phase control by latching or clutching (see for example, Bacelli et al. (2016) and Hals et al. (2011a)). As commercialization of wave energy will eventually be realized by farms of WECs, the development of control strategies for a WEC array is of high importance for improving the overall efficiency.
As waves scattered around each individual body, hydrodynamics of an array will differ from a group of isolated bodies, which affects power-capture capability and the control strategy that needs to be applied to the array. Evans (1979) and Falnes (1980) independently provided general formulas for the power-capture characteristics of an array of oscillating WECs using frequency-domain analysis, in the context of linear wave theory. The optimal power can be achieved by applying complex-conjugate control to the array. In their works, an interaction factor, named q factor, was introduced to represent the ratio of the power captured by an array to that captured by the sum from each member element acting in isolation. This approach was then applied in various studies of power absorption of different WEC arrays (Thomas and Evans 1981;Evans 1981;Falnes and Budal 1982;McIver 1994) with or without motion constraints on the devices. However, since the focus was on just the theoretical limit of the power generated by a given array, not many control strategies were developed for engineering applications other than resistive loading and complex-conjugate control, which have inherent difficulties in handling constraints and in achieving physically the required off-diagonal values in the damping matrix of the array. Further, in order to apply the theory for obtaining the optimal operating condition, hydrodynamic properties including wave-exciting forces, radiation damping coefficients, etc., for each individual device in the array configuration are needed. This can be a complex hydrodynamic problem itself with waves radiating and diffracting from one and another. Point-absorber theory has been used to approximate wave-interaction effects, which considers only the effect of radiation waves from one body on the others, but the diffracting waves, caused by the incident waves existing among the array, are neglected for convenience. More accurate methods include plane-wave theory (Simon 1982) where scattering waves were assumed to be plane waves; multiple scattering waves (Ohkusu 1974) where a multilevel iterative scheme was proposed to consider higher orders of scattering waves; and an algebraic exact interaction theory (Kagemoto and Yue 1986), where a direct matrix method was applied to find the hydrodynamics of the array using the diffraction properties of individual members. More recently, a generalized formulation for handling this issue was presented (Zhong and Yeung 2016). These methods can achieve different orders of interaction accuracy. However, for some, the computational cost can be prohibitively high for real-time implementation. The hydrodynamic interaction theory is not the focus here.
Recently, control of an array with physical and machinery constraints and wave-interaction effects taken into account is attracting more and more attention, which is attributed to the development of optimal control strategies for single WECs and the rapid increase in computing power. A model-based control strategy was developed and applied to unconstrained and constrained devices globally or independently (Barcelli and Ringwood 2013). The global controller optimized the performance of the whole system and operated each device in a coordinated manner, while the independent controller was used on an individual device with the knowledge of the surrounding wave field, but not of the actions of controllers on the other devices. Hydrodynamic coefficients were computed by the boundary element solver WAMIT which considers full wave interactions. It was shown that 10% or more performance improvement can be achieved for a WEC array by using the coordinated control, compared to the independent control. Li implemented decentralized, centralized, and distributed model predictive control (MPC) to an array of two WECs (Li and Belmont 2014a, b) with the point-absorber theory mentioned above to account for wave-interaction effects. For distributed MPC, each WEC was equipped with one controller and all the controllers were coordinated and run iteratively to achieve an optimal solution for the whole system. It was found by these authors that the distributed MPC can achieve almost the same amount of energy as that of the centralized MPC where the optimization was performed directly for the system; in the meantime, the former can reduce the computational burden to a large extent.
It is of interest to note that, unlike the case of a single device, the phase of the oscillating velocity of the individual cylinder in an array is not necessarily in phase with the wave-exciting force. The computational cost to find the optimal amplitude and phase of each oscillating body in an array and, most importantly, the control input to achieve the optimal conditions can be prohibitively demanding, especially with constraints being considered simultaneously. Furthermore, if a control strategy is to be implemented online, which better accounts for state changes and rejects disturbances compared to an offline control, the control strategy should be generated in a relatively short amount of time (compared with plant dynamics). MPC has shown its strong potential when applied to single WECs because of its capability in handling hard constraints on both states and the control input (see, e.g., Son et al. 2016). To overcome the challenge of the computational cost, a convex formulation of the optimization problem for the MPC was proposed (Zhong and Yeung 2017). The problem was recast as a quadratic programming (QP) problem, for which convexity can be guaranteed by tuning the weight of the penalty term on slew rates of the control input in the cost function; hence, the problem can be solved efficiently by a standard QP solver.
In this paper, such an MPC (Zhong and Yeung 2017) is applied to an array of WECs, in both regular and irregular waves. Point-absorber approximation will be adopted to simplify the consideration of wave-interaction effects. We will form a central controller to optimize the global performance of the energy absorption with constraints added on the control force and the oscillation amplitude of each device. Performance of the current MPC will be compared with a simple passive control strategy on the array. The interaction factor q of a WEC array under control will be studied, as well as the effects of constraints. In addition, reactive power required by the current MPC will be discussed. As the MPC requires that a wave-prediction unit be simultaneously incorporated to predict the dynamics of the system, we will assume in this work that there exists a wave predictor (Fusco and Ringwood 2012;Morris et al. 1998), which can estimate the wave elevation at a designated location for a certain future period. Figure 1 shows an artist conception of a farm of waveenergy devices which can act independently in a seaway or in a coordinated manner with possible communication among themselves. Ideally, we assume the incident-wave field is known and if a global controller is intelligent enough; the system may act cohesively to achieve the best global performance subject to some hardware or physical constraints. It is not difficult to envision each device may have its own controller communicating with a global controller. In what follows, we set up a framework so that the master controller would direct the members to act in the best interest of the "community."

The Equations of Motion
Without loss of generality, we consider the bodies in the array to be different in size. In the time domain, the equation of motion of the j th heaving point absorber in an array of N bodies can be written as follows: where m j and ζ j 3 (t) are mass and heave displacement of the absorber. Here  where K j is the hydrostatic stiffness coefficient, μ j 33 and λ j 33 are added-mass and damping coefficients at infinite frequency; note that λ j 33 (∞) = 0. The convolution integrals represent the fluid memory effects where K jl r (t) is the impulse response function with the superscripts jl indicating the effect of cylinder l on cylinder j . The K jl r can be obtained by the inverse Fourier transform of the damping coefficients in frequency domain (Wehausen 1971) as follows: where wave damping at infinite frequency vanishes, i.e., λ jl 33 (∞) = 0. For l = j , the λ jl 33 represents the effect of the motion of the cylinder k on the wave forces of cylinder j , which will be discussed in detail in Section 2.2. The wave-exciting force can be expressed as follows: for a given η(0, t), the incident-wave elevation at the vertical axis of the body. The impulse response function K e (t) can be evaluated by the inverse Fourier transform of the complex amplitude X 3 of the wave-exciting force induced by unit-amplitude waves of frequency σ as follows: As full knowledge of the incoming waves is presumed in this paper, we will compute the wave forces based on the frequency-domain analysis under linear wave theory, for regular waves in Section 4.1 and for irregular waves in Section 4.2. Substituting Eqs. 2a and 2b into Eq. 1, we obtain the following: where a corrective linear-viscous damping coefficient λ vis , obtained from experiments (Tom and Yeung 2014), is included to account for viscosity-based dissipation.

Point-Absorber Approximation
To model the wave-interaction effect among an array, we adopt the point-absorber approximation (Evans 1979;Falnes 1980) which simplifies the hydrodynamic interaction problem but considers the effect of the radiating waves from one device on the other devices. As indicated by McIver (1994), the basic assumption of this approximation is that the devices are small enough, relative to the wavelength of the incident waves, for the structure to be a weak scatterer. Based on this, all scattered waves were neglected in the computation of interaction effects. Hence, the waveexciting force on an individual cylinder is the same as the one on an isolated body except for a possible phase change governed by the coordinate origin of the cylinder. The damping coefficient λ jj 33 representing the hydrodynamic force induced by the motion of the j th cylinder itself will also not be affected by interacting waves. However, the "interaction" damping coefficient λ jl 33 with j = l, representing the force on the body j caused by the motion of the body k, will need to be evaluated. Following the derivation by Evans (1979), we can obtain the interaction damping coefficients λ jl 33 for the j th cylinder located at a global (cylindrical) coordinates (r j , α j ) induced by the lth cylinder located at (r l , α l ) as follows: where λ 33 is the radiation damping coefficient for the isolated cylinder, J 0 is the zero-th-order Bessel function, k is the wave number of the incoming waves, and d jl = r 2 j + r 2 l − 2r j r l cos(α m − α n ) is the distance between the j th and lth body with j, l = 1, 2, . . . , N. Note that λ jl 33 = λ 33 if j = l, and λ jl 33 = λ lj 33 . Substitution of Eq. 7 into Eq. 3 yields the retardation function K jl r (t). In addition, as derived by Evans (1979), the complex amplitude X j 3 of the wave-exciting force on cylinder j with a phase shift because of the location of the cylinder j can be obtained by the following: where β is the incident angle of the waves measured from the x-axis and X 3 for the isolated cylinder, which can be obtained from frequency-domain hydrodynamic analysis (Wehausen 1971).

State-Space Model
To form the state-space model based on the equation of motion (6), we denote the convolution of K jl r andζ 3 (t) by y jl r , and approximate it using the following linear subsystem: (l = j ) as the state vector, we can write the equation of motion for the j th cylinder as follows: and N j r = N l=1,l =j n jl r is the total number of states for the "interaction" radiation subsystems. The centralized model for a system of N WECs can be written aṡ The state-space model (12) is discretized using zeroorder hold (Hals et al. 2011b), where inputs are assumed to be piecewise constant over the sampling time T s . As a result, we obtain the following discrete-time model: where the subscript k indicates the time instant k.

Model-Predictive Control Formulation
With the state-space model constructed in Section 2, we can apply the procedure of the MPC formulation (Zhong and Yeung 2017) to build a MPC controller for multiple WECs.

The Optimization Problem
The goal is to maximize the total energy E extracted by all of the PTO systems over a predicted time horizon T h , where With the use of the discrete-time model (13a) and (13b) E can be written as E = T s This is equivalent to setting the cost function to be as follows: and minimize J 0 over the control inputs. For safety and long-term operation, constraints on the heaving motion of the absorber z and the machinery force u need to be included: where k = 0, 1, . . . , N h − 1 and j = 1, 2, . . . , N. However, it can be shown that the constrained optimization problem formed by Eqs. 16, 17a and 17b will result in a nonconvex QP (Li and Belmont 2014a), for which the global optimal solution is not guaranteed.
In this paper, we construct the cost function to be as follows: where u k = u k − u k−1 for k = 0, 1, 2, . . . , N h − 1, and || · || 2 represents the l 2 norm. u 0 := u 0 − u −1 where u −1 is defined as the control input obtained in the previous time horizon, and u −1 = 0 initially. The convexity of the current cost function J in Eq. 18 can be guaranteed when r is tuned to be larger than a certain value and the Hermitian matrix in the QP will be positive definite. This approach was previously used for a cost function with a penalty term on the control input for reducing the power consumption (Li and Belmont 2014a). The current controller, in addition to restricting consumed power, can also achieve a smaller slew rate without adding constraints on the slew rate, which can significantly increase computational costs. Nevertheless, it should be pointed out that a larger weight r on penalizing the slew rate can result in a larger loss in the energy extraction. For cases simulated here, a small r (r << 1) is sufficient to guarantee the convexity of the cost function.
With the smallest r, the impact of the term on energy extraction is negligible.

Quadratic Programming Formulation
We will cast the problem into a QP. The following notations will be used: which represent the series of inputs u and w at time instants k = 0, 1, . . . , N h − 1, where we recall that u and w are column vectors denoting the machinery forces and the wave-exciting forces on all of the WECs, as defined in Eqs. 11a and 11b. The optimization problem is composed of the cost function (18) and constraints (17a) and (17b) with the dynamics in Eqs. 13a, 13b, and 13c satisfied.
To form the QP, we propagate the state equation (13a) and use (13b) to obtainŷ(k + i|k), which denotes the state y k+i estimated at time k and stacks the output states to form the vector Y N h −1 0 := ŷ(k|k), can then be expressed in a matrix form as follows: where In a similar manner, we can express the vector of stateŝ z(k + i|k), i = 0, 1, . . . , N h − 1, denoted by Z N h −1 0 as follows: where we substitute C z for C in Eq. 21 to obtain S x,z and S u,z . In matrix form, the slew input vector U := [ u 0 ; u 1 ; u 2 ; . . . ; u N h −1 ] can be written as follows: with u −1 being the control input applied to the system currently and u −1 = 0 initially. In the MPC, u −1 is the input obtained from the previous receding-horizon computation.
The cost function J in Eq. 18 can then be written in a compact form: where R = r × I (NN h ) with I being the identity matrix. Substituting Eqs. 20 and 23 in Eq. 25 and neglecting the constant term, we can write J in a quadratic form as follows: with S R := (T ) T RT and S g = (T ) T Rg u . It can be seen that without penalizing the slew rate, the cost function will have a corresponding Hessian matrix H 0 = S u + S T u where S u in Eq. 21 has zero diagonal elements, which may lead to an indefinite H 0 , and hence a nonconvex problem. The addition of S R , aside from controlling the slew rate, adds positive values on the diagonal of the matrix H per Eq. 27a. By tuning r so that the minimal eigenvalue of H is larger than 0, we can have a positive definite H , and hence a convex J .
The constraints (17a) and (17b) can be written in a component-wise inequality. Thus, where it can be shown, with the use of Eq. 22, that with U max = f m,max × [1, . In summary, the QP with the cost function (26) subject to the constraints (28) can be written as follows:

MPC Scheme-Implementation
The MPC scheme as presented is implemented as follows: at time instant k, with the knowledge of estimated states and the predicted wave profile, the QP in Eqs. 30a and 30b is solved, and a sequence of optimal control inputs U * 0 = arg min(J ) is generated; the first set of the control inputs in the sequence is applied; the system dynamics is then moved forward to the next time step, and the optimization procedure is repeated. In this paper, we will assume the knowledge of a certain future period of the incoming waves is available.

Results of Simulations
We applied the MPC formulated in Section 3 to an array of the dual coaxial-cylinder WEC (Son et al. 2016). A nonlinear MPC (NMPC) was developed and used by these authors for the single WEC, where the PTO damping function was used as the optimizer. The performance of the current formulation of MPC on a single WEC has been compared to that of the NMPC (Zhong and Yeung 2017). It was shown that by tuning the penalty weight r, the current MPC strategy can achieve a better performance in both regular and irregular waves in terms of its ability to broaden bandwidth of the capture width, while simultaneously saves computational time and power because of the convex formulation of the optimization problem. In this work, results for two cylinders will be presented in head seas and in beam seas for both regular waves and irregular waves. The performance of the two-body array is compared with that obtained for a single cylinder. It is common to define a q factor as follows: where P u is the useful power obtained by multiple absorbers, P u,s is the power absorbed by a single or isolated absorber, N being the number of devices. For comparison, we will also present results of the WEC using a passive controller, i.e., without considering constraints, where each WEC was equipped with a simple linear damper of nonnegative damping coefficient, B g ≥ 0.
The PTO force in time domain can be obtained by the following: As in Falnes (2002), the optimal linear damping of the j th cylinder in regular waves can be obtained by the following: where λ 33 and μ 33 are hydrodynamic coefficients for the j th cylinder in isolation. For irregular waves, B g will be tuned to correspond to the peak frequency of the wave spectrum.
A dual coaxial-cylinder WEC, as shown in Fig. 2, was developed (Son and Yeung 2017). The (outer) cylinder or floater has a radius of 0.254 m and a draft of 0.635 m, and the WEC device has the total mass of 103.63 kg (including the added mass), and hydrostatic stiffness coefficient of 1 285.6 N/m. Water depth is 1.5 m. Additionally, the equivalent linear-viscous damping factor λ vis = 1.7λ 33 = 9.68 N · s/m was obtained from experiments (Son et al. 2016). The hydrodynamic coefficients for the single device with the coaxial geometry were obtained by an in-house code (Chau and Yeung 2012). A fifth-order state-space model was used to approximate its radiation subsystem to achieve an accuracy of 0.999. This impulse response K r (t) and its fifth-order approximation used is shown in Fig. 4a, as well as the damping coefficient λ 33 plotted versus the dimensional wave frequency σ .
For the array, we assumed two identical devices placed in an arrangement as shown in Fig. 3. To ensure the reasonable validity of the point-absorber theory, the center-to-center spacing L is set to L = 10a with a being the radius of the (outer) cylinder. Regarding the modeling of wave radiation forces, as λ jj 33 = λ 33 , the radiation subsystem for K jj r can be the same as the one for the isolated cylinder. A tenth-order state-space model was used to approximate the radiation subsystem for K jl r (l = j ) representing the force induced by interacting waves, where the damping coefficients were computed based on Eq. 7. The impulse response K jl r (t) and the corresponding damping coefficient λ jl 33 are shown in Fig. 4b.
As for the MPC controller, the oscillation amplitude of the absorber is constrained to be no larger than the maximal response amplitude obtained by the NMPC (Son and Yeung 2017): ζ 3,max = 5.1A 0 . In addition, the maximal machinery force is computed by f m,max = B g,max ζ 3,max σ , where σ is the angular frequency of incident waves and the maximal damping B g,max = 150 N·s/m (Son and Yeung 2017). We set the penalty weight r to be r = r/T 2 w , where the scaling factor T 2 w with T w being the period of incoming waves was introduced to reduce the effect of the time interval on the penalty of the slew rate; and r was chosen such that the cost function (18) becomes convex. The quadprog solver in MATLAB was used to solve the QP for the MPC.

Illustrative Results in a Regular Wave Environment
Under linear theory, the wave-exciting force induced by a regular wave of frequency σ on the cylinder j can be obtained by the following: where A 0 is the amplitude of an incoming wave, and X j 3 obtained by Eq. 8 is the complex amplitude of the wave-exciting force on the j th cylinder induced by a unitamplitude wave of frequency σ . For simulations, A 0 = 0.0254 m. For the discrete-time model in Eqs. 13a and 13b, the wave-exiciting force on cylinder j is obtained by the following: where w j k denotes the f j e at time instant k with k = 0, 1, . . . , N h −1. Simulations were carried out for 20T w with 80 time steps per wave period.
The nondimensional capture width C w and the maximal response amplitude operator (RAO) of absorbers in the array, plotted over the nondimensional frequency σ , are shown in Figs. 5 and 6 and compared with those obtained by the (unconstrained) passive control. The parameters were defined by as follows: with g being the gravitational acceleration, P wave the wave power, and P u (t) = 1 , indicating the energy extracted by the WEC array, which is then averaged over the N bodies. For a regular wave of amplitude A 0 and group velocity v g , the wave power can be obtained by P wave = 1 2 ρgA 2 0 v g . Results are shown for different values of penalty weight r with the minimal value r min = 10 −1.3 /T 2 w for guaranteeing the convexity of the cost function; also presented are the result for an array of two with unconstrained MPC and that of a single WEC controlled by the current MPC for reference.
Several observations can be made. Constraints largely affect the performance of the WECs at low frequency where large motion is needed to generate the optimal power. At higher frequencies, with the minimal penalty on the slew  . 7 The ratio of the reactive power to the power flowing from the absorbers to the PTO units (so-called "active power") for a two-WEC array in head seas (β = 0 • ) rate of the control input, the constrained MPC is able to produce approximately the same amount of the power as the unconstrained MPC. With an increase of the penalty weight r, less power was generated to have a milder changing rate of the input, and a smaller RAO was achieved. In general, the MPC outperformed the passive controller in terms of increasing the energy absorption, especially at frequencies far from the resonance frequency σ = 0.5668. However, the current MPC approach taking the machinery force as the control input will require power to motor the WEC system, i.e., the so-called "reactive power," while the passive control of the damping will not consume power. Figure 7 shows the ratio of the power flowing from the PTO units to the absorbers (reactive power P react ) to the power flowing from the absorbers to the PTO units (so-called "active power" P act ). The net power is the useful power, given by P u = P act − P react . It can be seen that less reactive power is required with larger r. Further, the closer the frequency is to the resonance frequency, the smaller the reactive power is; at the resonance frequency, no reactive power is needed. Similar results were found for a single WEC (Zhong and Yeung 2017). As shown in Fig. 5, at resonance frequency, all of the controllers for the WEC array achieve a similar amount of energy production which, on average, is smaller than the one generated by an isolated WEC. This is attributed to the destructive wave effect. In Fig. 8, the interaction factor q defined in Eq. 31 was shown for WEC arrays controlled by the constrained and the unconstrained MPC, and the complex-conjugate (CC) control in Evans (1979) and Falnes (1980) for comparison. It can be seen that values of the factor q of different control strategies follow a similar trend; at resonance frequency, the cylinders in beam seas (β = 90 • ) experienced constructive wave  (gray, dash-dotted) and f e on the right, for the two devices in regular waves of period T w = 1.8s (σ = 0.56). Constraints on ζ 3 are shown by dotted lines, while the constraints on F g are not active in this case effects with q > 1 while the ones in head seas (β = 0 • ) experiences the destructive effects. This suggests that the MPC is able to achieve the optimal condition of the array with constraints taken into account. Furthermore, it appears the effect of interacting waves on the energycapture ability of a WEC array compared to isolated WECs, is indepedent of the controller, with the assumption that the devices are weak scatterers. Hence, it is possible to predict the energy captured by a WEC array with a certain controller by obtaining the energy produced by a single isolated WEC with the same controller, and also the q factor from CC control, which can be easily computed. This way will be more efficient than running simulations for the whole array, and will be beneficial for the preliminary design of a large WEC array. It is also noticed that adding constraints will reduce wave effects as it restricts the heaving motion of devices which radiates waves. Time histories of the two-body array are shown in the Fig. 9 for waves of period T w = 1.8 s around the resonance frequency. It can be seen that the phases of the heaving velocity of the devices match that of corresponding wave-exciting forces, a similar result to the case of a single cylinder. However, it was pointed out (Thomas and Evans 1981) that, the phase matching at the optimal condition was not obtained for an array with k 0 L < 3 where L is the spacing between devices. Here, in order to reduce the effect of scattering waves so that the point-absorber approximation stays valid, we choose a relatively large spacing L = 10a and have k 0 L = 3.13. We believe the phase difference between the heaving velocity and the wave force will be more significant with the decrease in spacing or an increase in the number of devices.

Performance in Irregular Waves
The proposed MPC was then applied to study the performance of a WEC in irregular wave conditions. A modified Pierson-Moskowitz spectrum for fully developed seas was used to construct the irregular wave profile: where H s , T p , A j , and δ j are the significant wave height, peak period, wave amplitude, and phase angle of the waveexciting force the j th-component, respectively, and ε j is a random phase angle of the wave components. The random phase angles are uniformly distributed between 0 and 2π and constant with time. The range of angular frequencies used in the construction of Eqs. 37a, 37b, 37c, 37d, and 37e was between 0.1 and 8.75 rad/s spaced at 0.05 rad/s. The simulation runs for H s = 7.62 cm (3.0 in) and varying T p , carried out for 1 000 peak-wave periods with 40 time steps per peak period and the time horizon T h = 2T p . In this application, the RAO is constrained to be half of the significant wave height, i.e., ζ 3,max = H s /2, which was the maximal RAO used in the NMPC work (Son and Yeung 2017). The maximal machinery force is computed in the

Fig. 10
Averaged useful power in irregular waves same way as described in Section 4.1 with σ = 2π/T p . And the penalty weight r = 10 −1.1 /T 2 p is used. Results of the averaged useful power absorbed by the array in head seas (β = 0 • ) and beam seas (β = 90 • ) are shown in Fig. 10 and compared with the single-device case. The ratio of the reactive power to the "active" power is shown in Fig. 11. For both the two-device and singledevice cases, more useful power is extracted when peak wave conditions are closer to the resonance frequency T p = 1.784 s. It is observed that in general a smaller reactive power ratio is required for WECs in irregular waves than those in regular waves. This is expected as the composition of waves of different frequencies reduces the efforts required from the PTO unit to drive the absorber in phase with waves, compared to the case where the frequency of the monochromatic waves is substantially different from the resonance frequency. This fact is favorable for

Fig. 12
Interaction factor q for an array of two WECs in head seas and beam seas of irregular waves applications in real-sea environments. q factor is plotted over the nondimensional frequency for the peak waves in Fig. 12, and also compared with that obtained by using CC control. It is interesting to note that in irregular waves, the interaction factor of the WEC array does not appear to follow the trend of q obtained in the regular waves. In addition, with effects coming from waves of various frequencies, both the constructive and the destructive effects of radiating waves were much reduced.

Conclusions
We have extended the MPC for a single WEC (Zhong and Yeung 2017), for which a convex QP was formed for the optimization problem, and applied it to an array of heaving point absorbers so that an optimal control law can be obtained to maximize the global performance of the array with constraints on both the motion of devices and machinery forces satisfied. A penalty term introduced in the cost function of the MPC can ensure the convexity of the problem and penalize the slew rates of the machinery forces so that the reactive power required by the MPC can be adjusted. The wave-interaction effects among the array was taken into account in the dynamic model by adopting the point-absorber approximation.
In general, the current MPC outperformed the passive control strategy in terms of a broadening of the capturewidth bandwidth. By tuning the penalty weight r, the reactive power required by the system can be adjusted. The convex formulation allows efficient computations for an array of large number of elements, though we applied the method on only two bodies for demonstration. In regular waves, the q factor obtained by the current MPC follows the same trend of the one obtained by complex-conjugate control. This suggests that the current MPC provides a control law that optimizes the global performance of the array with physical constraints accommodated simultaneously. It should be pointed out that adding constraints on the motion reduces wave-interaction effects as it reduces the radiating waves. At resonance frequency, phase matching between the oscillation velocity and the wave-exciting force was obtained by both cylinders in the array. However, this may not be the case for arrays with less spacing distance or more elements, for which the effect of scattering waves would no longer be negligible. In irregular waves, the wave-interaction effects seems not be predictable by using the results obtained in the regular waves. With effects coming from waves of various frequencies, the q factor of the array in irregular waves of a certain (peak) frequency differs from the one obtained in regular waves of the same frequency. It is also because of the presence of multiple frequencies of wave that the reactive power required by the current MPC in irregular waves is in general lower than that in regular waves.
For future study, a more accurate model for considering wave-interaction effects (Zhong and Yeung 2016) will be adopted in the dynamics model. This theory provides an efficient method to obtain hydrodynamic properties of an arbitrary array of devices. Effects of the scattering waves on the performance of the WEC array and the control law can then be investigated and understood. We will also be looking into the energy-capture capability of arrays with different configurations under control.