Hi, I’m new to FEniCSx and I need help. I am trying to solve optimal control using FEniCSx but my code is not converging i.e (f is not converging to f_star), and the constraint problem solution (u) is zero everywhere. I need help to figure out what might be the problem.
The problem is to minimize the objective function:
The goal is to minimize the objective function:
$$ E(f) = \int_{\Omega} (u - u_{\text{star}})^2 + k |\nabla f|^2 , dx $$
subject to the constraint:
$$ -\Delta u = f \quad \text{and} \quad u = 0 \text{ on } \partial \Omega $$
Here, ( u_{\text{star}} ) is computed by solving the following “star problem” with ( f_{\text{star}} = 1 ):
$$ -\Delta u_{\text{star}} = 1 \quad \text{and} \quad u_{\text{star}} = 0 \text{ on } \partial \Omega $$
The corresponding adjoint equation is:
$$ -\Delta \lambda = -2(u - u_{\text{star}}) \quad \text{with} \quad \lambda = 0 \text{ on } \partial \Omega $$
Here is the code:
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx import fem, mesh, plot,io
import ufl
from dolfinx.fem.petsc import LinearProblem
from ufl import dx,grad,inner
from petsc4py.PETSc import ScalarType
from pathlib import Path
Max_iteration=10000 # Number of iteration
Tol=1e-3 # set the tolorance
k=1e-4 #penalization term
nx=ny=60
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(-2, -2), (2, 2)], [nx, ny], mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))
mu = Constant(domain, ScalarType(1e-5)) # fixed stepsize
# boundary condition
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)
#now we compute the constrain problem -laplacian(u)=f with u=0 on the bpundary
u=ufl.TrialFunction(V)
v=ufl.TestFunction(V)
f_n=fem.Function(V) #initial guess
f_n.x.array[:]=0.0 # setting the initial starting point which is also a source term for the constrain pb
# define the constraint problem
a=inner(grad(u), grad(v))*dx # bilinear form
L= f_n*v*dx # The linear form
# solve the constraint problem
u_sol=fem.Function(V)
problem= fem.petsc.LinearProblem(a,L, bcs=[bc])
u_sol=problem.solve()
#Test my code with this # equatuion -laplacian(u_star)=f_star
u_star=ufl.TrialFunction(V)
f_star = fem.Constant(domain, ScalarType(1.0)) ## Constant function of 1 across the domain
a_star= inner (grad(u_star), grad(v))*dx # Bilinear form
L_star = f_star*v*dx # Linear form
u_starsol=fem.Function(V)
problem_star=fem.petsc.LinearProblem( a_star, L_star, bcs=[bc])
u_starsol=problem_star.solve() #solve u_star problem
# next we define the adjoint problem
lambdaa=ufl.TrialFunction(V)
a_adj= inner(grad(lambdaa), grad(v))*dx # Bilinear form
L_adj= -2*inner(u_sol- u_starsol,v)*dx #Linear form right handside
lambdaa_sol=fem.Function(V)
problem_adj=fem.petsc.LinearProblem(a_adj, L_adj, bcs=[bc])
lambdaa_sol=problem_adj.solve()
# objective function
def compute_objective(u_sol, u_starsol, f_n, k):
# Term 1: (u - u_star)^2
term1 = ufl.inner(u_sol - u_starsol, u_sol - u_starsol) * ufl.dx
# Term 2: k * |grad(f)|^2
term2 = k *ufl.inner(grad(f_n), grad(f_n)) * ufl.dx
# Total objective function E(f)
E_form = fem.form(term1 + term2)
E_f= fem.assemble_scalar(E_form)
return E_f
#steepest descent loop
for n in range (Max_iteration):
# we now define the the steepest descent problem to compute f_n+1. replace f_n+1 by f
f = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a_sd = f * v * dx # bilinear form for the steepest descent problem
L_sd = f_n * v * dx + 2 * mu * k * inner(grad(f_n), grad(v)) * dx - lambdaa_sol*v*dx # right hand side steepest descent problem
f_sol=fem.Function(V)
problem_sd=fem.petsc.LinearProblem(a_sd, L_sd, bcs=[bc])
f_sol=problem_sd.solve()
# Compute the objective function E(f) for this iteration
E_f = compute_objective(u_sol, u_starsol, f_n, k)
print(f'Step {n}: Objective function E(f) = {E_f:e}')
print(f"Step {n}: f_n min={f_n.x.array.min():e}, max={f_n.x.array.max():e}")
print(f"u_sol: min={u_sol.x.array.min():e}, max={u_sol.x.array.max():e}")
print(f"u_starsol: min={u_starsol.x.array.min():e}, max={u_starsol.x.array.max():e}")
print(f"lambdaa_sol: min={lambdaa_sol.x.array.min():e}, max={lambdaa_sol.x.array.max():e}")
# Calculate the L2 error (difference) between f_sol and f_star
difference = np.sqrt(fem.assemble_scalar(fem.form(inner(f_sol- f_star, f_sol - f_star) * dx)))
print(f'Step {n}: L2 error= {difference:e}, min={f_sol.x.array.min():e} max= {f_sol.x.array.max():e}')
# Check for convergence
if difference < Tol:
print('Convergence is achieved.')
break
else:
# Update the previous solution
f_n.x.array[:] = f_sol.x.array[:]
