Optimal Control of Nonlinear, Nonautonomous, Energy Harvesting Systems Applied to Point Absorber Wave Energy Converters

1. Introduction

Optimal control theory provides a powerful framework for determining a dynamic system’s control law that extremizes an objective function [1]. It has been applied in many fields, including engineering, economics, and biology, to name a few. This study focuses on its application to energy maximization for a general class of second-order, nonlinear, time-varying dynamic systems.

A common yet challenging class of optimal control problems are the Linear Optimal Problems (LOPs), also known as singular problems. They occur when the inputs appear linearly in the system’s state equation and cost function, even though the state equations may be nonlinear in the states. This means that the control law cannot be formed directly by applying the stationary condition of optimality [2]. Several researchers have investigated the existence and solution of singular problems in optimal control [3,4,5,6,7]. Scardina [8] investigated the necessary conditions for optimality for singular problems using Pontryagin’s Minimum Principle (PMP). He exploited the fact that the value of the Hamiltonian, the central quantity in PMP that combines the system’s state, control inputs, and co-state variables, should be constant and reach the absolute minimum for an optimal control function while satisfying the system’s dynamics, co-state equations and any induced boundary conditions [6,7] to develop general solutions to singular optimal control problems.

Kelly et al. and Robbins [9,10,11] further explored the optimality conditions and introduced the generalized Legendre-Clebsch condition. It provides additional requirements for optimality, such as when the Hamiltonian is linear in the control input which often leads to a singularity. The generalized Legendre-Clebsch condition has been extensively studied in the context of optimal control, and its application has provided valuable insights into the existence and properties of optimal control solutions [11,12].

Creating an optimal solution for a nonlinear, nonautonomous system often requires finding both the control law and initial conditions such that the system states stay on the optimal path [13,14]. This is the approach taken in this paper. Limit cycles represent bounded periodic motion in a system, and their establishment is crucial for achieving sustained and efficient energy absorption for a point absorber wave energy converter subjected to periodic waves [15]. Due to the complexity of finding an energy-maximizing path of a nonlinear, nonautonomous system, sometimes it is acceptable to linearize the system first [16]. In contrast, the approach below exploits the model’s nonlinear, time-varying nature and applies to a wide range of dynamic systems.

In recent years, optimal control has gained significant attention in energy extraction applications, including designing and controlling wave energy converters (WECs). Wave energy converters are devices that transform water wave energy into valuable forms of energy, such as electricity, through a Power Take-Off system (PTO). The authors in [17] present a review and classification of wave energy converter technologies. This paper focuses on point absorber WECs, which are most commonly modeled using a linear approach based on Cummins equation [18] and are accurate for small motion about an equilibrium. However, energy capture can be increased for large motions where the model contains nonnegligible, nonlinear terms [19]. For heaving point absorbers, an important nonlinearity is the Froude-Krylov (FK) force [20] found by integrating the incident wave pressure over the buoy’s wetted surface area generating both buoyancy and hydrodynamic expressions.

In the pursuit of maximizing the energy output of Wave Energy Converters (WECs), diverse methodologies have been explored. It is widely acknowledged that the shape of a WEC plays a pivotal role in its dynamic response and, consequently, its energy harvesting potential. Therefore, various approaches have incorporated strategies for shape optimization, to increase energy extraction from ocean waves [21,22,23]. In parallel, strategies have focused on optimizing the geometry of WEC arrays [24]. Moreover, a WEC’s performance is highly dependent on the control law used. Thus, significant attention has been directed toward control system development. This domain has witnessed the proposal of diverse control strategies, including sliding mode control [25], numerical model predictive controls [26,27], and latching control techniques [28,29], each offering unique avenues for improving energy conversion efficiency.

Optimal strategies are often proposed with a comprehensive review by Ringwood et al. [30,31]. A well-studied energy-optimal strategy for point absorbers, whose dynamic response can be accurately modeled by a linear differential equation, is called complex conjugate control (CCC) [32]. Its origins date to the 1840s with Moritz von Jacobi’s maximum power transfer law for resistive circuits where the load and source resistance should be equal. This was later expanded to impedance matching, where the load should be equal to the circuit’s Thevinen equivalent source impedance for maximum power transfer [33]. Impedance matching is also used in many engineering fields, including the vibration community, to remove energy from vibrating structures [34,35]. This control strategy results in optimal energy extraction; however, it exaggerates the motion of the device, and it is restricted to linear WEC models. It should be noted that the resulting large motions are inconsistent with the model’s small motion assumption [20]. Therefore, this study develops an energy-maximizing control strategy for nonlinear, nonautonomous second-order periodic systems that applies to nonlinear WEC models as derived in [36].

Establishing limit cycles in WECs ensures the motion remains periodic and bounded, facilitating efficient energy conversion. The derivation of optimal control laws that lead to limit cycles has been investigated in different contexts. Wilson et al. [37] derived feedback linearizing optimal control law for nonlinear WECs, resulting in closed-loop limit cycles. In [38], an optimal control law was obtained by assuming a linear model for WECs exposed to regular waves, where they obtained an explicit expression for the optimal control force maximizing the absorbed energy. The control law was shown to be equivalent to the CCC in this specific linear case. However, they didn’t explore limit cycles to determine the optimal initial condition but instead relied on a switching function to keep the motion bounded.

This article builds upon the existing literature by presenting a general optimal control law for energy maximization in second-order, time-varying, nonlinear periodic dynamic systems. It extends the analysis of singular control problems and explores the existence of limit cycles to find the optimal initial conditions. The approach is applied using linear and nonlinear models motivated by recent WEC applications.

The remainder of the paper is organized as follows. Section 2 shows the derivation of the general optimal control law, while Section 3 applies the control law to both a linear and nonlinear WEC point absorber model. Section 4 provides some conclusions and ideas for further study.

2. General Optimal Control Law

After introducing the form of the second-order dynamic system, the optimal control law is derived, including examining the second-order optimality conditions.

2.1. Control Law Derivation

We consider the general, second-order, time-varying differential equation

$$M\ddot{\zeta }=F_e(\zeta ,\dot{\zeta },t)-F_c$$

(1)

where (1) is often used to model energy harvesting systems. M is the system’s total mass, $\zeta$ is position, and $F_e$ is the net external force acting on the system which could be either linear or nonlinear depending on the system explored. $F_e$ is exclusive of the externally applied control force, $u=F_c$.

In the remainder, we will develop a control law for u that maximizes energy extraction. For the WEC examples considered later, $F_e$ contains time-varying, measurable disturbances due to the waves exploited to extract energy.

A state variable form of (1), $\dot{\overrightarrow{x}}=\overrightarrow{f}(\overrightarrow{x},F_e\left(\overrightarrow{x}\right),u)$, is shown below

$$\begin{array}{cc}\hfill {\dot{x}}_1& =x_2\hfill \\ \hfill {\dot{x}}_2& ={\displaystyle \frac{F_e\left(\overrightarrow{x}\right)-u}{M}}\hfill \\ \hfill {\dot{x}}_3& =1\hfill \\ \hfill \end{array}$$

(2)

where the states are defined as $x_1=\zeta$, $x_2=\dot{\zeta }$ and $x_3=t$. Time is introduced as a state to manage the nonautonomous nature of (1). It is a classical technique used in control theory to use the mathematical tools developed for autonomous systems to analyze and solve nonautonomous systems [13].

Our optimal problem will be framed in the context of a free end-time problem. For a free end-time problem, the initial state $x\left(t_o\right)$ is given, but the final time $t_f$ is not fixed.

The energy E is the integral of power

$$E={\int }_{t_o}^{t_f}u{x}_{2}dt$$

(3)

In (3) we neglect any dissipating electric power terms, concentrating instead on the mechanical energy. This approach is deliberate, as our primary focus is on the mechanical aspects of the system rather than the power electronics. Thus our objective function is framed to maximize the mechanical energy extraction, or conversely, minimize its negative as shown below:

$$\begin{array}{cc}\hfill \underset{u\in [t_o,t_f]}{min}:& -E(x_2,u)\hfill \\ \hfill & \dot{\overrightarrow{x}}\left(t\right)=\overrightarrow{f}(\overrightarrow{x},F_e\left(\overrightarrow{x}\right),u)\hfill \\ \hfill & x\left(t_o\right)=x_o\hfill \end{array}$$

where

$$-E(x_2,u)=-{\int }_{t_o}^{t_f}u{x}_{2}dt={\int }_{t_o}^{t_f}\varphi (x_2,u)dt$$

(4)

and $\varphi (x_2,u)=-u{x}_2$.

The Hamiltonian is [1,13]

$$H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=\varphi (x_2,u)+{\overrightarrow{\lambda }}^T\overrightarrow{f}(\overrightarrow{x},F_e\left(\overrightarrow{x}\right),u)$$

(5)

where ${\overrightarrow{\lambda }}^T=\left[{\lambda }_1{\lambda }_2{\lambda }_3\right]$ is the adjoint vector also known as co-states.

Substituting (2) and our definition for $\varphi (x_2,u)$ into (5) gives the specific form of the Hamiltonian.

$$H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=-u{x}_2+{\lambda }_1{x}_2+{\displaystyle \frac{{\lambda }_2}{M}}(F_e\left(\overrightarrow{x}\right)-u)+{\lambda }_3$$

