FR 3049989 A1 20171013 FR 1653109 A 20160408 1 L'invention concerne le domaine des dispositifs pour convertir l'énergie des vagues en énergie électrique. En particulier, l'invention concerne le domaine de la détermination de la force d'excitation par la houle incidente sur un moyen mobile, notamment pour la commande d'un système houlomoteur. Les ressources énergétiques renouvelables connaissent un fort intérêt depuis quelques années. Ces ressources sont propres, gratuites et inépuisables, autant d'atouts majeurs dans un monde rattrapé par la diminution inexorable des ressources fossiles disponibles et prenant conscience de la nécessité de préserver la planète. Parmi ces ressources, l'énergie houlomotrice, source relativement méconnue au milieu de celles largement médiatisées comme l'éolien ou le solaire, contribue à la diversification indispensable de l'exploitation des énergies renouvelables. Les dispositifs, couramment appelés appareils "houlomoteurs", sont particulièrement intéressants, car ils permettent de produire de l'électricité à partir de cette source d'énergie renouvelable (l'énergie potentielle et cinétique des vagues) sans émission de gaz à effet de serre. Ils sont particulièrement bien adaptés pour fournir de l'électricité à des sites insulaires isolés. Par exemple, les demandes de brevet FR 2876751, et WO 2009/081042 décrivent des appareils pour capter l'énergie produite par le flot marin. Ces dispositifs sont composés d'un support flottant sur lequel est disposé un pendule monté mobile par rapport au support flottant. Le mouvement relatif du pendule par rapport au support flottant est utilisé pour produire de l'énergie électrique au moyen d'une machine de conversion de l'énergie (par exemple une machine électrique). La machine de conversion fonctionne en tant que générateur et en tant que moteur. En effet, pour fournir un couple ou une force qui entraîne le mobile, on fournit une puissance à la machine de conversion afin de mettre le mobile en résonance avec les vagues (mode moteur). Par contre, pour produire un couple ou force qui résiste au mouvement du mobile, on récupère une puissance via la machine de conversion (mode générateur). Dans l'intention d'améliorer le rendement et donc la rentabilité des dispositifs pour convertir l'énergie des vagues en énergie électrique (systèmes houlomoteurs), il est intéressant d'estimer la force d'excitation exercée par la houle incidente sur le système houlomoteur. La force appliquée par la houle aux dispositifs pour convertir l'énergie des vagues en énergie électrique ou d'autres formes d'énergie exploitable, couramment appelés houlomoteurs, n'est pas directement mesurable en conditions de fonctionnement normales dans un contexte industriel. Pour les convertisseurs ponctuels oscillants, il est possible de réaliser des essais spécifiques, consistant à bloquer la partie flottante de la machine et à 3049989 2 mesurer l'effort nécessaire pour contrer l'action de la vague (c'est-à-dire pour maintenir le flotteur immobile) à travers un ou plusieurs capteurs (de force ou de couple) placés sur la machine de conversion, appelée aussi système de prise de puissance (PTO), par exemple un générateur électrique assorti d'un dispositif permettant de transmettre le mouvement 5 oscillatoire du moyen mobile. Selon la terminologie courante, on appelle la quantité ainsi mesurée, force d'excitation de la houle (ou vague) incidente, en la distinguant notamment de la force (et de la houle) générée par le mouvement même du flotteur (force de radiation). Durant le fonctionnement normal du flotteur, en revanche, les mêmes capteurs d'effort mesureront seulement la force appliquée par le PTO, et non pas la force d'excitation de la 10 houle incidente. Il serait possible, en principe, de calculer l'ensemble des forces appliquées par la houle en intégrant des mesures de pression provenant de capteurs distribués sur toute la surface, mais il s'agit là d'une solution coûteuse et peu robuste qui peut difficilement être envisagée dans un contexte industriel. 15 Il y a peu de travaux scientifiques dans le domaine de l'estimation en temps-réel de la force (ou du couple) d'excitation de la vague. Par exemple, on peut citer le document Peter Kracht, Sebastian Perez-Becker, Jean-Baptiste Richard, and Boris Fischer. "Performance Improvement of a Point Absorber Wave Energy Converter by Application of an ObserverBased Control: Results from Wave Tank Testing". Dans IEEE Transactions on Industry 20 Applications, juillet 2015, qui décrit un algorithme d'estimation de la force de la houle, celui-ci est basé sur banc d'oscillateurs harmoniques indépendants et un observateur de Luenberger. Les résultats obtenus avec cette méthode montrent que les estimations sont retardées (déphasées) de manière significative par rapport à la vraie force d'excitation, ce qui est gênant notamment pour une utilisation dans le cadre du contrôle du système 25 houlomoteur. Plus généralement, cette approche considère une plage de fréquence de fréquence fixe (choisie a priori) pour le spectre de la houle. Pour que la méthode puisse marcher dans des conditions réalistes, où le spectre est variant dans le temps, il faut prendre en compte un nombre très élevé de fréquences, ce qui rend l'approche très lourde en termes de calculs. Dans le document Bradley A Ling. "Real-Time Estimation and Prediction of Wave 30 Excitation Forces for Wave Energy Control Applications",publié dans ASME 2015 34th International Conference on Ocean, Offshore and Arctic Engineering, une approche basée sur un filtre de Kalman étendu est proposée pour reconstruire le spectre en considérant la force d'excitation de la vague comme une (seule) sinusoïde dont les paramètres (amplitude, fréquence et phase) varient dans le temps. Toutefois, la méthode ne peut être efficace que 35 pour des vagues dans une bande de fréquences très étroite. 3049989 3 En outre, la demande de brevet FR2973448 (US9261070) décrit une méthode de commande pour convertisseurs ponctuels oscillants dont une des étapes est l'estimation (du spectre) de la force de la houle sur le flotteur (ou du couple sur la partie mobile du convertisseur), à partir d'un ensemble de sinusoïdes et d'un observateur de Luenberger. La 5 méthode s'apparente à celle publiée dans les documents précédemment cités, et présente a priori les mêmes désavantages. L'invention propose d'améliorer le fonctionnement d'un appareil houlomoteur en déterminant en temps réel les efforts exercés par la houle incidente sur le moyen mobile, afin 10 d'adopter en conséquence les meilleurs réglages des stratégies de récupération de l'énergie électrique. Le procédé selon l'invention est basé sur la construction d'un modèle de la force de radiation exercée sur le moyen mobile, et d'un modèle de la dynamique du système houlomoteur. L'invention utilise seulement des mesures de la cinématique du flotteur (position, vitesse et éventuellement accélération) et de la force appliquée par la machine de 15 conversion, mesures normalement disponibles sur un système houlomoteur, car utilisées pour son contrôle et pour sa supervision. Ainsi, grâce aux modèles utilisés, il est possible d'estimer la force exercée par la houle incidente pour toute plage de fréquences des vagues, tout en gardant un temps de calcul adapté à une détermination et un contrôle en temps réel. 20 Le procédé selon l'invention L'invention concerne un procédé de détermination de la force d'excitation exercée par la houle incidente sur un moyen mobile d'un système houlomoteur, ledit système houlomoteur convertissant l'énergie de la houle en énergie au moyen dudit moyen mobile coopérant avec une machine de conversion, ledit moyen mobile effectuant un mouvement 25 par rapport à ladite machine de conversion sous l'action de la houle. Pour ce procédé, on réalise les étapes suivantes : a) on mesure la position et la vitesse dudit moyen mobile ; b) on mesure la force F, exercée par ladite machine de conversion sur ledit moyen mobile ; 30 c) on construit un modèle de la force Frad de radiation exercée sur ledit moyen mobile, ledit modèle de la force Frad de radiation reliant ladite force de radiation à la vitesse dudit moyen mobile ; d) on construit un modèle de dynamique dudit système houlomoteur qui relie ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile, à ladite 35 position dudit moyen mobile, à ladite vitesse du moyen mobile, à ladite force F, 3049989 4 exercée par ladite machine de conversion sur ledit moyen mobile et à ladite force de radiation Frad exercée sur ledit moyen mobile ; et e) on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile au moyen dudit modèle de dynamique, dudit modèle de la force de 5 radiation, desdites position et vitesse mesurées et de ladite force F, mesurée exercée par ladite machine de conversion sur ledit moyen mobile. Selon un mode de réalisation, on construit ledit modèle de dynamique dudit système houlomoteur au moyen d'une équation du type : M2(t) = Fex(t) + Fhd(t)+F'd(t) Fu(t) avec M la masse du moyen mobile, 2 l'accélération dudit moyen mobile, Fhd la force de 10 rappel hydrostatique exercée sur ledit moyen mobile. Selon une alternative, on détermine ladite force Fhd de rappel hydrostatique au moyen d'une formule du type : Fhd(t) = -Kz(t) avec z la position dudit moyen mobile et K le coefficient de raideur hydrostatique. 15 Alternativement, on détermine ladite force Fhd de rappel hydrostatique par une fonction affine par morceau du type : Fhd = -Fiz(t) - GL pour au moins deux intervalles Sil de variation de ladite position z, avec F, et G, des coefficients constants de ladite fonction affine pour chaque intervalle ni. Conformément à une mise en oeuvre, on construit ledit modèle de la force Frad de 20 radiation par une équation du type Frad(t) = - Fr(t) avec Fr(t) = fot h(t - z)2(z)dr l'amortissement de radiation, 2 l'accélération dudit moyen mobile, M' la masse ajoutée à la fréquence infiniment élevée, 2 la vitesse dudit moyen mobile, h la réponse impulsionnelle qui relie la vitesse du moyen mobile à l'amortissement de radiation. Avantageusement, on mesure l'accélération dudit moyen mobile. 25 Selon un premier mode de réalisation, on détermine l'amortissement de radiation Fr au moyen d'une réponse temporelle approximée issue d'une réalisation d'état, du type : Fr(t) = Cr f(t) e(t-r)Ar Bri(r) ctr + Dr2(t) avec Ar, Br, Cr, Dr des matrices de la réalisation d'état. De préférence, on détermine ladite force d'excitation Fex exercée par la houle incidente 30 sur ledit moyen mobile (2) au moyen d'une équation du type : Fex(k)=(M + Mc0)2(k) + Eit -01 Crd Ark d-1-1 rsrd 2(0 Drdi(k) FiZ(k) + G + Fu(k) avec k le pas de temps dans le domaine discret, Ard,Brd,Crd,Drd les matrices de la réalisation d'état dans le domaine discret. Selon un deuxième mode de réalisation, on détermine ladite force d'excitation Fex 35 exercée par la houle incidente sur ledit moyen mobile au moyen d'un observateur d'état basé 3049989 5 sur un filtre de Kalman linéaire construit à partir d'un modèle de marche aléatoire de ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile. De manière avantageuse, on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile en mettant en oeuvre les étapes suivantes : 5 i) on initialise k=0, le vecteur d'état 2,(010)=.Fc,(0) et l'état de la matrice de covariance, P(010) = P0, ; ii) à chaque instant k, on acquiert lesdites mesures de position et de vitesse dudit moyen mobile y(k) = [z(k) z(k)] et la mesure de ladite force exercée par ladite machine de conversion sur ledit moyen mobile u(k) = Fu(k) ; et 10 iii) à chaque instant k, on détermine la force d'excitation exercée par la houle incidente sur ledit moyen mobile (2) 0(k) au moyen des équations suivantes : tei(k)= [0 1] fca(klk) avec : fca(klk -1) = Aajea(k -1Ik - 1) + BaFa(k -1) P(klk -1) = A,P(k -11k -1)ATa+ Q , K(k) = P(kIk -1)CI(CaP(klk -1)Cj; + R)-1 .2a(klk)= 2a(klk - 1) + K(k)(Y(k) - Ccv2a(k I k - 1) - DaFu(k)) P(klk) = (I - K(k)Ca(k))P(klk - 1) avec K(k) le gain de l'observateur, xa(k) le vecteur d'état inconnu xa(k) = 15 Lw(k)J, Aa,Ba,Ca,Da des matrices de la réalisation d'état, P la matrice de covariance du vecteur d'état, Q et R des matrices de calibration. Selon un troisième mode de réalisation, on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile au moyen d'un estimateur d'état à entrée inconnue basé sur une approche d'horizon glissant. 20 Conformément à une caractéristique, on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile en mettant en oeuvre les étapes suivantes à chaque instant k : i) on acquiert lesdites mesures de position et de vitesse dudit moyen mobile y(k) = [z(k) z(k)] et la mesure de ladite force exercée par ladite machine de 25 conversion sur ledit moyen mobile u(k) = Fu(k) ; et ii) on détermine la force d'excitation exercée par la houle incidente sur ledit moyen mobile w(kik) au moyen des équations suivantes : 1-4)(k 1k) = «N) avec «N) obtenu par résolution d'un problème de programmation quadratique suivant : 30 i [ s [ ,1 ming) T (P + WS Pl; 2c3' vT' I_Rws1 '1 avec P, R et S des matrices - 2 tel que He ,g'.+ HdUd 3049989 6 formées à partir de matrices de pondération de la fonction coût, la matrice qui représente l'ensemble des valeurs des variables inconnues sur la fenêtre d'estimation de dimension N, IFS la matrice, calculée à partir de la représentation d'état, qui correspond à l'effet de ces variables sur les valeurs des sorties sur 5 cette même fenêtre, H et Hd des matrices et ,g'. un vecteur obtenus à partir de contraintes imposées sur l'état et l'entrée inconnue, Jeo = 2(k - NIk - 1) et 00 = 0(k - NIk - 1) avec x le vecteur d'état et w l'entrée ud(k - N) ,i(k - N + 1) inconnue de l'estimateur d'état, Ud = u- . [ , ud étant l'entrée connue lid(k) de l'estimateur d'état avec ud(t)=Fhd(t) - Fu(t). 10 Selon une mise en oeuvre, on commande ledit système houlomoteur en fonction de ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile, déterminée. Présentation succincte des figures 15 D'autres caractéristiques et avantages du procédé selon l'invention, apparaîtront à la lecture de la description ci-après d'exemples non limitatifs de réalisations, en se référant aux figures annexées et décrites ci-après. La figure 1 illustre un système houlomoteur selon un mode de réalisation de l'invention. La figure 2 illustre une approximation de la force de rappel hydrostatique en fonction de 20 la position selon un exemple. La figure 3 illustre l'amplitude de la houle pour la mise en oeuvre de l'exemple illustratif. Les figures 4a à 4c illustrent la force d'excitation exercée par la houle incidente sur le moyen mobile en fonction du temps, pour une houle de type W1 conformément à la figure 3, obtenue par mesure et respectivement par un des modes de réalisation du procédé selon 25 l'invention. Les figures 5a à 5c illustrent la force d'excitation exercée par la houle incidente sur le moyen mobile en fonction du temps, pour une houle de type W3 conformément à la figure 3, obtenue par mesure respectivement par une des modes de réalisation du procédé selon l'invention. 30 Description détaillée de l'invention La présente invention concerne la détermination de la force d'excitation exercée par la houle incidente sur un moyen mobile d'un système houlomoteur, appelée aussi force 3049989 7 d'excitation de la houle ou de la vague. Un système houlomoteur est un système qui convertit l'énergie de la houle en énergie récupérable, en particulier en énergie électrique. Un système houlomoteur comporte généralement un moyen mobile, également appelé pendule ou flotteur, qui a un mouvement oscillatoire sous l'action de la houle. Le moyen 5 mobile coopère avec une machine de conversion, appelée aussi système de prise de puissance (PTO de l'anglais « power take off »), qui comporte dans la plupart des cas un générateur électrique couplé à un dispositif permettant d'adapter la transmission du mouvement oscillatoire, pour convertir le mouvement du moyen mobile en énergie récupérable. Dans certains cas, la machine de conversion peut agir en tant que moteur en 10 générant un effort sur le moyen mobile. En effet, pour récupèrer une puissance via la machine de conversion, on produit un couple ou force qui résiste au mouvement du mobile (mode générateur). Par contre, si la machine de conversion le permet, on peut fournir une puissance à la machine de conversion pour fournir un couple ou une force qui entraîne le mobile afin de l'aider à se mettre en résonance avec les vagues (mode moteur). 15 Le procédé selon l'invention est adapté à tout type de système houlomoteur avec au moins un moyen mobile, par exemple ceux décrits dans la demande de brevet FR2973448 (US9261070). Le procédé de commande selon l'invention peut être également appliqué un système houlomoteur appartenant à la catégorie de systèmes houlomoteurs avec colonnes d'eau oscillantes (OWC de l'anglais Oscillating Water Column). 20 La figure 1 illustre, de manière schématique, un exemple non limitatif de système houlomoteur 1. Le système houlomoteur 1 comporte un moyen mobile 2, qui flotte (la houle est représentée de manière schématique par la surface hachurée). Le moyen mobile 2 coopère, au moyen d'un bras de levier 4, avec une machine de conversion basée sur un générateur électrique 3, qui est, dans ce cas, un générateur linéaire. Le bras de levier 4 et le 25 générateur électrique sont articulés par rapport à une installation 6. Les mouvements du bras de levier sont représentés schématiquement par des flèches. Dans la suite de la description et pour les revendications, les termes vagues, flots marins, et houle sont considérés comme équivalents. De plus, dans la description et pour les revendications, le terme force désigne un effort 30 ou un couple. De la même manière, les termes position, vitesse et accélérations désignent à la fois des valeurs « linéaires » ou « angulaires », les valeurs « linéaires » pouvant être associées aux efforts, et les valeurs « angulaires » pouvant être associées aux couples. En outre, pour une meilleure compréhension, les différents modèles sont représentés en une seule dimension, toutefois, le procédé selon l'invention est adapté à des modèles à 35 plusieurs dimensions, pour des systèmes ayant plusieurs degrés de liberté dans leur mouvement. 3049989 8 Le procédé selon l'invention comporte les étapes suivantes : 1) mesure de la position et de la vitesse du moyen mobile 2) mesure de la force exercée par la machine de conversion (ou PTO) 3) construction du modèle de force de radiation 5 4) construction du modèle de dynamique 5) détermination de la force d'excitation de la houle incidente 6) (étape facultative) commande du système houlomoteur 1) Mesure de la position et de la vitesse du moyen mobile 10 Lors de cette étape, la position et la vitesse du moyen mobile sont mesurées. La position correspond au mouvement (par exemple distance ou angle) par rapport à la position d'équilibre du moyen mobile. Ces mesures peuvent être réalisées au moyen de capteurs, généralement implémentés sur un système houlomoteur pour son contrôle et/ou sa supervision. 15 Selon une mise en oeuvre de l'invention, lors de cette étape, on peut également mesurer ou estimer l'accélération du moyen mobile, car cette mesure ou estimation peut être utilisée dans les modèles mis en oeuvre par le procédé selon l'invention. Par exemple, l'accélération peut être mesurée au moyen d'un accéléromètre placé sur le moyen mobile. 20 2) Mesure de la force exercée par la machine de conversion (PTO) Lors de cette étape, on mesure la force (l'effort ou le cas échéant le couple) exercée par la machine de conversion PTO sur le moyen mobile. Cette mesure peut être réalisée au moyen d'un capteur, qui peut être un capteur d'effort ou un capteur de couple. Ce type de capteur est souvent installé, ou peut être facilement installé dans les systèmes 25 houlomoteurs, pour leur contrôle et/ou leur supervision. Alternativement, la mesure peut être remplacée par une estimation réalisée à partir de la consigne de force (ou de couple) envoyée au PTO. Pour l'exemple de système houlomoteur illustré sur la figure 1, un capteur de couple 5 peut être disposé entre le bras de levier 4 et le générateur électrique 3. 30 3) Construction d'un modèle de force de radiation Lors de cette étape, on construit un modèle de la force de radiation exercée sur le moyen mobile. Selon la théorie linéaire des vagues (comme décrit par exemple dans le document Faines J, Kurniawan A. "Fundamental formulae for wave-energy conversion". R. 35 Soc. open sci. 2: 140305, 2005, http://dx.doi.org/10.1098/rsos.140305), la force de radiation résulte de l'oscillation d'un corps immergé (donc dépend du mouvement du moyen mobile), tandis que la force d'excitation, résultant de la présence même d'un corps dans l'eau, ne 3049989 9 dépend pas du mouvement du corps immergé, mais de la vague incidente. En l'absence de houle incidente, la force de radiation amortit l'oscillation résiduelle du corps immergé jusqu'à l'arrêter. Il est important de noter que, bien que la théorie linéaire permette de relier la force d'excitation à l'élévation de la vague incidente à travers un modèle linéaire (dans le domaine 5 fréquentiel ou temporel), en pratique on ne peut pas s'en servir pour calculer la force d'excitation en ligne, même si c'était possible de mesurer l'élévation de la vague au centre de gravité du flotteur comme le requiert la théorie. En effet, la relation linéaire entre élévation de la vague et force d'excitation est non-causale, ce qui veut dire qu'on ne peut pas calculer la force d'excitation à un instant donné sans connaître l'élévation de la vague dans les 10 instants futurs (le calcul peut en revanche être fait hors ligne, une fois la vague passée). Dans un contexte de contrôle en temps-réel, la force d'excitation ne peut pas donc être considérée que comme une force exogène totalement inconnue agissant sur le flotteur. En revanche, toujours selon la théorie linéaire des vagues, la force de radiation est reliée au mouvement du flotteur, et plus précisément à son accélération et à sa vitesse par un modèle 15 linéaire (dans le domaine fréquentiel ou temporel) causal. Elle peut donc être calculée en ligne en utilisant les mesures d'accélération et vitesse courantes (et passées pour la vitesse). Selon une mise en oeuvre de l'invention, on construit ledit modèle de la force Frad de radiation par une équation du type : 20 Frad(t) = -Mc02(t) Fr(t) avec Fr(t) = fot h(t - 1-)2(r)ch- la composante de ladite force Frad de radiation qui dépend de la vitesse (courante et passée) du moyen mobile, que l'on peut appeler amortissement de radiation , 2 l'accélération du moyen mobile, 25 Mo, la masse ajoutée à la fréquence infiniment élevée, que l'on peut obtenir par des codes de codes de calcul BEM (« Boundary Element Method » signifiant méthode des éléments de frontières), tels que WAMIT (WAMIT, USA), ou Nemoh (Ecole Centrale de Nantes, France), à partir de la géométrie du moyen mobile, 2 la vitesse du moyen mobile, 30 h la réponse impulsionnelle qui relie la vitesse du moyen mobile à l'amortissement de radiation, dont les coefficients s'obtiennent à partir des paramètres hydrodynamiques du moyen mobile calculés par les mêmes codes de calcul BEM. La construction de ce modèle permet une détermination à tout instant de la force de 35 radiation, avec un temps de calcul limité. Ainsi, la détermination de la force exercée par la houle peut être déterminée à tout instant avec un temps de calcul court. 3049989 10 4) Construction du modèle de dynamique du système houlomoteur Lors de cette étape, on construit un modèle de dynamique du système houlomoteur. On appelle modèle de dynamique un modèle qui relie la force d'excitation exercée par la 5 houle incidente sur le moyen mobile, la force de radiation exercée sur le moyen mobile, la force de rappel hydrostatique exercée sur le moyen mobile, la force exercée par la machine de conversion sur le moyen mobile, à la position et la vitesse du moyen mobile. Ce type de modèle permet d'obtenir des résultats représentatifs du comportement du système houlomoteur, si les mouvements ne sont pas trop grands. 10 Avantageusement, le modèle de dynamique est obtenu par application du principe fondamental de la dynamique (deuxième loi de Newton) au niveau du moyen mobile. Selon une mise en oeuvre de l'invention, pour lequel on considère les efforts, on peut construire le modèle de dynamique du système houlomoteur par une équation du type : 15 M2(t) = F'(t) + Fhd(t)+F'd(t) - Fu(t) avec M la masse du moyen mobile, 2 l'accélération du moyen mobile, F' la force d'excitation exercée par la houle incidente sur le moyen mobile, Frad la force de radiation exercée sur le moyen mobile, 20 Fhd la force de rappel hydrostatique exercée sur le moyen mobile, et Fu la force exercée par la machine de conversion sur le moyen mobile. Ce modèle traduit un mouvement de translation vertical (typique des flotteurs ayant un mouvement de pillonement). Ce modèle est dérivé de la théorie linéaire des vagues. Conformément à une première variante de réalisation, la force de rappel hydrostatique 25 exercée sur le moyen mobile peut être approximée par une fonction linéaire de la position z définie par rapport à la position d'équilibre. Dans ce cas, la force de rappel hydrostatique peut s'écrire par une fonction du type : Fhd(t) = -Kz(t) avec z la position dudit moyen mobile définie par rapport à sa position d'équilibre et K le coefficient de raideur hydrostatique. Ainsi, la force de rappel hydrostatique peut être calculée à partir d'un modèle 30 simple si l'on dispose de la mesure de la position z. Cette fonction est particulièrement adaptée pour les petits déplacements z. Conformément à une deuxième variante de réalisation, la force de rappel hydrostatique exercée sur le moyen mobile peut être approximée par une fonction affine par morceaux de la position z. La fonction peut être linéaire sur au moins deux, de préférence trois intervalles 35 de la position z. Dans ce cas, la force de rappel hydrostatique peut s'écrire par une fonction du type : Fhd = -Fiz(t)- Gi pour au moins deux intervalles Sil de variation de la position z, 3049989 11 avec F, et G, des coefficients constants de la fonction affine pour chaque intervalle ni. Cette fonction est particulièrement adaptée pour les mouvements z plus larges. La figure 2 illustre cette approximation. La figure 2 illustre en trait continu, de manière schématique et non limitative, les variations de la force de rappel hydrostatique Fhd en fonction de la position z. 5 Sur cette figure, sont représentés en traits pointillés les fonctions linéaires de chaque morceau 0.1, 0.2 et 0.3. On remarque que cette fonction affine par morceaux permet une meilleure approximation de l'effort de rappel hydrostatique sur l'ensemble de la plage de variation de la position z. 10 Selon une mise en oeuvre de l'invention, pour lequel on considère les couples, on peut construire le modèle de dynamique du système houlomoteur par une équation du type : JÔ(t) = Mex(t) + Mhyd(t)+Mrad(t) Mu(t) avec J l'inertie du moyen mobile, Ô l'accélération angulaire du moyen mobile, 15 M' le couple d'excitation exercé par la houle incidente sur le moyen mobile, Mrad le couple de radiation exercé sur le moyen mobile, Mhd le couple de rappel hydrostatique exercée sur le moyen mobile, et Mu le couple exercé par la machine de conversion sur le moyen mobile. Ce modèle traduit un mouvement de rotation autour d'un axe horizontal (typique des 20 flotteurs ayant un mouvement de tangage). Ce modèle est dérivé de la théorie linéaire des vagues. Ce modèle est le « miroir » du modèle pour lequel on considère les efforts : les différents termes du modèle sont de même nature. De la même manière que pour la force de rappel hydrostatique, le couple de rappel hydrostatique peut être approximé par une fonction linéaire, ou par une fonction affine par 25 morceaux. Dans la suite de la description, seul le cas relatif aux efforts sera illustré, mais le cas des couples peut se déduire par transposition des équations dans un référentiel angulaire. 30 5) Détermination de la force d'excitation exercée par la houle incidente Lors de cette étape, on détermine, en temps réel, la force d'excitation exercée par la houle incidente sur le moyen mobile au moyen : des mesures de position et de vitesse (et éventuellement d'accélération) du moyen mobile, 35 de la mesure de la force exercée par la machine de conversion PTO sur le moyen mobile, 3049989 12 - du modèle de la force de radiation, et - du modèle de dynamique du système houlomoteur. Selon un premier mode de réalisation de l'invention (décrit en détail dans la suite), on 5 peut déterminer la force d'excitation exercée par la houle incidente sur le moyen mobile sur la base du calcul d'une approximation de la force de radiation. Ce mode de réalisation présente l'avantage de permettre la prise en compte d'une représentation de la force de rappel hydrostatique par une fonction affine par morceaux. Selon un deuxième mode de réalisation de l'invention (décrit en détail dans la suite), 10 on peut déterminer la force d'excitation exercée par la houle incidente sur le moyen mobile au moyen d'un observateur d'état basé sur un filtre de Kalman linéaire construit à partir d'un modèle de marche aléatoire de la force d'excitation exercée par la houle incidente sur le moyen mobile. Ce deuxième mode de réalisation présente l'avantage de prendre en compte en compte des incertitudes. 15 Selon un troisième mode de réalisation de l'invention (décrit en détail dans la suite), on peut déterminer la force d'excitation exercée par la houle incidente sur le moyen mobile au moyen d'un estimateur d'état d'entrée inconnue basé sur une approche à horizon glissant. Ce mode de réalisation présente l'avantage de prendre en compte des incertitudes. De plus, il permet la prise en compte d'une représentation de la force de rappel hydrostatique par une 20 fonction affine par morceaux. D'autres modes de réalisation peuvent être envisagés. Premier mode de réalisation Selon le premier mode de réalisation de l'invention, on peut déterminer la force 25 d'excitation exercée par la houle incidente sur le moyen mobile sur la base du calcul d'une approximation de la force de radiation. Pour ce mode de réalisation, on mesure également l'accélération du moyen mobile. Préalablement à cette étape, on connaît donc : les mesures de position z, de vitesse 2 et d'accélération 2 du moyen mobile, 30 la mesure de la force Fu exercée par la machine de conversion PTO sur le moyen mobile, le modèle de la force de radiation F'd(t) = -M,02(t) - Fr(t) avec Fr(t) = f(t)h(t - 1-)2(r)dr, et le modèle de dynamique du système houlomoteur 35 M2(t) = F'(t) + Fhd(t)+F'd(t) - Fu(t) 3049989 13 Ainsi si l'on sait calculer Fr on peut alors obtenir F'd, et finalement F'(t), seule inconnue du modèle de la dynamique du système houlomoteur. Mais il est n'est pas aisé de réaliser le calcul du produit de convolution dans le domaine temporel dans l'expression précédente pour Fr(t), à cause du grand nombre 5 d'opérations mathématiques à réaliser et de données à stocker. Pour l'éviter, on peut la considérer comme un système linéaire avec Fr(t) en sortie et i(r) en entrée. L'équation résultante dans le domaine fréquentiel (de Laplace), obtenue en appliquant la méthode de Prony, est : Fr(s) = wr(s)2(s) où Wr(s) est une fonction de transfert. Il est possible de mettre l'équation précédente, 10 qui était dans le domaine de Laplace, sous une forme d'état équivalente, par exemple : 1.icr (t) = Arxr(t) + Bri(t) Fr(t) = Crxr(t) + Drz(t) où xr est un état interne n'ayant pas de signification physique particulière et (Ar,Br,Cr,Dr) sont les matrices de la réalisation d'état. On rappelle qu'en théorie des systèmes (et en automatique), une représentation d'état permet de modéliser un système dynamique sous une forme matricielle, en utilisant des variables d'état. Cette représentation 15 peut être linéaire ou non, continue ou discrète. La représentation permet de déterminer l'état interne et les sorties du système à n'importe quel instant futur si l'on connaît l'état à l'instant initial et le comportement des variables d'entrée qui influent sur le système. La forme d'état permet d'exprimer l'évolution dans le temps de Fr comme il suit : t Fr(t) = CretArxr(0) + Cr f ét-r)ArBri(r) dr + Dri(t) 0 mais cette expression ne peut pas être utilisée directement car la valeur initiale de 20 l'état interne xr(0) est inconnue. On peut toutefois approximer l'évolution temporelle Fr comme il suit : t Fr(t) = Cr f e(t-r)ArBri(r)dr + Dri(t) 0 puisque CretArxr(0) 0 pour t suffisamment large, étant donné que la matrice Ar est asymptotiquement stable pour ce type de systèmes. Ce qui correspond à négliger la contribution de l'état initial à l'évolution temporelle de Fr, sachant qu'elle va s'estomper très 25 rapidement par rapport à la contribution liée à la vitesse z. Le premier mode de réalisation consiste donc à calculer Fex à partir de z, 2, 2 et Fu (connus) à partir d'une équation de la forme : Pex(t)=(1 Mco)2(t) + Cr fot e t-2 ArBr - ) dr + Dri(t)+ Kz(t) + Fu(t) si l'on utilise un modèle linéaire pour la force de rappel hydrostatique, ou 3049989 14 P'(0=(i + m,0)2(t) + Cr fot e t-2 ArBr - ) dr + D r (t) Fi2(t) + G Fu(t) pour chaque intervalle ni auquel z appartient, si on utilise un modèle affine par morceaux pour la force de rappel hydrostatique. Puisque dans la pratique, les mesures de z, 2, 2 et Fu sont des signaux 5 échantillonnés, on peut développer l'algorithme dans le temps discret k, en discrétisant comme il suit l'équation d'état introduite précédemment pour représenter l'évolution de Fr : rxr (k + 1) = Ardxr(k) Brdi(k) tFr(k) = Crdxr(k) Drdi(k) Fr peut être ainsi donné par une équation de la forme : k-1 Fr(k) = CrdArdxr(0) Ard1-1Brdi(1) Drdi(k) 1=0 et son approximation par la formule: k-1 Fr (k) crd Ardi-i i(1) Drdi(k) 1=0 Ce qui donne (pour le cas où la force de rappel hydrostatique est considérée affine par 10 morceaux) l'expression suivante : Pex(k)=(M M')2(k) e=01 Crd Ard 1 1Brd2(1) Drdi (k) F i2(k) + G Fu(k) pour chaque intervalle ni auquel z(k) appartient. L'aspect original de ce mode de réalisation est l'approximation de la réponse temporelle de Fr, l'amortissement de radiation (la composante de la force de radiation 15 dépendante de la vitesse du flotteur), calculée en passant par une représentation d'état équivalente et en négligeant la contribution de l'état initial de cette représentation. Deuxième mode de réalisation Selon un deuxième mode de réalisation de l'invention, on peut déterminer la force 20 d'excitation exercée par la houle incidente sur le moyen mobile au moyen d'un observateur d'état basé sur un filtre de Kalman linéaire construit à partir d'un modèle de marche aléatoire de la force d'excitation exercée par la houle incidente sur le moyen mobile. On rappelle qu'un observateur d'état, ou un estimateur d'état, est, en automatique et en théorie des systèmes, une extension d'un modèle représenté sous forme de représentation d'état. 25 Lorsque l'état du système n'est pas mesurable, on construit un observateur qui permet de reconstruire l'état à partir d'un modèle. Pour ce mode de réalisation, préalablement à cette étape, on connaît donc : les mesures de position z, et de vitesse 2 du moyen mobile, la mesure de la force Fu exercée par la machine de conversion PTO sur le moyen 30 mobile, 3049989 15 le modèle de la force de radiation F'd(t) = - M ,,,2 (t) - Fr(t) avec Fr(t) = f(t) h(t - T)2(t)clt , et le modèle de dynamique du système houlomoteur M 2 (t) = F'(t) + Fhd(t)+F'd(t) - Fu(t) 5 Dans cette approche, le problème d'estimation de la force d'excitation de la houle est transformé en un problème classique d'estimation d'état (qui peut être résolu avec un filtre de Kalman linéaire), en exprimant la dynamique de la force d'excitation de la houle par un modèle de marche aléatoire. L'avantage principal de cette méthode est la prise en compte d'incertitudes permettant de prendre en compte les bruits de mesure et les dynamiques non 10 modélisées. En remplaçant dans l'équation qui décrit le mouvement du flotteur M 2 (t) = F'(t) + Fhd(t)+F'd(t) - Fu(t) les expressions pour la force de rappel hydrostatique (avec le modèle linéaire) et le modèle de la force de radiation Fhd(t) - -Kz(t) Frad(t) = -M'2(t) - Fr(t) K 1 1 x2 (t) = Mco Xi (t) M + Mco Fr (t) + M + M, (F'(t) - Fu (t)) M + Comme pour le premier mode de réalisation, la dynamique de l'amortissement de radiation Fr(t) peut être décrite par une représentation d'état équivalente au produit de convolution entre la réponse impulsionnelle de l'amortissement de radiation et la vitesse du 20 moyen mobile. t Fr(t) = 1 h (t - 1-)2 (T)clt o soit (Fr(t) = Arxr(t) + Br2(t) (Fr(t) = Crxr(t) + Dri (0 on obtient (M + M')2(t) + K z (t) + Fr(t) = F'(t) - Fu(t) 15 Cette équation peut être mise sous forme d'état, en définissant f xi(t) = z(t) tx2(t) = 2(t) ce qui donne {±i(t) = x2 (t) OU 3049989 16 (fîr(t) = Arxr(t) + Brx2(t) Fr (t)= Crxr(t)+ Drx2(t) où xr est un état interne (non mesurable) n'ayant pas de signification physique particulière, et (Ar, Br, Cr, Dr) sont les matrices de la réalisation d'état. En combinant les deux représentations d'état, on obtient : X1(0 = x2(0 { K Dr Cr 1 jC2 (t) _ Xi (t) X2 (t) Xr (t) + (Fex(t) - Fu(t)) M + MG° M + M,x, M + M,x, M + M,x, jcr(t) = Brx2(t) + Arxr(t) ou, de manière équivalente r±(t) = Acx(t) + Bc(Fex(t) - Fu(t)) 67(0 = Coc(t) 5 où x(t) = [xl(t) x2(t) xT(t)r Y = [xi(t)1 [x2 (t)] et 0 1 0 0 K Dr Cr 1=o 1 0 1 0 01 [ Ac -[ m + m Co M + M,x, M + M co' Bc [m + m co, Cc 0 Br Ar 0 On doit estimer Fex(t) en utilisant y(t) (le vecteur qui contient les mesures de la position et de la vitesse du moyen mobile) et Fu(t), la mesure de la force du PTO. Pour ce faire, on discrétise d'abord le système d'état ci-dessus (on rappelle que les 10 mesures sont échantillonnées), avec la méthode de Tustin. Ce qui donne rx(k + 1) = Adx(k) + B d(Fex(k) - Fu(k)) (y(k) = Cdx(k) + B d(Fex(k) - Fu(k)) On note qu'il est important d'effectuer la discrétisation avec la méthode de Tustin, car elle fait apparaître un terme direct dans l'équation de sortie qui permet de relier y(k) à Fex(k), et donc d'estimer Fex(k) en utilisant y(k). Cela ne serait pas possible avec d'autres méthodes de discrétisation, qui relieraient y(k) à Fex(k - 1) et non pas à Fex(k). 15 En définissant w(k) = Fex(k) , u(k) = Fu(k), il résulte le modèle global de la dynamique du flotteur : rx(k + 1) = Adx(k) + Bdw(k) - Bdu(k) t y(k) = Cdx(k) + D dw(k) - D du(k) Pour estimer la force d'excitation w(k), on la considère comme un état, dont on définit l'évolution en introduisant un modèle mathématique qui relie w(k) à w(k - 1). 20 Pour ce mode de réalisation, on utilise un modèle de marche aléatoire pour décrire l'évolution de la force d'excitation : 3049989 17 w(k + 1) = w(k) + Em(k) OU Em(k), qui décrit la variation de w(k) , est considéré comme un nombre aléatoire. En d'autres termes, ce modèle suppose qu'à chaque instant k, la force d'excitation s'éloigne d'un pas (d'une quantité) aléatoire de sa valeur précédente, et que ces pas sont distribués indépendamment et identiquement en taille. 5 On couple alors ce modèle de marche aléatoire de la force d'excitation à un modèle de la dynamique du flotteur : {y(k)+ 1) = Adx(k) + Bdw(k) - Bdu(k) + Ex(k) ty(k) = Cdx(k) + Ddw(k) - Ddu(k) + ,u(k) plus robuste vis-à-vis des erreurs de modélisation, grâce à l'introduction de Ex(k) qui représente les dynamiques non modélisées (frottements du PTO, non-linéarité 10 hydrostatique, ...) et 1,1(k) qui décrit les bruits de mesure sur la position et la vitesse du flotteur. En combinant le modèle de marche aléatoire de la force d'excitation de la houle incidente et la dynamique du flotteur, on obtient le système d'état augmenté : {x(k + 1) = Adx(k)+Bdw(k)- Bdu(k)+ Ex(k) w(k + 1) = w(k) + Em(k) y(k) = Cdx(k) + Ddw(k) - Ddu(k) + ,u(k) 15 ou, de manière équivalente, {rxd (k + 1) = Aaxa(k) + Bau(k) + E(k) y(k) = Caxa(k) + Dau(k) + ,u(k) OÙ x(k) Ex(k) xa(k) - [w(k)1, e(k) - [Em(k)] et Aa 1 = [A 0 d Bd] a =[ Bal Ca = [Cd DA Da = -Dd 0 ' Ainsi, le problème d'estimation de w(k), devient un problème d'estimation d'état. 20 Une manière d'estimer le vecteur d'état inconnu xa(k), en prenant en compte des informations sur E(k) et ,u(k), est d'appliquer l'algorithme du filtre de Kalman linéaire, comme montré par la suite. Une alternative peut être d'utiliser l'approche d'estimation à horizon glissant (la même utilisée dans le troisième mode de réalisation pour estimer la force d'excitation, vue cette fois-ci comme une entrée inconnue). 25 L'algorithme du filtre de Kalman linéaire (KF) fournit la solution du problème d'optimisation suivant : { min (xa(0) - .,a(0))TPci1(xa(0) - .,a(0)) + 1 EU - 1)TQ-1EU - 1) + PU )7' R-1 PU) xa(k) j = 1 k 3049989 18 où Po, Q, et R sont des matrices de calibration de dimensions appropriées, .ca(0) est la valeur moyenne de l'état initial xa(0). En considérant xa(0) comme un vecteur aléatoire, non corrélé à E(k) et 1,1(k), de moyenne .,a(0) avec matrice de covariance P0, et E(k) et 1,1(k) comme deux processus 5 aléatoires de type bruit blanc avec moyenne nulle avec matrices de covariance Q et R respectivement, on peut résoudre le problème d'optimisation avec l'algorithme du filtre de Kalman linéaire. À chaque instant k, l'algorithme du filtre de Kalman linéaire calcule la solution à ce problème en passant par deux étapes. 10 La première étape est la mise à jour temporelle des estimations : t)?a (k 1 k - 1) = Au% (k - 11k - 1) + BaFu(k - 1) P(klk - 1) = AaP(k - 11k - 1)ATa + Q où .2a(k lk - 1) et P(kIk - 1) sont respectivement l'estimation de xa(k) et sa matrice de covariance obtenues en utilisant les mesures à partir de l'instant k - 1 (soit y(k - 1), y(k 2), ... et Fa(k - 1), Fa(k - 2) , ... ), et .2a(k1k) et P (kik) sont l'estimation de xa(k) et sa matrice de covariance obtenues en utilisant les mesures à partir de l'instant k (soit y(k),y(k - 1), ... 15 et Fu(k), Fu(k - 1) , ... ). La seconde étape est la mise à jour des mesures : K(k) = P(klk - 1)Cj; (CaP(kik - 1)Cj; + R)-1 { .2a(k1k) = 2a(klk - 1) + K(k)(y(k) - Ca.5ea(klk - 1) - DaMa(k)) P(k1k) = (I - K(k)Ca(k))P(kik - 1) Avec K(k) le gain de l'observateur, qui est un vecteur de paramètres qui pondère les erreurs entre les estimations et les mesures des sorties du système (position et vitesse). 20 Dans le cas du filtre de Kalman, ce n'est pas un paramètre de calibration, mais il est issu de l'estimation de la matrice de covariance. En résumant, pour le deuxième mode de réalisation, on peut mettre en oeuvre les étapes suivantes : - On initialise k = 0, .2a(010)=.%(0), P(010) = Po 25 - À chaque instant k : - on utilise : - les mesures de la position et de la vitesse du moyen mobile y(k) = [z(k) 2(k)] et de la force exercée par le PTO sur le moyen mobile u(k) = Fa(k) 30 - les résultats des estimations du pas précédent 2a(k - 11k - 1), P(k - 11k - 1) - les paramètres Q, R (matrices de covariance) 3049989 19 - on détermine la force d'excitation Fex exercée par la houle incidente sur le moyen mobile, noté pour ce mode de réalisation 0 (k) , en mettant en oeuvre les étapes suivantes : - on applique les deux étapes de l'algorithme du filtre de Kalman linéaire 5 pour obtenir 2a(klk), P (kik), ainsi on estime l'état complet avec sa matrice de covariance : rPa(k lk - 1) = Aaîca(k - 11k - 1) + BaFu(k - 1) o (kik - 1) = A aP (k - 11k - 1) ATa + Q , K (k) = P (kik - 1)Cj; (C aP (kik - 1)Cj; + R)-1 o )? a(k1k) = .5 e a(kIk - 1) + K (k)(Y (k) - C a)? a(k ik - 1) - D aFu(k)) P (kik) = (I - K (k)C a(k))P (kik - 1) (k avec xa(k) le vecteur d'état inconnu x a(k) = [wk))], Aa, Ba, Ca, Da des 10 matrices de la réalisation d'état, P la matrice de covariance du vecteur d'état, Q et R des matrices de calibration. - et on extrait la dernière composante du vecteur d'état estimé, c'est-à-dire la force d'excitation Fex exercée par la houle incidente sur le moyen mobile par une équation de la forme : 0(k) = [0 1] 2a(klk) 15 Troisième mode de réalisation Selon le troisième mode de réalisation de l'invention, on peut déterminer la force d'excitation exercée par la houle incidente sur le moyen mobile au moyen d'un estimateur d'état à entrée inconnue basé sur une approche d'horizon glissant. L'estimation à horizon 20 glissant est un algorithme qui détermine récursivement l'état et l'entrée inconnue en utilisant des données appartenant à une fenêtre temporelle finie. Pour ce mode de réalisation, préalablement à cette étape, on connaît donc : - les mesures de position z, de vitesse 2 du moyen mobile, - la mesure de la force Fu exercée par la machine de conversion PTO sur le moyen 25 mobile, - le modèle de la force de radiation Frac/ (t) = - M ,02 (0 - Fr (t) avec Fr (t) = f(t) h(t - 1-)2(r)dr , et - le modèle de dynamique du système houlomoteur 2(t) = Fex (t) + Fhd(O+Frad(0 - Fu (t) 30 Le deuxième mode de réalisation est complètement linéaire, et ne permet pas de prendre en compte directement une description non linéaire (plus précise) de la force de rappel hydrostatique. 3049989 20 Ce troisième mode de réalisation permet de s'affranchir de cette limitation. À la différence du deuxième mode de réalisation, la force d'excitation n'est plus considérée comme un état du modèle global qui représente la dynamique du flotteur, dont l'évolution est décrite par un modèle de marche aléatoire, mais est considérée comme une entrée inconnue 5 d'un modèle global (légèrement différent) représentant la dynamique du flotteur. En combinant l'équation qui décrit le mouvement du flotteur selon la théorie linéaire des vagues M2(t) = F'(t) + Fhd(t)+Frad(t) - Fu(t) et l'expression pour la force de radiation Frad(t) = -M'2(t) - Fr(t) on obtient (M + M,0)2(t) + Fr(t) = w(t) + ud(t) 10 où w(t) = Fex(t) et ud(t)=Fhd(t) - Fu(t). On considère ud(t) comme une entrée de contrôle du système décrit, connue car Fu(t) est mesurable et Fhd(t) peut se calculer à partir d'un modèle (linéaire ou non linéaire) de la force de rappel hydrostatique, tandis que w(t) est une entrée (une perturbation) inconnue. Cette équation peut être mise sous forme d'état, en définissant {x2(t) = z(t) tx2(t) = 2(t) 15 ce qui donne Comme pour le deuxième mode de réalisation, en partant de t Fr(t) = 1 h(t - 1-)2(t)clt 0 on peut exprimer l'amortissement de radiation Fr(t) comme la sortie d'un système linéaire sous forme de représentation d'état dont l'entrée est x2(t) = 2(t) soit (Fr(t) = Arxr(t) + Bri(t) (Fr(t) = Crxr(t) + Dri(t) 20 soit tîr(t) = Arxr (0 + Brx2(t) Fr(t) = Crxr(t) + Drx2(t) {±1(t) = x2(t) 1 1 jc2 (t) = M -F M,0 Fr(t) + M -F M (w(t) - ud(t)) ,0 3049989 21 5 où xr est un état interne (non mesurable) n'ayant pas de signification physique particulière et (Ar, Br, Cr, Dr) sont les matrices de la réalisation d'état. En combinant les deux représentations d'état, on obtient : {Xi (t) = x2(t) D, c° Cr 1 X2 (0 - X2 (0 Xr + (W(0 + (0) M M M c° M M Mc° 4(0 = BrX2 + ArXr ou, de manière équivalente r±(t) = Acx(t) + Bc(w(t) + ud(t)) 67(0 = Ccx(t) où x(t) = [x1(t) x2(0 x7:(01T [xi (01 [x2 (01 Y= et 0 1 0 Dr Cr Ac -[(:) m + m co M + M co' 0 Br Ar = [1 0 01 10 1 01 0 1 On cherche à estimer w(t) en utilisant y(t) le vecteur qui contient les mesures de la 10 position et de la vitesse du flotteur) et ud(t). On rappelle que ud(t)=Fhd(t) - Fu(t) et qu'on peut utiliser, par exemple, l'approximation affine par morceaux, pour calculer la force de rappel hydrostatique : Fhd = -Fiz(t) - GL pour chaque intervalle Sil, auquel z appartient. Pour ce faire, on peut d'abord discrétiser le système d'état ci-dessus (on rappelle que 15 les mesures sont échantillonnées), avec la méthode de Tustin. Ce qui donne rx(k + 1) = Adx(k) + Bdw(k) + Bdud(k) ty(k) = Cdx(k) + D dw(k) + D dud(k) Puis rendre le modèle plus réaliste en rajoutant des incertitudes (des bruits) sur l'état et la sortie : x(k + 1) = Adx(k) + B dw (k) + Bdud(k) + E(k) y(k) = Cdx(k) + D dw(k)+D dud (k) +µ(k) où E(k) et 1,1(k) sont des bruits blancs de moyenne nulle. Pour estimer l'entrée inconnue w(k) (et accessoirement l'état x(k)) on utilise la 20 méthode d'estimation à horizon glissant (en anglais « moving-horizon estimation », MHE). Le problème consiste à estimer, à chaque instant de temps k N le vecteur d'état x(k) et l'entrée inconnue w(k) = Fex(k), sur la base : 3049989 22 - d'estimations a priori : 2(k - Nik - 1), 0(k - Nik - 1) (respectivement les estimations de x(k - N) et w(k - N)) en utilisant les mesures à partir du temps k - 1 ; - d'informations sur les sorties : y(k - N), y(k - N + 1), ...,y(k) ; et 5 - d'informations sur l'entrée connue : ud(k - N), ud(k - N + 1), ...,ud(k) Soient : 20 = 2(k - Nik -1) = 0(k - Nik - 1) et : xj = 2(k - N + 11k) wj = f4)(k - N + j1k) = e(ic - N + j1k) 10 pour j = 0,1, ...,N. Le problème d'estimation à horizon glissant est formulé comme : min /(xo, wo, w1, --- --- 'EN) xo,www1,----E1-----EN avec la fonction objectif (coût) suivante N-1 J(.) = 11x0 - + iiwo -14,011Q,71 + II 11Q' + j=0 N-1 +- N + - y'(k - N + j)I1R2 -1 +111141i 141]-1112 j=0 j=0 où y'(k - N + j) = Cdxj + Ddwj + Ddud(k - N + j) 15 La notation MOL indique le produit vTMv pour un vecteur de réels U et une matrice symétrique de réels M. Il y a cinq termes dans la fonction coût : - Le premier terme 11x0 - 20 11p2,71 et le deuxième terme 11w0 - i o IIQ2,7i sommés ensemble forment le « coût d'arrivée », qui résume l'information préalable au 20 temps k - N et fait partie des données du problème d'estimation. - Le troisième terme et le quatrième terme pénalisent les bruits de procédé et de mesure respectivement avec les matrices de pondération Q et R. - Le dernier terme est un terme de régularisation qui sert à prendre en compte le fait que la force d'excitation de la houle est un signal régulier (lisse). 25 Ce dernier terme est important pour déterminer la force exercée par la houle de manière optimale. Une autre manière de prendre en compte la nature « lisse » de la force d'excitation de la houle consisterait à lui imposer des contraintes de variation. 3049989 23 En utilisant la démarche classique de l'estimation et de la commande à horizon glissant, la minimisation de la fonction coût peut être transformée en un problème de programmation quadratique QP (on rappelle qu'un problème de programmation quadratique est un problème d'optimisation dans lequel on minimise/maximise une fonction quadratique 5 sur un polyèdre convexe) en propageant : pc(k + 1) = Adx(k) + Bdw(k)+Bdud(k) + E(k) y(k) = Cdx(k) + Ddw(k)+Ddud(k) +µ(k) dans le temps. Ce qui donne : x1 = Adxo + Bdwo + Bdud(k - N) + 6.0 x2 = Adxi + Bdwi + Bdud(k - N + 1) + = Aâxo + AdBdwo + Bdwi + AdBdud(k - N) + Bdud(k - N + 1) + Ade() + xn = Ad X0 +Arw + --- + 1-B d 0 - B dw N -1 + Ad 1Bdud- N) + - - - +Bdud(k - 1) + 6.0 + - - - + ou, en forme compacte, X = (I)xxo + xW + xUd rxE où X = xi x.2 ,w = wi. w.2 ,u d = [ud(k - N) ,E = [ XN WN ud(k - N + 1) . ud(k) 1 6.2 EN 10 et Ad Ad ,tu _ [Bd AdB d 4_, = 0 OC) o. "I'x [ A1,11 -.-- x Ar 1 B d 0 0 Bd Bd Air 2 B d _ = 'lx I 0 0 Ad I 01 AN-1 /-1 d AN- 2 rl / En appliquant la même démarche à l'équation de sortie, et en introduisant les notations y (k)T1T , (1) = Y= [ y(k ccil - N Y CdAâ m + 1)T y(k - N + 2)T - ' 'r Y - rBd Dd [ C d Al c'ii C dAdB d CdBd C d Air 1 Bd Cd Air 2 B d o o o o Y C dB d Dd Cd 0 0 CdAd Cd 0 C d Ar 1 C d Ar 2 Cd on aboutit à un problème QP de la forme : {min (P + SR s)- 2 [01; -T VT] [RStps] X0 tel que He + HdUd où 3049989 24 {" = [WT 4 ET] V=Y-W U d Ws = [Wy (1)y FA avec une matrice qui représente l'ensemble des variables inconnues à estimer (les valeurs de l'entrée inconnue WT et de l'incertitude sur l'état ET pour tous les N pas de la fenêtre d'estimation considérée, et l'état initial 4 au temps k - N), tifs une matrice qui englobe les effets de ces mêmes variables inconnues sur les valeurs sorties sur la fenêtre 5 d'estimation Y, Wy une matrice qui correspond aux effets de l'entrée connue Ud sur Y, Wy Wy et tifs se calculant à partir de la représentation d'état définie auparavant, et P, R et S sont des matrices formées à partir des matrices de pondération de la fonction coût Q,R,Q(71-,Rdi et du paramètre de régularisation À comme il suit : Q= Q-1 0 0 Q-1R= R-1 0 0 [ o Q-1 o o R-1 --- 0 o o 0 l R-1 -Qd1-À À 0 0 ... 0 0 0 -À À 0 --- 0 0 0 0 Al = 0 0 -À À --- 0 0 ,P=[Al 0 /V 0,S Q(71 0 = [ 0 0 Q 0 Pd ' 0 0 0 0 --- À 0 0 0 0 0 --- -À À- Pour terminer .1--- 0 t x E Px,Px = tx : Hx gxl, gx> 0 à leur tour définies par les matrices Hw et Hx et les vecteurs gw et gx.Cette inégalité permet de prendre en compte des informations sur les amplitudes minimales et maximales de la force d'excitation de la houle incidente, de la position et de la vitesse du flotteur. 15 Le problème d'estimation à horizon glissant est un problème de programmation quadratique sous contraintes que l'on peut résoudre avec des solveurs tels que FORCES 0 (Embotech, Suisse) ou CVXGEN 0 (CVXGEN LLC, USA) par exemple. Soit Ç la solution du problème de programmation dynamique obtenue avec l'un de ces solveurs. L'estimation des forces d'excitation, à partir des mesures au temps k, est donnée 20 par : 0(k - N + j1k) = C(j),Vj = 1,2, ...,N En résumant, pour le troisième mode de réalisation, on peut mettre en oeuvre les étapes décrites ci-dessous. 0°1 3049989 25 On initialise à zéro les composants des vecteurs des mesures et estimations passées correspondants au N-1 pas de temps précédant l'instant actuel k = 0 À chaque instant k, 1. On acquiert : 5 - les mesures de la position et la vitesse du moyen mobile y(k) = [z(k) 2(k)] et de la force exercée par le PTO sur le moyen mobile u(k) = Fu(k), - les résultats des estimations du pas de temps précédent 2d(k - NIk - 1), 0(k - A 1 lk - 1) 10 - les paramètres Q, R, Q(71- , R(71- (matrices de pondération) et À (régularisation) 2. On calcule 0(1clic) qui correspond à l'estimation de la force exercée par la houle sur le moyen mobile, en mettant en oeuvre les étapes suivantes : i. on résout le problème de programmation quadratique OP suivant pour obtenir la solution Ç {min P(P + 'FT It'Fs) - 2 [4' g 177-i[Rtid } tel que He g '. + H utIu 15 ii. on estime la force d'excitation exercée par la houle incidente sur le moyen mobile w(k) par une équation de la forme : 0(1clic) = Ç(N). 6) Commande du système houlomoteur 20 Il s'agit d'une étape facultative. Lors de cette étape, on contrôle le système houlomoteur, de manière à prendre en compte la force d'excitation exercée par la houle incidente. Ainsi, il est possible de piloter le système houlomoteur, afin d'optimiser l'énergie récupérée. La commande peut consister en un contrôle du moyen mobile du système 25 houlomoteur, par exemple au moyen d'une machine de conversion électrique, pneumatique ou hydraulique, nommé système PTO (de l'anglais « power take-off »). Ce système PTO influence le mouvement du moyen mobile et permet de transférer l'énergie mécanique au réseau électrique, pneumatique ou hydraulique. La commande prédictive par modèle (MPC) est un exemple de méthode de commande de systèmes houlomoteurs. 30 Exemple illustratif Afin de montrer l'intérêt du procédé selon l'invention, un test du procédé selon l'invention a été réalisé en bassin. Le système houlomoteur correspond au système houlomoteur illustré en figure 1. Le système houlomoteur est équipé d'un capteur de couple 3049989 26 pour mesurer le couple exercé par la houle sur le moyen mobile. Préalablement, on détermine le couple exercée par la houle pour le cas où le contrôle est non actif (c'est-à-dire lorsque la machine de conversion ne génère pas de couple au système), et où le flotteur est maintenu à une position d'équilibre. Ces caractéristiques permettent de comparer le couple 5 d'excitation exercé par la houle incidente sur le moyen mobile pour ce premier test, à d'autres tests lorsque le flotteur est mobile, et la machine électrique activée. Pour la présente analyse, une concaténation de quatre types de houles différentes est considérée. Ces quatre types de houle W1, W2, W3, W4 (de l'amplitude la plus faible à l'amplitude la plus élevée) sont représentés en figure 3, qui illustre l'amplitude A en fonction 10 du temps T. Les figures 4a à 4c représentent une comparaison de la force (dans ce cas, le couple) de référence d'excitation exercée par la houle incidente sur le moyen mobile, avec la force (dans ce cas, le couple) d'excitation exercée par la houle incidente sur le moyen mobile déterminée selon un mode de réalisation de l'invention. Les figures 4a à 4c correspondent à 15 une houle de type W1 représentée sur la figure 3. La référence est notée REF, et est représentée par des traits pointillés. Les forces de la houle incidente déterminées par le procédé selon l'invention sont représentées par un trait continu. La figure 4a correspond au premier mode de réalisation de l'invention MR1, la figure 4b correspond au deuxième mode de réalisation de l'invention MR2, et la figure 4c correspond au troisième mode de réalisation 20 de l'invention MR3. On remarque sur ces figures 4a à 4c, que, pour les trois modes de réalisation, la courbe représentative de la force déterminée par le procédé selon l'invention est quasiment superposée à la courbe représentant de la force de référence. Ainsi, le procédé selon l'invention permet de déterminer de manière précise la force d'excitation de la houle incidente sur le moyen mobile. 25 Les figures 5a à 5c représentent une comparaison de la force (dans ce cas, le couple) de référence d'excitation exercée par la houle incidente sur le moyen mobile, avec la force (dans ce cas, le couple) d'excitation exercée par la houle incidente sur le moyen mobile déterminée selon un mode de réalisation de l'invention. Les figures 5a à 5c correspondent à une houle de type W3 représentée sur la figure 3. La référence est notée REF et est 30 représentée par des traits pointillés. Les forces de la houle déterminées par le procédé selon l'invention sont représentées par un trait continu. La figure 5a correspond au premier mode de réalisation de l'invention MR1, la figure 5b correspond au deuxième mode de réalisation de l'invention MR2, et la figure 5c correspond au troisième mode de réalisation de l'invention MR3. On remarque sur ces figures 5a à 5c, que, pour les trois modes de réalisation, la 35 courbe représentative de la force déterminée par le procédé selon l'invention est sensiblement superposée à la courbe représentant de la force de référence. Ainsi, le procédé selon l'invention permet de déterminer de manière précise la force d'excitation de la 3049989 27 houle incidente exercée sur le moyen mobile. On peut également observer que l'erreur entre la force de référence et la force déterminée est inférieure pour ces courbes que pour les courbes des figures 5a à 5c. REVENDICATIONS1. Procédé de détermination de la force d'excitation exercée par la houle incidente sur un moyen mobile (2) d'un système houlomoteur (1), ledit système houlomoteur (1) convertissant l'énergie de la houle en énergie au moyen dudit moyen mobile (2) coopérant avec une machine de conversion (3), ledit moyen mobile (2) effectuant un mouvement par rapport à ladite machine de conversion (3) sous l'action de la houle, caractérisé en ce qu'on réalise les étapes suivantes : a) on mesure la position et la vitesse dudit moyen mobile (2) ; b) on mesure la force F, exercée par ladite machine de conversion (3) sur ledit moyen mobile (2) ; c) on construit un modèle de la force Frad de radiation exercée sur ledit moyen mobile (2), ledit modèle de la force Frad de radiation reliant ladite force de radiation à la vitesse dudit moyen mobile (2) ; d) on construit un modèle de dynamique dudit système houlomoteur qui relie ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2), à ladite position dudit moyen mobile (2), à ladite vitesse du moyen mobile (2), à ladite force F, exercée par ladite machine de conversion (3) sur ledit moyen mobile (2) et à ladite force de radiation Frad exercée sur ledit moyen mobile (2) ; et e) on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2) au moyen dudit modèle de dynamique, dudit modèle de la force de radiation, desdites position et vitesse mesurées et de ladite force F, mesurée exercée par ladite machine de conversion (3) sur ledit moyen mobile (2). 2. Procédé selon la revendication 1, dans lequel on construit ledit modèle de dynamique dudit système houlomoteur au moyen d'une équation du type : M2(t) = Fex(t) + Fhd(t)+F'd(t) - Fu(t) avec M la masse du moyen mobile, 2 l'accélération dudit moyen mobile (2), Fhd la force de rappel hydrostatique exercée sur ledit moyen mobile (2). 3. Procédé selon la revendication 2, dans lequel on détermine ladite force Fhd de rappel hydrostatique au moyen d'une formule du type : Fhd(t) = -Kz(t) avec z la position dudit moyen mobile et K le coefficient de raideur hydrostatique. 4. Procédé selon la revendication 2, dans lequel on détermine ladite force Fhd de rappel hydrostatique par une fonction affine par morceau du type : Fhd = -Fiz(t) - Gi pour au moins deux intervalles Sil de variation de ladite position z, avec F, et G, des coefficients constants de ladite fonction affine pour chaque intervalle ni. 3049989 29 5. Procédé selon l'une des revendications précédentes, dans lequel on construit ledit modèle de la force Frac de radiation par une équation du type : Frad(t) = - M '2(t) - Fr(t) avec Fr(t) = f ot h(t - T)2(t)clt l'amortissement de radiation, 5 2 l'accélération dudit moyen mobile (2), M' la masse ajoutée à la fréquence infiniment élevée, 2 la vitesse dudit moyen mobile (2), h la réponse impulsionnelle qui relie la vitesse du moyen mobile à l'amortissement de radiation. 6. Procédé selon la revendication 5, dans lequel on mesure l'accélération dudit moyen mobile (2). 10 7. Procédé selon la revendication 6, dans lequel on détermine l'amortissement de radiation Fr au moyen d'une réponse temporelle approximée issue d'une réalisation d'état, du type : Fr(t) = Cr f ot e(t-T)ArBri(T) Dri(t) avec Ar, Br, Cr, Dr des matrices de la réalisation d'état. 15 8. Procédé selon les revendications 2 et 7, dans lequel on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2) au moyen d'une équation du type : Fex(k)=(M M')2(k) Eit101 Crd Arkd-1-1Brdi(1) Drdi(k) FiZ(k) Gi Fu(k) avec k le pas de temps dans le domaine discret, Ard, B rd, Crd, Drd les 20 matrices de la réalisation d'état dans le domaine discret. 9. Procédé selon la revendication 5, dans lequel on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2) au moyen d'un observateur d'état basé sur un filtre de Kalman linéaire construit à partir d'un modèle de marche 25 aléatoire de ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2). 10. Procédé selon la revendication 9, dans lequel on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2) en mettant en oeuvre les étapes 30 suivantes : i) on initialise k=0, le vecteur d'état 2,(010)=Fc,(0) et l'état de la matrice de covariance, P(010) = P0, ; 3049989 30 ii) à chaque instant k, on acquiert lesdites mesures de position et de vitesse dudit moyen mobile y(k) = [z(k) z(k)] et la mesure de ladite force exercée par ladite machine de conversion sur ledit moyen mobile u(k) = Fu(k) ; et iii) à chaque instant k, on détermine la force d'excitation exercée par la houle 5 incidente sur ledit moyen mobile (2) 0(k) au moyen des équations suivantes : 0(k) = [0 1] fca(k1k) avec : 1fca(klk -1) = Au5eu(k -11k -1) + BaFu(k -1) P(klk -1) = AaP(k -11k - flATa + Q K(k) = P(kIk - 1)Cj; (CaP(klk - 1)Cj; + R)-1 .2a(k1k) = 2a(klk - 1) + K(k)(Y(k) Ca.2a(k ik - 1) DaFu(k)) P(klk) = (I - K(k)Ca(k))P(klk - 1) avec K(k) le gain de l'observateur, xa(k) le vecteur d'état inconnu xa(k) = [w )1 x( Aco Ba, Ca, Da des matrices de la réalisation d'état, P la matrice de () 10 covariance du vecteur d'état, Q et R des matrices de calibration. 11. Procédé selon la revendication 5, dans lequel on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2) au moyen d'un estimateur d'état à entrée inconnue basé sur une approche d'horizon glissant. 15 12. Procédé selon la revendication 11, dans lequel on détermine ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile (2) en mettant en oeuvre les étapes suivantes à chaque instant k : i) on acquiert lesdites mesures de position et de vitesse dudit moyen mobile 20 y(k) = [z(k) z(k)] et la mesure de ladite force exercée par ladite machine de conversion sur ledit moyen mobile u(k) = Fu(k) ; et ii) on détermine la force d'excitation exercée par la houle incidente sur ledit moyen mobile ei(klk) au moyen des équations suivantes : 1-4)(k lk) = «N) avec «N) obtenu par résolution d'un problème de programmation 25 quadratique suivant : { S i el mir4. 1( (P + WIRWsg - 2 [01; -T X0 V 1 T '[RWsiS i avec P, R et S des matrices tel que He ,g'.+ HdUd formées à partir de matrices de pondération de la fonction coût, la matrice qui représente l'ensemble des valeurs des variables inconnues sur la fenêtre d'estimation de dimension N, W la matrice, calculée à partir de la représentation 30 d'état, qui correspond à l'effet de ces variables sur les valeurs des sorties sur cette même fenêtre, H et Hd des matrices et g un vecteur obtenus à partir de 3049989 31 contraintes imposées sur l'état et l'entrée inconnue, Jeo = 2(k - NIk - 1) et 00 = 0(k - NIk - 1) avec x le vecteur d'état et w l'entrée ud(k - N) ud(k - N + 1) inconnue de l'estimateur d'état, Ud = u- [ , ud étant l'entrée connue lid(k) de l'estimateur d'état avec ud(t)=Fhd(t) - Fu(t). 13. Procédé selon l'une des revendications précédentes, dans lequel on commande ledit système houlomoteur en fonction de ladite force d'excitation Fex exercée par la houle incidente sur ledit moyen mobile, déterminée (2). 5 10