# Optimal Control#

Thanks to the function BFGS it is possible to solve complex nonlinear optimization problem within FreeFem++. For example consider the following inverse problem

\begin{eqnarray} \min_{b, c, d\in R}J &=& \int_E(u-u_d)^2\\ -\nabla(\kappa(b, c, d)\cdot\nabla u) &=& 0\\ u|_\Gamma &=& u_\Gamma \end{eqnarray}

where the desired state $u_d$, the boundary data $u_\Gamma$ and the observation set $E\subset\Omega$ are all given. Furthermore let us assume that

\kappa(x)=1+bI_B(x)+cI_C(x)+dI_D(x)~~~\forall x\in\Omega

where $B,C,D$ are separated subsets of $\Omega$.

To solve this problem by the quasi-Newton BFGS method we need the derivatives of $J$ with respect to $b,c,d$. We self explanatory notations, if $\delta b,\delta c,\delta d$ are variations of $b,c,d$ we have

\begin{eqnarray} \delta J &\approx& 2\int_E(u-u_d)\delta u\\ -\nabla(\kappa\cdot\nabla\delta u) &\approx& \nabla(\delta\kappa\cdot\nabla u)\\ \delta u|_\Gamma &=& 0 \end{eqnarray}

Obviously $J'_b$ is equal to $\delta J$ when $\delta b=1,\delta c=0,\delta d=0$, and so on for $J'_c$ and $J'_d$.

All this is implemented in the following program :

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 // Mesh border aa(t=0, 2*pi){x=5*cos(t); y=5*sin(t);}; border bb(t=0, 2*pi){x=cos(t); y=sin(t);}; border cc(t=0, 2*pi){x=-3+cos(t); y=sin(t);}; border dd(t=0, 2*pi){x=cos(t); y =-3+sin(t);}; mesh th = buildmesh(aa(70) + bb(35) + cc(35) + dd(35)); // Fespace fespace Vh(th, P1); Vh Ib=((x^2+y^2)<1.0001), Ic=(((x+3)^2+ y^2)<1.0001), Id=((x^2+(y+3)^2)<1.0001), Ie=(((x-1)^2+ y^2)<=4), ud, u, uh, du; // Problem real[int] z(3); problem A(u, uh) = int2d(th)( (1+z[0]*Ib+z[1]*Ic+z[2]*Id)*(dx(u)*dx(uh) + dy(u)*dy(uh)) ) + on(aa, u=x^3-y^3) ; // Solve z[0]=2; z[1]=3; z[2]=4; A; ud = u; ofstream f("J.txt"); func real J(real[int] & Z){ for (int i = 0; i < z.n; i++) z[i] =Z[i]; A; real s = int2d(th)(Ie*(u-ud)^2); f << s << " "; return s; } // Problem BFGS real[int] dz(3), dJdz(3); problem B (du, uh) = int2d(th)( (1+z[0]*Ib+z[1]*Ic+z[2]*Id)*(dx(du)*dx(uh) + dy(du)*dy(uh)) ) + int2d(th)( (dz[0]*Ib+dz[1]*Ic+dz[2]*Id)*(dx(u)*dx(uh) + dy(u)*dy(uh)) ) +on(aa, du=0) ; func real[int] DJ(real[int] &Z){ for(int i = 0; i < z.n; i++){ for(int j = 0; j < dz.n; j++) dz[j] = 0; dz[i] = 1; B; dJdz[i] = 2*int2d(th)(Ie*(u-ud)*du); } return dJdz; } real[int] Z(3); for(int j = 0; j < z.n; j++) Z[j]=1; BFGS(J, DJ, Z, eps=1.e-6, nbiter=15, nbiterline=20); cout << "BFGS: J(z) = " << J(Z) << endl; for(int j = 0; j < z.n; j++) cout << z[j] << endl; // Plot plot(ud, value=1, ps="u.eps"); 

In this example the sets $B,C,D,E$ are circles of boundaries $bb,cc,dd,ee$ and the domain $\Omega$ is the circle of boundary $aa$.

The desired state $u_d$ is the solution of the PDE for $b=2,c=3,d=4$. The unknowns are packed into array $z$.

Note

It is necessary to recopy $Z$ into $z$ because one is a local variable while the other one is global.

The program found $b=2.00125,c=3.00109,d=4.00551$.

figure 1 shows $u$ at convergence and the successive function evaluations of $J$.

Fig. 1: On top the level lines of $u$. At the bottom the successive evaluations of $J$ by BFGS (5 values above 500 have been removed for readability)

Note that an adjoint state could have been used. Define $p$ by

\begin{eqnarray} -\nabla\cdot(\kappa\nabla p) &=& 2I_E(u-u_d)\\ p|_\Gamma &=& 0 \end{eqnarray}

Consequently

\begin{eqnarray} \delta J = -\int_{\Omega}(\nabla\cdot(\kappa\nabla p))\delta u\nonumber\\ = \int_\Omega(\kappa\nabla p\cdot\nabla\delta u) =-\int_\Omega(\delta\kappa\nabla p\cdot\nabla u) \end{eqnarray}

Then the derivatives are found by setting $\delta b=1, \delta c=\delta d=0$ and so on:

\begin{eqnarray} J'_b=-\int_B \nabla p\cdot\nabla u\\ J'_c=-\int_C \nabla p\cdot\nabla u\\ J'_d=-\int_D \nabla p\cdot\nabla u \end{eqnarray}

Note

As BFGS stores an $M\times M$ matrix where $M$ is the number of unknowns, it is dangerously expensive to use this method when the unknown $x$ is a Finite Element Function. One should use another optimizer such as the NonLinear Conjugate Gradient NLCG (also a key word of FreeFem++).