Optimal Control of the Poisson Equation with OpenFOAM


PDF version | Download Openfoam Code… | DOI

In this tutorial, we show how to use the C++ library OpenFOAM (Open Field Operation and Manipulation) in order to solve control problems for partial differential equations (PDE).

OpenFOAM is a free open source toolbox that allows the user to code customized solvers for Continuum Mechanics, with special attention to the field of Computational Fluid Dynamics. One can either use one of the multiple solvers already programmed for steady-state, transient or multiphase problems in fluid mechanics, or code his own solver depending on particular needs.

In this work, we will solve the optimal control problem

(1)   \begin{equation*}  \min _{u \in L^2 \left( \Omega \right)} \mathcal{J}\left( u\right) = \min _{u \in L^2 \left( \Omega \right)} \frac{1}{2} \int_{\Omega} \left( y - y_d \right) ^2 \mathrm{d} \Omega + \frac{\beta}{2} \int_{\Omega} u ^2 \mathrm{d} \Omega, \end{equation*}

where u is the control variable, y the state variable and y_d a target function. The minimization problem is subject to the elliptic partial differential equation

(2)   \begin{align*} -\Delta y &  = f + u && \mathrm{in} \enspace \Omega , \\ \end{align*}

(3)   \begin{align*} y & = 0 && \mathrm{on} \enspace \Gamma . \end{align*}

We use some classical gradient descent methods based on the adjoint methodology, namely the steepest descent and conjugate gradient methods. The corresponding adjoint system for 2 writes as,

(4)   \begin{align*} - \Delta \lambda & = y - y_d && \mathrm{in} \enspace \Omega,  \\ \end{align*}

(5)   \begin{align*} \lambda & = 0 \ && \mathrm{on} \enspace \Gamma.  \end{align*}

In the steepest descent gradient method, we compute the directional derivative of the cost function in 1,

(6)   \begin{equation*} J^\prime \left( u \right) \delta u = \int_{\Omega} \left( \lambda + \beta u \right) \delta u \mathrm{d} \Omega, \end{equation*}

where \lambda is solution to the elliptic problem in 4, and use the gradient to update the control variable as

(7)   \begin{equation*} u^{\left( n + 1 \right)} = u^{\left( n \right)} - \gamma \left( \lambda^{\left( n \right)} + \beta u^{\left( n \right)} \right), \end{equation*}

for some value of \gamma.

One of the main advantages of OpenFOAM is its friendly syntax to describe PDE, as can be seen in the lines 10 and 13 in the code below.

// Initialize L2 norm of control update
scalar pL2 = 10*tol;
 
// Get the volume of mesh cells
scalarField volField = mesh.V();
 
while (runTime.loop() && ::sqrt(pL2) > tol)
{
    // Solve primal equation
    solve(fvm::laplacian(k, y) + u);
 
    // Solve ajoint equation
    solve(fvm::laplacian(k, lambda) + y - yd);
 
    // Update control
    p = lambda + beta*u;
    u = u - gamma*p;
 
    // Control update norm
    pL2 = ::sqrt(gSum(volField*p.internalField()*p.internalField()));
 
    // Compute cost function
    J = 0.5*gSum(volField*(Foam::pow(y.internalField()-yd.internalField(),2) \
	    + beta*Foam::pow(u.internalField(),2)));
 
    // Display information
    Info << "Iteration " << runTime.timeName() << " - " \
	    << "Cost value " << J << " - " \
	    << "Control variation L2 norm " << ::sqrt( pL2 ) << endl;
 
    // Save current iteration results
    runTime.write();
}

In order to improve the convergence to the optimal control, the conjugate gradient method is applied also to the problem. First, the state variable must be separated in two terms as

(8)   \begin{equation*} y = y_u + y_{f}, \end{equation*}

where y_u solves the state equation with zero Dirichlet boundary conditions,

(9)   \begin{align*} -\Delta y_u &  = u && \mathrm{in} \enspace \Omega , \\ \end{align*}

(10)   \begin{align*} y_u & = 0 && \mathrm{on} \enspace \Gamma , \end{align*}

and y_{f} is the control-free solution to the state equation,

(11)   \begin{align*} -\Delta y_{f} &  = f && \mathrm{in} \enspace \Omega , \\ \end{align*}

(12)   \begin{align*} y_{f} & = 0 && \mathrm{on} \enspace \Gamma . \end{align*}

Now, using the above separation of the state variable, one part depending on the control and the other part on a known source term, the cost functional in 1 can be expressed as

(13)   \begin{equation*}  J \left( u \right) = \frac{1}{2} \left( y_u + y_f - y_d, y_u + y_f - y_d \right)_{L^2\left( \Omega \right)} + \frac{\beta}{2} \left( u , u \right) _{L^2\left( \Omega \right)}. \end{equation*}

We define a linear operator \Lambda: L^2\left( \Omega \right) \rightarrow L^2\left( \Omega \right) that takes a control u and returns the solution to problem 910, so that y_u = \Lambda u. We introduce as well its adjoint operator \Lambda^*: L^2\left( \Omega \right) \rightarrow L^2\left( \Omega \right) that takes the source term in 4 and solves the Poisson equation, so that \lambda = \Lambda^* \left( y - y_d \right).

The directional derivative of the functional 13 then reads as

(14)   \begin{equation*} J^\prime\left( u \right) \delta u = \left( \underbrace{ \left(  \Lambda^* \Lambda + \beta I  \right)}_{A_{cg}} u - \underbrace{ \Lambda^* \left( y_d - y_f \right)}_{b_{cg}}, \delta u \right) _{L^2\left( \Omega \right)}. \end{equation*}

After having identified A_{cg} and b_{cg} we can use the conjugate gradient method to reach the optimal control faster.

Algorithm 1 Optimal control with Conjugate Gradient Method

Require: y_D, u^{\left(0\right)}, \beta, y_d, tol

  1. n \gets 0
  2. compute the control-free solution, y_f
  3. b \gets \Lambda^*\left( y_d - y_f \right)
  4. z \gets \Lambda u
  5. g \gets \Lambda^*z + \beta u - b
  6. h \gets ||g||^2_{L^2\left(\left[0,T\right]\right)}
  7. h_a \gets h
  8. r \gets -g
  9. While ||r||_{L^2\left(\left[0,T\right]\right)} > tol do
  10.     z \gets \Lambda r
  11.     w \gets \Lambda^*z + \beta r
  12.     \alpha \gets \frac{h}{\left(r,w\right)_{L^2\left(\left[0,T\right]\right)}}
  13.     u \gets u + \alpha r
  14.     g \gets g + \alpha w
  15.     h_a \gets h
  16.     h \gets ||g||^2_{L^2\left(\left[0,T\right]\right)}
  17.     \gamma \gets \frac{h}{h_a}
  18.     r \gets -g + \gamma r
  19.     n \gets n + 1
// Compute b: right-hand side of Ax = b
solve(fvm::laplacian(k, y));
solve(fvm::laplacian(k, lambda) + yd - y);
volScalarField b = lambda;
 
// Compute A*f0
// Primal equation
solve(fvm::laplacian(k, y0) + u);
// Adjoint equation
solve(fvm::laplacian(k, lambda) + y0);
 
// Compute g0: initial gradient
volScalarField g = lambda + beta*u - b;
scalar gL2 = gSum( volField * g.internalField() * g.internalField() );
scalar gL2a = 0;
 
// Compute initial residual and its norm
volScalarField r = -g;
scalar rL2 = gL2;
 
volScalarField w = g;
 
scalar alpha = 0;
scalar gamma = 0;
 
while (runTime.loop() && ::sqrt( rL2 ) > tol)
{
    r.dimensions().reset( u.dimensions() );
 
    // Compute w = A*r
    solve(fvm::laplacian(k, y0) + r);
    solve(fvm::laplacian(k, lambda) + y0);
    w = lambda + beta*r;
 
    // Update alpha
    alpha = gL2 / gSum( volField * r.internalField() * w.internalField() );
 
    // Update control
    u = u + alpha*r;
 
    // Update cost gradient and its norm
    g = g + alpha*w;
    gL2a = gL2;
    gL2 = gSum( volField * g.internalField() * g.internalField() );
    gamma = gL2/gL2a;
 
    // Update residual and its norm
    r.dimensions().reset( g.dimensions() );
    r = -g + gamma*r;
    rL2 = gSum( volField * r.internalField() * r.internalField() );
 
    solve(fvm::laplacian(k, y) + u);
 
    J = 0.5*gSum(volField*(Foam::pow(y.internalField()-yd.internalField(),2) \
     + beta*Foam::pow(u.internalField(),2)));
 
    Info << "Iteration no. " << runTime.timeName() << " - " \
	    << "Cost value " << J << " - " \
	    << "Residual " << ::sqrt( rL2 ) << endl;
 
    runTime.write();
}
 

Both solvers, laplaceAdjointFoam for the steepest descent method with \gamma = 10 and laplaceCGAdjointFoam for the conjugate gradient method, have been tested in a square domain [0, 1] \times [0, 1] with zero Dirichlet boundary conditions and \beta = 10^{-3},10^{-4},10^{-5},10^{-6}. The target function is y_d = xy \sin \left( \pi x \right) \sin \left( \pi y \right).

The problem setup is called case in OpenFOAM and it is independent from the solver itself. A typical case in OpenFOAM has three folders: 0, where the problem fields are stored along with their boundary conditions; constant, that includes the mesh data and the physical properties; and system, with parameters regarding the numerical solution and problem output such as tolerances, linear solvers or write interval. In this example, the case folder is called laplaceAdjointFoamCase.

Figure 1.a: Functional value against the number of iterations of the gradient method. Steepest descent method. Steepest descent method
Figure 1.a: Functional value against the number of iterations of the gradient method. Steepest descent method. Steepest descent method
Figure 1.b: Functional value against the number of iterations of the gradient method. Steepest descent method. Conjugate gradient method.
Figure 1.b: Functional value against the number of iterations of the gradient method. Steepest descent method. Conjugate gradient method.
Figure 2.a
Figure 2.a: || y - y_d ||_{L^2 \left( \Omega \right)} against the number of iterations of the gradient method. Steepest descent method.
Figure 2.b
Figure 2.b: || y - y_d ||_{L^2 \left( \Omega \right)} against the number of iterations of the gradient method. Conjugate gradient method.
Bibliography

[1] The OpenFOAM Foundation, http://openfoam.org.
[2] F. Tröltzsch. Optimal control of partial differential equations: theory, methods, and applications. American Mathematical Soc., 2010.

Click here to review “Steepest Descent Method in OpenFOAM” on Github

Click here to review “Conjugate Gradient Method in OpenFOAM” on Github

Authors: Víctor Hernández-Santamaría, José Vicente Lorenzo, Enrique Zuazua
May, 2018