Extend overconstrained connector semantics to handle dynamic branches

[casella]

This is a second, improved attempt to address this issue, following #6238

The overconstrained connector semantics defined in Section ​9.4 of Modelica 3.4 considers a fixed connection graph, that can be handled statically at compile time.

One of the cited examples for overconstrained types in connectors is the phase information in AC power systems. The current semantics is able to represent synchronous systems whose topology remains unchanged during simulation. The phase reference is generated by one component of the synchronous system (normally a power generator) and distributed throughout the entire synchronous system via the overconstrained connector variable. In fact, it is possible to have multiple synchronous systems in the same Modelica model, that correspond to structurally disconnected connections sub-graphs, e.g. two national grids connected by a DC link, but in any case their topology is fixed at runtime.

When modelling AC trasmission systems, is is possible that in case of severe perturbations, some key circuit breakers are switched open (i.e., their conductance brough to zero), effectively splitting a single synchronous system into multiple synchronous islands. Note that this requires no structural changes in the grid equations, only to set some admittance values to zero.

When this happens, the two (or more) ensuing islands can settle into new steady states with different frequencies. Hence, if a single, whole-system-wide reference is still used, the phase angle of currents and voltages of the now disconnected island(s) will continue to rotate permanently with a frequency that is the difference between the local island frequency and the frequency of the island where the original root node was picked. This is very inconvenient, because it prevents variable step-size solvers to increase the step size, once the system settles to the new steady state.

This problem could be avoided by allowing to dynamically break and reconnect the breakable branches corresponding to connect() statements in the connection graph of ​Section 9.4. It would then be possible to break the synchronous connections established by transmission lines when their admittance is brough to zero, thus modelling effects of circuit breakers on the synchronous sub-system topology.

As a consequence, two or more disjoint connection graphs would be formed at the time of the breaker openings, each corresponding to a new synchronous island. The new graph topology should be analyzed at this point, picking a new root node for each newly formed island in the grid. Then, instead of having a single phase reference for the entire system, which would no longer be adequate, one would now have two or more independent phase references, one for each island, which would ensure that the variables of each island reach a steady state, thus avoiding the persistent sinusoidal oscillations found in the case of a statically determined connection topology.

From the point of view of the language specification, this extension only requires to add one optional Boolean input to the connect() statement

  connect(A.R, B.R, activeBranch = true);

where activeBranch tells whether a breakable branch should be added to the connection graph (possibly to be broken at a later stage to avoid loops), or not.

If activeBranch is set to a constant value of true (possibly by default), the behaviour is exactly the same as in Modelica 3.4. In case instead it is set to a Boolean variable, the breakable branch will only be added to the graph when activeBranch = true.

From an implementation point of view, this means that the graph connection analysis no longer can be performed statically at compile time, but needs to be deferred to run time. According to ​Section 9.4.2, everything except overconstrained connector variables is always handled in the same way, i.e., following ​Section 9.2, so the runtime analysis will be limited to overconstrained variables in connectors. Everthing else, including determining connection sets, building connection equations for flow variables, and handling stream connector variables, is unchanged and can still be managed statically at compile time.

In general this can be quite a big deal, because Modelica cannot handle structural dynamics at runtime. However, the special case of scalar overconstrained connector variables with zero-dimensional equalityConstraint() function output turns out to be very easy to handle, and is enough to address the issue with the synchronous AC grids. This is demonstrated by the attached examplary case.

The attached package contains components to build conceptual models of phasor-based AC power grid models. The components are overly simplified and only retain the minimal features that are needed to analyze this language semantics extension.

In this case, the overconstrained variable is actually not the phase angle but rather the angular velocity of the phasor reference frame, which is sufficient to ensure that all state variables are constant when the synchronous islands are all operating in steady state.

System 1 is a basic example with two generators G1 and G2, each connected to a local load, connected by an inductive transmission line T. The system starts in steady state with balanced loads, then at time = 1 the active power consumption of load L2 is reduced by 20%. This starts a transient, that eventually dies out once the primary frequency controls stabilize the system at a new, higher frequency. Note that at that point, all state variables are constant, allowing implicit solvers to take large time steps.

System 2 is similar to System1, except that line T is split in the series connection of lines T1a and T1b, running in parallel, with line T2. This allows to demonstrate the overconstrained connector mechanism, since T1a and T1b form a loop that must be broken to avoid getting an overconstrained system of equations. The overall impedance is the same, so the power and frequency transients are the same as in System1.

