Energy-maximising moment-based constrained optimal control of ocean wave energy farms

Funding information Italian Ministry for Research, Grant/Award Number: 2017YKXYXJ; Royal Society, Grant/Award Number: IEC/R1/180018; Science Foundation Ireland, Grant/Award Number: 13/IA/1886; European Union’s Horizon 2020 Research and Innovation Programme, Grant/Award Number: 739551 (KIOS CoE) Abstract Successful commercialisation of wave energy technology inherently incorporates the concept of an array of wave energy converters (WECs). These devices, which constantly interact via hydrodynamic effects, require optimised control that can guarantee maximum energy extraction from incoming ocean waves while ensuring, at the same time, that any physical limitations associated with device and actuator systems are being consistently respected. This paper presents a moment-based energy-maximising optimal control framework for WECs arrays subject to state and input constraints. The authors develop a framework under which the objective function (and system variables) can be mapped to a finite-dimensional tractable quadratic program (QP), which can be efficiently solved using state-of-the-art solvers. Moreover, the authors show that this QP is always concave, i.e. existence and uniqueness of a globally optimal solution is guaranteed under this momentbased framework. The performance of the proposed strategy is demonstrated through a case study, where (state and input constrained) energy-maximisation for a WEC farm composed of CorPower-like WEC devices is considered.


