Sundials solution greatly deviates for large perturbation, why?

New to Julia
,
1

Hello all,

First of, I am sorry for the massive code text.

The problem: in the following for small enough perturbations in P_in (ie jump = 1e-2), the IDA solver matches perfectly my solver. However, for anything larger, certain variables take values that violate the input domain of certain functions. I am unsure whether that it due to the problem being too stiff, inappropriate step size or anything else. Actually, I have no idea why this happens.

using Plots
using LinearAlgebra, ForwardDiff
using DifferentialEquations, Sundials

import Base: convert
convert(::Type{Float64}, x::ForwardDiff.Dual) = ForwardDiff.value(x)

font_size=15
line_width = 1

setprecision(256) 

mutable struct ElectrolyzerParameters
    nx::Int64
    ny::Int64
    nu::Int64
    nd::Int64
    ndadd::Int64
    flag::Bool
    scale_factor::Float64
    epsilon::Float64

    cp_lye::Float64       # Specific heat capacity of the lye solution [kJ/(kg·K)]
    U_tn::Float64         # Thermoneutral voltage [V]
    U_rev::Float64        # Reversible voltage [V]
    nc::Int               # Number of cells in the stack [-]
    z::Int                # Number of electrons transferred per reaction [-]
    F::Float64            # Faraday's constant [C/mol]
    A::Float64            # Electrode area [m²]
    CP_el::Float64        # Heat capacity of the stack [kJ/C]
    r1::Float64           # Ohmic resistance coefficient 1 [Ohm·m²]
    r2::Float64           # Ohmic resistance coefficient 2 [Ohm·m²·C⁻¹]
    s::Float64            # Activation overvoltage coefficient [V]
    t1::Float64           # Activation overvoltage coefficient [m²·A⁻¹]
    t2::Float64           # Activation overvoltage coefficient [m²·C·A⁻¹]
    t3::Float64           # Activation overvoltage coefficient [m²·C²·A⁻¹]
    f1::Float64           # Faraday efficiency parameter [A²·m^{-3} formerly in m·A^2·cm^-4]
    f2::Float64           # Faraday efficiency parameter
    A_c::Float64          # Active area of heat transferred [m²]
    h_c::Float64          # Heat transfer coefficient [W/(m²·K)]
    
    f_H2O_1_bar::Float64  # Nominal water outflow from separator 1 [kg/s or mol/s]
    Kp_H2O_1::Float64     # Proportional gain for water flow control in separator 1 [ (kg/s)/mol or similar ]
    n_H2O_1_set::Float64
    f_H2O_2_bar::Float64  # Nominal water outflow from separator 2 [kg/s or mol/s]
    Kp_H2O_2::Float64     # Proportional gain for water flow control in separator 2 [ (kg/s)/mol or similar ]
    n_H2O_2_set::Float64
    f_O2_bar::Float64     # Nominal oxygen outflow from separator 1 [mol/s]
    Kp_O2::Float64        # Proportional gain for oxygen flow control [ (mol/s)/mol or similar ]
    n_O2_set::Float64
    f_H2_bar::Float64     # Nominal hydrogen outflow from separator 2 [mol/s]
    Kp_H2::Float64        # Proportional gain for hydrogen flow control [ (mol/s)/mol or similar ]
    n_H2_set::Float64
    f_H2O_mu_bar::Float64
    Kp_H2O_mu::Float64

    sigma_T_stack_d::Float64
    sigma_n_H2O_1_d::Float64
    sigma_n_H2O_2_d::Float64
    sigma_T_1_d::Float64
    sigma_T_2_d::Float64

    u::Array{Float64,1}
    d::Array{Float64,1}
end

function AELmodel_DAE_julia1!(out, dx, x, p, t)
    # Unpack Differential & Algebraic States
    T_stack, n_H2O_1, n_H2O_2, T_1, T_2, n_O2, n_H2 = x[1:p.nx]
    n_H2O_1 *= p.scale_factor
    n_H2O_2 *= p.scale_factor
    n_O2 *= p.scale_factor
    n_H2 *= p.scale_factor
    # Unpack Algebraic States
    U_cell, I_cell, f_H2O_mix, T_mix, T_hex, f_H2O_in, T_in = x[p.nx+1:end]
    I_cell *= p.scale_factor
    f_H2O_mix *= p.scale_factor
    f_H2O_in *= p.scale_factor
    # Unpack Inputσ
    Q_cool, P_in = p.u
    Q_cool *= p.scale_factor
    P_in *= p.scale_factor
    # Unpack Disturbances
    T_amp, T_mu = p.d

    # P-controllers
    f_H2O_1 = p.f_H2O_1_bar + p.Kp_H2O_1 * ( p.n_H2O_1_set - n_H2O_1 )
    f_H2O_2 = p.f_H2O_2_bar + p.Kp_H2O_2 * ( p.n_H2O_2_set - n_H2O_2 )
    f_O2 = p.f_O2_bar + p.Kp_O2 * ( p.n_O2_set - n_O2 )
    f_H2 = p.f_H2_bar + p.Kp_H2 * ( p.n_H2_set - n_H2 )
    f_H2O_mu = p.f_H2O_mu_bar + p.Kp_H2O_mu * ( p.n_H2O_2_set - n_H2O_2 )

    # Faraday efficiency
    eta_F = (I_cell / p.A)^2 * p.f2 / (p.f1 + (I_cell / p.A)^2)

    # Reaction rate
    r = p.nc * eta_F * I_cell / (p.z * p.F)

    Q_gen = p.nc * I_cell * (U_cell - p.U_tn)
    Q_liq = f_H2O_in * p.cp_lye * ( T_in - T_stack )
    Q_amb = p.A_c * p.h_c * ( T_stack - T_amp )
    n_H2O_1_in = f_H2O_in / 2 + r
    n_H2O_2_in = f_H2O_in / 2 - 2 * r

    U_ohm = ( p.r1 + p.r2 * T_stack ) * I_cell / p.A
    U_act = p.s * log10( ( p.t1 + p.t2 / T_stack + p.t3 / T_stack^2 ) * I_cell / p.A + 1 )

    # Differential equations
    out[1] = (Q_gen + Q_liq - Q_amb) / p.CP_el - dx[1]
    out[2] = (n_H2O_1_in - f_H2O_1 - dx[2] ) / p.scale_factor
    out[3] = (n_H2O_2_in - f_H2O_2 - dx[3] ) / p.scale_factor
    out[4] = n_H2O_1_in * (T_stack - T_1) / n_H2O_1 - dx[4]
    out[5] = n_H2O_2_in * (T_stack - T_2) / n_H2O_2 - dx[5]
    out[6] = (r / 2 - f_O2 - dx[6] ) / p.scale_factor
    out[7] = (r - f_H2 - dx[7] ) / p.scale_factor

    # Algebraic equations
    out[8] = U_cell - ( p.U_rev + U_ohm + U_act )
    out[9] = (P_in - p.nc * U_cell * I_cell) / p.scale_factor
    out[10] = (f_H2O_mix - (f_H2O_1 + f_H2O_2)) / p.scale_factor
    out[11] = T_mix * f_H2O_mix * p.cp_lye - p.cp_lye * ( T_1 * f_H2O_1 + T_2 * f_H2O_2 )
    out[12] = T_hex - ( T_mix - Q_cool / ( f_H2O_mix * p.cp_lye ) )
    out[13] = (f_H2O_in - ( f_H2O_mix + f_H2O_mu )) / p.scale_factor
    out[14] = f_H2O_in * T_in * p.cp_lye - p.cp_lye * ( f_H2O_mix * T_hex + f_H2O_mu * T_mu )
    
end

