3 problems in one model - sym. differentiation, code gen. and simulation

[jakesson]

Problem present in omc 1.4.1

The model below results in several errors

1) The index-reduction algorithm fails since the differentiation rule for the
log function (and several others) is missing in the function
Derive.differentiateExpTime. Adding the corresponding rules enables OpenModelica
to generate simulation code.

2) The simulation code does not compile. The problem is resolved by adding a
line to modelica.h:

typedef modelica_real log_rettype;

It should be checked whether more typedefs are missing. With this modification
of modelica.h, the model compiles

3) The model is not correctly simulated. Some variables are ok, e.g. Mt, whereas
others are constant, e.g. p, which is not correct.

In the model, there are differential equations of the type

Es = rhos*us*Vs;
us = hs - p/rhos;
der(Es) = qs*hs - qc*hs;

Clearly, all factors defining Es are time varying. Is this a problem?

model CTest

parameter Real V=5;

parameter Real m=1130;
parameter Real A=15.5;
parameter Real alpha=2000;
parameter Real Cp=500;
Real p( start=3.80e5);
Real Tm(start=135,stateSelect=StateSelect.always);
Real qs;
Real qw;
Real Qp;
Real Qm;
Real Vs;
Real Vw;
Real rhos;
Real rhow;
Real Ms;
Real Mw(start=114);
Real Mt(start=120);
Real Es(start=5e7);
Real Ew;
Real Em(start=5.2e8);
Real qc;
Real Ts;
Real us;
Real uw;
Real hs;
Real hw;

equation

0.25=qs;
-0.3=qw;
10000=Qp;

Ms = rhos*Vs;

Mw = rhow*Vw;
der(Mt) = qs + qw;
Mt = Ms + Mw;
Es = rhos*us*Vs;
us = hs - p/rhos;
der(Es) = qs*hs - qc*hs;
Ew = rhow*uw*Vw;
uw = hw - p/rhow;
der(Ew) = qc*hs + qw*hw - Qm;
Em = m*Cp*Tm;
der(Em) = Qm + Qp;
Qm = alpha*A*(Ts - Tm);
V = Vs + Vw;
Ts = 0.1723*(log(p))^{3 - 3.388*(log(p))}2 +

37.71*log(p) - 148.7;

hs = 1e3*(-0.07402*(log(p))^{4 + 2.887*(log(p))}

3 - 39.58*(log(p))^{2 + 260*log(p) + 1824);}

hw = 1e3*(0.8842*(log(p))^{3 - 18.77*(log(p))}2

+ 200*log(p) - 748.5);

rhos = (0.005048*p + 64.26)/1000;
rhow = -0.3136*(log(p))^{3 + 6.792*(log(p))}2 -

52.43*log(p) + 1141;

annotation (uses(Modelica(version="2.2")));

end CTest;

[petar]

Johan, could you please sent me the changes you made to Derive.mo and I will
look at it and check it in to the svn. /Peter

[petar]

I've fixed the state selection problem. Now Tm is selected as state even though
it does not appear differentiated. But it still does not simulate, since the
non-linear equation can not be solved. Perhaps the analytic jacobian is required
for this, I am not sure.
Johan, do you know if the jacobian is required?

[jfrenkel]

The model now simulates, comparison with dymola result in no differences, also state selection is equal, the output is

No. of Equations: 22
No. of Variables: 22
##########################################################
Statistics
##########################################################
Number of independent Subsystems: 1
selected States: p, Tm, Vw
Number of States:                 3
Toplevel Inputs:                  0

Single Equations:  17
Array Equations:   0
Algorithms:        0
Complex Equations: 0

Equationsystems with constant Jacobian:     0 {}
Equationsystems with time varying Jacobian: 1 {14}
Equationsystems with nonlinear Jacobian:    0 {}
Equationsystems without analytic Jacobian:  0 {}

mixed Equationsystems with Single Equation:       0 {}
mixed Equationsystems with Array Equation:        0 {}
mixed Equationsystems with Algorithm:             0 {}
mixed Equationsystems with Complex Equation:      0 {}
mixed Equationsystems with constant Jacobian:     0 {}
mixed Equationsystems with time varying Jacobian: 0 {}
mixed Equationsystems with nonlinear Jacobian:    0 {}
mixed Equationsystems without analytic Jacobian:  0 {}
mixed Equationsystems without analytic Jacobian:  0 {}
mixed Equationsystems with linear Tearing System:  0 {}
mixed Equationsystems with nonlinear Tearing System:  0 {}
torn linear Equationsystems:    0 {}
torn nonlinear Equationsystems:  0 {}
##########################################################
SingleEquation: 14
1:  Ts:VARIABLE() .CTest, .Real type: Real 
1/1 (1): Ts = -148.7 + 0.1723 * log(p) ^ 3.0 + -3.388 * log(p) ^ 2.0 + 37.71 * log(p)
SingleEquation: 15
1:  Qm:VARIABLE() .CTest, .Real type: Real 
1/1 (1): Qm = alpha * A * (Ts - Tm)
SingleEquation: 16
1:  $DER.Em:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
1/1 (1): $DER.Em = 10000.0 + Qm
SingleEquation: 17
1:  Tm:STATE(1)(start = 135.0 stateSelect=StateSelect.alwas ) .CTest, .Real type: Real 
1/1 (1): $DER.Em = m * Cp * der(Tm)
SingleEquation: 6
1:  Em:DUMMY_STATE(start = 520000000.0 ) .CTest, .Real type: Real 
1/1 (1): Em = m * Cp * Tm
SingleEquation: 5
1:  Vs:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): V = Vs + Vw
SingleEquation: 4
1:  hs:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): hs = 1824000.0 + 1000.0 * (-0.07402 * log(p) ^ 4.0 + 2.887 * log(p) ^ 3.0 + -39.58 * log(p) ^ 2.0 + 260.0 * log(p))
SingleEquation: 3
1:  hw:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): hw = -748500.0 + 1000.0 * (0.8842 * log(p) ^ 3.0 + -18.77 * log(p) ^ 2.0 + 200.0 * log(p))
SingleEquation: 2
1:  rhos:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): rhos = 0.06426000000000001 + 0.000005048000000000001 * p
SingleEquation: 9
1:  us:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): us = hs + (-p) / rhos
SingleEquation: 10
1:  Es:DUMMY_STATE(start = 50000000.0 ) .CTest, .Real type: Real 
1/1 (1): Es = rhos * us * Vs
SingleEquation: 12
1:  Ms:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): Ms = rhos * Vs
SingleEquation: 1
1:  rhow:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): rhow = 1141.0 + -0.3136 * log(p) ^ 3.0 + 6.792 * log(p) ^ 2.0 + -52.43 * log(p)
SingleEquation: 7
1:  uw:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): uw = hw + (-p) / rhow
SingleEquation: 8
1:  Ew:DUMMY_STATE() .CTest, .Real type: Real 
1/1 (1): Ew = rhow * uw * Vw
SingleEquation: 11
1:  Mw:DUMMY_STATE(start = 114.0 ) .CTest, .Real type: Real 
1/1 (1): Mw = rhow * Vw
SingleEquation: 13
1:  Mt:DUMMY_STATE(start = 120.0 ) .CTest, .Real type: Real 
1/1 (1): Mt = Ms + Mw
Equationsystem Jacobian Time varying:
1:  $DER.Ms:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
2:  $DER.Mw:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
3:  $DER.rhow:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
4:  $DER.us:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
5:  $DER.hs:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
6:  $DER.uw:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
7:  $DER.hw:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
8:  p:STATE(1)(start = 380000.0 ) .CTest, .Real type: Real 
9:  $DER.rhos:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
10:  $DER.Ew:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
11:  qc:VARIABLE() .CTest, .Real type: Real 
12:  $DER.Es:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
13:  $DER.Vs:DUMMY_DER(fixed = false ) .CTest, .Real type: Real 
14:  Vw:STATE(1)() .CTest, .Real type: Real 
1/1 (1): $DER.Ew = rhow * (uw * der(Vw) + $DER.uw * Vw) + $DER.rhow * uw * Vw
2/2 (1): 0.0 = $DER.Vs + der(Vw)
3/3 (1): $DER.Es = rhos * (us * $DER.Vs + $DER.us * Vs) + $DER.rhos * us * Vs
4/4 (1): $DER.Es = 0.25 * hs - qc * hs
5/5 (1): $DER.Ew = qc * hs + -0.3 * hw - Qm
6/6 (1): $DER.Ms = rhos * $DER.Vs + $DER.rhos * Vs
7/7 (1): $DER.rhos = 0.000005048000000000001 * der(p)
8/8 (1): $DER.hw = 1000.0 * (2.6526 * log(p) ^ 2.0 * der(p) / p + -37.54 * log(p) * der(p) / p + 200.0 * der(p) / p)
9/9 (1): $DER.uw = $DER.hw + ((-der(p)) * rhow - (-p) * $DER.rhow) / rhow ^ 2.0
10/10 (1): $DER.hs = 1000.0 * (-0.29608 * log(p) ^ 3.0 * der(p) / p + 8.661 * log(p) ^ 2.0 * der(p) / p + -79.16 * log(p) * der(p) / p + 260.0 * der(p) / p)
11/11 (1): $DER.us = $DER.hs + ((-der(p)) * rhos - (-p) * $DER.rhos) / rhos ^ 2.0
12/12 (1): $DER.rhow = -0.9408 * log(p) ^ 2.0 * der(p) / p + 13.584 * log(p) * der(p) / p + -52.43 * der(p) / p
13/13 (1): $DER.Mw = rhow * der(Vw) + $DER.rhow * Vw
14/14 (1): -0.04999999999999999 = $DER.Ms + $DER.Mw
Jac:
{1,3}:(-uw) * Vw,
{1,6}:(-rhow) * Vw,
{1,10}:1.0,
{1,14}:(-rhow) * uw,
{2,13}:-1.0,
{2,14}:-1.0,
{3,4}:(-rhos) * Vs,
{3,9}:(-us) * Vs,
{3,12}:1.0,
{3,13}:(-rhos) * us,
{4,11}:hs,
{4,12}:1.0,
{5,10}:1.0,
{5,11}:-hs,
{6,1}:1.0,
{6,9}:-Vs,
{6,13}:-rhos,
{7,8}:-0.000005048000000000001,
{7,9}:1.0,
{8,7}:1.0,
{8,8}:-1000.0 * (2.6526 * log(p) ^ 2.0 / p + -37.54 * log(p) / p + 200.0 / p),
{9,3}:(-p) / rhow ^ 2.0,
{9,6}:1.0,
{9,7}:-1.0,
{9,8}:1.0 / rhow,
{10,5}:1.0,
{10,8}:-1000.0 * (-0.29608 * log(p) ^ 3.0 / p + 8.661 * log(p) ^ 2.0 / p + -79.16 * log(p) / p + 260.0 / p),
{11,4}:1.0,
{11,5}:-1.0,
{11,8}:1.0 / rhos,
{11,9}:(-p) / rhos ^ 2.0,
{12,3}:1.0,
{12,8}:0.9408 * log(p) ^ 2.0 / p + -13.584 * log(p) / p + 52.43 / p,
{13,2}:1.0,
{13,3}:-Vw,
{13,14}:-rhow,
{14,1}:-1.0,
{14,2}:-1.0