# Large DAE problem in Julia using DifferentialEquations
**URL:** https://discourse.julialang.org/t/large-dae-problem-in-julia-using-differentialequations/97555
**Category:** New to Julia
**Tags:** jump, differentialequation, sundials
**Created:** [April 16, 2023, 6:43pm UTC](https://discourse.julialang.org/t/large-dae-problem-in-julia-using-differentialequations/97555 "2023-04-16T18:43:10Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@hannemgr](https://discourse.julialang.org/u/hannemgr)
#### Post date: [April 16, 2023, 6:43pm UTC](https://discourse.julialang.org/t/large-dae-problem-in-julia-using-differentialequations/97555/1 "2023-04-16T18:43:10Z")
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)](https://docs.sciml.ai/DiffEqDocs/latest/tutorials/dae_example/)).
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:
```julia
# 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:
```julia
# 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!
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [April 16, 2023, 7:19pm UTC](https://discourse.julialang.org/t/large-dae-problem-in-julia-using-differentialequations/97555/2 "2023-04-16T19:19:09Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 16, 2023, 9:07pm UTC](https://discourse.julialang.org/t/large-dae-problem-in-julia-using-differentialequations/97555/3 "2023-04-16T21:07:39Z")
> [@hannemgr](#):
>
> 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
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.
> [@hannemgr](#):
>
> ```julia
> 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;
>
> ```
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:
> **[Automated Index Reduction of DAEs · ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/examples/modelingtoolkitize_index_reduction/)**
>
> Documentation for ModelingToolkit.jl.
for an example, and see if that changes the generated equations.