# How to compute the gradient of pde constrained optimal control problem with Neuman boundary control

Hello,
I’m currently try to solve an optimal control problem including the poisson equation with Neumann boundary conditions on one part of the boundary and Neumann control on another part of the boundary. The entire problem is formulated as follows:

\min J(u, f) = \frac{1}{2} \int_{\Omega} (u - u_d)^2 dx + \frac{\alpha}{2} \int_{\Gamma_1} q^2 \, ds
s.t.
\begin{align*} -\Delta u = f \, \, & \text{in } \Omega\\ \partial_n u = q \, \, & \text{on } \Gamma_1\\ \partial_n u = 0 \, \, & \text{on } \Gamma_2\\ \partial_n u = 0 \, \, & \text{on } \Gamma_3\\ \partial_n u = 0 \, \, & \text{on } \Gamma_4 \end{align*}

\begin{align*} & \Omega & \text{unit square}\\ & \Gamma_1 & \text{bottom boundary of the square}\\ & \Gamma_2 & \text{left boundary of the square}\\ & \Gamma_3 & \text{top boundary of the square}\\ & \Gamma_4 & \text{right boundary of the square}\\ & q \in L^2(\Gamma_2) & \text{control variable}\\ & u \in H^1(\Omega) & \text{state variable}\\ & \alpha > 0 & \text{penalization parameter}\\ & u_d & \text{desired state}\\ & f & \text{forcing term} \end{align*}

My approach to solve this problem is compute the gradient \partial_q J of the functional J in terms of adjoint variable \lambda for which I have to solve the adjoint problem:

\begin{align*} -\Delta \lambda = u-u_d \, \, & \text{in } \Omega\\ \partial_n \lambda = \, \, & \text{on } \Gamma_1\\ u = 0 \, \, & \text{on } \Gamma_2\\ u = 0 \, \, & \text{on } \Gamma_3\\ u = 0 \, \, & \text{on } \Gamma_4 \end{align*}

with the solution \lambda the gradient can be computed as follows

\partial_q J = \alpha q + \int_{\Gamma_1} \lambda \, ds .

Finally, to compute the the optimal values for q I’m using the simple steepest descent algorithm (using constant stepsize). The code is attached below.

The whole framework I try to implement using fenicsx, but a few problems have arisen that cannot be solved by myself.

The biggest problem for me is the calculation of the gradient. I think my solution is not quite right. Perhaps someone here has a hint to determine the gradient in such a way that only the boundary values for q are affected.

Thanks a lot!
Denis

from mpi4py import MPI
import numpy as np
from petsc4py import PETSc
from dolfinx import fem, mesh
import ufl

# define the mesh as unit square with n mesh elements
n = 64
domain = mesh.create_unit_square(MPI.COMM_WORLD, n, n, mesh.CellType.quadrilateral)

# define constant
alpha = 1e-6

# define the function space as piecewise continous lagrange functions
V = fem.FunctionSpace(domain,("CG", 1))

# define the function space for the control
W = fem.FunctionSpace(domain,("DG",0))

# identifying the facets contained in each boundary and create custom integration measure
boundaries = [(1, lambda x: np.isclose(x, 0)),
(2, lambda x: np.isclose(x, 1)),
(3, lambda x: np.isclose(x, 0)),
(4, lambda x: np.isclose(x, 1))]

# loop through each boundary and create MeshTags identifying the facets for each boundary
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# redefine the boundary measure
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

# define boundary facets
dirichlet_facets = []
dirichlet_facets.append(facet_tag.find(1))
dirichlet_facets.append(facet_tag.find(2))
dirichlet_facets.append(facet_tag.find(4))
dirichlet_facets = np.hstack(dirichlet_facets).astype(np.int32)

# get the indices of the deegre of freedom on the boundary
tdim = domain.topology.dim
fdim = tdim-1
domain.topology.create_connectivity(fdim,tdim)

# define the dirichlet boundary condition
boundary_dofs = fem.locate_dofs_topological(V, fdim, dirichlet_facets)
bc = fem.dirichletbc(PETSc.ScalarType(0.0), boundary_dofs,V)

# define the trial and test function
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# define  the desired  profile for the cost functional
ud    = fem.Function(V)
ud.interpolate(lambda x: 10.0-10.0*x)

# cost function
def objective(uh, ud, q, alpha):
Jt1     = alpha/2 * q ** 2 * ds(3)
Jt2     = 0.5*ufl.inner(uh-ud,uh-ud)*ufl.dx
J_form  = fem.form(Jt1 + Jt2)
J       = fem.assemble_scalar(J_form)
return J

# control variable
q  = fem.Function(W, name='control')
q.interpolate(lambda x: x)

# source term
f = fem.Function(V)

# solution variables for poisson and adjoint state
uh     = fem.Function(V)
lambdh = fem.Function(V)

# definition of the forward problem
# define linear variational problem
L_forward  = f*v*ufl.dx + ufl.inner(q,v)*ds(3)

# Assemble system
A_forward = fem.petsc.assemble_matrix(fem.form(a_forward))
A_forward.assemble()
b_forward =fem.petsc.create_vector(fem.form(L_forward))

# Create Krylov solver
solver = PETSc.KSP().create(A_forward.getComm())
solver.setOperators(A_forward)

# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
A_forward.setNullSpace(nullspace)

# orthogonalize b with respect to the nullspace ensures that
# b does not contain any component in the nullspace
nullspace.remove(b_forward)

# define the adjoint problem lambda as nonlinear problem
lambd = ufl.TrialFunction(V)

# define the adjoint problem as linear problem

lambdg = fem.Function(W)

J = objective(uh,ud,q,alpha)

while(J > 0.1):
# update the right hand side reusing the initial vector
with b_forward.localForm() as b_loc:
b_loc.set(0)
fem.petsc.assemble_vector(b_forward,fem.form(L_forward))

# solve the forward problem
solver.solve(b_forward,uh.vector)
uh.x.scatter_forward()

# update the right hand side reusing the initial vector
loc_ba.set(0)

# apply dirichlet bcs  to the vector

# solve the poisson eq
lambdh.x.scatter_forward()

# interpolate adjoint state on the DG function space
lambdg.interpolate(lambdh)

dJdf = alpha*q.vector + lambdg.vector

# steepest descent with constant beta as stepsize
beta = 0.1
delta_q = -beta*dJdf

q_new = q.vector + delta_q

# determine new value for the design variables
q.interpolate(q_new)

L2_error = fem.form(ufl.inner(uh - ud, uh - ud) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = np.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))

J = objective(uh,ud,q,alpha)

print(J)


This problem is not unique if you only have Neumann conditions.
See: Issues Solving Pure Neumann Problems in dolfinx - #3 by jkrokowski

I also dont think your Adjoint problem is correct, as you have changed your homogeneous Neumann conditions into Dirichlet conditions.

Thank you very much for the quick reply and the advice. I will look at them and correct them. 