\documentclass[a4paper,12pt]{article}
\usepackage{latexsym}
\usepackage{amsmath}
% amsfonts permet le \mathbb(C) attantion, peu de lettres disponibles
\usepackage{amsfonts}
\usepackage{eepic,psfig,pstricks,multirow}
\usepackage{floatflt}
\usepackage{fancyheadings}
\usepackage{footnpag}
\usepackage{rotating}
\usepackage{euscript}
\usepackage{tabularx,xspace}
%\usepackage{minitoc}
\usepackage[frenchb]{babel}
\usepackage[latin1]{inputenc}

\topmargin -1cm
\addtolength{\textheight}{4cm}
\tolerance=8000
\pretolerance=4000
\textwidth 15cm
\oddsidemargin 0.5cm

%
\newcommand{\bra}[1]{<{#1}\mid}
\newcommand{\ket}[1]{\mid{#1}>}
\newcommand{\mat}[1]{\underline{\underline{\mathcal{#1}}}}
\newcommand{\diag}[1]{{#1}^{d}}
\newcommand{\mean}[1]{\overline{{#1}}}
\newcommand{\scal}[2]{{<\, #1 \mid #2>}}
\newcommand{\ext}[2]{{\mid\, #1 >< #2\mid}}
\newcommand{\scalm}[3]{{<\, #1 {#2}^{\dagger} \, #2 #3>}}
\newcommand{\qrayl}[3]{\frac{\scalm{#1}{#2}{#3}}{\scal{#1}{#3}}}
\newcommand{\trac}[1]{{\mathcal{T}r (#1)}}


\begin{document}
%
\def\todec{{\overline {\overline {\Upsilon  }}}_{,dec} }
\def\tox{{\overline {\overline {\Upsilon  }}}_{,x} }
\def\dt{{\delta t}}
\def\edtaux{{exp (-{\delta t  \over \taux } ) }}
\def\deta{{\delta \eta } }
\def\Deta{{\Delta \eta } }
\def\Dt{{\delta t } }
\def\dfi{{\delta \varphi } }
\def\Dfi{{\Delta \varphi } }
\def\dts{{\delta t\over 2 } }
\def\dtaf{{\overline {\overline {\partial _\eta f }}} }
\def\d(#1)G{{\overline {\overline {\partial_#1 G}}} }
\def\dA{{\overline {\overline {\partial A}}} }
\def\dtA{\partial_t A}
\def\VD{{\mathcal V \mathcal D}}
%     \xspace est un "gentil" espace qui s'ajuste
\newcommand{\cad}{c'est-\`a-dire\xspace}
\newcommand{\Minik}{Mini$\_$ker\xspace}
% essai pour subsection a la JAM, cad italique a. titre
\makeatletter
% d'abord, numero = alphabetique
\renewcommand{\thesubsubsection}{\alph{subsubsection}.}
% style de la forme subsection
%\renewcommand{\subsubsection}{\@startsubsection
% {subsubsection}%                      %the name
% {2}%                               % the level
% {0mm}%                             % the indent
% {-\baselineskip}%                  % the before skip
% {0.5\baselineskip}%                % the after skip
% {\normalfont\large\itshape}}% % the style
%%  {\normalfont\normalsize\itshape}}% % the style
\makeatother

\title{\Minik, adjoint et sensibilités, un résumé\footnote{Al1, Décembre 2003; révisions
Juillet 2004 et Janvier 2005; section adjoint-propagateur et quadruple comparaison, de Juin 2005}\date{Juin 2005}}
\maketitle
%
{\it{\bf Avertissement} Ce texte prend appui sur les articles fondateurs de Cacuci, en
omettant ses aspects de rigueur mathématique, et le lecteur est prié de se rapporter aux
articles originaux pour ce qui concerne les ensembles supports et la possibilité de définir
l'opérateur adjoint. On a aussi simplifié les notations en perdant la distinction entre les
divers produits scalaires nécessaires, car il est facile d'en remonter ls signification
suivant les opérandes. On a donc indistinctement utilisé la
notation ``bra-ket'' de Dirac. Par rapport à Cacuci, on a étendu son
vecteur de paramètres aux transferts du TEF, et on a introduit un vecteur explicite de commande
pour rendre les applications plus souples.}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le problème et sa résolution par l'adjoint}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
On a un modèle écrit sous TEF, dépendant d'une loi de commande $h(t)$, et une trajectoire
$\mathcal{T}$, ayant - une fois - initialisé une loi de commande. On a un critère à maximiser :
%
\begin{equation}
J = \psi(\eta(T) ,h(T)) + \int_0 ^T {l(\eta(\tau),h(\tau)}\, d\tau
\end{equation}
%
A une trajectoire proche, provoquée par exemple par une modification $u(t)$ de $h$,
correspond une variation de la fonctionnelle $J$ :
%
\begin{equation}
dJ = \psi_x(\eta(T) ,h(T))\, x(T) + \int_0 ^T \{l_x(\eta(\tau),h(\tau))\, x(\tau) +
l_u(\eta(\tau),h(\tau))\, u(\tau) \}\, d\tau
\label{dj}
\end{equation}
%
on a noté $x$ la variation $\Deta$ de $\eta$ et on n'a pas précisé les produits scalaires
implicites. Ce qu'on veut, c'est éliminer $x$ de cette variation, de manière à ce qu'on
obtienne $\bra{d_u J}$, le gradient qui permet d'attaquer les algorithmes divers de descente
pour résoudre notre problème de commande optimale.

Ce problème est traité classiquement par maximisation de $J$ sous la contrainte des équations
du SLTC, avec l'introduction d'un vecteur de Lagrange $\lambda$, et du Hamiltonien du système.
Nous suivons ici l'approche de Cacucci, d'abord sans introduire explicitement les transferts,
pour ne pas alourdir la présentation. Dans la section~\ref{stef}, on étendra au TEF.

Le système d'équations du SLTC est :
%
\begin{eqnarray}
\partial_t x(t) &=&  A(t)\, x(t) + B(t)\, u(t)
\label{sltc}
\end{eqnarray}
%
On choisit un produit scalaire - qui en général comprendra une intégration temporelle -
et on introduit une variable adjointe $v(t)$, vecteur ligne, par laquelle on multiplie
scalairement~(\ref{sltc}) en isolant la partie commande à droite de l'égalité :
%
\begin{eqnarray}
\scal{v(t)}{\partial_t x(t)} - \bra{v(t)} A(t) \ket{x(t)} &=& \bra{v(t)}B(t) \ket{u(t)}
\label{scal}
\end{eqnarray}
%
{\it Précision mathématique\footnote{cette partie a été corrigée par stepH.} :} le
système~(\ref{scal}) est équivalent à~(\ref{sltc})
lorsqu'il est pris $\forall v(t)$ (Lax-Milgram). Ainsi, pour toute fonction $v(t)$,
l'égalité~(\ref{scal}) est vérifiée. On est donc en mesure de définir un nouveau problème
dont $v(t)$ est solution : ce système restera indépendant de $x(t)$
%
\begin{eqnarray}
-\scal{\partial_t v(t)}{x(t)} - \bra{v(t)} A(t) \ket{x(t)} &=& \scal{l_x(\eta(t),h(t))}{x(t)}
\label{scaladj}
\end{eqnarray}
%
cette fois ci pour tout $x$. Ce problème définit le système adjoint à~(\ref{sltc}) pour
la fonction de coût $J$. 

En effectuant une intégration par partie dans le premier produit scalaire\footnote{c'est là
que l'on restreind la généralité de la vision de Cacuci qui envisage des produits scalaires
plus généraux.}, on obtient :
%
\begin{eqnarray}
\scal{v(t)}{\partial_t x(t)} &-& \bra{v(t)} A(t) \ket{x(t)} \nonumber\\
 = \scal{l_x(\eta(t)}{h(t)),x(t)} &-& \scal{v(T)}{x(T)} + \scal{v(0)}{x(0)}
\label{scaladjb}
\end{eqnarray}
%
pour tout $v$, ce qui permet d'identifier le terme de droite de~(\ref{scaladj}) avec celui
de~(\ref{scaladjb}). De cette façon, on élimine le terme dépendant de $x(t)$ dans la variation
$dJ$ par une expression dépendant de $v(t)$ obtenue indépendamment de $x(t)$. Il suffit de
résoudre le système adjoint suivant :
%
\begin{eqnarray}
- \bra{\partial_t v(t)} &=& \bra{v(t)} A(t) + \bra{l_x(\eta(t),h(t))}
\label{adj}
\end{eqnarray}
%
les termes de borne de l'intégration par parties sont choisis en fonction de la forme de $dJ$
obtenue par remplacement dans l'intégrande du terme dépendant de $x$ :
%
\begin{eqnarray}
dJ &=& \scal{\psi_x(\eta(T) ,h(T))}{ x(T)} - \scal{v(T)}{x(T)} \nonumber\\
&\,&+ \int_0 ^T \{\bra{v(\tau)} B(\tau) \ket{u(\tau)} + \scal{l_h(\eta(\tau),h(\tau))}{ u(\tau)} \} d\tau \nonumber\\
&\,& + \scal{v(0)}{x(0)}
\label{djhtef}\end{eqnarray}
%
ainsi, en prenant comme conditions terminales $v(T) = \psi_x(T)$, on élimine la dépendance
en $x(T)$. $v(0)$ est déterminé par la résolution de l'adjoint, et $x(0)$ est connu.

In fine, on obtient le gradient voulu, qui ne dépend que de la correction $u(t)$ à appliquer
à la loi de commande, soit :
%
\begin{equation}
dJ = \scal{v(0)}{x(0)}+\int_0 ^T \{\bra{v(\tau)} B(\tau) \ket{u(\tau)}
       + \scal{l_h[\eta(\tau),h(\tau)]}{\, u(\tau)} \}\, d\tau
\label{duj}
\end{equation}
d'où on déduit le gradient recherché :
%
\begin{equation}
\nabla_u J(\tau) = \bra{v(\tau)} B(\tau) + \bra{l_h[\eta(\tau),h(\tau)]}
\label{gruj}
\end{equation}
Pour ceux qui n'auraient pas pris le temps de tout suivre, cette forme est en effet pertinente
en ce que les variables adjointes répercutent tout l'effet du changement de la commande sur
le segment de
trajectoire postérieur à $t$ (de $t$ à $T$), c'est le fondement du théorème du maximum de Pontryagine.

Le théorème du {\bf maximum de Pontryagine} dit en effet que pour la trajectoire optimale, pour tout $t$,
le Hamiltonien du système est maximum en la commande $h$ pour $\eta$ et $v$
déterminés par le système aller et retour\footnote{le cas où l'instant terminal $T$
est aussi à optimiser revient à ajouter la condition $H(T) = 0$, cf annexe} :
%
\begin{equation}
H(\eta,v,h) = \scal{v\,}{ g(\eta ,h)} + l(\eta,h)
\end{equation}
%
à chaque instant, ce qui peut être utile pour déterminer $u(t)$, par exemple avec une méthode
de « programmation dynamique ».
Globalement donc, avec l'adjoint, on a simplifié un problème de recherche d'extremum. Au départ,
on travaille dans un espace de dimension celle du vecteur commande fois le nombre de pas
d'intégration ({\it npas}), puisque la variation de $u(t)$ à l'instant $t$ se répercute
sur les instants ultérieurs.
Après passage de l'adjoint, on travaille indépendamment dans {\it npas} espaces, car alors,
en supposant que l'on soit proche de l'optimum, on peut optimiser $u(t)$ dans l'intervalle
$[t,t+\delta t]$ sans remettre en question l'optimalité de $u(t)$ en dehors de cet
intervalle.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Extension au TEF}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{stef}
Avec des signes relativement non conformes aux divers formulaires éven\-tuellement
entre eux contradictoires, le SLTC à la TEF peut s'écrire\footnote{ je ne trouve
pas de raison pour ne pas se restreindre à $B^u$ et $D^u$ diagonales, car on a introduit
ce vecteur pour simplifier l'utilisation de \Minik pour la commande optimale. Les cas plus
complexes peuvent être traités par l'intermédiaires de commandes-transferts.} :
%
\begin{equation}
\begin{bmatrix} \partial_t \!\cdot\, -A  & -B \\ -C^{\dagger} & 1-
D\end{bmatrix} \quad 
\begin{bmatrix} \Deta  \\ \Dfi \end{bmatrix} \quad = \quad
\begin{bmatrix} B^u u^{\eta}\\ D^u u^{\varphi} \end{bmatrix}
\label{sltc1} 
\end{equation}
%
Les matrices prenant leur définition la plus simple des Jacobiennes du TEF, \cad que si :
%
\begin{equation}
\left \lbrace
\begin{array}{l}
\partial_t \eta (t) =  g(\eta(t),\varphi(t),h^{\eta}(t)) \nonumber\\
\varphi(t) = f(\eta(t),\varphi(t),h^{\varphi}(t))
\end{array}
\right.
\label{tefbas}
\end{equation}
%
on a défini :
\begin{eqnarray}
A = \partial_{\eta} g &;& B = \partial_{\varphi} g; \nonumber\\
C^{\dagger} = \partial_{\eta} f &;& D = \partial_{\varphi}f ;\nonumber\\
B^h = \partial_{h} g &;& D^h = \partial_{h} f
\label{tefmat}
\end{eqnarray}
%
Le système direct sous TEF correspondant est ainsi :
%
\begin{eqnarray}
\partial_t \deta (\tau) &=& g_0 (t) \, + A(t)\, \deta (\tau) + B(t)\, \dfi (\tau)+ B^h (t)\, \delta h^{\eta}(\tau) \\
\dfi (\tau) &=& C^{\dagger}(t)\,\deta (\tau)\, + D (t)\,\dfi(\tau)\, + D^h (t) \delta h^{\varphi} (\tau)
\label{tef1}
\end{eqnarray}
%
avec $t\leq \tau < t+\delta t$, $\delta t , \delta t$ étant le pas de temps local
de l'intégration numérique.
et sa forme matricielle intégrée pour la trajectoire est :
%
\begin{equation}
\begin{bmatrix} I -\dts A  & -\dts B \\ -C^{\dagger} & 1-
D\end{bmatrix} \quad 
\begin{bmatrix} \deta  \\ \dfi \end{bmatrix} \quad = \quad
\begin{bmatrix} \dt g_0 + B^u u^{\eta}(t+\dts)\\ D^u u^{\varphi}(t+\dts) \end{bmatrix}
\label{tefint} 
\end{equation}
%
C'est la forme que j'ai utilisée pour le \Minik\footnote{contrairement à la version originelle
dans laquelle on avait une
distribution des signes et des $\dt$ un peu aléatoire.}, avec stockage des Jacobiennes,
parce qu'en pratique on en a de multiples usages (calcul de la matrice d'avance d'état,
analyse en Borel,\dots).

En ce qui concerne le SLTC et son adjoint,
on y a étendu l'application de la commande aux transferts par souci de généralité, que
l'on invoque aussi dans la fonction de coût :
%
\begin{equation}
J = \psi[\eta(T),\varphi(T) ,h(T)] + \int_0 ^T {l[\eta(\tau),\varphi(T),h(\tau)]}\, d\tau
\end{equation}
%
où la fonction vectorielle $h$ englobe les commandes sur les deux types de variables du TEF.
Le modèle adjoint (à~\ref{sltc1}) s'écrit donc de la manière suivante :
%
\begin{equation}
\begin{bmatrix} -\partial_t \!\cdot\, -A^{\dagger}  & -C \\ -B^{\dagger} & 1-D^{\dagger}\end{bmatrix} \quad 
\begin{bmatrix} v  \\ w \end{bmatrix} \quad = \quad
\begin{bmatrix} l_x\\ l_y \end{bmatrix}
\label{fadj} 
\end{equation}
%
avec comme conditions initiales : $v(T) = \psi_x(T)\quad \{ w(T)$ est déterminée par $v(T)\}$.
Ce qui fait que le système adjoint est par définition linéarisé.

Le gradient d'intérêt se décompose en deux parties $\,\nabla_{ux} J \, ,\, \nabla_{uy} J$ :
%
\begin{eqnarray}
dJ &=& \scal{v(0)}{x(0)} \nonumber\\
    &+&\int_0 ^T \{\bra{v(\tau)} B^u (\tau) \ket{u^x (\tau)}
    + \scal{l_{ux}(\eta(\tau),h^x(\tau))}{ u^x(\tau)} \}\, d\tau \nonumber \\
    &+&\!\!\!\!\int_0 ^T \{\bra{w(\tau)} D^u (\tau) \ket{u^y(\tau)}
       + \scal{l_{uy}(\eta(\tau),h^y(\tau))}{ u^y(\tau)} \}\, d\tau
\label{duuj}
\end{eqnarray}
avec des notations évidentes. Et ainsi finalement, on obtient :
\begin{eqnarray}
\nabla_{ux}J(\tau,T) &=& \bra{v(\tau)} B^u (\tau) \quad\,+\, \bra{l_{ux}(\eta(\tau),h^x(\tau))} \nonumber \\
\nabla_{uy}J(\tau,T) &=& \bra{w(\tau)} D^u (\tau) \quad+\, \bra{l_{uy}(\eta(\tau),h^y(\tau))}
\label{duuj2}
\end{eqnarray}
Le Hamiltonien correspondant peut être défini en ne conservant que les termes dépendants de $h$ :
%
\begin{eqnarray}
H(\eta,\varphi,v,w,h^x,h^y) &=& \scal{v\,}{ g(\eta,\varphi ,h^x)} \nonumber\\ 
+\quad \scal{w\,}{ f(\eta,\varphi ,h^y)} &+& l(\eta,\varphi,h)
\end{eqnarray}
à chaque instant.
%
%\vspace{0.2cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Intégration du système adjoint}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{adjf}
Le système~(\ref{fadj}) s'intègre comme pour le direct de $t$ à $t+\dt$ pour donner $\delta v$,
avec la source $l_x$ constante, comme pour les matrices, ce qui donne :
%
\begin{equation}
\begin{bmatrix} -I -\dts A^{\dagger}  & -\dts C \\ \qquad -B^{\dagger} & 1-D^{\dagger}\end{bmatrix} \quad 
\begin{bmatrix} \delta v  \\ \delta w \end{bmatrix} \quad = \quad
\begin{bmatrix} \dt ( l_x + A^{\dagger}v + C w)\\ l_y \end{bmatrix}
\label{fadjint} 
\end{equation}
\footnote{Le cas du 0 de la 2\textsuperscript{ième} ligne et du lien avec le pb du dfdpi a été résolu et amène à la définition
du gain non-$\mathcal{L}$, cf papier Propagateurs etc.}
Le plus simple consiste à prendre $\dt$ négatif pour obtenir l'incrément $\delta v$ vers
l'arrière. On résout ainsi ce système simplement par le même appel à la routine {\bf ker}
du \Minik en inversant les arguments $B$ et $C^{\dagger}$ et avec toutes les matrices
transposées; de même, le vecteur {\bf gamma}\footnote{pardon pour les lecteurs étrangers
au TEF, $\Gamma$ est le vecteur connu à droite de l'équation de cellule du système~(\ref{tefint}).}
 doit être remplacé. Ici, le remplaçant de $g_0$ est sous
forme linéaire $= A^{\dagger} v + C w$, et on a un terme source $l_x$.
On a besoin, comme pour le direct, de l'équivalent de l'équation des transferts qui n'est autre
que la deuxième ligne de~(\ref{fadj}) :
\begin{equation}
(1-D^{\dagger}) \ket{w} = B^{\dagger} \ket{v} +\, l_y
\label{fadjw} 
\end{equation}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ ... et \Minik}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{mini}
\Minik constitue une maquette logicielle destinée à développer et éprouver les méthodes exploitant
le SLTC rendues accessibles en simulation par l'utilisation du TEF.
La refondation du code originel établi par stb est basé sur la concision : ne pas avoir à écrire
chaque formule utile et dérivable plus d'une seule fois.
D'où la programmation suivante, qui se veut simple et de bon goût~: entre deux mises à jour de la
commande, on fait systématiquement un aller-retour dans \Minik.
\vspace{0.1cm}

{\bf En avant,} $u(t)$ est donnée (loi ou fichier). On calcule la trajectoire de manière
classique, en ajoutant comme diagnostic un calcul de $J$, pour vérifier l'approche
à l'optimum\footnote{cet aller permet aussi éventuellement de calculer des sensibilités en avant.}.
On stocke $\eta,\varphi,h$. Par rapport à une trajectoire standard, on a donc
des  $f\_set$ supplémentaires\footnote{cf manuels \Minik.}, pour la densité de coût $l(t)$,
et $\Psi(T)$.
\vspace{0.2cm}

{\bf Au retour,} à partir de T, initialiser $v$, remplacer le calcul
aller par restauration de $\eta,\varphi,h$ et déterminer l'incrément temporel (qui a pu varier en cours
de simulation aller). Même calcul des matrices\footnote{en restaurant $\varphi(t)$,
on n'a plus besoin d'itérer sur la D-loop.}, vecteurs donnés par~(\ref{fadjint}). Calcul de $w$
par résolution de~(\ref{fadjw}).
Appel de ker (le même qu'à l'aller) mais avec $\dt$ négatif, calcul et stockage des deux
vecteurs $d_u J(t)$ donnés par les formules~(\ref{duuj2}) et du Hamiltonien (cf fichiers
générés ---.data).
\vspace{0.2cm}

{\bf Phase hors \Minik - } il reste à corriger la commande $u(t)$ : à vos souhaits !
Il est clair que chacun aura à se faire la main sur des exos différents. Il y a le cas
important des coûts quadratiques, qu'il faudra développer. Voir quels algorithmes de
descente (par gradients conjugués ou quasi lorsqu'on peut développer la fonction de coût
au deuxième ordre seulement),
% ou via équation de Ricatti\footnote{c'est juste
% pour vous en mettre plein la vue et pouvoir utiliser les dérivées partielles de matrices qu'on a
% dù mettre au point pour les fameux transferts de Léontief avec Venance.},
comment aborder les
problèmes avec contrainte sur la commande etc. Une piste intéressante par son aspect
généraliste pourrait combiner l'adjoint et la programmation dynamique, cette dernière chose
que nous avions pratiquée avec Harold Vasselin pour optimiser la gestion de stocks
de froid et de chaud et ainsi climatiser à moindre frais un immeuble de bureaux au Québec
(cf sa thèse). Bref, y a de quoi alimenter quelques Zircons.

{\it Juillet 2004 : la version couplée à Minuit a été délivrée, cf exo ``pousse''
et le manuel du \Minik. Et en Juin 2005, on a aussi Lorenz 63, et les essais
de stepH avec NEDyM}. La conclusion de nombreux efforts de stepH avec Minuit n'est guère encouragente,
et il s'avère que le pilotage de Minuit lors de son balayage d'un espace paramétrique accidenté (le
modèle NEDyM comporte par construction de nombreuses poches de chaos) est difficile.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\Minik et LLSQ}\footnote{Linear Least SQuare fit, programme mettant en \oe uvre la méthode de
Householder de Paris Sud Informatique.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{llsq}
La difficulté pratique d'utiliser Minuit dans des cas complexes et l'intérêt
de la formule de Newton (cf annexe) nous amènent à proposer d'introduire dans \Minik le calcul des sensibilités
des variables à une liste de paramètres à déterminer. On a déjà montré comment il était possible
de calculer ces sensibilités par rapport à un paramètre choisi parmi les variables Fortran des modèles.
On propose ainsi une extension de ces calculs à la liste des {\bf free-parameter:---;}~, qui seraient inclus
dans la version 1.02.
\vspace{0.2cm}

%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{la Collab}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{{\bf retour au texte plus ancien}}Il existe déjà une certaine pratique dans la secte, on pourrait en profiter
pour faire un point sur les exos / difficultés / particularités rencontrées.

Sur ce qui a été fait dans ZOOM : avec un modèle explicite des sensibilités
(en avant, donc, \cad en doublant les équations du modèle), on gère le calcul d'un écart quadratique entre modèle / mesures
dans ZSTEER, et on utilise la méthode de Newton (via OLLSQ) pour déterminer la mise à jour
des paramètres (on relance donc ZOOM dans l'itération de l'assimilation
des données).
Olivier Fudym a ainsi fitté un modèle à deux paramètres par la méthode flash.
S'agissait d'un  modèle de séchage de grains (ou de bois ?). JLJ (et STB ?) on dû
faire un fit plus pédagogique avec une source sinusoïdale.
J'ignore ce qui a été fait avec Minuit.

C'est StepH qui veut joindre son modèle de climat-carbone avec celui de Vincent
Gitz qui travaille en «inter-temporel», \cad que le modèle détermine un coefficient
de propension des ménages à l'investissement au fil du temps censé maximiser
leur consommation intégrée sur un horizon fixé.
% Va faire ça dans \Minik en
% utilisant l'adjoint codé à Tarbes cet été avec stb d'après les notes de jyg.
Il s'agit donc non pas d'un fit paramétrique, mais de la recherche d'une
commande optimale, d'où ce texte. L'adjoint de \Minik a été validé par
comparaison avec les sensibilités automatiques à la perturbation d'une condition
initiale, et les deux calculs précédents permettent bien de retrouver la matrice
de propagation du système déterminée indépendamment (cf papier SLTC).
Ceci valide donc la partie matricielle,
les termes sources provenant de la commande et de l'intégrande de la fonction de coût
ont été validé avec la recherche d'optimum d'une loi de commande pour le cas
d'un exercice menant à une solution analytique (exo ``pousse'').
On obtient une précision de l'ordre du \% seulement, et ce problème numérique est
bien connu des fans du fit : d'où la version 1.01 de Mini$\_$ker en double précision, particulièrement
utile pour ce cas d'application. 
%\vspace{0.5cm}

\newpage
%\vspace{2cm}
\appendix{\bf ANNEXE : Le cas du FIT}
\label{annex}
%%------------
% 
\vspace{1cm}

Rappelons la méthode de Newton\footnote{que certains dénomment méthode de Gauss}
déjà largement utilisée avec ZOOM (cf Fudym, et JLJ et al.)
pour référence au problème de l'utilisation de l'adjoint dans le même contexte.
On dispose d'une série de mesures sur une séquence temporelle, rassemblées dans le
vecteur $\mu$. Par ailleurs, et aux mêmes instants pour simplifier, un modèle
fournit les mêmes observables $M$. Alors, au minimum de l'écart quadratique\footnote{on a
supposé encore pour simplifier que les mesures étaient indépendantes et de même incertitude.}
%
\begin{equation}
Q^2 = <(M-\mu)^2 >
\end{equation}
%
une variation $d\varpi$ des paramètres est liée à la variation du modèle par :
%
\begin{equation}
\ket{d\varpi}\, = -\left[\frac{\partial M}{\partial \varpi}^{\dagger}\frac{\partial M}{\partial \varpi}\right]^{-1}
       \frac{\partial M}{\partial \varpi}^{\dagger} \ket{M-\mu}
\label{newton}\end{equation}
%
On va établir cette formule car cela nous indiquera le lien avec l'adjoint.
%
\begin{equation}
Q^2(\varpi+d\varpi) = <(M-\mu)^2> + 2\scal{(M-\mu)}{\frac{\partial M}{\partial \varpi}d\varpi}
    + \scalm{d\varpi}{\frac{\partial M}{\partial \varpi}}{d\varpi}
\label{dpi}\end{equation}
%
par le fait que la matrice du terme du second degré en $d\varpi$ est symétrique définie
positive, le minimum est atteint pour $d\varpi$ donné par~\ref{newton}.

En effet, avec les termes du développement au second ordre de $Q^2$ :
%
\begin{equation}
Q^2(\varpi+d\varpi) = Q^2 + \scal{\nabla_{\varpi}Q^2}{d\varpi} + \frac{1}{2}\bra{d\varpi}\frac{\partial^2 Q^2}{\partial \varpi^2}\ket{d\varpi}
\label{d2pi}\end{equation}
%
on a explicitement :
%
\begin{eqnarray}
%\right\{
\nabla_{\varpi}Q^2 &=&2\bra{M-\mu}\frac{\partial M}{\partial \varpi}   \nonumber\\
\frac{\partial^2 Q^2}{\partial \varpi^2} &=& 2\frac{\partial M}{\partial \varpi}^{\dagger}\frac{\partial M}{\partial \varpi}
%\left.
\end{eqnarray}
%
d'où $d\varpi$ au minimum.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Utilisation de Minuit\footnote{\copyright CERN.}}
%-----------------------------------------------------------
Ce que fournit l'adjoint, c'est le premier de ces termes ($\nabla_{\varpi}Q^2$).
On constate ainsi dans ce cas
le travail qu'il reste à faire pour la recherche de l'extremum, avec pour ce cas
utilisation d'une méthode de gradients conjugués.

O.Fudym a pour une raison que j'ai oubliée,
fait appel à MINUIT\footnote{logiciel super-sophistiqué du CERN, méthode à
métrique variable, premières descentes par symplexe puis gradients conjugués.}.
Comme le modèle est appelé à partir de MINUIT, Olivier a du faire des acrobaties
avec des shells pour relancer ZOOM avec les paramètres récursivement modifiés.

La solution boite noire par utilisation de Minuit\footnote{option {\bf MinuiK} ou {\bf Monitor.}} transforme le programme {\bf Principal}
de \Minik en subroutine pilotée par {\bf DSDQ}, le programme gérant la recherche
d'optimum, avec une subroutine à l'interface {\bf FMINI}. Un manuel de Minuit
est disponible en ligne au CERN\footnote{voir aussi manuels \Minik sur la page web de {\bf ZOOM}}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Utilisation de la formule de Newton\footnote{version 1.02}}
%-----------------------------------------------------------
 Lors du traitement de ce problème de détermination paramétrique avec ZOOM, on a jusqu'à présent construit
explicitement le modèle de $\frac{\partial M}{\partial \varpi}$ dans chaque processeur de cellule et transfert,
stocké tous les éléments le long de la trajectoire  et
appliqué la formule de Newton~\ref{newton} par
appel à la routine OLLSQ\footnote{qui fournit également la matrice de covariance des
paramètres. Cette routine est issue d'un agorithme de Golub et programmée par l'équipe du PSI (Paris Sud Informatique
de l'université d'Orsay) Jean Rauel et Germain Bonne pour l'UNIVAC -- LLSQ et MLLSQ. Le nom de OLLSQ correspond
à la version jyg, avec la convention ``raw-wise'' adoptée pour cause de compatibilité avec la CERLIB et ZOOM.}.
Avec un appel à LLSQ après ZSTEER en fin de trajectoire, on détermine
les nouveaux paramètres et on boucle au début de ZINIT.

La version 1.02 de \Minik va permettre d'utiliser la formule de Newton beaucoup plus directement. On prévoit
en effet que la {\bf free-parameters'} liste provoque le calcul de la matrice $\frac{\partial M}{\partial \varpi}$
à chaque pas de temps du calcul.

Lors d'une rencontre avec un instant ($i$) de mesure fournissant $ket{\mu_i}$,
signalé par {\bf ZObs True}, on incrémente le calcul de $d\varpi$ (variable {\bf ddpi}, fichier de sortie
{\bf ddpi.data}) par la formule~(\ref{newton}) que l'on décline sur chaque instant de mesure ($i$) avec
l'hypothèse que la matrice d'erreur de mesure $\sigma_i$ est indépendante d'une prise de mesure
à l'autre :
%
\begin{equation}
\boxed{
\ket{d\varpi}\,=\sum_i dd\varpi_i\, = -\sum_i \left[\frac{\partial M}{\partial \varpi}^{\dagger}\sigma^{-1}_i\frac{\partial M}{\partial \varpi}\right]^{-1}_i
       \frac{\partial M}{\partial \varpi}^{\dagger} \sigma^{-1}_i\ket{M_i-\mu_i}
}\label{newtonbis}\end{equation}
%
On peut éventuellement prévoir de ne commencer ces calculs que lorsque  la dimension du vecteur des mesures dépasse
suffisammment le nombre de paramètres. Dans la méthode de Householder, on minimise en pratique la norme
euclidienne de $\displaystyle \frac{\partial M}{\partial \varpi} \sigma^{-\frac{1}{2}}_i + \ket{M_i-\mu_i}$
(cf~\ref{dpi}).
Les cas les plus habituellement traités possèdent une simple erreur de mesure $\sigma^{\frac{1}{2}}$ diagonale.
\vspace{1cm}

\hrule
\vspace{1cm}
\end{document}
\end

