Как исправить «IDASolve: сбой корректора» при решении DAE с scikits.odes и Sundial? - PullRequest
Новая
       49

Как исправить «IDASolve: сбой корректора» при решении DAE с scikits.odes и Sundial?

0 голосов
/ 24 июня 2019

Я пытаюсь решить систему DAE (2 ОДУ и 1 алгебраическое уравнение), используя решатель IDA из Солнечных часов (https://computation.llnl.gov/projects/sundials/ida), через пакет Python scikits.odes (https://scikits -одес.readthedocs.io ).

Я использую scikits.odes 2.4.0, Солнечные часы 3.1.1 и Python 3.6 64bit.

Вот код: 

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.

Есть ли у вас какие-либо идеи о том, откуда она может взяться?

Ответы [ 2 ]

0 голосов
/ 01 июля 2019

LutzL прав насчет начальных условий. Благодаря также ответу ACHindmarsh на форуме SUNDIALS (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', и теперь решателю удается интегрировать мою систему.

0 голосов
/ 28 июня 2019

Непосредственной причиной должно быть то, что начальный вектор не является непротиворечивым состоянием, поскольку он нарушает алгебраическую часть 

0 == m1 @ y3 + v1

, поскольку y1=[0,0] и v1=[0.3, 9]*1e-4 ненулевые.

Кроме того, насколько я вижу, ваша система имеет индекс 2, для этого требуется специализированный решатель DAE.DASPK, который используется в Sundials / IDA, обычно решает только DAE индекса-1.В зависимости от версии, определенные специальные классы DAE индекса-2 также могут быть решены.Из документации оболочки 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

Если бы кто-то хотел включить полное состояние в ОДУ, ему нужно было бы еще раз дифференцировать все уравнения (кроме, возможно, третьего) (таким образом, индекс дифференциации 2).Затем состояние DAE получается из состояния системы ODE (содержащей некоторые производные компоненты) посредством проекции.

...