[FEATURE] Add homotopy where relevant

[pierre-eliep-met]

Is your feature request related to a problem? Please describe.
For model convergence, we know that adding homotopy is a good thing, but where is it best to add it without having to go through everything ?

I am thinking of the heat exchange equations, maybe th W variable in flow model because it adds non linearity absolutely everywhere, maybe in T = state.T

what else ?

[AndreaBartolini]

homotopy is not easy to be implemented and requires an accurate preliminary design... we can organize a call to discuss about it.

[casella]

I'd say it is very easy to implement homotopies, less easy to implement ones that actually work :) I'm also in favour of a live discussion. Some possibly useful strategies are listed in this paper.

[pierre-eliep-met]

From Francesco's last email added here for traceability :

• replacing quadratic and fluid-property dependent pressure loss correlations with pressure losses depending linearly on the flow rate, with coefficients based on nominal/design conditions
• replacing heat transfer coefficients obtained from correlations with constant nominal/design values
• replacing flow rates with nominal values in energy balance equations; however, this can lead to singular problems in some cases, i.e., when the flow rate must be modulated to obtain certain outlet temperatures)
• for very complicated plants, it can be useful to add decouplers that fix the outlet enthalpy at lambda = 0 to some nominal/design value. This avoids the first Newton iteration to get too far from the solution, to points where the solution cannot be recovered

[pierre-eliep-met]

@casella, @AndreaBartolini, just to properly understand homotopy, if I replace :

T_sensor.T = T; // T is an output

By

T_sensor.T = homotopy(T, T_0); // T is an output, T_0 is a parameter

Does it add an equation to the system ? For the first homotopy iteration, is it equivalent to solving the same system but with T_sensor.T = T_0 equation instead of T_sensor.T = T, which means 1 additional equation to the system ?

[AndreaBartolini]

As far as I know the homotopy oprator doesn't "add" equations but "substitute" the "actual" model with a linear combination of "simplified" and "actual" models during the initialization phase only.

From the paper cited above:
homotopy(actual=a*b, simplified=c/d)

is expanded to
λ (a · b) + (1 − λ ) (c/d)

During the initialization phase the λ coefficient goes from 0 to 1 in a number of steps that can be set by the following simulation flag in OpenModelica

-ils=value or -ils value
Value specifies the number of steps for homotopy method (required: -iim=symbolic). The value is an Integer with default value 1.

[pierre-eliep-met]

But then, in your example, if a and b are outputs of the model, and c and d are inputs, then at λ=0, it changes the number of equations, no ?
Isn't it the problem stated by @casella during our meeting, when we fix inlet temperature of heat exchangers in series ?

[pierre-eliep-met]

I ran a first test, with homotopy as done in branch test-homotopy, and got the following results, which are a bit frustrating : homotopy is implemented for h_in and DP friction equations, and it did not help at all, even for the metroscopia NPP models, while I am sure of the values of the init parameters (finishing with _0) used in homotopy equations

Model	Without homotopy (main branch)	h_in homotopy (dymola)	DP homotopy (dymola)	Both homotopy (dymola)	Both homotopy (OpenModelica)
MetroscopiaNPPP_reverse	0 errors, CV ✅	0 errors, CV ✅	0 errors, CV ✅	0 errors, CV ✅	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.
MetroscopiaNPP_direct	5 errors, then CV ✅	1 error, then CV ✅	5 errors then CV ✅	2 errors, then CV ✅	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.
MetroscopiaNPP_faulty	0 errors, CV ✅	0 errors, CV ✅	0 errors, CV ✅	0 errors, CV ✅	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.
MetroscopiaNPP_internal_ramp_direct	5 errors, then CV ✅	1 error, then CV ✅	5 errors, then CV ✅	2 errors, then CV ✅	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.
ParallelTurboFWP_direct	No CV ❌	No CV ❌	No CV ❌	No CV ❌	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.
MetroscopiaCCGT_causality_reverse	4 errors then CV ✅	8 errors then CV ✅	4 errors then CV ✅	8 errors then CV ✅	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.
MetroscopiaCCGT_causality_direct	No CV ❌	No CV ❌	No CV ❌ could not solve simplified initialization for homotopy method	No CV ❌ could not solve simplified initialization for homotopy method	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.
MetroscopiaCCGT_faulty (no start values)	No CV ❌	No CV ❌	No CV ❌ could not solve simplified initialization for homotopy method	No CV ❌ could not solve simplified initialization for homotopy method	❌ Failed to solve the initialization problem with global homotopy with equidistant step size.

