Control óptimo de una ecuación de Poisson mediante el algoritmo del gradiente con paso óptimo y proyección

Resumen: En este ejemplo se considera un problema de control relativo a la ecuación de Poisson.

Problema considerado

$\Omega\subset\mathbb{R}^N$ $\omega\subset\Omega$ $N=2$ $N=3$ .

$u\,\colon\,\omega \to \mathbb{R}$ solución del problema (de control) siguiente

\begin{matrix} (1) & {\begin{cases} Minimizar J (u) = \frac{a}{2} \int_{Ω} | y - y_{d} |^{2} d x + \frac{b}{2} \int_{ω} | u |^{2} d x \\ sujeto a u \in U_{a d} \subset L^{2} (ω) \end{cases} \end{matrix}

$a,b>0$ $y_d\in L^2(\Omega)$ $y$ es solución del problema

\begin{matrix} (2) & {\begin{cases} - Δ y = u 1_{ω} en Ω \\ y = 0 sobre \partial Ω \end{cases} \end{matrix}

${\cal U}_{ad}$ $u$ control $\eqref{eq.poisson1}$ ecuación de estado $y$ estado asociado a $u$ .

$u$ ${\cal U}_{ad}$ $y$ $\eqref{eq.poisson1}$ $y_d$ sin tener que realizar para ello un esfuerzo excesivo.

$\hat u$ control óptimo $\hat y$ $\hat p$ estado adjunto asociado a $\hat u$ , solución del problema siguiente, denominado sistema de optimalidad:

\begin{matrix} (3) & {\begin{cases} {\begin{cases} - Δ \hat{y} = \hat{u} 1_{ω} en Ω \\ \hat{y} = 0 sobre \partial Ω \end{cases} \\ {\begin{cases} - Δ \hat{p} = \hat{y} - y_{d} en Ω \\ \hat{p} = 0 sobre \partial Ω \end{cases} \\ \hat{u} = P_{a d} (- \frac{a}{b} \hat{p} |_{ω}) \end{cases} \end{matrix}

$\mathbf{P}_{ad}$ $L^2(\omega)$ ${\cal U}_{ad}$ . Tomaremos aquí, para facilitar la resolución,

\begin{matrix} (4) & U_{a d} = {u \in L^{2} (ω) : \underset{―}{u} \leq u \leq \overset{―}{u} c.p.d. en Ω} . \end{matrix}

Con esta elección, la operación de proyección se hace punto a punto:

Algoritmo del gradiente con paso óptimo

$\eqref{eq.minimizar}$ $\eqref{eq.poisson.optimality}$ . Más precisamente, recurriremos a un método de descenso, concretamente el algoritmo del gradiente con paso óptimo:

$u^0\in{\cal U}_{ad}$

$n\geq 0$ $u^n\in{\cal U}_{ad}$ ,

$d^n = \nabla J(u^n)$

$\rho^n$ , solución de

\begin{matrix} (5) & {\begin{cases} Min H (ρ) := J (u^{n} - ρ d^{n}) \\ sujeto a ρ \geq 0 \end{cases} \end{matrix}

$u^{n+1} = \mathbf{P}_{ad} \Big( u^n-\rho^n\, d^n \Big)$

(B.4) Si

\begin{matrix} (6) & ∥ u^{n + 1} - u^{n} ∥_{L^{2} (ω)} \leq ε (o mejor \frac{∥ u^{n + 1} - u^{n} ∥_{L^{2} (ω)}}{∥ u^{n + 1} ∥_{L^{2} (ω)}} \leq ε), \end{matrix}

$u^{n+1}$ como control óptimo. En caso contrario, repetir el paso (B).

Cálculo del paso óptimo

$J\,\colon\, {\cal U}_{ad} \to \mathbb{R}$ $C^1$ ${\cal U}_{ad}$ y que su gradiente viene dado por

\begin{matrix} (7) & \nabla J (u) = a p |_{ω} + b u, \end{matrix}

$p$ $u$ .

$u$ $d$ $\rho$ $H\,\colon\, \mathbb{R}^+ \to \mathbb{R}$ $H(\rho) = J(u - \rho d)$ $y$ $u$ $p$ $u$ $\overline{y}$ $d$ $u-\rho d$ $y-\rho\,\overline{y}$ y

\begin{matrix} (8) & \begin{aligned} H (ρ) & = J (u - ρ d) = \frac{a}{2} \int_{Ω} | y - ρ \overset{―}{y} - y_{d} |^{2} d x + \frac{b}{2} \int_{ω} | u - ρ d |^{2} d x \\ = \frac{a}{2} \int_{Ω} [| y - y_{d} |^{2} - 2 ρ \overset{―}{y} (y - y_{d}) + ρ^{2} | \overset{―}{y} |^{2}] d x + \frac{b}{2} \int_{ω} [| u |^{2} - 2 ρ u d + ρ^{2} | d |^{2}] d x \\ = ρ^{2} [\frac{a}{2} \int_{Ω} | \overset{―}{y} |^{2} d x + \frac{b}{2} \int_{ω} | d |^{2} d x] - ρ [a \int_{Ω} \overset{―}{y} (y - y_{d}) d x + b \int_{ω} d u d x] + \dots \end{aligned} \end{matrix}

$\rho$ .

$H'(\rho)=0$ si y solo si

\begin{matrix} (9) & ρ = \frac{a \int_{Ω} \overset{―}{y} (y - y_{d}) d x + b \int_{ω} d u d x}{[a \int_{Ω} | \overset{―}{y} |^{2} d x + b \int_{ω} | d |^{2} d x]} \end{matrix}

$p$ $u$ , se tiene que

\begin{matrix} (10) & \int_{Ω} \nabla p \cdot \nabla v d x = \int_{Ω} (y - y_{d}) v d x \forall v \in H_{0}^{1} (Ω) \end{matrix}

$v = \overline{y}$ se tiene

\begin{matrix} (11) & \int_{Ω} (y - y_{d}) \overset{―}{y} d x = \int_{Ω} \nabla p \cdot \nabla \overset{―}{y} d x = \int_{Ω} (d 1_{ω}) p d x = \int_{ω} (d 1_{ω}) p d x \end{matrix}

Así pues, se tiene que

\begin{matrix} (12) & \begin{aligned} a \int_{Ω} \overset{―}{y} (y - y_{d}) d x + b \int_{ω} d u d x & = a \int_{ω} (d 1_{ω}) p d x + b \int_{ω} d u d x \\ = \int_{ω} d (a p + b u) d x = \int_{ω} | d |^{2} d x \end{aligned} \end{matrix}

Finalmente, el paso óptimo viene dado por la expresión

\begin{matrix} (13) & ρ = \frac{\int_{ω} | d |^{2} d x}{a \int_{Ω} | \overset{―}{y} |^{2} d x + b \int_{ω} | d |^{2} d x} . \end{matrix}

Código

Recordamos que las formulaciones débiles correspondientes a las ecuaciones de estado y de estado adjunto conducen a la misma forma bilineal. Por tanto, bastará definir un único problema e ir cambiando adecuadamente los segundos miembros en cada caso.


//   ControlCond_FP.edp
//
//   OPTIMAL CONTROL OF A CAPACITOR
//   ALGORITHM: GRADIENT OPTIMAL STEP + PROJECTION
//
//   |  Min J(u) = (a/2)\int_\Omega |y-y_d|^2 + (b/2) \int_\omega |u|^2
//   |  subject to u\in U_{ad} = { u\inL^2(\omega)/ umin <= u <= umax }
//   where
//   |  -\Delta y = u 1_\omega  in \Omega
//   |    y = 0      over \partial\Omega
//
//   Optimality system:  (u,y,p) is the solution to
//   | | -\Delta y = u 1_\omega  in \Omega
//   | |   y = 0     over \partial\Omega
//   | 
//   | | -\Delta p = y-y_d   in \Omega
//   | |   p = 0     over \partial\Omega
//   |
//   | | u = Projection_ad( (-a/b) p|\omega )

verbosity = 0;

//
//  Problem data
//
real umin=0., umax=100.;
real a = 0.5, b = 1.-a, abc = a/b;
real cyd = 1.;                             // coeff. for desired state
func yd = cyd*(x*x-y*y);                   // desired state

//
//  Data for computations
//
real errFP = 1.e-5;                        // Tolerance for stop test
int  itermax = 500;                        // Maximum number of iterations

//
//  Borders
//
int LC = 12,  LR = 9;                       // Labels
// outer boundary
border C(t = 0, 2*pi){x=5*cos(t); y=5*sin(t); label = LC;}
// inner boundary
border Rec1(t = 0., 1.){ x=2.-t;  y=3.;       label = LR;}
border Rec2(t = 0., 1.){ x=1.;    y=3.-6.*t;  label = LR;}
border Rec3(t = 0., 1.){ x=1.+t;  y=-3.;      label = LR;}
border Rec4(t = 0., 1.){ x=2.;    y=-3.+6.*t; label = LR;}

plot (C(50)+Rec1(5)+Rec2(30)+Rec3(5)+Rec4(30), dim = 2, cmm = "Borders", wait = true);

//
// MESH
//
mesh Th = buildmesh(C(50)+Rec1(5)+Rec2(30)+Rec3(5)+Rec4(30));
plot(Th, cmm = "Mesh", wait = true);
        
int Lomega = Th(1.9,2.9).region;            // region number of \omega
func rom  = (region==Lomega);               // caract. function of \omega

//
//  F.E. space
//
fespace Vh(Th,P1);
Vh sol, phi;                            // unknown and test func.
Vh yst, pst;                            // state, adjoint state
Vh u = 0., ydh = yd;                    // control, desired state
Vh F;                                   // RHS of pde's; to be defined at each iteration
plot(ydh, value = true, fill = true, cmm = "Desired state",    wait = true);
plot(u,   value = true, fill = true, cmm = "Starting control", wait = true);

//
//  Problem
//
problem Poisson(sol,phi,init = 1) = 
          int2d(Th)(dx(sol)*dx(phi)+dy(sol)*dy(phi))
        - int2d(Th)(F*phi)
        + on(LC, sol = 0.)
;

// -------------------------------------
//  Projection function on U_{ad}
//    umin \leq u \leq umax
// -------------------------------------
//  puu = max(min(uu,umax),umin)
// -------------------------------------
func real[int] PUad(real[int] & uu)
{
  int nu = uu.n;
  real[int] puu(nu);
  for (int i = 0; i < nu; i++){
    puu[i] = min( max(uu[i], umin), umax);
  }
  return puu;
}
// -------------------------------------

//
//  Computations
//
real error, erroru;
real[int] Uerror(1000);
Vh uold = u, udif;
Vh d, w;                          // gradient and its associated state
real A, B;                        // auxiliary variables
real rho;                         // optimal step

int k;
for (k = 0; k < itermax; k++)
{
  F = u;
  Poisson;
  yst = sol;
  plot(yst, value=true, fill=true, cmm="State, iteration "+ k, wait=true);
  
  F = yst - ydh;
  Poisson;
  pst = sol;
  plot(pst, value=true, fill=true, cmm="Adjoint state, iteration "+ k, wait=true);
  
  d   = a*rom*pst + b*u;
  
  //  Optimal step compùtation
  F = d;
  Poisson;
  w = sol;
  B = int2d(Th)(d*d);
  A = int2d(Th)(w*w);
  A = a*A + b*B;
  rho = B/A;
  u = u - rho*d;
  u[] = PUad(u[]);
  u   = rom*u;
  plot(u, value=true, fill=true, cmm="Control, iteration "+ k, wait=true);

  udif = u-uold;
  error = udif[].l2;                         // Other norms: l1, linfty
  erroru = u[].l2;
  error = error/erroru;
  Uerror[k] = error;
  cout << " Iteration = " << k << ", error = " << error << ", rho = " << rho << endl;
  
  if (error <= errFP) break;
  
  uold = u;
} 

if (k == itermax) {
  cout << " Attention:: too many iterates " << endl;
  cout << " Stop test not attended " << endl;
  k = k-1;
}

//
//  plot of error
//
real[int] absx(0:k);
plot( [absx, Uerror(0:k)], cmm = "End of iterates. Evolution of the error. Final error= "+error, wait = true);

Ejercicio 1

Modificar este programa para utilizar el método del gradiente con paso fijo.

Ejercicio 2

Escribir un programa para resolver, mediante el algoritmo del gradiente con pado óptimo, el siguiente problema:

\begin{matrix} (14) & {\begin{cases} Minimizar J (u) = \frac{a}{2} \int_{Ω} | y - y_{d} |^{2} d x + \frac{b}{2} \int_{ω} | u |^{2} d x \\ sujeto a u \in U_{a d} \subset L^{2} (ω) \end{cases} \end{matrix}

${\cal U}_{ad}$ $\eqref{eq.uad}$ $y$ la solución de

\begin{matrix} (15) & {\begin{cases} - Δ y + A_{1} \partial_{1} y + A_{2} \partial_{2} y + B y = u 1_{ω} en Ω \\ y = 0 sobre \partial Ω \end{cases} \end{matrix}

$A_1, \, A_2, \, B \,\colon\, \Omega \to \mathbb{R}$ son dadas y se supone que poseen las propiedades necesarias para que la forma bilineal correspondiente a la formulación débil sea coerciva.

El sistema de optimalidad para este problema es el siguiente:

\begin{matrix} (16) & {\begin{cases} {\begin{cases} - Δ \hat{y} + A_{1} \partial_{1} \hat{y} + A_{2} \partial_{2} \hat{y} + B \hat{y} = \hat{u} 1_{ω} en Ω \\ \hat{y} = 0 sobre \partial Ω \end{cases} \\ {\begin{cases} - Δ \hat{p} - \partial_{1} (A_{1} \hat{p}) - \partial_{2} (A_{2} \hat{p}) + B \hat{p} = \hat{y} - y_{d} en Ω \\ \hat{p} = 0 sobre \partial Ω \end{cases} \\ \hat{u} = P_{a d} (- \frac{a}{b} \hat{p} |_{ω}) \end{cases} \end{matrix}

Ejercicio 3

Escribir un programa para resolver, mediante el algoritmo del gradiente con paso fijo, el siguiente problema:

\begin{matrix} (17) & {\begin{cases} Minimizar J (u) = \frac{a}{2} \int_{Ω} | y - y_{d} |^{2} d x + \frac{b}{2} \int_{ω} | u |^{2} d x \\ sujeto a u \in U_{a d} \subset L^{2} (ω) \end{cases} \end{matrix}

${\cal U}_{ad}$ $\eqref{eq.uad}$ $(y, z)$ es la solución de

\begin{matrix} (18) & {\begin{cases} - Δ y + A y + B z = u 1_{ω} en Ω \\ - Δ z + C y + D z = 0 en Ω \\ y = 0, \partial_{n} z = 0 sobre \partial Ω \end{cases} \end{matrix}

$A, \, B, \, C, \, D \colon\, \Omega \to \mathbb{R}$ son dadas y se supone que poseen las propiedades necesarias para que la forma bilineal de la correspondiente a la formulación débil sea coerciva.

El sistema de optimalidad para este problema es ahora:

\begin{matrix} (19) & {\begin{cases} {\begin{cases} - Δ \hat{y} + A \hat{y} + B \hat{z} = \hat{u} 1_{ω} en Ω \\ - Δ \hat{z} + C \hat{y} + D \hat{z} = 0 en Ω \\ \hat{y} = 0, \partial_{n} \hat{z} = 0 sobre \partial Ω \end{cases} \\ {\begin{cases} - Δ \hat{p} + A \hat{p} + C \hat{q} = \hat{y} - y_{d} en Ω \\ - Δ \hat{q} + B \hat{p} + D \hat{q} = 0 en Ω \\ \hat{p} = 0, \partial_{n} \hat{q} = 0 sobre \partial Ω \end{cases} \\ \hat{u} = P_{a d} (- \frac{a}{b} \hat{p} |_{ω}) \end{cases} \end{matrix}

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla