Observer-Based Unknown Input Estimator of Wave Excitation Force for a Wave Energy Converter

Several energy maximization control approaches for point-absorber wave-energy converter (PAWEC) systems require knowledge of the wave excitation force (WEF) that is not measurable during the PAWEC operation. Many WEF estimators have been proposed based on stochastic PAWEC modeling using the Kalman filter (KF), extended KF (EKF), or receding-horizon estimation. Alternatively, a deterministic WEF estimator is proposed here based on the fast unknown input estimation (FUIE) concept. The WEF is estimated as an unknown input obviating the requirement to represent its dynamics. The proposed observer-based unknown input estimator (OBUIE) inherits the capability of estimating fast-changing signals from the FUIE, which is important when considering irregular wave conditions. Unlike preceding methods, the OBUIE is designed based on a PAWEC model, including the nonlinear viscous drag force. It has been shown that the nonlinear viscous drag force is essential for accurate PAWEC model description, within the energy maximization control role. The performance of the proposed estimator is evaluated in terms of PAWEC conversion efficiency in a single degree-of-freedom PAWEC device operating in regular and irregular waves. Simulation results are obtained using MATLAB to evaluate the estimator under different control methods and subject to parametric uncertainty.

prediction. However, WEF is not directly measurable during 48 PAWEC operation and hence it should be calculated/estimated from other available or redundant PAWEC measurements. 50 Several studies have proposed WEF estimators based on the Kalman Filter (KF) [5], [6], [7], [8] and the Extended 52 Kalman Filter (EKF) [9]. In general, good estimation of WEF is obtained using the KF-based approaches utilizing 54 measurement of PAWEC output(s) and input. However, the KF-based WEF estimators necessitate the WEF dynamics to 56 be represented in the estimator model. The WEF dynamics are usually approximated by a set of Harmonic Oscillators with 58 a range of frequencies that should be specified based on the incoming wave, which is unknown. The number of Harmonic 60 Oscillator frequencies and their values affect the observer accuracy as highlighted by [6] [8]. In addition, the model order 62 is noticeably increased when using KF-Based methods due to augmentation of extra states for the Harmonic Oscillator 64 with the original system states. This adds complexity especially when considering complex PAWEC cases, e.g. multi 66 Degree-Of-Freedom (DOF) PAWECs or PAWEC arrays. The representation of WEF dynamics using harmonic oscillators 68 is unnecessarily complex and other approaches are proposed [10]. These authors describe a receding-horizon approach to 70 estimate the WEF utilizing an iterative solution of a quadratic programming problem shown to give accurate estimation of 72 WEF validated against experimental data. [10] also proposed an estimator using a KF coupled with a Random-Walk WEF 74 model to overcome the complexity of harmonic oscillator WEF representation. The Random-Walk approach shows a 76 comparable performance to the receding-horizon approach but with less complexity and computation burden. 78 As in most scientific studies two approaches to physical modelling usually exist, namely stochastic and deterministic 80 methods. In the context of WEF estimation the former includes the KF as well as the receding horizon approaches discussed 82 above. Here deterministic modelling methods are all based on consideration of non-linear dynamics and uncertainty. For 84 WEF estimation this requires an understanding of the nonlinear viscous force effect (see Section II for description) which 86 is shown to be important [11] for accurate PAWEC model utilized for energy maximization control. 88 Motivated by the importance of viscous force consideration for model-based energy maximization control and the 90 fact that most of WEF estimator proposed so far have not explicitly consider this important force, this paper proposes a 92 deterministic WEF estimator utilizing the fast unknown input estimation concept. It is worth mentioning that deterministic 94 methods [12], [13] and the stochastic receding-horizon method [10] are based on unknown input estimation as well, but they 96 all assume linear model without the viscous force term. Since a vertical cylindrical buoy is considered, f b (t) is proportional to the displacement p(t) [15] where K m is the hydrostatic stiffness. K m = 710 for the considered PAWEC [15]. The viscous force is given by [17] f where v w (t) is the water surface velocity near the heaving buoy, assumed measured, and K v is a constant that depends 140 on the buoy geometry and buoy water interaction, for the considered PAWEC K v = 50.96 [11]. The PAWEC used here 142 is considered a slender structure according to [17]. So that the viscous force is essential in the accurate hydrodynamic 144 description, see [11] for more explanation. The time domain radiation force is given by: where M ∞ is the added mass at infinite frequency, M ∞ = 6.5Kg for the considered PAWEC [11]. K r (t) is the so called impulse response function, or the kernel function, of the radiation force. The term in (4) represents a convolution integral between K r (t) and 146 v(t) and can be approximated by a finite order state-space subsystem [18], [19] which makes the system model more 148 suitable for control design. Hence, [14] uses the Matlab function imp2ss to directly deduce the following state-space 150 model from the impulse response function K r (t) obtained using the boundary element method software NEMOH: A 3 rd order model is chosen in [14] leading to: Using (1)-(6) the PAWEC hydrodynamics can be represented in state space form as: f ul (t) represents the unknown and nonlinear part of the model. Section III describes an observer-based technique to estimate f ul (t).

