Problem with control constraint optimization with NonlinearVariationalSolver

Question

Hi,
I am considering the simple optimization problem
$0.5 \cdot \| y - y_d \|^2 + \alpha \|u^2\|^2$
such that $-\Delta y = u$ in $\Omega$ and zero Dirchlet boundary conditions. This can be easy solved with the use of the Lagrange-Funktion

from dolfin import *
mesh = UnitSquareMesh(100,100)
U = FunctionSpace(mesh, "CG", 1)
Z = MixedFunctionSpace([U, U, U])

z = Function(Z)
(y, p, u) = split(z)

F = inner(grad(y) , grad(p))*dx - u*p*dx

yd = Expression(...)

J = inner(y-yd, y-yd)*dx + Constant(alpha)*inner(u,u)*dx

bc_y = DirichletBZ(Z.sub(0), 0.0, "on_boundary")
bc_p = DirichletBZ(Z.sub(1), 0.0, "on_boundary")
bcs = [bc_y, bc_p]

L = J + F
kkt = derivative(L, z, TestFunction(Z))

solve(kkt == 0, z, bcs)

plot(y)
plot(u)
interactive()

I now want to add the control constaints $-1 \leq u \leq 1$. According to several threads here I should use PETScSNESSolver. I am not sure how to change the Form of F properly to NonVariationalForm F(u,v):

Second: How do i impose box constraints only to u (as part of z) )

from dolfin import *
mesh = UnitSquareMesh(100,100)
U = FunctionSpace(mesh, "CG", 1)
Z = MixedFunctionSpace([U, U, U])
v = TestFunction(U)

z = Function(Z)
(y, p, u) = split(z)

F = inner(grad(y) , grad(p))*v*dx - u*p*v*dx

yd = Expression(...)

J = inner(y-yd, y-yd)*v*dx + Constant(alpha)*inner(u,u)*v*dx

bc_y = DirichletBZ(Z.sub(0), 0.0, "on_boundary")
bc_p = DirichletBZ(Z.sub(1), 0.0, "on_boundary")
bcs = [bc_y, bc_p]

L = J + F
kkt = derivative(L, z, TestFunction(Z))
d_kkt = derivative(kkt, z, TestFunction(Z))

# How to impose correctly? This is obviously wrong?
constraint_u = Constant(-1)
constraint_l = Constant(1)
umin = interpolate(constraint_l, Z.sub(2))
umax = interpolate(constraint_u, Z.sub(2))

snes_solver_parameters = {"nonlinear_solver": "snes",
                      "snes_solver": {"linear_solver": "lu",
                                      "maximum_iterations": 20,
                                      "report": True,
                                      "error_on_nonconvergence": False}}

problem = NonlinearVariationalProblem(kkt, u, bc, d_kkt)
problem.set_bounds(umin, umax)

solver  = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
info(solver.parameters, True)

(iter, converged) = solver.solve()

plot(y)
plot(u)
interactive()

Of course the code above is not working - the definition of umin and umax is wrong:

Error: Unable to create function
Error: Connot be created from subspace. Consider collapsing the function space

How do i have to change the code to make it run? Or: What am i doing wrong.

Best regards

Answer 1

You need to pass bounds for the entire mixed space function. Since you only want to impose bounds for the state component, you can set the remaining bounds to $\pm \infty$. You can FunctionAssigner for this as suggested by Kent.

See complete example below.

from dolfin import *
import numpy

N = 64
alpha = 1e-4
d = Constant(0.9)

y_lower = Constant(0.0) # lower bound for state
y_upper = Constant(1.0) # upper bound 

mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CG", 1)
W = VectorFunctionSpace(mesh, "CG", 1, dim = 3)

w = Function(W)
u, y, p = split(w) # = (control, state, adjoint)

# Lagrangian
L = (y - d)**2 * dx + Constant(alpha) * u**2 * dx \
  + inner(grad(y), grad(p)) * dx - u * p * dx

# First derivative
J = derivative(L, w, TestFunction(W))

# second derivative
H = derivative(J, w, TrialFunction(W))

# bcs
bc_y = DirichletBC(W.sub(1), 0.0, "on_boundary")
bc_p = DirichletBC(W.sub(2), 0.0, "on_boundary")
bcs = [bc_y, bc_p]


snes_solver_parameters = {"nonlinear_solver": "snes",
                          "snes_solver": {"linear_solver": "lu",
                                          "maximum_iterations": 20,
                                          "report": True,
                                          "error_on_nonconvergence": False}}

problem = NonlinearVariationalProblem(J, w, bcs, H)

# handle the bounds
lower = Function(W)
upper = Function(W)

ninfty = Function(V); ninfty.vector()[:] = -numpy.infty
pinfty = Function(V); pinfty.vector()[:] =  numpy.infty

fa = FunctionAssigner(W, [V, V, V])
fa.assign(lower, [ninfty, interpolate(y_lower, V), ninfty])
fa.assign(upper, [pinfty, interpolate(y_upper, V), pinfty])

# set bounds and solve
problem.set_bounds(lower, upper)

solver  = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
info(solver.parameters, True)

(iter, converged) = solver.solve()

u, y, p = w.split()

parameters["plotting_backend"] = "matplotlib"
from matplotlib import pyplot
pyplot.figure()
ax = plot(u)
pyplot.title("Optimized control")
pyplot.colorbar(ax)

pyplot.figure()
ax = plot(y)
pyplot.title("Optimized state")
pyplot.colorbar(ax)

pyplot.show()

Comments

Thank you!

FunctionAssigner was exactly what i was looking for!

Still I have another question - related to this topic here, so I don't open a new thread. In my Lagrange-Function I have some parameter: The value alpha and a function lambda. In the end i want to implement an iterative method, so just this two parameter are changing.

In C++ and UFL i would handle this via Coefficient, and assign J.lmbda = lmbda, H.lmbda = lmbda. Is there a way to copy this in Python? I don't want to call FFC JIT in every iteration to compute the derivates and construct the solver.

The best way would be if I am able to set up the solver just one time in the beginning, and in each iteration i just assign the parameters and call solve. Is this somehow possible?

The following code is working, but very inefficient

k = 0
while k<100:
    w = Function(W)
    u, y, p = split(w) # = (control, state, adjoint)

    # Lagrangian
    L = 0.5*(y - d)**2 * dx + Constant(alpha) * ( 0.5*u**2 * dx - u*lmbda*dx ) \
        + inner(grad(y), grad(p)) * dx - u * p * dx

    # First derivative
    J = derivative(L, w, TestFunction(W))

    # second derivative
    H = derivative(J, w, TrialFunction(W))


    problem = NonlinearVariationalProblem(J, w, bcs, H)
    # set bounds and solve
    problem.set_bounds(lower, upper)

    solver  = NonlinearVariationalSolver(problem)
    solver.parameters.update(snes_solver_parameters)

    (iter, converged) = solver.solve()

    u, y, p = w.split()

    # update lambda
    lmbda = 1.0/alpha*p + lmbda

        k  = k+1

---

Try making the scalar parameter a Constant in the form, i.e.

# Lagrangian
L = 0.5*(y - d)**2 * dx \
    + Constant(alpha) * (0.5 * u**2 * dx - u * Constant(lmbda)*dx) \
    + inner(grad(y), grad(p)) * dx - u * p * dx

---

Hi,

lmbda is a Function, not a constant, so this is not working. Im not sure how to avoid recompilation.

Answer 2

You should use FunctionAssigner.