function AELmodel_DAE(t, x, u, d, p)
    # Unpack Differential States
    T_stack, n_H2O_1, n_H2O_2, T_1, T_2, n_O2, n_H2 = x[1:p.nx]
    n_H2O_1 *= p.scale_factor
    n_H2O_2 *= p.scale_factor
    n_O2 *= p.scale_factor
    n_H2 *= p.scale_factor
    # Unpack Algebraic States
    U_cell, I_cell, f_H2O_mix, T_mix, T_hex, f_H2O_in, T_in = x[p.nx+1:end]
    I_cell *= p.scale_factor
    f_H2O_mix *= p.scale_factor
    f_H2O_in *= p.scale_factor
    # Unpack Inputs
    Q_cool, P_in = u
    Q_cool *= p.scale_factor
    # f_H2O_mu *= p.scale_factor
    P_in *= p.scale_factor
    # Unpack Disturbances
    T_amp, T_mu = d

    # P-controllers
    f_H2O_1 = p.f_H2O_1_bar + p.Kp_H2O_1 * ( p.n_H2O_1_set - n_H2O_1 )
    f_H2O_2 = p.f_H2O_2_bar + p.Kp_H2O_2 * ( p.n_H2O_2_set - n_H2O_2 )
    f_O2 = p.f_O2_bar + p.Kp_O2 * ( p.n_O2_set - n_O2 )
    f_H2 = p.f_H2_bar + p.Kp_H2 * ( p.n_H2_set - n_H2 )
    f_H2O_mu = p.f_H2O_mu_bar + p.Kp_H2O_mu * ( p.n_H2O_2_set - n_H2O_2 )

    # Faraday efficiency
    eta_F = (I_cell / p.A)^2 * p.f2 / (p.f1 + (I_cell / p.A)^2)
    deta_FdI = 2 * p.f1 * p.f2 * I_cell / (p.A * ( p.f1 + ( I_cell / p.A )^2 ) )^2

    # Reaction rate
    r = p.nc * eta_F * I_cell / (p.z * p.F)
    drdI = p.nc * ( deta_FdI * I_cell + eta_F ) / (p.z * p.F)

    Q_gen = p.nc * I_cell * (U_cell - p.U_tn)
    Q_liq = f_H2O_in * p.cp_lye * ( T_in - T_stack )
    Q_amb = p.A_c * p.h_c * ( T_stack - T_amp )
    n_H2O_1_in = f_H2O_in / 2 + r
    n_H2O_2_in = f_H2O_in / 2 - 2 * r

    U_ohm = ( p.r1 + p.r2 * T_stack ) * I_cell / p.A
    dUohmdTstack = p.r2 * I_cell / p.A
    dUohmdI = ( p.r1 + p.r2 * T_stack ) / p.A
    U_act = p.s * log10( ( p.t1 + p.t2 / T_stack + p.t3 / T_stack^2 ) * I_cell / p.A + 1 )
    dUactdTstack = - p.s * I_cell * ( p.t2 / T_stack^2 + 2 * p.t3 / T_stack^3 ) / ( ( ( p.t1 + p.t2 / T_stack + p.t3 / T_stack^2 ) * I_cell / p.A + 1 ) * p.A )
    dUactdI = p.s * ( p.t1 + p.t2 / T_stack + p.t3 / T_stack^2 ) / ( ( ( p.t1 + p.t2 / T_stack + p.t3 / T_stack^2 ) * I_cell / p.A + 1 ) * p.A )

    # Differential equations
    f1 = (Q_gen + Q_liq - Q_amb) / p.CP_el
    f2 = n_H2O_1_in - f_H2O_1
    f3 = n_H2O_2_in - f_H2O_2
    f4 = n_H2O_1_in * (T_stack - T_1) / (n_H2O_1)
    f5 = n_H2O_2_in * (T_stack - T_2) / (n_H2O_2)
    f6 = r / 2 - f_O2
    f7 = r - f_H2

    # Algebraic equations
    g1 = U_cell - ( p.U_rev + U_ohm + U_act )
    g2 = P_in - p.nc * U_cell * I_cell
    g3 = f_H2O_mix - (f_H2O_1 + f_H2O_2)
    g4 = T_mix * f_H2O_mix * p.cp_lye - p.cp_lye * ( T_1 * f_H2O_1 + T_2 * f_H2O_2 )
    g5 = T_hex - ( T_mix - Q_cool / ( f_H2O_mix * p.cp_lye ) )
    g6 = f_H2O_in - ( f_H2O_mix + f_H2O_mu )
    g7 = f_H2O_in * T_in * p.cp_lye - p.cp_lye * ( f_H2O_mix * T_hex + f_H2O_mu * T_mu )
    
    F = [f1; f2/p.scale_factor; f3/p.scale_factor; f4; f5; f6/p.scale_factor; f7/p.scale_factor; g1; g2/p.scale_factor; g3/p.scale_factor; g4; g5; g6/p.scale_factor; g7]

    if p.flag
        J = zeros(p.nx+p.ny, p.nx+p.ny)
        J[1,1] = - ( f_H2O_in * p.cp_lye + p.A_c * p.h_c ) / p.CP_el
        J[1,8] = p.nc * I_cell / p.CP_el
        J[1,9] = p.scale_factor * p.nc * (U_cell - p.U_tn) / p.CP_el
        J[1,13] = p.scale_factor * p.cp_lye * ( T_in - T_stack ) / p.CP_el
        J[1, 14] = f_H2O_in * p.cp_lye / p.CP_el
        
        J[2,2] = p.Kp_H2O_1
        J[2,9] = drdI
        J[2,13] = 0.5

        J[3,3] = p.Kp_H2O_2
        J[3,9] = - 2.0 * drdI
        J[3,13] = 0.5

        J[4,1] = n_H2O_1_in / (n_H2O_1)
        J[4,2] = - n_H2O_1_in * ( T_stack - T_1 ) / (n_H2O_1)^2
        J[4,4] = - n_H2O_1_in / (n_H2O_1)
        J[4,9] = drdI * ( T_stack - T_1 ) / (n_H2O_1)
        J[4,13] = 0.5 * ( T_stack - T_1 ) / (n_H2O_1)
        
        J[5,1] = n_H2O_2_in / (n_H2O_2)
        J[5,3] = - n_H2O_2_in * ( T_stack - T_2 ) / (n_H2O_2)^2
        J[5,5] = - n_H2O_2_in / (n_H2O_2)
        J[5,9] = - 2.0 * drdI * ( T_stack - T_2 ) / (n_H2O_2)
        J[5,13] = 0.5 * ( T_stack - T_2 ) / (n_H2O_2)

        J[6,6] = p.Kp_O2
        J[6,9] = 0.5 * drdI

        J[7,7] = p.Kp_H2
        J[7,9] = drdI

        J[8,1] = - dUohmdTstack - dUactdTstack
        J[8,8] = 1.0
        J[8,9] = - p.scale_factor * (dUohmdI + dUactdI)

        J[9,8] = - p.nc * I_cell / p.scale_factor
        J[9,9] = - p.nc * U_cell

        J[10,2] = p.Kp_H2O_1
        J[10,3] = p.Kp_H2O_2
        J[10,10] = 1.0

        J[11,2] = p.scale_factor * p.cp_lye * T_1 * p.Kp_H2O_1
        J[11,3] = p.scale_factor * p.cp_lye * T_2 * p.Kp_H2O_2
        J[11,4] = - p.cp_lye * f_H2O_1
        J[11,5] = - p.cp_lye * f_H2O_2
        J[11,10] = p.scale_factor * p.cp_lye * T_mix
        J[11,11] = p.cp_lye * f_H2O_mix

        J[12,10] = -  p.scale_factor * Q_cool / ( p.cp_lye * (f_H2O_mix)^2 )
        J[12,11] = - 1.0
        J[12,12] = 1.0

        J[13,3] = p.Kp_H2O_mu
        J[13,10] = - 1.0
        J[13,13] = 1.0

        J[14,3] = p.scale_factor * p.cp_lye * T_mu * p.Kp_H2O_mu
        J[14,10] = - p.scale_factor * p.cp_lye * T_hex
        J[14,12] = - p.cp_lye * f_H2O_mix
        J[14,13] =  p.scale_factor * p.cp_lye * T_in
        J[14,14] = p.cp_lye * f_H2O_in
        
    else
        J = nothing
    end

    return F, J
