Solving index 2 DAEs: ValueError: array must not contain infs or NaNs

[nathan-t4]

I am trying to use your package to solve the following DAE with index 2, but the solver is running into ValueError: array must not contain infs or NaNs.

Is this error because the solver cannot resolve index 2 DAEs yet? I've checked that the initial conditions are consistent by checking F(t0, y0, yp0) = 0.

import numpy as np
import matplotlib.pyplot as plt
import jax
import jax.numpy as jnp
import time
from functools import partial
from scipy_dae.integrate import solve_dae, consistent_initial_conditions
from scipy.optimize._numdiff import approx_derivative

def F(t, y, yp, u, A, system_data):
    """Define implicit system of differential algebraic equations."""
    R, L, C, splits = system_data
    AC, AR, AL, AI, AV = A

    y = np.nan_to_num(y)
    yp = np.nan_to_num(yp)
    qc, phi, e, jv = np.split(y, splits)
    qp, phip, ep, jvp = np.split(yp, splits)
    i, v = u(t)
    g = lambda e : (AR.T @ e) / R 

    def H(y):
        qc, phi, e, jv = np.split(y, splits)
        return 0.5 * ((L**(-1)) * phi.T @ phi + (C**(-1)) * qc.T @ qc)

    # dH = jax.grad(H)(y)
    # dH = np.split(dH, splits)
    dH1 = phi / L
    dH0 = qc / C

    F0 = AC @ qp + AL @ dH1 + AV @ jv + AR @ g(e) + AI @ i 
    F1 = phip - AL.T @ e
    F2 = AC.T @ e - dH0                   # algebraic equation
    F3 = AV.T @ e - v                     # algebraic equation
    return np.concatenate((F0, F1, F2, F3))

def jac(t, y, yp, u, A, system_data, f=None):
    n = len(y)
    z = np.concatenate((y, yp))

    def fun_composite(t, z):
        y, yp = z[:n], z[n:]
        return F(t, y, yp, u, A, system_data)

    J = approx_derivative(lambda z: fun_composite(t, z), 
                            z, method="2-point", f0=f)
    J = J.reshape((n, 2 * n))
    Jy, Jyp = J[:, :n], J[:, n:]
    return Jy, Jyp

def f(t, z):
    y, yp = z[:7], z[7:]
    return np.concatenate((yp, F(t, y, yp)))

def get_microgrid_params():
    AC = np.array([[1, 0, 0, -1]]).T
    AR = np.array([[0, -1, 1, 0]]).T
    AL = np.array([[0, 0, -1, 1]]).T
    AV = np.array([[-1, 1, 0, 0]]).T
    AI = np.array([[0, 0, 0, 0]]).T
    splits = np.array([len(AC.T), 
                        len(AC.T) + len(AL.T), 
                        len(AC.T) + len(AL.T) + len(AC)])
    R = np.array(1.0)
    L = np.array(1.0)
    C = np.array(1.0)
    A = (AC, AR, AL, AV, AI)
    system_data = (R, L, C, splits)
    return (A, system_data)

def get_microgrid_policy(t):
    return np.array([[0, np.sin(t)]]).T # [i, v]

# time span
t0 = 0
t1 = 1e4
t_span = (t0, t1)
t_eval = np.linspace(0, 1, num=100)

# solver options
method = "Radau"
# method = "BDF" # alternative solver
atol = rtol = 1e-4

# system parameters
A, system_data = get_microgrid_params()
AC, AR, AL, AV, AI = A
len_y = len(AC.T) + len(AL.T) + len(AC) + len(AV.T)
# initial conditions
u = get_microgrid_policy
func = partial(F, A=A, system_data=system_data, u=u)
jac = partial(jac, A=A, system_data=system_data, u=u)
y0 = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
yp0 = np.array([0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0])
print(func(0.0, y0, yp0))
# y0, yp0, fnorm = consistent_initial_conditions(func, jac, t0, y0, yp0)
# solve DAE
start = time.time()
sol = solve_dae(func, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, t_eval=t_eval, jac=jac, stages=1)
end = time.time()
print(f'elapsed time: {end - start}')
t = sol.t
y = sol.y

# visualization
fig, ax = plt.subplots()
ax.set_xlabel("t")
ax.plot(t, y[0], label="q")
ax.plot(t, y[1], label="phi")
ax.plot(t, y[2], label="e")
ax.plot(t, y[3], label="jv")
ax.legend()
ax.grid()
plt.show()

[JonasBreuling]

Thank you for sharing your experience with the package. I reviewed your problem and noticed a few issues:

• The Jacobian matrix for your problem seems to be singular, as indicated by the LU decomposition. Debugging the integration methods shows that the initial step is never accepted because the system cannot be solved successfully.
• This singularity is why the call to consistent_initial_conditions reports ValueError: Index greater than one.
• Using s=1 in the Radau solver leads to the backward Euler method, which is not recommended due to its very low order.

