modelica / ModelicaSpecification

Specification of the Modelica Language
https://specification.modelica.org
Creative Commons Attribution Share Alike 4.0 International
104 stars 40 forks source link

What is the semantics of delay(<expr>,0)? #3287

Open casella opened 2 years ago

casella commented 2 years ago

Section 3.7.4 defines the delay() operator as follows:

delay(𝑒π‘₯π‘π‘Ÿ, π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’, π‘‘π‘’π‘™π‘Žπ‘¦π‘€π‘Žπ‘₯) delay(𝑒π‘₯π‘π‘Ÿ, π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’)

Evaluates to 𝑒π‘₯π‘π‘Ÿ(time - π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’) for πšπš’πš–πšŽ>πšπš’πš–πšŽ.πšœπšπšŠπš›πš+π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’ and 𝑒π‘₯π‘π‘Ÿ(time.start) for πšπš’πš–πšŽβ‰€πšπš’πš–πšŽ.πšœπšπšŠπš›πš+π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’. The arguments, i.e., 𝑒π‘₯π‘π‘Ÿ, π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’ and π‘‘π‘’π‘™π‘Žπ‘¦π‘€π‘Žπ‘₯, need to be subtypes of Real. π‘‘π‘’π‘™π‘Žπ‘¦π‘€π‘Žπ‘₯ needs to be additionally a parameter expression. The following relation shall hold: 0β‰€π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’β‰€π‘‘π‘’π‘™π‘Žπ‘¦π‘€π‘Žπ‘₯, otherwise an error occurs. If π‘‘π‘’π‘™π‘Žπ‘¦π‘€π‘Žπ‘₯ is not supplied in the argument list, π‘‘π‘’π‘™π‘Žπ‘¦π‘‡π‘–π‘šπ‘’ needs to be a parameter expression.

According to the text, delayTime can actually be zero, and in that case, delay(<expr(time)>, 0) evaluates to <expr(time)>. This has consequences on the structural properties of the system. For example, consider this MWE:

model M
  Real x,y;
  parameter Real T;
equation
  x + y = sin(time);
  x -  delay(y, T) = 0;
end M;

If T > 0, x can be computed from the second equation, fetching the delayed value of y from a buffer, then the first equation can be solved for y.

If T = 0, according to a strict interpretation of the specification text, during simulation delay(y, T) should be the same as y, so x and y would be part of the same strong component and would need to be solved simultaneously.

This would actually make T a structural parameter.

The non-normative text of 3.7.4.1 muddies the waters instead of clarifying this point, as it says:

In principle, delay could break algebraic loops.

I guess this would mean that if I write delay(<expr>, 0), this would be computed using the last computed values of the variables depends upon, without considering the structural dependency of the output on the input when computing the incidence matrix, BLT, etc. This would be similar to the (in)famous "memory block" of Simulink. It would be dangerous to use, but it could make some sense, as long as the influence of the output of delay() on the system dynamics is weak in some sense (more on this later).

But then, the text continues with:

For simplicity, this is not supported

What does this mean? Not supported by whom? Does this mean "tools should not support it". And then, what exactly is not allowed or should not be considered? This is quite unclear.

because the minimum delay time has to be given as additional argument to be fixed at compile time.

What is "the minimum delay time"? The normative text mentions an optional maximum delay time input, but there is no mention of a "minimum delay time". Furthermore, the specification explicitly allows delayTime to be zero, so it really seems there is no such thing. I thought this may be a leftover of previous versions of the MLS, but it turns out this text has remained unchanged since Modelica 1.1 (Modelica 1.0 did not contain an explicit specification of the delay() operator, just mentioned the need of it).

Furthermore, the maximum step size of the integrator is limited by this minimum delay time in order to avoid extrapolation in the delay buffer.]

Ditto, what is "this minimum delay time" is not clear.

I think the non-normative part should be rewritten, once we've clarified the following issues:

  1. is delayTime = 0 allowed? I think the answer should be yes, since this has been in the specification forever, and there are libraries that make use of this feature (see further down), so disallowing this would break backwards compatibility
  2. what is the semantics of the operator in that case? Should it be interpreted strictly (i.e. delay(<expr>,0) could be substituted by <expr> in the Modelica code, producing exactly the same results), or should it be interpreted like the "memory block" of Simulink?

I am definitely in favour of a strict interpretation, i.e. delay(<expr>,0) is identical to <expr>, for the following reason. Suppose one wants to implement a FOPTD model, with transfer function G(s) = exp(-tau*s)/(1 + sT):

model FOPTD
  parameter SI.Time tau = 0 "Time delay";
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  T*der(y) + y = delay(u, tau);
end FOPTD;

When instantiating the model, by default tau = 0, so you get a standard first order low-pass PT1 filter. It would be weird that in this case a "memory" block was placed between the input and the output, since the model is strictly proper, the output is a state, so there is no need to add any decoupling with an infinitesimal delay. In fact, this could cause serious harm to the accuracy of the simulation results, since the ODE solver will be missing important structural information, for no apparent reason.

If there is agreement on this interpretation, we should rewrite the non-normative text accordingly.

casella commented 2 years ago

Related to #2953, though the scope is a bit different.

casella commented 2 years ago

Adding @AnHeuermann and @kabdelhad to the loop.

henrikt-ma commented 2 years ago

This was resolved in #3245.

casella commented 2 years ago

OK, I was not aware of that. It definitely makes sense in those cases where delayTime is a time-varying expression. On the other hand, it introduces a significant non-backwards compatibility, that could maybe be avoided by specifiying the meaning of delay(<expr>, <param_expr>) to be <expr> if param_expr = 0.

Does this mean that the implementation of FOPTD should now be something like this:

model FOPTD
  parameter SI.Time tau = 0 "Time delay";
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  T*der(y) + y = if tau > 0 then delay(u, tau) else u;
end FOPTD;

or maybe something like this, if you want to still be able to change a nonzero value of tau at runtime without recompiling:

model FOPTD
  parameter SI.Time tau = 0 "Time delay";
  parameter Boolean positive_tau = tau > 0 annotation(Evaluate = true);
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  T*der(y) + y = if positive_tau then delay(u, tau) else u;
end FOPTD;

All things considered, this makes a lot of sense, though the original formulation I proposed is of course a lot more compact.

Comments?

henrikt-ma commented 2 years ago

OK, I was not aware of that. It definitely makes sense in those cases where delayTime is a time-varying expression. On the other hand, it introduces a significant non-backwards compatibility, that could maybe be avoided by specifiying the meaning of delay(<expr>, <param_expr>) to be <expr> if param_expr = 0.

That would force unnecessary evaluation of parameters just to be able to determine whether the delay is zero at translation time. The unnecessary evaluation of parameters could be avoided by making a special case for a constant expression of value 0, but I doubt that the added complexity of such a specification would be worth it. I mean, in order to be useful the delayTime is unlikely to just be the literal 0, and the question is whether constant expression of value 0 captures any cases of interest, or if one would need to go for an even more obscure formulation.

henrikt-ma commented 2 years ago

Comments?

Your example points to a problem we haven't considered; this expression is illegal for tau = 0:

if tau > 0 then delay(u, tau) else u

That is, we need to ensure there is a way to guard against invalid use of delay(…, 0). (The solution cannot be ST.Time tau(min = Constants.small) to make the if-condition a constant expression, as the purpose is to also handle the structurally different case of tau = 0.)

I couldn't see an easy way to describe what would count as a structurally disabled delay expression (where the constraint on delayTime would be lifted), and this leads be back to a small variation on your idea, @casella:

model FOPTD_constant
  constant SI.Time tau(min = 0) = 0 "Time delay"; // Constant in order to allow zero delay.
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  T * der(y) + y = delay(u, tau);
end FOPTD_constant;

That is, I now realize that it is useful to have a special rule for delayTime being a constant expression of value 0.

With a more complicated rule we could also support this:

model FOPTD_evaluated
  parameter SI.Time tau(min = 0) = 0 "Time delay" annotation(Evaluate = true); // Evaluated in order to allow zero delay.
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  T * der(y) + y = if tau > 0 then delay(u, tau) else u;
end FOPTD_evaluated;

The rule would then looks something like this:

Except in structurally disabled branches of \lstinlin!if!-expressions and \lstinline!if!-equations, the following relation shall hold: $0 < \mathit{delayTime} \leq \mathit{delayMax}$, otherwise an error occurs. Here, a \emph{structurally disabled} branch means a branch that at translation-time can be determined to never be active, and the presence of \lstinline!delay($\ldots$, $p$)! in a branch should not trigger evaluation of $p$ by itself.

casella commented 2 years ago

Your example points to a problem we haven't considered; this expression is illegal for tau = 0:

if tau > 0 then delay(u, tau) else u

Would this be OK?

model FOPTD
  parameter SI.Time tau = 0 "Time delay";
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  if tau > 0 then
    T*der(y) + y = if tau > 0 then delay(u, tau) else u;
  else
    T*der(y) + y = u;
end FOPTD;

And what about this?

model FOPTD
  parameter SI.Time tau = 0 "Time delay";
  parameter Boolean positive_tau = tau > 0 annotation(Evaluate = true);
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  if positive_tau then
    T*der(y) + y = delay(u, tau) 
  else
    T*der(y) + y = u;
  end if;
end FOPTD;
casella commented 2 years ago

That is, I now realize that it is useful to have a special rule for delayTime being a constant expression of value 0.

"constant expression" won't be very useful, because you can't change easily that upon instantiation :)

"parameter expression" would be ok, but then we have to specify that when the expression has value zero, the operator returns exactly its first input, without any kind of decoupling or infinitesimal delay implied.

With a more complicated rule we could also support this:

model FOPTD_evaluated
  parameter SI.Time tau(min = 0) = 0 "Time delay" annotation(Evaluate = true); // Evaluated in order to allow zero delay.
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  T * der(y) + y = if tau > 0 then delay(u, tau) else u;
end FOPTD_evaluated;

The rule would then looks something like this:

Except in structurally disabled branches of \lstinlin!if!-expressions and \lstinline!if!-equations, the following relation shall hold: 0<delayTime≀delayMax, otherwise an error occurs. Here, a \emph{structurally disabled} branch means a branch that at translation-time can be determined to never be active, and the presence of \lstinline!delay($\ldots$, p)! in a branch should not trigger evaluation of p by itself.

LGTM.

henrikt-ma commented 2 years ago

That is, I now realize that it is useful to have a special rule for delayTime being a constant expression of value 0.

"constant expression" won't be very useful, because you can't change easily that upon instantiation :)

I didn't say final constant, just constant:

FOPTD_constant myFoptd(tau = 1); /* Not a problem. */

But yes, it is very frustrating that we still can't set Dialog(visisble = true) for tau after having decided on a solution back in November 2018 (#2211).

"parameter expression" would be ok, but then we have to specify that when the expression has value zero, the operator returns exactly its first input, without any kind of decoupling or infinitesimal delay implied.

I find it misleading to make it a parameter expression when we know that the semantics require it to be evaluated. I think we should look for a solution that avoids evaluation of parameters, and it seems like we might be able to work something out…

henrikt-ma commented 2 years ago

Would this be OK?

model FOPTD
  parameter SI.Time tau = 0 "Time delay";
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  if tau > 0 then
    T*der(y) + y = if tau > 0 then delay(u, tau) else u;
  else
    T*der(y) + y = u;
end FOPTD;

