Invalid smooth() statement in SlewRateLimiter

[casella]

Description of the problem

The core equation of the Modelica.Blocks.Nonlinear.SlewRateLimiter block reads:

ModelicaStandardLibrary/Modelica/Blocks/Nonlinear.mo

Lines 206 to 210 in e653d4f

	if strict then
	der(y) = smooth(1, (if noEvent(val<Falling) then Falling else if noEvent(val>Rising) then Rising else val));
	else
	der(y) = if val<Falling then Falling else if val>Rising then Rising else val;
	end if;

This is likely inspired by the code of the Limiter block

ModelicaStandardLibrary/Modelica/Blocks/Nonlinear.mo

Lines 25 to 39 in e653d4f

	if strict then
	if homotopyType == Types.LimiterHomotopy.NoHomotopy then
	y = smooth(0, noEvent(if u > uMax then uMax else if u < uMin then uMin else u));
	else
	y = homotopy(actual = smooth(0, noEvent(if u > uMax then uMax else if u < uMin then uMin else u)),
	simplified=simplifiedExpr);
	end if;
	else
	if homotopyType == Types.LimiterHomotopy.NoHomotopy then
	y = smooth(0,if u > uMax then uMax else if u < uMin then uMin else u);
	else
	y = homotopy(actual = smooth(0,if u > uMax then uMax else if u < uMin then uMin else u),
	simplified=simplifiedExpr);
	end if;
	end if;

In the case of the Limiter block, the strict = true version makes sense when the block output should never exceed the upper and lower limit, e.g. to avoid division-by-zero, sqrt-neg-number or log-neg-number issues. Regular event detection necessarily implies computing out-of-bounds values during the solver iterations to find the root of the zero-crossing functions, so it is necessary to disable event detection to avoid that; luckily, since the noEvent expression is continuous, this can be done without too much damage for the ODE integrator.

The SlewRateLimiter block also has a strict = true version that disables event generation and introduces the smooth() operator. Unfortunately, in this case the expression on the RHS of the differential equation is not continuously differentiable, as it is incorrectly declared by the smooth() operator; in fact, when the input undergoes large enough step changes, that expression is discontinuous, as demonstrated by this simple MWE:

model TestSlewRate
  Modelica.Blocks.Nonlinear.SlewRateLimiter slewRateLimiter(Rising = 300, Td = 0.03);
  Modelica.Blocks.Sources.Step step(height = 1e5, startTime = 100);
equation
  connect(step.y, slewRateLimiter.u);
annotation(
    uses(Modelica(version = "4.1.0")),
  experiment(StartTime = 0, StopTime = 500, Tolerance = 1e-06, Interval = 10));
end TestSlewRate;

As F. Cellier very clearly discusses in Chapter 9.2 of the book Continuous System Modelling, relying on error estimation and step-size adaption algorithms to handle discontinuous variables without explicit event handling is not a good idea and can lead to completely wrong results. On the other hand, if the output of the derivative block has a derivative that slightly exceeds the given limits during the iterations to find the roots of the zc-function generating the event, nothing harmful happens.

Suggested solution

Based on the analyis reported above I would deprecate the strict = true behaviour, by moving the strict parameter to a Deprecated tab and changing the code to

der(y) = if val<Falling then Falling else if val>Rising then Rising else val; 
if strict then 
   assert("strict = true is deprecated and has no effect in this version of the library; the strict parameter will be removed in future versions of the library", AssertionLevel.warning);
end if; 

[casella]

OpenModelica ticket #13341 shows one concrete case where the application of strict = true leads to wrong, tool-dependent results. Taking events into account leads to correct results.

Unfortunately I was not able to replicate that behaviour on a simpler MWE.

[HansOlsson]

