Models Buildings.ThermalZones.ISO13790.Examples.xxx fails verification due to wrong expansion of the scalar product inside a function

[AndreaBartolini]

The following models fail the verification vs the Dymola reference files:

• Buildings.ThermalZones.ISO13790.Examples.FreeFloating
• Buildings.ThermalZones.ISO13790.Examples.FreeFloatingHVAC
• Buildings.ThermalZones.ISO13790.Examples.HeatingCooling
• Buildings.ThermalZones.ISO13790.Examples.HeatingCoolingHVAC

The root cause is the same for all the models, the analysis is made with reference to the model Buildings.ThermalZones.ISO13790.Examples.FreeFloatingHVAC.

The variable that fails the verification is the temperature zonHVAC.TAir, as shown below:

the failure is caused by the wrong expansion of the scalar product y :=a*xp inside the function Buildings.Utilities.Math.Functions.polynomial, that implements the polynomial interpolation:
y = a1 + a2*x + a3*x^2 + a4*x^3 ...

The function implementation is the following:

function polynomial "Polynomial function"
  extends Modelica.Icons.Function;
 input Real x "Independent variable";
 input Real a[:] "Coefficients";
 output Real y "Result";
protected
 parameter Integer n = size(a, 1)-1;
 Real xp[n+1] "Powers of x";
algorithm
  xp[1] :=1;
  for i in 1:n loop
     xp[i+1] :=xp[i]*x;
  end for;
  y :=a*xp;
end polynomial;

In all the above listed examples the array of the coefficients a contains only one element:
a = {1}

so the expected output is y=1, constant for the entire transient, while the output produced by omc is y=0, constant for the entire transient.

If inside the function the scalar product y :=a*xp; is substituted by its expansion y := sum(a[i]*xp[i] for i in 1:size(a,1)); then the output becomes the right one ( y=1, constant for the entire transient) and the transient overlaps the reference.

I tried to analyze the two cases by means of the transformational debugger, if the scalar product y :=a*xp; is used then inside the transformational debugger only the initial assignment is present:

initial (assign) zon5R1C.win.polynomial[1].y := 0.0

and in principle it is what is expected because the output is constant for the entire transient, what is wrong is the value 0.0 instead of 1.0

If the expansion y := sum(a[i]*xp[i] for i in 1:size(a,1)); is used then the function is not simplified and in the transformational debugger we found the following regular assignement:

regular (assign) zonHVAC.win.polynomial[1].y := Buildings.Utilities.Math.Functions.polynomial(zonHVAC.win.HDirTil[1].inc, zonHVAC.win.polynomial[1].a)

that calls the function every simulation steps.

Starting from here, it seems that the problem should be in some simplification made by the BE. I tried to compile the model by using -d=optdaedump but after 20min the translation was still in progress and I stopped the compiler.

Building master branch a18ff66b4239fc9816cc319441f79b26ae6c4c40
OpenModelica v1.26.0-dev-59-g5f10fe7efc (64-bit)
WIN 11 PRO (64 bit)

[casella]

@phannebohm, this is on you 😃

[phannebohm]

Off topic but may I suggest you use Horner's method for evaluating polynomials:

function horner "Polynomial evaluation"
  extends Modelica.Icons.Function;
  input Real x "Independent variable";
  input Real a[:] "Coefficients";
  output Real y "Result";
algorithm
  y := 0;
  for i in size(a, 1):-1:1 loop
    y := y*x + a[i];
  end for;
end horner;

It uses less memory and is AFAIK numerically more stable than the naive approach. Or is there some reason not to use that?

[phannebohm]

The function horner as above works for the case a = {1} in the old and new backend. The function polynomial as of now only works with the new backend. I will investigate why the old backend produces these bad results. Thanks @AndreaBartolini for the thorough analysis so far!

[phannebohm]

I found the problem:
The preopt module evalFunc saw that xp[1] is set to a constant so it added the replacement rule xp[1] -> 1. The replacement then saw a * xp and tried to scalarize xp into its elements like {xp[1], xp[2], ...} in order to apply the constant replacement. However because the size of xp (which is n+1) is not constant, it actually "scalarized" xp into the empty array {}. Then the simplify module saw a * {} and simplified that to 0.0.

I opened #13963 to not replace arrays if their size is unknown. That should fix these models.