(6)

The optimal control, $u^{*}$, and the resulting optimal trajectory in both states and co-states must certain optimality conditions.

$$\begin{array}{cc}\hfill & {\partial }_{\lambda }H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=\dot{\overrightarrow{x}}\hfill \\ \hfill & {\partial }_{x}H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=-\dot{\overrightarrow{\lambda }}\hfill \\ \hfill & {\partial }_{u}H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=0\hfill \end{array}$$

(7)

The three equations obtained from ${\partial }_{\lambda }H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=\dot{\overrightarrow{x}}$ of (7) are a restatement of the system state equations of (2) and the co-state differential equations, obtained from ${\partial }_{x}H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=-\dot{\overrightarrow{\lambda }}$ of (7) are known as Euler-Lagrange equations. The expansion of ${\partial }_{x}H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=-\dot{\overrightarrow{\lambda }}$ yields

$$\begin{array}{cc}\hfill {\dot{\lambda }}_1& =-{\displaystyle \frac{{\lambda }_2}{M}}{\partial }_{x_1}{F}_e\left(\overrightarrow{x}\right)\hfill \\ \hfill {\dot{\lambda }}_2& =u-{\lambda }_1-{\displaystyle \frac{{\lambda }_2}{M}}{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)\hfill \\ \hfill {\dot{\lambda }}_3& =-{\displaystyle \frac{{\lambda }_2}{M}}{\partial }_{x_3}{F}_e\left(\overrightarrow{x}\right)\hfill \\ \hfill \end{array}$$

(8)

The stationary condition ${\partial }_{u}H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u)=0$ of (7) is given as

$$-x_2-{\displaystyle \frac{{\lambda }_2}{M}}=0$$

(9)

One more condition is necessary since a free final time problem is considered. This condition is known as the transversality condition which provide additional boundary conditions at the terminal time $t_f$.

$$H(\overrightarrow{x_f},\overrightarrow{{\lambda }_f},F_e\left(\overrightarrow{x}\right),u^{*})=0$$

(10)

Since H is linear in u, (9) does not directly express an optimal u it is called a singular problem which could be challenging to solve analytically. Fortunately, our equations could mathematically be manipulated to derive a closed-form solution for the control force u in terms of the states $\overrightarrow{x}$ and $F_e\left(\overrightarrow{x}\right)$, as shown in the remainder of this section.

The stationary condition of (9) relates $x_2$ and ${\lambda }_2$

$${\lambda }_2=-M{x}_2$$

(11)

which is differentiated and combined with ${\dot{\lambda }}_2$ of (8) to obtain:

$$M{\dot{x}}_2=-u+{\lambda }_1+{\displaystyle \frac{{\lambda }_2}{M}}{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)$$

(12)

Substituting ${\dot{x}}_2$ of (2) and (11) into (12) and simplifying, we obtain the algebraic equation for ${\lambda }_1$.

$${\lambda }_1=F_e\left(\overrightarrow{x}\right)+x_2{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)$$

(13)

After obtaining the algebraic expressions of the first two co-states, we can obtain the algebraic expression of the third co-state ${\lambda }_3$. Although ${\lambda }_3$ is not used explicitly in the optimal control law, it is important to obtain all the co-state expressions to prove optimality later by showing that they are satisfied when applying the assumed optimal u. The value of the Hamiltonian of the optimal control function is constant across time and has the absolute minimum [8]. This could be further seen mathematically by taking the derivative of the Hamiltonian with respect to time. Since we have transformed our nonautonomous system to an autonomous system, time is now treated as a state [13].

$${\displaystyle \frac{dH}{dt}}={\partial }_{x}H{\partial }_{t}x+{\partial }_{\lambda }H{\partial }_t\lambda +{\partial }_{u}H{\partial }_{t}u$$

(14)

Substituting (7) in (14) yields to ${\displaystyle \frac{dH}{dt}}=0$. Thus, let the Hamiltonian of (6) equal an arbitrary constant $H_o$ as shown below

$$H(\overrightarrow{x},\overrightarrow{\lambda },F_e\left(\overrightarrow{x}\right),u^{*})=H_o$$

(15)

Substituting the algebraic equations of the two other co-states of (11) and (13) into (15), we get (16).

$$x_2^2{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)+{\lambda }_3=H_o.$$

(16)

If $u^{*}$ is the optimal control and $\overrightarrow{x}$ is the vector of corresponding optimal states, then there exists a co-state such that the Hamiltonian $H=0$ [7]. Moreover, this satisfies our transversality condition of (10). Thus, imposing $H_o=0$ , yields the algebraic expression of ${\lambda }_3$ as follows:

$${\lambda }_3=-x_2^2{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)$$

(17)

The algebraic expressions of the first two co-states can now be combined with (2) and (8) to find an algebraic equation for $u^{*}$, which is the optimal control force. Differentiating (13) with respect to time ($x_3$), we obtain:

$$\begin{array}{cc}\hfill {\dot{\lambda }}_1=& {\dot{x}}_1{\partial }_{x_1}{F}_e\left(\overrightarrow{x}\right)+2{\dot{x}}_2{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)\hfill \\ \hfill & +{\dot{x}}_3{\partial }_{x_3}{F}_e\left(\overrightarrow{x}\right)+x_2\left({\dot{x}}_1{\partial }_{x_1}{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)\right.\hfill \\ \hfill & \left.+{\dot{x}}_2{\partial }_{x_2}{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)+{\dot{x}}_3{\partial }_{x_3}{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)\right)\hfill \end{array}$$

(18)

Substituting $\dot{{\lambda }_1}$ of (8) and $\overrightarrow{\dot{x}}$ of (2) in (18) and simplifying, the optimal control force becomes:

$$u^{*}=F_e\left(\overrightarrow{x}\right)+M\left({\displaystyle \frac{x_2^2{\partial }_{x_1}{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)+x_2{\partial }_{x_3}{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)+{\partial }_{x_3}{F}_e\left(\overrightarrow{x}\right)}{2{\partial }_{x_2}{F}_e\left(\overrightarrow{x}\right)+x_2{\partial }_{x_2}^2{F}_e\left(\overrightarrow{x}\right)}}\right)$$

(19)

which can also be written in terms of the state definitions as:

$$F_c^{*}=F_e(\zeta ,\dot{\zeta },t)+M\left({\displaystyle \frac{{\dot{\zeta }}^2{\partial }_{\zeta }{\partial }_{\dot{\zeta }}{F}_e(\zeta ,\dot{\zeta },t)+\dot{\zeta }{\partial }_t{\partial }_{\dot{\zeta }}{F}_e(\zeta ,\dot{\zeta },t)+{\partial }_t{F}_e(\zeta ,\dot{\zeta },t)}{2{\partial }_{\dot{\zeta }}{F}_e(\zeta ,\dot{\zeta },t)+\dot{\zeta }{\partial }_{\dot{\zeta }}^2{F}_e(\zeta ,\dot{\zeta },t)}}\right)$$

(20)

Interestingly, the control law of (20) still applies for varying mass systems. Consider a system of a varying mass $M(\zeta ,\dot{\zeta },t)$ such that its governing equation is:

$$M\ddot{\zeta }+{\displaystyle \frac{dM}{dt}}\dot{\zeta }=F_e(\zeta ,\dot{\zeta },t)-u$$

(21)

Writing (21) in the form of (1), we get

$$\ddot{\zeta }=F_m(\zeta ,\dot{\zeta },t)$$

(22)

where

$$F_m={\displaystyle \frac{Fe-\frac{dM}{dt}\dot{\zeta }}{M}}$$

(23)

Applying the control law of (20) such that $M=1$ and $F_e=F_m$, we obtain the optimal control law of a varying mass system as follows:

$$F_c^{*}=F_m(\zeta ,\dot{\zeta },t)+{\displaystyle \frac{{\dot{\zeta }}^2{\partial }_{\zeta }{\partial }_{\dot{\zeta }}{F}_m(\zeta ,\dot{\zeta },t)+\dot{\zeta }{\partial }_t{\partial }_{\dot{\zeta }}{F}_m(\zeta ,\dot{\zeta },t)+{\partial }_t{F}_m(\zeta ,\dot{\zeta },t)}{2{\partial }_{\dot{\zeta }}{F}_m(\zeta ,\dot{\zeta },t)+\dot{\zeta }{\partial }_{\dot{\zeta }}^2{F}_m(\zeta ,\dot{\zeta },t)}}$$

(24)

which has an identical form to (20).

2.2. Second Order Optimality Condition

In addition to satisfying the first-order Euler-Lagrange equations of (7), Kelly et al. and Robbins [9,10,11] showed the existence of another necessary condition known as the generalized Legendre-Clebsch condition when the optimal solution satisfies (25).

$${(-1)}^k{\partial }_u\left[{\displaystyle \frac{d^{2k}}{d{t}^{2k}}}{\partial }_{u}H(\overrightarrow{x},\overrightarrow{\lambda },u)\right]\ge 0,k=1,2,3\dots$$

(25)

where k denotes the number of equations used for multiple control inputs. For our single input case, $k=1$. Writing t in its state form, $x_3$, (25) can be simplified to (26).

