Merge pull request #2561 from casella/LimiterHomotopy

Add more options for homotopy-based initialization of limiters

Modelica/Blocks/Continuous.mo

@@ -918,10 +918,6 @@ to compute u by an algebraic equation.
       annotation(Dialog(enable=withFeedForward));
     parameter .Modelica.Blocks.Types.InitPID initType= .Modelica.Blocks.Types.InitPID.DoNotUse_InitialIntegratorState
       "Type of initialization (1: no init, 2: steady state, 3: initial state, 4: initial output)"
-       annotation(Evaluate=true,
-        Dialog(group="Initialization"));
-    parameter Boolean limitsAtInit = true
-      "= false, if limits are ignored during initialization"
       annotation(Evaluate=true, Dialog(group="Initialization"));
     parameter Real xi_start=0
       "Initial or guess value for integrator output (= integrator state)"
@@ -936,8 +932,14 @@ to compute u by an algebraic equation.
     parameter Real y_start=0 "Initial value of output"
       annotation(Dialog(enable=initType == .Modelica.Blocks.Types.InitPID.InitialOutput, group=
             "Initialization"));
+    parameter Modelica.Blocks.Types.LimiterHomotopy homotopyType = Modelica.Blocks.Types.LimiterHomotopy.Linear
+      "Simplified model for homotopy-based initialization"
+      annotation (Evaluate=true, Dialog(group="Initialization"));
     parameter Boolean strict=false "= true, if strict limits with noEvent(..)"
       annotation (Evaluate=true, choices(checkBox=true), Dialog(tab="Advanced"));
+    parameter Boolean limitsAtInit=true
+      "Has no longer an effect and is only kept for backwards compatibility (the implementation uses now the homotopy operator)"
+      annotation (Dialog(tab="Dummy"),Evaluate=true, choices(checkBox=true));
     constant Modelica.SIunits.Time unitTime=1 annotation (HideResult=true);
     Modelica.Blocks.Interfaces.RealInput u_ff if withFeedForward
       "Optional connector of feed-forward input signal"
@@ -984,7 +986,8 @@ to compute u by an algebraic equation.
       uMax=yMax,
       uMin=yMin,
       strict=strict,
-      limitsAtInit=limitsAtInit)
+      limitsAtInit=limitsAtInit,
+      homotopyType=homotopyType)
       annotation (Placement(transformation(extent={{70,-10},{90,10}})));
   protected
     parameter Boolean with_I = controllerType==SimpleController.PI or
@@ -1242,13 +1245,23 @@ to compute u_m by an algebraic equation.
 </p>
 
 <p>
-If parameter <strong>limitAtInit</strong> = <strong>false</strong>, the limits at the
-output of this controller block are removed from the initialization problem which
-leads to a much simpler equation system. After initialization has been
-performed, it is checked via an assert whether the output is in the
-defined limits. For backward compatibility reasons
-<strong>limitAtInit</strong> = <strong>true</strong>. In most cases it is best
-to use <strong>limitAtInit</strong> = <strong>false</strong>.
+When initializing in steady-state, homotopy-based initialization can help the convergence of the solver,
+by using a simplified model a the beginning of the solution process. Different options are available. 
+<ul>
+<li><strong>homotopyType=Linear</strong> (default): the limitations are removed from the simplified model,
+making it linear. Use this if you know that the controller will not be saturated at steady state.</li>
+<li><strong>homotopyType=UpperLimit</strong>: if it is known a priori the controller will be stuck at the upper
+limit yMax, this option assumes y = yMax as a simplified model.</li>
+<li><strong>homotopyType=LowerLimit</strong>: if it is known a priori the controller will be stuck at the lower
+limit yMin, this option assumes y = yMin as a simplified model.</li>
+<li><strong>homotopyType=NoHomotopy</strong>: this option does not apply any simplification and keeps the
+limiter active throughout the homotopy transformation. Use this if it is unknown whether the controller
+is saturated or not at initialization and if the limitations on the output must be enforced throughout 
+the entire homotopy transformation.</li>
+</ul>
+</p>
+<p>
+The parameter <strong>limitAtInit</strong> is obsolete since MSL 3.2.2 and only kept for backwards compatibility.
 </p>
 </html>"));
   end LimPID;

