Removing structurally non-observable variables

[casella]

In many cases, when dealing with performance-critical models, there are some parts that are not of interest for some simulations, so their computation should be skipped. One way to handle this, which is used e.g. for visualization features in the Modelica MultiBody library, is via conditional components, but this is cumbersome and needs to be programmed explicitly.

A better idea is to declare what are the outputs of interest for the simulation, e.g. via annotations or making them top-level outputs of the model, performing a structural observability analysis, and only keep the parts of the model that are strictly necessary to compute those outputs. In the case of electrical system, one could get rid of several auxiliary variables (e.g. angles and moduli of voltage and current phasors) that are not of interest and take some time to compute because of the arctan and square root functions.

A similar problem is already solved by the backend to separate equations that are not needed to compute the derivatives and the zero-crossing functions,, which are only solved when the step is accepted. However, at that point all other variables are computed.

An algorithm to perform this operation is available in this paper. It can be easily implemented as an extension of existing matching algorithms in the OpenModelica backend. This could also be used in more sophisticated cases, e.g. to remove off-plane dynamic equations (and states!) from a planar model built with the 3D MultiBody library, provided that the sensors are only capturing in-plane motion and forces.

[phannebohm]

While starting to assess the implementation effort for this, I found the flag -d=evalOutputOnly by @vwaurich, which apparently already does what is asked for in this ticket. I made a couple changes so that only equations for top-evel outputs are computed, and subsystems without top-level outputs are removed completely.

@casella is that what you were looking for? I only tested very small MWEs. Could you please perform some more elaborate tests yourself?

[phannebohm]

An algorithm to perform this operation is available in this paper.

BTW, the old implementation for -d=evalOutputOnly is not as smart as the algorithm in your paper. It acts after the usual causalization and runs a fill algorithm on the "task graph" of the system to remove unneeded equations. Only sorting needed parts might be a bit more performant and I don't see any drawbacks. @kabdelhak, we should directly do that in the new backend.

[casella]

While starting to assess the implementation effort for this, I found the flag -d=evalOutputOnly by @vwaurich, which apparently already does what is asked for in this ticket. I made a couple changes so that only equations for top-evel outputs are computed

This is one good idea, particularly for FMI. The only problem is that it requires models to be built specifically for this purpose. You can't make a variable of an existing model an output, without changing the existing model's source code. Of course you can always inherit and add output aliases.

I wonder if we could also have an option to put the list of "outputs" as a regexp in a compiler flag, so we don't need to change the source code for that.

and subsystems without top-level outputs are removed completely.

I'm not sure this is always a good idea. Unless the sub-system has protected variables, other sub-systems can access it via dot-notation, and also with other top-level connectors that are not outputs. In this case, you would break them. Do I miss something?

Could you please perform some more elaborate tests yourself?

Will do next year. If you are interested, you can start with two cases which are trivial to set up, and can be scaled up to arbitrarily large scale.

ScalableTestGrids.Models contains power grid models where phase angles are computed at all buses and transformer ports. These variables require computing arctangents (which are expensive) and are not required to advance the simulation. You can inherit from them and add output Real omega = GEN_1_1.omega;. That's it. Computing this output requires to solve the entire power flow, but will leave out those phase outputs, as well as a lot of numerous other auxiliary outputs, e.g. per unit quantities and stuff like that, which is useful for monitoring purposes but unnecessary to carry on the simulation.

ScalableTestSuite.Mechanical.Strings.ScaledExperiments contains models of strings suspended from a slider to which an impulse is applied at some point in time, generating travelling waves. The motion happens entirely in the x-y plane. My conjecture is that if you extend them and add output Real p[2] = bodybox[N+1].frame_b.r_0[1:2];, i.e., only observe the x-y coordinates of the loose end of the string, you should be able to get rid of all the equations relative to out-of-plane motion. To do so, it may be necessary to set -d=evaluateAllParameters.

[casella]

BTW, the old implementation for -d=evalOutputOnly is not as smart as the algorithm in your paper.

Our algorithm is actually some kind of "structural observability analysis".

It acts after the usual causalization and runs a fill algorithm on the "task graph" of the system to remove unneeded equations. Only sorting needed parts might be a bit more performant and I don't see any drawbacks.

Neither do I. In fact, I'm still not sure that the existing algorithm obtains the same result. Does it also remove states and differential equations, if those are not influencing the top-level outputs?

@kabdelhak, we should directly do that in the new backend.

Let me know how it works :)

[phannebohm]

and subsystems without top-level outputs are removed completely.

I'm not sure this is always a good idea. Unless the sub-system has protected variables, other sub-systems can access it via dot-notation, and also with other top-level connectors that are not outputs. In this case, you would break them. Do I miss something?

Pardon, I'm talking about sub-systems in the mathematical sense. If a variable influences another they are both part of the same system. Otherwise you would be right of course. The removal of whole sub-systems was done in the beginning and then reverted shortly after, so there might be something I overlooked, but from my understanding this is valid.

