Practical Optimal Control of a Wave-Energy Converter in Regular Wave Environments

The governing equation of motion of a single point absorber in the time domain, also known as the Cummings equation [18], can be written as

where $\ddot{z}(t)$ and $\dot{z}(t)$ respectively represent the acceleration $(\text{m}/\text{s}^{2})$ and velocity $(\text{m}/\text{s})$ of the body at time $t$ in the heave mode, $M$ (kg) is the mass of the body, $m$ (kg) is the added mass (representing the mass of fluid displaced by the body during operation), and $K(\cdot)$ $(\text{kg}/\text{s}^{2})$ is the memory function of the radiation force that represents the hydrodynamic damping of the system. The forces considered for the current analysis include the excitation force on the body, $F_{\text{ex}}$ (N); the hydrostatic force, $F_{\text{h}}$ (N); the power take-off (PTO) force, $F_{\text{PTO}}$ (N); and, the end stop force $F_{\text{es}}$ (N) ensuring some stopping conditions. We assume linear wave theory and ignore viscous drag on the system [19].

The excitation force can be written as

where $f_{i}$ ($1/\text{s}$) are wave frequencies, $\phi_{i}$ (rad) are random phases, $\bar{F}_{\text{ex}}$ ($\text{N}/\text{m}$) represent complex vectors of wave excitation force per meter of wave amplitude in the frequency domain, $\Delta f$ ($1/\text{s}$) is an appropriate frequency step, and $S(f)$ ($\text{m}^{2}/\text{Hz}$) represents the spectrum of the incident wave field [20]. For single regular wave conditions investigated in this paper, the excitation force $F_{\text{ex}}$ is linear in terms of $w(t)$.

The hydrostatic force can be written as

where $\rho$ ($\text{kg}/\text{m}^{3}$) is the density of the water, $g$ ($\text{m}/\text{s}^{2}$) is the acceleration due to gravity and $S_{w}$ ($\text{m}^{2}$) is the cross section area of the WEC at the water plane. Since the parameters in equation (3) are constant, it can be rewritten as

Due to the mechanical complexity of the PTO system, it is common to limit operating stroke length of the WEC by introducing physical bounds, known as end stops [21]. Numerically, the end stops are modeled as a restoring force with large restoring coefficient [8]. The end stop force activates only when the relative motion between the WEC bodies exceeds the limits; for the current simulations, a limit of $-3$ m to $3$ m is used. The end stop force acting on the WEC is given as

where $K_{es}$ ($\text{N}/\text{m}$) is the restoring coefficient, $z_{es}$ (m) is the stroke limit and $u$ is the step function. In most operating conditions, the stroke length is within the limit.

In the present analysis, the point absorber is attached to the seabed with a long mooring line so that the mooring force acting on the system is negligible. The power take-off (PTO) force is ideal in the sense that it can exert any desired force within some limits.

We estimate the system coefficients for the control system from a numerical simulation of the point absorber. In our numerical simulation, the hydrodynamic coefficients, such as added mass, hydrodynamic damping coefficient and excitation forces, are calculated in the frequency domain using WAMIT [22], which is a potential flow solver based on a boundary element method (BEM). The linear wave theory assumption is adopted in the BEM solver, i.e., the wave amplitude is assumed to be very small compared with the wavelength. Therefore, the body-boundary conditions and free surface are linearized, and the gradient of velocity potential, $\Phi$, satisfies the Laplace equation ($\nabla^{2}\Phi=0$) in the fluid domain. The added-mass and damping coefficients calculated using velocity potential are given by

where $m_{ij}$ is the added-mass coefficient and $\zeta_{ij}$ is the hydrodynamic damping coefficient.

The excitation force ($X_{i}$) on the WEC due to the incident wave potential ($\phi_{0}$) is given by

The multi-body dynamics of the WEC is modeled using WEC-Sim [23], which is a MATLAB-based, open-source code. In WEC-Sim, the complex interaction between the incident waves, device motion, and power take-off mechanism is modeled in the time domain. The power performance and device motion of the WEC are simulated using the radiation and diffraction method. In the present analysis, the hydrodynamic bodies (float and torus) of the WEC are surface meshed and constrained to operate in one degree of freedom (heave, in the present case). With the hydrodynamic coefficients calculated from WAMIT, WEC-Sim calculates the time-domain motion and power performance of the WEC. In Section 3.6, we describe how we estimated the system coefficients, including the hydrodynamic system, control, and wave coefficients.

The objective of WEC control is the maximization of the total energy extracted throughout a given horizon $T$, which can be simply described as

