PDE-constrained optimisation

PDE-constrained optimisation problems are problems of the form:

\[ \begin{align}\begin{aligned}\min_{u,m} J(u, m)\\\mathrm{subject~to}\\F(u, m) = 0\\l_u \le m \le l_b\\g(m) \le 0\end{aligned}\end{align} \]

where \(m\) contains the optimisation variables, \(J\) is a real valued objective functional, \(F(u, m) = 0\) is the PDE with solution \(u\). The bounds and inequality constraints can be used to restrict the feasible optimisation variables.

For an introduction to the mathematics, see the chapter in the mathematical background.

The reduced functional

While it is possible to solve the optimisation problem above directly, we often prefer to form the so-called reduced problem. Given that for every \(m\) the PDE yields a unique solution \(u\), we can define a solution operator \(u(m)\). Substituting this operator into the optimisation problem yields the reduced problem:

\[ \begin{align}\begin{aligned}\min_{m} J(u(m), m)\\\mathrm{subject~to}\\l_u \le m \le l_b\\g(m) \le 0\end{aligned}\end{align} \]

The advantage of solving this formulation is that the PDE-constraint is exactly satisfied at each optimisation iteration. In particular, the optimisation loop can be terminated as soon as the functional is sufficiently reduced by the optimisation algorithm, without any feasibility iterations.

The functional in the reduced form can be seen as a function that only depends on the optimisation variable m, that is:

\[\tilde J(m) := J(u(m), m)\]

The definition of this reduced functional \(\tilde J\) is the first step of solving an optimisation problem with dolfin-adjoint. It is created with:

reduced_functional = ReducedFunctional(J, m)

where J is an AdjFloat and m is a Control.

Important: ReducedFunctional works by replaying the simulation record of dolfin-adjoint. Therefore, make sure that you execute the forward model once before using it.

Solving the optimisation problem

Once the reduced functional is defined, we are only one step away from solving the optimisation problem:

m_opt = minimize(reduced_functional)

or if a maximization problem is to be solved:

m_opt = maximize(reduced_functional)

By default, the optimisation problem is solved using limited memory BFGS method with bound support.

Features

Important: Please make sure that you have scipy >= 0.11 installed. Older scipy versions are only partly supported and require different arguments. You can check your scipy version with

import scipy
print scipy.__version__

Choosing the optimisation algorithm

The optimisation module currently supports following optimisation algorithms:

  • CG: The nonlinear conjugate gradient algorithm.

  • BFGS: The Broyden–Fletcher–Goldfarb–Shanno (BFGS) method.

  • L-BFGS-B: A limited memory BFGS implementation with bound support.

  • SLSQP: The sequential least squares quadratic programming algorithm.

  • TNC: The truncated Newton algorithm with bound support.

  • Nelder-Mead: The Simplex algorithm (gradient-free).

  • Newton-CG: The truncated Newton algorithm.

  • Anneal: The simulated annealing method (gradient-free).

  • COBYLA: Constrained optimization by linear approximation.

  • Powell: The Powell’s method (gradient-free).

More details about the algorithms can be found on the scipy.optimize web page.

This list can be generated by calling:

print_optimization_methods()

By default, the framework uses the L-BFGS-B method. A different algorithm can be selected by adding the method argument to minimize or maximize and providing one of the names from the list above, e.g.:

m_opt = minimize(reduced_functional, method = 'SLSQP')

Callbacks

Often one wants to add a callback function that is executed after every optimisation iteration, for example to save or plot the functional or parameter values. The optimisation framework provides two ways how this can be achieved.

Option 1

One can attach callbacks functions to ReducedFunctional object which are executed whenever the functional is evaluated. There are separate callbacks for functional evaluation and functional gradient evaluation.

The following code example prints the functional value, functional gradient and the associated scalar parameter:

def eval_cb(j, m):
  print "j = %f, m = %f." % (j, float(m))

def derivative_cb(j, dj, m):
  print "j = %f, dj = %f, m = %f." % (j, dj, float(m))

