Linearization model

[MohammadHadiAlizadeh]

Hello
When I use the Linearize command for a hybrid system( contains if- else statement) I get the following error:

>> a=omc.linearize()
Index in position 2 is invalid. Array indices must be positive integers or logical values.

Error in OMMatlab/getLinearMatrixValues (line 843)
                    tmpMatrix(r,c)=val;

Error in OMMatlab/getLinearMatrix (line 818)
            tmpMatrix_A=getLinearMatrixValues(obj,matrix_A);

Error in OMMatlab/linearize (line 777)
                    result=getLinearMatrix(obj);

[MohammadHadiAlizadeh]

@arun3688 can you please shed some light or do possible fixes? Thanks in advance

[arun3688]

@MohammadHadiAlizadeh could you make a minimal example so that i can reproduce the error

[MohammadHadiAlizadeh]

@arun3688 I think the main problem is because of the model complexity. Unfortunately, I'm not allowed to share the model here. But I will try to generate the error in a test model, But I'm not sure whether it will cause the same error.

[MohammadHadiAlizadeh]

@arun3688 I have no idea why I'm getting such a tiny number and this error. Everything is OK in the minimal models, even hybrid ones! But when I removed the second mode of the problem to simplify it and did the linearization on the single mode problem, the following issue propagated.

[C:/Users/AT/AppData/Local/Temp/tp38c6880a_1e96_4a66_afbb_97651737c732/linearized_model.mo:179:727-179:749:writable] Error: 1.118765557412997e-312 cannot be represented by a double on this machine
[C:/Users/AT/AppData/Local/Temp/tp38c6880a_1e96_4a66_afbb_97651737c732/linearized_model.mo:179:1288-179:1310:writable] Error: 4.081919022638463e-315 cannot be represented by a double on this machine

Output argument "result" (and maybe others) not assigned during call to "OMMatlab/linearize".

Error in Linearize (line 5)
a=omc.linearize()

[MohammadHadiAlizadeh]

@arun3688 @casella This is the linearized model created during the omc.linearize() command; I'm attaching the file because it is too long to write inside the comment:
linearized_model.zip

As you can see, there are some numbers in the order of 1e-300 inside Matrix A that cause the error:

.
.
.

**Error: 1.118765e-312 cannot be represented by a double on this machine**
[C:/Users/AT/AppData/Local/Temp/tp38c6880a_1e96_4a66_afbb_97651737c732/linearized_model.mo:179:1288-179:1310:writable] Error: 4.081919022638463e-315 cannot be represented by a double on this machine

Output argument "result" (and maybe others) not assigned during call to "OMMatlab/linearize".

Error in Linearize (line 5)
a=omc.linearize()

One solution is to round the numbers like 1e-300 to 0 inside the linearized model during the matrix call, but I'm not sure how to do it with OMMatlab.

I'd be grateful if you could help to fix this issue.

[arun3688]

@MohammadHadiAlizadeh I will see whether i can reproduce the issue with the linearized model provided and try to fix the problem

[MohammadHadiAlizadeh]

Thanks @arun3688 .

When I try to load the linearized model, as the following:

import OMMatlab.*
 omc= OMMatlab()
 omc.ModelicaSystem("linearized_model.mo","linearized_model");

The same Error rises as when I tried to use omc.linearize() for the main model!!

.
.
.
 Error: 1.48219693752374e-323 cannot be represented by a double on this machine
[F:/Jul/linearized_model.mo:179:2422-179:2443:writable] Error: 1.48219693752374e-323 cannot be represented by a double on this machine

[casella]

This is really weird, I never saw such small numbers before. @arun3688 any idea why they show up?

[arun3688]

@casella I really don't know, the original model is not provided so it is hard to say, but i will investigate the issue

[MohammadHadiAlizadeh]

@arun3688 @casella I developed an example with very similar dynamics to the actual model. You can reproduce the error by using it. The files are attached:
simplified.zip

Also, this is the Matlab command I'm using for linearization

clear all
import OMMatlab.*
omc= OMMatlab()
omc.ModelicaSystem("Eshgh.mo","Eshgh",["mid.mo"]);

a=omc.linearize()

Thanks for your help

[casella]

@arun3688, as I understand it this linearized model is obtained by numerical differentiation. So, there are many elements that are supposed to be zero, but end up being computed as the difference of two very close numbers, hence resulting in numbers such as 2.159042138773611e-81 or 1.153519257740235e-315. We could try to avoid that in the first place by using structural information about the Jacobian, since the backend knows a-priori what is the sparsity pattern of the matrix, but I guess this would require some significant work, and maybe we don't bother.

Searching for the error message brought me here in the parser code:

      double d = strtod(chars,&endptr);
      if (!(*endptr == 0 && errno==0)) {
        c_add_source_message(NULL,2,ErrorType_syntax, ErrorLevel_error, "\%s cannot be represented by a double on this machine", (const char **)&chars, 1, $start->line, $start->charPosition+1, LT(1)->line, LT(1)->charPosition+1, ModelicaParser_readonly, ModelicaParser_filename_C_testsuiteFriendly);
        ModelicaParser_lexerError = ANTLR3_TRUE;
      }

As I understand, e.g. by looking here, the way underflow is handled by strtod() changed in C++11, i.e., when underflow happens, errno is not returned to be zero, but rather ERANGE.

I guess we need to update that part of the code to cope with that correctly.

Adding @sjoelund, @adrpo and @perost for comments.

[arun3688]

@MohammadHadiAlizadeh @casella , I can reproduce the issue with your model and i see that the issue is not related with OMMatlab , but with OpenModelica itself, It is not possible to build the generated linearized_model.mo in OpenModelica and this needs to be fixed in the run time in linearize.cpp

[MohammadHadiAlizadeh]

@arun3688 Thanks for investigating the Error.
@casella, How could the error be fixed? I need the linearized model as the prediction model for MPC. Can you please lend a helping hand?

[casella]

@MohammadHadiAlizadeh, as far as I understand, the problem is with the Modelica parser, which actually generates the error because of underflow not being handled correctly in C11. I would just fix that, e.g.

double d = strtod(chars,&endptr);
      if (!(*endptr == 0 && (errno==0 || (errno == ERANGE && d == 0)))) {
        c_add_source_message(NULL,2,ErrorType_syntax, ErrorLevel_error, "\%s cannot be represented by a double on this machine", (const char **)&chars, 1, $start->line, $start->charPosition+1, LT(1)->line, LT(1)->charPosition+1, ModelicaParser_readonly, ModelicaParser_filename_C_testsuiteFriendly);
        ModelicaParser_lexerError = ANTLR3_TRUE;
      }`

However, I'm not too much into changing the code myself, @arun3688 could you take care of that?

We could also improve linearize.cpp to avoid printing out those small numbers in the first place, but that would take more time, maybe we can do it later. At the end of the day it doesn't matter (the result is practically the same), we just need to avoid aborting the linear dump in case of underflow. @arun3688 do I miss something?

[casella]

@MohammadHadiAlizadeh in the meantime you can try to turn on the compiler flag --generateSymbolicJacobian. This doesn't always work, but if it does, it should generate proper 0s in all those places.

[MohammadHadiAlizadeh]

@casella , thanks for your continuous support. I used the flag but I'm still getting the same error. Can you please check out the attached code? I'm afraid I'm escaping something or entering the flag improperly.
simplified.zip

This is the Matlab code I'm using in:

clear all
import OMMatlab.*

omc=OMMatlab();

omc.ModelicaSystem("Eshgh.mo","Eshgh",["mid.mo"],"--generateSymbolicJacobian=true");

a=omc.linearize()

[arun3688]

@MohammadHadiAlizadeh It should be
omc.ModelicaSystem("Eshgh.mo","Eshgh",["mid.mo"], "--generateSymbolicJacobian");

[arun3688]

@casella I think we should directly fix in linearize.cpp, because that is function where the linearized model is generated and the error occurs when the generated linearized model is loaded and the parser fails, So the linearize should be fixed and the question is when to make the approximation eg. e-30

[MohammadHadiAlizadeh]

@MohammadHadiAlizadeh It should be omc.ModelicaSystem("Eshgh.mo","Eshgh",["mid.mo"], "--generateSymbolicJacobian");

@arun3688, thanks, but still getting the same error... :|

[arun3688]

@MohammadHadiAlizadeh that is true, it doesn't always work

[MohammadHadiAlizadeh]

@casella I think we should directly fix in linearize.cpp, because that is function where the linearized model is generated and the error occurs when the generated linearized model is loaded and the parser fails, So the linearize should be fixed and the question is when to make the approximation eg. e-30

@arun3688, Maybe the most standard effort is to make the threshold a user-specified parameter so it can be adjusted according to the specific application. It might need to develop a new flag.
However, 1e-30 is too low. 1e-10 is effectively zero too.

[casella]

@MohammadHadiAlizadeh introducing arbitrary "low" thresholds for numbers that should actually be zero is not a good idea. The elements on the diagonal of the A matrix are dimensionally the inverse of a time constant, so for applications with time constants of less than a year they should be larger than 3e-8, otherwise they can be approximated as zero. However, off-diagonal terms have dimensions that depend on the unit of the state variables involved; specifically, element A[i,j] has the the unit of x[i] divided by the unit of x[j] divided by "s". So it can be very large or very small, if you have pressures of high-pressure systems in Pa or powers of large power plants in W.

I wouldn't really do that, too dangerous.

@arun3688 there is no reason why the parser should fail when reading, e.g., 2.34e-325. The reason it currently fails is that at some point strtod changed the implementation and started to generate an error different from zero in case of underflow, when it actually returns zero, which is the expected result. So, we need to fix the parser first, so it doesnt fail in case of underflow, but just return a value of zero, which is exactly what we need 😅.

[sjoelund]

I'm not sure we should change the parser to silently convert small non-zero numbers to zero. Values that are non-zero are usually very different from values that are actually zero... Might be good to standardize the behaviour in the Modelica specification.

[casella]

@sjoelund, we need to change the parser so that it can convert to zeros numbers that cause an underflow of Real. The obvious reason is that the best approximation we have in IEEE 754 arithmetics of, say, 2e-350 is actually 0.0. Generating an error in this case is just stupid IMHO.

[casella]

BTW, this is what strtod() used to do until C11 😅

[adrpo]

We should not convert silently, we should do a warning and say we will translate it to zero.

There are however two cases here:

• too small // give a warning and transform to 0
• too large // give an error

How do we differentiate between the two?
Searching for e- in the number seems weird and people might actually not use the e notation.

[adrpo]

SO to the rescue :)
https://stackoverflow.com/a/30607006/1263889

[adrpo]

For this model:

model HelloWorld
  parameter Real x_too_small = 2e-350;
  parameter Real x_too_large = 2e350;
end HelloWorld;

Another tool gives:

Translation of HelloWorld:
The DAE has 0 scalar unknowns and 0 scalar equations.
The variable x_too_large=INF.0 is not in range [-2.22507e-308,1.79769e+308].

which is basically in line with what I want to do.

[adrpo]

With PR: OpenModelica/OpenModelica#9687 we get:

adrpo33@ida-0030 MINGW64 /c/home/adrpo33/dev/OMTesting/bugs/ommatlab_17
# omc HelloWorld.mo
Error processing file: HelloWorld.mo
Failed to parse file: HelloWorld.mo!
[C:/home/adrpo33/dev/OMTesting/bugs/ommatlab_17/HelloWorld.mo:2:32-2:38:writable] Warning: Underflow: 2e-350 cannot be represented by a double on this machine. It will be converted to 0.0.
[C:/home/adrpo33/dev/OMTesting/bugs/ommatlab_17/HelloWorld.mo:3:32-3:37:writable] Error: Overflow: 2e350 cannot be represented by a double on this machine

So we will convert underflow to zero and refuse to parse on overflow.

[adrpo]

@MohammadHadiAlizadeh: the parser is now changed to handle underflow as zero.
I see you are using Windows, I will start a windows nightly build so you can test it.

[MohammadHadiAlizadeh]

@adrpo , thanks a lot. I will be looking forward to hearing from you about the nightly build version.

[MohammadHadiAlizadeh]

@adrpo I installed the nightly build version and runned linearization, and I got some error like this:

[C:/Users/AT/AppData/Local/Temp/tpda5d3c2e_252b_4dd9_91cb_d61fd8e58ece/linearized_model.mo:89:1063-89:1085:writable] Warning: Underflow: 4.940656458412465e-324 cannot be represented by a double on this machine. It will be converted to 0.0.
[C:/Users/AT/AppData/Local/Temp/tpda5d3c2e_252b_4dd9_91cb_d61fd8e58ece/linearized_model.mo:121:704-121:726:writable] Error: Overflow: 2.121995644721841e-317 cannot be represented by a double on this machine
[C:/Users/AT/AppData/Local/Temp/tpda5d3c2e_252b_4dd9_91cb_d61fd8e58ece/linearized_model.mo:121:1337-121:1359:writable] Error: Overflow: 7.120236347223051e-310 cannot be represented by a double on this machine

Why does it treat some of the small numbers as an overflow?

[casella]

@adrpo some comments on my side.

I fully agree overflow should generate an error. No question about it.

Regarding underflow, I don't think we should generate a warning. Nobody in their right mind will ever write values like 3e-350 himself in a model (e.g. for a parameter value), there is no physical model for which such small value will be significantly different from zero for all practical purposes. In 99.9% of the cases, these underflows will be the result of some internal computations of OMC (in this case, a numerical differentiation) that ought to be zero in theory, but get some very small value because of rounding errors. The problem is, normal users have no idea about these internal computations, and most importantly they have no means to control it.

Basically, we'd be scaring users with warnings about a complete non-issue, that they have no means to control. Why should we do that? BTW, also the other tool doesn't do that, for a good reason 😅

So, I would remove the warning in case of underflow.

[casella]

Why does it treat some of the small numbers as an overflow?

I guess this is a bug in the fix 😄

[casella]

SO to the rescue :) https://stackoverflow.com/a/30607006/1263889

This is a slightly different issue, which I'd say is related to #8865. We have to fix that once an for all in 1.21.0.

[adrpo]

Hopefully I did the proper check now :)
OpenModelica/OpenModelica#9693

[adrpo]

Ok. I somehow still think we should keep the warning as the user might just have made a mistake and they should be informed so they can fix it.
Of course when this comes from linearization that should be fixed when we generate the MO file so we don't trigger underflow as that would confuse people.

[perost]

I think there's still something wrong with the check, because e.g. -4.940656458412465e-324 is actually a valid double and can be handled by strtod. But the parser gives a warning and converts it to 0 anyway.

[casella]

I think there's still something wrong with the check, because e.g. -4.940656458412465e-324 is actually a valid double and can be handled by strtod. But the parser gives a warning and converts it to 0 anyway.

Yes, I also noticed that. But in any case, I think that underflow of doubles is really not a big deal.

[perost]

So it seems numbers like -4.940656458412465e-324 are defined as subnormal, which is a special case for when the exponent is 0. strtod is allowed (but not required) to give an error in that case. So we could just ignore the error from strtod if the returned number is smaller than DBL_MIN but larger than 0.

[MohammadHadiAlizadeh]

@adrpo Thanks for the fix, now I can open and run the linearized_model.mo with OMs latest nightly build version. But still OMMatlab raises some error that I'm going to explain in the next comment.

[MohammadHadiAlizadeh]

@arun3688 , This is the error with OMMatlab, can you please have a look at that?
Actually, the linearize command generates the linearized_model.mo but it can not transfer the matrices to Matlab.

corresponding files: simplified.zip
Error:

Index in position 2 is invalid. Array indices must be positive integers or logical values.

Error in OMMatlab/getLinearMatrixValues (line 843)
                    tmpMatrix(r,c)=val;

Error in OMMatlab/getLinearMatrix (line 818)
            tmpMatrix_A=getLinearMatrixValues(obj,matrix_A);

Error in OMMatlab/linearize (line 777)
                    result=getLinearMatrix(obj);

Error in Main_Matlab_Linearizer (line 8)
a=omc.linearize()

[adrpo]

Unfortunately I have no Matlab to check this. @arun3688, please have a look.

[MohammadHadiAlizadeh]

I changed the OMMatlab code, and now it is working #18

[perost]

OpenModelica/OpenModelica#9696 fixes the parser so that we no longer give a warning about subnormal Real values. Mathematically speaking it probably doesn't matter at all, but we shouldn't give warnings about numbers that can naturally occur in the models.

[MohammadHadiAlizadeh]

I have another question :))
Based on OM linearize command documentation, "all variables that appear differentiated and that are selected as states at this time instant are treated as states x of the model." So, How can I define a variable as a state variable whose derivative is not provided in the model and expect linearize to treat it as a state variable? As you might know, some of the variables in most of the actual physical systems may not have the derivative but have to be treated as the state variables(e.g., Enthalpies).

Sth, like x2 in the following model, needs to be known and linearized as a state variable as well as x. I know x2 is a function of x1 and is accounted for der(x) calculation when the perturbation is applied, but is there any way to consider it as an independent state?

model test1

parameter Real a=2;

input Real u=0.1;
Real x;
output Real y1;
Real x2;

initial equation

x=1;

equation

der(x)= x^2+2*u+x2^3;

y1=2*x;
x2=5*x^2;

end test1;

@casella @adrpo @perost do you have any idea?
Thanks

[casella]

@MohammadHadiAlizadeh I'm sorry but I'm not sure I understand your question. In the example you provide, x2 is an algebraic variable in the DAE formulation. When the tool brings the DAEs into explicit ODE form, it takes them into account if they are needed to compute the state derivatives or the outputs: assuming the DAE is index 1, the algebraic equations are solved for the algebraic variables, assuming the states to be known. So, conceptually speaking, the ODE form will be

der(x) = x^2 + 2*u + (5*x^2)^3

and that is what gets linearized:

der(dx) = (2*x + 5^3*6*x^5)*dx + 2*du

this is not done by substitution of expressions as I showed you here, but is rather done in an efficient mixed symbolic and numerical fashion.

If the system has index > 1, e.g.

der(x1) = y;
der(x2) = -y;
x1 = x2;

then it is first brought to index 1 using Pantelides' algorithm and the dummy derivative algorithm; then, linearization is applied.

You can read about this stuff extensively in Cellier, Kofman "Continuous Time System Simulation", Springer 2006, that's an excellent introductory textbook.

[MohammadHadiAlizadeh]

Thank you. I'm closing the issue.