Mise en équation d'un problème de géométrie dans l'espace pour un héliostat.

[pierrepoulpe]

Bonjour à tous,

Je bute sur un problème qui ne rentre pas dans le cadre des cas typiques étudiés à l'école. J'aurai besoin d'un peu d'aide pour bien poser le problème.

Tout d'abord, je vous explique le contexte : j'ai construit un héliostat (un dispositif qui suit le soleil) pour positionner un miroir de telle sorte que les rayons du soleil soient toujours réfléchi dans la même direction. Ça diffère un peu d'un héliostat qui sert à exposer de façon optimale un panneau photovoltaïque : dans mon cas j'ai besoin de beaucoup plus de précision. Un écart d'un degré, et je suis 40 cm à côté de ma cible à une distance de 24m, alors que dans le cas d'un panneau à exposer ca représente une perte de surface incidente 0,01%....

Cet héliostat est composé de 3 parties :
- un châssis fixé sur un toit
- un 'chariot' qui est articulé au châssis selon un axe approximativement 'est ouest' et qui sert à contrôler la 'hauteur'
- un miroir (monté sur un cadre) articulé au chariot selon un axe perpendiculaire au premier, appelons le azimut (même si c'est faux)

Je sais mesurer avec précision un angle relatif entre les parties mobiles et leurs parties porteuses respectives. J'utilise un codeur incrémental sur lequel est monté un pignon entrainé par une courroie synchrone collée sur une poulie sur la partie mobile. Une butée électrique me permet d'étalonner le dispositif à chaque démarrage et donc d'avoir des valeurs répétables.

Je sais déterminer avec précision pour quel angle sur l'axe hauteur j'ai mon axe azimut qui est vertical (avec un niveau à bulle).
Je sais déterminer avec précision pour quel angle sur l'axe azimut j'ai la normale de mon miroir perpendiculaire à l'axe hauteur.

Pour effectuer tout mes calculs, je raisonne dans un espace orthonormé avec x qui correspond à l'axe est-ouest, y nord-sud, z verticale. (nord s'entend géographique, le magnétique ne m'est d'aucun secours)
Dans cet espace je sais calculer avec précision la position apparente du soleil à partir de la date et de l'heure.

Ce que je suis incapable de déterminer avec précision, avec les outils à ma disposition, c'est le vecteur héliostat-cible (qui est approximativement un vecteur horizontal orienté au sud), et le vecteur correspondant à l'axe hauteur de mon héliostat (qui est approximativement un vecteur horizontal orienté ouest)

Ce que je pensais faire :
- Me placer au niveau de ma cible, et orienter mon miroir jusqu'à ce que la normale de mon miroir soit confondue avec mon vecteur cible : c'est simple, c'est quand je me vois dans le miroir! (je peux utiliser une longue vue, ou un pointeur laser) => je connais alors les angles (1) sur mon héliostat correspondant à mon vecteur cible.

- Je règle manuellement la position de l'héliostat pour que les rayons du soleil pointent vers ma cible. Je note la position du soleil dans mon repère géographique à l'instant T, et les angles (2) de mon héliostat.
Connaissant les angles (1) et (2) de mon héliostat, je peux déterminer l'angle entre la normale de mon miroir à l'instant T, et ma cible. En multipliant par 2, j'ai l'angle entre ma cible et le soleil à l'instant T, nommons le alpha.

- Si je répète quelques heures après l'opération précédente (à T'), j'obtiens une autre paire d'angles (3), et de je détermine un angle alpha'.

- Connaissant deux vecteurs (Vsoleil et Vsoleil'), les deux angles formés (alpha et alpha') avec un 3ème vecteur inconnu (appellons le Vcible), il me semble que je dois être capable de déterminer ce 3ème vecteur?
Ma première question : comment je résous ce premier problème? je ne sais pas comment le poser.

- Connaissant deux vecteurs (Vsoleil et Vcible), les angles de rotation appliqués autour de deux vecteurs inconnus, donc on sait néanmoins qu'il sont perpendiculaires, est on capable de déterminer ces deux vecteurs inconnus? et comment?



