Data-Driven Model Predictive Control for Wave Energy Converters Using Gaussian Process

Abstract

The energy harvested by an ocean wave energy converter (WEC) can be enhanced by a well-designed wave-by-wave control strategy. One of such superior control methods is model predictive control (MPC), which is a nonlinear constrained optimization control strategy. A limitation of the classical MPC algorithm is its requirement of an accurate WEC dynamic model for real-time implementation. This article overcomes this challenge by proposing a data-driven MPC scheme for wave energy converters. The data-based WEC model is developed by a Gaussian process (encompassing mean predictions and symmetric uncertainties) for a more accurate description of nonlinear and unmodeled system dynamics. A cross-entropy solver for data-driven MPC is employed for rapid, high-performance results, which samples trajectories from Gaussian distributions based on the concept of the symmetry principle. The proposed strategy is verified numerically by simulations which demonstrate its superior performance over a classical complex-conjugate controller.

1. Introduction

Wave energy is a widespread and promising renewable source of energy from which electricity can be generated without any carbon dioxide emissions. Nonetheless, appropriate control strategies are critical to optimize wave energy absorption efficiency and guarantee operational safety at minimal hardware expense. This is especially important for point-absorber wave energy converters (PAWECs). PAWECs can generally be described as oscillators excited by waves, whose horizontal dimensions are much smaller than the prevailing wavelength [1].

Many control strategies have been proposed to achieve a practically implementable optimal power objective. The classical WEC controller design methods are mainly developed based on the principle that the optimal energy output can be achieved when the resonance frequency of the WEC adaptively matches the dominant frequencies of the incident waves under control. Classical WEC controllers include impedance matching [2], complex-conjugate [3], latching [4], declutching [5] controllers, etc. These methods are incapable of handling WEC physical constraints effectively. To maximize energy absorption while maintaining safe operation, WEC control is essentially a constrained optimal control problem subject to persistent disturbances [6], which can be solved by model predictive control (MPC) or MPC-like control algorithms [7]. MPC is an online optimization technique which generates an optimal control sequence by solving a constrained optimization problem at each sampling instant. Since MPC intrinsically handles multiple objectives and constraints, it is suitable for WEC control. Nonetheless, there are still some challenges with the practical application of MPC of WECs. The majority of MPC-based WEC controllers utilize the widespread linear WEC model based on Cummins’ equation [8]. This is developed with the assumption that the device has only small movements around the equilibrium position (the stationary water) [9]. Nonetheless, this linear model is impractical since the WEC is expected to have large movements to maintain resonance with the wave excitation force. The large WEC motion can result in strong nonlinear dynamics such as viscous force on the device, which is verified to be non-negligible in the MPC design [10]. To this end, nonlinear MPC (NMPC) control schemes are introduced to deal with the nonlinear effects of WEC dynamics. To achieve the alleviation of computational burden and better control performance, nonlinear buoyancy force effects are incorporated into online system optimization by the pseudospectral approach which exploits the flatness property in the constructed NMPC scheme [11]. The nonlinear WEC model, considering mooring forces, is linearized and optimized in the NMPC framework [12], which indicates that the nonlinear mooring force alone cannot fully realize the potential of NMPC. However, most of these nonlinear models account for some weak nonlinear effects, which just are approximated by the linear approach and fail to represent the “unmodeled dynamics” of the system well. Some strong nonlinear hydrodynamics, especially the viscous damping, mooring force and Froude–Krylov forces should be well-described to achieve a more accurate WEC model. Therefore, it is important to propose an appropriate WEC model including nonlinear dynamics for better control performance.

Data-driven methods have the capability to capture the key nonlinear wave–WEC interaction hydrodynamics and do not require physical information, i.e., structure parameters, as they are solely based on system measurements such as body position, velocity, wave elevation, etc. Recently, there has been an explosion of data-driven modeling approaches using machine learning techniques, which presents good generalization capabilities as well as statistical inference, making them good at mapping highly nonlinear physical relationships [13]. Due to the complex dynamics and high nonlinearity of WECs, it is desired to introduce machine learning methods into WEC modelling for better control performance. The work in [14] proposed an artificial neural network for wave prediction, which is further combined with a linear real-time MPC controller. The work in [15] presented a model-free deep reinforcement learning controller for WEC that outperforms the model-based control in terms of wave energy production, but it is rather computationally expensive. The authors of [16] presented a novel data-driven approach for updating reactive control parameters based on a Multifidelity Gaussian process model for a WEC device. The machine learning applications for WEC control are focused on wave force forecasting [14], wave excitation force estimation [17], complex wave hydrodynamic approximation [18], direct controller optimization using a reinforcement learning algorithm [15,19], etc. However, the use of machine learning techniques to describe overall WEC dynamics and to study corresponding control methods is relatively underdeveloped in the literature. This study fills this gap. Gaussian Process (GP) is a powerful kernel-based learning method and presents strong potential for analyzing implicit patterns between a series of training datasets. This approach is good at dealing with complex nonlinear regression problems using nonparametric models. The GP method provides the advantages of modeling flexibility, uncertainty estimation as well as learning smoothness and noise parameters from a training dataset [20,21]. Therefore, the GP technique is employed in this study to model the overall WEC dynamics. To optimize the performance of this complex system, the data-driven model predictive control is introduced and formulated for the proposed GP-based WEC model that maximizes the wave energy absorption. Meanwhile, the physical limitations on the body movements and control efforts are considered in designing this data-driven MPC.

