Repository: python-control/python-control

Issue:
stability_margins Fails for Many FRD Systems at Stab Margins
There is a fairly obvious copy-and-paste error in the stab margins section for FRD systems in stability_margins. The calculation for the wstab array checks the indexes only for `i+1 < len(sys.omega)`, copied from the phase margin test, and then goes on to access `sys.omega[i-1:i+2]`. Clearly this should be checking both that `i+2 < len(sys.omega)` and that `i-1 >= 0`.

So to correct this, line 244 in margins.py becomes something like:
```
              for i in widx if i+2 < len(sys.omega) and i-1 >= 0 and
```


Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/31746921/badge)](https://coveralls.io/builds/31746921)

Coverage increased (+0.01%) to 84.209% when pulling **13cabb7ba9b27a082c0dfe93d9677a73d8a5ebed on bnavigator:fix-407** into **a9d869e0215b3dfc630ebf67ca1cce8b988a2a53 on python-control:master**.

bnavigator: First commit is just the adjustment of the check within list the construction, but then I could not let it go and cleaned up the whole code block. Now see if I can come up with a meaningful unit test...

Patch:
diff --git a/.travis.yml b/.travis.yml
index d0201031..022e48c6 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -19,7 +19,7 @@ python:
 
 # Test against multiple version of SciPy, with and without slycot
 #
-# Because there were significant changes in SciPy between v0 and v1, we 
+# Because there were significant changes in SciPy between v0 and v1, we
 # test against both of these using the Travis CI environment capability
 #
 # We also want to test with and without slycot
@@ -84,7 +84,6 @@ before_install:
       sudo apt-get update -qq;
       sudo apt-get install liblapack-dev libblas-dev;
       sudo apt-get install gfortran;
-      sudo apt-get install cmake;
     fi
   # use miniconda to install numpy/scipy, to avoid lengthy build from source
   - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then
@@ -104,7 +103,7 @@ before_install:
   # Install scikit-build for the build process if slycot is being used
   - if [[ "$SLYCOT" = "source" ]]; then
       conda install openblas;
-      conda install -c conda-forge scikit-build;
+      conda install -c conda-forge cmake scikit-build;
     fi
   # Make sure to look in the right place for python libraries (for slycot)
   - export LIBRARY_PATH="$HOME/miniconda/envs/test-environment/lib"
diff --git a/control/frdata.py b/control/frdata.py
index 14705947..c57cf09b 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -400,17 +400,30 @@ second has %i." % (self.outputs, other.outputs))
 
     # Method for generating the frequency response of the system
     def freqresp(self, omega):
-        """Evaluate a transfer function at a list of angular frequencies.
-
-        mag, phase, omega = self.freqresp(omega)
-
-        reports the value of the magnitude, phase, and angular frequency of
-        the transfer function matrix evaluated at s = i * omega, where omega
-        is a list of angular frequencies, and is a sorted version of the input
-        omega.
-
+        """Evaluate the frequency response at a list of angular frequencies.
+
+        Reports the value of the magnitude, phase, and angular frequency of
+        the requency response evaluated at omega, where omega is a list of
+        angular frequencies, and is a sorted version of the input omega.
+
+        Parameters
+        ----------
+        omega : array_like
+            A list of frequencies in radians/sec at which the system should be
+            evaluated. The list can be either a python list or a numpy array
+            and will be sorted before evaluation.
+
+        Returns
+        -------
+        mag : (self.outputs, self.inputs, len(omega)) ndarray
+            The magnitude (absolute value, not dB or log10) of the system
+            frequency response.
+        phase : (self.outputs, self.inputs, len(omega)) ndarray
+            The wrapped phase in radians of the system frequency response.
+        omega : ndarray or list or tuple
+            The list of sorted frequencies at which the response was
+            evaluated.
         """
-
         # Preallocate outputs.
         numfreq = len(omega)
         mag = empty((self.outputs, self.inputs, numfreq))
diff --git a/control/iosys.py b/control/iosys.py
index 1b29b5b0..1fe10346 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -1645,8 +1645,10 @@ def find_eqpt(sys, x0, u0=[], y0=None, t=0, params={},
         # and were processed above.
 
         # Get the states and inputs that were not listed as fixed
-        state_vars = np.delete(np.array(range(nstates)), ix)
-        input_vars = np.delete(np.array(range(ninputs)), iu)
+        state_vars = (range(nstates) if not len(ix)
+                      else np.delete(np.array(range(nstates)), ix))
+        input_vars = (range(ninputs) if not len(iu)
+                      else np.delete(np.array(range(ninputs)), iu))
 
         # Set the outputs and derivs that will serve as constraints
         output_vars = np.array(iy)
diff --git a/control/lti.py b/control/lti.py
index c9a58f9c..8db14794 100644
--- a/control/lti.py
+++ b/control/lti.py
@@ -55,8 +55,8 @@ class LTI:
         Parameters
         ----------
         strict: bool, optional
-            If strict is True, make sure that timebase is not None.  Default 
-            is False. 
+            If strict is True, make sure that timebase is not None.  Default
+            is False.
         """
 
         # If no timebase is given, answer depends on strict flag
@@ -75,8 +75,8 @@ class LTI:
         sys : LTI system
             System to be checked
         strict: bool, optional
-            If strict is True, make sure that timebase is not None.  Default 
-            is False. 
+            If strict is True, make sure that timebase is not None.  Default
+            is False.
         """
         # If no timebase is given, answer depends on strict flag
         if self.dt is None:
@@ -421,6 +421,7 @@ def evalfr(sys, x):
         return sys.horner(x)[0][0]
     return sys.horner(x)
 
+
 def freqresp(sys, omega):
     """
     Frequency response of an LTI system at multiple angular frequencies.
@@ -430,13 +431,20 @@ def freqresp(sys, omega):
     sys: StateSpace or TransferFunction
         Linear system
     omega: array_like
-        List of frequencies
+        A list of frequencies in radians/sec at which the system should be
+        evaluated. The list can be either a python list or a numpy array
+        and will be sorted before evaluation.
 
     Returns
     -------
-    mag: ndarray
-    phase: ndarray
-    omega: list, tuple, or ndarray
+    mag : (self.outputs, self.inputs, len(omega)) ndarray
+        The magnitude (absolute value, not dB or log10) of the system
+        frequency response.
+    phase : (self.outputs, self.inputs, len(omega)) ndarray
+        The wrapped phase in radians of the system frequency response.
+    omega : ndarray or list or tuple
+        The list of sorted frequencies at which the response was
+        evaluated.
 
     See Also
     --------
@@ -472,9 +480,9 @@ def freqresp(sys, omega):
         #>>> # frequency response from the 1st input to the 2nd output, for
         #>>> # s = 0.1i, i, 10i.
     """
-
     return sys.freqresp(omega)
 
+
 def dcgain(sys):
     """Return the zero-frequency (or DC) gain of the given system
 
diff --git a/control/margins.py b/control/margins.py
index 193b6c59..7bdcf6ca 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -1,12 +1,12 @@
-"""margin.py
+"""margins.py
 
 Functions for computing stability margins and related functions.
 
 Routines in this module:
 
-margin.stability_margins
-margin.phase_crossover_frequencies
-margin.margin
+margins.stability_margins
+margins.phase_crossover_frequencies
+margins.margin
 """
 
 # Python 3 compatibility (needs to go here)
@@ -211,41 +211,40 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
 
     else:
         # a bit coarse, have the interpolated frd evaluated again
-        def mod(w):
-            """to give the function to calculate |G(jw)| = 1"""
+        def _mod(w):
+            """Calculate |G(jw)| - 1"""
             return np.abs(sys._evalfr(w)[0][0]) - 1
 
-        def arg(w):
-            """function to calculate the phase angle at -180 deg"""
+        def _arg(w):
+            """Calculate the phase angle at -180 deg"""
             return np.angle(-sys._evalfr(w)[0][0])
 
-        def dstab(w):
-            """function to calculate the distance from -1 point"""
+        def _dstab(w):
+            """Calculate the distance from -1 point"""
             return np.abs(sys._evalfr(w)[0][0] + 1.)
 
         # Find all crossings, note that this depends on omega having
         # a correct range
-        widx = np.where(np.diff(np.sign(mod(sys.omega))))[0]
+        widx = np.where(np.diff(np.sign(_mod(sys.omega))))[0]
         wc = np.array(
-            [ sp.optimize.brentq(mod, sys.omega[i], sys.omega[i+1])
-              for i in widx if i+1 < len(sys.omega)])
+            [sp.optimize.brentq(_mod, sys.omega[i], sys.omega[i+1])
+             for i in widx])
 
         # find the phase crossings ang(H(jw) == -180
-        widx = np.where(np.diff(np.sign(arg(sys.omega))))[0]
+        widx = np.where(np.diff(np.sign(_arg(sys.omega))))[0]
         widx = widx[np.real(sys._evalfr(sys.omega[widx])[0][0]) <= 0]
         w_180 = np.array(
-            [ sp.optimize.brentq(arg, sys.omega[i], sys.omega[i+1])
-              for i in widx if i+1 < len(sys.omega) ])
+            [sp.optimize.brentq(_arg, sys.omega[i], sys.omega[i+1])
+             for i in widx])
 
         # find all stab margins?
-        widx = np.where(np.diff(np.sign(np.diff(dstab(sys.omega)))))[0]
-        wstab = np.array([ sp.optimize.minimize_scalar(
-                  dstab, bracket=(sys.omega[i], sys.omega[i+1])).x
-              for i in widx if i+1 < len(sys.omega) and
-              np.diff(np.diff(dstab(sys.omega[i-1:i+2])))[0] > 0 ])
-        wstab = wstab[(wstab >= sys.omega[0]) *
-                      (wstab <= sys.omega[-1])]
-
+        widx, = np.where(np.diff(np.sign(np.diff(_dstab(sys.omega)))) > 0)
+        wstab = np.array(
+            [sp.optimize.minimize_scalar(_dstab,
+                                         bracket=(sys.omega[i], sys.omega[i+1])
+                                         ).x
+             for i in widx])
+        wstab = wstab[(wstab >= sys.omega[0]) * (wstab <= sys.omega[-1])]
 
     # margins, as iterables, converted frdata and xferfcn calculations to
     # vector for this
@@ -254,13 +253,13 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
         GM = 1.0/gain_w_180
     SM = np.abs(sys._evalfr(wstab)[0][0]+1)
     PM = np.remainder(np.angle(sys._evalfr(wc)[0][0], deg=True), 360.0) - 180.0
-    
+
     if returnall:
         return GM, PM, SM, w_180, wc, wstab
     else:
         if GM.shape[0] and not np.isinf(GM).all():
             with np.errstate(all='ignore'):
-                gmidx = np.where(np.abs(np.log(GM)) == 
+                gmidx = np.where(np.abs(np.log(GM)) ==
                                  np.min(np.abs(np.log(GM))))
         else:
             gmidx = -1
@@ -276,7 +275,6 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
 
 
 # Contributed by Steffen Waldherr <waldherr@ist.uni-stuttgart.de>
-#! TODO - need to add test functions
 def phase_crossover_frequencies(sys):
     """Compute frequencies and gains at intersections with real axis
     in Nyquist plot.
diff --git a/control/statesp.py b/control/statesp.py
index 1779dfbf..d176e98c 100644
--- a/control/statesp.py
+++ b/control/statesp.py
@@ -462,11 +462,8 @@ but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs))
                                     self.B)) + self.D
         return array(resp)
 
-    # Method for generating the frequency response of the system
     def freqresp(self, omega):
-        """Evaluate the system's transfer func. at a list of freqs, omega.
-
-        mag, phase, omega = self.freqresp(omega)
+        """Evaluate the system's transfer function at a list of frequencies
 
         Reports the frequency response of the system,
 
@@ -479,26 +476,22 @@ but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs))
 
         Parameters
         ----------
-        omega : array
+        omega : array_like
             A list of frequencies in radians/sec at which the system should be
             evaluated. The list can be either a python list or a numpy array
             and will be sorted before evaluation.
 
         Returns
         -------
-        mag : float
+        mag : (self.outputs, self.inputs, len(omega)) ndarray
             The magnitude (absolute value, not dB or log10) of the system
             frequency response.
-
-        phase : float
+        phase : (self.outputs, self.inputs, len(omega)) ndarray
             The wrapped phase in radians of the system frequency response.
-
-        omega : array
+        omega : ndarray
             The list of sorted frequencies at which the response was
             evaluated.
-
         """
-
         # In case omega is passed in as a list, rather than a proper array.
         omega = np.asarray(omega)
 
diff --git a/control/xferfcn.py b/control/xferfcn.py
index cb351de0..96f3b5ca 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -639,19 +639,36 @@ class TransferFunction(LTI):
 
         return out
 
-    # Method for generating the frequency response of the system
     def freqresp(self, omega):
-        """Evaluate a transfer function at a list of angular frequencies.
+        """Evaluate the transfer function at a list of angular frequencies.
 
-        mag, phase, omega = self.freqresp(omega)
+        Reports the frequency response of the system,
 
-        reports the value of the magnitude, phase, and angular frequency of
-        the transfer function matrix evaluated at s = i * omega, where omega
-        is a list of angular frequencies, and is a sorted version of the input
-        omega.
+             G(j*omega) = mag*exp(j*phase)
 
-        """
+        for continuous time. For discrete time systems, the response is
+        evaluated around the unit circle such that
+
+             G(exp(j*omega*dt)) = mag*exp(j*phase).
+
+        Parameters
+        ----------
+        omega : array_like
+            A list of frequencies in radians/sec at which the system should be
+            evaluated. The list can be either a python list or a numpy array
+            and will be sorted before evaluation.
 
+        Returns
+        -------
+        mag : (self.outputs, self.inputs, len(omega)) ndarray
+            The magnitude (absolute value, not dB or log10) of the system
+            frequency response.
+        phase : (self.outputs, self.inputs, len(omega)) ndarray
+            The wrapped phase in radians of the system frequency response.
+        omega : ndarray or list or tuple
+            The list of sorted frequencies at which the response was
+            evaluated.
+        """
         # Preallocate outputs.
         numfreq = len(omega)
         mag = empty((self.outputs, self.inputs, numfreq))
diff --git a/examples/robust_mimo.py b/examples/robust_mimo.py
index 402d9148..d4e1335e 100644
--- a/examples/robust_mimo.py
+++ b/examples/robust_mimo.py
@@ -44,7 +44,7 @@ def triv_sigma(g, w):
     w - frequencies, length m
     s - (m,n) array of singular values of g(1j*w)"""
     m, p, _ = g.freqresp(w)
-    sjw = (m*np.exp(1j*p*np.pi/180)).transpose(2, 0, 1)
+    sjw = (m*np.exp(1j*p)).transpose(2, 0, 1)
     sv = np.linalg.svd(sjw, compute_uv=False)
     return sv
 


Test patch:
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
index 7e803a9e..9d59a197 100644
--- a/control/tests/freqresp_test.py
+++ b/control/tests/freqresp_test.py
@@ -7,14 +7,16 @@
 # including bode plots.
 
 import unittest
+import matplotlib.pyplot as plt
 import numpy as np
+from numpy.testing import assert_array_almost_equal
+
 import control as ctrl
 from control.statesp import StateSpace
 from control.xferfcn import TransferFunction
 from control.matlab import ss, tf, bode, rss
 from control.exception import slycot_check
-from control.tests.margin_test import assert_array_almost_equal
-import matplotlib.pyplot as plt
+
 
 class TestFreqresp(unittest.TestCase):
    def setUp(self):
@@ -51,7 +53,7 @@ class TestFreqresp(unittest.TestCase):
       for ax in plt.gcf().axes:
          # Make sure there are 2 lines in each subplot
          assert len(ax.get_lines()) == 2
-      
+
       # Generate two plots as a list; should be on the same axes
       plt.figure(2); plt.clf();
       ctrl.bode_plot([ctrl.tf([1], [1,2,1]), ctrl.tf([5], [1, 1])])
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
old mode 100644
new mode 100755
index 85404b44..2f60c7bc
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -1,315 +1,312 @@
-#!/usr/bin/env python
-#
-# margin_test.py - test suite for stability margin commands
-# RMM, 15 Jul 2011
+#!/usr/bin/env pytest
+"""
+margin_test.py - test suite for stability margin commands
 
+RMM, 15 Jul 2011
+BG, 30 Juin 2020 -- convert to pytest, gh-425
+"""
 from __future__ import print_function
-import unittest
+
 import numpy as np
-from control.xferfcn import TransferFunction
-from control.frdata import FRD
+from numpy import inf, nan
+from numpy.testing import assert_allclose
+import pytest
+
+from control.frdata import FrequencyResponseData
+from control.margins import margin, phase_crossover_frequencies, \
+                            stability_margins
 from control.statesp import StateSpace
-from control.margins import *
-
-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.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
-    except TypeError as e:
-        print("Error", e, "with", x, "and", y)
-        #raise e
-    np.testing.assert_array_almost_equal(x, y, ndigit)
-
-class TestMargin(unittest.TestCase):
-    """These are tests for the margin commands in margin.py."""
-
-    def setUp(self):
-        # system, gain margin, gm freq, phase margin, pm freq
-        s = TransferFunction([1, 0], [1])
-        self.tsys = (
-        (TransferFunction([1, 2], [1, 2, 3]),
-         [], [], [], []),
+from control.xferfcn import TransferFunction
+
+
+s = TransferFunction([1, 0], [1])
+
+# (system, stability_margins(sys), stability_margins(sys, returnall=True))
+tsys = [(TransferFunction([1, 2], [1, 2, 3]),
+         (inf, inf, inf, nan, nan, nan),
+         ([], [], [], [], [], [])),
         (TransferFunction([1], [1, 2, 3, 4]),
-        [2.001], [1.7321], [], []),
+         (2., inf, 0.4170, 1.7321, nan, 1.6620),
+         ([2.],     [], [1.2500, 0.4170], [1.7321], [], [0.1690, 1.6620])),
         (StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]],
