Control óptimo de una ecuación de Poisson mediante un algoritmo de tipo punto fijo

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 al control $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 $\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}

$z\in L^2(\omega)$ , entonces

\begin{matrix} (5) & \begin{matrix} P_{a d} (z) (x) = {\begin{cases} \overset{―}{u} & si z (x) > \overset{―}{u} \\ z (x) & si \underset{―}{u} \leq z (x) \leq \overset{―}{u} \\ \underset{―}{u} & si z (x) < \underset{―}{u} \end{cases} para x c.p.d. en ω \end{matrix} \end{matrix}

Algoritmo de resolución

Teniendo en cuenta la caracterización del control óptimo, se puede considerar que éste es un punto fijo de la aplicación

\begin{matrix} u \in U_{a d} \to y (estado) \to p (estado adjunto) \to u^{'} \in U_{a d} \end{matrix}

$\eqref{eq.minimizar}$ :

$u^0\in{\cal U}_{ad}$ $\varepsilon > 0$
$n\geq 0$ $u^n \in {\cal U}_{ad}$ $y^n$ solución de
$\begin{matrix} (6) & {\begin{cases} - Δ y^{n} = u^{n} 1_{ω} en Ω \\ y^{n} = 0 sobre \partial Ω \end{cases} \end{matrix}$
$p^n$ solución de
$\begin{matrix} (7) & {\begin{cases} - Δ p^{n} = y^{n} - y_{d} en Ω \\ p = 0 sobre \partial Ω \end{cases} \end{matrix}$
(2.3) Tomar
$\begin{matrix} (8) & u^{n + 1} = P_{a d} (- \frac{a}{b} p^{n} |_{ω}) \end{matrix}$
(2.4) Si
$\begin{matrix} (9) & ∥ u^{n + 1} - u^{n} ∥_{L^{2} (ω)} \leq ε (o bien \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 2.

Código

$\eqref{eq.poisson.yn}$ $\eqref{eq.poisson.pn}$ conducen a la misma forma bilineal. Por lo tanto, bastará con definir un único problema e ir cambiando adecuadamente los segundos miembros correspondientes a cada caso.


//   ControlCond_FP.edp
//
//   OPTIMAL CONTROL OF A CAPACITOR
//   FIXED-POINT ALGORITHM
//
//   |  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.2, 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 = 100;                        // 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;

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);
  
  u   = -abc*pst;
  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 << 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

Escribir un programa para resolver el siguiente problema:

\begin{matrix} (10) & {\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} (11) & {\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}$ y se supone que poseen las propiedades necesarias para que la correspondiente forma bilineal de la formulación débil sea coerciva.

El sistema de optimalidad para este problema es:

\begin{matrix} (12) & {\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 2

Escribir un programa para resolver el siguiente problema:

\begin{matrix} (13) & {\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)$ la solución de

\begin{matrix} (14) & {\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}$ y se supone que poseen las propiedades necesarias para que la correspondiente forma bilineal de la formulación débil sea coerciva.

sistema de optimalidad $(\hat u, (\hat y,\hat z), (\hat p, \hat q))$ verifican

\begin{matrix} (15) & {\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