INTRODUCTION
Among the available renewable energy sources, ocean wave energy, once economically viable, can make a valuable contribution towards a sustainable, global, energy mix. Ocean waves represent a massive and untapped source of clean energy: the wave energy resource has been estimated (worldwide) to be around 3. 7 [TW] and about 32, 000 [TWh/year] in [1] and [2], which would cover ≈20% of current global energy consumption. Despite being a vast resource, wave energy conversion technology has not yet reached commercial reality. The main reason for the lack of proliferation of wave energy can be attributed to the fact that harnessing the irregular reciprocating motion of the sea is not as straightforward as, for example, extracting energy from the wind. This is clearly reflected in the current high installa-tion, operation, maintenance and decommissioning costs, which are hindering wave energy converters (WECs) in reaching economic viability.
As a direct consequence, the roadmap towards successful commercialisation of wave energy systems inherently embodies the concept of WEC arrays (sometimes called farms), which incorporate several devices in a common sea area [3]. This can effectively reduce the associated levelised cost of energy (LCoE) through an economy of scale, and hence any realistic effort to commercialise a novel device must include both a single WEC and a WEC array development process.
Though the economy of scale facilitated through the development of arrays can effectively reduce the LCoE, it is well-established that WECs require of an optimised control which can consistently ensure maximum energy extraction from incoming ocean waves (see, for instance, [4,5]). Furthermore, any realistic attempt to realise such a process must ensure that any physical limitations of both device and power take-off (PTO) system (i.e. actuator) are being consistently respected, hence minimising the risk of component damage.
Controlling arrays of wave energy devices is intrinsically challenging: Devices composing a WEC farm are normally installed in close proximity, mostly motivated by practical considerations, such as sharing of electrical infrastructure and mooring systems and space limitations [3]. This directly complicates both modelling and real-time energy-maximising control design when compared to the case of a single device, since the motion of each WEC is directly affected by waves (i.e. radiation effects) generated by adjacent devices. In the following, we present a brief review of the main optimal control strategies developed in the literature for the WEC array case, considering both early (simplified) theoretical results and advanced optimisation-based control strategies.
A pioneering, but unconstrained approach, to the control of WEC farms can be found in the early study [6], where the conditions for optimal energy absorption are presented for the case of regular wave excitation. Rapidly following this publication, [7] incorporates constraints in the motion (amplitude and velocity) of the device, presenting the first result on constrained control of WEC farms with potential practical implications. Contemporary studies in this subject include 'advanced' control strategies, such as model predictive control (MPC) [8][9][10], design based on evolutionary algorithms [11] and mean weighted residual methods [12,13]. These control strategies are briefly discussed in the following paragraph.
The study performed in [8] is an extension of the MPC strategy for a single WEC developed in [14], and presents a constrained distributed MPC formulation, where the objective function is modified with a weighting term to ensure convexity of the problem (and hence, uniqueness of the solution). The authors essentially divide the WEC array into smaller sub-farms weakly coupled to the adjacent devices, in an attempt to reduce the computational burden required by the strategy. We note that this strong computational burden can be mainly attributed to the parametrisation of radiation forces inherently required by MPC, as detailed in [5]. A similar MPC formulation to that of [8], though not distributed, can be found in both [9] and [10], where the energy-maximising objective function is also modified by the introduction of a convexification term. Using a different approach, the authors of [11] utilise a differential evolution algorithm (see, for example, [15]) to compute an optimal control law based on the current sea state, though the strategy does not consider either state (displacement and velocities of the devices in the array) nor input (PTO force) constraints, hence challenging the practicality of the proposed solution. Finally, a spectralbased controller can be found in [12,13]. Though the strategy presented in [12,13] considers the original energy-maximising objective function, we note that the authors do not explicitly guarantee existence and uniqueness of the optimal control solution, hence successful convergence towards a unique optimal law is not specifically ensured, but only demonstrated numerically with specific examples.
A novel energy-maximising optimal control framework for the case of a single WEC was recently proposed in [16]. Such a framework is based on the system-theoretic concept of moments (discussed herein in Section 2), and maps the original energymaximising optimal control problem for WECs into a concave Quadratic Program (QP), systematically guaranteeing a unique solution for the target energy-maximising control objective, subject to both state and input constraints. Subsequently, a computationally efficient calculation of the moment-based optimal control law is achieved by means of the state-of-the-art QP solvers, such as those extensively described in, for instance, [17]. Despite the fact that [16] effectively accomplishes the energymaximising control objective, subject to state and input constraints, the mathematical formalism proposed considers only the single-input, single-output (SISO) case, precluding the application of the strategy to the case in which an array of WECs is involved.
Following the array roadmap for a successful WEC commercialisation, we present, in this paper, an energy-maximising control framework for WEC arrays, exploiting the system-theoretic notion of moments presented in [18]. The hydrodynamic interactions between bodies (or devices) are fully exploited to compute the optimal control law, therefore optimally maximising the energy extraction for a WEC array from a given wave field, subject to both state and input constraints. To this end, the paper provides the following contributions: • We propose a method to map the energy-maximising objective function for WEC arrays (and system variables) to a finite-dimensional QP problem. In particular, unlike most of the model based energy-maximising control strategies reported for both single WECs and WEC farms, we show that this moment-based strategy does not require an a priori parametric approximation of the radiation force (convolution) term, but actually provides an analytical description of the convolution operation in terms of moments. This, together with the QP characteristic of the objective function in the moment-domain, renders this moment-based strategy highly efficient in computational terms, appealing for realtime applications. • We show that the resulting QP problem is always concave, i.e. we guarantee existence and uniqueness of a globally optimal energy-maximising solution for the WEC array, under state and input constraints. This systematically guarantees globally optimal performance of the moment-based strategy. • Finally, we present an extensive case study based on a WEC farm composed of CorPower-like devices [19], where the performance of the proposed controller is assessed, under state and input constraints.
The remainder of this paper is organised as follows. Section 2 discusses key concepts behind the moment-based framework for both SISO and multi-input multi-output (MIMO) systems. Section 3 formally introduces the energy-maximising problem for WEC farms, while Section 4 details the moment-based analysis of the constrained optimal control formulation. Finally, Section 5 discusses the application case described in the previous paragraph, while Section 6 encompasses the main conclusions of this study.
We note that a preliminary study on this subject has been presented in [20]. The present paper formalises the theoretical results presented in [20], while also significantly extending the application case to a full-scale WEC array composed of heaving CorPower-like devices.