where $u(t)$ is the control force at time $t$.In the context of the governing equation described in the previous section, $u(t)$ represents the PTO force at time $t$. Based on (9), the control force is considered to have a linear relationship with the states of the system.

Applying a forward Euler discretization to the system dynamics in (9), we have the following discrete-time optimal control problem:

subject to:

where $\Delta t_{k}$ is the discretization interval for the $k$th time period, $k\in\{0,1,\dots,N\}$. In the discretized system dynamics represented by (17), $A\in\mathbb{R}^{2\times 2}$ is the coefficient matrix for the discretized system involving the velocity and the position of the device, $b\in\mathbb{R}^{2}$ is the coefficient vector for the control force $u(t)$, and $c\in\mathbb{R}^{2}$ is the coefficient vector for the wave position $w(t)$. Equation (12) enforces the initial condition of the system given an initial point, whereas (18) limits the control force applied into the system by $\gamma$, denoting the maximum physical limit on the force, in Newtons. We refer to this constraint as a “hard bound.”

Observe that $\Delta t_{k}$ can be different for different $k$, based on the specific discretization method used. In the formulation above, the objective in (11) is approximated via the trapezoidal method. If an equidistant discretization method is considered, the objective in (11) would be given by

which is common in the literature; see, e.g., [11, 15].

The optimal solution to the formulation above may involve control forces that are undesirable from the perspective of the safety and longevity of the device. One concern is that the force might alternate quickly between a large positive quantity and a large negative quantity, creating a hammering effect that can cause fatigue over the long term. Another concern is that the force might stay at a large positive or a large negative quantity for a long time and apply constant pressure to the mechanical components, especially to the joint parts of the device. This, too, may result in undue wear-and-tear on the device.

A solution to these potential causes of damage is to impose a “soft bound” on the control force, wherein one penalizes the objective function for any control value that exceeds a certain fraction of the hard bound constraint (18). In particular, we incorporate into the objective the term

where $\rho\geq 0$ is the per unit penalty for exceeding the soft bound, $(\cdot)^{+}=\max\{0,\cdot\}$, and $\eta\in(0,1]$ is the fraction of the hard bound at which to impose the soft bound.

Finally, it is very common to have a physical limit on the displacement of the torus relative to the float. As is done in several works [14, 11, 16, 12], we consider a bound on the maximum displacement of the device (which is in fact the relative displacement between the torus and the float), given as

where $\delta$ is the constant of displacement in meters. Another interpretation of (21) is that it prevents the torus from going underwater completely, or rising above the sea level and slamming onto the sea surface.

Many WEC control models in the literature do not account for the cost or energy required to apply the control. In other words, they maximize the energy captured by the WEC using formulations similar to (11)–(18), but they ignore the energy that must be expended to exert this control force. As a result, the net energy captured by the device may be much less than predicted by the model—it may even be negative.

As a simple example, suppose that there is no incident wave, i.e., $F_{\text{ex}}(t)\equiv 0$. For certain values of the remaining parameters, and with initial conditions, the optimal solution to (11)–(18) (plotted in Figure 2) has non-zero control force and a strictly positive optimal objective function value $\approx 10^{5}$ (i.e., a strictly positive energy extraction). This suggests that we can extract energy even if there is no wave. Of course, this is erroneous, and the discrepancy comes from the model’s failure to account for the energy expended to exert the control force. In this simple example, the motion of the WEC comes solely from the control force (since there is no incident wave), and the energy required to exert that force must exceed the energy extracted from the motion. Hence, this solution involves a net loss of energy, even though the model suggests a gain. Any optimization model that does not incorporate the cost of control will suffer from this same sort of error.

Although the “cost” of the control force has not been included explicitly in WEC control models, one could argue that it is accounted for implicitly in some control strategies. For example, relatively simpler control methods such as latching control [24] or phase control [25] restrict the control policy to have a certain form, rather than optimizing over all possible controllers; since these control methods are common, one might infer that they do not incur excessive control costs. However, we argue that these control methods are too simplistic, since they do not promote device safety and offer less flexibility to be extended to irregular wave environments. Therefore, advanced control models such as the one proposed here are more promising. For such models, including the cost of control as an explicit term in the model is essential.

Of course, it may be difficult to derive an explicit expression for the control cost. Such expressions may be highly complex, and may depend on the specific power conversion system. Therefore, as a proxy, we assume that the instantaneous cost of applying a force $u$ has the functional form $\mathcal{O}(u^{2})$, i.e., that it equals $\lambda_{1}u^{2}$ for a constant $\lambda_{1}$. We chose a form by calibrating various nonlinear models using data from [26], which considers magnetic damping as a controller for a specific system. More details on our derivation of this functional form are given in Appendix A. Of course, for systems that use different PTO systems and control strategies than those in [26], the control cost may have a form that is different from $\mathcal{O}(u^{2})$. For such systems, one can use the approach we outline in Appendix A, or another approach, to estimate the functional form, and then insert that form into our optimization model. That is, our framework is flexible enough to handle any form of the instantaneous control cost, and our contribution lies in the modeling and solution approach, rather than in the specific form used for this cost.

The optimal control strategy resulting from the basic formulation in Section 3.1 is likely to be non-smooth since any arbitrary control policy satisfying the bounds is allowed. For example, in bang-bang control policies, the control is always at its bounds. Some works in the literature, e.g., [27], try to exploit this type of structure and develop efficient procedures for finding the optimal solution within the family of bang-bang policies. However, control policies, like bang-bang, that have sharp jumps in the control might not be desirable for two main reasons:

1. 1.

   It may not be possible to apply a non-smooth control policy, since most power take-off systems are best suited for smooth control forces.

2. 2.

   Non-smooth control would create additional wear and tear on the components of the WEC.

Thus, we would like the control sequence to be smooth in our formulation. One straightforward way to encourage smoothness is to add the term

to the objective function, where $\lambda_{2}$ is a parameter to adjust the importance of smoothness. By approximating $\Delta u(t_{k})$ from the differences in $u(t_{k})$ and applying the trapezoidal rule, one obtains the term

From another perspective, this term can also be considered as the cost of changing the control force, with a cost coefficient of $\lambda_{2}$.

Combining all the elements discussed in Sections 3.1–3.4, we have the complete formulation of the discrete-time optimal control problem:

subject to

The objective function (23) maximizes the total energy extraction minus the cost of applying the control force and penalties for violations of the soft bounds on the control limits. Constraints (24) and (25) are used in conjunction with the second term in the objective as an equivalent linear representation of the soft bound condition (20). Specifically, the auxiliary variables $\alpha_{k}$ describe the amount by which the force exceeds the soft bound, which are penalized in the objective with a constant per unit cost $\rho$.

We make the following remarks about this complete model.

• 1.

  It is assumed that the unit cost violation term $\rho$ is a constant. In order to achieve different bound structures, one can use a function $\rho(\cdot)$ instead of specifying only one cost coefficient $\rho$, where the function depends on the amount of violation.

• 2.

  Power smoothness of the operation can be handled internally in the mechanical/electrical power take-off component [28, 29].

• 3.

  The same mechanical/electrical component is used for both taking power from the system and feeding power into the system to apply the control.

• 4.

  Since the parameter $\lambda_{1}$ gives the cost of applying the control force, it also captures the conversion inefficiency of the system. For instance, if we assume a conversion efficiency ratio $0<\psi<1$, then we have a system that can only capture a proportion $\psi$ of the total energy injected into the system by the controller, and we can extract only a proportion $\psi$ of the total energy harvested. In this case, the inefficiency in the conversion system can be incorporated into $\lambda_{1}$ by properly adjusting this coefficient.

We simulate the wave–structure interaction of the WEC using WEC-Sim. The time-domain simulation of a wave–structure interaction is generated by the radiation and diffraction method. In this method, the hydrodynamic coefficients are calculated with a frequency-domain boundary element method and the system dynamics is solved in the time domain. We then estimate the coefficients $A$, $b$, and $c$ in the linear equation (17) by simulating the device under different wave conditions and control policies. Specifically, we use the following three-step procedure. The estimations in each step are done using least-squares regression.

1. 1.

   First, we consider a free-decay experiment in which the device is given an initial condition with no wave input. From this experiment, we estimate the hydrodynamic coefficient matrix $A$.

2. 2.

   Second, we consider a free-floating experiment in which the device is floating with some wave input. Using the previously estimated coefficient matrix $A$, we estimate the wave coefficient vector $c$ in this step.

3. 3.

   Finally, we conduct an experiment with the same input wave and some control policy. Using the previously estimated coefficients $A$ and $c$, we estimate the control coefficient vector $b$. A caveat here is that we cannot simply use the control policy used in the simulation to estimate the coefficients, since the system representation in (17) is only an approximation. That said, we know that in regular waves, the control policy has a sinusoidal form with the same period as the wave input. Thus, to estimate $b$, we consider sinusoidal control policies with a fixed amplitude (as the force limit), with the same period as the wave input, and with all possible shifts within the discretization level; we then choose the coefficients that resulted in the least error in the regression.

