A Comparative Study of Model Predictive Control and Optimal Causal Control for Heaving Point Absorbers

: Efforts by various researchers in recent years to design simple causal control laws that can be applied to WEC devices suggest that these controllers can yield similar levels of energy output as those of more complex non-causal controllers. However, most studies were established without adequately considering device and power conversion system constraints which are relevant design drivers from a cost and economic point of view. It is therefore imperative to understand the benefits of MPC compared to causal control from a performance and constraint handling perspective. In this paper, we compare linear MPC to a casual controller that incorporates constraint handling to benchmark its performance on a one DoF heaving point absorber in a range of wave conditions. Our analysis demonstrates that MPC provides significant performance advantages compared to an optimized causal controller, particularly if significant constraints on device motion and/or forces are imposed. We further demonstrate that distinct control performance regions can be established that correlate well with classical point absorber and volumetric limits of the wave energy conversion device.


Introduction
As the field of wave energy conversion transitions from traditional passive control techniques to advanced optimal control strategies for power maximization, it becomes increasingly important to understand the requirements, capabilities, limitations, and benefits of each method to choose the best optimization strategy for a given application.
The control system affects power capture, structural loads, and power-takeoff (PTO) design. To achieve true economic optimality in a wave energy conversion system design, the control system needs to be able to optimize performance given the constraints imposed by the wave energy converter (WEC) system. This is due to the fact that economic outcomes are optimized if a given component or sub-system is continually utilized at its rated operational condition during a majority of its operational life. These rated operational conditions provide constraints within which the device needs to operate, and the WEC control system is responsible for enforcing these constraints because exceeding them will result in mechanical breakdown. As such, the constraint handling of a controller is incredibly important to maximize the economic competitiveness of a WEC device. Constrained optimal control as offered by a Model Predictive Control (MPC) algorithm framework, which is also referred to as non-causal control, is a key tool in this optimization process.
Control systems for wave energy conversion are broadly classified as either causal or non-causal controllers. A causal controller uses information from sensors that monitor the outputs of a dynamical system, in our case the WEC. The sensor information is used as feedback by the controller to follow a desired command signal. Non-causal controllers such as MPC leverage wave prediction to generate control commands in a feed-forward manner rather than using feedback from sensors that measure system output variables.
Combination of non-causal (feed-forward) and causal (feedback) control is also possible in a hybrid manner and can be used to correct for modeling and/or wave prediction errors in non-causal controllers.
Within the literature there are two well-published methods for implementing causal control that incorporate constraint handling. The first approach tries to approximate complex conjugate control (CCC) ( [1][2][3][4]) with a causal feedback law [5]. This approach does well in the unconstrained case where no motion limits are imposed. To allow for constraints, these methods were subsequently augmented to incorporate an MPC controller that works for the sole purpose of enforcing constraints [5]. Performance of the CCC approximating causal control methods is sub-optimal when constraints are imposed, because the impedance matching value obtained from the CCC approximation provides in-sufficient damping for the constrained case.
The second approach for designing a causal controller comes from recognizing that the causal WEC control problems fall into the Linear Quadratic Gaussian (LQG) paradigm ( [6][7][8]). A non-standard LQG optimization problem is established with the objective of maximizing power. The problem treats the input wave excitation force as a stationary stochastic disturbance with a known spectrum. A wave-shaping filter is designed to model the spectrum and this filter model is incorporated in the optimization along with the system dynamics model of the WEC. The LQG optimal controller derived in this manner does not rely on approximations to CCC and will therefore achieve optimality in a broader set of sea states. Stroke protection techniques can be explicitly designed to protect against end-stop violations [7,9]. The control parameters are optimized offline and implemented using a gain schedule to adapt to changing wave conditions. A block-diagram of a casual controller with non-linear stroke protection is shown in Figure 1. Non-causal performance optimization of WECs is typically implemented using a receding horizon MPC framework. To guarantee optimality, MPC needs a forecast of oncoming waves 20 to 30 s into the future to plan optimal control commands. Wave forecasting can be conducted by placing probes up-wave, measuring the free-surface elevation in the wave field [10], or using a wave-sensing radar [11,12] to characterize the wave-field propagation. These up-wave measurements are marched forward in space and time to predict the wave elevation at the target WEC location. The predicted wave excitation force input, in addition to the states of a system dynamics model, are used to solve a quadratic optimization problem to find the optimal control trajectory. The states of the device are marched forward in time and the optimization process is repeated at each step with new prediction information. Figure 2 shows the basic MPC control architecture. Considerable research has been undertaken on expanding the capabilities of MPC to meet the requirements of an entire range of device topologies and PTO configurations ( [13][14][15][16][17][18][19]). Some of the common constraints that can be handled effectively within the MPC optimization include: (1) PTO force constraints; (2) WEC position, velocity, and acceleration constraints; (3) discrete (ON/OFF) control of PTO forces; and (4) power-flow constraints. MPC can accommodate non-linear device dynamics, such as viscous drag and non-linear models of powertrain losses, which need to be considered in the optimization problem to guarantee optimality [14].

