Code direct

[alesiagr]

fausse alerte, erreur trouvée, je remets les résultats d'ici peu

[jbcaillau]

@alesiagr dispo à partir de 17:30 pour un point ?

[alesiagr]

Yes! 17h30 me va :)

[alesiagr]

Code direct sans l'initial guess:

La contrainte n'est pas touché:

Code avec initial guess:

Contrainte n'est pas touchée non plus:

[alesiagr]

Join Zoom Meeting
https://univ-cotedazur.zoom.us/j/89665456865?pwd=ZcTr0bQ1irw8rpVt5rvxtqbvFUloXL.1

Meeting ID: 896 6545 6865
Passcode: 311867

[alesiagr]

J'écris plutôt ici :)

J'ai réussi à faire converger un arc au lieu de la trajectoire entière:

Ensuite j'essaie d'augmenter l'arc (= diminuer le temps initial), en prenant en compte la solution précédente comme initial guess. Le soldeur dit "Optimal solution found":

Mais par contre la trajectoire est très discretisé. Quand j'essaie d'augmenter grid_size (101 au lieu de 100), ça diverge...

Je suppose que la solution trouvée n'en est pas vraiment une? Ou est-ce que c'est juste le pb qui est très sensible aux CI? Tu penses que je devrais faire quoi? Plutôt réessayer une autre manière?

[alesiagr]

Si j'essaie de le faire progressivement (je résous un arc court, ensuite petit à petit j'augmente le grid_size, ensuite je re-résous un arc plus grand), j'obtiens ça:

[jbcaillau]

@alesiagr bravo !!! ça converge, tu as un point de départ donc maintenant ça va le faire, quitte à passer par une homotopie (sur $t_f$, les conditions aux limites, etc.)

OK pour 18:00 sur ton zoom

[alesiagr]

J'y suis! :)

[jbcaillau]

Point du 3/7/2024

@alesiagr

• essayer d'abord en enlevant la contrainte $r^2(t) \geq \varepsilon$
• comprendre la structure : la contrainte de boîte $|\beta| \leq \lambda \pi/2$ active un arc bang, qui semble activer la contrainte d'état (touch point)
• passer en 2D vrai

x = (r₁, r₂, v₁, v₂) ∈ R^4, state

• essayer de régulariser la contrainte thermique et le terme du potentiel dans le coût selon $r_1^2 + r_2^2 + \varepsilon^2$
• ploter les coûts pour initial guess + solution

[alesiagr]

Merci beaucoup @jbcaillau! Juste 2 choses qui manquent:

• rajouter le point (0,0) sur le plot pour visualiser la singularité (s'il y en a une)

• cas échéant rajouter la pénalisation dans le coût avec log

[jbcaillau]

Merci beaucoup @jbcaillau! Juste 2 choses qui manquent:

• cas échéant rajouter la pénalisation dans le coût avec log

@alesiagr yes. sur le deuxième point, c'est soit ça, soit ça la régularisation en $r^2 + \varepsilon^2$ dans contrainte et potentiel.

[alesiagr]

C'est intéressant de comparer les coûts en effet. On "perd"de l'énergie pour en gagner plus à la fin. (je sais pas trop encore à quoi c'est lié).

J'essaie tjrs de faire converger...

[jbcaillau]

@ldellelc wants to give some advice 🙂

[alesiagr]

I am all ears

[ldellelc]

energie plus petite = demi-grand axe plus petit.. peut etre que cela est fait pour aller plus proche au soleil justemment:

"reculer pour mieux sauter!" :)

Pour voir si c'est le cas essaie de plotter la norme de r en fonction du temps (pour les deux trajectoires.. mais rescale l'axe y pour mieux observer ce qu'il se passe en proximité du passage au perihelion) ainsi que l'altitude du perihelion en fonction du temps aussi

[jbcaillau]

@alesiagr

Point du 8/7/2024

• voir ici premiers exemples multi-grille
• raffiner au niveau du perihélie
• régulariser si nécessaire la singularité en $1/r$

$$\frac{1}{\sqrt{r_1^2 + r_2^2 + \varepsilon^2}}$$

par exemple avec ta propre fonction norme :

my_norm(x; eps=1e-9) = sqrt(sum(x.^2) + eps^2)

• passer en cartésiennes sur l'angle en remplaçant selon $u_1 = \cos\beta$, $u_2 = \sin\beta$, en ajoutant les contraintes $u_1^2(t) + u_2^2(t) \leq 1$ et, pour aider le solveur (voir si ça change ou non la convergence), $-[1, 1] \leq u(t) \leq [1, 1]$.

[alesiagr]

@jbcaillau Voici le code simple pour JuliaCon: https://github.com/mctao-inria/sail-therm/blob/main/interstellar_sail2D_JuliaCon.jl

1. Un petit exemple de tir sur une petite portion de l'arc
2. Exemple de continuation sur le T0 pour retrouver la trajectoire

[jbcaillau]

@jbcaillau Voici le code simple pour JuliaCon: https://github.com/mctao-inria/sail-therm/blob/main/interstellar_sail2D_JuliaCon.jl

1. Un petit exemple de tir sur une petite portion de l'arc
2. Exemple de continuation sur le T0 pour retrouver la trajectoire

@alesiagr Super, je teste et reviens vers toi.

[alesiagr]

@jbcaillau J'AI LA TRAJECTOIRE ENTIÈRE!!!!!!!

bon, c'était l'ancien run, donc avec la contrainte non convexifiée (sum(u^2) == 1), et avec la time_grid toujours uniforme.

Donc le code est très lent, mais je pense une dois la question de time_grid réglée, je pourrais le rendre plus efficace. De même pour la convexification de la contrainte.

[alesiagr]

@ldellelc @jbcaillau c'est intéressant la structure, on refait un tour supprimé pour regagner de l'énergie à la fin. De plus on remarque bien les bangs, mais aussi (probablement) des touch-points de la contrainte.

[jbcaillau]

👍🏽 @alesiagr intéressant ce tour supplémentaire. voyons si tu peux améliorer le calcul / la vitesse en passant en cartésiennes sur le contrôle, et en raffinant la grille au fur et à mesure. en tout cas, c'est forcément assez délicat car la trajectoire trouvée est très différente de l'initial guess.

[alesiagr]

Ce changement ce passe tout à la fin (quand je fais la continuation sur t0, en allant vers t0 = 0, ce qui est noté par le point vers dans la trajectoire), donc c'est possible que ça soit un min local. En tout cas je pense c'est mieux de refaire le run avec la grid non uniforme, car le pb est très sensible.

[jbcaillau]

@alesiagr quelle est la valeur de l'énergie juste avant que le tour ne soit rajouté, et juste après ? si c'est (beaucoup) moins bon juste après, c'est probablement un min local (ou un cas de bifurcation plus sérieux, mais bon)

[alesiagr]

c'est très proche. mais on perd de l'énergie, donc je me dis qu'on reste bloqué dans un local min

1.0993742095119743 juste avant le rajout du 2 arc
1.048803277238019 juste après

[alesiagr]

Mais, encore une fois, je pense ça ira mieux une fois la grid raffinée. Parce que j'ai vraiment la traj très linéarisée en passant proche du soleil, ce qui est un soucis majeur à mon avis

[alesiagr]

Hello @jbcaillau . La time-grid variable fonctionne bien, je fais les simulations. Mon inquiétude est que j'obtiens des résultats où les adjoints sont fortement oscillants. Je voulais donc te demander ce que t'en penses :

[jbcaillau]

@alesiagr une question concernant le comportement oscillatoire des approximations des adjoints (par les multiplicateurs de Lagrange), sachant que c'est vraisemblablement dû à la contrainte d'état : est-ce que c'est comme ça sur toutes tes simus ? (je me souviens plus de tes premiers plots.)

[alesiagr]

@jbcaillau tu as raison: non. C'est depuis que j'ai la time_grid non-uniforme on dirait

[alesiagr]

@jbcaillau de plus, avant j'avais bien le saut dans les adjoint exactement au même endroit que le potentiel
touch-point de la contrainte

[jbcaillau]

@alesiagr essaie stp de régénérer une simu de ce style, même avec un petit bout d'arc

[alesiagr]

@jbcaillau
Je viens de le faire: tjrs oscillatoire, meme avec time grid uniforme.

Je pense que le pb est pas le time_grid, mais plutot u = [cos β , sin β ]

Regarde une des simus que j'avais quand u = β :

[alesiagr]

Je viens de refaire les sims avec u = β : c'est encore oscillatoire....
Je ne comprends pas

[alesiagr]

OK. Donc tout à la fin de JuliaCon, juste avant que je commence à faire du grid nonuniforme, c'était pas oscillatoire, bien avec u = [cos β, sin β ]

[alesiagr]

Le soucis était que j'avais pas assez de points et tout à la fin je restais coincée dans le local minimum avec des passages au perihelion supplémentaires.

[alesiagr]

Du coup c'est juste après la mis à jour de OptimalControl.jl 😅 Est-ce possible que c'est lié?

[jbcaillau]

Du coup c'est juste après la mis à jour de OptimalControl.jl 😅 Est-ce possible que c'est lié?

A priori non :

• j'ai refait le test de Goddard, adjoints réguliers comme avant la mise à jour ¹
• y a un autre cas test avec des adjoints très oscillatoires, c'est aussi avec une contrainte sur l'état
  https://github.com/control-toolbox/BVP-DAE/blob/main/julia/goddard.ipynb

Footnotes

1. pour être précis, mise à jour de ADNLPModels.jl dont on se sert pour interfaces le solveur opti (Ipopt) ↩

[alesiagr]

Ok, je vois.
C'est étrange cette histoire 🧐

[alesiagr]

Je vois pour le cas Goddard oscillatoire... Et c'est quoi BVP-DAE?

[jbcaillau]

Je vois pour le cas Goddard oscillatoire... Et c'est quoi BVP-DAE?

Une discrétisation du problème aux deux bouts par collocation plutôt que tir, en laissant la maximisation du Hamiltonien implicite (ça te rappellera qqchose 🙂)

[alesiagr]

Ah oui, je vois :)

Bon j'essaie quand même de réfléchir sur cette histoire, mais ça me surprend que les adjoint ont décidé d'un coup de devenir oscillants...

Je vois pour le cas Goddard oscillatoire... Et c'est quoi BVP-DAE?

Une discrétisation du problème aux deux bouts par collocation plutôt que tir, en laissant la maximisation du Hamiltonien implicite (ça te rappellera qqchose 🙂)

[alesiagr]

@jbcaillau on vient de faire une réunion avec notre stagiaire Vincent, qui travaille également sur Julia avec OptimalControl.jl. Il s'avère qu'il a le même problème que moi: les adjoints sont devenus fortement oscillatoires après la mise à jour de la toolbox.

[jbcaillau]

@jbcaillau on vient de faire une réunion avec notre stagiaire Vincent, qui travaille également sur Julia avec OptimalControl.jl. Il s'avère qu'il a le même problème que moi: les adjoints sont devenus fortement oscillatoires après la mise à jour de la toolbox.

Check control-toolbox/CTDirect.jl#184

[jbcaillau]

@alesiagr ça donne quoi les nouveaux résultats ? sachant que le bug de récupération des adjoints (qui n'affectait pas le calcul de l'optimimum, c'était un post-traitement) ajoutait non seulement des oscillations, mais décalait aussi les valeurs (ce qui explique que le saut de l'adjoint ne coïncidait plus avec le supposé touch point).

[alesiagr]

@jbcaillau je n'ai pas encore fait la simulation entière, mais ça donne bien selon les résultats préliminaires!

[alesiagr]

On dirait même que c'est pas vraiment un saut dans les adjoint. En vrai si on "sauter" la partie de l'adonis à perihelion, on dirait que sa dynamique est bien continue

[jbcaillau]

heu... adonis ? tu es sûre 😅 ?

[jbcaillau]

@alesiagr essaie de bouger la vapeur de la contrainte pour voir si ça change