Nonlinear Optimal Control for a PMLSG-VSC Wave Energy Conversion Unit

Gerasimos Rigatos ^1,*, Pierluigi Siano ², Mohammed AL-Numay ³, Masoud Abbaszadeh ⁴, Gennaro Cuccurullo ⁵

1. Unit of Industrial Automation, Industrial Systems Institute, 26504, Rion Patras Greece

2. Department of Management & Innovation Systems, University of Salerno, Fisciano (SA), 84084, Italy

3. Department of Electrical Engineering, King Saud University, Riyadh, 11421, Saudi Arabia

4. Department ECS Engineering, Rensselaer Polytechnic Institute, 12065, Troy NY, USA

5. Department of Industrial Engineering, University of Salerno, Fisciano, 84084, Italy

* Correspondence: Gerasimos Rigatos

Academic Editor: Nagi Abdussamie

Special Issue: Marine Renewable Energy

Received: September 16, 2023 | Accepted: February 04, 2024 | Published: February 21, 2024

Journal of Energy and Power Technology 2024, Volume 6, Issue 1, doi:10.21926/jept.2401006

Recommended citation: Rigatos G, Siano P, AL-Numay M, Abbaszadeh M, Cuccurullo G. Nonlinear Optimal Control for a PMLSG-VSC Wave Energy Conversion Unit. Journal of Energy and Power Technology 2024; 6(1): 006; doi:10.21926/jept.2401006.

© 2024 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

Electric power generation from wave power is a promising solution for covering the continuously growing energy needs in an environmentally friendly and without-pollution manner [1,2,3,4]. However, forces caused by waves vary in amplitude and frequency, and the voltage outputs produced by power generators receiving such actuation are not directly usable by the primary electricity grid [5,6,7,8]. To align the voltage outputs of power generators with the power grid, Voltage Source Converters (VSCs) must be connected to them [9,10,11,12]. Wave energy conversion utilizing Permanent Magnet Linear Synchronous Generators (PMLSGs) and AC/DC Voltage Source Converters (VSCs) enables efficient harnessing of wave energy and the generation of substantial electric power, which can be integrated into the electricity grid [13,14,15,16]. Controlling the motion of such a generator under wave-induced excitation is a non-trivial task that necessitates sophisticated nonlinear control algorithms [17,18,19,20]. Besides, controlling AC/DC voltage source converters to synchronize wave power units with the main electricity grid is a delicate procedure requiring advanced nonlinear control methods [21,22,23].

Additionally, control of the integrated wave power generation systems, which comprise PMLSGs serially connected to VSCs, is a subtle task that should compensate for the complicated nonlinear dynamics and multivariable structure of these energy conversion units [24,25,26,27]. Among recent results in PMLSG and VSC-based wave power conversion systems, one can distinguish methods for maximizing the amount of wave power finally converted to exploitable electric power [28,29,30,31,32]. Several other results exist on the stabilization of WECs and on the VSC-based integration of these power units to the electricity grid [33,34,35,36,37].

In this article, a novel nonlinear optimal control method is developed for the dynamic model a wave energy conversion unit that includes a Permanent Magnet Linear Synchronous Generator (PMLSG) serially connected to a three-phase AC/DC Voltage Source Converter (VSC) [38,39,40]. The dynamic model of the integrated PMLSG-VSC-based power unit is presented in state-space form, and differential flatness properties are proven about it. Differential flatness is also an implicit proof of the system’s controllability and its input-output linearizability [39]. Indeed, through successive differentiations of the power unit’s flat outputs, this system is transformed into the input-output linearized form, and a stabilizing feedback controller is developed about it using the eigenvalues assignment technique. Applying the nonlinear optimal control schemes mentioned above can help avoid the complicated and elongated state-space model transformations characteristic of this control method. In the article’s nonlinear optimal control approach, the dynamic model of the PMSLG-VSC-based wave energy conversion system undergoes first approximate linearization through Taylor series expansion [41,42,43,44]. This approximate linearization process occurs at each sampling instance around a time-varying operating point defined by the present value of the system’s state vector and the last sampled value of the control inputs vector. The linearization primarily relies on computing the system’s Jacobian matrices and updating them at each sampling interval. The modeling error arising from the truncation of higher-order terms in the Taylor series is regarded as a perturbation that is asymptotically compensated for by the robustness of the control algorithm.

A stabilizing H-infinity feedback controller is designed for the approximately linearized PMSLG-VSC-based wave energy conversion unit model. This controller solves the optimal control problem for the wave energy conversion system under model uncertainty and external perturbations. It represents a min-max differential game between (i) the control inputs, which try to minimize a quadratic function of the state vector’s tracking error, and (ii) the model uncertainty or perturbation terms, which try to maximize this cost function. To compute the controller’s gains, an algebraic Riccati equation is solved repetitively at each time step of the control method [45,46,47]. The global stability properties of this control scheme are demonstrated through Lyapunov analysis. Initially, it is established that the control loop satisfies the H-infinity tracking performance criterion, indicating significant robustness to model uncertainties and external perturbations [48,49,50].

Moreover, global asymptotic stability for this control approach is shown to hold under moderate conditions. Thus, the nonlinear optimal control method achieves fast and accurate tracking of reference setpoints by the state variables of the wave energy conversion system under moderate variations of the control inputs [38,39,40]. Finally, the H-infinity Kalman Filter is used as a robust state estimator to implement state estimation-based control without the need to measure the entire state vector of the power unit.

The nonlinear optimal control method has several advantages compared to other nonlinear control schemes one could have considered for the dynamic model of the PMLSG-VSC wave energy conversion unit: (1) unlike global linearization-based control approaches, such as Lie algebra-based control and differential flatness theory-based control, the optimal control approach does not rely on complicated transformations (diffeomorphisms) of the system’s state variables. Besides, the computed control inputs are applied directly to the initial nonlinear model of the PMSMG-VSC wave energy conversion unit and not to its linearized equivalent. The inverse transformations met in global linearization-based control are avoided, and consequently, one does not come against the related singularity problems. (2) Unlike Model Predictive Control (MPC) and Nonlinear Model Predictive Control (NMPC), the proposed control method is of proven global stability. It is recognized that MPC is a linear control approach, and if applied to the nonlinear dynamics of the wave energy conversion unit, the stability of the control loop may be compromised. Besides, in NMPC, the convergence of its iterative search for an optimum depends on initialization and parameter value selection. Consequently, the global stability of this control method cannot always be assured. (3) Unlike sliding-mode control and backstepping control, the proposed optimal control method does not require the state-space description of the system to be found in a specific form. Regarding sliding-mode control, when the controlled system is not in the input-output linearized form, defining the sliding surface can be an intuitive procedure. Regarding backstepping control, it is known that it cannot be directly applied to a dynamical system if the related state-space model is not in the triangular (backstepping integral) form. (4) unlike PID control, the proposed nonlinear optimal control method is of proven global stability; the selection of the controller’s parameters does not rely on a heuristic tuning procedure, and the stability of the control loop is assured in the case of changes of operating points, (5) unlike multiple local models-based control the nonlinear optimal control method uses only one linearization point and needs the solution of only one Riccati equation to compute the stabilizing feedback gains of the controller. Consequently, regarding computation load, the proposed control method for the wave energy conversion system’s dynamics is much more efficient.

The article’s structure is as follows: In Section 2, the dynamic model of the PMLSG-VSC-based wave energy conversion system is analyzed, and its state-space description is formulated. In Section 3, the differential flatness properties of the PMLSG-VSC wave energy conversion system are proven, and a feedback controller is designed about it using transformation in the input-output linearized form through differentiation of the flat outputs. Section 4 defines a nonlinear optimal controller for the PMLSG-VSC wave energy conversion system through approximate linearization of its state-space model using Taylor series expansion. In Section 5, the global stability properties of the nonlinear optimal control method are proven through Lyapunov stability analysis. Besides, the H-infinity Kalman Filter is used to implement state estimation-based control. In Section 6, the performance of the nonlinear optimal control method for the PMSLG-VSC-based wave energy conversion system is tested through simulation experiments. Finally, in Section 7, concluding remarks are stated.

2. Dynamic Model of the PMLSG-VSC Wave Energy Conversion Unit

An overview of previously analyzed published results on the control and optimization problem of wave energy conversion systems is provided in Table 1. The diagram of the wave energy conversion unit, consisting of a Permanent Magnet Linear Synchronous Generator (PMLSG) serially connected to an AC/DC three-phase Voltage Source Converter (VSC), is illustrated in Figure 1.

Table 1 Overview of research results on control and optimization of WECS.

Figure 1 Diagram of the wave energy conversion unit that comprises a Permanent Magnet Linear Synchronous Generator (PMLSG) serially connected to an AC/DC three-phase voltage source converter (VSC).

The motion of the moving part (mover) of the PMSLG is given by:

\[ M\ddot{z}\ = \ -{b_g}\dot{z}-{b_w}\dot{z}-{k_s}z+F_e+F_{br} \tag{1} \]

where ${b_g}\dot{z}$ is the friction force that resists the motion of the moving part of the tubular PMLSG, ${b_w}\dot{z}$ is the hydrodynamic resistance force that resists the motion of the buoy, ${k_s}z$ is a spring-type force exerted on the buoy, $F_e$ is the excitation force due to the incident wave on the buoy, $F_{br}$ is a brake force which can decelerate the mover. Additional forces that can be included in the dynamic model of the motion of the PMLSG’s mover are the Archimedean (lift) force, denoted by $F_a$, applied to the partially submerged buoy, and the radiated forces, denoted by $F_r$ acting on the buoy due to the waves generated by the buoy’s oscillations [50].

The electromagnetic force $F_{em}$ which is applied to the moving part of the PMLSG because of the generator’s magnetic field, is given by:

\[ F_{em}={3 \over 2}{\pi \over \tau_p}[{i_d}{i_q}(L_d-L_q)-\Psi_{PM}{i_q}] \tag{2} \]

where $\tau_p$ is the poles’ pitch, $i_d$ is the d-axis component of the current in the primary (mover) part of the generator, $i_q$ is the q-axis component of the current in the primary (mover) part of the generator, $L_d$ is the d-axis inductance at the primary (mover) part of the generator, $L_q$ is the q-axis inductance at the primary (mover) part of the generator, $\Psi_{PM}$ is the magnetic flux due to permanent magnets.

The dynamics of the electrical part of the Permanent Magnet Linear Synchronous Generator are obtained by decomposing the electric circuit of the generator in the dq rotating frame. The following differential equations give them [9,50]:

\[ v_d\ = \ -R{i_d}-{{d{\psi_d}} \over {dt}}+{v_t}{\pi \over \tau_p}\psi_q \tag{3} \]

\[ v_q\ = \ -R{i_q}-{{d{\psi_q}} \over {dt}}-{v_t}{\pi \over \tau_p}\psi_d \tag{4} \]

where $v_t=\dot{z}$ is the linear velocity of the primary part (mover) of the PMSLG along the z-axis, while about the magnetic flux, one has the following equations.

\[ \psi_d={L_d}{i_d}-\psi_{PM} \tag{5} \]

\[ \psi_q={L_q}{i_q} \tag{6} \]

Considering that the external forces applied to the PMLSG’s mover are: (i) The excitation force $F_e$, which arises from the incidence of waves on the buoy; (ii) The brake force $F_{br}$ which can be exerted on the mover by an electrohydraulic actuator. The dynamics of the PMLSG are then given by the following equations [9,50]:

\[ \ddot{z}\ = \ -{{(b_g+b_w)} \over M}\dot{z}-{k_v \over M}z+{3 \over {2M}}{\pi \over \tau_p}[{i_d}{i_q}(L_d-L_q)-\Psi_{PM}{i_q}]+{F_e \over M}+{F_{br} \over M} \tag{7} \]

\[ {{d{i_d}} \over {dt}}\ = \ -{R \over L_d}i_d+{\dot{z} \over L_d}{\pi \over \tau_p}({L_q}{i_q})-{v_d \over L_d} \tag{8} \]

\[ {{d{i_d}} \over {dt}}\ = \ -{R \over L_q}i_q-{\dot{z} \over L_q}{\pi \over \tau_p}({L_d}{i_d}-\Psi_{PM})-{v_q \over L_q} \tag{9} \]

Next, by defining the state variables $x_1=z$, $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_q$, and the control input variables $u_1=F_{br}$, $u_2=v_d$, $u_3=v_q$ the dynamic model of the PMLSG becomes.

\[ \begin{array} \ \dot{x}_1=x_2 \\ \\ \dot{x}_2\ = \ -{{(b_g+b_w)} \over M}x_2-{k_v \over M}x_1+{3 \over {2M}}{\pi \over \tau_p}[{x_3}{x_4}(L_d-L_q)-\Psi_{PM}{x_4}]+{F_e \over M}+{ {1\over M} {u_1}} \\ \\ \dot{x}_3\ = \ -{R \over L_d}x_3+{x_2 \over L_d}{\pi \over \tau_p}({L_q}{x_4})-{ {1\over L_d}{u_2}} \\ \\ \dot{x}_4\ = \ -{R \over L_q}x_4-{x_2 \over L_q}{\pi \over \tau_p}({L_d}{x_3}-\Psi_{PM})-{ {1\over L_q}u_3} \end{array} \tag{10} \]

The diagram of the AC/DC three-phase voltage source converter is provided in Figure 2. The following equations describe the dynamics of the AC/DC converter.

Figure 2 Diagram of the three-phase AC/DC voltage source converter.

