Pde constraint optimization solver

Hello,
I’ve added my working code below, it’s a simple code but it’s not working correctly.
I want to do pde constraint optimization, I minimize the following functional, first solve my model for a prescribed displacement and then solve it for F_i, desired force, to minimize it. (I want my force-disp graph to pass through the desired force points after the optimization).
my functional:

I left the code running below, but the result never provides correctly. I also tried different solvers.
I would appreciate it if you could help me because I couldn’t get it to work somehow.
Thank you in advance!

from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
import numpy as np

try:
	import cyipopt
except ImportError:
	print("""This example depends on IPOPT and pyipopt. \
  When compiling IPOPT, make sure to link against HSL, as it \
  is a necessity for practical problems.""")
	raise

# Set log level
set_log_level(LogLevel.WARNING)

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True,
               "eliminate_zeros": True,
               "precompute_basis_const": True,
               "precompute_ip_const": True}
l1 = 1e-1
l2 = 1e-9
normal = Constant((0, 1))
chi0 = 0.9
delta_t = 1e-1
E1, nu1 = 30.5, 0.4  # base material
E2, nu2 = 1e-10, 0.0  # void material
mu1, lmbda1 = Constant(E1 / (2 * (1 + nu1))), Constant(E1 * nu1 / ((1 + nu1) * (1 - 2 * nu1)))
mu2, lmbda2 = Constant(E2 / (2 * (1 + nu2))), Constant(E2 * nu2 / ((1 + nu2) * (1 - 2 * nu2)))
max_disp = [-0.1]
f_desired = [-0.11]
store_u = np.linspace(0, 1.0 * min(max_disp), 100)
displacement_by_step = store_u
eps = Constant(1e-6)
c1, c2 = 10, 1e-4

# Create mesh and define function space
n = 32
mesh = UnitSquareMesh(n, n)
W = FunctionSpace(mesh, "CG", 2)
V_vector = VectorFunctionSpace(mesh, "CG", 2)
W_tensor = TensorFunctionSpace(mesh, "CG", 2)
force_FS = VectorFunctionSpace(mesh, "Real", 0)
df = TestFunction(force_FS)

def bcs_(disp):
    
    def boundary_bot(x, on_boundary):
        tol = 1E-7
        return on_boundary and near(x[1], 0.0, tol)

    bc_bot = DirichletBC(V_vector, [0.0, 0.0], boundary_bot)
    
    def boundary_top(x, on_boundary):
        tol = 1E-7
        return on_boundary and (near(x[1], 1.0, tol))
    
    bc_top = DirichletBC(V_vector, disp, boundary_top)
 
    bcs = [bc_bot, bc_top]
    boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundary_subdomains.set_all(0)
    AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)
    dss = ds(subdomain_data=boundary_subdomains)
    return bcs, dss

def H(a, l):
    return 0.5 * (1 + a / (l + (a ** 2) ** 0.5 ))

# initial chi expression
chi0 = Constant(0.9)
chi = Function(W)
chi.interpolate(chi0)
file = File(
    "./Results" + "/chi_initial.pvd");
file << chi

# Stored strain energy density (compressible neo-Hookean model)
def psi(mu, Ic, J_, lmbda):
    return (mu / 2) * (Ic - 3) - mu * ln(J_) + (lmbda / 2) * (ln(J_)) ** 2

def kinematics(u):
	v = TestFunction(V_vector)  # Test function
	d = u.geometric_dimension()
	I = Identity(d)  # Identity tensor
	F = variable(I + grad(u))  # Deformation gradient
	C = F.T * F  # Right Cauchy-Green tensor
	# Invariants of deformation tensors
	Ic = tr(C)
	J_ = det(F)
	return Ic, J_, v, F