Above all, the present study is instructive for developing a data-driven MPC approach relying on black-box identification of a GP model from the WEC input/output datasets. The contributions are illustrated as follows:

• A novel data-driven WEC model using machine learning techniques and targeting the control perspective is proposed, promising to advance state-of-the-art WEC modelling. The PAWEC system dynamics are learned by the Gaussian Process model, which aims to capture the nonlinear system characteristics with mean value and uncertainties.
• Developing a new data-driven MPC scheme based on the GP model for efficient and real-time implementation in the actual operation of WEC. The cross-entropy technique is introduced to deal with the trajectory optimization for fast, sample-efficient and high performance.
• The investigation of the performance of the data-driven MPC compared with the classical complex-conjugate controller is expected to fill the gap in the literature.
• The developed GP-based MPC scheme is validated in a small-sized and single-type WEC in this study, which can generally be applied to any WECs across different deployment prototypes (e.g., sizes, shapes) and other energy-maximizing control problems.

The arrangement of the paper sections is as follows. Section 2 presents the modeling of point absorbers, including classical linear state-space modeling and Gaussian-Process-based modeling methods. Section 3 introduces the complex-conjugate controller and data-driven MPC scheme using the CEM solver. The detailed simulation results are presented in Section 4 and finally the conclusions in Section 5.

2. Classical WaveStar PAWEC Modelling and Gaussian Process Regression

2.1. Classical WaveStar PAWEC Modelling

In this work, a specific type of WEC called point absorber (i.e., WaveStar WEC [22]) is employed to demonstrate the efficacy of the proposed strategy. The WaveStar WEC is a multiabsorber concept utilizing several hemispherical floats attached to a single platform. Figure 1 illustrates the WaveStar in the WEC-Sim simulator considering a point absorber. Assuming linear wavetheory and small body motions, the PAWEC response can be achieved from the superposition of the inertial, hydrostatic, viscous, radiation, PTO control and excitation forces. The floater–wave dynamic response (motions and forces) in WEC-Sim is achieved by calculating the motion equation for each body around the gravity center according to the Cummins’ equation [8,23]:

$$\begin{array}{cc}\hfill (m+A_{\infty })\ddot{X}\left(t\right)& =-{\int }_0^t{K}_r(t-\tau )\dot{X}\left(\tau \right)d\tau \hfill \\ \hfill & +F_{ext}\left(t\right)+F_{vis}\left(t\right)+F_{hs}\left(t\right)+F_{pto}\left(t\right)\hfill \end{array}$$

(1)

where m represents the buoy mass, $A_{\infty }$ represents added buoy mass at infinite wave frequency, $X\left(t\right)$ represents the buoy displacement vector, $K_r$ represents the radiation impulse response dynamic, $F_{ext}\left(t\right)$ denotes the wave excitation force acting on the floater, $F_{hs}\left(t\right)$ presents the hydrostatic restoring force, $F_{pto}\left(t\right)$ represents the PTO force and $F_{vis}\left(t\right)$ is the quadratic viscous drag term computed using Morison’s equation. Other effects, including mechanical friction and mooring force, are neglected.

Moreover, the WaveStar moment around the point A in Figure 1 can be expressed using a Newtonian rotation equation [1]:

$$J\ddot{\theta }\left(t\right)=M_r\left(t\right)+M_{ext}\left(t\right)+M_{vis}\left(t\right)+M_{hs}\left(t\right)+M_{pto}\left(t\right),$$

(2)

where J denotes the moment of the rotating body’s inertia, $\ddot{\theta }\left(t\right)$ is the time derivative of the buoy angular velocity $\dot{\theta }\left(t\right)$, $M_r\left(t\right)$ is the moment of the radiation force and similar explanations follow for the other moments.

Generally, a finite state-space subsystem of the 2nd order is adopted to approximate the radiation moment [24], illustrated in Equation (3). In Equation (3), $(A_r;B_r;C_r;D_r)$ are the state-space matrices for the variable $x_r\left(t\right)$, which is an internal state without specific physical meaning and is achieved by Prony’s method to express the inconvenient product $M_r\left(t\right)={\int }_0^t{K}_r(t-\tau )\dot{X}\left(\tau \right)d\tau$ [24].

