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

Abstract

Pursuing sustainable energy solutions has prompted researchers to focus on optimizing energy extraction from renewable sources. Control laws that optimize energy extraction require accurate modeling, often resulting in time-varying, nonlinear differential equations. An energy-maximizing optimal control law is derived for time-varying, nonlinear, second-order, energy harvesting systems. We demonstrate that sustaining periodic motion under this control law when subjected to periodic disturbances necessitates identifying appropriate initial conditions, inducing the system to follow a limit cycle. The general optimal solution is applied to two point absorber wave energy converter models: a linear model where the analytical derivation of initial conditions suffices and a nonlinear model demanding a numerical approach. A stable limit cycle is obtained for the latter when the initial conditions lie within an ellipse centered at the origin of the phase plane. This work advances energy-maximizing optimal control solutions for nonautonomous nonlinear systems with application to point absorbers. The results also shed light on the significance of initial conditions in achieving physically realizable periodic motion for periodic energy harvesting systems.

1. Introduction

Optimal control theory provides a powerful framework for determining a dynamic system’s control law that extremizes an objective function [1]. Applied across various fields—engineering, economics, and biology, to name a few—this framework is essential for energy maximization in second-order, nonlinear, time-varying dynamic systems.

One challenging class of optimal control problems is Linear Optimal Problems (LOPs), also known as singular problems. These arise when control inputs appear linearly in the system’s state equation and cost function, even though state equations may be nonlinear. This structure complicates forming the control law directly using the stationary condition of optimality [2]. Several researchers have investigated singular optimal control problems, focusing on solution existence and formulation [3,4,5,6,7]. Scardina [8] explored the necessary conditions for optimality in singular problems using Pontryagin’s Minimum Principle (PMP), leveraging the constancy of the Hamiltonian to achieve optimality while satisfying system dynamics, co-state equations, and boundary conditions [6,7].

Kelly et al. and Robbins [9,10,11] introduced the generalized Legendre–Clebsch condition to address optimality in singular problems, particularly where the Hamiltonian is linear in the control inputs. Extensive study of this condition has provided valuable insights into the existence and properties of optimal solutions [11,12].

Finding optimal solutions for nonlinear, nonautonomous systems often requires both a control law and initial conditions that keep system states on an optimal path [7,13]. Limit cycles, representing bounded periodic motion, are crucial for achieving efficient, sustained energy absorption in applications like point absorber wave energy converters subject to periodic waves [14]. Given the complexity of nonlinear, nonautonomous systems, linearization may be acceptable in some cases [15], though this study leverages the nonlinear, time-varying nature of the model.

In recent years, optimal control has drawn significant interest for energy extraction, especially in the design and control of Wave Energy Converters (WECs). These devices convert ocean wave energy into electricity through a Power Take-Off system (PTO). The authors in [16] provide a review and classification of WEC technologies. This paper focuses on point absorber WECs, typically modeled linearly using Cummins’ equation [17], suitable for small motion around an equilibrium. However, energy capture can increase significantly with large motions that introduce nonlinearities [18], such as the Froude–Krylov force [19] from integrating incident wave pressure over the buoy’s wetted surface.

Diverse methods have been explored to optimize WEC energy output. The shape of a WEC significantly impacts its response and energy harvesting capacity, with various strategies proposed for shape optimization [20,21,22]. In parallel, optimizing WEC array geometry has been investigated to exploit WEC-to-WEC interactions while minimizing disturbances [23]. Control laws also critically influence WEC performance, and control system development has been an active research area. Proposed control strategies include sliding mode control [24,25], model predictive control [26,27,28], latching techniques [29,30], and various other maximizing schemes [31,32,33,34], each offering unique improvements in energy efficiency.