In System 3, line T2 also includes a circuit breaker mechanism, that brings its susceptance B to zero when the T2.open signal becomes true. From that point onwards, there is no longer any exchange of power along the lines, so the system is effectively split into two synchronous islands, one containing G1 and L1, the other one containing G2 and L2. The frequencies of the two drift apart, causing the phasors and state variables of G2 and L2, that still use G1.port.omegaRef as a reference angular velocity, to oscillate forever even when the new steady state is reached.

In order to avoid this, the connection graph needs to be split when the breaker embedded in T2 is opened. This behaviour is achieved in System4 by using the modified TransmissionLineVariableBranch component, that contains two public and one protected ports

  ACPort port_a;
  ACPort port_b;
protected
  ACPort port_b_int;

and the following equations

  port_a.i = Complex(0,B_act)*(port_a.v - port_b_int.v);
  when open then
    closed = false;
    B_act = 0;
  elsewhen close then
    closed = true;
    B_act = B;
  end when;
  port_a.omegaRef = port_b_int.omegaRef;
  Connections.branch(port_a.omegaRef, port_b_int.omegaRef);
  connect(port_b_int, port_b, closed);

where the breakable branch between the protected port port_b_int and the public port_b is only activated when closed = true.

Handling the case of Real overdetermined connector variables with empty constraint function is relatively straighforward.

During flattening, the front end should indentify all the structurally connected subgraphs (i.e., ignoring the values of activeBranch), then identify those of them that only contain static breakable branches (for which activeBranch = true by default or by means of constant-evaluated binding) and those that contain at least one dynamically breakable branch.

The former should be handled as in Modelica 3.4. For the latter, all the relevant information about overconstrained types and variables, and about nodes and branches of the corresponding subgraphs, should be passed to the backend using some suitable formalism. Subgraphs that are structurally unconnected, regardless of activeBranch status of their breakable branches should be identified.

I guess we should try to re-use the front-end code in the backend as much as possible, and even in the runtime code (e.g. for graph analysis) as long as the MetaModelica code of the frontend can be translated into C/C++ code.

The generated code should then perform the following operations at initialization and each time any activeBranch connection attribute changes:

• re-build the connection graphs of all the structurally unconnected sub-graphs to which the branches that changed their activation status belong, according to the rules set forth in ​Section 9.4, only accounting for branches with activeBranch = true
• in case a previously de-activated branch is re-established, check that the corresponding connection equation involving the overconstrained connector variables is actually verified within some tolerance, otherwise aborting the simulation; this would correspond to the requirement that two islands can only be reconnected if the are operating exactly at the same frequency or phase. This check could be defined by adding one more function reconnectable() to the overconstrained type definition, which returns a boolean indicating if the reconnection is feasible (see below).
• identify disjoint sub-graphs within each structurally connected, re-built sub-graph, and select one root node for each of them, according the standard rules set forth in ​Section 9.4.2
• break all loops in those sub-graphs
• set each overconstrained connector variable of each sub-graph to be equal to the value of the corresponding root node variable (this can be done trivially via pointers!)

If a frequency reference is used, the reconnectable function could be implemented as

function reconnectable
  input ReferenceFrequency f1;
  input ReferenceFrequency f2;
  output Boolean canReconnect;
protected
  constant Real eps = 1e-6;
algorithm
  canReconnect := abs(f1 - f2) < eps;
end reconnectable;

For phase references, one should account for 2*pi ambiguity

function reconnectable
  input ReferenceAngle theta1;
  input ReferenceAngle theta2;
  output Boolean canReconnect;
protected
  constant Real eps = 1e-6;
  SI.Angle delta;
algorithm
  delta := (theta1 - theta2)- div(theta1 - thetha2, 2*pi)*2*pi;
  canReconnect := abs(delta) < eps;
end reconnectable;

I suppose this feature could be implemented rather easily in OpenModelica, as a first prototype implementation of this language extension. It could be tested in an extended version of the attached package, including some more examples with a bit more nodes and edges in the connection graphs, and finally tested on an extended version of the ​PowerGrids library, that already envisioned such a feature.

If the experiments with the prototype are successful, this could form the basis of an MCP to introduce this feature into the next version of the Modelica Specification.

[casella]

