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).

Claims

1. A method for controlling a wave-energy system for generating exploitable energy, the system including a floating moving part which moves in response to waves incident on the floating moving part and which drives a power take-off system comprising one of an electric, pneumatic or hydraulic conversion machine with an excitation force Fex produced in response to the waves incident on the floating moving part, a moveable member pivotally attached to the floating moving part which is coupled to the power take-off system for influencing the moving part by transferring energy from the power take-off system to the floating moving part, and allowing transfer of mechanical energy from the floating moving part, to an electrical, pneumatic or hydraulic network, the method comprising: a) measuring position and speed of the floating moving part in relation to an equilibrium position thereof; b) measuring a force Fu applied by the power take-off system to the floating moving part; c) constructing a linear model of radiation force Frad relating acceleration and speed of the floating moving part; d) constructing a model of dynamics of the wave-energy system which relates the excitation force Fex to the position of the floating moving part, to the speed of the floating moving part, and to the force Fu applied by the power take-off system of the floating moving part; e) determining the excitation force Fex produced in response to the incident waves to the floating moving part by use of the model of dynamics, the model of the radiation force Frad, a measured position and speed of the floating moving part in step a, and the force Fu applied by the power take-off to the floating moving part; and f) controlling the operation of the wave energy system by controlling the power take-off system dependent on the determined excitation force Fex applied by the waves incident on the floating moving part to adjust recovery of exploitable energy as a function of the determined excitation force Fex applied by the waves incident on the floating moving part.

2. The method as claimed in claim 1, wherein the model of dynamics of the wave-energy system is constructed by use of an equation: M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t) where M is the mass of the floating moving part, {umlaut over (z)} is acceleration of the floating moving part, and F.sub.hd is a hydrostatic restoring force applied to the floating moving part.

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

4. The method as claimed in claim 3, wherein radiation damping F.sub.r is determined by use of an approximated time response obtained from a state realization: F.sub.r(t)=C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t) where A.sub.r, B.sub.r, C.sub.r, and D.sub.r are state realization matrices.

5. The method as claimed in claim 2, wherein the hydrostatic restoring force F.sub.hd is determined by a piecewise affine function: F.sub.hd=F.sub.iz(t)G.sub.i for at least two variation intervals .sub.i of a position z, where F.sub.i and G.sub.i are constant coefficients of the piecewise affine function for each interval .sub.i.

6. The method as claimed in claim 5, wherein the excitation force Fex applied by the waves incident to the floating moving part is determined by use of an equation: {circumflex over (F)}.sub.ex(k)=(M+M.sub.){umlaut over (z)}(k)+.sub.l=0.sup.k-1C.sub.rd A.sub.rd.sup.k-1-lB.sub.rd(l)+D.sub.rd(k)+F.sub.iz(k)+G.sub.i+F.sub.u(k) where k is a time step in a discrete domain, and A.sub.rd, B.sub.rd, C.sub.rd, and D.sub.rd are state realization matrices in the discrete domain.

7. The method as claimed in claim 2, wherein acceleration of the floating moving part is measured.

8. The method as claimed in claim 2, wherein excitation force Fex applied by the waves incident on the floating moving part is determined by an equation: {circumflex over (F)}.sub.ex(k)=(M+M.sub.){umlaut over (z)}(k)+.sub.l=0.sup.k-1C.sub.rdA.sub.rd.sup.k-1-lB.sub.rd(l)+D.sub.rd(k)+F.sub.iz(k)+G.sub.i+F.sub.u(k) where k is time step in a discrete domain, and A.sub.rd, B.sub.rd, C.sub.rd, and D.sub.rd are state realization matrices in a discrete domain.

9. The method as claimed in claim 1, wherein the model of the radiation force F.sub.rad is constructed by an equation: F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.0.sup.th(t)()d is radiation damping, {umlaut over (z)} is acceleration of the floating moving part, and M.sub. is mass added to an infinitely high frequency, is speed of the floating moving part, and h is a pulse response relating speed of floating moving part to radiation damping.

10. The method as claimed in claim 9, wherein the acceleration of floating moving part is measured.

11. The method as claimed in claim 10, wherein the radiation damping F.sub.r is determined by use of an approximated time response obtained from a state realization: F.sub.r(t)=C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t) where A.sub.r, B.sub.r, C.sub.r, and D.sub.r are state realization matrices.