def optimization(chi):
    max_error = 0
    count = 0
    for i in max_disp:
        u0 = Function(V_vector)
        u = Function(V_vector)
        step_size = np.linspace(0, i, 100, endpoint=True)
        for j in step_size:
            disp = Constant([0.0, j])
            bcs, dss = bcs_(disp)
            Ic, J_, v, F = kinematics(u)
            # Stored strain energy density (compressible neo-Hookean model)
            Psi = H((chi - 0.5), l2) * psi(mu1, Ic, J_, lmbda1) + (1 - H((chi - 0.5), l2)) * psi(mu2, Ic, J_, lmbda2)
            # Total potential energy
            Pi = Psi * dx
            # Compute first variation of Pi (directional derivative about u in the direction of v)
            dPi = derivative(Pi, u, v)
            file2 = File("./Results" + "/stress-opt.pvd");
            u.assign(u0)
            solve(dPi == 0, u, bcs,
                  form_compiler_parameters=ffc_options)
            u0.assign(u)
            # Piola-Kirchhoff Stress
            stress = diff(Psi, F)
            stress_calculation = (project(stress,  W_tensor))
            file2 << stress_calculation
            t = dot(stress, normal)
            t = assemble(t[1] * dss(1))
        max_error += ((t - f_desired[count]) ** 2)
        print("desired f", f_desired[count])
        print("calculated f", t)
        count += 1
    J = c1 * max_error + assemble((((1/c2) * (chi * (1 - chi))  ** 2) + 0.5 * c2 * (inner(grad(chi), grad(chi)))) * dx)
    print("max_error", max_error)
    print("J", J)

    controls = File("./Results" + "/control_iterations_guess.pvd")
    chi_viz = Function(W, name="ControlVisualisation")

    def eval_cb(j, chi):
        chi_viz.assign(chi)
        controls << chi_viz
    control = Control(chi)
    rf = ReducedFunctional(J, control, eval_cb_post=eval_cb)
    problem = MinimizationProblem(rf)  # , constraints=const
    parameters_ipopt = {"acceptable_tol": 1.0e-10, "maximum_iterations": 100}
    solver = IPOPTSolver(problem, parameters=parameters_ipopt)  # , parameters=parameters_ipopt
    f_opt = solver.solve()

    file = File(
        "./Results" + "/Xnew.pvd");
    file << f_opt

    plot(f_opt, title="f_opt")
    plt.savefig('./Results' + '/Xnew' + '--' + str(i) + '.png')

    print(f"Initial J: {rf(chi)}")
    print(f"min opt {min(f_opt.vector().get_local())}, max opt {max(f_opt.vector().get_local())}")
    print(f"Optimal J {rf(f_opt)}")

    return f_opt

opt_chi = optimization(chi)
print("Optimization Done")


def forward(opt_chi, Traclist, u0):
    u = Function(V_vector)
    Ic, J_, v, F = kinematics(u)
    # Stored strain energy density (compressible neo-Hookean model)
    Psi = (H((opt_chi - 0.5), l2)) * psi(mu1, Ic, J_, lmbda1) + (1 - H((opt_chi - 0.5), l2)) * psi(mu2, Ic, J_, lmbda2)
    # Total potential energy
    Pi = Psi * dx
    # Compute first variation of Pi (directional derivative about u in the direction of v)
    dPi = derivative(Pi, u, v)
    # Piola-Kirchhoff Stress
    stress = diff(Psi, F)
    # Solve variational problem
    disp = Constant([0.0, 0.0])
    bcs, dss = bcs_(disp)
    recu = []
    for d in displacement_by_step:
        print(d)
        disp.assign(Constant([0.0, d]))
        u.assign(u0)
        solve(dPi == 0, u, bcs)
        u0.assign(u)
        recu.append(u(0, 1))
        file1 << u
        stress_calculation = (project((stress), W_tensor))
        traction = dot(stress, normal)
        storedtraction = assemble((traction[1]) * dss(1))
        Traclist.append(storedtraction)
    return store_u, Traclist

Traclist = []
u0 = Function(V_vector)
store_u, Traclist = forward(H((opt_chi - 0.5), l2), Traclist, u0)
Traclist2 = []
u0 = Function(V_vector)
store_u2, Traclist2 = forward(H((chi - 0.5), l2), Traclist2, u0)

file = File(
    "./Results" + "/chi_last.pvd");
file << project(H((opt_chi - 0.5), l2), W)


# Plot
K = -1
plt.figure()
displacement_by_step = [x * K for x in displacement_by_step]
Traclist = [y * K for y in Traclist]
plt.plot(displacement_by_step, Traclist, label="After- opt")
Traclist2 = [y * K for y in Traclist2]
plt.plot(displacement_by_step, Traclist2, label="Before- opt")
max_disp2 = [x * K for x in max_disp]
f_desired2 = [y * K for y in f_desired]
plt.plot(max_disp2, f_desired2, 'o', label="Desired Points")
plt.legend()
plt.xlabel("Displacement")
plt.ylabel("Force")
plt.title("Graph for the displacement" + str(max(max_disp)))
plt.savefig('./Results' + '/Graph.png')

I use a different chi expression from an image.

Hi, do I need to minimize my code? Or I can provide more information @dokken