FreeFEM Documentation on GitHub

stars - forks

# 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{split}\begin{array}{rcl} \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{array}\end{split}$

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)\quad\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{split}\begin{array}{rcl} \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{array}\end{split}$

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

Fig. 40 and Fig. 41 show $$u$$ at convergence and the successive function evaluations of $$J$$.

Fig. 40 Level line of $$u$$.

Fig. 41 Successive evaluations of $$J$$ by BFGS (5 values above 500 have been removed for readability)

Optimal control

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

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

Consequently:

$\begin{split}\begin{array}{rcl} \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{array}\end{split}$

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

$\begin{split}\begin{array}{rcl} 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{array}\end{split}$

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).

Table of content