-            [[1., 0.]], [[0.]]),
-        [], [], [147.0743], [2.5483]),
-        ((8.75*(4*s**2+0.4*s+1))/((100*s+1)*(s**2+0.22*s+1)) * 
-         1./(s**2/(10.**2)+2*0.04*s/10.+1), 
-        [2.2716], [10.0053], [97.5941, -157.7904, 134.7359],
-        [0.0850, 0.9373, 1.0919]))
-        
-
-        """
-        sys1 = tf([1, 2], [1, 2, 3]);
-        sys2 = tf([1], [1, 2, 3, 4]);
-        sys3 = ss([1, 4; 3, 2], [1; -4], ...
-            [1, 0], [0])
-        s = tf('s')
-        sys4 = (8.75*(4*s^2+0.4*s+1))/((100*s+1)*(s^2+0.22*s+1)) * ...
-                                      1.0/(s^2/(10.0^2)+2*0.04*s/10.0+1);
-        """
-
-        self.sys1 = TransferFunction([1, 2], [1, 2, 3])
-        # alternative
-        # sys1 = tf([1, 2], [1, 2, 3])
-        self.sys2 = TransferFunction([1], [1, 2, 3, 4])
-        self.sys3 = StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]],
-            [[1., 0.]], [[0.]])
-        s = TransferFunction([1, 0], [1])
-        self.sys4 = (8.75*(4*s**2+0.4*s+1))/((100*s+1)*(s**2+0.22*s+1)) * \
-                                      1./(s**2/(10.**2)+2*0.04*s/10.+1)
-        self.stability_margins4 = \
-          [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,
-            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=(
-                float('Inf'), 44.4594, float('NaN'), 0.7907)),
-            dict(sys='type3', K=1.0, digits=3, result=(
-                0.0626, 37.1748, 0.1119, 0.7951)),
-            )
-        
-        
-        # from "A note on the Gain and Phase Margin Concepts
-        # Journal of Control and Systems Engineering, Yazdan Bavafi-Toosi,
-        # Dec 2015, vol 3 iss 1, pp 51-59
-        #
-        # A cornucopia of tricky systems for phase / gain margin
-        # Still have to convert more to tests + fix margin to handle
-        # also these torture cases
-        """
-        % matlab compatible
-        s = tf('s');
-        h21 = 0.002*(s+0.02)*(s+0.05)*(s+5)*(s+10)/( ...
-                (s-0.0005)*(s+0.0001)*(s+0.01)*(s+0.2)*(s+1)*(s+100)^2 );
-        h23 = ((s+0.1)^2 + 1)*(s-0.1)/( ...
-                ((s+0.1)^2+4)*(s+1) );
-        h25a = s/(s^2+2*s+2)^4; h25b = h25a*100;
-        h26a = ((s-0.1)^2 + 1)/( ...
-                (s + 0.1)*((s-0.2)^2 + 4) ) ;
-        h26b = ((s-0.1)^2 + 1)/( ...
-                (s - 0.3)*((s-0.2)^2 + 4) );
-        """
-        self.yazdan = {
-            'example21' :
-            0.002*(s+0.02)*(s+0.05)*(s+5)*(s+10)/(
-                (s-0.0005)*(s+0.0001)*(s+0.01)*(s+0.2)*(s+1)*(s+100)**2 ),
-            
-            'example23' :
-            ((s+0.1)**2 + 1)*(s-0.1)/(
-                ((s+0.1)**2+4)*(s+1) ),
-            
-            'example25a' :
-            s/(s**2+2*s+2)**4,
-            
-            'example26a' :
-            ((s-0.1)**2 + 1)/(
-                (s + 0.1)*((s-0.2)**2 + 4) ),
-            
-            'example26b': ((s-0.1)**2 + 1)/(
-                (s - 0.3)*((s-0.2)**2 + 4) )
-        }
-        self.yazdan['example24'] = self.yazdan['example21']*20000
-        self.yazdan['example25b'] = self.yazdan['example25a']*100
-        self.yazdan['example22'] = self.yazdan['example21']*(s**2 - 2*s + 401)
-        self.ymargin = (
-            dict(sys='example21', K=1.0, digits=2, result=(
-                0.0100,  -14.5640,         0,    0.0022)),
-            dict(sys='example21', K=1000.0, digits=2, result=(
-                0.1793, 22.5215, 0.0243, 0.0630)),
-            dict(sys='example21', K=5000.0, digits=4, result=(
-                4.5596, 21.2101, 0.4385, 0.1868)),
-                    )
-
-        self.yallmargin = (
-            dict(sys='example21', K=1.0, result=(
-                [0.01, 179.2931, 2.2798e+4, 1.5946e+07, 7.2477e+08],
-                [0, 0.0243, 0.4385, 6.8640, 84.9323],
-                [-14.5640],
-                [0.0022]))
-        )
-            
-            
-    def test_stability_margins(self):
-        omega = np.logspace(-2, 2, 2000)
-        for sys,rgm,rwgm,rpm,rwpm in self.tsys:
-            print(sys)
-            out = np.array(stability_margins(sys))
-            gm, pm, sm, wg, wp, ws = out
-            outf = np.array(stability_margins(FRD(sys, omega)))
-            print(out,'\n', outf)
-            #print(out != np.array(None))
-            assert_array_almost_equal(
-                out, outf, 2)
-        # final one with fixed values
-        assert_array_almost_equal(
-            [gm, pm, sm, wg, wp, ws],
-            self.stability_margins4, 3)
-
-    def test_margin(self):
-        gm, pm, wg, wp = margin(self.sys4)
-        assert_array_almost_equal(
-            [gm, pm, wg, wp],
-            self.stability_margins4[:2] + self.stability_margins4[3:5], 3)
-
-        
-    def test_stability_margins_all(self):
-        for sys,rgm,rwgm,rpm,rwpm in self.tsys:
-            out = stability_margins(sys, returnall=True)
-            gm, pm, sm, wg, wp, ws = out
-            print(sys)
-            for res,comp in zip(out, (rgm,rpm,[],rwgm,rwpm,[])):
-                if comp:
-                    print(res, '\n', comp)
-                    assert_array_almost_equal(
-                        res, comp, 2)
-
-    def test_phase_crossover_frequencies(self):
-        omega, gain = phase_crossover_frequencies(self.sys2)
-        assert_array_almost_equal(omega, [1.73205,  0.])
-        assert_array_almost_equal(gain, [-0.5,  0.25])
-
-        tf = TransferFunction([1],[1,1])
-        omega, gain = phase_crossover_frequencies(tf)
-        assert_array_almost_equal(omega, [0.])
-        assert_array_almost_equal(gain, [1.])
-
-        # testing MIMO, only (0,0) element is considered
-        tf = TransferFunction([[[1],[2]],[[3],[4]]],
-                              [[[1, 2, 3, 4],[1,1]],[[1,1],[1,1]]])
-        omega, gain = phase_crossover_frequencies(tf)
-        assert_array_almost_equal(omega, [1.73205081,  0.])
-        assert_array_almost_equal(gain, [-0.5,  0.25])
-
-    def test_mag_phase_omega(self):
-        # test for bug reported in gh-58
-        sys = TransferFunction(15, [1, 6, 11, 6])
-        out = stability_margins(sys)
-        omega = np.logspace(-2,2,1000)
-        mag, phase, omega = sys.freqresp(omega)
-        #print( mag, phase, omega)
-        out2 = stability_margins((mag, phase*180/np.pi, omega))
-        ind = [0,1,3,4]   # indices of gm, pm, wg, wp -- ignore sm
-        marg1 = np.array(out)[ind]
-        marg2 = np.array(out2)[ind]
-        assert_array_almost_equal(marg1, marg2, 4)
-
-    def test_frd(self):
-        f = np.array([0.005, 0.010, 0.020, 0.030, 0.040,
-              0.050, 0.060, 0.070, 0.080, 0.090,
-              0.100, 0.200, 0.300, 0.400, 0.500,
-              0.750, 1.000, 1.250, 1.500, 1.750,
-              2.000, 2.250, 2.500, 2.750, 3.000,
-              3.250, 3.500, 3.750, 4.000, 4.250,
-              4.500, 4.750, 5.000, 6.000, 7.000,
-              8.000, 9.000, 10.000 ])
-        gain = np.array([  0.0,   0.0,   0.0,   0.0,   0.0,
-                        0.0,   0.0,   0.0,   0.0,   0.0,
-                   0.0,   0.1,   0.2,   0.3,   0.5,
-                   0.5,  -0.4,  -2.3,  -4.8,  -7.3,
-                  -9.6, -11.7, -13.6, -15.3, -16.9,
-                 -18.3, -19.6, -20.8, -22.0, -23.1,
-                 -24.1, -25.0, -25.9, -29.1, -31.9,
-                 -34.2, -36.2, -38.1 ])
-        phase = np.array([    0,    -1,    -2,    -3,    -4,
-                     -5,    -6,    -7,    -8,    -9,
-                    -10,   -19,   -29,   -40,   -51,
-                    -81,  -114,  -144,  -168,  -187,
-                   -202,  -214,  -224,  -233,  -240,
-                   -247,  -253,  -259,  -264,  -269,
-                   -273,  -277,  -280,  -292,  -301,
-                   -307,  -313,  -317 ])
-        # calculate response as complex number
-        resp = 10**(gain / 20) * np.exp(1j * phase / (180./np.pi))
-        # frequency response data
-        fresp = FRD(resp, f*2*np.pi, smooth=True)
-        s=TransferFunction([1,0],[1])
-        G=1./(s**2)
-        K=1.
-        C=K*(1+1.9*s)
-        TFopen=fresp*C*G
-        gm, pm, sm, wg, wp, ws = stability_margins(TFopen)
-        assert_array_almost_equal(
-            [pm], [44.55], 2)
-
-    def test_nocross(self):
-        # what happens when no gain/phase crossover?
-        s = TransferFunction([1, 0], [1])
-        h1 = 1/(1+s)
-        h2 = 3*(10+s)/(2+s)
-        h3 = 0.01*(10-s)/(2+s)/(1+s)
-        gm, pm, wm, wg, wp, ws = stability_margins(h1)
-        assert_array_almost_equal(
-            [gm, pm, wg, wp],
-            [float('Inf'), float('Inf'), float('NaN'), float('NaN')]) 
-        gm, pm, wm, wg, wp, ws = stability_margins(h2)
-        self.assertEqual(pm, float('Inf'))
-        gm, pm, wm, wg, wp, ws = stability_margins(h3)
-        self.assertTrue(np.isnan(wp))
-        omega = np.logspace(-2,2, 100)
-        out1b = stability_margins(FRD(h1, omega))
-        out2b = stability_margins(FRD(h2, omega))
-        out3b = stability_margins(FRD(h3, omega))
-
-    def test_zmore_margin(self):
-        print("""
-        warning, Matlab gives different values (0 and 0) for gain
-        margin of the following system:
-        {type2!s}
-        python-control gives inf
-        difficult to argue which is right? Special case or different
-        approach?
-
-        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:]:
-            res = margin(sdict[test['sys']]*test['K'])
-            print("more margin {}\n".format(sdict[test['sys']]),
-                  res, '\n', test['result'])
-            assert_array_almost_equal(
-                res, test['result'], test['digits'])
-        sdict = self.yazdan
-        for test in self.ymargin:
-            res = margin(sdict[test['sys']]*test['K'])
-            print("more margin {}\n".format(sdict[test['sys']]),
-                  res, '\n', test['result'])
-            assert_array_almost_equal(
-                res, test['result'], test['digits'])
-
-
-if __name__ == "__main__":
-    unittest.main()
+                    [[1., 0.]],           [[0.]]),
+         (inf, 147.0743, inf, nan, 2.5483,   nan),
+         ([], [147.0743], [], [], [2.5483],   [])),
+        ((8.75*(4*s**2+0.4*s+1)) / ((100*s+1)*(s**2+0.22*s+1))
+         / (s**2/(10.**2)+2*0.04*s/10.+1),
+         (2.2716,  97.5941, 0.5591, 10.0053, 0.0850, 9.9918),
+         ([2.2716],  [97.5941, -157.7844, 134.7359], [1.0381, 0.5591],
+          [10.0053], [0.0850,     0.9373,   1.0919], [0.4064, 9.9918])),
+        (1/(1+s),  # no gain/phase crossovers
+         (inf, inf, inf, nan, nan, nan),
+         ([], [], [], [], [], [])),
+        (3*(10+s)/(2+s),  # no gain/phase crossovers
+         (inf, inf, inf, nan, nan, nan),
+         ([], [], [], [], [], [])),
+        (0.01*(10-s)/(2+s)/(1+s),  # no phase crossovers
+         (300.0, inf, 0.9917, 5.6569, nan, 2.3171),
+         ([300.0], [], [0.9917], [5.6569], [], 2.3171))]
+
+
+def compare_allmargins(actual, desired, **kwargs):
+    """Compare all elements of stability_margins(returnall=True) result"""
+    assert len(actual) == len(desired)
+    for a, d in zip(actual, desired):
+        assert_allclose(a, d, **kwargs)
+
+
+@pytest.mark.parametrize("sys, refout, refoutall", tsys)
+def test_stability_margins(sys, refout, refoutall):
+    """Test stability_margins() function"""
+    out = stability_margins(sys)
+    assert_allclose(out, refout, atol=1.5e-2)
+    out = stability_margins(sys, returnall=True)
+    compare_allmargins(out, refoutall, atol=1.5e-2)
+
+
+@pytest.mark.parametrize("sys, refout, refoutall", tsys)
+def test_stability_margins_omega(sys, refout, refoutall):
+    """Test stability_margins() with interpolated frequencies"""
+    omega = np.logspace(-2, 2, 2000)
+    out = stability_margins(FrequencyResponseData(sys, omega))
+    assert_allclose(out, refout, atol=1.5e-3)
+
+
+@pytest.mark.parametrize("sys, refout, refoutall", tsys)
+def test_stability_margins_3input(sys, refout, refoutall):
+    """Test stability_margins() function with mag, phase, omega input"""
+    omega = np.logspace(-2, 2, 2000)
+    mag, phase, omega_ = sys.freqresp(omega)
+    out = stability_margins((mag, phase*180/np.pi, omega_))
+    assert_allclose(out, refout, atol=1.5e-3)
+
+
+@pytest.mark.parametrize("sys, refout, refoutall", tsys)
+def test_margin_sys(sys, refout, refoutall):
+    """Test margin() function with system input"""
+    out = margin(sys)
+    assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
+
+
+@pytest.mark.parametrize("sys, refout, refoutall", tsys)
+def test_margin_3input(sys, refout, refoutall):
+    """Test margin() function with mag, phase, omega input"""
+    omega = np.logspace(-2, 2, 2000)
+    mag, phase, omega_ = sys.freqresp(omega)
+    out = margin((mag, phase*180/np.pi, omega_))
+    assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
+
+
+def test_phase_crossover_frequencies():
+    """Test phase_crossover_frequencies() function"""
+    omega, gain = phase_crossover_frequencies(tsys[1][0])
+    assert_allclose(omega, [1.73205,  0.], atol=1.5e-3)
+    assert_allclose(gain, [-0.5,  0.25], atol=1.5e-3)
+
+    tf = TransferFunction([1], [1, 1])
+    omega, gain = phase_crossover_frequencies(tf)
+    assert_allclose(omega, [0.], atol=1.5e-3)
+    assert_allclose(gain, [1.], atol=1.5e-3)
+
+    # testing MIMO, only (0,0) element is considered
+    tf = TransferFunction([[[1], [2]],
+                           [[3], [4]]],
+                          [[[1, 2, 3, 4], [1, 1]],
+                           [[1, 1], [1, 1]]])
+    omega, gain = phase_crossover_frequencies(tf)
+    assert_allclose(omega, [1.73205,  0.], atol=1.5e-3)
+    assert_allclose(gain, [-0.5,  0.25], atol=1.5e-3)
+
+
+def test_mag_phase_omega():
+    """Test for bug reported in gh-58"""
+    sys = TransferFunction(15, [1, 6, 11, 6])
+    out = stability_margins(sys)
+    omega = np.logspace(-2, 2, 1000)
+    mag, phase, omega = sys.freqresp(omega)
+    out2 = stability_margins((mag, phase*180/np.pi, omega))
+    ind = [0, 1, 3, 4]   # indices of gm, pm, wg, wp -- ignore sm
+    marg1 = np.array(out)[ind]
+    marg2 = np.array(out2)[ind]
+    assert_allclose(marg1, marg2, atol=1.5e-3)
+
+
+def test_frd():
+    """Test FrequencyResonseData margins"""
+    f = np.array([0.005, 0.010, 0.020, 0.030, 0.040,
+                  0.050, 0.060, 0.070, 0.080, 0.090,
+                  0.100, 0.200, 0.300, 0.400, 0.500,
+                  0.750, 1.000, 1.250, 1.500, 1.750,
+                  2.000, 2.250, 2.500, 2.750, 3.000,
+                  3.250, 3.500, 3.750, 4.000, 4.250,
+                  4.500, 4.750, 5.000, 6.000, 7.000,
+                  8.000, 9.000, 10.000])
+    gain = np.array([  0.0,   0.0,   0.0,   0.0,   0.0,
+                       0.0,   0.0,   0.0,   0.0,   0.0,
+                       0.0,   0.1,   0.2,   0.3,   0.5,
+                       0.5,  -0.4,  -2.3,  -4.8,  -7.3,
+                      -9.6, -11.7, -13.6, -15.3, -16.9,
+                     -18.3, -19.6, -20.8, -22.0, -23.1,
+                     -24.1, -25.0, -25.9, -29.1, -31.9,
+                     -34.2, -36.2, -38.1])
+    phase = np.array([  0,    -1,    -2,    -3,    -4,
+                       -5,    -6,    -7,    -8,    -9,
+                      -10,   -19,   -29,   -40,   -51,
+                      -81,  -114,  -144,  -168,  -187,
+                     -202,  -214,  -224,  -233,  -240,
+                     -247,  -253,  -259,  -264,  -269,
+                     -273,  -277,  -280,  -292,  -301,
+                     -307,  -313,  -317])
+    # calculate response as complex number
+    resp = 10**(gain / 20) * np.exp(1j * phase / (180./np.pi))
+    # frequency response data
+    fresp = FrequencyResponseData(resp, f*2*np.pi, smooth=True)
+    s = TransferFunction([1, 0], [1])
+    G = 1./(s**2)
+    K = 1.
+    C = K*(1+1.9*s)
+    TFopen = fresp*C*G
+    gm, pm, sm, wg, wp, ws = stability_margins(TFopen)
+    assert_allclose([pm], [44.55], atol=.01)
+
+
+def test_frd_indexing():
+    """Test FRD edge cases
+
+    Make sure frd objects with non benign data do not raise exceptions when
+    the stability criteria evaluate at the first or last frequency point
+    bug reported in gh-407
+    """
+    # frequency points just a little under 1. and over 2.
+    w = np.linspace(.99, 2.01, 11)
+
+    # Note: stability_margins will convert the frd with smooth=True
+
+    # gain margins
+    # p crosses -180 at w[0]=1. and w[-1]=2.
+    m = 0.6
+    p = -180*(2*w-1)
+    d = m*np.exp(1J*np.pi/180*p)
+    frd_gm = FrequencyResponseData(d, w)
+    gm, _, _, wg, _, _ = stability_margins(frd_gm, returnall=True)
+    assert_allclose(gm, [1/m, 1/m], atol=0.01)
+    assert_allclose(wg, [1., 2.], atol=0.01)
+
+    # phase margins
+    # m crosses 1 at w[0]=1. and w[-1]=2.
+    m = -(2*w-3)**4 + 2
+    p = -90.
+    d = m*np.exp(1J*np.pi/180*p)
+    frd_pm = FrequencyResponseData(d, w)
+    _, pm, _, _, wp, _ = stability_margins(frd_pm, returnall=True)
+    assert_allclose(pm, [90., 90.], atol=0.01)
+    assert_allclose(wp, [1., 2.], atol=0.01)
+
+    # stability margins
+    # minimum abs(d+1)=1-m at w[1]=1. and w[-2]=2., in nyquist plot
+    w = np.arange(.9, 2.1, 0.1)
+    m = 0.6
+    p = -180*(2*w-1)
+    d = m*np.exp(1J*np.pi/180*p)
+    frd_sm = FrequencyResponseData(d, w)
+    _, _, sm, _, _, ws = stability_margins(frd_sm, returnall=True)
+    assert_allclose(sm, [1-m, 1-m], atol=0.01)
+    assert_allclose(ws, [1., 2.], atol=0.01)
+
+
+"""
+NOTE:
+Matlab gives gain margin 0 for system `type2`, python-control gives inf
+Difficult to argue which is right? Special case or different approach?
+
+Edge cases, like `type0` 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.
+"""
+tzmore_sys = {
+    '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)}
+tzmore_margin = [
+    dict(sys='typem1', K=2.0, atol=1.5e-3, result=(
+        float('Inf'), -120.0007, float('NaN'), 0.5774)),
+    dict(sys='type0', K=0.8, atol=1.5e-3, result=(
+        10.0014, float('inf'), 1.7322, float('nan'))),
+    dict(sys='type0', K=2.0, atol=1e-2, result=(
+        4.000,  67.6058,  1.7322,   0.7663)),
+    dict(sys='type1', K=1.0, atol=1e-4, result=(
+        float('Inf'), 144.9032, float('NaN'), 0.3162)),
+    dict(sys='type2', K=1.0, atol=1e-4, result=(
+        float('Inf'), 44.4594, float('NaN'), 0.7907)),
+    dict(sys='type3', K=1.0, atol=1.5e-3, result=(
+        0.0626, 37.1748, 0.1119, 0.7951)),
+    ]
+tzmore_stability_margins = []
+
+"""
+from "A note on the Gain and Phase Margin Concepts
+Journal of Control and Systems Engineering, Yazdan Bavafi-Toosi,
+Dec 2015, vol 3 iss 1, pp 51-59
+
+A cornucopia of tricky systems for phase / gain margin
+TODO: still have to convert more to tests + fix margin to handle
+also these torture cases
+"""
+yazdan = {
+    'example21':
+    0.002*(s+0.02)*(s+0.05)*(s+5)*(s+10)/(
+        (s-0.0005)*(s+0.0001)*(s+0.01)*(s+0.2)*(s+1)*(s+100)**2),
+    'example23':
+    ((s+0.1)**2 + 1)*(s-0.1)/(
+        ((s+0.1)**2+4)*(s+1)),
+    'example25a':
+    s/(s**2+2*s+2)**4,
+    'example26a':
+    ((s-0.1)**2 + 1)/(
+        (s + 0.1)*((s-0.2)**2 + 4)),
+    'example26b': ((s-0.1)**2 + 1)/(
+        (s - 0.3)*((s-0.2)**2 + 4))
+}
+yazdan['example24'] = yazdan['example21']*20000
+yazdan['example25b'] = yazdan['example25a']*100
+yazdan['example22'] = yazdan['example21']*(s**2 - 2*s + 401)
+ymargin = [
+    dict(sys='example21', K=1.0, atol=1e-2,
+         result=(0.0100, -14.5640,  0, 0.0022)),
+    dict(sys='example21', K=1000.0, atol=1e-2,
+         result=(0.1793, 22.5215, 0.0243, 0.0630)),
+    dict(sys='example21', K=5000.0, atol=1.5e-3,
+         result=(4.5596, 21.2101, 0.4385, 0.1868)),
+    ]
+ystability_margins = [
+    dict(sys='example21', K=1.0, rtol=1e-3, atol=1e-3,
+         result=([0.01, 179.2931, 2.2798e+4, 1.5946e+07, 7.2477e+08],
+                 [-14.5640],
+                 [0.2496],
+                 [0, 0.0243, 0.4385, 6.8640, 84.9323],
+                 [0.0022],
+                 [0.0022])),
+    ]
+
+tzmore_sys.update(yazdan)
+tzmore_margin += ymargin
+tzmore_stability_margins += ystability_margins
+
+
+@pytest.mark.parametrize('tmargin', tzmore_margin)
+def test_zmore_margin(tmargin):
+    """Test margins for more tricky systems"""
+    res = margin(tzmore_sys[tmargin['sys']]*tmargin['K'])
+    assert_allclose(res, tmargin['result'], atol=tmargin['atol'])
+
+
+@pytest.mark.parametrize('tmarginall', tzmore_stability_margins)
+def test_zmore_stability_margins(tmarginall):
+    """Test stability_margins for more tricky systems with returnall"""
+    res = stability_margins(tzmore_sys[tmarginall['sys']]*tmarginall['K'],
+                            returnall=True)
+    compare_allmargins(res, tmarginall['result'],
+                       atol=tmarginall['atol'],
+                       rtol=tmarginall['rtol'])
diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py
index 647ddd20..966f700d 100644
--- a/control/tests/rlocus_test.py
+++ b/control/tests/rlocus_test.py
@@ -4,13 +4,14 @@
 # RMM, 1 Jul 2011
 
 import unittest
+import matplotlib.pyplot as plt
 import numpy as np
+from numpy.testing import assert_array_almost_equal
+
 from control.rlocus import root_locus, _RLClickDispatcher
 from control.xferfcn import TransferFunction
 from control.statesp import StateSpace
 from control.bdalg import feedback
-import matplotlib.pyplot as plt
-from control.tests.margin_test import assert_array_almost_equal
 
 
 class TestRootLocus(unittest.TestCase):
diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py
index f2cdf910..e0012a37 100644
--- a/control/tests/sisotool_test.py
+++ b/control/tests/sisotool_test.py
@@ -1,34 +1,42 @@
 import unittest
+import matplotlib.pyplot as plt
 import numpy as np
+from numpy.testing import assert_array_almost_equal
+
 from control.sisotool import sisotool
-from control.tests.margin_test import assert_array_almost_equal
 from control.rlocus import _RLClickDispatcher
 from control.xferfcn import TransferFunction
-import matplotlib.pyplot as plt
+
 
 class TestSisotool(unittest.TestCase):
     """These are tests for the sisotool in sisotool.py."""
 
     def setUp(self):
         # One random SISO system.
-        self.system = TransferFunction([1000],[1,25,100,0])
+        self.system = TransferFunction([1000], [1, 25, 100, 0])
 
     def test_sisotool(self):
-        sisotool(self.system,Hz=False)
+        sisotool(self.system, Hz=False)
         fig = plt.gcf()
-        ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3]
+        ax_mag, ax_rlocus, ax_phase, ax_step = fig.axes[:4]
 
         # Check the initial root locus plot points
-        initial_point_0 = (np.array([-22.53155977]),np.array([0.]))
+        initial_point_0 = (np.array([-22.53155977]), np.array([0.]))
         initial_point_1 = (np.array([-1.23422011]), np.array([-6.54667031]))
         initial_point_2 = (np.array([-1.23422011]), np.array([06.54667031]))
-        assert_array_almost_equal(ax_rlocus.lines[0].get_data(),initial_point_0)
-        assert_array_almost_equal(ax_rlocus.lines[1].get_data(),initial_point_1)
-        assert_array_almost_equal(ax_rlocus.lines[2].get_data(),initial_point_2)
+        assert_array_almost_equal(ax_rlocus.lines[0].get_data(),
+                                  initial_point_0, 4)
+        assert_array_almost_equal(ax_rlocus.lines[1].get_data(),
+                                  initial_point_1, 4)
+        assert_array_almost_equal(ax_rlocus.lines[2].get_data(),
+                                  initial_point_2, 4)
 
         # Check the step response before moving the point
-        step_response_original = np.array([ 0.,          0.02233651,  0.13118374,  0.33078542,  0.5907113,   0.87041549, 1.13038536,  1.33851053,  1.47374666,  1.52757114])
-        assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_original)
+        step_response_original = np.array(
+            [0., 0.02233651, 0.13118374, 0.33078542, 0.5907113, 0.87041549,
+             1.13038536, 1.33851053, 1.47374666, 1.52757114])
+        assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],
+                                  step_response_original, 4)
 
         bode_plot_params = {
             'omega': None,
@@ -43,24 +51,37 @@ class TestSisotool(unittest.TestCase):
         }
 
         # Move the rootlocus to another point
-        event = type('test', (object,), {'xdata': 2.31206868287,'ydata':15.5983051046, 'inaxes':ax_rlocus.axes})()
-        _RLClickDispatcher(event=event, sys=self.system, fig=fig,ax_rlocus=ax_rlocus,sisotool=True, plotstr='-' ,bode_plot_params=bode_plot_params, tvect=None)
+        event = type('test', (object,), {'xdata': 2.31206868287,
+                                         'ydata': 15.5983051046,
+                                         'inaxes': ax_rlocus.axes})()
+        _RLClickDispatcher(event=event, sys=self.system, fig=fig,
+                           ax_rlocus=ax_rlocus, sisotool=True, plotstr='-',
+                           bode_plot_params=bode_plot_params, tvect=None)
 
         # Check the moved root locus plot points
         moved_point_0 = (np.array([-29.91742755]), np.array([0.]))
         moved_point_1 = (np.array([2.45871378]), np.array([-15.52647768]))
         moved_point_2 = (np.array([2.45871378]), np.array([15.52647768]))
-        assert_array_almost_equal(ax_rlocus.lines[-3].get_data(),moved_point_0)
-        assert_array_almost_equal(ax_rlocus.lines[-2].get_data(),moved_point_1)
-        assert_array_almost_equal(ax_rlocus.lines[-1].get_data(),moved_point_2)
+        assert_array_almost_equal(ax_rlocus.lines[-3].get_data(),
+                                  moved_point_0, 4)
+        assert_array_almost_equal(ax_rlocus.lines[-2].get_data(),
+                                  moved_point_1, 4)
+        assert_array_almost_equal(ax_rlocus.lines[-1].get_data(),
+                                  moved_point_2, 4)
 
         # Check if the bode_mag line has moved
-        bode_mag_moved = np.array([ 111.83321224,   92.29238035,   76.02822315,   62.46884113,   51.14108703, 41.6554004,    33.69409534,   27.00237344,   21.38086717,   16.67791585])
-        assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],bode_mag_moved)
+        bode_mag_moved = np.array(
+            [111.83321224, 92.29238035, 76.02822315, 62.46884113, 51.14108703,
+             41.6554004, 33.69409534, 27.00237344, 21.38086717, 16.67791585])
+        assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],
+                                  bode_mag_moved, 4)
 
         # Check if the step response has changed
-        step_response_moved = np.array([[ 0.,          0.02458187,  0.16529784 , 0.46602716 , 0.91012035 , 1.43364313, 1.93996334 , 2.3190105  , 2.47041552 , 2.32724853] ])
-        assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_moved)
+        step_response_moved = np.array(
+            [0., 0.02458187, 0.16529784, 0.46602716, 0.91012035, 1.43364313,
+             1.93996334, 2.3190105, 2.47041552, 2.32724853])
+        assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],
+                                  step_response_moved, 4)
 
 
 if __name__ == "__main__":