State variables with equality constraints

[MarkusOlssonModelon]

Is the following model allowed?

model Example
    constant Real eps = ...;
    Real x(start = 0, fixed = true);
    Real y(start = 0, fixed = true);
equation
    der(x) = time;
    der(y) = time;
    assert(noEvent(abs(x - y) < eps), "some message");
end Example;

Many solvers will estimate the derivatives in the jacobian by evaluating the equations with perturbations applied to the state variables. Dependent on the choices of eps and the perturbation, this can cause the assert to fail.

[henrikt-ma]

I don't think that the assert should be evaluated at such intermediate points of an integration method. However, while it seems easy enough to avoid for the model given here, I can imagine that you could come up with a variation of it where the assert is in the body of a function which needs to be evaluated in the vicinity of the true solution trajectory. In such a situation, I'd blame the model for lack of robustness.

[HansOlsson]

Many solvers will estimate the derivatives in the jacobian by evaluating the equations with perturbations applied to the state variables. Dependent on the choices of eps and the perturbation, this can cause the assert to fail.

It's always difficult to reason about models without underlying reason for them: in this case the reason for the assert. Is it for regression-checking - or is there some actual modeling reason?

[MarkusOlssonModelon]

It's always difficult to reason about models without underlying reason for them: in this case the reason for the assert. Is it for regression-checking - or is there some actual modeling reason?

The original model had an overconstrained connection graph with multiple roots. (I know that can be fixed by using Connections.potentialRoot)

[HansOlsson]

It's always difficult to reason about models without underlying reason for them: in this case the reason for the assert. Is it for regression-checking - or is there some actual modeling reason?

The original model had an overconstrained connection graph with multiple roots. (I know that can be fixed by using Connections.potentialRoot)

There seems to be some steps missing here. Do you mean that the graph was handled by introducing these asserts to ensure some condition for the non-selected roots (or for the equalityConstraints)?

[MarkusOlssonModelon]

It's always difficult to reason about models without underlying reason for them: in this case the reason for the assert. Is it for regression-checking - or is there some actual modeling reason?

The original model had an overconstrained connection graph with multiple roots. (I know that can be fixed by using Connections.potentialRoot)

There seems to be some steps missing here. Do you mean that the graph was handled by introducing these asserts to ensure some condition for the non-selected roots (or for the equalityConstraints)?

Here is a minimal example:

model Example
    type Angle
        extends Real;
        function equalityConstraint
            input Angle theta1;
            input Angle theta2;
            output Real residue[0];
        algorithm
            assert(abs(theta1 - theta2) < 1.e-10 , "Consistent angles");
        end equalityConstraint;
    end Angle;

    connector Plug
        Angle theta;
    end Plug;

    model Source
        Plug plug(theta(start = 0, fixed = true));
    equation
        Connections.root(plug.theta);
        der(plug.theta) = time;
    end Source;

    Source s1;
    Source s2;
equation
    connect(s1.plug, s2.plug);
end Example;

[HansOlsson]

Ah, overconstrained connection graphs in that domain - not multibody.

[HansOlsson]

The problem is that if they are different states with the same derivative we cannot guarantee that they proceed exactly in sync (even ignoring the Jacobian-computation). Note that this is used as an example in the specification.

An alternative would be to rewrite: der(p.theta) = 2*Modelica.Constants.pi*50 // 50 Hz; as: p.theta = 2*Modelica.Constants.pi*50*time+startPhase; // 50 Hz

[HansOlsson]

Design meeting: Two parts:

1. Use the alternative rewriting in example and push it to MSL. (And check if works in user-examples).
2. Should assertions be valid during Jacobian-evaluation or only on the trajectory?

Arguments for only evaluating assertions on the trajectory:

• Since example like this will otherwise break.
• Check that solution is correct with very tight tolerances - without being limited by finite difference approximations

Arguments for evaluating the assertions also for Jacobian and during intermediate steps in iterative solvers:

• Equations should be valid in a neighborhood of the trajectory
• Hard to define which assertions to skip (code after assertions in algorithms including calls of external functions)

[HansOlsson]

Phone meeting: Cannot easily disable for Jacobians, since for FMUs it is just model calls. (Unless providing directional derivatives.) Propose PR for spec and MSL. ( @HansOlsson ) Check if Buildings library, and git annotate. (Just remove from spec)

[eshmoylova]

Here is an example from the Buildings library: ACACTransformerFull. The equations with der on lines 87-88:

  theRef = PhaseSystem_p.thetaRef(terminal_p.theta);
  omega = der(theRef);

define der(terminal_p.theta[1]). theta is defined as a ReferenceAngle with equalityConstraint, essentially of the form given in PartialPhaseSystem after some inheritance with modification on constants.:

  type ReferenceAngle "Reference angle for connector"
    extends Modelica.Units.SI.Angle;

    function equalityConstraint "Assert that angles are equal"
      extends Modelica.Icons.Function;
      input ReferenceAngle theta1[:];
      input ReferenceAngle theta2[:];
      output Real residue[0];
    algorithm
      for i in 1:size(theta1, 1) loop
        assert(abs(theta1[i] - theta2[i]) < Modelica.Constants.eps,
          "Angles theta1 and theta2 are not equal over the connection.");
      end for;
    end equalityConstraint;
    annotation (Documentation(info="<html>
This type defines the voltage angle (used by the phasorial approach) for a specific connector that extends
<a href=\"modelica://Buildings.Electrical.PhaseSystems.PartialPhaseSystem\">
Buildings.Electrical.PhaseSystems.PartialPhaseSystem</a>.
</html>"));
  end ReferenceAngle;

Unlike in the example in the specification , however, it looks like der is used to define omega as a helper variable to compute some other quantities, instead of an equation to be solved for theta. But I am not familiar enough with the application to be certain about it.

[HansOlsson]

Looking more it seems that Buildings-library doesn't have the corresponding problem.

Consider Buildings.Electrical.AC.OnePhase.Sources.FixedVoltage where PhaseSystem.thetaRef(terminal.theta) = 2*Modelica.Constants.pi*f*time; so there's no state and thus no problem of equality constraints for state-variables. That also corresponds to the proposed idea for MSL.