I made some reflections, and I think I made a mistake in the original design described in this ticket. There, I assumed that the break-up of islanded sub-systems should take place by optionally forcing the break-up of breakable branches, i.e., connection equations. I now understand this was not a good idea.

The problem with that idea is that breaking connections not only influences the set of equations related to the overconstrained connector variables (which is OK), but it also changes the connection sets, and hence the equations for flow and effort variables. This introduces additional structural variability in the system, which would be hard to manage.

The good news is, this is really not necessary. The break-up of synchronous sub-systems doesn't really need to undo connections. It just needs the admittance of the branch to become zero, so that there is no longer any power flow through it, and the two sub-systems at each side of the branch can rotate at different speeds unconstrained. This does not require to change the structure of the equations, it is just a consequence of the current through the line becoming zero because the admittance is now zero, so that there is no longer any power flow depending on the relative angle, that keeps the two parts of the grid in sync with each other.

So, the right way to handle dynamic connection sets is to inform the compiler that when the admittance is brought to zero, we are actually breaking the branch corresponding to the transmission line, i.e. the unbreakable branch defined by the Connections.branch statement. In this way, the runtime can update the connected spanning trees and come up with better reference values for the overconstrained variables in the connectors.

I was indeed misled by the concept of breakable/unbreakable branches, so I didn't want to break unbreakable branches. In fact, "breakable" in this context means "breakable by a topology analysis algorithm that avoids redundant equations by breaking loops into spanning trees". If we open a switch by setting the admittance to zero, we know the branch has to be broken there, we don't need an algorithm to determine that.

Of course this means we also have to explicitly remove the equation that relates the overconstrained variables on the two connectors of the line, in case it gets opened, otherwise we may end up with redundant equations. This requires an exception to the rule set forth in Section 8.3.4, because we'd need to have a conditional equation which is only activated when the branch is active. This exception is justified because we are not running the risk to not have the right number of equations, since this is handled nicely by the graph topology algorithm.

In the next comment I'll lay out the revised proposal.

[casella]

The overconstrained connector semantics defined in Section ​9.4 of Modelica 3.5 considers a fixed connection graph, that can be handled statically at compile time.

One of the cited examples for overconstrained types in connectors is the phase information in AC power systems. The current semantics is able to represent synchronous systems whose topology remains unchanged during simulation. The phase reference is generated by one component of the synchronous system (either an infinite bus, or a large synchronous generator for islanded systems) and distributed throughout the entire synchronous system via the overconstrained connector variables. In fact, it is possible to have multiple synchronous systems in the same Modelica model, that correspond to structurally disconnected connections sub-graphs, e.g. two national grids connected by an undersea DC link, but in any case their topology is fixed at runtime.

When modelling AC trasmission systems, is is possible that in case of severe perturbations, some key circuit breakers are switched open (i.e., their admittance brought to zero), effectively splitting a single synchronous system into multiple synchronous islands. Note that this requires no structural changes in the grid equations, only to set some admittance values to zero.

When this happens, the two (or more) ensuing islands can settle into new steady states with different frequencies. Hence, if a single, whole-system-wide reference is still used, the phase angle of currents and voltages of the now disconnected island(s) will continue to rotate permanently with a frequency that is the difference between the local island frequency and the frequency of the island where the original root node was picked. This is very inconvenient, because it prevents variable step-size solvers to increase the step size, once the system settles to the new steady state.

This problem could be avoided by allowing to dynamically add or remove the unbreakable branches corresponding to Connections.branch() statements in the connection graph of ​Section 9.4, based on the status of the corresponding breaker components, and to add or remove the equations that show up in the same if-equation branches where the Connections.branch() statements are declared. It would then be possible to break up the original synchronous connections established by transmission lines when their admittance is brought to zero, thus modelling effects of circuit breakers on the synchronous sub-system topology.

As a consequence, two or more disjoint connection graphs would be formed at the time of the breaker openings, each corresponding to a new synchronous island. The new graph topology should be analyzed at this point, picking a new root node for each newly formed island in the grid. Then, instead of having a single phase reference for the entire system, which would no longer be adequate, one would now have two or more independent phase references, one for each island, which would ensure that the variables of each island reach a steady state, thus avoiding the persistent sinusoidal oscillations found in the case of a statically determined connection topology.

From the point of view of the language specification, this extension requires to relax the rule in ​Section 8.3.4, by adding a statement such as the following one at the end of the section:

The only exception to this rule is when the if branch only contains Connection.branch() statements and equations exclusively involving the overconstrained variables of the connectors declared in the Connection.branch() statement. The rationale is that in this case, the overconstrained connector semantics takes care of always getting a balanced system.

and by adding the definition of the reconnectable() function to the overconstrained type definitions given in Section ​9.4.1.

From an implementation point of view, this means that parts of the graph connection analysis no longer can be performed statically at compile time, but need to be deferred to run time, because some unbreakable branches may be optionally activated or deactivated at runtime.

According to ​Section 9.4.2, everything except overconstrained connector variables is always handled in the same way, i.e., following ​Section 9.2, so the runtime analysis will be limited to overconstrained variables in connectors. Everthing else, including determining connection sets, building connection equations for flow variables, and handling stream connector variables, is unchanged and can still be managed statically at compile time.

In general, handling overconstrained connector variables dynamically at runtime could be quite a big deal, because Modelica tools normally cannot handle structural dynamics at runtime. However, the special case of scalar overconstrained connector variables with zero-dimensional equalityConstraint() function output turns out to be very easy to implement, and is enough to address the issue with the synchronous AC grids. This is demonstrated by the attached examplary case.

The attached package DynamicOverconstrainedConnectors2.mo contains components to build conceptual models of phasor-based AC power grid models, and will be used to demonstrate the proposed approach. The components are overly simplified and only retain the minimal features that are needed to analyze this language semantics extension.

In this case, the overconstrained variable is actually not the phase angle but rather the angular velocity of the phasor reference frame, which is sufficient to ensure that all state variables are constant when the synchronous islands are all operating in steady state, and avoids introducing extra unnecessary states to relate phase and frequency of different disconnected sub-grids.

System 1 is a basic example with two generators G1 and G2, each connected to a local load, connected by an inductive transmission line T. The system starts in steady state with balanced loads, then at time = 1 the active power consumption of load L2 is reduced by 20%. This starts a transient, that eventually dies out once the primary frequency controls stabilize the system at a new, higher frequency. Note that at that point, all state variables are constant, allowing implicit solvers to take large time steps.

System 2 is similar to System1, except that line T is split in the series connection of lines T1a and T1b, running in parallel, with line T2. This allows to demonstrate the overconstrained connector mechanism, since T1a and T1b form a loop that must be broken to avoid getting an overconstrained system of equations. The overall impedance is the same, so the power and frequency transients are the same as in System1.

In System 3, line T2 also includes a circuit breaker mechanism, that brings its susceptance B to zero when the T2.open signal becomes true. From that point onwards, there is no longer any exchange of power along that line, so the system is effectively split into two synchronous islands, one containing G1 and L1, the other one containing G2 and L2. The frequencies of the two drift apart, causing the phasors and state variables of G2 and L2, that still use G1.port.omegaRef as a reference angular velocity, to oscillate forever even when the new steady state is reached.

In order to avoid this, the connection graph needs to be split when the breaker embedded in T2 is opened. This behaviour is achieved in System4 by using the TransmissionLineVariableBranch model:

  type ReferenceAngularSpeed "Overconstrained type for per-unit reference angular speed"
    extends SI.PerUnit;
    function equalityConstraint
      input ReferenceAngularSpeed omega1;
      input ReferenceAngularSpeed omega2;
      output SI.PerUnit residue[0] "No constraints";
    end equalityConstraint;
  end ReferenceAngularSpeed;

  connector ACPort "Port for per unit 3-phase AC systems"
    SI.ComplexPerUnit v "Per unit voltage phasor referred to phase";
    flow SI.ComplexPerUnit i "Per unit current phasor referred to phase";
    ReferenceAngularSpeed omegaRef "Reference per-unit angular speed";
  end ACPort;

  model TransmissionLineVariableBranch "Purely inductive transmission line model with time-varying connection branch"
    parameter SI.PerUnit B = -5.0 "Line series per unit susceptance";
    discrete SI.PerUnit B_act "Actual value of per unit susceptance including breaker status";
    Boolean closed "State of line breaker";
    Boolean open = false "Command to open the line breaker";
    Boolean close = false "Command to close the line breaker";
    ACPort port_a;
    ACPort port_b;
  initial equation
    closed = true;
    B_act = B;
  equation
    port_a.i + port_b.i = Complex(0);
    port_a.i = Complex(0,B_act)*(port_a.v - port_b.v);
    when open then
      closed = false;
      B_act = 0;
    elsewhen close then
      closed = true;
      B_act = B;
    end when;
    // Equations affected by the proposed change to the specification:
    if closed then
      port_a.omegaRef = port_b.omegaRef;
      Connections.branch(port_a.omegaRef, port_b.omegaRef);
    end if;
  end TransmissionLineVariableBranch;

Handling the case of Real overdetermined connector variables with empty constraint function is relatively straighforward.

During flattening, the front end should indentify all the potentially connected subgraphs (i.e., ignoring the fact that some Connections.branch() statements are inside if-equations), then identify those of them that only contain static Connections.branch() statements, and those that contain at least one unbreakable branch declared within a conditional Connections.branch() statement.

The former subgraphs should be handled by the front-end as usual in Modelica 3.5. For the latter, the front-end will process the connection sets and the flow, effort and stream variable equations as usual; however, it should pass to the back-end the relevant information about overconstrained types and variables, and about nodes and branches of the corresponding subgraphs. Subgraphs that are structurally disconnected, regardless of the activation of some unbreakable branches, should be identified as such.

I guess we should try to re-use the front-end code in the backend as much as possible, and even in the runtime code (e.g. for graph analysis) as long as the MetaModelica code of the frontend can be translated into C/C++ code.

The generated code should then perform the following operations during initialization and each time any condition triggering a conditional equation containing Connections.branch statements changes:

• re-build the connection graphs of all the structurally unconnected sub-graphs to which the branches that changed their activation status belong, according to the rules set forth in ​Section 9.4, only accounting for branches that are actually active
• in case a previously de-activated branch is re-established, check that the corresponding connection equation involving the overconstrained connector variables is actually verified within some tolerance, otherwise abort the simulation; this would correspond to the requirement that two islands can only be reconnected if the are operating exactly at the same frequency or phase. This check can be defined by adding one more function reconnectable() to the overconstrained type definition, which returns a boolean indicating if the reconnection is feasible (see below).
• identify disjoint sub-graphs that are internally connected, and select one root node for each of them, according the standard rules set forth in ​Section 9.4.2
• break all loops in those sub-graphs
• set each overconstrained connector variable of each sub-graph to be equal to the value of the corresponding root node variable (this can be done trivially via pointers!)

Regarding the last two points, the back-end processing in the case of overconstrained variables being declared as equal in the unbreakable conditional branches, and of empty constraint equation, is trivial. At the end of the day, for each connected sub-graph, all the overconstrained connector variables are equal to the one of the selected root node. I guess generating code for this is also no big deal.

If a frequency reference is used as overconstrained connector variable (e.g. in the PowerGrids library, the reconnectable function could be implemented as

function reconnectable
  input ReferenceFrequency f1;
  input ReferenceFrequency f2;
  output Boolean canReconnect;
protected
  constant Real eps = 1e-6;
algorithm
  canReconnect := abs(f1 - f2) < eps;
end reconnectable;

For phase references, as in the PowerSystems library, one should also account for 2*pi ambiguity

function reconnectable
  input ReferenceAngle theta1;
  input ReferenceAngle theta2;
  output Boolean canReconnect;
protected
  constant Real eps = 1e-6;
  SI.Angle delta;
algorithm
  delta := (theta1 - theta2)- div(theta1 - thetha2, 2*pi)*2*pi;
  canReconnect := abs(delta) < eps;
end reconnectable;

I suppose this feature could be implemented rather easily in OpenModelica, as a first prototype implementation of this language extension. It could be tested in an extended version of the attached DynamicOverconstrainedConnectors2.mo package, including some more examples with a bit more nodes and edges in the connection graphs, and finally tested on an extended version of the ​PowerGrids library, that could be very easily extended to include this feature, by just changing a few lines of code.

Another application of this extension could be incompressible fluid systems. These systems require a component setting the pressure for each connected subsystem, corresponding to a pressurizer or expansion tank. When closing some strategic valves, it is possible to split the circuit in two or more unconnected parts, which should of course have at least one expansion tank each. At that point, new pressure-setting components (i.e. root nodes) should be identified, one for each newly formed subsystem.

I don't yet have the code ready for this system, but I can produce it overnight if we proceed with this proposal.

If the experiments with the prototype are successful, this text could form the basis of an MCP to introduce this feature into the next version of the Modelica Specification.

[casella]

Ticket retargeted after milestone closed