$$\begin{array}{c}\hfill \begin{array}{cc}& {\dot{x}}_r\left(t\right)=A_r{x}_r\left(t\right)+B_r\dot{\theta }\left(t\right),\hfill \\ & M_r\left(t\right)=C_r{x}_r\left(t\right)+D_r\dot{\theta }\left(t\right).\hfill \end{array}\end{array}$$

(3)

where

$$\begin{array}{cc}\hfill & A_r=\left[\begin{array}{cc}-13.59& -13.35\\ 8.0& 0\end{array}\right],\hfill \\ \hfill & B_r=\left[\begin{array}{c}8\\ 0\end{array}\right],\hfill \\ \hfill & C_r=\left[\begin{array}{cc}4.73& 0.50\end{array}\right],\hfill \\ \hfill & D_r=-0.1586.\hfill \end{array}$$

Note that the wave excitation force $F_{ext}\left(t\right)$ and its relation to $M_{ext}\left(t\right)$ are as follows, where $L_A$ denotes the distance from gravity center to point A in Figure 1 [25]:

$$M_{ext}\left(t\right)=F_{ext}\left(t\right)L_A$$

(4)

It is assumed that the small displacement $M_{hs}\left(t\right)$ has an approximate proportional relationship with the angular displacement $\theta \left(t\right)$, where $K_{hs}=92.33{Nmrad}^{-1}$ represents the hydrostatic factor [22]:

$$M_{hs}\left(t\right)=-K_{hs}\theta \left(t\right)$$

(5)

Moreover, the linear term can be used to approximate the nonlinear viscous moment $M_{vis}\left(t\right)$ with a linear damping coefficient $K_v=1.8{Nmrad}^{-1}{s}^{-1}$ [22]:

$$M_{vis}\left(t\right)=-K_v\dot{\theta }\left(t\right)$$

(6)

Using Equations (2)–(5), the classical WaveStar PAWEC model can then be summarized in a linear state-space expression, illustrated as:

$$\begin{array}{cc}\hfill & \dot{x}\left(t\right)=Ax\left(t\right)+B(u\left(t\right)+F_{ext}\left(t\right)),\hfill \\ \hfill & y\left(t\right)=Cx\left(t\right),\hfill \\ \hfill & x\left(t\right)={\left[\theta \left(t\right)\dot{\theta }\left(t\right)x_r\left(t\right)\right]}^T,\hfill \end{array}$$

(7)

where

$$\begin{array}{cc}\hfill & A=\left[\begin{array}{ccc}0& 1& {\mathbf{0}}_{1\times 2}\\ -\frac{K_{hs}}{J_t}& -\frac{K_v}{J_t}& -\frac{C_r}{J_t}\\ {\mathbf{0}}_{2\times 1}& B_r& A_r\end{array}\right],\hfill \\ \hfill & B=\frac{1}{J_t}\left[\begin{array}{c}0\\ 1\\ {\mathbf{0}}_{2\times 1}\end{array}\right],\hfill \\ \hfill & C=\left[\begin{array}{ccc}1& 0& {\mathbf{0}}_{1\times 2}\\ 0& 1& {\mathbf{0}}_{1\times 2}\end{array}\right].\hfill \end{array}$$

(8)

where $u\left(t\right)=F_{pto}\left(t\right)=\left(M_{pto}\left(t\right)/L_A\right)$ denotes the control input, $J_t=J+J_{\infty }=1.36$ kgm² means the equivalent moment of inertia and ${\mathbf{0}}_{m\times n}$ is a zero matrix with m rows and n columns.

The electromechanical conversion efficiency $\eta$ of the linear generator adopted here is 0.7. In this case, the absorbed wave power by the WaveStar is depicted as [24]:

