METHOD FOR CONTROLLING A WAVE-ENERGY SYSTEM BY DETERMINING THE EXCITATION FORCE APPLIED BY WAVES INCIDENT UPON A MOVING PART OF THE SAID SYSTEM

Abstract

The invention relates to the real-time determination of the forces applied by waves incident upon the moving part (2) of a wave-energy system (1). The method according to the invention is based on the construction of a model of the radiation force applied to the moving part (2) and a model of the dynamics of the wave-energy system (1). The invention uses only measurements of the kinematics of the moving part (2) and the force applied by the conversion machine (3) to the moving part (2).

Description

The invention relates to the sector of devices for converting the energy of waves into electrical energy. In particular, the invention relates to the sector dealing with determination of the excitation force of waves incident upon a moving part, in particular for controlling a wave-energy system.

For some years there has been a growing interest in renewable energy resources. These resources are clean, free and limitless, these being major advantages in a world where there is a relentless decrease in the fossil fuel reserves available and a growing awareness of the need to preserve the planet. Among these resources, wave power energy, which is a relatively unknown energy source compared to other widely advertised energy sources such as wind power and solar energy, helps ensure an indispensable diversification of the renewable energy sector. The devices commonly known as “wave-energy converters” are of particular interest because they allow electricity to be produced from this renewable energy source (potential and kinetic energy of waves) without the emission of greenhouse gases. They are particularly suitable for providing electricity in isolated island areas.

For example, the patent applications FR 2876751 and WO 2009/081042 describe apparatus for capturing the energy produced by the sea tide. These devices are composed of a floating support on which a pendulum is mounted so as to be movable in relation to the floating support. The relative movement of the pendulum in relation to the floating support is used to produce the electrical energy by means of an energy conversion machine (for example an electrical machine). The conversion machine operates both as a generator and as a motor. In fact, in order to provide a torque or a force which drives the moving element, power is supplied to the conversion machine in order to place the moving element in synchronism with the waves (motor mode). On the other hand, in order to produce a torque or force which resists the movement of the moving element, power is recovered via the conversion machine (generator mode).

In order to improve the performance and therefore the profitability of the devices for converting the energy of waves into electrical energy (wave-energy systems), it is of interest to estimate the excitation force applied by the wave incident upon the wave-energy system.

The force applied by waves to the devices for converting the energy of the waves into electrical energy or other forms of exploitable energy, commonly known as wave-energy converters, cannot be directly measured under normal operating conditions in an industrial context. In the case of oscillating point converters, it is possible to carry out specific tests consisting in keeping still the floating part of the machine and measuring the force needed to oppose the action of the wave (namely to keep the float immobile) by means of one or more (force or torque) sensors positioned on the conversion machine, also referred to as a “power take-off” (PTO) system, for example an electric generator combined with a device allowing transmission of the oscillating movement of the moving part. In present-day terminology the quantity thus measured is called the excitation force of the incident swell (or wave), this being distinguished in particular from the force (and the swell) generated by the movement itself of the float (radiation force). During normal operation of the float, on the other hand, the same force sensors will measure only the force applied by the PTO and not the excitation force of the incident wave. It would be possible, in principle, to calculate all the forces applied by the wave by incorporating pressure measurements obtained from sensors distributed over the whole surface, but this represents a costly solution and unreliable solution which may be difficult to implement in an industrial context.

Few scientific studies exist in the field of real-time estimation of the excitation force (or torque) of waves. Among those which exist there is, for example, the document by Peter Kracht, Sebastian Perez-Becker, Jean-Baptiste Richard and Boris Fischer: “Performance Improvement of a Point Absorber Wave-Energy system by Application of an Observer-Based Control: Results from Wave Tank Testing”. In IEEE Transactions on Industry Applications, July 2015, which describes an algorithm for estimating the force of the wave, the latter is based on a bank of independent harmonic oscillators and a Luenberger observer. The results obtained with this method show that the estimations are significantly delayed (out-of-phase) in relation to the real excitation force, this posing problems in particular for use in connection with control of the wave-energy system. More generally, this approach considers a fixed frequency range (chosen a priori) for the wave spectrum. In order for the method to be able to work in realistic conditions, where the spectrum varies over time, a very large number of frequencies must be taken into account, and this makes the approach very complex in terms of calculations. In the document by Bradley A Ling: “Real-Time Estimation and Prediction of Wave Excitation Forces for Wave Energy Control Applications”, published in ASME 2015 34th International Conference on Ocean, Offshore and Arctic Engineering, an approach based on an extensive Kalman filter is proposed in order to reconstruct the spectrum where the excitation force of the wave is regarded as being a (single) sinusoid, the parameters (amplitude, frequency and phase) vary over time. However, the method may be effectively used only for waves in a very narrow frequency range.

Moreover, the patent application FR2973448 (U.S. Pat. No. 9,261,070) describes a control method for oscillating point converters, one of the steps of which involves the estimation (of the spectrum) of the force of the wave on the float (or of the torque on the moving part of the converter), using a set of sinusoids and a Luenberger observer. The method resembles that published in the aforementioned documents and has a priori the same drawbacks.

The invention aims to improve the operation of a wave-energy system by determining in real time the forces applied by waves incident upon the moving part so as to adjust consequently in an optimum manner the electrical energy recovery operations. The method according to the invention is based on the construction of a model of the radiation force applied to the moving part and a model of the dynamics of the wave-energy system. The invention uses only measurements of the kinematics of the float (position, speed and if necessary acceleration) and the force applied by the conversion machine, said measurements being normally available on a wave-energy system, because they are used for control and management thereof. Thus, owing to the models used, it is possible to estimate the force applied by the incident wave for the entire range of frequencies of the waves, while maintaining a calculation time suitable for real-time determination and control.

THE METHOD ACCORDING TO THE INVENTION

The invention relates to a method for controlling a wave-energy system by determining the excitation force applied by waves incident upon a moving part of a wave-energy system, the said wave-energy system converting the energy of the waves into energy by means of the said moving part cooperating with a conversion machine, the said moving part performing a movement in relation to the said conversion machine as a result of the action of the waves. For this method, the following steps are performed:

• • a) measurement of the position and speed of the said moving part;
  • b) measurement of the force F_u applied by the said conversion machine to the said moving part;
  • c) construction of a model of the radiation force F_rad applied to the said moving part, the said model of the radiation force F_rad relating the said radiation force to the speed of the said moving part;
  • d) construction of a model of the dynamics of the said wave-energy system which relates the said excitation force F_ex applied by waves incident upon the said moving part, to the said position of the said moving part, to the said speed of the moving part, to the said force F_ex applied by the said conversion machine to the said moving part, and to the said radiation force F_rad applied to the said moving part; and
  • e) determination of the said excitation force F_ex applied by waves incident upon the said moving part by means of the said dynamics model, the said model of the radiation force, the said position and speed measured and the said measured force F_u applied by the said conversion machine to the said moving part; and
  • f) controlling the said wave-energy system depending on the said determined excitation force F_ex applied by waves incident upon the said moving part.

According to one embodiment, the said dynamics model of the said wind-energy system is constructed by means of an equation of the type: M{umlaut over (z)}(t)=F_ex(t)+F_hd(t)+F_rad(t)−F_u(t) where M is the mass of the moving part, {umlaut over (z)} is the acceleration of the said moving part, and F_hd is the hydrostatic restoring force applied to the said moving part.

According to an alternative, the said hydrostatic restoring force F_hd is determined by means of a formula of the type:

F_hd(t)=−Kz(t) where z is the position of the said moving part and K is the hydrostatic stiffness coefficient.

Alternatively, the said hydrostatic restoring force F_hd is determined by means of a piecewise affine function of the type: F_hd=−F_iz(t)−G_i for at least two variation intervals al of the said position z, where F_i and G_i are constant coefficients of the said affine function for each interval Ω_i.

According to a mode of implementation the said model of the radiation force F_rad is constructed by an equation of the type: F_rad(t)=−M_∞{umlaut over (z)}(t)−F_r(t) where F_r(t)=∫₀^th(t−τ)ż(τ)dτ is the radiation damping, {umlaut over (z)} is the acceleration of the said moving part, M_∞ is the mass added to the infinitely high frequency, ż is the speed of the said moving part, h is the pulse response which relates the speed of the moving part to the radiation damping.

Advantageously, the acceleration of the said moving part is measured.

According to a first embodiment, the radiation damping F_r is determined by means of an appointed time response resulting from a state realization, of the type: F_r(t)=C_r∫₀^te^(t-τ)ArB_rż(τ)dτ+D_rż(t) where A_r, B_r, C_r, D_r are matrices of the state realization.

Preferably, the said excitation force F_ex applied by waves incident upon the said moving part (2) is determined by means of an equation of the type: {circumflex over (F)}_ex(k)=(M+M_∞){umlaut over (z)}(k)+Σ_l=0^k-1C_rdA_rd^k-1-lB_rdż(l)+D_rdż(k)+F_iz(k)+G_i+F_u(k) where k is the time step in the discrete domain, and A_rd, B_rd, C_rd, D_rd are the matrices of the realization state in the discrete domain.

According to a second embodiment, the said excitation force F_ex applied by waves incident upon the said moving part is determined by means of a state observer based on a linear Kalman filter constructed using a random-walk model of the said excitation model F_ex applied by waves incident upon the said moving part.

Advantageously, the said excitation force F_ex applied by waves incident upon the said moving part is determined by performing the following steps:

• • i) initialization of k=0, the state vector {circumflex over (x)}_a(0|0)=x_a(0) and the state of the covariance matrix, P(0|0)=P₀;
  • ii) at each instant k, acquisition of the said position and speed measurements of the said moving part y(k)=[z(k) ż(k)] and measurement of the said force applied by the said conversion machine to the said moving part u(k)=F_u(k); and
  • iii) at each instant k, determination of the excitation force applied by waves incident upon the said moving part (2) ŵ(k) by means of the following equations: ŵ(k)=[0 1] {circumflex over (x)}_a(k|k) where:

$\{\begin{array}{c}{\hat{x}}_a\left(k|k-1\right)=A_a{\hat{x}}_a\left(k-1|k-1\right)+B_aF_u\left(k-1\right)\\ P\left(k|k-1\right)=A_aP\left(k-1|k-1\right)A_a^T+Q\end{array},K\left(k\right)=P\left(k|k-1\right){C_a^T\left(C_aP\left(k|k-1\right)C_a^T+R\right)}^{-1}{\hat{x}}_a\left(k|k\right)={\hat{x}}_a\left(k|k-1\right)+K\left(k\right)\left(y\left(k\right)-C_a{\hat{x}}_a\left(k|k-1\right)-D_aF_u\left(k\right)\right)P\left(k|k\right)=\left(I-K\left(k\right)C_a\left(k\right)\right)P\left(k|k-1\right)$

• • • where K(k) is the gain of the observer, x_a(k) is the unknown state vector

$x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],$

A_a, B_a, C_a, D_a are state realization matrices, P is the covariance matrix of the state vector, and Q and R are calibration matrices.

According to a third embodiment, the said excitation force F_ex applied by waves incident upon the said moving part is determined by means of an unknown input state estimator based on a moving-horizon approach.

According to a characteristic feature, the said excitation force F_ex applied by waves incident upon the said moving part is determined performing the following steps at each instant k:

• • i) acquisition of the said position and speed measurements of the said moving part y(k)=[z(k) ż(k)] and measurement of the said force applied by the said conversion machine to the said moving part u(k)=F_u(k); and
  • ii) determination of the excitation force applied by waves incident upon the said moving part ŵ(k|k) by means of the following equations: ŵ(k|k)=ξ(N) where ξ(N) is obtained by means of resolution of a quadratic programming problem:

$\hspace{1em}\{\begin{array}{c}\underset{\xi }{\min }\left\{{\xi }^T\left(P+{\Psi }_s^TR\Psi \right)\xi -2\left[{\hat{w}}_0^T{\hat{x}}_0^TV^T\right]\left[\begin{array}{c}S\\ R{\Psi }_s\end{array}\right]\xi \right\}\\ \mathrm{such}\mathrm{as}H_{\xi }\xi \le g_{\xi }+H_dU_d\end{array}$

where P, R and S are matrices formed from cost function weighting matrices, ξ is the matrix which represents all the values of the unknown values in the estimation window of size N, Ψ_s is the matrix, calculated using the state-space representation, which corresponds to the effect of these variables on the values of the outputs in this same window, H_ξ and H_d are matrices and g_ξ a vector obtained from constraints imposed on the state and the unknown input, {circumflex over (x)}₀={circumflex over (x)}(k−N|k−1) and ŵ₀=ŵ(k−N|k−1) where x is the state vector and w is the unknown input of the state estimator,

$U_d=\left[\begin{array}{c}{u}_d\left(k-N\right)\\ u_d\left(k-N+1\right)\\ ⋮\\ u_d\left(k\right)\end{array}\right],$

u_d being the known input of the state estimator where u_d(t)=F_hd(t)−F_u(t).

According to a mode of implementation, operation of the said wave-energy system is controlled depending on the said determined excitation force F_ex applied by waves incident upon the said moving part.

BRIEF DESCRIPTION OF THE FIGURES

Further characteristic features and advantages of the method according to the invention will become clear from a reading of the description below of non-limiting examples of embodiment, with reference to the attached figures described below.

FIG. 1 shows wave-energy system according to an embodiment of the invention.

FIG. 2 shows an approximation of the hydrostatic restoring force as a function of the position according to an example.

FIG. 3 shows the amplitude of waves for implementation of the illustrative example.

FIGS. 4a to 4c show the excitation force applied by waves incident upon the moving part as a function of time, for a wave of the type W1 according to FIG. 3 obtained by measurement and respectively by one of the embodiments of the method according to the invention.

FIGS. 5a to 5c show the excitation force applied by waves incident upon the moving part as a function of time, for a wave of the type W3 according to FIG. 3, obtained by measurement respectively by one of the embodiments of the method according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

The present invention relates to the determination of the excitation force applied by waves incident upon a moving part of a wave-energy system, called also excitation force of the swell or wave. A wave-energy system is a system which converts the energy of the wave into recoverable energy, in particular electrical energy. A wave-energy system generally comprises a moving part, which is also called pendulum or float, which has an oscillating movement when acted on by the wave. The moving part cooperates with a conversion machine, called also “power take-off” (PTO) system, which in most cases comprises an electric generator connected to a device allowing the transmission of the oscillating movement to be adapted, so as to convert the movement of the moving part into recoverable energy. In certain cases, the conversion machine may act as a motor generating a force on the moving part. In fact, in order to recover power via a conversion machine, a torque or force which resists the movement of the moving element is produced (generator mode). On the other hand, if the conversion machine allows it, power may be supplied to the conversion machine so as to provide a torque or a force which drives the moving element so as to help it be placed in synchronism with the waves (motor mode). The method according to the invention relates to the control of a wave-energy system depending on the excitation force determined,

The method according to the invention is suitable for all types of wave-energy system with at least one moving part, for example those described in the patent application FR2973448 (U.S. Pat. No. 9,261,070). The control method according to the invention may also be applied to a wave-energy system belonging to the category of wave-energy systems with oscillating water columns (OWC).

FIG. 1 shows, in schematic form, a non-limiting example of a wave-energy system 1. The wave-energy system 1 comprises a floating moving part 2 (the wave is represented in schematic form by the hatched area). The moving part 2 cooperates, by means of a lever arm 4, together with a conversion machined based on the electric generator 3 which, in this case, is a linear generator. The lever arm 4 and the electric generator are hingeably mounted relative to an installation 6. The movements of the lever arm are indicated in schematic form by means of arrows.

In the continuation of the description and the claims, the terms “waves”, “sea tides” and “swell” are considered to be equivalent.

Moreover, in the description and the claims, the term “force” indicates a force or a torque. Similarly, the terms “position”, “speed” and “accelerations” all indicate values which are “linear” or “angular”, whereby the “linear” values may be associated with forces, and the “angular” values may be associated with “torques”.

Furthermore, for better comprehension, the different models are shown in a single dimension; however, the method according to the invention is suitable for models with various dimensions, designed for systems which have several degrees of freedom in their movement.

The method according to the invention comprises the following steps:

1) measurement of the position and speed of the moving part;

2) measurement of the force applied by the conversion machine (or PTO)

3) construction of the radiation force model

4) construction of the dynamics model

5) determination of the excitation force of the incident wave

6) controlling of the wave-energy system

1) Measurement of the Position and Speed of the Moving Part:

During this step, the position and the speed of the moving part are measured. The position corresponds to the movement (for example distance or angle) in relation to the equilibrium position of the moving part. These measurements may be carried out by means of sensors, which are generally provided on a wave-energy system for control and/or management thereof.

According to a mode of implementation of the invention, during this step, it is also possible to measure or estimate the acceleration of the moving part, since this measurement or estimation may be used in the models implemented by the method according to the invention. For example, the acceleration may be measured by means of an accelerometer mounted on the moving part.

2) Measurement of the Force Applied by the Conversion Machine (PTO)

During this step, the force (or if applicable the torque) applied by the conversion machine (PTO) to the moving part is measured. This measurement may be performed by means of a sensor, which may be a force sensor or a torque sensor. This type of sensor is often installed, or may be easily installed in wave-energy systems, for control and/or management thereof. Alternatively, the measurement may be replaced by an estimation carried out based on a set force (or torque) value sent to the PTO.

By way of example of the wave-energy system shown in FIG. 1, a torque sensor 5 may be arranged between the lever arm 4 and the electric generator 3.

3) Construction of a Radiation Force Model

During this step, a model of the radiation force applied to the moving part is constructed. According to the linear wave theory (as described for example in the document by Falnes J, Kumiawan A. “Fundamental formulae for wave-energy conversion”. R. Soc. open sci. 2: 140305, 2005, http: //dx. doi. org/10. 1098/rsos. 140305), the radiation force is produced by the oscillation of an immersed body (and therefore depends on the movement of the moving part), while the excitation force, resulting from the presence itself of a body in the water, does not depend on the movement of the immersed body, but on the incident wave. If there is no incident wave, the radiation force dampens the residual oscillation of the immersed body until it causes it to stop. It is important to note that, although the linear theory allows the excitation force to be related to the height of the incident wave by means of a linear model (in the frequency or time domain), in practice it may not be used to calculate the excitation force linearly, even if it was possible to measure the height of the wave at the center of gravity of the float as required by the theory. In fact, the linear relation between height of the wave and excitation force is not random, this meaning that the excitation force cannot be calculated at a given instant without knowing the height of the wave in future instants (the calculation may on the other hand be performed non-linearly once the wave has passed). In the context of real-time control, the excitation force may therefore not be regarded as being a totally unknown exogenous force acting on the float. On the other hand, still according to the linear wave theory, the radiation force is related to the movement of the float and more precisely to its acceleration and its speed by a random model which is linear (in the frequency or time domain). It may therefore be calculated linearly using the current acceleration and speed measurements (and past speed measurements).

According to a mode of implementation of the invention, the said model of the radiation force F_rad is constructed by an equation of the type:

F _rad(t)=−M_∞{umlaut over (z)}(t)−F_r(t)

where F_r(t)=∫₀^th(t−τ)ż(τ)dτ is the component of the said radiation force F_rad which depends on the (current and past) speed of the moving part, which may be called radiation damping,

{umlaut over (z)} is the acceleration of the moving part,

M_∞ is the mass added to the infinitely high frequency, which may be obtained by means of codes of BEM (Boundary Element Method) calculation codes, such as WAMIT® (WAMIT, USA), or Nemoh® (Ecole Centrale de Nantes, France), based on the geometry of the moving part,

ż is the speed of the moving part,

h is the pulse response which relates the speed of the moving part to the radiation damping, the coefficients of which are obtained from the hydrodynamic parameters of the moving part calculated by the same BEM calculation codes.

The construction of this model allows the determination at any instant of the radiation force, with a limited calculation time. Thus the determination of the force applied by the wave may be determined at any instant with a short calculation time.

4) Construction of the Dynamics Model of the Wave-Energy System

During this step, a dynamics model of the wave-energy system is constructed. “Dynamics model” refers to a model which relates the excitation force applied by waves incident upon the moving part, the radiation force applied to the moving part, the hydrostatic restoring force applied to the moving part and the force applied by the conversion on the moving part, to the position and speed of the moving part. With this type of model it is possible obtain results representing the behavior of the wave-energy system, if the movements are not too big.

Advantageously, the dynamics model is obtained by applying the fundamental principle of dynamics (Newton's second law) to the moving part.

According to a mode of implementation of the invention, where the forces are considered, the dynamics model of the wave-energy system may be constructed by an equation of the type:

M{umlaut over (z)}(t)=F_ex(t)+F_hd(t)+F_rad(t)−F_u(t)

where M is the mass of the moving part,

{umlaut over (z)} is the acceleration of the moving part,

F_ex is the excitation force applied by waves incident upon the moving part,

F_rad is the radiation force applied to the moving part,

F_hd is the hydrostatic restoring force applied to the moving part, and

F_u is the force applied by the conversion machine to the moving part.

This model translates a vertical translation movement (typically of floats which have a pitching movement). This model is derived from the linear wave theory.

In accordance with a first variation of embodiment, the hydrostatic restoring force applied to the moving part may be approximated using a linear function of the position z defined in relation to the equilibrium position. In this case, the hydrostatic restoring force may be written by a function of the type: F_hd(t)=−Kz(t) where z is the position of the said moving part defined in relation to its equilibrium position and K is the hydrostatic stiffness coefficient. Thus the hydrostatic restoring force may be calculated using a simple model if the measurement of the position z is available. This function is particularly suitable for small displacements z.

According to a second variation of embodiment, the hydraulic restoring force applied to the moving part may be approximated using a piecewise affine function of the position z. The function may be linear for at least two—preferably three—intervals of the position z. In this case, the hydrostatic restoring force may be written by a function of the type: F_hd=−F_iZ(t)−G_i for at least two variation intervals Ω_i of the position z, where F_i and G_i are constant coefficients of the affine function for each interval Ω_i. This function is particularly suitable for the larger movements z. FIG. 2 illustrates this approximation. FIG. 2 shows in continuous lines, in a schematic and non-limiting manner, the variations of the hydrostatic restoring force F_hd as a function of the position z. This figure shows in broken lines the linear functions of each piece Ω₁, Ω₂ and Ω₃. It can be noted that this piecewise affine function allows a better approximation of the hydrostatic restoring force over the whole of the variation range of the position z.

According to a mode of implementation of the invention, where the torques are considered, the dynamics model of the wave-energy system may be constructed by an equation of the type:

J{umlaut over (θ)}(t)=M_ex(t)+M_hyd(t)+M_rad(t)−M_u(t)

where J is the inertia of the moving part,

{umlaut over (θ)} is the angular acceleration of the moving part,

M_ex is the excitation torque applied by waves incident upon the moving part,

M_rad is the radiation torque applied to the moving part,

M_hd is the hydrostatic restoring torque applied to the moving part, and

M_u is the torque applied by the conversion machine to the moving part.

This model converts a rotational movement about a horizontal axis (typical of floats which have a pitching movement). This model is derived from the linear wave theory. This model is a “mirror” of the model where the forces are considered: the various terms of the model are of the same nature.

In the same way as for the hydrostatic restoring force, the hydrostatic restoring torque may be approximated by a linear function or by a piecewise affine function.

In the description below, only the example relating to forces will be described, but the example relating to torques may be deduced by means of transposition of the equations in an angular reference system.

5) Determination of the Excitation Force Applied by Incident Waves

During this step, the excitation force applied by waves incident upon the moving part is determined by means of:

• • measurements of the position and speed (and optionally acceleration) of the moving part,
  • measurement of the force applied by the conversion machine (PTO) to the moving part,
  • the model of the radiation force, and
  • the dynamics model of the wave-energy system.

According to a first embodiment of the invention (described in detail below), it is possible to determine the excitation force applied by waves incident upon the moving part using the calculation of an approximation of the radiation force. This embodiment has the advantage that a representation of the hydrostatic restoring force by a piecewise affine function may be taken into account.

According to a second embodiment of the invention (described in detail below), the excitation force applied by waves incident upon the moving part may be determined by means of a state observer based on a linear Kalman filter constructed using a random-walk model of the excitation force applied by waves incident upon the moving part. This second embodiment has the advantage that uncertainty factors are taken into account into account.

According to a third embodiment of the invention (described in detail below), the excitation force applied by waves incident upon the moving part may be determined by means of an unknown input state estimator based on a moving-horizon approach. This embodiment has the advantage that uncertainty factors are taken into account. Moreover a representation of the hydrostatic restoring force by a piecewise affine function may be taken into account.

Other embodiments may be envisaged.

First Embodiment

According to the first embodiment of the invention, it is possible to determine the excitation force applied by waves incident upon the moving part using the calculation of an approximation of the radiation force.

For this embodiment, the acceleration of the moving part is also measured. Prior to this step, the following are therefore known:

• • the measurements of the position z, speed ż and acceleration {umlaut over (z)} of the moving part,
  • the measurement of the force F_u applied by the conversion machine (PTO) to the moving part,
  • the model of the radiation force F_rad(t)=−M_∞{umlaut over (z)}(t)−F_r(t) where F_r(t)=∫₀^th(t−τ)ż(τ)dτ, and
  • the dynamics model of the wave-energy system M{umlaut over (z)}(t)=F_ex(t)+F_hd(t)+F_rad(t)−F_u(t).

Thus, knowing how to calculate F_r it is therefore possible to obtain F_rad, and finally F_ex(t), the only unknown factor of the dynamics model of the wave-energy system.

However, it is not easy to perform the calculation of the convolution product within the time domain in the preceding expression for F_r(t), owing to the large number of mathematical operations which must be performed and data to be stored. In order to avoid this, it is possible to regard it as a linear system with F_r(t) at the output and ż(τ) at the input. The resultant equation in the (Laplace) frequency domain, obtained by applying Prony's method, is:

F _r(s)=W_r(s)ż(s)

where W_r(s) is a transfer function. The preceding equation, which was in the Laplace domain, may be expressed in an equivalent state form, for example:

$\hspace{1em}\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_rx_r\left(t\right)+B_r\stackrel{.}{z}\left(t\right)\\ F_r\left(t\right)=C_rx_r\left(t\right)+D_rz\left(t\right)\end{array}$

where x_r is an internal state which does not have any particular physical significance and (A_r, B_r, C_r, D_r) are state realization matrices. It is pointed that in system theory (and in automatic control engineering), a state-space representation allows a dynamic system to be modelled in matrix form using state variables. This representation may be linear or non-linear, continuous or discrete. The representation allows determination of the internal state and the outputs of the system at any future instant knowing the state at the initial instant and the behavior of the input variables which influence the system.

The state form allows the evolution in time of F_r to be expressed as follows:

F _r(t)=C_re^tArx_r(0)+C_r∫₀^te^(t-τ)ArB_rż(τ)dτ+D_rż(t)

but this expression may not be used directly because the initial value of the internal state x_r(0) is unknown. However, it is possible to approximate the time evolution F_r as follows:

F _r(t)=C_r∫₀^te^(t-τ)ArB_rż(τ)dτ+D_rż(t)

since C_re^tArx_r(0)≈0 where t is sufficiently wide, given that the matrix A_r is asymptotically stable for this type of system. This corresponds to ignoring the contribution of the initial state to the time evolution of F_r, knowing that it will recede very rapidly compared to the contribution linked to the speed ż.

The first embodiment therefore consists in calculating F_ex from z, ż, {umlaut over (z)} and F_u (which are known) using an equation with the form:

{circumflex over (F)} _ex(t)=(M+M_∞){umlaut over (z)}(t)+C_r∫₀^te^(t-τ)ArB_rż(τ)dτ+D_rż(t)+Kz(t)+F_u(t)

if a linear model is used for the hydrostatic restoring force, or

{circumflex over (F)} _ex(t)=(M+M_∞){umlaut over (z)}(t)+C_r∫₀^te^(t-τ)ArB_rż(τ)dτ+D_rz(t)+F_iz(t)+G_i+F_u(t)

for each interval Ω_i to which z belongs, if a piecewise affine model is used for the hydrostatic restoring force.

Since in practice the measurements of z, ż, {umlaut over (z)} and F_u are sampled signals, the algorithm may be developed within the discrete time k, discretizing as follows the state equation introduced beforehand in order to represent the evolution of F_r:

$\hspace{1em}\{\begin{array}{c}{x}_r\left(k+1\right)=A_{\mathrm{rd}}x_r\left(k\right)+B_{\mathrm{rd}}\stackrel{.}{z}\left(k\right)\\ F_r\left(k\right)=C_{rd}x_r\left(k\right)+D_{rd}\stackrel{.}{z}\left(k\right)\end{array}$

F_r may thus be expressed by an equation with the form:

$F_r\left(k\right)=C_{rd}A_{rd}^kx_r\left(0\right)+\sum _{l=0}^{k-1}C_{rd}A_{rd}^{k-1-l}B_{rd}\stackrel{.}{z}\left(l\right)+D_{rd}\stackrel{.}{z}\left(k\right)$

and its approximation expressed by the formula:

$F_r\left(k\right)=\sum _{l=0}^{k-1}C_{rd}A_{rd}^{k-1-l}B_{rd}\stackrel{.}{z}\left(l\right)+D_{rd}\stackrel{.}{z}\left(k\right)$

This gives (in the case where the hydrostatic restoring force is considered to be a piecewise affine function) the following expression:

{circumflex over (F)} _ex(k)=(M+M_∞){umlaut over (z)}(k)+Σ_l=0^k-1C_rdA_rd^k-1-lB_rdż(l)+D_rdż(k)+F_iz(k)+G_i+F_u(k)

for each interval Ω_i to which z(k) belongs.

The original aspect of this embodiment is the approximation of the time response of F_r, the radiation damping (the component of the radiation force dependent on the speed of the float), calculated by means of an equivalent state-space representation and ignoring the contribution of the initial state of this representation.

Second Embodiment

According to a second embodiment of the invention, the excitation force applied by waves incident upon the moving part may be determined by means of a state observer based on a linear Kalman filter constructed using a random-walk model of the excitation force applied by waves incident upon the moving part. It is pointed out that a state observer, or state estimator, is in automatic control engineering and system theory an extension of a model represented in the form of a state-space representation. When the system state cannot be measured, an observer allowing the state to be reconstructed from a model is constructed.

For this embodiment, prior to this step, the following are therefore known:

• • the measurements of the position z, and speed Z of the moving part,
  • the measurement of the force F_u applied by the conversion machine (PTO) to the moving part,
  • the model of the radiation force F_rad(t)=−M_∞{umlaut over (z)}(t)−F_r(t) where F_r(t)=∫₀^th(t−τ)ż(τ)dτ, and
  • the dynamics model of the wave-energy system M{umlaut over (z)}(t)=F_ex(t)+F_hd (t)+F_rad(t)−F_u(t).

In this approach, the problem of estimation of the excitation force of the waves is transformed into a classic state estimation problem (which may be resolved with a linear Kalman filter), where the dynamics of the excitation force of the waves are expressed by a random-walk model. The main advantage of this method is the taking into consideration of uncertainty factors which allows measurement noise and non-modelled dynamics to be taken into account.

By replacing in the equation which describes the movement of the float

M{umlaut over (z)}(t)=F_ex(t)+F_hd(t)+F_rad(t)−F_u(t)

the expressions for the hydrostatic restoring force (with the linear model) and the model of the radiation force

F _hd(t)=−Kz(t)

F _rad(t)=−M_∞{umlaut over (z)}(t)−F_r(t)

the following is obtained

(M+M_∞){umlaut over (z)}(t)+Kz(t)+F_r(t)=F_ex(t)−F_u(t)

This equation may be expressed in state form, by defining

$\hspace{1em}\{\begin{array}{c}{x}_1\left(t\right)=z\left(t\right)\\ x_2\left(t\right)=\stackrel{.}{z}\left(t\right)\end{array}$

which gives

$\hspace{1em}\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=-\frac{K}{M+M_{\infty }}x_1\left(t\right)-\frac{1}{M+M_{\infty }}F_r\left(t\right)+\frac{1}{M+M_{\infty }}\left(F_{\mathrm{ex}}\left(t\right)-F_u\left(t\right)\right)\end{array}$

In the same way as for the first embodiment, the dynamics of radiation damping F_r(t) may be described by a state-space representation equivalent to the convolution product of the pulse response of the radiation damping and the speed of the moving part.

F _r(t)=∫₀^th(t−τ)ż(τ)dτ

namely

$\hspace{1em}\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_rx_r\left(t\right)+B_r\stackrel{.}{z}\left(t\right)\\ F_r\left(t\right)=C_rx_r\left(t\right)+D_r\stackrel{.}{z}\left(t\right)\end{array}\mathrm{or}\hspace{1em}\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_rx_r\left(t\right)+B_rx_2\left(t\right)\\ F_r\left(t\right)=C_rx_r\left(t\right)+D_rx_2\left(t\right)\end{array}$

where x_r is an internal (non-measurable) state which does not have any particular physical significance, and (A_r, B_r, C_r, D_r) are matrices of the state realization.

Combining the two state-space representations, the following are obtained:

$\hspace{1em}\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=\begin{array}{c}-\frac{K}{M+M_{\infty }}x_1\left(t\right)-\frac{D_r}{M+M_{\infty }}x_2\left(t\right)-\frac{C_r}{M+M_{\infty }}x_r\left(t\right)+\\ \frac{1}{M+M_{\infty }}\left(F_{\mathrm{ex}}\left(t\right)-F_u\left(t\right)\right)\end{array}\\ {\stackrel{.}{x}}_r\left(t\right)=B_rx_2\left(t\right)+A_rx_r\left(t\right)\end{array}$

or, in an equivalent manner

$\{\begin{array}{c}\stackrel{.}{x}\left(t\right)=A_cx\left(t\right)+B_c\left(F_{\mathrm{ex}}\left(t\right)-F_u\left(t\right)\right)\\ y\left(t\right)=C_cx\left(t\right)\end{array}\mathrm{where}x\left(t\right)={\left[x_1\left(t\right)x_2\left(t\right)x_r^T\left(t\right)\right]}^Ty=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]\mathrm{and}A_c=\left[\begin{array}{ccc}0& 1& 0\\ -\frac{K}{M+M_{\infty }}& -\frac{D_r}{M+M_{\infty }}& -\frac{C_r}{M+M_{\infty }}\\ 0& B_r& A_r\end{array}\right],B_c=\left[\begin{array}{c}0\\ \frac{1}{M+M_{\infty }}\\ 0\end{array}\right],C_c=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

It is necessary to estimate F_ex(t) using y(t) (the vector which contains the position and the speed measurements of the moving part) and F_u(t), the measurement of the force of the PTO.

In order to do this, firstly the state system above is discretized (it should be remembered that the measurements are sampled), using Tustin's method. This gives

$\hspace{1em}\{\begin{array}{c}x\left(k+1\right)=A_dx\left(k\right)+B_d\left(F_{\mathrm{ex}}\left(k\right)-F_u\left(k\right)\right)\\ y\left(k\right)=C_dx\left(k\right)+B_d\left(F_{\mathrm{ex}}\left(k\right)-F_u\left(k\right)\right)\end{array}$

It should be noted that it is important to perform discretization using Tustin's method because it allows a direct term to appear in the output equation which allows y(k) to be related to F_ex(k), and therefore F_ex(k) to be estimated using y(k). This would not be possible with other discretization methods, which would relate y(k) to F_ex(k−1) and not to F_ex(k).

By defining

w(k)=F_ex(k),u(k)=F_u(k),

the global model of the float dynamics is obtained:

$\hspace{1em}\{\begin{array}{c}x\left(k+1\right)=A_dx\left(k\right)+B_dw\left(k\right)-B_du\left(k\right)\\ y\left(k\right)=C_dx\left(k\right)+D_dw\left(k\right)-D_du\left(k\right)\end{array}$

In order to estimate the excitation force w(k), it is regarded as a state, the evolution of is defined by introducing a mathematical model which relates w(k) to w(k−1).

For this embodiment, a random-walk model is used to describe the evolution of the excitation force:

w(k+1)=w(k)+ϵ_m(k)

where ϵ_m(k), which describes the variation of w(k), is regarded as a random number. In other words this model assumes that at each instant k the excitation force deviates by one random step (quantity) from its preceding value and that these steps are distributed independently and identically in terms of size.

This random-walk model of the excitation force is then combined with a dynamics model of the float:

$\hspace{1em}\{\begin{array}{c}x\left(k+1\right)=A_dx\left(k\right)+B_dw\left(k\right)-B_du\left(k\right)+{ϵ}_x\left(k\right)\\ y\left(k\right)=C_dx\left(k\right)+D_dw\left(k\right)-D_du\left(k\right)+\mu \left(k\right)\end{array}$

which is more resistant to modelling errors, owing to the introduction of ϵ_x(k) which represents the unmodelled dynamics (friction of the PTO, hydrostatic non-linearity, . . . ) and μ(k) which describes the measurement noise affecting the position and speed of the float.

By combining the random-walk model of the excitation force of the incident waves and the dynamics of the float. the increased state system is obtained:

$\hspace{1em}\{\begin{array}{c}x\left(k+1\right)=A_dx\left(k\right)+B_dw\left(k\right)-B_du\left(k\right)+{ϵ}_x\left(k\right)\\ w\left(k+1\right)=w\left(k\right)+{ϵ}_m\left(k\right)\\ y\left(k\right)=C_dx\left(k\right)+D_dw\left(k\right)-D_du\left(k\right)+\mu \left(k\right)\end{array}$

or, in an equivalent manner

$\{\begin{array}{c}{x}_a\left(k+1\right)=A_ax_a\left(k\right)+B_au\left(k\right)+ϵ\left(k\right)\\ y\left(k\right)=C_ax_a\left(k\right)+D_au\left(k\right)+\mu \left(k\right)\end{array}\mathrm{where}x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],ϵ\left(k\right)=\left[\begin{array}{c}{ϵ}_x\left(k\right)\\ {ϵ}_m\left(k\right)\end{array}\right]\mathrm{and}A_a=\left[\begin{array}{cc}{A}_d& B_d\\ 0& 1\end{array}\right],B_a=\left[\begin{array}{c}-B_d\\ 0\end{array}\right],C_a=\left[\begin{array}{cc}{C}_d& D_d\end{array}\right],D_a=-D_d$

Thus the problem of estimating w(k) becomes a state estimation problem.

One way of estimating the unknown state vector x_a(k), taking into account information about ϵ(k) and μ(k), is to apply the linear Kalman filter algorithm, as shown below. An alternative may be to use the moving-horizon estimation approach (the same approach used in the third embodiment to estimate the excitation force, viewed this time as an unknown input).

The linear Kalman filter (KF) algorithm provides a solution for the following optimization problem:

$\underset{x_a\left(k\right)}{\min }=\left\{{\left(x_a\left(0\right)-{\stackrel{\_}{x}}_a\left(0\right)\right)}^TP_0^{-1}\left(x_a\left(0\right)-{\stackrel{\_}{x}}_a\left(0\right)\right)+\sum _{j=1}^k{ϵ\left(j-1\right)}^TQ^{-1}ϵ\left(j-1\right)+{\mu \left(j\right)}^TR^{-1}\mu \left(j\right)\right\}$

where P₀, Q and R are calibration matrices of suitable dimensions, {circumflex over (x)}_a(0) is the mean value of the initial state x_a(0).

Considering x_a(0) to be a random vector, not related to ϵ(k) and μ(k), of mean {circumflex over (x)}_a(0) with covariance matrix P₀, and ϵ(k) and μ(k) as two random processes of the mean-zero white noise type with covariance matrices Q and R respectively, the problem of optimization may be solved using the linear Kalman filter algorithm.

At each instant k, the linear Kalman filter algorithm calculates the solution to this problem by means of two steps.

The first step is the temporal updating of the estimations:

$\hspace{1em}\{\begin{array}{c}{\hat{x}}_a\left(k|k-1\right)=A_a{\hat{x}}_a\left(k-1|k-1\right)+B_aF_u\left(k-1\right)\\ P\left(k|k-1\right)=A_aP\left(k-1|k-1\right)A_a^T+Q\end{array}$

where {circumflex over (x)}_a(k|k−1) and P(k|k−1) are respectively the estimation of x_a(k) and its covariance matrix obtained using the measurements from the instant k−1 (namely y(k−1), y(k−2), . . . and F_u(k−1), F_u(k−2), . . . ), and {circumflex over (x)}_a(k|k) and P(k|k) are the estimation of x_a(k) and its covariance matrix obtained using the measurements from the instant k (namely y(k), y(k−1), . . . and F_u(k), F_u(k−1), . . . ).

The second step is the updating of the measurements:

$\hspace{1em}\{\begin{array}{c}K\left(k\right)=P\left(k|k-1\right){C_a^T\left(C_aP\left(k|k-1\right)C_a^T+R\right)}^{-1}\\ {\hat{x}}_a\left(k|k\right)={\hat{x}}_a\left(k|k-1\right)+K\left(k\right)\left(y\left(k\right)-C_a{\hat{x}}_a\left(k|k-1\right)-D_aM_u\left(k\right)\right)\\ P\left(k|k\right)=\left(I-K\left(k\right)C_a\left(k\right)\right)P\left(k|k-1\right)\end{array}$

with K(k) the gain of the observer, which is a parameter vector which weights the errors between the estimations and the measurements of the outputs of the system (position and speed). In the case of the Kalman filter, it is not a calibration parameter, but it is obtained from the estimation of the covariance matrix.

To summarize, for the second embodiment, the following steps may be performed:

• • Initialization of k=0, {circumflex over (x)}_a(0|0)={circumflex over (x)}_a(0), P(0|0)=P₀
  • At each instant k:
    • the following are used:
      • the measurements of the position and speed of the moving part y(k)=[z(k) ż(k)] and the force applied by the PTO to the moving part u(k)=F_u (k)
      • the results of the estimations of the preceding step {circumflex over (x)}_a(k−1|k−1), P(k−1|k−1)
      • the parameters Q, R (covariance matrices)
    • determination of the excitation force F_ex applied by the waves incident upon the moving part, known for this embodiment ŵ(k), by performing the following steps:
      • the two steps of the linear Kalman filter algorithm are applied to obtain {circumflex over (x)}_a(k|k), P(k|k), thus the complete state with its covariance matrix are estimated:

$\hspace{1em}\{\begin{array}{c}{\hat{x}}_a\left(k|k-1\right)=A_a{\hat{x}}_a\left(k-1|k-1\right)+B_aF_u\left(k-1\right)\\ P\left(k|k-1\right)=A_aP\left(k-1|k-1\right)A_a^T+Q\end{array},\begin{array}{c}K\left(k\right)=P\left(k|k-1\right){C_a^T\left(C_aP\left(k|k-1\right)C_a^T+R\right)}^{-1}\\ {\hat{x}}_a\left(k|k\right)={\hat{x}}_a\left(k|k-1\right)+K\left(k\right)\left(y\left(k\right)-C_a{\hat{x}}_a\left(k|k-1\right)-D_aF_u\left(k\right)\right)\\ P\left(k|k\right)=\left(I-K\left(k\right)C_a\left(k\right)\right)P\left(k|k-1\right)\end{array}$

• • • • where x_a(k) is the unknown state vector

$x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],$

• • • •  A_a, B_a, C_a, D_a are state realization matrices, P is the covariance matrix of the state vector, and Q and R are calibration matrices.
      • and extraction of the last component of the last component of the estimated state vector, namely the excitation force F_ex applied by the waves incident upon the moving part, by means of an equation with the form: ŵ(k)=[0 1] {circumflex over (x)}_a(k|k)

Third Embodiment

According to the third embodiment of the invention, the excitation force applied by waves incident upon the moving part may be determined by means of an unknown input state estimator based on a moving-horizon approach. The moving-horizon estimation is an algorithm which determines recursively the state and the unknown input using data belong to a finite time window.

For this embodiment, prior to this step, the following are therefore known:

• • the measurements of the position z, and speed ż of the moving part,
  • the measurement of the force F_u applied by the conversion machine (PTO) to the moving part,
  • the model of the radiation force F_rad(t)=−M_∞{umlaut over (z)}(t)−F_r(t) where F_r(t)=∫₀^th(t−τ)ż(τ)dτ, and
  • the dynamics model of the wave-energy system {umlaut over (z)}(t)=F_ex(t)+F_hd(t)+F_rad(t)−F_u(t).

The second embodiment is completely linear and does not allow a non-linear (more precise) description of the hydrostatic restoring force to be taken into account directly.

This third embodiment allows this limitation to be overcome. Unlike the second embodiment, the excitation force is no longer regarded as a global state model which represents the dynamics of the float, the evolution of which is described by a random-walk model, but is regarded as an unknown input of a (slightly different) global model representing the dynamics of the float.

By combining the equation which describes the movement of the float according to the linear wave theory

M{umlaut over (z)}(t)=F_ex(t)+F(t)+F_rad(t)−F_u(t)

and the expression for the radiation force

F _rad(t)=−M_∞{umlaut over (z)}(t)−F_r(t)

the following is obtained

(M+M_∞){umlaut over (z)}(t)+F_r(t)=w(t)+u_d(t)

where w(t)=F_ex(t) and u_d(t)=F_hd(t)−F_u(t). u_d(t) is regarded as a control input of the system described, known because F_u(t) is measurable and F_hd(t) may be calculated from a from a (linear or non-linear) model of the hydrostatic restoring force, while w(t) is an unknown input (disturbance).

This equation may be expressed in state form, by defining

$\hspace{1em}\{\begin{array}{c}{x}_1\left(t\right)=z\left(t\right)\\ x_2\left(t\right)=\stackrel{.}{z}\left(t\right)\end{array}$

which gives

$\hspace{1em}\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=-\frac{1}{M+M_{\infty }}F_r\left(t\right)+\frac{1}{M+M_{\infty }}\left(w\left(t\right)-u_d\left(t\right)\right)\end{array}$

In the same way as for the second embodiment, using

F _r(t)=∫₀^th(t−τ)ż(τ)dτ

the radiation damping F_r(t) may be expressed as the output of a linear system in the form of a state-space representation, the input of which is x₂(t)=ż(t)

either

$\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_rx_r\left(t\right)+B_r\stackrel{.}{z}\left(t\right)\\ F_r\left(t\right)=C_rx_r\left(t\right)+D_r\stackrel{.}{z}\left(t\right)\end{array}\mathrm{or}\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_rx_r\left(t\right)+B_rx_2\left(t\right)\\ F_r\left(t\right)=C_rx_r\left(t\right)+D_rx_2\left(t\right)\end{array}$

where x_r is a (non-measurable) internal state not having any particular physical significance and (A_r, B_r, C_r, D_r) are the state realization matrices By combining the two state-space representations, the following is obtained:

$\hspace{1em}\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=-\frac{D_r}{M+M_{\infty }}x_2\left(t\right)-\frac{C_r}{M+M_{\infty }}x_r\left(t\right)+\frac{1}{M+M_{\infty }}\left(w\left(t\right)-u_d\left(t\right)\right)\\ {\stackrel{.}{x}}_r\left(t\right)=B_rx_2\left(t\right)+A_rx_r\left(t\right)\end{array}$

or, in an equivalent manner

$\hspace{1em}\{\begin{array}{c}\stackrel{.}{x}\left(t\right)=A_cx\left(t\right)+B_c\left(w\left(t\right)+u_d\left(t\right)\right)\\ y\left(t\right)=C_cx\left(t\right)\end{array}\mathrm{where}x\left(t\right)={\left[\begin{array}{ccc}{x}_1\left(t\right)& x_2\left(t\right)& x_r^T\left(t\right)\end{array}\right]}^Ty=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]\mathrm{and}A_c=\left[\begin{array}{ccc}0& 1& 0\\ 0& -\frac{D_r}{M+M_{\infty }}& -\frac{C_r}{M+M_{\infty }}\\ 0& B_r& A_r\end{array}\right],B_c=\left[\begin{array}{c}0\\ \frac{1}{M+M_{\infty }}\\ 0\end{array}\right],C_c=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

It is attempted to estimate w(t) using y(t) (the vector which contains the measurements of the position and the speed of the float) and u_d(t). It is pointed out that u_d(t)=F_hd(t)−F_u(t) and that it is possible to use, for example, the piecewise affine approximation in order to calculate the hydrostatic restoring force:

F _hd =−F _i z(t)−G_i for each interval Ω_i, to which z belongs.

To do this, it is first possible to discretize the abovementioned state system (it should be remembered the measurements are sampled), using Tustin's method. This gives:

$\hspace{1em}\{\begin{array}{c}x\left(k+1\right)=A_dx\left(k\right)+B_dw\left(k\right)+B_du_d\left(k\right)\\ y\left(k\right)=C_dx\left(k\right)+D_dw\left(k\right)+D_du_d\left(k\right)\end{array}$

To make the model more realistic by adding uncertainty factors (noise) to the state and the output:

$\hspace{1em}\{\begin{array}{c}x\left(k+1\right)=A_dx\left(k\right)+B_dw\left(k\right)+B_du_d\left(k\right)+ϵ\left(k\right)\\ y\left(k\right)=C_dx\left(k\right)+D_dw\left(k\right)+D_du_d\left(k\right)+\mu \left(k\right)\end{array}$

where ϵ(k) and μ(k) are mean-zero white noise.

In order to estimate the unknown input w(k) (and secondarily the state x(k)) the moving-horizon estimation (MHE) method is used. The problem consists in estimating, at each time instant k≥N, the state vector x(k) and the unknown input w(k)=F_ex(k), on the basis of:

• • a priori estimations: {circumflex over (x)}(k−N|k−1), ŵ(k−N|k−1) (the estimations of x(k−N) and w(k−N) respectively) using the measurements from the time k−1;
  • information about the outputs: y(k−N), y(k−N+1), . . . , y(k); and
  • information about the known input: u_d(k−N), u_d(k−N+1), . . . , u_d(k)

Namely:

{circumflex over (x)} ₀ ={circumflex over (x)}(k−N|k−1)

ŵ ₀ =ŵ(k−N|k−1)

and:

x _j ={circumflex over (x)}(k−N+j|k)

w _j =ŵ(k−N+j|k)

ϵ_j={circumflex over (ϵ)}(k−N+j|k)

for j=0, 1, . . . , N.

The moving-horizon estimation problem is formulated as:

$\underset{x_0,w_0,w_1,\dots ,{ϵ}_1,\dots ,{ϵ}_N}{\min }J\left(x_0,w_0,w_1,\dots ,{ϵ}_1,\dots ,{ϵ}_N\right)$

with the following objective function (cost):

$J=(·)={x_0-{\hat{x}}_0}_{P_0^{-1}}^2+{w_0-{\hat{w}}_0}_{Q_0^{-1}}^2+\sum _{j=0}^{N-1}{{ϵ}_j}_{Q^{-1}}^2++\sum _{j=0}^N{y\left(k-N+j\right)-\hat{y}\left(k-N+j\right)}_{R^{-1}}^2+\sum _{j=0}^{N-1}{w_j-w_{j-1}}^2$

where

ŷ(k−N+j)=C_dx_j+D_dw_j+D_du_d(k−N+j)

The notation ∥v∥_M² indicates the product v^TMv for a vector of reals v and a symmetrical matrix of reals M.

There are five terms in the cost function:

• • The first term

${x_0-{\hat{x}}_0}_{P_0^{-1}}^2$

and a second term

${w_0-{\hat{w}}_0}_{Q_0^{-1}}^2$

added together form the “arrival cost” which summarizes information prior to the time k−N and forms part of the data of the estimation problem.

• • The third term and the fourth term penalize the method and measurement noise, respectively, with the weighting matrices Q and R.
  • The last term is a regularization term which has the function of taking into account the fact that the excitation force of the waves is a regular (smooth) signal.

This latter term is important for determining the force applied by waves in an optimum manner. Another way of taking into account the “smooth” nature of the excitation force of waves would consist in imposing variation constraints thereon.

Using the classic approach of a moving-horizon estimation and control, the minimization of the cost function may be transformed into a quadratic programming QP problem (a quadratic programming problem is an optimization problem in which a quadratic function is minimized/maximized on a convex polyhedron) by propagating:

$\hspace{1em}\{\begin{array}{c}x\left(k+1\right)=A_dx\left(k\right)+B_dw\left(k\right)+B_du_d\left(k\right)+ϵ\left(k\right)\\ y\left(k\right)=C_dx\left(k\right)+D_dw\left(k\right)+D_du_d\left(k\right)+\mu \left(k\right)\end{array}$

• • over time. This gives:

$x_1=A_dx_0+B_dw_0+B_du_d\left(k-N\right)+{ϵ}_0$ $\begin{array}{c}x_2=A_dx_1+B_dw_1+B_du_d\left(k-N+1\right)+{ϵ}_1\\ =A_d^2x_0+A_dB_dw_0+B_dw_1+A_dB_du_d\left(k-N\right)+\\ B_du_d\left(k-N+1\right)+A_d{ϵ}_0+{ϵ}_1\end{array}$ $⋮$ $x_n=A_d^Nx_0+A_d^{N-1}B_dw_0+\dots +B_dw_{N-1}+A_d^{N-1}B_du_d\left(k-N\right)+\dots +B_du_d\left(k-1\right)+A_d^{N-1}{ϵ}_0+\dots +{ϵ}_1$

• • or, in compact form,

$X={\Phi }_xx_0+{\Psi }_xW+{\Psi }_xU_d+{\Gamma }_xE$ $\mathrm{where}$ $X=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ A_N\end{array}\right],W=\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ W_N\end{array}\right],U_d=\left[\begin{array}{c}{u}_d\left(k-N\right)\\ u_d\left(k-N+1\right)\\ ⋮\\ u_d\left(k\right)\end{array}\right],E=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_N\end{array}\right]$ $\mathrm{and}$ ${\Phi }_x=\left[\begin{array}{c}{A}_d\\ A_d^2\\ ⋮\\ A_d^N\end{array}\right],{\Psi }_x=\left[\begin{array}{ccccc}{B}_d& 0& \dots & 0& 0\\ A_dB_d& B_d& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ A_d^{N-1}B_d& A_d^{N-2}B_d& \dots & B_d& 0\end{array}\right],{\Gamma }_x=\left[\begin{array}{cccc}I& 0& \dots & 0\\ A_d& I& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ A_d^{N-1}& A_d^{N-2}& \dots & I\end{array}\right]$

Applying the same approach to the output equation and introducing the notations

$Y={\left[{y\left(k-N+1\right)}^T{y\left(k-N+2\right)}^T{y\left(k\right)}^T\right]}^T,{\Phi }_y=\left[\begin{array}{c}{C}_dA_d\\ C_dA_d^2\\ ⋮\\ C_dA_d^N\end{array}\right],{\Psi }_y=\left[\begin{array}{ccccc}{C}_dB_d& D_d& \dots & 0& 0\\ C_dA_dB_d& C_dB_d& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ C_dA_d^{N-1}B_d& C_dA_d^{N-2}B_d& \dots & C_dB_d& D_d\end{array}\right],{\Gamma }_y=\left[\begin{array}{cccc}{C}_d& 0& \dots & 0\\ C_dA_d& C_d& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ C_dA_d^{N-1}& C_dA_d^{N-2}& \dots & C_d\end{array}\right]$

results in a QP problem with the form:

$\{\begin{array}{c}\underset{\xi }{\min }\left\{{\xi }^T\left(P+{\Psi }_s^TR{\Psi }_s\right)\xi -2\left[{\hat{w}}_0^T{\hat{x}}_0^TV^T\right]\left[\begin{array}{c}S\\ R{\Psi }_s\end{array}\right]\xi \right\}\\ \mathrm{such}\mathrm{as}H_{\xi }\xi \le g_{\xi }+H_dU_d\end{array}\mathrm{where}\{\begin{array}{c}\xi =\left[W^Tx_0^TE^T\right]\\ V=Y-{\Psi }_yU_d\\ {\Psi }_s=\left[{\Psi }_y{\Phi }_y{\Gamma }_y\right]\end{array}$

where ξ is a matrix which represents all the unknown variables to be estimated (the values of the unknown input W^T and the uncertainty about the state E^T for all the N steps of the estimation window considered, and the initial state x₀^T at the time k−N), Ψ_s is a matrix which includes the effects of these same unknown variables on the values output in the estimation window Y, Ψ_y, is a matrix which corresponds to the effects of the known input U_d on Y, Ψ_y Ψ_y and Ψ_s calculated based on the state-space representation defined above, and P, R and S are matrices formed from weighting matrices of the cost function Q, R, Q₀⁻¹, R₀⁻¹ and the regularization parameter λ as follows:

$Q=\left[\begin{array}{cccc}{Q}^{-1}& 0& \dots & 0\\ 0& Q^{-1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & Q^{-1}\end{array}\right],R=\left[\begin{array}{cccc}{R}^{-1}& 0& \dots & 0\\ 0& R^{-1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & R^{-1}\end{array}\right],{\Lambda }_1=\left[\begin{array}{ccccccc}{Q}_0^{-1}-\lambda & \lambda & 0& 0& \dots & 0& 0\\ 0& -\lambda & \lambda & 0& \dots & 0& 0\\ 0& 0& -\lambda & \lambda & \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \dots & \lambda & 0\\ 0& 0& 0& 0& \dots & -\lambda & \lambda \end{array}\right],P=\left[\begin{array}{ccc}{\Lambda }_1& 0& 0\\ 0& P_0^{-1}& 0\\ 0& 0& Q\end{array}\right],S=\left[\begin{array}{cccc}{Q}_0^{-1}& 0& \dots & 0\\ 0& P_0^{-1}& \dots & 0\end{array}\right]$

To conclude

H _ξ ξ≤g _ξ +H _d U _d

is a linear inequality which, with the aid of the matrices H_ξ and H_d and the vector g_ξ, retranscribes in compact form polytopic constraints on the unknown input and state:

$\hspace{1em}\{\begin{array}{c}w\in P_w,P_w=\left\{w\text{:}H_w\le g_w\right\},g_w>0\\ x\in P_x,P_x=\left\{x\text{:}H_x\le g_x\right\},g_x>0\end{array}$

which in turn are defined by the matrices H_w and H_x and the vectors g_w and g_x. With this inequality it is possible to take into account information about the minimum and maximum amplitudes of the excitation force of the incident waves, and the position and speed of the float.

The moving-horizon estimation problem is a constrained quadratic programming problem which may be resolved with solvers such as FORCES® (Embotech, Suisse) or CVXGEN® (CVXGEN LLC, USA) for example.

Let it be assumed that ξ* is the solution of the dynamic programming problem obtained with one of these solvers. The estimation of the excitation forces, using measurements at the time k, is given by:

ŵ(k−N+j|k)=ξ*(j),∀i=1,2, . . . ,N

To summarize, for the third embodiment, the steps described below may be performed.

The components of the vectors of the past estimations and measurements corresponding to the N−1 time step preceding the present instant k=0 are initialized to zero.

At each instant k,

• • 1. The following are acquired:
    • the measurements of the position and the speed of the moving part y(k)=[z(k) ż(k)] and of the force applied by the PTO to the moving part u(k)=F_u (k),
    • the results of the estimations of the preceding time step {circumflex over (x)}_d(k−N|k−1), ŵ(k−N|k−1)
    • the parameters Q, R, Q₀⁻¹, R₀⁻¹ (weighting matrices) and λ (regularization)
  • 2. The calculation ŵ(k|k) is carried out, this corresponding to the estimation of the force applied by waves to the moving part, by performing the following steps:
    • i. The following quadratic programming problem QP is resolved to obtain the solution ξ*

$\hspace{1em}\{\begin{array}{c}\underset{\xi }{\min }\left\{{\xi }^T\left(P+{\Psi }_s^TR{\Psi }_s\right)\xi -2\left[{\hat{w}}_0^T{\hat{x}}_0^TV^T\right]\left[\begin{array}{c}S\\ R{\Psi }_s\end{array}\right]\xi \right\}\\ \mathrm{such}\mathrm{as}H_{\xi }\xi \le g_{\xi }+H_dU_d\end{array}$

• • • ii. the excitation force applied by waves incident upon the moving part w(k) is estimated by means of an equation with the form:

ŵ(k|k)=ξ*(N).

6) Controlling of the Wave-Energy System

During this step the wave-energy system is controlled so as to take into account the excitation force applied by the incident wave. This it is possible to drive the wave-energy system so as to optimize the recovered energy.

Controlling may consist of control of the moving part of the wave-energy system, for example by means of an electric, pneumatic or hydraulic conversion machine, called PTO (power take-off system). This PTO system influences the movement of the moving part and allows the mechanical energy to be transferred to the electrical, pneumatic or hydraulic network. Model predictive control (MPC) is an example of a method for controlling wave-energy systems.

Illustrative Example

In order to illustrate the advantages of the method according to the invention, the method according to the invention was tested in a tank. The wave-energy system corresponds to the wave-energy system shown in FIG. 1. The wave-energy system is equipped with a torque sensor for measuring the torque applied by the wave to the moving part. Beforehand, the torque applied by the wave is determined in the case where control is not active (namely the conversion machine does not generate any torque on the system) and the float is kept in an equilibrium position. These characteristics allows the excitation torque applied by the wave incident upon the moving part for this first test to be compared with other tests when the float is moving and the electric machine is activated.

For the present analysis a series of four types of different waves is considered. These four types of wave W1, W2, W3, W4 (from smallest to largest amplitude) are shown in FIG. 3 which illustrates the amplitude A as a function of the time T.

FIGS. 4a to 4c show a comparison of the excitation reference force (in this case the torque) applied by the wave incident upon the moving part, with the excitation force (in this case the torque) applied by the wave incident upon the moving part determined according to an embodiment of the invention. FIGS. 4a to 4c correspond to a wave of the type W1 shown in FIG. 3. The reference wave is indicated by REF and shown as a broken line. The forces of the incident wave determined by the method according to the invention are shown as a continuous line. FIG. 4a corresponds to the first embodiment of the invention MR1, FIG. 4b corresponds to the second embodiment of the invention MR2 and FIG. 4c corresponds to the third embodiment of the invention MR3. It can noted in these FIGS. 4a to 4c that, in all three embodiments, the graph representing the force determined by the method according to the invention is more or less superimposed on the graph representing the reference force. Thus, with the method according to the invention, it is possible to determine precisely the excitation force of the wave incident upon the moving part.

FIGS. 5a to 5c show a comparison of the excitation reference force (in this case the torque) applied by waves incident upon the moving part, with the excitation force (in this case the torque) applied by waves incident upon the moving part determined according to an embodiment of the invention. FIGS. 5a to 5c correspond to a wave of the type W3 shown in FIG. 3. The reference wave is indicated by REF and shown as a broken line. The forces of the wave determined by the method according to the invention are shown as a continuous line. FIG. 5a corresponds to the first embodiment of the invention MR1, FIG. 5b corresponds to the second embodiment of the invention MR2 and FIG. 5c corresponds to the third embodiment of the invention MR3. It can noted in these FIGS. 5a to 5c that, in all three embodiments, the graph representing the force determined by the method according to the invention is more or less superimposed on the graph representing the reference force. Thus, with the method according to the invention, it is possible to determine precisely the excitation force of the wave incident upon the moving part. It is also possible to note that the error between the reference force and the force which is determined is smaller for these graphs than for the graphs shown in FIGS. 4a to 4c.

Claims

1. A method for controlling a wave-energy system by determining the excitation force applied by waves incident upon a moving part of the wave-energy system, the wave-energy system converting the energy of the waves into energy by means of the moving part cooperating with a conversion machine, the moving part performing a movement in relation to the conversion machine as a result of the action of the waves, wherein the following steps are performed: a) measurement of the position and speed of the moving part;b) measurement of the force Fu applied by the conversion machine to the moving part;c) construction of a model of the radiation force Frad applied to the moving part, the model of the radiation force Frad relating the radiation force to the speed of the moving part;d) construction of a model of the dynamics of the wave-energy system which relates the excitation force Fex applied by the wave incident upon the moving part, to the position of the moving part, to the speed of the moving part, to the force Fu applied by the conversion machine to the moving part, and to the radiation force Frad applied to the moving part;e) determination of the excitation force Fex applied by the wave incident upon the moving part by means of the dynamics model, the model of the radiation force, the position and speed measured and the measured force Fu applied by the conversion machine to the moving part; andf) controlling of the wave-energy system depending on the determined excitation force Fex applied by the wave incident upon the moving part.

2. The method as claimed in claim 1, wherein the dynamics model of the wave-energy system is constructed by means of an equation of the type: M{umlaut over (z)}(t)=Fex(t)+Fhd(t)+Frad(t)−Fu(t) where M is the mass of the moving part, {umlaut over (z)} is the acceleration of the moving part, Fhd is the hydrostatic restoring force applied to the moving part.

3. The method as claimed in claim 2, wherein the hydrostatic restoring force Fhd is determined by means of a formula of the type: Fhd(t)=−Kz(t) where z is the position of the moving part and K is the hydrostatic stiffness coefficient.

4. The method as claimed in claim 2, wherein the hydrostatic restoring force Fhd is determined by a piecewise affine function of the type: Fhd=−Fiz(t)−Gi for at least two variation intervals Ωi of the position z, where Fi and Gi are constant coefficients of the affine function for each interval Ωi.

5. The method as claimed in claim 1, wherein the model of the radiation force Frad is constructed by an equation of the type: Frad(t)=−M∞{umlaut over (z)}(t)−Fr(t) where Fr(t)=∫oth(t−τ)ż(τ)dτ is the radiation damping, {umlaut over (z)} is the acceleration of the moving part, M∞ is the mass added to the infinitely high frequency, ż is the speed of the moving part, h is the pulse response which relates the speed of the moving part to the radiation damping.

6. The method as claimed in claim 5, wherein the acceleration of the moving part is measured.

7. The method as claimed in claim 6, wherein the radiation damping Fr is determined by means of an approximated time response obtained from a state realization, of the type: Fr(t)=Cr∫0te(t-τ)ArBrż(τ)dτ+Drż(t) where Ar, Br, Cr, Dr are state realization matrices.

8. The method as claimed in claim 2, wherein the excitation force Fex applied by the wave incident upon the moving part is determined by means of an equation of the type: {circumflex over (F)}ex(k)=(M+M∞){umlaut over (z)}(k)+Σl=0k-1CrdArdk-1-lBrdż(l)+Drdż(k)+Fiz(k)+Gi+Fu(k) where k is the time step in the discrete domain, and Ard, Brd, Crd, Drd are the state realization matrices in the discrete domain.

9. The method as claimed in claim 5, wherein the excitation force Fex applied by the wave incident upon the moving part is determined by means of a state observer based on a linear Kalman filter constructed from a random-walk model of the excitation force Fex applied by the said wave incident upon the said moving part.

10. The method as claimed in claim 9, wherein the excitation force Fex applied by the wave incident upon the moving part is determined by performing the following steps: i) initialization of k=0, the state vector {circumflex over (x)}a(0|0)=xa(0) and the state of the covariance matrix, P(0|0)=P0;ii) at each instant k, acquisition of the position and speed measurements of the moving part y(k)=[z(k) ż(k)] and the measurement of the force exerted by the conversion machine on the moving part u(k)=Fu(k); andiii) at each instant k, determination of the excitation force applied by the wave incident upon the moving part (2) ŵ(k) by means of the following equations: ŵ(k)=[0 1]{circumflex over (x)}a(k|k) where:

11. The method as claimed in claim 5, wherein the excitation force FP applied by the wave incident upon the moving part is determined by means of an unknown input state estimator based on a moving-horizon approach.

12. The method as claimed in claim 11, wherein the excitation force Fex applied by the wave incident upon the moving part is determined by performing the following steps at each instant k: iv) acquisition of the position and speed measurements of the moving part y(k)=[z(k) ż(k)] and the measurement of the force applied by the conversion machine to the moving part u(k)=Fu(k); andv) determination of the excitation force applied by the wave incident upon the moving part ŵ(k|k) by means of the following equations:ŵ(k|k)=ξ(N) where ξ(N) is obtained by resolution of a following quadratic programming problem: