Variability of expressions inside `initial` sections.

[qlambert-pro]

In section 3.8.3, the variability of expressions inside initial algorithm and initial equation are defined to be discrete-time. Shouldn't they rather be parametric? Since this reduces the level of checking, making the change should not create any backward compatibilities.

[gwr69]

One might also add when initial() as Section 3.8.3 includes

Expressions in the body of a when-clause

as being discrete-time.

This seems to be handled rather inconsistently in Modelica tools?

The following code should clearly be in violation of Section 3.8.3 as far as I understand the specs:

model ParameterModel
  Real x(start = 1.0, fixed = true);
  parameter Real k(fixed = false);
initial equation
  k = 2*x;
equation
  x = -k*der(x);  
end ParameterModel;

[henrikt-ma]

Since this reduces the level of checking, making the change should not create any backward compatibilities.

The problems wouldn't come in the from of backwards incompatibility, but in the form of opening up for unhandled cases of equations taking advantage of the more permissive variability rules. Still, I agree that the specification should be better at explaining why the expressions are discrete-time rather than parameter. If there is no good reason for it, I suppose the sensible thing to do would be to actually change them to be parameter.

[henrikt-ma]

The following code should clearly be in violation of Section 3.8.3 as far as I understand the specs:

Are you referring to the master branch of the specification, or the 3.4 release? In 3.4, I don't think there is any variability requirements on your initial equation. On master, on the other hand, I think we are missing this case of an initial equation, as your equation wouldn't be considered contributing to the solution of k. This, however, looks like exactly the kind of problem that this issue is all about — it would be resolved elegantly if the variability of x was parameter in the initial equation.

[henrikt-ma]

One might also add when initial() as Section 3.8.3 includes

Expressions in the body of a when-clause

as being discrete-time.

I think this one isn't that clear, since initial() isn't in general the only condition triggering a when-clause, but reducing everything to parameter would only make sense for the clause when used during initialization. To me, it would seem better to keep the discrete-time variability for all when-clauses instead of making an exception for the ones triggered by initial() alone.

[gwr69]

I believe I looked atthe lonked master version given as html.

The problem I am facing is solving for some initial value that needs to be passed as an initialValue parameter for a component (where it sets the start attribute). How to do that with discrete-time variability then?

[henrikt-ma]

The problem I am facing is solving for some initial value that needs to be passed as an initialValue parameter for a component (where it sets the start attribute). How to do that with discrete-time variability then?

I'm afraid there is no good answer to that, but I invide others to prove me wrong. I think you present a good use case that should increase the interest in finding a good resolution to this issue.

[HansOlsson]

Let's get back to the issue.

In section 3.8.3, the variability of expressions inside initial algorithm and initial equation are defined to be discrete-time. Shouldn't they rather be parametric? Since this reduces the level of checking, making the change should not create any backward compatibilities.

I'm unsure what change it would actually cause, can you explain that in more detail? Is this issue based on some specific problem?

I can understand the historic reasons for making them discrete-time - I believe "initial equation" started out as "when initial()". Thinking more I can see one additional reason: we generally don't have systems of equations for parameters: Thus parameter-bindings must acyclical ( https://specification.modelica.org/master/class-predefined-types-and-declarations.html#acyclic-bindings-of-constants-and-parameters ), whereas initial equations can include cyclical equations. I know it's not a strong argument.

On the other hand I can see some weirdness related to parameters with fixed=false, but I'm not sure this change would be the best solution.

[qlambert-pro]

The specific problem is the one that Guido presented. The change would be to give any expression inside of the body of an initial section variability parametric. It wouldn't pause a problem as far as variability checking goes but as Henrik pointed out, if it did create problems they would be created at initialization.

[gwr69]

@HansOlsson Is the nature of solving for initial values parametric or discrete-time? I am probably in over my head, but I see at least some significant subset of such problems that I would consider rather parametric. How to declare that?

[HansOlsson]

@HansOlsson Is the nature of solving for initial values parametric or discrete-time? I am probably in over my head, but I see at least some significant subset of such problems that I would consider rather parametric. How to declare that?

Do you mean something like:

  parameter Real p(fixed=false);
  Real x(start=p);                         
initial equation
  p+p^3=9;

I can see a problem here, but to me this seems different than variability and thus I'm not sure that just changing variability will solve that. I recall that we do something extra for cases related to this - but I will have to investigate more.

[sjoelund]

OpenModelica does not do anything special for that case from what I can see. It's just an equation system containing a parameter as unknown.

[qlambert-pro]

I think there seems to be a good rule of thumb between when an expression can be evaluated and its variability, constant expressions are evaluated around flattening, parameter expressions are evaluated at initialization, discrete-time expressions are evaluated at event and the rest are continuous-time. As a result, and although I understand the historic argument, it feels weird that, initial section expressions would be discrete-time rather than parametric.

[henrikt-ma]

OpenModelica does not do anything special for that case from what I can see. It's just an equation system containing a parameter as unknown.

I agree. There is nothing strange about solving for p in that initial equation, as it has the maximum variability in its equation. Then, using it for x.start doesn't seem problematic in any way.

I can see a problem here, but to me this seems different than variability and thus I'm not sure that just changing variability will solve that. I recall that we do something extra for cases related to this - but I will have to investigate more.

Right, we also treat the initialization of x a bit different in this case compared to more simple cases. However, I'd consider the more simple cases just optimized variants of the general case seen here.

[henrikt-ma]

Thinking more I can see one additional reason: we generally don't have systems of equations for parameters: Thus parameter-bindings must acyclical

As I see it, this is a special feature that comes with using a declaration equation for a parameter. For parameters without a declaration equation I have no expectation of them not ending up in a system of equations.

[HansOlsson]

Well, if you want something more complicated consider:

  parameter Real p3=4;
  parameter Real p2(fixed=false, start=p3);
  parameter Real p(fixed=false, start=p2);
  Real x(start=p, fixed=true);
initial equation 
  p2+p2^3=7+x;
  p+p^3=9+x;
equation 
  der(x)=-x;

In this case Dymola first solve for p, x so that x.start=p ignoring the start-value for p, and then for p2 (using x).

Note: Added fixed=true for x. And the logic is that if a guess-value cannot be computed we ignore it, instead of having a system of equations just for a guess-value.

[qlambert-pro]

What do people think of the following example:

model ParameterModel
  Real x(start = 1.0, fixed = true);
  parameter Real k(fixed = false);
initial algorithm
  k = 2*x;
equation
  x = -k*der(x);  
end ParameterModel;

[gwr69]

Should probably be:

model ParameterModel
  Real x(start = 1.0, fixed = true);
  parameter Real k(fixed = false);
initial algorithm
  k := 2*x;
equation 
  x = -k*der(x);  
end ParameterModel;

and since x has fixed = true and start = 1.0 wouldn't this be really be equivalent to

model ParameterModel
  Real x(start = 1.0, fixed = true);
  parameter Real k(fixed = false);
initial equation
  k = 2*x;
equation 
  x = -k*der(x);  
end ParameterModel;

and also to

model ParameterModel
  Real x(start = 1.0, fixed = true);
  parameter Real k(fixed = false);
equation
  when initial() then
     k = 2*x;
  end when; 
  x = -k*der(x);  
end ParameterModel;

[HansOlsson]

Should probably be:

model ParameterModel
  Real x(start = 1.0, fixed = true);
  parameter Real k(fixed = false);
initial algorithm
  k := 2*x;
equation 
  x = -k*der(x);  
end ParameterModel;

As far as I can tell that is the previous model

and since x has fixed = true and start = 1.0 wouldn't this be really be equivalent to

model ParameterModel
 Real x(start = 1.0, fixed = true);
 parameter Real k(fixed = false);
initial equation
 k = 2*x;
equation 
 x = -k*der(x);  
end ParameterModel;

Yes, that would be the same.

and also to

model ParameterModel
 Real x(start = 1.0, fixed = true);
 parameter Real k(fixed = false);
equation
 when initial() then
    k = 2*x;
 end when; 
 x = -k*der(x);  
end ParameterModel;

No, that model shouldn't work as there is nothing special with "when initial()" and you need to solve "k" from "initial equation"s, and there are no initial equations here.

[HansOlsson]

However, what is still not clear to me is the following: it seems it would allow more models. But do we have any actual models that are currently invalid that would be valid if we change the semantics in this way?

[qlambert-pro]

