# Large DAE problem in Julia using DifferentialEquations

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!

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.

From whom . 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:

https://docs.sciml.ai/ModelingToolkit/stable/examples/modelingtoolkitize_index_reduction/

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

1 Like