Robust optimization of control parameters for WEC arrays using stochastic methods

This work presents a new computational optimization framework for the robust control of parks of Wave Energy Converters (WEC) in irregular waves. The power of WEC parks is maximized with respect to the individual control damping and stiffness coefficients of each device. The results are robust with respect to the incident wave direction, which is treated as a random variable. Hydrodynamic properties are computed using the linear potential model, and the dynamics of the system is computed in the frequency domain. A slamming constraint is enforced to ensure that the results are physically realistic. We show that the stochastic optimization problem is well posed. Two optimization approaches for dealing with stochasticity are then considered: stochastic approximation and sample average approximation. The outcomes of the above mentioned methods in terms of accuracy and computational time are presented. The results of the optimization for complex and realistic array configurations of possible engineering interest are then discussed. Results of extensive numerical experiments demonstrate the efficiency of the proposed computational framework.


Introduction
Ocean waves and tides are a renewable source of energy with large predicted growth over the next years and decades. Proposals of the European Commission [2] regard as "realistic and achievable" a target of 1 GW installed capacity by 2030 and of 40 GW by 2050 for ocean energy. 1 A massive deployment of wave energy is expected to drive costs down and make this choice competitive, especially for islands and offshore applications [1]. Since single devices for wave energy conversion have small production capabilities and a limited possibility of scaling up compared to other technologies such as wind turbines [6], it is crucial from the economical point of view to install them in arrays. The relatively small distance between devices in arrays implies the appearance of mutual interactions. The result of these interactions is known as the park effect [5], and it has to be carefully taken into account in the design process.
In order to better utilize the energy resource and improve survivability in extreme conditions, control strategies of devices in the array need to be appropriately designed. In particular, hydrodynamic interactions in the controlled system should be favourable in terms of extracted power. The most common control techniques for wave energy converters are presented in [20], reviewing constant and time-varying damping control, reactive control, latching and model predictive control. Classical results obtained through frequency domain analysis for regular waves with control parameters kept constant in time are presented in [11]. Strategies for irregular waves (in time domain) and related results are discussed in [17]. In an array, the same control parameters may be assigned to all devices, or each device may be assigned its own. The latter approach is referred to as individual optimization and it has been performed in [7] and [28] using SQP (Sequential Quadratic Programming), optimizing for a single wave direction and then checking the performance for different directions. Another aspect of array design and optimization, that will not be treated in this work, is the layout of devices. Layout optimization can be performed either before or together with control optimization. For more details, see the reviews [33,15].
In this work, we present a numerical strategy for the robust optimization (maximization) of the power production of arrays with reactive control: damping and stiffness coefficients of the generators are optimized. We allow each body to have its own control coefficients and assume that the corresponding parameters are constant in time. As shown in [10], damping is actually a nonlinear function of the velocity for real electrical generators. Our assumptions, though, allow the use of a linear framework, and the results can still be a guide for the design process. We also allow control stiffnesses to assume negative values, meaning that the control action must be able to counteract the hydrostatic restoring force, that is the force acting on a floating body when it is displaced from its buoyancy equilibrium position [11,Sec. 5.9]. Even though a negative stiffness could seem counterintuitive or unphysical, some technologies are in fact capable of achieving this effect: see, e.g., [30,34]. For this reason, we do not exclude apriori this case from the set of admissible solutions. Moreover, the optimal solutions must fulfill the slamming constraint, meaning that devices must not leave the water, and in some cases a constraint of positive stiffness is added, to assess the benefit of using negative stiffness technologies.
To obtain good performances in real scenarios, the result of optimization needs to be robust against the uncertainty in the incident waves direction. To the best of our knowledge, this problem has never been addressed in the literature. We formulate a stochastic (robust) optimal control problem governed by linear hydrodynamic equations, describing the physics of the system, and characterized by control and state constraints. The numerical solution of this optimization problem requires an adequate treatment of different issues. A reduced approach based on the adjoint equation allows us to work in the space of solutions of the hydrodynamic equations and consider the control as the only optimization variable. Stochasticity is addressed by employing two different approaches suitable for stochastic optimization problems: SAA (Sample Average Approximation) and SA (Stochastic Approximation) [27]. The slamming (state) constraint is treated by using a penalty approach, while the stiffness (control) constraint is enforced by projection. Therefore, our proposed numerical computation framework is obtained by two nested iterative processes. The outer iteration is essentially a quadratic penalty method increasing at each iteration the weight of the penalty term for the slamming constraint. At each penalty iteration a stochastic optimal control problem is solved by a projection method based either on SA or SAA. The performances obtained by individually optimizing the coefficients of each device are compared with those obtained by setting the coefficients of all devices equal to the tuning parameters of the isolated case. We show that sets of parameters corresponding to feasible solutions for isolated devices can lead to unfeasible (incompatible with the constraints) solutions when the devices are installed in arrays, and that the optimization process leads to sensible increases of the power output in most cases.
The paper is organized as follows. In Section 2, we detail the hydrodynamic and mechanical models used and discuss the constraints. In Section 3, we study the wellposedness of the problem and present optimization methods and algorithms. Sections 4 and 5 contain numerical experiments performed with the proposed algorithms. In Section 6 the obtained results are commented and conclusions are drawn.