$${\partial }_u\left[{\displaystyle \frac{d^2}{d{x_3}^2}}{\partial }_{u}H(\overrightarrow{x},\overrightarrow{\lambda },u)\right]\le 0$$

(26)

Applying the test of (26) to the expression of (9) and doing the necessary substitutions and simplifications, yields:

$${\partial }_u\left[{\displaystyle \frac{d^2}{d{t}^2}}{\partial }_{u}H\right]={\displaystyle \frac{2{\partial }_{\dot{\zeta }}{F}_e\left(\overrightarrow{x}\right)+\dot{\zeta }{\partial }_{\dot{\zeta }}^2{F}_e\left(\overrightarrow{x}\right)}{M^2}}$$

(27)

According to (27), the Legendre-Clebsch condition applies when $2{\partial }_{\dot{\zeta }}{F}_e(\zeta ,\dot{\zeta },t)+\dot{\zeta }{\partial }_{\dot{\zeta }}^2{F}_e(\zeta ,\dot{\zeta },t)\le 0$. The sign of this term depends on the presence of $\dot{\zeta }$ in $F_e(\zeta ,\dot{\zeta },t)$ as well as its value and thus on the specifics of the application. For some optimal control problems, the Legendre-Clebsch condition is not satisfied for fixed time intervals. The control law may still be optimal over regions where it applies. Modifying the initial conditions can induce distinct motion patterns within the system, potentially guiding it toward a bounded optimal trajectory. This condition will be explored in detail for the control law derived in Section 3 for point absorber wave energy converters.

3. Application to Point Absorber Wave Energy Converters (WEC)

The control law of (24) is applied to two, time-varying, WEC model cases, one linear and the other nonlinear. Moreover, our solution is assumed periodic since a periodic input disturbance is considered [39], the waves, thus our boundary conditions will be written as follows:

$$\begin{array}{cc}\hfill x\left(t_o\right)& =x_o\hfill \\ \hfill x\left(t_f\right)& =x(t_o+nT)\hfill \\ \hfill \end{array}$$

(28)

where $x_o$ is to be determined and n is the number of periods.

For the linear case, (24) reduces to a well-known optimal solution along with initial condition requirements. The nonlinear case is more interesting, resulting in a nonlinear, time-varying control law. Initial condition regions that produce the required stable limit cycle motion are also explored. Since the system is time-varying, Poincare maps are used to ensure the motion is quasi-periodic.

3.1. Dynamic Model

Point absorber WECs exhibiting small motions about an equilibrium are often modeled using Cummins equation shown in (29), for the heave-only case [18]

$$M\ddot{\zeta }+{\int }_{-\infty }^t{h}_r(t-\tau )\dot{\zeta }\left(\tau \right)d\tau +k\zeta =f_f$$

(29)

where $\zeta$ is the buoy displacement relative to its equilibrium, and the mass M includes physical and added mass components, $M=m+m_a$. The added mass, $m_a$, is due to the water displaced by the oscillating buoy. The convolution term is the radiation damping, where $h_r$ is the radiation impulse response function. The linear buoyancy force coefficient is denoted as k, and $f_f$ contains all additional external forces acting on the buoy, such as the wave excitation and PTO forces. For regular waves, the convolution term becomes a single, constant coefficient b, and a special case of Cummin’s equation is given in (30).

$$M\ddot{\zeta }+b\dot{\zeta }+k\zeta =f_f$$

(30)

Next, we will generalize (30) to include nonlinear hydrodynamic terms while also showing the PTO force, u, explicitly

$$M\ddot{\zeta }+b\dot{\zeta }=F_h(\zeta ,t)-F_c$$

(31)

where $F_h(\zeta ,t)$ can be any continuous nonlinear, time-varying function, for example, due to nonlinear Froude-Krylov force expressions for non-cylindrical buoys undergoing large motions [40]. Comparing (31) to (1) we note that $F_e(\zeta ,\dot{\zeta },t)$ is

$$F_e(\zeta ,\dot{\zeta },t)=F_h(\zeta ,t)-b\dot{\zeta }$$

(32)

While we’ve kept the linear damping expression from Cummin’s Equation for regular waves, $b\dot{\zeta }$ can be replaced by a more general, perhaps nonlinear, damping expression.

When examining simulation results, we’ll use the physical parameters shown in Table 1 from [36] for the linear and nonlinear examples below. Though the optimal control law of (20) applies to any continuous, nonlinear, time-varying second-order model, however, to ease the comparison with previous literature [37], we will restrict our attention to regular waves with elevation

$$\eta \left(t\right)=Asin\left(\omega t\right)$$

(33)

where A and $\omega$ are the wave amplitude and frequency, respectively.

3.2. Optimal Control Law