The intitial algorithm example is currently invalid, since the left-hand side of the statement has a lower variability than the right-hand side.

[gwr69]

I was simply pointing at the fact, that within an initial algorithm section we would write statements not equations so it should be k := 2*x and not k = 2*x.

UPDATE: Ah, indeed, this is not a use case for when initial()—it has to be either initial algorithm or initial equation.

[HansOlsson]

The intitial algorithm example is currently invalid, since the left-hand side of the statement has a lower variability than the right-hand side.

Ok, that is something to consider. But in the future - please post such statements in the post itself, so that we don't have to guess what the issue is.

Regarding that issue:

https://specification.modelica.org/master/equations.html#initialization-initial-equation-and-initial-algorithm States that fixed=false parameters are unknown during the initialization (and thus can be solved at all), so we need to either clarify that they can also be left-hand-side variables in "initial algorithm" (which makes sense as they are unknown during the initialization); or make expressions in "initial equation/algorithm" parameter expressions and we might still need the same rule saying that only unknown variables during initialization may appear as lhs-variables in "initial algorithm"s.

I can see one additional issue with making expressions in "initial equation/algorithm" parameter expressions: the variability of pre becomes very messy, as pre(p) is a discrete-time expression even for a parameter p.

[henrikt-ma]

or make expressions in "initial equation/algorithm" parameter expressions and we might still need the same rule saying that only unknown variables during initialization may appear as lhs-variables in "initial algorithm"s.

I would have hoped that such a "lhs-variable rule" wouldn't be needed, since the "known" variables should be known due to being solved by some other equation. Hence, general handling of overdetermined equations should be sufficient.

[henrikt-ma]

I can see one additional issue with making expressions in "initial equation/algorithm" parameter expressions: the variability of pre becomes very messy, as pre(p) is a discrete-time expression even for a parameter p.

But isn't this simply because we have specified pre in the wrong way? Since pre requires the operand to be discrete-time, can't pre just follow normal variability rules and preserve the variability of the operand?

[HansOlsson]

or make expressions in "initial equation/algorithm" parameter expressions and we might still need the same rule saying that only unknown variables during initialization may appear as lhs-variables in "initial algorithm"s.

I would have hoped that such a "lhs-variable rule" wouldn't be needed, since the "known" variables should be known due to being solved by some other equation. Hence, general handling of overdetermined equations should be sufficient.

If we don't have a rule about lhs-variables we allow

  parameter Real p=2;
  parameter Real q(fixed=false);
initial algorithm 
 p:=q+1; // It's an algorithm that seems to assign to p, but we then solve it for q

I noticed that Dymola currently accepts this, but it really looks weird to me.

[HansOlsson]

But isn't this simply because we have specified pre in the wrong way? Since pre requires the operand to be discrete-time, can't pre just follow normal variability rules and preserve the variability of the operand?

Possibly. In that case pre(p) will be the same as p for parameters (in addition to becoming a parameter expression).

Sounds good to me.

[henrikt-ma]

Well, if you want something more complicated consider:

  parameter Real p3=4;
  parameter Real p2(fixed=false, start=p3);
  parameter Real p(fixed=false, start=p2);
  Real x(start=p, fixed=true);
initial equation 
  p2+p2^3=7+x;
  p+p^3=9+x;
equation 
  der(x)=-x;

In this case Dymola first solve for p, x so that x.start=p ignoring the start-value for p, and then for p2 (using x).

Note: Added fixed=true for x. And the logic is that if a guess-value cannot be computed we ignore it, instead of having a system of equations just for a guess-value.

I actually find this off-topic, but couldn't resist trying it out anyway. SystemModeler is actually rejecting this currently due to an internal error check that might be too harsh. With a small modification, I can get some sort of initialization, but I don't like what I see:

• p and x are solved first in a system with x as iteration variable. Hence, the start = p2 of p ends up being unused, but it it were used, it would have been some sort of uninitialized read.
• We fail to report that x isn't initialized to its start before solving the nonlinear system of equations.
• Finally we solve for p2. By changing the equation to a quadratic polynomial, p2 + p2^2 = 7 + x, I can verify that changing p4 actually affects which of the two solutions we get — nothing strange at this point.

