A Fast Model Predictive Control Framework for Multi-Float and Multi-Mode-Motion Wave Energy Converters

Recently, multi-float and multi-mode-motion wave energy converters (M-WECs) have been developed to improve energy conversion capability. Although model predictive control (MPC) can be very effective to solve the constrained energy maximization control problem of point absorber WECs, the increased complexity of the M-WEC hydrodynamics can bring significant challenges due to computational demand. This brief proposes a novel computational-efficient fast MPC (FMPC) design method for the M-WECs requiring complex linear hydrodynamic models. The controller design objective is to maximize the energy conversion with some available wave forecasting information and to satisfy state and control input constraints to ensure safe operation. The main advantage of the proposed FMPC is the reduced computational burden with a negligible impact on performance. A demonstrative numerical simulation based on a 1:50 laboratory-scale M-WEC design, M4, for which linear hydrodynamics has been verified experimentally, is presented to verify the efficacy of the proposed control method in terms of both computational load and energy output.


z s
States associated with radiation force dynamics. P out Converted power considering losses. E out Converted energy considering losses. H s Significant wave height. T p Peak period. n p Wave prediction horizon. n c Fast MPC control horizon.

I. INTRODUCTION
O CEAN waves provide a substantial renewable energy source, with the potential capacity of 2 terawatts (TWs) worldwide [1]. Despite the tremendous potential, wave energy technology is still relatively immature for commercialization compared with other renewable energy sources, e.g., wind energy and solar energy, mainly due to its high unit electricity cost [2]. To bring down the cost, it is necessary to improve the energy conversion efficiency. This can be achieved by the combination of a good hydrodynamic design and a well-designed associated wave energy converter (WEC) controller.
Various types of point absorber design have been developed, and their associated control policies have been extensively studied [3], [4], [5]. However, wave energy capture of point absorbers in regular waves is limited to wave power (power transmitted per meter crest width) times wavelength divided by 2π in heave, twice this for surge and pitch and three times for heave and pitch and/or surge in combination [6]. For a two beam raft system, the limit is 4/3π times wavelength [7].
To improve energy conversion, in particular to make capacity similar to wind turbines, multi-float and multi-mode-motion WECs (M-WECs) have been developed, such that energy absorbed from different modes of motion from each float combines constructively [8], [9]. A reliable and efficient controller is essential to exploit the potential of an M-WEC. Early conventional control policies based on the impedance matching principle are unsuitable for M-WECs, because they cannot optimally handle the hydrodynamic complexity and the safety and machine constraints that may come from design configurations and the limited capacity of power takeoff (PTO) mechanisms [4]. Those constraints are essential to guarantee the safe operation of the WEC PTO, by limiting power or torque [10], [11].
The M-WEC controller design objective is to maximize the energy output while subject to safety constraints and power and/or torque limitations. Recent studies [4], [12] reveal that the WEC control problem is essentially a constrained optimal control problem (OCP), and with short-term forecasting of incoming waves, the energy conversion efficiency can be remarkably improved [12]. Being able to handle constraints optimally, many model predictive control (MPC)-based algorithms have been proposed for point absorber WECs [12], [13], [14]. However, the computational burden can increase significantly with the increase in the WEC modeling complexity. On the one hand, the performance of the WEC MPC algorithms heavily relies on the accuracy of the control-oriented model of an M-WEC. If unaccounted for in the WEC control, the modeling inaccuracies can significantly degrade the control performance [15], [16]. On the other hand, it is challenging to develop a simple model of an M-WEC that represents the WEC dynamics in complex real sea conditions with adequate fidelity. In particular, it is shown in [17] that the model order of an M-WEC that represents a class of attenuator type is substantial, generally at least in the order of hundreds of states. The existing WEC MPC algorithms, when applied to the M-WEC requiring a complicated model, are unable to be computed in real time and are, thus, ineffective.
This brief aims to solve the problem by proposing a computational-efficient fast MPC (FMPC) method for M-WECs requiring high-dimensional linear models. It is important to note that this discussion on control relates to linear hydrodynamics, and in practice, some point absorbers and rafts show nonlinear effects, particularly in Froude-Krylov and drag forces, e.g., [18]. In contrast, the M-WEC M4 developed by The University of Manchester, Manchester, U.K., has shown linear hydrodynamic response characteristics in both operational and surprisingly extreme waves as demonstrated by direct analysis of experimental data [19] and linear diffraction modeling [9], [20]. Note that this linear prediction does not extend to mooring force [20] though not of concern here. A picture of M-WEC M4 [9] with the three-float configuration is shown in Fig. 1. Linear hydrodynamic modeling makes the control problem more tractable, and nonlinear effects generally reduce power capture. Since control to increase power capture generally increases response, it is important that linear modeling also predicts large response well.
The main idea of the proposed FMPC is to move part of the control calculation offline to reduce the online computational burden. This is achieved by first solving the unconstrained OCP offline with an analytic solution and then designing an MPC to cope with constraints. The concept of augmenting an offline-designed causal control law with constraint-handling mechanism has been verify for the control of power absorber WEC in our previous work [21]. However, the proposed FMPC method for M-WEC significantly differs from [21] in the following aspects. First, a noncausal control law is solved offline for the unconstrained OCP, enabling the explicit incorporation of wave forecasting information to notably improve the energy conversion efficiency. Second, the MPC is designed with a much shorter control horizon to reduce the computational burden while retaining the ability to incorporate a more extended wave prediction horizon to avoid performance degradation.
This brief is organized as follows. The control problem is formulated in Section II. In Section III, first, an MPC strategy for M-WECs is designed following conventional approaches as a comparison basis. Then, the novel FMPC is developed. Some implementation issues are also addressed. Numerical simulations based on a laboratory-scale M-WEC M4 is given in Section IV to demonstrate the efficacy. Finally, this brief is concluded in Section V.
Notation: R n (R a×b ) denotes the space of all real n-dimensional vectors (a-by-b-dimensional matrices). I [a,b] denotes a set of integers from a to b. For column vectors a and b, [a, b] denotes a column vector [a T b T ] T . n p and n c denote the disturbance prediction horizon and control horizon, respectively. Subscripts k and i denote the time steps and prediction steps, respectively. I n denotes the n-by-n identity matrix. 0 a×b denotes an a-by-b matrix whose all elements are 0. "s.t." is the abbreviation of "subject to."

II. CONTROL PROBLEM SETUP
To demonstrate the efficacy of the proposed FMPC, we present our study based on the M-WEC M4 with 1-1-1 configuration shown in Figs. 1 and 2. Here, we briefly explain the working principle of this attenuator type M-WEC. The persistent wave excitation provides kinetic energy to the hinge arms, which leads to the flexing motions about the hinge. The motions then drive a PTO mechanism located at the hinge to generate electricity. Since an attenuator type M-WEC operates parallel to the wave direction and effectively rides the waves and the roll motion of the device can be minimized by some mechanical designs, we study the vertical planar motion. For more details, please refer to [8], [9], and [17]. The model of the M-WEC is in the generalized coordinates q ∈ R 4 , which are selected as follows: Here, x O denotes the hinge position in the surge direction; z O denotes the hinge position in the heave direction; θ L denotes the angular position of the left arm, and θ R denotes the Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
angular position of the right arm. Here, the directions of the generalized coordinates are defined in Fig. 2. Following the similar modeling procedure as [17], we obtain the following: where M is the 4 × 4 mass and inertia matrix; f e,q (t), f rd,q (t), and f rs,q (t) are the wave excitation force, radiation damping force, and hydrostatic restoring force in generalized coordinates, respectively. f pto,q (t) is the PTO torque, which acts as the control input. Drag force is negligible for the M4 WEC supported by computational fluid dynamics result [22].
Mooring has a negligible influence on response and energy conversion. This has been demonstrated experimentally by tethering from a buoy or from a beam above water [23] and by modeling including or excluding a buoy mooring [20]. A mooring model is, thus, omitted here.
where z s is the state associated with radiation force dynamics, control input u(t) := f pto,q (t) ∈ R, and disturbance input w(t) := f e,q (t) ∈ R n w , the final state-space representation of the control-oriented model for M-WEC can be written as follows:ẋ Here, the continuous-time state-space matrices (A c , B uc , B wc ) adopt the same value as [17]. n x and n w are the dimension of state and disturbance, respectively. Please refer to [17, Sec. II] for the details. To facilitate the MPC design, the continuous-time control-oriented model is discretized via zero-order holder (ZOH) and a sampling time t s , which leads to where (A, B u , B w ) are the corresponding discrete-time state-space matrices. With the direction of M-WEC modeling defined in Fig. 2, the instant converted power, with the consideration of the quadratic dissipation effect in the PTO [24], can be expressed by the following: The energy captured, using the control in a ZOH convention, for a time interval between time steps k 1 and k 2 , can be expressed by the following: where To guarantee the safe operation of the M-WEC, and to reduce the maintenance cost, we define the following safety constraints on the angular positions of the left and right arms: where θ L ,max and θ R,max are the position limits of left and right arms, respectively. Due to the capacities of the PTO mechanism, an important constraint on control input torque provided by the PTO is imposed, which can be expressed as follows: where u max is the PTO torque limit. At each time instant k, an n p -step long wave forecasting information w k is assumed to be available, which is represented as follows: where w k , w k+1 , . . . , w k+n p −1 are the wave excitation force [with respect to the generalized coordinates (1)] at time instant k, k + 1, . . . , k + n p − 1, respectively. In summary, the design objective of the M-WEC MPC is to maximize the energy output (5) while satisfying the safety constraints (6). In Section III, we will focus on proposing numerical tractable algorithms to solve this control problem.

A. Conventional MPC Formulation
Based on the availability of an n p -step wave prediction w k , an MPC law can be derived by recursively solving the following constrained OCP: where R, S x , S w , C L , and C R are defined in (5), and then apply the first element of optimal solution u * k as the control input to the system.
Remark 1: The result of [25] shows that for a passive system, (8) leads to convex quadratic programming (QP). The WEC systems typically satisfy this passivity assumption. This is because if such an assumption is violated, it will imply that the WEC can convert more energy than it absorbs from wave excitation.
Remark 2: Although the QP can be formulated following the standard procedure as in [10], the computational burden is still relatively high for real-time implementation due to the increased complexity of the M-WEC model. As a result, only a reduced control horizon can be used to avoid the computational complexity problems, leading to performance degradation. Thus, a trade-off can be inevitable between computational burden and performance with the conventional MPC design approach for M-WECs with complex linear dynamics. In Section III-B, a computational-efficient algorithm will be developed to tackle OCP (8) to avoid this trade-off.

B. Novel FMPC Formulation
The design procedure of the proposed FMPC can be described in two steps.
1) Find the analytic solution of OCP (8) without considering the constraints, which leads to a linear time-variant control law whose coefficients are determined offline. 2) Design an MPC to track the calculated unconstrained optimal control law, considering constraints (6). First, we consider the following unconstrained OCP. Theorem 1: The unconstrained OCP of (8) can be solved by µ k := [μ k , μ k+1 , . . . , μ k+n p −1 ] for i ∈ I [0,n p −1] , where the time-variant coefficients K x,i , K w,i , and K s,i are determined by the following: and the coefficients P i+1 and s i+1 are calculated from backward iterations with final conditions P n p = 0 n x ×n x and s n p = 0 n x ×1 , provided that R + B T u P i+1 B u ≥ 0 for all i ∈ I [0,n p −1] . Here, the coefficients i := (A + B u K x,i ) T and i := i P i+1 B w + K T x,i S w . Proof: Define the optimal cost-to-go function (from prediction step i to n p ) as follows: for i ∈ I [0,n p ] , with stage cost (x k+i , u k+i ) := (1/2)Ru 2 k+i + u k (S x x k+i + S w w k+i ). Since the WEC model (4) is linear time-invariant and the stage cost is quadratic, following [26] and [27], a guess of the value of the optimal cost-to-go where P i ∈ R n x ×x x , s i ∈ R n x ×1 , and a i ∈ R are the coefficients to be determined. By comparing (13) with (14), we have the final conditions P n p = 0 n x ×n x , s n p = 0 n x ×1 , and a n p = 0.
With the Bellman optimality principle [28], we have the following: The second-order condition of (15) gives From the first-order condition, the optimal law μ satisfies which leads to (10) and (11).
To determine the value of P i+1 in (11) and s i+1 in (10), we substitute optimal control action μ k+i and optimal cost-togo function v(x, i ) in (15) with (10) and (14), respectively. Since (15) satisfies for all x k+i , both the coefficients of terms that are quadratically depending on and linearly depending on x k+i are 0, which leads to (12). The proof is completed.
Remark 3: Theorem 1 solves the unconstrained OCP (9) with an analytic form solution. The calculation of unconstrained optimal control law can be separated into two parts: 1) offline calculation of the coefficients P i+1 , K x,i , K w,i , and K s,i via (12b) and 2) online calculation of s i+1 via (12b). To reduce the computational burden, we reinvestigate the online calculation of s i+1 .
Corollary 1: The calculation of the unconstrained optimal control law (10) can be further simplified as follows: for i ∈ I [0,n p −1] . Here, w k is an n p -step long wave prediction at time instant k; the feedback coefficients K x,i take the same value as calculated from (11), and the feedforward coefficients K d,i are calculated from where K x,i , K w,i , and K s,i are calculated from (11), and i and i are calculated from (12), respectively. Proof: From (12b) and final condition s n p = 0, we have s n p −1 = n p −1 w k+n p −1 s n p −2 = n p −2 s n p −1 + n p −2 w k+n p −2 = n p −2 n p −1 w k+n p −1 + n p −2 w k+n p −2 . . .
Remark 4: Corollary 1 provides a simplified calculation of unconstrained optimal control law μ, where the coefficients K x,i and K d,i are computed offline. Theorem 1 and Corollary 1 provide a simple linear controller that can explicitly use the wave preview information to improve the energy output. However, the proposed linear controller can lose its efficacy when the safety constraints violation becomes a prominent issue. To explicitly handling constraints, we consider the following FMPC framework: (20) where μ k+i is the unconstrained optimal control law calculated from (17), and v k+i is an online optimization variable to be determined. The working principle of the FMPC is to use the first term μ k+i to maximize the energy and the additional term v k+i to cope with constraints. In summary, the FMPC recursively solves the following OCP at each time instant k: where , with n c as the control horizon chosen to be very short to reduce the computational burden. Coefficients K x,i , K d,i , u max , θ L ,max , and θ R,max take the same value as in Corollary 1. With the solution of OCP (21), the control action to be applied at time instant k is determined by the following: where v * k is the first element of v * k . Remark 5: The proposed FMPC differs from the MPC designed following the conventional approach, because the following hold.
1) In the FMPC formulation (21), the control horizon n c can be chosen to be much shorter than the wave prediction horizon, i.e., n c n p , to substantially reduce the computational burden, while using a much prolonged wave prediction to improve control performance. 2) In the conventional MPC formulation (8), the control horizon has be chosen the same as the wave prediction horizon to use the information of wave prediction, i.e., n c = n p . The predicted state trajectory x k := [x k , x k+1 , . . . , x k+n c −1 ] can be calculated from for . . , w k+n p −1 ], and A K ,i := A + BK x,i ; the coefficients L i ∈ R n x ×n c and C i ∈ R n x ×n p n w are defined by The predicted input trajectory u k := [u k , u k+1 , . . . , u k+n c −1 ] can be calculated from The safety constraints on state and input (6) can be rewritten as follows: where here, we design a Luenberger observer to retrieve the full information of states, which takes the following form: Here,x denotes the estimated state information that will be passed to the QP of FMPC; y k = Cx k is the measured output.
In this case, we assume that the surge and heave positions and the angular positions of the arms are directly measurable, i.e., C := [I 4 0 4×(n x −4) ], where (A, C) is assumed to be observable. The observer gain L is properly designed with ρ(A − LC) < 0, such that the real state information x k can be estimated withx k with an acceptable error e k = x k −x k . For more details of designing the observer gain, please refer to [26]. Design Procedure: The proposed FMPC can be designed following the steps below.

