使用 scikits.odes 和日晷求解 DAE 时如何解决“IDASolve：校正器收敛失败”？

Question

我正在尝试使用 Sundials (https://computation.llnl.gov/projects/sundials/ida), through the Python package scikits.odes (https://scikits-odes.readthedocs.io) 的求解器 IDA 求解 DAE 系统（2 个 ODE 和 1 个代数方程）。

我正在使用 scikits.odes 2.4.0、日晷 3.1.1 和 Python 3.6 64 位。

这是代码：

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)

当我运行它时，出现以下错误：

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

你知道它可能来自哪里吗？

Answer 1

直接原因应该是初始向量不一致，违反了代数部分

0 == m1 @ y3 + v1

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

除此之外，据我所知，您的系统有索引 2，这需要专门的 DAE 求解器。 Sundials/IDA中使用的DASPK，一般只求解index-1 DAE。根据版本的不同，index-2 DAE 的某些特殊类也可以解决。从 R wrapper 文档中可以了解到，如果已知变量的最大微分阶数，则可以获得索引最多为 3 的解。我不知道这个 python 包装器是否为此准备。

没有 DAE 求解器的解决方案通过手动索引约简

系统

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

可以转化为 ODE 内核和状态重建程序。 ODE 的简化状态由分量 [ b, c3 ] 和方程

组成

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

并用于状态重建（从 b,c3,c3' 开始，使用 ODE 函数求导数）

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

如果想在 ODE 中包含完整状态，则需要再次微分所有方程（可能第三个方程除外）（因此微分指数为 2）。然后通过投影从ODE系统（包含一些导数组件）的状态获得DAE的状态。

Answer 2

LutzL 关于初始条件的说法是正确的。还要感谢 SUNDIALS 论坛 (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) 中 A.C.Hindmarsh 的回答，我发现求解器 IDA 的包装器可以采用一个选项 compute_initcond 来指定我们想要的初始条件求解器自行计算。因此，我将此选项设置为 'y0'，求解器现在成功集成了我的系统。

使用 scikits.odes 和日晷求解 DAE 时如何解决“IDASolve：校正器收敛失败”？

How to fix ‘IDASolve : the corrector convergence failed’ when solving a DAE with scikits.odes and Sundial?

python

ode

scikits

sundials

没有 DAE 求解器的解决方案通过手动索引约简