The Device Model
We consider a one degree-of-freedom heaving point absorber as an example for our comparison study. A general schematic of the heaving buoy WEC is shown in Figure 3 and key dimensions are provided in Table 1. This device has well-established theoretical limits on power absorption which can be used to identify the design trade-off space for controls optimization.  The device geometry was modeled using a commercially available Boundary Element code called WAMIT [20]. The added mass, radiation damping, and excitation force kernels were obtained as a function of frequency by post-processing the WAMIT outputs. Figure 4 shows the frequency dependent parameters plotted against the wave period.
z is the velocity, and ..
z is the acceleration of the buoy in heave. m is the displaced mass of the buoy, F R is the radiation force, F e is the wave-excitation force, F m is the PTO force, F s is the hydrostatic restoring force given by F s = −k z, where k is the buoyancy stiffness, defined as k = ρgS, where ρ is the water density, g the gravitational constant, and S the water-plane area. F d is the drag force related to viscous effects in real fluid. The viscous damping coefficient was assumed to be linear according to the The radiation force f R (t), can be expressed as: where m ∞ is the infinite frequency added mass and h r (t) is the impulse response function of the radiation force. The excitation force F e (t) is expressed as: where h e (t) is the excitation force impulse response function and η(t) is the wave elevation at the device location. The impulse response function relating the wave elevation to the excitation force affecting the device is non-causal. The main reason is that the chosen input, i.e., the wave elevation at the device location, is not the direct cause of the output, i.e., the interaction force between the wavefield and the device. To recast the system dynamics into state-space form, the radiation term can be modelled as an embedded sub-system through the following state-space realization: This leads to the following state-space model: .

Model Predictive Control (MPC)
A wave energy converter with a linearized device dynamics model can be optimized using a linear MPC framework. In this type of problem, the non-linearities such as viscous drag are approximated using linear relationships. Under the assumption of no loss in the power generation process, optimizing the device average power absorbed P a at a given instant t 0 over a defined control horizon T h can be achieved by determining the optimal control trajectory u * (t) that maximizes the following cost function: where v(t) is the instantaneous device velocity. The minus sign is due to the convention of considering absorbed energy with a negative sign. Assuming a fixed time step, we discretize the integral and set up the optimization problem as a minimization by changing the sign in Equation (9). The new optimization objective, therefore, is given by: where N is the number of time intervals over the control horizon T h and S v is a linear operator extracting the WEC velocity from the state-space vector. The control force and state vector, however, are not independent variables, and are constrained by the system dynamics equation of the WEC, which in discrete time is defined as: To preserve mechanical and structural integrity, motion constraints and force limits are imposed. These constraints limit the maximum actuation force and the WEC device velocity and vertical displacement, i.e., where S p is a linear operator extracting the WEC displacement from the state vector. The cost function, together with the constraints, represents a linear MPC problem in its standard formulation. Let us define the following consolidated vectors based on the sequence of states and control commands. This will help in simplifying the notation and problem setup: Let us define Λ v , Λ p as block diagonal matrices having the velocity extraction matrix S v and the position extraction matrix Λ p , repeated N times respectively on the main block diagonal. Furthermore, let I be the N × N identity matrix and let ξ be a N × 1 vector of ones. In this manner, the cost function can then be expressed as follows: The inequality constraints can be reformulated using this vector notation as follows: where By recursively applying the discrete-time dynamics equations, it is possible to express X as a function of the control vector U, the excitation force vector F e , and the initial condition x 0 : where This allows us to rewrite the MPC problem as the following cost function and constraint equations: min The maximization of absorbed power requires the solution of a constrained convex optimization problem, for which well-consolidated routines, such as interior-point-convex or active-set methods, are available in literature [21,22]. Positive definiteness of the Hessian is, in general, guaranteed for the optimization of a point absorber device, unless the time step chosen for the conversion of the continuous time model into discrete time is too large to represent the actual behavior of the WEC device. At each timestep, an MPC problem needs to be solved, and the first value of the optimal solution vector U * is applied to the system. The system is marched forward in time by integrating the system dynamics using a standard Runge-Kutta (RK) scheme. The above MPC formulation also assumes that a perfect excitation force prediction is available 20 to 30 s into the future using up-wave measurement probes.

Causal Feedback Control
The first step in optimal causal feedback control design is to determine an analytically tractable stochastic disturbance model to represent the incident wave force f e (t). This step is nontrivial for two reasons. Firstly, the power spectral density typically used to characterize the free surface elevation η(t), denoted S η (ω), does not constitute a rational spectrum (i.e., it is not a ratio of rational polynomials of frequency ω). Secondly, the transfer function from the free surface elevation to the incident wave force, equal to involves the solution to a partial differential equation (due to fluid-structure interaction), and consequently it also does not evaluate to a rational transfer function. Consequently, the resultant power spectral density of the incident wave force, denoted as S e (ω) = S η (ω)|H e (jω)| 2 is an irrational function of frequency. To make the causal control design problem tractable, this spectrum must be approximated, as where the transfer function G(s) is a stable, minimum-phase, strictly-proper, rational transfer function. Equivalently, we seek to find matrices A e , B e , and C e such that in the above equation, With such a formulation for G(s), we may equivalently represent the incident wave force as filtered white noise; i.e., .
x e (t) = A e x e (t) + B e w(t) where w(t) is a white noise process with unit spectral intensity.
There are a number of ways to find appropriate matrices A e , B e , and C e for the noise filter model described above. Here, we make use of the subspace-based spectral factorization approach. Because this approach is applied exactly as described in [23], we will forgo the details here in the interest of brevity.
With the accomplishment of the finite-dimensional approximation of the incident wave force, we then append the internal states x e (t) to the state x(t) of the physical system; i.e., we define x a (t) = x T (t) x T e (t) T to arrive at an augmented state space model: .
x a (t) = A a x a (t) + B a u(t) + E a w(t) y(t) = C a x a (t) + D a u(t) where matrices A a , B a , E a , C a , and D a are appropriately defined. Our objective is to design a causal feedback law that maps y into u, to maximize the objective where E{.} denotes the expectation in stationarity, and R is a nonnegative parameter. For the case in which R = 0, the above expression is simply the mean power absorption. In design, R is typically chosen to be greater than zero, for two reasons. Firstly, it can be used to (approximately) quantify the power dissipation in the power train of the WEC, thus changing the performance objective from absorbed power to generated power. Secondly, it can be used as a tuning parameter to balance the power generation objective against the need to keep the mean-square control force magnitude below a desired threshold.
If there were no constraints on the response of the WEC, the above problem could be solved in closed form. More specifically, it distills to a sign-indefinite Linear Quadratic Gaussian (LQG) control problem, which is guaranteed to have a unique, stabilizing solution assuming that the open-loop mapping from u to . z is passive [6]. However, for the WEC under consideration here, there are constraints on both the displacement z and the control force u. In the presence of these constraints, the closed-form LQG control solution cannot be used without modification. To address these constraints, in this paper we implement the methodology outlined in detail in [9]. Because the methodology used here is identical to that one, we only provide an overview here of the basic steps in the nonlinear control design procedure.
The first step is to design a linear feedback controller of the form . ξ = A c ξ + B c y u = C c ξ in which A c , B c , and C c are optimized to maximize P a , subject to two constraints. The first of these constraints is a relaxed version of the WEC displacement constraint, in which we require only that the mean-square displacement be below a bound related to p max , i.e., where σ < 1 is a design parameter, chosen here to be 0.25. The second constraint is the restriction that the open-loop controller must be stable; i.e., that A c be a Hurwitz matrix. Convex techniques for accomplishing this constrained optimization are covered extensively in [9]. The resultant linear controller is guaranteed to result in positive power generation, be robust to saturations in the input u, and maintain mean-square displacements below p 2 max . However, it does not guarantee that the displacement limit is satisfied at all times. The second step is to augment the above linear control design with a nonlinear strokeprotection feedback loop that ensures that the displacement limit is satisfied at all times (rather than just in the mean-square sense). This stroke protection loop, and its relationship to the linear control design, is illustrated in Figure 1. It accepts as inputs the position and velocity, and outputs a supplemental nonlinear force q. This quantity q is sent back to the linear controller through an augmented input matrix, i.e., via the equations: where E c is a matrix of design parameters. The design of the feedback function that maps z, . z into q, and the design of E c , are described in fully in [9]. Here it suffices to say that q can be viewed as being similar to the restoring force of a "virtual hardening spring," providing no force when the displacement is far from its limit, but increasing to infinity when the displacement approaches its limit. The participation of velocity . z in the function can be thought of as providing damping. Meanwhile, the matrix E c is designed such that the open-loop transfer function from q to . z is passive, which guarantees unconditional stability of the closed-loop system. The controller resulting from the above design framework is guaranteed to protect the displacement from reaching its limit for the case in which the force u is unbounded. When saturation limits are imposed on u, of the form then the displacement protection is not strictly guaranteed. Indeed, when both the force and displacement are constrained, there may not exist a controller that can honor both. However, proper tuning of the various design parameters in the control design can be undertaken to assure that displacement constraints are satisfied in all but extremely rare cases in stationary response.

Point Absorber Limit
Consider a heaving axisymmetric body oscillating without constraints in resonance with an incoming regular wave of period T (s) and height H (m). Let J E (W/m) denote the wave power level, k w (1/m) denote the wavenumber, and λ (m) the wavelength. In this case, the point absorber effect in wave-energy extraction [24] imposes a fundamental limit on the average absorbed power P a . As described in [24,25], this limit can be derived from the relationship between absorption width d a = P a /J E (m) and wavelength λ, as shown below: Or equivalently: The last equality in Equation (28) is valid for deep water conditions. This is referred to as the Point Absorber Limit. The theoretical limit may be reached if the average absorbed power equals half the average excitation power, which happens when the radiated power is equal to the absorbed power. Furthermore, it is known that the Point Absorber Limit can only be reached up to a certain wave period and height. This depends on the motion amplitude constraints on the absorber. Beyond, only a lower relative power absorption can be realized, as defined by a device's volumetric limit.

Volumetric Limit
The heaving body sweeps a finite volume during its oscillation cycle based on physical limits. This volumetric constraint imposes a limit on the average power absorbed from each oscillation cycle. As shown in [25], the average power P a in the motion-constrained case is limited by the expression: This is referred to as the volumetric limit. With advanced control, we expect to achieve the active theoretical limit based on the device's operational region.

Controller Design Trade-Off Space
The controller design trade-off space for a WEC is defined by the Point Absorber and Volumetric Limit. For any given wave period, we will call the lesser of the two the "active theoretical limit" or the "theoretical limit" for simplicity. The controller performance cannot exceed the theoretical limit, and this defines the design trade-off space for controller optimization. As shown in Figure 5, the Point Absorber Limit is "active" for shorter wave periods. We will call this set of wave periods Region I. The Volumetric Limit is "active" for longer wave periods and we will call this set of wave periods Region II. The cross-over period, where the Point Absorber and Volumetric Limits are equal, demarcates the two regions.
Equations (28) and (29) show that both upper limits are wave height dependent and describe the WEC device performance in sinusoidal waves. To understand theoretical upper limits for an irregular wave train, we compute the wave height for a sinusoidal wave with the equivalent average power density of the irregular wave train as follows. The power in a monochromatic wave of height H and period T is given by: The equivalent wave height (H) of a monochromatic wave with the same power as the polychromatic wave can be found by equating the two. Assuming the period of the monochromatic wave to be the same as the peak period of the polychromatic sea state, we obtain the wave height of the equivalent sinusoidal wave as: This approach allows us to classify regions in the scatter diagram based on point absorber theory, and to draw meaningful conclusions regarding the importance of MPCbased control. For this purpose, we divided the entire set of input wave conditions into three distinct regions based on theoretical limits. Figure 6 illustrates the partitioning of sea states, where Region I corresponds to wave conditions where performance is dominated by the Point Absorber Limit. Region II corresponds to wave conditions where performance is dominated by the Volumetric Limit. Finally, Region III is associated with large waves where controllers are expected to maintain rated power production and minimize structural loads on the system. Region III is not important from a performance optimization perspective, but it is important that a control law active in this region will continue to produce power at the rated capacity, while minimizing structural loads. Although slightly counter-intuitive, structural loads in Region III can be smaller than those in Regions I and II because we are effectively de-tuning the system to the incident wave conditions. Applying the above region-based approach to a small and a large heaving point absorber provides the following regions. To understand their significance in context, we simply superimposed the regions onto the frequency distribution of the DoE's reference resource in Humboldt Bay, California. As shown in Figure 7, the reference geometry used for this paper has only a few sea states, where causal control has the potential to provide similar performance as that of the MPC-based controller. This is in stark contrast to the SANDIA Wavebot, which has a much larger volumetric displacement and as a result performs well in most sea-states using causal control (see Figure 8).   The comparison in Table 2 shows a buoy with a much larger volumetric displacement, similar to that tested by SANDIA [26]. It shows that fewer regions are present where volumetric limits dominate and, therefore, MPC becomes less important.

Benchmarking Performance of Causal Control Design Methods
To provide a fundamental understanding of the unconstrained controls' performance, we compared Optimal Causal Control, which is based on the LQG paradigm [7], with SAN-DIA's causal controller designed based on the CCC approximation principle-Proportional Integral (PI) and Feedback Resonator (MPC-FBR). A set of MATLAB and Simulink files containing methods for implementing PI and FBR controllers is available online through the Marine Hydro-Kinetic Data Repository [26]. Both controllers were simulated to compare performance in different wave conditions using no wave prediction information. MPC was also simulated with the same inputs to serve as an upper limit for benchmarking performance of the two causal controller design methods.
The results of this initial trade-off study demonstrated that the algorithm provided by SANDIA provided significantly lower power capture across all sea states, but especially in longer-period waves, where the impedance-matching condition of the underlying controls law is physically impossible due to the device's volumetric constraints. Because we found the LQG method outperformed the CCC controller, the remainder of comparisons in this paper were carried out between the LQG causal control method and a linear MPC controller. Time-domain simulations were carried out using Optimal Causal Control using the LQG method outlined in [9] and linear Model Predictive Control for a complete set of wave conditions. A discrete time-step of 0.01 s and a simulation length of 3600 s was utilized. A Pierson-Moskowitz (PM) wave spectrum with the following formulation was chosen to generate the wave train using the WAFO Toolbox [27]. A plot of the input PM spectrum is shown in Figure 9 S where ω p = 2π T p and ω n = ω ω p .