I'm not very familiar with modeling electrical circuits, but I assume that there is an error in your equations. If you can share any books, weblinks, or papers discussing your problem, I can investigate further. Is this a well-known system I can study somewhere?

I'm not sure if this is helpful, but the one-transistor amplifier problem can be solved without any issues. See the code below.

import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize._numdiff import approx_derivative
from scipy_dae.integrate import solve_dae, consistent_initial_conditions

# Problem parameters
Ub = 6
R0 = 1000
R = 9000
alpha = 0.99
beta = 1e-6
Uf = 0.026
Ue = lambda t: 0.4 * np.sin(200 * np.pi * t)  

# The constant, singular mass matrix
c = 1e-6 * np.arange(1, 4)
M = np.zeros((5, 5))
M[0, 0] = -c[0]
M[0, 1] =  c[0]
M[1, 0] =  c[0]
M[1, 1] = -c[0]
M[2, 2] = -c[1]
M[3, 3] = -c[2]
M[3, 4] =  c[2]
M[4, 3] =  c[2]
M[4, 4] = -c[2]

# Hairer and Wanner's RADAU5 requires consistent initial conditions
# which they take to be
u0 = np.zeros(5)
u0[1] = Ub / 2
u0[2] = Ub / 2
u0[3] = Ub

# Perturb the algebraic variables to test initialization.
u0[3] = u0[3] + 0.1
u0[4] = u0[4] + 0.1

def F(t, u, up):
    f23 = beta * (np.exp((u[1] - u[2]) / Uf) - 1)
    dudt = np.array([
        -(Ue(t) - u[0]) / R0,
        -(Ub / R - u[1] * 2 / R - (1 - alpha) * f23),
        -(f23 - u[2] / R),
        -((Ub - u[3]) / R - alpha * f23),
        u[4] / R,
    ])
    return M @ up - dudt

def jac(t, y, yp, f):
    n = len(y)
    z = np.concatenate((y, yp))

    def fun_composite(t, z):
        y, yp = z[:n], z[n:]
        return F(t, y, yp)

    J = approx_derivative(lambda z: fun_composite(t, z), 
                          z, method="2-point", f0=f)
    J = J.reshape((n, 2 * n))
    Jy, Jyp = J[:, :n], J[:, n:]
    return Jy, Jyp

# time span
t0 = 0
t1 = 0.05
t_span = (t0, t1)

method = "Radau"
# method = "BDF" # alternative solver

# consistent initial conditions
up0 = np.zeros_like(u0)
print(f"u0: {u0}")
print(f"up0: {up0}")
u0, up0, fnorm = consistent_initial_conditions(F, jac, t0, u0, up0)
print(f"u0: {u0}")
print(f"up0: {up0}")
print(f"fnorm: {fnorm}")

# f0 = F(t0, u0, up0)
# Jy0, Jyp0 = jac(t0, u0, up0, f0)
# J0 = Jy0 + Jyp0
# print(f"Jy0:\n{Jy0}")
# print(f"Jyp0:\n{Jyp0}")
# print(f"J0:\n{Jyp0}")
# print(f"rank(Jy0):  {np.linalg.matrix_rank(Jy0)}")
# print(f"rank(Jyp0): {np.linalg.matrix_rank(Jyp0)}")
# print(f"rank(J0):   {np.linalg.matrix_rank(J0)}")
# print(f"J0.shape: {J0.shape}")
# exit()

# solver options
atol = rtol = 1e-6

# dae solution
start = time.time()
sol = solve_dae(F, t_span, u0, up0, atol=atol, rtol=rtol, method=method)
end = time.time()
print(f"elapsed time: {end - start}")
t = sol.t
u = sol.y
success = sol.success
status = sol.status
message = sol.message
print(f"success: {success}")
print(f"status: {status}")
print(f"message: {message}")
print(f"nfev: {sol.nfev}")
print(f"njev: {sol.njev}")
print(f"nlu: {sol.nlu}")

# visualization
fig, ax = plt.subplots()
ax.set_xlabel("t")
ax.set_title("One transistor amplifier problem")
ax.plot(t, Ue(t), label="input voltage Ue(t)")
ax.plot(t, u[4], label="output voltage U5(t)")
ax.grid()
ax.legend()
plt.show()

[nathan-t4]

Thanks Jonas! Any suggestions on how I can solve a general DAE with a singular Jacobian?

I am considering equation (4a) in this paper. These equations are essentially Kirchhoff's laws.

[JonasBreuling]

I've spotted a bug, see the working version below. You unpacked the matrices wrong

    # AC, AR, AL, AI, AV = A # <= this was the error!
    AC, AR, AL, AV, AI = A

Since AI is set to zero, this leads to a model with zero AV, which causes problems since jv is not involved in the residual.