Applying the optimal solution of Equation (20) for the $F_e\left(\overrightarrow{x}\right)$ of Equation (32) gives the optimal control law of Equation (34).

$$F_c^{*}=F_h(\zeta ,t)-b\dot{\zeta }-{\displaystyle \frac{M}{2b}}{\partial }_t{F}_h(\zeta ,t)$$

(34)

where the co-states of Equation (13), 11 and 17 become

$$\begin{array}{cc}\hfill {\lambda }_1& =F_h(\zeta ,t)-2b\dot{\zeta }\hfill \\ \hfill {\lambda }_2& =-M\dot{\zeta }\hfill \\ \hfill {\lambda }_3& =b{\dot{\zeta }}^2\hfill \\ \hfill \end{array}$$

(35)

The second-order optimality condition, Equation (27), is shown in Equation (36).

$${\displaystyle \frac{2{\partial }_{\dot{\zeta }}{F}_e(\zeta ,\dot{\zeta },t)+\dot{\zeta }{\partial }_{\dot{\zeta }}^2{F}_e(\zeta ,\dot{\zeta },t)}{M^2}}=-{\displaystyle \frac{2b}{M^2}}\le 0$$

(36)

which is always satisfied since both b and M are positive. The optimal state trajectories can now be derived using the control law of Equation (34). Substituting Equation (34) into Equation (31) yields the optimal acceleration

$${\ddot{\zeta }}^{*}={\displaystyle \frac{{\partial }_t{F}_h}{2b}}$$

(37)

Equation (37) can be integrated in closed form, if $F_h$ permits, or numerically to obtain the optimal velocity, displacement, and co-state trajectories by substitution into Equation (35).

3.3. Linear Point Absorber WEC

Consider the model of a linear point absorber below that has been used previously to derive the complex conjugate control (CCC) law [41]

$$M\ddot{\zeta }+b\dot{\zeta }+k\zeta =F\left(t\right)-F_c$$

(38)

where $F\left(t\right)$ is the excitation force. To apply the optimal control law of (34), we only need to identify the correct expression for $F_h$, and in this case

$$F_h(\zeta ,t)=F\left(t\right)-k\zeta$$

(39)

which results in

$$F_c^{*}=F\left(t\right)-b\dot{\zeta }-k\zeta -{\displaystyle \frac{M}{2b}}{\partial }_{t}F\left(t\right)$$

(40)

This result was previously obtained by Zou et al. using Euler-Lagrange equations applied directly to a linear point absorber model. [38] which was in accordance with the results obtained using CCC [41].

Figure 1. Point absorber WEC with a cylindrical shaped buoy.

Figure 1. Point absorber WEC with a cylindrical shaped buoy.

However, one often-cited practical limitation of CCC is that the buoy motion becomes larger than the model assumptions used in its derivation [25,31,32,42]. Thus, switching functions were used to keep the linear system bounded. In our work, we have managed to achieve that without the need for a switching function but by obtaining the proper initial condition of (28) to drive our system to a stable limit cycle. To examine this, we will consider the optimality requirement that the buoy’s motion at least resides on a stable limit cycle.

Using (37) the optimal acceleration is

$${\ddot{\zeta }}^{*}={\displaystyle \frac{\dot{F}\left(t\right)}{2b}}$$

(41)

which can be integrated twice to obtain its velocity and position

$$\begin{array}{cc}\hfill {\dot{\zeta }}^{*}& ={\displaystyle \frac{F\left(t\right)}{2b}}+C_1\hfill \\ \hfill {\zeta }^{*}& =\int {\displaystyle \frac{F\left(t\right)}{2b}}dt+C_{1}t+C_0\hfill \end{array}$$

(42)

where $C_1$ and $C_0$ are integration factors.

To achieve a stable limit cycle, bounded motion, $C_1=0$, resulting in the optimal buoy velocity of (43), the same result obtained from CCC [19].

$${\dot{\zeta }}^{*}={\displaystyle \frac{F\left(t\right)}{2b}}$$

(43)

Achieving $C_1=0$ depends on the excitation force expression, $F\left(t\right)$, and the buoy initial conditions defined when the control law is activated. To examine this, consider a particular wave excitation force expression, $F\left(t\right)$, used for a cylindrical point absorber buoy in [37].

$$F\left(t\right)=k\eta \left(t\right)=kAsin\left(\omega t\right)$$

(44)

and its corresponding model

$$M\ddot{\zeta }+b\dot{\zeta }+k\left[\zeta -\eta \left(t\right)\right]=-F_c$$

(45)

where $\eta \left(t\right)$ is the wave elevation.

From (43), at time $t_o=0$,