[...] I'm still not sure that the existing algorithm obtains the same result. Does it also remove states and differential equations, if those are not influencing the top-level outputs?

It claims to do so and I checked a few cases which confirm that. For example this one

model MWE
  output Real out;
  Real x(start=0, fixed=true);
  Real y, w;
equation
  der(x) = -x;
  y = sin(time);
  out = x^2;
  w + x = time;
end MWE;

removes y because it is a separate sub-system (in the mathematical sense) and w because it is not needed to compute out, so it basically reduces to

der(x) = -x;
out = x^2;

but if I change out = x^2; to out = y^2;, then y is kept and x and der(x) are removed.

[casella]

Pardon, I'm talking about sub-systems in the mathematical sense.

Aha 😅

If a variable influences another they are both part of the same system. Otherwise you would be right of course. The removal of whole sub-systems was done in the beginning and then reverted shortly after, so there might be something I overlooked, but from my understanding this is valid.

I'm not sure exactly what you mean by sub-system in this context, but I assume what you did is ok. Maybe it will be described formally in some paper sooner or later?

[...] I'm still not sure that the existing algorithm obtains the same result. Does it also remove states and differential equations, if those are not influencing the top-level outputs?

It claims to do so and I checked a few cases which confirm that. For example this one

model MWE
  output Real out;
  Real x(start=0, fixed=true);
  Real y, w;
equation
  der(x) = -x;
  y = sin(time);
  out = x^2;
  w + x = time;
end MWE;

removes y because it is a separate sub-system (in the mathematical sense) and w because it is not needed to compute out, so it basically reduces to

der(x) = -x;
out = x^2;

but if I change out = x^2; to out = y^2;, then y is kept and x and der(x) are removed.

OK, very good.

[casella]

I'll test it ASAP

[casella]

I ran some tests, see the attached test package TestNonObservable.mo.txt

In ODE mode, the feature works fine. Unfortunately, DAE mode doesn't have it. MWE1 and MWE2 demonstrate that clearly on a small test case. The problem is, RTE needs DAE mode for their models. Small models in principle should also run in ODE mode, but when I tried it on an example of ScalableTestGrids, it failed at the start of the simulation.

We basically have two options:

1. port the existing -d=evalOutputOnly feature to daeMode in the old backend,
2. implement my algorithm in the new backend for both odeMode and daeMode,

[casella]

@kabdelhak there is still some issue here. I tried MW2 and it worked, then I realized there were no differential equations, so euler was actualy ran. not daeMode. I tried MW3 with an output involving a state, and indeed it worked. Here is the updated package

TestNonObservable.mo.txt

However, when I tested TestType1, it doesn't still work fine. TestType1Full_daeMode is the base test, which works fine. Before your commit, TestType1Reduced_daeMode was actually not reduced at all, but it worked fine, while TestType1Reduced_odeMode failed at initialization. After your commit, also TestType1Reduced_daeMode fails at initialization, in the same way and probably for the same reason the odeMode model already failed previously.

I made a quick check with the debugger. Apparently, the application of the reduction algorithm changes the structure of the strong components, eventually causing the convergence failure. Looking at the initial-lambda0 set of equations, which is the first to be solved during initialization (and which fails), the TestType1Full_daeMode model has one strong component with these iteration variables:

and a separate one with these iteration variables

while the TestType1Reduced_daeMode has a single strong component with variables coming from the two previously shown sets:

This makes no sense, removing some equations that are not needed should not introduce couplings between variables that were originally computed in two separate strong components.

Can you please check that? Thanks!

[kabdelhak]

Well for this i would need to debug the whole module, let's see how long that takes. But at least it seems like it is broken in general and not only for DAEMode.

[casella]

Well for this i would need to debug the whole module

I guess so. I could try to provide a smaller MWE if that helps you, please let me know.

But at least it seems like it is broken in general and not only for DAEMode.

Absolutely.

[kabdelhak]

I guess so. I could try to provide a smaller MWE if that helps you, please let me know.

yes please, this one is hard to debug.

Also you i think you are right, removing non output dependent stuff should not impact algebraic loops. Initialization on the other hand might be different since we only look at the simulation system while removing stuff. It might be that fixed start attribute variables have been removed which impacts the way the initial system has to be solved. I am not entirely sure that it is wrong what is happening here.

[casella]

Also you i think you are right, removing non output dependent stuff should not impact algebraic loops. Initialization on the other hand might be different since we only look at the simulation system while removing stuff.

Your're right. In fact this may cause some initial equations to become redundant, if they are not given as fixed = true attributes on removed variables, but involve actual initial equations. Not entirely sure how to handle this, we should think about it. But in any case this is not relevant for RTE, we are only removing algebraic equations on the leafs of the dependency graph, not states.

It might be that fixed start attribute variables have been removed which impacts the way the initial system has to be solved.

I can't see how removing initial equations can add dependencies between formerly independent strong components. You'd have to add equations to get that outcome.

I am not entirely sure that it is wrong what is happening here.

Me neither 😅

[kabdelhak]

While trying to come up with a model to explain my point of structural changes to the initial system due to removed variables, i stumbled into a minimal model that shows the problem pretty well, it seems like.

model test_remove
  Real x;
  output Real y,z;
initial equation
  x = y;
equation
  x = sin(time);
  der(z) = y - 3;
  y = z^2 - 10;
end test_remove;

My point being is that in the simulation system x can be removed, but from the initial system it cannot, since y depends on x in that case. We don't seem to check for that. Simulating this model works fine, but adding -d=evalOutputOnly returns

test_remove_06inz.c:20:56: error: use of undeclared identifier 'x'
  (data->localData[0]->realVars[3] /* y variable */) = x;
                                                       ^
1 error generated.
make: *** [<builtin>: test_remove_06inz.o] Error 1
make: *** Waiting for unfinished jobs....
Compilation process failed. Exited with code 2.

which is a direct of x being removed but being needed for initialization. The module is quite complicated, uses some tasking graphs which i would need to understand first. I don't know if it's easy to just also check the initial system.

[kabdelhak]

As an addition to the previous comment i want to state that i could also remove all initial equations containing any removed variable. Is there a problem in always doing that?

[casella]

While trying to come up with a model to explain my point of structural changes to the initial system due to removed variables, i stumbled into a minimal model that shows the problem pretty well, it seems like.

In principle this is a legitimate model, though it seems to me that writing initial equations that depend on variables that not interesting for the output seems quite unlikely in real applications.

As an addition to the previous comment I want to state that i could also remove all initial equations containing any removed variable. Is there a problem in always doing that?

I'm not sure. A safer approach would probably be to add all the variables that show up in initial equations and variables with fixed = true to the set of outputs, so the problem is avoided a-priori. However, this would make the algorithm unnecessarily cumbersome in most cases.

A better approach would probably be to initially ignore initial equations and fixed = true variables when reducing the system. Then you can check if any eliminated variables shows up in initial equations. If none does (I guess that's >95% of the cases) you are done, otherwise you add those who do to the set of outputs and repeat.

We have to be careful that these checks do not unfavourably scale for large systems. A brute-force implementation would scale as (N*M) where N is the number of eliminated variables and M is the number of initial equations. Not sure if you can make it more efficient than that. On the other hand, the number of initial equations in most real-life models is a fraction of the total number of equations. Also, this optimization could be done on an array-preserving representation of the system, which would probably reduce the complexity in a drastic way.

[casella]

Thanks @kabdelhak, do all tests in the attached package work fine now?

[phannebohm]

Oops, MWE2 still fails. It contains a state variable that is unneeded for the outputs.

Our approach in PR #10247 converts all unneeded variables to parameters so that they are still accessible for initialization. (A second step could be to remove all unneeded parameters from the initial system, while keeping those needed in the simulation system.)

But since the old backend does not properly distinguish between states and their derivatives, we get in trouble when changing a state to parameter.

Also we didn't consider the case that a derivative is needed but the corresponding state isn't. In general all four combinations of state/derivative and needed/unneeded are possible, but that can only be handled by the new backend. In the old backend we will either keep both or scrap both.

[casella]

Good. I'll try it out on a larger grid example, and if it works we're in business. Thanks!

[casella]

@phannebohm, @kabdelhak, I checked the result of the power system MWE thoroughly. The method works fine, results match the ones of the full model, and the gain in efficiency is substantial, as the number of simulation equations is slashed by half and the overall simulation time is reduced by 33%.

However, the strategy of #10247 has a questionable outcome, from an end-user point of view: the variables removed from the simulation problem and re-added as parameters to the initial equations are still visible in the simulation results, where they show up as parameters with a constant value.

This is ok in the MWE I gave you, which is a steady-state simulation. However, during transient simulations, those variables would in general be time-varying. If they can be displayed in the simulation results, but are then shown to remain constant throughout the simulation, that is plain wrong, and very confusing too.

My suggestion for the time being is to set the visibility attribute of all those variables-turned-parameters to hidden, so they are not saved in the .mat file and not displayed at all in the simulation results, that is unnecessarily cluttered by this mostly uninteresting stuff. In the power system MWE, that would mean getting rid of half of the variables in the variable browser.

Can you do that?

If in the future a method is implemented in the NB to get rid of those that are not necessary to solve the initialization problem of the reduced model (see #10262), then we could think of actually displaying the surviving ones, but then their names should be modified, e.g. adding the $initial suffix, so that, e.g. port.VPu becomes port.VPu$initial. In this way it will still be possible to visualize variables that were relevant for the initialization problem, but it will be clear that they are constant because they are only used during initialization.

[casella]

Great, thanks! I'll test it ASAP.

[casella]

See results here. The impact on simulation time is quite good, some 40% reduction in simulation time. Unfortunately there are serious scaling issues, see #10413.