$$P\left(t\right)=-{\eta }_{pto}{F}_{pto}\left(t\right)\dot{X}\left(t\right),{\eta }_{pto}=\left\{\begin{array}{ccc}\eta \hfill & if\hfill & -F_{pto}\left(t\right)\dot{X}\left(t\right)>0\hfill \\ 1/\eta \hfill & if\hfill & -F_{pto}\left(t\right)\dot{X}\left(t\right)\le 0\hfill \end{array}\right.$$

(9)

2.2. Gaussian-Process-Based Modeling Method

Gaussian process is a data-driven and nonparametric Bayesian approach to obtain regression functions from the provided sample dataset. The GP regression provides the advantages of modeling flexibility, uncertainty estimation as well as learning smoothness and noise parameters from a training dataset [21]. Generally, a GP represents the function distribution based on the training dataset. A GP is a stochastic process including an infinite parameter set, any finite subsets of which are jointly Gaussian-distributed. Then, the priori statistics of a GP model $f\left(x\right)$ can be described by the mean function $m\left(x\right)$ and the related covariance function $k(x,x^{*})$:

$$\begin{array}{cc}\hfill f\left(x\right)\sim & \mathcal{G}\mathcal{P}(m\left(x\right),k(x,x^{*}))\hfill \\ \hfill m\left(x\right)& =\mathbb{E}\left[f\right(x\left)\right]\hfill \\ \hfill k(x,x^{*})& =cov(f\left(x\right),f\left(x^{*}\right))\hfill \end{array}$$

(10)

where $x\in {\mathbb{R}}^D$ denotes the dynamic input vector, D is the input dimension and $f\left(x\right)$ and $f\left(x^{*}\right)$ are random Gaussian scalar variables indexed by inputs x and $x^{*}$. Typically, $k(x,x^{*})$ is called a kernel function parametrized by some specific hyperparameters. The covariance function exhibits symmetry and non-negative determinism, which describes the process behavior and defines the proximity between arbitrary random points of the Gaussian function [26]. The regression problem here is equivalent to inferring the features of the unknown PAWEC model from the input and output time series.

Assume a training dataset $\mathcal{D}=(X,Y)$ in a GP regression problem, where the matrix $X=[x_1,x_2,\dots ,x_n]$ presents the input vector and vector $Y=[y_1,y_2,\dots ,y_n]$ denotes the corresponding scalar outputs. The GP model can be regarded as a prior within the range of plausible dynamic functions. After being updated through the training dataset, the GP model will be a posterior over functions. The achieved GP posterior model is defined as:

$$y_i=f\left(x_i\right)+{ϵ}_i{ϵ}_i\sim \mathcal{N}(0,{\sigma }^2)$$

(11)

where y, $f\left(x\right)$, and $ϵ$ represent the observed outputs, GP model values (including uncertainty estimates) and a white noise (errors) with a mean of 0 and variance of ${\sigma }^2$, respectively.

With n-many training data points in the training input matrix $x\in {\mathbb{R}}^{N\ast D}$ and observation target vector $y\in {\mathbb{R}}^N$, any arbitrary collection of function values $f\left(x_i\right)$ is jointly Gaussian-distributed:

$${[f\left(x_1\right),f\left(x_2\right),\cdots ,f\left(x_N\right)]}^T\sim \mathcal{N}(\mu ,K+{\sigma }^{2}I)$$

(12)

where covariance function $K=K_{ij}=k(x_i,x_j)$ with $N*N$ entries and mean function $\mu ={\mu }_i=m\left(x_i\right)$ with N entries.

The squared exponential (SE) is the adopted covariance function in this study, which is referred to as the Radial Basis Function (RBF) [20]. Any continuous function can be approximated when the number of RBFs is large enough. This is expressed as [16]:

$$k(x_i,x_j)=h^{2}exp\left[-{\left(\frac{x_i-x_j}{\lambda }\right)}^2\right]$$

(13)

where two related hyperparameters h and $\lambda$ decide the model input and output scales. Large values of $\lambda$ enhance the smoothness of the selected kernel function and vice versa.

The hyperparameter values h and $\lambda$ in the SE covariance can be achieved by the optimization of the log marginal likelihood function [20]) as follows:

$$logp\left(y\right|\theta )=-\frac{1}{2}log\left|K\right|-\frac{1}{2}{y}^T{K}^{-1}y-\frac{n}{2}log\left(2\pi \right)$$

(14)

Standard gradient nonconvex optimization strategies such as Rdrop, CG or Broyden–Fletcher–Goldfarb–Shanno (BFGS) are normally adopted in the hyperparameter optimization. It is noteworthy that highly appropriate initial hyperparameter values before the training process can promote faster convergence of optimization and avoid convergence to an unsatisfactory regional optimum.

The new target prediction $f^{*}$ for a given input $X^{*}$ from the posterior can be achieved after the training process using the extended joint distribution in Equation (15). The mean can be interpreted as a deterministic prediction, while the variance can be viewed as a measure of the uncertainty for that deterministic prediction.

$$\left[\begin{array}{c}{f}^{*}\\ Y\end{array}\right]\sim \left(\left[\begin{array}{c}m\left(X^{*}\right)\\ m\left(X\right)\end{array}\right],\left[\begin{array}{cc}k(X^{*},X^{*})& k(X^{*},X)\\ k(X,X^{*})& K+{\sigma }^{2}I\end{array}\right]\right)$$

(15)

where $k(X^{*},X)=k{(X,X^{*})}^T=[k(X_1,X^{*}),\cdots ,k(X_N,X^{*})]$. According to the Joint Gaussian Distribution Theorem [20], the forecasting result for the new target becomes:

