Force-Limited Control of Wave Energy Converters using a Describing Function Linearization
Rebecca McCabe Maha Haji



 


[footnoteinfo]This material is based on work supported by National Science Foundation Graduate Research Fellowship Grant No. DGE–2139899.

First]Rebecca McCabe 
Second]Maha N. Haji 
*Sibley School of Mechanical and Aerospace Engineering, 
Cornell University, 
   Ithaca, NY 14853 USA 
   [First]e-mail: rgm222@cornell.edu
   [Second]e-mail: maha@cornell.edu

§ ABSTRACT                 
Actuator saturation is a common nonlinearity. In wave energy conversion, force saturation conveniently limits drivetrain size and cost with minimal impact on energy generation. However, such nonlinear dynamics typically demand numerical simulation, which increases computational cost and diminishes intuition. This paper instead uses describing functions to approximate a force saturation nonlinearity as a linear impedance mismatch. In the frequency domain, the impact of controller impedance mismatch (such as force limit, finite bandwidth, or parameter error) on electrical power production is shown analytically and graphically for a generic nondimensionalized single degree of freedom wave energy converter in regular waves. Results are visualized with Smith charts. Notably, systems with a specific ratio of reactive to real mechanical impedance are least sensitive to force limits, a criteria which conflicts with resonance and bandwidth considerations. The describing function method shows promise to enable future studies such as large-scale design optimization and co-design.



Wave energy converters, constrained control, systems with saturation, nonlinear and optimal marine system control, describing functions, impedance mismatch, linearization.








§ INTRODUCTION





 §.§ Motivation

Ocean wave energy converters (WECs) are an immature yet promising source of renewable energy to decarbonize coastal grids and offshore systems. Maximum power transfer requires impedance matching of WEC controls, powertrain, and hydrodynamics, but plant uncertainty, controller bandwidth, and physical constraints prevent perfect matching. Actuator force limits are especially relevant given waves' high-force, low-speed nature and the scaling of device cost with force. WEC controllers must maximize power while obeying force limits, although the resulting nonlinearity requires computationally-costly numerical optimization. There is a gap of rapid and intuitive methods suitable for early tradeoff analysis. To this end, the present paper demonstrates linear analytic treatment of force saturation limits and other sources of impedance mismatch.



 §.§ Literature Review

Prior studies on WEC powertrain constraints investigate position, velocity, force, rate of change of force, and power flow direction. Constrained numerical optimization is typically used. For example, <cit.> review model predictive control while <cit.> apply pseudo-spectral methods.


A minority of work tackles the problem analytically. <cit.> use the Pontryagin principle to show that a bang-singular arc-bang controller is optimal, revealing that saturating the unsaturated optimal solution can still be optimal in certain cases. 
They derive analytical piecewise expressions for the optimal control force, but finding the corresponding power still requires numerical simulation.
<cit.> present a geometric tool to analyze simultaneous force and position constraints. They provide analytical relationships between power and root-mean-square signals, which relate to upper and lower bounds on constraint values, but the bounds are not tight. <cit.> introduce another geometric tool which accounts only for position constraints and focuses on hydrodynamics over powertrain. The tool is visually and mathematically similar to a Smith chart. To the authors' knowledge, no prior work uses describing functions to address WEC constraints, though <cit.> use them to model a WEC with nonlinear damping, and <cit.> use the same technique to model drag, calling it a Fourier approximation. Outside of wave energy, <cit.> apply describing functions to model velocity and torque saturation on an impedance-controlled robot and experimentally validate the results.












 §.§ Paper Outline and Contribution

Section <ref> presents the paper's first contribution: application of linear theory to analyze the relationship between impedance mismatch, signal amplitude, and power for a WEC in regular waves. Explicit analytical expressions are derived, and results are visualized on a Smith chart. As the second contribution, section <ref> suggests that constraints typically solved numerically, such as force limits, be linearized with describing functions in order to apply the previous section's results. The approximation is derived and discussed in the context of impedance mismatch. 