Modelica/Blocks/Nonlinear.mo

@@ -7,21 +7,38 @@ package Nonlinear
       block Limiter "Limit the range of a signal"
         parameter Real uMax(start=1) "Upper limits of input signals";
         parameter Real uMin= -uMax "Lower limits of input signals";
-        parameter Boolean strict=false
-      "= true, if strict limits with noEvent(..)"
+        parameter Boolean strict=false "= true, if strict limits with noEvent(..)"
           annotation (Evaluate=true, choices(checkBox=true), Dialog(tab="Advanced"));
+        parameter Types.LimiterHomotopy homotopyType = Modelica.Blocks.Types.LimiterHomotopy.Linear "Simplified model for homotopy-based initialization"
+          annotation (Evaluate=true, Dialog(group="Initialization"));
         parameter Boolean limitsAtInit=true
-      "Has no longer an effect and is only kept for backwards compatibility (the implementation uses now the homotopy operator)"
+          "Has no longer an effect and is only kept for backwards compatibility (the implementation uses now the homotopy operator)"
           annotation (Dialog(tab="Dummy"),Evaluate=true, choices(checkBox=true));
         extends Interfaces.SISO;
+  protected
+        Real simplifiedExpr "Simplified expression for homotopy-based initialization";
 
       equation
         assert(uMax >= uMin, "Limiter: Limits must be consistent. However, uMax (=" + String(uMax) +
                              ") < uMin (=" + String(uMin) + ")");
+        simplifiedExpr = (if homotopyType == Types.LimiterHomotopy.Linear then u
+                          else if homotopyType == Types.LimiterHomotopy.UpperLimit then uMax
+                          else if homotopyType == Types.LimiterHomotopy.LowerLimit then uMin
+                          else 0);
         if strict then
-           y = homotopy(actual = smooth(0, noEvent(if u > uMax then uMax else if u < uMin then uMin else u)), simplified=u);
+          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
-           y = homotopy(actual = smooth(0,if u > uMax then uMax else if u < uMin then uMin else u), simplified=u);
+          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;
         annotation (
           Documentation(info="<html>
@@ -31,6 +48,20 @@ as long as the input is within the specified upper and lower
 limits. If this is not the case, the corresponding limits are passed
 as output.
 </p>
+<p>
+The parameter <code>homotopyType</code> in the Advanced tab specifies the
+simplified behaviour if homotopy-based initialization is used:
+<ul>
+<li><code>NoHomotopy</code>: the actual expression with limits is used</li>
+<li><code>Linear</code>: a linear behaviour y = u is assumed (default option)</li>
+<li><code>UpperLimit</code>: it is assumed that the output is stuck at the upper limit u = uMax</li>
+<li><code>LowerLimit</code>: it is assumed that the output is stuck at the lower limit u = uMin</li>
+</ul>
+</p>
+<p>
+If it is known a priori in which region the input signal will be located, this option can help 
+a lot by removing one strong nonlinearity from the initialization problem. 
+</p>
 </html>"), Icon(coordinateSystem(
           preserveAspectRatio=true,
           extent={{-100,-100},{100,100}}), graphics={
@@ -99,7 +130,11 @@ as output.
   block VariableLimiter "Limit the range of a signal with variable limits"
     extends Interfaces.SISO;
     parameter Boolean strict=false "= true, if strict limits with noEvent(..)"
-      annotation (Evaluate=true, choices(checkBox=true));
+      annotation (Evaluate=true, choices(checkBox=true), Dialog(tab="Advanced"));
+    parameter Types.VariableLimiterHomotopy homotopyType = Modelica.Blocks.Types.VariableLimiterHomotopy.Linear "Simplified model for homotopy-based initialization"
+      annotation (Evaluate=true, Dialog(group="Initialization"));
+    parameter Real ySimplified = 0 "Fixed value of output in simplified model"
+      annotation (Dialog(tab="Advanced", enable=homotopyType == Modelica.Blocks.Types.VariableLimiterHomotopy.Fixed));
     parameter Boolean limitsAtInit=true
       "Has no longer an effect and is only kept for backwards compatibility (the implementation uses now the homotopy operator)"
       annotation (Dialog(tab="Dummy"),Evaluate=true, choices(checkBox=true));
@@ -109,13 +144,27 @@ as output.
     Interfaces.RealInput limit2
       "Connector of Real input signal used as minimum of input u"
       annotation (Placement(transformation(extent={{-140,-100},{-100,-60}})));
+  protected
+    Real simplifiedExpr "Simplified expression for homotopy-based initialization";
   equation
     assert(limit1 >= limit2, "Input signals are not consistent: limit1 < limit2");
-
+    simplifiedExpr = (if homotopyType == Types.VariableLimiterHomotopy.Linear then u
+                      else if homotopyType == Types.VariableLimiterHomotopy.Fixed then ySimplified
+                      else 0);
     if strict then
-       y = homotopy(actual = smooth(0, noEvent(if u > limit1 then limit1 else if u < limit2 then limit2 else u)), simplified=u);
+      if homotopyType == Types.VariableLimiterHomotopy.NoHomotopy then
+        y = smooth(0, noEvent(if u > limit1 then limit1 else if u < limit2 then limit2 else u));
+      else
+        y = homotopy(actual = smooth(0, noEvent(if u > limit1 then limit1 else if u < limit2 then limit2 else u)),
+                     simplified=simplifiedExpr);
+      end if;
     else
-       y = homotopy(actual = smooth(0,if u > limit1 then limit1 else if u < limit2 then limit2 else u), simplified=u);
+      if homotopyType == Types.VariableLimiterHomotopy.NoHomotopy then
+        y = smooth(0,if u > limit1 then limit1 else if u < limit2 then limit2 else u);
+      else
+        y = homotopy(actual = smooth(0,if u > limit1 then limit1 else if u < limit2 then limit2 else u),
+                     simplified=simplifiedExpr);
+      end if;
     end if;
 
     annotation (
@@ -127,6 +176,19 @@ limits specified by the two additional inputs limit1 and
 limit2. If this is not the case, the corresponding limit
 is passed as output.
 </p>
+<p>
+The parameter <code>homotopyType</code> in the Advanced tab specifies the
+simplified behaviour if homotopy-based initialization is used:
+<ul>
+<li><code>NoHomotopy</code>: the actual expression with limits is used</li>
+<li><code>Linear</code>: a linear behaviour y = u is assumed (default option)</li>
+<li><code>Fixed</code>: it is assumed that the output is fixed at the value <code>ySimplified</code></li>
+</ul>
+</p>
+<p>
+If it is known a priori in which region the input signal will be located, this option can help 
+a lot by removing one strong nonlinearity from the initialization problem. 
+</p>
 </html>"), Icon(coordinateSystem(preserveAspectRatio=true, extent={{-100,-100},{100,
               100}}), graphics={
           Line(points={{0,-90},{0,68}}, color={192,192,192}),

Modelica/Blocks/Types.mo

@@ -15,7 +15,7 @@ package Types
         "Table points are interpolated (by Steffen splines) such that the monotonicity is preserved and the first derivative is continuous")
     "Enumeration defining the smoothness of table interpolation";
 
-  type Extrapolation = enumeration(
+    type Extrapolation = enumeration(
       HoldLastPoint
         "Hold the first/last table point outside of the table scope",
       LastTwoPoints
@@ -24,13 +24,13 @@ package Types
       NoExtrapolation "Extrapolation triggers an error")
     "Enumeration defining the extrapolation of table interpolation";
 
-  type TimeEvents = enumeration(
+    type TimeEvents = enumeration(
       Always "Always generate time events at interval boundaries",
       AtDiscontinuities "Generate time events at discontinuities (defined by duplicated sample points)",
       NoTimeEvents "No time events at interval boundaries")
     "Enumeration defining the time event handling of time table interpolation";
 
-  type Init = enumeration(
+    type Init = enumeration(
       NoInit
         "No initialization (start values are used as guess values with fixed=false)",
       SteadyState
@@ -53,7 +53,7 @@ package Types
   </dl>
 </html>"));
 
-  type InitPID = enumeration(
+    type InitPID = enumeration(
       NoInit
         "No initialization (start values are used as guess values with fixed=false)",
       SteadyState
@@ -91,7 +91,7 @@ initialization definition.
   </dl>
 </html>"));
 
-  type SimpleController = enumeration(
+   type SimpleController = enumeration(
       P "P controller",
       PI "PI controller",
       PD "PD controller",
@@ -121,6 +121,19 @@ initialization definition.
       Cosine "Cosine regularization")
     "Enumeration defining the regularization around zero";
 
+  type LimiterHomotopy = enumeration(
+      NoHomotopy "Homotopy is not used",
+      Linear "Simplified model without limits",
+      UpperLimit "Simplified model fixed at upper limit",
+      LowerLimit "Simplified model fixed at lower limit")
+    "Enumeration defining use of homotopy in limiter components" annotation (Evaluate=true);
+
+  type VariableLimiterHomotopy = enumeration(
+      NoHomotopy "Simplified model = actual model",
+      Linear "Simplified model: y = u",
+      Fixed "Simplified model: y = ySimplified")
+    "Enumeration defining use of homotopy in variable limiter components" annotation (Evaluate=true);
+
   class ExternalCombiTimeTable
     "External object of 1-dim. table where first column is time"
     extends ExternalObject;

Modelica/Electrical/Analog/Examples/OpAmps.mo

@@ -566,7 +566,8 @@ package OpAmps "Examples with operational amplifiers"
     Modelica.Electrical.Analog.Ideal.IdealizedOpAmpLimted opAmp(
       Vps=Vps,
       Vns=Vns,
-      out(i(start=0)))
+      out(i(start=0)),
+      homotopyType=Modelica.Blocks.Types.LimiterHomotopy.LowerLimit)
       annotation (Placement(transformation(extent={{0,-10},{20,10}})));
     Modelica.Electrical.Analog.Basic.Ground ground
       annotation (Placement(transformation(extent={{-20,-100},{0,-80}})));
@@ -639,7 +640,8 @@ package OpAmps "Examples with operational amplifiers"
     Modelica.Electrical.Analog.Ideal.IdealizedOpAmpLimted opAmp(
       Vps=Vps,
       Vns=Vns,
-      out(i(start=0)))
+      out(i(start=0)),
+      homotopyType=Modelica.Blocks.Types.LimiterHomotopy.UpperLimit)
       annotation (Placement(transformation(extent={{0,10},{20,-10}})));
     Modelica.Electrical.Analog.Basic.Ground ground
       annotation (Placement(transformation(extent={{-20,-100},{0,-80}})));
@@ -707,8 +709,9 @@ package OpAmps "Examples with operational amplifiers"
     parameter SI.Resistance R2=1000 "Resistance 2 for adjusting the Schmitt trigger voltage level";
     parameter SI.Resistance R=1000 "Arbitrary resistance";
     parameter SI.Capacitance C=1/f/(2*R*log(1 + 2*R1/R2)) "Calculated capacitance to reach the desired frequency f";
-    Modelica.Electrical.Analog.Ideal.IdealizedOpAmpLimted opAmp(Vps=Vps, Vns=
-          Vns) annotation (Placement(transformation(extent={{0,-10},{20,10}})));
+    Modelica.Electrical.Analog.Ideal.IdealizedOpAmpLimted opAmp(
+      Vps=Vps,
+      Vns=Vns) annotation (Placement(transformation(extent={{0,-10},{20,10}})));
     Modelica.Electrical.Analog.Basic.Ground ground
       annotation (Placement(transformation(extent={{-20,-80},{0,-60}})));
     Modelica.Electrical.Analog.Sensors.VoltageSensor vOut annotation (Placement(
@@ -776,8 +779,11 @@ package OpAmps "Examples with operational amplifiers"
     parameter SI.Frequency f=10 "Desired frequency";
     parameter SI.Resistance R=1000 "Arbitrary resistance of integrator part";
     parameter SI.Capacitance C=Vps/VAmp/(4*f*R) "Calculated capacitance of integrator part to reach f";
-    Modelica.Electrical.Analog.Ideal.IdealizedOpAmpLimted opAmp1(Vps=Vps, Vns=
-          Vns)
+    Modelica.Electrical.Analog.Ideal.IdealizedOpAmpLimted opAmp1(
+      Vps=Vps,
+      Vns=Vns,
+      strict=false,
+      homotopyType=Modelica.Blocks.Types.LimiterHomotopy.UpperLimit)
       annotation (Placement(transformation(extent={{-60,10},{-40,-10}})));
     Modelica.Electrical.Analog.Basic.Resistor r2(R=R2, i(start=Vps/R2))
       annotation (Placement(transformation(
@@ -793,7 +799,8 @@ package OpAmps "Examples with operational amplifiers"
     Modelica.Electrical.Analog.Ideal.IdealizedOpAmpLimted opAmp2(
       Vps=Vps,
       Vns=Vns,
-      v_in(start=0))
+      v_in(start=0),
+      strict=false)
       annotation (Placement(transformation(extent={{30,-10},{50,10}})));
     Modelica.Electrical.Analog.Basic.Capacitor c(C=C, v(fixed=true, start=0))
       annotation (Placement(transformation(extent={{50,10},{30,30}})));

Modelica/Electrical/Analog/Ideal.mo

@@ -675,6 +675,10 @@ If the input voltage is vin larger than 0, the output voltage is out.v = VMax.
       annotation (Dialog(enable=not useSupply));
     parameter SI.Voltage Vns=-15 "Negative supply voltage"
       annotation (Dialog(enable=not useSupply));
+    parameter Boolean strict=true "= true, if strict limits with noEvent(..)"
+      annotation (Evaluate=true, choices(checkBox=true), Dialog(tab="Advanced"));
+    parameter Modelica.Blocks.Types.LimiterHomotopy homotopyType = Modelica.Blocks.Types.LimiterHomotopy.Linear "Simplified model for homotopy-based initialization"
+      annotation (Evaluate=true, Dialog(group="Initialization"));
     SI.Voltage vps "Positive supply voltage";
     SI.Voltage vns "Negative supply voltage";
     SI.Voltage v_in=in_p.v - in_n.v "Input voltage difference";
@@ -700,14 +704,34 @@ If the input voltage is vin larger than 0, the output voltage is out.v = VMax.
     Modelica.Electrical.Analog.Interfaces.NegativePin s_n(final i=-i_s, final v=
          vns) if useSupply "Optional negative supply pin" annotation (Placement(
           transformation(extent={{-10,-110},{10,-90}})));
+  protected
+    SI.Voltage simplifiedExpr "Simplified expression for homotopy-based initialization";
   equation
     if not useSupply then
       vps = Vps;
       vns = Vns;
     end if;
     in_p.i = 0;
     in_n.i = 0;
-    v_out = homotopy(actual = smooth(0, if V0*v_in<vns then vns else if V0*v_in>vps then vps else V0*v_in), simplified=V0*v_in);
+    simplifiedExpr = (if homotopyType == Modelica.Blocks.Types.LimiterHomotopy.Linear then V0*v_in
+                      else if homotopyType == Modelica.Blocks.Types.LimiterHomotopy.UpperLimit then vps
+                      else if homotopyType == Modelica.Blocks.Types.LimiterHomotopy.LowerLimit then vns
+                      else 0);
+    if strict then
+      if homotopyType == Modelica.Blocks.Types.LimiterHomotopy.NoHomotopy then
+        v_out = smooth(0, noEvent(if V0*v_in>vps then vps else if V0*v_in<vns then vns else V0*v_in));
+      else
+        v_out = homotopy(actual = smooth(0, noEvent(if V0*v_in>vps then vps else if V0*v_in<vns then vns else V0*v_in)),
+                         simplified=simplifiedExpr);
+      end if;
+    else
+      if homotopyType == Modelica.Blocks.Types.LimiterHomotopy.NoHomotopy then
+        v_out = smooth(0, if V0*v_in>vps then vps else if V0*v_in<vns then vns else V0*v_in);
+      else
+        v_out = homotopy(actual = smooth(0, if V0*v_in>vps then vps else if V0*v_in<vns then vns else V0*v_in),
+                         simplified=simplifiedExpr);
+      end if;
+    end if;
     annotation (defaultComponentName="opAmp",
       Icon(coordinateSystem(preserveAspectRatio=false, extent={{-100,-100},{100,
               100}}), graphics={