Algorithm 1 Implementation of the Proposed FMPC
The implementation of the FMPC is summarized in Algorithm 1, and the implementation framework is shown in Fig. 3.

IV. NUMERICAL SIMULATIONS
In this section, numerical simulations are carried out here based on a laboratory prototype of M-WEC M4 to verify the efficacy of the proposed FMPC method. The laboratory scale is about 1:50. Note that for Froude scaling, which is standard for hydrodynamic modeling, a length scale l gives a time scale of √ l, a force scale of l 3 , a power scale of l 3.5 , and so on. As reported in [17], the control-oriented model of M-WEC M4 has the model order of 136, which stimulates the demand of the proposed FMPC method. The simulations adopt similar settings as [17], where the sampling time is t s = 0.009 s. The PTO torque limit is (31a) and the PTO parasitic loss is modeled by 3 × 10 −3 u 2 k . For safe operational purposes, the maximal angular position limits are θ L ,max = 12 degs, θ R,max = 12 degs. (31b) The simulation sections are organized in the following order. First, after demonstration of the design procedure of the proposed method, we show that the safety constraints on state and control input (6) can be ensured. Second, the comparative simulation results are provided between FMPC and MPC designed following Section III-A and [10], referred as "CMPC." The results verify that the FMPC can achieve the maximal energy production and significantly reduce the computational load.