Si vous avez besoin de photos pour visualiser, j'en ferai.


Merci pour votre aide.
-----

[Dlzlogic]

Bonjour,
Je résume à ma façon, pour être sûr d'avoir bien compris.
Vous avez 2 objets fixes :
1- l'héliostat qui possède deux axes de rotation. Vous connaissez la position du soleil à partir de ce point.
2- une cible. Le but est que la cible reçoive à tout moment l'image du soleil réfléchie par le miroir de l'héliostat.

La question est donc : "quels doivent être les deux angles, vertical et horizontal pour obtenir cela ?".
Imaginons le soleil exactement au sud. Vous aurez pris soin de positionner la cible exactement au sud de votre héliostat.
A cet instant la correction de l'angle horizontale est nulle. La bissectrice de l'angle Soleil-Héliostat-Cible est orthogonale au miroir.
Cet angle SHC dépend de la hauteur du soleil (un angle variable) et de l'angle entre H et C, fixe et qu'il faut connaitre avec une bonne précision.
Lorsque le soleil n'est pas au sud, le raisonnement concernant les angles est le même, au lieu de faire le calcul dans le plan vertical, on le fait dans le plan horizontal.

A la relecture, je me demande si votre question n'est pas justement de calculer, ou mesurer la position relative de la cible par rapport à l'héliostat.
Pour l'avoir exactement au sud, puisque vous connaissez la position du soleil à midi ( heure solaire locale), le plus simple me parait être de faire une trace à ce moment précis. Pour la différence d'altitude, avec un théodolite, c'est pas très difficile, encore faut-il en avoir un. Ce serait un bon exercice pour des élèves topographes.

[pierrepoulpe]

Bonjour, merci pour votre réponse, et désolé pour mon retour tardif.

1- l'héliostat qui possède deux axes de rotation.

Tout à fait

Vous connaissez la position du soleil à partir de ce point.

Je connais la position du soleil à partir de ce point, mais dans un référentiel que je vais appeler 'géographique'. Mais les axes de rotation de mon héliostat ne son pas alignés sur les axes de mon référentiel géographique. Même si le premier axe de mon héliostat est à peu près horizontal et orienté est-ouest, il est à considérer comme quelconque et inconnu.

2- une cible. Le but est que la cible reçoive à tout moment l'image du soleil réfléchie par le miroir de l'héliostat.

Exactement

La question est donc : "quels doivent être les deux angles, vertical et horizontal pour obtenir cela ?".
Imaginons le soleil exactement au sud. Vous aurez pris soin de positionner la cible exactement au sud de votre héliostat.

La question est effectivement de calculer les deux angles à appliquer sur les axes de l'héliostat pour obtenir cela. Mais comme dit plus haut, l'héliostat n'est qu'approximativement aligné sur un axe est-ouest, pour rajouter un peu de complexité, ma cible n'est pas plein sud par rapport à mon héliostat.

Je pourrais travailler le temps de cet 'étalonnage' avec une cible plein sud, encore faut il déterminer le sud avec précision. Une boussole, même avec un viseur type course d'orientation, et même en incluant la dérive magnétique, ne me permettra de dépasser le degré de précision.