But I'm still wondering about these equations. Why is AI zero? From where did you take the involved parameters for all the matrices? Do you know if this system is of index one or two?

And finally, the working code:

import numpy as np
import matplotlib.pyplot as plt
import time
from functools import partial
from scipy_dae.integrate import solve_dae, consistent_initial_conditions
from scipy.optimize._numdiff import approx_derivative

def F(t, y, yp, u, A, system_data):
    """Define implicit system of differential algebraic equations."""
    R, L, C, splits = system_data
    # AC, AR, AL, AI, AV = A # <= this was the error!
    AC, AR, AL, AV, AI = A

    y = np.nan_to_num(y)
    yp = np.nan_to_num(yp)
    qc, phi, e, jv = np.split(y, splits)
    q_t, phi_t, e_t, jv_t = np.split(yp, splits)
    i, v = u(t)

    F0 = AC @ q_t + AL @ (phi / L) + AV @ jv + AR @ (AR.T @ e) / R + AI @ i
    F1 = phi_t - AL.T @ e
    F2 = AC.T @ e - qc / C                # algebraic equation
    F3 = AV.T @ e - v                     # algebraic equation
    return np.concatenate((F0, F1, F2, F3))

def jac(t, y, yp, u, A, system_data, f=None):
    n = len(y)
    z = np.concatenate((y, yp))

    def fun_composite(t, z):
        y, yp = z[:n], z[n:]
        return F(t, y, yp, u, A, system_data)

    J = approx_derivative(lambda z: fun_composite(t, z), 
                          z, method="3-point", f0=f)
    J = J.reshape((n, 2 * n))
    Jy, Jyp = J[:, :n], J[:, n:]
    return Jy, Jyp

def get_microgrid_params():
    AC = np.array([[1, 0, 0, -1]]).T
    AR = np.array([[0, -1, 1, 0]]).T
    AL = np.array([[0, 0, -1, 1]]).T
    AV = np.array([[-1, 1, 0, 0]]).T
    AI = np.array([[0, 0, 0, 0]]).T
    splits = np.array([len(AC.T), 
                       len(AC.T) + len(AL.T), 
                       len(AC.T) + len(AL.T) + len(AC)])
    R = np.array([1.0])
    L = np.array([1.0])
    C = np.array([1.0])
    A = (AC, AR, AL, AV, AI)
    system_data = (R, L, C, splits)
    return (A, system_data)

def get_microgrid_policy(t):
    return np.array([[0, np.sin(t)]]).T # [i, v]

# time span
t0 = 0
t1 = 2e1
t_span = (t0, t1)
t_eval = np.linspace(t0, t1, num=int(1e3))

# solver options
method = "Radau"
# method = "BDF" # alternative solver
atol = rtol = 1e-4

# system parameters
A, system_data = get_microgrid_params()
AC, AR, AL, AV, AI = A
len_y = len(AC.T) + len(AL.T) + len(AC) + len(AV.T)

# initial conditions
u = get_microgrid_policy
func = partial(F, A=A, system_data=system_data, u=u)
jac = partial(jac, A=A, system_data=system_data, u=u)
y0 = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
yp0 = np.array([0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0])
f0 = func(t0, y0, yp0)
print(f"f0: {f0}")
Jy0, Jyp0 = jac(t0, y0, yp0)
J0 = Jy0 + Jyp0
print(f"Jy0:\n{Jy0}")
print(f"Jyp0:\n{Jyp0}")
print(f"J0:\n{Jyp0}")
print(f"rank(Jy0):  {np.linalg.matrix_rank(Jy0)}")
print(f"rank(Jyp0): {np.linalg.matrix_rank(Jyp0)}")
print(f"rank(J0):   {np.linalg.matrix_rank(J0)}")
print(f"J0.shape: {J0.shape}")
# y0, yp0, fnorm = consistent_initial_conditions(func, lambda t, y, yp, f: jac(t, y, yp), t0, y0, yp0)

# solve DAE
start = time.time()
sol = solve_dae(func, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, jac=jac, t_eval=t_eval)
end = time.time()
print(f'elapsed time: {end - start}')
t = sol.t
y = sol.y

# visualization
fig, ax = plt.subplots(2, 2)

ax[0, 0].set_xlabel("t")
ax[0, 0].plot(t, y[0], label="q")
ax[0, 0].legend()
ax[0, 0].grid()

ax[0, 1].set_xlabel("t")
ax[0, 1].plot(t, y[1], label="phi")
ax[0, 1].legend()
ax[0, 1].grid()

ax[1, 0].set_xlabel("t")
ax[1, 0].plot(t, y[2], label="e")
ax[1, 0].legend()
ax[1, 0].grid()

ax[1, 1].set_xlabel("t")
ax[1, 1].plot(t, y[3], label="jv")
ax[1, 1].legend()
ax[1, 1].grid()

plt.show()