Inefficient translation of models with vectorized components

[Rüdiger Franke]

This issue follows up #5110. It sketches an important use case for arrays: vectorized model components. The following model was created using OMEdit (stripped graphical annotations here):

model Vectorized
  parameter Integer n = 3;
  Modelica.Blocks.Continuous.Integrator integrator1[n];
  Modelica.Blocks.Interfaces.RealInput u[n];
  Modelica.Blocks.Interfaces.RealOutput y[n];
equation
  connect(u, integrator1.u);
  connect(integrator1.y, y);
  annotation(
    uses(Modelica(version = "3.2.2")));
end Vectorized;

omc -d=newInst,-nfScalarize only declares u and y as arrays and rolls out everything else:

class Vectorized
  parameter Integer n = 3;
  parameter Real integrator1[1].k(unit = "1") = 1.0 "Integrator gain";
  parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput) integrator1[1].initType = Modelica.Blocks.Types.Init.InitialState "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real integrator1[1].y_start = 0.0 "Initial or guess value of output (= state)";
  Real integrator1[1].u "Connector of Real input signal";
  Real integrator1[1].y(start = integrator1[1].y_start) "Connector of Real output signal";
  parameter Real integrator1[2].k(unit = "1") = 1.0 "Integrator gain";
  parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput) integrator1[2].initType = Modelica.Blocks.Types.Init.InitialState "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real integrator1[2].y_start = 0.0 "Initial or guess value of output (= state)";
  Real integrator1[2].u "Connector of Real input signal";
  Real integrator1[2].y(start = integrator1[2].y_start) "Connector of Real output signal";
  parameter Real integrator1[3].k(unit = "1") = 1.0 "Integrator gain";
  parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput) integrator1[3].initType = Modelica.Blocks.Types.Init.InitialState "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real integrator1[3].y_start = 0.0 "Initial or guess value of output (= state)";
  Real integrator1[3].u "Connector of Real input signal";
  Real integrator1[3].y(start = integrator1[3].y_start) "Connector of Real output signal";
  input Real[3] u;
  output Real[3] y;
initial equation
  integrator1[3].y = integrator1[3].y_start;
  integrator1[2].y = integrator1[2].y_start;
  integrator1[1].y = integrator1[1].y_start;
equation
  u[1] = integrator1[1].u;
  u[2] = integrator1[2].u;
  u[3] = integrator1[3].u;
  integrator1[1].y = y[1];
  integrator1[2].y = y[2];
  integrator1[3].y = y[3];
  der(integrator1[1].y) = integrator1[1].k * integrator1[1].u;
  der(integrator1[2].y) = integrator1[2].k * integrator1[2].u;
  der(integrator1[3].y) = integrator1[3].k * integrator1[3].u;
end Vectorized;

The desired result should read instead:

class VectorizedDesired
  parameter Integer n = 3;
  parameter Real integrator1.k[n](unit = "1") = 1.0 "Integrator gain";
  parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput) integrator1.initType[n] = Modelica.Blocks.Types.Init.InitialState, "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real integrator1.y_start[n] = 0.0 "Initial or guess value of output (= state)";
  Real integrator1.u[n] "Connector of Real input signal";
  Real integrator1.y[n](start = integrator1.y_start) "Connector of Real output signal";
  input Real[3] u;
  output Real[3] y;
initial equation
  integrator1.y = integrator1.y_start;
equation
  u = integrator1.u;
  integrator1.y = y;
  der(integrator1.y) = integrator1.k .* integrator1.u;
end VectorizedDesired;

A couple of interesting questions arises:

1. Should array dimensions be "normalized" by collecting them at the end, i.e. integrator1.u[n] instead of integrator1[n].u?
2. Can n be kept as structural parameter during the translation?
3. Is each required for Real integrator1.y[n](start = 0.0). The current implementation (see #5110) works fine without each.
4. Are initializations of arrays with scalars allowed, i.e. parameter Real integrator1.y_start[n] = 0.0? This would simplify the generation of efficient code, i.e. a for loop assigning 0.0 to each array element instead of creation of a temporary array resulting from fill(0.0, n).
5. Note: the operator changes from * to .*! This appears quite tricky if possible at all in the general case.
6. Shall the frontend better introduce a for loop for all equations of a vectorized component and the backend keep that loop, analyzing only one iteration of it?

   ...
   equation
     for i in 1:n loop
       u[i] = integrator1.u[i];
       integrator1.y[i] = y[i];
       der(integrator1.y[i]) = integrator1.k[i] * integrator1.u[i];
     end for;
   end VectorizedDesired;
   

[Adrian Pop]

I guess this should be partially possible. We expand the complex components using the type array dimensions here:
​https://github.com/OpenModelica/OMCompiler/blob/master/Compiler/NFFrontEnd/NFFlatten.mo#L507
so that should probably be the starting point.

There are several possibilities to handle this, one would be to have some sort of for loops for component declarations and equations:

class VectorizedDesired
  parameter Integer n = 3;
  for in in 1:n
    parameter Real integrator1[i].k(unit = "1") = 1.0 "Integrator gain";
    parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput) integrator1[i].initType = Modelica.Blocks.Types.Init.InitialState, "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
    parameter Real integrator1[i].y_start = 0.0 "Initial or guess value of output (= state)";
    Real integrator1[i].u "Connector of Real input signal";
    Real integrator1[i].y(start = integrator1.y_start) "Connector of Real output signal";
  end for;
  input Real[3] u;
  output Real[3] y;
initial equation
  for i in 1:n loop
    integrator1[i].y = integrator1[i].y_start;
  end for;
equation
  for i in 1:n loop
    u[i] = integrator1[i].u;
    integrator1[i].y = y[i];
    der(integrator1[i].y) = integrator1[i].k * integrator1[i].u;
  end for;
end VectorizedDesired;

I think that moving the array dimensions at the end of the hierarchical component reference is not that good because:

• the component might contain array components of basic type
• the handling of bindings would become more complicated
• one would need interception of component references containing indexes, i.e. integrator1[2].u -> integrator1.u[2]

At least if we use for loops we clearly know what the structure was and where the actual array dimensions are. At the end we might be able to "normalize" the components (collect array dimensions at the end of the hierarchical component reference) but we can keep the array dimension order (which one is first and such) because we know if from the iterator index.

Another solution (similar to yours but more correct Modelica) would be to move all the dimension to the type in the correct order (we should know the order and the dimensions from the typing).

class VectorizedDesired
  parameter Integer n = 3;
  parameter Real[n] integrator1.k(unit = "1") = 1.0 "Integrator gain";
  parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput)[n] integrator1.initType = Modelica.Blocks.Types.Init.InitialState, "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real[n] integrator1.y_start = 0.0 "Initial or guess value of output (= state)";
  Real[n] integrator1.u "Connector of Real input signal";
  Real[n] integrator1.y(start = integrator1.y_start) "Connector of Real output signal";
  input Real[3] u;
  output Real[3] y;
initial equation
  integrator1.y = integrator1.y_start;
equation
  u = integrator1.u;
  integrator1.y = y;
  der(integrator1.y) = integrator1.k * integrator1.u;
end VectorizedDesired;

We'll see what Per thinks as well about this and we'll decide how to proceed.

[Adrian Pop]

With PR: ​https://github.com/OpenModelica/OMCompiler/pull/2682 we get:

//adrpo33@ida-0030 MINGW64 /c/home/adrpo33/dev/OMTesting/bugs/5144
//$ ~/dev/OpenModelica/build/bin/omc -d=newInst,-nfScalarize v.mos
class Vectorized
  parameter Integer n = 3;
  Real[3] integrator1.y(start = integrator1.y_start) "Connector of Real output signal";
  Real[3] integrator1.u "Connector of Real input signal";
  parameter Real[3] integrator1.y_start = 0.0 "Initial or guess value of output (= state)";
  parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput)[3] integrator1.initType = Modelica.Blocks.Types.Init.InitialState "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real[3] integrator1.k(unit = "1") = 1.0 "Integrator gain";
  input Real[3] u;
  output Real[3] y;
initial equation
  integrator1.y = integrator1.y_start;
equation
  u[1] = integrator1[1].u;
  u[2] = integrator1[2].u;
  u[3] = integrator1[3].u;
  integrator1[1].y = y[1];
  integrator1[2].y = y[2];
  integrator1[3].y = y[3];
  der(integrator1.y) = integrator1.k * integrator1.u;
end Vectorized;

I tried to get rid of connect expansion but that is not that easy. If I disable the expansion here: ​https://github.com/OpenModelica/OMCompiler/blob/master/Compiler/NFFrontEnd/NFConnections.mo#L146
the model fails to flatten. I guess connection sets handling doesn't like arrays.

[Rüdiger Franke]

The vectorized equation

  der(integrator1.y) = integrator1.k * integrator1.u;

should use element wise multiplication:

  der(integrator1.y) = integrator1.k .* integrator1.u;

There are for sure more issues in the general case. This is why I thought that a for loop might work out better for vectorized components.

[Per Östlund]

You get some declarations that are incorrect now, like (notice the scalar bindings):

parameter Real[3] integrator1.k(unit = "1") = 1.0 "Integrator gain";

I guess generating fill calls might work, i.e.:

parameter Real[3] integrator1.k(unit = fill("1", 3)) = fill(1.0, 3) "Integrator gain";

Replying to adrpo:

I tried to get rid of connect expansion but that is not that easy. If I disable the expansion here: ​https://github.com/OpenModelica/OMCompiler/blob/master/Compiler/NFFrontEnd/NFConnections.mo#L146
the model fails to flatten. I guess connection sets handling doesn't like arrays.

Yes, the connection handling only expects scalars with constant subscripts. Handling arrays would be much more complicated, since you can have cases such as:

model M
  C c1[3], c2[3], c3;
equation
  connect(c1, c2);    // Connect arrays c1 and c2.
  connect(c1[1], c3); // Connect only part of c1 with c3.
end M;

The first connect would create the set {c1, c2} if we don't expand it, while the second connect would require us to split that set so that we get {c1[1], c2[1], c3}, {c1[2], c2[2]} and {c1[3], c2[3]} (you could also use e.g. c1[2:3] instead of expanding it fully, but that would be a nightmare to handle efficiently). This would of course also mean that we'd have to look connectors up twice, first without subscripts, and then with subscripts if not found.

[Adrian Pop]

Yes, PR: ​https://github.com/OpenModelica/OMCompiler/pull/2682 fixes just the types of the variables, nothing else. It needs more work to fix equations and bindings.

[Rüdiger Franke]

Replying to perost:

You get some declarations that are incorrect now, like (notice the scalar bindings):

parameter Real[3] integrator1.k(unit = "1") = 1.0 "Integrator gain";

I guess generating fill calls might work, i.e.:

parameter Real[3] integrator1.k(unit = fill("1", 3)) = fill(1.0, 3) "Integrator gain";

It would be important to propagate the information that each array element has the same value to the code generation, something like:

parameter Real[3] integrator1.k(each unit = "1", each value = 1.0);

An array structure with fill-operator should also work, but would need to be implemented in the OpenModelica data structures:

parameter Real[3] integrator1.k(each unit = "1") = fill(1.0, 3);

Avoid today's expansions, like:

parameter Real[3] integrator1.k(unit = {"1", "1", "1"}) = {1.0, 1.0, 1.0};

because this requires the unnecessary treatment of potentially large arrays in the simulation code.

[Rüdiger Franke]

As for loops do work in the backend to some extend now (see #5110), it would be interesting to see them introduced by the frontend for arrays of components ...

[Rüdiger Franke]

The following model works with omc -d=newInst,-nfScalarize --simCodeTarget=Cpp --std=3.3. Note that

• I added a clock as the backend is more gentle to such models and discrete-time models are my primary target anyway (could be achieved with a block Modelica_Synchronous.RealSignals.Sampler.AssignClockVectorized)
• used for loops that appear easier (contain a pattern of a scalar equation -- no problems with dot products, vectorized if conditions and the like)

Couldn't the new frontend create such for loops with -nfScalarize for arrays of component models?

model VectorizedFor
  parameter Integer n = 3;
  Real[3] integrator1_y(start = integrator1_y_start) "Connector of Real output signal";
  Real[3] integrator1_u "Connector of Real input signal";
  parameter Real[3] integrator1_y_start = fill(0.0, n) "Initial or guess value of output (= state)";
  //parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput)[3] integrator1_initType = Modelica.Blocks.Types.Init.InitialState "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real[3] integrator1_k(each unit = "1") = fill(1.0, n) "Integrator gain";
  input Real[3] u(each start = 1);
  output Real[3] y;
  Real[3] uAssignClock;
  Clock clock = Clock(Clock(0.1), solverMethod="ImplicitEuler");
initial equation
  integrator1_y = integrator1_y_start;
equation
  for i in 1:n loop
    uAssignClock[i] = sample(u[i], clock);
  end for;
  for i in 1:n loop
    uAssignClock[i] = integrator1_u[i];
    integrator1_y[i] = y[i];
  end for;
  for i in 1:n loop
    der(integrator1_y[i]) = integrator1_k[i] * integrator1_u[i];
  end for;
end VectorizedFor;

[Rüdiger Franke]

PR2712 (​https://github.com/OpenModelica/OMCompiler/pull/2712) builds on PR2682 above. It keeps arrays in connections (works for signals) and vectorizes equations with for loops per default. The following example:

model VectorizedClocked
  constant Integer n = 10;
  Modelica_Synchronous.ClockSignals.Clocks.PeriodicRealClock periodicClock1(period = 0.1, useSolver=true, solverMethod="ImplicitEuler");
  Modelica.Blocks.Interfaces.RealInput u[n](start = 1:n);
  Modelica.Blocks.Continuous.Integrator integrator1[n](y_start = 1:n);
  Modelica.Blocks.Interfaces.RealOutput y[n];
  Modelica_Synchronous.RealSignals.Sampler.AssignClockVectorized assignClock1(n = n);
equation
  connect(y, integrator1.y);
  connect(integrator1.u, assignClock1.y);
  connect(u, assignClock1.u);
  connect(periodicClock1.y, assignClock1.clock);
  annotation(
    uses(Modelica(version = "3.2.2"), Modelica_Synchronous(version = "0.92.1")));
end VectorizedClocked;

produces with omc --std=3.3 -d=newInst,-nfScalarize:

class VectorizedClocked
  constant Integer n = 10;
  parameter Real periodicClock1.period(quantity = "Time", unit = "s") = 0.1 "Period of clock (defined as Real number)";
  parameter Boolean periodicClock1.useSolver = true "= true, if solverMethod shall be explicitly defined";
  parameter String periodicClock1.solverMethod = "ImplicitEuler" "Integration method used for discretized continuous-time partitions";
  Clock periodicClock1.y;
  input Real[10] u(start = /*Real[n]*/(1:10));
  Real[10] integrator1.y(start = integrator1.y_start) "Connector of Real output signal";
  Real[10] integrator1.u "Connector of Real input signal";
  parameter Real[10] integrator1.y_start = /*Real[n]*/(1:10) "Initial or guess value of output (= state)";
  parameter enumeration(NoInit, SteadyState, InitialState, InitialOutput)[10] integrator1.initType = Modelica.Blocks.Types.Init.InitialState "Type of initialization (1: no init, 2: steady state, 3,4: initial output)";
  parameter Real[10] integrator1.k(unit = "1") = 1.0 "Integrator gain";
  output Real[10] y;
  parameter Integer assignClock1.n(min = 1) = 10 "Size of input signal vector u (= size of output signal vector y)";
  Real[10] assignClock1.u "Connector of clocked, Real input signal";
  Real[10] assignClock1.y "Connector of clocked, Real output signal";
  Clock assignClock1.clock;
initial equation
  integrator1.y = integrator1.y_start;
equation
  y = integrator1.y;
  integrator1.u = assignClock1.y;
  u = assignClock1.u;
  periodicClock1.y = assignClock1.clock;
  periodicClock1.y = Clock(Clock(0.1), "ImplicitEuler");
  for $i in 1:10 loop
    der(integrator1.y[$i]) = integrator1.k[$i] * integrator1.u[$i];
  end for;
  when assignClock1.clock then
    assignClock1.y = assignClock1.u;
  end when;
end VectorizedClocked;

The vectorized model simulates with --preOptModules-=inlineArrayEqn --postOptModules-=inlineArrayEqn,solveSimpleEquations --simCodeTarget=Cpp.

[Rüdiger Franke]

After some more commits, the vectorized model simulates without backend tweaks; just use

-d=newInst,-nfScalarize --simCodeTarget=Cpp

Generated for-equations place subscripts at the right places now, i.e. integrator1[$i].u instead of integrator1.u[$i] in the above example to avoid misunderstandings (see ​https://github.com/OpenModelica/OMCompiler/commit/789a1ecd90789f47c06317f042dc3af4d39e6b87).

Three new cases have been added to the testsuite (​https://github.com/OpenModelica/OpenModelica-testsuite/tree/master/openmodelica/cppruntime):

• testVectorizedBlocks.mos -- above example
• testVectorizedSolarSystem.mos -- physical connectors with flow variables
• testVectorizedPowerSystem.mos -- array equations that need to be solved after inlining

The latter case still phases problems with the backend module removeSimpleEquations. Need to additionally use:

--preOptModules-=removeSimpleEquations --postOptModules-=removeSimpleEquations

But this is out of the scope of this ticket.