Bug in CVODES or its interface

[zegkljan]

I was interested in the Gravity Turn Maneuver (found here). However, it is written for CasADi 2.4.2. I tried to modify it to work with 3.1.1 (see the code at the very end of this post). However, when I run it (in python 3.5.2), the solver seems to be stuck with one CPU core working at 100%. When I set the print_level option of IPOPT to 6, I get the following output:

> python3 gturn3.py 

List of options:

                                    Name   Value                # times used
                           linear_solver = mumps                     2
                             print_level = 6                         2
                                     tol = 1e-05                     2

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

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

Scaling parameter for objective function = 1.000000e+00
objective scaling factor = 1
No x scaling provided
No c scaling provided
No d scaling provided
Moved initial values of x sufficiently inside the bounds.
Initial values of s sufficiently inside the bounds.
MUMPS used permuting_scaling 5 and pivot_order 0.
           scaling will be 77.
Number of doubles for MUMPS to hold factorization (INFO(9)) = 37484
Number of integers for MUMPS to hold factorization (INFO(10)) = 28536
Factorization successful.
Least square estimates max(y_c) = 2.070943e+03, max(y_d) = 0.000000e+00
Total number of variables............................:     1799
                     variables with only lower bounds:      898
                variables with lower and upper bounds:      901
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1500
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


**************************************************
*** Update HessianMatrix for Iteration 0:
**************************************************

and then I'm waiting forever.

However, using CasADi 2.4.2, it works fine and I get progress reports right away and after about 80 seconds it finishes and I get the results. Could this be me failing to properly convert the code to work with 3.1.1 or is it a bug? This is my very first time touching CasADi (in fact I came across the gravity turn rather than CasADi and I just wanted to get it working with the newest versions of everything...) so I won't be surprised if it is me. However, right now I have no idea. I can post additional info/data if instructed.

Contents of my gturn3.py file:

# ----------------------------------------------------------------
# Gravity Turn Maneuver with direct multiple shooting using CVodes
# (c) Mirko Hahn
# ----------------------------------------------------------------
import sys
sys.path.append('./casadi-3.1.1')

from casadi import *
 
# Artificial model parameters
N = 300                   # Number of shooting intervals
vel_eps = 1e-6            # Initial velocity (km/s)
 
# Vehicle parameters
m0   = 11.3               # Launch mass (t)
m1   = 1.3                # Dry mass (t)
g0   = 9.81e-3            # Gravitational acceleration at altitude zero (km/s^2)
r0   = 6.0e2              # Radius at altitude zero (km)
Isp  = 300.0              # Specific impulse (s)
Fmax = 600.0e-3           # Maximum thrust (MN)
 
# Atmospheric parameters
cd  = 0.021                        # Drag coefficients
A   = 1.0                          # Reference area (m^2)
H   = 5.6                          # Scale height (km)
rho = (1.0 * 1.2230948554874)      # Density at altitude zero
 
# Target orbit parameters
h_obj = 75                # Target altitude (km)
v_obj = 2.287             # Target velocity (km/s)
q_obj = 0.5 * pi          # Target angle to vertical (rad)

# Create symbolic variables
x = SX.sym('[m, v, q, h, d]')      # Vehicle state
u = SX.sym('u')                    # Vehicle controls
T = SX.sym('T')                    # Time horizon (s)
 
# Introduce symbolic expressions for important composite terms
Fthrust = Fmax * u
Fdrag   = 0.5e3 * A * cd * rho * exp(-x[3] / H) * x[1]**2
r       = x[3] + r0
g       = g0 * (r0 / r)**2
vhor    = x[1] * sin(x[2])
vver    = x[1] * cos(x[2])
 
# Build symbolic expressions for ODE right hand side
mdot = -(Fmax / (Isp * g0)) * u
vdot = (Fthrust - Fdrag) / x[0] - g * cos(x[2])
hdot = vver
ddot = vhor / r
qdot = g * sin(x[2]) / x[1] - ddot
 
# Build the DAE function
ode = [
    mdot,
    vdot,
    qdot,
    hdot,
    ddot
]
quad = u
dae = {'x': x, 'p': vertcat(u, T), 'ode': T * vertcat(*ode), 'quad': T * quad}
I = integrator("I", "cvodes", dae, {'t0': 0.0, 'tf': 1.0 / N})
 
# Specify upper and lower bounds as well as initial values for DAE parameters,
# states and controls
p_min  = [120.0]
p_max  = [600.0]
p_init = [120.0]
 
u_min  = [0.0]
u_max  = [1.0]
u_init = [0.5]
 
x0_min  = [m0, vel_eps,       0.0, 0.0, 0.0]
x0_max  = [m0, vel_eps,  0.5 * pi, 0.0, 0.0]
x0_init = [m0, vel_eps, 0.05 * pi, 0.0, 0.0]
 
xf_min  = [m1, v_obj, q_obj, h_obj, 0.0]
xf_max  = [m0, v_obj, q_obj, h_obj, inf]
xf_init = [m1, v_obj, q_obj, h_obj, 0.0]
 
x_min  = [m1, vel_eps, 0.0, 0.0, 0.0]
x_max  = [m0,     inf,  pi, inf, inf]
x_init = [0.5 * (m0 + m1), 0.5 * v_obj, 0.5 * q_obj, 0.5 * h_obj, 0.0]
 
# Useful variable block sizes
np = 1                         # Number of parameters
nx = x.size1()                 # Number of states
nu = u.size1()                 # Number of controls
ns = nx + nu                   # Number of variables per shooting interval
 
# Introduce symbolic variables and disassemble them into blocks
V = MX.sym('X', N * ns + nx + np)
P = V[0]
X = [V[np+i*ns:np+i*ns+nx] for i in range(0, N+1)]
U = [V[np+i*ns+nx:np+(i+1)*ns] for i in range(0, N)]
 
# Nonlinear constraints and Lagrange objective
G = []
F = 0.0
 
# Build DMS structure
x0 = p_init + x0_init
for i in range(0, N):
    Y = I(x0=X[i], p=vertcat(U[i], P))
    G = G + [Y['xf'] - X[i+1]]
    F = F + Y['qf']
 
    frac = float(i+1) / N
    x0 = x0 + u_init + [x0_init[i] + frac * (xf_init[i] - x0_init[i]) for i in range(0, nx)]
 
# Lower and upper bounds for solver
lbg = 0.0
ubg = 0.0
lbx = p_min + x0_min + u_min + (N-1) * (x_min + u_min) + xf_min
ubx = p_max + x0_max + u_max + (N-1) * (x_max + u_max) + xf_max
 
# Solve the problem using IPOPT
nlp = {'x': V, 'f': m0 - X[-1][0], 'g': vertcat(*G)}
S = nlpsol("S", "ipopt", nlp, {'ipopt': {'tol': 1e-5, 'print_level': 6, 'linear_solver': 'mumps'}})
r = S(x0=x0,
      lbx=lbx,
      ubx=ubx,
      lbg=lbg,
      ubg=ubg)
 
# Extract state sequences and parameters from result
x = r['x']
f = r['f']
T = x[0]
 
t = [i * (T / N) for i in range(0, N+1)]
m = x[np  ::ns].get()
v = x[np+1::ns].get()
q = x[np+2::ns].get()
h = x[np+3::ns].get()
d = x[np+4::ns].get()
u = x[np+nx::ns].get() + [0.0]

[jgillis]

Nice catch!
it's definitely number-crunching, but seems to be stuck in the Newton root-finder.
passing "nonlinear_solver_iteration": "functional" to the options seems to resolve the situation, while "fsens_err_con": True, "quad_err_con": True makes the convergence a bit faster.

@jaeandersson interested to have a look?

[jaeandersson]

Can you try replacing "cvodes" with "idas"?

[zegkljan]

@jgillis I can confirm it works for me as well.
@jaeandersson That works too.

Now I'm getting errors regarding the get() calls at the end of my script, which is just 2.4.2 vs 3.1.1 API thing, so it's back to me. Thanks for your help!

[jaeandersson]

@jaeandersson That works too.

So it works with idas but not with cvodes? I've seen that before.

[zegkljan]

@jaeandersson Sorry for the brevity. Exactly as you say, it works with idas but not with cvodes.

[jaeandersson]

Looks like a bug in CVODES then, or in CasADi's interface to CVODES. I'd use IDAS instead until we figure out what it is. CVODES and IDAS are both from the SUNDIALS suite and use the same variable step-size, variable-order BDF method by default.

[jaeandersson]

I think I've fixed it now... the test suite gets the following error now though:

======================================================================
FAIL: test_hess6 (integration.Integrationtests)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/jaeandersson/dev/casadi/test/python/integration.py", line 737, in test_hess6
    self.assertAlmostEqual(H_out[0][0],(q0*tend**6*exp(tend**3/(3*p)))/(9*p**4)+(2*q0*tend**3*exp(tend**3/(3*p)))/(3*p**3),9,"Evaluation output mismatch")
  File "/Users/jaeandersson/dev/casadi/test/python/helpers.py", line 224, in assertAlmostEqual
    unittest.TestCase.assertAlmostEqual(self,first/n,second/n,places=places,msg=msg + "  scaled by %d" %n)
AssertionError: Evaluation output mismatch 1.1014858783948605e+02 <-> 1.1014858783785127e+02  scaled by 1

----------------------------------------------------------------------

False positive, @jgillis ?

[jaeandersson]

I changed significant digits to 8...