$$\begin{array}{c}\hfill \mu \left(f^{*}\right)=m\left(X^{*}\right)+k(X^{*},X){[K+{\sigma }^{2}I]}^{-1}(Y-m\left(X\right))\\ \hfill \sigma \left(f^{*}\right)=k(X^{*},X^{*})-k(X^{*},X){[K+{\sigma }^{2}I]}^{-1}k(X,X^{*})\end{array}$$

(16)

3. Control for Optimal Power Extraction from WEC

3.1. Complex-Conjugate Control

Complex-conjugate control is a classical method to extract the optimal mechanical power of PAWECs, which achieves energy absorption maximization through impedance matching. As required, the PTO can be assumed to work as an additional spring or inertia to compensate for the inherent WEC reactance. The WEC floater can move in phase with the wave excitation force (WEF) by adding a PTO force. Therefore, optimal power absorption can be achieved when the complex conjugate of the device’s mechanical impedance is matched to the PTO impedance, or when the PTO’s reactance compensates the device’s reactance. For convenience of control, the WEC-Sim floater–wave dynamic (1) is typically expressed in the frequency domain as [27]:

$$\left(i\omega (m+A_{\omega })+B_v+R\left(\omega \right)+\frac{S}{i\omega }\right)i\omega X\left(\omega \right)=F_{ext}\left(\omega \right)+F_{pto}\left(\omega \right)$$

(17)

where $X\left(\omega \right)$ denotes the Fourier transform of the corresponding time-domain signal, $B_v$ and $R\left(\omega \right)$ are the PTO system viscous and radiation damping, respectively, and S represents the hydrostatic restoring coefficient matrix.

Furthermore, the WEC system intrinsic impedance is defined as:

$$Z_i\left(\omega \right)=i\omega (m+A_{\omega })+B_v+R\left(\omega \right)+\frac{S}{i\omega }$$

(18)

According to linear potential flow theory [1], the optimal PTO impedance corresponds to the complex-conjugate of the intrinsic mechanical impedance of the WEC system in the frequency domain, i.e.,

$$\begin{array}{cc}\hfill Z_{pto}\left(\omega \right)& =Z_i^{*}\left(\omega \right)\hfill \\ \hfill & =-i\omega (m+A_{\omega })+B_v+R\left(\omega \right)-\frac{S}{i\omega }\hfill \end{array}$$

(19)

In the wave complex-conjugate control, the system PTO force can be derived by adding the mass, damping and spring terms, illustrated as:

$$\begin{array}{cc}\hfill F_{pto}\left(t\right)& =M_{pto}\ddot{X}\left(t\right)+B_{pto}\dot{X}\left(t\right)+C_{pto}X\left(t\right)\hfill \\ \hfill & =-(m+A)\ddot{X}\left(t\right)+(B_v+R)\dot{X}\left(t\right)-SX\left(t\right)\hfill \end{array}$$

(20)

The PTO spring coefficient $C_{pto}$ changes with the PTO geometry. Given the PTO geometry, $C_{pto}$ turns out to be a constant value. However, the values of $M_{pto}$ and $B_{pto}$ need to be regulated with different $\omega$ values. Assuming $B_{pto}=B_v+R$, Equation (20) becomes:

$$M_{pto}\ddot{X}\left(t\right)+C_{pto}X\left(t\right)=-(m+A)\ddot{X}\left(t\right)-SX\left(t\right)$$

(21)

There are two unknown coefficients in Equation (21), namely $M_{pto}$ and $C_{pto}$. When a random certain value is provided for one of these two unknown coefficients, the other coefficient can be computed according to Equation (21). When $M_{pto}=0$ in Equation (21), it is referred to as the complex-conjugate control [28]:

$$F_{pto}\left(t\right)=B_{pto}\dot{X}\left(t\right)+C_{pto}X\left(t\right)$$

(22)

The values of PTO damping and spring coefficients remain unchanged for a single wave frequency. The complex-conjugate control is known to be noncausal, which means that future knowledge of the wave excitation force is required for physical implementation [29]. The WEF forecasting can be approximated by correlated wave measurements using different methods [30], but this is not the focus of this study. Therefore, ideal future WEF information in the time domain is assumed in this paper. To summarize, Figure 2 provides the illustration of the computation of WEC device motions and corresponding complex-conjugate controller in a block diagram.

3.2. Data-Driven MPC Design with Cross-Entropy Optimization

3.2.1. Cross-Entropy Optimization

To determine the control actions, the MPC solved by the learning model searches the optimal action sequence for the system and executes the first action of the obtained sequence, discarding the remaining control actions. Normally, this search is repeated after each step in the real environment to resolve errors in model predictions and enable interaction with the environment. In this study, the planning step is achieved with the Cross-Entropy Method (CEM), a fast, sample-efficient and high-performing algorithm. The CEM algorithm is a global derivative-free optimization technique and a sampling-based trajectory optimization method applied to model predictive control, which presents symmetry [31]. It provides several appeals including the possibility to optimize complex black-box functions, higher robustness, no need for gradient knowledge and lower sensitivity to local optima.