reduced_functional = ReducedFunctional(J, Control(nu),
                                       eval_cb = eval_cb,
                                       derivative_cb = derivative_cb)

In most gradient-based optimisation methods, the gradient is evaluated at the beginning of a new optimisation iteration. Hence, if one wants to plot the progress of the optimisation, the derivative callback is the natural choice.

Option 2

Alternatively, one can attach a callback function to the minimize (or maximize) routine (see below). However, this callback takes only the parameter value as an argument and therefore this method is not suitable if one wants to plot the functional values during the optimisation.

def iter_cb(m):
  print "m = ", m

m_opt = minimize(reduced_functional, method = 'SLSQP', callback = iter_cb)

Advanced optimisation options

Each optimisation algorithm supports different features and hence has different configuration options. To be able to access these options, any arguments that are unknown to minimize or maximize will be passed to the optimisation algorithm.

The most relevant options that can be used with all supported optimisation methods are:

  • tol: Tolerance for termination. For detailed control, use solver-specific options.

  • options: A dictionary of solver options. All methods accept the following generic options:

    • maxiter: Maximum number of iterations to perform.

    • disp: Set to True to print convergence messages.

    • gtol: The iteration loops stops if the gradient norm drops below this tolerance.

more info For method-specific options, see scipy’s function show_options(‘minimize’, method).

For example:

m_opt = minimize(reduced_functional, method = 'SLSQP', tol = 1e-10, options = {'disp': True})

Multiple parameters

The optimisation module can handle multiple optimisation parameters. Simply pass a list of parameters to ReducedFunctional:

reduced_functional = ReducedFunctional(J, [m1, m2, ...])
m_opt = minimize(reduced_functional)

Bounds

If the optimisation algorithm supports bounds of the form \(b_u < m < b_u\) this functionality can be used by adding the bounds argument to minimize or maximize.

reduced_functional = ReducedFunctional(J, m)
m_opt = minimize(reduced_functional, bounds = (m_lb, m_ub))

where m_lb and m_ub are objects of the same type than the parameter that contain the lower and the upper bound values.

If the bounds are constants, a set of floats can be passed alternatively, e.g.:

reduced_functional = ReducedFunctional(J, m)
m_opt = minimize(reduced_functional, bounds = (0.0, 1.0))

In the case where multiple parameters are optimised, the bound parameter must consist of a list whose elements contains the bounds for each parameter, i.e.

reduced_functional = ReducedFunctional(J, [m1, m2, ...])
m_opt = minimize(reduced_functional, bounds = [(m1_lb, m1_ub), (m2_lb, m2_ub), ...])

where each of the m1_lb, m1_ub, m2_lb, … are objects of the same type as the parameter.

Example

The following example shows the code for solving the optimal control of the heat equation:

""" Solves the optimal control problem for the heat equation """
from dolfin import *
from dolfin_adjoint import *

# Setup
n = 200
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), n, n)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V, name="State")
m = Function(V, name="Control")
v = TestFunction(V)

# Run the forward model once to create the simulation record
F = (inner(grad(u), grad(v)) - m*v)*dx
bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u, bc)

# The functional of interest is the normed difference between desired
# and simulated temperature profile
x = SpatialCoordinate(mesh)
u_desired = exp(-1/(1-x[0]*x[0])-1/(1-x[1]*x[1]))
J = assemble((0.5*inner(u-u_desired, u-u_desired))*dx)

# Run the optimisation
reduced_functional = ReducedFunctional(J, Control(m))
# Make sure you have scipy >= 0.11 installed
m_opt = minimize(reduced_functional, method = "L-BFGS-B",
                 tol=2e-08, bounds = (-1, 1), options = {"disp": True})

more info Download the optimisation code.

This prints the following output that contains various information, such as the final functional value:

$ python optimal_control.py
...
N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
40401      6      8      6     0     0   8.869D-09   1.601D-05
F =   1.6009498922520026E-005

CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL