Parametric representation of arrays of wave energy converters for motion simulation and unknown input estimation_ A moment-based approach

When it comes to parameterisation of dynamical models for arrays of Wave Energy Converters (WECs), the mostused approach, within the wave energy literature, provides a state-space representation whose order (dimension) increases quadratically with the number of devices composing the WEC array. This represents a major drawback for key WEC design elements, such as motion simulation and unknown input estimation, being the latter essential to effectively maximise energy extraction from ocean waves. We present herein a multi-input, multioutput (MIMO) parameterisation strategy based on a system-theoretic interpretation of moments. The state-space representations computed with this moment-based approach exactly match the steady-state behavior of the target WEC system at specific (user-selected) interpolation points, providing efficient low dimensional models that can accurately represent the input-output dynamics of WEC arrays. Moreover, we show that there exists an intrinsic connection between wave excitation force estimation strategies and the moment-based parameterisation method proposed in this paper. We exploit this mathematical correlation to provide low order models that deliver the same degree of wave excitation force estimation accuracy to that obtained by implementing the currently-used parameterisation methods, with mild computational requirements. The performance of the strategy is analysed in terms of a case study, considering a WEC array composed of state-of-the-art CorPower-like devices, for both WEC motion simulation and wave excitation force estimation scenarios.


Introduction
Among the different renewable energy sources, wave energy has one of the highest power densities [1]. Despite being a vast resource, the cost involved in generating power from ocean waves is still prohibitive, compared to state-of-the-art wind or solar technologies [1]. To be precise, the current high installation, operation, maintenance, and decommissioning costs are hindering these novel wave energy extraction technologies from reaching a commercialisation stage. As a direct consequence of this, the roadmap towards successful commercialisation of Wave Energy Converters (WECs) naturally embodies the development of so-called WEC arrays or farms, which effectively incorporates several devices in a common sea area [2]. The devices composing a WEC farm are commonly installed in close proximity, mainly motivated by the underlying practical considerations, such as space limit, cable deployment, electricity delivery, and general maintenance.
Given that each WEC composing an array represents not only a wave absorber but also a wave generator [3,Chapter 8], the motion of each WEC is directly affected by the waves generated by adjacent devices (i.e. radiation effects). This feature further complicates the modelling task when compared to the case of a single device, often producing mathematical representations that are not tractable for key components facilitating the maximisation of energy absorption from ocean waves, such as real-time optimal controllers and wave excitation force estimators [4]. In light of this, at this development stage, modelling accurate but simple dynamical models is crucial for both the optimisation of the different components of each isolated WEC and layout of the array (i.e. optimal position of devices as a function of their hydrodynamic interactions), and the development of wave excitation force estimation strategies, which are essential in the design of optimal controllers that can maximise power conversion from ocean waves [4,5].
Boundary Element Methods (BEMs), readily available in up-to-date codes such as WAMIT or NEMOH [6], are the most commonly used methods to characterise the hydrodynamics of WECs. Nevertheless, the so-called hydrodynamic coefficients, obtained with these methods, are computed in the frequency-domain, characterising only the steady-state motion of the device under analysis. This represents a drawback for wave excitation force estimation strategies, which inherently require a time-domain representation of the WEC dynamics.
The motion of each floating device composing the WEC array can be expressed in the time-domain in terms of the well-known Cummins' equation [7], whose parameters can be computed from the frequencydomain hydrodynamic coefficients obtained by BEM solvers, following the theoretical framework presented in [8]. The resulting time-domain formulation is a Volterra integro-differential equation [9] which contains a convolution term accounting for the memory effects associated with radiation forces acting on the device. If we take into account that the WEC array is composed of N devices, and that WECs (hydrodynamically) interact between each other due to radiation effects, the total number of convolution operation rises to N 2 . In other words, a non-parametric impulse response function is required to model each cross-coupling term arising from radiation effects between devices.
The mere existence of these non-parametric convolution terms represent a drawback for most wave energy applications, including motion simulation and wave excitation force estimation theory perspectives. From a pure simulation point of view, it is well-known that the explicit computation of the convolution operator is computationally inefficient, often worsened by the necessity of a small (time) discretisation step to obtain accurate numerical integration. 1 Concerning modern state estimation techniques, they are virtually always based on the availability of a state-space description (i.e. a set of first-order differential equations) of the system under analysis [4].
Motivated by these aforementioned drawbacks, researchers often seek a parametric approximation of this radiation force subsystem in terms of a linear time-invariant representation, making explicit use of the hydrodynamic parameters computed with BEM-based codes. To be precise, the prevailing methodology is to approximate each of the radiation force-related convolution terms as a single-input single-output (SISO) system, even though the problem is inherently of a multipleinput multiple-output (MIMO) nature, as a consequence of the characteristics of the WEC array. Examples that consider this approach can be found in [5,11,12], among others. This so-called herein "multi-SISO" approach often leads to an unnecessarily high order parameterisation of the radiation force dynamics, leading to computationally inefficient representations. In fact, in this "multi-SISO" method each convolution operator is approximated with a parametric form of order n , r yielding an input-output representation of the WEC array of order + N n N 2 r 2 . Though computationally more efficient than solving N 2 convolutions, the final dimension (order) of the WEC array model increases quadratically with N, potentially rendering any wave excitation force estimation strategy intractable for real-time energy-maximising control applications [4].
As discussed earlier in this section, the estimation of wave excitation force is a key element for the optimal extraction of energy from ocean waves. Explicitly, the determination of an optimal control force that maximises energy absorption always requires 2 knowledge of the instantaneous and future wave excitation force e acting on each WEC composing the array [4,13]. For a fixed body, e can be measured by integrating the pressure over the submerged body surface. This is certainly not the case when considering WECs, given that the wave excitation force is an unmeasurable quantity for a moving body, since the integrated pressure over the submerged WEC no longer represents e solely. To overcome this issue, several unknown input estimation strategies have been proposed within the wave energy literature. Particularly, the strategy presented in [5] has been proven to outperform stateof-the-art unknown input (observer) strategies applied in the wave energy field, both in terms of estimation accuracy and computational time [14]. This technique is based on optimal state estimation theory, and presents a combination of Kalman filtering [15] and the internal model principle 3 of control theory [16]. Note that, as highlighted previously, this estimation strategy requires a parameterisation of the input-output dynamics of the array in terms of a state-space representation, to successfully design a Kalman observer.
Following the roadmap towards successful WEC commercialisation, and the discussion offered in the previous paragraphs, the objective of this study is twofold. The model parameterisation problem for wave energy converter arrays is addressed first, proposing an extension of the SISO moment-matching-based identification framework developed in [17], for MIMO systems. Moment-based methods interpolate a certain number of points on the complex plane, termed moments, which are directly related to the frequency response of the target dynamical system under analysis. In fact, the transfer function of the approximating model obtained by this moment-based approach exactly matches the steady-state behavior of the target system at these specific interpolation points, which are user-selected. Secondly, and taking explicit advantage of the (frequency) interpolation feature of our moment-based strategy, the existence of an intrinsic connection between the wave excitation force estimation strategy presented in [5], and the moment-based parameterisation method proposed in this paper, is shown. This mathematical correlation can be exploited to compute exceptionally low-order models which provide the same degree of wave excitation force estimation accuracy as the currently implemented parameterisation methods, but with a significant improvement in terms of computational requirements. This has strong practical implications, being particularly appealing for real-time (combined) estimation and optimal control of wave energy farms, hence contributing to the roadmap towards successful commercialisation of WEC technologies. Note that a preliminary study on this subject was presented in [18]. This paper extends the results of [18], while also showing and exploiting the synergy between moments and the unknown input wave excitation force estimation problem. In addition, the application case is extended to a full-scale WEC array constituted by state-of-the-art heaving CorPower-like [19] devices.
The remainder of this paper is organised as follows. To keep this paper reasonably self-contained, Section 2 recalls key concepts behind the moment-matching framework for both SISO and MIMO systems. Section 3 presents the dynamics of a WEC array written in a suitable structure for the theoretical elements of this paper, while Section 4 presents a moment-based analysis of such a system. Section 5 proposes a moment-matching-based methodology to compute finite-order parametric models for the force-to-motion dynamics of the WEC array. Section 6 deals with a case study, for both motion simulation and wave excitation force estimation, explicitly showing the appealing characteristics of the moment-based parameterisation proposed in this paper. Finally, Section 7 encompasses the main conclusions of this study.

Notation and preliminaries
Standard notation is considered through this study, with 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 a negative real part. The notation q indicates the set of all positive natural numbers up to q, i.e. = q {1, 2, , } q . The symbol 0 stands for any zero element, dimensioned according to the context. The symbol n denotes the identity matrix of the space × n n . The spectrum of a matrix × A , n n i.e. the set of its eigenvalues, is denoted as λ(A . The symbol × e ij q q q denotes a matrix with 1 in the ij-entry and 0 elsewhere. Finally, the symbol n n denotes a vector with odd components equal to 1 and even components 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. Definition 1. [20] (Kronecker sum) The Kronecker sum between two matrices P 1 and P 2 , with × P n n 1 and × P , k k 2 is defined (and denoted) as where p j n with j , m the vector valued operator vec is defined as  (2)

Moment-based theory
With the final aim of keeping this paper reasonably self-contained, we briefly recall some of the key concepts behind the moment-matching framework, as developed in studies such as [21][22][23], for both SISO and MIMO systems. In particular, we make a special emphasis on the formal definition of moment, using a system-theoretic approach. Note that the formalism behind moments has been originally proposed within the field of model order reduction, being adapted for the wave energy parameterisation case in [17].

Moments for SISO systems
Consider a finite-dimensional, SISO, continuous-time system described, for , by the state-space model 4 and assume that (3) is minimal (i.e controllable and observable). The idea of the moment-based model order reduction technique is based on interpolating the transfer function of the original system (and the derivatives of this) and the transfer function of the reduced-order model (and the derivatives of this) at these specific interpolation points s i .
One of the key results behind this framework is developed in [21], where the moments of system (3) are shown to be in a one-to-one relation with the steady-state response (provided it exists) of the output of the interconnection between a signal generator (sometimes referred to as exogenous system [25]) and system (3).

Lemma 1 ([21,22]). Consider system (3) and the autonomous signal generator
Then, there exists a one-to-one relation between the moments η 0 (s 1 ), η 0 (s 2 ), … , η 0 (s ν ), with s i ∈ λ(S) for all i , and the steady-state response CΠξ of the output y of the interconnection of system (3) with the signal generator (5) (as in Fig. 1). In fact, the moments are uniquely determined by the matrix CΠ.

Remark 3.
From now on, we refer to the matrix CΠ ≡ Y as the momentdomain equivalent of y.

Moments for MIMO systems
Consider a finite-dimensional, MIMO, continuous-time system described, for , by the state-space model with 6 x t ( ) , and assume that (7) is minimal. We now recall an adaptation of Lemma 1 for the MIMO case, which shows an explicit connection with the signal generator defined for the SISO case (5), in terms of the operator S.

Lemma 2 ([13]). Consider system (7) and the autonomous multiple-output signal generator
, 0 and assume that the triple of matrices L S ( , , (0)) q is minimal. Let × n q be the (unique) solution of the Sylvester equation Then, there exists a one-to-one relation between the moments η 0 (s 1 ), η 0 (s 2 ), … , η 0 (s ν ), with s i ∈ λ(S) for all i , and the steady-state response CΠΞ of the output y of the interconnection of system (7) with the signal generator (8).
Herein, we are interested in retaining the exact same steady-state response for the WEC array, in spite of the consequent increase in model order with respect to the tangential approach. The rationality behind this argument stems from the fact that retaining an exact steady-state response in the parametric representation, greatly increases the accuracy of unknown input estimation procedures, as further discussed in Section 6.

Remark 5.
Analogously to the SISO case, the moment for system (7) is computed in terms of the unique solution of a Sylvester equation, i.e. Eq. (9).
Based on this steady-state characterisation of moments, we recall the following key definition from [21,22], adapted for the MIMO case. Definition 4. [21,22] Consider the signal generator (8). The system described by the equations q q is a model of system (7) at S if system (10) has the same moments at S as system (7).

Lemma 3 ([21,22]).
Consider system (7) and the signal generator (8). Then, the system (10) is a model of system (7) where = Y C is the moment-domain equivalent of the output of system (7) computed from (9) and P is the unique solution of the Sylvester equation Remark 6. The steady-state output of the reduced order model (10) exactly matches the steady-state output of the system resulting from the interconnection of system (7) and the signal generator (8).

Governing motion equations for a WEC array
This section aims to introduce the basics behind linear modelling of arrays of wave energy converters in the time-domain, according to Cummins' equation [7]. The modelling assumptions considered herein are consistent across a wide variety of WEC applications presented in the literature, such as those concerning both WEC array estimation [5,11], and control studies [12,27,28], amongst others (the reader is referred to [4] for a comprehensive list of related studies and the modelling hypotheses considered).

Time-domain formulation
The equation of motion for an array of N wave energy converters can be expressed in time-domain according to Newton's second law, obtaining the hydrodynamic formulation [29], [ 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 [7], is ( ), μ ∞ > 0 represents the added-mass 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 where the (total) external input U + : N is composed of both the wave excitation force, and the control input, i.e. U = e . The internal stability of Eq. (15), for the WEC case, has been analysed and guaranteed for any physically meaningful values of the parameters and the radiation mapping K involved, see [29].

Remark 7.
We note that, if required by the application, a direct parameterisation of the control input can be used within this same framework, i.e. one can parameterise [17,30] as can be obtained using different optimal or suboptimal control strategies (see, for example, the strategies reported in [4]).

Remark 8.
Though excluded in this paper (motivated by a large number of studies considering the same modelling assumptions discussed at the beginning of this section), linearised viscous forces [31] can also be included straightforwardly in Eq. (15). In addition, we note that nonlinear viscous effects can be accommodated within the parameterisation framework described in Section 5, by following [32].

Frequency-domain formulation
As discussed in Section 1, BEM solvers are widely considered in the literature to characterise (hydrodynamically) WEC devices, mainly due to their efficient computational performance, when compared to timedomain methods [33]. Given that these hydrodynamic codes provide a frequency-response charaterisation of the device under analysis, it is natural to regard the motion of the WEC array using a frequency-domain description. Applying the Fourier transform to (15), and considering the velocity of each device as measurable outputs i.e. the mapping , the following representation where × H: denotes the force-to-velocity frequency response mapping of the WEC array, holds. Note that the expression force-tovelocity (or more generally, force-to-motion) is used here to denote the frequency response of the WEC considering the total external force as the input to the system. H(jω) can be computed as a function of a well- Fig. 1. Schematic of the interconnection between the system (3) and the signal generator (5) (adapted from [21]). 7 where B(ω) and A(ω) represent the radiation damping, and the radiation added mass matrix of the device, respectively. These frequencydependent device parameters are virtually always calculated using hydrodynamic codes at a finite set of uniformly-spaced frequency where ω l and ω u represents the lower and upper bound of the range, respectively. The selection of such a frequency range depends explicitly on the application, as discussed in [17].

Moment-based WEC array formulation
To consider the theoretical results recalled in Section 2 for this WEC array case, the equation of motion presented in (15) needs to be rewritten in a more suitable structure. The following state-space representation, for the WEC array dynamics, is proposed: i i i . The function u + : , N assumed to be the input of system (19), is defined as Under this assumption, the matrices in (19) can be written in compact form as follows: where ij is the ijth element of the inverse generalised mass matrix + M µ ( ) 1 . Within the moment-based framework, the input mapping U (i.e. the total external force), is expressed as an autonomous multiple-output implicit form signal generator as where the dimension of S is as in (8) and write the matrix S in a real block-diagonal form as Under this selection of matrices, the moments of system (19), driven by the signal generator (23), can be computed by solving a specific Sylvester equation (see Theorem 2). Such an equation can be specialised for the WEC array case as where × N N 2 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 (25) straightforwardly as = C . Aiming to keep this paper reasonably selfcontained, we recall the following lemma from [13], which defines the quantity .

Lemma 4 ([13]
). The moment-domain equivalent of the convolution integral in (14) can be computed as where A ij (ω) is the added-mass matrix, B ij (ω) is the radiation damping matrix of the device at each specific frequency induced by the eigenvalues of S, and μ ∞ij is the ijth entry of the matrix μ ∞ .
Following Lemma 4, we now present a compact result to directly compute , by solving (25), with defined as in (26).

Lemma 5.
Suppose system (19) is asymptotically stable in the Lyapunov sense. Then, the moment-domain equivalent of the output y φ of system (19), i.e. , can be uniquely computed as where the matrix R × N N is defined as Proof. The proof of this lemma follows from a direct application of the vec operator (as in Definition 2) to Eq. (25) and, hence, it is omitted. □

Parametric models achieving moment-matching
In this section, we pose an algorithm to compute an approximated parametric time-domain model of the input-output dynamics (force-tomotion) of the WEC array, based on the moment-domain representation of the N-device WEC array dynamics derived in Section 4. In particular, we consider the theoretical structure behind the family of systems (10) achieving moment-matching, in synergy with some of the main notions behind frequency-domain subspace-based identification methods, as proposed in [34] and briefly recalled in the subsequent paragraph.
Following the well-established theory on subspace-based methods reported in [34], both the dynamic and output matrix from system (19) can be approximated in terms of the singular value decomposition of the Hankel matrix Ĥ , constructed 9 from the input-output frequencydomain data of the WEC array (18) computed at the finite set of 9 We refer the reader to [34] for further details on the computation of the Hankel matrix Ĥ .
where the singular value decomposition of the Hankel matrix is written as Remark 11. Note that, with the knowledge of the frequency sampling associated to the set Ω, the continuous-time equivalent matrix Â can be computed directly from Â d using, for instance, the bilinear (Tustin) mapping, as discussed in [34].

Remark 12.
If Â , d computed as in (31), has unstable eigenvalues, one can always project such a set into the complex unit circle, following the procedure described in [34]. This guarantees the computation of a Jury matrix Â , d if required.
We are now ready to pose an algorithm to compute a momentmatching based model of the WEC array, exploiting the result of Proposition 5 and the system-theoretic structure of (10).
Algorithm 1 is based on the idea of building the parametric model F matching the f (user-defined) frequencies of the set F , exploiting the system structure of (10) and solving for an equality-constrained optimisation problem. Summarising, this optimisation process aims for the computation of the input matrix G opt that minimises the difference between the target frequency response and that of (10) (in terms of the matrix Euclidean norm) while guaranteeing the moment-matching conditions in the obtained parametric model. The optimisation problem of Algorithm 1 is a constrained least-squares problem and can be solved using computationally efficient state-of-the-art solvers as those detailed in, for example, Boyd and Vandenberghe [35].
Remark 13. The model computed with Algorithm 1 has dimension Nν with = f 2 , where N is the number of devices in the WEC array and f is the number of interpolation points (frequencies) selected in the set F . This is a consequence of the fact that, for each frequency ω i and each device, both ± jω i are chosen as eigenvalues of the real-valued matrix S N . Remark 14. The notation F refers to an approximated parametric time-domain model of the force-to-velocity dynamics of the device under analysis in (19), by matching the frequencies selected in the set F using the procedure described in Algorithm 1.

Case study: array of CorPower-like devices
This section presents an application case to illustrate and analyse the proposed strategy, based on the square WEC array layout studied in [5], and depicted herein in Fig. 2. This particular layout is composed of = N 4 energy converters, where each of the four devices composing the WEC farm is a state-of-the-art full-scale CorPower 10 -like device oscillating in heave (translational motion). Note that this type of device is often considered as a case study due to its intrinsic geometrical complexity (see, for example, Giorgi and Ringwood [36]), and is illustrated herein in Fig. 3, along with its corresponding physical dimensions specified in metres. These dimensions are based on the experimental study performed in [37]. From now on, and without any loss of 1 Select a set of f interpolation points (frequencies) F = {ω p } f p=1 ⊂ R + to achieve moment-matching. 2 Compute the matrix I N ⊗ S following (24) and select any L e such that the pair (L e , I N ⊗ S ) is observable. 3 Calculate the moment-domain equivalent of the output of system (19) V using the result of Lemma 5. 4 Compute the matricesÂ Nν andĈ Nν as described in equation (31). 5 Consider the familiy of systems (10) and set F = F ϕ =Â Nν and Q = Q ϕ =Ĉ Nν . 6 Consider the frequency response of (10) as a function of G ϕ i.e.
Using the frequency set Ω = {ω i } M i=1 , compute the input matrix G opt ϕ with the following optimisation-based procedure: 7 Compute a Nν-dimensional WEC array time-domain model achieving moment-matching at S from (10) as  Fig. 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 A(ω), respectively. 11 Note that, since the devices composing the WEC farm are identical (i.e. CorPower-like devices), the corresponding hydrodynamic characteristics (including interactions arising due to radiation effects) present symmetrical behavior, in accordance to the layout depicted in Fig. 2. That said, only three elements of the matrices × B A { ( ), ( )} 4 4 are required to completely characterise the hydrodynamic parameters of the complete array. These are plotted in Fig. 4, along with the corresponding symmetry pattern 12  For both the numerical device motion, and wave excitation force estimation analysis of this section, we consider irregular waves as inputs, with a Spectral Density Function (SDF) computed according to the  Note that there is a one-to-one relation between the colors of the lines and the corresponding symmetry pattern depicted in the figure. 11 For this study, the hydrodynamic characteristics of the WEC array are computed using the BEM software NEMOH [6]. 12 The reader is referred to [ Fig. 5). The total duration of each generated wave and, hence, each simulation, is set to 200 [s]. Finally, without any loss of generality, the incident wave direction β is set to = 0 (see Fig. 2).
As anticipated in Section 5, we consider the proposed momentmatching-based strategy for both WEC motion simulation (in the frequency-and time-domains) and wave excitation force estimation (timedomain). In order to obtain statistically meaningful and consistent results for the time-domain scenarios, and since the waves are generated from sets of random amplitudes [40], it is found that the mean of 15 simulations is necessary to obtain a 95% confidence interval with a half-width of 0.2% of the mean, computed as in [5].
From now on, we denote the frequency-domain model corresponding to the WEC array H(jω) as the target frequency-domain response. In addition, we use the notation H ij (jω) for the ijth element of the matrix H(jω). More precisely, H : ij 0 is the frequency response mapping between the wave excitation force acting on the ith device f , ei and the output velocity of the jth device x j .

WEC array motion simulation
Following the theoretical results developed in Section 4 and the method proposed in Algorithm 1, we now proceed to the explicit computation of a moment-based parametric model F to approximate the force-to-velocity dynamics described in Eq. (19). To fulfill this objective, we make explicit use of the target frequency response operator H(jω) computed using BEM codes. We recall that a key feature of this moment-based strategy is that the user is allowed to select a set of frequencies F to interpolate, i.e. a set of frequency points where the approximating model F exactly matches the steady-state response of in (19). As discussed in [17] for the SISO case (i.e. single device in the array), a sensible selection of the interpolation set can be achieved by analysing the gain of the target frequency response, and selecting points that characterise dynamically important features of the WEC. By way of example, a sensible selection always includes the resonant frequency of the device under analysis. We note that this is the frequency where the maximum amplification occurs, i.e. the frequency characterising the H -norm of the WEC system. For this MIMO case, it is well-known that the system gain depends intrinsically on the corresponding input direction [41], hence the selection of these dynamically relevant points cannot be done by simply inspecting each H ij (jω) individually. As a matter of fact, the gain of a MIMO system is defined in terms of the singular values of H(jω), plotted, for our CorPower-array case, in Fig. 6 (dashed-blue).
Based on the previous discussion, and to test the proposed strategy in the motion (velocity) simulation case, we consider different frequency interpolation points sets F , as follows:  precise, the presented moment-based strategy is able to preserve the H -norm of the target system by simply including the corresponding frequency in the interpolation set. The set F 2 sim additionally includes the frequency point where the second maximum amplification peak occurs. Finally, the sets F 3 sim and F 4 sim expand F 2 sim by including a low-, and a low-and high-frequency component, respectively.
We begin the assesment of the moment-based parameterisation approach in terms of WEC array motion simulation, by illustrating the performance of the approximating models F 1 sim and F 4 sim . Fig. 6 presents the so-called sigma-values plot, for both the target response H (jω) (dashed-blue), and the approximating frequency response mappings (solid-green) F H 1 sim (top) and F H 4 sim (bottom). The interpolation points selected for the computation of each approximating parametric structure are denoted by an empty black circle. Note that, as expected by the theoretical foundations of this moment-based strategy, the approximating models have exactly the same gain as the target model H (jω) of the WEC array, for each element of the corresponding interpolation set F sim . In addition, it can be readily appreciated that, by a sensible selection of the interpolation frequency set F , 1 sim the model F , 1 sim i.e. a parametric description computed using a single interpolation point, already provides an accurate frequency-domain description when compared with the target steady-state response of the WEC array under study. Though considering F 1 sim as an interpolation set provides quite accurate results, the decrease in the overall approximation error from system F 1 sim to F 4 sim can be clearly appreciated (further discussed in Table 1).
As a conclusive graphical illustration of the frequency-domain performance for the models computed via this moment-based strategy, Fig. 7 presents the Bode diagrams for the mappings H 11 (jω), H 12 (jω) and  H 14 (jω). Note that, due to the underlying symmetry of the WEC array illustrated in Fig. 4, these mappings are indeed sufficient to completely characterise the frequency response function H(jω). Analogously to what is presented in Fig. 6, Fig. 7 uses the same colors and lines to characterise the target response, and the approximating models F H 1 sim (top) and F H 4 sim (bottom). Once again, it can be appreciated that, by a sensible selection of the set F , 1 the moment-based approximating model, computed using a single interpolation point, already presents accurate results, though this performance is indeed improved by the set F 4 . Note that both parametric models have exactly the same response (both in magnitude and phase) as the target WEC array H(jω), for each corresponding interpolation set F sim .
To provide a precise measure of the performance of the momentbased parametric approximations of the WEC array considered in this case study, Table 1  We note that the first row of Table 1 includes what we call the "multi-SISO" approach, which essentially constitutes a parametric model of the WEC array model described in (19), obtained by approximating each individual convolution operator (arising due to hydrodynamic interaction between devices) with a SISO system, i.e. in a "decoupled" fashion. As discussed in Section 1, this is the predominant approach in the wave energy literature, and hence is included in Table 1 for the sake of comparison. In this paper, the strategy used to compute an approximation of each of these convolution terms separately, is the frequency-domain method presented in [42]. The dimension of each approximating model is set 14 to 4, which results in a full state-space description of the WEC array of dimension (order) 72.
It is straightforward to acknowledge, from Table 1, that the frequency-domain performance of the moment-matching based models is always over 95% accurate, being able to successfully represent the target WEC array even in the case where a single interpolation frequency is (sensibly) chosen. In fact, it is noteworthy that the momentmatching based model F 1 sim provides almost the same frequency response performance as the "multi-SISO" approach, but with ≈ 44% improvement in computational time (due to the substantial difference in model order). This performance is progressively improved by using the interpolation sets F , 2 sim F 3 sim and F , 4 sim with a mild increase in computational demand. In fact, the moment-based model with interpolation set F 2 sim already provides a better result than those of the "multi-SISO" approach, both in terms of accuracy and computational time. This is indeed directly associated with a sensible selection of the interpolation frequencies, as previously discussed in this same section.

Remark 15.
It is certainly possible to improve the computational requirements of the "multi-SISO" approach, by subsequently applying a model order reduction technique. Naturally, this has some drawbacks: Firstly, the quality of the resulting (reduced) model directly depends on the first parameterisation stage, so that, for the non-parametric WEC case, the main focus should always be put in producing accurate models from the corresponding non-parametric data. Secondly, as expected from any model reduction technique, there is an intrinsic trade-off between improving the computational requirements of the model, and the final approximation quality.
The increase of approximation quality, when considering the different interpolation sets in time-domain simulations, is consistent with the previous frequency-domain results, though it can be appreciated that the approximating model  Finally, to briefly illustrate the transient response of the approximating models computed with the presented strategy, Fig. 8 (15). It can be readily noted that the velocity computed with the target and approximating model perfectly overlap throughout both the transient period (approximate length of 20 [s]), and the steady-state regime. We note that the transient response obtained with the "multi-SISO" approach behaves similarly to that of F , 4 sim and it is excluded from the plot for the sake of visibility.

Wave excitation force estimation
The estimation strategy considered in this paper is that utilised in [5], where the estimation of wave excitation force for a WEC array is discussed. This strategy belongs to the class of unknown input observers, where the system's input (wave excitation force e acting on the device) is estimated using only velocity measurements of the WEC array (i.e. the vector t ( )) based on a direct application of the internal model principle [16]. To this end, a Kalman Filter (KF) [15] is used, in conjunction with a harmonic oscillator model, to describe the dynamics of the excitation force. In fact, we note that the dynamical model used to describe the excitation force is exactly defined by the signal generator proposed in (23). In other words, there is a natural synergy between the moment-matching-based parameterisation strategy proposed in this paper, and the unknown input estimation problem inherently present in wave energy systems. 15 This intrinsic connection, together with a summary of the unknown input estimation strategy presented in [5] (for the sake of completeness), are discussed in the subsequent paragraph.
With regard to the unknown input observer design, we begin by defining the so-called augmented system (see [16]), given by the following continuous-time system, described, for , in state-space form as where ϵ(t) and ζ(t) represent the process and measurement (white) noises, with associated covariance matrices Q and R , respectively. The extended matrices, and state-vector, involved in (34), are defined as where t ( ) N 2 contains the system and signal generator state-vectors, i.e. Ξ e (t) and Θ φ (t), respectively. Following [5], we define an optimal continuous-time KF observer as where the optimal gain K K + t t : , ( ), can be computed [15] as, Finally, the estimated wave excitation force can be directly computed in terms of the estimated state-vector as As can be directly appreciated from Eq. (35), F a explicitly includes the dynamic matrix of the signal generator (23), which is utilised in the moment-based parameterisation method discussed in Section 5. We now show that this mathematical connection can be exploited to improve the performance of the unknown input observer substantially, by a wise selection of the set of interpolation points F . To be precise, if we compute a parametric representation by moment-matching the inputoutput dynamics of the WEC array F , interpolating the same set of frequencies used to describe the internal model of the wave excitation force in the observer defined in (36) (i.e. the set λ(S)), the performance of the unknown input estimator improves when compared to the case where the currently most-used parameterisation approach (i.e. "multi-SISO") is considered to describe the WEC array motion. This is associated with the fact that the model based on the moment-matching strategy presented herein has exactly the same steady-state response as the target non-parametric model of the WEC array computed with BEM solvers at the key frequencies utilised to describe the excitation input. We note that this improvement in performance is given both in terms of estimation quality and computational effort required by the observer, as demonstrated in the remainder of this section.
To illustrate the advantages of using the moment-based parameterisation method proposed in this paper within the unknown input estimation problem, we design a KF observer with a matrix S completely characterised by the sets (see Eq. (24) 16 Using the frequency sets defined in Eq. (39), we compute the moment-based statespace models F 1 est and F 2 est following the procedure described in Algorithm 1. Considering these two moment-based models, specifically designed to correlate with the optimal observer of Eq. (36), and the momentmatching models F 1 sim and F , 2 sim computed for the WEC array motion simulation case of Section 6.1, we design a KF using each of the parametric representations listed above, and assess the performance of the different unknown input observers in terms of the following indicators. Note that, analogously to Section 6.1, we also include in the upcoming discussion a KF designed using the state-space model arising from the "multi-SISO" model.  Fig. 5). T-Gain % improvement in normalised run-time (i.e. the ratio between the time required to compute the estimated wave excitation force, and the length of the simulation itself) with respect to the slowest model (normalised run-time indicated in table between parenthesis).
We begin the analysis of Table 2 by making an explicit performance comparison between the moment-matching-based models F 1 sim and F , 1 est for the unknown input estimation problem. Note that the dimension (order) of both models is exactly the same, i.e. the same number of interpolation points are used to compute both parametric representations. Despite sharing the same model complexity, the performance results are significantly different, as discussed in the following. Whilst F 1 sim provides a much better overall approximation in frequency-domain (with almost 96% of accuracy) than F 1 est , the performance of the estimator designed using F 1 sim is quite poor, with a NRMSE e T of ≈ 24%. Instead, having a frequency-domain NRMSE F of almost 32%, the KF designed used F 1 est performs substantially better in terms of wave excitation force estimation quality, with less than 10% of estimation error. This can be attributed to the fact that F 1 est has exactly the same steady-state response as the WEC array at the frequency points selected to represent the internal model of the excitation force embedded in the optimal observer (36), where the majority of the spectral content of the input signal is contained (see Fig. 5). Fig. 9 graphically illustrates the performance of the KF designed using , where it can be appreciated that both time traces, i.e. target (dashed-blue) and estimated (solid-red) wave excitation forces, for each of the four devices composing the analysed WEC array, are qualitatively identical. Fig. 10 shows a Bode diagram analogous to that of Fig. 7, where the frequency-domain description and performance of the parametric models  Table 2, the overall frequency-domain fit for the selected frequency range of It is noteworthy to highlight that exploiting this inherent relation between moment-matching parameterisation and the internal model principle, used to estimate the wave excitation force, using a single interpolation point F , 1 est provides similar estimation accuracy results as the case where the "multi-SISO" approach is considered, but with ≈ 81% improvement in computational requirements. In other words, the estimator using the moment-based model F 1 est computes 81% faster than the currently most-used method, for the same degree of wave excitation force estimation accuracy, and is therefore especially suited for real-time applications.
Finally, Table 2 also provides results for the KF observers designed using the moment-based models F 2 sim and F , 2 est where the perfomance (frequency-domain fitting and wave excitation force estimation quality) is subsequently improved, for both cases. Note that the situation described in the previous paragraph is repeated, i.e. the KF designed using F 2 est provides better estimation performance due to the particular selection of the interpolation points (same as those carefully selected to describe the internal model of the input in (36)).

Conclusions
In this paper, we present a MIMO moment-based parameterisation strategy for wave energy applications, making specific emphasis in motion simulation and wave excitation force estimation. This approach allows for the computation of state-space representation characterising the input-output dynamics of WEC arrays which exactly matches the target steady-state behavior of the array at a set of user-selected frequencies. We explicitly discuss how the selection of these interpolation points have a key role in the performance of the computed models, both for WEC array motion simulation and wave excitation force estimation applications. In fact, we show that there exists an intrinsic mathematical relation between the unknown input wave excitation force estimation problem and the moment-based strategy presented herein, which we exploit by a sensible selection of the set of interpolation points, in synergy with the internal model principle utilised to estimate e . Finally, we analyse the performance of the strategy using a state-ofthe-art WEC array composed of CorPower-like devices, showing that we can outperform the "multi-SISO approach" (i.e. currently used method) both in terms of accuracy and computational requirements, hence providing models that are especially suited to design real-time energymaximisiation strategies, contributing to the roadmap towards successful commercialisation of WEC devices.