A well-known energy-optimal strategy for linear point absorbers is Complex Conjugate Control (CCC) [35]. Originating with Jacobi’s 19th-century maximum power transfer law, CCC has become a widely studied method for maximizing energy capture. Impedance matching, an extension of Jacobi’s law, is also applied in engineering fields, including vibration energy removal [36,37,38]. While CCC achieves optimal energy extraction, it can exaggerate device motion and is limited to linear WEC models, often conflicting with small motion assumptions [19]. This study addresses these limitations by developing an energy-maximizing control strategy for nonlinear, nonautonomous, second-order periodic systems, as derived in [39].

Establishing limit cycles for WECs ensures periodic and bounded motion, essential for efficient energy conversion. Previous studies [40] derived feedback linearizing optimal control laws for nonlinear WECs, producing closed-loop limit cycles, while others assumed linear models under regular waves to obtain explicit control expressions, equivalent to CCC in linear cases [41]. However, these studies often relied on switching functions to bound motion, without exploring limit cycles for optimal initial conditions.

This article builds upon existing literature by presenting a general optimal control law for energy maximization in second-order, nonlinear, time-varying systems. Our approach extends singular control analysis and investigates limit cycles to determine optimal initial conditions, applying both linear and nonlinear models inspired by recent WEC advancements. A significant finding is that our derived control law naturally reduces to the classical CCC approach in the linear case, demonstrating consistency with this trusted method. This serves as an indirect validation and suggests that its extension to nonlinear systems, as conducted here, should also be effective. By maintaining CCC’s core principles while adapting them to nonlinear dynamics, our method offers a robust, theoretically sound strategy for energy maximization in both linear and nonlinear applications.

The paper is organized as follows: Section 2 derives the general optimal control law, Section 3 applies it to linear and nonlinear WEC point absorber models, Section 4 discusses the numerical results from Section 3, and Section 5 provides conclusions and directions for future research.

2. General Optimal Control Law

The goal is to develop a control strategy designed to hypothetically maximize energy extraction for energy harvesting systems. Our approach applies optimal control theory, a classical theory, where the aim is to derive a control law—a guiding rule—that directs the system to follow an ideal trajectory over time. This trajectory, chosen to maximize the harvested energy, is derived based on the system’s dynamic behavior and the external forces it encounters. The control law defines how the system’s motion should be continuously adjusted to remain on this optimal path, ultimately maximizing energy output while ensuring stable and periodic operation. To demonstrate this method, we apply it to a wave energy converter models.

The details of deriving the optimal control law will be covered in the section below. After introducing the form of the second-order dynamic system, the optimal control law is derived, including examining the second-order optimality conditions. We begin by deriving the control law for a broad class of nonlinear, time-varying systems and then apply it to specific examples to showcase its capabilities.

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)