end

function ImplicitEulerFixedStepSizeDAE(funJac, tspan, dt, x0, n_x, n_y, args...;
    tol=1e-8, maxit=1e2)
    # Compute the time step.
    t0, tf = tspan
    N = Int64(round((tf - t0) / dt))
    n = length(x0)

    # Allocate memory
    X = zeros(Float64, n, N+1)
    T = zeros(Float64, N+1)
    res_all = Vector{Vector{Float64}}(undef, N)

    # Set initial conditions.
    T[1] = t0
    X[:, 1] = x0

    # Time-stepping loop.
    for k in 1:N
        # Evaluate funJac at the current time and state.
        F, _ = funJac(T[k], X[:, k], args...)
        # Extract the differential part f (first n_x components)
        f = F[1:n_x]
        T[k+1] = T[k] + dt
        # Use an explicit Euler step for the differential part as initial guess;
        # the algebraic part is left unchanged.
        xinit = copy(X[1:n_x+n_y, k])
        xinit[1:n_x] .= xinit[1:n_x] .+ dt * f
        # Use Newton's method to solve the implicit step.
        X[1:n_x+n_y, k+1], res_hist = NewtonsMethodDAE(funJac, T[k], X[:, k], dt, xinit, n_x, n_y, args...; tol, maxit)
        res_all[k] = res_hist
    end

    X = X'
    return T, X, res_all
end

function NewtonsMethodDAE(funJac, tk, Xk, dt, Xinit, n_x, n_y, args...; tol=1e-8, maxit=1e2)
    t = tk + dt
    X = copy(Xinit)
    # Evaluate the function and its Jacobian at (t, X)
    F, J = funJac(t, X, args...)
    # Partition F into f (differential part) and g (algebraic part)
    f = F[1:n_x]
    g = F[n_x+1:n_x + n_y]
    # Residual function
    R = vcat( X[1:n_x] - Xk[1:n_x] - dt * f,  g )
    res_hist = []

    k = 0
    while k < maxit && norm(R, Inf) > tol
        k += 1
        # Partition the Jacobian J:
        #   A = ∂f/∂x, B = ∂f/∂y, C = ∂g/∂x, D = ∂g/∂y.
        A = J[1:n_x, 1:n_x]
        B = J[1:n_x, n_x+1:n_x + n_y]
        C = J[n_x+1:n_x + n_y, 1:n_x]
        D = J[n_x+1:n_x + n_y, n_x+1:n_x + n_y]
        # Newton matrix for the DAE:
        # [ I - dt*A    - dt*B ]
        # [    C          D    ]
        J_newton = [ I(n_x)-dt*A   -dt*B;
        C                D ]
        # Solve for the update
        dx = J_newton \ R
        X -= dx
        # Reevaluate F, J and R at updated X
        F, J = funJac(t, X, args...)
        f = F[1:n_x]
        g = F[n_x+1:n_x + n_y]
        R = vcat( X[1:n_x] - Xk[1:n_x] - dt * f,  g )
        push!(res_hist, norm(R, Inf))
    end
    return X, res_hist
end

mm_H2O = 18.015         # Molar Mass of water            [g/mol]
mm_O2 = 16 * 2          # Molar Mass of oxygen           [g/mol]
mm_H2 = 1.008 * 2       # Molar Mass of hydrogen         [g/mol]

nx = 7
ny = 7
nu = 2
nd = 2
ndadd = 5
flag = true
scale_factor = 1e0
epsilon = 1e-5

cp_lye = 3.101*mm_H2O                     # Specific heat capacity of the lye solution          [J/(mol·K)]
U_tn = 1.482                              # Thermoneutral voltage                               [V]
U_rev = 1.229                             # Reversible voltage                                  [V]
nc = 230                                  # Number of cells in the stack                        [-]
z = 2                                     # Number of electrons transferred per reaction        [-]
F = 96485.33                              # Faraday's constant                                  [C/mol]
A = 2.6                                   # Electrode area                                      [m²]
CP_el = 51322.1*1e3                       # Heat capacity of the stack                          [J/C]
r1 = 2.18e-4                              # Ohmic resistance coefficient 1                      [Ohm·m²]
r2 = -4.25e-7                             # Ohmic resistance coefficient 2                      [Ohm·m²·C⁻¹]
s = 117.93e-3                             # Activation overvoltage coefficient                  [V]
t1 = -145.29e-3                           # Activation overvoltage coefficient                  [m²·A⁻¹]
t2 = 11.794                               # Activation overvoltage coefficient                  [m²·C·A⁻¹]
t3 = 395.68                               # Activation overvoltage coefficient                  [m²·C²·A⁻¹]
f1 = 120.0                                # Faraday efficiency parameter                        [m·A^2·cm^-4]
f2 = 0.98                                 # Faraday efficiency parameter                        [?]
A_c = 131.56                              # Active area of heat transferred                     [m²]
h_c = 5.5                                 # Heat transfer coefficient                           [W/(m²·K)]

# Set up fake p-controlers, f_out = f_out_set + K ( n_set - n )
f_H2O_1_bar = 6.6e4/2
Kp_H2O_1 =  1e0
n_H2O_1_set =  1e5/2

f_H2O_2_bar = 6.6e4/2
Kp_H2O_2 =  1e0
n_H2O_2_set = 1e5/2

f_O2_bar = 3.0/2
Kp_O2 = 1e0  
n_O2_set = 1e5/2

f_H2_bar = 6.0/2
Kp_H2 = 1e0  
n_H2_set = 1e5/2

f_H2O_mu_bar = 5e0
Kp_H2O_mu = 1e0

sigma_T_stack_d = 1.0
sigma_n_H2O_1_d = 1.0
sigma_n_H2O_2_d = 1.0
sigma_T_1_d = 1.0
sigma_T_2_d = 1.0

u = [5e5; 2e3]
d = [2e1; 2e1]

p = ElectrolyzerParameters(nx, ny, nu, nd, ndadd, flag, scale_factor, epsilon,
    cp_lye, U_tn, U_rev, nc, z, F, A, CP_el, r1, r2, s, t1, t2, t3, f1, f2, A_c, h_c, 
    f_H2O_1_bar, Kp_H2O_1, n_H2O_1_set, f_H2O_2_bar, Kp_H2O_2, n_H2O_2_set, f_O2_bar, Kp_O2, n_O2_set,f_H2_bar, Kp_H2, n_H2_set, f_H2O_mu_bar, Kp_H2O_mu,
    sigma_T_stack_d, sigma_n_H2O_1_d, sigma_n_H2O_2_d, sigma_T_1_d, sigma_T_2_d, u, d
)

t0 = 0.0
tf = 60.0 * 60.0 * 1
Ts = 15
tspan = range(t0, tf, Int64(round(tf/Ts)))
N = length(tspan)
dt = 5.0

p.flag = true

d = repeat(d0,1,N)

jump = 1e-1
P_ins = u_ss[2] * [1*ones(Int64(1*N/3),1); (1+jump)*ones(Int64(2*N/3),1)]
jump = 0e-1
Q_cools = u_ss[1] * [1*ones(Int64(1*N/3),1); (1+jump)*ones(Int64(2*N/3),1)]
u = vcat(Q_cools', P_ins')

X0_1 = copy(X0)
T = zeros(0,1)
X = zeros(0, length(X0_1))

X0_2 = copy(X0)
dX0,_ = AELmodel_DAE(t0, X0_2, u[:,1], d[:,1], p)
T2 = zeros(0,1)
X2 = zeros(0, length(X0_2))
differential_vars = [true for i in 1:p.nx] 
append!(differential_vars, [false for i in 1:p.ny])

for i in 1:N-1
    # Custom  DAESolver
    Tk, Xk, Rk = ImplicitEulerFixedStepSizeDAE(AELmodel_DAE, (tspan[i], tspan[i+1]), dt, X0_1, p.nx, p.ny, u[:,i], d[:,i], p)
    T = [T; Tk]
    X = [X; Xk]
    X0_1 = Xk[end,:]

    # Julia DAE Solver
    p.u = u[:,i]
    p.d = d[:,i]
    prob = DAEProblem(AELmodel_DAE_julia1!, dX0, X0_3, (tspan[i], tspan[i+1]), p; differential_vars=differential_vars)
    sol = solve(prob, IDA()) #, dt=dt, adaptive=false)
    tmp = transpose(Array(sol))
    T3 = [T3; sol.t[end]]
    X3 = [X3; tmp[end,:]']
    X0_3 = X3[end,:]
    dX0,_ = AELmodel_DAE(t0, X0_3, u[:,i+1], d[:,i+1], p)
end


DomainError with -1.0181106338562343e7:
log10 was called with a negative real argument but will only return a complex result if called with a complex argument. Try log10(Complex(x)).
DomainError detected in the user `f` function. This occurs when the domain of a function is violated.
For example, `log(-1.0)` is undefined because `log` of a real number is defined to only output real
numbers, but `log` of a negative number is complex valued and therefore Julia throws a DomainError
by default. Cases to be aware of include:

* `log(x)`, `sqrt(x)`, `cbrt(x)`, etc. where `x<0`
* `x^y` for `x<0` floating point `y` (example: `(-1.0)^(1/2) == im`)

Within the context of SciML, this error can occur within the solver process even if the domain constraint
would not be violated in the solution due to adaptivity. For example, an ODE solver or optimization
routine may check a step at `new_u` which violates the domain constraint, and if violated reject the
step and use a smaller `dt`. However, the throwing of this error will have halted the solving process.

Thus the recommended fix is to replace this function with the equivalent ones from NaNMath.jl
(https://github.com/JuliaMath/NaNMath.jl) which returns a NaN instead of an error. The solver will then
effectively use the NaN within the error control routines to reject the out of bounds step. Additionally,
one could perform a domain transformation on the variables so that such an issue does not occur in the
definition of `f`.

For more information, check out the following FAQ page:
https://docs.sciml.ai/Optimization/stable/API/FAQ/#The-Solver-Seems-to-Violate-Constraints-During-the-Optimization,-Causing-DomainErrors,-What-Can-I-Do-About-That?


Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math .\math.jl:33
  [2] _log
    @ .\special\log.jl:295 [inlined]
  [3] log10(x::Float64)
    @ Base.Math .\special\log.jl:262
  [4] AELmodel_DAE_julia1!(out::Vector{Float64}, dx::Vector{Float64}, x::Vector{Float64}, p::ElectrolyzerParameters, t::Float64)
    @ Main c:\Users\dsous\Documents\DTU\Thesis\John's\Code\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X24sZmlsZQ==.jl:106
  [5] DAEFunction
    @ C:\Users\dsous\.julia\packages\SciMLBase\sYmAV\src\scimlfunctions.jl:2492 [inlined]
  [6] idasolfun(t::Float64, u::Ptr{Sundials._generic_N_Vector}, du::Ptr{Sundials._generic_N_Vector}, resid::Ptr{Sundials._generic_N_Vector}, funjac::Sundials.FunJac{1, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Nothing, Nothing, ElectrolyzerParameters, Nothing, Nothing, Nothing, Nothing, Vector{Float64}})
    @ Sundials C:\Users\dsous\.julia\packages\Sundials\2hgkX\src\common_interface\function_types.jl:93
  [7] IDACalcIC
    @ C:\Users\dsous\.julia\packages\Sundials\2hgkX\lib\libsundials_api.jl:3312 [inlined]
  [8] initialize_dae!(integrator::Sundials.IDAIntegrator{1, ElectrolyzerParameters, DAESolution{Float64, 2, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, DAEProblem{Vector{Float64}, Vector{Float64}, Tuple{Float64, Float64}, true, ElectrolyzerParameters, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, @Kwargs{}, Vector{Bool}}, IDA{:Dense, Nothing, Nothing}, SciMLBase.HermiteInterpolation{Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}, SciMLBase.DEStats, Nothing, Nothing}, IDA{:Dense, Nothing, Nothing}, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Sundials.FunJac{1, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Nothing, Nothing, ElectrolyzerParameters, Nothing, Nothing, Nothing, Nothing, Vector{Float64}}, Nothing, Sundials.DEOptions{DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Vector{Float64}, Vector{Float64}, Nothing, CallbackSet{Tuple{}, Tuple{}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing, Sundials.IDADefaultInit}, initializealg::Sundials.IDADefaultInit)
    @ Sundials C:\Users\dsous\.julia\packages\Sundials\2hgkX\src\common_interface\integrator_utils.jl:212
  [9] __init(prob::DAEProblem{Vector{Float64}, Vector{Float64}, Tuple{Float64, Float64}, true, ElectrolyzerParameters, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, @Kwargs{}, Vector{Bool}}, alg::IDA{:Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}; verbose::Bool, dt::Nothing, dtmax::Float64, save_on::Bool, save_start::Bool, callback::Nothing, abstol::Float64, reltol::Float64, saveat::Vector{Float64}, tstops::Vector{Float64}, d_discontinuities::Vector{Float64}, maxiters::Int64, timeseries_errors::Bool, dense_errors::Bool, save_everystep::Bool, save_idxs::Nothing, dense::Bool, save_timeseries::Nothing, save_end::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, advance_to_tstop::Bool, stop_at_next_tstop::Bool, userdata::Nothing, initializealg::Sundials.IDADefaultInit, kwargs::@Kwargs{})
    @ Sundials C:\Users\dsous\.julia\packages\Sundials\2hgkX\src\common_interface\solve.jl:1316
 [10] __init
    @ C:\Users\dsous\.julia\packages\Sundials\2hgkX\src\common_interface\solve.jl:981 [inlined]
 [11] __solve(prob::DAEProblem{Vector{Float64}, Vector{Float64}, Tuple{Float64, Float64}, true, ElectrolyzerParameters, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, @Kwargs{}, Vector{Bool}}, alg::IDA{:Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; calculate_error::Bool, kwargs::@Kwargs{})
    @ Sundials C:\Users\dsous\.julia\packages\Sundials\2hgkX\src\common_interface\solve.jl:15
 [12] __solve (repeats 5 times)
    @ C:\Users\dsous\.julia\packages\Sundials\2hgkX\src\common_interface\solve.jl:3 [inlined]
 [13] #solve_call#35
    @ C:\Users\dsous\.julia\packages\DiffEqBase\PbBEl\src\solve.jl:635 [inlined]
 [14] solve_call
    @ C:\Users\dsous\.julia\packages\DiffEqBase\PbBEl\src\solve.jl:592 [inlined]
 [15] #solve_up#44
    @ C:\Users\dsous\.julia\packages\DiffEqBase\PbBEl\src\solve.jl:1128 [inlined]
 [16] solve_up
    @ C:\Users\dsous\.julia\packages\DiffEqBase\PbBEl\src\solve.jl:1106 [inlined]
 [17] solve(prob::DAEProblem{Vector{Float64}, Vector{Float64}, Tuple{Float64, Float64}, true, ElectrolyzerParameters, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, @Kwargs{}, Vector{Bool}}, args::IDA{:Dense, Nothing, Nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\dsous\.julia\packages\DiffEqBase\PbBEl\src\solve.jl:1043
 [18] solve(prob::DAEProblem{Vector{Float64}, Vector{Float64}, Tuple{Float64, Float64}, true, ElectrolyzerParameters, DAEFunction{true, SciMLBase.FullSpecialize, typeof(AELmodel_DAE_julia1!), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, @Kwargs{}, Vector{Bool}}, args::IDA{:Dense, Nothing, Nothing})
    @ DiffEqBase C:\Users\dsous\.julia\packages\DiffEqBase\PbBEl\src\solve.jl:1033
 [19] top-level scope
    @ c:\Users\dsous\Documents\DTU\Thesis\John's\Code\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X24sZmlsZQ==.jl:450

the system is initiated from a steady state btw under some constraints (T_stack < 80, Q_cool < 1e6)

Thank you very much in advance, I look forward to all replies.

2

Did you try using NaNMath.log10 instead like the error message suggests?

Within the context of SciML, this error can occur within the solver process even if the domain constraint
would not be violated in the solution due to adaptivity. For example, an ODE solver or optimization
routine may check a step at new_u which violates the domain constraint, and if violated reject the
step and use a smaller dt. However, the throwing of this error will have halted the solving process.

Thus the recommended fix is to replace this function with the equivalent ones from NaNMath.jl
(GitHub - JuliaMath/NaNMath.jl: Julia math built-ins which return NaN and accumulator functions which ignore NaN) which returns a NaN instead of an error. The solver will then
effectively use the NaN within the error control routines to reject the out of bounds step. Additionally,
one could perform a domain transformation on the variables so that such an issue does not occur in the
definition of f.

3

Yes I did. Completely nonsensical results along with error:

[IDAS ERROR]  IDACalcIC
  Newton/Linesearch algorithm failed to converge.

along each time interval after the perturbation occurs. For small perturbations though it works as intended.

4

You’re not really allowed to do this (Style Guide · The Julia Language). If you previously got an error that went away after adding this definition, you’re now almost certainly getting incorrect derivatives somewhere inside the DifferentialEquations solvers.

I don’t see any BigFloats in your code, so setprecision will have no effect.

5

I removed the definition and the error remains. Also, maybe I am wrong, but if the derivative computations were incorrect, shouldn’t the solver fail to solve the system for small perturbations as well?

6

Not necessarily. Newton’s method could still converge with an incorrect Jacobian, depending on how close it is. You lose quadratic convergence though.

Sundials doesn’t use autodiff so it’s unrelated, but yes that piece of code is very bad and I recommend not overloading autodiff internals.

As for your core problem though, have you tried writing it in mass matrix form for Rodas5P or like methods?

The issue is that time stepping needs to propose steps and then reject them. If you error in your f function that would stop adaptivity from working because you will just error before it adapts. Instead what you want is to signal with a NaN and then let it reject the step. That said, I’m a bit surprised that Sundials fails with NaN signals, so trying Rodas5P or other native Julia methods should be fine with it.

7

Just to confirm: Did you run your code in a fresh Julia session after removing the convert method, or in the same Julia session you were working in before? And did you get the exact same error, or a (slightly) different one?

Edit: don’t worry about this given that @ChrisRackauckas said autodiff is a red herring here.

8

Regarding IDA, i fixed the time step dt to be equal to the sampling time Ts (15 seconds, non adaptive), increased the setpoints of the fake p-controllers and rescaled the units from kilounits to units and it worked for large perturbations. I don’t know why but it need. However, it does not work for a smaller time step and weirdly enough does not work for smaller perturbations (same error but smaller number).

I also tested Rodas5P, but no matter the step size and the size of the perturbations (and NaNMaths.log10) i get this:

┌ Warning: Instability detected. Aborting
└ @ SciMLBase C:\Users\dsous\.julia\packages\SciMLBase\sYmAV\src\integrator_interface.jl:631

for an adaptive step, i get warning like this:

┌ Warning: At t=51.87840874722878, dt was forced below floating point epsilon 7.105427357601002e-15, and step error estimate = 1.2216632914214376. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase C:\Users\dsous\.julia\packages\SciMLBase\sYmAV\src\integrator_interface.jl:623

If I dont use NaNMaths, i get the domain violation error.

9

Check the condition number of the Jacobian at that time, is the DAE index 1?

10

The condition number is massive through out the whole simulation

Condition Number: 1.9549581435838647e11

the DAE has index 1, that is, unless I am mistaken, I can differentiate the algebraic equations once to reveal and isolate the algebraic variables.

11

Usually this is a sign that your DAE is changing index. Do you have any algebraic equations which are vanishing, i.e. something going to zero and getting a 0 = 0 condition?

12

The algebraic states are always non-zero, no trivial solutions. The only thing I can imagine going wrong is that the model formulation is incorrect and the implicit euler works by “accident”. If I linearize the whole DAE system with respect to states, inputs and disturbances, the system matrices are incredibly ill-conditioned:

=======Analytical Scaled=======
Matrix: df/dx = Ax1
  Rank: 7
  Size: (7, 7)
  Condition Number: 27.630560054881187
  Eigenvalues: [-0.6606, -0.6600, -0.0719, 1.0000, 1.0000, 1.0000, 1.0000]

Matrix: df/dy = Ay1
  Rank: 3
  Size: (7, 7)
  Condition Number: Inf
  Eigenvalues: [-0.0165, 0.0000, 0.0000, 0.0000, 0.0000, 0.0177, 0.0215]

Matrix: df/du = B1
  Rank: 0
  Size: (7, 2)
  Condition Number: Inf
  Eigenvalues: Not applicable (matrix is not square)

Matrix: df/dd = E1
  Rank: 1
  Size: (7, 2)
  Condition Number: Inf
  Eigenvalues: Not applicable (matrix is not square)

Matrix: dg/dx = Ax2
  Rank: 4
  Size: (7, 7)
  Condition Number: Inf
  Eigenvalues: [-1844501.3994, 0.0000, 0.0000, 0.0000, 0.0000, 0.0070, 1.0000]

Matrix: dg/dy = Ay2
  Rank: 7
  Size: (7, 7)
  Condition Number: 2.574868839963498e7
  Eigenvalues: [-417.2775, 1.0000, 1.0000, 1.0000, 1.2509, 3688064.0036, 3688376.9353]

Matrix: dg/du = B2
  Rank: 2
  Size: (7, 2)
  Condition Number: 3.6880640035708477e6
  Eigenvalues: Not applicable (matrix is not square)

Matrix: dg/dd = E2
  Rank: 1
  Size: (7, 2)
  Condition Number: Inf
  Eigenvalues: Not applicable (matrix is not square)

Matrix: Ac = Ax1 - Ay1 * ( Ay2 \ Ax2 )
  Rank: 7
  Size: (7, 7)
  Condition Number: 16797.777105168134
  Eigenvalues: [-0.7322, -0.6603, -0.5000, -0.0001, 1.0000, 1.0000, 1.0000]

Matrix: Bc = B1 - Ay1 * ( Ay2 \ B2 )
  Rank: 2
  Size: (7, 2)
  Condition Number: 706.4701174170237
  Eigenvalues: Not applicable (matrix is not square)

Matrix: Cc = E1 - Ay1 * ( Ay2 \ E2 )
  Rank: 2
  Size: (7, 2)
  Condition Number: 2.1139228150591085e15
  Eigenvalues: Not applicable (matrix is not square)

i am not really sure anymore if that is relevant information.

13

Did you try moving to mass matrix form FBDF and trying FBDF(linsolve = QRFactorization(LinearAlgebra.ColumnNorm())) ?