I would like to not rely on this being a valid guard against evaluating the delay for tau = 0. I find it subtle enough that unbalanced if-equation branches lead to evaluation of the conditions. (I would appreciate the convenience of always being able to fill in the delay buffers of all delay expressions that haven't been completely eliminated from the integration problem.)

henrikt-ma commented 2 years ago

And what about this?

model FOPTD
  parameter SI.Time tau = 0 "Time delay";
  parameter Boolean positive_tau = tau > 0 annotation(Evaluate = true);
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  if positive_tau then
    T*der(y) + y = delay(u, tau) 
  else
    T*der(y) + y = u;
  end if;
end FOPTD;

This would be OK according to the structurally disabled branches proposal, but wouldn't have any advantage over

parameter SI.Time tau = 0 "Time delay" annotation(Evaluate = true);
…
equation
  if tau > 0 then
    T*der(y) + y = delay(u, tau); /* Structurally disabled when tau = 0 */
casella commented 2 years ago

I didn't say final constant, just constant:


FOPTD_constant myFoptd(tau = 1); /* Not a problem. */

The problem is, when I double-click on the block to set its parameters in a GUI, tau won't show up in the parameter window, making this solution pretty useless from a practical point of view.

But yes, it is very frustrating that we still can't set Dialog(visisble = true) for tau after having decided on a solution back in November 2018 (#2211).

I understand #2211 discusses parameters, not constants. Do I miss something?

I find it misleading to make it a parameter expression when we know that the semantics require it to be evaluated. I think we should look for a solution that avoids evaluation of parameters, and it seems like we might be able to work something out…

I know this is less than ideal, but we have exactly the same problem, e.g., with parameters setting the number of elements in an array. I think we have to live with it. Or postpone the solution of this issue more or less indefinitely πŸ˜….

Besides, evaluating the parameter when it has zero value at runtime is the simplest way to implement the proposed behaviour, but is by no means always necessary. For example, if the delay() operator is not breaking algebraic loops, one could devise a smart enough runtime that can handle zero and nonzero values correctly. No big deal.

casella commented 2 years ago

The specification could say something like:

if delayTime is a constant or parameter expression, then delay(<expr>, delayTime) should behave exactly as <expr> if the value of delayTime is zero. Otherwise, delayTime should be greater than zero, and the tool can assume that its output does not instantaneously depend on <expr>.

How this gets implemented is besides the scope of the specification. Yes, one simple implementation is to evaluate it, which prevents changing it at runtime. More sophisticated implementations could be devised that allow to go from zero to non-zero values without recompiling. This is purely a tool issue.

BTW, this specification would be backwards compatible in most cases that are easily handled, and only prevent tricky cases like the one that was put forth in #3245, which would most likely fail in all tools anyway. The current formulation after #3245 is a lot more non-backwards compatible, probably for no compelling reason.

henrikt-ma commented 2 years ago

It is also possible to avoid the need to evaluate tau (this is my personal favorite):

model FOPTD_unevaluated
  constant Boolean zeroDelay = false;
  parameter SI.Time tau(min = Constants.small) = 0 "Time delay, if not zeroDelay" annotation(Dialog(enable = not zeroDelay));
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  /* The delay expression below should be structurally eliminated if zeroDelay. */
  T * der(y) + y = if zeroDelay then u else delay(u, tau); 
end FOPTD_unevaluated;
casella commented 2 years ago

But you'll have to make zeroDelay a parameter with Evaluate = true, otherwise you won't get it in the GUI, as mentioned above πŸ˜…

henrikt-ma commented 2 years ago

But yes, it is very frustrating that we still can't set Dialog(visisble = true) for tau after having decided on a solution back in November 2018 (#2211).

I understand #2211 discusses parameters, not constants. Do I miss something?

Yes, that it is also very useful for constants, especially since tools have a tendency to not make constants visible by default. (It may be useful with other variability prefixes as well, for instance to insist on showing a binding equation even though it is final and cannot be modified.)

henrikt-ma commented 2 years ago

The current formulation after #3245 is a lot more non-backwards compatible, probably for no compelling reason.

But it is only a hypothetical backwards incompatibility unless there was actually a tool that had a working implementation for a zero delay. If you are aware of such a tool, that would be very important input to the discussion!

henrikt-ma commented 2 years ago

How this gets implemented is besides the scope of the specification. Yes, one simple implementation is to evaluate it, which prevents changing it at runtime. More sophisticated implementations could be devised that allow to go from zero to non-zero values without recompiling. This is purely a tool issue.

There are many things we could define and say it's purely a tool issue whether it works or not, but the more constructs we have that tools treat differently, the more damage to the famous Modelica model portability between tools. Therefore, if we can define semantics for delay that are easy enough to support fully in all tools, this will be beneficial for model portability.

(On a side note, I believe that System Modeler is actually in rather good shape for being able to go for the "more sophisticated implementation", but I'm not sure users would actually appreciate the complexity this adds to the simulation execution; exponential growth of the number of cases to consider, longer model translation time, bigger size of simulation executables, more difficult to debug simulation issues, etc.)

HansOlsson commented 2 years ago

In general, if predictate(x) then foo(x) else bar(x) should only evaluate foo(x) is predicate(x) is true.

For delay it is a different, and I would say that we should not allow delay at all in conditional code, unless the condition is structural. Basically, delay (and spatialDistribution) implicitly requires that we store the expression (note that it takes an expression as first argument not a variable) - and to store that expression we need to evaluate it, which breaks the normal logic that conditional code isn't executed.

Compare if p>0 then sqrt(p)*x else 0 with if p>0 then delay(sqrt(p)*x, 2) else 0.

henrikt-ma commented 2 years ago

But you'll have to make zeroDelay a parameter with Evaluate = true, otherwise you won't get it in the GUI, as mentioned above πŸ˜…

I don't need to, because I have Dialog(__Wolfram_show = true) while waiting for #2211. :)

henrikt-ma commented 2 years ago

For delay it is a different, and I would say that we should not allow delay at all in conditional code, unless the condition is structural.

You are right, that's a good observation. I suppose the best way to phrase it might actually be something similar to what we have for unbalanced if-equations, and explain non-normatively that this gives the tools the freedom to evaluate the conditions as an easy way of ensuring that guarded expressions don't get incorrectly evaluated. Such a formulation allows a tool to leave conditions unevaluated if it can guarantee that evaluation of the expression to be delayed is safe.

To make life harder for tool vendors one could ease up the requirement so that branching conditions don't need to be evaluable if the expression to be delayed is just a component reference (covering a small but important subset of what could be evaluated safely).

casella commented 2 years ago

The current formulation after #3245 is a lot more non-backwards compatible, probably for no compelling reason.

But it is only a hypothetical backwards incompatibility unless there was actually a tool that had a working implementation for a zero delay. If you are aware of such a tool, that would be very important input to the discussion!

Dymola handles it without a hiccup. OMC didn't (now it does, according to what I am proposing here), and that's why we are here discussing this issue, to make sure it is not a tool issue πŸ˜ƒ .

henrikt-ma commented 2 years ago

Dymola handles it without a hiccup. OMC didn't (now it does, according to what I am proposing here), and that's why we are here discussing this issue, to make sure it is not a tool issue πŸ˜ƒ .

OK. In what sense does it work; with evaluation of parameters, by dynamically changing structure at initialization time, or in some other sense?

To make it concrete, how is this model handled?

model PossiblyZeroDelay
  parameter Real delayTime(min = 0) = 1.0;
  Real x(start = 1.0);
  Real y = delay(x, delayTime);
equation
  x * y = exp(-time);
  annotation(experiment(StopTime = 5.0));
end PossiblyZeroDelay;

Is x expected to be determined by a nonlinear system of equations in both x and y, just in order to capture the structure in case delayTime is zero? (I guess that the example can be modified to get even more drastic consequences on the equation system structure.)

casella commented 2 years ago

@henrikt-ma I have no idea how Dymola handles this internally.

What I know is that we are dealing with a library that uses blocks similar to my FOPTD block, which most of the time have tau = 0, which apparently gave no problem in Dymola. When we started using the library in OMC, we were getting weird numerical issues with Jacobians in some cases, which led the solver to take millions of steps instead of just a few tens of thousands for a simulation. Until I just removed the delay(expr,0) operator, and at that point everything started to work fine, exactly as in Dymola. So, we solved the problem by tentatively evaluating the parameter expression at compile time, and removing the delay operator in case delayTime evaluates to zero.

Anyway, from a library developer's perspective I see these very clear requirements

How this is actually implemented, and to which extent you can change the delayTime parameter at runtime without recompiling the model, these are all tool issues IMHO. I would just like to make sure that these requirements are well understood by tool developers, so that I can be confident that I will actually get no performance penalty when delayTime = 0, regardless of the tool I'm using. I don't think the text of the specification should get into actual implementation details, that's beyond the scope of the language spec.

casella commented 2 years ago

To make it concrete, how is this model handled?

model PossiblyZeroDelay
  parameter Real delayTime(min = 0) = 1.0;
  Real x(start = 1.0);
  Real y = delay(x, delayTime);
equation
  x * y = exp(-time);
  annotation(experiment(StopTime = 5.0));
end PossiblyZeroDelay;

If delayTime > 0, I expect that y is a know variable at each time step, because it refers to past computed values that can be retrieved from a suitable buffer, so the equation is solved for x in this way: x := exp(-time)/y.

If delayTime = 0, then x = y, so I guess it could be solved as x := sqrt(exp(-time)), or possibly via tearing and Newton iterations, with

y := x;
res[1] = x*y - exp(-time);

Whether you can change the value of delayTime at runtime without recompiling the model is just a tool issue. For example, a tool could tentatively evalutate delayTime at compile time; if it is zero, it will constant-evaluate the parameter, remove the delay operator and will not allow to change the parameter at runtime; if it is > 0, it may allow to change it at runtime, issuing a runtime error if it is changed to zero.

Or, it could generate runtime code that handles both cases at runtime, without the need of recompiling anything.

I don't think we should bother about these details in the language specification.

henrikt-ma commented 2 years ago

How this is actually implemented, and to which extent you can change the delayTime parameter at runtime without recompiling the model, these are all tool issues IMHO. I would just like to make sure that these requirements are well understood by tool developers, so that I can be confident that I will actually get no performance penalty when delayTime = 0, regardless of the tool I'm using. I don't think the text of the specification should get into actual implementation details, that's beyond the scope of the language spec.

I don't find that a viable principle for Modelica language design decisions. To be the technology we want it to be, we have to design the language so that it becomes reasonable to implement it in tools. Suggesting that tools could handle this with a runtime error in case the parameter is set to zero at initialization is an indication to me that it might not be a reasonable implementation effort to do it properly.

Also, I believe that unpredictable evaluation of parameters is one of the central reasons for poor portability of models between tools. To tentatively evaluate and make a decision based on the value sounds like a new level of non-portability.

For these reasons, I still think that the realistic way of specifying this so that tools have a fair chance of implementing it correctly (without the need to resort to always evaluating the delayTime), is to make the change of structure explicit in the model; require that delayTime is positive, and use a separate parameter/constant for selecting between the structurally different cases. In addition to making the change of structure explicit, it also makes it clear that there is no need to evaluate delayTime parameter expression.

Thinking more about how to apply this design, I'm changing the role and renaming the structure-controlling parameter (which I now make Evaluate = true in order to avoid discussions about visibility of constants in GUIs):

block FOPTD_unevaluated2
  parameter Boolean enableDelay = false annotation(Evaluate = true);
  parameter SI.Time tau(min = Constants.small) = 1.0 "Time delay, if enableDelay" annotation(Dialog(enable = enableDelay));
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  /* The delay expression below should be structurally eliminated if not enableDelay. */
  T * der(y) + y = if enableDelay then delay(u, tau) else u;
end FOPTD_unevaluated2;
henrikt-ma commented 2 years ago

Here's a more complete example that shows the drawback of allowing a delayTime parameter to take the value zero.

Consider a big closed-loop system where both plant and controller have direct feed-through, and the sensor is modeled with a small time delay in the form of the FixedDelay block. In essence, the block looks like this:

block FixedDelay "Delay block with fixed DelayTime"
  extends Modelica.Blocks.Interfaces.SISO;
  parameter Modelica.SIunits.Time delayTime(start = 1);
equation
  y = delay(u, delayTime);
end FixedDelay;

Now, if an implementation does the ambitious thing and doesn't fail with an feature-not-implemented-fully error if the delayTime is set to zero after translation, this means that the translated model must include the complicated equation systems for a closed loop with direct feed-though, in addition to the simpler equation systems of the "broken" loop. This has several disadvantages for a user who is only interested in the "broken" loop behavior:

I don't think we can expect users to understand that they would need to specify a positive min for the delay parameter to avoid this:

FixedDelay sensorDelay(delayTime(min = Modelica.Constants.small) = 0.01);

For these reasons, it is better to have rules that require a user to make an active translation-time choice between the zero and non-zero delay behaviors.

casella commented 2 years ago

I don't find that a viable principle for Modelica language design decisions. To be the technology we want it to be, we have to design the language so that it becomes reasonable to implement it in tools.

No doubt.

Suggesting that tools could handle this with a runtime error in case the parameter is set to zero at initialization is an indication to me that it might not be a reasonable implementation effort to do it properly.

I'm not following you here. A trivial implementation will cause a runtime error if you change the value from non-zero to zero at runtime. The error will just say: please set the parameter to zero and recompile. No big deal. A less trivial implementation may take care of this without the need of recompiling. Notice that the language specification doesn't event refer to a "compilation" and a "runtime" stages. One could very well build an interpreted Modelica tool.

For these reasons, I still think that the realistic way of specifying this so that tools have a fair chance of implementing it correctly (without the need to resort to always evaluating the delayTime), is to make the change of structure explicit in the model; require that delayTime is positive, and use a separate parameter/constant for selecting between the structurally different cases. In addition to making the change of structure explicit, it also makes it clear that there is no need to evaluate delayTime parameter expression.

This is a perfectly legitimate opinion. As in the case of our never ending discussion on units of literals, it has the disadvantage of being non-backwards compatible. After 20+ years of Modelica code development, I am strongly in favour of preserving backwards compatibility, unless there are really compelling reasons not to do so. The reason is that a lot of effort went into developing and testing libraries, and breaking them causes an untold amount of waste of time. We shouldn't do that lightly. If we were starting to develop a new language from scratch, I would certainly follow your advice. But we aren't.

Precisely as in that case, what I'm trying to do here is to find a backwards compatible solution. In any case, I believe that some backwards-compatible improvement is much better than no improvement at all.

Thinking more about how to apply this design, I'm changing the role and renaming the structure-controlling parameter (which I now make Evaluate = true in order to avoid discussions about visibility of constants in GUIs):

block FOPTD_unevaluated2
  parameter Boolean enableDelay = false annotation(Evaluate = true);
  parameter SI.Time tau(min = Constants.small) = 1.0 "Time delay, if enableDelay" annotation(Dialog(enable = enableDelay));
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  /* The delay expression below should be structurally eliminated if not enableDelay. */
  T * der(y) + y = if enableDelay then delay(u, tau) else u;
end FOPTD_unevaluated2;

Writing this kind of code would also be possible with my proposal. The problem is that by making this design compulsory we are deprecating libraries that could be perfectly functional in all tools if we made the semantics of zero delay explicit. Currently, they are perfectly functional, but only in some tools that made this choice implicitly. I would like to make it explicit so they become portable. As in the case of units of literals, I don't want to solve all the problems of Modelica, just trying to make a very modest advance of the state of the art, without breaking things πŸ˜….

casella commented 2 years ago

Here's a more complete example that shows the drawback of allowing a delayTime parameter to take the value zero.

Consider a big closed-loop system where both plant and controller have direct feed-through, and the sensor is modeled with a small time delay in the form of the FixedDelay block. In essence, the block looks like this:

block FixedDelay "Delay block with fixed DelayTime"
  extends Modelica.Blocks.Interfaces.SISO;
  parameter Modelica.SIunits.Time delayTime(start = 1);
equation
  y = delay(u, delayTime);
end FixedDelay;

Now, if an implementation does the ambitious thing and doesn't fail with an feature-not-implemented-fully error if the delayTime is set to zero after translation, this means that the translated model must include the complicated equation systems for a closed loop with direct feed-though, in addition to the simpler equation systems of the "broken" loop. This has several disadvantages for a user who is only interested in the "broken" loop behavior:

* Longer time to translate model.

* Bigger size of simulator.

* More difficult to understand equation system structure of translated model.

* Risk of encountering warnings and errors that only affect the direct-feedthrough loop.

Well, first of all it is well-know among control engineers that a closed-loop model with direct feed-through in the loop transfer function is not very meaningful, because it is neglecting high-frequency dynamics that has a significant chance of being crucial to determine the closed-loop stability. Hence, in your example, the case with zero delay would not be very interesting. It is interesting from a purely academic point of view for a mathematician, but not for a control engineer.

To the contrary, I believe that in all well-posed models, this problem of having a conditional algebraic loops depending on delay being zero or non-zero is purely theoretical, but will hardly ever present itself in practice. So I won't be much worried by that. The only structural issue that must absolutely be considered regards the Jacobian of the RHS of the ODEs. Consider this MWE:

model M
  Real x1, x2;
  parameter Real tau(min = 0);
equation
  der(x1) = x2;
  der(x2) = delay(x1, tau);
end M;

The Jacobian is [0, 1; 1, 0] if tau = 0, otherwise it is [0, 1, 0, 0]. This seems not such a big deal to support. I just wanted to make this more explicit in the specification, rather than leaving it to the good will of tool implementors. With the supreme goals of a) backwards compatibility and b) tool portability.

