This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Problem with control constraint optimization with NonlinearVariationalSolver

+1 vote

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

asked Apr 22, 2016 by Kokett FEniCS Novice (300 points)

2 Answers

+3 votes
 
Best answer

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()
answered Apr 26, 2016 by Magne Nordaas FEniCS Expert (13,820 points)
selected Apr 27, 2016 by Kokett

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.

0 votes

You should use FunctionAssigner.

answered Apr 25, 2016 by Kent-Andre Mardal FEniCS Expert (14,380 points)
...