12. The method as claimed in claim 11, wherein excitation force Fex applied by the waves incident on the floating moving part is determined by an equation: {circumflex over (F)}.sub.ex(k)=(M+M.sub.){umlaut over (z)}(k)+.sub.l=0.sup.k-1C.sub.rdA.sub.rd.sup.k-1-lB.sub.rd(l)+D.sub.rd(k)+F.sub.iz(k)+G.sub.i+F.sub.u(k) where k is time step in a discrete domain, and A.sub.rd, B.sub.rd, C.sub.rd, and D.sub.rd are state realization matrices in a discrete domain.

13. The method as claimed in claim 9, wherein the excitation force Fex applied by the waves incident on the floating moving part is determined by use of a state observer based on a linear Kalman filter constructed from a random-walk model of the excitation force Fex applied by the waves incident on the floating moving part; and wherein the excitation force Fex applied by the waves incident on the floating moving part is determined by performing steps of: i) initialization of k=0, state vector {circumflex over (x)}.sub.a(0|0)=x.sub.a(0) and state of the covariance matrix, P(0|0)=P.sub.0; ii) at each instant k, acquisition of the position and speed measurements of floating moving part y(k)=[z(k) (k)] and measurement of force exerted by the power take-off on the floating moving part u(k)=F.sub.u(k); and iii) at each instant k, determination of the excitation force Fex applied by incident waves on the floating moving part (2) (k) by use of the equations:
(k)=[01]{circumflex over (x)}.sub.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},\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(1-K\left(k\right)C_a\left(k\right)\right)P\left(k|k-1\right)\end{array}$ where K(k) is gain of an observer, x.sub.a(k) is an unknown state vector, $x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],$ A.sub.a, B.sub.a, C.sub.a, and D.sub.a are state realization matrices, P is a covariance matrix of the unknown state vector, and Q and R are calibration matrices.

14. The method as claimed in claim 9, wherein excitation force Fex applied by the waves incident on the floating moving part is determined by use of an unknown input state estimator based on a moving-horizon method.

15. The method as claimed in claim 14, wherein the excitation force Fex applied by the incident waves on the floating moving part is determined by performing the following steps at each time instant k: i) acquisition of position and speed measurements of the floating moving part and measurement of the force Fex y(k)=[z(k) (k)] and measurement of the force applied by the conversion machine to the floating moving part u(k)=F.sub.u(k); and ii) determination of the excitation force Fex applied by the incident waves on the floating moving part (k|k) by use of the equations: (k|k)=(N) where (N) is obtained by resolution of a following quadratic programming problem: $\{\begin{array}{c}\underset{}{\min }\left\{{}^T\left(P+{}_s^{T}R{}_s\right)-2\left[{\hat{w}}_0^T{\hat{x}}_0^T{V}^T\right]\left[\begin{array}{c}S\\ R{}_s\end{array}\right]\right\}\\ \mathrm{such}\mathrm{as}{H}_{}{g}_{}+H_d{U}_d\end{array}$ where P, R and S are matrices formed from matrices for weighting of a cost function, is a matrix which represents all the values of unknown variables in an estimation window of size N, .sub.s is a matrix, calculated from a state-space representation, corresponding to an effect of variables on output values an identical window, H.sub. and H.sub.d are matrices and g.sub. is a vector obtained from constraints imposed on the state representation and an unknown input, {circumflex over (x)}.sub.0={circumflex over (x)}(kN|k1) and .sub.0=(kN|k1) where x is a state vector and w is an 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)\\ \mathrm{.Math.}\\ u_d\left(k\right)\end{array}\right],$ u.sub.d being the known input of state estimator with u.sub.d(t)=F.sub.hd(t)F.sub.u(t).

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) 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.

(2) FIG. 1 shows a wave-energy system according to an embodiment of the invention.

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

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

(5) 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.

(6) 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

(7) The present invention relates to the determination of the excitation force applied by waves incident upon a moving part of a wave-energy system, also called an 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 a 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 electrical generator connected to a device allowing the transmission of the oscillating movement to be adapted, 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, power may be supplied to the conversion machine to provide a torque or a force which drives the moving element to help the moving platform become synchronized with the waves (motor mode). The method according to the invention relates to the control of a wave-energy system depending on the determined excitation force.

(8) The method according to the invention is suitable for all types of wave-energy systems with at least one moving part, for example those described in the French patent application FR2973448 which corresponds to 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).

(9) 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). Movement of the moving part 2 in response to waves communicates motion via lever arm 4 to a conversion machine 3 which is an electrical generator that as illustrated is a linear generator. The lever arm 4 and the electrical generator 3 are pivotably mounted relative to an installation 6. The movements of the lever arm 4 are indicated in schematic form by use of arrows.

(10) In the description of the invention, the terms waves, sea tides and swell are considered to be equivalent.

(11) 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.

(12) 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.

(13) The method according to the invention comprises the following steps: 1) measuring of the position and speed of the moving part; 2) measuring of the force applied by the conversion machine (or PTO) 3) constructing a radiation force model 4) constructing a model of dynamics 5) determining an excitation force of an incident wave 6) controlling the wave-energy system

(14) 1) Measurement of the Position and Speed of the Moving Part

(15) 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 sensors, which are generally provided on a wave-energy system for at least one of control and management thereof.

(16) 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.

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

(18) 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 at least one of control and management thereof. Alternatively, the measurement may be replaced by an estimation carried out based on a set force (or torque) value applied to the PTO.

(19) 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.

(20) 3) Construction of a Radiation Force Model

(21) 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, Kurniawan 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 use 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, which means 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).

(22) According to a mode of implementation of the invention, the model of the radiation force Frad is constructed by an equation of the type:
F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t);

(23) where F.sub.r(t)=.sub.0.sup.th(t)()d is the component of the radiation force F.sub.rad which depends on the (current and past) speed of the moving part, which is called radiation damping;

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

(25) M.sub. 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;

(26) is the speed of the moving part;

(27) 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.

(28) 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.

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

(30) During this step, a model of dynamics 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 to obtain results representing the behavior of the wave-energy system, if the movements are not too big.

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

(32) 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.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t)

(33) where M is the mass of the moving part;

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

(35) F.sub.ex is the excitation force applied by waves incident upon the moving part;

(36) F.sub.rad is the radiation force applied to the moving part;

(37) F.sub.hd is the hydrostatic restoring force applied to the moving part; and

(38) F.sub.u is the force applied by the conversion machine to the moving part.

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

(40) 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.sub.hd(t)=Kz(t) where z is the position of the 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.

(41) 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 twopreferably threeintervals of the position z. In this case, the hydrostatic restoring force may be written by a function of the type: F.sub.hd=F.sub.iz(t)G.sub.i for at least two variation intervals .sub.i of the position z, where Fi and Gi are constant coefficients of the affine function for each interval .sub.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 Fhd as a function of the position z. This figure shows in broken lines the linear functions of each piece .sub.1, .sub.2 and .sub.3. 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.

(42) According to a mode of implementation of the invention, where the torques are considered, the model of dynamics of the wave-energy system may be constructed by an equation of the type:
J{umlaut over ()}=M.sub.ex(t)+M.sub.hyd(t)+M.sub.rad(t)M.sub.u(t)

(43) where J is the inertia of the moving part;

(44) {umlaut over ()} is the angular acceleration of the moving part;

(45) M.sub.ex is the excitation torque applied by waves incident upon the moving part;

(46) M.sub.rad is the radiation torque applied to the moving part;

(47) M.sub.hd is the hydrostatic restoring torque applied to the moving part; and

(48) M.sub.u is the torque applied by the conversion machine to the moving part.

(49) 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.

(50) 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.

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

(52) 5) Determination of the Excitation Force Applied by Incident Waves

(53) During this step, the excitation force applied by waves incident upon the moving part is determined by use 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.

(54) 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.

(55) 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 use 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.

(56) 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 use 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.

(57) Other embodiments may be envisaged.

First Embodiment

(58) 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.

(59) 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.sub.u applied by the conversion machine (PTO) to the moving part, the model of the radiation force F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.0.sup.th(t)()d, and the model of dynamics of the wave-energy system M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t).

(60) Thus, from knowing how to calculate F.sub.r it is therefore possible to obtain F.sub.rad, and finally F.sub.ex(t), which is the only unknown factor of the dynamics model of the wave-energy system.

(61) However, it is not easy to perform the calculation of the convolution product within the time domain in the preceding expression for F.sub.r(t), because of 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.sub.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.sub.r(s)=W.sub.r(s)(s)

(62) where W.sub.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:

(63) $\{\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}$

(64) where x.sub.r is an internal state which does not have any particular physical significance and (A.sub.r, B.sub.r, C.sub.r, and D.sub.r) are state realization matrices. It is to be noted 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.

(65) The state form allows the evolution in time of F.sub.r to be expressed as follows:
F.sub.r(t)=C.sub.re.sup.tA.sup.rx.sub.r(0)+C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t)

(66) but this expression may not be used directly because the initial value of the internal state x.sub.r(0) is unknown. However, it is possible to approximate the time evolution F.sub.r as follows:
F.sub.r(t)=C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t)

(67) since C.sub.re.sup.tA.sup.rx.sub.r(0)0 where t is sufficiently great, given that the matrix A.sub.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.sub.r, knowing that it will recede very rapidly compared to the contribution linked to the speed .

(68) The first embodiment therefore calculates F.sub.ex from z, , {umlaut over (z)} and F.sub.u (which are known) using an equation with the form:
{circumflex over (F)}.sub.ex(t)=(M+M.sub.){umlaut over (z)}(t)+C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r(t)d+D.sub.r(t)+Kz(t)+F.sub.u(t)

(69) if a linear model is used for the hydrostatic restoring force, or
{circumflex over (F)}.sub.ex(t)=(M+M.sub.){umlaut over (z)}(t)+C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t)+F.sub.iz(t)+G.sub.i+F.sub.u(t)

(70) for each interval .sub.i to which z belongs, if a piecewise affine model is used for the hydrostatic restoring force.

(71) Since in practice the measurements of z, , {umlaut over (z)} and F.sub.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.sub.r:

(72) $\{\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_{rd}{x}_r\left(k\right)+D_{rd}\stackrel{.}{z}\left(k\right)\end{array}$

(73) F.sub.r may thus be expressed by an equation with the form:

(74) $F_r\left(k\right)=C_{rd}{A}_{rd}^k{x}_r\left(0\right)+\underset{l=0}{\overset{k-1}{.Math.}}{C}_{rd}{A}_{rd}^{k-1-l}{B}_{rd}\stackrel{.}{z}\left(l\right)+D_{rd}\stackrel{.}{z}\left(k\right)$

(75) and its approximation expressed by the formula:

(76) $F_r\left(k\right)=\underset{l=0}{\overset{k-1}{.Math.}}{C}_{rd}{A}_{rd}^{k-1-l}{B}_{rd}\stackrel{.}{z}\left(l\right)+D_{rd}\stackrel{.}{z}\left(k\right)$

(77) This gives (in the case where the hydrostatic restoring force is considered to be a piecewise affine function) the following expression:
{circumflex over (F)}.sub.ex(k)=(M+M.sub.){umlaut over (z)}(k)+.sub.t=0.sup.k-1C.sub.rdA.sub.rd.sup.k-1-iB.sub.rd(l)+D.sub.rd(k)F.sub.iz(k)G.sub.iF.sub.u(k)

(78) for each interval .sub.i to which z(k) belongs.

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

Second Embodiment

(80) According to a second embodiment of the invention, the excitation force applied by waves incident upon the moving part may be determined by use 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.

(81) For this embodiment, prior to this step, the following are therefore known: the measurements of the position z, and speed 2 of the moving part, the measurement of the force F.sub.u applied by the conversion machine (PTO) to the moving part, the model of the radiation force F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.0.sup.th(t)()d, and the model of dynamics of the wave-energy system M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t).

(82) 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.

(83) By replacing in the equation which describes the movement of the float
M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t)

(84) the expressions for the hydrostatic restoring force (with the linear model) and the model of the radiation force
F.sub.hd(t)=Kz(t)
F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t)

(85) the following is obtained
(M+M.sub.){umlaut over (z)}(t)+Kz(t)+F.sub.r(t)=F.sub.ex(t)F.sub.u(t)

(86) This equation may be expressed in state form, by defining

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

(88) which gives

(89) 0$\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=-\frac{K}{M+M_{}}{x}_1\left(t\right)-\frac{1}{M+M_{}}{F}_r\left(t\right)+\frac{1}{M+M_{}}\left(F_{\mathrm{ex}}\left(t\right)-F_u\left(t\right)\right)\end{array}$

(90) In the same way as for the first embodiment, the dynamics of radiation damping F.sub.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.sub.r(t)=.sub.0.sup.th(t)()d

(91) namely

(92) $\{\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}$