Note that all the coefficients are estimated independently to ensure that they accurately represent the system. An advantage of this procedure is that the hydrodynamic coefficient matrix $A$ is estimated only once since the device remains the same. However, if a different wave input is considered, one needs to start over from step 2. If the new wave input has the same period, it is straightforward to scale the wave coefficients $c$ accordingly. However, it is not so clear how to adjust the coefficients when a wave input with a different period is considered. This problem may be a possible future direction.

Typical wave conditions for the east coast of the United States have an average wave period of 4–6 seconds and significant wave height around 6 meters [30]. We use these numbers to guide our choice of wave inputs in our estimation, which are summarized in Table 1. Note that these conditions are off-resonance for the device $(\text{wave period}\neq 10s)$.

The estimated coefficients for the system in which the input wave has 6-meter significant wave height and 4-second wave period are:

In order to visually inspect the dynamics of the hydrodynamic system using the coefficients from (26), we rolled forward the equations in (17). Figure 3 plots the resulting trajectories using using a sinusoidal control policy that respects the bounds on the control variables. (The control policy used has a similar form to the policy used in the third phase of the estimation). From the figure, it is clear that the system appears to be stable, with trajectories that conform to expectations. Moreover, by comparing the real and computed trajectories in Figure 4, it is clear that the estimated coefficients yield dynamics that accurately represent the system.

The coefficients for 6-meter significant wave height and 5-second wave period are:

Finally, the coefficients for 6-meter, 6-second waves are:

The values of the objective function coefficients $\lambda_{1}$ and $\lambda_{2}$, which are intended to capture the cost of applying and changing the control force, are difficult to quantify precisely, but it is possible to identify a range of plausible values for them. Note that if we set $\lambda_{1}$ or $\lambda_{2}$ to too small a value, then we are saying that the control force can be applied nearly for free; this will result in the control force altering the device motion dramatically, not following the regular period of the waves, as described in Section 3.3. For sufficiently large values of the coefficients, the control force is so expensive that we will not apply it at all, and the device will simply follow the wave. With a regular wave input, the velocity of the device should have the same period as the wave; therefore, our suggestion is to choose the smallest values of $\lambda_{1}$ and $\lambda_{2}$ that result in the period of the velocity being sufficiently close to that of the wave.

For these experiments, the optimization time horizon was set as 7 times the wave period, $\Delta t_{k}=0.01$, $\gamma=10^{6}$, $\delta=3$, and $\eta=1$ (meaning that the soft bound equals the hard bound). We then considered two sets of pairs $(\lambda_{1},\lambda_{2})$ as follows.

1. (a)

   $\lambda_{2}=0$ with $\lambda_{1}\in\{10^{-10},10^{-9},\dots,10^{-4}\}$.

2. (b)

   $\lambda_{1}=0$ with $\lambda_{2}\in\{10^{-10},10^{-9},\dots,10^{-4}\}$.

Using these pairs of values, we solved the optimization problem and calculated the average realized period of the velocity trajectory. Figure 5 shows the calculated average period of the velocity trajectory for a wave input with a period of 4 seconds and a significant wave height of 6 meters.

The figure bears out the suggestion above that the calculated velocity trajectory has a shorter period than the waves if the coefficients are small, and approaches the period of the waves as the coefficients increase. We choose the smallest $\lambda_{1}$ and $\lambda_{2}$ values such that the average period of the velocity is within $5\%$ of the wave period. For this instance, we chose $\lambda_{1}=10^{-7}$ and $\lambda_{2}=10^{-9}$.

As a sanity check, Figure 6 provides optimal objective values and energy absorption over different $(\lambda_{1},\lambda_{2})$ values. The figure shows that with larger cost coefficients, the energy absorption and the objective function go to 0, suggesting that the optimal control trajectory goes to 0, which makes sense for large $(\lambda_{1},\lambda_{2})$ values. Overall, these results help us to verify that $\lambda_{1}=10^{-7}$ and $\lambda_{2}=10^{-9}$ are large enough to keep the period of the device in sync with the wave, yet small enough that the device is capable of absorbing energy in the model.

When we conducted similar experiments considering the other wave inputs, we found that slightly different $(\lambda_{1},\lambda_{2})$ values satisfied our criterion. This is due to the fact that by changing the period of the wave input, the optimal control sequence has a different period, thus the cost for applying control might require a different coefficient to ensure that the period of the device and the wave are in sync. Since, in a realistic setting, the parameters $\lambda_{1}$ and $\lambda_{2}$ are likely to be fixed (not tunable for different wave conditions), as a final choice of $\lambda_{1}$ and $\lambda_{2}$ we choose the largest values among all trials with different wave inputs. This led us to choose $\lambda_{1}^{*}=10^{-6}$ and $\lambda_{2}^{*}=10^{-6}$ for our remaining experiments.