Notation and preliminaries
Standard notation is considered throughout this study, any exceptions detailed in this section. ℝ + (ℝ − ) denotes the set of non-negative (non-positive) real numbers. ℂ 0 denotes the set of pure-imaginary complex numbers and ℂ <0 denotes the set of complex numbers with negative real part. The symbol 0 stands for any zero element, dimensioned according to the context. The notation ℕ q indicates the set of all positive natural numbers up to q, i.e. ℕ q = {1, 2, … , q}. The symbol n denotes the identity matrix of order n, while the notation 1 n×m is used to denote an n × m Hadamard identity matrix (i.e. an n × m matrix with all entries equal to 1). The spectrum of a matrix A ∈ ℝ n×n , i.e. the set of its eigenvalues, is denoted by (A). The superscript ⊺ denotes the transposition operator. The symbol ⨁ denotes the direct sum of n matrices, i.e.
The symmetric-part of a matrix A ∈ ℝ n×n is defined (and denoted) as ℋ{A} = (A + A ⊺ )∕2. The complex conjugate of a matrix A ∈ ℂ n×m is denoted as A. The Kronecker product between two matrices M 1 ∈ ℝ n×m and M 2 ∈ ℝ p×q is denoted by M 1 ⊗ M 2 ∈ ℝ np×mq , while the Kronecker delta function is denoted as i j , for {i, j } ⊂ ℕ. The convolution between two functions f and g over the set Ω ⊂ ℝ, i.e. ∫ Ω f ( )g(t − )d is denoted as f * g. The set of all real-valued square integrable functions is denoted as L 2 (ℝ), i.e. L 2 (ℝ) = { f ∶ ℝ → ℂ| ∫ ℝ | f (x) 2 |dx < +∞}. Let f and g be two functions belonging to L 2 ( ), where  ⊂ ℝ is closed. Then, the inner-product between f and g is given by ⟨ f , g⟩ = ∫  f ( )g( )d . The symbol e q i j ∈ ℝ q×q denotes a matrix with 1 in the i j entry and 0 elsewhere. Likewise, the symbol e q i ∈ ℝ q denotes a vector with 1 in the i entry and 0 elsewhere. The symbol n ∈ ℝ n denotes a vector with odd entries equal to 1 and even entries equal to 0.
In the remainder of this section the formal definitions of two important operators are presented, since their definition in the literature can often be ambiguous.

PRELIMINARIES ON MOMENT-BASED THEORY
We recall and extend, in this section, some of the fundamental concepts behind the so-called moment-based framework, as developed in key studies such as [18,22].

Definition 3.
The matrix C Π, with Π solution of the Sylvester equation (6), is the moment of system (4) at the signal generator (5).

MIMO case
Consider a finite-dimensional, MIMO, continuous-time system described, for t ∈ ℝ + , by the following model, given in statespace form asẋ C ∈ ℝ q×n and assume that (7) is minimal, i.e. controllable and observable.
We are now ready to present an adaptation of Lemma 1 for the MIMO case, i.e. system (7).

and S is as in Lemma 1. Consider the autonomous multiple-output signal generatoṙ
Ξ(t ) ∈ ℝ q and assume that the pair (( q ⊗ S ), Ξ(0)) is excitable [23]. Then, there is a unique matrix Π ∈ ℝ n×q which solves the 3 We focus on square systems, motivated by the WEC control application. 4 Although we assume the same dynamic matrix S for all u i to simplify the notation, each input can be driven by an independent signal generator, i.e.̇i = S i i , u i = L i i .

Sylvester equation
and the steady-state response of the output of the interconnected system is y ss (t ) = C Π Ξ(t ).
Proof. The proof follows the same arguments as in the SISO case considered in [18] and, hence, is omitted for brevity. □ Remark 2. As in the SISO case presented in Section 2, the moment for system (7) is computed in terms of the unique solution of a Sylvester equation, i.e. Equation (10).
Remark 3. We note that previous literature in momentmatching, for the MIMO case, utilises the so-called tangential interpolation framework (see, for instance, [25,26]), where moments are defined along specific directions in ℂ q . Though this tangential approach does not require 'inflation' of the matrix S as a function of the number of inputs (as in (9)), the one-to-one relationship between moments and the steadystate output response of system (7) is lost. Herein there is special interest in retaining such a relation, in spite of the consequent increase in the order of the associated signal generator (9).

ENERGY-MAXIMISING CONTROL FORMULATION
As discussed in Section 1, this paper proposes a moment-based energy-maximising control framework for wave energy farms. In the remainder of this paper, we consider the case of an array composed of N devices each with a single degree-of-freedom (DoF), to simplify the notation. Multiple DoF devices can be incorporated within this framework in a straightforward fashion with minor modelling modifications (see, for example, [27,Chapter 8]).
The energy-maximising control problem for a wave energy farm composed of N WEC devices can be informally posed as follows: compute the optimal control input (i.e. PTO force) acting on each body u i such that the time-averaged energy absorbed by the wave energy farm is maximised over a time interval  = [0, T ] ⊂ ℝ + . This energy-maximising criterion can be posed in terms of an objective function, by noting that the total useful energy converted by the PTO of each WEC in the array can be computed as whereẋ i and P denote the velocity of the ith device and the total instantaneous power of the WEC array, respectively.

Equations of motion of an array of WECs
We introduce, in this section, the fundamentals of time-domain linear modelling of arrays of WECs. Note that the assumptions considered herein are consistent across a wide variety of WEC control and estimation studies presented in the literature, such as those utilised by both WEC array estimation [28,29], and control studies [8][9][10], among others (see [5]).
The equation of motion for an array of N WECs can be expressed in the time-domain according to Newton's second law, obtaining the following hydrodynamic formulation [30] [27, Chapter 8]: where M = ⨁ N i=1 m i is the mass matrix of the buoy array with m i the mass of the ith device, and each element of the vectors the excitation force f e i , the hydrostatic restoring force f h i and the radiation force f r i acting on the ith device (i ∈ ℕ N ) of the array, respectively. 5 The control variable  (t ) is composed of the PTO forces u i (i.e. control inputs) exerted on each device.
Continuing with the description of Equation (12), the linearised hydrostatic force  h can be written as −S h , where S h = ⨁ N i=1 s h i and each s h i > 0 denotes the hydrostatic stiffness of the ith WEC of the array. The radiation force  r is modelled from linear potential theory and, using the well-known Cummins' equation [31], is where 6 ∞ = lim →+∞Ã ( ), ∞ > 0 represents the addedmass matrix at infinite frequency and , contains the (causal) radiation impulse response of each device (if i = j ) and each interaction due to the radiated waves created by the motion of other devices (if i ≠ j ). Finally, the equation of motion of the WEC array can be expressed as The dynamical system described in terms of the Volterra integro-differential equation (14), for the WEC array case, is internally stable (in the Lyapunov sense) and strictly passive with respect to the output (velocity), for any physically meaningful values of the parameters and the radiation mapping K ∶ ℝ ⟶ ℝ N ×N involved, see [27,30]. 5 Note that diffraction effects are included, within the linear potential flow framework adopted in our manuscript, in the description of the wave excitation force input. 6 See [30] for the definition ofÃ( ).

Optimal control formulation and motion constraints
As discussed throughout Section 1, any realistic attempt to design an energy-maximising optimal controller for WECs should consider both state (displacement and velocity) and input (PTO force) constraints, given that the associated unconstrained optimal solution is often unrealistic in terms of body motion and PTO force requirements (see, for instance, [5,30]).
In particular, we consider constraints on the displacement x i and velocityẋ i of each WEC composing the array, simultaneously with constraints on each PTO force u i , which can be compactly written, for all i ∈ ℕ N , as 7 With the definition of the objective function in Equation (11), the governing dynamics of the WEC farm in (14), and the set of state and input constraints defined in (15), the energymaximising optimal control problem can be posed as WEC array dynamics (14), state and input constraints (15).

MOMENT-BASED WEC ARRAY FORMULATION
To consider the moment-based theoretical framework outlined in Section 2 on this WEC array case, the equation of motion in (14) needs to be re-written in a 'suitable' structure. We propose the state-space representation: The function ∶ ℝ + → ℝ N , assumed to be the input to the system (17), is defined as 7 Note that there is no loss of generality in assuming that the maximum allowed values are the same for the N devices composing the array. This is considered to simplify the notation.
Under this specific representation, the matrices in (17) can be written, in compact form, as with each A i j ∈ ℝ 2×2 , B i j ∈ ℝ 2 defined as where  i j is the i jth element of the inverse generalised mass matrix, i.e. (M + ∞ ) −1 .
Within the moment-based framework, presented in Section 2, each ith entry of the vectors  e and  are expressed as the output of the signal generatorṡ where the dimension of S , L e i and L u i are as in (5), i (t ) ∈ ℝ and the pairs (L e i , S ) and (L u i , S ) are observable. Given the characteristics of (S ), we consider the finite set ℱ = { p } f p=1 ⊂ ℝ + and write the matrix S in block-diagonal form as where = 2 f , f ≥ 0 integer. Note that the specific structure of the matrix S in (22) is inherently motivated by standard assumptions within the field of numerical generation of ocean waves. This is further discussed in Remark 9. Finally, both the control force and excitation force vectors are expressed as the solution of the autonomous multiple-output signal generator aṡ where, without loss of generality, the initial condition of the signal generator is chosen as Ξ e (0) = N .
Remark 4. To simplify the notation used throughout the upcoming results, and to explicitly focus this manuscript on the formu-lation of a non-linear moment-based controller, it is assumed that the moment-domain equivalent L e , characterising the wave excitation  e as in Equation (17), is known, i.e. full (instantaneous and future) knowledge of  e is available over the time interval  ⊂ ℝ + . This is without loss of generality, since estimation and forecasting algorithms for  e (which are often required due to the inherent difficulty behind measuring wave excitation forces in a moving body [32]) can be incorporated straightforwardly, by following the adaptation of the moment-based representation of  e for the receding-horizon control method presented in [32, Section IV-A], without further modifications.
Remark 5. Though beyond the scope of this study, we note that sensitivity and robustness of moment-based optimal control solutions with respect to potential errors arising in the wave excitation force estimation and forecasting processes has been analysed, in the SISO case, in [32,Section V]. In addition, we refer the interested reader to [33], which derives a framework to assess the sensitivity of a general class of energy-maximising WEC controllers to wave excitation force prediction errors.
Under this selection of matrices, the moments of system (17), driven by the autonomous signal generator (23), can be computed by solving the following Sylvester equation (see Lemma 2) where Π ∈ ℝ 2N ×N and  is the moment-domain equivalent of the radiation matrix convolution term. The moment-domain equivalent of the velocity can be expressed in terms of the solution of (24) straightforwardly as  = C Π . Nonetheless, the term  depends on Π , hence we cannot yet solve (24). We now define the quantity  and then provide an explicit solution for (24). (13) can be computed as

Proposition 1. The moment-domain equivalent of the convolution integral in
where each ℛ i j ∈ ℝ × is a block-diagonal matrix defined as whereÃ( ) i j is the added-mass matrix,B( ) i j is the radiation damping matrix 8 of the device at each specific frequency induced by the eigenvalues of S, and ∞i j is the i j th entry of the matrix ∞ .
Proof. The proof can be carried out analogously to the SISO case in [16] and, hence, is omitted. □ Remark 6. As briefly discussed in Section 1 and explicitly expressed in Proposition 1, this moment-based strategy does not require an a priori parametric approximation of the radiation force (convolution) term, but actually provides an analytical description of the convolution operation in moment-domain. This characteristic, together with the QP formulation presented in this section, renders this moment-based strategy highly efficient in computational terms and, hence, appealing for realtime applications.
With the analytical definition of the moment-domain equivalent of the radiation force convolution term in (25), we state the following two propositions, which address the uniqueness of the solution of the Sylvester equation (24) and the explicit computation of the moment equivalent .

Proposition 2. The solution of the Sylvester equation (24) is unique if and only if
where the matrix G p is defined as Proof. See the Appendix for the proof. □
Proof. Recall that  = C Π . Then, Equation (31) follows directly from (A.1) (see Section 6). □ Remark 7. Equation (28) always holds for the WEC array case: it follows from the internal stability of (17) (see Section 3.1) that (G p ) ⊂ ℂ <0 for all p ∈ ℕ f , with the matrices G p as in Proposition 2.
Remark 8. Note that, given the structure of the matrices L u and L e in (23), the moment-domain equivalent  can always be i ∈ ℝ denotes the moment-domain equivalent of the velocity of the ith device.
Propositions 2 and 3 explicitly show how to compute the (unique) moment-domain equivalent of the output of system (17), i.e. the velocities of the WEC array elements. With this last result, we address the formulation of (16) using a moment-based approach. Specifically, we show that the energy-based objective function  can be greatly simplified under the proposed moment-based framework. (11) and the representation for u i as in (21). Define the set ℱ considered to compute S in (22)

Proposition 4. Suppose (28) holds, and consider the expression for the instantaneous power P in
Then, the absorbed power  over the time period  = [0, T ], with T = 2 ∕ 0 , can be computed as where  i denotes the moment-domain equivalent of the velocity of the ith device.
Proof. See the Appendix for the proof. □ Remark 9. The selection of the set ℱ in Proposition 4 follows from a standard assumption in the numerical generation of ocean waves: the so-called free-surface elevation, fully characterising the input wave, can be described as a finite sum of f harmonics of a (sufficiently small) fundamental frequency 0 (see, for instance, [34]).
Note that Proposition 4 explicitly shows that, under the presented moment-based strategy, the objective function of (11) can be computed as the sum of N inner-product operations in ℝ 1×N . Moreover, under the presented moment-based strategy, the (unconstrained 9 ) optimisation problem associated with (16) has a strictly concave QP formulation, as detailed in the following proposition. (16). Then, under the same assumptions of Proposition 4, the optimal control law  opt = L opt u Ξ can be computed in the moment-domain as