[AndreaBartolini]

But then, in your example, if a and b are outputs of the model, and c and d are inputs, then at λ=0, it changes the number of equations, no ? Isn't it the problem stated by @casella during our meeting, when we fix inlet temperature of heat exchangers in series ?

input and output are relative concepts in Modelica, they are used to tell to the compiler that some variables shall be supplied via bounding equations. The important thing is that the overall model is balanced after you have fixed the boundary conditions.
The number of equations doesn't change, due to both actual and simplified model are essentially equations, so they don't define new variables.

[AndreaBartolini]

I ran a first test, with homotopy as done in branch test-homotopy, and got the following results, which are a bit frustrating : homotopy is implemented for h_in and DP friction equations, and it did not help at all, even for the metroscopia NPP models, while I am sure of the values of the init parameters (finishing with _0) used in homotopy equations
...

To deal with homotopy is not easy, no surprise if it doesn't work at the first tentative... let us some days to investigate your examples...

[pierre-eliep-met]

To deal with homotopy is not easy, no surprise if it doesn't work at the first tentative... let us some days to investigate your examples...

No problem, tell me when you will have the time to look at it !

[casella]

@pierre-eliep-met we first analyzed MetroscopiaNPPP_reverse. The cause of failure in this case is articulated.

The lambda = 0 homotopy process starts and then fails with

Error in region computation of IF97 steam tables(p = 1.94e+06, h = 1564.22)

Unfortunately for some reason the function calling stack is not present in the error message. However, a quick search in the flat model shows that 1.94e6 is the value of HPT_P_out_sensor.P. The corresponding enthalpy HPT_P_out_sensor.h is computed by nonlinear system 2161. Turning on the doomsday device LOG_NLS_V, we got the log of the iterations of this nonlinear system:

############ Solve nonlinear system 2161 at time 0 ############
initial variable values:
x0 [4-dim]
nls status
variables
[ 1]                HP_extract.h_in  =       -47044.389   nom =          2500000   min =           -1e+10   max =            1e+10
[ 2]             HPT_P_out_sensor.h  =        1564.2201   nom =           500000   min =           -1e+10   max =            1e+10
[ 3]                P_cond_sensor.h  =        110604.54   nom =           500000   min =           -1e+10   max =            1e+10
[ 4]                turbines_eta_is  =       -235.11112   nom =              0.8   min =                0   max =                1
******************************************************
SYSTEM SOLVED
indRow: [4-dim]
indCol: [5-dim]
vector x (solution): [5-dim]
regular initial point!!!
******************************************************
NEWTON SOLVER STARTED! equation number:  2161
maximum number of function evaluation:  400
nls status
variables
[ 1]                HP_extract.h_in  =          2500000   nom =          2500000   min =           -1e+10   max =            1e+10
[ 2]             HPT_P_out_sensor.h  =           500000   nom =           500000   min =           -1e+10   max =            1e+10
[ 3]                P_cond_sensor.h  =           500000   nom =           500000   min =           -1e+10   max =            1e+10
[ 4]                turbines_eta_is  =              0.8   nom =              0.8   min =                0   max =                1
Iteration: 1
------------------------------------------------------
newton step
variables
[ 1]                HP_extract.h_in  =        183498.24   step =       -2316501.8   old =          2500000
[ 2]             HPT_P_out_sensor.h  =        36480.428   step =       -463519.57   old =           500000
[ 3]                P_cond_sensor.h  =         118490.4   step =        -381509.6   old =           500000
[ 4]                turbines_eta_is  =       -254.07235   step =       -254.87235   old =              0.8
Need to damp, grad_f =   -3.6864308190e+13
Need to damp, error_f =    4.2932684630e+06
Need to damp this!! lambda1 =    1.0000000000e+00
Need to damp, error_f1 =    3.3290545492e+05
Need to damp, forced error =    1.4745723276e+13
function values: [4-dim]
scaled function values: [4-dim]
error measurements:
delta_x        =   2.3930274669e+06
SOLVING NON-LINEAR SYSTEM USING MIXED SOLVER (Newton/Homotopy solver)
delta_x_scaled =   3.1859404838e+02
newtonXTol          =   1.0000000000e-12
error_f        =   3.3290545492e+05
error_f_scaled =   1.4074727195e-01
newtonFTol          =   1.0000000000e-12
Iteration: 2
newton step
variables
[ 1]                HP_extract.h_in  =        -48715.09   step =       -232213.33   old =        183498.24
[ 2]             HPT_P_out_sensor.h  =        1423.9318   step =       -35056.496   old =        36480.428
[ 3]                P_cond_sensor.h  =        110598.25   step =       -7892.1416   old =         118490.4
[ 4]                turbines_eta_is  =       -234.89776   step =        19.174586   old =       -254.07235
Need to damp, grad_f =   -2.2165208383e+11
Need to damp, error_f =    3.3290545492e+05
Need to damp this!! lambda1 =    1.0000000000e+00
EQUATION NUMBER: 2161
Need to damp, error_f1 =    2.3280773363e+03
Need to damp, forced error =    8.8660833531e+10
function values: [4-dim]
scaled function values: [4-dim]
error measurements:
delta_x        =   2.3497717542e+05
delta_x_scaled =   2.3968520439e+01
newtonXTol          =   1.0000000000e-12
error_f        =   2.3280773363e+03
error_f_scaled =   1.0286670902e-03
TIME:   0.0000000000e+00
newtonFTol          =   1.0000000000e-12
Iteration: 3
newton step
variables
[ 1]                HP_extract.h_in  =       -47044.522   step =         1670.568   old =        -48715.09
[ 2]             HPT_P_out_sensor.h  =        1564.2183   step =        140.28657   old =        1423.9318
[ 3]                P_cond_sensor.h  =        110604.54   step =         6.284006   old =        110598.25
[ 4]                turbines_eta_is  =        -235.1111   step =      -0.21333293   old =       -234.89776
Need to damp, grad_f =   -1.0839888168e+07
Need to damp, error_f =    2.3280773363e+03
Need to damp this!! lambda1 =    1.0000000000e+00
Need to damp, error_f1 =    1.7858313325e-01
Need to damp, forced error =    4.3359552671e+06
function values: [4-dim]
scaled function values: [4-dim]
number of function calls (so far!):  0
error measurements:
delta_x        =   1.6764597914e+03
delta_x_scaled =   2.6666714724e-01
newtonXTol          =   1.0000000000e-12
error_f        =   1.7858313325e-01
error_f_scaled =   7.2644846311e-08
newtonFTol          =   1.0000000000e-12
Iteration: 4
newton step
variables
[ 1]                HP_extract.h_in  =       -47044.389   step =        0.1331652   old =       -47044.522
[ 2]             HPT_P_out_sensor.h  =        1564.2201   step =     0.0017559175   old =        1564.2183
[ 3]                P_cond_sensor.h  =        110604.54   step =   -0.00046142134   old =        110604.54
[ 4]                turbines_eta_is  =       -235.11112   step =   -1.7816597e-05   old =        -235.1111
Need to damp, grad_f =   -6.3783870960e-02
System values [4-dim]
Need to damp, error_f =    1.7858313325e-01
Need to damp this!! lambda1 =    1.0000000000e+00
Need to damp, error_f1 =    1.1591491242e-09
Need to damp, forced error =    2.5513548384e-02
function values: [4-dim]
scaled function values: [4-dim]
error measurements:
delta_x        =   1.3317758140e-01
delta_x_scaled =   2.2270810824e-05
newtonXTol          =   1.0000000000e-12
Nominal values [4-dim]
error_f        =   1.1591491242e-09
error_f_scaled =   4.5587456277e-16
newtonFTol          =   1.0000000000e-12
Iteration: 5
newton step
Need to damp, grad_f =   -2.6872533842e-18
Need to damp, error_f =    1.1591491242e-09
Need to damp this!! lambda1 =    1.0000000000e+00
Need to damp, error_f1 =    4.6053133342e-11
Need to damp, forced error =    1.0749013537e-18
Scaling values [5-dim]
function values: [4-dim]
scaled function values: [4-dim]
error measurements:
delta_x        =   8.9294872765e-10
delta_x_scaled =   1.5945419522e-13
newtonXTol          =   1.0000000000e-12
error_f        =   4.6053133342e-11
error_f_scaled =   3.4606398682e-17
newtonFTol          =   1.0000000000e-12
NEWTON SOLVER DID CONVERGE TO A SOLUTION!!!

