Large DAE problem in Julia using DifferentialEquations

New to Julia
, ,
1

Hi,
I am currently working on a project which requires the solution of a large-scale DAE problem (differential algebraic equations). I have tried some example test codes such as Robertson equation which has worked fine (Link: Implicitly-Defined Differential-Algebraic Equations (DAEs)).

However, my problem is on a larger scale with 137 variables and equations. In order to find suitable initial conditions, I used JuMP optimization on a steady-state model of the problem/project and set the differential variables du0 as all zeros.

I cannot share the entire code due to the confidentiality of the project. So I made a small working example. (The “par” variables are defined in a separate dictionary, and are set as constants.) The example’s steady-state model:

# User packages 
using JuMP, Ipopt, Plots, PGFPlotsX, LaTeXStrings

include("Parameters.jl") # parmeters

# tolerance of Ipopt
tol  = 1e-8; #1e-8; # sjekk tolerance
iter = 8000;

function SS_Model(par, tol, iter)
    SSEl = Model(Ipopt.Optimizer)
    set_optimizer_attribute(SSEl, "max_iter", iter)
    set_optimizer_attribute(SSEl, "tol", tol)

    ############# Variables ###############################
    #input
    @variable(SSEl, 0 <= p1 <= 1, start = 0.1);
    @variable(SSEl, 0 <= p2 <= 1, start = 0.1);

    @variable(SSEl, 0.5 <= V_a_rec <= 1000, start = 6);
    @variable(SSEl, 0.5 <= V_c_rec <= 1000, start = 6);
    @variable(SSEl, 0 <= n_a_out, start = 1);
    @variable(SSEl, 0 <= n_c_out, start = 1);

    ############# Expressions #############################

    @NLexpression(SSEl, n_c_rec, (V_c_rec*10^3)/22.4); # m^3 to mol ideal gas
    @NLexpression(SSEl, n_a_rec, (V_a_rec*10^3)/22.4);

    ############# Model (written as constraints) ##########

    @NLconstraint(SSEl, n_a_rec - 100*p1*sqrt(par[:psto] - par[:psep]) == 0);
    @NLconstraint(SSEl, n_c_rec - 100*p2*sqrt(par[:psto] - par[:psep]) == 0);
    
    @NLconstraint(SSEl, n_a_rec - n_a_out == 0);
    @NLconstraint(SSEl, n_c_rec - n_c_out == 0);

    #objective
    @NLobjective(SSEl, Max, n_a_rec + n_c_rec);
    return SSEl
end

SSEl = SS_Model(par, tol, iter)

optimize!(SSEl)

ssx = value.(all_variables(SSEl))
ssv = name.(all_variables(SSEl))

dict = Dict(ssv[1] => ssx[1];)

for i in 2:length(ssx)
    push!(dict, ssv[i] => ssx[i]);
end

p = ssx[1:2]
x0 = ssx[3:6]
dx0 = zeros(4);

The example’s DAE model:

# User packages 
using Sundials
import DifferentialEquations

include("ssmodel.jl")

function f(out, dx, x, p, t)

    #Variables#########
    V_a_rec = x[1];
    V_c_rec = x[2];
    n_a_out = x[3];
    n_c_out = x[4];
    
    d_na = dx[1];
    d_nc = dx[2];

    #Expressions#######
    n_a_rec = (V_a_rec*10^3)/22.4;
    n_c_rec = (V_c_rec*10^3)/22.4;

    #Model#############
    out[1] = n_a_rec - 100*p[1]*sqrt(par[:psto] - par[:psep]);
    out[2] = n_c_rec - 100*p[2]*sqrt(par[:psto] - par[:psep]);

    out[3] = n_a_rec - n_a_out - d_na; 
    out[4] = n_c_rec - n_c_out - d_nc;
    
end

tspan = (0.0, 1.0);

differential_vars = [false, false, true, true]

prob = DAEProblem(f, dx0, x0, tspan, p, differential_vars = differential_vars)
sol = solve(prob, IDA())

The problem that occurs when I try to solve the large-scale problem that I have formulated is:

  1. DomainErrors where the solver sets the states to negative variables.
  2. I get the error message:

[IDAS ERROR] IDACalcIC
The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover.

retcode: InitialFailure

Now, I am wondering if someone has had a similar problem. Additionally, I have heard that IDA or DifferentialEquations.jl might not be able to solve problems on such a large scale. I would really appreciate all answers. Thanks in advance!

2

if recommend trying dfbdf. in my experience it is a lot more stable. you also might want to try formulating this as a mass matrix problem.

3

From whom :sweat_smile: . It’s being used in some examples with about 10 million equations on a few thousand cores that I know of, so it’s at least tested to a decent scale. Just choose the right method for the right size.

Your algebraic equations are not a function of the algebraic variables, so your DAE is not Index 1. IDA requires Index 1 DAEs. I recommend using ModelingToolkit to handle this kind of thing better. See:

for an example, and see if that changes the generated equations.

1 Like