Apparently, the problems we had in OMC could be traced to the fact that we were not generating the right Jacobian because we missed some dependencies because of delay, and we didn't consider that the delay was actually not there. So, this is not a purely theoretical issue, but actually happened in practical industrial usage.

HansOlsson commented 1 year ago

Here's a more complete example that shows the drawback of allowing a delayTime parameter to take the value zero. Consider a big closed-loop system where both plant and controller have direct feed-through, and the sensor is modeled with a small time delay in the form of the FixedDelay block. In essence, the block looks like this:

block FixedDelay "Delay block with fixed DelayTime"
  extends Modelica.Blocks.Interfaces.SISO;
  parameter Modelica.SIunits.Time delayTime(start = 1);
equation
  y = delay(u, delayTime);
end FixedDelay;

Now, if an implementation does the ambitious thing and doesn't fail with an feature-not-implemented-fully error if the delayTime is set to zero after translation, this means that the translated model must include the complicated equation systems for a closed loop with direct feed-though, in addition to the simpler equation systems of the "broken" loop. This has several disadvantages for a user who is only interested in the "broken" loop behavior:

* Longer time to translate model.

* Bigger size of simulator.

* More difficult to understand equation system structure of translated model.

* Risk of encountering warnings and errors that only affect the direct-feedthrough loop.