Hydrodynamic modeling
To model the hydrodynamic behaviour of WEC (Wave Energy Converter) arrays, we use linear potential wave theory [9]. In particular, we make the kinematic assumption of irrotational flow (zero vorticity), which implies the existence of a scalar potential φ whose gradient is the flow velocity. This assumption is well verified for external flows, outside boundary layers [18]. Under the assumption that the displacements of the free surface and of the surfaces of floating bodies are small, boundary conditions can be linearized and enforced at the reference surfaces. This avoids the difficulty of dealing with moving boundaries and permits solving the problem in the frequency domain. Finally, we assume that viscous effects can be neglected.
For a single frequency ω, the potential is written as φ( andφ is the solution of the hydrodynamic problem Figure 1: Domain of the hydrodynamic problem (slice) where Ω ⊂ R 3 is the fluid domain, bounded by the sea bottom surface Γ b , the reference free surface Γ s and the device immersed surfaces Γ d, , and unbounded on the lateral sides (see Fig. 1). The complex amplitude of the velocity of the -th device is denoted byv . We assume that the bodies are restricted to vertical motion (heave). Hence, we havev whereζ is the complex amplitude of the vertical position of the -th device. The free surface elevationη can be computed using the kinematic conditionη = iφω/g. For the -th device, the dynamic problem in the z direction is where m is the mass, k is the combined hydrostatic and mechanical stiffness, c and s are the equivalent damping and stiffness coefficients of the electrical generator, respectively, andF is the vertical component of the hydrodynamic force. The latter is obtained from the potentialφ:F n z being the vertical component of the outer surface normal of the body. This is a consequence of the unsteady, linearized Bernoulli equation, which holds thanks to the assumption of inviscid flow.
The equations above show that there is a two-way coupling between the hydrodynamic problem and the mechanical problem. The two may be solved monolithically. However, for the purpose of control design, it is more common and insightful to define intermediate hydrodynamic quantities, namely added mass, radiation damping and excitation force [11]. In this way, once the hydrodynamic problem has been solved to determine such quantities, the dynamic problem can be treated independently. To apply this approach, the potential is split into the sum of three harmonic functions: an incident potentialφ 0 , a diffraction potentialφ d and a radiation potentialφ r . Given a wave of direction θ and height H, the complex amplitude of the incident potentialφ 0 is given by [9,Sec. 3.5]: This function satisfies the Laplace equation in Ω and the boundary conditions on the bottom and the free surface, but not the boundary condition on device surfaces. We can then rewrite (2.1d) as Such condition can be satisfied by requiring that the diffraction and radiation potentials satisfy The radiation potential φ r can then be decomposed as a sum of the effects of radiation from each body,φ so that each ϕ satisfies a condition of unit velocity amplitude on the -th body and an homogeneous condition on all the others. Diffraction and radiation potentials are also required to satisfy a radiation condition at infinity, which is needed for energy conservation. The full potentialφ =φ 0 +φ d +φ r is then computed by solving a single diffraction problem to obtainφ d and a radiation problem for each body to get ϕ , (2.5) Using the decomposition of the potential, equation (2.3) can be rewritten aŝ where the first term is the excitation forceF e, , and the m-th term of the summation is the radiation impedance R m . This is further split into real (in phase with velocity) and imaginary (out of phase) parts: where A m is called added mass and B m radiation damping. Summarising, using the above decomposition of the potential together with (2.6) and (2.7), the dynamical system (2.2) can be recast as a complex linear system for each frequency: where the impedance matrix Z(ω q ) is defined as and δ m is the Kronecker symbol. The dependence on frequency ω q and on direction θ has been made explicit. Notice that only excitation forces depend on the wave direction. The only terms coupling the motions of different devices are the added mass and radiation damping coefficients.

Modeling of constraints
The model described in Section 2.1 needs to be complemented with constraints to guarantee that the solutions are physical and that the control forces are feasible. A description of some possible constraints is reported in [7], namely slamming, stroke and force constraints. The slamming constraint guarantees that the bodies do not leave the water, and it is needed because the linear potential model is built on the assumption of small displacements. The stroke constraint limits the amplitude of body motions; it is related to the dimensions of the power take-off system. Finally, the force constraint ensures that the control force can be realized by an electromechanical system. In this work, since we aim for generality, we focus on the slamming constraint, that is a purely hydrodynamic one. Conversely, force and stroke constraints depend strongly on the design details of a specific power take-off system. Once a system is specified, the other two constraints can be imposed with the same methodologies presented here.
The slamming constraint is a state constraint expressed in time domain as where d is the draft of the -th device and η is the wave height at the center of the -th device. Let us now formulate this constraint in frequency domain. First, we consider the approximation introduced in [7], which consists in neglecting diffraction and radiation effects and taking the wave height from the incident wave whose complex amplitude iŝ For monochromatic waves, constraint (2.9) implies ζ −η 0 ≤ d . (2.10) For irregular waves, it is not possible to obtain an explicit formula such as (2.10). The constraint can instead be imposed in a statistical sense by limiting the fraction of time above threshold t at , meaning the fraction of time for which the constraint is violated, or the fraction of peaks above threshold Q, where constraint violations are considered as discrete events. It has been observed that for deep sea waves the free surface elevation η(t) at a fixed point is well described as a Gaussian process [26]. Moreover, if the input of a linear time invariant system is a Gaussian process, then its output is also a Gaussian process. We thus work under the approximation that w(t) is a Gaussian process, which implies that the constraint can be enforced on the root mean square value w rms [19].
Such quantity can be computed from frequency domain data: The constraint is then written for each body as w ,rms ≤ αd , (2.12) where α is a constant chosen to control the exceedance probability. In particular, the fraction of time above threshold is where Φ(·) is the normal cumulative distribution function. The fraction of peaks above threshold is, for a narrow banded process [19], (2.13) The process we consider does not have a narrow band. However, it can be shown that, for a generic Gaussian process, (2.13) is a cautionary estimate (upper bound) of the exceedance probability. As an example, setting α = 1/2 and assuming that the constraint (2.12) is active, i.e. w ,rms = αd , yields t at, = 4.55% and Q p, = 13.5%. For convenience, we rewrite (2.12) using (2.11) as The two constraints are equivalent in the sense that they have the same active set: In addition to the slamming constraint, a constraint of positive control stiffness is enforced in some of the simulations to quantify the advantage of using negative stiffness technologies.

Forces stochasticity modeling
As noted in section 2.1, the wave climate influences the system dynamics through the excitation forceF e . Since our aim is to obtain control parameters that are robust with respect to the wave direction, we define a realization as a condition with irregular unidirectional waves, where the direction is sampled from a probability distribution. By irregular we mean that waves are not monochromatic, but rather a superposition of harmonics. Frequency and amplitude of each harmonic are the same for all realizations. The wave climate is described by the directional spectrum S(f, θ), which is the distribution of wave power with respect to frequency f and direction θ. Its general form is [13] S(f, θ) = S(f )D(θ|f ), with short (high-frequency) waves generally showing a larger directional spreading than long waves. In this work, we consider the simpler approximate form S(f, θ) = S(f )D(θ).
For the frequency spectrum we use the two-parameter Piersov-Moskovitz form [11], whose expression is and where H s is the significant wave height and f p is the peak frequency. The spectrum is discretized using the deterministic spectral amplitude method [9] by first defining N f intervals of center f q and width ∆f q . Intervals are chosen so that they all correspond to the same power fraction (meaning they have in general different widths). The incident wave potential is then a superposition of potentials of the form (2.4), where the q-th component has height H q = 2 2S(f q )∆f q . If necessary, the wave profile in time domain can be reconstructed as with φ i ∼ U(0, 2π). In the limit of N f → ∞, η(x, t) is a Gaussian process [32]. For the directional part, choosing a spectrum that has a closed-form inverse cumulative distribution makes sampling straightforward. A suitable choice is Donelan's spreading function [12]: with the corresponding cumulative distribution functioñ where θ 0 is the dominating direction and β is a scale parameter: for increasing β, the distribution is more sharply peaked (see Fig. 2). To sample a wave direction, we use the inverse transformation technique [16, Sec. 9.2]: we first obtain a sample ξ from a uniformly distributed variable X ∼ U(0, 1), and then apply the transformation θ =D −1 (ξ).

Computational aspects and Haskind's relation
The hydrodynamic problem is solved with the boundary element method (BEM), based on boundary integral reformulations of (2.5). In particular, the source distribution method is used [22,Sec. 4.11]: an unknown source distribution σ is defined on the device surfaces Γ d, , = 1, . . . , N b . The potential in any point x of the domain is obtained as the effect of such surface distribution through a convolution integral with a Green function G, which satisfies the boundary conditions on ∂Ω \ ∪ Γ d, : (2.16) The boundary conditions on the body surface are then enforced: where the b.c. depends on the considered problem in (2.5) The two equations (2.16), (2.17) are transformed into linear algebraic equations by discretization: Once the first system is solved, the potential is obtained from the second. We now discuss the computational effort required to construct a realization of the excitation force. In principle, one would need to draw a sample θ and then, for all frequencies, compute the corresponding incident wave potential using (2.4) and use it as a boundary condition for a diffraction problem. This requires storing or recomputing N f dense, and possibly large, matrices and solving each system once. This would make the process very computationally demanding. A much cheaper approach would be to run a batch of simulations for different values of the wave direction before the optimization process and then recover the force by interpolation for each realization. This strategy, however, could result in large inaccuracies because of the difficulty and ambiguity of interpolating complex functions [29]. Another possibility is the use of Haskind's relation, an exact integral relation which allows computing the diffraction force for any wave direction as a combination of the potential of incident waves and radiation potentials [11,Sec. 5.4]. The recomputation of the incident wave field is relatively cheap, while radiation potentials would need to be computed anyway to obtain the added mass and radiation damping matrices. For these reasons, this approach is intermediate in computational cost between full diffraction simulations and interpolation. The use of Haskind's relation is the choice that has been made here. Its expression iŝ The integral is discretized as a sum on the mesh elements.

Design of damping and stiffness by optimization
In this section, we formalize the optimization problem and present numerical methods for its solution.

Problem statement
The cost function is the average power with negative sign, which can be written as The motion of each body is realized in the time domain as a sum of harmonics. Because of the orthogonality of such functions, the mean square of the signal depends only on the amplitudes of the harmonic components, phases being irrelevant. If incident waves are realized as the sum of N f monochromatic waves, as explained in Section 2.3, the cost is obtained by summing the contributions of all frequency components: is a diagonal matrix with C = c , and we denote by the superscript θ all quantities depending on the wave angle. Here,ζ θ q ∈ C N b is the vector of complex amplitudes of all bodies at the q-th frequency, whileζ θ ∈ C N b N f is the collection of allζ θ q for q = 1, . . . , N f . For a single wave angle, the problem of power maximization can then be stated as The admissible set of controls U ad is defined by where ε is a positive constant, and γ = 0 if the positive stiffness constraint is applied, γ = −∞ otherwise. Next, we discuss sufficient conditions for Z q (u) to be invertible.
Splitting into real and imaginary parts and solving formally for x r and x i yields Since Z is invertible, the expressions in (3.3) are well defined for any f and thus Z is invertible.
Consider now the impedance matrix defined in (2.8). Matrices A and B are symmetric, and because of energy conservation B is also positive semidefinite [11, Sec. 5.2, 6.5]. The other involved matrices are diagonal and thus symmetric. Then, the real part of the impedance matrix Z r = −ω 2 (M + A) + K + S is symmetric. Regarding the imaginary part Z i = −ω(B + C), we have that B is symmetric positive semidefinite. If we further assume that all bodies have positive generator damping, implying that the mean extracted power is positive, then C is diagonal with positive values. Thus Z i is symmetric positive definite and Z r is symmetric, satisfying the hypotheses of Lemma 3.1.
Since the matrix Z q (u) is invertible for any admissible control vector u, we can formally eliminate the constraint Z q (u)ζ θ =F θ and rewrite (3.2) in the reduced form The following properties hold for the case of a single wave direction θ.
Proof. Lemma 3.1 ensures that Z q (u) is invertible for all u ∈ U ad . Then,F θ q = 0 impliesζ θ q = 0. From this and using again the assumption that u ∈ U ad , so that in particular c l > 0, it follows that ζ θ q H Cζ θ q > 0. From the definition of J (3.1) we have The following asymptotics hold for u 2 → ∞ and ∀q, ∀θ: In particular, lim u 2 →∞ J(u) = 0.
Proof. From the definition (2.8) of the impedance matrix, one directly finds The equivalence between the Frobenius norm and the 2- For the second statement, consider where B 0 ( C) is a ball of radius C centered in the origin.
Proof. Lemmas 3.2 and 3.3 imply that In particular, if J has a minimum over such set, then the latter is also a minimum over U ad . Proof. Invertibility of the impedance matrix Z q (u) guarantees that the reduced cost J is a continuous function of u. Lemma 3.4 implies that the problem can be equivalently posed on the compact set U ad ∩ B 0 C . The existence of a minimum in such set is guaranteed by the Weierstrass theorem.
is a compact set and h is continuous for all , h −1 (0) is compact. Then, U ad is also compact and the same reasonings of Theorem 3.1 can be applied.
In the stochastic case, (3.4) turns into where E θ denotes the expected value with respect to θ.
admits a solution for all C > 0.
Proof. Lemma 3.2 guarantees that J(u, θ) ≤ 0 ∀θ, ∀u ∈ U ad . Then, we have where the last equality follows from Lemma 3.3. Then, the problem can be equivalently recast in a compact set U ad ∩ B 0 C , for some suitable positive constant C. Since such set is compact and J is a continuous function of both its arguments, the Heine-Cantor theorem implies J is also uniformly continuous. We now use uniform continuity and the fact that D(θ) is a probability distribution, so that it has unit integral. We have where δ does not depend on θ. Then, if u − v < δ, the following inequality proves that j(u) is also continuous. Then, existence of a solution for the unconstrained problem follows from the same reasonings applied in the proofs of Lemma 3.4 and Theorem 3.1.
For the constrained problem, we can prove, as done above for j, that E θ [h (u)] is continuous, so its inverse maps {0} into a compact set. Then, the remainder of the proof is the same as for Theorem 3.2.

Numerical optimization strategy
In this section, we present a computational framework for the solution of the stochastic optimization problem (3.2) modelled in Sections 2 and 3. Problem (3.2) has three different constraints that must be appropriately treated: a system of complex algebraic equations (the state problem), the slamming constraint (a constraint on the state) and the positive stiffness constraint (a constraint on the controls). First, the slamming (state) constraint is treated by using a penalty method, as described in Section 3.2.1. Our framework is based on a reduced approach working in the space of solutions to the state equation and considers the control u as the only optimization variable. To do so, the gradient of the (reduced) cost j (and its penalized counterpart Q µ ) is obtained by a Lagrangian approach leading to an adjoint problem, as show in Section 3.2.2. Each iteration of the penalty method requires the solution of a stochastic problem with control constraint (positive stiffness constraint). This is achieved by either an SAA approach (Section 3.2.3), or by a Robust SA (Section 3.2.4). In both cases, the control constraint is enforced by projection onto the admissible set at each iteration. Implementation details and a comparison between SAA and Robust SA are given in Section 3.2.5.

Penalty method
The penalized cost function used to enforce the slamming constraint is defined as and the penalty method consists in solving the sequence of problems min u∈U ad Each of these problems is solved using either sample average approximation or robust stochastic approximation; these methods are described in Sections 3.2.3 and 3.2.4, respectively. The penalized cost function can be rewritten as an expected value thanks to linearity, and, again because of linearity, also the gradient of Q µ may be computed as the expected value of the quantity in brackets. The optimization procedure is detailed in Algorithm 1.
Algorithm 1 Stochastic quadratic penalty algorithm Compute u j+1 ∈ U ad by minimizing Q j+1 , using either SA or SAA with initial guess u j and tolerance τ j 5: Convergence results for Algorithm 1 can be obtained from [23,Th. 17.1,17.2]. It is proved that, if at each iteration the exact minimum of Q is found and if µ j → ∞, then every limit point of the sequence u j is a global solution of the original problem. It is further proved that if at each iteration the solution of the optimization problem is computed approximately with vanishing tolerance k τ , as in Algorithm 1, either a limit point of the sequence u j is infeasible and corresponds to a stationary point of the penalty term, or it is a KKT point for the original problem. The above results require differentiability of the penalty term, which is satisfied in our case, as ∇h = ∇[g ] 2 + = 2[g ] + ∇g . Since convergence is only guaranteed in the limit, the method is terminated when constraint violations are below a finite tolerance, that is when the obtained iterate is in, or close enough to, the feasible set.

Gradient computation and first-order optimality condition
Algorithm 1 requires finding the minimum of the penalized cost function Q µ . In order to do this, the computation of its gradient is needed. The SA method needs a stochastic gradient at each iterate, that is, a gradient computed for a single realization of the wave angle θ. The SAA method, instead, needs an estimate of the true gradient, which is the expected value of the stochastic gradient and thus can be computed as a suitable average of stochastic gradients. For a single given angle θ, the gradient of the reduced cost can be obtained by using the Lagrangian function defined as The problem is set in the space of complex vectors over the field of real numbers, with scalar product (u, v) = Re[u H v], thus the presence of the real part in the Lagrangian. The adjoint equation is obtained by differentiation with respect to the state vector, that is independent of the frequency. The components of the stochastic gradient are obtained by differentiating with respect to the controls: The first-order optimality condition requires [23] u where P U ad is the projector onto U ad : (3.7)

SAA
The sample average approximation (SAA) method is based on writing the expected value in Problem (3.5) as an integral, and approximating such integral with a suitable quadrature rule: The same approximation is applied to the gradient. One then obtains a deterministic problem that can be solved using standard optimization algorithms. In this work, we use the gradient descent method with step length determined according to the Armijo rule [23,Chap. 3]. The quadrature rule is defined by the choice of nodes θ i and weights w i . In the following, we consider the Monte Carlo [27] and Gauss-Legendre [25,31] methods. For the Monte Carlo method, the nodes θ i are obtained by sampling the wave direction, following the procedure described in Section 2.3, and the choice of weights corresponds to averaging, w i = 1/N . In the case of Gauss-Legendre, instead, nodes are defined as the zeros of Legendre polynomials. Numerical libraries provide nodes and weights θ i , w i on the interval [−1, 1], which are then mapped to the desired interval [a, b] using In our case, the interval [a, b] is chosen as [θ min , θ max ]. Here, θ min and θ max define an effective interval of wave angles, which depends on the parameters β, θ 0 , such that the probability is considered negligible outside such interval. Notice that each weight for Gauss-Legendre integration includes the value of the probability distribution function in the corresponding node. This is not the case for Monte Carlo, where the probability distribution is accounted for in the sampling phase.  The order of convergence of Gauss-Legendre, instead, is estimated at a value of 2.53, while repeating the test for the unconstrained problem leads to a spectral convergence of the Gauss-Legendre method, as shown in Fig. 4 (right). In particular, we obtain a behaviour of the type k exp(−mN ), where m ≈ 0.658. The loss of spectral convergence in the constrained case is explained by the presence of the max operator in the slamming constraint; such operator is present in the right hand side of the adjoint equation.

Robust SA
The stochastic approximation method (SA) moves, at each iteration, in the direction of a gradient computed from a single sample of the random variable. We call this vector the stochastic gradient. Since there is no guarantee that the stochastic gradient defines a descent direction, the line-search stepsize strategies suitable for deterministic problems cannot be applied. Instead, the sequence of stepsizes is fixed once the initial step is chosen. A suitable initial step t 0 can be estimated by performing a line search using a single sample of the stochastic gradient.
Some stepsize strategies for SA method are discussed in [21]. We choose the robust variable stepsize strategy, that is, the sequence of steps is A convergence estimate is available, based on the assumptions that that Q is convex and that the admissible set U ad is convex where N is the number of iterations, c Q is a constant independent on N andũ N K is a weighted average defined asũ Such average is taken as the solution of the optimization problem. As a stopping indicator, we use the norm of an average of the stochastic gradients on the most recent iterations: where W is the window size. This is a modification of a criterion proposed in [24] and it can be interpreted as an estimate of the gradient of the corresponding deterministic problem. Such quantity appears in the convergence proof of the penalty method mentioned above, suggesting that our choice is reasonable.

Algorithm 2 Robust Stochastic Approximation (variable stepsize)
1: Set initial step t, window size W 2: while k ≤ it out max or ind > τ in do Findζ q by solving the state problem with (u k ,F θ q )

7:
Find y q by solving the adjoint problem with (u k ,F θ q ) 8: end for 9: Compute the stochastic gradient G from (3.6)

Implementation and comparison between SAA and robust SA
Hydrodynamic data is computed using the code Capytaine [4], version 1.3. The Haskind relation is applied to compute diffraction properties. Such routine and the numerical optimization codes are written in Python, using Numba for performance improvement.
The performances of stochastic approximation and sample average approximation using Monte Carlo and Gauss-Legendre integration are compared in Fig. 5. A constrained case with 15 bodies and fixed penalization parameter has been considered, using the same initial guess and initial step. The reference solution was obtained by a computation with Gauss-Legendre using 320 nodes. The relative error u−u ref / u ref is plot against execution time and number of iterations. The test is run on an average laptop. We can observe that the SAA methods require a much smaller number of iterations with respect to SA, but, since for the former many problems need to be solved at each iteration, the methods are comparable in terms of computational time. In terms of number of iterations, the histories of Monte Carlo and Gauss-Legendre are almost identical up to a certain value of the error, after which Gauss-Legendre exhibits better accuracy. However, in terms of computational times, the histories of the two methods are similar only when the same number of nodes is used. The error of stochastic approximation, instead, initially decreases at a rate comparable to (in terms of number of iterations) or faster than (in terms of computational time) the other methods, and then starts oscillating, as expected from an SA approach (see, e.g., [21] and references therein). Moreover, the results of Fig. 5 provide useful insights into the choice of the method to use in a design phase. If one seeks indicative results in relatively short times, then an SA approach is indeed a reasonable choice: on the one hand, it converges faster in the initial phase when the starting point is far from the optimal solution; on the other hand, it does not require the choice of a number of nodes or of an effective interval of wave angles where the probability distribution is correctly represented. Instead, if one is interested in highly accurate results, then Gauss-Legendre is the method to choose. However, this strategy requires an appropriate choice of an effective interval of wave angles, which is not needed by a Monte Carlo approach. Moreover, if the number of stochastic parameters becomes very large, then the Gauss-Legendre method could become unfeasible, and an SA strategy is more appropriate.

Numerical experiments
The numerical framework presented in the previous sections was used to compute the optimal control parameters of arrays of two different models of WECs (see Table 1). In particular, the results presented in this section are obtained by simulations run using SA. The data for the first case is taken from [14]. The second case is derived from the first by rescaling the generator mass proportionally with respect to the volume. In the first case, the resonance frequency of the uncontrolled body is higher than the peak spectrum frequency. The second case is tuned in such a way that the uncontrolled body frequency is lower than the peak spectrum frequency.
In both cases, the region of the spectrum containing 99.9% of the power is discretized into 30 bins of equal power, so that each realization is obtained as a superposition of 30 harmonics.
To quantify the performance of the obtained configurations, we define the interaction factor q as the ratio between the total average power produced by the array and the sum of the powers of the isolated devices:

Single-body optimization
We first consider optimization for a single device. For case 1, if negative stiffnesses are allowed, the optimal power is 7636 W, achieved by setting the control parameters to c = 31820.7 N · s/m, s = −27022.2 N/m. The outer convergence history is shown in Figure 6 (left), superimposed on a contour map of the absorbed power. The region above the red dashed line is feasible with respect to the slamming constraint, and the black dot is the constrained optimum obtained by exhaustive search over a rectangular grid. If instead only positive values of stiffness are accepted, the optimal power is 5886 W, and the corresponding controls are c = 61774.8 N · s/m, s = 0. Figure 6 (right) shows the statistical meaning of the slamming constraint, which has been formulated in Section 2.2 in the frequency domain. The fraction of peaks exceeding the bound can be controlled by varying α in (2.12). For case 2, the optimal power is 5712 W, obtained with controls c = 733196 N · s/m, s = 1411344 N/m, without imposing constraints on the sign of s. The results of this section are used as initial guesses for the multi-body optimization runs: all devices are initially assigned the same control parameters. Table 2 shows the results of numerical experiments with 2 bodies of type 1. Configuration θ 0 = 0°corresponds to bodies aligned along the mean wave vector; θ 0 = 90°c orresponds to bodies aligned perpendicular to the mean wave vector. In both cases, the resulting interaction factor q is larger when the spreading in wave direction is small, that is when β is large (see (2.15)). In most cases, the slamming constraint is violated for the initial guess. This results in a reduction of power between the initial guess and the optimized result, the latter being compatible with the constraint up to the tolerance. The results for control parameters in the symmetric case θ 0 = 90°are not exactly symmetric: this is due to the use of a stochastic optimization method. The asymmetry may be reduced by increasing the window size used for averaging the iterates. Figure 7 reports the convergence history for test 1. The two jumps correspond to a switch to a higher value of the penalization parameter. After a jump, the cost starts decreasing again, but stops at a larger value than the one before the jump, which would now be unattainable because of the stronger enforcement of the constraint.

Multi-body optimization: devices on two concentric arcs
We now consider an array of 15 devices in an arrangement with 2 concentric arcs, derived from the ones considered in [14], for case 1. The geometry of the array and the results without constraints on the sign of the stiffness are shown in Fig. 8 and 9. For all the following tests, the dominant wave direction is zero: the incident wave vector is aligned with the x axis. We observe that the bodies with the smallest stiffnesses are the ones located in the downwave arc, at the center. They are thus the bodies with the lowest β θ 0 Initial violation ∆P (%)   tuning frequency, and they receive the largest power increase. Conversely, the devices at the outer edges of the array have the largest stiffnesses and they are subject to a moderate power decrease. Regarding dampings, bodies in the upwave arc have larger values than the ones in the downwave arc, corresponding to wider response peaks: the bodies at the front, roughly speaking, are tuned so that they are able to absorbe power in a wider range of frequencies. The global result is a power increase of 3.44% compared  Table 2) to the initial guess, and an interaction factor of 0.966. If the constraint of positive stiffnesses is enforced, the results of Fig. 10 are obtained. Only a very small power increase of 0.31% is obtained, and an interaction factor of 0.980. Contrary to the previous case, bodies in the downwave arc have larger control dampings than the ones in the upwave arc. All optimal control stiffnesses are zero; this is explained by the large, positive frequency distance between the response peaks of the uncontrolled bodies and the frequencies where most of wave power is distributed.
For the same kind of geometry, with bodies of type 2, the results are reported in Figures 11 and 12. All optimal stiffnesses are positive, without imposing any constraint. As in the corresponding arrangement for case 1, bodies in the upwave arc have generally larger stiffnesses than the ones in the downwave arc, but in this case there is the exception of the body located at y = 0 in the downwave row. There is no evident behaviour for dampings between the two arcs. The power increase obtained by optimization is of 3.99%; the interaction factor is 1.144. Notice that the slamming constraint is never active, due to the relatively large draft of the bodies. Due to this, positive interactions in the array are promoted and they are not limited by the constraint.   Figure 15. This time, the largest power increases are for the devices at the center of the array, which are also the ones with the largest control stiffnesses. The power increase is of 0.41%, and the obtained interaction factor is 0.816. A 4 × 4 square arrangement has also been simulated for the case 2 considered in the previous section. The results are shown in Figures 16 and 17. Devices in the most downwave row are subject to a power reduction, while all the others see a power increase. Bodies in the third row are assigned negative stiffnesses. The global power is increased by 9.16% and the interaction factor is 0.835.
The process is repeated with the constraint of positive stiffnesses: see Fig. 18, 19. In terms of total power, the result is around 76.3 kW in both cases, but in this second