As a further sanity check, as well as to illustrate the sensitivity of the optimal solution with respect to the values of $\lambda_{1}$ and $\lambda_{2}$, we evaluate the error that arises from using incorrect cost coefficients in the model. In particular, suppose that the true cost coefficients are $\lambda^{*}=(\lambda_{1}^{*},\lambda_{2}^{*})=(10^{-6},10^{-6})$. We then define the following quantities:

In Figure 7, we plot the values

as a function of $\lambda_{1}$ with $\lambda_{2}=\lambda_{2}^{*}$ (top plot) and also as a function of $\lambda_{2}$ with $\lambda_{1}=\lambda_{1}^{*}$ (bottom plot). When $\lambda_{1}=\lambda_{1}^{*}$ (resp. $\lambda_{2}=\lambda_{2}^{*}$), the ratio is equal to 1 in the top plot (resp. bottom plot). Otherwise, the ratio shows the objective value lost by solving the problem with the “wrong” cost coefficients. We see from this experiment that it is arguably much worse to underestimate the cost coefficients than to overestimate them. In particular, if the problem is solved with a cost coefficient that is too large, then the plotted ratio is less than 1, but at least still positive. On the other hand, if the problem is solved with a cost coefficient that is too small, then the plotted ratio can be negative. These experiments support our inclusion of these costs in the objective function, and support our idea of choosing the largest coefficients among all trials.

In order to decide on appropriate values for $\eta$ (the soft bound ratio) and $\rho$ (the penalty term for exceeding the bound), we model and solve the optimization problem (23)–(25) under multiple values of $\eta$ and $\rho$. In particular, we set the problem horizon as one full wave period, $\Delta t_{k}=0.01$, $\gamma=10^{6}$, and $\delta=3$, and $\lambda$ values set as in the previous section. We tested all combinations of $\eta$ values in $\{\frac{1}{50},\frac{2}{50},\dots,\frac{49}{50},\frac{50}{50}\}$ and $\rho$ values in $\{\frac{1}{50},\frac{2}{50},\dots,\frac{49}{50},\frac{50}{50}\}$, resulting in a total of $50\times 50=2500$ instances. The whole procedure was relatively fast since we limit the horizon to one full wave period, and since IPOPT can solve the model (23)–(25) quickly. It is worth mentioning that all the instances are successfully solved to local optimality. Note that each of the instances is solved multiple times with different initial points, and in all of these cases we get the same optimal solution, suggesting that the local optimum may be globally optimal.

Figure 8 plots the optimal objective function value under different safety parameters. As can be seen from the figure, the optimal objective value decreases as $\rho$ increases and as $\eta$ decreases. Moreover, the optimal objective function value is flat for $\eta$ greater than $\approx 0.2$, since below that range, the optimal control trajectory is already bounded away from $\gamma$. This effect will be demonstrated in more detail in the next section. Also, it can be said that the effect of $\eta$ is more direct, and that it inherently affects the impact of $\rho$ on the objective function.

Of course, the choice of the safety parameters primarily depends on the application and is subject to the modeler’s preferences. We will further explore the impact of $\eta$ and $\rho$ on the optimal solution in the next section.

Figure 9 plots the total CPU time (user CPU + system CPU) returned by the solver IPOPT integrated in the modelling software AMPL. The total CPU time required to solve the problem is relatively small throughout the parameter grid, except when the fraction of soft bound ($\eta$) gets small and becomes an active constraint. In this region, the problem becomes dramatically harder as $\eta\rightarrow 0$, with CPU times roughly quadrupled compared to the other instances.

In this section, we use the receding horizon approach described in Section 4, using the previously determined $\lambda$ values, initial time $T_{0}=200$, horizon $T=\text{wave period}$, time step $\Delta t=0.01$, controller update horizon $T_{c}=\frac{1}{10}\times(\text{wave period})$, and number of receding periods $K=100$. We also test each of the following $\eta$ and $\rho$ values:

1. 1.

   $\eta=1$ (soft bound equals hard bound),

2. 2.

   $\eta=0.1$ and $\rho=1$,

3. 3.

   $\eta=0.1$ and $\rho=0.1$.

In order to ensure that the control sub-sequence resulting from each problem at the end of a receding period is connected to the next one, we set the initial point of the control trajectory in each receding period as the last entry of the implemented control trajectory obtained from the solution of the previous problem.