$$C_1={\dot{\zeta }}^{*}\left(0\right)-{\displaystyle \frac{F\left(0\right)}{2b}}$$

(46)

To make $C_1=0$ the initial velocity is

$${\dot{\zeta }}^{*}\left(0\right)={\displaystyle \frac{F\left(0\right)}{2b}}=0$$

(47)

Thus, for any initial condition of the form $[{\zeta }^{*}\left(0\right),0]$ where ${\zeta }^{*}\left(0\right)$ is any convenient initial position, a limit cycle is obtained having bounded periodic motion. However, it is always preferred to have the motion symmetric about $\zeta =0$. This requires that $C_0=0$ of (43).

From (43) with $C_1=0$ at $t_o=0$

$$C_0={\zeta }^{*}\left(0\right)-\int {\displaystyle \frac{F\left(t\right)}{2b}}{|}_{t=0}$$

(48)

Setting $C_0=0$, the initial condition of the position ${\zeta }^{*}\left(0\right)$ is

$${\zeta }^{*}\left(0\right)=\int {\displaystyle \frac{F\left(t\right)}{2b}}{|}_{t=0}=-{\displaystyle \frac{kA}{2b\omega }}$$

(49)

The optimal, centered, periodic solution of (45), with its required initial condition $x_o$, is

$$\begin{array}{cc}\hfill {\zeta }^{*}& =\int {\displaystyle \frac{F\left(t\right)}{2b}}=-{\displaystyle \frac{kA}{2bw}}cos\left(\omega t\right)\hfill \\ \hfill {\dot{\zeta }}^{*}& ={\displaystyle \frac{F\left(t\right)}{2b}}={\displaystyle \frac{kA}{2b}}sin\left(\omega t\right)\hfill \\ \hfill x_o& =[{\zeta }^{*}\left(0\right),{\dot{\zeta }}^{*}\left(0\right)]=\left[-{\displaystyle \frac{kA}{2b\omega }},0\right]\hfill \\ \hfill \end{array}$$

(50)

Figure 2 shows simulation results using the model parameters of Table 1 for four different initial conditions indicated by black dots. When the initial condition lacks any initial velocity, [-1.75, 0], [-0.1594, 0], and [1.25, 0], the buoy exhibits bounded, periodic motion indicated by the closed, red contours in the phase space where [-0.1594, 0] produces the desired centered motion. The presence of initial velocity causes a gradual displacement growth over time, as shown by the blue trajectory of Figure 2 for initial conditions of [-1.4, 1]. It’s important to note that the buoy oscillation frequency aligns with the input frequency for linear systems, resulting in the repetition of limit cycles every T seconds.

Figure 3 shows the steady-state simulation results obtained after applying the optimal control law with the centering initial conditions, [-0.1594, 0]. As expected, the periodic motion is bounded and centered about $\zeta =0$. Furthermore, the energy extracted has a net positive slope.

In summary, to extract maximum energy using a linear point absorber, it is essential to apply the optimal control law of Equation (40) while ensuring appropriate initial conditions to achieve a bounded periodic motion. Alternatively, a reference tracking controller can be used to follow the bounded analytical periodic optimal states. This can be easily achieved when $F\left(t\right)$ is analytically integrable.

3.4. Nonlinear Point Absorber WEC

Nonlinear WEC models introduce challenges, such as obtaining analytical solutions to compute initial conditions to achieve limit cycles. This arises when ${\partial }_t{F}_h(\zeta ,t)$, used in the optimal acceleration trajectory of (37), is not integrable.

However, not all nonlinear systems exhibit this challenge due to the form of $F_h$. Consider a nonlinear system such that the $F_h$ of (31) is

$$F_h(\zeta ,t)=g\left(\zeta \right)+f\left(t\right)$$

(51)

where $g\left(\zeta \right)$ is a nonlinear function and $f\left(t\right)$ is any integrable function of time. For example, the Duffing Oscillator has this form where the $a_i$ are constant coefficients.

$$\ddot{\zeta }+a_1\dot{\zeta }+a_2\zeta +a_3{\zeta }^3=a_{4}cos\omega t-u$$

(52)

Based on (52)

$$\begin{array}{cc}\hfill g\left(\zeta \right)& =-a_2\zeta -a_3{\zeta }^3\hfill \\ \hfill f\left(t\right)& =a_{4}cos\omega t\hfill \end{array}$$

(53)

the optimal control law is

$$u^{*}=g\left(\zeta \right)+f\left(t\right)-a_1\dot{\zeta }-{\displaystyle \frac{{\partial }_{t}f\left(t\right)}{2a_1}}$$

(54)

and the rest becomes identical to the linear system of (45).

Unfortunately, this is a special case, and most models of nonlinear WECs possess nonlinearities of the form $\sum a_{ij}{f}^i\left(t\right)g^j\left(\zeta \right)$ where i and j are integers.

Consider the model of a nonlinear point absorber WEC with an hourglass-shape buoy [37]

$$M\ddot{\zeta }+b\dot{\zeta }+k{\left[\zeta -\eta \left(t\right)\right]}^3=-F_c$$

(55)

with $F_h=-k{\left[\zeta -\eta \left(t\right)\right]}^3$ that is due solely to the buoyancy force. From (34) the optimal control law is (55)

$$F_c^{*}=k{\left[\eta -\zeta \right]}^3-b\dot{\zeta }-{\displaystyle \frac{3Mk}{2b}}\dot{\eta }{(\zeta -\eta )}^2$$

(56)

with an optimal acceleration profile, from (37), of

$${\ddot{\zeta }}^{*}={\displaystyle \frac{3k}{2b}}\dot{\eta }{(\zeta -\eta )}^2$$

(57)

Since (57) can not be integrated in closed form, obtaining expressions for the initial conditions that give limit cycle motion is impossible. Thus, we must resort to numerical techniques.

To numerically explore the initial condition space, the system of (55) and the control law of (56) were simulated using a third-order, explicit Runge-Kutta solver with a time step of 0.01 s with the parameters values of Table 1. The results for two different wave amplitude cases are shown in Figure 5, where the simulation duration was 50 s. The left figure is for $A=1.3m$, and on the right, $A=0.65m$.

The blue dashed ellipses indicate the initial conditions resulting in limit cycles. The ellipse corresponding to the higher wave amplitude case has major and minor axes of $0.3m/s$ and $0.1m$, respectively, in contrast to the smaller wave amplitude case, whose major and minor axes are $0.5m/s$ and $0.3m$. The initial condition region increases with decreasing wave amplitude in these two cases. Other initial conditions may exist outside the ellipses that lead to a limit cycle. However, none were found using an extensive grid search of the phase space.

The red trajectories of Figure 5 show buoy motions when the initial conditions are inside the ellipses, indicated by black dots. While these trajectories provide anecdotal information regarding position and velocity extremums they cannot be used to determine limit cycle existence.

Poincare maps were used to confirm limit cycle behavior with the simulation time increased to 500 seconds. Given an initial condition, it is instructive to view the phase planes of Figure 5 with time being along a third orthogonal axis. The red trajectories then move around the phase plane in the direction of increasing time. Poincare maps were made from sampling points on this three-dimensional trajectory at the excitation period of $T=1$ seconds and are shown in Figure 6 for the two initial condition cases of Figure 5. Limit cycle motion is confirmed since the Poincare Maps are closed contours.

Figure 7 shows the dynamics of the nonlinear point absorber for a wave of amplitude A upon the implementation of the optimal control law and choosing proper initial conditions, [0.05, 0.1], within the corresponding blue ellipse of Figure 5. As with the phase plane plots of Figure 7, periodicity is difficult to infer from this relatively short time duration. Once again, the energy time history has a net positive slope.

To further validate the optimality of the control law for both the linear and nonlinear cases, the value of the Hamiltonian is plotted in Figure 8 and is zero for all time. This is a consequence of the PMP and ensures the simulated co-states match the solutions of the co-state differential equations of Equation (35) arising from the Euler-Lagrange equations.

Overall, the presented figures and analyses demonstrate the promising performance of the optimal control approach, providing valuable insights into the dynamic behavior and energy production capabilities of the considered WEC systems.

4. Conclusion and Future Work

Optimal control analysis can be used to maximize the energy of a renewable source even when its second-order model is nonlinear and nonautonomous. To achieve bounded motion the initial condition region resulting in a stable limit cycle should be included in the optimal control analysis. The size and location of the region will, in general, be a function of the wave amplitude and frequency.

A relatively small set of initial conditions was found for the nonlinear hourglass-shaped buoy. It may be possible to find buoy shapes yielding larger initial condition regions extending the practical application of the optimal control law. Consequently, using the control law during shape optimization is one of the next steps to be exploited, finding the form of the nonlinearity that leads to a more extensive set of initial conditions. Finding a more systematic way to obtain initial conditions to drive the system to a limit cycle should also be explored. It should be noted that while the examples presented here used regular waves, the optimal control law form accommodates irregular waves as well.

Initial conditions and wave parameters can result in intriguing limit cycle behavior that could be exploited for operational reasons. One case is illustrated in Figure 9 for $A=0.611$ and optimal initial conditions of [0.2, 0]. The buoy response progresses between the five closed contours of the Poincare map. Shaping the location, number, and size of the contours could be used to manage PTO bounds and motion constraints.

5. References