This is a very rare case in which a bad initial guess (HPT_P_out_sensor.h = 500000, P_cond_sensor.h = 500000) doesn't cause the Newton solver to fail, but rather to converge to a completely bogus solution:

[ 1]                HP_extract.h_in  =       -47044.389   nom =          2500000   min =           -1e+10   max =            1e+10
[ 2]             HPT_P_out_sensor.h  =        1564.2201   nom =           500000   min =           -1e+10   max =            1e+10
[ 3]                P_cond_sensor.h  =        110604.54   nom =           500000   min =           -1e+10   max =            1e+10
[ 4]                turbines_eta_is  =       -235.11112   nom =              0.8   min =                0   max =                1
error_f_scaled =   3.4606398682e-17

i.e., by computing completely unreasonable values for some of the unknowns (e.g. an isentropic efficiency of -235), the solver miraculously manages to get the residual to 2.3e-17. The problem is, the solver is completely agnostic regarding the physical meaning of these unknowns, so it accepts the solution without a hiccup. But then, the enthalpy value 1564 J/kg at 19.4 bar is too low (probably ice around 0.3 °C or something) so the solution process stops. In case this was accepted, the value HP_extract.h_in = -47044.389 would have broken the calculation later on anyway, not to mention the negative eta_is.

There are two fundamental issues here. The first, as specified in the old homotopy paper is that the specific enthalpy vs. temperature relationship is not too nonlinear, provided the initial guess is in the right phase region, otherwise all hell breaks loose with the Newton solver. The problem here is that 500,000 J/kg is likely too low a value for the slightly humid vapour at the outlet of the HPT, since water at such enthalpy and 19.4 bars is heavily subcooled liquid at 120 °C. Ditto for p_cond_sensor.h.

So, it is not important to have precise start values for those enthalpies, but you have to be careful that they fall in the right region to avoid these problems.

The second issue is that the specific enthalpy type that you are using for the water/steam medium is SI.specificEnthalpy, where you changed the start and nominal attribute. Hence, it does not have any min/max attribute, as the medium-specific WaterIF97.SpecificEnthalpy type instead does. If a specific enthalpy type with min = 40000 (corresponding to liquid water at about 10 °C) was used, then this bogus solution would have generated an error or at least a warning with a proper message (too low specific enthalpy) instead of continuing until a function call to specificEnthalpy_ph broke in a later computation. So, either you use the WaterIF97-specific types, or you should add appropriate min/max attributes to your custom-defined types.

This also applies to the eta_is variable, that should probably have something like min = 0.5, max = 1.0 attributes, so that a solution with eta = -235 is flagged as inappropriate.

The general idea is always the same: the solver is agnostic, the modeller can make it less so, by adding extra information that the modeller is aware of, e.g. in terms of min/max bounds, order of magnitude (nominal), etc.

Last, but not least, the fact that eta_is can become an iteration variable gives room to extremely nonlinear and wild behaviour if Newton does not behave. One suggestion could be to use homotopy for turbine efficiencies as well, using the nominal value for the simplified model. This will make the initial system a lot less nonlinear, and a lot less prone on converging on totally bogus solutions.

Please try to implement these improvements and then re-run this test case, we'll see if there are other problems that were masked by these ones, or if you get a solution straight out of the box.

[casella]

We also checked MetroscopiaNPP_direct. It fails during initialization at lambda = 0 when trying to solve system 1917, which only has one iteration variable, namely steam_dryer.P. The provided start value is 1.94e6, and if you turn on LOG_NLS_V, you can see that after one iteration this variable changes wildly, actually wants to become negative. We reduced the start value to 1.64e6 and now after one iteration it became 7.2e6.

