Equality Constraint For Optimal Control Problem

General
1

Dear FEniCS Community,

I am currently solving optimal control problems (OCP) with PDEs as constraints in dolfin_adjoint. The problem is solved by scipy.
Now I want to add an equality constraint to one of the cells of the controls, i.e. the one cell which is targeted should have the the given value at the end.
The code so far uses the class EqualityConstraint to set the value of one cell to 1.
In this example the cell with index 0 should be set to 1.

Code 
from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
from matplotlib import cm

# Create mesh
n = 4
mesh = UnitSquareMesh(n, n)

# Define function spaces and functions
Y = FunctionSpace(mesh, 'CG', 1)
U = FunctionSpace(mesh, 'DG', 0)
y_ = TrialFunction(Y)
v_ = TestFunction(Y)
f = interpolate(Expression("x[0]*x[1]", degree=3), U)

# Bonudary condition
bc = DirichletBC(Y, 0.0, "on_boundary")

# Weak formulation and solving Poisson
A = dot(grad(y_), grad(v_)) * dx
L = f * v_ * dx
y = Function(Y)
solve(A == L, y, bc)

# Set up OCP
w = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=3)
d = 1 / (2 * pi ** 2)
y_d = Expression("d*w", d=d, w=w, degree=3)
lam = Constant(1e-5)
J = assemble( (0.5*inner(y-y_d, y-y_d))*dx + 0.5*lam*f**2*dx )

# Equality constraint for exactly 1 cell
class EqCon(EqualityConstraint):
    def __init__(self, target_cell_index):
        self.target_cell_index = target_cell_index  # Index of target cell

    def value(self, u):
        cell_value = u.vector().get_local()[self.target_cell_index] # get value at target cell
        return cell_value - 1.0 # return the difference to 1 (should equal 0 for the constraint)
    
    def evaluate(self, u):
        return self.value(u)

# Solve OCP
eq_con = EqCon(0)
control = Control(f)
rf = ReducedFunctional(J, control)
u_opt = minimize(rf, bounds = (0.0, 1.0),
                 constraints = eq_con,
                 tol=1e-15,
                 method="L-BFGS-B",
                 options={"maxiter": 20})

# Plot control
print("u_opt:\n", u_opt.vector()[:])
p = plot(u_opt, title="control", vmin=0.0, vmax=1.0, cmap=cm.gist_yarg)
plt.colorbar(p)
plt.gca().set_xticklabels([])
plt.gca().set_yticklabels([])
plt.savefig("controlExample.png", dpi=300, bbox_inches="tight", pad_inches=0)
plt.close()

But the plot for the optimal control looks like this:
|

What I get
What I get1348×1196 15.2 KB
|
What I want
What I want1348×1196 15.3 KB
|

The left control is what I get and the right picture is what I want (for the right picture I just changed the value of the 0th control to 1 afterwards just for showcase).
Is it even possible to set constraints like this?

Thanks for any help or suggestion!

2

Have you considered the demo: dolfin-adjoint/docs/source/documentation/poisson-topology/poisson-topology.rst at 5f020398c59173d879f9dbf3fda36e8acafc4afb · dolfin-adjoint/dolfin-adjoint · GitHub