Speed up NDSolve for DAE

Question

I'm handling with a large DAE system with some events (around 400 equations and 20 events) and lots of them are nonlinear. Let's take the equations from this post as an example.

(* effective heat capacity of building *)
Cwirk = 50 25 3;
(* import outdoor temperature *)
li = Import[
   "http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/\
725958TYA.CSV"];
(* interpolate outdoor temperature *)
tae = Interpolation[
  Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]]
eq = {
  (* equation for building *)
  Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], 
  (* heating capacity of floor heating system *)
  Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, 
  (* the water outlet temperature of floor heating *)
  tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 vF[tau] 0.22)] == 0,
  (* a simple P-controller for the flow rate *) 
  vF[tau] - Clip[(20 + 20 (20 - tt[tau])), {10^(-10), 100}] == 0}
ic = FindRoot[eq[[All, 1]] == 0 /. tau -> 1, 
   Transpose[{{Q[1], tt[1], tr[1], vF[1]}, {3000, 20, 20, 20}}]] /. 
  Rule -> Equal
workingPrecision = 8;
AbsoluteTiming[
 sol = NDSolve[
    SetPrecision[Join[eq, ic], workingPrecision], {Q, tt, tr, 
     vF}, {tau, 1, 24 365}, WorkingPrecision -> workingPrecision, 
    MaxSteps -> Infinity];]
(*{31.2564,Null}*)

Question 1:

Solving these 4 equations already takes around 30 seconds for me. For the full 400 equations, my poor laptop needs about 3 hours, which is little bit too long. Is there any way to accelerate the NDSolve for DAE or at least this example?

Question 2:

I've found some posts in this forum state that vectorize the equations will speed up the NDSolve. Can DAEs also be written in vector form and will it help for the speed?

As to question 2, vectorizing speeds up NDSolve mainly because it speeds up the pre-processing stage of NDSolve, AFAIK. (Related: https://mathematica.stackexchange.com/a/158519/1871) So personally I don't think that'll help for your problem. — xzczd, Nov 25 '19 at 06:01
@xzczd thanks for the comment. What if I compile all the equations? Will it help? — 407PZ, Nov 25 '19 at 09:10
I doubt, because reading the document of e.g. StateSpace method, the pre-processing of DAE solver seems to heavily depend on… Er… structure analysis of DAE system. But I confess I never tried compiling DAE system, perhaps it's worth exploring. — xzczd, Nov 25 '19 at 09:39
You guys are gonna love this ... I spent like an hour trying to run "NDSolve" in parallel by breaking it up into six different ranges, thinking it would be faster to solve a bunch of small ranges rather than a large range. I finally realized that all the kernels were solving from the initial condition so there was no improvement in speed. Moral of the story: I don't think this is the ideal candidate for a parallel computation solution. — NeilTheSeal, Nov 26 '19 at 06:22
@NeilTheSeal It's not that surprising considering how the numeric DAE solver works. (It has to start from the initial value, no matter how you set the domain of definition. ) Another possibility for parallelization is to use the option mentioned here, but still, I doubt if this will help, because OP's system only involves 400 equations, which is still small. — xzczd, Nov 27 '19 at 04:02

score 5 · Answer 1 · answered Nov 26 '19 at 05:53

For this set of equations, at least, eliminating the algebraic variables is straightforward and decreases the runtime by a factor of three.

Cwirk = 50 25 3;
li = Import["http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/\725958TYA.CSV"];
tae = Interpolation[Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]];
eq = 100 (28 - tt[tau] - (28 - tt[tau]) 
  Exp[-0.9/(7/6 Clip[(20 + 20 (20 - tt[tau])), {10^(-10), 100}] 0.22)]) 
  Clip[(20 + 20 (20 - tt[tau])), {10^(-10), 100}] 7/6 - 200 (tt[tau] - tae[tau]) 
  == Cwirk tt'[tau];
ic = FindRoot[eq[[1]] == 0 /. tau -> 1, {tt[1], 20}] /. Rule -> Equal;
AbsoluteTiming[sol = NDSolveValue[{eq, ic}, tt, {tau, 1, 24 365}, MaxSteps -> Infinity];]

which has a runtime of about 6.7 seconds, as compared with 21 seconds on my computer for the code in the question. The algebraic variables then can be computed by back substitution, if desired.

Even 6.7 seconds may seem large. It is so, because the source term. tae, is very noisy. To the extent that it can be smoothed, further runtime improvements can be expected. Smoothing the source term also would reduce the vulnerability of the numerical computation to crashing when small details, like WorkingPrecision, are changed modestly.

For completeness, the solution can be plotted as

Plot[sol[tau], {tau, 1, 24 365}, PlotRange -> All, AxesLabel -> {tau, "tt"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, 15}]

which appears to be the same as that in the post cited in the question.

Speed up NDSolve for DAE

1 Answers1