This means that steam_dryer.P has a very strongly nonlinear effect on the equations, otherwise the behaviour of the first Newton iteration would not change so much for relatively small changes of the initial guess.

If you look into the nonlinear system of equations, the simplified energy balance with Q_0 instead of Q shows up:

(torn) steam_dryer.steam_phase.W := steam_dryer.steam_phase.Q_0 * (steam_dryer.steam_phase.h_out - steam_dryer.steam_phase.h_in)

This homotopy can be a bit dangerous in some cases, so I would suggest to try avoid using it and see what happens with this test case.

[casella]

More analysis on MetroscopiaNPP_direct. This is the failing nonlinear system:

The iteration variable steam_dryer.P. Based on that, the saturated liquid and vapour enthalpies are computed, then the specific enthalpy of 99% dry saturated steam leaving, then the transferred power steam_phase.W, via an energy balance on the steam side, which has the flow rate fixed to steam_dryer.steam_phase.Q_0.

The iteration variable steam_dryer.P needs then to be adjusted by Newton iterations to match the energy balance on the liquid side. Since the two flow rates are fixed, the only way to do that is to adjust the saturation properties, but I guess this is a bit too overconstrained, and results in a heavily nonlinear system of equations. I'm not exactly sure how this steam dryer operates in real life, but I guess there could be some control loop acting on the liquid water flow - if you fix the flow rate to Q_0, this loop is opened and the simplified model cannot work properly.

I would then suggest to switch off using Q_0 in the energy balance for this component. More in general, I would switch that homotopy off by default, because it can be dangerous in several non-trivial contexts. It may be useful in some cases, but you should turn it on carefully.

[pierre-eliep-met]

The problem here is that 500,000 J/kg is likely too low a value for the slightly humid vapour at the outlet of the HPT, since water at such enthalpy and 19.4 bars is heavily subcooled liquid at 120 °C. Ditto for p_cond_sensor.h.

In the turbine, all enthalpies have an appropriate start value. 5e5 J/kg is the default value in the enthalpy unit, but it is redeclared in the StodolaTurbine's definition :

  extends MetroscopeModelingLibrary.Partial.BaseClasses.FlowModel(
    Q_0=1500,
    P_in_0=60e5,
    P_out_0=55e5,
    h_in_0 = 2.7e6,
    h_out_0= 2.6e6,
    ...
  )
  Utilities.Units.SpecificEnthalpy h_real(start=h_out_0); // Enthalpy after real decompression
  Utilities.Units.SpecificEnthalpy h_is(start=h_out_0/0.8); // Enthalpy after isentropic decompression

For the eta_is it has unit Yield (which is the same as Fraction), who's bounded between 0 and 1...

So I don't really understand this start value or bound issue, since we already implemented that

[AndreaBartolini]

The problem is not in the turbine but in the pressure sensor downstream the turbine (HPT_P_out_sensor).
Probably it is necesary to update the start values also in that model (or most probaly into its instantiation).

[casella]

@pierre-eliep-met, @valentind-met we analyzed this issue more in depth.

One of the problems here is probably in the handling of start attributes for alias variables. When generating code from a Modelica modular model, it is often the case that many variables end up in the flat model that are equal to each other (up to a minus sign), e.g. because of connector variables.

Modelica tools, including OMC, use alias elimination algorithms to remove the trivial equations and replace all instances of variables from the alias set with just one representative variable. In most cases, one of these variables has a carefully selected start attribute, which is meant to be used in order to achieve convergence during initialization, while other ones (e.g. connector variables, or sensor variables in this case) have generic default start attributes. The question then is, how to retain the carefully selected value among the set of alias variables.

This is a well-known problem, that has plagued the community of Modelicans doing thermo-fluid system models since the early 2000's. The problem was first addressed in 2012 by adding section 8.6.2 to the language specification, which proposed a criterion for priority ranking for start values in alias sets. That criterion was revised in October 2021 (see modelica/ModelicaSpecification#3008) and is now included in the 3.6 version of the language specification, which is in the final approval and voting phase.

This criterion is implemented in Dymola, where apparently it does its job quite well. Unfortunatetly, it is still not implemented in OpenModelica, see OpenModelica/OpenModelica#5813. We will make some further enquiries, and if it turns out that start attribute priority is actually the problem preventing the successful simulation of non-trivial Metroscopia model in OMC, then it will become essential to provide to the implementation of the latest rules there. At that point, you may discuss some arrangement with the OSMC to get that done, e.g. a DFD contract. I guess the effort should be limited, around 1 person-month.

[pierre-eliep-met]

Thanks for the details. Yes I had the same idea of the problem. To try to tackle that we have defined all the start values parameters in all components of the metroscopia models, so all variables should now have an appropriate start value, and there should now be appropriate values for all equations of simplified models through homotopy.

I'll test that asap and keep you updated, but I must say I don't really have the time at the moment, I hope I will next week or the one after.

[pierre-eliep-met]

Actually I took the time to rerun simulations on metroscopia NPP direct and reverse which now have start values parameters everywhere, and it works better in dymola (less errors for direct model) but nothing has changed in OMEdit...

my changes are on branch test-homotopy if you can have a look @casella

[pierre-eliep-met]

And even the MetroscopiaNPP_direct_with_StartValues does not simulate, while all variable start values are good, and all init parameters used in homotopy simplified models are close to solution...

[casella]

In fact, there is another issue with the MetroscopiaNPP_direct model. If you compile it with a recent nightly build (as we do, so we get the latest improvements of the software), we get the following notification:

The linear system: 
1 : HPT_P_out_sensor.Q + HP_extract_P_sensor.Q - HPT_P_in_sensor.Q = 0.0
2 : superheater_bleed_P_sensor.Q + HPT_P_in_sensor.Q - P_steam_sensor.Q = 0.0
3 : superheater_bleed_P_sensor.Q - 1e-05 = -1.001 + superheater_drains_pipe.Q - -1.001
4 : HP_heater.Q_hot_in + (-superheater_drains_pipe.Q) - HP_extract_P_sensor.Q = 0.0
[
  1.0 , 0.0 , 0.0 , -1.0 ;
  0.0 , 0.0 , 1.0 , 1.0 ;
  0.0 , -1.0 , 1.0 , 0.0 ;
  -1.0 , -1.0 , 0.0 , 0.0
]
  *
[
  HP_extract_P_sensor.Q ;
  superheater_drains_pipe.Q ;
  superheater_bleed_P_sensor.Q ;
  HPT_P_in_sensor.Q
]
  =
[
  -HPT_P_out_sensor.Q ;
  P_steam_sensor.Q ;
  1e-05 ;
  -HP_heater.Q_hot_in
]
 might be structurally or numerically singular for variable HPT_P_in_sensor.Q since U(4,4) = 0.0. It might be hard to solve. Compilation continues anyway.

Indeed, the LHS of eq. 3 is the sum of the LHS of equations 1., 2., and 4., so the system is indeed structurally singular. Equations 1., 2., and 4. clearly stem from mass balance equations at connection points, as it is also reported by the declarative debugger, e.g.:

Unfortunately, we couldn't really figure out where eq. 3. comes from, because the debugger doesn't show any significant symbolic transformation chain, nor it points to any line in the original Modelica code. We need to investigate this issue further; Dymola does not report any structural singularity, so I guess something goes wrong in OMC in this respect.

There are also two further issues with the system structure, which may be related to the previous point. The backend further reports:

[12] 13:17:13 Translation Warning
Assuming fixed start value for the following 7 variables:
         loopBreaker.loop_flow_error:VARIABLE(flow=false fixed = true ) MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse, MetroscopeModelingLibrary.WaterSteam.Pipes.LoopBreaker type: Real
         HP_reheater_drains_control_valve.Cv:VARIABLE(start = 10000.0 unit = "m4/(s.N5)" fixed = true )  "Cv"MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse, MetroscopeModelingLibrary.WaterSteam.Pipes.ControlValve type: Real
         HP_heater.HX_subcooling.epsilon:VARIABLE(min = 0.0 max = 1.0 start = HP_heater.HX_subcooling.epsilon_0 unit = "1" fixed = true nominal = 0.5 ) MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse, MetroscopeModelingLibrary.WaterSteam.HeatExchangers.Reheater, MetroscopeModelingLibrary.Power.HeatExchange.NTUHeatExchange type: Real
         steam_dryer_liq_out_pipe.Qv:VARIABLE(min = 0.0 start = 0.1 unit = "m3/s" fixed = true nominal = 0.1 )  "Mean volumetric flow rate"MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse, MetroscopeModelingLibrary.WaterSteam.Pipes.PressureCut type: Real
         LP_reheater_drains_control_valve.Cv:VARIABLE(start = 10000.0 unit = "m4/(s.N5)" fixed = true )  "Cv"MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse, MetroscopeModelingLibrary.WaterSteam.Pipes.ControlValve type: Real
         cold_source.Qv_out:VARIABLE(max = 0.0 start = -1.0 unit = "m3/s" fixed = true nominal = -0.1 ) MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse, MetroscopeModelingLibrary.WaterSteam.BoundaryConditions.Source type: Real
         superheater.HX_condensing.epsilon:VARIABLE(min = 0.0 max = 1.0 start = superheater.HX_condensing.epsilon_0 unit = "1" fixed = true nominal = 0.5 ) MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse, MetroscopeModelingLibrary.WaterSteam.HeatExchangers.Superheater, MetroscopeModelingLibrary.Power.HeatExchange.NTUHeatExchange type: Real

This is quite strange: adding fixed = true attributes is something that is done if the initialization problem is underdetermined, but as the Metroscopia models are purely algebraic and without any fixed = false parameters and no initial equation, there is no reason why this should happen.

Even more strange is what is reported afterwards:

[13] 13:17:14 Translation Warning
The initial conditions are over specified. The following 7 initial equations are redundant, so they are removed from the initialization_lambda0 system:
         1e-05 * CW_P_in_sensor.P = CW_P_in
         HP_heater_T_drains_sensor.T = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(HP_reheater_drains_control_valve.P_in, HP_heater_T_drains_sensor.h, 0, 0).T
         1e-05 * Q_feedwater_sensor.P = HP_heater_P_out
         1e-05 * deaerator_inlet_pipe.P_in = LP_heater_P_out
         superheater_T_out_sensor.T = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(LPT1.P_in, superheater_T_out_sensor.h, 0, 0).T
         superheater_drains_P_sensor.h = Modelica.Media.Water.IF97_Utilities.BaseIF97.Regions.boilingcurve_p(MetroscopeModelingLibrary.Examples.Nuclear.MetroscopiaNPP.MetroscopiaNPP_reverse.superheater.WaterSteamMedium.setSat_p(superheater.hot_side_deheating.P_in).psat).h
         steam_dryer.steam_phase.h_out = 0.99 * steam_dryer.h_vap_sat + 0.01000000000000001 * steam_dryer.h_liq_sat.

after adding 7 extra initial equations, the backend removes 7 other equations, that are declared to be "redundant" but really don't seem redundant at all to me.

Small wonder the model cannot be initialized even with all start values in place 😅.

Can we share this model with @kabdelhak, which is currently responsible for the OMC backend development, so he can have a look and figure out what is going wrong?

Thanks!

[pierre-eliep-met]

No problem to share that with them, but they also have to signe the license contract, as well as you had done before (I'm going to send it to you by email), sorry for the inconvenience, I'm looking forward to full open source publication...

[casella]

@pierre-eliep-met, @AndreaBartolini, I think I have a clue on what the root cause of this problem is, see OpenModelica/OpenModelica#10497. I hope we can solve the issue in OMC with the simple MWEs reported there, without the need of actually sharing the Metroscopia models.

[pierre-eliep-met]

@casella and @AndreaBartolini I guess for now there has not been any change on the issue linked above ? Just to know if we can investigate some new things or still need to wait.
For your information, I also tested the same homotopy methods with a complete NPP model and it did not improve convergence, though it was a rapid test, I didn't spend a lot of time to investigate

[AndreaBartolini]

The Osmc are still investigating the problem. It is better to wait that the OpenModelica/OpenModelica#10497 issue will be solved.