A. Simulation Set I: Time Simulation
The FMPC is designed following the design procedure detailed in Section III-C. The wave prediction length of 4.5 s (n p = 500) is used. To reduce the computational burden, the FMPC is designed with a short control horizon (n c = 5).
The first time simulation is presented based on a 100-s Joint North Sea Wave Project (JONSWAP) wave profile [29] with a significant wave height H s = 0.03 m and a peak period T p = 1.2 s. Fig. 4 shows the control input responses using FMPC. We can see that when constraints are active, the online optimized auxiliary variable v provides compensation, such that the input constraints violations are avoided. The maximal amplitude of θ L or θ R is 10.5 • , which indicates that the state constraints are also effectively handled for the 100-s simulation time. This demonstrates the working principle of the online constraint-handling mechanism of the proposed FMPC.
To demonstrate the universal efficacy of the proposed FMPC in different sea conditions, the comprehensive simulation results are provided based on the JONSWAP wave spectra, with a fixed significant wave height H s = 0.03 m and peak periods ranging from T p = 0.4 to 2.0 s, each with 0.2-s interval. Fig. 5 shows the time percentages when the auxiliary control input part v = 0, i.e., constraints are active, in different   sea conditions. Here, a large percentage indicates strong constraints handling effect from FMPC. For the sea conditions with peak periods between T p = 1.2 s and T p = 2.0 s, constraints are active for 20%-30% of the simulation time, indicating that the constraints handling mechanism of FMPC contributes significantly to ensuring the safe operation of the M-WEC in those sea states. Fig. 6 shows the average energy conversion powers and corresponding capture width ratios (CWRs) in different sea conditions. The CWR is defined as average power absorbed divided by wave power per meter crest width divided by wavelength associated with frequency at the centroid of wave spectrum. The sea states with peak period ranging from 0.8 to 1.8 s have much higher average powers and CWRs, representing most effective operating sea conditions for M-WEC M4. By comparing Fig. 6 with Fig. 5, we find that the effect of constraints handling with FMPC is aligned with the average power output, which means that the constraints handling mechanism of FMPC is more effective in the sea states where higher wave power is captured.

B. Simulation Set II: Computational Time Analysis
In simulation set II, the proposed FMPC is compared with CMPC to further demonstrate the advantages of the proposed FMPC in reducing computational loads. The simulations are run on a standard PC with Intel Core i5 Processor at 2.6 GHz and 16-GB memory, with the same sampling rate t s = 0.009 s as used in simulation set I and [17]. The simulation presented in Fig. 7 is based on the same 100-s JONSWAP wave profile with significant wave height H s = 0.03 m and peak period T p = 1.2 s as used in simulation set I.
Since the CMPC can be designed with different prediction horizons, the red line in Fig. 7 (top) shows the average time spent computing a single step and average power output for CMPC with different prediction horizons. It can be concluded that with longer prediction, longer wave prediction information is used, which results in a better energy conversion efficacy. However, the increase of prediction length in CMPC will bring computational issues. When the prediction exceeds 1.8 s, the CMPC cannot yield its solution within a time interval. Moreover, when the prediction horizon is longer than 4.5 s, the PC cannot carry out the corresponding optimization due to computational space issues.
The blue lines in Fig. 7 show the average time spent on computing a single step and average power output when using the proposed FMPC with different prediction horizons. The control horizon is set as n = 5. We can see that the average time spent in calculating the MPC does not have a significant change, while the average power output increases notably with the use of longer wave prediction information. Overall, compared with CMPC, the computational load of FMPC using long wave prediction information is dramatically reduced, which means with the proposed FMPC, we can use longer prediction to get better energy conversion efficiency. With the affordable computational load, CMPC can only use 1.8 s of wave prediction, while the FMPC can make the best use of the 6.8-s wave prediction, which leads to over 22% of energy conversion efficiency, in the ideal case.
However, wave prediction inevitably comes with inaccuracies. To test the robustness of FMPC subject to prediction error, we deliberately add wave prediction sequences with ramped error increasing with time, as shown in Fig. 8. This Top: snapshot of imperfect wave excitation predictions, in the surge direction. Bottom: average power output of FMPC using 100% accurate prediction (gray); imperfect prediction with underestimated magnitudes (red); imperfect prediction with overestimated magnitudes (blue). is in alignment with the feature of the state-of-the-art wave prediction techniques, which provide more accurate prediction for shorter horizons than for longer horizons [30]. We can see that the performances of the three cases are similar when the prediction discrepancy is minor and efficiency is slightly degraded when the prediction error is significant. In this case, a 4.8-s (four times the peak period equivalently) prediction length is an appropriate choice for FMPC, representing a 21% performance improvement compared with CMPC using a 1.8-s wave excitation prediction.
V. CONCLUSION In this brief, a novel computational-efficient FMPC is proposed, which tackles a critical problem of optimal controller design for M-WECs. A prerequisite is that the hydrodynamics may be assumed to be linear. After solving the unconstrained energy maximization control problem offline with an analytic solution, we design an MPC to cope with constraints, with a short control horizon to reduce the computational burden yet a long wave prediction horizon to improve efficiency.
Numerical simulations based on M-WEC M4 verify that the FMPC can effectively reduce the computational load compared with the MPC designed following the conventional approach. The reduced computational requirements make the proposed FMPC suitable for M-WECs that require high-dimensional linear models. Our future work will be focused on the real-time implementation of FMPC on a prototype of M-WEC M4 in wave basins.