Pilote automatique & PID

[Hitman458]

Hello,

Je me tente à la programmation d'un pilote automatique, j'ai déjà tapé du code pour tester le PID mais..... ça ne marche pas -________-"

Code : Tout sélectionner

erreur = consigne - mesure
commande = kp * erreur + ki * integrale_erreur + kd * derivee_erreur

Pour régler le PID, on commence par le correcteur proportionnel, donc on met ki=0 et kd=0, pour simplifer :

Code : Tout sélectionner

commande = kp *erreur

et c'est ce qu'on trouve sur le wiki de ksp, plus exactement :

Code : Tout sélectionner

Vector3d error = (Vector3d)vessel.transform.InverseTransformDirection((Vector3)dir);
yaw_command = k * error.x;
pitch_command = -k * error.z;

Que ce passe-t-il concrètement dans KSP ? Mon satellite oscille infiniment autour de l'orientation demandée. L'amplitude de l'oscillation reste la même sans jamais converger vers la valeur demandée. (Bon, si on active le SAS ça marche car le SAS freine l'oscillation, mais je considère que c'est tricher, j'aimerai que ça fonctionne sans le SAS).

Après avoir regardé mon satellite faire du gangnam style pendant des heures et après de multiples tentatives de réglage de kp, ki et kd, j'ai soudain eu une illumination : c'est logique que ça ne marche pas avec cette équation. Les commandes yaw et pitch sont en quelque sorte des intégrateurs. Plus yaw/pitch va être grand, plus l'accélération va être grande, or on ne veut pas régler une accélération, on veut régler une position.

J'ai regardé le code source de MechJeb et c'est bien compliqué, ils ont bien un PID dans leur code mais c'est un PID trafiqué, avant d'injecter la commande, il font tout un tas de manipulation de vecteur mais leur code n'est pas très bien commenté. En gros ils pondèrent la commande avec un vecteur d'inertie et un vecteur de moment et tout un tas d'autres coefficients... ce n'est plus vraiment un PID là en fait

J'ai pas envie de copier coller du code que je ne comprend pas, et ça perdrais tout l'intérêt de ma démarche. J'aimerais réussir à faire un PID simple sans user de formules magiques.

Quelques liens pour mieux comprendre :
1) Comment faire un PID
2) Le code du PID de MechJeb
3) Le code qui récupère la valeur du PID et qui fait de la sorcellerie dessus avant d'injecter les valeurs dans les commandes de KSP

Depuis que j'ai eu mon bac, j'ai bouffé du PID chaque jour à toutes les sauces pendant 6 ans sauf que.... on ne nous a jamais appris à en programmer un algorithme, c'est un peu frustrant . Surtout qu'au final, quand je vois ce que je trouve sur internet concernant la programmation d'un PID, on est très loin des délires de transformés de Laplace, d'études fréquentielles, de réponses indicielles, de Dirac,... bref de tout ce que j'ai vu en cours

Si j'arrive à faire fonctionner ce PID, j'envisage de faire un petit tuto sur le développement d'addon dans KSP

Quelqu'un aurait une idée d'algo qui fonctionnerait ?

[Stellaris]

Héba, tu te lances dans du lourd !
Z'auraient pu le commenter leur code, la source de MJ est quand même la caverne d'Ali Baba en terme d'infos !
Désolé de ne pas pouvoir t'aider, mais tu peux peut être poster un message sur le fofo US : Ils ont beaucoup plus d''expérience en terme de modding !

[Hitman458]

Je viens d'essayer de traduire le code de MechJeb en français, une sorte de reverse engineering de code xD
Maintenant je sais ce qu'ils font mais certains passages relèvent du mystère Je suis remonté jusqu'aux premiers commits de leur repository pour essayer de voir l'évolution de l'algo J’espérais qu'ils soient partis d'un algo simple au tout début mais en fait non dès le premier commit de code on arrive déjà à peu près à l'algo final qui est utilisé aujourd'hui....

Je continue de creuser leur code, je publierais peut être le résultat de mes recherches

[Hitman458]

Comment fonctionne le pilote automatique de MechJeb en théorie :

A) 1er étape, calculer l'erreur :
1) on calcule la précision : une formule magique, y a rien a comprendre c'est de la sorcellerie, c'est du min et max avec le moment_disponible + la_poussé_dispo * la_poussé * 20 / |moment_d'inertie|
2) attitudeReference : on détermine le repère par rapport auquel on veut manœuvrer.
3) attitudeTarget : on détermine l'attitude cible qu'on veut atteindre.
4) target: on calcule la cible à atteindre dans le repère de référence.
5) delta : c'est l'erreur entre la l'attitude actuelle et la l'attitude cible.
6) deltaEuler : delta converti en angles d'Euler, au lieu d'avoir des angles compris entre 0° et 360°, on a des angles compris entre -180° et 180°.

On entre maintenant dans la partie la plus mystérieuse :

7) torque : on calcule le moment disponible, moment = moment_disponible + moment_de_la_poussé * la_poussé_des_réacteurs
8) inertia : on calcul l'inertie, inertia = sign(moment_angulaire) * moment_angulaire^2 * moment * moment_d'inertie oO j'ai p'tet vu ça en prépa mais je m'en rappel pas
9) err : c'est le deltaEuler (calculé en A.6) converti en radian
10) on ajoute l'inertie à cette erreur
11) erreur = erreur * moment_inertie * moment^(-1) (le ^(1) c'est la matrice inverse)

Conclusion : Mechjeb utilise les moments pour pondérer l'erreur en position... pourquoi ? je sais pas Mais c'est cette sorcellerie là qui devrait empêcher le satellite d'osciller infiniment sans jamais converger vers la consigne. J'aimerais bien comprendre ce schmilblick.

B) 2iem étape, le PID :
là c'est pas très compliqué.
1) On calcule l'intégrale de l'erreur
2) On calcule la dérivé de l'erreur.
3) resultat = kp * erreur + ki * intégrale de l'erreur + kd * dérivée de l'erreur
4) on contraint resultat à rester entre une valeur min et une valeur max

C) 3iem étape, injecter le résultat du PID dans les commandes de KSP :
1) on calcul la limite d'application du lacet et du tangage : drive_limit = |erreur| * constante / precision. (précision calculé en A.1)
2) on contraint le restultat du PID à rester entre -drive_limit et +drive_limit
3) une autre formule magique, resultat = resultat_précédent + (resultat - resultat_précédent) * une_valeur
4) on recontraint le nouveau resultat à rester entre -drive_limit et +drive_limit (et là surprise certains paramètres transmis à la fonction ne sont pas utilisés).

J'appelle ça une usine à gaz °_°

[Hitman458]

Re, j'ai réussi à faire un PID qui fonctionne, je ne sais pas pourquoi cela ne marchait pas avant mais maintenant c'est bon. (J'espère que le fait que ça marche maintenant n'a rien à voir avec la dernière MAJ).

Du coup maintenant ça devient plus facile pour faire des alignements et des sychronisations d'orbites xD

https://github.com/kerbalspaceprogram-fr/Autopilot

[Mahzel]

Pour ton calcul de vecteurs (celui ou tu as mis ???? Oo ...... MechJeb2) j'ai trouvé ça sur le wiki KSP

Code : Tout sélectionner

//unit vectors in the up (normal to planet surface), east, and north (parallel to planet surface) directions
Vector3d eastUnit = vessel.mainBody.getRFrmVel(position).normalized; //uses the rotation of the body's frame to determine "east"
Vector3d upUnit = (position - vessel.mainBody.position).normalized;
Vector3d northUnit = Vector3d.Cross(upUnit, eastUnit); //north = up cross east

Je comprends rien au code de MJ sur ce sujet, par contre celui là je le comprends. Je suppose juste que MJ calcule le "nord"... Calcul un peu trop poussé pour ma pomme.
Tant que ça marche comme tu l'as fait celà-dit, c'est l'essentiel. Je suis intéressé si tu fait un test avec cependant (mes tests sur le sujet viendront pas avant quelques temps :p)

[Hitman458]

Oui oui tu as raison, la version du wiki est même la première que j'avais implémentée et elle est simple et logique ˆˆ

Par contre le résultat était faut...... le sens du vecteur était inversé.

upUnit, y a rien à dire c'est ok
northUnit idem.

J'en ai déduit que eastUnit pointait vers l'ouest. Je voulais pas rajouter un moins comme un barbare sans comprendre.
C'est getRFrmVel le problème.

The reference frame of the planet's surface rotates with respect to an inertial reference frame. For example, a ship landed on the planet's surface is actually moving with respect to an inertial frame, because the planet is rotating. This function gives the difference between a velocity in the rotating frame and the same velocity in the inertial frame. You can use it to convert an inertial frame velocity (like the one returned by Orbit.GetVel()) to a rotating frame velocity by (rotating frame velocity) = (inertial frame velocity) - body.GetRFrmVel(position).

J'ai supposé une chose mais je ne sais pas si c'est ça :
Si le vaisseau tourne dans le même sens que la rotation de la planète, alors getRFrmVel pointe vers l'Est
Si le vaisseau tourne dans le sens inverse de la rotation de la planète, alors getRFrmVel pointe vers l'Ouest.

P'tet que j'avais fait une erreur, je ferais encore un essai avec la version du wiki car je la préfère.

Je vois un peu l'idée du calcul de MechJeb mais il faut être ultra farfelu pour oser imaginer une équation comme ça.

[Hitman458]

Ah!!!! j'ai été trop trop naïf, c'est une des premières choses qui m'était venue en tête en plus mais je l'avais pas vérifié : c'est le produit vectoriel qui est faux.

Vect(A,B) = -Vect(B,A) et pour avoir un repère orthonormé direct il faut faire nécessairement north = east ^ up et non pas north = up ^ east comme c'est écrit dans le wiki.

Faudrait faire corriger le wiki.

Code : Tout sélectionner

Vector3d northUnit = Vector3d.Cross (eastUnit, upUnit);

Maintenant c'est good, je vais pouvoir remplacer la ligne de MechJeb par ça.

edit :
Par contre j'aimerai quand même avoir la confirmation que getRFrmVel pointe toujours vers l'Est et ne dépend pas du sens de rotation de la planète. La plupart des planètes tournent dans le sens direct (anti-horraire), mais il n'y en a pas une qui tourne dans le sens indirect (horaire) ?

[Hitman458]

Je suis entrain de rajouter plus de fonctionnalités a mon pilote automatique (décollage automatique avec auto staging).

Mais pour aller plus loin, j'ai un petit soucis

Quelqu'un sait comment récupérer les Maneuver Nodes en C# ? Et d'en créer via le code ?

Concrètement, j'ai besoin de calculer du deltaV, du temps de burn, le temps restant avant le passage au Node.

Il est possible de reprogrammer le code de l'algo/calcul, mais j'aimerai utiliser le système natif de KSP pour m'épargner cette peine.

edit :: pour créer un ManeuverNode :

Code : Tout sélectionner

double ut = Planetarium.GetUniversalTime();
Vector3d nodeDV =  new Vector3d(radialPlus, normalMinus, prograde)  ;
ManeuverNode mn = vessel.patchedConicSolver.AddManeuverNode(ut);
mn.OnGizmoUpdated(nodeDV, UT);