I'd like to blame all those failures on a problem that isn't well posed, but unfortunately I don't find support in the specification for rejecting this model as invalid. So yes, this is definitely something more complicated than the other examples in this issue, but I don't think we should try to work out the example here and now. Instead, I invide those who are interested about complicated initialization problems to join the ongoing Flat Modelica discussions on the matter, where a common understanding of what full Modelica initialization can and cannot do is necessary in order to work out a more elementary initialization model for Flat Modelica:

• https://github.com/modelica/ModelicaSpecification/pull/2604 (concrete proposal under development)
• https://github.com/modelica/ModelicaSpecification/pull/2561 (wider discussion on initialization and other topics)

[henrikt-ma]

If we don't have a rule about lhs-variables we allow

  parameter Real p=2;
  parameter Real q(fixed=false);
initial algorithm 
 p:=q+1; // It's an algorithm that seems to assign to p, but we then solve it for q

I noticed that Dymola currently accepts this, but it really looks weird to me.

I think Dymola does the right thing. This is the kind of surprises that come with the power of Modelica.

[HansOlsson]

Trying to get back on track and summarizing the proposed changes, and thinking more:

• Change expressions in "initial equation" and "initial algorithm" to be parametric.
• Remove exception for: pre, edge, change, der, initial, terminal, (but not sample, right?) that makes them non-parametric even with parameter-input; since otherwise some initial equations will still non be parametric.
• Or: keep exception but say that "initial equation" and "initial algorithm" parametric trumps that; i.e. so that pre is still a discrete-expression regardless of sub-expressions, but any expression in initial equations is a parameter expression.

The first variant (change of operators) seems like a clean-up, but apart from that I don't see any obvious benefits - and if we keep the exceptions we don't even get a clean-up - and to me it just becomes messier.

The change of operator will make previously illegal models legal Consider:

  parameter Real p=2;
  parameter Real q=pre(p); // Currently illegal, since pre is not parametric, but would be legal with change

I just view that as a curiosity. But more importantly it will break the semantics, as if initial() then ... as used in many models will not make sense if initial is a parameter-expression (as it actually changes during the simulation). And even if odd you can currently use initial() in initial equations; so if we keep it as a discrete expression but change "initial equation" and "initial algorithm" we will break those models.

I might have missed something, especially some benefit of making expressions in "initial equation" parameter expressions if we keep the exceptions for the operators, if so please enlighten me. Otherwise I think we should just keep it as is, and focus on more impactful changes.

The major reason I can see for having them as discrete expressions is that we use "discrete expressions" as a sort of generic catch-all when we don't want anything odd (no special crossing function treatment); and thus also use it that for functions. I admit that it is weak, but at least it's something.

I do see some problems with variability and initial equations, in particular solving for start-value (as above) and that when solving them we mix them with normal equations - which is a bit odd; but I don't see how this proposal solves it.

[henrikt-ma]

I think it would be hard to express the variability rules without exceptions in one form or another for built-in operators like der, initial, and terminal. The specification might get some added clarity by introducing the concept variability-preserving operator for all the other built-in operators (by listing all the operators that are not variability-preserving directly under 3.8 Variability of Expressions), so that the main variability rules could be stated in an exception-free manner, for example:

A function call or a variability-preserving operator with parameter subexpressions is a parameter expression.

The current exceptions for pre, edge and change appear to be a different story; they are variability-preserving as far as I can tell.

Regarding what trumps what, I think it is important for the clarity of the variability rules that one should always look for the first (lowest variability) matching rule, and never need to look further to find the variability of an expression. Hence, if the rule for initial equation and initial algorithm was placed under parameter expressions, this should be enough.

[henrikt-ma]

But more importantly it will break the semantics, as if initial() then ... as used in many models will not make sense if initial is a parameter-expression (as it actually changes during the simulation). And even if odd you can currently use initial() in initial equations; so if we keep it as a discrete expression but change "initial equation" and "initial algorithm" we will break those models.

Yes, this really sounds problematic. I think it would be enlightening to take a closer look at such an example, do you have one to share?

[henrikt-ma]

I might have missed something, especially some benefit of making expressions in "initial equation" parameter expressions if we keep the exceptions for the operators, if so please enlighten me. Otherwise I think we should just keep it as is, and focus on more impactful changes.

This is an impactful change as it allows non-fixed parameters to be solved from equations also involving (the initial value of) time-varying variables. Two concrete variants of this need were given above:

• Compute a parameter, say initialValue, to be used as start attribute of another variable. Here, initialValue depends on the initialization of some other time-varying variable, say z. Even if initialValue is just a simple expression in z (for example, initialValue := z), this expression cannot be used directly as a start modifier, as start has parameter variability. Currently, assigning z to initialValue in an initial algorithm isn't possible due to variability violation.
• Currently, the initial equation z = p; isn't considered contributing to the solvability of p in the perfect matching rule.

I'd also like to add that the current situation is problematic because it doesn't agree with what I believe is common intuition about the variability of things in the initialization problem.

[HansOlsson]

I might have missed something, especially some benefit of making expressions in "initial equation" parameter expressions if we keep the exceptions for the operators, if so please enlighten me. Otherwise I think we should just keep it as is, and focus on more impactful changes.

This is an impactful change as it allows non-fixed parameters to be solved from equations also involving (the initial value of) time-varying variables. Two concrete variants of this need were given above:

• Compute a parameter, say initialValue, to be used as start attribute of another variable. Here, initialValue depends on the initialization of some other time-varying variable, say z. Even if initialValue is just a simple expression in z (for example, initialValue := z), this expression cannot be used directly as a start modifier, as start has parameter variability. Currently, assigning z to initialValue in an initial algorithm isn't possible due to variability violation.

This seems like the following:

  parameter Real initialValue(fixed=false);
  parameter Real init2=2*initialValue; // Or was this an issue?
  Real z;
  Real x(start=initialValue+init2);
initial algorithm 
  initialValue:=z; // Is this the potential problem?
initial equation 
  z+der(z)=1;
equation 
  der(z)=-2*z;
  der(x)=-x;

Dymola currently allows it, and I agree that this should legal - if that is unclear we have to make some change. There is no corresponding issue for the corresponding "initial equation".

• Currently, the initial equation z = p; isn't considered contributing to the solvability of p in the perfect matching rule.

I don't see the issue here, obviously we need a separate perfect matching rule for the initialization problem.

I'd also like to add that the current situation is problematic because it doesn't agree with what I believe is common intuition about the variability of things in the initialization problem.

[henrikt-ma]

This seems like the following:

initial algorithm 
  initialValue:=z; // Is this the potential problem?

Dymola currently allows it, and I agree that this should legal - if that is unclear we have to make some change.

Yes, this violates the variability rule for assignment statements, and would be solved elegantly by also giving 'z' parameter variability in the 'initial algorithm'.

[henrikt-ma]

There is no corresponding issue for the corresponding "initial equation".

• Currently, the initial equation z = p; isn't considered contributing to the solvability of p in the perfect matching rule.

I don't see the issue here, obviously we need a separate perfect matching rule for the initialization problem.

I suppose you mean that the perfect matching rule needs to be applied once for the initialization problem, and once for the integration problem, but that the rule is the same in both cases? If so, then the problem is that the perfect matching rule (when applied to the initialization problem) only considers z solvable in z = p, as the solution for p would have variability that isn't compatible with p. Again, this would be solved elegantly by giving z parameter variability in the initial equation.

[HansOlsson]

This seems like the following:

initial algorithm 
  initialValue:=z; // Is this the potential problem?

Dymola currently allows it, and I agree that this should legal - if that is unclear we have to make some change.

Yes, this violates the variability rule for assignment statements, and would be solved elegantly by also giving 'z' parameter variability in the 'initial algorithm'.

Or, one can say that it is already legal, as 'initialValue' is in an initial algorithm, and thus has 'discrete-time variability'.

[HansOlsson]

I suppose you mean that the perfect matching rule needs to be applied once for the initialization problem, and once for the integration problem, but that the rule is the same in both cases?

Well, same idea of perfect matching, but the initial problem has an extended set of variables and equations to match; which should be clarified.

[henrikt-ma]

This seems like the following:

initial algorithm 
  initialValue:=z; // Is this the potential problem?

Dymola currently allows it, and I agree that this should legal - if that is unclear we have to make some change.

Yes, this violates the variability rule for assignment statements, and would be solved elegantly by also giving 'z' parameter variability in the 'initial algorithm'.

Or, one can say that it is already legal, as 'initialValue' is in an initial algorithm, and thus has 'discrete-time variability'.