These contributions provide three main benefits. First, they offer an intuitive process for rapid tradeoff analysis early in the design process, with the opportunity to apply standard linear frequency-domain tools. Second, they are computationally efficient enough to integrate with more expensive techniques like design optimization or control co-design, either directly or as an initial guess for solvers of higher fidelity dynamics. Third, the results are analytical and differentiable, providing gradients for sensitivities or to accelerate convergence of outer optimizations.



§ PEAK LIMITING IN THE LINEAR CASE




 §.§ General Impedance-Mismatched System

The analysis starts with a generic linear system modeled as a Thévenin equivalent circuit with AC voltage source V_th and complex source impedance Z_th, shown in Fig. <ref>.


The load impedance Z_L is to be selected, with the conflicting goals of maximizing average power transfer P_L and minimizing the peak amplitude of load current |I_L| or voltage |V_L|. Maximum power transfer occurs when there is impedance matching, meaning Z_L = Z_th^* where * indicates complex conjugate. The load average power, peak voltage, and peak current at this matched point, denoted P_L^m, |V_L^m|, and |I_L^m| respectively, are found as:

    P_L^m = |V_th|^2/8 (Z_th), 
        |V_L^m| = |V_th| |Z_th|/2 (Z_th), 
        |I_L^m| = |V_th|/2 (Z_th)

where  means the real part and  the imaginary part.

To consider all possibilities of the unmatched case, we set Z_L = z Z_th^* for arbitrary complex number z. The space of z can be visualized using a Smith chart, where (z) is on a curved horizontal axis and (z) is on a curved vertical axis. The axes are curved such that the chart can be simultaneously read as a standard polar plot of the complex reflection coefficient Γ, which is a transformation of z defined as Γ = z-1/z+1. The impedance-matched case of z=1, Γ = 0 is found at the center of the plot, the minimum voltage at z=0, Γ = -1 on the left, and the minimum current at z →∞, Γ = 1 on the right.

The average power, peak voltage, and peak current in the unmatched case can be found using standard circuit techniques and expressed as fractions of their matched counterparts. <cit.> derives the power ratio:

    P_L/P_L^m = 1 - |Γ|^2

This relationship is visualized on the Smith chart in Fig. <ref>. As the impedance ratio z gets further away from the impedance-matched condition z=1 at the center of the circle, the power lowers quadratically.





The corresponding voltage and current ratios are:

    |V_L|/|V_L^m|,|I_L|/|I_L^m| = √(|Γ|^2 +  2 ϵ(Γ) + 1/α^2 |Γ|^2 + 2 α(Γ) + 1)

where ϵ indicates sign: ϵ=1 for voltage ratio and ϵ=-1 for current ratio, and α = (Z_th)/(Z_th) is a parameter related to the phase of the source impedance. These relationships are visualized on the Smith charts in Fig. <ref> (a) and (b), where contours show the ratios as a function of z (and thus Γ) for various values of α. Only positive α (inductive Z_th) is shown for brevity. The contours for negative α (capacitive Z_th) can be found by reflecting the graphs over the horizontal axis ((Γ) = 0) due to symmetry. On the Smith charts, points where the voltage and current ratios exceed 1 are shaded. These points are undesirable because they absorb less power than the baseline z=1, Γ=0 matched case while having higher peaks. The optimal contours (lowest voltage and current ratios for a given power ratio) are traced out with dashed lines. Because the optimal voltage reduction path requires decreasing the impedance and the optimal current reduction path requires increasing the impedance, it is not possible to follow both paths simultaneously. For sufficiently high values of |α|, it is possible to reduce both current and voltage simultaneously (i.e. avoid the shaded region of both Smith charts in Fig. <ref> (a) and (b)), although for α=0, decreasing voltage implies increasing current and vice versa. This tradeoff is explored further in section (c) of Fig. <ref>, a pareto front showing the nondominated combinations of the voltage, current, and power ratios.

The optimal paths (·)^opt of Fig. <ref> are derived by setting the derivative of the voltage and current ratios (<ref>) with respect to ∠ Γ equal to zero, for ∠ Γ= tan^-1 ((Γ)/(Γ)):

    ∠ Γ ^opt = 2 tan^-1[ α^2 |Γ|^2 + 1/σ + ϵα (1 + |Γ|)^2] + ϵcos^-1[ -2α| Γ|/σ]

where σ = √((α^2|Γ|^2 + 1)^2 + α^2 (|Γ|^2 + 1)^2) and sign indicator ϵ=1 obtains the angle at minimum voltage and ϵ=-1 minimum current. <cit.> obtained a similar result for a two-port system. Further manipulation of (<ref>) shows that the minimum-voltage and minimum-current angles are supplementary: ∠ Γ^opt,V+ ∠ Γ^opt,I= π. Using this relation and plugging into (<ref>) reveals that the ratios are equal at their respective optima: |V_L|/|V_L^m|^opt,V= |I_L|/|I_L^m|^opt,I. Intuitively, voltage and current reductions are symmetric.

The optimal contours are aggregated in Fig. <ref> to show the tradeoff between power and voltage/current. Interestingly, as |α| grows (i.e., as Z_th becomes less resistive and more reactive), there is less of a power penalty for a given voltage or current reduction. In other words, power in the pure reactive Z_th case is least sensitive to voltage and current limits. For fixed (Z_th),  (<ref>) shows the baseline matched power is independent of α, suggesting that plant design should maximize |α|.


Note that most of the optimal contours of Fig. <ref> (a) and (b) require a z with nonzero imaginary part, i.e. a Z_L that is more or less reactive than the impedance-matched case, rather than merely scaled up or down. This highlights the difference between “constrained optimal control" and “optimal constrained control." The former refers to the unconstrained optimal controller that has been scaled or saturated until it meets the constraint, while the latter refers to the controller that is optimal for the constrained problem, which is distinct under the present assumption of linear control. Scaling down or saturating the control signal computed with the unconstrained optimal impedance yields a signal with a fundamental amplitude identical to one computed with a proportionally scaled linear control impedance. This enforces (z)=0, which is evidently non-optimal for all but α=0. This is distinct from the complex z in the optimal profiles of (<ref>) and Fig. <ref>, which would be considered optimal constrained control. Interestingly, <cit.> show that the optimal constrained and constrained optimal controllers are identical if nonlinear control is allowed. In summary, classical theory reveals the effect of impedance mismatch on power, current, and voltage, informing the choice of a load impedance to balance power generation and peak limiting constraints.



 §.§ Wave Energy Converters



Applying the preceding analysis to a WEC requires a Thévenin equivalent circuit for WEC dynamics. This study assumes a single degree of freedom floating body coupled to a power take-off (PTO) with a drivetrain and a linear or rotational synchronous surface permanent magnet electric generator. Substituting a hydraulic or other impedance is straightforward provided the system remains linear. 

The generator model is non-ideal and the objective is electrical power. The frequency domain WEC dynamics are:

    ((m+A)s^2 + B_h s + K_h) X + F_P = F_e       Body
       G τ_PTO = F_PTO, G s X = Ω      Gear ratio
       τ_PTO = (B_d + K_d/s) Ω + τ_gen      PTO
       τ_gen = K_t I, V = I(R + s L) - K_t Ω      Generator
       V = (B_c + K_c/s) I       Controller
       P_elec = 0.5 (I^* V)       Power

with Laplace variable s, mass m, added mass A, hydrodynamic damping B_h, hydrostatic stiffness K_h, WEC position X, power take-off and wave excitation forces F_P and F_e, effective gear ratio G, drivetrain mechanical stiffness and damping K_d and B_d, generator torque and rotation speed τ_gen and Ω, generator torque constant K_t, generator q-axis current and voltage I and V, generator resistance and inductance R and L, controller stiffness and damping K_c and B_c, and average electrical power P_elec. The controller acts between V and I rather than τ_gen and Ω as is more typical, providing equivalent dynamics and making the relevant Thévenin equivalent more convenient to define. B_d and K_d can capture mooring and drag forces. Frequency dependence of hydrodynamic coefficients A, B_h, and F_e is omitted since this study assumes regular waves. Fig. <ref> shows a block diagram of the dynamics.


Since electrical, not mechanical, power is the focus, choose the controller as the load impedance: Z_L = Z_c = B_c + K_c/s. Thus the Thévenin equivalent load voltage/current are the generator q-axis voltage/current: V_L = V, I_L = I. Thévenin parameters from <cit.> are then:

    Z_th = Z_w + K_t^2 G^2/Z_m,     |V_th| = K_t G F_e/|Z_m|

with winding impedance Z_w = R + s L and mechanical impedance Z_m = B_h + G^2B_d + (m+A)s + K_h + G^2 K_d/s.
Substituting (<ref>) into (<ref>) and applying the Haskind relation between B_h and F_e yields the matched power:

    P_L^m = G_0 J/k𝒟/1+ℛ/𝒟 (1+α_m^2),  G_0 = 
             1    heave
    
             2    surge/pitch

with incident energy density J, wavenumber k, and gain G_0. <cit.> show this applies for any WEC shape. Nondimensional quantities ℛ=R B_h/(K_t^2 G^2), 𝒟=B_h/(B_h+G^2B_d),  and α_m = (Z_m)/(Z_m) are the normalized resistance and ratios of hydrodynamic to total damping and reactive to real impedance, respectively.

This matched power is not possible when there is controller impedance mismatch, either intentional with the intent to obey a given amplitude limit, or unintentional due to parameter uncertainty or controller bandwidth limitations in broadband waves. The mismatched power P_L is found by multiplying the matched power P_L^m  in (<ref>) by the power ratio in (<ref>). The required reflection coefficient Γ to meet constraints is found via (<ref>) or graphically via Fig. <ref> and <ref>. Doing so requires an expression for α, found via (<ref>):

    α = (Z_th)/(Z_th) =  ℒℛ (1+α_m^2) - 𝒟α_m/ℛ (1+α_m^2)+𝒟

with ℒ=ω L/ R. Typically ℒ≈ 0 since the wave period 2 π / ω far exceeds the winding electrical time constant L/R.





Section <ref> allows choice of Γ to limit q-axis variables V and I, but limiting other quantities may be desired. Limiting V has little direct utility. Instead, a limit on phase voltage V_s represents the voltage of a vector drive which affects the generator's torque-speed curve; a limit on position X represents actuator stroke or kinematic constraints; and a limit on apparent power S sizes the PTO including reactive power. These amplitudes are:

    X = F_e/s/Z_m + K_t^2 G^2/Z_w - z Z_th^*,    V_s^2 = V^2 + (L p Ω I)^2 
    S_max,min/P_L^m = P_L/P_L^m±|V_L|/|V_L^m||I_L|/|I_L^m|√(1+α^2)

where p is the number of machine poles. Smith plots similar to Fig. <ref> could be made to visualize the effect of these limits, but more parameters besides α must be swept to visualize the design space. Meanwhile, a generator force limit of |F_gen| ≤ F_max is achieved with a q-axis current limit of |I| ≤ I_max = F_max/K_t G, which is well-captured by the earlier parameterization. For brevity, the rest of this paper focuses on current (force) limits, which drive PTO cost and are of primary design interest. Thus, Γ is selected to set current ratio |I_L|/|I_L^m| = I_max/|I_L^m|, thereby enforcing the limit.

