IRIA
laboria

laboratoire de recherche
en informatique
et automatique

# Contrôle sous-optimal d'un réacteur chimique par une méthode de perturbation singulière

Maurice Robin

<figure><img src="image_1.png" /></figure>

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

Rapport de Recherche n°2

Janvier 1973
CONTROLE SOUS-OPTIMAL D'UN RÉACTEUR CHIMIQUE
PAR UNE MÉTHODE DE PERTURBATION SINGULIERE

Maurice Robin

RÉSUMÉ

On utilise une méthode de perturbations singulières pour contrôler
un système d'équations différentielles constituant un modèle simple de
réacteur chimique. On compare les résultats avec ceux obtenus par un
algorithme de gradient classique.
P L A N

INTRODUCTION. 1

I - PROBLEME DE CONTROLE D'UN REACTEUR. 3

1.1. Description du modèle. 3

1.2. Introduction d'un petit paramètre. 6

1.3. Critère d'optimisation. 7

1.4. Essais numériques. 7

1.5. Résultats. 8

II - EXPOSE FORMEL DE LA METHODE. 24

2.1. Formulation du problème. 24

2.2. Développement. 26

2.3. Un cas simple. 29

2.4. Méthode numérique. 30

CONCLUSIONS. 32

BIBLIOGRAPHIE. 33
## INTRODUCTION

Les modèles de systèmes de production contiennent souvent des variables dont l'évolution est très rapide par rapport à celle d'autres variables du système. Si le modèle est linéaire ceci se traduit par des valeurs propres d'ordres de grandeur très différents.

Les méthodes explicites d'intégration ne sont alors stables que pour des pas de discrétisation très faibles ce qui conduit très rapidement à des temps de calcul prohibitifs pour contrôler le système.

Ce type de difficultés nous a conduit à essayer (sur un exemple très simple mais assez caractéristique) une méthode issue de la théorie des perturbations sin-gulières (5), (10), (11), (4) et déjà utilisée par O'MALLAY (1) pour le problème de contrôle linéaire-quadratique. SANNUTI et KOKOTOVIC (6) ont envi-sagé un cas légèrement différent. (cf § 3).

Les problèmes analogues pour les équations aux dérivées partielles sont traités dans LIONS (12). La méthode formelle sera exposée au §2 en supposant que le problème a été ramené à la forme :

$$ (I) \quad \left\{ \begin{array}{l} \displaystyle \frac{dx}{dt} = f(x,y,u,t,\epsilon) & x(0) = \xi \\[1em] \displaystyle \epsilon \frac{dy}{dt} = g(x,y,u,t,\epsilon) & y(0) = \eta \end{array} \right. $$

où est un "petit paramètre" scalaire (on peut étendre le problème à plusieurs petits paramètres cf (14)) où l'on doit minimiser par rapport à u une fonctionnelle $J(u, \epsilon)$
- La définition du petit paramètre sur un problème particulier peut poser quelques problèmes qui seront évoqués au § 1.

- On a comparé cette méthode avec un algorithme de gradient où la trajectoire et le gradient sont calculés par intégration des équations directes et adjointes.
# I - PROBLEME DE CONTROLE D'UN REACTEUR CHIMIQUE.

## 1.1. Description du modèle - hypothèses physiques.

- On considère un réacteur chimique qui transforme une matière première de concentration donnée en constituant A en un produit B qui est stocké.

- On fait les hypothèses suivantes :

- On suppose que la température, les concentrations, et l'activité du catalyseur sont les mêmes dans tout le réacteur (ceci permet d'obtenir un modèle en équations différentielles).

- La réaction est exothermique, du premier ordre, et accélérée par un catalyseur caractérisé uniquement par un indice d'activité ($z_1$) ; la vitesse de réaction lui est proportionnelle.

- Le réacteur est refroidi par la circulation d'un fluide de température fixée ($a_{33}$), et on contrôle le coefficient d'échange de chaleur $u_2$ entre le refroidisseur et le réacteur (qui est une fonction non explicitée ici du débit de ce refroidisseur cf (13)).

- On contrôle également le débit d'entrée ($u_1$) dans le réacteur et dont la concentration $C_0$ est prise comme unité.

La production de B est donc $u_1(t) \cdot (1-z_2(t))$ si $z_2(t)$ désigne la concentration de A dans le catalyseur.

Les équations résultent alors de l'écriture des bilans matières et énergie pour le réacteur et le stock ; $z_i$ ($i=1,4$ désignant respectivement l'activité du catalyseur, la concentration de A, la température du réacteur et le niveau de stock de produit B :
$$
\begin{equation} \tag{1.1}
\left\{
\begin{aligned}
\frac{dz_1}{dt} &= - a_{11} (1 + a_{12} z_2) z_1 \exp \left(-\frac{a_{13}}{z_2}\right) \\
\frac{dz_2}{dt} &= a_{21} u_1(t) (1-z_2) - a_{22} z_1 z_2 \exp \left(-\frac{a_{23}}{z_2}\right) \\
\frac{dz_3}{dt} &= a_{31} u_1(t) (1-z_3) + a_{32} u_2(t) (a_{33} - z_3) + \\
&\qquad + a_{34} (1+a_{35} z_3) z_1 z_2 \exp \left(-\frac{a_{36}}{z_3}\right) \\
\frac{dz_4}{dt} &= u_1(t) (1-z_2) - d(t)
\end{aligned}
\right.
\end{equation}
$$

$z_i(0)$ donné

$t \in (0,T)$

Notations de l'exemple et données numériques utilisées.

$z_1(x_1)$ activité du catalyseur (sans dimension),

$z_2(y_1)$ concentration du produit A (normalisée-sans dimension),

$z_3(y_2)$ température dans le réacteur (normalisée),

$z_4(x_2)$ niveau du stock (masse),

$u_1(t)$ débit d'entrée (masse par unité de temps),

$u_2(t)$ coefficient de transfert de chaleur entre le réacteur et le refroidisseur,

$d(t)$ fonction donnée, représente le taux de "vente" de produit B.
Valeurs numériques :

$$ 
\begin{array}{l}
a_{13} = a_{23} = a_{37} = 10. \\
a_{11} = 10^4 \\
a_{12} = a_{21} = 0.1 \\
a_{31} = 0.6 \\
a_{32} = 0.25 \\
a_{33} = 0.8 \\
a_{35} = 0.1 \\
d(t) = 20 \quad \forall t
\end{array} 
\quad 10^5 \le a_{22} = a_{34} \le 10^7 \\
\text{(suivant les exemples)}
$$

valeurs initiales :

$$ x_1(0) = 1. \qquad y_1(0) = 1. $$

$$ x_2(0) = 50. \qquad y_2(0) = 1. $$

Dans tous les exemples on calcule la solution sur (0,1) avec le pas h=0.01.
## 1.2. Introduction d'un petit paramètre.

- Avec les données numériques citées plus loin, il faut un pas de $10^{-4}$ pour qu'une méthode de RUNGE-KUTTA (explicite) soit stable.

Dans ce cas présent, la "raideur" du système vient de ce que la vitesse de réaction (essentiellement $a_{22} \exp(-\frac{a_{23}}{z_3})$) est toujours grande devant la vitesse de "désactivation" du catalyseur ($a_{11} \exp(-\frac{a_{13}}{z_3})$). Ceci amène à essayer pour $\epsilon$ le rapport des coefficients $a_{11}$ et $a_{12}$ dans le cas où $a_{23}$ et $a_{13}$ sont du même ordre de grandeur.

Compte tenu des ordres de grandeurs des autres coefficients (cf valeurs numériques) on a adopté le modèle suivant :

$$ (1.2) \quad \left\{ \begin{aligned} \frac{dx_1}{dt} &= -a_{11} (1+a_{12}y_1) x_1 \exp(-\frac{a_{13}}{y_2}) \\ \frac{dx_2}{dt} &= u_1(t)(1-y_1) - d(t) \\ \epsilon \frac{dy_1}{dt} &= \epsilon a_{21} u_1(1-y_1) - a'_{22} y_1 x_1 \exp(-\frac{a_{23}}{y_2}) \\ \epsilon \frac{dy_2}{dt} &= \epsilon a_{31} u_1(1-y_2) + a'_{32} u_2 (a_{33} - y_2) + a'_{34} (1+a_{35}y_2) y_1 x_1 \exp(-\frac{a_{36}}{y_2}) \end{aligned} \right. $$ où en règle générale $a' = \epsilon a$

## Extensions possibles du modèle.

L'exemple ci-dessus est facilement généralisable au cas de plusieurs réactions simultanées faisant intervenir plusieurs produits, éventuellement plusieurs réacteurs. Ceci n'introduit que des complications techniques et la même méthode reste applicable.
### 1.3. **Problème de contrôle.**

On a construit un problème test en se donnant un contrôle de référence $u_d(t)$ et en intégrant les équations (1.4) pour obtenir $z_d(t)$. On prend alors un critère d'optimisation du type minimisation de l'écart par rapport à l'objectif $(u_d, z_d)$, soit

$$ (1.3) \qquad J = \int_0^T \left\{ \sum_i \alpha_i |x_i(t) - x_{ji}(t)|^2 + \sum_j \beta_j |y_j(t) - y_{dj}(t)|^2 + \sum_\ell \gamma_\ell |u_\ell(t) - u_{d\ell}(t)|^2 \right\} dt $$

le coût optimal est nul bien entendu.

### 1.4. **Essais numériques.**

On a traité numériquement deux problèmes différents :

#### Problème 1

Les coefficients du modèle initial (1.1) étant fixés, comment varie l'appro- ximation quand on adopte des définitions différentes pour le petit paramètre ?

Ici la trajectoire optimale est la même quel que soit $\varepsilon$. On compare égale- ment le résultat à celui obtenu par l'algorithme de gradient.

#### Problème 2

On envisage des systèmes plus ou moins "raides" en faisant varier les coef- ficients du modèle. On prend $\varepsilon = a_{11}/a_{22}$ (la trajectoire optimale est alors différente dans chaque cas).

Dans les deux problèmes on arrête le développement au terme en $\varepsilon$ (un développement d'ordre supérieur semble quasiment inutilisable cf § 2.4).
Exemple, supposons que l'on veuille ramener le système

$$ \frac{dx}{dt} = x - y \qquad \frac{dy}{dt} = -x - 100y $$

à la forme (I)

On peut choisir $\epsilon = \frac{1}{100}$ et prendre

$$ \epsilon \frac{dy}{dt} = -\epsilon x - y $$

Mais on peut tout aussi bien prendre $\epsilon = \frac{1}{50}$

$$ \text{d'où } \epsilon \frac{dy}{dt} = -\epsilon x - 2y $$

Dans le cas non linéaire, l'approximation d'ordre 2 en $\epsilon$ sera différente suivant les définitions du petit paramètre (ce n'est pas vrai dans le cas linéaire).

Le problème 2 serait ici d'envisager

$$ \frac{dx}{dt} = x - y $$

$$ \frac{dy}{dt} = -x - ay \quad \text{avec } a=10, a=100, a=1000 $$

et toujours $\epsilon = \frac{1}{a}$

## 1.5. Résultats

Les tableaux suivants donnent les temps de calculs des deux méthodes et la valeur de la fonction économique pour la méthode sous optimale.

On donne également les courbes de la trajectoire et du contrôle pour chaque exemple.

### - Remarques sur ces résultats

1) - Il est difficile de comparer la méthode sous-optimale à la méthode du gradient :
Il faut tenir compte de la précision obtenue sur le coût, et du temps de calcul. Or le temps de calcul d'une méthode de gradient dépend beaucoup de la solution initiale, de la procédure d'ajustement du ρ etc.

On ne peut pas savoir ici quel temps de calcul de l'algorithme de gradient est nécessaire pour atteindre une précision égale à celle de la méthode sous optimale.

2) - Les domaines d'application des deux méthodes sont différents comme le montre le cas $ε = 10^{-3}$ du problème 2.

Le temps de calcul d'une itération rend la méthode de gradient très coûteuse alors que la précision de la méthode sous optimale est améliorée considérablement.

### Note

On a utilisé un algorithme de gradient à pas divisé avec valeurs pour le test d'arrêt (variation relative de la fonction économique)

test 1. variation relative < $10^{-2}$

test 2 variation relative < $10^{-4}$

Il n'y a que les résultats du test 1 qui apparaissent sur les courbes.
$$
\underline{\text{Problème 1}} \qquad (a_{22} = a_{24} = 10^6)
$$

Méthode approchée

<table><thead><tr><th>ε</th><th>temps de calcul (mn)</th><th>valeur de la fonctionnelle</th></tr></thead><tbody><tr><td>10<sup>-4</sup></td><td>24</td><td>15.7</td></tr><tr><td>10<sup>-3</sup></td><td>9</td><td>13.6</td></tr><tr><td>5 10<sup>-3</sup></td><td>7</td><td>14.6</td></tr><tr><td>10<sup>-2</sup></td><td>4.5</td><td>15.9</td></tr><tr><td>5 10<sup>-2</sup></td><td>2.5</td><td>20.8</td></tr><tr><td>10<sup>-1</sup></td><td>1.5</td><td>36.3</td></tr></tbody></table>

Méthode de gradient

<table><thead><tr><th>test</th><th>temps 1 itération mn</th><th>valeur. coût</th><th>temps total mn</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1.4</td><td>10</td></tr><tr><td>2</td><td>1</td><td>0.5</td><td>21</td></tr></tbody></table>

Note : les temps sont donnés pour CII 10 070
## - PROBLEME 2 -

$$ ( \varepsilon = a_{11}/a_{22} ) $$

Méthode approchée

Méthode de gradient

<table><thead><tr><th>ε</th><th>temps calcul mn</th><th>valeur du coût</th><th>temps l'itération mn</th><th>pas d'intégration</th><th>temps total mn</th><th>valeur du coût</th></tr></thead><tbody><tr><td>10<sup>-3</sup></td><td>11</td><td>1.6</td><td>5</td><td>10<sup>-5</sup></td><td></td><td></td></tr><tr><td>10<sup>-2</sup></td><td>4.5</td><td>15.9</td><td>1.</td><td>10<sup>-4</sup></td><td>21</td><td>0.5</td></tr><tr><td>10<sup>-1</sup></td><td>2.</td><td>64.</td><td>0.2</td><td>10<sup>-3</sup></td><td>8</td><td>0.05</td></tr></tbody></table>
<figure><img src="image_2.png" /></figure>

FIGURE 1

CAS 1 CONTROLE UI

ε = 10⁻¹

ε = 10⁻² o 5 10⁻²

ε = 10⁻³ +

Trait plein : Contrôle optimal

-x Méthode de gradient TEST 1
<figure><img src="image_3.png" /></figure>

FIGURE 2 -
<figure><img src="image_4.png" /><figcaption>FIGURE 3</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 4</figcaption></figure>

CAS 1

ETAT X2 (Stock)

- + ε = 10⁻³

- - ε = 10⁻²

- 0 ε = 5 · 10⁻²

-- ε = 10⁻¹

x gradient TEST 1
<figure><img src="image_6.png" /></figure>

FIGURE 5
<figure><img src="image_7.png" /><figcaption>- FIGURE 6 -</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 7</figcaption></figure>
- FIGURE 8 -

<figure><img src="image_9.png" /></figure>

CAS 2
ETAT XI (Catalyseur)

$$ \epsilon = 10^{-3} $$
<figure><img src="image_10.png" /><figcaption>FIGURE 9</figcaption></figure>
<figure><img src="image_11.png" /></figure>

FIGURE 10

CAS 2
ETAT Y2
ε = 10⁻³
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>
## II - EXPOSE FORMEL DE LA METHODE.

### 2.1. Formulation du problème.

- On supposera que le problème a été ramené à la forme suivante

$$ (2.1) \qquad \frac{dx(t)}{dt} = f(x(t),y(t),u(t),t,\varepsilon) \quad x(0) = \xi $$

$$ (2.2) \qquad \varepsilon \frac{dy(t)}{dt} = g(x(t),y(t),u(t),t,\varepsilon) \quad y(0) = \eta $$

$$ \text{où } x(t) \in \mathbb{R}^n, y(t) \in \mathbb{R}^m, u(t) \in \mathbb{R}^p, \quad t \in (0,T) $$

$\varepsilon$ est un "petit" paramètre scalaire.

On peut minimiser

$$ (2.3) \qquad J_{\varepsilon} = \int_{0}^{T} L(x,y,u,t,\varepsilon)dt $$

sous les contraintes (1.1)-(1.2).

On notera $P_\varepsilon$ le problème précédent et $P_0$ le problème

Minimiser

$$ (2.4) \qquad \left\{ \begin{array}{l} \displaystyle \int_{0}^{T} L(x,y,u,t,0)dt \\[1em] \frac{dx}{dt} = f(x,y,u,t,0) \quad x(0) = \xi \\[1em] g(x,y,u,t,0) = 0 \end{array} \right. $$

On fera alors les hypothèses suivantes :

H-1 $P_0$ et $P_\varepsilon$ ont une solution unique

H-2 Si $u_0$ est la solution de $P_0$, $x_0$, $y_0$ la trajectoire correspondante, on suppose que $f,g,L$ sont continûment différentiables $R+1$ fois en $(x,y,u,t,\varepsilon)$ au voisinage de $(x_0(t),y_0(t),t,0)$ $\forall t \in (0,T)$
H-3 $\frac{\partial g}{\partial y}(x_o(t), y_o(t), t, 0)$ est inversible $\forall t \in (0, T)$

H-4 $\frac{\partial g}{\partial y}(x_o(t), y_o(t), t, 0)$ a des valeurs propres à parties réelles strictement négatives $(R_l(\lambda_i(t)) \le -\mu < 0) \quad \forall t \in (0, T)$

Ecrivons le système adjoint correspondant à $P_\epsilon$

$$ (2.5) \qquad - \frac{dy}{dt} = f_x^* \varphi + g_x^* \psi - \frac{\partial L}{\partial x} \qquad \varphi(\tau) = 0 $$

$$ (2.6) \qquad - \epsilon \frac{d\psi}{dt} = f_y^* \varphi + g_y^* \psi - \frac{\partial L}{\partial y} \qquad \psi(\tau) = 0 $$

Les conditions nécessaires d'optimalités donnent

$$ (2.7) \qquad f_u^* \varphi + g_u^* \psi - \frac{\partial L}{\partial u} = 0 $$

(* signifie : transposée)

Le problème que l'on veut résoudre sera alors :

Trouver une solution de (2.1), (2.2), (2.5), (2.6), (2.7) soit $x(t, \epsilon)$, $y(t, \epsilon)$, $\varphi(t, \epsilon)$, $\psi(t, \epsilon)$, $u(t, \epsilon)$ sous la forme (cf. (4) par exemple)

$$ (2.8) \qquad \left\{ \begin{aligned} x(t, \epsilon) &= \hat{x}(t, \epsilon) + X(\tau, \epsilon) + \chi(\sigma, \epsilon) \\ y(t, \epsilon) &= \hat{y}(t, \epsilon) + Y(\tau, \epsilon) + \tilde{\chi}(\sigma, \epsilon) \\ \varphi(t, \epsilon) &= \hat{\varphi}(t, \epsilon) + \Phi(\tau, \epsilon) + \Theta(\sigma, \epsilon) \\ \psi(t, \epsilon) &= \hat{\psi}(t, \epsilon) + \Psi(\tau, \epsilon) + \Lambda(\sigma, \epsilon) \end{aligned} \right. $$

Avec $\tau = t/\epsilon$ $\sigma = (T-t)/\epsilon$

où $X, Y, \Phi, \Psi$ (resp. $\tilde{x}$ ...) tendent vers 0 quand $\tau$ (resp. $\sigma$) tend vers l'infini.

Plus précisément on cherche les termes des développements en $\epsilon$ de $\hat{x}$, $X$, $\tilde{x}$...etc
## 2.2. Développements formels.

Il est immédiat que le premier terme des développements de $\hat{x}, \hat{y}, \hat{\phi}, \hat{\psi}, \hat{u}$ est obtenu en résolvant $P_o$.

à l'ordre 1 en $\epsilon$, on obtient (pour $\hat{x}...etc$)

$$ (2.9) \quad \left\{ \begin{aligned} \frac{dx_1}{dt} &= f_x x_1 + f_y y_1 + f_u u_1 + f_\epsilon \\ \frac{dy_o}{dt} &= g_x x_1 + g_y y_1 + g_u u_1 + g_\epsilon \end{aligned} \right. $$

et en posant $H = f^*\varphi + g^*\psi - L$

$$ (2.10) \quad \left\{ \begin{aligned} -\frac{d\varphi_1}{dt} &= H_{xx} x_1 + H_{xy} y_1 + H_{x\psi} \psi_1 + H_{xu} u_1 + H_{x\epsilon} \\ -\frac{d\psi_0}{dt} &= H_{yx} x_1 + \dots \\ H_{ux} x_1 + H_{uy} y_1 + H_{u\psi} \psi_1 + \dots &= + H_{y\epsilon} \\ H_{u\epsilon} &= 0 \end{aligned} \right. $$

où toutes les dérivées sont prises en $x_o(t), y_o(t), \psi_o(t), \psi_o(t), u_o(t)$ déterminés par $P_o$.

On remarque que $P_o$ peut être résolu seul, mais que ce n'est plus le cas de $P_1$ car on ne connaît pas les conditions initiales sur $x_1$ et finales sur $\varphi_1$.

On va voir que ces conditions nécessitent la détermination préalable des termes correcteurs $y_o, \psi_o, y_o, \Lambda_o, x_1, \Phi_1$.

L'identification des développements en $t=0$ et $t=T$ donne en effet (on a déjà pris $x_o(0) = \xi$)
$$
(2.11) \qquad X_o(0) = 0
$$

$$
(2.12) \qquad x_t(0) + x_r(0) = 0 \quad r \ge 1
$$

$$
(2.13) \quad y_r(0) + Y_r(0) = \eta \delta_{ro} \quad r \ge 0
$$

$$
(2.14) \qquad \varphi_o(T) = 0
$$

$$
(2.15) \qquad \varphi_r(T) + \theta_r(0) = 0 \quad r \ge 1
$$

$$
(2.16) \qquad \Theta_o(0) = 0
$$

$$
(2.17) \qquad \psi_r(T) + \Lambda_r(0) = 0 \quad r \ge 0
$$

$X, Y, \Theta, \Lambda$ (resp. $X, Y, \Phi, \Psi$) sont négligés en $\frac{T}{\varepsilon}$

$P_o$ déterminera $y_o(0)$ et $\psi_o(T)$ d'où par (2.13) et (2.17)

$$
(2.18) \quad \left\{
\begin{aligned}
Y_o(0) &= \eta - y_o(0) \\
\Lambda_o(0) &= - \psi_o(T)
\end{aligned}
\right.
$$

Il faut maintenant pouvoir calculer $y_o(\tau)$, $\Lambda_o(\sigma)$ puis $X_l(\tau)$ et $\Theta_l(\sigma)$ de
façon à obtenir $x_l(0)$ et $\varphi_l(T)$ par (2.12) et (2.15).

Le premier terme correctif en $\tau$ est obtenu en développant en $\epsilon$ les équations
correspondant à $x - \hat{x}$, $y - \hat{y}$ etc... après avoir fait le changement de va-
riable $\tau = t/\epsilon$.

On obtient

$$
\frac{dX_o(\tau)}{d\tau} = 0
$$

$$
\frac{d \Phi_o(\tau)}{d \tau} = 0
$$
$$
(2.19) \quad - \frac{dY_o}{d\tau} = g(x_o(0) + X_o(\tau), y_o(0) + Y_o(\tau), u_o(0) + U_o(\tau), 0, 0)
$$

$$
(2.20) \quad - \frac{d\Psi_o}{d\tau} = f_y^* \cdot (\varphi_o(0) + \Phi_o(\tau)) + \tilde{g}_y^* (\psi_o(0) + \psi_o(\tau)) - \frac{\partial L}{\partial y}
$$

$$
(2.21) \quad f_u^* (\varphi_o(0) + \Phi_o(\tau)) + g_u^* (\psi_o(0) + \psi_o(\tau)) - \frac{\partial L}{\partial u} = 0
$$

où toutes les dérivées sont prises en x₀(0) + X₀(τ), y₀(0) + Y₀(τ)...

$$
\text{On a bien entendu } X_o(\tau) \equiv 0 \quad \Phi_o(\tau) \equiv 0 \quad \left( \lim_{T \to \infty} X_o(\tau) = 0 \right)
$$

On a $Y_o(0) - (2.10)$ et en pratique on prend $\psi_o(\frac{T}{\varepsilon}) = 0$. On peut alors
intégrer

$$
(2.22) \quad \frac{dX_1}{d\tau} = f(x_o(0), y_o(0) + Y_o(\tau), u_o(0) + U_o(\tau), 0, 0) - f(x_o(0), y_o(0), u_o(0), 0, 0)
$$

$$
(2.23) \text{d'où } x_1(0) = -x_1(0) = \int_0^\infty \left(\frac{dX_1}{d\tau}\right) d\tau
$$

Le calcul est identique pour le premier terme correctif en σ.

La détermination de $\mathcal{Y}_o$, $\Lambda_o$ permettant d'intégrer

$$
\begin{equation} \tag{2.24}
\left\{
\begin{aligned}
\frac{d}{d\sigma} \left( \frac{1}{\sigma} \right) &= g_x(x_o(T), y_o(T) + \mathcal{Y}_o(\sigma), u_o(T) + \mathcal{U}_o(\sigma), T, 0) \cdot (\psi_o(T) + \Lambda_o(\sigma)) - \\
&\quad - \frac{\partial L}{\partial x}(x_o(T), y_o(T) + \mathcal{Y}_o(\sigma), \dots) - g_x(x_o(T), y_o(T), u_o(T), T, 0) \psi_o(T) + \\
&\quad + \frac{\partial L}{\partial y}(x_o(T), y_o(T), u_o(T), T, 0) \\
\text{et } \psi_1(T) &= - \frac{1}{\sigma}(0) = \int_0^\infty \left( \frac{d}{d\sigma} \left( \frac{1}{\sigma} \right) \right) d\sigma
\end{aligned}
\right.
\end{equation}
$$
2.3. **Un cas simple** (cf (6))

Supposons f,y linéaires en y, et en u, L indépendant de y et quadratique en u.

i.e.

$$
\frac{dx}{dt} = A(x,t,\epsilon) y + B(x,t,\epsilon)u + C(x,t,\epsilon)
$$

$$
(2.25) \quad \epsilon \frac{dy}{dt} = D(x,t,\epsilon) y + E(x,t,\epsilon)u + F(x,t,\epsilon)
$$

$$
J = \int_{0}^{T} \left[ L(x,t,\epsilon) + \frac{1}{2} u^{T} R(x,t,\epsilon) u \right] dt
$$

(Toujours avec les hypothèses H1 à H4 qui ici suffisent)

On note

$$
A_o = A(x_o(0), 0, 0)
$$

$$
D_o = D(x_o(0), 0, 0)
$$

$$
D_T = D(x_o(T), T, 0) \quad \text{etc...}
$$

On obtient :

$$
(2.26) \quad \frac{dy_o}{d\tau} (\tau) = D_o(y_o(0) + Y_o(\tau)) + E_o(u_o(0) + U_o(\tau)) + F_o
$$

$$
(2.27) \quad - \frac{d\psi^o}{d\tau}(\tau) = D_o^* \psi_o
$$

$$
(2.28) \quad U_o(\tau) = R_o^{-1} E_o \psi_o
$$

$$
(2.29) \text{d'autre part} \quad \psi_o(T) = D_o^{*-1} A_o^* \varphi_o(T) = 0
$$

$$
(2.30) \text{d'où } \Lambda_o(0) = 0 \quad \text{(par (2.17))}
$$
$$
(2.31) \qquad \frac{d\mathcal{H}(\sigma)}{d\sigma} = - D_T \mathcal{Y}_o(\sigma) - E_T \mathcal{U}_o(\sigma)
$$

$$
(2.32) \qquad \frac{d\Lambda_o}{d\sigma}(\sigma) = D_T^* \Lambda_o(\sigma)
$$

$$
(2.33) \qquad \mathcal{U}_o(\sigma) = R_T^{-1} E_T^* \Lambda_o(\sigma)
$$

$$
(3.8) \text{ et } (3.6) \text{ donnent } \Lambda_o(\sigma) \equiv 0 \text{ donc } \mathcal{U}_o(\sigma) \equiv 0
$$

H4 et (3.7) impliquent que $\mathcal{Y}_o(\sigma) \equiv 0$ est la seule solution qui vérifie

$$
\lim_{\sigma \to \infty} \mathcal{Y}_o(\sigma) = 0
$$

Comme on a également

$$
(2.34) \qquad \frac{d}{d\sigma} \mathcal{Y}_1(\sigma) = (D'_{xT} \mathcal{Y}_o(\sigma) + E'_{xT} \mathcal{U}_o(\sigma)) \quad \psi_o(T)
$$

On a $@_1(\sigma) = 0 \quad \forall \sigma$ donc $\varphi_1(T) = 0$

2.4. **Méthode numérique.**

On s'arrêtera au premier ordre en ε, ce qui s'est montré justifié dans
les applications numériques. ( par contre,on ne peut en général se contenter
du premier terme (solution de P₀) qui n'a souvent aucune signification con-
crète).

La procédure qui ressort du § 2 se résume ainsi.

1) Résoudre le problème P<sub>o</sub>

C'est un problème de contrôle de dimension réduite par rapport à P<sub>ε</sub> et
qui a l'avantage de ne porter que sur les équations les moins "raides".

2) Résoudre les problèmes (à deux points) en $Y_o, \Psi_o$ et $\mathcal{Y}_o, \Lambda_o$

3) Intégrer les équations en $X_l$ et $\mathcal{Y}_l$ d'où les conditions $x_l(0), \varphi_l(T)$.

4) Résoudre le problème à deux points en $x_l, \varphi_l$.
5) On a alors

$$
\left\{
\begin{aligned}
x(t) &= x_o(t) + \epsilon x_1(t) \\
y(t) &= y_o(t) + Y_o(\tau) + \mathcal{Y}_o(\sigma) + \epsilon y_1(t) \\
u(t) &= u_o(t) + U_o(\tau) + \mathcal{U}_o(\sigma) + \epsilon u_1(t)
\end{aligned}
\right.
$$

On voit facilement les limitations de l'application de cette méthode : notamment l'étape 2 qui consiste à résoudre deux problèmes à deux points sur 0,T/ε
Il est évident que ceci ne sera praticable que quand le problème permet des simplifications ou des approximations supplémentaires.

Considérations pratiques

Pour l'exemple du § 1 :

1 - La résolution de $g(x_o, y_o, u^o, t, 0) = 0$ est i immédiate ($y_1^o = 0, y_1^o = a_{33}$)
Les équations en $x_o, y_o$ sont intégrées par la méthode de la tangente améliorée avec un pas $h=0.01$ sur $(0,1)$

2 - Le calcul de $y_o, \psi_o$ est effectué par itération : le membre d'itérations est en général très faible car $y_o$ et $\psi_o$ tendent très rapidement vers 0 (rappelons qu'on n'utilise que les valeurs $y_o(\frac{nh}{\epsilon})$.

La détermination de $x_l(0), \psi_l(T)$ (qui comprend le calcul précédent) est l'étape la plus coûteuse en temps de calcul car il nécessite l'intégration de systèmes différentiels sur $(0, \frac{T}{\epsilon})$ avec un pas de l'ordre de $\frac{h}{\epsilon}$.

3 - Le problème en $x_l, \psi_l$ est linéaire (car le critère est quadratique en $u$) et se ramène donc à l'intégration de systèmes différentiels (méthode de la tangente améliorée, pas $h$).

4 - Dans l'algorithme de gradient, les équations sont intégrées par une méthode de Runge-Kutta explicite d'ordre 4.
- **CONCLUSIONS** -

Il ressort des résultats obtenus que la méthode présentée ne donne une bonne approximation de la solution que si le système est assez "raide" et ceci avec un temps de calcul raisonnable. Il ne convient sans doute d'utiliser cette méthode que quand le temps d'intégration du système par une méthode classique laisse prévoir que l'algorithme de gradient sera très coûteux.

Il faut signaler également deux difficultés de mise en oeuvre de cette méthode.

- Il n'ya pas de procédure systématique pour déterminer le petit paramètre, sauf peut-être une étude numérique préalable comme celle envisagée dans MIRANKER (3).

- Les calculs préalables à la programmation deviennent rapidement important quand la dimension du système augmente.
- BIBLIOGRAPHIE -

(1) O'MALLEY R.E. SIAM J. Control Vol. 10 n°3 Août 1972.

(2) O'MALLEY R.E. SIAM J. Applied Math. Vol. 17 n°3 Mai 1969.

(3) MIRANKER Conférence I.R.I.A. 1972

(4) HOPPENSTEADT F. Communication in pure and Appl. Math XXIV 1971. p807-840.

(5) WASOV Asymptotic Expansion for ordinary differential equations Wiley.

(6) SANNUTJ et KOKOTOVIC Joint Automatic control conference 1969.

(7) FLATTO-LEVINSON J. rational Mech. Anal. 4 (1955) p.943-950

(8) LEVIN-LEVINSON J. Rational Mech. Anal. 3 (1954) p.247-270

(9) HOPPENSTEADT F. J. of Diff. Equations Vol 5 1969 p 106

(10) VASIL'EVA A.B. Russian Math Surveys 1963-18-3 p13

(11) VASIL'EVA A.B. U.R.S.S. Computational Math and Math. physics. 1963 3-4 p823-863.

(12) LIONS J.L. Livre à paraître

(13) GOULD Chemical process control - Addison-Wesley. 1969.

(14) HOPPENSTEADT F. J. of Differential Equations Vol. 5 1969. pp 106-116.

(15) O'MALLEY R.E. Aut. Rat. Mech. and Annal Vol 40. 1971. n°3 p 209.