As a minimum it should be smooth(0, as it is continuous, but not continuously differentiable, when the inputs change continuously.
However, it is odd that smooth is only applied to the strict-part.

[HansOlsson]

Thinking more:

• I agree that the "strict" noEvent doesn't really make sense here, and even if smooth(0, might be possible, it seems best without it.
• The MWE does not show any issue. Smooth is only a local guarantee that the expression itself is smooth based on the variables used in it, not that the variables used in the expression actually vary smoothly, or as stated in the specification: "Note that smooth does not guarantee a smooth output if any of the occurring variables change discontinuously."
• I tried to track the reason for the problem in this model:
  • Initial condition of SlewRateLimiter #2140 changed the initialization, but no real change for this equation.
  • SlewRateLimiter parameter strict has no effect #1798 made it slightly better, as it previously used min/max (which doesn't generate events) in both branches! In retrospect I should have noticed the other errors as well.
  • Limiter (and LimPID) does not treat limits literally #987 added strict
  • Add new RateLimiter block to MSL #529 and the corresponding commit 1a53fdd added the original version, but doesn't explain those design choices.

[henrikt-ma]

• I agree that the "strict" noEvent doesn't really make sense here, and even if smooth(0, might be possible, it seems best without it.

This looks like precisely the kind of expression where smooth(0, …) is often used. I don't see any harm in expressing the continuity guarantee, so why would it be best without it?

[henrikt-ma]

Deprecated tab and changing the code to

if strict then
assert("strict = true is deprecated and has no effect in this version of the library; the strict parameter will be removed in future versions of the library", AssertionLevel.warning);
end if;

Just noting that this is probably what you mean:

initial equation
  assert(not strict, "Strict mode is deprecated and has no effect in this version of the library; the strict parameter will be removed in future versions of the library.", AssertionLevel.warning);

Note that the assertion string should be formulated under the assumption that the tool will automatically present the violated condition, see https://github.com/modelica/ModelicaSpecification/pull/3641/changes.

[henrikt-ma]

• The MWE does not show any issue. Smooth is only a local guarantee that the expression itself is smooth based on the variables used in it, not that the variables used in the expression actually vary smoothly, or as stated in the specification: "Note that smooth does not guarantee a smooth output if any of the occurring variables change discontinuously."

Just so that we are on the same page, the following is an MWE, right?

model TestSlewRate
  Modelica.Blocks.Sources.Sine u(amplitude = 1.0, f = 1.0);
  Modelica.Blocks.Nonlinear.SlewRateLimiter slewRateLimiter(Rising = 1, Falling = -1, Td = 1);
  Modelica.Blocks.Continuous.Der d1;
  Modelica.Blocks.Continuous.Der d2;
equation
  connect(u.y, slewRateLimiter.u);
  connect(slewRateLimiter.y, d1.u);
  connect(d1.y, d2.u);
  annotation(
    Documentation(
      figures = {
        Figure(
          title = "Smoothness check",
          identifier = "busy-thales",
          preferred = true,
          plots = {
            Plot(
              identifier = "exprs",
              curves = {
                Curve(y = der(slewRateLimiter.val))
              }
            ),
            Plot(
              identifier = "smooth-der",
              curves = {
                Curve(y = d2.y)
              }
            )
          },
          caption = "%(plot:exprs) In the slewRateLimiter, the only Real expression inside smooth(1, …) with non-zero derivative of time is %(variable:slewRateLimiter.val). It is seen that this derivative is continuous.

%(plot:smooth-der) Since the Real expressions inside smooth(1, …) have continuous derivatives of time, the smoothness guarantee states that %(variable:der(slewRateLimiter.y)) shall have continuous derivative of time. That is, %(variable:slewRateLimiter.y) should be two times continuously differentiable with respect to time, and %(variable:d2.y) is this second order derivative. It is seen that the smoothness guarantee is violated.
"
        )
      }
    ),
    uses(
      Modelica(
        version = "4.1.0"
      )
    ),
    experiment(
      StopTime = 3.0
    )
  );
end TestSlewRate;

[casella]

• The MWE does not show any issue. Smooth is only a local guarantee that the expression itself is smooth based on the variables used in it, not that the variables used in the expression actually vary smoothly, or as stated in the specification: "Note that smooth does not guarantee a smooth output if any of the occurring variables change discontinuously."

@HansOlsson thanks for pointing this out. Indeed, smooth(0, ...) could be kept, as it reflects an objective property of the conditional expression. What is really bad is noEvent, since its argument could change discontinuously, as shown by the MWE.

The MWE does not show any issue

Yes. Unfortunately I was not able to reproduce the issue that plagued #13341, i.e. the output of a filtered step not having the maximum allowed slope throughout the entire transient. The point of the MWE is just to show that the noEvent expression can be brutally discontinuous, so it's not a good idea to suppress events there, in general.

[casella]

Just so that we are on the same page, the following is an MWE, right?

I'd say so. I would probably set Td = 0.1, so that the output that you get is actually a slew-limited sine. Td = 1 is too slow for a 6.28 rad/s sinusoid. But that doesn't really matter w.r.t. the continuity issue: the expression inside smooth is continuous in val, but it's not differentiable when val hits the Rising and Falling limits. In other words, the expression is continuous but it has corner points where it is not differentiable.