and even more margin fixing and testing, now correct for systems with

repagh · repagh · commit 8d22279da0d5 · 2017-06-10T14:14:51.000+02:00
differentiating s, and return NAN for phase crossover frequency when
no gain margin found
diff --git a/control/margins.py b/control/margins.py
@@ -87,6 +87,14 @@ def _polysqr(pol):
 # RvP, July 8, 2015, augmented to calculate all phase/gain crossings with
 #                    frd data. Correct to return smallest phase
 #                    margin, smallest gain margin and their frequencies
+# RvP, Jun 10, 2017, modified the inclusion of roots found for phase
+#                    crossing to include all >= 0, made subsequent calc
+#                    insensitive to div by 0
+#                    also changed the selection of which crossings to
+#                    return on basis of "A note on the Gain and Phase
+#                    Margin Concepts" Journal of Control and Systems
+#                    Engineering, Yazdan Bavafi-Toosi, Dec 2015, vol 3
+#                    issue 1, pp 51-59, closer to Matlab behavior??
 def stability_margins(sysdata, returnall=False, epsw=0.0):
     """Calculate stability margins and associated crossover frequencies.
 
@@ -104,7 +112,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
         minimum stability margins.  For frequency data or FRD systems, only one
         margin is found and returned.
     epsw: float, optional
-        Frequencies below this value (default 1e-8) are considered static gain,
+        Frequencies below this value (default 0.0) are considered static gain,
         and not returned as margin.
 
     Returns
@@ -263,17 +271,19 @@ def dstab(w):
     if returnall:
         return GM, PM, SM, w_180, wc, wstab
     else:
-        if GM.shape[0]:
+        if GM.shape[0] and not np.isinf(GM).all():
             with np.errstate(all='ignore'):
                 gmidx = np.where(np.abs(np.log(GM)) == 
                                  np.min(np.abs(np.log(GM))))
+        else:
+            gmidx = -1
         if PM.shape[0]:
             pmidx = np.where(np.abs(PM) == np.amin(np.abs(PM)))[0]
         return (
-            (not GM.shape[0] and float('inf')) or GM[gmidx][0],
+            (not gmidx != -1 and float('inf')) or GM[gmidx][0],
             (not PM.shape[0] and float('inf')) or PM[pmidx][0],
             (not SM.shape[0] and float('inf')) or np.amin(SM),
-            (not w_180.shape[0] and float('nan')) or w_180[gmidx][0],
+            (not gmidx != -1 and float('nan')) or w_180[gmidx][0],
             (not wc.shape[0] and float('nan')) or wc[pmidx][0],
             (not wstab.shape[0] and float('nan')) or wstab[SM==np.amin(SM)][0])
 
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
@@ -16,9 +16,9 @@ def assert_array_almost_equal(x, y, ndigit=4):
     x = np.array(x)
     y = np.array(y)
     try:
-        if np.isfinite(x).any() and (np.isfinite(x) == np.isfinite(y)).all() and \
-            (x[np.logical_not(np.isfinite(x))] ==
-            y[np.logical_not(np.isfinite(y))]).all():
+        if np.isfinite(x).any() and \
+           np.equal(np.isfinite(x), np.isfinite(y)).all() and \
+           np.equal(np.isnan(x), np.isnan(y)).all():
             np.testing.assert_array_almost_equal(
                 x[np.isfinite(x)], y[np.isfinite(y)], ndigit)
             return
@@ -69,15 +69,25 @@ def setUp(self):
           [2.2716, 97.5941, 0.5591, 10.0053, 0.0850, 9.9918]
 
         """
+        hm1 = s/(s+1);
+        h0 = 1/(s+1)^3;
         h1 = (s + 0.1)/s/(s+1);
         h2 = (s + 0.1)/s^2/(s+1);
         h3 = (s + 0.1)*(s+0.1)/s^3/(s+1);
         """
         self.types = {
+            'typem1': s/(s+1),
+            'type0': 1/(s+1)**3,
             'type1': (s + 0.1)/s/(s+1),
             'type2': (s + 0.1)/s**2/(s+1),
             'type3': (s + 0.1)*(s+0.1)/s**3/(s+1) }
-        self.tmargin = ( self.types, 
+        self.tmargin = ( self.types,
+            dict(sys='typem1', K=2.0, digits=3, result=(
+                float('Inf'), -120.0007, float('NaN'), 0.5774)),
+            dict(sys='type0', K = 0.8, digits=3, result=(
+                10.0014, float('inf'), 1.7322, float('nan'))),
+            dict(sys='type0', K = 2.0, digits=2, result=(
+                4.000,  67.6058,  1.7322,   0.7663)),
             dict(sys='type1', K=1.0, digits=4, result=(
                 float('Inf'), 144.9032, float('NaN'), 0.3162)),
             dict(sys='type2', K=1.0, digits=4, result=(
@@ -89,10 +99,10 @@ def setUp(self):
         
         # from "A note on the Gain and Phase Margin Concepts
         # Journal of Control and Systems Engineering, Yazdan Bavafi-Toosi,
-        # Dec 205, vol 3 iss 1, pp 51-59
+        # Dec 2015, vol 3 iss 1, pp 51-59
         #
         # A cornucopia of tricky systems for phase / gain margin
-        # Still have to convert this to tests + fix margin to handle
+        # Still have to convert more to tests + fix margin to handle
         # also these torture cases
         """
         % matlab compatible
@@ -273,11 +283,17 @@ def test_zmore_margin(self):
         print("""
         warning, Matlab gives different values (0 and 0) for gain
         margin of the following system:
-        {:s}
+        {type2!s}
         python-control gives inf
         difficult to argue which is right? Special case or different
         approach?
-        """.format(str(self.types['type2'])))
+
+        edge cases, like
+        {type0!s}
+        which approaches a gain of 1 for w -> 0, are also not identically
+        indicated, Matlab gives phase margin -180, at w = 0. for higher or
+        lower gains, results match
+        """.format(**self.types))
                 
         sdict = self.tmargin[0]
         for test in self.tmargin[1:]:
