Dw-02953 RR221

IRIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

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

CONTRÔLE OPTIMAL
DE SYSTÈMES GOUVERNÉS
PAR DES PROBLÈMES
AUX VALEURS PROPRES

Fulbert MIGNOT
Christian SAGUEZ
Jean-Paul VAN DE WIELE

Rapport de Recherche N° 221

Mars 1977

42p
CONTRÔLE OPTIMAL DE SYSTEMES GOUVERNÉS
PAR DES PROBLEMES AUX VALEURS PROPRES

F. Mignot*, C. Saguez**, J.P. Van de Wiele***

Résumé :

Nous étudions dans cet article le contrôle optimal de systèmes gouvernés par des problèmes aux valeurs propres. Nous montrons l’existence d’un contrôle optimal et nous obtenons des conditions d’optimalité du premier ordre. Nous étudions aussi deux exemples numériques dont nous présentons la solution complète. Nous utilisons la méthode des éléments finis pour l’approximation et une méthode de gradient pour la minimisation. Enfin, nous exposons dans le détail les algorithmes numériques.

Abstract :

In this paper, we study the optimal control of systems governed by an eigenvalue problem. We prove the existence of an optimal control and we obtain first order optimality conditions. We present also two numerical examples and we give the complete solution of both problems, involving approximation by the finite element method and minimization by a gradient method. Finally numerical algorithms are detailed.

<footnote>* U.E.R. Mathématiques, Lille I.</footnote>

<footnote>**IRIA, Rocquencourt.</footnote>

<footnote>***IRIA, Rocquencourt (Détaché C.N.R.S.).</footnote>
## Introduction

Dans cet article, nous étudions les problèmes de contrôle optimal de systèmes dont l'état est la première fonction propre d'un opérateur elliptique. De tels systèmes interviennent dans la modélisation du fonctionnement des réacteurs nucléaires.

La première partie est un exposé, sur un exemple type, de la méthode de résolution d'un tel problème. Après avoir montré l'existence d'un contrôle optimal, on calcule, au moyen d'un état adjoint, le gradient de la fonctionnelle à minimiser. Ces résultats, qui s'appliquent aussi aux problèmes discrets associés permettent une résolution pratique.

Les deuxième et troisième parties sont constituées de deux applications numériques dans lesquelles on met en oeuvre la méthode des éléments finis et la méthode de la descente suivant la projection du gradient.

Enfin la dernière partie est consacrée à la présentation des algorithmes utilisés dans le calcul des valeurs et fonctions propres ainsi que dans la résolution des systèmes linéaires singuliers.
# I) Résultats théoriques; Exemple type

## 1) Hypothèses et Notations

On désigne par $\Omega$ un ouvert borné de $\mathbb{R}^n$ et par $V$ un espace de Hilbert muni du produit scalaire $a(.,..)$ tel que $V \subset L^2(\Omega)$, l'injection étant compacte et dense. Si $V'$ désigne le dual de $V$ et si $L^2(\Omega)$ est identifié à son dual, on a :

$$V \subset\subset L^2(\Omega) \subset\subset V'$$

l'isomorphisme $A$ entre $V$ et $V'$ associé au produit scalaire $a(.,..)$ ($\mathbb{V}u, v \in V$, $a(u,v) = \langle Au, v \rangle$) définit dans $L^2(\Omega)$ un opérateur auto-adjoint, noté encore $A$, non borné, d'inverse compact avec

$$D(A) = \{u \mid u \in V, Au \in L^2(\Omega)\}$$

## 2) Position du problème :

Pour $\rho \in L^\infty(\Omega)$, $\rho(x) \ge \alpha > 0$, on désigne par $\mu_\rho$ la plus petite valeur propre du problème :

$$(1-1) \qquad \text{Au} = \mu_{\rho}u$$

On note $E_\rho$ le sous-espace propre associé à $\mu_\rho$ et on suppose $\mu_\rho$ simple, c'est à dire que $F_\rho$ est de dimension 1.

L'ensemble des contrôles admissibles est le convexe $C$

$$(1-2) \quad C = \{\rho \mid \rho \in L^\infty(\Omega); 0 < \alpha(x) \le \rho(x) \le \beta(x) < +\infty\}$$

$\alpha(x)$, $\beta(x)$ étant données dans $L^\infty(\Omega)$.

Soient $v \in V$ et $a \in \mathbb{R} - \{0\}$. On suppose que pour tout $\rho$ de $C$, $E_\rho$ n'est pas contenu dans le noyau de $v$. Alors l'état $(u_\rho, \mu_\rho)$ du système est solution du problème :

$$(1-3) \qquad \begin{cases} \text{Au}_\rho = \mu_\rho \rho u_\rho, \\ \langle v, u_\rho \rangle = a. \end{cases}$$
On définit la fonction coût $J(\rho)$ par :

$$ (1-4) \qquad J(\rho) = \frac{1}{2} \int_{\Omega} (u_{\rho}(x) - \bar{u}_{\rho})^2 dx \quad \text{avec} $$

$$ (1-5) \qquad \bar{u}_{\rho} = \frac{1}{\text{mes}(\Omega)} \int_{\Omega} u_{\rho}(x) dx. $$

On considère alors le problème de contrôle suivant :

$$ (1-6) \qquad \begin{cases} \text{Trouver } \rho_o \in C \text{ tel que} \\ J(\rho_o) = \inf_{\rho \in C} J(\rho) \end{cases} $$

**Remarque 1-1 :** Afin de simplifier, on ne fait pas intervenir la valeur pro-
pre $\mu_\rho$ dans la fonctionnelle $J$, mais les démonstrations s'appliquent
au cas de fonctionnelles du type

$$ J_1(\rho) = \frac{1}{2} \int_{\Omega} (u_{\rho}(x) - \bar{u}_{\rho})^2 dx + f(\mu_{\rho}) $$

ou du type

$$ J_2(\rho) = \frac{1}{2} \int_{\Omega} (u_{\rho}(x) - z_d)^2 dx \quad f(\mu_{\rho}) $$

avec $f$ fonction positive et $z_d$ fonction de $L^2(\Omega)$.

### 3) Un résultat d'existence :

On a le :

**Théorème 1-1 :** *Le problème (1-6) admet un contôle optimal* $\rho_o$

C'est une conséquence immédiate du :

**Lemme 1-1 :** *L'application* $\rho \to (u_\rho, \mu_\rho)$ *est continue de* $C \subset L^\infty(\Omega)$ *muni de la topologie faible étoile dans* $V \times R$.

En effet soit $\rho_n$ une suite de $C$ convergeant vers $\rho$
dans $L^\infty(\Omega)$ faible étoile, alors la suite $\mu_n$ est bornée car d'après le principe du minimax de COURANT-HILBERT [1].
$$
\mu_{\rho_n} = \mu_n = \min_{u \in V} \frac{a(u, u)}{\int_{\rho_n}^2 dx}
$$

De même la suite $v_n = \frac{u_n}{\|u_n\|_{L^2}}$ ($u_\rho = u_n$), est bornée dans V. Alors
la formulation variationnelle de Au = $\mu\rho u$ montre qu'il existe une
sous-suite de ($u_n, v_n$) convergeant dans $\mathbb{R} \times V$ faible ($\mathbb{R} \times L^2(\Omega)$
fort) vers ($\mu, v$) solution de :

$$
(1-7) \qquad \int_0^\infty \frac{Av}{v^2} dx = 1.
$$

De plus le principe du minimax entraine que $\mu = \mu_\rho$. Or $u_n = \frac{a}{\langle v, v_n \rangle} v_n$, com-
me $E_\rho \notin \text{Kerv} \langle v, v \rangle \neq 0$, la sous-suite de $(u_n, U_n)$ converge vers $(\mu_\rho, u_\rho)$ l'état
du système associé au contrôle $\rho$. Comme il y a unicité de la limite, toute la
suite $(u_n, U_n)$ converge vers $(\mu_\rho, U_\rho)$.

La convergence forte de $u_n$ vers $u_\rho$ dans $V$ résulte des inéga-
lités suivantes obtenues à partir des relations (1-3) valables pour
$(\rho_n, u_n)$ et $(\rho, u_\rho)$ :

$$
\begin{gather*}
a(u_{\rho} - u_{n}, u_{\rho} - u_{n}) = (\mu_{\rho}u - u_{n}\rho_{n}u_{n}, u - u_{n})_{L^{2}(\Omega)} \\
\|u - u_{n}\|_{V}^{2} \leq |u - u_{n}|_{L^{2}}^{2} \cdot |\mu_{\rho}u - u_{n}\rho_{n}u_{n}|_{L^{2}}^{2}
\end{gather*}
$$

<u>Remarque 12.</u> L'ensemble des contrôles optimaux est un compact de
$L^\infty(\Omega)$ faible étoile.

<u>Remarque 1.3.</u> Il n'y a pas en général unicité car à deux fonctions
ρ homothétiques correspond la même fonction propre.

4) <u>Conditions Nécessaires du premier ordre.</u>

Nous utiliserons les propriétés de différentiabilité de l'application
ρ → (u_ρ, μ_ρ) pour transformer la condition d'optimalité :
$$
(1-8) \qquad J'(\rho_o)(\rho - \rho_o) \ge 0 \qquad \forall \rho \ge C
$$

On pose $K = \ker v$ et on munit $D(A)$ de la norme du graphe. On a le
résultat suivant :

**Théorème 1-2** : l'application $\rho \to (u_\rho, \mu_\rho)$ de $L^\infty(\Omega)$ dans
$(D(A) \cap v^{-1}(a)) \times \mathbb{R}$ est différentiable au sens de Fréchet en tout point
de C. Sa différentielle au point $\rho_o = (u_o, \mu_o)$ :
$(\rho - \rho_o) \to (u, M)$ est définie implicitement par :

$$
(1-9) \quad \left\{
\begin{array}{@{}l@{}}
  \text{A}U - \mu_o \rho_o U - \mu_o [(\rho - \rho_o)u_o - \displaystyle\frac{\int_{\Omega} (\rho - \rho_o) u_o^2 dx}{\rho_o u_o^2 dx} \cdot \rho_o u_o] = 0 \\
  \text{M} = \mu_o \displaystyle\frac{\int_{\Omega} (\rho - \rho_o) u_o^2 dx}{\rho_o u_o^2 dx}
\end{array}
\right.
$$

Démonstration :

L'application F définie par :

$$
(1-10) \quad (u, \mu, \rho) \to Au - \mu\rho u = F(u, \mu, \rho)
$$

de $(D(A) \cap v^{-1}(a)) \times \mathbb{R} \times L^\infty(\Omega)$ dans $L^2(\Omega)$ est $C^\infty$. Montrons que
si $\rho_o \in C$, $\frac{\partial F}{\partial(u,\mu)}(u_o, \mu_o, \rho_o) \in \text{Isom}((D(A) \cap K) \times \mathbb{R}, L^2(\Omega))$.

En effet $\frac{\partial F}{\partial(u,\mu)}(u_o, \mu_o, \rho_o)$ est l'application

$$
(1-11) \quad \begin{cases} (U,M) \to (AU - \mu_o \rho_o U) - M\rho_o u_o \\ (D(A) \cap K) \times \mathbb{R} \to L^2(\Omega) \end{cases}
$$

ce qui s'écrit matriciellement (avec $L^2(\Omega) = (\mathbb{R}u_o)^\perp \oplus \mathbb{R}u_o$)
($\perp$ désigne l'orthogonalité dans $L^2(\Omega)$),
$$
(1-12) \quad \left( \begin{matrix} (\mathbb{R} u_o)^L \\ \oplus \\ (\mathbb{R} u_o) \end{matrix} \right) = \left( \begin{matrix} (A - \mu_o \rho_o) & - \rho_o u_o + \left( \frac{\int \rho_o u_o^2 dx}{\int u_o^2 dx} \right) u_o \\ 0 & - \frac{\int \rho_o u_o^2 dx}{\int u_o^2 dx} \end{matrix} \right) \quad \left( \begin{matrix} D(A) \cap K \\ \oplus \\ \mathbb{R} u_o \end{matrix} \right)
$$

Grâce à (1-12), il suffit de prouver que $(A - \mu_o \rho_o) \in \operatorname{Isom}((D(A) \cap K, (\mathbb{R}u_o))^\perp)$,
or d'une part d'après l'alternative de Freedholm, $(A - \mu_o \rho_o)(D(A)) = (\mathbb{R}u_o)^\perp$,
car $\operatorname{Ker}(A - \mu_o \rho_o) = (\mathbb{R}u_o)$ et d'autre part $u_o \notin D(A) \cap K$. Donc
$A - \mu_o \rho_o / D(A) \cap K$ est bien un isomorphisme.

Alors d'après le théorème des fonctions implicites il existe

d'un voisinage $\mathcal{l}$ de $u_o$ dans $D(A) \cap v^{-1}(a)$, d'un voisinage $\mathcal{v}$ de
$\mu_o$ dans $\mathbb{R}$ et un voisinage $\mathcal{v}^-$ de $\rho_o$ dans $L^\infty(\Omega)$ tel que l'ense-
ble des zéros de $F(u,\mu,\rho)$ appartenant à $\mathcal{l} \times W \times \mathcal{v}^-$ est le graphe
d'une unique application $\Phi, \phi \in \mathcal{C}^\infty(\mathcal{v}, u \times W)$.

$$
(1-13) \quad \rho \to (\phi_1(\rho), \phi_2(\rho)).
$$

De plus $\Phi_1(\rho) = u_\rho$ et $\Phi_2(\rho) = \mu_\rho$. En effet l'application
$\rho \to (u_\rho, \mu_\rho)$ étant continue dans $V \times \mathbb{R}$, il existe $\tilde{V} \subset V$ tel que
pour $\rho \in \tilde{V}$, $(u_\rho, \mu_\rho) \in \mathcal{U} \times \mathcal{W}$ et vérifie $F(u_\rho, \mu_\rho, \rho) = 0$.

Donc d'après l'unicité $u_\rho = \Phi_1(\rho)$ et $\mu_\rho = \Phi_2(\rho)$ $\forall \rho \in \tilde{V}$.

Enfin la différentielle de l'application $\rho \to (u_\rho, \mu_\rho)$ est définie im-
plicitement par

$$
(1-14) \quad \mathrm{AU} - \mu_o \rho_o U - M\rho_o u_o = \mu_o (\rho - \rho_o) u_o
$$

ce qui par projection sur $(\mathbb{R}u_o)$ et $(\mathbb{R}u_o)^\perp$ donne (1-9).
Soit $\rho_o$ un contôle optimal de (1-6). Les équations aux variations (1-9) nous permettent de simplifier (1-8) qui s'écrit encore

$$ (1-15) \quad \int_{\Omega} (u_o - \bar{u}_o) U dx \ge 0 \quad \forall \rho \in C. $$

On introduit l'état adjoint $q$ solution de :

$$ \begin{cases} (1-16) & Aq - \mu_o \rho_o q = (u - \bar{u}_o) - \left( \int (u_o - \bar{u}_o) u_o dx \right)^{\frac{\nu}{a}} \\ (1-17) & q \in K \end{cases} $$

(1-16)(1-17) admet une solution unique. En effet le second membre de (1-16) est orthogonal à $u_o$ dans la dualité $\langle V, V' \rangle$ et $(A - \mu_o \rho_o)$ est un isomorphisme entre $K$ et $(u_o)^\perp_{V'} = \{v' \in V', \langle u_o, v' \rangle_{V,V'} = 0\}$ car $u_o \notin K$.

En utilisant $q$ dans (1-15) on obtient :

$$ \begin{align*} \int_{\Omega} (u_o - \bar{u}_o) U dx &= \langle Aq - \mu_o \rho_o q + (u_o - \bar{u}_o, u_o) \rangle_{a,U} \\ &= \langle q, AU - \mu_o \rho_o U \rangle + \frac{1}{a} \langle u_o - \bar{u}_o, u_o \rangle \langle v, U \rangle \\ &= \langle q, Au - \mu_o \rho_o U \rangle \\ &= \mu_o \int_{\Omega} q [u_o(\rho - \rho_o) - \frac{\int_{\Omega} u_o^2 (\rho - \rho_o) dx}{\int_{\Omega} \rho_o u_o^2 dx} \rho_o U_o] dx. \end{align*} $$

La condition d'optimalité (1-8) est donc équivalente à $(\mu_o)$ étant positif):

$$ (1-18) \quad \int_{\Omega} \left( qu_o - \frac{\int_{\rho_o} q u_o \rho_o dx}{\int_{\rho_o} u_o^2 dx} u_o^2 - \rho_o - \rho_o' \right) dx \ge 0 \quad \forall \rho_o \in C $$

On a alors le théorème :

**Théorème 1-3** : Si $\rho_o$ est un contrôle optimal du problème (1-6) $(u_o, \mu_o, \rho_o, q)$ est solution du système (1-3)(1-16)(1-17)(1-18).
**Remarque 1-4** : la condition (1-17) peut être remplacée par $q \in (\mathbb{R}u_o)^\perp$, on obtient un nouvel état adjoint p solution de (1-16) et (1-17)' : $p \in (\mathbb{R}u_o)^\perp$. La nouvelle condition d'optimalité s'obtient en remplaçant q par p dans 1.18.

**Exemple**

Soient $\Omega$ un ouvert borné de $\mathbb{R}^n$, $\Gamma = \partial\Omega$, $\sigma \in L^\infty(\Omega)$, $\sigma \ge 0$.

On prend
$$V = H^1(\Omega)$$
$$a(u,v) = \int_{\Omega} (\operatorname{grad} u \cdot \operatorname{grad} v + uv) + \int_{\Gamma} \sigma uv \, d\Gamma$$

Pour $\rho \in L(\Omega)$, $\rho \ge 0$, soit $\mu_\rho$ la plus petite valeur propre de l'opérateur $Au = \mu_\rho \rho u$ et $u_\rho$ la solution de :

$$\begin{cases} -\Delta u + u = \mu_{\rho} \rho u \quad \text{dans } \Omega \\ \frac{\partial u}{\partial n} + \sigma u = 0 & \text{sur } \Gamma \\ \int_{\Gamma} u(x) d\Gamma = 1 \end{cases}$$

Soit $o$ un contrôle optimal du problème

$$J(\rho_o) = \min_{\alpha \le \rho \le \beta} J(\rho) \quad \text{avec } J(\rho) = \int_{\Omega} (u_{\rho} - \bar{u}_{\rho})^2 dx$$

On introduit l'état adjoint $q$ ($u_{\rho_o} = u_o$, $\mu_{\rho_o} = \mu_o$)

$$\begin{cases} - \Delta q + q - \mu_o \rho_o q = u_o - \bar{u}_o \quad \text{dans } \Omega \\ \frac{\partial q}{\partial n} + \sigma q = - \int (u_o - \bar{u}_o) u_o \quad \text{sur } \Gamma \\ \int_{\Gamma} q(x) dx = 1 \end{cases}$$

Alors $\rho_o$ vérifie l'inéquation :
$$\forall \rho \in C \left\{ \left( qu_o - \frac{\int_{\Omega} qu_o \, dx}{\rho_o u_o^2} \right) (\rho - \rho_o) dx \ge 0 \right.$$
**Remarque 1-5** Les démonstrations précédentes s'appliquent aussi au cas d'une fonction propre simple quelconque. Le cas de fonctionnelles d'une ou plusieurs valeurs propres a été traité par GARABEDIAN-SCHIFFER [1]. Le contrôle dans les coefficients de plus haut degré (coefficients de A) conduit aux mêmes difficultés que celles rencontrées dans l'étude des problèmes elliptiques. Un contre-exemple est donné dans MURAT [1]. L'étude des variations de valeurs et fonctions propres en fonction du domaine ont été étudiées par de nombreux auteurs (problèmes isopérimétriques) HADAMARD [1], COURANT-HILBERT [1], HERSH-PAYNE SCHIFFER [1].

Nous nous proposons maintenant d'étudier la résolution numérique de problèmes de ce type. Deux exemples seront traités. Bien qu'ils ne soient pas exactement du modèle présenté dans la première partie, les principes de résolution sont les mêmes.

## II Exemple I

### 1) Position du problème

On désigne par :

$\Omega$ un ouvert borné de $\mathbb{R}^n$

$V = H_0^1(\Omega)$ muni du produit scalaire $a(u,v) = \int \operatorname{grad} u \operatorname{grad} v \,dx,$

On a la suite d'injection

$$V \hookrightarrow L^2(\Omega) \xrightarrow{\iota} V' = H^{-1}(\Omega).$$

Le convexe $C$, $(\alpha, \beta \in L^\infty(\Omega))$

$$ (2-1) \qquad C = \{\rho, \rho \in L^\infty(\Omega), 0 < \alpha(x) \le \rho(x) \le \beta(x) \text{ pp}\} $$

constitue l'ensemble des contrôles admissibles.
L'état du système est le couple $(u_\rho, \mu_\rho)$ solution unique (COURANT-HILBERT [1]) de

$$
(2-2) \qquad - \Delta u_{\rho} + \rho u_{\rho} = \mu_{\rho} \quad u_{\rho}
$$

$$
(2-3) \qquad u_{\rho} \ge 0
$$

$$
(2-4) \quad \int_{\Omega} u_{\rho}^{2}(x) dx = 1
$$

On définit la fonction coût suivante :

$$
(2-5) \quad J(\rho) = \frac{1}{2} \int_{\Omega} |u|_{\rho}^{2}(x) - z_{d}(x) |^{2} dx, \quad z_{d} \text{ donné dans } L^{2}(\Omega)
$$

On a alors le problème de contrôle optimal

$$
(2-6) \quad \begin{cases} \text{Trouver } \rho_o \in C \text{ tel que} \\ J(\rho_o) = \inf_{\rho \in C} J(\rho) \end{cases}
$$

2) Enoncé des résultats :

Théorème 2-1 : Sous les hypothèses précédentes, le problème (2-6) admet une solution.

Ce théorème est la conséquence du

Lemme 2-1 : L'application $\rho \mapsto (u_\rho, \mu_\rho)$ est continue de $L^\infty(\Omega)$ faible étoile dans $V \times \mathbb{R}$.

La preuve est analogue à celle du lemme 1-1.

En général la solution de (2-6) n'est pas unique (deux contrôles dif-
férent d'une constante définissent la même fonction propre), mais dans
le cas où $z_d \ge 0$ on a la propriété suivante :
**Théorème 2-2.** On se place sous les hypothèses du théorème 2-1, avec

$$z_d \ge 0,$$

(Soit $\mu \in \mathbb{R}$ et $C_\mu = \{\mu \in \mathbb{C}, \mu_\rho = \mu\}$. Alors il existe une unique solution $\rho_1$ du problème :

$$ (2-7) \quad \left\{ \begin{array}{l} J(\rho_1) = \inf_{\rho \in C_\mu} J(\rho), \\ \rho_1 \in C_\mu \end{array} \right. $$

Preuve :

i) **existence** : l'application $\rho \to \mu_\rho$ étant continue de $L^\infty(\Omega)$ faible étoile dans $\mathbb{R}$, $C_\mu$ est un compact de $L^\infty(\Omega)$ faible étoile; donc d'après le lemme 2-1 le problème (2-7) admet une solution.

ii) **Unicité** : Soient $\rho_1$ et $\rho_2$ deux solutions de (2-7), $u_1 = u_{\rho_1}$ et $u_2 = u_{\rho_2}$, supposons que $\rho_1 \neq \rho_2$ donc $u_1 \neq u_2$. Nous allons construire $\rho_o, \rho_o \in C_\mu$ tel que $J(\rho_o) < J(\rho_1) = J(\rho_2)$. Soit

$$ (2-8) \quad \rho_o(x) = \frac{\rho_1 u_1 + \rho_2 u_2}{u_1 + u_2} $$

Alors $\rho_o \in C_\mu$ en effet $\alpha \le \rho_o(x) \le \beta$ donc $\rho_o \in C$

et $(u = \frac{u_1 + u_2}{\|u_1 + u_2\|}, \mu)$ est solution de (2-2)(2-3)(2-4) avec $\rho = \rho_o$.

L'égalité $J(\rho_1) = J(\rho_2)$ implique,

$$ (2-9) \quad \int_{\Omega} z_d(x) u_1(x) dx = \int_{\Omega} z_d(x) u_2(x) dx, $$

donc
$$
\begin{align*}
2 J(\rho_o) &= \int_{\Omega} \left(z_d - \frac{u_1 + u_2}{\|u_1 + u_2\|}\right)^2 dx \\
2 J(\rho_o) &= 1 + \int_{\Omega} z_d^2 dx - 2 \int_{\Omega} z_d^2 \frac{u_1 + u_2}{\|u_1 + u_2\|} dx \\
(2-10) \qquad 2 J(\rho_o) &= 1 + \int_{\Omega} z_d^2 dx - \frac{4}{\|u_1 + u_2\|_L^2} dx
\end{align*}
$$

Or $u_1$ et $u_2$ n'étant pas colinéaires $\|u_1 + u_2\|_L^2 < 2$, donc on obtient à partir de (2.10) :

$$
(2-11) \quad 2 J(\rho_o) < 1 + \int_{\Omega} z_d^2 dx - 2 \int_{\Omega} z_d u_1 dx = 2 J(\rho_1)
$$

ce qui démontre le théorème.

Un raisonnement analogue conduit au résultat :

Corollaire 2-1 : Il existe un unique $\rho_o$ solution du problème

$$
(2-12) \quad \left\{
\begin{array}{ll}
J(\rho_o) = \displaystyle\ Inf_{\rho \in C} J(\rho) & \\
\rho_o \in C & \\
\mu_{\rho_o} = \displaystyle\ Inf_{\rho} \mu_{\rho} & \\
\{\rho, \rho \in C; J(\rho) = J(\rho_o)\} &
\end{array}
\right.
$$

3) Une méthode de régularisation :

Afin d'obtenir une borne convergence des algorithmes, nous rempla-
çons le problème (2-12) par une suite de problèmes perturbés ($P_\epsilon$).

Soit $\epsilon > 0$. Posons :

$$
(2-13) \quad J_{\epsilon}(\rho) = J(\rho) + \epsilon \mu_{\rho}
$$

On démontre comme dans le théorème 2-1, qu'il existe $\rho_\epsilon$ solution du problème $(P_\epsilon)$ suivant :
$$
(2-14) \quad \left\{
\begin{array}{l}
(\rho_{\varepsilon}) \qquad J(\rho_{\varepsilon}) + \varepsilon \mu_{\rho_{\varepsilon}} = \operatorname{Inf}_{\rho \in C} (J_{\varepsilon}(\rho)) \\[1em]
\qquad \qquad \rho_{\varepsilon} \in C.
\end{array}
\right.
$$

On a alors le :

Théorème 2-3 : la suite $\rho_\epsilon$ converge vers $\rho_o$ solution de (2-12) dans $L^\infty$ faible étoile.

Preuve

Posons $\mu_{\rho_{\epsilon}} = \mu_{\epsilon}$ et $\mu_{\rho_{o}} = \mu_{o}$

D'après les définitions de $\rho_o$ et $\rho_\epsilon$,

$$
J(\rho_o) \le J(\rho_\epsilon) \le J(\rho_\epsilon) + \epsilon \mu_\epsilon \le J(\rho_o) + \epsilon \mu_o
$$

Donc

$$
(2-15) \quad 0 \le J(\rho_{\epsilon}) - J(\rho_{o}) \le \epsilon \mu_{o}
$$

Donc $J(\rho_{\epsilon})$ tend vers $J(\rho_{o})$ quand $\epsilon$ tend vers zéro et la suite $\rho_{\epsilon}$ est une suite minimisante du problème (2-12).

D'autre part on déduit de

$$
\epsilon(\mu_{\epsilon} - \mu_{o}) \le J(\rho_{o}) - J(\rho_{\epsilon}) \le 0
$$

que $\mu_{\epsilon} \le \mu_{o}$.

Soit $\rho_o^\sim$ une valeur d'adhérence dans $L^\infty(\Omega)$ faible étoile de la suite
$\rho_\epsilon$. D'une part on a $J(\rho_o^\sim) = J(\rho_o)$ et d'autre part l'application $\rho \to \mu_\rho$
étant continue (lemme 2-1) il vient $\mu_\rho^\sim \le \mu_o$. Il en résulte, grâce
au corollaire 2-1 que $\mu_\rho^\sim = \mu_o$ et $\rho_o^\sim = \rho_o$

Conditions Nécessaires du premier ordre.

semblable
Un raisonnement à celui de la première partie conduit au
résultat suivant : (La sphère unité de $L^2(\Omega)$ est notée $S^1$ et
l'orthogonalité dans $L^2(\Omega)$ , ).
**Théorème 2-4** : L'application $\rho \to (u_\rho, \mu_\rho)$ de $L^\infty(\Omega)$ dans $(V \cap S^1) \times \mathbb{R}$ est différentiable au sens de Fréchet en tout point, $\rho$ de $L^\infty(\Omega)$.

Sa différentielle au point $\rho_o$ ($\mu_o = \mu_{\rho_o}$ ; $u_o = u_{\rho_o}$) est une application de $L^\infty(\Omega)$ dans $(V \cap (\mathbb{R}u_o)) \times \mathbb{R}$ :

$$ (\rho - \rho_o) \to (U, \Lambda) $$

définie implicitement par :

$$ (2-16) \qquad -\Delta U + \rho_o U - \lambda_o U = - (\rho - \rho_o) \bar{u}_o + \left( \int_{\Omega} (\rho - \rho_o) u_o^2 dx \right) u_o $$

$$ (2-17) \qquad \Lambda = \int_{\Omega} (\rho - \rho_o) u_o^2 dx $$

Introduisons alors l'état adjoint $q_\epsilon$, qui est l'unique solution de :

$$ (2-18) \qquad -\Delta q_\epsilon + \rho_\epsilon q_\epsilon - \mu_\epsilon q_\epsilon = z_\epsilon - \left( \int_\Omega z_\epsilon u_\epsilon dx \right) u_\epsilon $$

$$ (2-19) \qquad \int_{\Omega} q_{\epsilon} u_{\epsilon} dx = 0 $$

La condition d'extrémalité

$$ (2-20) \qquad J'(\rho_\epsilon)(\rho - \rho_\epsilon) \geq 0 \quad \forall \rho_\epsilon \in C $$

qui, explicitée, s'écrit

$$ (2-21) \qquad J'(\rho_\epsilon) (\rho - \rho_\epsilon) + \epsilon \frac{\partial \mu_\epsilon}{\partial \rho_\epsilon} (\rho - \rho_\epsilon) \geq 0, \forall \rho_\epsilon \in C, $$

devient grâce à (2-18) et (2-19) :

$$ (2-22) \qquad \frac{1}{2} \int_{\Omega} (\rho - \rho_{\epsilon}) u_{\epsilon} q_{\epsilon} dx + \epsilon \int_{\Omega} (\rho - \rho_{\epsilon}) u_{\epsilon}^2 dx \geq 0, \quad \forall \rho_{\epsilon} \in C $$

**Remarque 2-1** Dans le cas du problème (2-6), on définit un état adjoint par le système (2-18) (2-19) en remplaçant $\rho_\epsilon$, $\mu_\epsilon$, $u_\epsilon$ par $\rho_o$, $\mu_o$, $u_o$. Ceci conduit à la condition d'optimalité

$$ (2-23) \qquad \int_{\Omega} u_o q_o (\rho - \rho_o) dx \geq 0, \quad \forall \rho_o \in C $$
Dans le cas où $\alpha(x)$ et $\beta(x)$ sont des constantes nous avons la :

**Proposition 2-1.** Si $z_d$ n'est pas vecteur propre d'un opérateur $(\Delta - \rho I)$

($\rho \in \mathbb{C}$), tout contrôle optimal de (2-6) est quasi bang-bang c'est-à-dire :

$$ \inf \text{ ess } \rho(x) = \alpha \quad \sup \text{ ess } \rho(x) = \beta $$

### Preuve

Supposons par exemple que $\inf \text{ ess } \rho_o(x) > \alpha$, il existe donc $\eta$, $\eta > 0$ tel que $(\rho_o - \eta)$ soit intérieur à C et $(\rho_o - \eta)$ est encore contrôle optimal de (2-6). La relation d'optimalité (2-23) implique :

$$ u_o(x) \cdot q_o(x) = 0 $$

comme $u_o(x) > 0$ pp. il en résulte que $q_o(x) = 0$ p.p. donc $z_d \neq (\int_{\Omega} z_d u_o dx) u_o$ ce qui contredit l'hypothèse faite sur $z_d$.

### 4) Approximation du problème régularisé

Dans les applications numériques nous avons choisi pour $\Omega$ un domaine rectangulaire de $\mathbb{R}^2$ et utilisé des éléments finis triangulaires

#### 4-1) Approximation de V et $L^\infty(\Omega)$

On approche V par un espace de dimension finie $V_h$ de fonctions continues linéaires sur chaque triangle et nulles au bord
Soit $(w_i)_{i=1}^m$ une base de $V_h$.

$L^\infty(\Omega)$ est approché par l'espace $W_h$ des fonctions constantes sur chaque triangle ou plus finement par l'espace des fonctions constantes sur chacun des domaines $D_i$ définis de la façon suivante :

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

$D_i$ "entoure" le noeud interne $M_i$ et sa frontière est obtenue en joignant les milieux des côtés des triangles afférents à $M_i$ aux centres de gravité correspondants :

Pour une étude détaillée de la méthode des éléments finis nous renvoyons à RAVIART [1].

## 4.2. - Le problème approché ($P_{\epsilon h}$)

On considère le problème aux valeurs propres suivant :

$$ \text{pour } \rho_h \in W_h \text{ trouver } (y_h, \lambda_h) \in V_h \times \mathbb{R} \text{ tel que} \tag{2.24} $$
$$ \forall v_h \in V_h, a(y_h, v_h) + \int_{\Omega} \rho_h y_h v_h \, dx = \lambda_h \int_{\Omega} y_h v_h \, dx. $$

(2.24) équivaut au système matriciel :

$$ \text{(2.25)} \quad \left\{ \begin{array}{l} (\tilde{y}_h, \lambda_h) \in \mathbb{R}^m \times \mathbb{R} \\ A_h \tilde{y}_h = \lambda_h B_h \tilde{y}_h, \quad \text{où } m = \dim V_h, \quad y_h = \sum_{i=1}^m \tilde{y}_{h_i} w_i, \\ \text{et } A_h, B_h \text{ sont des matrices symétriques définies positives.} \end{array} \right. $$

Donc le problème aux valeurs propres généralisé (2.25) admet une plus petite valeur propre strictement positive $\mu_h$,
$$
(2.26) \qquad \mu_h = \min_{\substack{v_h \\ v_h \in V_h}} \frac{a(v_h, v_h) + \int_{\Omega} \rho_h v_h^2 dx}{\int_{\Omega} v_h^2 dx}
$$

On peut choisir l'élément fini de référence et la numérotation des
noeuds de la triangulation de façon à ce que $A_h$ soit irréductible et ait
tous ses éléments non diagonaux négatifs ou nuls. Dans ce cas $A_h^{-1}$ est positive
(VARGA [1]). $B_h$ est toujours positive.

Par suite d'après un théorème de Perron, la matrice positive $A_h^{-1}$ $B_h$
admet un vecteur propre positif unique à une constante multiplicative près ;
la valeur propre correspondante est la plus grande valeur propre, c'est-à-dire
$\mu_h^{-1}$, et est simple :

Soit $\hat{u}_h = (\hat{u}_{h,i})_{i=1,...,m}$ l'unique vecteur propre positif tel que
$$
u_h = \sum_{i=1}^{m} \hat{u}_{h,i} w_i \in V_h \quad \text{ait une norme } L^2(\Omega) \text{ égale à 1.}
$$

Le couple $(u_h, \mu_h) \in V_h \times \mathbb{R}$ est l'état du système approché et
$\tilde{u}_h$ est solution du système :

$$
(2.27) \qquad \left\{
\begin{array}{ll}
\hat{u}_h \in \mathbb{R}^m \\
A_h \cdot \hat{u}_h = \mu_h B_h \quad \hat{u}_h \\
\hat{u}_h \ge 0 \\
\| \hat{u}_h \| = \left\| \sum_{i=1}^m \hat{u}_{h,i} w_i \right\|_{L^2(\Omega)} = 1
\end{array}
\right.
$$

L'ensemble des contrôles admissibles est :

$$
(2.28) \qquad C_h = C \cap W_h
$$

La fonctionnelle coût approché s'écrit :
$$
(2.29) \qquad J_{\epsilon,h}(\rho_h) = \frac{1}{2} \int_{\Omega} (u_{\rho_h} - z_d)^2 dx + \epsilon \mu_{\rho_h}.
$$

On obtient alors le problème de contrôle ($P_{\epsilon,h}$)

$$
(2.30) \quad \left\{
\begin{array}{ll}
\text{Trouver } \rho_{\epsilon,h} \in C_h \text{ tel que :} & \\
\qquad J_{\epsilon,h}(\rho_{\epsilon,h}) = \displaystyle\min_{\rho_h \in C_h} J_{\epsilon,h}(\rho_h) &
\end{array}
\right.
$$

On obtient l'existence d'un contrôle optimal, comme au théorème 2.1. En général il n'y a pas unicité. On a le résultat de convergence suivant :

Théorème 2.5. - Soit $\rho_{\epsilon,h}$ une solution de (2.30). Lorsque $h$ tend vers zéro, la suite $(\rho_{\epsilon,h})$ converge vers une solution $\rho_{\epsilon}$ du problème (2.14) dans $L^{\infty}(\Omega)$ faible étoile.

La démonstration s'articule comme suit.

Soient $\mu_{\epsilon,h}$ et $u_{\epsilon,h}$ la valeur propre et le vecteur propre correspondant à $\rho_{\epsilon,h}$. On établit des majorations a priori des $(\rho_{\epsilon,h})$ $(u_{\epsilon,h})$ et $(\mu_{\epsilon,h})$. Alors on peut extraire des sous-suites convergentes :

$$
\left\{
\begin{array}{l}
\rho_{\epsilon h} \to \rho_{\epsilon}^* \text{ dans } L^{\infty}(\Omega) \text{ faible étoile} \\
u_{\epsilon,h} \to u_{\epsilon}^* \text{ dans } V \text{ faible et } L^2(\Omega) \text{ fort} \\
\mu_{\epsilon,h} \to \mu^* \text{ dans } \mathbb{R}^+
\end{array}
\right.
$$

Puis, en utilisant des résultats de densité, on montre que

$$
u_{\epsilon}^{*} = u_{\rho_{\epsilon}}^{*} \quad \text{et} \quad \mu_{\epsilon} = \mu_{\rho_{\epsilon}}^{*}
$$

On en déduit
$$ \lim_{h \to 0} J_{\epsilon h}(\rho_{\epsilon h}) = J_{\epsilon}(\rho_{\epsilon})^* $$

Enfin la densité de l'espace $U = W_h$ dans $L^\infty(\Omega)$ implique :
$$
J_{\varepsilon}(\rho^{*}) = \underset{\rho \in C}{\operatorname{Min}} J_{\varepsilon}(\rho) .
$$

V - CALCUL DU GRADIENT DE Jεh.-

Comme précédemment, $\perp$ désigne l'orthogonalité dans $L^2(\Omega)$. On désigne par $(u_{h_o}, \mu_{h_o})$ l'état du système approché associé à $\rho_{h_o}$. Par une méthode analogue au théorème 1.2. nous démontrons le :

**Théorème 2.6.-** La fonctionnelle $J_{\epsilon h}$ définie de $W_h$ dans $R$ est différentiable et sa dérivée au point $\rho_{h_o}$ est donnée par :

$$
(2.31) \quad J'_{\epsilon h}(\rho_{h_o}) (\rho - \rho_{h_o}) = \int_{\Omega} u_{h_o} q_{h_o} (\rho - \rho_{h_o}) \, dx + \epsilon \int_{\Omega} u_{h_o}^2 (\rho - \rho_{h_o}) \, dx,
$$

où $q_{h_o}$, l'état adjoint discret, est l'unique solution du système :

$$
(2.32) \quad \left\{
\begin{array}{@{}l@{}}
  a(q_{h_o}, v_h) + \int_{\Omega} \rho_{h_o} q_{h_o} v_h \, dx - \mu_{h_o} \int_{\Omega} q_{h_o} v_h \, dx = \int_{\Omega} z_{d,h_o}^1 v_h \, dx \\
  \qquad \forall v_h \in V_h .
\end{array}
\right.
$$

$$
(2.33) \quad \left\{ \begin{array}{l} q_{h_o} \in V_h \cap (R u_{h_o})^{\perp}, \\ \phantom{q_{h_o} \in } v_h \in V_h. \end{array} \right.
$$

avec

$$
(2.34) \quad \int_{\Omega} z_{d,h_o}^1 v_h \, dx = \int_{\Omega} \left( z_d - \int_{\Omega} z_d u_{h_o} \, dx \right) v_h \, dx .
$$

Sous forme matricielle (2.32) - (2.33) s'écrit :

$$
(2.35) \quad \left\{
\begin{array}{l}
A_{h_o} \hat{q}_{h_o} - \mu_{h_o} B_{h_o} \hat{q}_{h_o} = \left( \int_{\Omega} z_{d,h_o}^1 \cdot w_i \right)_{i=1}^m \\
q_{h_o} = \sum_{i=1}^m \hat{q}_{h_o i} w_i \perp u_{h_o} = \sum_{i=1}^s \hat{u}_{h_o i} w_i
\end{array}
\right.
$$

$$
(2.36) \quad \left\{
\begin{array}{l}
q_{h_o} = \sum_{i=1}^{m} \hat{q}_{h_o i} w_i \perp u_{h_o} = \sum_{i=1}^{s} \hat{u}_{h_o i} w_i \\
\phantom{q_{h_o} = } \qquad A_{h_o} \hat{q}_{h_o} = \sum_{i=1}^{m} \hat{q}_{h_o i} w_i \perp B_{h_o} \hat{u}_{h_o} = \sum_{i=1}^{s} \hat{u}_{h_o i} w_i
\end{array}
\right.
$$

où $\hat{q}_{h_o}$ et $\hat{u}_{h_o}$ sont les vecteurs des composantes de $q_{h_o}$ et $u_{h_o}$ dans
la base $(w_i)_{i=1}^m$ et $A_{h_o}$ et $B_{h_o}$ les matrices symétriques définies positives
introduites précédemment.

## VI - RESULTATS NUMERIQUES.

Pour le calcul numérique effectif nous nous sommes limités au cas d'un ouvert rectangulaire $[0, d_1[x]] \cup [0, d_2]$ et nous avons utilisé des éléments finis triangulaires de degré 1.

Les triangulations sont du type de celle de la fig. 1 qui indique également comment sont numérotés des noeuds. On note i le numéro du noeud $M_i$.

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

### 6.1. - Exemple test :

Supposons qu'il existe une fonction $\rho_t$ de C telle que :

$$ z_d = u(\rho_t) \quad \text{et} \quad \left\{ \begin{array}{ll} \text{Inf} & \text{ess } \rho_t(x) = \alpha \\ x \in \Omega \\ \\ \text{Sup} & \text{ess } \rho_t(x) = \beta \\ x \in \Omega \end{array} \right. $$

Alors il est facile de voir que $\rho_t$ est la solution unique du
problème (2.6.) :

Pour

$$ \rho_t = \beta = 10 \text{ sur } ]0, \frac{d_1}{2}[\times]0, d_2[ $$

$$ \rho_t = \alpha = 0 \text{ sur } ]\frac{d_1}{2}, d_1[\times]0, d_2[ $$

$$ d_1 = d_2 = 2.76 $$

nous avons les résultats suivants :

$$ \begin{array}{l} \text{contrôle initial: } \rho_o = 0 \text{ sur } \Omega \quad \left\{ \begin{array}{l} \text{Après une itération} \\ \text{de gradient} \end{array} \right. \\ \text{coût initial : } J_o = 0.32 \qquad \text{Contrôle final : } \rho_f \\ \hfill (\text{cf. tableau T1}) \end{array} $$

<table><caption>Tableau T1</caption><thead><tr><th>i</th><th>ρ<sub>f</sub></th><th>ρ<sub>t</sub></th><th>i</th><th>ρ<sub>f</sub></th><th>ρ<sub>t</sub></th></tr></thead><tbody><tr><td>1</td><td>6.80</td><td>10</td><td>11</td><td>10</td><td>10</td></tr><tr><td>2</td><td>10</td><td>10</td><td>12</td><td>10</td><td>10</td></tr><tr><td>3</td><td>10</td><td>10</td><td>13</td><td>10</td><td>10</td></tr><tr><td>4</td><td>10</td><td>10</td><td>14</td><td>10</td><td>10</td></tr><tr><td>5</td><td>0</td><td>10</td><td>15</td><td>10</td><td>10</td></tr><tr><td>6</td><td>0</td><td>0</td><td>16</td><td>0</td><td>0</td></tr><tr><td>7</td><td>0</td><td>0</td><td>17</td><td>0</td><td>0</td></tr><tr><td>8</td><td>0</td><td>0</td><td>18</td><td>0</td><td>0</td></tr><tr><td>9</td><td>0</td><td>0</td><td>19</td><td>0</td><td>0</td></tr><tr><td>10</td><td>0</td><td>0</td><td>20</td><td>0</td><td>0</td></tr></tbody></table>

Nous rappelons que i est le numéro du noeud n_i de la triangulation.
## 6.2. - "Aplatissement optimal".-

On prend pour $z_d$ la fonction constante $\gamma$ tel que :

$$ \int_{\Omega} z_{d}^{2} dx = 1 $$

$$ d_{1} = d_{2} = 2.76 $$

$$ \alpha = 0, \beta = 30 $$

$$ \text{Contrôle initial : } \rho_{o} = 0 $$

$$ \text{Coût initial : } J_{o} = 0.38 $$

Les variations du coût en fonction du nombre d'itérations de gradient sont représentées dans le tableau suivant :

<table><caption>Tableau T2</caption><thead><tr><td>Nombre d'itérations de gradient</td><td>Coût</td><td>Temps de calcul<br/>10 070</td></tr></thead><tbody><tr><td>1</td><td>0.24</td><td></td></tr><tr><td>2</td><td>0.22</td><td></td></tr><tr><td>5</td><td>0.17</td><td></td></tr><tr><td>8</td><td>0.166</td><td></td></tr><tr><td>10</td><td>0.165</td><td></td></tr><tr><td>12</td><td>0.163</td><td>3 heures</td></tr></tbody></table>

Après 12 itérations de gradient le contrôle $\rho_f$ a l'allure suivante :
<figure><img src="image_4.png" /></figure>

Les numéros indiqués en abscisse et en ordonnée sont les numéros des noeuds de la triangulation.

Les figures 3 et 4 ci-dessous représentent une section parallèle à l'axe des x et une section selon la diagonale secondaire de Ω des fonctions $z_d, u_h(\rho_o)$ et $u_h(\rho_f)$.

Les valeurs propres correspondantes sont :

$$ \lambda_h(\rho_o) = 2.67 \qquad \lambda_h(\rho_f) = 12.14 $$

Ces résultats ont été obtenus avec une triangulation comportant 400 noeuds intérieurs.
<figure><img src="image_5.png" /><figcaption>Fig. 4 : Diagonale secondaire de Ω</figcaption></figure>
### III - DEUXIEME EXEMPLE - APPLICATION A UN PROBLEME DE CONTROLE PONCTUEL.-

#### 1) Position du problème.-

Soient $\Omega = ]0, R[ \text{ un ouvert de } R \ (R<+\infty)$

$E$ un compact de $\Omega$

$b$ un point de $E$

On considère le système suivant :

$$ (3.1.) \qquad \begin{cases} - M^2 \Delta u + \frac{2M^2}{\lambda} u(b) \delta(x-b) + u = \mu \Phi(x) u \\ \frac{\partial u}{\partial n}(0) = 0 \end{cases} $$

$$ (3.2.) \qquad \begin{cases} \frac{\partial u}{\partial n}(0) = 0 \\ \frac{\partial u}{\partial n} \end{cases} $$

$$ (3.3.) \qquad \begin{cases} u(R) = 0, \text{ où } \Phi(x) \text{ est une fonction positive, continue par morceaux.} \\ u(0) = 1. \end{cases} $$

Le point $b$ joue le rôle de contrôle.

L'état du système sera alors la fonction propre associée à la plus petite valeur propre positive $\mu$ telle que :

$$ (3.4.) \qquad u(0) = 1. $$

On définit la fonctionnelle :

$$ (3.5.) \qquad J(b) = \int_{\Omega} (u-\bar{u})^2 x \, dx \quad \text{avec} \quad \bar{u} = \frac{1}{R^2} \int_{\Omega} u \, x \, dx $$

On considère le problème de contrôle :

$$ (3.6.) \qquad \begin{cases} \text{Trouver } \bar{b} \in E \text{ tel que :} \\ J(\bar{b}) = \begin{cases} \text{Min } J(b) \\ b \in E \end{cases} \end{cases} $$

**Remarque 3.1.-** Ce problème peut s'interpréter comme la recherche du positionnement optimal de barres absorbantes dans un réacteur nucléaire, $u$ représentant le flux de neutrons.

On démontre que l'application qui a b' associe $u_b$ solution
de (3.1) (3.2) (3.3) (3.4) est continue de $\cdot E$ dans $V = \{\phi \mid \phi \in H^{-1}(\Omega), \phi(R) = 0\}$
On en déduit la proposition :

**Proposition 3.1.** - Le problème (3.6.) admet au moins une solution.

**Orientation** - Nous allons appliquer la théorie présentée au I directement à un problème approché. Notons que la méthode du I ne s'applique pas au problème continu, car le terme $u(b) \delta(x-b)$ ne possède pas de bonnes propriétés de différentiabilité.

2) Le problème approché.

Nous utilisons des éléments finis d'ordre 1. Soit $(w_i)_{i=1,...,n}$ la base d'éléments finis considérés. Si on définit les matrices suivantes :

$$A_{ij} = \int_{\Omega} (M^2 \operatorname{grad} w_i \operatorname{grad} w_j + w_i w_j) + \frac{2M^2}{\lambda} w_i^{(b)} w_j^{(b)}$$

$$B_{ij} = \int_{\Omega} \Phi(x) w_i w_j$$

le problème approché s'écrit :

$$(3.7.) \qquad \left\{ \begin{array}{l} \displaystyle Au_h = \mu B u_h \\[1em] \displaystyle u_h \quad \text{étant la plus petite valeur propre positive.} \end{array} \right.$$

$$(3.8.) \qquad u_h = \sum_{i=1}^{m} u_h^i w_i \quad \text{et} \quad u_h^l = l$$

$$(3.9.) \qquad J_h^{(b)} = \int_{\Omega} \left[ \sum_{i=1}^{n} u_h^i v_i(x) \right]^2 x \, dx$$

$$\text{avec } v_i(x) = w_i(x) - \frac{2}{R^2} \int_{\Omega} w_i(x) \, x \, dx.$$

On note

$$a_{ij} = \int_{\Omega} v_i(x) v_j(x) \, x \, dx \quad \text{alors} \quad J_h^{(b)} = \int_{\Omega} u_h^i u_h^j \, a_{ij}.$$

On a alors le problème de contrôle approché :
$$
(3.10.) \quad \left\{
\begin{array}{l}
\text{Trouver } \bar{b}_h \in E \text{ tel que} \\
J_h(\bar{b}_h) = \min_{b \in E} J_h(b) .
\end{array}
\right.
$$

On démontre les résultats suivants :

*Proposition 3.2.* - Le problème (3.10.) admet une solution (*b*<sub>*h*</sub>, *u*<sub>*h*</sub>, *h*) . Et quand *h* tend vers 0 on a les convergences suivantes :

$$
b_h \to b \text{ dans } E
$$

$$
\mu_h \to \mu \text{ dans } \mathbb{R}^+
$$

$$
u_h \to u \text{ dans } V = \{\phi \mid \phi \in H^1(\Omega) ; \phi(R) = 0\} \text{ faible}
$$

$$
u_h \to u \text{ dans } L^2(\Omega) \text{ fort}
$$

(*b*, *µ*, *g*) étant une solution du problème (3.6.)

Schéma de la démonstration . L'existence de (b<sub>h</sub>, u<sub>h</sub>, μ<sub>h</sub>) est immédiate

Par utilisation du théorème du Min - Max, on a :

$$
\mu_h = \min_{\substack{n \\ \{u=\sum_i u_i \cdot w_i\}}} \frac{\langle \text{Au}, u \rangle}{\mathbb{R}^n}
$$

On en déduit que $|\mu_h| \le K_1$

$$
\text{et} \quad \| |u_h| \|_V^2 \le K_2
$$

($K_1$ et $K_2$ étant deux constantes indépendantes de $h$).

On extrait alors des sous-suites $b_h, y_h, u_h$ telles que :
$$b_h \rightarrow b \text{ dans } E$$

$$\mu_h \rightarrow \mu \text{ dans } \mathbb{R}^+$$

$$u_h \rightarrow u \text{ dans } V \text{ faible}$$

En utilisant la formulation faible du problème, on montre que $(b, \mu, u)$ vérifient (3.1) ... (3.4).

On en déduit enfin que $b$ est un contrôle optimal.

Comme dans le cadre de la théorie développée au paragraphe I, on définit alors pour le problème approché un état adjoint $p$ solution de

$$(3.11) \qquad (A-\mu B) p = - A y + \frac{y(Ay,y)}{(y,y)_{\mathbb{R}^n}}$$

$$(3.12) \qquad (p, By)_{\mathbb{R}^n} = 0$$

Remarque 3.2.- On peut obtenir formellement cet état adjoint en introduisant le Lagrangien suivant :

$$L(u,\mu,b,p,\lambda) = \int_0^t u A u + (p, (A-\mu B)u)_{\mathbb{R}^n} + \sum_i \lambda_i u_i \\ \text{avec } \lambda_i \le 0.$$

Compte tenu de (3.11) et (3.12) la fonctionnelle $b \mapsto J_h(b)$ a pour dérivée :

$$(3.13) \qquad J'_h(b) = \frac{1}{2} \left( p, \frac{\partial A}{\partial b} y \right)_{\mathbb{R}^n}.$$

Pour résoudre numériquement ce problème nous avons utilisé la méthode du gradient. La méthode de résolution du système adjoint est présentée au paragraphe IV.

3) Généralisation.- Le problème lié au positionnement de plusieurs barres conduit au système :
$$ (3.14) \qquad \begin{cases} - M^2 \Delta u + \sum_i \frac{2M^2}{\lambda_i} u(b_i) \delta(x-b_i) + u = \mu \Phi(x) u \\ \frac{\partial u}{\partial n} = 0 \end{cases} $$

$$ (3.15) \qquad \frac{\partial u}{\partial n} = 0 $$

$$ (3.16) \qquad u(R) = 0 $$

On garde la même fonction coût qui dépend donc des $b_i$.

Dans ce cas nous avons appliqué pour la résolution numérique une méthode de relaxation barre par barre, plus précisément, on suppose qu'il y a p barres $b_i$, et on note $G_i$ le gradient de $J$ par rapport à $b_i$. Indiquons l'organigramme du calcul.

1) Initialisations $b_i^1$, $n=0$.

2) $n = n + 1$ ; $i = 0$.

3) $i = i + 1$, on effectue un pas de gradient pour la barre $b_i$

$$ b_i^{n+1} = b_i^n - \rho_n G_i^n $$

$\rho_n$ étant choisi tel que $b_{i-1}^{n+1} < b_i^{n+1} < b_{i+1}^n$.

4) Si $i = p$ on continue, sinon retour en (3).

5) Test final de convergence sur la variation relative du critère, si le test n'est pas vérifié retour en (2).

Sinon fin.

Nous présentons les résultats numériques dans le cas de 1-2-3 et 4 barres à positionner. La fonction $\Phi(x)$ est prise constante par morceaux

$$ \Phi(x) = k_i \quad \text{pour } x \in ]R_{i-1}, R_i[ $$

dans les exemples traités $R_0 = 0$ ; $R_1 = 0,80 \text{ m}$ , $R_2 = 1,13 \text{ m}$ ; $R_3 = 1,38 \text{ m}$ et $E = [R_1, R_2]$, $k_1 = 1$ ; $k_2 = 1,2$ ; $k_3 = 1,2$.
Pour illustrer la convergence de l'algorithme, nous donnons les résultats à chaque itération dans le cas du positionnement de deux barres (dans les 2 premières colonnes sont indiquées les abscisses (en mètres) des barres, la troisième colonne correspond au coût.

<table><thead><tr><th>1ère barre</th><th>2ème barre</th><th>J</th></tr></thead><tbody><tr><td>0,8300</td><td>1,100</td><td>5,2302</td></tr><tr><td>0,8300</td><td>1,047</td><td>4,1968</td></tr><tr><td>0,8916</td><td>1,047</td><td>1,6952</td></tr><tr><td>0,8916</td><td>1,066</td><td>1,6892</td></tr><tr><td>0,9141</td><td>1,050</td><td>1,2759</td></tr><tr><td>0,9143</td><td>1,050</td><td>1,2747</td></tr><tr><td>0,9143</td><td>1,057</td><td>1,2560</td></tr><tr><td>0,9155</td><td>1,057</td><td>1,2487</td></tr></tbody></table>
Positionnement d'une barre

<figure><img src="image_6.png" /></figure>
Positionnement de deux barres

<figure><img src="image_7.png" /></figure>
Positionnement de trois barres

<figure><img src="image_8.png" /></figure>
Positionnement de quatre barres

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

flux de neutrons

Données

$\lambda=0,5$

$M^2=0,005$

$k_1=1$

$k_2=1,2$

$k_3=1,2$
IV - METHODES NUMERIQUES.-

Dans cette partie nous développons les aspects proprement numériques
de ce type de problème. On utilise un algorithme du gradient avec projection.

1) - *Algorithme du gradient avec projection.* -

Étape 1 : Initialisations :

n = 0
choix d'un contrôle admissible x<sub>1</sub>
résolution de l'équation d'état pour le contrôle x<sub>1</sub>.
calcul de la fonction coût pour le contrôle x<sub>1</sub>.

Étape 2 : n = n + 1

calcul de l'état adjoint (résolution d'un système
singulier d'équations linéaires).

Étape 3 : calcul du gradient de la fonction coût soit G_n.

Étape 4 : $\hat{x}_{n+1} = x_n - \rho_n G_n$
($\rho_n$ étant choisi de façon sous optimale)

Étape 5 : Projection du contrôle $\hat{x}_{n+1}$ sur le convexe des contrôles admissibles, soit $x_{n+1}$.

Étape 6 : Calcul de l'état du système associé à $x_{n+1}$
calcul de la fonction coût.

Étape 7 : Test final sur la variation relative du coût.

Si le test n'est pas vérifié retour en 2.

Remarque 4.1.- Il est facile de modifier l'algorithmesci-dessus
pour y inclure la régularisation introduite au paragraphe II.

Remarque 4.2.- On ne sait démontrer aucune propriété de convexité
pour les fonctionnelles à minimiser. On ne peut donc démontrer que la convergence que la convergence vers un minimum local.

Nous allons maintenant préciser les méthodes de résolution du système direct et du système adjoint.

## 2) Calcul des plus petits éléments propres du problème généralisé

$$AX = \mu Bx.$$

Dans ce paragraphe comme dans le suivant $||.||$ représente la norme euclidienne dans $\mathbb{R}^m$ ; (.,.) est le produit scalaire associé.

A et B sont des matrices, à coefficients réels symétriques définies positives de dimensions (mxm). Soit S(A,B) le spectre généralisé du couple (A,B) ; c'est-à-dire l'ensemble des $\mu \in \mathbb{C}$ tels qu'il existe $x \in \mathbb{R}^m$ : $Ax = \mu Bx$.

S(A,B) se compose de m éléments réels strictement positifs.

$$0 < \mu_1 \le \mu_2 \le \dots \le \mu_n.$$

On cherche à calculer $\mu_l$ et un vecteur propre associé $x_l$ (Dans le cas qui nous occupe $\mu_l$ est simple).

Nous nous contentons d'exposer les principes essentiels de la méthode suivie ici. Cette méthode est exposée dans BATEH [1].

## Méthode de la puissance inverse.

Lorsqu'on dispose d'une bonne approximation $\mu_o$ de la valeur propre $\mu$, un moyen pour calculer un vecteur propre associé est la méthode de la puissance inverse avec shift.

Un vecteur initial normalisé $x_o$ étant choisi, on calcule $\bar{x}_l$ solution de :

$$(4.1.) \qquad (A - \mu_o B) \bar{x}_l = B x_o$$
(le choix de $x_o$ est discuté dans WILKINSON [1]).

On pose :

$$
(4.2.) \qquad x_1 = \frac{\bar{x}_1}{||x_1||}
$$

$$
(4.3.) \quad \mu_l = \mu_o + \frac{((A - \mu_o)^B x_l, x_l)}{(B^T x_l, x_l)}
$$

D'une façon générale, connaissant $\bar{x}_k$ on calcule $\bar{x}_{k+1}$ solution
de :

$$
(4.4.) \qquad (A - \mu_o)^B \bar{x}_{k+1} = B \bar{x}_k
$$

Puis

$$
(4.5.) \qquad x_{k+1} = \frac{\bar{x}_{k+1}}{\|\bar{x}_{k+1}\|}
$$

$$
(4.6.) \quad \mu_{k+1} = \mu_o + \frac{((A - \mu_o)^B x_{k+1}, x_{k+1})}{(B^T x_{k+1}, x_{k+1})}
$$

(la suite $(\mu_k, x_k)$ converge vers $(\mu_1, x)$ où $x$ est un vecteur propre normalisé associé à $\mu$ si l'approximation $\mu_o$ est plus proche de $\mu$ que toute autre valeur propre.

Donc ici il suffit de trouver une approximation de $\mu_1$ plus petite
que $\mu_1$. De plus comme $\mu_1$ est simple, on convergera bien, au signe près,
vers l'unique vecteur propre cherché.

Méthode de la sècante. Factorisation L D L^t :

Comme $\mu_1$ est le plus petit zéro du polynôme caractéristique :

$$
(4.7.) \qquad P(\mu) = \det(A - \mu B)
$$
l'idée est de calculer une approximation de $\mu_1$ par la méthode de la sécante. On calcule $P(\lambda)$ en factorisant $(A-\lambda B)$ sous la forme $L D L^t$ où $L$ est une matrice triangulaire inférieure avec des 1 sur la diagonale principale et $D = (d_{ij})$ une matrice diagonale. On a alors :

$$ (4.8.) \qquad P(\mu) = \det D = \prod_{i=1}^{m} d_{ii} $$

Cette factorisation bien que coûteuse est intéressante d'un double point de vue.

i) Si on choisit efficacement les noeuds de la triangulation les matrices $A$ et $b$ sont bandes. La matrice $L$ possède alors exactement la même structure ce qui permet une grande économie de place mémoire.

ii) On peut contrôler le processus de façon à ne jamais "dépasser" $\mu_1$. En effet le nombre d'éléments négatifs ou nuls sur $D$ est exactement égal au nombre de valeurs propres comptées avec leur ordre de multiplicité, inférieure ou égale à $\mu$.

**3) Résolution numérique d'un système creux et singulier d'équations linéaires**
$Cy = z.$-

i) **Systèmes linéaires et problèmes aux moindres carrés.-**

Soit $C$ une matrice réelle symétrique singulière de dimension $(mxm)$.

$z$ est un $m$-vecteur réel. On cherche à résoudre le système :

$$ (4.9.) \qquad \begin{cases} Cy = z \\ y \in \mathbb{R}^m \cap \text{Im}(C) \end{cases} $$

$\text{Im}(C)$ est l'espace image de l'application linéaire défini par $C$.
Le système (4.9.) (4.10.) admet au plus une solution. Toutefois il n'admet évidemment pas de solution si z n'est pas dans Im(C). Dans le cas qui nous occupe, cette condition, bien que théoriquement vraie, n'est pas numériquement assurée. En effet z est orthogonal à un vecteur x qui n'est qu'approximativement dans le noyau de C = A - μB, puisque μ et x sont des approximations des premiers éléments propres du couple (A,B). C'est pourquoi il est préférable d'un point de vue numérique de définir y comme l'unique solution du problème aux moindres carrés sans contraintes théoriquement équivalent au système (4.9) (4.10) :

$$ (4.11.) \qquad ||Cy' - z|| = \min_{t \in \mathbb{R}^m} ||Ct - z|| $$

$$ (4.12.) \qquad y' \in \mathbb{R}^m \cap \text{Im}(C) $$

ce qui est équivalent à

$$ (4.13.) \qquad \begin{cases} \tilde{C} \tilde{y} = \begin{bmatrix} I & C \\ C & 0 \end{bmatrix} \begin{bmatrix} r \\ y' \end{bmatrix} = \begin{bmatrix} z \\ 0 \end{bmatrix} = \tilde{z} \\ \\ \tilde{y} \in \mathbb{R}^{2m} \cap \text{Im}(\tilde{C}) \end{cases} $$

$\tilde{C}$ est une matrice symétrique, singulière de dimension (2m, 2m) et (4.13.)

(4.14.) admet toujours une solution.

Pour résoudre alors ce problème nous utilisons une méthode itérative du type gradient conjugué.
BIBLIOGRAPHIE.

K.J. BATHE [1] - Solution Methods for large generalized eigenvalue Problems in structural engineering. UC SEM 71-2 Structural Engineering Laboratory. University of California, Berkeley, Nov. 71

M. CHICCO [1] - some properties of first eigenvalue and first eigen function of linear second elliptic partial differential equations in divergence forms
(Bulletins U.M.I. t.5. 1972 pp. 245-254).

COURANT-HILBERT [1] - Methods of mathematical physics
(Vol. 1 Interscience Publishers N.Y. 1953).

G. DANDEU [1] - L'algorithme de Le Foll
(Journées d'études sur les équations intégrales
Lab. d'études Math. des signaux. Université de Nice
France - Dec. 1969.)

PR. GARABEDIAN - M. SHIFFER [1] - Variational problems in the theory of elliptic partial differential equations.
J. of. Rat. Mach. and Analysis (1953) 2 pp. 137-171.

HADAMARD [1] - Mémoire sur le problème d'analyse relatif à l'équilibre des plaques élastiques encastrées.
(oeuvres complètes Ed du C.N.R.S. Paris 1968. p. 12).

F. MIGNOT [1] - Contrôle de fonction propre.
(Note au CRAS Série A. t. 280 1975 pp. 333-335).

[2] - Contrôle dans les inéquations elliptiques
(à paraître dans Journal of Functional Analysis).

F. MURAT [1] - Contre-exemples pour divers problèmes où le contrôle intervient dans les coefficients
(à paraître dans Annali di Matematica pura et applicata
F. MURAT - J. SIMON [1] - Sur le contrôle par un domaine géométrique
(Publication du Laboratoire d'Analyse Numérique
Université Paris VI - 1976).

C.C. PAIGE - M.A. SAUNDERS [1] - Solution of sparse indefinite systems of
linear equations. Siam Journal on numerical
analysis - Vol. 12 - n° 4 - Sept. 1975 p. 617-629.

P.A. RAVIART [1] - Méthode des éléments finis
(cours DEA Analyse numérique Paris VI).

SAGUEZ [1] - Contrôle ponctuel et contrôle en nombres entiers
de systèmes distribués (thèse Docteur Ingénieur
Paris VI - 1974).

G.W. STEWART [1] - The numerical treatment of large eigenvalue problems.
Proceedings of IFIP-Congress 74.
Information processing 74 Stockholm. Vol. 3
mathematical aspects of information processing
p. 666-672.

J.P. VAN DE WIELE [1] - Résolution numérique d'un problème de contrôle
optimal de valeurs propres et vecteurs propres
(Thèse 3e cycle Paris VI 1974).

R.S. VARGA [1] - Matrix Iterative Analysis
(Prentice Hall New Jersey 1962).

J.H. WILKINSON [1] - The algebraic eigenvalue problem
(clarendon Press - Oxford 1965).