Well, first of all it is well-know among control engineers that a closed-loop model with direct feed-through in the loop transfer function is not very meaningful, because it is neglecting high-frequency dynamics that has a significant chance of being crucial to determine the closed-loop stability. Hence, in your example, the case with zero delay would not be very interesting. It is interesting from a purely academic point of view for a mathematician, but not for a control engineer.

To the contrary, I believe that in all well-posed models, this problem of having a conditional algebraic loops depending on delay being zero or non-zero is purely theoretical, but will hardly ever present itself in practice. So I won't be much worried by that. The only structural issue that must absolutely be considered regards the Jacobian of the RHS of the ODEs. Consider this MWE:

model M
  Real x1, x2;
  parameter Real tau(min = 0);
equation
  der(x1) = x2;
  der(x2) = delay(x1, tau);
end M;

The Jacobian is [0, 1; 1, 0] if tau = 0, otherwise it is [0, 1, 0, 0].

Theoretically yes, but if you approximate the delay (as is often necessary) it may be [0, 1; min(1, tau), 0].

Apparently, the problems we had in OMC could be traced to the fact that we were not generating the right Jacobian because we missed some dependencies because of delay, and we didn't consider that the delay was actually not there. So, this is not a purely theoretical issue, but actually happened in practical industrial usage.

Hmm... It seems you did some special handling of delay to ignore it.

casella commented 1 year ago

Theoretically yes, but if you approximate the delay (as is often necessary) it may be [0, 1; min(1, tau), 0].

I'm sorry, I can see why it should be so. According to the specification, Section 3.7.4.1, delay is approximated by interpolating past buffered values.

If, e.g., tau = 0.9, why should be the partial derivative of der(x2) w.r.t. x1 be 0.9? If any change in x1 is only seen on the derivative of x2 after 0.9 seconds, IMHO that derivative should definitely be zero.

My point is, if there is a delay, the coupling between the ODEs that have delays in-between takes place exclusively through the delay buffers, so the implicit solver is not aware of that. However, if there is zero delay, it better should.

Apparently, the problems we had in OMC could be traced to the fact that we were not generating the right Jacobian because we missed some dependencies because of delay, and we didn't consider that the delay was actually not there. So, this is not a purely theoretical issue, but actually happened in practical industrial usage.

Hmm... It seems you did some special handling of delay to ignore it.

See the above comment: if tau > 0, delay(x1,tau) is just like an exogenous input, coming out of the buffer. It loses any explicit relationship with state variable x1, so the element of the Jacobian for an implicit solver should be zero. Of course if the delay is very short, the solver will eventually need to decrease the step size below that to get meaningful results, but that's inevitable if you have short but non-zero delays in your model.

Do I miss something?

HansOlsson commented 1 year ago

Theoretically yes, but if you approximate the delay (as is often necessary) it may be [0, 1; min(1, tau), 0].

I'm sorry, I can see why it should be so. According to the specification, Section 3.7.4.1, delay is approximated by interpolating past buffered values.

But if the step-size exceeds the delay-time that is problematic and limiting the step-size to a very short delay-time as indicated in the specification will be inefficient. You can either extrapolate instead (which has some numeric issues) or use an implicit method where you interpolate the previous iterate for buffered value at the end of the step.

casella commented 1 year ago

I see what you mean. One could, e.g., use J = [0, 1; 1 0] when the step size h > tau and J = [0, 1; 0, 0] when h < tau.

The truth is, variable-step size integrators with error control assume that the system is finite-dimensional and there are no "hidden states" (e.g. delay lines or other types of memory) embedded in the RHS of the ODE. As soon as we have that, the error estimator's outputs must be taken with a grain (or two) of salt, I guess. Every time we simulate systems with delay, we are stretching solvers beyond the theory and there is some heuristics involved.

In any case, how do you deal efficiently with systems with small delay is a quality of implementation issue. Allow me to come back to the original topic of this issue before we drift off too much (I kind of see a recurring pattern here πŸ˜…).

There are three issues in the non-normative text:

After this very long discussion, do we have some ideas to rewrite the non-normative text to make it clearer and consistent with the normative text, without changing the normative part, which would introduce non-backwards compatible changes?

As usual, the goal of this ticket is not to change the specification in a fundamental, non-backwards compatible way (that would be another ticket), but just to improve parts which are inconsistent or unclear.

I can give it a try.

henrikt-ma commented 1 year ago
  • it does not make clear what is the semantics of delay when delayTime = 0. In fact, IMHO it muddies the waters by talking about a minimum delay time, whereby the normative part states 0 <= delayTime, so it explicitly allows zero delay.

As of https://github.com/modelica/ModelicaSpecification/issues/3245 it should be clear that delayTime = 0 is not allowed, and the reason for this is that allowing it makes things so very complicated. It would only be counter-productive to specify that this is allowed if we cannot reasonably expect tools to handle this case well, and I still haven't seen a good proposal for how to handle it in the general case of a parametric delayTime. I also don't think that a special rule for a constant delayTime of zero would be that useful in practice.

What I've seen in the discussion above, however, is that we're lacking some rules for how delay may be combined with conditional constructs, and if we give such rules one could consider allowing any value of delayTime as long as the delay expression is statically eliminated. (It could be handy to allow leaving a delayTime parameter at a default of zero as long as the delay is statically eliminated, meaning that the user will be forced to select a positive value only if the delay is statically activated.)

casella commented 1 year ago

As of #3245 it should be clear that delayTime = 0 is not allowed, and the reason for this is that allowing it makes things so very complicated.

