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?

4

2 回答 2

0

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

0 == m1 @ y3 + v1

asy1=[0,0]并且v1=[0.3, 9]*1e-4是非零的。

除此之外，据我所知，您的系统具有索引 2，这需要专门的 DAE 求解器。Sundials/IDA 中使用的 DASPK，一般只解决 index-1 DAE。根据版本，某些特殊类别的 index-2 DAE 也可以得到解决。从 R 包装器文档中可以了解到，如果变量的最大微分阶数已知，则可以获得最多为 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 的状态。

于 2019-06-28T13:51:05.527 回答
0

LutzL 关于初始条件是正确的。还要感谢 SUNDIALS 论坛中 ACHindmarsh 的回答（http://sundials.2283335.n4.nabble.com/IDA-ERROR-IDASolve-the-corrector-convergence-failed-repeatedly-or-with-h-hmin- td4655649.html），我更深入地查看了 Scikits.Odes 的文档（https://github.com/bmcage/odes/blob/master/scikits/odes/sundials/ida.pyx），我发现求解器的包装器 IDA 可以采用一个选项compute_initcond来指定我们希望求解器自行计算的初始条件。因此，我将此选项设置为'y0'，求解器现在成功集成了我的系统。

于 2019-07-01T14:49:46.090 回答