III. OBSERVER-BASED UNKNOWN INPUT ESTIMATOR
In this Section an approach for the estimation of f e (t) is presented, inspired by the FUIE [20], [21] strategy. This technique is used to estimate f ul (t) as the unknown input signal in (8).

A. System Model
Consider the model (8), where the state vector  The PAWEC model (8) with matrices (9) satisfies conditions 1-3 directly (checked using Matlab). Moreover, the force 176 summation f ul (t) satisfies condition 4, as the wave motion is continuous.

178
Remark 1: To facilitate the OBUIE design it is assumed that both the PAWEC position and velocity are available, and this 180 explains the choice of matrix C as in (9). This could appear as a limitation of the proposed observer. However, it is important 182 to measure the PAWEC displacement, for example using a linear variable differential transducer (LVDT). Nevertheless, 184 the velocity measurement is usually expensive to obtain. Here the velocity is estimated by a soft sensing using a band-pass 186 filtered version of the low frequency position signal, since the band-pass filter is a band-limited differentiator with derivative 188 action being most accurate at low frequencies [12]. This approach also has the advantage of filtering higher frequency 190 noise effects. An interesting alternative is to use an inertial measurement unit which can be a cost effective way to obtain 192 v(t), see [22], for PAWEC practical implementation. However, in a real application both would be valuable to achieve suitable 194 redundancy to enhance the PAWEC reliability in the face of inertial measurement unit malfunction.

B. Observer-based WEF estimator
The following observer is proposed for the system (8) wherex(t) ∈ R 5 is the observer state estimate,f ul (t) is the estimate of the unknown input,ŷ(t) ∈ R 2 is the observer output, and e y (t) :=ŷ (t) − y(t) . The observer gain L 1 is designed such thatĀ = A−L 1 C is stable. Define the following error signals: e x (t):=x (t) −x(t) and e f (t) :=f ul (t)−f ul (t), then the error dynamics are described bẏ The following learning law, inspired from [21], is used to get f ul (t):f where Γ > 0 is a learning rate, chosen by the user, and L 2 198 is a design gain. Lemma 1 [23]: Given a scalar α > 0 and a symmetric positive definite matrix G, the following inequality holds: Theorem 1: Under conditions 1-4 above and given positive scalars µ and σ, if there exist symmetric positive definite matrices P ∈ R 5×5 and G ∈ R 1 , such that the following constraints hold: and where then the observer (11) and the update law (13) give the uniformly ultimately bounded estimation errors e x (t) and e f (t).
Proof: To prove the stability of the error system (12), and hence the estimation convergence, assume the following Lyapunov function then proceed to prove thatV (x,f ul , t) < 0. Using (12) and 200 (13) the derivative of (16) can be written aṡ Using Lemma 1, it follows that wheref ≥ ḟ ul (t) satisfies condition 4. Substituting (18) 202 into (17), it follows that: and with Since B ul is full column rank, when Π < 0, then (19) can

206
Based on Lyapunov stability theory the estimation errors e x and e f are ultimately bounded, i.e. they converge to a small } around the origin. Actually, it is not possible to achieve asymptotic conver-210 gence, i.e. e x and e f convergence to 0, unless the excitation force is constant or piecewise constant, which imply its 212 first derivative is 0 valued. It is then possible to prove the asymptotic convergence using the Barbalat's Lemma. As the 214 WEF cannot be assumed constant or even piecewise constant, it is reasonable to ensure the estimation errors converge to a 216 small set Ψ. Practically Ψ is close to zero and can be further tightened towards the origin by increasing the learning rate 218 Γ which in turn decreases ∆. The small size of Ψ gives confidence of the effectiveness of the estimator to get accurate 220 enough results. This is the case as we can see from the simulation results.

222
Remark 2: The equality constraint (14) is proposed to simplify the proof and ensure that the optimization problem can be feasible. However, it is difficult to solve equality constraints such as (14) using the Matlab LMI toolbox. Following the idea presented in [24], (14) is replaced by where η is a positive constant chosen by the designer to be small enough to enable a good "approximation" of the equality 224 (14) with the inequality (24). This procedure is used by a number of investigators for cases of mixed equality/inequality 226 problems [20], [21].
After obtaining an estimate of f ul (t) using the OBUIE, the 228 WEF is calculated by subtracting (3) from (10). This requires knowledge of the water surface velocity, v w (t), which can be 230 obtained by differentiation of the water surface elevation measurement. The differentiation operation can be approximated 232 using a filtered differentiation of the water surface elevation measurement, approximated by a suitable band-pass filter 234 (see [12] and Remark 1). An alternative application approach would be to obtain v w (t) using an inertial measurement unit. 236 The solution of the LMI constraints set (24) and (15) is obtained using the Matlab LMI toolbox, with parameters 238 µ = 2, η = 1 × 10 −4 and σ = 100. The resulting observer gains are: The learning rate Γ > 0 is chosen by the designer. Section IV discusses how to select a suitable value for Γ and 244 how the choice can influence the OBUIE performance as well as the PAWEC energy conversion efficiency.

246
Remark 3: The values of the scalar constants η, µ and σ affect the performance of the designed observer as well as the 248 feasibility of the LMI solution. These are chosen using trial and error, with an additional tuning requirement. However, 250 there are some simple guidelines for design tuning: i) the variable η should be as small as possible to ensure accurate 252 approximation in (24). Hence, η is set as a small positive constant less than unity. η can then be increased gradually 254 towards unity until the LMI feasibility is reached. ii) µ and σ are set with small values first, for example unity, and then 256 gradually increased to tune the design gains L 1 and L 2 . It has been found that smaller values of µ and σ cause L 1 258 and L 2 to have entries with small magnitudes, which means that the observer would be more stable to disturbances, e.g., 260 measurement noise, but it may become slow and cannot capture the high frequency content of the unknown input 262 signal. On the other hand, large values of µ and σ result in a high-gain observer which can have fast enough dynamics 264 to capture high frequency changes in the unknown input. However, a high-gain observer may be sensitive to disturbance 266 effects. In summary, the designer should iterate among values of µ and σ starting with small values seeking a balance in the 268 designed observer between sensitivity to disturbing effects and fast response. So, the methodology here is to design 270 (off-line) the values of the OBUIE gain matrices L 1 and L 2 to give acceptable robust performance over all considered sea 272 states. The estimator tuning aspect is then left to an on-line speed tuning feature via the single learning rate parameter Γ. 274

