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:

- DomainErrors where the solver sets the states to negative variables.
- 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!