Proposition 5. Consider the (unconstrained) energy-maximising optimal control problem
Proof. See the Appendix for the proof. □ Proposition 6. The QP formulation in (34) is strictly concave for any physically meaningful values of the parameters of (17).
Proof. See the Appendix for the proof. □ Remark 10. Propositions 5 and 6 have a strong impact on the practicality of the moment-based solution proposed in this study: the original optimal control formulation in (16) can be transformed into a QP program, which always has a unique (global) maximiser due to the fact that strict concavity is guaranteed. Hence, we can use well-known and highly efficient stateof-the-art quadratic programming solvers [17].

State and input constraints in moment-domain
Following the moment-based framework for a single device proposed in [16], we map the set of motion constraints using their respective moment-domain equivalents 10 , i.e.
⊂ ℝ + be a set of uniformly spaced time instants. We propose to enforce the set of constraints defined in at the set  : Defining the matrices Λ ∈ ℝ NN c ×N 2 and Δ ∈ ℝ 2NN c ×N 2 as and substituting  using (31), the motion constrained energymaximising optimal control law can be written in momentdomain as an inequality-constrained QP problem, i.e.: Note that the moment-domain equivalent of the position x i (t ) can be expressed [35] as The moment-based constrained QP formulation of (37) offers a globally optimal solution for the energy-maximising problem for WECs under state and input constraints. The performance of such a formulation is analysed in Section 5, both in terms of energy absorption, and constraint satisfaction.

APPLICATION TO A WEC ARRAY
We present, in this section, an application case to illustrate the proposed moment-based strategy, based on the regularpolytope-type WEC array layout depicted in Figure 2, composed of N = 5 converters. Each of the five devices composing this WEC farm is a full-scale CorPower-like device oscillating in heave (translational motion). Such a device is illustrated in Figure 3, with its corresponding physical dimensions specified in metres. Moreover, to fully characterise the wave farm, Figure 4 presents the hydrodynamic characteristics of the WEC array considered in this application case in terms of its corresponding radiation damping and radiation added-mass matrices, i.e.B( ) andÃ( ), respectively. Note that, due to the fact that the devices composing the WEC farm are identical (i.e. CorPower-like devices), the corresponding hydrodynamic characteristics (including interactions due to radiation effects) present symmetrical behaviour, in accordance with the layout depicted in Figure 2. That said, only three elements of the matrices {B( ),Ã( )} ⊂ ℝ 5×5 are required to completely characterise the hydrodynamic parameters of the farm. These are plotted in Figure 4, along with the corresponding symmetry pattern 11 for both matricesÃ andB. We note that, consistent with the main literature in WEC array control (see Section 1), we consider, for the remainder of this section, the linear model   (14) (with the parameters corresponding to the CorPower-like WEC array presented in Figure 4), for both design and performance evaluation of the proposed momentbased control strategy. Given that our primary objective is the development of a novel and efficient optimal control strategy for WEC arrays with guarantees of globally optimal solutions, the application case presented in this section is intended as a proof of concept for our novel controller, rather than an assessment in a 'realistic' environment.
Remark 11. For the next simulation results, the normalised runtime (i.e. ratio between the time required to compute the energymaximising optimal control input for the duration of the simulation, and the length of the simulation itself) of the proposed WEC array moment-based controller is always less than a second (for a MATLAB-based ; more customised coding can likely reduce this by an order of magnitude), being consistent with the typical sampling time of a full-scale WEC (see, for instance, [16]). Note that real-time application of the proposed controller can be effectively performed in a receding-horizon fashion, by directly following the moment-based method described in [32, Section IV] (see also Remark 4).
The performance assessment of the presented momentbased strategy initially considers the case of regular waves, taking into consideration both state and input constraints. We remind the reader that, as discussed in Section 3.2, the considering motion constraints is required due to the fact that the unconstrained energy-maximising optimal solution often requires unrealistic values for the physical variables of the WEC system [30]. Constraining the motion of the device, however, can consequently lead to a decrease in the total absorbed energy. This motivates us to analyse the effect of enforcing the set of state and input constraints, defined in Section 3.2, in terms of total power absorbed by the WEC farm presented in 2, when using the moment-based strategy proposed in this paper. In Then, we consider the following power absorption ratio as a performance indicator: where  con,R A|C T is the total power absorption for a regular wave of period T with either displacement ( con,R A ) or control force ( con,R C ) constrained to Figure 5 illustrates the results obtained for R P with varying wave period T , and both displacement and control force constraint factors R A (left column) and R C (right column), respectively. Furthermore, the results presented herein are for three different wave heights, i.e. H ∈ ℕ 3 m. A key element to high-light from Figure 5 is that the proposed moment-based strategy is able to maintain a constant performance with respect to H , giving extremely similar power absorption ratio results for the full set of analysed wave heights. Focusing on the left column of Figure 5, where the displacement of the device is constrained following Equation (39), it is noteworthy that with a constraint of 40% of the optimal unconstrained motion the energy-maximising moment-based strategy is capable of extracting ≈80% of the unconstrained optimal result for the totality of the analysed periods, being almost 90% for some values of T . Similar behaviour can be appreciated in the right column of Figure 5, where now the maximum PTO force is constrained within the optimal energy-maximising control computation, as in Equation (40). Note that the deterioration in power performance becomes higher in the case in which the PTO force (control input) is constrained, while a milder effect can be appreciated in the case of displacement constraints. We note that this is consistent with previous results, such as those reported in [7] (simplified theoretical analysis) and [12,13] (numerical assessment).
Completing the results for regular excitation, Figure 6 illustrates the WEC array motions under optimally controlled conditions (left column), along with each corresponding momentbased energy-maximising control laws (right column). The input wave is considered to have a wave height H = 2 m and a period T = 8 s. The state and input constraints, for each device composing the array, are set as follows: ⋄ Maximum allowed displacement X max = 2 m; ⋄ Maximum allowed velocity V max = 2 m/s; ⋄ Maximum control force U max = 1 × 10 6 N.
More precisely, the left column of Figure 6 shows displacement (solid black), velocity (dashed black) and wave excitation force input (dotted grey), for each device composing the array, from device 1 (first row) to device 5 (last row). The constraint limits for displacement and velocity are denoted with a dash-dotted red line.
We note that there are some key features that can be appreciated in the left column of Figure 6, which we detail in the following. To begin with, it is straightforward to notice that the state constraints are being consistently respected for the totality of the devices composing the WEC array, illustrating the capability of the moment-based strategy to maximise energy absorption while respecting the physical limitations of each device. Moreover, we note that even in this fully constrained case, the velocity of the device under optimal control conditions remains in-phase with the wave excitation force, agreeing with the well-known theoretical results for unconstrained energy-maximisation of (single) WECs [30]. The right column of Figure 6 presents the control inputs for each device computed with the momentbased strategy (solid black), used to elicit each corresponding motion results, along with each wave excitation force (dotted grey). Once again, it can be appreciated that the PTO force constraints (dash-dotted red) are being respected consistently, showing the ability of the strategy to handle both state and input constraints simultaneously. Finally, we also note that the optimal  [N]

