NDSolve DAE with Constraints

The 10^-10 in Clip is indeed causing trouble, but not in Exp[-0.9/(7/6 vF[tau] 0.22)], because the mconly warning still pops up even after the technique mentioned in this post is used.

So, what's the root of the problem? I'm not sure, nevertheless, the problem can be resolved by changing 10^-10 to something a bit larger:

eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], 
   Q[tau] == 100 (28 - tr[tau]) vF[tau] 7/6, 
   tr[tau] == tt[tau] + (28 - tt[tau]) (10^-10 + Exp[-0.9/(7/6 vF[tau] 0.22)]), 
   vF[tau] - Clip[(20 + 20 (20 - tt[tau])), {0.1(*See here!*), 100}] == 0, 
   tt[1] == 15, tr[1] == 15, 
   vF[1] == 100, Q[1] == 0};
sol0 = NDSolveValue[eq, tt, {tau, 0, 24}, MaxSteps -> Infinity]

 {sol0, sol2}[24] // Through
 (* {20.2045, 20.2045} *)

The result at tau == 24 is the same as that of sol2, so I think this solution is acceptable?

The block of code in the question following "However, if I eliminate some equations with Rule, it runs smoothly" works, because the transformation changes the differential-algebraic system to a purely differential system, for which NDSolve uses more robust solution methods. The same can be accomplished more generally by expressing the variables appearing only algebraically as functions of tt. Here, the functions are

Qf[tt_?NumericQ] := 100 (28 - trf[tt]) vFf[tt] 7/6;
trf[tt_?NumericQ] := tt + (28 - tt) Exp[-(270/77)/vFf[tt]];
vFf[tt_?NumericQ] := Clip[(20 + 20 (20 - tt)), {10^-10, 100}];

and the system, now a single ODE, is solved by

Cwirk = 50 25 3;
eq = {Qf[tt[tau]] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau], tt[1] == 15};
sol1 = Flatten[NDSolveValue[eq, {tt}, {tau, 0, 24}]];
ul = sol1[[1]]["Domain"][[1, 2]];
Plot[sol1[[1]][tau], {tau, 0, ul}, AxesLabel -> {tau, tt}, 
    LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]

In general, the variables appearing only algebraically can be converted into functions, although some care must be taken, if they are not single-valued.

As an aside, the original system of equations in the question, when rewritten as

Cwirk = 50 25 3;
eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau],
      vF[tau] == Clip[(20 + 20 (20 - tt[tau])), {3 10^-3, 100}] , 
      tr[tau] == tt[tau] + (28 - tt[tau]) Exp[-0.9/(7/6 vF[tau] 0.22)],
      Q[tau] == 100 (28 - tr[tau]) vF[tau] 7/6,
      tt[1] == 15, tr[1] == 15, vF[1] == 100, Q[1] == 0};
sol0 = Flatten[NDSolveValue[eq, {tt, tr, vF, Q}, {tau, 0, 24}]];

also runs to completion, because the lower limit of Clip now is set to 3 10^-3, which I would think would be good enough for most purposes. (The exponential in the question can assume a minimum value as small as about 10^-305.) Incidentally, the warning messages, NDSolveValue::ivres returned by this latter block of code indicate that the NDSolve has corrected the initial conditions for consistency:

Table[sol0[[i]][1], {i, 4}]
(* {15., 27.5521, 100., 5226.02} *)

Using these values instead of those in the question eliminates the warning messages.

Addendum: Result for imported li

The first block of code in this answer also can handle the imported version of li, despite the fact that the resulting tae is very noisy.

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]}]];
Plot[tae[tau], {tau, 1, 24 365}, AxesLabel -> {tau, "tae"}, PlotRange -> All, 
    LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]

Then, the first block of code, integrated over {tau, 1, 24 365}, yields

Incidentally, the second block of code with the lower bound of Clip set to 5 10^-2 reproduces the final plot above.

I found out that setting the WorkingPrecision to 8 simply solves this problem.

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

Cwirk = 50 25 3;
eq = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau],
      Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0, 
      tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 vF[tau] 0.22)] == 0,
      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} *)

But somehow the Rule method doesn't work anymore on my MMA 12.0 mac.

rs = {Q -> 100 (28 - tr) vF 7/6, 
   tr -> tt[tau] + (28 - tt[tau]) Exp[-0.9/(7/6 vF 0.22)], 
   vF -> Clip[20 + 20 (20 - tt[tau]), {10^(-10), 100}]};
