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 $$ \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 $$ 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 $$ \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 $$ 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

 // file optimcontrol.edp
 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 Vh(th,P1);
 Vh Ib=((x^2+y^2)<1.0001),
    Ic=(((x+3)^2+ y^2)<1.0001),
    Ie=(((x-1)^2+ y^2)<=4),
 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);
 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;

 real[int] dz(3), dJdz(3);

 problem B(du,uh)

 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;
  cout << "BFGS: J(z) = " << J(Z) <<  endl;
  for(int j=0;j<z.n;j++) cout<<z[j]<<endl;

In this example the sets \(B,C,D,E\) are circles of boundaries \(bb,cc,dd,ee\) are 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\). Notice that 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 [figcontrol] shows \(u\) at convergence and the successive function evaluations of \(J\). [htbp]

[figcontrol]On the left the level lines of \(u\). On the right 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 $$ -\nabla\cdot(\kappa\nabla p)=2I_E(u-u_d),   p|_\Gamma=0 $$ Consequently

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

Then the derivatives are found by setting \(\delta b=1, \delta c=\delta d=0\) and so on: $$ 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 $$ Remark 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++). See the file algo.edp in the examples folder.