This would require a non-standard reading of section 3.8, depending on some sort of unclear conflict resolution between 3.8.2 (the variable is declared parameter) and 3.8.3 (appearing inside initial algorithm). I believe the standard reading is that the variability is the one in the first section with a matching rule, which is 3.8.2 (parameter) in this case.

[henrikt-ma]

I suppose you mean that the perfect matching rule needs to be applied once for the initialization problem, and once for the integration problem, but that the rule is the same in both cases?

Well, same idea of perfect matching, but the initial problem has an extended set of variables and equations to match; which should be clarified.

Sure, it certainly wouldn't hurt.

[henrikt-ma]

Another thing to keep in mind here is how one should or shouldn't be able to use external objects in the initialization of parameters. For example, in MSL 3.2.2, the CombiTimeTable contained the following initial algorithm, where tableID is an external object, and tableOnFileRead, t_minScaled and t_maxScaled parameters:

initial algorithm
  if tableOnFile then
    tableOnFileRead := readTableData(tableID, false, verboseRead);
  else
    tableOnFileRead := 1.;
  end if;
  t_minScaled := getTableTimeTmin(tableID, tableOnFileRead);
  t_maxScaled := getTableTimeTmax(tableID, tableOnFileRead);

Note that the tableID argument to the three external function calls make these calls have discrete-time variability (since the continuous-time variability of tableID is capped to discrete-time inside the initial algorithm).

Was it because of variability violation that this was changed in MSL 3.2.3, so that this can be taken as evidence that the variability error here is really the intention of the specification? If so, we need to be careful when resolving this issue, so that we don't make it legal. One way of doing this would be to say that inside initial sections, the variability of all component references except external objects is at most parameter.

This also points at a slightly different way of looking at the bigger picture in this issue: Instead of fiddling with the variability of any expression inside initial sections, it might be better to target just the variables appearing in expressions, and derive other expression variabilities in the normal way. If nothing else, this paves the way for the more restrictive use of external objects in initial sections.

[HansOlsson]

Another thing to keep in mind here is how one should or shouldn't be able to use external objects in the initialization of parameters. For example, in MSL 3.2.2, the CombiTimeTable contained the following initial algorithm, where tableID is an external object, and tableOnFileRead, t_minScaled and t_maxScaled parameters:

That change wasn't due to a problem with variability, but because the algorithm had a side-effect (for tableOnFileRead, and removing that means that t_m**Scaled can be initialized without an initial algorithm) which caused problems in terms of being fragile and not handling "continue": https://github.com/modelica/ModelicaStandardLibrary/issues/1941 https://github.com/modelica/ModelicaStandardLibrary/issues/1899

[HansOlsson]

Two possibilities:

• Raise variability in initial algorithm (and initial equation), so ignore variability-check in initial algorithm (and initial equation). Somewhat odd interpretation of current text.
• Lower variability in initial algorithm (and initial equation). Change of current text. Are they equivalent?

Quentin test-implement. MSL and MSL-test-suite.

[qlambert-pro]

Can you specify a little bit more what you mean by "ignore variability-check in initial algorithm" In particular, how would you account for the following model:

model NonStructuralForRange
  Real[10] x;
  Integer n;
initial equation
  for k in 1 : n loop
    x[k] = k;
  end for;
equation
  der(x) = -x;
end NonStructuralForRange;

[qlambert-pro]

My implementation of the second strategy, introduced a variability level between discrete-time and parametric. One could picture it as strictly parametric variability as opposed to potentially structurally parametric variability (aka parametric variability in the specs). I considered expressions inside initial sections to have strictly parametric variability.

Lowering the variability in initial algorithm didn't lead to any failures in MSL-3.2.3 and MSL-4.0.0, neither did deactivating variability checking inside initial sections.

The second strategy allows the rejection of the previous model on variability basis, the first strategy does not.

Let me know if I should test something else.

[henrikt-ma]

My implementation of the second strategy, introduced a variability level between discrete-time and parametric. One could picture it as strictly parametric variability as opposed to potentially structurally parametric variability (aka parametric variability in the specs). I considered expressions inside initial sections to have strictly parametric variability.

Formalization of this variability level is proposed in #2754, so we should probably try to finish that one before returning to this issue.