In the CEM-based MPC design, CEM functions at each time step to optimize the h-step planning problem over action sequences [32]. CEM iteration is referred to as a step in the inner loop of CEM that implements the optimization of the sampling distribution. The outer loop step represents progress in the real environment by performing a control action. Typically, the next step is to consider the planning problem one time step later. The CEM algorithm can be viewed as an evolution scheme which optimizes the system cost function by finding appropriate “individuals”. A sample of individuals is drawn from the population distribution as control actions for MPC. These individuals are then evaluated and sorted according to the cost function and a fixed number of “elite” candidates are selected as the elite set. This elite set determines the population parameters for the next iteration. In the standard CEM case, the population is modeled by a Gaussian probability distribution based on the concept of the symmetry principle with mean $\mu$ and a diagonal covariance matrix $diag\left({\delta }^2\right)$, where $\mu ,\delta \in {\mathbb{R}}^n$. By fitting $\mu$ and $\delta$ to the elite set, the sampling distribution is centered around the near-optimal solutions with high performance. After multiple selection processes, the best solution close to the global optimum is obtained.

To achieve strong performance with few samples, reduce computational cost and make it suitable for real-time trajectory optimization in MPC, several improvements can be employed to CEM [33]. Ideally, CEM action samples should yield trajectories that explore the state space to a large extent. The momentum term $\alpha$ [34] is employed to refit the distribution between the CEM iterations, aiming to estimate parameters of the sampling distribution using only a small elite set. Here, ${\mu }_t^{i+1}=\alpha {\mu }^i+(1-\alpha ){\mu }_{elitese{t}_i}/K$, where $\alpha \in [0,1]$, i denotes the index of CEM iterations and K is the size of the elite set. Actions are sampled from the unmodified normal Gaussian distributions and the results are clipped to lie inside the permitted action interval, which allows to sample maximal actions more frequently. Once the inner optimization loop is completed, the optimized Gaussian distribution and the elite set produced at each inner CEM iteration are stored, while a small fraction of them are added to the pool for the next iteration, rather than discarding the elite set during each CEM iteration. Moreover, a small fraction of the elite set of the last CEM iteration is stored and added with each stochastic action to employ it in the next step. Here, the elite reuse fraction $\xi =0.3$ is introduced due to the fact that when an optimum is closed, the standard deviation automatically reduces, thus narrowing down the search space and fine-tuning the final solution. It turns out to be sufficient to sample fewer action sequences as the CEM iterations progress. Therefore, an exponential decrease in the sample size of a constant coefficient $\gamma$ is introduced. The population size of iteration i is $N_i=max(N{\gamma }^{-i},2K)$, where $max$ ensures that the population size is at least twice the size of the elite set. In this CEM algorithm, the first action of the best seen action sequence is performed. The overall CEM algorithm is summarized in Algorithm 1.

Algorithm 1: The CEM optimization algorithm.

3.2.2. Data-Driven MPC Formulation

In this section, the GP-based MPC controller aims to find optimal control actions along the prediction horizon according to the given performance objectives and physical constraints. For the WEC optimization problem, the data-driven MPC minimizes the discrete form of the finite-horizon cost function of the nonlinear energy output criterion over a certain prediction horizon $T_H$, aiming to maximize the energy output:

$$J\left(t\right)=-{\int }_t^{t+T_H}P\left(t\right)dt=-{\int }_t^{t+T_H}{\eta }_{pto}{F}_{pto}\left(t\right)\dot{X}\left(t\right)dt$$

(23)

with

$$\dot{X}\left(t\right)=L_A\dot{\theta }\left(t\right)$$

(24)

where moment arm $L_A$ is the known and variable parameter. The mean value and covariance of $\dot{X\left(t\right)}$ can be derived from the GP model of $\dot{\theta }\left(t\right)$ according to Equations (24) and (25). The dynamic PAWEC system is illustrated by the GP regression:

$$\Delta x\sim \mathcal{G}\mathcal{P}(\Delta F_{ext},\Delta F_{pto})$$

(25)

where $\Delta x=[\Delta \theta ,\Delta \dot{\theta }]$ denotes the output increments of the GP-based PAWEC model and $\Delta F_{ext}$ and $\Delta F_{pto}$ represent the input increments of the GP-based PAWEC model. The MPC makes full use of the probability knowledge of the GP-based WEC model. The data-driven GP-MPC scheme is a sampling-based algorithm, and the optimal PTO force $F_{pto}(k+i|k),i=1,2,\cdots ,H$ is provided by the sampling-based CEM trajectory optimization algorithm. To achieve noncasual control, H samples of wave excitation force prediction $F_{ext}$ are assumed to be available. In this sense, $\Delta F_{ext}$ and $\Delta F_{pto}$ can be achieved by the difference between the previous moment and the next moment.