(93) where x.sub.r is an internal (non-measurable) state which does not have any particular physical significance, and (A.sub.r, B.sub.r, C.sub.r, D.sub.r) are matrices of the state realization.

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

(95) $\{\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_{}}{x}_1\left(t\right)-\frac{D_r}{M+M_{}}{x}_2\left(t\right)-\frac{C_r}{M+M_{}}{x}_r\left(t\right)+\\ \frac{1}{M+M_{}}\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}$

(96) or, in an equivalent manner

(97) $\{\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_{}}& -\frac{D_r}{M+M_{}}& -\frac{C_r}{M+M_{}}\\ 0& B_r& A_r\end{array}\right],B_c=\left[\begin{array}{c}0\\ \frac{1}{M+M_{}}\\ 0\end{array}\right],C_c=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

(98) It is necessary to estimate F.sub.ex(t) using y(t) (the vector which contains the position and the speed measurements of the moving part) and F.sub.u(t), the measurement of the force of the PTO.

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

(100) $\{\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}$

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

(102) By defining
w(k)=F.sub.ex(k),u(k)=F.sub.u(k),

(103) the global model of the float dynamics is obtained:

(104) $\{\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}$

(105) 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(k1).

(106) For this embodiment, a random-walk model is used to describe the evolution of the excitation force:
w(k+1)=w(k)+.sub.m(k)

(107) where .sub.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.

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

(109) $\{\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)+\left(k\right)\end{array}$

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

(111) 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:

(112) $\{\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)+\left(k\right)\end{array}$

(113) or, in an equivalent manner

(114) $\{\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)+\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$

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

(116) One way of estimating the unknown state vector x.sub.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).

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

(118) $\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)+\underset{j=1}{\overset{k}{.Math.}}{\left(j-1\right)}^T{Q}^{-1}\left(j-1\right)+{\left(j\right)}^T{R}^{-1}\left(j\right)\right\}$

(119) where P.sub.0, Q and R are calibration matrices of suitable dimensions, x.sub.a(0) is the mean value of the initial state x.sub.a(0).

(120) Considering x.sub.a(0) to be a random vector, not related to (k) and (k), of mean x.sub.a(0) with covariance matrix P.sub.0, 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.

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

(122) The first step is the temporal updating of the estimations:

(123) 0$\{\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}$

(124) where {circumflex over (x)}.sub.a(k|k1) and P(k|k1) are respectively the estimation of x.sub.a(k) and its covariance matrix obtained using the measurements from the instant k1 (namely y(k1), y(k2), . . . and F.sub.u(k1), F.sub.u(k2), . . . ), and {circumflex over (x)}.sub.a(k|k) and P(k|k) are the estimation of x.sub.a(k) and its covariance matrix obtained using the measurements from the instant k (namely y(k), y(k1), . . . and F.sub.u(k), F.sub.u(k1), . . . ).

(125) The second step is the updating of the measurements:

(126) $\{\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}$

(127) K(k) is 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.

(128) To summarize, for the second embodiment, the following steps may be performed: Initialization of k=0, {circumflex over (x)}.sub.a(0|0)={circumflex over (x)}.sub.a(0), P(0|0)=P.sub.0 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.sub.u(k) the results of the estimations of the preceding step {circumflex over (x)}.sub.a(k1|k1), P(k1|k1) the parameters Q, R (covariance matrices) determination of the excitation force F.sub.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)}.sub.a(k|k), P(k|k), thus the complete state with its covariance matrix are estimated:

(129) $\{\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.sub.a(k) is the unknown state vector

(130) $x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],$ A.sub.a, B.sub.a, C.sub.a, D.sub.a are state realization matrices with P being the covariance matrix of the state vector and Q and R being calibration matrices, and extraction of the last component of the last component of the estimated state vector, namely the excitation force F.sub.ex applied by the waves incident upon the moving part, by use of an equation with the form: (k)=[0 1]{circumflex over (x)}.sub.a(k|k)

Third Embodiment

(131) According to the third embodiment of the invention, the excitation force applied by waves incident upon the moving part may be determined by use 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.

(132) For this embodiment, prior to this step, the following are therefore known: the measurements of the position z, and speed 2 of the moving part, the measurement of the force F.sub.u applied by the conversion machine (PTO) to the moving part, the model of the radiation force F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.0.sup.t h(t)()d, and the dynamics model of the wave-energy system {umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t).

(133) 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.