Besides facilitating control, the formulation also informs plant design. By inspection of (<ref>), the maximum-power plant design in the matched (unconstrained) case is (ℛ,𝒟,α_m)^opt=(0,1,0), intuitively minimizing loss due to resistance, friction, and reactive power. α_m relates to the uncontrolled mechanical damping ratio ζ and natural frequency ω_n as α_m = (ω^2-ω_n^2)/(2 ζωω_n), so the α_m=0 plant resonates (ω=ω_n) passively at the wave frequency.

Amplitude limits make design more complicated. Substituting (ℛ,𝒟,α_m)^opt into (<ref>) gives α=0, but section <ref> established that high |α| is desirable to minimize the effect of amplitude limits on power. Specifically, if ℒ=0, maximum |α| requires α_m^2=1+𝒟/ℛ instead of α_m=0. Meanwhile, the well-known Bode-Fano limit implies that for good broadband matching, |α| must instead be minimized, implying either α_m=0 or |α_m|→∞. Therefore, the three goals of maximizing matched power (min|α_m|), minimizing the effect of constraints (max|α|), and maximizing bandwidth (min|α|) conflict. For the best tradeoff, the plant α_m must maximize P_L across the wave spectrum, the subject of future control co-design work.




§ SATURATION NONLINEARITIES



Section <ref> considered a linear impedance mismatch Z_L = z Z_th^*. In linear control, the current waveform is sinusoidal, so the fundamental amplitude equals the peak amplitude and never exceeds current limit I_max. 
Nonlinear control allows non-sinusoidal current, increasing the fundamental. 




 §.§ Describing Functions


Consider a nonlinear saturation control law of the form

    I_L,sat(t) = I_max sat(I_temp(t)/I_max)

where sat is the unit saturation function and I_temp(t) is the time domain output of a linear impedance controller I_temp = V_L,sat/Z_C, stored temporarily in controller memory and not physically realized. 
Fig. <ref> illustrates this nonlinearity. Unlike the unsaturated case (a), the orange saturation blocks in cases (b)-(d) create harmonics. If I_temp(t) is sinusoidal, I_L,sat(t) is a saturated sine, shown in (e). Technically, I_temp is nonsinusoidal because I_L,sat harmonics propagate through linear blocks, shown with multiple orange lines in (b). However, the second-order low-pass plant dynamics Z_th substantially reduce harmonics in downstream quantities V_S,sat, V_L,sat, and I_temp. Thus, (c)-(d) neglect harmonics of I_temp. Cases (c) and (d) differ in the location of the approximation: (c) depicts the sinusoidal-input describing function method from <cit.>, taking the fundamental of the saturated signal, and (d) is the higher order sinusoidal-input describing function from <cit.>, preserving harmonics of V_L,sat. The higher order method is useful if the low-pass assumption is poor, which may be true for broadband WECs. The saturated-sine current signal is decomposed into a sum of harmonics:

    I_L,sat(t) ≈∑_n |I_L,sat,n| sin(n ω t + ψ)

where ψ is the same phase as I_temp(t), since saturation does not alter phase. The nth harmonic amplitude |I_L,sat,n| is defined as f_sat,n |I_temp|. The saturation factor f_sat,n is found via Fourier analysis:

    f_sat,n =
    
    1    ℐ≥ 1, n = 1
    
    0    ℐ≥ 1,   n ≠ 1 
    2/π(ℐ√(1 - ℐ^2) + θ)    ℐ <1,   n = 1 
     
    0    ℐ< 1,  n = 2,4... 
    4/πn √(1 - ℐ^2)sinθ - ℐcosθ/n (n^2 - 1)   ℐ< 1 ,  n = 3,5...

where ℐ=I_max/|I_temp| and θ = n sin^-1ℐ. Fig. <ref> depicts the first 7 harmonics. As the signal saturates, the fundamental decreases and higher harmonics emerge.


At each harmonic frequency n, the saturation block is approximated by gain f_sat,n, so the equivalent load impedance is Z_L = Z_C / f_sat,n rather than the unsaturated Z_L = Z_C. This makes the nondimensional impedance for the nth harmonic z_n = Z_C  / (f_sat,n Z_th^*). 
Now the response is the sum of linear responses at different frequencies:

    P_L = ∑_nP_L,n(z_n)

In classic describing functions, only the fundamental current is used: I_L,sat≈ I_L,sat,1. This system contains only one harmonic, allowing the use of the tools from section <ref>, though it is still only quasi-linear because f_sat,1 depends on amplitude. The fundamental of the saturated current is a factor f_sat,1/ℐ above the peak, thus nonlinear control allows higher current ratios for the same current limit. The resulting power increase is found from Fig. <ref> by moving on the x-axis a factor of f_sat,1/ℐ higher than a given starting point. If |I_temp|≫ I_max, then ℐ→0 and f_sat,1/ℐ→4/π, yielding a square wave with fundamental 4/π higher than is possible in linear control.



To compute ℐ otherwise, the dynamics are used to find complex-valued I_temp:

    I_L,sat,1/I_temp = V_th-V_L,sat/I_tempZ_th = V_th-Z_C I_temp/I_tempZ_th

Controller impedance Z_C is still undecided. A procedure aiming to maximize power might use equations (<ref>) and (<ref>) to select Z_C so fundamental |I_L,sat,1| approaches its unconstrained value |I_L^m|. However, invoking Pontryagin's principle avoids this effort. <cit.> show analytically that steady-state optimal nonlinear control of a force-limited WEC is simply unconstrained optimal control with saturation, so Z_C=Z_th^*. (<ref>) is then solved using amplitude condition f_sat,1 = |I_L,sat,1/I_temp| and phase condition (I_L,sat,1/I_temp) = 0. Combining with (<ref>) and (<ref>) gives a transcendental equation. If an analytical solution is desired, arcsine can be algebraically approximated.










 §.§ Limitations and Extensions

The sinusoidal-input describing function assumes a saturated-sine current waveform. This is reasonable for WECs (a) in regular waves, (b) with low-pass dynamics |Z_th(nω)| ≪ |Z_th(ω)| for n≥3. Critically, (a) fails in a broadband ocean environment. Accuracy in irregular waves remains to be tested. <cit.> suggest that sinusoidal-input describing functions still apply for some non-sinusoidal excitation, but <cit.> show that systems with multiple sinusoids as input require other methods. Even for regular waves, the filter assumption (b) could be violated for broadband WECs, such as small WECs, or for narrow-band WECs analyzed far below the resonant frequency. This should be assessed with representative frequency-dependent hydrodynamics. Meanwhile, adding multiple constraints is unproblematic since arbitrarily many amplitudes may be analyzed/limited in the linear system, as long as all nonlinearity is filtered by subsequent low-pass dynamics. However, bump-stop position saturation is one problematic example where the nonlinearity occurs after the dynamics. Future work includes an error analysis compared to fully nonlinear simulations.






























§ CONCLUSION

This work presented an analytical method to handle WEC generator force (current) constraints. First, linear theory was reviewed, and key relationships for an impedance-mismatched Thévenin-equivalent circuit were visualized on Smith charts and Pareto fronts. Then, specific wave energy dynamics were introduced, highlighting the implications of nondimensional resistance, damping, and mechanical impedance on power and sensitivity to amplitude limits. Next, describing function theory for the saturation nonlinearity was introduced, using the plant's low-pass nature to justify neglecting higher harmonics.

A variety of future work is possible. The effect of irregular waves and frequency-dependent dynamics on describing function accuracy should be investigated. The linear framework for voltage and current limits could easily extend to other constraints like position, velocity, phase voltage, or apparent power, and multiple degrees of freedom could be considered. When combined with cost estimates, this offers a powerful strategy for techno-economic tradeoff analysis and design optimization. All in all, while analytical approximations do not take the place of nonlinear numerical optimization, they provide the intuition and computational speed that may unlock improved WEC designs with substantial climate impact.

The code for this work is available open-source at <https://github.com/symbiotic-engineering/IFAC_CAMS_2024/>.


















