0

I'm trying to solve a system of DAE (2 ODE and 1 algebraic equation) using the solver IDA from Sundials (https://computation.llnl.gov/projects/sundials/ida), through the Python package scikits.odes (https://scikits-odes.readthedocs.io).

I'm using scikits.odes 2.4.0, Sundials 3.1.1 and Python 3.6 64bit.

Here is the code :

import numpy as np
from scikits.odes.dae import dae

SOLVER = 'ida'
extra_options = {'old_api': False, 'algebraic_vars_idx': [0, 1]}

##### Parameters #####

# vectors
v1 = np.array([3.e-05, 9.e-04])
v2 = np.array([-0.01])
v3 = np.array([100])

# matrices
m1 = np.array([[-1, 1, -1], [0, -1, 0]])
m2 = np.array([[1, 0, 0]])
m3 = np.array([[0, 0, 1]])
m4 = np.array([[0., 0., 0., 0., 0., 0.],
               [0., 0., 0., 0., 0., 0.],
               [0., 0., 2000., 0., 0., 0.],
               [0., 0., 0., 13e3, 0., 0.],
               [0., 0., 0., 0., 13e3, 0.],
               [0., 0., 0., 0., 0., 13e3]])

##### Equations #####

def f(_, y):
    y1 = y[:2]
    y2 = y[2:3]
    y3 = y[3:]
    return np.hstack((m1 @ y3 + v1,
                      - m2 @ y3 - v2,
                      - 2e2 * np.abs(y3*1000) ** 0.852 * y3 + m1.T @ y1 + m2.T @ y2 + m3.T @ v3))

def right_hand_side(_, y, ydot, residue):

    residue[:] = f(_, y) - m4 @ ydot

    return 0

##### Initial conditions and time grid #####

tout = np.array([0.,  300.])

y_initial = np.array([0., 0., 0., 0., 0., 0.])

ydot_initial = np.array([0., 0., 0., 0., 0., 0.])

##### Solver #####

dae_solver = dae(SOLVER, right_hand_side, **extra_options)
output = dae_solver.solve(tout, y_initial, ydot_initial)
print(output.values.y)

When I run it, I get the following error :

[IDA ERROR]  IDASolve
  At t = 0 and h = 2.86102e-07, the corrector convergence failed repeatedly or with |h| = hmin.

Do you have any idea of from where it could come from?

Camille C
  • 27
  • 7

2 Answers2

0

The immediate cause should be that the initial vector is not a consistent state, as it violates the algebraic part

0 == m1 @ y3 + v1

as y1=[0,0] and v1=[0.3, 9]*1e-4 is non-zero.

Apart from that, as far as I can see, your system has index 2, this requires a specialized DAE solver. DASPK, which is used in Sundials/IDA, generally only solves index-1 DAE. Depending on the version, certain special classes of index-2 DAE can also be solved. From the R wrapper documentation one can learn that if the maximal differentiation orders of the variables are known, that solutions for index up to 3 can be obtained. I do not know if this python wrapper is prepared for that.

Solution without DAE solver via manual index reduction

The system

0 = -c1+c2-c3 + v11
0 =    -c2    + v12

m*b' = -c1 - v2

M*c1' = f(c1) - a1     + b 
M*c2' = f(c2) + a1-a2   
M*c3' = f(c3) - a1     + v3

can be transformed into an ODE kernel and a state reconstruction procedure. The reduced state for the ODE consists of the components [ b, c3 ] with the equations

b'  = -(c3 + v2)/m
c3' = 0.5*(f(c3)-f(v11+v12-c3)+v3-b)/M

and for the state reconstruction (starting from b,c3,c3' using the ODE function for derivatives)

c1 = v11+v12-c3
c2 = v22
a1 = f(c1)+b+M*c3'
a2 = f(c2)+a1

If one wanted to include the full state in the ODE one would need to differentiate all equations (except possibly the third) once more (thus differentiation index 2). Then the state of the DAE is obtained from the state of the ODE system (containing some derivative components) via projection.

Community
  • 1
  • 1
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • your solution without DAE solver is very interesting. It not only transforms the DAE to an ODE (making it easier to solve) ; it also reduces the number of unknowns in the system, and thus reduces the memory needed. Thanks! – Camille C Jul 01 '19 at 15:03
0

LutzL is right about initial conditions. Thanks also to an answer of A.C.Hindmarsh in the SUNDIALS forum (http://sundials.2283335.n4.nabble.com/IDA-ERROR-IDASolve-the-corrector-convergence-failed-repeatedly-or-with-h-hmin-td4655649.html), I took a deeper look in the documentation of Scikits.Odes (https://github.com/bmcage/odes/blob/master/scikits/odes/sundials/ida.pyx), and I found that the wrapper to solver IDA can take an option compute_initcond to specifies the initial condition that we want the solver to compute by itself. Thus I set this option to 'y0' and the solver now succeeds to integrate my system.

Camille C
  • 27
  • 7