Equation (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 Equation (1), $\dot{\overrightarrow{x}}=\overrightarrow{f}(\overrightarrow{x},F_e(\overrightarrow{x}),u)$, is shown below

$$\begin{array}{cc}\hfill {\dot{x}}_1& =x_2\hfill \\ \hfill {\dot{x}}_2& ={\displaystyle \frac{F_e(\overrightarrow{x})-u}{M}}\hfill \\ \hfill {\dot{x}}_3& =1\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 Equation (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 [7].

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(t_o)$, 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 Equation (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 on the power electronics. Thus, our objective function is framed to maximize the mechanical energy extraction, or conversely, minimize its negative. It is important to note that, while limiting the Power Take-Off (PTO) force is a practical consideration in many applications, our optimization problem does not impose specific constraints on the PTO force or system motion during the derivation of the optimal control law. The primary objective here is to derive and analyze the control law that maximizes energy extraction without predefined limitations on force or motion, allowing for a pure examination of the optimal trajectory and its energy-maximizing potential. Consequently, the optimization problem is focused on achieving this ideal trajectory and evaluating its implications under unrestricted conditions. Our optimization problem is thus:

$$\begin{array}{cc}\hfill \underset{u\in [t_o,t_f]}{min}:& -E(x_2,u)\hfill \\ \hfill & \dot{\overrightarrow{x}}(t)=\overrightarrow{f}(\overrightarrow{x},F_e(\overrightarrow{x}),u)\hfill \\ \hfill & x(t_o)=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,7]

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

(5)

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

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

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

(6)

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

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

(7)

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

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

(8)

The stationary condition ${\partial }_{u}H(\overrightarrow{x},\overrightarrow{\lambda },F_e(\overrightarrow{x}),u)=0$ of Equation (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 provides additional boundary conditions at the terminal time, $t_f$ [7].

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

(10)

Since H is linear in u, Equation (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(\overrightarrow{x})$, as shown in the remainder of this section.

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

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

(11)

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

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

(12)

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

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

(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 [7].

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

(14)

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

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

(15)

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

$$x_2^2{\partial }_{x_2}{F}_e(\overrightarrow{x})+{\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 Equation (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(\overrightarrow{x})$$

(17)

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

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

(18)

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

$$u^{*}=F_e(\overrightarrow{x})+M\left({\displaystyle \frac{x_2^2{\partial }_{x_1}{\partial }_{x_2}{F}_e(\overrightarrow{x})+x_2{\partial }_{x_3}{\partial }_{x_2}{F}_e(\overrightarrow{x})+{\partial }_{x_3}{F}_e(\overrightarrow{x})}{2{\partial }_{x_2}{F}_e(\overrightarrow{x})+x_2{\partial }_{x_2}^2{F}_e(\overrightarrow{x})}}\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 Equation (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 Equation (21) in the form of Equation (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 Equation (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 Equation (20).

2.2. Second-Order Optimality Condition

In addition to satisfying the first-order Euler–Lagrange equations of Equation (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 Equation (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$, Equation (25) can be simplified to Equation (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 Equation (26) to the expression of Equation (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(\overrightarrow{x})+\dot{\zeta }{\partial }_{\dot{\zeta }}^2{F}_e(\overrightarrow{x})}{M^2}}$$

(27)

According to Equation (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 Equation (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 [42], the waves; thus, our boundary conditions will be written as follows:

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

(28)

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

For the linear case, Equation (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, as shown in Equation (29), for the heave-only case [17].

$$M\ddot{\zeta }+{\int }_{-\infty }^t{h}_r(t-\tau )\dot{\zeta }(\tau )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 Equation (30).

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

(30)

Next, we will generalize Equation (30) to include nonlinear hydrodynamic terms while also showing the PTO force, $F_c$, 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 [43]. Comparing Equation (31) to Equation (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)

When examining simulation results, we will use the physical parameters shown in

Table 1. WEC parameters.

Feature	Symbol	Value	Units
Mass	M	109,626	kg
Radiation Damping	b	20,000	N/(m/s)
Hydrostatic Stiffness	k	30,819	N/m
Wave Amplitude	A	$1.3$	m
Wave Frequency	$\omega$	$2\pi$	rad
Wave Period	T	1	s

from [39] for the linear and nonlinear examples below. However, though the optimal control law of Equation (20) applies to any continuous, nonlinear, time-varying, second-order model, to ease the comparison with previous literature [40] we will restrict our attention to regular waves with elevation

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

(33)

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

In the case of irregular waves, our control law still applies by assuming

$$\eta (t)=\sum _i^n{A}_{i}sin({\omega }_{i}t)$$

(34)

where n is the number of wave components and $A_i$ and ${\omega }_i$ are the corresponding wave amplitude and frequency of every component. Moreover, $b\dot{\zeta }$ can be replaced by a more general, perhaps nonlinear, damping expression that better fits irregular waves.

3.2. Optimal Control Law

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

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

(35)

where the co-states of Equations (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 \end{array}$$

(36)

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

$${\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$$

(37)

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 (35). Substituting Equation (35) into Equation (31) yields the optimal acceleration

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

(38)

Equation (38) 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 (36).

3.3. Linear Point Absorber WEC

Consider the model of a linear point absorber below for the buoy of Figure 1, which has been used previously to derive the Complex Conjugate Control (CCC) law [18].

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

(39)

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

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

(40)

which results in

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

(41)

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

However, one often-cited practical limitation of CCC is that the buoy motion becomes larger than the model assumptions used in its derivation [24,35,44,45]. 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 by obtaining the proper initial condition of Equation (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 Equation (38), the optimal acceleration is

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

(42)

which can be integrated twice to obtain its velocity and position

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

(43)

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 Equation (44), the same result obtained from CCC [18].

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

(44)

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

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

(45)

and its corresponding model

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

(46)

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

From Equation (44), at time $t_o=0$,

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

(47)

To make $C_1=0$, the initial velocity is

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

(48)

Thus, for any initial condition of the form $[{\zeta }^{*}(0),0]$, where ${\zeta }^{*}(0)$ 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 Equation (44).

From Equation (44), with $C_1=0$ at $t_o=0$,

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

(49)

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

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

(50)

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

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

(51)

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 ∂_tF_h(ζ,t), used in the optimal acceleration trajectory of Equation (38), 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 Equation (31) is

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

(52)

where $g(\zeta )$ is a nonlinear function and $f(t)$ 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$$

(53)

Based on Equation (53)

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

(54)

the optimal control law is

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

(55)

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

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

Consider the model of a nonlinear point absorber WEC with an hourglass-shaped buoy [40],

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

(56)

with $F_h=-k{\left[\zeta -\eta (t)\right]}^3$, which is due solely to the buoyancy force. See also Figure 2. From Equation (35), the optimal control law is Equation (56),

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

(57)

with an optimal acceleration profile, from Equation (38), of

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

(58)

Since Equation (58) cannot 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, as illustrated in Section 4.2.

4. Numerical Simulations and Discussion

All simulations were simulated using 2023b MATLAB’s third-order, explicit Runge–Kutta solver with a time step of 0.01 s and with the parameter values of Table 1.

4.1. Linear WEC

Figure 3 shows the 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 3 for initial conditions of [−1.4, 1]. It is 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 4 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$ due to the choices of the proper initial conditions that lead the system to exhibit a limit cycle. A decent amount of energy is harnessed while undergoing a realistic motion.

In summary, to extract maximum energy using a linear point absorber, it is essential to apply the optimal control law of Equation (41) 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(t)$ is analytically integrable.

4.2. Nonlinear WEC

As previously discussed, identifying initial conditions that guide the system toward a stable limit cycle analytically is infeasible; therefore, a numerical approach was used to explore the initial condition space. Figure 5 illustrates the resulting phase portraits for two different wave amplitudes, $A=1.3\mathrm{m}$ (left) and $A=0.65\mathrm{m}$ (right), over a 50-s simulation.

The blue dashed ellipses represent the regions of initial conditions that produce limit cycles. For the larger wave amplitude ($A=1.3\mathrm{m}$), the ellipse has smaller major and minor axes (approximately $0.3\mathrm{m}/\mathrm{s}$ and $0.1\mathrm{m}$). This contrasts with the smaller wave amplitude case ($A=0.65\mathrm{m}$), where the region of initial conditions expands, with major and minor axes of approximately $0.5\mathrm{m}/\mathrm{s}$ and $0.3\mathrm{m}$. This trend suggests that as wave amplitude decreases, the initial condition space required to drive the system into a stable limit cycle becomes larger, indicating a higher degree of flexibility in selecting initial conditions for smaller amplitudes.

The red trajectories in each plot represent the system’s motion over time when the initial conditions fall within the corresponding blue ellipses (indicated by black dots). These trajectories display complex patterns of velocity and position oscillations that depend heavily on the initial energy input from the chosen starting points. In both amplitude cases, the trajectories converge within bounded regions, which implies that the system can achieve stable, periodic motion when initial conditions are carefully chosen within the identified ellipses.

This figure also highlights the effect of wave amplitude on the nature of limit cycles. For higher amplitudes, the system exhibits more pronounced oscillations in both velocity and position, as seen in the denser and broader phase-space coverage in the left plot. This observation is consistent with the increased energy associated with larger wave amplitudes, resulting in greater system displacement and velocity variations. Conversely, the smaller amplitude case on the right shows more compact oscillations in phase space, reflecting the reduced energy input.

Lastly, while these initial condition ellipses capture the main regions that lead to limit cycles, the figure does not exclude the possibility of additional regions outside these ellipses that may also yield stable cycles. However, an extensive grid search of the phase space did not reveal any such additional regions. This focused analysis reinforces that, within the explored parameter space, the blue ellipses capture the essential regions for achieving bounded periodic motion.

Poincare maps were employed to confirm the presence of limit cycle behavior, with the simulation time extended to 500 s. By sampling the system’s position and velocity at intervals corresponding to the excitation period, $T=1$ second, we obtain discrete points in the phase space, capturing the system’s state at each cycle. Figure 6 presents these Poincare maps for the two initial condition cases from Figure 5. The closed contours observed in the maps confirm the existence of stable limit cycles, as the system consistently revisits these states, indicating periodic behavior. This approach reinforces that, once the system is driven to a limit cycle, it remains on this bounded, repeating trajectory indefinitely, provided that external conditions and parameters remain unchanged.

In Figure 7, we examine the dynamic behavior of the nonlinear point absorber Wave Energy Converter (WEC) over a shorter period of 50 s for a wave amplitude of A, using initial conditions within the previously identified stable region ([0.05, 0.1]).

The energy time history shows a net positive trend, demonstrating the energy harvesting potential of the system under the derived optimal control law. Despite fluctuations, the overall energy increases steadily, indicating that the control strategy effectively extracts energy from the incoming waves over time.

The applied force (PTO force) oscillates symmetrically around zero, with peak values that remain within reasonable bounds for the chosen initial conditions and wave amplitude. This symmetry suggests that the control law balances the extraction and dissipation of energy, maintaining the system within its operational limits while achieving efficient energy capture.

Both velocity and position exhibit clear periodic oscillations, aligning with the Poincare maps that confirmed the limit cycle behavior. The velocity oscillations are bounded, with maximum and minimum values stabilizing over time, demonstrating the control law’s capacity to achieve sustainable, repeatable motion. Similarly, the position oscillations show that the WEC remains within a predictable range, ensuring that it operates within the physical constraints of the system.

These insights highlight the successful application of the optimal control law, as it achieves stable, periodic motion while having positive energy extraction. Together, Figure 6 and Figure 7 provide a comprehensive view of the system’s behavior under the optimal control law, emphasizing the effectiveness of the derived initial conditions and control strategy.

To further validate the optimality of the control law for both linear and nonlinear cases, the value of the Hamiltonian is plotted in Figure 8. This figure displays the Hamiltonian as a function of position, velocity, and time, highlighting its conservation over time when using the derived optimal control law, which confirms that the control law derived not only ensures positive energy harvesting but also maintains Hamiltonian conservation across both linear and nonlinear systems. Thus, the simulated co-states align with the solutions of the co-state differential equations of Equation (36), validating the efficacy of the control strategy.

Initial conditions and wave parameters can lead to diverse and intriguing limit cycle behaviors that may be advantageous for specific operational purposes. One example is illustrated in Figure 9, where the wave amplitude $A=0.611$ and the initial conditions are set to [0.2, 0]. The phase space trajectory (left) reveals a complex pattern where the buoy response oscillates between distinct regions, indicating a multi-contour limit cycle.

The corresponding Poincare map (right) highlights this behavior further, showing five closed contours. These contours represent the system’s state at consistent intervals. The presence of multiple closed contours implies that the buoy’s motion does not settle into a single, fixed limit cycle but instead alternates among several periodic states. This type of behavior could be particularly useful in managing the Power Take-Off (PTO) force and other operational constraints by allowing the system to distribute its motion across multiple states, potentially reducing peak loads.

By strategically adjusting the location, number, and size of these contours through the selection of wave parameters and initial conditions, operators could exert finer control over the system’s motion constraints. This flexibility can be used to achieve desired energy harvesting performance while ensuring that the buoy’s motion remains within predefined bounds, thus mitigating potential wear and tear on the PTO mechanism.

In summary, the multi-contour limit cycle observed here showcases an operationally beneficial phenomenon, where adjusting system parameters allows for controlled oscillations across multiple stable regions, which enables adaptable operational behaviors.

While this control law was derived to maximize energy extraction, we do not claim it to be the ultimate or most optimal method for all conditions. Rather, it represents a promising approach with considerable potential, particularly given its flexibility and effectiveness across both linear and nonlinear systems. Notably, in the linear case, this control law reduces to the well-established Complex Conjugate Control (CCC) strategy—a key finding that validates the robustness of our approach and links it to proven energy-maximizing methods. Extending this control law to nonlinear systems has further demonstrated its utility, showing that it can successfully maintain stability and maximize energy capture even in more complex, nonlinear dynamics. This generality and adaptability make it a valuable tool, potentially applicable to a wide range of energy-harvesting systems beyond the specific cases tested here, and a strong foundation for further optimization or refinement.

A study by Wilson et al. [40] sought to extend the Complex Conjugate Control (CCC) method to nonlinear systems, specifically applying it to the same nonlinear system analyzed in our work. Using identical wave conditions and buoy parameters, we observed that our derived control law outperformed their optimal resistive damping feedback controller, generating significantly more energy over the same period. This comparison underscores the effectiveness of our approach in nonlinear scenarios, where traditional extensions of CCC may be limited. By moving beyond the resistive feedback framework, our control law demonstrates a more comprehensive energy-harvesting capability, making it a compelling alternative in complex, nonlinear energy-harvesting applications. This highlights both the robustness and practical advantage of our method, especially when applied to systems that involve nonlinear dynamics.

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.

5. Conclusions and Future Work

This study has successfully developed and validated an optimal control approach for maximizing energy extraction in nonlinear, nonautonomous, second-order systems, exemplified through two types of WEC model. By deriving an optimal control law tailored for such complex systems, we have provided a robust framework capable of handling both linear and nonlinear dynamics without the need for switching functions or approximations. This control law has proven effective in sustaining stable, periodic motion, which is essential for efficient and consistent energy conversion.

The findings underscore the significance of carefully selected initial conditions that lead to bounded, periodic trajectories, as these result in limit cycles that allow the WEC to operate within its physical constraints while optimizing energy capture. Importantly, the initial condition space was observed to vary with wave amplitude, indicating the need for adaptive strategies in dynamic, real-world ocean environments.

Applying this control law to both linear and nonlinear WEC models, we observed that it reduces to the established complex conjugate control method for linear cases, thereby confirming the law’s consistency with well-recognized energy-maximizing strategies. Extending this control to nonlinear systems, we noted superior energy harvesting capabilities when compared with traditional resistive feedback methods, particularly under high wave amplitudes. This reveals the control law’s flexibility and its potential to enhance energy capture across a range of nonlinear conditions.

Future research will expand on this foundation by exploring two main areas. First, integrating this control law within buoy shape optimization could yield configurations that maximize initial condition flexibility and stability. Second, a systematic approach to selecting initial conditions for achieving stable limit cycles will be developed to ensure practical implementation in variable and unpredictable marine environments. Real-time experiments could be carried out with frictionless testbeds [46] to observe the performance of the controller.

Furthermore, while this work focused on regular wave patterns, the control law’s inherent structure accommodates irregular waves, which offers exciting possibilities for deploying it in real-world ocean energy applications where wave patterns are often unpredictable. This adaptability enhances the practical value of the proposed control strategy, which holds promise not only for wave energy conversion but also potentially for other renewable energy systems with these irregularities.