Optimal Control of the Poisson equation

Optimal Control of the Poisson equation#

Section author: Jørgen S. Dokken (dokken@simula.no).

Original implementation in dolfin-adjoint was by Simon W. Funke.

This demo solves the mother problem of PDE-constrained optimization: the optimal control of the Possion equation. Physically, this problem can be interpreted as finding the best heating/cooling of a cooktop to achieve a desired temperature profile.

This example introduces the basics of how to solve optimization problems with DOLFINx-adjoint.

Problem definition#

Mathematically, the goal is to minimize the following tracking type functional:

\[ \min_{f \in Q} J(u) = \frac{1}{2} \int_{\Omega} (u - d)^2 ~\mathrm{d}x + \frac{\alpha}{2}\int_{\Omega} f^2~\mathrm{d} x \]

subject to the Poisson equation with Dirichlet boundary conditions:

\[\begin{split} \begin{align} - \kappa \Delta u &= f && \text{in } \Omega, \\ u &= 0 && \text{on } \partial\Omega, \\ a &\leq f \leq b && \end{align} \end{split}\]

where \(\Omega\) is the domain of interest, \(u: \Omega \mapsto \mathbb{R}\) is the unknown temperature, \(\kappa\in\mathbb{R}\) is the thermal diffusivity, \(f: \Omega \mapsto \mathbb{R}\) is the unknown control function acting as a source term, \(d:\Omega\mapsto \mathbb{R}\) is the desired temperature profile, and \(\alpha\in[0,\infty)\) is a Tikhonov regularization parameter, and \(a,b\in\mathbb{R}\) are the lower and upper bounds on the control function \(f\). Note that \(f(x)>0\) corresponds to heating, while \(f(x)<0\) corresponds to cooling.

It can be shown that this problem is well-posed and has a unique solution, see for instance Section 1.5 [Ulb09] or [Troltzsch10].

Implementation#

We start by import the necessary modules for this demo, which includes mpi4py, dolfinx and dolfinx_adjoint.

import os
import sys

from mpi4py import MPI

import dolfinx

Next we import Moola, which is a Python package containing a collection of optimization solvers specifically designed for PDE-constrained optimization problems.

import moola
import numpy as np
import pandas
import pyadjoint
import pyvista
import ufl
from moola.adaptors import DolfinxPrimalVector  # noqa: E402

import dolfinx_adjoint

We configure Pyvista for rendering

Hide code cell source 
pyvista.set_jupyter_backend("html")
if sys.platform == "linux" and (os.getenv("CI") or pyvista.OFF_SCREEN):
    pyvista.start_xvfb(0.05)

Next we create a regular mesh of a unit square, which will be our domain \(\Omega\). Some optimization algorithms suffer from bad performance when the mesh is non-uniform (i.e. when the mesh is partially refined). To demonstrate that Moola does not have this issue, we will refine our mesh in the center of the domain. We use the DOLFINx refine function to do this.

n = 16
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n)


def refinement_region(x, tol=1e-14):
    """Define a conditional to determine the refinement region.

    Args:
        x: Coordinates of mesh nodes, shape (num_nodes, 3)
        tol: Tolerance for being within the region.
    """
    clause_x = abs(x[0] - 0.5) < 0.25 + tol
    clause_y = abs(x[1] - 0.5) < 0.25 + tol
    return clause_x & clause_y


mesh.topology.create_connectivity(1, mesh.topology.dim)
edges_to_refine = dolfinx.mesh.locate_entities(mesh, 1, refinement_region)
refined_mesh_data = dolfinx.mesh.refine(mesh, edges_to_refine)
refined_mesh = refined_mesh_data[0]
tdim = refined_mesh.topology.dim
del mesh

Next we use Pyvista to plot the mesh

grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(refined_mesh))
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, color="lightgrey")
plotter.view_xy()
plotter.show()

Then we define the discrete function spaces \(V\) and \(Q\) for the state and control variable, respectively

V = dolfinx.fem.functionspace(refined_mesh, ("Lagrange", 1))  # type: ignore[arg-type]
Q = dolfinx.fem.functionspace(refined_mesh, ("Discontinuous Lagrange", 0))  # type: ignore[arg-type]

The optimization algorithm will use the value of the control function \(f\) as an initial guess for the optimization. A zero intial guess seems to be too simple: For example the L-BFGS algorithm will find the optimial control with just two iterations. To make it more interesting, we choose a non-zero initial guess.

f = dolfinx_adjoint.Function(Q, name="Control")
f.interpolate(lambda x: x[0] + x[1])  # Set intial guess

Note

As opposed to standard DOLFINx code, we use dolfinx_adjoint.Function to create the control function. This is so that we can track it throughout the program on the computational tape.

We also create a state variable that we will store the solution to the Poisson equation in.

uh = dolfinx_adjoint.Function(V, name="State")

Next, we define the variational formulation of the Poisson equation.

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
kappa = dolfinx.fem.Constant(refined_mesh, dolfinx.default_scalar_type(1.0))  # Thermal diffusivity
F = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx - f * v * ufl.dx
a, L = ufl.system(F)

We create the Dirichlet BC

refined_mesh.topology.create_connectivity(tdim - 1, tdim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(refined_mesh.topology)
exterior_dofs = dolfinx.fem.locate_dofs_topological(V, tdim - 1, exterior_facets)
zero = dolfinx.fem.Constant(refined_mesh, dolfinx.default_scalar_type(0.0))
bc = dolfinx.fem.dirichletbc(zero, exterior_dofs, V)

Next, we define a dolfinx_adjoint.LinearProblem instance, which overloads the dolfinx.fem.petsc.LinearProblem class.

Note

When creating the :py:func:dolfinx_adjoint.LinearProblem, we can specify the solver options that are passed on to the underlying PETSc Krylov subspace solver. This is also the place to pass in solver options for the first and second order adjoint equations and the tangent linear model (TLM) equation.

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
}
problem = dolfinx_adjoint.LinearProblem(
    a,
    L,
    u=uh,
    bcs=[bc],
    petsc_options=petsc_options,
    adjoint_petsc_options=petsc_options,
    tlm_petsc_options=petsc_options,  # type: ignore
)
problem.solve()

(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False, float64, []), (2,)), 1), Basix element (P, triangle, 1, gll_warped, unset, False, float64, [])), 1),
 4,
 1)

Note

Note that we can pass in solver options for the adjoint equation via the keyword argument adjoint_petsc_options. As we a solving a linear, symmetric problem, i,e, a self-adjoint problem, we use the same options for both the forward and adjoint problems.

Before we can start the optimization, we need to specity the control variable and the functional of interest. For this example we use \(d(x,y)=\frac{1}{2\pi^2} \sin(\pi x)\sin( \pi y)\) as the desired temperature profile

x, y = ufl.SpatialCoordinate(refined_mesh)
d = 1 / (2 * ufl.pi**2) * ufl.sin(ufl.pi * x) * ufl.sin(ufl.pi * y)

The functional is written out in ufl and assembled with dolfinx_adjoint.assemble_scalar

alpha = dolfinx.fem.Constant(refined_mesh, dolfinx.default_scalar_type(1.0e-6))  # Tikhonov regularization parameter
alpha.name = "alpha"  # type: ignore
J_symbolic = 0.5 * ufl.inner(uh - d, uh - d) * ufl.dx + 0.5 * alpha * ufl.inner(f, f) * ufl.dx
J = dolfinx_adjoint.assemble_scalar(J_symbolic)

The next step is to formulate the so-called reduced optimization problem. The idea is that the solution \(u\) cna be considered as a functin of \(f\): given a value of \(f\), we can solve the Poisson equation and obtain the associated solution \(uh\). Bu denoting this solution function as \(u(f)\) we can rewrite the original optimization problem as a reduced problem:

\[ \min_{f \in Q} \hat J(u(f), f) = \frac{1}{2} \int_{\Omega} (u(f) - d)^2 ~\mathrm{d}x + \frac{\alpha}{2}\int_{\Omega} f^2~\mathrm{d} x \]

Note that no PDE-constraint is required anymore, since it is implicitly contained in the solution function.

dolfinx-adjoint can automatically reduce the optimization problem by using pyadjoint.ReducedFunctional. This object solves the forward PDE using the pyadjoint.Tape each time the functional is evaluated. Additionally, it derives and solves the adjoint equation each time the function gradient should be evaluated.

control = pyadjoint.Control(f)
Jhat = pyadjoint.ReducedFunctional(J, control)

Now that all ingredients are in place, we can perform the optimization. For this, we employ the moola.MoolaOptimizationProblem to generate a problem that is compatible with the Moola framework.

optimization_problem = pyadjoint.MoolaOptimizationProblem(Jhat)
f_moola = DolfinxPrimalVector(f)

Then, we wrap the control function into a Moola object, and create a `moola.BFGS`` solver for solving the optimisation problem

optimization_options = {"jtol": 0, "gtol": 1e-9, "Hinit": "default", "maxiter": 100, "mem_lim": 10, "rjtol": 0}
solver = moola.BFGS(optimization_problem, f_moola, options=optimization_options)
solution = solver.solve()

L-BFGS method.
------------------------------
Line search:		 strong_wolfe
Maximum iterations:	 100

iteration = 0:	objective = 0.0001438611042603543:	grad_norm = 0.0007841745828175678:
iteration = 1:	objective = 1.9506963766488212e-05:	grad_norm = 0.00010356235912243852:	delta_J = 0.00012435414049386608:
iteration = 2:	objective = 1.5553208955156504e-05:	grad_norm = 8.861605017903993e-05:	delta_J = 3.953754811331708e-06:
iteration = 3:	objective = 3.010046175058102e-06:	grad_norm = 2.127212976822824e-05:	delta_J = 1.2543162780098401e-05:
iteration = 4:	objective = 1.8673662092088006e-06:	grad_norm = 1.6644826983888767e-05:	delta_J = 1.1426799658493015e-06:
iteration = 5:	objective = 5.923951821332709e-07:	grad_norm = 3.7066721705558705e-06:	delta_J = 1.2749710270755296e-06:

iteration = 6:	objective = 4.807412579873319e-07:	grad_norm = 3.0650713236135173e-06:	delta_J = 1.1165392414593898e-07:
iteration = 7:	objective = 2.716701154912723e-07:	grad_norm = 1.571287104275792e-06:	delta_J = 2.0907114249605962e-07:
iteration = 8:	objective = 2.2159945224985335e-07:	grad_norm = 1.7388814343553936e-06:	delta_J = 5.007066324141897e-08:
iteration = 9:	objective = 2.1602912928896518e-07:	grad_norm = 2.3357449403932515e-06:	delta_J = 5.570322960888171e-09:
iteration = 10:	objective = 1.9256190281008928e-07:	grad_norm = 6.633042295564967e-07:	delta_J = 2.3467226478875907e-08:
iteration = 11:	objective = 1.8784461159695448e-07:	grad_norm = 6.087498354746665e-07:	delta_J = 4.717291213134796e-09:

iteration = 12:	objective = 1.680561798243559e-07:	grad_norm = 8.355672513272449e-07:	delta_J = 1.978843177259857e-08:
iteration = 13:	objective = 1.5481165825312897e-07:	grad_norm = 6.872473582798471e-07:	delta_J = 1.3244521571226943e-08:
iteration = 14:	objective = 1.528615337803931e-07:	grad_norm = 2.864038115600805e-06:	delta_J = 1.950124472735878e-09:
iteration = 15:	objective = 1.468339890218053e-07:	grad_norm = 1.378267872870736e-06:	delta_J = 6.027544758587798e-09:
iteration = 16:	objective = 1.4335633126598514e-07:	grad_norm = 5.834782031192003e-07:	delta_J = 3.477657755820151e-09:
iteration = 17:	objective = 1.4132525789113624e-07:	grad_norm = 5.47136468016973e-07:	delta_J = 2.031073374848899e-09:

iteration = 18:	objective = 1.3960848396497156e-07:	grad_norm = 6.228882454658727e-07:	delta_J = 1.7167739261646828e-09:
iteration = 19:	objective = 1.380943122427587e-07:	grad_norm = 8.02205722480128e-07:	delta_J = 1.5141717222128564e-09:
iteration = 20:	objective = 1.364603681173245e-07:	grad_norm = 2.5449188043197694e-07:	delta_J = 1.633944125434206e-09:
iteration = 21:	objective = 1.3605446849105251e-07:	grad_norm = 2.1553063403861824e-07:	delta_J = 4.058996262719824e-10:
iteration = 22:	objective = 1.351337279205677e-07:	grad_norm = 3.163602607334195e-07:	delta_J = 9.207405704848115e-10:
iteration = 23:	objective = 1.3422908553490524e-07:	grad_norm = 3.1904647256772685e-07:	delta_J = 9.046423856624576e-10:

iteration = 24:	objective = 1.3219755582242469e-07:	grad_norm = 3.7506658183916456e-07:	delta_J = 2.0315297124805592e-09:
iteration = 25:	objective = 1.3213364308320176e-07:	grad_norm = 4.965521684046833e-07:	delta_J = 6.391273922292941e-11:
iteration = 26:	objective = 1.3155681101402125e-07:	grad_norm = 3.28171956977112e-07:	delta_J = 5.768320691805079e-10:
iteration = 27:	objective = 1.3139973933208658e-07:	grad_norm = 2.1893301562644314e-07:	delta_J = 1.570716819346706e-10:
iteration = 28:	objective = 1.3126226106474632e-07:	grad_norm = 9.218779469645262e-08:	delta_J = 1.3747826734025758e-10:
iteration = 29:	objective = 1.3120963760344938e-07:	grad_norm = 8.819698309416155e-08:	delta_J = 5.262346129693842e-11:

iteration = 30:	objective = 1.3098285560734137e-07:	grad_norm = 1.2100235731855799e-07:	delta_J = 2.2678199610801196e-10:
iteration = 31:	objective = 1.3090066538677157e-07:	grad_norm = 3.365722843898437e-07:	delta_J = 8.21902205698028e-11:
iteration = 32:	objective = 1.3076102666413506e-07:	grad_norm = 1.0341572170501043e-07:	delta_J = 1.3963872263650775e-10:
iteration = 33:	objective = 1.307291799823549e-07:	grad_norm = 7.01378292429097e-08:	delta_J = 3.1846681780161484e-11:
iteration = 34:	objective = 1.3070219143643606e-07:	grad_norm = 8.346469635797048e-08:	delta_J = 2.698854591883611e-11:
iteration = 35:	objective = 1.3063924344172124e-07:	grad_norm = 1.1102154613322705e-07:	delta_J = 6.294799471482073e-11:

iteration = 36:	objective = 1.305479278007062e-07:	grad_norm = 1.6117402682475145e-07:	delta_J = 9.131564101502802e-11:
iteration = 37:	objective = 1.305293453898813e-07:	grad_norm = 1.563985445535282e-07:	delta_J = 1.8582410824916552e-11:
iteration = 38:	objective = 1.3048738774860097e-07:	grad_norm = 6.199825449747471e-08:	delta_J = 4.195764128033038e-11:
iteration = 39:	objective = 1.30471962422973e-07:	grad_norm = 5.492851850807347e-08:	delta_J = 1.5425325627969854e-11:
iteration = 40:	objective = 1.3045415345909568e-07:	grad_norm = 9.7532148198221e-08:	delta_J = 1.7808963877316356e-11:
iteration = 41:	objective = 1.3042428205376165e-07:	grad_norm = 1.0226149327873368e-07:	delta_J = 2.98714053340275e-11:

iteration = 42:	objective = 1.3038041381139885e-07:	grad_norm = 7.674191296743462e-08:	delta_J = 4.386824236280233e-11:
iteration = 43:	objective = 1.303686414706444e-07:	grad_norm = 1.7662661080999413e-07:	delta_J = 1.177234075445317e-11:
iteration = 44:	objective = 1.3034820375359734e-07:	grad_norm = 3.0605291075622174e-08:	delta_J = 2.0437717047052647e-11:
iteration = 45:	objective = 1.3034620020005397e-07:	grad_norm = 2.102669965274069e-08:	delta_J = 2.00355354337631e-12:
iteration = 46:	objective = 1.303426717251546e-07:	grad_norm = 2.0118977954415333e-08:	delta_J = 3.5284748993579865e-12:
iteration = 47:	objective = 1.3033463323389997e-07:	grad_norm = 3.4766649545739576e-08:	delta_J = 8.038491254639429e-12:

iteration = 48:	objective = 1.303325586797203e-07:	grad_norm = 9.666734851078971e-08:	delta_J = 2.074554179676196e-12:
iteration = 49:	objective = 1.3032902402892783e-07:	grad_norm = 2.8841940342072736e-08:	delta_J = 3.534650792468103e-12:
iteration = 50:	objective = 1.3032767869286106e-07:	grad_norm = 1.949724095199671e-08:	delta_J = 1.3453360667625192e-12:
iteration = 51:	objective = 1.303253152294287e-07:	grad_norm = 2.2470497933874727e-08:	delta_J = 2.36346343236743e-12:
iteration = 52:	objective = 1.303217438512714e-07:	grad_norm = 2.8691986403394645e-08:	delta_J = 3.5713781572870244e-12:
iteration = 53:	objective = 1.3031928211397128e-07:	grad_norm = 3.209497015469935e-08:	delta_J = 2.4617373001270514e-12:

iteration = 54:	objective = 1.303159133582194e-07:	grad_norm = 2.9859746257374074e-08:	delta_J = 3.3687557518741888e-12:
iteration = 55:	objective = 1.3031462029416804e-07:	grad_norm = 4.466505852309923e-08:	delta_J = 1.2930640513661203e-12:
iteration = 56:	objective = 1.3031287766765783e-07:	grad_norm = 2.3139872393237434e-08:	delta_J = 1.7426265102151242e-12:
iteration = 57:	objective = 1.3031143296861154e-07:	grad_norm = 1.7818542135281655e-08:	delta_J = 1.4446990462828769e-12:
iteration = 58:	objective = 1.3031084975971408e-07:	grad_norm = 1.8850707469602246e-08:	delta_J = 5.8320889745969e-13:
iteration = 59:	objective = 1.303091606552831e-07:	grad_norm = 1.2890759218292052e-08:	delta_J = 1.6891044309955524e-12:

iteration = 60:	objective = 1.3030864881490178e-07:	grad_norm = 3.537388275551271e-08:	delta_J = 5.118403813041634e-13:
iteration = 61:	objective = 1.3030769981982707e-07:	grad_norm = 1.610217865112296e-08:	delta_J = 9.489950747165902e-13:
iteration = 62:	objective = 1.3030711832916908e-07:	grad_norm = 1.4722151474578562e-08:	delta_J = 5.814906579860206e-13:
iteration = 63:	objective = 1.3030669143521868e-07:	grad_norm = 1.2430067558316044e-08:	delta_J = 4.268939504034354e-13:

Tolerance reached: grad_norm < gtol.

iteration = 64:	objective = 1.3030632301298418e-07:	grad_norm = 6.77000109351731e-09:	delta_J = 3.684222345006902e-13:

Next, we update the control function with the optimal value found by Moola and solve the forward problem to get the optimal state variable.

f_opt = solution["control"].data
f.x.array[:] = f_opt.x.array  # f_opt.x.array.copy()
problem.solve(annotate=False)

(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False, float64, []), (2,)), 1), Basix element (P, triangle, 1, gll_warped, unset, False, float64, [])), 1),
 4,
 1)

Error analysis#

For our desired temperature profile, there exist an analytic solution on the following form

\[\begin{split} \begin{align*} f_{analytic}(x,y) &= \frac{1}{1 + \alpha 4\pi^4} \sin(\pi x) \sin(\pi y)\\ u_{analytic}(x,y) &= \frac{1}{2\pi^2} f_{analytic}(x,y) \end{align*} \end{split}\]
# We use these, and compute the
f_analytic = 1 / (1 + alpha * 4 * pow(ufl.pi, 4)) * ufl.sin(ufl.pi * x) * ufl.sin(ufl.pi * y)
u_analytic = 1 / (2 * ufl.pi**2) * f_analytic

err_u = dolfinx_adjoint.error_norm(u_analytic, uh, norm_type="L2", annotate=False)
err_f = dolfinx_adjoint.error_norm(f_analytic, f, norm_type="L2", annotate=False)
print(f"Error in state variable: {err_u:.3e}")
print(f"Error in control variable: {err_f:.3e}")

Error in state variable: 7.602e-05
Error in control variable: 3.683e-02

We visualize the results using Pyvista.

Hide code cell source 
# Plotting space for source
plotting_space = dolfinx.fem.functionspace(refined_mesh, ("Discontinuous Lagrange", 1))  # type: ignore[arg-type]
u_plot = dolfinx.fem.Function(plotting_space)
u_plot.interpolate(uh)
f_plot = dolfinx.fem.Function(plotting_space)
f_plot.interpolate(f)

# Interpolate the analytic solutions to the plotting space
u_expr = dolfinx.fem.Expression(u_analytic, plotting_space.element.interpolation_points)
u_ex = dolfinx.fem.Function(plotting_space)
u_ex.interpolate(u_expr)
f_expr = dolfinx.fem.Expression(f_analytic, plotting_space.element.interpolation_points)
f_ex = dolfinx.fem.Function(plotting_space)
f_ex.interpolate(f_expr)

# Attach data to the plotting grid
plotting_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(plotting_space))  # type: ignore[arg-type]
plotting_grid.point_data["u_optimal"] = u_plot.x.array
plotting_grid.point_data["u_exact"] = u_ex.x.array
plotting_grid.point_data["f_optimal"] = f_plot.x.array
plotting_grid.point_data["f_exact"] = f_ex.x.array

# Create plotter for all variables
plotter = pyvista.Plotter(shape=(2, 2))
plotter.subplot(0, 0)
scale_factor = 25
warped_u = plotting_grid.warp_by_scalar("u_optimal", factor=scale_factor)
plotter.add_mesh(warped_u, color="lightgrey", scalars="u_optimal")
plotter.subplot(0, 1)
warped_u_ex = plotting_grid.warp_by_scalar("u_exact", factor=scale_factor)
plotter.add_mesh(warped_u_ex, color="lightgrey", scalars="u_exact")
plotter.subplot(1, 0)
scale_factor = 1
warped_f = plotting_grid.warp_by_scalar("f_optimal", factor=scale_factor)
plotter.add_mesh(warped_f, scalars="f_optimal", show_edges=True)
plotter.subplot(1, 1)
warped_f_ex = plotting_grid.warp_by_scalar("f_exact", factor=scale_factor)
plotter.add_mesh(warped_f_ex, scalars="f_exact", show_edges=True)
plotter.link_views((0, 1))
plotter.link_views((2, 3))
plotter.show()

Convergence analysis (mesh independence)#

It is highly desirable that the optimisation algorithm achieve mesh independence: i.e., that the required number of optimisation iterations is independent of the mesh resolution. Achieving mesh independence requires paying careful attention to the inner product structure of the function space in which the solution is sought.

We therefore perform a mesh convergence analysis with the BFGS algorithm and the Newton-CG algorithm. By expanding the dropdown below, you will see the same implementation as above, but wrapped in a function that can be called with different mesh resolutions.

Hide code cell source 
def solve_optimal_problem(N: int, use_newton: bool = False) -> dict[str, float | int]:
    """Solve the optimal control problem for a given mesh resolution.

    Args:
        N: Number of elements in each direction of the mesh.
        use_newton: Whether to use Newton-CG instead of BFGS for the optimization.
    Returns:
        A dictionary containing the results of the optimization problem.
    """
    # Reset tape
    pyadjoint.get_working_tape().clear_tape()
    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)
    mesh.topology.create_connectivity(1, mesh.topology.dim)
    edges_to_refine = dolfinx.mesh.locate_entities(mesh, 1, refinement_region)
    refined_mesh_data = dolfinx.mesh.refine(mesh, edges_to_refine)
    refined_mesh = refined_mesh_data[0]
    refined_mesh = mesh  # Use the original mesh for testing
    tdim = refined_mesh.topology.dim
    del mesh

    V = dolfinx.fem.functionspace(refined_mesh, ("Lagrange", 1))  # type: ignore[arg-type]
    Q = dolfinx.fem.functionspace(refined_mesh, ("Discontinuous Lagrange", 0))  # type: ignore[arg-type]

    f = dolfinx_adjoint.Function(Q, name="Control")
    f.interpolate(lambda x: x[0] + np.sin(2 * x[1]))  # Set intial guess

    uh = dolfinx_adjoint.Function(V, name="State")

    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    kappa = dolfinx.fem.Constant(refined_mesh, dolfinx.default_scalar_type(1.0))  # Thermal diffusivity
    F = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx - f * v * ufl.dx
    a, L = ufl.system(F)

    refined_mesh.topology.create_connectivity(tdim - 1, tdim)
    exterior_facets = dolfinx.mesh.exterior_facet_indices(refined_mesh.topology)
    exterior_dofs = dolfinx.fem.locate_dofs_topological(V, tdim - 1, exterior_facets)
    zero = dolfinx.fem.Constant(refined_mesh, dolfinx.default_scalar_type(0.0))
    bc = dolfinx.fem.dirichletbc(zero, exterior_dofs, V)

    problem = dolfinx_adjoint.LinearProblem(
        a, L, u=uh, bcs=[bc], petsc_options=petsc_options, adjoint_petsc_options=petsc_options
    )
    problem.solve()

    x, y = ufl.SpatialCoordinate(refined_mesh)
    d = 1 / (2 * ufl.pi**2) * ufl.sin(ufl.pi * x) * ufl.sin(ufl.pi * y)

    alpha = dolfinx.fem.Constant(refined_mesh, dolfinx.default_scalar_type(1.0e-6))  # Tikhonov regularization parameter
    J_symbolic = 0.5 * ufl.inner(uh - d, uh - d) * ufl.dx + 0.5 * alpha * ufl.inner(f, f) * ufl.dx
    J = dolfinx_adjoint.assemble_scalar(J_symbolic)

    control = pyadjoint.Control(f)
    Jhat = pyadjoint.ReducedFunctional(J, control)

    optimization_problem = pyadjoint.MoolaOptimizationProblem(Jhat)
    f_moola = DolfinxPrimalVector(f)
    if use_newton:
        solver = moola.NewtonCG(
            optimization_problem,
            f_moola,
            options={
                "gtol": 1e-9,
                "maxiter": 20,
                "display": 0,
                "ncg_hesstol": 0,
            },
        )
    else:
        opts = optimization_options.copy()
        opts["display"] = 0  # Turn down verbosity
        solver = moola.BFGS(optimization_problem, f_moola, options=opts)
    sol = solver.solve()
    f_opt = sol["control"].data
    f.x.array[:] = f_opt.x.array.copy()
    problem.solve(annotate=False)
    f_analytic = 1 / (1 + alpha * 4 * pow(ufl.pi, 4)) * ufl.sin(ufl.pi * x) * ufl.sin(ufl.pi * y)
    u_analytic = 1 / (2 * ufl.pi**2) * f_analytic
    err_u = dolfinx_adjoint.error_norm(u_analytic, uh, norm_type="L2", annotate=False)
    err_f = dolfinx_adjoint.error_norm(f_analytic, f, norm_type="L2", annotate=False)

    num_cells = refined_mesh.topology.index_map(tdim).size_local
    local_h = np.max(refined_mesh.h(tdim, np.arange(num_cells, dtype=np.int32)))
    global_h = refined_mesh.comm.allreduce(local_h, op=MPI.MAX)
    return {
        "Number of iterations": sol["iteration"],
        "h": global_h,
        "L2_u": err_u,
        "L2_f": err_f,
        "J": sol["objective"],
        "|dJ/df|": sol["grad_norm"],
    }

BFGS results#

We run the BFGS algorithm for different mesh resolutions and collect the results.

results_bfgs = []
for N in [16, 32, 64, 128]:
    results_bfgs.append(solve_optimal_problem(N, use_newton=False))

pandas.DataFrame(results_bfgs)
Number of iterations h L2_u L2_f J |dJ/df|
0 67 0.088388 0.000089 0.038032 1.321135e-07 1.249291e-08
1 85 0.044194 0.000022 0.017306 1.259477e-07 7.413259e-09
2 69 0.022097 0.000006 0.008543 1.251555e-07 5.729077e-09
3 63 0.011049 0.000003 0.005042 1.250073e-07 1.280819e-08

We observe that the number of iterations is independent of the mesh resolution, and that the errors in the state and control goes down with a rate of 1.

results_newton_cg = []
for N in [16, 32, 64, 128]:
    results_newton_cg.append(solve_optimal_problem(N, use_newton=True))

Newton-CG results#

We can also check the convergence of an algorithm using Hessian information, like Newton-CG. If we set the tolerance of the Conjugate Gradient solver to zero, the algorithm will converge in one iteration. However, if we set it to a small value, we should observe a low number of iterations, that is also independent of the mesh resolution.

newton_results = pandas.DataFrame(results_newton_cg)
newton_results
Number of iterations h L2_u L2_f J |dJ/df|
0 3 0.088388 0.000089 0.037917 1.321064e-07 2.855536e-11
1 3 0.044194 0.000022 0.017105 1.259415e-07 4.246664e-11
2 3 0.022097 0.000005 0.008276 1.251510e-07 1.109632e-10
3 3 0.011049 0.000001 0.004101 1.249983e-07 2.178553e-10

References#

[Troltzsch10]

Fredi Tröltzsch. Optimal control of partial differential equations: theory, methods, and applications. Volume 112. American Mathematical Soc., 2010.

[Ulb09]

Stefan Ulbrich. Analytical Background and Optimality Theory, pages 1–95. Springer Netherlands, Dordrecht, 2009. doi:10.1007/978-1-4020-8839-1_1.

Hide code cell source 
assert newton_results["Number of iterations"].max() == 3
Contents