Optimal control problem

Question

Hi, I'm new to using FEniCS and would like to use it for solving problems in optimal control.

More specifically I have been trying to implement the following simple heat diffusion problem on the unit square

min J(y,u) = 1/2 * Integral (y-y_d)^2 dx + alpha/2 * Integral u^2 dx
subj. to     -Laplace(y) = u     in the domain, and 
                       y = 0     on the boundary

where y is our state, u is the control and y_d is the desired state.

(For future reference: is it possible to use LaTeX on this Q&A forum?)

I have attempted to set this up following some of the demos as examples, however, I receive an error, and I can't seem to figure out what exactly the problem is and how I change the code correctly.

Here is my code. I suspect the error somehow resides in the fact that I have both y and u as TrialFunctions, however, I'm not sure how to get around this. I have done quite a bit of searching around, but have not been able to find any good examples for this kind of problem.

# File: fenicsLinVarSol.py 
from dolfin import *

# Create mesh and define space
mesh = UnitSquareMesh(5, 5, 'crossed')
V = FunctionSpace(mesh, 'Lagrange', 1) # CG - Continuous Galerkin (same)
W = V*V

# Define boundary conditions
y0 = Constant(0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, y0, boundary)

# Define variational problem
(y,u) = TrialFunction(W)
v = TestFunction(V)
f = Constant(0)

yd = Expression("2*sin(3.14*x[0])")
alpha = Constant(0.1)

a = inner(grad(y), grad(v))*dx - u*v*dx
L = f*v*dx

# Define function for the solution
w = Function(W)

# Define goal functional
J = ( ( w[0] - yd )**2 )*dx + alpha*( w[1]**2 )*dx

# Define error tolerance
tol = 1.e-3

problem = LinearVariationalProblem(a, L, w, bc)
solver = AdaptiveLinearVariationalSolver(problem, J)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
solver.parameters["error_control"]["dual_variational_solver"]["symmetric"] = True
solver.solve(tol)

solver.summary()

# Plot solution
plot(u.root_node(), title="Solution on initial mesh")
plot(u.lead_node(), title="Solution on final mesh")
interactive()

Here is what I receive from python when I run it.

~$ python Documents/python/fenicsLinVarSol.py 
Generating forms required for error control, this may take some time...
Replacement expressions must have the same shape as what they replace.
Traceback (most recent call last):
  File "Documents/python/fenicsLinVarSol.py", line 45, in <module>
solver = AdaptiveLinearVariationalSolver(problem, J)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/adaptivesolving.py", line 58, in __init__
ec = generate_error_control(self.problem, goal)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/adaptivesolving.py", line 145, in generate_error_control
(ufl_forms, is_linear) = generate_error_control_forms(problem, goal)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/adaptivesolving.py", line 215, in generate_error_control_forms
forms = generator.generate_all_error_control_forms()
  File "/usr/lib/python2.7/dist-packages/ffc/errorcontrol/errorcontrolgenerators.py", line 111, in generate_all_error_control_forms
(a_R_T, L_R_T) = self.cell_residual()
  File "/usr/lib/python2.7/dist-packages/ffc/errorcontrol/errorcontrolgenerators.py", line 163, in cell_residual
L_R_T = replace(self.weak_residual, {v_h: v_T})
  File "/usr/lib/python2.7/dist-packages/ufl/algorithms/replace.py", line 69, in replace
return map_integrand_dags(Replacer(mapping2), e)
  File "/usr/lib/python2.7/dist-packages/ufl/algorithms/replace.py", line 39, in __init__
"Replacement expressions must have the same shape as what they replace.")
  File "/usr/lib/python2.7/dist-packages/ufl/assertions.py", line 37, in ufl_assert
if not condition: error(*message)
  File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 151, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Replacement expressions must have the same shape as what they replace.

Answer 1

You have a funny way of setting up the system. I would expect the
use of the Lagrange multiplier method to set up a 3x3 system
and consider the optimality condition. You can do this
(roughly) by multiplying ´a´by the test function of the multiplier and add ´a´ to ´J´.
Then differentiate and get the linear problem.

Have a look eg. at the thesis of Ole Elvetun. He implemented many
problems similar to yours in FEniCS. See here:
http://folk.uio.no/kent-and/phdthesis_OLE.pdf
He looked at problems related to precondition rather than adaptivity
but the norm considerations are probably useful.

Comments

I understand that my code may be way off for what I want to solve. I have for the most part simply followed the demo here: http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/auto-adaptive-poisson/python/documentation.html

This demo was not exactly what I wanted, but it has been the only documentation of use of goal-functionals, as I have with my functional J(y,u), so I figured I'd try to modify it.

Is there somewhere I can read more about how to use the Lagrange multiplier method that you mention?

Thank you for your reply though, I will take a look at that thesis now! :)