5
$\begingroup$

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?

Improve this question
$\endgroup$
5
  • $\begingroup$ As to question 2, vectorizing speeds up NDSolve mainly because it speeds up the pre-processing stage of NDSolve, AFAIK. (Related: mathematica.stackexchange.com/a/158519/1871) So personally I don't think that'll help for your problem. $\endgroup$
    – xzczd ♦
    Commented Nov 25, 2019 at 6:01
  • $\begingroup$ @xzczd thanks for the comment. What if I compile all the equations? Will it help? $\endgroup$
    – 407PZ
    Commented Nov 25, 2019 at 9:10
  • $\begingroup$ 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. $\endgroup$
    – xzczd ♦
    Commented Nov 25, 2019 at 9:39
  • 1
    $\begingroup$ 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. $\endgroup$
    – NeilTheSeal
    Commented Nov 26, 2019 at 6:22
  • 1
    $\begingroup$ @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. $\endgroup$
    – xzczd ♦
    Commented Nov 27, 2019 at 4:02

1 Answer 1

Reset to default
5
$\begingroup$

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}]

enter image description here

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.