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], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 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[1])


# 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[0])

# 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
a_forward  = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
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)

Adj   = (ufl.dot(ufl.grad(lambd), ufl.grad(v)) - (uh-ud)*v) * ufl.dx

# define the adjoint problem as linear problem 
a_adjoint        = ufl.lhs(Adj)
L_adjoint        = ufl.rhs(Adj)

a_adjform   = fem.form(a_adjoint) 
L_adjform   = fem.form(L_adjoint) 

A_adjoint = fem.petsc.assemble_matrix(a_adjform, bcs=[bc])
A_adjoint.assemble()
b_adjoint = fem.petsc.create_vector(L_adjform)

solverAdj = PETSc.KSP().create(domain.comm) 
solverAdj.setOperators(A_adjoint)                 
solverAdj.setType(PETSc.KSP.Type.PREONLY)   
solverAdj.getPC().setType(PETSc.PC.Type.LU) 

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
    with b_adjoint.localForm() as loc_ba:
        loc_ba.set(0)
    fem.petsc.assemble_vector(b_adjoint, L_adjform)
    
    # apply dirichlet bcs  to the vector
    fem.petsc.apply_lifting(b_adjoint, [a_adjform], [[bc]])
    b_adjoint.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,   mode=PETSc.ScatterMode.REVERSE) 
    fem.petsc.set_bc(b_adjoint, [bc])
   
    # solve the poisson eq
    solverAdj.solve(b_adjoint, lambdh.vector)
    lambdh.x.scatter_forward()       
    
    # interpolate adjoint state on the DG function space
    lambdg.interpolate(lambdh)
  
    dJdf = alpha*q.vector + lambdg.vector
    dJdf.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
  
    # 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. :grinning: