smoothMin does not respect minimum

[Mathadon]

The following example shows that smoothMin (and smoothMax, smoothLimit) need not respect that the output should be the minimum of the two inputs. Furthermore a local maximum exists.

model test

  Annex60.Utilities.Math.SmoothMin smoothMin(deltaX=1)
    annotation (Placement(transformation(extent={{-40,0},{-20,20}})));
  Modelica.Blocks.Sources.Ramp ramp(height=10, duration=1)
    annotation (Placement(transformation(extent={{-98,8},{-78,28}})));
  Modelica.Blocks.Sources.Constant const(k=1)
    annotation (Placement(transformation(extent={{-100,-24},{-80,-4}})));
equation 
  connect(ramp.y, smoothMin.u1) annotation (Line(
      points={{-77,18},{-60,18},{-60,16},{-42,16}},
      color={0,0,127},
      smooth=Smooth.None));
  connect(const.y, smoothMin.u2) annotation (Line(
      points={{-79,-14},{-54,-14},{-54,4},{-42,4}},
      color={0,0,127},
      smooth=Smooth.None));
  annotation (uses(Annex60(version="0.1"), Modelica(version="3.2.1")), Diagram(
        coordinateSystem(preserveAspectRatio=false, extent={{-100,-100},{100,100}}),
        graphics));
end test;

The problem is that when u1 == u2 then y = (u1 + u2)/2. Hence whenever one input is constant and the other is increasing then there will be an overshoot starting right after the point where u1 == u2. This effect is stronger whenever the slope of either of the inputs is larger so it seems difficult to find a proper solution for this without making the function inputs complicated.

However some sort of a work around could be to take for instance when u1 == u2 then y = 0.975*(u1 + u2)/2and to rework the rest of the interpolation accordingly. If the user then chooses delta to be sufficiently small then the minimum should be respected and no local maximum should appear. This implementation however poses a problem when the two input signals are nearly equal but large. Therefore this implementation would require some nominal difference value to be defined, which again makes it more complicated.

Additionally the user guide could mention this problem and it could hint towards the use of the min function since it does not generate events anyway.

Note that this problem is for instance quite important in the implementation of PressureIndependentValve where a similar local maximum exists. This causes numerical problems since multiple solutions to the flow problem may exist.

How should we approach with this issue? I think that just fixing the documentation is sufficient.

[mwetter]

I suggest you update the documentation. A have not yet found a reliable way to avoid this undershoot or overshoot. We typically address this in a form such as in http://simulationresearch.lbl.gov/modelica/releases/latest/help/Buildings_BoundaryConditions_WeatherData_BaseClasses.html#Buildings.BoundaryConditions.WeatherData.BaseClasses.CheckRadiation (see parameter HMin.

[mwetter]

This is now on the master

[Mathadon]

I may have found (part of) the solution for this problem.

we could use:

xMax = (x1^6 + x2^6)^(1/6);
xMin = x1 + x2 - xMax ;

E.g. this example model

model test
  parameter Real power=4;
  Real y1=0.1;
  Real y2=(sin(time)+1)*0.1;
  Real yMax = (y1^power+y2^power)^(1/power);
  Real yMin= y1 + y2-yMax;
  Real smoothMin=IBPSA.Utilities.Math.Functions.smoothMin(y1,y2,0.1);
end test;

returns:

we just need to figure out a way to normalise the numbers properly and find a way to support negative numbers.

[Mathadon]

The 'softmax' formulation on wikipedia also overshoots but does work for negative numbers. The log formulation does not work for negative numbers.

[thorade]

Our library has some more smooth functions, they are here: https://github.com/UdK-VPT/BuildingSystems/tree/master/BuildingSystems/Utilities/SmoothFunctions Especially, there are softcut_upper and softcut_lower. These were implemented long time ago in Smile and then ported to Modelica. At some point we want to either remove them from our library or improve them and add them to IBPSA, see https://github.com/UdK-VPT/BuildingSystems/issues/10

[mwetter]

@Mathadon @thorade : Shall this be closed or is there any activity?

[Mathadon]

@mwetter, can be closed if it's a nuisance :)

[thorade]

For comparison, here is the output of the softcut_upper function:

model test
  parameter Real power=4;
  Real y1=0.1;
  Real y2=(sin(time)+1)*0.1;

  Real smoMin1=IBPSA.Utilities.Math.Functions.smoothMin(y1,y2,0.05);
  Real smoMin2=BuildingSystems.Utilities.SmoothFunctions.softcut_upper(y1,y2,0.05);
end test;

I'll do some more testing, and I do not yet understand how our function works, and whether it is C1 or C2... If it works as expected, I would clean it up and add documentation.

[Mathadon]

within BuildingSystems.Utilities.SmoothFunctions;
function softcut_upper "Softly cuts to upper limit"
  input Real x;
  input Real x_ulimit;
  input Real r;
  output Real y;
protected
  Real x_start;
  Real x_end;
  Real help;
  Real SQRT_TWO = 1.4142135624;
algorithm
    x_start := x_ulimit - r * (1.0-1.0/SQRT_TWO);
    x_end := x_start + r/SQRT_TWO;
    if x <= x_start then
      y := x;
    elseif x >= x_end then
      y := x_ulimit;
    else
      help := x_ulimit - x - r * (1.0 - SQRT_TWO);
      y := x_ulimit - r + sqrt(r*r - help*help);
    end if;
end softcut_upper;

Looks quite good! I need to check later whether this function is symmetric, wether it works for negative arguments and whether it is C2-differentiable.

[mwetter]

It looks good if it works for all arguments and if it is C2.

[thorade]

After rewriting the softcut_upper function, I have arrived at

function smoothMin2 "Softly cuts to upper limit"
  input Real x;
  input Real x_ulimit;
  input Real r;
  output Real y;

algorithm 
  y :=     if x <= x_ulimit - r + r/sqrt(2) then x
       elseif x >= x_ulimit - r + r*sqrt(2) then x_ulimit
       else x_ulimit - r + sqrt(r^2 - (x_ulimit - x - r + r*sqrt(2))^2);

  annotation (Inline=true, smoothOrder=1);
end smoothMin2;

To my limited understanding, the values and the first derivative are continuous, so it is as good or bad as the current implementation, except that it respects the specified limit and does not overshoot.

[Mathadon]

I think the current implementation is infinitely differentiable? See #736

[thorade]

But in the code it says smoothOrder=1!? https://github.com/ibpsa/modelica-ibpsa/blob/f58e63e1e281f2886d127f649ce1853b9678bfa2/IBPSA/Utilities/Math/Functions/smoothMin.mo#L14

[Mathadon]

That's a bug ;)

Well, actually, I think it's infinitely differentiable, but I haven't been able to convince Michael yet such that we can change that annotation =)