Mathematica NDSolve for DAE and System Modeling

Question

I'm a student for HVAC engineering and I've been using Mathematica for my study (mostly differential-algebraic equation system or Finite-Element-Analysis) and my hobby project for years. I really enjoy Mathematica, but I always have been struggling with NDSolve and its DAE solver, because it seems to be not very robust (see this post).

I have this feeling that NDSolve (or IDA) does not particularly favor discontinuous functions such as Clip, Min, Max. For most engineering problems and modeling method, I think that these functions are not avoidable. I've been browsing this forum for solutions quite a long time and there appeared some solutions for it, but they are mostly only for some specific case and not applicable for my case (e.g. a simple P-controller with a minimum of a close-to-zero value such as $MachineEpsilon will very likely trigger the NDSolve::ndsz error for unknown reasons, not sure if it's a bug or not).

Till recent, I found that since v11 Mathematica allows connection to Modelica with CreateSystemModel. I've heard that Modelica has a strong DAE solver and I would really like to test this Mathematica-Modelica feature. But the documentation for System Modeling seems to be not very detailed. For the first trial I would like to realize a simple house heating model from this post.

(* 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 = Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]
eq = {
  (* equation for building *)
  Q[tau] - 200 (tt[tau] - ta[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] - Max[Min[(20 + 20 (20 - tt[tau])), 100], 10^(-10)] == 0}

Here is what I've got so far.

model = CreateSystemModel[eq, tau]
data = CreateDataSystemModel[tae]
wholeModel = 
 ConnectSystemModelComponents[{model, data}, {"ta" \[Element] data}]

The connection seems to be unsuccessful and SystemModelSimulate returns nothing. It would be really appreciated if a simple example can be provided for the new System Model feature.

EDIT

Here is a simplified version of the above problem. It is about the thermal behavior of a house. The thermal resistance is R (e.g. $1/200 \frac{K}{W}$) and it has a heat capacity of Cwirk (e.g. $50 \cdot 25 \cdot 3 \frac{W h}{K}$). The variable t is the room temperature and ta is the outdoor ambient temperature. Q stands for the heating output of the in-house heating device (we can set Q to a constant such as $1000 W$). So the conservation equation is

$$ C_{Wirk} \cdot \frac{dt}{d\tau} = - \frac{t - t_a}{R} + \dot{Q} $$

The t_a is obviously an input to this model. We can use the weather data from the list tae.

li = Import[
   "http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/data/tmy3/\
725958TYA.CSV"];
(* interpolate outdoor temperature *)
tae = Transpose[{Range[8760], Drop[Drop[li, 1][[All, 32]], 1]}]

It's unclear for me how to connect the weather data to the equation/model. It would be very appreciated if any explanation or maybe an example can be provided.

Answer 1

Based on your simplified version, and using hard-coded data, create the model with the data source:

data = {{0, 6}, {1, 9}, {2, 1}, {3, 9}, {4, 7}, {5, 6}, {6, 1}, {7, 
    8}, {8, 0}, {9, 2}, {10, 9}};
datamodel = CreateDataSystemModel[data]

Define your parameters:

cwirk = 50 25 3;
q = 1000;
r = 1/200;

Create the model containing the equations. Note that we are defining the ta as a Real input:

eqmodel = 
 CreateSystemModel[{cwirk t'[tau] == (t[tau] - ta[tau])/r + q}, 
  tau, {ta \[Element] "RealInput"}]

Connect them. Note how we connect the output from the data model, data.y[1], to the input we created in the equation model, eqs.ta:

simmodel = ConnectSystemModelComponents[
   {"data" \[Element] datamodel, "eqs" \[Element] eqmodel},
   {"data.y[1]" -> "eqs.ta"}
]

And then you can simulate and plot:

sim = SystemModelSimulate[simmodel, 10]
SystemModelPlot[sim, {"eqs.t"}]