The data-driven MPC using GP-based PAWEC model can be illustrated as follows:

$$\begin{array}{cc}\hfill & arg\underset{F_{pto}(·)}{max}J\left(k\right)\hfill \\ \hfill with:& J\left(k\right)=-\sum _{i=1}^H{\eta }_{pto}{F}_{pto}(k+i|k)\dot{X}(k+i|k)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill s.t.& X_{k|k}=X_0\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & X_{min}\le X_{k+i|k}\pm 2{\sigma }_X\le X_{max}i=0,1,\cdots ,H\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & {F_{pto}}_{min}\le F_{pto}(k+i|k)\le {F_{pto}}_{max}\hfill \end{array}$$

(29)

Let $X_{k+i|k}$ represents the prediction of X at sample $k+i$ where the prediction is achieved at sampling time k, the PTO force $F_{pto}(k+i|k)$ is the optimal control input at sample $k+i$ obtained at sample k and H is the prediction horizon. Equation (27) indicates the initial knowledge of states which is assumed to be available. Equation (28) represents the buoy position constraints $X_{min}=-0.08 m$ and $X_{max}=0.08 m$. Given that this buoy position X prediction is achieved from the GP, including mean ${\mu }_X$ and covariance ${\sigma }_X$, the truncated GP $X_{k+i|k}\pm 2{\sigma }_X$ is employed to guarantee that the estimated uncertain X is within the safe region. Equation (29) represents the control input constraints that can be provided by the PTO mechanism ${F_{pto}}_{min}=-80 N$ and ${F_{pto}}_{max}=80 N$.

This GP model (Equation (25)) can realize the sequence-to-sequence state prediction for MPC optimization, avoiding the error accumulation. That is to say, ${[\Delta x_{k+1|k},\Delta x_{k+2|k},\cdots ,\Delta x_{k+H|k}]}^T$ can be obtained at the same k time from the GP regression (Equation (25)). Therefore, $\dot{\theta }(k+i|k),i=1,2,\cdots ,H$ can be achieved as follows:

$$\begin{array}{cc}\hfill x_{k+1|k}& =\Delta x_{k+1|k}+x_{k|k}\hfill \\ \hfill x_{k+2|k}& =\Delta x_{k+2|k}+x_{k+1|k}\hfill \\ \hfill ⋮\\ \hfill x_{k+H|k}& =\Delta x_{k+H|k}+x_{k+H-1|k}\hfill \end{array}$$

(30)

Moreover, $\dot{X}(k+i|k),i=1,2,\cdots ,H$ can be obtained through the relationship between $\dot{X}\left(t\right)$ and $\dot{\theta }\left(t\right)$ (Equation (24)). The proposed GP-based data-driven MPC strategy is summarized in Figure 3, where the stop criteria indicates that the iteration number reaches CEM iterations. In summary, the optimal PTO force can thus be achieved from the designed data-driven MPC strategy.

4. Simulations

The numerical simulations to demonstrate the effectiveness of the proposed data-driven MPC using the GP modeling method are provided in this section. The WaveStar model in the WEC-Sim simulator is employed here. The sampling time is 0.1 s. The 100 s record of buoy angular position is utilized as the training dataset for GP model fitting. With the help of this step, the final GP model ($posterior$) is achieved from the initial GP model ($prior$) through updating the relevant hyperparameters h and $\lambda$ in Equation (13). Samples of 20 s of buoy angular position and velocity (i.e., [100, 120] s) are estimated to test the performance of the resulting GP regression model, as illustrated in Figure 4, where the second subplot provides the zoomed-in regression result. Apart from the obtained mean values, uncertainty estimates of state points are also provided. It can be seen that the buoy angular position $\theta$ outputted by the GP regression model can follow the true value very closely. It is interesting to note that the uncertainty range is not large, implying that the mean value estimates are accurate.

Moreover, the GP regression results for a 350 s duration are compared with the outputs of the linear state-space model and the real WEC-Sim nonlinear model outputs, and part of the simulation results are as illustrated in Figure 5. It can be seen that the buoy angular position and velocity estimated by the GP model are very similar to the actual output values of WEC-Sim, presenting much better regression performance compared with that of the linear state-space model. In contrast, the state-space model presents worse regression results in terms of high-frequency components. This is because the GP model provides a more complete description of the PAWEC nonlinear effects and unknown properties. The results verify the feasibility and validity of the designed GP model for the WaveStar model. Moreover, the GP regression performance is evaluated by normalized mean squared error (NMSE). The smaller the NMSE is the better the performance. The NMSE results are illustrated in