(134) This third embodiment overcomes this limitation. 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. It is now regarded as an unknown input of a (slightly different) global model representing the dynamics of the float.

(135) By combining the equation which describes the movement of the float according to the linear wave theory
M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t)

(136) and the expression for the radiation force
F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t)

(137) The following is obtained
(M+M.sub.){umlaut over (z)}(t)+F.sub.r(t)=w(t)+u.sub.d(t)

(138) where w(t)=F.sub.ex(t) and u.sub.d(t)=F.sub.hd(t)F.sub.u(t). u.sub.d(t) is regarded as a control input of the system described, known because F.sub.u(t) is measurable and F.sub.hd(t) may be calculated from a (linear or non-linear) model of the hydrostatic restoring force, while w(t) is an unknown input (disturbance).

(139) This equation may be expressed in state form, by defining

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

(141) which gives:

(142) $\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=-\frac{1}{M+M_{}}{F}_r\left(t\right)+\frac{1}{M+M_{}}\left(w\left(t\right)-u_d\left(t\right)\right)\end{array}$

(143) In the same way as for the second embodiment, using
F.sub.r(t)=.sub.0.sup.th(t)()d

(144) the radiation damping F.sub.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.sub.2(t)=(t)

(145) either

(146) $\{\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}$

(147) where x.sub.r is a (non-measurable) internal state not having any particular physical significance and (A.sub.r, B.sub.r, C.sub.r, and D.sub.r) are the state realization matrices

(148) By combining the two state-space representations, the following is obtained:

(149) $\{\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_{}}{x}_2\left(t\right)-\frac{C_r}{M+M_{}}{x}_r\left(t\right)+\frac{1}{M+M_{}}\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}$

(150) or, in an equivalent manner

(151) $\{\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_{}}& -\frac{C_r}{M+M_{}}\\ 0& B_r& A_r\end{array}\right],B_c=\left[\begin{array}{c}0\\ \frac{1}{M+M_{}}\\ 0\end{array}\right],C_c=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

(152) 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.sub.d(t). It is pointed out that u.sub.d(t)=F.sub.hd(t)F.sub.u(t) and that it is possible to use, for example, the piecewise affine approximation in order to calculate the hydrostatic restoring force:
F.sub.hd=F.sub.iz(t)G.sub.i for each interval .sub.i, to which z belongs.

(153) 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:

(154) $\{\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}$

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

(156) 0$\{\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)+\left(k\right)\end{array}$

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

(158) 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 is estimating, at each time instant kN, the state vector x(k) and the unknown input w(k)=F.sub.ex(k), on the basis of: a priori estimations: {circumflex over (x)}(kN|k1), (kN|k1) (the estimations of x(kN) and w(kN) respectively) using the measurements from the time k1; information about the outputs: y(kN), y(kN+1), . . . , y(k); and information about the known input: u.sub.d(kN), u.sub.d(kN+1), . . . , u.sub.d(k)

(159) Namely:
{circumflex over (x)}.sub.0={circumflex over (x)}(kN|k1)
.sub.0(kN|k1)
and:
x.sub.j={circumflex over (x)}(kN+j|k)
w.sub.j=(kN+j|k)
.sub.j={circumflex over ()}(kN+j|k)
for j=0, 1, . . . , N.

(160) The moving-horizon estimation problem is formulated as:

(161) $\underset{x_0,w_0,w_1,\mathrm{.Math.},{}_1,\mathrm{.Math.},{}_N}{\min }J\left(x_0,w_0,w_1,\mathrm{.Math.},{}_1,\mathrm{.Math.},{}_N\right)$

(162) with the following objective function (cost):

(163) $J=(.Math.)={.Math.x_0-{\hat{x}}_0.Math.}_{P_0^{-1}}^2+{.Math.w_0-{\hat{w}}_0.Math.}_{Q_0^{-1}}^2+\underset{j=0}{\overset{N-1}{.Math.}}{.Math.{}_j.Math.}_{Q^{-1}}^2++\underset{j=0}{\overset{N}{.Math.}}{.Math.y\left(k-N+j\right)-\hat{y}\left(k-N+j\right).Math.}_{R^{-1}}^2+\underset{j=0}{\overset{N-1}{.Math.}}{.Math.w_j-w_{j-1}.Math.}^2$

(164) where
(kN+j)=C.sub.dx.sub.j+D.sub.dw.sub.j+D.sub.du.sub.d(kN+j)