Controller Performance with Unconstrained PTO Force
First, let us look at the unconstrained case where no PTO force limits are imposed on the controller. Motion constraints are still enforced because the displacement of the device relative to the equilibrium is limited to ±2 m. Figure 10 shows a plot of the ratio of average absorbed power using causal control normalized to the average absorbed power using MPC. For short period waves, the ratio of causal controller response closely follows that of MPC and the performance ratio is within 90%. This corresponds to Region I, where power absorption is dominated by the Point Absorber Limit. In long-period waves, the performance ratio starts to deteriorate for wave periods >10 s. This corresponds to Region II, where power absorption is dominated by the volumetric constraint. The results show that, unlike in the case of MPC, the stroke-limited power absorption capability of causal controllers is sub-optimal.

Comparison of Annual Average Energy Capture for Reference Site
A simulation study was carried out to benchmark performance of the two control approaches by calculating the overall annual average energy capture for a deep-water reference site. The heaving buoy WEC was simulated for all of the wave conditions in the scatter diagram shown in Figure 7 for the reference site at Humboldt Bay, CA. Both controllers were simulated with a stroke limit of ±2 m, and no limits on the peak PTO force were imposed.
Assuming a plant capacity factor of 30%, we obtained rated power of 473 kW using MPC and 375 kW with causal control (see Table 3). At this site, the causal controller captures only 79% of the MPC-based controller. Region II 81% 82% Figure 10 shows the ratio of average absorbed power between causal control and MPC at each sea state. Note that in the regions which are point absorber limited (Region I), the two controllers perform equally well. For long-period waves where performance is upper bound by the Volumetric Limit (Region II), as expected, MPC outperforms causal control. Overall, we note that the percentage of sea state occurrences in Region II is 67.8% and is responsible for close to 81% of the respective annual energy captured from all sea states.

Controller Performance with Constrained PTO Force
Controllers must often operate with limits imposed on the maximum PTO force due to hardware limitations. Therefore, it is important to determine how well the two control approaches handle force constraints. In this situation, the controller must attempt to maximize power absorption while simultaneously limiting the maximum forces used on the PTO. Not all controllers can achieve this dual objective effectively, which typically results in sub-optimal performance for causal controllers. Figure 15 shows a plot of average power capture vs. force constraint for both control approaches. Note that for the same choice of PTO force limit, the average power capture with MPC significantly exceeds causal control. Another way of interpreting these results is that MPC requires less PTO force to achieve the same power capture performance as that of causal control. Figures 16-19 show time domain comparison of the device response with causal control and MPC when the PTO force is limited to 0.5 MN.     In addition to the control commands issued by the linear causal controller, the causal controller also has a non-linear stroke protection block which adds to the requirements of the PTO machinery force. Causal controllers require careful tuning of the non-linear control parameters to reduce spikes in PTO machinery force. Tuning the stroke protection parameters helps to smooth the PTO force response but it comes at the cost of sacrificing the amount of power that is absorbed. Figure 17 shows a comparison of the PTO force response with causal control and MPC. Note that MPC is able to limit the PTO force effectively, and does not violate the force limit. In contrast, causal control occasionally violates the force limit.

Conclusions
In this paper, we presented results of a performance comparison study between causal and non-causal methods in a design trade-off space formed by theoretical upper limits on power absorption. While exploring this design trade-off space, we show that identification of operating regions based on device displaced volume and site-specific annual spectrum data are key drivers in identifying suitable control methods for a given WEC application. MPC and causal control performance is similar for Region I, where performance is constrained by Point Absorber Limits and assuming that the PTO is capable of full four-quadrant control. For operating Region II, where performance is limited by Volumetric Limits, MPC shows a significant performance advantage compared to the causal control method.
It should be noted that the cost of most in-ocean systems is directly proportional to their displaced volume. This means that the ratio of average power to volume (P/V) is a key indicator of structural efficiency of a WEC device. Small point absorbers tend to have better P/V ratios, which drives economically viable solutions towards smaller devices. As shown in this study, these smaller point absorbers will produce power predominantly in Region II, where volumetric constraints dominate and MPC provides significant performance advantages. Furthermore, constrained optimal control offered by MPC provides benefits beyond pure performance gains. The ability to directly handle a wide range of constraints as part of the optimization problem ensures constraints can be gently applied, thus reducing transient loads that are harmful to the PTO and device structural components, while maximizing power production.