Extension to multiple random parameters
The methods presented in the previous sections can be generalized to the case of multiple random variables. We present an application of practical interest involving two random parameters, based on a discussion presented in [8]. This allows us to discuss a case where our optimization framework computes an optimal solution that is, from a physical point of view, less intuitive than the ones presented in the previous sections. Consider a sea state formed by a superposition of wind waves and swell waves. Wind waves  Figure 18: Map of power before (left) and after (right) array optimization; case 2, 16 devices, positive stiffness constraint are generated locally by the action of the wind on the free surface, while swell waves originate from storms occurring far from the measurement point. The frequency spectrum of such sea state typically presents two peaks: a relatively narrow, low-frequency peak corresponding to the swell component, and a relatively wide, high-frequency peak corresponding to the wind component. Assuming linearity of the interaction between the two wave systems, the spectrum can be written as S(f, θ) = S w (f, θ) + S s (f, θ), the subscripts w and s being referred to the wind contribution and the swell contribution, respectively. Consistently with the approximation adopted for the case of a single random direction, we write S w (f, θ) = S w (f )D w (θ), and the analogous expression for swell. The frequency part is given by the Pierson-Moskovitz spectrum (2.14), while the directional part is given by the Donelan distribution function (2.15).
To formulate the robust optimization problem in this case, we assume that peak We define a realization as the sum of two irregular waves, each one unidirectional with directions θ w and θ s , sampled from D w (θ) and D s (θ), respectively. The excitation force acting on each device is, accordingly, the sum of the excitation forces due to the wind and swell contributions. We present a case with the 15-device configuration of type 2 (see Tab Fig. 21. Let us discuss the obtained results from a physical point of view. First, it is possible to see that the distributions of damping and stiffness parameters is not symmetric with respect to the x axis. This is due to the asymmetry in the wave directions. More interestingly, one can observe that the parameters are not clearly clustered: even though one may expect to obtain two clusters, dealing with the two

Discussion and conclusions
In this work, a new computational framework for the robust optimization of arrays of WEC was proposed and tested. In particular, we modelled an optimization problem for the maximization of the power of WEC parks with respect to the individual control damping and stiffness coefficients of each device. The results are robust with respect to the incident wave direction, which is treated as a random variable. The numerical solution to this stochastic optimization problem is obtained as a combination of penalty iterations and stochastic approaches, namely robust SA and SAA based on Monte Carlo or Gauss-Legendre quadrature integration.
The results of our numerical simulations are in agreement with the ones presented (for deterministic problems and unidirectional waves) in classical references in the field of marine engineering, like, e.g., [5], where it was observed that downwave rows are in general less productive than upwave rows. This behaviour is made even more evident by optimization, especially in square arrangements of bodies, where the difference in power between the first and the last row generally increases. Upwave bodies extract more power than they would in isolated conditions, so they benefit from the interactions with downwave bodies. These positive interactions might be impeded or reduced by the slamming constraint. In particular, it is possible that bodies tuned to have admissible oscillations individually (i.e., in isolated conditions) will violate the constraint when installed in an array. In our tests, we observed that when the initial guess is feasible, then the proposed optimization framework leads to performance increases of a few percentage points. When the initial guess is unfeasible, the method leads to a condition that is feasible and that might have a smaller power than the unfeasible initial guess. In most of the explored cases, allowing the stiffness to have negative values provided larger power that enforcing positive stiffness. In this case, technological feasibility considerations should be made about whether such solution is realizable. Finally, by proposing an example with a double-peaked wave spectrum, we showed that our computational framework can be used to control and optimize more complex configurations, so that it is capable to provide useful indications for the preliminary control design of realistic WEC arrays.