eq2 = Q - 200 (tt[tau] - tae[tau]) == tt'[tau];
AbsoluteTiming[
 sol2 = Flatten[
   NDSolveValue[{eq2 //. rs, tt[1] == 18.887256814210946}, 
    tt, {tau, 1, 24 365}, MaxSteps -> Infinity]]]

NDSolveValue::ndcf: Repeated convergence test failure at tau == 3497.9989465540943`; unable to continue

EDIT

Here is another solution for this case. I always thought that NDSolve does not particularly like discontinuous functions such as Clip, Min, Max. It is possible to replace them with some events. Below a discrete variable mode is used to indicate if the function is below zero or above 100. If below zero, then the mode will be set to -1. If above 100, it will be set to 1. Otherwise, it is 0. Using the mode variable, Clip can be avoided.

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

eq1 = {
  Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau],  
  Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0,
  tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/7/6 vF[tau] 0.22] == 0,
  vF[tau] - DiscreteDelta[mode[tau] - 1] 100 - 
    DiscreteDelta[mode[tau]] (20 + 20 (20 - tt[tau])) - 
    DiscreteDelta[mode[tau] + 1] $MachineEpsilon == 0}
ic1 = Join[
   FindRoot[eq1[[All, 1]] == 0 /. {tau -> 1, mode[_] -> 0}, 
    Transpose[{{Q[1], tt[1], tr[1], vF[1]}, {3000, 20, 20, 
       20}}]], {mode[1] -> 0}] /. Rule -> Equal
evts1 = {WhenEvent[(20 + 20 (20 - tt[tau])) < $MachineEpsilon, 
   mode[tau] -> -1], 
  WhenEvent[(20 + 20 (20 - tt[tau])) > 100, mode[tau] -> 1], 
  WhenEvent[(20 + 20 (20 - tt[tau])) < 100, mode[tau] -> 0], 
  WhenEvent[(20 + 20 (20 - tt[tau])) > $MachineEpsilon, 
   mode[tau] -> 0]}
AbsoluteTiming[
 sol1 = NDSolve[
     Join[eq1, evts1, ic1], {Q, tt, tr, vF, mode}, {tau, 1, 24 365}, 
     DiscreteVariables -> mode][[1]];]
(* {73.7385, Null} *)

It's worth mentioning that this method is 2x slower than simply setting WorkingPrecision to 8.

I would like to share some new findings to this post after nearly 2 years.

As described in the question, running the original code gives me the NDSolve::mconly error.

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]}]]
Cwirk = 50 25 3;
eq0 = {Q[tau] - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau],
       Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0,
       tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 (vF[tau]) 0.22)] == 0,
       vF[tau] - Clip[20 + 20 (20 - tt[tau]), {$MachineEpsilon, 100}] == 0}
ic0 = FindRoot[eq0[[All, 1]] == 0 /. tau -> 1, 
               Transpose[{{Q[1], tt[1], tr[1], vF[1]}, {3000, 20, 20, 20}}], 
               MaxIterations -> 10000] /. Rule -> Equal

sol0 = Monitor[NDSolve[Rationalize[Join[eq0, ic0]], {Q, tt, tr, vF}, 
                       {tau, 1, 24 365}, 
                       EvaluationMonitor :> (monitor = Row[{"t = ", tau}])], 
               monitor][[1]]
(* NDSolve::mconly error *)

But if I modify the third equation tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 (vF[tau]) 0.22)] == 0 into tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 Max[$MachineEpsilon, vF[tau]] 0.22)] == 0, NDSolve returns a different error message NDSolve::ndsz. This is quite interesting. My guess is that NDSolve thought it might be a good idea to try vF[tau] == 0 in the equation at some step which might result in NDSolve::mconly error for IDA solver.

The next thing I tried is to eliminate the second equation Q[tau] - 100 (28 - tr[tau]) vF[tau] 7/6 == 0 and surprisingly it works.

eq2 = {100 (28 - tr[tau]) vF[tau] 7/6 - 200 (tt[tau] - tae[tau]) == Cwirk tt'[tau],
       tr[tau] - tt[tau] - (28 - tt[tau]) Exp[-0.9/(7/6 Max[vF[tau], $MachineEpsilon] 0.22)] == 0,
       vF[tau] - Clip[20 + 20 (20 - tt[tau]), {$MachineEpsilon, 100}] == 0}
ic2 = FindRoot[eq2[[All, 1]] == 0 /. tau -> 1, 
               Transpose[{{tt[1], tr[1], vF[1]}, {20, 20, 20}}], 
               MaxIterations -> 10000] /. Rule -> Equal
AbsoluteTiming[sol2 = Monitor[NDSolve[Rationalize[Join[eq2, ic2]], {tt, tr, vF}, {tau, 1, 24 365}, EvaluationMonitor :> (monitor = Row[{"t = ", tau}])], monitor][[1]];]
(* 64.5397 *)

This leads me to take a deeper look into the second equation. Actually for the equation solving, the second equation is not necessary. We can solve the tt and calculate the Q afterwards. After some brief investigation, it turns out that the second equation is quite stiff. If we insert the forth equation and the third equation in to the second, we will have a function $Q(tt(tau))$. If we plot the function $Q(tt(tau))$, we can easily see where the problem is.

rl1 = {Q[tau] -> 100 (28 - tr[tau]) vF[tau] 7/6, 
       tr[tau] -> tt[tau] + (28 - tt[tau]) Exp[-0.9/(7/6 Max[$MachineEpsilon, vF[tau]] 0.22)],
       vF[tau] -> Clip[20 + 20 (20 - tt[tau]), {$MachineEpsilon, 100}]}
Plot[(rl1[[1, 2]] //. Drop[rl1, 1]) /. {tt[tau] -> x}, {x, 15, 25}, AxesLabel -> {"tt", "Q"}]

The curve becomes very steep when the $tt$ increases to 21 and after 21 the curve becomes flat. If we extract the values for the last steps in sol0, we will see that the last values for $tt$ is around 21.

After all these back and forth and trying for really couple years, I come to a conclusion concerning using NDSolve, which I want to tell myself repeatedly and share with everyone:

Always check the equation first before using NDSolve, even they have physical meaning.