Table 1. NMSE comparison results of the WaveStar model state regression.

Comparison Items	GP Regression Model	State-Space Model
NMSE of angular position ($\theta$)	0.2971	0.9951
NMSE of angular velocity ($\dot{\theta }$)	0.3609	0.9607

. The NMSE results of the GP-based model are smaller than those of the linear state-space model, which validates the effectiveness of the proposed GP-based model.

The classical complex-conjugate controller is designed according to the authors’ previous work [35]. The prediction horizon is $H=20$. CEM iteration is set to 5. The size of elite set is 20 and the sample decay is $\gamma =1.25$. The last time control output is set to the initial mean ${\mu }_{init}$ of the elite set, while the initial covariance ${\delta }_{elitese{t}_i}$ in the range $[1,12]$ increases as the prediction horizon increases to guarantee the continuity of control output. The irregular wave scenario with JONSWAP spectrum is adopted here to validate the proposed strategy, with a significant wave height of 0.1042 m, peak wave cycle of 1.836 s and peak enhancement parameter of 3.3, as illustrated in Figure 6.

The MPC controller is activated after 10 s of the wave fluctuation to avoid numerical instability during the initial transient. Under these settings, the control action output under the data-driven MPC during one CEM iteration is compared with the corresponding output of the complex-conjugate controller, as illustrated in Figure 7. The red line denotes the complex-conjugate controller output and the green line represents the 20 control actions from the elite set of the CEM-based MPC solver, while the other lines represent a subset of normal sample control actions. It can be seen that the PTO force generated by the complex-conjugate controller exceeds the control input constraints (i.e., [−80, 80] N) for a period of time (i.e., horizon steps [8, 15]), which violates the physical limits. Moreover, a small fraction of sample control actions by the data-driven MPC also exceed the floater physical constraints but are discarded. The control actions presenting better energy output constitute the final elite set. It is emphasized that the 20 control actions in the elite set are within the control limits. Furthermore, the control action $F_{pto}$, the corresponding system buoy position X and velocity $\dot{X}$ under the proposed data-driven MPC strategy are illustrated in Figure 8, where the black line denotes the best elite sequence of the elite set. It is worth noting that the optimal PTO force and the buoy velocity should be maximized under the premise of the same phase, thus realizing energy output optimization.

The buoy position during the simulation time is depicted in Figure 9, where the second subplot provides a zoomed-in figure between [300, 350] s. It can be seen that the buoy position X is within the physical limits [−0.08, 0.08] m. Moreover, the buoy position X computed by the data-driven MPC is always larger than that of the complex-conjugate controller, implying that the floater motion is largely exaggerated by the proposed MPC scheme, thereby achieving energy maximization.

Figure 10a,b provides the extracted instantaneous energy output and the accumulated power output of the two above-mentioned controllers, respectively. It can be seen that the proposed GP-based MPC can yield more instantaneous power. Especially, the proposed method, which can obtain 14% more energy than the classical complex-conjugate controller under irregular wave conditions for the considered simulation duration. The results reveal that with the help of the proposed GP-based MPC scheme the PAWEC can achieve more energy output, which verifies its effectiveness.

Moreover, when the proposed data-driven MPC strategy is applied in real practice, an initial energy maximization controller is required to generate the original training datasets for the GP-based WEC model regression. The relevant sensors are required to collect the measurement inputs/outputs. The measurement noise should be reprocessed before being injected to the GP model training. Once the training datasets are ready, the GP-based WEC model can be trained offline. The regression performance should be validated and compared with the real observations before the GP model is applied to the data-driven MPC.

5. Conclusions

This paper proposes a data-driven MPC strategy for PAWEC energy maximization that relies on a GP-based probabilistic model from black-box identification using input/output training datasets. The proposed GP-based PAWEC model can describe the system nonlinear parts and the unmodeled dynamics more accurately than classical linear state-space models, providing the system mean value and corresponding symmetric uncertainty. The resulting nonlinear MPC problem is solved and optimized by the CEM algorithm, making full use of the uncertainty information of the probabilistic GP-WEC model. This novel GP-based WEC model targets the control perspective and is combined with the data-driven MPC, which is expected to advance the state-of-the-art WEC energy-maximizing control scheme. The proposed method achieves better energy optimization performance than the conventional complex-conjugate controller; it can also prevent system constraint violations. The experimental results under the irregular wave conditions verify that the proposed strategy has the merits of better energy-maximization performance and better constraint compliance. The proposed GP-based model can include PTO dynamics in the future, completing a wave-to-wire WEC model, to improve the model-based energy-maximizing controller performance. Although the point absorber is employed as a case study, the proposed GP-based MPC using the CEM optimization framework can be applied to the control of other WECs and other similar energy-maximizing control problems.