Boundary conditions for constrained optimization

Question

I want to solve a PDE-constrained optimization problem $\min f(y,u)$ s.t. $c(y,u)=0$ with homogeneous Dirichlet boundary conditions. In this context it seems like I make some error imposing the boundary conditions.
So this is what I observe:
Everything works fine when I implement the Lagrangian as
$$ L(y,u,p) = f(y,u) + pc(y,u) + 10^9\int_\Gamma y(s)p(s)ds,$$
i. e. if I incorporate Dirichlet boundary conditions with the help of a penalty term.

Now, if I implement the Lagrangian as
$$ L(y,u,p) = f(y,u) + pc(y,u)$$ and explicitly set boundary conditions via
DirichletBC bcy{V[0], y0, boundary}; // boundary condition for state variable
DirichletBC bcp{V[2], p0, boundary}; // boundary condition for adjoint variable
my solver does not converge any more. I assume that in this case calling bc.apply will affect
the diagonal matrix blocks instead of the blocks related to state and adjoint equation.

Am I right about this? If so is there a way to correctly set the boundary conditions for my needs?

Best Regards

P.S.: At least with my Firefox the Preview seems not to recognize MathJax.

Answer 1

DirichletBC.apply will set the relevant diagonal entries to 1 and set off-diagonal entries on the same rows (across all blocks) to zero. This will correctly enforce the boundary conditions for the system of equations. You can use assemble_system instead of assemble if you wish to preserve symmetry.

Here is a complete python example implementing a simple PDE-constrained optimization problem. The problem is solved in one Newton iteration.

from dolfin import *

mesh = UnitSquareMesh(16,16)
a = Constant(1e-4)

V = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, V, V])

bcs = [DirichletBC(W.sub(1), 0., "on_boundary"),
       DirichletBC(W.sub(2), 0., "on_boundary")]

w = Function(W); u, y, p = split(w)

f = (1-y)**2 * dx + a * u**2 * dx
L = f + inner(grad(y), grad(p)) * dx - u * p * dx

dL = derivative(L, w)
H  = derivative(dL, w)

u, y, p = w.split()
solve(H == -dL, w, bcs)

Comments

Thanks for your answer. Do I understand correctly that this yields different matrices describing the differential operator depending on where it is used, since the diagonal of the complete system is considered?
If it is used in a PDE context, then the diagonal entries, associated with dofs on the boundary, of the stiffness matrix are set to one.
If it used in an optimal control context, then the differential operator is not on the diagonal and thus its diagonal entries, related to the boundary, in its discretization will be set to zero.
Thus to clarify my initial question: The application of the boundary conditions affects all blocks, but the stiffness matrix gets modified differently depending on its appearance in a PDE or optimal control context?

---

That is correct.