I was not aware of #3245 when we started this discussion (I missed that, my fault), and I don't think it was a good decision, because it breaks backwards compatibility. We should not do that light-heartedly.

Restricting delayTime to be > 0 is perfectly fine if delayTime is a time-varying variable or expression, because that is a nightmare to implement and it's not worth going through it, I agree with that. In fact, since nobody implemented it, there is no library that needs this feature, so the problem of breaking non-backwards compatibility by introducing this restriction does not exist in practice.

On the other hand, restricting delayTime to be > 0 is overly restrictive if it is a parameter or constant expression, which can be handled, and it will break backwards compatibility on existing libraries that use this feature, which IMHO is really bad, because it just annoys people unnecessarily, making them waste time to "fix" code that was previously perfectly fine. I believe one of the strong points of Modelica has been its remarkable stability and backwards compatibility over 25 years, which promoted code reusability. Other modelling and simulation languages do not have this good track record.

There was not much discussion in #3245, just @HansOlsson writing "That makes sense". Maybe it does, but I'd like to have a bit more feedback from the language group (and even from the MA, which ultimately votes the new releases of MLS) before we introduce non-backwards compatible changes in the language, that make existing libraries illegal.

It would only be counter-productive to specify that this is allowed if we cannot reasonably expect tools to handle this case well

I absolutely agree. The key is to introduce the restriction that delayTime > 0 only for time-varying delayTime, not for the case of delayTime being a parameter (or constant) expression.

There is a very straightforward way to handle this.

If delayTime is a parameter or constant expression, like in the FOPTD model

model FOPTD
  parameter SI.Time tau(min = 0) = 0 "Time delay";
  parameter SI.Time T "Time constant";
  input Real u;
  output Real y;
equation
  T*der(y) + y = delay(u, tau);
end FOPTD;

the tool tentatively evaluates it at compile time, and if it is zero it just removes the delay operator outright. If it is not, then the delay operator is retained, and changing delayTime to zero at runtime could trigger a runtime error, if the runtime cannot handle that.

That's it.

We already do it in OMC and it works very well. In fact, there is also something that works equally well in Dymola, though I have no idea how is it implemented. So we have two working implementations, which is the normal requirement for a Modelica change proposal.

As far as I understand, this proposal provides exactly the same functionality of #3245 in all cases where delayTime > 0. Additionally, it allows to handle the FOPDT model (which is legal in MLS 3.5!) also for the zero delay case, which is not possible with #3245. And this is possible with a straightforward implementation, that could be mentioned in the non-normative text.

So, IMHO, the proposal to allow delayTime >= 0 when delayTime it is a parameter or constant expression is dominant (in the Pareto sense) with respect to #3245: it covers more cases, and it doesn't do any harm in all cases that are already handled by #3245. Do I miss something?

I also don't think that a special rule for a constant delayTime of zero would be that useful in practice.

This is your personal opinion. We are working with a library that needs this feature, I find it perfectly reasonable according to the currently released MLS 3.5, and I think it won't be good if it is made illegal by #3245 when MLS 3.6 is eventually released. There are probably other libraries in the same situation, since the need of describing dynamic behaviour with optional delay is something rather common in systems and control.

HansOlsson commented 1 year ago

It seems we need to revisit this. To me there are several different parts here.

The last point is important, since to me some proposed checks for delayTime=0 do not make sense, e.g., having delayTime in an unused branch. However, if we fully allow delayTime=0 we get cases like:

model M
  Real x,y;
  parameter Real delayTime=0;
equation
  delay(x, delayTime)=y;
  der(y)=1-y;
end M;

Where the model suggest that we implicitly compute x from the zero-delay, which doesn't make sense to me. Obviously, if it were x=delay(y, delayTime); it could be made legal - but it will cause some other problems.

henrikt-ma commented 1 year ago

the tool tentatively evaluates it at compile time, and if it is zero it just removes the delay operator outright. If it is not, then the delay operator is retained, and changing delayTime to zero at runtime could trigger a runtime error, if the runtime cannot handle that.

The problem I see is that it could be considered a poor implementation if the runtime cannot handle it; it becomes a "quality of implementation" thing that pushes implementations to handle the most general and difficult situation in order to be considered of better quality. However, I don't think that our users are actually interested in a solution to this complex situation (reasons were given above). Second, I am not aware of any other construct in the specification that would rely on the trick of tentatively evaluating a parameter expression just to see if it has a value that would make the evaluation permanent, and it is not a kind of expected behavior I'd like to introduce.

I also don't think that a special rule for a constant delayTime of zero would be that useful in practice.

This is your personal opinion.

Let's call it a belief. Having a solution based on constant variability would make a lot more sense to me then something based around value-based parameter evaluation. I'm just afraid it wouldn't solve any interesting problems.

We are working with a library that needs this feature, I find it perfectly reasonable according to the currently released MLS 3.5, and I think it won't be good if it is made illegal by #3245 when MLS 3.6 is eventually released. There are probably other libraries in the same situation, since the need of describing dynamic behaviour with optional delay is something rather common in systems and control.

Yes, we all understand that we need models that can switch between having a delay and not having a delay. As has been illustrated in examples above, this could be achieved using explicit switching of the system structure instead of switching dynamically based on a parameter value.

Let me repeat that I don't think we're introducing any real backwards compatibility as long as there aren't any tools that can present a sensible strategy for dealing with the old specification. I suppose there is no tool that can correctly handle a time-varying that becomes zero, so which is a good example of why we should add at least restrictions without worrying about backward compatibility.

I feel sorry for the library developers, as always when there is code that has some sort of dependency on a flaw in the specification, but this won't be the last time libraries would need to be updated after things have been either fixed or "clarified" in the specification. Why couldn't explicit modeling of the structural change be a solution for the library that needs this feature?

casella commented 1 year ago

the tool tentatively evaluates it at compile time, and if it is zero it just removes the delay operator outright. If it is not, then the delay operator is retained, and changing delayTime to zero at runtime could trigger a runtime error, if the runtime cannot handle that.

The problem I see is that it could be considered a poor implementation if the runtime cannot handle it; it becomes a "quality of implementation" thing that pushes implementations to handle the most general and difficult situation in order to be considered of better quality.

@henrikt-ma, nobody will think bad of an implementation where you need to recompile the model when you change the parameter from zero from nonzero. The fact that you may have more or less good implementations is not a good argument to get rid of a language feature. Instead, people may be bothered if we declare their hitherto legal and working code illegal. Which is the net outcome of #3245 as it stands, in the case of parameter-varying delayTime.

However, I don't think that our users are actually interested in a solution to this complex situation (reasons were given above).

If your users are not interested, just give them the trivial implementation, since they won't use it anyway. Why do you want to bother other users?

Second, I am not aware of any other construct in the specification that would rely on the trick of tentatively evaluating a parameter expression just to see if it has a value that would make the evaluation permanent, and it is not a kind of expected behavior I'd like to introduce.

I don't see why this should be a problem. In any case, it's not mentioned in the normative part, we could as well skip this in the non-normative text. I understood you were worried that support for zero parameter-varying delayTime could not implemented, so I wanted to make it clear that it can be implemented, very easily.

We are working with a library that needs this feature, I find it perfectly reasonable according to the currently released MLS 3.5, and I think it won't be good if it is made illegal by #3245 when MLS 3.6 is eventually released. There are probably other libraries in the same situation, since the need of describing dynamic behaviour with optional delay is something rather common in systems and control.

Yes, we all understand that we need models that can switch between having a delay and not having a delay. As has been illustrated in examples above, this could be achieved using explicit switching of the system structure instead of switching dynamically based on a parameter value.

But you'll need to change the source code! It seems to me that you place zero value on backwards compatibility. I have a very different opinion: Modelica has been around for 25 years, there is lot of code out there, and we should try not to break it unless it is absolutely necessary.

Other modelling languages that I won't mention here follow a different path, continuously introducing non-backwards-compatible changes marketed as "improvements". With the result that you need to put extra work to just keep your old code alive. As a end-user, I find it very annoying and also a bit arrogant, and I stopped using those languages a long time ago also for this reason. I'd like Modelica to keep a different course.

I understand that the dream of every compiler and language designer is to design the best possible language without constraints, but Modelica is mature enough that we need to make compromises, because of the value of existing code. This is how I see it in general, also beyond this specific case.

Let me repeat that I don't think we're introducing any real backwards compatibility as long as there aren't any tools that can present a sensible strategy for dealing with the old specification. I suppose there is no tool that can correctly handle a time-varying that becomes zero, so which is a good example of why we should add at least restrictions without worrying about backward compatibility.

@henrikt-ma, as I wrote here, I'm totally fine with this restriction πŸ˜ƒ. BTW, there is no implementation that can handle time-varying delays becoming zero, so we are not breaking any backwards compatibility by introducing this restriction, because no library can be using a feature that is not supported by any tool, even if theoretically allowed by the specification.

My proposal is to limit this restriction to the case of time-varying delayTime, while allowing delayTime >= 0 for parameter-varying expressions. This case is allowed since Modelica 2.0, is used by existing libraries and is currently supported by at least two tools. I would find it slightly odd that now we decide to make it illegal because some implementors find it difficult to support it in a way that they personally find satisfying.

casella commented 1 year ago

The last point is important, since to me some proposed checks for delayTime=0 do not make sense

I'm sorry, what proposed checks are you referring to? Can you make concrete examples? We had some many different statements and code fragments in this endless ticket that I'm a bit lost.

However, if we fully allow delayTime=0 we get cases like:

model M
  Real x,y;
  parameter Real delayTime=0;
equation
  delay(x, delayTime)=y;
  der(y)=1-y;
end M;

Where the model suggest that we implicitly compute x from the zero-delay, which doesn't make sense to me.

I don't see any problem with that. When delayTime = 0 the model can be solved, otherwise it can't, so an error will be generated.

It's the same as writing x^2 + p = 0, this can only be solved reliably if p < 0, but it's not a good reason to forbid using x^2 in equations.

It's a fact of life that one can write nonsense models in Modelica, in which cases the models will be rejected by the tool at some stage.

Obviously, if it were x=delay(y, delayTime); it could be made legal - but it will cause some other problems.

Sorry, I don't get this. x=delay(y, delayTime); with delayTime = 0 has always been legal from MLS 2.0, when delay was introduced with delayTime >= 0, to the latest release MLS 3.5. As I understand, this works fine in Dymola and people have been using this feature for a long time. What problems are you referring to?

HansOlsson commented 1 year ago

The last point is important, since to me some proposed checks for delayTime=0 do not make sense

I'm sorry, what proposed checks are you referring to? Can you make concrete examples? We had some many different statements and code fragments in this endless ticket that I'm a bit lost.

The proposed check that caused even if delayTime>0 then delay(x, delayTime) else x to fail.

However, if we fully allow delayTime=0 we get cases like:

model M
  Real x,y;
  parameter Real delayTime=0;
equation
  delay(x, delayTime)=y;
  der(y)=1-y;
end M;

Where the model suggest that we implicitly compute x from the zero-delay, which doesn't make sense to me.

I don't see any problem with that. When delayTime = 0 the model can be solved, otherwise it can't, so an error will be generated.

Obviously, it will fail if delayTime > 0 - but should we allow it for delayTime = 0, and defer checking which case it is until simulation?

However, it seems clear that we should for now revert the change of forbidding delayTime=0, and have a larger discussion about it.

HansOlsson commented 1 year ago

Phone meeting: Just revert (for the moment will return later): Favor: Elena, Gerd, Ben, Markus, Hans Against: Abstain: Henrik, Martin, Quentin, Luigi

casella commented 1 year ago

Obviously, it will fail if delayTime > 0 - but should we allow it for delayTime = 0, and defer checking which case it is until simulation?

Why not?

However, it seems clear that we should for now revert the change of forbidding delayTime=0, and have a larger discussion about it.

I'm also in favour. I'm sorry but I could not join the phone meeting, I need to work asynchronously.