FIGURE 6
Results for regular wave excitation. The left column of Figure 6 shows displacement (solid black), velocity (dashed black) and wave excitation force input (dotted grey), for each device composing the array, for devices 1 (top) to 5 (bottom). The right column of Figure 6 presents (in the same order) the corresponding control inputs for each device computed with the moment-based strategy (solid black), used to elicit the corresponding motion results, along with the wave excitation force input (dotted grey). The dash-dotted red lines represent constraint values control force is shifted by ≈ ∕2 [rad] with respect to the wave excitation force input, also agreeing with the theoretical (unconstrained maximum) power absorption conditions [30] for an isolated WEC device. We now presents results under irregular wave excitation, randomly generated using a Joint North Sea Wave Project (JON-SWAP) spectral density function [36], following the methodology presented in [34]. In particular, we consider a JONSWAP spectrum analogous to the regular excitation case presented in Figure 6, i.e. with peak period T p = 8 s and significant wave height H s = 2 m. The peak enhancement factor is set to = 3.3 (see [36]). Both the state and input constraints for each device are also set to the exact same values as those for the regular excitation case of Figure 6, i.e. X max = 2 m, V max = 2 m/s and U max = 1 × 10 6 N. Figure 7 presents motion (left column) and energymaximising control input (right column) results for this irregular wave case. We note that this figure is analogous to Fig-[

FIGURE 7
Results for irregular wave excitation. The left column of Figure 7 shows displacement (solid black), velocity (dashed black) and wave excitation force input (dotted grey), for each device composing the array, for devices 1 (top) to 5 (bottom). The right column of Figure 7 presents (in the same order) the corresponding control inputs for each device computed with the moment-based strategy (solid black), used to elicit the corresponding motion results, along with the wave excitation force input (dotted grey). The dash-dotted red lines represent constraint values ure 6, and the same indexing to variables and devices is used. We begin by noting that, as can be appreciated from Figure 7, the moment-based strategy is able to maximise energy absorption, while systematically respecting both state and input constraints for the case of this irregular wave, according to the control design objective, and hence providing a strong practical result in a realistic sea description. We also note that the velocity and excitation force of each device presents the in-phase 12 optimal energy absorption condition for regular unconstrained motion. In fact, this behaviour is consistent with what has been reported previously in the energy-maximising moment-based strategy presented in [16], for the isolated (single) WEC case.
To finalise the results under irregular wave excitation, Figure 8 shows performance results in terms of absorbed energy, for the energy-maximising moment-based strategy presented in Absorbed energy for the moment-based energy-maximising control strategy presented in this paper this paper, constrained in both displacement (X max = 2 m) and velocity (V max = 2 m/s). The time-length of the simulation is set to T = 120 s. Note that, since the waves are generated from sets of random amplitudes [34], it is found that a mean of ≈30 simulations (per sea state) is necessary to obtain statistically consistent performance results, in terms of absorbed energy.

CONCLUSIONS
This study formally introduces a moment-based energymaximising technique for WEC farms, providing a mathematical framework for array optimal control design with strong practical value, thus helping in the roadmap towards successful commercialisation of WEC technologies. We demonstrate that the optimal control problem for WEC arrays can be substantially simplified using moment-based theory, explicitly mapping the objective function to a strictly concave finite-dimensional QP problem, by means of moments. In addition, the paper details how to systematically handle both state and input constraints, by making explicit use of the advantages inherently present in the moment-domain formulation. The combination between energy-maximisation, successful simultaneous state and input constraint handling and computational efficiency (due to the nature of the objective function in moment-domain) has strong practical advantages, providing an optimal control framework that can maximise energy absorption from incoming sea waves, respect intrinsic physical limitations, and compute in realtime. Finally, this paper demonstrates the use of the proposed method by means of a full-scale WEC array composed by five CorPower-like devices, explicitly assessing the performance of the moment-based strategy for both regular and irregular wave excitation, under both state and input constraints.

Proof of Proposition 2
A direct application of the vec operator to Equation (24) (and considering Property 1 and the bilinearity and associativity property of the Kronecker product) yields the equivalent linear system of equations where the matrix Φ ∈ ℝ 2N ×2N is defined as It is straightforward to conclude from (A.1) that the solution of the Sylvester equation (24) is unique if and only if 13 0 ∉ (Φ ). As a consequence of the block-structure of each of the matrices involved in (A.2), we can always write the matrix Φ in a blockdiagonal structure, i.e. Φ = ⨁ f p=1 Φ p . Therefore, the matrix Φ is invertible if and only if each block Φ p is invertible.
After some algebraic manipulations of (A.2), each block composing Φ can be expressed as , (A.3) where ir p and im p are defined in (30). Consider now the matrix .4) and the similarity transformation W Φ p = W Φ p W −1 which yields Since the spectrum of a matrix remains invariant under a coordinate transformation, we can conclude from (A.5) that Φ p is invertible if and only if 14 Finally, the claim follows repeating the same analysis for each of the blocks of Φ p with p ∈ ℕ f .

Proof of Proposition 4
If (28) holds, the objective function  in (11), which is defined over the time period  , can be expressed in terms of  and L u as Proof of Proposition 5 We first note that, since  , where considering the bilinearity and associativity property of the Kronecker product, the matrix M can be equivalently written as from which the claim follows.

Proof of Proposition 6
Note that the quadratic program (QP) defined in (34) is strictly concave if and only if ℋ{ N ⊗ Φ ℛ ⊺ } = N ⊗ ℋ{Φ ℛ ⊺ } is positive-definite. Let (t ) =  e (t ) −  (t ) = [id 1 (t ), … , id N (t )]. Recall that, since system (17) is strictly passive, the relation [39] holds for any time interval [t 1 , t 2 ] ∈ ℝ. Assume, without any loss of generality, that the time interval is set to [0, T ]. Then, it follows from the representation of the input (t ) as in (23), i.e. = L id Ξ, Proposition 4, and the passivity condition in Equation (A.15), that the relation 16) holds. Performing the same analysis as in Proposition 5, the proof of our claim follows noting that which holds if and only if ℋ{Φ ℛ ⊺ } is positive-definite.