\[ L^{'}{i}_d^{'}\ = \ -R^{'}i_d^{'}+L^{'}\omega_{dq}i_q^{'}+v_d^{'}-{V_{dc} \over 2}\eta_1 \tag{11} \]

\[ L^{'}{i}_q^{'}\ = \ -L^{'}\omega_{dq}i_d^{'}-R^{'}i_q^{'}+v_q^{'}-{V_{dc} \over 2}\eta_2 \tag{12} \]

\[ C_{dc}\dot{V}_{dc}\ = \ -{1 \over R_c}V_{dc}+{3 \over 4}i_d^{'}\eta_1+{3 \over 4}i_q^{'}\eta_2 \tag{13} \]

where $i_d^{'}$, $i_q^{'}$ are the line currents $i_a$, $i_b$, $i_c$ after transformation in the dq reference frame, $v_d^{'}$, $v_q^{'}$ are the line voltages $v_a$, $v_b$, $v_c$ after transformation in the dq reference frame, $V_{dc}$ is the voltage output of the converter, $\eta_1$ and $\eta_2$ stand for external control inputs. The resistance $R^{'}$ models the line losses, while the resistance $R_c$ represents the converter losses. The parameter $L^{'}$ denotes inductance, and the coefficient $\omega_{dq}$ represents the grid’s frequency.

Because of connecting the AC/DC converter to the PMLSG holds that $\eta_1={k_1}i_d^{'}$, $\eta_2={k_2}i_q^{'}$, where $k_1$, $k_2$ are constants. By defining the state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_7=V_{dc}$, the dynamic model of the AC/DC converter becomes [38]:

\[ \begin{array} \ \dot{x}_5\ = \ -{R^{'} \over L^{'}}x_5+\omega_{dq}x_6+{v_d^{'} \over L^{'}}-{{k_1} \over {2L^{'}}}{x_3}{x_7} \\ \\ \dot{x}_6\ = \ -\omega_{dq}x_5-{R^{'} \over L^{'}}x_6+{v_d^{'} \over L^{'}}-{{k_2} \over {2L^{'}}}{x_4}{x_7} \\ \\ \dot{x}_7\ = \ -{1 \over {R_cC_{dc}}}x_7+{{3{k_1}} \over {4C_{dc}}}{x_3}{x_5}+{{3{k_2}} \over {4C_{dc}}}{x_4}{x_6} \end{array} \tag{14} \]

Combining the dynamic model of the PMLSG given in Eq. (10) with the dynamic model of the VSC given in Eq. (14), one obtains the following state-space model for the wave energy conversion unit:

\[ \begin{array} \ \dot{x}_1=x_2 \\ \\ \dot{x}_2\ = \ -{{(b_g \ + \ b_w)} \over M}x_2-{k_v \over M}x_1+{3 \over {2M}}{\pi \over \tau_p}[{x_3}{x_4}(L_d-L_q)-\Psi_{PM}{x_4}]+{F_e \over M}+{1 \over M}u_1 \\ \\ \dot{x}_3\ = \ -{R \over L_d}x_3+{x_2 \over L_d}{\pi \over \tau_p}({L_q}{x_4})-{1 \over L_d}u_2 \\ \\ \dot{x}_4\ = \ -{R \over L_q}x_4-{x_2 \over L_q}{\pi \over \tau_p}({L_d}{x_3}-\Psi_{PM})-{1 \over L_q}u_3\\ \\ \dot{x}_5\ = \ -{R^{'} \over L^{'}}x_5+\omega_{dq}x_6+{v_d^{'} \over L^{'}}-{{k_1} \over {2L^{'}}}{x_3}{x_7} \\ \\ \dot{x}_6\ = \ -\omega_{dq}x_5-{R^{'} \over L^{'}}x_6+{v_d^{'} \over L^{'}}-{{k_2} \over {2L^{'}}}{x_4}{x_7} \\ \\ \dot{x}_7\ = \ -{1 \over {R_cC_{dc}}}x_7+{{3{k_1}} \over {4C_{dc}}}{x_3}{x_5}+{{3{k_2}} \over {4C_{dc}}}{x_4}{x_6} \end{array} \tag{15} \]

Thus, the state-space model of the wave energy conversion system can be expressed in the following matrix form:

\[\begin{pmatrix} \dot{x}_1 \\ \\ \dot{x}_2 \\ \\ \dot{x}_3 \\ \\ \dot{x}_4 \\ \\ \dot{x}_5 \\ \\ \dot{x}_6 \\ \\ \dot{x}_7 \end{pmatrix} = \begin{pmatrix} x_2 \\ \\ -{{(b_g \ + \ b_w)} \over M}x_2-{k_v \over M}x_1+{3 \over {2M}}{\pi \over \tau_p}[{x_3}{x_4}(L_d-L_q)-\Psi_{PM}{x_4}]+{F_e \over M} \\ \\ -{R \over L_d}x_3+{x_2 \over L_d}{\pi \over \tau_p}({L_q}{x_4}) \\ \\ -{R \over L_q}x_4-{x_2 \over L_q}{\pi \over \tau_p}({L_d}{x_3}-\Psi_{PM}) \\ \\ -{R^{'} \over L^{'}}x_5+\omega_{dq}x_6+{v_d^{'} \over L^{'}}-{{k_1} \over {2L^{'}}}{x_3}{x_7} \\ \\ -\omega_{dq}x_5-{R^{'} \over L^{'}}x_6+{v_d^{'} \over L^{'}}-{{k_2} \over {2L^{'}}}{x_4}{x_7} \\ \\ -{1 \over {R_cC_{dc}}}x_7+{{3{k_1}} \over {4C_{dc}}}{x_3}{x_5}+{{3{k_2}} \over {4C_{dc}}}{x_4}{x_6} \end{pmatrix} + \begin{pmatrix} 0 & 0 & 0 \\ & & & \\ {1 \over M} & 0 & 0 \\ & & & \\ 0 & -{1 \over L_d} & 0 \\ & & & \\ 0 & 0 & -{1 \over L_q} \\ & & & \\ 0 & 0 & 0 \\ & & & \\ 0 & 0 & 0 \\ & & & \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} u_1 \\ \\ u_2 \\ \\ u_3 \end{pmatrix} \tag{16} \]

Finally, the dynamic model of the wave energy conversion system can also be represented in the concise nonlinear affine-in-the-input state-space form:

\[ \dot{x}=f(x)+g(x)u \tag{17} \]

where $x\ {\in} \ R^{7{\times}1}$, $f(x)\ {\in} \ R^{7{\times}1}$, $g(x)\ {\in} \ R^{7{\times}3}$, and $u\ {\in} \ R^{3{\times}1}$.

3. Differential Flatness of the PMLSG-VSC Wave Energy Conversion Unit

3.1 Proof of Differential Flatness of the Wave Energy Conversion System

Differential flatness for a dynamical system is confirmed if: (i) All its state variables can be expressed as differential functions of a subset of selected state vector elements that constitute the flat outputs vector; (ii) The flat outputs are differentially independent, meaning they are not connected through a relation in the form of an homogeneous differential equation [39]. It will be demonstrated that the dynamic model of the PMLSG-VSC wave energy conversion unit is differentially flat with the flat outputs vector:

\[ Y=[y_1,y_2,y_3]^T=[x_1, \ {{3L^{'}} \over 4}(x_5^2+x_6^2)+{C_{dc} \over 2}x_7^2, \ x_6]^T \tag{18} \]

By differentiating the flat output $y_2$ and after some intermediate operations, one obtains:

\[ \dot{y}_2\ = \ -{{3R^{'}} \over 2}(x_5^2+x_6^2)-{1 \over R_c^{'}}{x_7^2}+{{3L^{'}} \over 2}{v_d}x_5+{{3L^{'}} \over 2}{v_q}x_6 \tag{19} \]

Besides from the definition of the flat output $y_2={{3L^{'}} \over 4}(x_5^2+x_6^2)+{C_{dc} \over 2}x_7^2$ and after solving for $x_7^2$ one has:

\[ x_7^2={{2{y_2}-{3L^{'} \over 2}(x_5^2+x_6^2)} \over C_{dc}} \tag{20} \]

By substituting Eq. (20) into Eq. (19) and using that the third flat output is $y_3=x_2$, one obtains:

\[ \dot{y}_2\ = \ -{{3R^{'}} \over 2}x_5^2-{{3R^{'}} \over 2}y_3^2-{2 \over {{R_c}{C_{dc}}}}y_2+{{3L^{'}} \over {2{R_c}{C_{dc}}}}x_5^2+ {{3L^{'}} \over {2{R_c}{C_{dc}}}}y_2^2+{{3{L^{'}}} \over 2}{v_d}{x_5}+{{3{L^{'}}} \over 2}{v_q}{y_3} \tag{21} \]

By regrouping terms in Eq. (21), one obtains a binomial for $x_5$:

\[ [-{{3R^{'} \over 2}+{{3L^{'}} \over {2{R_c}{C_{dc}}}}}]x_5^2+{{3L^{'}} \over 2}{v_d}x_5+ [-{{3R^{'} \over 2}}{y_3^2}-{2 \over {{R_c}{C_{dc}}}}y_2+{{3L^{'}} \over {2{R_c}{C_{dc}}}}y_3^2+{{3{L^{'}}} \over 2}{v_q}{y_3}]=0 \tag{22} \]

By defining coefficients $a=[-{{3R^{'} \over 2}+{{3L^{'}} \over {2{R_c}{C_{dc}}}}}]$, $b={{3L^{'}} \over 2}{v_d}$ and $c=[-{{3R^{'} \over 2}}{y_3^2}-{2 \over {{R_c}{C_{dc}}}}y_2+{{3L^{'}} \over {2{R_c}{C_{dc}}}}y_3^2+{{3{L^{'}}} \over 2}{v_q}{y_3}]$ the binomial of Eq. (22) is written as:

\[ a{x_5^2}+b{x_5}+c=0 \tag{23} \]

By solving Eq. (23) for $x_5$ and by retaining the maximum of the two solutions, one obtains:

\[ x_5={{-b+\sqrt{|b^2-4ac|}} \over {2a}}{\Rightarrow}x_5=h_5(Y,\dot{Y}) \tag{24} \]

which signifies that the state variable $x_5$ is a differential function of the flat outputs vector. Next, returning to Eq. (20):

\[ x_7=\sqrt{{{2{y_2}-{3L^{'} \over 2}(h_5(Y,\ddot{Y})^2+y_3^2)} \over C_{dc}}}{\Rightarrow}x_7=h_7(Y,\dot{Y}) \tag{25} \]

which signifies that state variable x₇ is also a differential function of the flat outputs vector. From the 5th row of the state-space model, one solves for $x_3$. It holds that:

\[ x_3\ = \ -2L^{'}\cdot{{\dot{x}_5+{R^{'} \over L^{'}}x_5-\omega_{dq}x_6-v_d} \over {{K_1}{x_2}}} \tag{26} \]

where $x_5$, $x_6$, $x_7$ are differential functions of the flat outputs vector $Y$. This gives:

\[ x_3=h_3(Y,\dot{Y}) \tag{27} \]

which signifies that $x_3$ is a differential function of the flat outputs of the system. Additionally, from the 6th row of the state-space model, one solves for $x_4$. This gives:

\[ x_4\ = \ -2L^{'}\cdot{{\dot{x}_6+\omega_{dq}x_5+{R^{'} \over L^{'}}x_6-v_q} \over {{K_2}{v_q}}} \tag{28} \]

where $x_5$, $x_6$, $x_7$ are differential functions of the flat outputs vector $Y$. This gives:

\[ x_4=h_4(Y,\dot{Y}) \tag{29} \]

which signifies that $x_4$ is a differential function of the flat outputs of the system. Furthermore, from the 1st equation of the state-space model, one has:

\[ x_2=\dot{x}_1{\Rightarrow}x_2=h_2(Y,\dot{Y}) \tag{30} \]

which signifies that $x_2$ is a differential function of the flat outputs vector. As a result, all state variables $x_i$, $i=1,2,\cdots,7$ are differential functions of the flat outputs vector. Furthermore, by solving for the control input $u_1$ from the second row of the state-space model, we obtain:

\[ u_1=M\dot{x}_2+{K_3}{x_1}+b{x_2}+{3 \over 2}{\pi \over \tau_p}[{x_3}{x_4}(L_d-L_q)-\Psi_{PM}{x_4}]-F_e \tag{31} \]

where the excitation force $F_e$ is regarded as independent of $x_i$. Therefore, $u_1$ is a differential function of the flat outputs vector or:

\[ u_1=h_{u_1}(Y,\dot{Y}) \tag{32} \]

Equivalently, from the 3rd row of the state-space model, one solves for control input $u_2$. It holds that:

\[ u_2\ = \ -[{L_d}\dot{x}_3+R{x_3}+{x_2}{\pi \over \tau_p}({L_q}{x_4})]{\Rightarrow}u_2=h_{u_2}(Y,\dot{Y}) \tag{33} \]

Finally, from the 4th row of the state-space model, one solves for control input $u_3$. It holds that:

\[ u_3\ = \ -[{L_q}\dot{x}_4+R{x_4}+{x_2}{\pi \over \tau_p}({L_d}{x_3})]{\Rightarrow}u_3=h_{u_2}(Y,\dot{Y}) \tag{34} \]

which signifies that $u_3$ is a differential function of the flat outputs vector. As a result of the above, all state variables and the control inputs are differential functions of the flat outputs vector, and the wave energy conversion system is differentially flat.

The differential flatness of the system is an implicit proof of the system’s controllability. It also demonstrates that the system is input-output linearizable through successive differentiations of its flat outputs. Finally, it allows for solving the setpoints definition problem. First, one selects setpoints in an unconstrained manner for the flat outputs of the system. Based on these setpoints, one can also obtain reference values for the rest of the state vector elements of the system using the differential relations, which connect them to the flat outputs.

3.2 Input-Output Linearization and Eigenvalues Assignment-Based Control

Establishing differential flatness properties for the dynamic model of the wave energy conversion system model allows for transforming the associated state-space model into an input-output linearized form and subsequently into the canonical Brunovsky form. Linearization is achieved through successive differentiations of the flat outputs. The state-space model of the system, as given in Eq. (15), is shown to comprise the following state equations:

\[ \begin{array} \ \dot{x}_1=f_1(x) \\ \dot{x}_2=f_2(x)+g_2(x)u_1 \\ \dot{x}_3=f_3(x)+g_3(x)u_2 \\ \dot{x}_4=f_4(x)+g_4(x)u_3 \\ \dot{x}_5=f_5(x) \\ \dot{x}_6=f_6(x) \\ \dot{x}_7=f_7(x) \end{array} \tag{35} \]

where $f_1(x)=x_2$, $g_1(x)=0$, $f_2(x)\ = \ -{{(b_g+b_w)} \over M}x_2-{k_v \over M}x_1+{3 \over {2M}}{\pi \over \tau_p}[{x_3}{x_4}(L_d-L_q)-\Psi_{PM}{x_4}]+{F_e \over M}$, $g_2(x)={1 \over M}$, $f_3(x)\ = \ -{R \over L_d}x_3+{x_2 \over L_d}{\pi \over \tau_p}({L_q}{x_4})$, $g_3(x)\ = \ -{1 \over L_d}$, $f_4(x)\ = \ -{R \over L_q}x_4-{x_2 \over L_q}{\pi \over \tau_p}({L_d}{x_3}-\Psi_{PM})$, $g_4(x)\ = \ -{1 \over L_d}$, $f_5(x)\ = \ -{R^{'} \over L^{'}}x_5+\omega_{dq}x_6+{v_d^{'} \over L^{'}}-{{k_1} \over {2L^{'}}}{x_3}{x_7}$, $g_5(x)=0$, $f_6(x)\ = \ -\omega_{dq}x_5-{R^{'} \over L^{'}}x_6+{v_d^{'} \over L^{'}}-{{k_2} \over {2L^{'}}}{x_4}{x_7}$, $g_6(x)=0$, $f_7(x)\ = \ -{1 \over {R_cC_{dc}}}x_7+{{3{k_1}} \over {4C_{dc}}}{x_3}{x_5}+{{3{k_2}} \over {4C_{dc}}}{x_4}{x_6}$, $g_7(x)=0$.

By differentiating twice the first flat output of the system $y_1=x_1$ one obtains

\[ \ddot{y}_1=f_a(x)+g_{a1}(x)u_1+g_{a2}(x)u_2+g_{a3}(x)u_3 \tag{36} \]

where about functions $f_a(x)$, $g_{a1}(x)$, $g_{a2}(x)$, $g_{a3}(x)$ one has that:

\[ f_{a1}(x)\ = \ -{{(b_g+b_w)} \over M}x_2-{k_v \over M}x_1+{3 \over {2M}}{\pi \over \tau_p}[{x_3}{x_4}(L_d-L_q)-\Psi_{PM}{x_4}]+{F_e \over M} \tag{37} \]

\[ g_{a1}(x)={1 \over M} g_{a2}(x)=0 \ g_{a3}(x)=0 \tag{38} \]

By differentiating three times flat output $y_2={{3L^{'}} \over 4}(x_5^2+x_6^2)+{C_{dc} \over 2}x_7^2$, one obtains:

\[ y_2^{(3)}=f_b(x)+g_{b1}(x)u_1+g_{b2}(x)u_2+g_{b3}(x)u_3 \tag{39} \]

where about functions $f_b(x)$, $g_{b1}(x)$, $g_{b2}(x)$, $g_{b3}(x)$ one has that:

\[ f_{b}(x)=f_{b1}(x)+f_{b2}(x)+f_{b3}(x) \tag{40} \]

with functions to be $f_{b1}(x)$, $f_{b2}(x)$, $f_{b3}(x)$, given by:

\[ \begin{array} \ f_{b1}(x)= -3R^{'}\dot{x}_5[-{R^{'} \over L^{'}}x_5+\omega_{dq}{x_6}-{k_1 \over {2L^{'}}}{x_3}{x_7}+v_d] \\ \qquad -3R^{'}{x_5}[-{R^{'} \over L^{'}}\dot{x}_5+\omega_{dq}{\dot{x}_6}-{k_1 \over {2L^{'}}}{x_3}{\dot{x}_7}] \\ \qquad \qquad -3R^{'}\dot{x}_6[-\omega_{dq}{x_5}-{R^{'} \over L^{'}}x_6-{k_2 \over {2L^{'}}}{x_4}{x_7}+v_q] \\ \qquad -3R^{'}\dot{x}_6[-\omega_{dq}{\dot{x}_5}-{R^{'} \over L^{'}}\dot{x}_6-{k_2 \over {2L^{'}}}{x_4}{\dot{x}_7}] \\ \qquad \qquad \ \ -{1 \over R_c}\dot{x}_7[-{1 \over {{R_c}{C_{dc}}}}x_7+{{3{k_1}} \over {4C_{dc}}}{x_3}{x_5}+{{3{k_2}} \over {4C_{dc}}}{x_4}{x_6}] \\ \qquad \qquad \ \ -{1 \over R_c}{x_7}[-{1 \over {{R_c}{C_{dc}}}}\dot{x}_7+{{3{k_1}} \over {4C_{dc}}}{x_3}{\dot{x}_5}+{{3{k_2}} \over {4C_{dc}}}{x_4}{\dot{x}_6}] \\ \quad {{3L^{'}} \over 2}{v_d}[-{R^{'} \over L^{'}}\dot{x}_5+\omega_{dq}\dot{x}_6-{k_1 \over {2L^{'}}}{\dot{x}_7}x_3] \\ \quad {{3L^{'}} \over 2}{v_q}[-\omega_{dq}\dot{x}_5-{R^{'} \over L^{'}}\dot{x}_6-{k_2 \over {2L^{'}}}{\dot{x}_7}x_4] \end{array} \tag{41} \]

\[ f_{b2}(x)=[3R^{'}x_5{k_1 \over {2L^{'}}}x_7-{1 \over R_c}{x_7}{{3{k_1}} \over {4C_{dc}}}x_5-{{3L^{'}} \over 2}{v_d}{k_1 \over 2L^{'}}x_7]f_3(x) \tag{42} \]

\[ f_{b3}(x)=[3R^{'}x_6{k_2 \over {2L^{'}}}x_7-{1 \over R_c}{x_7}{{3{k_2}} \over {4C_{dc}}}x_6-{{3L^{'}} \over 2}{v_q}{k_2 \over 2L^{'}}x_7]f_4(x) \tag{43} \]

while about functions $g_{b1}(x)$, $g_{b2}(x)$, $g_{b3}(x)$ one has that:

\[ \begin{array} \ g_{b1}(x)=0 \\ g_{b2}(x)=[3R^{'}x_5{k_1 \over {2L^{'}}}x_7-{1 \over R_c}{x_7}{{3{k_1}} \over {4C_{dc}}}x_5-{{3L^{'}} \over 2}{v_d}{k_1 \over 2L^{'}}x_7]g_3(x) \\ g_{b3}(x)=[3R^{'}x_6{k_2 \over {2L^{'}}}x_7-{1 \over R_c}{x_7}{{3{k_2}} \over {4C_{dc}}}x_6-{{3L^{'}} \over 2}{v_q}{k_2 \over 2L^{'}}x_7]g_4(x) \end{array} \tag{44} \]

Moreover, by differentiating twice flat output $y_3$, one obtains:

\[ \ddot{y}_3=f_c(x)+g_{c1}(x)u_1+g_{c2}(x)u_2+g_{c3}(x)u_3 \tag{45} \]

where about function $f_c(x)$, it holds that:

\[ f_{c}(x)\ = \ -\omega_{dq}\dot{x}_5-{R^{'} \over L^{'}}\dot{x}_6-{{k_1} \over {2L^{'}}}{\dot{x}_7}x_4 - {{k_1} \over {2L^{'}}}{x_7}[-{R \over L_q}x_4-{x_2 \over L_q}{\pi \over \tau_p}({L_d}{x_3})]+{x_2 \over L_q}{\pi \over \tau_p}\Psi_{PM} \tag{46} \]

while about functions $g_{c1}(x)$, $g_{c2}(x)$, $g_{c3}(x)$, one has that:

\[ \begin{array} \ g_{c1}(x)=0 \\ g_{c2}(x)=0 \\ g_{c3}(x)={k_1 \over 2L^{'}}x_7{1 \over L_q} \end{array} \tag{47} \]

In aggregate, the input-output linearized dynamics of the wave energy conversion system, which comprises a PMLSG and a VSC, is given by:

\[ \begin{equation} \begin{pmatrix} \ddot{y}_1 \\ y_2^{(3)} \\ \ddot{y}_3 \end{pmatrix}= \begin{pmatrix} f_a(x) \\ f_b(x) \\ f_c(x) \end{pmatrix}+ \begin{pmatrix} g_{a1}(x) & g_{a2}(x) & g_{a3}(x) \\ g_{b1}(x) & g_{b2}(x) & g_{b3}(x) \\ g_{c1}(x) & g_{c2}(x) & g_{c3}(x) \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix}\end{equation} \tag{48} \]

The previous state-space model in matrix form is also outlined in the following equation:

\[ \dot{\tilde{Y}}=\tilde{F}+\tilde{G}\tilde{u} \tag{49} \]

Moreover, by defining the virtual control inputs:

\[ \begin{array} \ v_1=g_{a1}(x)u_1+g_{a2}(x)u_2+g_{a3}(x)u_3 \\ v_2=g_{b1}(x)u_1+g_{b2}(x)u_2+g_{b3}(x)u_3 \\ v_3=g_{c1}(x)u_1+g_{c2}(x)u_2+g_{c3}(x)u_3 \end{array} \tag{50} \]

the following equations give the input-output linearized model of the wave energy conversion system:

\[ \ddot{y}_1=v_1 \quad y_2^{(3)}=v_2 \quad \ddot{y}_3=v_3 \tag{51} \]

By defining the state variables $z_1=y_1$, $z_2=\dot{y}_1$, $z_3=y_2$, $z_4=\dot{y}_2$, $z_5=\ddot{y}_2$, $z_6=y_3$, $z_7=\dot{y}_3$ one arrives also at a canonical Brunovsky form of the wave energy conversion system which comprises the following equations: $\dot{z}_1=z_2$, $\dot{z_2}=v_1$, $\dot{z}_3=z_4$, $\dot{z}_4=z_5$, $\dot{z}_5=v_2$, $\dot{z}_6=z_7$, and $\dot{z}_7=v_3$. In matrix form, the canonical Brunovsky state-space description of the system becomes:

\[ \begin{equation}\begin{pmatrix} \dot{z}_1 \\ \dot{z}_2 \\ \dot{z}_3 \\ \dot{z}_4 \\ \dot{z}_5 \\ \dot{z}_6 \\ \dot{z}_7 \end{pmatrix}= \begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} z_1 \\ z_2 \\ z_3 \\ z_4 \\ z_5 \\ z_6 \\ z_7 \end{pmatrix}+ \begin{pmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix}\end{equation} \tag{52} \]

with measurement equation:

\[ \begin{equation}\begin{pmatrix} z_1^m \\ z_2^m \\ z_3^m \end{pmatrix}= \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} z_1 \\ z_2 \\ z_3 \\ z_4 \\ z_5 \\ z_6 \\ z_7 \end{pmatrix}\end{equation} \tag{53} \]

The stabilizing feedback control for the wave energy conversion system is:

\[ \begin{array} \ v_1=\dot{z}_{1,d}-K_{1,1}(\dot{z}_1-\dot{z}_{1,d})-K_{2,1}(z_1-z_{1,d}) \\ v_2=z_{3,d}^{(3)}-K_{1,2}(\ddot{z}_3-\ddot{z}_{3,d})-K_{2,2}(\dot{z}_3-\dot{z}_{3,d})-K_{3,2}(z_3-z_{3,d}) \\ v_3=\ddot{z}_{6,d}-K_{1,3}(\dot{z}_6-\dot{z}_{6,d})-K_{2,3}(z_6-z_{6,d}) \end{array} \tag{54} \]

By substituting the feedback control law of Eq. (54) into the state-space model of Eq. (52), one obtains the following closed-loop system dynamics:

\[ \begin{array} \ (\ddot{z}_1-\ddot{z}_{1,d})+K_{1,1}(\dot{z}_1-\dot{z}_{1,d})+K_{2,1}(z_1-z_{1,d})=0 \\ (z_3^{(3)}-z_{3,d}^{(3)})+K_{1,2}(\ddot{z}_3-\ddot{z}_{3,d})+K_{2,2}(\dot{z}_3-\dot{z}_{3,d})+K_{3,2}(z_3-z_{3,d})=0 \\ (\ddot{z}_6-\ddot{z}_{6,d})+K_{1,3}(\dot{z}_6-\dot{z}_{6,d})+K_{2,3}(z_6-z_{6,d})=0 \end{array} \tag{55} \]

Next, by defining the tracking error variables $e_1=z_1-z_{1,d}$, $e_3=z_3-z_{3,d}$, $e_6=z_6-z_{6,d}$, one obtains that the tracking error dynamics of the closed-loop system is given by:

\[ \begin{array} \ \ddot{e}_1+K_{1,1}\dot{e}_1+K_{2,1}e_1=0 \\ e_3^{(3)}+K_{1,2}\ddot{e}_3+K_{2,2}\dot{e}_3+K_{3,2}e_3=0 \\ \ddot{e}_6+K_{1,3}\dot{e}_6+K_{2,3}e_6=0 \end{array} \tag{56} \]

Besides, by selecting the feedback gains $(K_{1,1},K_{2,1})$, $(K_{1,2},K_{2,2},K_{3,2})$ and $(K_{1,3},K_{2,3})$ so as the characteristic polynomials which are associated with the differential equations of Eq. (56) to be Hurwitz stable, one has that:

\[ \begin{array} \ lim_{t{\rightarrow}\infty}e_1=0 \ {\Rightarrow} \ lim_{t{\rightarrow}\infty}z_1=z_{1,d} \ {\Rightarrow} \ lim_{t{\rightarrow}\infty}y_1=y_{1,d} \\ lim_{t{\rightarrow}\infty}e_3=0 \ {\Rightarrow} \ lim_{t{\rightarrow}\infty}z_3=z_{3,d} \ {\Rightarrow} \ lim_{t{\rightarrow}\infty}y_2=y_{2,d} \\ lim_{t{\rightarrow}\infty}e_6=0 \ {\Rightarrow} \ lim_{t{\rightarrow}\infty}z_6=z_{6,d} \ {\Rightarrow} \ lim_{t{\rightarrow}\infty}y_3=y_{3,d} \end{array} \tag{57} \]

After computing the control inputs vector $\tilde{v}=[v_1,v_2,v_3]^T$, one can also determine the actual control inputs vector $\tilde{u}$ of Eq. (49), which should be applied to the initial nonlinear state-space model of the wave energy conversion system. The inverse transformation gives this:

\[ \tilde{u}=\tilde{G}^{-1}(\tilde{v}-\tilde{F}) \tag{58} \]

By demonstrating that the flat outputs of the system converge to their associated setpoints it is also ensured that all state variables$x_i(t)$ converge to their reference values. However, as shown above, control of the PMLSG and VSC-based wave energy conversion system through input-output linearization and eigenvalues assignment is a complicated procedure requiring complex forward and inverse transformations. A nonlinear optimal control approach is applied to overcome this complexity.

4. Nonlinear Optimal Control for the Wave Energy Conversion System

4.1 Approximate Linearization of the Dynamics of the Wave Energy Conversion System

The dynamic model of the wave energy conversion system undergoes approximate linearization around the temporary operating point $(x^{*},u^{*})$, which is updated at each sampling instant and is defined by the present value of the system’s state vector $x^{*}$ and by the last sampled value of the control inputs vector $u^{*}$. The linearization process relies on the first-order Taylor series expansion and the computation of the corresponding Jacobian matrices. The modeling error, resulting from the truncation of higher-order terms in the Taylor series expansion, is regarded as a perturbation that is asymptotically compensated for by the robustness of the control algorithm [38,40].

Through the proposed linearization process, the wave energy conversion system, being initially in the nonlinear state-space form $\dot{x}=f(x)+g(x)u$, is turned into the equivalent linear state-space description:

\[ \dot{x}=Ax+Bu+\tilde{d} \tag{59} \]

where $A$, $B$ are the Jacobian matrices of the system, and $\tilde{d}$ is the disturbances vector, which may comprise (i) the modeling error due to the truncation of higher-order terms in the Taylor series expansion, (ii) exogenous perturbations, (iii) sensor measurement noise of any distribution. The Jacobian matrices of the system are given by:

\[ A={\nabla_x}[f(x)+g(x)u]\mid_{(x^{*},u^{*})}{\Rightarrow}A={\nabla_x}f(x)\mid_{(x^{*},u^{*})} \tag{60} \]

\[ B={\nabla_u}[f(x)+g(x)u]\mid_{(x^{*},u^{*})}{\Rightarrow}B=g(x)\mid_{(x^{*},u^{*})} \tag{61} \]

The linearization approach adopted for implementing the nonlinear optimal control scheme yields a highly accurate system dynamics model. Consider, for instance, the following affine-in-the-input state-space model.

\[ \begin{eqnarray} && \qquad \qquad \qquad \qquad\dot{x}=f(x)+g(x)u{\Rightarrow}\\ \dot{x}&=&[f(x^{*})+{\nabla_x}f(x)\mid_{x^{*}}(x-x^{*})]+[g(x^{*})+{\nabla_x}g(x)\mid_{x^{*}}(x-x^{*})]u^{*} \\ && +g(x^{*})u^{*}+g(x^{*})(u-u^{*})+\tilde{d}_1{\Rightarrow}\\ \dot{x}&=&[{\nabla_x}f(x)\mid_{x^{*}}+{\nabla_x}g(x)\mid_{x^{*}}u^{*}]x+g(x^{*})u \\ &&-[{\nabla_x}f(x)\mid_{x^{*}}+{\nabla_x}g(x)\mid_{x^{*}}u^{*}]x^{*}+f(x^{*})+g(x^{*})u^{*}+\tilde{d}_1 \end{eqnarray} \tag{62} \]

where $\tilde{d}_1$ is the modeling error due to the truncation of higher order terms in the Taylor series expansion of $f(x)$ and $g(x)$. Next, by defining $A=[{\nabla_x}f(x)\mid_{x^{*}}+{\nabla_x}g(x)\mid_{x^{*}}u^{*}]$, $B=g(x^{*})$, one obtains:

\[ \dot{x}=Ax+Bu-Ax^{*}+f(x^{*})+g(x^{*})u^{*}+\tilde{d}_1 \tag{63} \]

Moreover, by denoting $\tilde{d}\ = \ -Ax^{*}+f(x^{*})+g(x^{*})u^{*}+\tilde{d}_1$ about the cumulative modeling error term in the Taylor series expansion, one has

\[ \dot{x}=Ax+Bu+\tilde{d} \tag{64} \]

which is the approximately linearized model of the dynamics of the system of Eq. (59). The term $f(x^{*})+g(x^{*})u^{*}$ is the derivative of the state vector at $(x^{*},u^{*})$, which is almost annihilated by $-Ax^{*}$.

Next, the computation of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$ proceeds as follows:

First row of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$: ${{{\partial}f_1} \over {{\partial}x_1}}=0$, ${{{\partial}f_1} \over {{\partial}x_2}}=1$, ${{{\partial}f_1} \over {{\partial}x_3}}=0$, ${{{\partial}f_1} \over {{\partial}x_4}}=0$, ${{{\partial}f_1} \over {{\partial}x_5}}=0$, ${{{\partial}f_1} \over {{\partial}x_6}}=0$, ${{{\partial}f_1} \over {{\partial}x_7}}=0$

Second row of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$: ${{{\partial}f_2} \over {{\partial}x_1}}\ = \ -{k_3 \over M}$, ${{{\partial}f_2} \over {{\partial}x_2}}\ = \ -{b \over M}$, ${{{\partial}f_2} \over {{\partial}x_3}}={3 \over {2M}}{\pi \over \tau_p}{x_4}(L_d-L_q)$, ${{{\partial}f_2} \over {{\partial}x_4}}={3 \over {2M}}{\pi \over \tau_p}[x_3(L_d-L_q)-\Psi_{PM}]$, ${{{\partial}f_2} \over {{\partial}x_5}}=0$, ${{{\partial}f_2} \over {{\partial}x_6}}=0$, ${{{\partial}f_2} \over {{\partial}x_7}}=0$

Third row of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$: ${{{\partial}f_3} \over {{\partial}x_1}}=0$, ${{{\partial}f_3} \over {{\partial}x_2}}\ = \ -{1 \over L_d}{\pi \over \tau_p}L_q$, ${{{\partial}f_3} \over {{\partial}x_3}}\ = \ -{R \over L_d}$, ${{{\partial}f_3} \over {{\partial}x_4}}={x_2 \over L_d}{\pi \over \tau_p}L_q$, ${{{\partial}f_3} \over {{\partial}x_5}}=0$, ${{{\partial}f_3} \over {{\partial}x_6}}=0$, ${{{\partial}f_3} \over {{\partial}x_7}}=0$

Fourth row of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$: ${{{\partial}f_4} \over {{\partial}x_1}}=0$, ${{{\partial}f_4} \over {{\partial}x_2}}\ = \ -{1 \over L_d}{\pi \over \tau_p}(L_d{x_3}-\Psi_{PM})$, ${{{\partial}f_4} \over {{\partial}x_3}}\ = \ -{x_2 \over L_d}{\pi \over \tau_p}L_d$, ${{{\partial}f_4} \over {{\partial}x_4}}\ = \ -{R \over L_d}$, ${{{\partial}f_4} \over {{\partial}x_5}}=0$, ${{{\partial}f_4} \over {{\partial}x_6}}=0$, ${{{\partial}f_4} \over {{\partial}x_7}}=0$

Fifth row of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$: ${{{\partial}f_5} \over {{\partial}x_1}}=0$, ${{{\partial}f_5} \over {{\partial}x_2}}=0$, ${{{\partial}f_5} \over {{\partial}x_3}}\ = \ -{{{k_1}{x_7}} \over {2L^{'}}}$, ${{{\partial}f_5} \over {{\partial}x_4}}=0$, ${{{\partial}f_5} \over {{\partial}x_5}}\ = \ -{R^{'} \over L^{'}}$, ${{{\partial}f_5} \over {{\partial}x_6}}=\omega_{dq}$, ${{{\partial}f_5} \over {{\partial}x_7}}\ = \ -{{{k_1}{x_3}} \over {2L^{'}}}$

Sixth row of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$: ${{{\partial}f_6} \over {{\partial}x_1}}=0$, ${{{\partial}f_6} \over {{\partial}x_2}}=0$, ${{{\partial}f_6} \over {{\partial}x_3}}=0$, ${{{\partial}f_6} \over {{\partial}x_4}}\ = \ -{{{k_2}{x_7}} \over {2L^{'}}}$, ${{{\partial}f_6} \over {{\partial}x_5}}\ = \ -\omega_{dq}$, ${{{\partial}f_6} \over {{\partial}x_6}}\ = \ -{R^{'} \over L^{'}}$, ${{{\partial}f_6} \over {{\partial}x_7}}\ = \ -{{{k_2}{x_4}} \over {2L^{'}}}$

Seventh row of the Jacobian matrix ${\nabla_x}f(x)\mid_{(x^{*},u^{*})}$: ${{{\partial}f_7} \over {{\partial}x_1}}=0$, ${{{\partial}f_7} \over {{\partial}x_2}}=0$, ${{{\partial}f_7} \over {{\partial}x_3}}={3 \over {4C_{dc}}}{k_1}{x_5}$, ${{{\partial}f_7} \over {{\partial}x_4}}={3 \over {4C_{dc}}}{k_2}{x_6}$, ${{{\partial}f_7} \over {{\partial}x_5}}= {3 \over {4C_{dc}}}{k_1}{x_3}$, ${{{\partial}f_7} \over {{\partial}x_6}}={3 \over {4C_{dc}}}{k_2}{x_4}$, ${{{\partial}f_7} \over {{\partial}x_7}}\ = \ -{1 \over {{R_c}C_{dc}}}$.

4.2 Equivalent Linearized Dynamics of the Wave Energy Conversion System

After linearization around its current operating point, the dynamic model for the PMLSG-VSC-based wave energy conversion system is expressed as:

\[ \dot{x}=Ax+Bu+d_1 \tag{65} \]

Parameter $d_1$ stands for the linearization error in the PMLSG-VSC-based wave energy conversion system that was given previously in Eq. (59). The reference setpoints for the state vector of the aforementioned dynamic model are denoted by ${{x_d}}=[x_1^{d},\cdots,x_7^{d}]$. After applying the control input $u^{*}$, this setpoint vector is tracked. At every time-instant, the control input $u^{*}$ is assumed to differ from the control input u appearing in Eq. (65) by an amount equal to ${\Delta}u$, that is, $u^{*}=u+{\Delta}u$.

\[ \dot{x}_d=Ax_d+Bu^{*}+d_2 \tag{66} \]

The dynamics of the controlled system described in Eq. (65) can also be expressed as:

\[ \dot{x}=Ax+Bu+Bu^{*}-Bu^{*}+d_1 \tag{67} \]

Moreover, by denoting $d_3\ = \ -Bu^{*}+d_1$ as an aggregate disturbance term, one obtains the following form:

\[ \dot{x}=Ax+Bu+Bu^{*}+d_3 \tag{68} \]

By subtracting Eq. (66) from Eq. (68), one has:

\[ \dot{x}-\dot{x}_d=A(x-x_d)+Bu+d_3-d_2 \tag{69} \]

By denoting the tracking error as $e=x-x_d$ and the aggregate disturbance term as $L\tilde{d}=d_3-d_2$, the tracking error dynamics becomes:

\[ \dot{\pmb{e}}=Ae+Bu+L\tilde{d} \tag{70} \]

where $L$ is the disturbance inputs gain matrix. After applying an H-infinity feedback control scheme, the above-linearized form of the PMLSG-VSC-based wave energy conversion system can be efficiently controlled.

4.3 The Nonlinear H-Infinity Control

The initial nonlinear model of the PMLSG-VSC-based wave energy conversion system is in the form:

\[ \dot{x}=f(x,u) \ \ x\ {\in} \ R^n, \ u\ {\in} \ R^m \tag{71} \]

The model of the PMLSG-VSC-based wave energy conversion system is linearized at each iteration of the control algorithm around its current operating point ${(x^{*},u^{*})}=(x(t),u(t-T_s))$. The linearized equivalent of the system is described by:

\[ \dot{x}=Ax+Bu+L\tilde{d} \ \ x\ {\in} \ R^n, \ u\ {\in} \ R^m, \ \tilde{d}\ {\in} \ R^q \tag{72} \]

where matrices A and B are obtained from the computation of the previously defined Jacobians and vector $\tilde{d}$ denotes disturbance terms due to linearization errors. At the same time, $L$ is a disturbance inputs gain matrix. The problem of disturbance rejection for the linearized model is described by:

\[ \begin{array} \ \dot{x}=Ax+Bu+L\tilde{d} \\ y=Cx \end{array} \tag{73} \]

where $x\ {\in} \ R^n$, $u\ {\in} \ R^m$, $\tilde{d}\ {\in} \ R^q$ and $y\ {\in} \ R^p$ cannot be handled efficiently if the classical $LQR$ control scheme is applied. This is because of the existence of the perturbation term $\tilde{d}$. The disturbance term $\tilde{d}$ apart from modeling (parametric) uncertainty and external perturbation terms, can also represent noise terms of any distribution.

In the $H_{\infty}$ control approach, a feedback control scheme is designed for setpoints tracking by the system’s state vector and simultaneous disturbance rejection, considering that the disturbance affects the system in the worst possible manner. The effects of disturbances are incorporated in the following quadratic cost function:

\[ J(t)={1 \over 2}{\int_0^T}[{y^T}(t)y(t)+r{u^T}(t)u(t)-{\rho^2}{\tilde{d}^T}(t)\tilde{d}(t)]dt, \ \ r,{\rho}>0 \tag{74} \]

The significance of the negative sign in the cost function’s term that is associated with the perturbation variable $\tilde{d}(t)$ is that the disturbance tries to maximize the cost function $J(t)$ while the control signal $u(t)$ tries to minimize it. The physical meaning of the abovementioned relation is that the control signal and the disturbances compete within a min-max differential game. This problem of min-max optimization can be written as ${min_{u}}{max_{\tilde{d}}}J(u,\tilde{d})$.

The objective of the optimization procedure is to compute a control signal $u(t)$ which can compensate for the worst possible external disturbance imposed on the PMLSG-VSC-based wave energy conversion system. It can be observed that the solution to the min-max optimization problem is directly dependent on the value of the parameter $\rho$. This implies an upper bound exists on the magnitude of disturbances the control signal can nullify.

4.4 Computation of the Feedback Control Gains

For the linearized system described by Eq. (73), the cost function defined in Eq. (74) is characterized, where the coefficient r determines the penalization of the control input, and the weight coefficient ρ determines the reward for the effects of disturbances. It is assumed that (i) The energy that is transferred from the disturbances signal $\tilde{d}(t)$ is bounded, that is ${\int_0^{\infty}}{\tilde{d}^T(t)}\tilde{d}(t){dt}<\infty$, (ii) matrices $[A, B]$ and $[A, L]$ are stabilizable, (iii) matrix $[A, C]$ is detectable. In the case of a tracking problem, the optimal feedback control law is given by:

\[ u(t)\ = \ -Ke(t) \tag{75} \]

with $e=x-x_d$ as the tracking error and $K={1 \over r}{B^T}P$ where $P$ is a positive definite symmetric matrix. As demonstrated in Section 5, matrix $P$ is obtained from the solution of the Riccati equation:

\[ {A^T}P+PA+Q-P({2 \over r}B{B^T}-{1 \over {\rho^2}}L{L^T})P=0 \tag{76} \]

where $Q$ is a positive semi-definite symmetric matrix, the worst-case disturbance is given by:

\[ \tilde{d}(t)={1 \over \rho^2}{L^T}Pe(t) \tag{77} \]

The solution of the H-infinity feedback control problem for the PMLSG-VSC-based wave energy conversion system and the computation of the worst-case disturbance that the related controller can sustain comes from the superposition of Bellman’s optimality principle when considering that this renewable energy unit is affected by two separate inputs (i) the control input $u$ (ii) the cumulative disturbance input $\tilde{d}(t)$. Solving the optimal control problem for $u$, that is, for the minimum variation (optimal) control input that achieves elimination of the state vector’s tracking error, gives $u\ = \ -{1 \over r}{B^T}Pe$. Equivalently, solving the optimal control problem for $\tilde{d}$, for the worst-case disturbance that the control loop can sustain gives $\tilde{d}={1 \over \rho^2}{L^T}Pe$.

The diagram of the control loop for the PMLSG-VSC-based wave energy conversion system is illustrated in Figure 3.

Figure 3 Diagram of the control scheme for the PMLSG-VSC-based wave energy conversion system.

5. Lyapunov Stability Analysis

5.1 Stability Proof

Through Lyapunov stability analysis it will be demonstrated that the proposed nonlinear control scheme guarantees $H_{\infty}$ tracking performance for the PMLSG-VSC-based wave energy conversion system and achieves asymptotic convergence to the reference setpoints in the presence of bounded disturbance terms. The dynamics of the tracking error for the PMLSG-VSC-based wave energy conversion system are expressed in the form [38,40].

\[ \dot{\pmb{e}}=Ae+Bu+L\tilde{d} \tag{78} \]

where in the PMLSG-VSC based wave energy conversion unit’s case $L\ {\in} \ R^{7{\times}7}$ is the disturbance inputs gain matrix. The variable $\tilde{d}$ represents model uncertainties and external disturbances in the wave energy conversion system’s model. The following Lyapunov equation is considered.

\[ V={1 \over 2}{e^T}Pe \tag{79} \]

where $e=x-x_d$ is the tracking error. By differentiating concerning time, one obtains:

\[ \begin{array} \ \dot{V}={1 \over 2}{\dot{e}^T}Pe+{1 \over 2}{e^T}P\dot{e}{\Rightarrow}\\ \dot{V}={1 \over 2}{[Ae+Bu+L\tilde{d}]^T}Pe+{1 \over 2}{e^T}P[Ae+Bu+L\tilde{d}]{\Rightarrow} \\ \end{array} \tag{80} \]

\[ \begin{array} \ \dot{V}={1 \over 2}[{e^T}{A^T}+{u^T}{B^T}+{\tilde{d}^T}{L^T}]Pe +{1 \over 2}{e^T}P[Ae+Bu+L\tilde{d}]{\Rightarrow} \\ \end{array} \tag{81} \]

\[ \dot{V}={1 \over 2}{e^T}{A^T}Pe+{1 \over 2}{u^T}{B^T}Pe+{1 \over 2}{\tilde{d}^T}{L^T}Pe+ {1 \over 2}{e^T}PAe+{1 \over 2}{e^T}PBu+{1 \over 2}{e^T}PL\tilde{d} \tag{82} \]

The previous equation is rewritten as:

\[ \dot{V}={1 \over 2}{e^T}({A^T}P+PA)e+({1 \over 2}{u^T}{B^T}Pe+{1 \over 2}{e^T}PBu)+({1 \over 2}{\tilde{d}^T}{L^T}Pe+{1 \over 2}{e^T}PL\tilde{d}) \tag{83} \]

Assumption: For given positive definite matrix $Q$ and coefficients $r$ and $\rho$ there exists a positive definite matrix $P$, which is the solution of the following matrix equation.

\[ {A^T}P+PA\ = \ -Q+P({2 \over r}B{B^T}-{1 \over \rho^2}L{L^T})P \tag{84} \]

Moreover, the following feedback control law is applied to the system.

\[ u\ = \ -{1 \over r}{B^T}Pe \tag{85} \]

By substituting Eq. (84) and Eq. (85) in $\dot{V}$, one obtains:

\[ \dot{V}={1 \over 2}{e^T}[-Q+P({2 \over r}B{B^T}-{1 \over {\rho^2}}L{L^T})P]e+{e^T}PB(-{1 \over r}{B^T}Pe)+{e^T}PL\tilde{d}{\Rightarrow} \tag{86} \]

\[ \dot{V}\ = \ -{1 \over 2}{e^T}Qe+({1 \over r}{e^T}PB{B^T}Pe-{1 \over {2\rho^2}}{e^T}PL{L^T}Pe)-{1 \over r}{e^T}PB{B^T}Pe+{e^T}PL\tilde{d} \tag{87} \]

which after intermediate operations gives:

\[ \dot{V}\ = \ -{1 \over 2}{e^T}Qe-{1 \over {2\rho^2}}{e^T}PL{L^T}Pe+{e^T}PL\tilde{d} \tag{88} \]

or, equivalently:

\[ \dot{V}\ = \ -{1 \over 2}{e^T}Qe-{1 \over {2\rho^2}}{e^T}PL{L^T}Pe+{1 \over 2}{e^T}PL\tilde{d}+{1 \over 2}{\tilde{d}^T}{L^T}Pe \tag{89} \]

Lemma: The following inequality holds [38,47].

\[ {1 \over 2}{e^T}L\tilde{d}+{1 \over 2}\tilde{d}{L^T}Pe-{1 \over {2\rho^2}}{e^T}PL{L^T}Pe\ {\leq}\ {1 \over 2}{\rho^2}{\tilde{d}^T}\tilde{d} \tag{90} \]

Proof: The binomial $({\rho}{\alpha}-{1 \over \rho}b)^2$ is considered. Expanding the left part of the above inequality one gets.

\[ \begin{array} \ {\rho^2}{\alpha^2}+{1 \over {\rho^2}}{b^2}-2ab \geq 0 \Rightarrow {1 \over 2}{\rho^2}{\alpha^2}+{1 \over {2\rho^2}}{b^2}-ab \geq 0 \Rightarrow\\ ab-{1 \over {2\rho^2}}{b^2} \leq {1 \over 2}{\rho^2}{\alpha^2} \Rightarrow {1 \over 2}ab+{1 \over 2}ab-{1 \over {2\rho^2}}{b^2} \leq {1 \over 2}{\rho^2}{\alpha^2} \end{array} \tag{91} \]

The following substitutions are carried out: $a=\tilde{d}$ and $b={e^T}{P}L$ and the previous relation becomes:

\[ {1 \over 2}{\tilde{d}^T}{L^T}Pe+{1 \over 2}{e^T}PL\tilde{d}-{1 \over {2\rho^2}}{e^T}PL{L^T}Pe\ {\leq}\ {1 \over 2}{\rho^2}\tilde{d}^T\tilde{d} \tag{92} \]

Eq. (92) is substituted in Eq. (89) and the inequality is enforced, thus giving [38,47].

\[ \dot{V}\ {\leq}\ -{1 \over 2}{e^T}Qe+{1 \over 2}{\rho^2}{\tilde{d}^T}\tilde{d} \tag{93} \]

Eq. (93) shows that the $H_{\infty}$ tracking performance criterion is satisfied. The integration of $\dot{V}$ from 0 to $T$ gives:

\[ \begin{array} \ {\int_0^T}\dot{V}(t)dt\ {\leq}\ -{1 \over 2}{\int_0^T}{||e||_Q^2}{dt}+{1 \over 2}{\rho^2}{\int_0^T}{||\tilde{d}||^2}{dt}{\Rightarrow}\\ 2V(T)+{\int_0^T}{||e||_Q^2}{dt}\ {\leq}\ 2V(0)+{\rho^2}{\int_0^T}{||\tilde{d}||^2}dt \end{array} \tag{94} \]

Moreover, if there exists a positive constant $M_d>0$ such that:

\[ \int_0^{\infty}{||\tilde{d}||^2}dt \leq M_d \tag{95} \]

then one gets:

\[ {\int_0^{\infty}}{||e||_Q^2}dt \leq 2V(0)+{\rho^2}{M_d} \tag{96} \]

Thus, the integral ${\int_0^{\infty}}{||e||_Q^2}dt$ is bounded. Moreover, $V(T)$ is bounded, and from the definition of the Lyapunov function $V$ in Eq. (79), it becomes clear that $e(t)$ will also be bounded since $e(t) \ \in \ \Omega_e=\{e|{e^T}Pe\ {\leq}\ 2V(0)+{\rho^2}{M_d}\}$. According to the above and using Barbalat’s Lemma, one obtains $lim_{t \rightarrow \infty}{e(t)}=0$.

After following the stages of the stability proof, one arrives at Eq. (93), which shows that the H-infinity tracking performance criterion holds. By selecting the attenuation coefficient $\rho$ to be sufficiently small and, in particular, to satisfy $\rho^2<||e||^2_Q / ||\tilde{d}||^2$ one has that the first derivative of the Lyapunov function is upper bounded by 0. This condition holds at each sampling instance; consequently, global stability for the control loop can be concluded.

5.2 Robust State Estimation with the Use of the $H_{\mathrm{\infty }}$ Kalman Filter

The control loop has to be implemented using information provided by a few sensors and by processing only a small number of state variables. There is no need to measure the speed of the moving part of the PMSLG. A filtering scheme is proposed to reconstruct the missing information about the state vector of the wave energy conversion system. Subsequently, based on this scheme, state estimation-based control can be applied [38,47]. By denoting as $A(k)$, $B(k)$, and $C(k)$ the discrete-time equivalents of matrices $A$, $B$, and $C$ of the linearized state-space model of the system, the recursion of the $H_{\infty}$ Kalman Filter for the model of the PMLSG-VSC wave energy conversion unit, can be formulated in terms of a measurement update and a time update part.

Measurement update:

\[ \begin{array} \ D(k)=[I-{\theta}W(k)P^{-}(k)+{C^T}(k)R(k)^{-1}C(k)P^{-}(k)]^{-1}\\ K(k)=P^{-}(k)D(k){C^T}(k)R(k)^{-1}\\ \hat{x}(k)=\hat{x}^{-}(k)+K(k)[y(k)-C\hat{x}^{-}(k)] \end{array} \tag{97} \]

Time update:

\[ \begin{array} \ \hat{x}^{-}(k+1)=A(k)x(k)+B(k)u(k) \\ P^{-}(k+1)=A(k)P^{-}(k)D(k)A^T(k)+Q(k) \end{array} \tag{98} \]

where it is assumed that parameter $\theta$ is sufficiently small to assure that the covariance matrix $P{^{-}(k)}^{-1}-{\theta}W(k)+C^T(k)R(k)^{-1}C(k)$ will be positive definite. When $\theta=0$ the $H_{\infty}$ Kalman Filter becomes equivalent to the standard Kalman Filter. One can measure only a part of the state vector of the PMLSG-VSC-based wave energy conversion system and estimate it by filtering the rest of the state vector elements. For instance, one can measure state variables associated with the PMLSG $x_1=z$, $x_3=i_d$, $x_4=i_q$, and state variables associated with the VSC $x_5=i_d^{'}$, $x_6=i_q^{'}$, and $x_7=V_{dc}$ and can estimate through filtering the speed of the PMLSG’s mover $x_2=\dot{z}$.

6. Simulation Tests

Simulation experiments further confirm the global stability properties of the proposed nonlinear optimal control method for the PMLSG-VSC-based wave energy conversion system and demonstrate the precise tracking of setpoints achieved for the state variables of the wave energy conversion system. To implement the article’s nonlinear optimal control approach, the algebraic Riccati equation of Eq. (84) has to be repetitively solved at each time step of the control algorithm.

The sampling period has been chosen to be $T_s=0.01$ sec. The parameters of this Riccati equation are gains $r$, $\rho$ and matrix $Q$. The method is computationally efficient because even at a PC with moderate processing power (multi-core i7 Intel processor at 2.8 GHz), one achieves the solution of the Riccati equation within a time-interval which is much shorter than the sampling period. Figures 4 to 19 present the state variables of the wave energy conversion unit in blue, highlight the associated setpoints in red, and show the estimates generated by the H-infinity Kalman Filter in green. It can be noted that the proposed nonlinear optimal control scheme achieves fast and accurate tracking of setpoints by the state variables of the PMLSG-VSC-based wave energy conversion system under moderate variations of the control inputs.

The parameters that define the transient performance of the nonlinear optimal control method are gains $r$, $\rho$, and $Q$. For relatively small values of $r$, the steady-state tracking error of the state variables is eliminated. Additionally, for relatively large values of the elements of the diagonal matrix $Q$, faster convergence of the state variables to their reference values is achieved. Moreover, the smallest value of the attenuation coefficient $\rho$, for which one still obtains a valid solution of the method’s Riccati equation, is the one that provides the control loop with maximum robustness.

Figure 4 Tracking of setpoint 1 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$, of the VSC to their reference setpoints.

Figure 5 Tracking of setpoint 1 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

Figure 6 Tracking of setpoint 2 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$ of the VSC to their reference setpoints.

Figure 7 Tracking of setpoint 2 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

Figure 8 Tracking of setpoint 3 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$ of the VSC to their reference setpoints.

Figure 9 Tracking of setpoint 3 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

Figure 10 Tracking of setpoint 4 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$ of the VSC to their reference setpoints.

Figure 11 Tracking of setpoint 4 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

Figure 12 Tracking of setpoint 5 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$ of the VSC to their reference setpoints.

Figure 13 Tracking of setpoint 5 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

Figure 14 Tracking of setpoint 6 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$ of the VSC to their reference setpoints.

Figure 15 Tracking of setpoint 6 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

Figure 16 Tracking of setpoint 7 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$ of the VSC to their reference setpoints.

Figure 17 Tracking of setpoint 7 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

Figure 18 Tracking of setpoint 8 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) convergence of state variables $x_2=\dot{z}$, $x_3=i_d$, $x_4=i_d$ of the PMLSG to their reference setpoints (red line: setpoint, blue line: real value, green line: estimated value), (b) convergence of state variables $x_5=i_d^{'}$, $x_6=i_q^{'}$, $x_{7}=V_{dc}$ of the VSC to their reference setpoints.

Figure 19 Tracking of setpoint 8 for the PMLSG-VSC based wave energy conversion system under nonlinear optimal control (a) control inputs $u_1$ to $u_3$ applied to the wave energy conversion unit, (b) tracking error $e_1$, $e_2$, $e_3$, for state variables $x_2=\dot{z}$, $x_3=i_d$, $x_7=V_{dc}$ of the wave energy conversion unit.

To elaborate on the tracking performance and on the robustness of the proposed nonlinear optimal control method for the PMLSG and VSC-based wave energy conversion unit the following Tables are given: (i) Table 2 which provides information about the accuracy of tracking of the reference setpoints by the state variables of the wave power conversion unit’s state-space model, (ii) Table 3 which provides information about the robustness of the control method to parametric changes in the model of the wave power conversion unit’s dynamics (change ${\Delta}a\%$ in the mass $M$ of the moving part of the PMLSG), (iii) Table 4 which provides information about the precision in state variables’ estimation that the H-infinity Kalman Filter achieves, (iv) Table 5 which provides the convergence times of the wave power conversion unit’s state variables to the associated setpoints.

Table 2 Tracking RMSE for the WEC system in the disturbance-free case.

Table 3 Tracking RMSE for the WEC system in the case of disturbances.

Table 4 RMSE × 10⁻³ for the estimation performed by the H-infinity KF.

Table 5 Convergence time (sec) for the WEC system’s state variables.

The nonlinear optimal control method, analyzed in the present article, is novel compared to past attempts to solve the nonlinear optimal control problem for renewable energy systems. Unlike past approaches, in the new nonlinear optimal control method, linearization is performed around a temporary operating point, which is defined by the present value of the system’s state vector and by the last sampled value of the control inputs vector and not at points that belong to the desirable trajectory (setpoints). Besides, the Riccati equation used for computing the feedback gains of the controller is new, and so is the global stability proof for this control method. Compared to NMPC (Nonlinear Model Predictive Control), which is a popular approach for treating the optimal control problem in electrical machines and drives, the new nonlinear optimal (H-infinity) control scheme is of proven global stability and the convergence of its iterative search for the optimum does not depend on initial conditions and ad-hoc parametrization of the controller. Notably, the nonlinear optimal control method is more applicable to a broader class of dynamic systems than approaches based on the solution of State Dependent Riccati Equations (SDRE). The application of SDRE approaches is limited to dynamic systems that can be transformed into the Linear Parameter Varying (LPV) form. Besides, the nonlinear optimal control method performs better than nonlinear optimal control schemes, which use the approximation of the solution of the Hamilton-Jacobi-Bellman equation by Galerkin series expansions. The stability properties of the Galerkin series expansion-based optimal control approaches are still unproven. Finally, it is worth noting that the results presented in the article have implications for industrial applications and hold the potential for utilization in various electric power generation schemes.

7. Conclusions

Wave energy conversion is a promising approach for raising the electricity rates produced from renewable energy sources. To this end, the present article has proposed a novel nonlinear optimal control method for the dynamic model of wave energy conversion systems, which comprises a Permanent Magnet Linear Synchronous Generator (PMLSG) serially connected with a three-phase AC/DC Voltage Source Converter (VSC). It has been proven that the dynamic model of the PMLSG-VSC-based wave energy unit is differential flat, which stands for an implicit proof of the system’s controllability. Moreover, through successive differentiations of the system’s flat outputs, the associated dynamic model has been written in the input-output linearized form, and a stabilizing feedback controller has been designed about it using the eigenvalues assignment concept. In order to circumvent the complex and protracted state-space model transformations associated with the control above method, the article introduces a nonlinear optimal control method for the PMLSG-VSC-based power unit. The latter method does not need complex transformations of the state-space description or any changes in state variables. Conversely, it guarantees rapid and precise tracking of reference setpoints by all wave energy conversion system state variables, even under moderate variations of the wave energy conversion unit.

In the nonlinear optimal control method, the dynamic model of the PMLSG-VSC-based wave energy conversion unit undergoes approximate linearization with the use of first-order Taylor series expansion and through the computation of the associated Jacobian matrices. The linearization process takes place at each sampling instance around a temporary operating point, defined by the present value of the system’s state vector and by the last sampled value of the control inputs vector. A stabilizing H-infinity feedback controller is designed for the approximately linearized system model. To select the controller’s feedback gains, an algebraic Riccati equation has to be repetitively solved at each time step of the control algorithm. The global stability properties of the control scheme are proven through Lyapunov stability analysis. Initially, it is shown that the control loop meets the H-infinity tracking performance criterion, indicating its robustness to model uncertainty and external perturbations. Furthermore, it is established that global asymptotic stability conditions are also satisfied under moderate conditions. The H-infinity Kalman Filter has been used as a robust state estimator to implement sensorless feedback control without measuring the entire state vector of the PMLSC-VSC-based energy conversion system. The satisfactory performance of the nonlinear optimal control approach for the PMLSG-VSC power unit has been further confirmed through simulation experiments.

Acknowledgments

(a) G. Rigatos has been partially supported by Grant Ref. 301022 “Nonlinear optimal and flatness-based control methods for complex dynamical systems” of the Unit of Industrial Automation, Industrial Systems Institute, Greece (b) Pierluigi Siano and Mohammed AL-Numay acknowledge financial support from the Researchers Supporting Project Number (RSP2023R150), King Saud University, Riyadh, Saudi Arabia.

Author Contributions

The authors’ contribution in this manuscript is given by the order of appearance of their names in the article’s list of authors. G. Rigatos had the initial concept of the article and has contributed to the methodology, to the formulation of the analytical part, to the development of software code and to the validation of results. P. Siano, M. Al-Numay, M. Abbaszadeh and G. Cuccurullo have contributed to the article’s methodology, to the confirmation of the analytical part and to validation of results.

Competing Interests

The authors have declared that no competing interests exist.