Quant à se servir des rayons du soleil à midi pile (au soleil), je ne peux pas me servir de l'héliostat car je ne sais toujours pas l'aligner 'pile au sud', et en utilisant l'ombre d'un objet suffisamment haut pour avoir un peu de précision (j'ai 24m entre mon héliostat et ma cible), il faut que j'y réfléchisse, mais j'ai peur d'accumuler des erreurs, et surtout de ne pas être capable d'estimer ces erreurs cumulées.

A cet instant la correction de l'angle horizontale est nulle. La bissectrice de l'angle Soleil-Héliostat-Cible est orthogonale au miroir.
Cet angle SHC dépend de la hauteur du soleil (un angle variable) et de l'angle entre H et C, fixe et qu'il faut connaitre avec une bonne précision.
Lorsque le soleil n'est pas au sud, le raisonnement concernant les angles est le même, au lieu de faire le calcul dans le plan vertical, on le fait dans le plan horizontal.

A la relecture, je me demande si votre question n'est pas justement de calculer, ou mesurer la position relative de la cible par rapport à l'héliostat.
Pour l'avoir exactement au sud, puisque vous connaissez la position du soleil à midi ( heure solaire locale), le plus simple me parait être de faire une trace à ce moment précis. Pour la différence d'altitude, avec un théodolite, c'est pas très difficile, encore faut-il en avoir un. Ce serait un bon exercice pour des élèves topographes.

En fait mon code fait tous ses calculs avec des vecteurs, et j'ai déjà prévu d'avoir l'axe premier de mon héliostat sous la forme d'un vecteur quelconque dans le référentiel géographique. Pour l'instant, il est égal à {1 0 0} en première approximation. Mais je n'ai aucun soucis pour calculer les 2 angles de l'héliostat connaissant le vecteur soleil Vsol, le vecteur axe de rotation Vaxehorizontal , et le vecteur cible Vcible:

Je somme et normalise Vsol + Vcible, ce qui me donne la bissectrice, et donc la normale que doit adopter mon miroir : Vnormale
Je fait le produit vectoriel (+normalisation) de Vnormale * Vaxehorizontal, ce qui me donne la normale(=perpendiculaire) au plan passant par ces deux vecteurs : c'est mon second axe de rotation Vaxevertical (qui n'est pas toujours vertical, c'est un abus de langage)
Je fait le produit vectoriel de Vaxeverticale * Vaxehorizontale, ce qui me donne la projection de ma normale miroir sur le plan perpendiculaire à l'axe horizontal Vnormaleproj
Je calcul l'angle à appliquer sur le second axe avec le produit scalaire entre la normale miroir, et la normale projetée sur le plan perpendiculaire au premier axe. Vnormaleproj . Vaxehorizontale
Et je calcule l'angle du premier axe avec le produit scalaire entre la verticale Vvert et la normale projetée Vnormaleproj.

Voici la fonction principale, c'est en C (arduino), la librairie de fonction vectorielle est maison, le 3 ème paramètre est le résultat de la fonction.
Les 2 dernières lignes convertissent les angles en degré en nombre de pulsation à compter sur un codeur incrémental entrainé par les partie mobiles en rotation.

void CalculCible() {  
  dt = now();

  dtSolar = SolarTime(dt); // heure en UTC
  
  ElevationSoleil = Hauteur(dtSolar); //radians
  AzimutSoleil = Azimut(dtSolar); //radians
  // z+ haut, z- bas, x- est x+ ouest, y+ sud, y- nord
  soleil.x = cos(ElevationSoleil) * sin(AzimutSoleil);
  soleil.y = cos(ElevationSoleil) * cos(AzimutSoleil);
  soleil.z = sin(ElevationSoleil);

  SommeVectoriel(&soleil,&(cible),&normaleMiroir);

  ProdVectoriel(&normaleMiroir,&axeH,&axeV);
  ProdVectoriel(&axeH,&axeV,&projAxeH);
  AhCible = acos(ProdScalaire(&normaleMiroir,&projAxeH))  * radToDeg; //degré
  if (acos(ProdScalaire(&axeH,&normaleMiroir))  * radToDeg > 90) 
    AhCible = AhCible * -1;

  ElevationCible = 90 - (acos(ProdScalaire(&verticale,&projAxeH))  * radToDeg); //degré

  regul[V].cible = ElevationCible * regul[V].degToPulse + regul[V].pulseOffset;  
  regul[H].cible = AhCible * regul[H].degToPulse + regul[H].pulseOffset;
}

Pour résumer, mon problème est de déterminer Vaxehorizontale (axeH dans mon code), connaissant pour 2 instants donnés, le vecteur soleil, ElevationCible, AhCible (cible étant ici à entendre au sens cible de ma boucle de régulation d'angle, rien à voir avec la cible de renvoi de lumière, désolé pour la source de confusion)

PS : sans théodolite que je ne possède pas, je saurai déterminer simplement et précisément la différence d'altitude entre héliostat et cible : 30m de tuyaux d'arrosage, deux bouteilles en plastique....

[Dlzlogic]

Bonjour,
Apparemment, vous avez cerné le problème, il ne reste plus(!) que l'orientation.
Par ailleurs vous avez bien compris qu'il est difficile de calculer la précision, c'est à dire l'erreur commise. Voila la méthode que je vous propose : vous adoptez au mieux, un position, en notant soigneusement tous les éléments, et vous mesurez une dizaine de résultats. Vous adoptez une seconde position, toujours en notant soigneusement tous les éléments, et de nouveau une dizaine de résultats. Avec 2, ça devrait suffire.
Pour chacune des position, et chacun des résultats, vous aurez une valeur, par exemple l'heure, et un écart sur la position idéale, distance en hauteur et en azimut, par rapport à la cible.
C'est une méthode d'on appelle l'étalonnage. Dans votre cas, c'est celle qui me parait la préférable.
Si vous n'arrivez pas à exploiter les deux listes, n'hésiter pas à me les envoyer.

[A voir en vidéo sur Futura]

[pierrepoulpe]

Merci pour cette réponse rapide.

1/ Je ne suis pas sûr d'avoir compris votre proposition. Pouvez vous préciser ce que vous entendez par position et que vous faites varier 2 fois? Par résultat je comprend l'écart entre le rayon réfléchit et la cible, ainsi que la date/heure.
Ce n'est pas le plus évident d'un point de vue pratique, il faut qu'il fasse suffisamment beau pour que je puisse relever une dizaine de points dans la journée avec des écarts de temps significatifs, et vu la complexité du bâti (nombreux décrochés de plusieurs mètres, terrain étroit sans accès chez les voisins) ca risque d'être compliqué.

2/ En parallèle je persiste à essayer de trouver une solution mathématique avec les données que j'ai en entrée.
Je vais essayer de reformuler mon problème :

Soit deux repères d'origine commune : Rgeo et Rheliostat
Soit 3 vecteurs Vsol, Vsol', Vcible, et 2 bissectrices Vsol-cible et Vsol'cible

Je connais Vsol et Vsol' dans Rgeo, et je peux donc calculer l'angle entre ces deux vecteurs.
Je connais Vcible, Vsol-cible et Vsol'cible dans Rheliostat. Je peux donc calculer les angles entre Vcible et respectivement Vsol-cible / Vsol'-cible. En multipliant par deux, j'obtiens respectivement les angles entre Vcible / Vsol et Vcible / Vsol'

Je souhaite connaitre Vcible dans Rgeo.
Pour résumer j'ai deux vecteurs connus, un vecteur inconnu, et les 3 angles entre ces 3 vecteurs.

Connaissant un vecteur, l'angle avec un second vecteur, intuitivement, il me semble que l'ensemble des solutions de ce second vecteur décrit un cône.
Si je fais entrer le second vecteur connu en jeu, j'obtiens un deuxième cône. Je sais que ma solution Vcible fais partie de ces deux cônes.

Mon problème ne se résume il pas à trouver les deux demi-droites d'intersection entre deux cônes de sommet commun? Si oui, comment cela se résout il?
J'ai bien conscience que j'obtiendrai deux solution possibles, ce n'est pas un problème pour mon cas pratique, je n'aurai aucun mal à reconnaitre la bonne solution.

[Dlzlogic]

Lorsqu'il s'agit de changement de repère en 3D, je n'utilise jamais les angles, c'est trop compliqué, mais j'utilise la formule suivante
X = TX + XX.x + XY.y + XZ.z
Y = TY + YX.x + YY.y + YZ.z
Z = TZ + ZX.x + ZY.y + ZZ.z
La méthode est de calculer les coordonnées de 4 points dans les deux référentiels et de résoudre le système à 12 inconnues.
Dans votre cas, il me semble que la difficulté est justement que vous ne connaissez pas la position relative des 2 systèmes, toutes vos valeurs sont approximatives, donc pas vraiment utilisables. C'est pourquoi je vous proposais cette solution d'étalonnage. Il n'est pas du tout indispensable que les 10 mesures soient faites la même journée, ni qu'elles soient faites à des intervalles régulier, la seule chose importante est d'en avoir assez pour avoir un bon résultat. Comme je n'ai aucune idée de la précision de l'ensemble, 2 fois dix me paraissent une bonne approche.
En tout cas, je reste persuadé qu'une méthode dont l'élément principal est la mesure me parait la meilleure.
Par ailleurs, pour résoudre mathématiquement cela, la seule solution est de travailler indépendamment dans le plan vertical et dans le plan horizontal, mais comme votre axe est "à peu près horizontal", votre cible "à peu près au sud", tout calcul rigoureux me parait risqué.

[pierrepoulpe]

Bonjour,

Ça fait quelque mois, mais j'aime bien donner des retours.

Premièrement, je me suis résolu à déterminer le vecteur héliostat->cible par la mesure :
Pour l'élévation, c'est relativement facile, un tuyau d'arrosage, deux bouteilles....
Pour l'azimut, 3 fil à plomb, un premier au niveau de l'observateur, et à une distance de quelques mètres, un assistant qui positionne les deux autres fils à plomb, l'un de telle sorte à être d'en l'alignement avec l'héliostat, l'autre avec le clocher de l'église à 1km de là. Geoportail me donnant directement l'azimut entre mon point d'observation et le clocher, qui par le plus grand des hasard est exactement au nord (à l'erreur de positionnement près), je me retrouve avec un simple problème de trigo....

Secondo, ca m'énervait de ne pas trouver de solution mathématique à mon problème. Et je suis du genre un peu têtu....
En plus, ca me permet, en répétant le calcul plus de deux fois, d'évaluer la dispersion et donc la précision.

Bref, après moult recherches, la solution était relativement simple, comme toujours.
Je connais deux vecteurs, je connais les deux angles qu'ils forment respectivement avec un troisième vecteur inconnu, et je sais que ces trois vecteurs sont normés (norme=1). angle ~= produit scalaire

Soit U et V mes vecteur connus, W mon vecteur recherché,
Soit
U.W = i
V.W = j
(et |W| = 1)
En développant le produit scalaire avec les trois composantes x y z, je me retrouve avec un système de 2 équations à 3 inconnues.

xu*xw + yu*yw + zu*zw = i
xv*xw + yv*yw + zv*zw = j

je passe xu/xv à droite
xw = i/xu - yu/xu * yw - zu/xu * zw
xw = j/xv - yv/xv * yw - zv/xv * zw

"Combinaison linéaire"* : J'additionne le deuxième équation avec la première multiplié par le facteur qui va bien pour faire disparaitre yw : -(yv*xu)/(xv*yu)
xw * (-(yv*xu)/(xv*yu) + 1) = -(yv*xu)/(xv*yu) * i/xu + j/xv
-(yv*xu)/(xv*yu) * -zu/xu * zw - zv/xv * zw
-(yv*xu)/(xv*yu) * -yu-xu * yw -yv/xv * yw
=> (-yv*xu*yu / xv*yu*xu - yv/xv) * yw
=> (yv/xv - yv/xv) * yw
=> 0 * yw


xw * (-(yv*xu)/(xv*yu) + 1) = -(yv*xu*i)/(xv*yu*xu) + j/xv + ((yv*xu*zu)/(xv*yu*xu) - zv/xv) * zw
= a = b = c
Pour lisibilité
a = (-(yv*xu)/(xv*yu) + 1)
b = -(yv*xu*i)/(xv*yu*xu) + j/xv
c = ((yv*xu*zu)/(xv*yu*xu) - zv/xv)

xw * a = b + c * zw

-b + a*xw = c*zw
-b/c + a/c*xw = zw

zw = -b/c + a/c*xw
zw = (a*xw - b) / c

Et de même :
d = (-(zv*xu)/(xv*zu) + 1)
e = -(zv*xu*i)/(xv*zu*xu) + j/xv
f = ((zv*xu*zu)/(xv*zu*xu) - yv/xv)
yw = (d*xw - e) / f


Ensuite je réintroduis le fait que la norme de W est de 1
xw² + yw² + zw² = 1
xw² + ((d*xw - e) / f)² + ((a*xw - b) / c) ² = 1
vw² + ((d*xw)² - 2*e*d*xw + e²) / f² + ((a*xw)² - 2*b*a*xw + b²) / c² = 1
vw² + d²/f²*xw² - (2*e*d)/f²*xw + e²/f² + a²/c²*xw² - (2*b*a)/c²*xw + b²/c² = 1
(1 + d²/f² + a²/c²) * xw² + (-2ed/f² - 2ba/c²) * xw + (e²/f² + b²/c² - 1) = 0

A = (1 + d²/f² + a²/c²)
B = (-2ed/f² - 2ba/c²)
C = (e²/f² + b²/c² - 1)

D = B² - 4AC

xw = -B - racine(D) / 2A
ou
xw = -B + racine(D) / 2A

Je retrouve sans soucis yw et zw, et il me reste plus qu'à éliminer la mauvaise solution (devrait être facile, comme j'ai une connaissance approximative de la bonne réponse)


* Je me suis largement insipiré de ça
http://www.google.fr/url?sa=t&rct=j&...57752919,d.d2k


Question : Peut on simplifier mon raisonnement? Ou tout du moins le rendre plus lisible?

Comment gérer le cas ou une des composantes de mes vecteurs en entrée serait nulle, et risquerait de faire des divisions par zéro? (les seul cas dans la vraie vie, ce serait quand le soleil est plein sud, plein ouest/est (possible printemps été), et au lever/coucher)

Pierre

[pierrepoulpe]

Et j'ai fait le test en python et ca marche, miracle. Désolé pour les noms de variables, il sont pas cohérents avec ma démonstration.

#!/usr/bin/env python
# -*- coding: utf8 -*-
import math

def Norme(v):
	return (v[0]**2+v[1]**2+v[2]**2)**0.5

def Normalize(v):
	norme = Norme(v)
	v[0] = v[0] / norme
	v[1] = v[1] / norme
	v[2] = v[2] / norme
	return 

def ProdScalaire(v1,v2):
	return (v1[0] * v2[0]) + (v1[1] * v2[1]) + (v1[2] * v2[2])

def Angle(v1,v2):
	return math.acos(ProdScalaire(v1,v2)) * 180 / 3.1415927

v1 = [-0.8,1,0.1]
v2 = [0.7,0.4,0.1]
v3 = [0.2,0.7,1.05]

Normalize(v1)
Normalize(v2)
Normalize(v3)

#print Angle(v1,v2)
#print Angle(v2,v3)
#print Angle(v1,v3)

i = ProdScalaire(v1,v3)
j = ProdScalaire(v2,v3)

xa = v1[0]
ya = v1[1]
za = v1[2]
xb = v2[0]
yb = v2[1]
zb = v2[2]
xz = v3[0]
yz = v3[1]
zz = v3[2]

print "xz : " + str(xz)
print "yz : " + str(yz)
print "zz : " + str(zz)
print


a = 1 - ((yb * xa) / (xb * ya))
b = (-1 * yb *xa * i) / (xb*ya*xa) + j/xb
c = (za * yb * xa) / (xa * xb * ya) - (zb / xb)

d = 1 - ((zb * xa) / (xb * za))
e = (-1 * zb *xa * i) / (xb*za*xa) + j/xb
f = (ya * zb * xa) / (xa * xb * za) - (yb / xb)

eqa = 1 + d**2/f**2 + a**2/c**2
eqb = (-2*e*d) / f**2 + (-2*b*a) / c**2
eqc = e**2/f**2 + b**2/c**2 - 1


eqd = eqb ** 2 - 4*eqa*eqc
xz1 = (-1 * eqb - eqd ** 0.5) / (2 * eqa)
xz2 = (-1 * eqb + eqd ** 0.5) / (2 * eqa)
yz1 = a/c*xz1 - b/c
yz2 = a/c*xz2 - b/c
zz1 = d/f*xz1 - e/f
zz2 = d/f*xz2 - e/f

print xz1
print xz2
print yz1
print yz2
print zz1
print zz2

[pierrepoulpe]

correction, le produit scalaire donne le cosinus, pas directement l'angle. Mais ca ne change rien à le recherche des inconnues.