(165) The notation v.sub.M.sup.2 indicates the product v.sup.TMv for a vector of reals v and a symmetrical matrix of reals M.

(166) There are five terms in the cost function: The first term

(167) ${.Math.x_0-{\hat{x}}_0.Math.}_{P_0^{-1}}^2$
and the second term

(168) ${.Math.w_0-{\hat{w}}_0.Math.}_{Q_0^{-1}}^2$
added together form the arrival cost which summarizes information prior to the time kN and forms part of the data of the estimation problem. The third term and the fourth term limit the method and measurement noise, respectively, with the weighting matrices Q and R. The last term is a regularization term which has a function accounting for the fact that the excitation force of the waves is a regular (smooth) signal.

(169) 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 be imposing variation constraints thereon.

(170) 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:

(171) $\{\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)+\left(k\right)\end{array}$

(172) over time. This gives:

(173) $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}$$\mathrm{.Math.}$$x_n=A_d^N{x}_0+A_d^{N-1}{B}_d{w}_0+\mathrm{.Math.}+B_d{w}_{N-1}+A_d^{N-1}{B}_d{u}_d\left(k-N\right)+\mathrm{.Math.}+B_d{u}_d\left(k-1\right)+A_d^{N-1}{}_0+\mathrm{.Math.}+{}_1$

(174) or, in compact form,

(175) $X={}_x{x}_0+{}_{x}W+{}_x{U}_d+{}_{x}E$$\mathrm{where}$$X=\left[\begin{array}{c}{x}_1\\ x_2\\ \mathrm{.Math.}\\ A_N\end{array}\right],W=\left[\begin{array}{c}{w}_1\\ w_2\\ \mathrm{.Math.}\\ W_N\end{array}\right],U_d=\left[\begin{array}{c}{u}_d\left(k-N\right)\\ u_d\left(k-N+1\right)\\ \mathrm{.Math.}\\ u_d\left(k\right)\end{array}\right],E=\left[\begin{array}{c}{}_1\\ {}_2\\ \mathrm{.Math.}\\ {}_N\end{array}\right]$$\mathrm{and}$${}_x=\left[\begin{array}{c}{A}_d\\ A_d^2\\ \mathrm{.Math.}\\ A_d^N\end{array}\right],{}_x=\left[\begin{array}{ccccc}{B}_d& 0& \mathrm{.Math.}& 0& 0\\ A_d{B}_d& B_d& \mathrm{.Math.}& 0& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}& \mathrm{.Math.}\\ A_d^{N-1}{B}_d& A_d^{N-2}{B}_d& \mathrm{.Math.}& B_d& 0\end{array}\right],{}_x=\left[\begin{array}{cccc}I& 0& \mathrm{.Math.}& 0\\ A_d& I& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ A_d^{N-1}& A_d^{N-2}& \mathrm{.Math.}& I\end{array}\right]$

(176) Applying the same approach to the output equation and introducing the notations

(177) $Y={\left[{y\left(k-N+1\right)}^T{y\left(k-N+2\right)}^T{y\left(k\right)}^T\right]}^T,{}_y=\left[\begin{array}{c}{C}_d{A}_d\\ C_d{A}_d^2\\ \mathrm{.Math.}\\ C_d{A}_d^N\end{array}\right],{}_y=\left[\begin{array}{ccccc}{C}_d{B}_d& D_d& \mathrm{.Math.}& 0& 0\\ C_d{A}_d{B}_d& C_d{B}_d& \mathrm{.Math.}& 0& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}& \mathrm{.Math.}\\ C_d{A}_d^{N-1}{B}_d& C_d{A}_d^{N-2}{B}_d& \mathrm{.Math.}& C_d{B}_d& D_d\end{array}\right],{}_y=\left[\begin{array}{cccc}{C}_d& 0& \mathrm{.Math.}& 0\\ C_d{A}_d& C_d& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ C_d{A}_d^{N-1}& C_d{A}_d^{N-2}& \mathrm{.Math.}& C_d\end{array}\right]$

(178) results in a QP problem with the form:

(179) $\{\begin{array}{c}\underset{}{\min }\left\{{}^T\left(P+{}_s^{T}R{}_s\right)-2\left[{\hat{w}}_0^T{\hat{x}}_0^T{V}^T\right]\left[\begin{array}{c}S\\ R{}_s\end{array}\right]\right\}\\ \mathrm{such}\mathrm{as}{H}_{}{g}_{}+H_d{U}_d\end{array}\mathrm{where}\{\begin{array}{c}=\left[W^T{x}_0^T{E}^T\right]\\ V=Y-{}_y{U}_d\\ {}_s=\left[{}_y{}_y{}_y\right]\end{array}$

(180) where is a matrix which represents all the unknown variables to be estimated (the values of the unknown input W.sup.T and the uncertainty about the state E.sup.T for all the N steps of the estimation window considered, and the initial state x.sub.0.sup.T at the time kN), .sub.s is a matrix which includes the effects of these same unknown variables on the values output in the estimation window Y, .sub.y is a matrix which corresponds to the effects of the known input U.sub.d on Y, .sub.y .sub.y and .sub.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.sub.0.sup.1, R.sub.0.sup.1 and the regularization parameter as follows:

(181) 0$Q=\left[\begin{array}{cccc}{Q}^{-1}& 0& \mathrm{.Math.}& 0\\ 0& Q^{-1}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ 0& 0& \mathrm{.Math.}& Q^{-1}\end{array}\right],R=\left[\begin{array}{cccc}{R}^{-1}& 0& \mathrm{.Math.}& 0\\ 0& R^{-1}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ 0& 0& \mathrm{.Math.}& R^{-1}\end{array}\right],{}_1=\left[\begin{array}{ccccccc}{Q}_0^{-1}-& & 0& 0& \mathrm{.Math.}& 0& 0\\ 0& -& & 0& \mathrm{.Math.}& 0& 0\\ 0& 0& -& & \mathrm{.Math.}& 0& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}& \mathrm{.Math.}\\ 0& 0& 0& 0& \mathrm{.Math.}& & 0\\ 0& 0& 0& 0& \mathrm{.Math.}& -& \end{array}\right],P=\left[\begin{array}{ccc}{}_1& 0& 0\\ 0& P_0^{-1}& 0\\ 0& 0& Q\end{array}\right],S=\left[\begin{array}{cccc}{Q}_0^{-1}& 0& \mathrm{.Math.}& 0\\ 0& P_0^{-1}& \mathrm{.Math.}& 0\end{array}\right]$

(182) To conclude
H.sub.g.sub.+H.sub.dU.sub.d

(183) is a linear inequality which, with the aid of the matrices H.sub. and H.sub.d and the vector g.sub., retranscribes in compact form polytopic constraints on the unknown input and state:

(184) $\{\begin{array}{c}w{P}_w,P_w=\left\{w\text{:}{H}_w{g}_w\right\},g_w>0\\ x{P}_x,P_x=\left\{x\text{:}{H}_x{g}_x\right\},g_x>0\end{array}$

(185) which in turn are defined by the matrices H.sub.w and H.sub.x and the vectors g.sub.w and g.sub.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.

(186) 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.

(187) 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:
(kN+j|k)=*(j),j=1,2, . . . ,N

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

(189) The components of the vectors of the past estimations and measurements corresponding to the N1 time step preceding the present instant k=0 are initialized to zero.

(190) 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.sub.u(k), the results of the estimations of the preceding time step {circumflex over (x)}.sub.d(kN|k1), (kN|k1) the parameters Q, R, Q.sub.0.sup.1, R.sub.0.sup.1 (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 *

(191) $\{\begin{array}{c}\underset{}{\min }\left\{{}^T\left(P+{}_s^{T}R{}_s\right)-2\left[{\hat{w}}_0^T{\hat{x}}_0^T{V}^T\right]\left[\begin{array}{c}S\\ R{}_s\end{array}\right]\right\}\\ \mathrm{such}\mathrm{as}{H}_{}{g}_{}+H_d{U}_d\end{array}$ ii. The excitation force applied by waves incident upon the moving part w(k) is estimated by use of an equation with the form:
(k|k)=*(N).

(192) 6) Controlling of the Wave-Energy System

(193) During this step, the wave-energy system is controlled to account for the excitation force applied by the incident wave. Thus it is possible to drive the wave-energy system to optimize the recovered energy.

(194) Controlling may be controlling the moving part of the wave-energy system, for example by use of an electric, pneumatic or hydraulic PTO conversion machine. The 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

(195) 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 allow 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 electrical machine is activated.

(196) 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.

(197) 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 should be 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.

(198) 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.