Using Ipopt/pyipopt to solve constrained problem (NOT adjoint)?

dolfin-adjoint
1

I have a “tricky” mixed problem that I can’t (easily) solve using PETSc/Tao/etc. I’m minimizing an energy functional, \Pi(\mathbf{q}), \mathbf{q}:=[\mathbf{u | a | b | c}] subject to some dirichlet boundary conditions on the displacement field \mathbf{u}, and additional inequality and bound constraints on the extra degrees of freedom a,b,c:

min \Pi(\mathbf{q}), \mathbf{q}:=[\mathbf{u | a | b | c}] \\ s.t. \\ u = u_D \forall x \in \Gamma_D \\a_i >0 \\ b_i >0 \\c_i >0 \\ c_i < \pi,\\ a_i >b_i \\

I got ok results with a DIY primal-dual interior point method (with dynamic relaxation), but I need something more robust.

I was wondering if there is an easy way to hack the IPOPT interface provided by dolfin adjoint. It seems like it should work, but I am not sure how to pass in my boundary conditions and bounds. Inspecting the pyipopt interface to IPOPT https://github.com/live-clones/dolfin-adjoint/blob/master/dolfin_adjoint/optimization/ipopt_solver.py, it seems like the bounds correspond to lb and ub. I’m not sure what nconstraints is (does it include bounds and constraints or just constraints?)…

For the modified hyperelasticity demo, how can I impose the Dirichlet boundary conditions /interface to IPOPT? I appreciate any tips. Thanks

from dolfin import *
from dolfin_adjoint import *

mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0", "0.0"), degree=1)
r = Expression(("scale*0.0",
                "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
                "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
                scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3, degree=1)

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]

du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, -0.5, 0.0))  # Body force per unit volume
T  = Constant((0.1,  0.0, 0.0))  # Traction force on the boundary

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

Ic = tr(C)
J  = det(F)

E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
F = derivative(Pi, u, v)
J = derivative(F, u, du)
# Solve using newton
solve(F == 0, u, bcs, J=J)  #works
file = File("newton.pvd");
file << u  

# solve using ipopt
Jhat = ReducedFunctional(assemble(Pi), Control(u))
problem = MinimizationProblem(Jhat,
                              bounds=None,  # default to +- inf
                              constraints=None  # bcs go here???
                              )
solver = IPOPTSolver(problem, parameters={"maximum_iterations": 50,
                     "acceptable_tol": 1.0e-8}) 
u = solver.solve()  # obviously does not work
file = File("ipopt.pvd");
file << u
2

As the IPOPT solver is based on CYIPOPT, you should check its documentation: https://pythonhosted.org/ipopt/tutorial.html#executing-the-solver