optimization problem by use IPOPTSolver,the range of guessed value is very close to the target value.

[jingbinzhao415]

Hi,

I am working on optimization of material parameters. The objective is to minimize the displacement or (volume) difference between experimental measurements and the output of the finite element simulation to obtain the hyperelastic material parameters which are the similar as my target values.

My problem is: my target values of material parameter are a =1.291 and b=9.726, the result can be converged when the ranges of given guessed values are very close to the target values. (for example, the range of a can be accept is from 0.5 to 1.5, if larger than 1.5, like 3, the solver will be crashed.)
Everything works fine until to the optimization part, I am not sure if my forward problem has problem or the setting of parameters in optimization part is wrong.
I will be very grateful for all your suggestions and any help!!

below is my small benchmark code.

`

import dolfin as df
import dolfin_adjoint as da
import json
try:
    from dolfin_adjoint import (
        Constant,
        DirichletBC,
        Expression,
        UnitCubeMesh,
        interpolate,
        Mesh,
    )
except ImportError:
    from dolfin import (
        UnitCubeMesh,
        Expression,
        Constant,
        DirichletBC,
        interpolate,
        Mesh,
    )

import pulse

from collections import deque
import time
pulse.iterate.logger.setLevel(10)

# Create mesh
N = 6
mesh = da.UnitCubeMesh(N, N, N)
# Create subdomains
class Free(df.SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > (1.0 - df.DOLFIN_EPS) and on_boundary


class Fixed(df.SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < df.DOLFIN_EPS and on_boundary

def foward(mesh,param_a,param_b):
    # Create a facet fuction in order to mark the subdomains
    ffun = df.MeshFunction("size_t", mesh, 2)
    ffun.set_all(0)

    # Mark the first subdomain with value 1
    fixed = Fixed()
    fixed_marker = 1
    fixed.mark(ffun, fixed_marker)

    # Mark the second subdomain with value 2
    free = Free()
    free_marker = 2
    free.mark(ffun, free_marker)

    # Create a cell function (but we are not using it)
    cfun = df.MeshFunction("size_t", mesh, 3)
    cfun.set_all(0)


    # Collect the functions containing the markers
    marker_functions = pulse.MarkerFunctions(ffun=ffun, cfun=cfun)
    # Create mictrotructure
    V_f = df.VectorFunctionSpace(mesh, "CG", 1)

    # Fibers
    f0 = interpolate(Expression(("1.0", "0.0", "0.0"), degree=1), V_f)
    # Sheets
    s0 = interpolate(Expression(("0.0", "1.0", "0.0"), degree=1), V_f)
    # Fiber-sheet normal
    n0 = interpolate(Expression(("0.0", "0.0", "1.0"), degree=1), V_f)


    # Collect the mictrotructure
    microstructure = pulse.Microstructure(f0=f0, s0=s0, n0=n0)
    # Create the geometry
    geometry = pulse.Geometry(
        mesh=mesh,
        marker_functions=marker_functions,
        microstructure=microstructure,
    )

    material_parameters = dict(
        a=param_a,
        a_f=5.0,
        b=param_b,
        b_f=5.0,
        a_s=0.0,
        b_s=0.0,
        a_fs=0.0,
        b_fs=0.0,
    )


    # Select model for active contraction
    active_model = pulse.ActiveModels.active_strain
    # active_model = "active_stress"

    # Set the activation
    activation = Constant(0.0)

    # Create material

    material = pulse.HolzapfelOgden(
        active_model=active_model,
        parameters=material_parameters,
        f0=geometry.f0,
        s0=geometry.s0,
        n0=geometry.n0,
        activation=activation,
        T_ref=1.0,
        )

    # Make Dirichlet boundary conditions
    def dirichlet_bc(W):
        V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
        return DirichletBC(V, Constant((0.0, 0.0, 0.0)), fixed)

    lvp = da.Constant(0.0)
    # Make Neumann boundary conditions
    neumann_bc = pulse.NeumannBC(traction=lvp, marker=free_marker)

    # Collect Boundary Conditions
    bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))

    # Create problem
    problem = pulse.MechanicsProblem(geometry, material, bcs)


    pulse.iterate.iterate(problem, lvp, da.Constant(2.0),max_nr_crash=100, max_iters=100)
    u, p = problem.state.split(deepcopy=True)
    return u

parameter_a = da.Constant(1.291)
parameter_b = da.Constant(9.726)
guessed_parameter_a = da.Constant(3.0)
guessed_parameter_b = da.Constant(15.0)
u_syn = foward(mesh,parameter_a,parameter_b)
V = df.VectorFunctionSpace(mesh, "CG", 2)
u_synthetic = da.project(u_syn, V)
u_crl = foward(mesh,guessed_parameter_a,guessed_parameter_b)
def cost_function(u_model, u_data):
    norm = lambda f: da.assemble(df.inner(f, f) * df.dx)
    return norm(u_model - u_data)

J = cost_function(u_crl,u_synthetic)
print(J)

controls = [da.Control(guessed_parameter_a), da.Control(guessed_parameter_b)]
cost_func_values = []
control_values = []
def eval_cb(j, m):
    """Callback function"""
    cost_func_values.append(j)
    control_values.append(m)

reduced_functional = da.ReducedFunctional(
    J,
    controls,
   eval_cb_post=eval_cb,
)


problem = da.MinimizationProblem(reduced_functional, bounds=[(0, 4.0),(0, 16.0)])
parameters = {
    #"tolerance": 1e-30,
    "limited_memory_initialization": "scalar2",
    "maximum_iterations": 20,
}
solver = da.IPOPTSolver(problem, parameters=parameters)

optimal_control = solver.solve()
for i in optimal_control:
    print(i.values())

`
error message:

List of options:

                                Name   Value                # times used
               hessian_approximation = limited-memory            7
       limited_memory_initialization = scalar2                   1
                            max_iter = 20                        1
                         print_level = 6                         2

This is Ipopt version 3.14.4, running with linear solver MUMPS 5.2.1.

Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0

Hessian approximation will be done in the space of all 2 x variables.

List of options:

                                Name   Value                # times used
               hessian_approximation = limited-memory            7
       limited_memory_initialization = scalar2                   1
                            max_iter = 20                        1
                         print_level = 6                         2

This is Ipopt version 3.14.4, running with linear solver MUMPS 5.2.1.

Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0

Hessian approximation will be done in the space of all 2 x variables.

Scaling parameter for objective function = 1.000000e+00
objective scaling factor = 1
No x scaling provided
No c scaling provided
No d scaling provided
Initial values of x sufficiently inside the bounds.
Initial values of s sufficiently inside the bounds.
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 2
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0

---

*** Summary of Iteration: 0:

---

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 0.00e+00 6.55e-04 0.0 0.00e+00 - 0.00e+00 0.00e+00 0

---

*** Beginning Iteration 0 from the following point:

---

Current barrier parameter mu = 1.0000000000000000e+00
Current fraction-to-the-boundary parameter tau = 0.0000000000000000e+00

||curr_x||_inf = 1.5000000000000000e+01
||curr_s||_inf = 0.0000000000000000e+00
||curr_y_c||_inf = 0.0000000000000000e+00
||curr_y_d||_inf = 0.0000000000000000e+00
||curr_z_L||_inf = 1.0000000000000000e+00
||curr_z_U||_inf = 1.0000000000000000e+00
||curr_v_L||_inf = 0.0000000000000000e+00
||curr_v_U||_inf = 0.0000000000000000e+00

***Current NLP Values for Iteration 0:

                               (scaled)                 (unscaled)

Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 6.5483061078119853e-04 6.5483061078119853e-04
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 1.5000000010000001e+01 1.5000000010000001e+01
Overall NLP error.......: 1.5000000010000001e+01 1.5000000010000001e+01

Number of Iterations....: 0

                               (scaled)                 (unscaled)

Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 6.5483061078119853e-04 6.5483061078119853e-04
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 1.5000000010000001e+01 1.5000000010000001e+01
Overall NLP error.......: 1.5000000010000001e+01 1.5000000010000001e+01

Number of objective function evaluations = 1
Number of objective gradient evaluations = 1
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 0
Total seconds in IPOPT = 350.550

EXIT: Stopping optimization at current point as requested by user.

RuntimeError Traceback (most recent call last)
Untitled-1.ipynb Cell 5 in <cell line: 24>()
17 parameters = {
18 #"tolerance": 1e-30,
19 "limited_memory_initialization": "scalar2",
20 "maximum_iterations": 20,
21 }
22 solver = da.IPOPTSolver(problem, parameters=parameters)
---> 24 optimal_control = solver.solve()
25 for i in optimal_control:
26 print(i.values())

File ~/.local/lib/python3.9/site-packages/pyadjoint/optimization/ipopt_solver.py:199, in IPOPTSolver.solve(self)
197 """Solve the optimization problem and return the optimized controls."""
198 guess = self.rfn.get_controls()
--> 199 results = self.ipopt_problem.solve(guess)
200 new_params = [control.copy_data() for control in self.rfn.controls]
201 self.rfn.set_local(new_params, results[0])

File ~/anaconda3/envs/fenicsjoint/lib/python3.9/site-packages/cyipopt/cython/ipopt_wrapper.pyx:642, in ipopt_wrapper.Problem.solve()

File ~/anaconda3/envs/fenicsjoint/lib/python3.9/site-packages/cyipopt/cython/ipopt_wrapper.pyx:676, in ipopt_wrapper.objective_cb()

File ~/.local/lib/python3.9/site-packages/pyadjoint/reduced_functional_numpy.py:36, in ReducedFunctionalNumPy.call(self, m_array)
...

---

*** DOLFIN version: 2019.1.0
*** Git changeset: 7f46aeb0b296da5bbb1fb0845822a72ab9b09c55

{
"name": "RuntimeError",
"message": "\n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n*** [email protected]\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a minimal running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error: Unable to solve nonlinear system with NewtonSolver.\n*** Reason: Newton solver did not converge because maximum number of iterations reached.\n*** Where: This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2019.1.0\n*** Git changeset: 7f46aeb0b296da5bbb1fb0845822a72ab9b09c55\n*** -------------------------------------------------------------------------\n",
"stack": "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)\n\u001b[1;32mUntitled-1.ipynb Cell 5\u001b[0m in \u001b[0;36m<cell line: 24>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 17\u001b[0m parameters \u001b[39m=\u001b[39m {\n\u001b[1;32m 18\u001b[0m \u001b[39m#"tolerance": 1e-30,\u001b[39;00m\n\u001b[1;32m 19\u001b[0m \u001b[39m"\u001b[39m\u001b[39mlimited_memory_initialization\u001b[39m\u001b[39m"\u001b[39m: \u001b[39m"\u001b[39m\u001b[39mscalar2\u001b[39m\u001b[39m"\u001b[39m,\n\u001b[1;32m 20\u001b[0m \u001b[39m"\u001b[39m\u001b[39mmaximum_iterations\u001b[39m\u001b[39m"\u001b[39m: \u001b[39m20\u001b[39m,\n\u001b[1;32m 21\u001b[0m }\n\u001b[1;32m 22\u001b[0m solver \u001b[39m=\u001b[39m da\u001b[39m.\u001b[39mIPOPTSolver(problem, parameters\u001b[39m=\u001b[39mparameters)\n\u001b[0;32m---> 24\u001b[0m optimal_control \u001b[39m=\u001b[39m solver\u001b[39m.\u001b[39;49msolve()\n\u001b[1;32m 25\u001b[0m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m optimal_control:\n\u001b[1;32m 26\u001b[0m \u001b[39mprint\u001b[39m(i\u001b[39m.\u001b[39mvalues())\n\nFile \u001b[0;32m~/.local/lib/python3.9/site-packages/pyadjoint/optimization/ipopt_solver.py:199\u001b[0m, in \u001b[0;36mIPOPTSolver.solve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 197\u001b[0m \u001b[39m"""Solve the optimization problem and return the optimized controls."""\u001b[39;00m\n\u001b[1;32m 198\u001b[0m guess \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrfn\u001b[39m.\u001b[39mget_controls()\n\u001b[0;32m--> 199\u001b[0m results \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mipopt_problem\u001b[39m.\u001b[39;49msolve(guess)\n\u001b[1;32m 200\u001b[0m new_params \u001b[39m=\u001b[39m [control\u001b[39m.\u001b[39mcopy_data() \u001b[39mfor\u001b[39;00m control \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrfn\u001b[39m.\u001b[39mcontrols]\n\u001b[1;32m 201\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrfn\u001b[39m.\u001b[39mset_local(new_params, results[\u001b[39m0\u001b[39m])\n\nFile \u001b[0;32m~/anaconda3/envs/fenicsjoint/lib/python3.9/site-packages/cyipopt/cython/ipopt_wrapper.pyx:642\u001b[0m, in \u001b[0;36mipopt_wrapper.Problem.solve\u001b[0;34m()\u001b[0m\n\nFile \u001b[0;32m~/anaconda3/envs/fenicsjoint/lib/python3.9/site-packages/cyipopt/cython/ipopt_wrapper.pyx:676\u001b[0m, in \u001b[0;36mipopt_wrapper.objective_cb\u001b[0;34m()\u001b[0m\n\nFile \u001b[0;32m~/.local/lib/python3.9/site-packages/pyadjoint/reduced_functional_numpy.py:36\u001b[0m, in \u001b[0;36mReducedFunctionalNumPy.call\u001b[0;34m(self, m_array)\u001b[0m\n\u001b[1;32m 31\u001b[0m \u001b[39m"""An implementation of the reduced functional evaluation\u001b[39;00m\n\u001b[1;32m 32\u001b[0m \u001b[39m that accepts the control values as an array of scalars\u001b[39;00m\n\u001b[1;32m 33\u001b[0m \n\u001b[1;32m 34\u001b[0m \u001b[39m"""\u001b[39;00m\n\u001b[1;32m 35\u001b[0m m_copies \u001b[39m=\u001b[39m [control\u001b[39m.\u001b[39mcopy_data() \u001b[39mfor\u001b[39;00m control \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcontrols]\n\u001b[0;32m---> 36\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mrf\u001b[39m.\u001b[39;49m\u001b[39m__call__\u001b[39;49m(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mset_local(m_copies, m_array))\n\nFile \u001b[0;32m~/.local/lib/python3.9/site-packages/pyadjoint/tape.py:69\u001b[0m, in \u001b[0;36mno_annotations..wrapper\u001b[0;34m(args, kwargs)\u001b[0m\n\u001b[1;32m 66\u001b[0m \u001b[39m@wraps\u001b[39m(function)\n\u001b[1;32m 67\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mwrapper\u001b[39m(\u001b[39m\u001b[39margs, \u001b[39m\u001b[39m\u001b[39m\u001b[39mkwargs):\n\u001b[1;32m 68\u001b[0m \u001b[39mwith\u001b[39;00m stop_annotating():\n\u001b[0;32m---> 69\u001b[0m \u001b[39mreturn\u001b[39;00m function(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\nFile \u001b[0;32m~/.local/lib/python3.9/site-packages/pyadjoint/reduced_functional.py:140\u001b[0m, in \u001b[0;36mReducedFunctional.call\u001b[0;34m(self, values)\u001b[0m\n\u001b[1;32m 136\u001b[0m \u001b[39mwith\u001b[39;00m stop_annotating():\n\u001b[1;32m 137\u001b[0m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtape\u001b[39m.\u001b[39m_bar(\u001b[39m"\u001b[39m\u001b[39mEvaluating functional\u001b[39m\u001b[39m"\u001b[39m)\u001b[39m.\u001b[39miter(\n\u001b[1;32m 138\u001b[0m \u001b[39mrange\u001b[39m(\u001b[39mlen\u001b[39m(blocks))\n\u001b[1;32m 139\u001b[0m ):\n\u001b[0;32m--> 140\u001b[0m blocks[i]\u001b[39m.\u001b[39;49mrecompute()\n\u001b[1;32m 142\u001b[0m func_value \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mscale \u001b[39m*\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfunctional\u001b[39m.\u001b[39mblock_variable\u001b[39m.\u001b[39mcheckpoint\n\u001b[1;32m 144\u001b[0m \u001b[39m# Call callback\u001b[39;00m\n\nFile \u001b[0;32m~/.local/lib/python3.9/site-packages/pyadjoint/block.py:350\u001b[0m, in \u001b[0;36mBlock.recompute\u001b[0;34m(self, markings)\u001b[0m\n\u001b[1;32m 347\u001b[0m prepared \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mprepare_recompute_component(inputs, relevant_outputs)\n\u001b[1;32m 349\u001b[0m \u001b[39mfor\u001b[39;00m idx, out \u001b[39min\u001b[39;00m relevant_outputs:\n\u001b[0;32m--> 350\u001b[0m output \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mrecompute_component(inputs,\n\u001b[1;32m 351\u001b[0m out,\n\u001b[1;32m 352\u001b[0m idx,\n\u001b[1;32m 353\u001b[0m prepared)\n\u001b[1;32m 354\u001b[0m \u001b[39mif\u001b[39;00m output \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m 355\u001b[0m out\u001b[39m.\u001b[39mcheckpoint \u001b[39m=\u001b[39m output\n\nFile \u001b[0;32m~/.local/lib/python3.9/site-packages/dolfin_adjoint_common/blocks/solving.py:469\u001b[0m, in \u001b[0;36mGenericSolveBlock.recompute_component\u001b[0;34m(self, inputs, block_variable, idx, prepared)\u001b[0m\n\u001b[1;32m 467\u001b[0m func \u001b[39m=\u001b[39m prepared[\u001b[39m2\u001b[39m]\n\u001b[1;32m 468\u001b[0m bcs \u001b[39m=\u001b[39m prepared[\u001b[39m3\u001b[39m]\n\u001b[0;32m--> 469\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_forward_solve(lhs, rhs, func, bcs)\n\nFile \u001b[0;32m~/.local/lib/python3.9/site-packages/fenics_adjoint/blocks/variational_solver.py:81\u001b[0m, in \u001b[0;36mNonlinearVariationalSolveBlock._forward_solve\u001b[0;34m(self, lhs, rhs, func, bcs, kwargs)\u001b[0m\n\u001b[1;32m 79\u001b[0m solver \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mbackend\u001b[39m.\u001b[39mNonlinearVariationalSolver(problem, \u001b[39m\u001b[39m\u001b[39mself\u001b[39m\u001b[39m.\u001b[39msolver_args, \u001b[39m\u001b[39m\u001b[39m*\u001b[39m\u001b[39mself\u001b[39m\u001b[39m.\u001b[39msolver_kwargs)\n\u001b[1;32m 80\u001b[0m solver\u001b[39m.\u001b[39mparameters\u001b[39m.\u001b[39mupdate(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39msolver_params)\n\u001b[0;32m---> 81\u001b[0m solver\u001b[39m.\u001b[39;49msolve(\u001b[39m*\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49msolve_args, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49msolve_kwargs)\n\u001b[1;32m 82\u001b[0m \u001b[39mreturn\u001b[39;00m func\n\n\u001b[0;31mRuntimeError\u001b[0m: \n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n*** [email protected]\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a minimal running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error: Unable to solve nonlinear system with NewtonSolver.\n*** Reason: Newton solver did not converge because maximum number of iterations reached.\n*** Where: This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2019.1.0\n*** Git changeset: 7f46aeb0b296da5bbb1fb0845822a72ab9b09c55\n*** -------------------------------------------------------------------------\n"
}