In this tutorial we solve the optimal control problem
$$\min J(y, u) = \frac{1}{2} \int_{\Omega} (y - y_d)^2 dx + \frac{\alpha}{2} \int_{\Omega} u^2 dx$$ s.t. $$\begin{cases} - \Delta y = f + u & \text{in } \Omega\\ y = 1 & \text{on } \partial\Omega \end{cases}$$
where $$\begin{align*} & \Omega & \text{unit square}\\ & u \in L^2(\Omega) & \text{control variable}\\ & y \in H^1(\Omega) & \text{state variable}\\ & \alpha > 0 & \text{penalization parameter}\\ & y_d & \text{a piecewise constant desired state}\\ & f & \text{forcing term} \end{align*}$$ using an adjoint formulation solved by a one shot approach.
The test case is from section 5.1 of
F. Negri, G. Rozza, A. Manzoni and A. Quarteroni. Reduced Basis Method for Parametrized Elliptic Optimal Control Problems. SIAM Journal on Scientific Computing, 35(5): A2316-A2340, 2013.
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
import multiphenicsx.fem
import multiphenicsx.fem.petsc
L1 = 1.0
L2 = 3.0
H = 1.0
mesh_size = 0.05
gmsh.initialize()
gmsh.model.add("mesh")
p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(L1, 0.0, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(L1 + L2, 0.0, 0.0, mesh_size)
p3 = gmsh.model.geo.addPoint(L1 + L2, H, 0.0, mesh_size)
p4 = gmsh.model.geo.addPoint(L1, H, 0.0, mesh_size)
p5 = gmsh.model.geo.addPoint(0.0, H, 0.0, mesh_size)
l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p4)
l2 = gmsh.model.geo.addLine(p4, p5)
l3 = gmsh.model.geo.addLine(p5, p0)
l4 = gmsh.model.geo.addLine(p1, p2)
l5 = gmsh.model.geo.addLine(p2, p3)
l6 = gmsh.model.geo.addLine(p3, p4)
line_loop_rectangle_left = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3])
line_loop_rectangle_right = gmsh.model.geo.addCurveLoop([l4, l5, l6, -l1])
rectangle_left = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_left])
rectangle_right = gmsh.model.geo.addPlaneSurface([line_loop_rectangle_right])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l0, l4, l5, l6, l2, l3], 1)
gmsh.model.addPhysicalGroup(2, [rectangle_left], 1)
gmsh.model.addPhysicalGroup(2, [rectangle_right], 2)
gmsh.model.mesh.generate(2)
Info : Meshing 1D... Info : [ 0%] Meshing curve 1 (Line) Info : [ 20%] Meshing curve 2 (Line) Info : [ 30%] Meshing curve 3 (Line) Info : [ 50%] Meshing curve 4 (Line) Info : [ 60%] Meshing curve 5 (Line) Info : [ 80%] Meshing curve 6 (Line) Info : [ 90%] Meshing curve 7 (Line) Info : Done meshing 1D (Wall 0.000685929s, CPU 0.001124s) Info : Meshing 2D... Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay) Info : [ 50%] Meshing surface 2 (Plane, Frontal-Delaunay) Info : Done meshing 2D (Wall 0.0393158s, CPU 0.039379s) Info : 1967 nodes 3958 elements
mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()
# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
boundaries_1 = boundaries.indices[boundaries.values == 1]
# Define associated measures
dx = ufl.Measure("dx")(subdomain_data=subdomains)
viskex.dolfinx.plot_mesh(mesh)
viskex.dolfinx.plot_mesh_tags(mesh, subdomains, "subdomains")
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")
Y = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
U = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
Q = Y.clone()
(y, u, p) = (ufl.TrialFunction(Y), ufl.TrialFunction(U), ufl.TrialFunction(Q))
(z, v, q) = (ufl.TestFunction(Y), ufl.TestFunction(U), ufl.TestFunction(Q))
alpha = 0.01
y_d_1 = 1.0
y_d_2 = 0.6
ff = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0))
bc0 = petsc4py.PETSc.ScalarType(0)
bc1 = petsc4py.PETSc.ScalarType(1)
a = [[ufl.inner(y, z) * dx, None, ufl.inner(ufl.grad(p), ufl.grad(z)) * dx],
[None, alpha * ufl.inner(u, v) * dx, - ufl.inner(p, v) * dx],
[ufl.inner(ufl.grad(y), ufl.grad(q)) * dx, - ufl.inner(u, q) * dx, None]]
f = [ufl.inner(y_d_1, z) * dx(1) + ufl.inner(y_d_2, z) * dx(2),
None,
ufl.inner(ff, q) * dx]
a[2][2] = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(p, q) * dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)), v) * dx
a_cpp = dolfinx.fem.form(a)
f_cpp = dolfinx.fem.form(f)
bdofs_Y_1 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_1)
bdofs_Q_1 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_1)
bc = [dolfinx.fem.dirichletbc(bc1, bdofs_Y_1, Y),
dolfinx.fem.dirichletbc(bc0, bdofs_Q_1, Q)]
(y, u, p) = (dolfinx.fem.Function(Y), dolfinx.fem.Function(U), dolfinx.fem.Function(Q))
J = (0.5 * ufl.inner(y - y_d_1, y - y_d_1) * dx(1) + 0.5 * ufl.inner(y - y_d_2, y - y_d_2) * dx(2)
+ 0.5 * alpha * ufl.inner(u, u) * dx)
J_cpp = dolfinx.fem.form(J)
# Extract state forms from the optimality conditions
a_state = ufl.replace(a[2][0], {q: z})
f_state = ufl.replace(f[2], {q: z})
a_state_cpp = dolfinx.fem.form(a_state)
f_state_cpp = dolfinx.fem.form(f_state)
bc_state = [bc[0]]
# Assemble the linear system for the state
A_state = dolfinx.fem.petsc.assemble_matrix(a_state_cpp, bcs=bc_state)
A_state.assemble()
F_state = dolfinx.fem.petsc.assemble_vector(f_state_cpp)
dolfinx.fem.apply_lifting(F_state, [a_state_cpp], [bc_state])
F_state.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(F_state, bc_state)
# Solve
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A_state)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_state, y.x.petsc_vec)
y.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()
<petsc4py.PETSc.KSP at 0x7f5aa82ff920>
J_uncontrolled = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(J_cpp), op=mpi4py.MPI.SUM)
print("Uncontrolled J =", J_uncontrolled)
assert np.isclose(J_uncontrolled, 0.24)
Uncontrolled J = 0.24000000000000202
viskex.dolfinx.plot_scalar_field(y, "uncontrolled state")
# Assemble the block linear system for the optimality conditions
A = dolfinx.fem.petsc.assemble_matrix_block(a_cpp, bcs=bc)
A.assemble()
F = dolfinx.fem.petsc.assemble_vector_block(f_cpp, a_cpp, bcs=bc)
# Solve
yup = dolfinx.fem.petsc.create_vector_block(f_cpp)
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F, yup)
yup.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()
<petsc4py.PETSc.KSP at 0x7f5aa83a6b60>
# Split the block solution in components
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(yup, [Y.dofmap, U.dofmap, Q.dofmap]) as yup_wrapper:
for yup_wrapper_local, component in zip(yup_wrapper, (y, u, p)):
with component.x.petsc_vec.localForm() as component_local:
component_local[:] = yup_wrapper_local
J_controlled = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(J_cpp), op=mpi4py.MPI.SUM)
print("Optimal J =", J_controlled)
assert np.isclose(J_controlled, 0.158485065)
Optimal J = 0.15848479530337567
viskex.dolfinx.plot_scalar_field(y, "state")
viskex.dolfinx.plot_scalar_field(u, "control")
viskex.dolfinx.plot_scalar_field(p, "adjoint")