Suivi stage I. Mimoun

[jbcaillau]

Point du 17/06/2024

@ibtissammim

• tester control-toolbox sous Julia, notamment les tutos
• s'approprier le problème de saturation de deux spins décrit Section 2.1 (voir (2)) et 4 de ce papier :

• commencer à coder ce problème en Julia pour le résoudre par une approche directe, comme fait au 6.1 dans la référence citées, pour reproduire les résultats Figures 9 et 10

@

[jbcaillau]

See also: control-toolbox/CTDirect.jl#119

[jbcaillau]

Point du 24/06/2024

@ibtissammim

• vérifier attendus sur rapport 2ème année ENSEEIHT (option big data & HPC)
• poursuivre résolution du pb de la réf. ci-dessus
• infos sur projets Julia, voir ici
• julia + vscode + git (under wsl) set

[ibtissammim]

Bonjour Monsieur,
Est ce que les constantes du problème sont bien :
ε = 0.1
Γ = 9.855e-2
γ = 3.65e-3 ?
Cordialement,
Ibtissam

[ibtissammim]

Bonjour Monsieur,

Le problème ne converge plus. En effet, j'ai constaté que j'avais mal défini les conditions initiales. Après avoir effectué les changements nécessaires, le problème ne converge plus et j'obtiens l'erreur suivante :
EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
knot-vectors must be unique and sorted in increasing order

Je vais mettre la nouvelle version du code sur Git.

Cordialement,
Ibtissam MIMOUN

[jbcaillau]

@ibtissammim Deux exemples avec "warm start" pour continuation sur condition initiale :
https://control-toolbox.org/CTDirect.jl/stable/tutorial.html#Initial-guess-options
https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-batch.html

[ibtissammim]

D'accord, merci beaucoup.

[jbcaillau]

Point du 1er juillet 2024

@ibtissammim

• mettre au propre dans un code commun les 3 continuations utilisées (sur condition initiale)
• comparer (structure des trajectoires dans la sphère de Bloch) les différents minima locaux trouvés
• utiliser aussi l'initialisation par le système à un seul spin
• réaliser un notebook markdown du type de celui-là (voir source ici)
• résolution par tir simple du pb de saturation / contraste, cf. Anciens codes #2 ; voir notamment Section 7.1 du papier pour calcul des contrôles singuliers

Misc

• install juliaup

[jbcaillau]

Point du 3 juillet 2024

@ibtissammim

• @jbcaillau passer à la CI avec Documenter.jl
• pb contraste : code direct
• pb contraste : code indirect (voir Suivi stage I. Mimoun #1 (comment))

[ibtissammim]

@jbcaillau
J'ai un petit soucis avec la fonction de tir. En effet, je vois pas d'où vient le H01 et H1 dans la section 7.1.

[jbcaillau]

@jbcaillau J'ai un petit soucis avec la fonction de tir. En effet, je vois pas d'où vient le H01 et H1 dans la section 7.1.

$$ H_0(x,p) = \langle p,F_0(x) \rangle,\quad H_{01} = \lbrace H_0,H_1 \rbrace $$

(crochet de Poisson, voir ici)

[ibtissammim]

@jbcaillau
C'est quoi les fonctions : pz1 et pz2 de la page 35? Merci

[jbcaillau]

@ibtissammim

Point du 8/7/2024

• terminer la doc
• pour le tir multiple, utiliser fsolve de MINPACK.jl, voir ici

[ibtissammim]

@jbcaillau
Bonjour monsieur,
J'ai réussi à faire marché la fonction de tir, avec MINPACK, mais la solution ne converge pas, comment changer le p0? est ce que je le fais aléatoirement?
Merci

[jbcaillau]

bonjour @ibtissammim ; merci pour les nouvelles. pour le p0, surtout pas ! normalement, tu dois avoir une bonne initial guess avec le direct et sol.costate(0). si ça ne suffit pas à faire converger, re-résout le pb avec une grille plus fine.

[ibtissammim]

Bonjour @jbcaillau , shoot! me retourne une valeur, que signifie-t-elle ? Et est ce que la norme du vecteur s de shoot! doit être à une certaine valeur ?
La methode indirecte converge mais en changeant la tolérance de fsolve qui est normalement $10^{-8}$ à $10^{-6}$

[jbcaillau]

Bonjour @jbcaillau , shoot! me retourne une valeur, que signifie-t-elle ? Et est ce que la norme du vecteur s de shoot! doit être à une certaine valeur ? La methode indirecte converge mais en changeant la tolérance de fsolve qui est normalement 10−8 à 10−6

@ibtissammim shoot! retourne la dernière valeur affectée, ici une composante du tir :

spin/src/tir_saturation.jl

Line 75 in f7def56

s[8] = (pf[2] + pf[4]) * γ - 1

C'est plus propre d'ajouter à la fin un return nothing.

OK pour les tolérances, c'est normal (on reverra ensemble).

[ibtissammim]

Bonjour @jbcaillau, Lorsque j'essaye de déployer la page web de tir_saturation.md que je viens d'écrire exactement comme bisaturation.md, j'ai cette erreur :

Alors que j'ai effectué des changements dans documentation.yml :

Et aussi des modifications dans make.jl:

Est que je dois changer un autre fichier ?
Merci

[jbcaillau]

Bonjour @ibtissammim il ne faut pas faire comme ça dans Documentation.yml : fais ton .md sans te poser de question, je fais les modifs après pour rendre le déploiement possible.

[ibtissammim]

D'accord @jbcaillau, merci.

[ibtissammim]

Bonjour @jbcaillau @ocots,

J'ai un petit souci avec l'utilisation du solveur linéaire MA57 avec HSL. Les valeurs du contrôle sont l'opposé de ce que je dois trouver. Je ne vois pas d'où vient cette erreur.

Merci!!

[ocots]

Ce n'est pas nécessairement une erreur car c'est symétrique. C'est tout de même légèrement surprenant que le solveur converge vers la solution symétrique.

[jbcaillau]

@ocots en effet, symétrie $(y, z, u) \mapsto (-y, z, -u)$

[jbcaillau]

@ibtissammim je viens de faire des modifs dans la doc (et dans les fichiers CI / doc), tout passe et c'est mergé dans la branche main #9

Par contre, la doc tir_saturation.md est à compléter (j'ai supprimé la fin qui ne passait pas, et semblait obsolète

NB. Il faut visiblement mettre des global pour les variables qui apparaissent dans les boucles (alors que ce n'est plus la peine en julia standard¹ / hors Documenter)

Footnotes

1. https://docs.julialang.org/en/v1/manual/variables-and-scoping/ ↩

[ibtissammim]

D'accord @jbcaillau

[jbcaillau]

@ibtissammim pour la suite, voir :

• le code BVP sur cet exemple cet exemple (Van der Pol, et tester d'abord la version primale, a priori moins précise mais plus simple)
• le papier qui explique la méthode

[ibtissammim]

Bonjour @jbcaillau,
A quoi correspond le paramètre mu? puisque je le vois pas définit dans la fonction init. Je suppose que cette fonction définit la dynamique, et donc puisque dans la dynamique du problème j'ai besoin de \epsilon, cette variable "z" doit le contenir. est ce que c'est correct?

Voici ce comment j'ai définit la fonction ode :

[ibtissammim]

Bonjour @jbcaillau, je suis bloquée et je ne comprends pas comment notre problème de spin peut être traité de la même manière que les autres problèmes.

[jbcaillau]

bonjour @ibtissammim rdv 16:30 ce mercredi pour un point zoom :
https://univ-cotedazur.zoom.us/my/caillau

Bonjour @jbcaillau, je suis bloquée et je ne comprends pas comment notre problème de spin peut être traité de la même manière que les autres problèmes.

[ibtissammim]

D'accord @jbcaillau, Merci :)

[jbcaillau]

@ibtissammim control-toolbox/BVP-DAE#1 (comment)

[PlMlsn]

@ibtissammim @jbcaillau Si je peux me permettre, je pense que le plus simple est de partir du problème Goddard primal-dual. D'une part c'est un problème à temps minimal comme celui que vous traitez. Par ailleurs, en réalité, la méthode primale-duale est assez simple à implémenter car elle est très systématique.

[PlMlsn]

Bonjour @jbcaillau, A quoi correspond le paramètre mu? puisque je le vois pas définit dans la fonction init. Je suppose que cette fonction définit la dynamique, et donc puisque dans la dynamique du problème j'ai besoin de \epsilon, cette variable "z" doit le contenir. est ce que c'est correct? Voici ce comment j'ai définit la fonction ode :

En général, je note z la variable contenant les solutions des équations algébriques. Plus précisément, il s'agit du contrôle u et de tous les multiplicateurs des contraintes. Dans ce cas, la variable mu est le multiplicateur associé à la contrainte d'état.

[jbcaillau]

Bonjour @jbcaillau, A quoi correspond le paramètre mu? puisque je le vois pas définit dans la fonction init. Je suppose que cette fonction définit la dynamique, et donc puisque dans la dynamique du problème j'ai besoin de \epsilon, cette variable "z" doit le contenir. est ce que c'est correct? Voici ce comment j'ai définit la fonction ode :

En général, je note z la variable contenant les solutions des équations algébriques. Plus précisément, il s'agit du contrôle u et de tous les multiplicateurs des contraintes. Dans ce cas, la variable mu est le multiplicateur associé à la contrainte d'état.

🙏🏽 merci @PlMlsn en effet, on a vu ça hier : $\mu$ pour $g(x) \leq 0$, $\eta$ pour $c(x,u) \leq 0$, de mémoire

[jbcaillau]

@ibtissammim @jbcaillau Si je peux me permettre, je pense que le plus simple est de partir du problème Goddard primal-dual. D'une part c'est un problème à temps minimal comme celui que vous traitez. Par ailleurs, en réalité, la méthode primale-duale est assez simple à implémenter car elle est très systématique.

@PlMlsn dans l'immédiat, il s'agit de tester ton code python sur le cas de contrôle quantique (pas de réimplémenter). merci pour le conseil !

[ibtissammim]

Bonjour @jbcaillau, A quoi correspond le paramètre mu? puisque je le vois pas définit dans la fonction init. Je suppose que cette fonction définit la dynamique, et donc puisque dans la dynamique du problème j'ai besoin de \epsilon, cette variable "z" doit le contenir. est ce que c'est correct? Voici ce comment j'ai définit la fonction ode :

En général, je note z la variable contenant les solutions des équations algébriques. Plus précisément, il s'agit du contrôle u et de tous les multiplicateurs des contraintes. Dans ce cas, la variable mu est le multiplicateur associé à la contrainte d'état.

Merci beaucoup pour votre réponse @PlMlsn.

[ibtissammim]

Bonjour @jbcaillau,

J'espère que vous allez bien. J'ai une question concernant l'initialisation du co-état. Dans le problème de Zermelo, qui est un problème de temps minimal, nous devons écrire une fonction initialize qui initialise les états et les co-états. Cependant, je n'ai aucune idée de la manière d'initialiser pp dans le cas du problème de spin.

Merci beaucoup.

[jbcaillau]

bonjour @ibtissammim ; tu dois pouvoir réutiliser la solution que tu as précédemment trouvée soit par le code direct, soit par le code indirect (tir). Par exemple, avec la solution du direct (sol = solve(ocp; ...)) tu récupères l'état, sol.state, et l'adjoint (co-état), sol.costate. Dans les deux cas, il s'agit de fonctions du temps que tu peux évaluer en n'importe quel t. Voir par ex. ce use case.

Je te propose un point visio demain (vendredi) 11:00, dis moi stp si OK pour toi :
https://univ-cotedazur.zoom.us/my/caillau

[ibtissammim]

Bonjour @jbcaillau, je suis désolée, je viens de voir votre message. Oui c'est ok pour moi.
Merci

[jbcaillau]

Point du 16/08/2024

• essayer un Zermelo sans contrainte d'obstacle (= sans contrainte sur l'état), encore plus semblable au pb quantique
• pour l'init, récupérer les valeurs du code Julia, par ex. en rescalant tout sur $[0,1]$ (scaling à faire sur $p$, je pense) :

ocp1 = @def begin
    s ∈ [0, 1], time
    y = (x1, x2, x3, x4, tf) ∈ R^5, state
    u ∈ R, control
    tf(s) ≥ 0
    -1 ≤ u(s) ≤ 1
    x(0) == [0, 1, 0, 1]
    x(1) == [0, 0, 0, 0]

    derivative(y)(s) == tf(s) * [(-Γ*x1(s) -u(s)*x2(s)),
             (γ*(1-x2(s)) +u(s)*x1(s)),
             (-Γ*x3(s) -(1-ϵ)* u(s)*x4(s)),
             (γ*(1-x4(s)) +(1-ϵ)*u(s)*x3(s)),
             0 * tf(s) ]
    tf(1) → min
end

• ajouter une init. linéaire pour $p_{t_f}$ allant de $0$ à $-1$ (convention max du hamiltonien)

[ibtissammim]

Bonjour @jbcaillau,
J'espère que vous allez bien. Je suis entrain de résoudre le problème ci-dessus, et j'ai une question : Comment dois-je initialiser tf ? est ce que tf(0) = tf que j'ai déjà obtenu avant ? parce que la solution de ce problème est très loin de ce qu'on doit obtenir normalement.
Merci

[ibtissammim]

Bonjour @jbcaillau,

J'espère que vous allez bien. J'ai réussi à résoudre le problème précédent. Cependant, j'ai une nouvelle question concernant le code Python. J'ai pu corriger le problème de dimension, mais maintenant le code ne converge plus. Pour les conditions initiales, je les ai définies en modifiant les 4 co-états, mais je ne suis pas certain que ce soit correct.

Merci.

[PlMlsn]

Bonjour @jbcaillau,

J'espère que vous allez bien. J'ai réussi à résoudre le problème précédent. Cependant, j'ai une nouvelle question concernant le code Python. J'ai pu corriger le problème de dimension, mais maintenant le code ne converge plus. Pour les conditions initiales, je les ai définies en modifiant les 4 co-états, mais je ne suis pas certain que ce soit correct.

Merci.

Si vous voulez que je jette un œil au code python pour voir ce qui ne marche pas, c'est tout à fait possible.

[ibtissammim]

Bonjour @PlMlsn,
Merci beaucoup pour votre réponse. Ce code Python se trouve dans le fichier BiSaturation.py, situé dans le dossier SpinBiSaturation de ce repo.

[PlMlsn]

Bonjour @PlMlsn, Merci beaucoup pour votre réponse. Ce code Python se trouve dans le fichier BiSaturation.py, situé dans le dossier SpinBiSaturation de ce repo.

J'ai parcouru le fichier et j'ai corrigé les erreurs de calcul de Jacobien. C'est un des points faibles de la méthode, cela nécessite de faire pas mal de calculs fastidieux.

Je ne sais pas si le code va désormais mieux converger car le problème a l'air très complexe mais je pense qu'il n'y a plus d'erreur dans les formules de dérivées.

Paul

[ibtissammim]

Je vais vérifier si cela fonctionne. Merci beaucoup pour votre aide, @PlMlsn.

[ibtissammim]

Bonjour @PlMlsn,
Il y a eu un petit changement par rapport à ce que j'avais auparavant. Cette fois-ci, les erreurs résiduelles s'affichent, mais après le message "Solving Complete", le code ne converge plus.
Cela pourrait-il être dû à l'initialisation du problème ? J'ai initialisé p5 à 0 et p1 à p4 avec linspace(0, -1, n).
Merci.

[PlMlsn]

Bonjour @PlMlsn, Il y a eu un petit changement par rapport à ce que j'avais auparavant. Cette fois-ci, les erreurs résiduelles s'affichent, mais après le message "Solving Complete", le code ne converge plus. Cela pourrait-il être dû à l'initialisation du problème ? J'ai initialisé p5 à 0 et p1 à p4 avec linspace(0, -1, n). Merci.

En fait quand l'erreur résiduelle augmente comme cela, c'est que le schéma de Newton pour résoudre les équations de collocations ne converge pas. Il faut regarder ce qu'il y a dans infos voir si la NLSE est bien proche de zéro. Si ce n'est pas le cas c'est que l'initialisation n'est pas suffisamment bonne.

[ibtissammim]

Bonjour @PlMlsn,
Merci beaucoup pour votre aide. Les résidus infos.NLSE_infos.bc_residual et infos.NLSE_infos.ode_residual dans infos sont proches de zéro. Y a-t-il d'autres variables à vérifier pour comprendre pourquoi le problème ne converge pas ?

[PlMlsn]

Bonjour @PlMlsn, Merci beaucoup pour votre aide. Les résidus infos.NLSE_infos.bc_residual et infos.NLSE_infos.ode_residual dans infos sont proches de zéro. Y a-t-il d'autres variables à vérifier pour comprendre pourquoi le problème ne converge pas ?

J'ai modifié l'objet infos pour afficher le max des résidus, en testant je vois que les erreurs ode_residual et ae_residual sont très loin de zéro. C'est bien le schéma de Newton qui ne converge pas.

[ibtissammim]

Point du 16/08/2024

* essayer un Zermelo sans contrainte d'obstacle (= sans contrainte sur l'état), encore plus semblable au pb quantique
* pour l'init, récupérer les valeurs du code Julia, par ex. en rescalant tout sur [0,1](scaling à faire sur p, je pense) :

ocp1 = @def begin
    s ∈ [0, 1], time
    y = (x1, x2, x3, x4, tf) ∈ R^5, state
    u ∈ R, control
    tf(s) ≥ 0
    -1 ≤ u(s) ≤ 1
    x(0) == [0, 1, 0, 1]
    x(1) == [0, 0, 0, 0]

    derivative(y)(s) == tf(s) * [(-Γ*x1(s) -u(s)*x2(s)),
             (γ*(1-x2(s)) +u(s)*x1(s)),
             (-Γ*x3(s) -(1-ϵ)* u(s)*x4(s)),
             (γ*(1-x4(s)) +(1-ϵ)*u(s)*x3(s)),
             0 * tf(s) ]
    tf(1) → min
end

* ajouter une init. linéaire pour ptf allant de 0 à  −1
   (convention max du hamiltonien)

Je n'ai pas bien compris cette indication : "ajouter une init. linéaire pour p_{tf} allant de 0 à − 1 " . Est ce que cela veut signifier que p1, ..p5 va de 0 à -1 ?

[jbcaillau]

@ibtissammim je passe au bureau vers 16:30. on fait un point rapide si tu es dispo.

[jbcaillau]

@ibtissammim et point final jeudi 16:00, a priori en visio :

https://univ-cotedazur.zoom.us/my/caillau

[ibtissammim]

Bonjour @jbcaillau,
J'espère que vous allez bien. Le problème ne converge toujours pas, même après avoir changé les conditions initiales (par contre l'erreur a diminué légèrement pour n = 501).

[jbcaillau]

@ibtissammim OK : point 16:00 sur zoom. à tout
https://univ-cotedazur.zoom.us/my/caillau

[ibtissammim]

Bonjour @PlMlsn,
J'espère que vous allez bien. Pour la fonction 'initialize', j'ai initialisé p, q et u avec les valeurs obtenues à partir de mon code Julia 'bisat_scaling.jl'. Bien que l'erreur pour n = 501 ait légèrement diminué, le problème ne converge toujours pas. Pourriez-vous m'aider à comprendre d'où vient cette erreur ?
Merci d'avance pour votre aide

[ibtissammim]

Bonjour @jbcaillau,

Voici un petit bilan de ce que j'ai fait pendant mon stage :

• Résolution du problème des deux spins par la méthode directe, ce qui a permis de trouver le minimum global ainsi que des minima locaux.

• Résolution du problème des deux spins par la méthode indirecte.

• Résolution du problème des deux spins avec un code Python.

Améliorations/points à traiter :

• Faire fonctionner le code Python (le problème ne converge pas pour le moment).

• Raffiner la solution obtenue avec la méthode indirecte (je n'ai rien changé dans le code, et la dernière fois que je l'ai exécuté, la solution était optimale, mais ce n'est plus le cas actuellement).

[jbcaillau]

@ibtissammim j’archive ça : merci pour le travail ! bon retour à toulouse 🙂