Unit error in Modelica.Clocked.Examples.Systems.Utilities.ComponentsThrottleControl.IntakeManifold

[qlambert-pro]

The following equation doesn't pass unit checking:

m_ao_der = -0.366 + 0.08979 * N * P_m - 0.0337 * N * P_m ^ 2 + 0.0001 * N ^ 2 * P_m;

N is expressed in "rad/s" and P_m in "bar" so - 0.0337 * N * P_m ^ 2 and 0.0001 * N ^ 2 * P_m don't have the same unit, for instance, despite being added.

[henrikt-ma]

I'll bring up one more, similar, unit error also affecting the Modelica.Clocked.Examples.Systems.EngineThrottleControl example.

In the ThrottleBody model,

ModelicaStandardLibrary/Modelica/Clocked/Examples/Systems/Utilities/ComponentsThrottleControl/ThrottleBody.mo

Line 28 in 7051091

f_Theta = 2.821 - 0.05231*Theta + 0.10299*Theta^2 - 0.00063*Theta^3; // different to paper 0.0063*Theta^3

As far as I understand, the only unit of angle that would allow this kind of equation would have been "rad", but both f_Theta and Theta have unit "deg".

Edit: Fixed typo "degC" → "deg".

[HansOlsson]

As far as I understand, the only unit of angle that would allow this kind of equation would have been "rad", but both f_Theta and Theta have unit "degC".

You mean "deg" - not "degC". To me that is one of the problematic cases: someone has fitted a 3rd degree polynomial between angle and angle (?):

• Obviously one could "fix" the input by removing unit from input-angle (divide by unit-angle in degrees), but that only documents that there's something odd going on and wouldn't help finding any issues. (The same for the IntakeManifold.)
• As for the output f_Theta I have no idea why it is an angle. Looking at the original paper it is just a throttle function f(Theta), and to me it just seems like an error to think it is an angle.

Interestingly it could be that f_Theta is seen as angle due to incorrect unit-analysis; someone saw f_Theta = 2.821 - 0.05231*Theta +...; and thought that thus it had the same unit as Theta...

Added: For the first someone it seems to be: Prabhakar, R. 1975,“Optimal and Suboptimal Control of Automotive Engine Efficiency and Emissions,” Ph.D. Thesis, Purdue University. West Lafayette, IN

[HansOlsson]

And more importantly. Why does it say: 0.00063*Theta^3; // different to paper 0.0063*Theta^3
Is it:

• A bug in the model, and someone just a wrote comment documenting it instead of correcting it?
• A typo in the paper (NONLINEAR ENGINE MODEL FOR DRIVETRAIN SYSTEM DEVELOPMENT) and someone corrected it in the model? What does the original paper say?
• An intentional change for unspecified reason?

[mestinso]

What is the current thinking/best practices for these types of equations (in the first post) and what does the spec say about it, if anything? It's not yet obvious there there is an issue here. Bear with me here as I explain my thinking and understanding...

1. It seems fairly obvious that this is just an empirical (multivariate) polynomial, which are very common in engineering. For example, the following form is very typical (full quadratic with two inputs):
   y = a + b*x1 + c*x1^2 + d*x2 + e*x2^2 + f*x1*x2
2. For such cases, it's extremely common that the units of x1, x2, and y are different from each other (same goes for the higher order terms). So, what does this mean? It means the coefficients don't just have a value, they also have units. In other words, [a]=[y], [b]=[y]/[x1], [c]=[y]/[x1]^2, etc....
3. Given this, there should be no surprise here that the units for the different terms don't match in the original issue here.
4. So from a practical engineering standpoint, there is nothing wrong here. Therefore, the real question is what is the appropriate way to handle these situations in Modelica.

I expect the answer to following question should help guide us: for the equation below, should 1.0 be assumed to have units of "1" or units of "" (which could potentially be deduced to have units of [y]/[x])?
y = 1.0*x;
I am not seeing anything in the Modelica spec that clarifies this point.

Furthermore, should there be any difference in assumption between cases 1 and 2 below?

model Case1
  Real y(unit="m");
  Real x(unit="K") = 1.0;
  Real a=1.0;
equation 
  y = a*x;
end Case1;

model Case2
  Real y(unit="m");
  Real x(unit="K") = 1.0;
equation 
  y = 1.0*x;
end Case2;

In case 1, a has units "" to start, and my understanding is that Modelica tools are expected to be able to deduce the units to be "m/K". Furthermore, I expect this isn't necessarily considered a bad practice and tools aren't going to complain.

And then in case 2, the question i posed earlier naturally arises: should the modelica tool assume 1.0 has units of "" or "1"? If the answer is that the units should be assumed to be "1", then I agree there is an issue here. But if the assumption is that the units are "" and could later be deduced, I see no such issue. At the very least I'd like to see this clearly specified in the modelica spec.

[henrikt-ma]

In the language group, we have worked quite a bit towards standardization in the area of unit checking last year. The topic is bit, but our starting point has been precisely to make it clear that the literal 1.0 below should have unit "1" in order to not ruing unit checking:

y = 1.0 * x;

See:

• Unit checking ModelicaSpecification#3255
• Handle units of literal constants by means of rules for the empty unit ModelicaSpecification#3257
• Added "four-line" proposal for units of literals ModelicaSpecification#3266

Please join us in these specification topics if you are interested, and let us try to keep these design discussions in the ModelicaSpecification repository.

[mestinso]

@henrikt-ma

Thanks. I suppose I was aware of the items you mentioned, but they are either "open" or in "draft". Shouldn't we wait until they are closed to declare items like this as an "issue". Or would you characterize the situation as the writing is on the wall, and values like 1.0 should be considered as having unit="1" and action can be taken now?

Also, can you clarify your current expectation for case 1: I assume/hope that one is not considered an issue? I think this is relevant for this issue because it might drive what good should look like/how to resolve the issue here.

For example, if I have the following equation: y = a + b*x1 + c*x1^2 + d*x2 + e*x2^2 + f*x1*x2 or something similar, personally what I have done is to declare a vector (which might come from a record) and use it like the following:

Real coefficients[6] = {[insert values here]};
equation
y = coefficients[1] + coefficients[2]*x1 + ...

The key is while I would normally declare the units for x1, x2, and y, I'm not going through and declaring the units for the 6 coefficients (leaving them as "").

[henrikt-ma]

Thanks. I suppose I was aware of the items you mentioned, but they are either "open" or in "draft". Shouldn't we wait until they are closed to declare items like this as an "issue". Or would you characterize the situation as the writing is on the wall, and values like 1.0 should be considered as having unit="1" and action can be taken now?

Yes, if you ask me. In addition to all ongoing language design discussions, I believe this principle is also going to be enforced one way or another also in Dymola.

Note that Modelica tools already have unit checking features even though there is no formal specification of what can be checked. Still, I'd say it's a good idea to pay attention to (all?!) unit problems reported by tools, as this will increase changes of finding errors, and reduce risk of having to fix these problems later when they become errors according to the standard.

[henrikt-ma]

Also, can you clarify your current expectation for case 1: I assume/hope that one is not considered an issue? I think this is relevant for this issue because it might drive what good should look like/how to resolve the issue here.

We're on a slippery slope here, with clear risk of ending up in design discussions in the wrong place. My expectation for 1 is that initial specification of unit checking will leave it undefined, allowing it to at least pass for some time. However, with a slight variation as constant Real a(unit = "") = 1.0;, a discussion at the most recent language design meeting in Lund makes me expect it to get rejected just as if the literal 1.0 had been used directly in the equation. It's tricky.

[henrikt-ma]

The key is while I would normally declare the units for x1, x2, and y, I'm not going through and declaring the units for the 6 coefficients (leaving them as "").

To make it work out unit-wise, I'd recommend creating helper variables with unit "1" for each of x1, x2, and y. The actual computation can then be performed cleanly in unit "1", and unit errors can be detected in the helper variable conversions:

protected
  constant Real x1_u(unit = "m") = 1;
  Real x1_1(unit = "1") = x1 / x1_u;
  constant Real x2_u(unit = "kg") = 1;
  Real x2_1(unit = "1") = x2 / x2_u;
  constant Real y_u(unit = "s") = 1;
  Real y_1(unit = "1") = y / y_u;
equation
  y_1 = coefficients[1] + coefficients[2]*x1_1 + ...

[HansOlsson]

The key is while I would normally declare the units for x1, x2, and y, I'm not going through and declaring the units for the 6 coefficients (leaving them as "").

To make it work out unit-wise, I'd recommend creating helper variables with unit "1" for each of x1, x2, and y. The actual computation can then be performed cleanly in unit "1", and unit errors can be detected in the helper variable conversions:

protected
  constant Real x1_u(unit = "m") = 1;
  Real x1_1(unit = "1") = x1 / x1_u;
  constant Real x2_u(unit = "kg") = 1;
  Real x2_1(unit = "1") = x2 / x2_u;
  constant Real y_u(unit = "s") = 1;
  Real y_1(unit = "1") = y / y_u;
equation
  y_1 = coefficients[1] + coefficients[2]*x1_1 + ...

I think we need to discuss this first.

On one hand it makes sense and we don't want to allow such sloppy unit-mixing everywhere, but:

• This rewrite just seems like tedious work that doesn't really improve the model quality.
• Having new variables (with one of several naming schemes) obscures the original formula.

One possibility would be to have some annotation on the equation (equations can have annotations; or one could even consider if for the entire model) stating that unit-checking is more forgiving here.

Obviously there are some details to investigate: how forgiving, what annotation, etc?

[henrikt-ma]

I agree that the language could be better at supporting this case. However, for now I suggest that we do it in the elaborate but clearly unit-correct way in the MSL, and initiate a language discussion about how to better support this case in the future.

[HansOlsson]

I agree that the language could be better at supporting this case. However, for now I suggest that we do it in the elaborate but clearly unit-correct way in the MSL, and initiate a language discussion about how to better support this case in the future.

I would prefer to do it in the other way - since we haven't completed the unit-discussion for the language and it seems logical that specifying the unit-checking goes together with a mechanism for escaping it.

In my view MSL is both a library that people can use, and a "template" for how to write models - and for correlations such as this I believe we need a readable good and consistent way to handle it.

Additionally there are risks with rewriting equations - especially to make them "unit-free" as it may hide typos (new and old).

[HansOlsson]

Just a note regarding f_Theta - as far as I can tell the paper Crossley and Cook doesn't give it a unit. (The give unit for mass "g", and pressure "bar", etc.)

Trying the proposed unit-rules indicates a problem in the model - even after removing the correlation (that we know is unit-sloppy - an alternative would be to implement support for disableUnitCheck-annotation, but simply removing the equation was faster).

However, changing f_Theta to have unit="g/s" clears up that unit-error, and correctly gives that tmp1, tmp2, g_Pm, pratio, gpratio all have unit="1" (as they are ratios).

In IntakeManifold a similar removal of m_ao_der = -0.366 + 0.08979*N*P_m - 0.0337*N*P_m^2 + 0.0001*N^2*P_m; doesn't indicate unit-errors, but that RTVmRatio has unit="bar/g".