IV. SIMULATION RESULTS
A. OBUIE performance assessment and tuning 276 In this part of the simulation study the OBUIE performance evaluation and the tuning of its learning rate Γ are illustrated 278 under irregular wave conditions. Three wave states taken from the PM spectrum are considered. The considered wave states 280 architecture used in the simulation study is shown in Fig. 1.
Note: The wave states S1, S2 and S3 above are suitable for 284 the 1/50th scale PAWEC considered here. Scaling them up can give the corresponding (real) sea state using the Froude 286 scale ratio 1/50 [14]. For example, S1 corresponds to a sea state with: H s = 0.07 × 50 = 3.5m, and ω p = 4/ √ 50 = 288 0.5657 rad/s. WEF estimation accuracy is evaluated using the so-called Goodness of Fit (GoF): wheref e is the mean value of f e (t   [25]. Figure 2 shows the comparison 294 of the estimated and exact WEF values with learning rates: 2, 8 and 20. It is assumed that the exact WEF is available 296 in simulation, although the actual WEF cannot be available in a real experiment [12]. Fig. 2 shows that the estimated WEF 298 accuracy improves by increasing Γ. GoF values, calculated using (25), at the three considered values of Γ are 59.5%, 300 93.9% and 97.9%, respectively. The GoF is investigated further by varying the learning rate 302 among a wide range: Γ ∈ (1, 300), and the result is presented in Fig. 3a. The three wave sates S1,S2 and S3 are considered. 304 At Γ > 100 the GoF is very close to 100% for the 3 wave states. By zooming in for Γ ≤ 40, it appears that for Γ ≥ 20 a 306 GoF ≥ 95% is obtained. In addition, wave states with smaller significant wave heights (e.g., 0.07m here) produce slightly 308 better GoF compared with those of larger significant wave heights (e.g., 0.15 m here). Figure 3b illustrates the relation 310 between the GoF and the peak frequency of incoming wave at Γ ∈ (10, 100). Fig. 3b shows that the GoF is high at lower 312 peak frequencies and gradually declines as the peak frequency increases. This suggests that the observer "speed" should be 314 high enough to catch the fast-changing WEF frequencies, particularly at high peak frequency sea states. To prove this, 316 Γ is gradually increased from 10 to 100. It is clear that increasing Γ rises the GoF over all frequency range. Fore 318 example, Γ ≥ 50 is enough to make the GoF >90 over the whole frequency range.

320
To illustrate the effect of measurement noise on the estimator accuracy, white noise signals each with zero mean 322 and variance 1 × 10 −7 are added to each of the position and velocity outputs. In Fig. 4 the GoF is plotted against Γ 324 for four cases: (a) noise-free outputs, (b) noisy outputs, (c) filtered noisy outputs, and filtered noisy outputs with filtered 326 estimation. In the noisy output case, the GoF increases with Γ till a certain point, around Γ=50, then decreases as Γ increase. 328 At Γ = 100, for example, a 8 % drop occurs in GoF due to noise. However, filtering y(t), e.g. using a low-pass filter, 330 improves the GoF. Furthermore, filtering of both y(t) andf ul gives rise to additional GoF enhancement.

332
The above analysis indicates that higher GoF is obtained by increasing Γ when the output is noise-free. However, high 334 Γ will increase the sensitivity of the OBUIE to external disturbances, e.g., measurement noise, hence it can deteriorate 336 the OBUIE accuracy. So that, a balance between sensitivity suitable filtering of measurement data and estimated WEF is shown to improve the GoF.

B. Effects of control strategy and parametric uncertainty
This part of the simulation investigates the effect of using different control strategies on the performance of the proposed OBUIE as well as on the PAWEC performance in terms of: power conversion efficiency and amplitudes of the required PTO force. The conversion efficiency, also known as the relative capture width [2], is defined as: (26) where P W ave is the incident wave power on the buoy: J = ρg 2 32π T H 2 is known as the capture width [16], and it represents the wave power per unit width (one meter) of incident wave, T = 2π/ω is the wave period, and H is the wave height. The average mechanical power captured by the buoy is given by: where T s is the simulation time.  Loading (RL) [25] strategy where the PTO force can have only 348 a damping term and energy can only flow in the direction from the heaving buoy to the PTO. (2) The more advanced ACC 350 method [25] which is a form of Reactive Control (RC) [16] where the PTO can have both damping and stiffness effects. 352 Under reactive control, energy can flow in both directions between the heaving buoy and the PTO mechanism. Relative 354 to the RL, the ACC is known to amplify both the PTO force and the system motion (i.e., displacement and velocity). This 356 results in the amplification of the effect of the nonlinear PAWEC dynamics (e.g., the nonlinear viscous drag force). 358 In this part of the simulation harmonic (regular) waves with frequency range: 2 − 7 rad/s and wave height 0.07 m are 360 considered.
In addition, the effect of parametric uncertainty that may 362 arise in the system due to modeling errors, component ageing problems or possibly extreme wave states, is investigated. 364 Hence, some parametric uncertainties are intentionally introduced in the values of PAWEC parameters K m and M t . To 366 simulate the stochastic nature of the parametric uncertainty, the simulated parameter variations are sampled from a normal 368 distribution with mean value equal to the nominal value and with a standard deviation of 10% of nominal value. See Table I 370 for the stochastic values of perturbed PAWEC parameters. The Monte-Carlo experiment is used to get consistent outcomes 372 about the stochastic uncertainty. This is achieved by running the simulation for 50 times (1000s each) with learning rate 374 Γ = 20 considering concurrently all the randomly sampled parameter variations described above. Then, mean values for: 376 GoF, conversion efficiency, and maximum of PTO force, are computed. So, the following simulation results include two 378 cases: CASE 1) refers to the nominal model parameters, and CASE 2) refers to model with parametric uncertainty.

380
The GoF is shown in Fig. 5 for both RL and ACC controllers with and without parametric uncertainty. The figure shows 382 that the GoF is high with both RL and ACC in the nominal case (CASE 1). So that, the OBUIE almost gives similar 384 performance with both control types. On the other hand, with the uncertain model (CASE 2) parametric perturbations cause 386 (in average) a noticeable decline in GoF compared with the nominal case. The decline in GoF under uncertainty is more 388 significant near resonance, i.e., around 5 rad/s, with minimum values of about 58% for RL and about 75% with ACC. In 390 the lower frequency range, ω ≤ 4 rad/s, the GoF is above 85% with both RL and ACC. In the higher frequency range, 392 ω ≥ 6 rad/s, the GoF is ≥ 75% with ACC and ≥ 85% for RL. With uncertainty, the RL shows a relatively better 394 performance than the ACC in the off-resonance regions. Figure 6a illustrates the conversion efficiency against wave 396 frequency. It is clear that the ACC gives higher efficiency than the RL with both nominal and perturbed models. With 398 the RL. Note that the peak efficiency points with both ACC and RL shift away from the resonance point, 5 rad/s for 402 the considered PAWEC. This is because of considering the nonlinear viscous drag force in the PAWEC model, unlike the 404 linear PAWEC model which always gives the peak efficiency at resonance frequency [11]. The efficiency drops due to 406 uncertainty with both RL and ACC. Uncertainty degrades the conversion efficiency near resonance points with about 5% and 408 14% for RL and ACC, respectively. The decline in efficiency with the ACC is relatively high compared with the RL.

410
The enhanced efficiency under ACC control comes at a price of an increase in required PTO force. This is highlighted in 412 Fig. 6b which shows that ACC requires PTO force up to 60 N near its peak efficiency point while it is always less than 20 N 414 with RL. The higher PTO force may place restrictions on the size/rating of the PTO mechanism, an issue that should be 416 considered early at the design stage.

C. WEF prediction 418
As highlighted in Section I, energy maximization control necessitates WEF prediction. So, WEF estimators are usually 420 accompanied with a way to achieve WEF forecasting as described by [7], [26], [27]. It has been shown that a simple 422 Auto-Regressive model trained with past data of WEF (or water surface elevation) can be used to accurately predict 424 future values up to 2 significant wave periods [26]. As this technique is mature, it is combined with the OBUIE to get 426 WEF prediction, similar to [6], [8]. An example is shown in Fig. 7 with an irregular wave sampled from the PM spectrum 428 with: H s = 0.15m and ω p = 5 rad/s. The Simulink sampling time is set to 0.02s and the learning rate is taken as Γ = 7.

430
Two prediction horizons are considered, t hrz1 = 1.25s and t hrz2 = 2.5s. It is evident from Fig. 7 that an Auto-Regressive 432 model of order 50 can accurately predict the WEF up to one SWP in the future, with a GoF ≈ 94%. However, doubling the 434 prediction horizon, i.e., t hrz2 = 2.5s, degrades the prediction accuracy to GoF ≈ 71%.

V. DISCUSSION
A. Applicability to other PAWEC geometries 438 The proposed WEF estimator handles the (nonlinear) effect of the nonlinear viscous drag force as the considered PAWEC 440 has a cylindrical geometry for which the nonlinear viscous drag force is significant. Other nonlinear effects arise with 442 different geometries. For example, the (nonlinear) Froude-Krylov (FK) force is important in PAWECs with non-uniform 444 cross section areas [28], [29], e.g., cone or sphere. Significant work is carried out by the studies [28], [29], [30] in order 446 to provide a fast but accurate enough representation of FK force in a wide range of realistic PAWECs that can be used 448 in Model-Based control systems design. Based on the authors' current knowledge, such models with FK force cannot yet be 450 used in model-based control design.
In keeping with several control and estimation methods, the 452 OBUIE is a model-based technique that depends on a certain model structure as in (8). Further investigation of (8) reveals 454 that this structure has two parts: i) a linear part: Ax(t)+Bu(t) and ii) (unknown + nonlinear) part: B ul f ul . The second part 456 includes f e (t) plus nonlinear term(s), e.g., a nonlinear viscous drag force. For the OBUIE applicability to a certain PAWEC 458 with any degrees of freedom or geometries, It is required (a) to represent the PAWEC dynamics model in the form of (8) 460 (b) to have suitable information to calculate the nonlinear term(s), i.e., decouple the nonlinear term(s) from the (unknown 462 + nonlinear) part. Then the OBUIE can be used to estimate the (unknown + nonlinear) part as a single accumulated unknown 464 input. Consequently, f e (t) is obtained via subtraction. In this study the OBUIE is used to estimate f ul , i.e., the (unknown + 466 nonlinear) part, then the knowledge of v(t) and v w (t) is used to calculate the nonlinear viscous drag force, i.e., the nonlinear 468 term. Finally, (3) is subtracted from (10) to compute the WEF.

B. Comparison with EKF and unknown input observer esti-470 mators
The proposed OBUIE is shown to handle accurate nonlinear 472 forces, e.g. viscous or friction, in the PAWEC under the availability of suitable information to decouple them from the 474 estimation of (accumulated) unknown input. To the best of the authors' knowledge, available WEF estimators have not 476 explicitly discussed the case of nonlinearity in the PAWEC model and in most cases reported linear models are assumed. 478 Except in [9] where an EKF-based WEF estimator was designed for a nonlinear model. But it requires the computation 480 of Jacobian matrices step by step to linearise the model, adding considerable computational burden. In addition, the EKF may 482 fail if the linearized model is far from the actual nonlinear process. This is actually a well-known potential limitation of 484 the EKF which requires the determination of a domain of convergence that is difficult to determine [31], especially for 486 a system as complex as wave motion with varying properties. This is in agreement with [6] where the authors reported that 488 EKF is limited to regular wave and it may diverge when applied to irregular wave scenarios. However, it is possible that 490 some available WEF estimators such as [10] and [12] could be extended to handle nonlinear terms in a similar fashion as 492 presented in this paper.
The OBUIE has some similarity with the unknown input 494 observer based estimators in [12] and [13], as both do not need representation of WEF dynamics in the estimation model. 496 However, the OBUIE has more flexibility than [12] and [13] as the user can tune a learning rate (Γ) to increase the speed of the observer and hence capture possibly fast-changing WEF dynamics, leading to better estimation accuracy. In contrast, 500 the speed of the WEF estimators proposed in [12], [13] is determined during the (off-line) design stage using an LMI-502 based strategy for pole-assignment. The downside is that the LMI design may fail to find a feasible solution.

VI. CONCLUSION AND FUTURE WORK
The paper presents the OBUIE technique to estimate the 506 WEF on a single DOF PAWEC. Unlike some WEF estimators available in the literature, the proposed OBUIE has a simple 508 structure as it does not need a representation of the WEF dynamics in the estimator. It is also able to handle the non-510 linear effect due to viscous force but requires the knowledge of the water surface velocity. Finally, the OBUIE does not 512 require heavy computation. The Matlab simulation results show the OBUIE validity to estimate the WEF corresponding 514 to basic (RL) and advanced (ACC) controllers, and subject to both regular and irregular wave conditions. Under parametric 516 perturbations, the OBUIE continues to deliver acceptable WEF estimation but with reduced PAWEC conversion efficiency.

518
This comes as a consequence of parametric uncertainty. The OBUIE is based on the cylindrical-shaped PAWEC dynamics.

520
Further analysis is ongoing to study the applicability of OBUIE on PAWECs with complex geometries and multi-DOF, 522 where more nonlinear effects may arise and be signification.