Bug in new mover implementation

[Mathadon]

@rubenbaetens showed me a bug where the new mover implementation is causing problems.

When prescribing pressure and dp_in=0 the non-linear solver is sometimes unable to solve r_N. This is probably caused by a singularity in the similarity law equations for dp=0 and m_flow=0. Will need to debug this further later. A possible work-around is not to set dp_in=0 but instead dp_in = small value or to disable the similarity law computation (however the parameter for this has a final value right now).

[mwetter]

@rubenbaetens @Mathadon Do you have a small test case that reproduces the problem. This will be needed to verify that there is singularity that the tools cannot handle, and to ensure that we corrected it.

[Mathadon]

I tried to reproduce a similar bug that I hot-fixed earlier but it does not re-appear. I will let you know if it re-appears. @rubenbaetens do you still have this problem?

[rubenbaetens]

i'll try to reproduce it, but i solved it by applying a m_flow_small i think so the bug doesn't appear in my models anymore and it might take some effort trying to reproduce it (in case it wasn't caused by other mistakes in my model).

On Mon, May 9, 2016 at 12:26 AM, Filip notifications@github.com wrote:

I tried to reproduce a similar bug that I hot-fixed earlier but it does not re-appear. I will let you know if it re-appears. @rubenbaetens https://github.com/rubenbaetens do you still have this problem?

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/iea-annex60/modelica-annex60/issues/458#issuecomment-217750024

[Mathadon]

I think the following model produces the bug I was hunting:

model test
  extends ControlledFlowMachineDynamic(fan3(redeclare 
        Buildings.Fluid.Movers.Data.Pumps.Wilo.Stratos32slash1to12 per));
end test;

for instance:

ERROR: Failed to solve non-linear system using Newton solver.
To get more information: Turn on Simulation/Setup/Debug/Nonlinear solver diagnostics/Details
Solution to systems of equations not found at time = 2E-06
   Nonlinear system of equations number = 9
   Infinity-norm of residue = 0.676148
   Iteration is not making good progress.
   Accumulated number of residue       calc.: 404
   Accumulated number of symbolic Jacobian calc.: 11
   Last values of solution vector:
fan3.y_actual = 0
   Last values of residual vector:
{ -0.676148 }

[Mathadon]

The algebraic loop that is failing is:

// Nonlinear system of equations number = 9.
// It depends on the following parameters: 
//   fan3.eff.VDelta_flow
//   fan3.eff.cBar[1]
//   fan3.eff.cBar[2]
//   fan3.eff.dpDelta
//   fan3.eff.preDer3[10]
//   fan3.eff.preDer3[11]
//   fan3.eff.preDer3[1]
//   fan3.eff.preDer3[2]
//   fan3.eff.preDer3[3]
//   fan3.eff.preDer3[4]
//   fan3.eff.preDer3[5]
//   fan3.eff.preDer3[6]
//   fan3.eff.preDer3[7]
//   fan3.eff.preDer3[8]
//   fan3.eff.preDer3[9]
// It depends on the following timevarying variables: 
//   fan3.VMachine_flow
//   fan3.senRelPre.p_rel
// Start values for iteration variables of non-linear system of 1 equations: 
//   fan3.y_actual

equation // Residual equations
0 = fan3.senRelPre.p_rel-Buildings.Fluid.Movers.BaseClasses.Characteristics.pressure
(
Buildings.Fluid.Movers.BaseClasses.Characteristics.flowParametersInternal(
n = 11, 
V_flow = {0.0, 2.11830535572E-05, 0.000167865707434, 0.000700939248601, 
  0.0012450039968, 0.00177258193445, 0.00227268185452, 0.00272332134293, 
  0.00312450039968, 0.00345423661071, 0.005142914914853692}, 
dp = {59303.20534016437, 59279.55363280863, 59115.77671583633, 59002.16827665496,
   57354.827841611084, 54451.37923921096, 50291.289062817115, 44873.49049855453,
   38337.46221942317, 32076.924141548094, 0.0}
), 
fan3.VMachine_flow, 
fan3.y_actual, 
fan3.eff.VDelta_flow, 
fan3.eff.dpDelta, 
0.005142914914853692, 
59303.20534016437, 
fan3.eff.preDer3, 
0.05, 
fan3.eff.cBar, 
2882.762320687327);

so it seems that Buildings.Fluid.Movers.BaseClasses.Characteristics.pressure needs to be reformulated such that it can be solved towards y_actual = 0 for dp=0.

[Mathadon]

One problem seems to be that Buildings.Fluid.Movers.BaseClasses.Characteristics.pressure has local minima, probably caused by the function structure:

if r_N >= delta then
     dp := performanceCurve(V_flow=V_flow, r_N=r_N, d=d,
                            per=per, dimD=dimD);
  elseif r_N <= delta/2 then
    dp := flowApproximationAtOrigin(r_N=r_N, V_flow=V_flow,
                                    VDelta_flow=  VDelta_flow, dpDelta=dpDelta,
                                    delta=delta, cBar=cBar);
  else
    dp := Modelica.Fluid.Utilities.regStep(x=r_N-0.75*delta,
                                           y1=performanceCurve(V_flow=V_flow, r_N=r_N, d=d,
                                                               per=per, dimD=dimD),
                                           y2=flowApproximationAtOrigin(r_N=r_N, V_flow=V_flow,
                                                   VDelta_flow=VDelta_flow, dpDelta=dpDelta,
                                                   delta=delta, cBar=cBar),
                                           x_small=delta/4);
  end if;

Following graph shows the residual error as a function of r_N when getting stuck with the following error message:

ERROR: Failed to solve non-linear system using Newton solver.
To get more information: Turn on Simulation/Setup/Debug/Nonlinear solver diagnostics/Details
Solution to systems of equations not found at time = 0.116
   Nonlinear system of equations number = 9
   Infinity-norm of residue = 3.0777
   Iteration is not making good progress.
   Accumulated number of residue       calc.: 3512
   Accumulated number of symbolic Jacobian calc.: 303
   Last values of solution vector:
fan3.y_actual = 0.0360794
   Last values of residual vector:
{ -3.0777 }

The real solution at this time step is r_N = +-0.029 and the solution at the previous time step r_N=0.0066.

[Mathadon]

Some more analysis:

Below is a schematic of one of the parallel branches of ControlledFlowMachineDynamic.

The problem seems to be the position of the pump mass flow rate measurement, combined with the compressibility of the fluid. Consider for instance the case where dp_in increases from 0 Pa to 1 Pa. At that point in time vol.p -bou1.p = 1 and vol.p-bou.p=0. Therefore the mass flow rate through res1 will be >0 and the mass flow rate through res will be equal to zero. So the pump model 'sees' a mass flow rate of zero for dp_in=1 because of the delay introduced by the mixing volume. This is a valid operating point of the pump, but is probably more difficult to solve numerically. Furthermore I can imagine that when using DASSL, a situation can occur when m_flow<0 when dp_in>0 due to iterations on the state value of vol.p, which is probably more problematic.

Switching the position of the mixing volume and the mass flow rate measurement seems to fix the problem. edit: there are still problems with this fix

[Mathadon]

In https://github.com/iea-annex60/modelica-annex60/tree/issue458_singularityInMovers I changed Annex60.Fluid.Movers.Validation.ControlledFlowMachine to demonstrate a bug.

When running simulateModel("Annex60.Fluid.Movers.Validation.ControlledFlowMachine", stopTime=600, method="dassl", resultFile="ControlledFlowMachine"); I get:

The odd result here is that for a declining r_N, dp rises again when entering the regularisation region around 0.05 < r_N<0.025 and even fluctuates in the region for r_N<0.025 when using flowApproximationAtOrigin(.).

[Mathadon]

I think I implemented a fix in d388262 . This removes the existing regularisation from the paper but seems to work well in the unit tests. In the unit tests the bump in Annex60.Fluid.Movers.Examples.MoverContinuous is also removed. (see below)

A numerical Jacobian is however introduced by the use of InverseXRegularized. @mwetter do you know how to get rid of this?

If we continue with this we should

• [x] remove all cBar etc computations
• [x] try to inline the code such that the c-code contains no function calls

[rubenbaetens]

I adopted these changes locally and it solved the issues i had in bigger models. I do notice the simulation is a bit slower though, but that's preferred in contrast to a crashing model.

[Mathadon]

@rubenbaetens the numeric jacobian is fixed in 5960479 and 7683e02 This should fix the computing time.

[Mathadon]

I'm still experiencing problems with the current problem formulation. @rubenbaetens you too?

edit: this was for a case where the pump was off, but where mass was being forced through it by a mass flow source

[rubenbaetens]

I don't experience any problem anymore (but didn't implement the solution for the Numerical Jacobian)

[rubenbaetens]

Could we move forward on this issue ? Since i adopted this implementation, i no longer experienced an error on this topic. Hence i would liek to introduce it to IDEAS as soon as possible to ease workflow.

[mwetter]

@mathadon: In preparation of Friday's call, I'd like to look at the proposed code and run a diff against the master to see the changes. Which branch is current, issue458_fix or issue458_singularityInMovers?

[Mathadon]

The code probably needs some work before we converge, but the latest version is on issue458_fix.

[mwetter]

@Mathadon Can you please have a look at https://github.com/iea-annex60/modelica-annex60/compare/issue458_fix#diff-49e9a61488342b0aaf0ad2168a815951R60 before the call. I am not comfortable with the current formulation because I don't understand how it behaves for very small r_N, but I provide an alternate idea.

[mwetter]

@Mathadon why does inverseXRegularized has the section

   annotation (smoothOrder=2,
+  derivative(order=2,

smoothOrder=2 is correct, but derivative(order=2, ... is wrong as der_inverseXRegularized is the first derivative, and hence it should be order=1, or simply be omitted.

• [x] Please clarify and/or fix if you agree.

@Mathadon, @rubenbaetens I refactored the movers, not as in the file I sent to Filip because this caused unstable behavior in Examples.PumpsSeries. The pressures of the two pumps oscillated with an amplitude of around 1E20 Pascal but opposite sign, so the model run, the equations were satisfied, but it created garbage results. Good thing we have so many test. The reason was that in the limit, the regularization was correct, but this is an asymptotic behavior only, and the model was not yet operating in this asymptotic region, hence the V_flow/max(small, r_N) term dominated the r_N^2 term. Therefore, I reformulated Characteristics/pressure.mo differently and the tests are now fine.

I still want to test it on Buildings before we merge it to the master. It would be good if you can test your problematic case in IDEAS with this version (commit 220d4c520057c93c4acdce9d25e9a91f2f7d6ad9). All regression tests in Annex60 have been updated, and the results correspond to Dymola 2017.

[Mathadon]

I set order=1 and checked the solution.

I'll try to run my model again with the new solution this afternoon.

[Mathadon]

Trouble in paradise :(

One of my non-linear algebraic loops converges up to

Scaled residual 6.549607005268222E-13
Old Residual 1.22845989090722E-13
Predicted relative decrease 0
Actual relative decrease -4.331559502875851
basement.tabsPumps.pmpNrd.y_actual = 4.87526E-08

It then starts a line search to find a sign change in the residual, but this is not found:

and dymola throws an error. The point near the origin of this plot is { 0, 1.22845989090722e-13, -1.22845989090722e-13}.

So even though the residual is very small the Newton solver does not consider the solution to be found, which may or may not be an issue of Dymola. Dymola may reason that if no sign change in the residual can be detected near the point of convergence, the solution is not a true solution? So what may help is to get a sign change in the residual function for r_N<0. This can be obtained by replacing

dp:=r_N^2*Annex60.Utilities.Math.Functions.cubicHermiteLinearExtrapolation(
              x=rat,
              x1=per.V_flow[i],
              x2=per.V_flow[i + 1],
              y1=per.dp[i],
              y2=per.dp[i + 1],
              y1d=d[i],
              y2d=d[i+1]);

by

dp:=sign(r_N)*r_N^2*Annex60.Utilities.Math.Functions.cubicHermiteLinearExtrapolation(
              x=rat,
              x1=per.V_flow[i],
              x2=per.V_flow[i + 1],
              y1=per.dp[i],
              y2=per.dp[i + 1],
              y1d=d[i],
              y2d=d[i+1]);

this way dp becomes negative for negative r_N, which should introduce a sign change in typical residual functions.

I also used this factor in my previous implementation, which worked fine. I tested this change and indeed this bug disappeared. However, I am now experiencing issues where a negative y_actual is considered to be a valid solution for VMachine_flow>0 and dp_in>0..

Any ideas?

[mwetter]

@Mathadon I would also have tried he same approach as you did with dp:=sign(r_N)*r_N^2*Annex60....

We should keep your new implementation as this will move the fan curves into the III. quadrant for r_N < 0. Can you check if the cubicHermiteLinearExtrapolation is very steep in the area where you have divergence. We could consider something like

if r_N >= 0 then
dp:=sign(r_N)*r_N^2*Annex60.Utilities.Math.Functions.cubicHermiteLinearExtrapolation(
              x=rat,
              x1=per.V_flow[i],
              x2=per.V_flow[i + 1],
              y1=per.dp[i],
              y2=per.dp[i + 1],
              y1d=d[i],
              y2d=d[i+1]);
else
  dp:= -r_N^2* ( dpMax - dpMax/V_flow_max * V_flow);
end if;

which would ensure that in the III. quadrant, the fan curve is not very steep but rather has the average slope of the fan curve.