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)