diff --git a/control/robust.py b/control/robust.py
index 54e970a5d..b69a8ba11 100644
--- a/control/robust.py
+++ b/control/robust.py
@@ -43,9 +43,13 @@
 import numpy as np
 from .exception import *
 from .statesp import StateSpace
+from .xferfcn import tf
 from .statefbk import *
+import numbers
 
-def h2syn(P,nmeas,ncon):
+__all__ = ['h2syn', 'hinfsyn', 'augw', 'mixsyn', 'makeweight']
+
+def h2syn(P, nmeas, ncon):
     """H_2 control synthesis for plant P.
 
     Parameters
@@ -72,14 +76,14 @@ def h2syn(P,nmeas,ncon):
     >>> K = h2syn(P,nmeas,ncon)
 
     """
-    #Check for ss system object, need a utility for this?
-
-    #TODO: Check for continous or discrete, only continuous supported right now
-        # if isCont():
-        #    dico = 'C'
-        # elif isDisc():
-        #    dico = 'D'
-        # else:
+    # Check for ss system object, need a utility for this?
+
+    # TODO: Check for continous or discrete, only continuous supported right now
+    # if isCont():
+    #    dico = 'C'
+    # elif isDisc():
+    #    dico = 'D'
+    # else:
     dico = 'C'
 
     try:
@@ -87,10 +91,10 @@ def h2syn(P,nmeas,ncon):
     except ImportError:
         raise ControlSlycot("can't find slycot subroutine sb10hd")
 
-    n = np.size(P.A,0)
-    m = np.size(P.B,1)
-    np_ = np.size(P.C,0)
-    out = sb10hd(n,m,np_,ncon,nmeas,P.A,P.B,P.C,P.D)
+    n = np.size(P.A, 0)
+    m = np.size(P.B, 1)
+    np_ = np.size(P.C, 0)
+    out = sb10hd(n, m, np_, ncon, nmeas, P.A, P.B, P.C, P.D)
     Ak = out[0]
     Bk = out[1]
     Ck = out[2]
@@ -100,7 +104,8 @@ def h2syn(P,nmeas,ncon):
 
     return K
 
-def hinfsyn(P,nmeas,ncon):
+
+def hinfsyn(P, nmeas, ncon):
     """H_{inf} control synthesis for plant P.
 
     Parameters
@@ -136,14 +141,14 @@ def hinfsyn(P,nmeas,ncon):
 
     """
 
-    #Check for ss system object, need a utility for this?
+    # Check for ss system object, need a utility for this?
 
-    #TODO: Check for continous or discrete, only continuous supported right now
-        # if isCont():
-        #    dico = 'C'
-        # elif isDisc():
-        #    dico = 'D'
-        # else:
+    # TODO: Check for continous or discrete, only continuous supported right now
+    # if isCont():
+    #    dico = 'C'
+    # elif isDisc():
+    #    dico = 'D'
+    # else:
     dico = 'C'
 
     try:
@@ -151,11 +156,11 @@ def hinfsyn(P,nmeas,ncon):
     except ImportError:
         raise ControlSlycot("can't find slycot subroutine sb10ad")
 
-    n = np.size(P.A,0)
-    m = np.size(P.B,1)
-    np_ = np.size(P.C,0)
+    n = np.size(P.A, 0)
+    m = np.size(P.B, 1)
+    np_ = np.size(P.C, 0)
     gamma = 1.e100
-    out = sb10ad(n,m,np_,ncon,nmeas,gamma,P.A,P.B,P.C,P.D)
+    out = sb10ad(n, m, np_, ncon, nmeas, gamma, P.A, P.B, P.C, P.D)
     gam = out[0]
     Ak = out[1]
     Bk = out[2]
@@ -202,21 +207,21 @@ def _size_as_needed(w, wname, n):
     """
     from . import append, ss
     if w is not None:
-        if not isinstance(w,StateSpace):
+        if not isinstance(w, StateSpace):
             w = ss(w)
-        if 1==w.inputs and 1==w.outputs:
-            w = append(*(w,)*n)
+        if 1 == w.inputs and 1 == w.outputs:
+            w = append(*(w,) * n)
         else:
             if w.inputs != n:
-                msg=("{}: weighting function has {} inputs, expected {}".
-                      format(wname,w.inputs,n))
+                msg = ("{}: weighting function has {} inputs, expected {}".
+                       format(wname, w.inputs, n))
                 raise ValueError(msg)
     else:
-        w = ss([],[],[],[])
+        w = ss([], [], [], [])
     return w
 
 
-def augw(g,w1=None,w2=None,w3=None):
+def augw(g, w1=None, w2=None, w3=None):
     """Augment plant for mixed sensitivity problem.
 
     Parameters
@@ -225,7 +230,6 @@ def augw(g,w1=None,w2=None,w3=None):
     w1: weighting on S; None, scalar, or k1-by-ny LTI object
     w2: weighting on KS; None, scalar, or k2-by-nu LTI object
     w3: weighting on T; None, scalar, or k3-by-ny LTI object
-    p: augmented plant; StateSpace object
 
     If a weighting is None, no augmentation is done for it.  At least
     one weighting must not be None.
@@ -254,12 +258,12 @@ def augw(g,w1=None,w2=None,w3=None):
     ny = g.outputs
     nu = g.inputs
 
-    w1,w2,w3 = [_size_as_needed(w,wname,n)
-                for w,wname,n in zip((w1,w2,w3),
-                                     ('w1','w2','w3'),
-                                     (ny,nu,ny))]
+    w1, w2, w3 = [_size_as_needed(w, wname, n)
+                  for w, wname, n in zip((w1, w2, w3),
+                                         ('w1', 'w2', 'w3'),
+                                         (ny, nu, ny))]
 
-    if not isinstance(g,StateSpace):
+    if not isinstance(g, StateSpace):
         g = ss(g)
 
     #       w         u
@@ -270,11 +274,11 @@ def augw(g,w1=None,w2=None,w3=None):
     #  v  [ I    |    -g    ]
 
     # error summer: inputs are -y and r=w
-    Ie = ss([],[],[],np.eye(ny))
+    Ie = ss([], [], [], np.eye(ny))
     # control: needed to "distribute" control input
-    Iu = ss([],[],[],np.eye(nu))
+    Iu = ss([], [], [], np.eye(nu))
 
-    sysall = append(w1,w2,w3,Ie,g,Iu)
+    sysall = append(w1, w2, w3, Ie, g, Iu)
 
     niw1 = w1.inputs
     niw2 = w2.inputs
@@ -284,43 +288,45 @@ def augw(g,w1=None,w2=None,w3=None):
     now2 = w2.outputs
     now3 = w3.outputs
 
-    q = np.zeros((niw1+niw2+niw3+ny+nu,2))
-    q[:,0] = np.arange(1,q.shape[0]+1)
+    q = np.zeros((niw1 + niw2 + niw3 + ny + nu, 2))
+    q[:, 0] = np.arange(1, q.shape[0] + 1)
 
     # Ie -> w1
-    q[:niw1,1] = np.arange(1+now1+now2+now3,
-                           1+now1+now2+now3+niw1)
+    q[:niw1, 1] = np.arange(1 + now1 + now2 + now3,
+                            1 + now1 + now2 + now3 + niw1)
 
     # Iu -> w2
-    q[niw1:niw1+niw2,1] = np.arange(1+now1+now2+now3+2*ny,
-                                    1+now1+now2+now3+2*ny+niw2)
+    q[niw1:niw1 + niw2, 1] = np.arange(1 + now1 + now2 + now3 + 2 * ny,
+                                       1 + now1 + now2 + now3 + 2 * ny + niw2)
 
     # y -> w3
-    q[niw1+niw2:niw1+niw2+niw3,1] = np.arange(1+now1+now2+now3+ny,
-                                              1+now1+now2+now3+ny+niw3)
+    q[niw1 + niw2:niw1 + niw2 + niw3, 1] = np.arange(1 + now1 + now2 + now3 + ny,
+                                                     1 + now1 + now2 + now3 + ny + niw3)
 
     # -y -> Iy; note the leading -
-    q[niw1+niw2+niw3:niw1+niw2+niw3+ny,1] = -np.arange(1+now1+now2+now3+ny,
-                                                       1+now1+now2+now3+2*ny)
+    q[niw1 + niw2 + niw3:niw1 + niw2 + niw3 + ny, 1] = -np.arange(1 + now1 + now2 + now3 + ny,
+                                                                  1 + now1 + now2 + now3 + 2 * ny)
 
     # Iu -> G
-    q[niw1+niw2+niw3+ny:niw1+niw2+niw3+ny+nu,1] = np.arange(1+now1+now2+now3+2*ny,
-                                                            1+now1+now2+now3+2*ny+nu)
+    q[niw1 + niw2 + niw3 + ny:niw1 + niw2 + niw3 + ny + nu, 1] = np.arange(
+        1 + now1 + now2 + now3 + 2 * ny,
+        1 + now1 + now2 + now3 + 2 * ny + nu)
 
     # input indices: to Ie and Iu
-    ii = np.hstack((np.arange(1+now1+now2+now3,
-                              1+now1+now2+now3+ny),
-                    np.arange(1+now1+now2+now3+ny+nu,
-                              1+now1+now2+now3+ny+nu+nu)))
+    ii = np.hstack((np.arange(1 + now1 + now2 + now3,
+                              1 + now1 + now2 + now3 + ny),
+                    np.arange(1 + now1 + now2 + now3 + ny + nu,
+                              1 + now1 + now2 + now3 + ny + nu + nu)))
 
     # output indices
-    oi = np.arange(1,1+now1+now2+now3+ny)
+    oi = np.arange(1, 1 + now1 + now2 + now3 + ny)
 
-    p = connect(sysall,q,ii,oi)
+    p = connect(sysall, q, ii, oi)
 
     return p
 
-def mixsyn(g,w1=None,w2=None,w3=None):
+
+def mixsyn(g, w1=None, w2=None, w3=None):
     """Mixed-sensitivity H-infinity synthesis.
 
     mixsyn(g,w1,w2,w3) -> k,cl,info
@@ -356,8 +362,60 @@ def mixsyn(g,w1=None,w2=None,w3=None):
     """
     nmeas = g.outputs
     ncon = g.inputs
-    p = augw(g,w1,w2,w3)
+    p = augw(g, w1, w2, w3)
+
+    k, cl, gamma, rcond = hinfsyn(p, nmeas, ncon)
+    info = gamma, rcond
+    return k, cl, info
+
 
-    k,cl,gamma,rcond=hinfsyn(p,nmeas,ncon)
-    info = gamma,rcond
-    return k,cl,info
+def makeweight(dcgain, wc, hfgain, order=1, Ts=None):
+    """ Compute a weight transfer function, noted W
+
+    Parameters
+    ----------
+        dcgain : double positive scalar
+            Low frequency gain of 1/W; should be < 1
+        wc : double strictly positive scalar
+            Design frequency (where |W| is approximately 1)
+        hfgain : double positive scalar
+            High frequency gain of 1/W; should be > 1
+        order : int strictly positive scalar - optional
+            Order of the output transfer function
+        Ts : double strictly positive scalar - optional
+            Sampling time in case W is discrete
+    W - SISO LTI object
+    """
+    if not isinstance(dcgain, numbers.Real):
+        raise ControlArgument("dcgain should be a real scalar number")
+    if dcgain < 0:
+        raise ValueError("dcgain should be a positive real scalar number")
+
+    if not isinstance(wc, numbers.Real):
+        raise ControlArgument("wc should be a real scalar number")
+    if wc <= 0:
+        raise ValueError("wc should be a strictly positive real scalar number")
+
+    if not isinstance(hfgain, numbers.Real):
+        raise ControlArgument("hfgain should be a real scalar number")
+    if hfgain < 0:
+        raise ValueError("hfgain should be a positive real scalar number")
+
+    if not isinstance(order, numbers.Integral):
+        raise ControlArgument("order should be an integral scalar number")
+    if order <= 0:
+        raise ValueError("order should be a strictly positive integral scalar number")
+
+    # TODO: Implement discrete behavior
+    if Ts is not None:
+        raise ControlNotImplemented("makeweight cannot return discrete weight function")
+    if Ts is not None and not isinstance(Ts, numbers.Real):
+        raise ControlArgument("Ts should be a real scalar number")
+    if Ts is not None and Ts <= 0:
+        raise ValueError("Ts should be a strictly positive real scalar number")
+
+    s = tf([1, 0], [1])
+
+    W = (s + wc * (dcgain ** (1. / order))) ** order / (s / (hfgain ** (1. / order)) + wc) ** order
+
+    return W
diff --git a/control/tests/robust_test.py b/control/tests/robust_test.py
index e948f8e2c..aa016dc40 100644
--- a/control/tests/robust_test.py
+++ b/control/tests/robust_test.py
@@ -2,13 +2,13 @@
 import numpy as np
 import control
 import control.robust
-from control.exception import slycot_check
+from control.exception import slycot_check, ControlArgument, ControlNotImplemented
 
 
 class TestHinf(unittest.TestCase):
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testHinfsyn(self):
-        "Test hinfsyn"
+    def test_hinfsyn(self):
+        """Test hinfsyn"""
         p = control.ss(-1, [1, 1], [[1], [1]], [[0, 1], [1, 0]])
         k, cl, gam, rcond = control.robust.hinfsyn(p, 1, 1)
         # from Octave, which also uses SB10AD:
@@ -26,10 +26,11 @@ def testHinfsyn(self):
 
     # TODO: add more interesting examples
 
+
 class TestH2(unittest.TestCase):
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testH2syn(self):
-        "Test h2syn"
+    def test_h2syn(self):
+        """Test h2syn"""
         p = control.ss(-1, [1, 1], [[1], [1]], [[0, 1], [1, 0]])
         k = control.robust.h2syn(p, 1, 1)
         # from Octave, which also uses SB10HD for H-2 synthesis:
@@ -44,303 +45,298 @@ def testH2syn(self):
 
 
 class TestAugw(unittest.TestCase):
-    "Test control.robust.augw"
+    """Test control.robust.augw"""
 
     # tolerance for system equality
     TOL = 1e-8
 
-    def siso_almost_equal(self,g,h):
+    def siso_almost_equal(self, g, h):
         """siso_almost_equal(g,h) -> None
         Raises AssertionError if g and h, two SISO LTI objects, are not almost equal"""
         from control import tf, minreal
-        gmh = tf(minreal(g-h,verbose=False))
-        if not (gmh.num[0][0]<self.TOL).all():
+        gmh = tf(minreal(g - h, verbose=False))
+        if not (gmh.num[0][0] < self.TOL).all():
             maxnum = max(abs(gmh.num[0][0]))
-            raise AssertionError('systems not approx equal; max num. coeff is {}\nsys 1:\n{}\nsys 2:\n{}'.format(maxnum,g,h))
-
+            raise AssertionError(
+                'systems not approx equal; max num. coeff is {}\nsys 1:\n{}\nsys 2:\n{}'.format(
+                    maxnum, g, h))
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testSisoW1(self):
-        "SISO plant with S weighting"
+    def test_augw_siso_w1(self):
+        """SISO plant with S weighting"""
         from control import augw, ss
-        g = ss([-1.],[1.],[1.],[1.])
-        w1 = ss([-2],[2.],[1.],[2.])
-        p = augw(g,w1)
-        self.assertEqual(2,p.outputs)
-        self.assertEqual(2,p.inputs)
+        g = ss([-1.], [1.], [1.], [1.])
+        w1 = ss([-2], [2.], [1.], [2.])
+        p = augw(g, w1)
+        self.assertEqual(2, p.outputs)
+        self.assertEqual(2, p.inputs)
         # w->z1 should be w1
-        self.siso_almost_equal(w1,p[0,0])
+        self.siso_almost_equal(w1, p[0, 0])
         # w->v should be 1
-        self.siso_almost_equal(ss([],[],[],[1]),p[1,0])
+        self.siso_almost_equal(ss([], [], [], [1]), p[1, 0])
         # u->z1 should be -w1*g
-        self.siso_almost_equal(-w1*g,p[0,1])
+        self.siso_almost_equal(-w1 * g, p[0, 1])
         # u->v should be -g
-        self.siso_almost_equal(-g,p[1,1])
-
+        self.siso_almost_equal(-g, p[1, 1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testSisoW2(self):
-        "SISO plant with KS weighting"
+    def test_augw_siso_w2(self):
+        """SISO plant with KS weighting"""
         from control import augw, ss
-        g = ss([-1.],[1.],[1.],[1.])
-        w2 = ss([-2],[1.],[1.],[2.])
-        p = augw(g,w2=w2)
-        self.assertEqual(2,p.outputs)
-        self.assertEqual(2,p.inputs)
+        g = ss([-1.], [1.], [1.], [1.])
+        w2 = ss([-2], [1.], [1.], [2.])
+        p = augw(g, w2=w2)
+        self.assertEqual(2, p.outputs)
+        self.assertEqual(2, p.inputs)
         # w->z2 should be 0
-        self.siso_almost_equal(ss([],[],[],0),p[0,0])
+        self.siso_almost_equal(ss([], [], [], 0), p[0, 0])
         # w->v should be 1
-        self.siso_almost_equal(ss([],[],[],[1]),p[1,0])
+        self.siso_almost_equal(ss([], [], [], [1]), p[1, 0])
         # u->z2 should be w2
-        self.siso_almost_equal(w2,p[0,1])
+        self.siso_almost_equal(w2, p[0, 1])
         # u->v should be -g
-        self.siso_almost_equal(-g,p[1,1])
-
+        self.siso_almost_equal(-g, p[1, 1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testSisoW3(self):
-        "SISO plant with T weighting"
+    def test_augw_siso_w3(self):
+        """SISO plant with T weighting"""
         from control import augw, ss
-        g = ss([-1.],[1.],[1.],[1.])
-        w3 = ss([-2],[1.],[1.],[2.])
-        p = augw(g,w3=w3)
-        self.assertEqual(2,p.outputs)
-        self.assertEqual(2,p.inputs)
+        g = ss([-1.], [1.], [1.], [1.])
+        w3 = ss([-2], [1.], [1.], [2.])
+        p = augw(g, w3=w3)
+        self.assertEqual(2, p.outputs)
+        self.assertEqual(2, p.inputs)
         # w->z3 should be 0
-        self.siso_almost_equal(ss([],[],[],0),p[0,0])
+        self.siso_almost_equal(ss([], [], [], 0), p[0, 0])
         # w->v should be 1
-        self.siso_almost_equal(ss([],[],[],[1]),p[1,0])
+        self.siso_almost_equal(ss([], [], [], [1]), p[1, 0])
         # u->z3 should be w3*g
-        self.siso_almost_equal(w3*g,p[0,1])
+        self.siso_almost_equal(w3 * g, p[0, 1])
         # u->v should be -g
-        self.siso_almost_equal(-g,p[1,1])
-
+        self.siso_almost_equal(-g, p[1, 1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testSisoW123(self):
-        "SISO plant with all weights"
+    def test_augw_siso_w123(self):
+        """SISO plant with all weights"""
         from control import augw, ss
-        g = ss([-1.],[1.],[1.],[1.])
-        w1 = ss([-2.],[2.],[1.],[2.])
-        w2 = ss([-3.],[3.],[1.],[3.])
-        w3 = ss([-4.],[4.],[1.],[4.])
-        p = augw(g,w1,w2,w3)
-        self.assertEqual(4,p.outputs)
-        self.assertEqual(2,p.inputs)
+        g = ss([-1.], [1.], [1.], [1.])
+        w1 = ss([-2.], [2.], [1.], [2.])
+        w2 = ss([-3.], [3.], [1.], [3.])
+        w3 = ss([-4.], [4.], [1.], [4.])
+        p = augw(g, w1, w2, w3)
+        self.assertEqual(4, p.outputs)
+        self.assertEqual(2, p.inputs)
         # w->z1 should be w1
-        self.siso_almost_equal(w1,p[0,0])
+        self.siso_almost_equal(w1, p[0, 0])
         # w->z2 should be 0
-        self.siso_almost_equal(0,p[1,0])
+        self.siso_almost_equal(0, p[1, 0])
         # w->z3 should be 0
-        self.siso_almost_equal(0,p[2,0])
+        self.siso_almost_equal(0, p[2, 0])
         # w->v should be 1
-        self.siso_almost_equal(ss([],[],[],[1]),p[3,0])
+        self.siso_almost_equal(ss([], [], [], [1]), p[3, 0])
         # u->z1 should be -w1*g
-        self.siso_almost_equal(-w1*g,p[0,1])
+        self.siso_almost_equal(-w1 * g, p[0, 1])
         # u->z2 should be w2
-        self.siso_almost_equal(w2,p[1,1])
+        self.siso_almost_equal(w2, p[1, 1])
         # u->z3 should be w3*g
-        self.siso_almost_equal(w3*g,p[2,1])
+        self.siso_almost_equal(w3 * g, p[2, 1])
         # u->v should be -g
-        self.siso_almost_equal(-g,p[3,1])
-
+        self.siso_almost_equal(-g, p[3, 1])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testMimoW1(self):
-        "MIMO plant with S weighting"
+    def test_augw_mimo_w1(self):
+        """MIMO plant with S weighting"""
         from control import augw, ss
-        g = ss([[-1.,-2],[-3,-4]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]])
-        w1 = ss([-2],[2.],[1.],[2.])
-        p = augw(g,w1)
-        self.assertEqual(4,p.outputs)
-        self.assertEqual(4,p.inputs)
+        g = ss([[-1., -2], [-3, -4]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]])
+        w1 = ss([-2], [2.], [1.], [2.])
+        p = augw(g, w1)
+        self.assertEqual(4, p.outputs)
+        self.assertEqual(4, p.inputs)
         # w->z1 should be diag(w1,w1)
-        self.siso_almost_equal(w1,p[0,0])
-        self.siso_almost_equal(0, p[0,1])
-        self.siso_almost_equal(0, p[1,0])
-        self.siso_almost_equal(w1,p[1,1])
+        self.siso_almost_equal(w1, p[0, 0])
+        self.siso_almost_equal(0, p[0, 1])
+        self.siso_almost_equal(0, p[1, 0])
+        self.siso_almost_equal(w1, p[1, 1])
         # w->v should be I
-        self.siso_almost_equal(1, p[2,0])
-        self.siso_almost_equal(0, p[2,1])
-        self.siso_almost_equal(0, p[3,0])
-        self.siso_almost_equal(1, p[3,1])
+        self.siso_almost_equal(1, p[2, 0])
+        self.siso_almost_equal(0, p[2, 1])
+        self.siso_almost_equal(0, p[3, 0])
+        self.siso_almost_equal(1, p[3, 1])
         # u->z1 should be -w1*g
-        self.siso_almost_equal(-w1*g[0,0],p[0,2])
-        self.siso_almost_equal(-w1*g[0,1],p[0,3])
-        self.siso_almost_equal(-w1*g[1,0],p[1,2])
-        self.siso_almost_equal(-w1*g[1,1],p[1,3])
+        self.siso_almost_equal(-w1 * g[0, 0], p[0, 2])
+        self.siso_almost_equal(-w1 * g[0, 1], p[0, 3])
+        self.siso_almost_equal(-w1 * g[1, 0], p[1, 2])
+        self.siso_almost_equal(-w1 * g[1, 1], p[1, 3])
         # # u->v should be -g
-        self.siso_almost_equal(-g[0,0],p[2,2])
-        self.siso_almost_equal(-g[0,1],p[2,3])
-        self.siso_almost_equal(-g[1,0],p[3,2])
-        self.siso_almost_equal(-g[1,1],p[3,3])
-
+        self.siso_almost_equal(-g[0, 0], p[2, 2])
+        self.siso_almost_equal(-g[0, 1], p[2, 3])
+        self.siso_almost_equal(-g[1, 0], p[3, 2])
+        self.siso_almost_equal(-g[1, 1], p[3, 3])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testMimoW2(self):
-        "MIMO plant with KS weighting"
+    def test_augw_mimo_w2(self):
+        """MIMO plant with KS weighting"""
         from control import augw, ss
-        g = ss([[-1.,-2],[-3,-4]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]])
-        w2 = ss([-2],[2.],[1.],[2.])
-        p = augw(g,w2=w2)
-        self.assertEqual(4,p.outputs)
-        self.assertEqual(4,p.inputs)
+        g = ss([[-1., -2], [-3, -4]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]])
+        w2 = ss([-2], [2.], [1.], [2.])
+        p = augw(g, w2=w2)
+        self.assertEqual(4, p.outputs)
+        self.assertEqual(4, p.inputs)
         # w->z2 should be 0
-        self.siso_almost_equal(0, p[0,0])
-        self.siso_almost_equal(0, p[0,1])
-        self.siso_almost_equal(0, p[1,0])
-        self.siso_almost_equal(0, p[1,1])
+        self.siso_almost_equal(0, p[0, 0])
+        self.siso_almost_equal(0, p[0, 1])
+        self.siso_almost_equal(0, p[1, 0])
+        self.siso_almost_equal(0, p[1, 1])
         # w->v should be I
-        self.siso_almost_equal(1, p[2,0])
-        self.siso_almost_equal(0, p[2,1])
-        self.siso_almost_equal(0, p[3,0])
-        self.siso_almost_equal(1, p[3,1])
+        self.siso_almost_equal(1, p[2, 0])
+        self.siso_almost_equal(0, p[2, 1])
+        self.siso_almost_equal(0, p[3, 0])
+        self.siso_almost_equal(1, p[3, 1])
         # u->z2 should be w2
-        self.siso_almost_equal(w2, p[0,2])
-        self.siso_almost_equal(0,  p[0,3])
-        self.siso_almost_equal(0,  p[1,2])
-        self.siso_almost_equal(w2, p[1,3])
+        self.siso_almost_equal(w2, p[0, 2])
+        self.siso_almost_equal(0, p[0, 3])
+        self.siso_almost_equal(0, p[1, 2])
+        self.siso_almost_equal(w2, p[1, 3])
         # # u->v should be -g
-        self.siso_almost_equal(-g[0,0], p[2,2])
-        self.siso_almost_equal(-g[0,1], p[2,3])
-        self.siso_almost_equal(-g[1,0], p[3,2])
-        self.siso_almost_equal(-g[1,1], p[3,3])
-
+        self.siso_almost_equal(-g[0, 0], p[2, 2])
+        self.siso_almost_equal(-g[0, 1], p[2, 3])
+        self.siso_almost_equal(-g[1, 0], p[3, 2])
+        self.siso_almost_equal(-g[1, 1], p[3, 3])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testMimoW3(self):
-        "MIMO plant with T weighting"
+    def test_augw_mimo_w3(self):
+        """MIMO plant with T weighting"""
         from control import augw, ss
-        g = ss([[-1.,-2],[-3,-4]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]])
-        w3 = ss([-2],[2.],[1.],[2.])
-        p = augw(g,w3=w3)
-        self.assertEqual(4,p.outputs)
-        self.assertEqual(4,p.inputs)
+        g = ss([[-1., -2], [-3, -4]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]])
+        w3 = ss([-2], [2.], [1.], [2.])
+        p = augw(g, w3=w3)
+        self.assertEqual(4, p.outputs)
+        self.assertEqual(4, p.inputs)
         # w->z3 should be 0
-        self.siso_almost_equal(0, p[0,0])
-        self.siso_almost_equal(0, p[0,1])
-        self.siso_almost_equal(0, p[1,0])
-        self.siso_almost_equal(0, p[1,1])
+        self.siso_almost_equal(0, p[0, 0])
+        self.siso_almost_equal(0, p[0, 1])
+        self.siso_almost_equal(0, p[1, 0])
+        self.siso_almost_equal(0, p[1, 1])
         # w->v should be I
-        self.siso_almost_equal(1, p[2,0])
-        self.siso_almost_equal(0, p[2,1])
-        self.siso_almost_equal(0, p[3,0])
-        self.siso_almost_equal(1, p[3,1])
+        self.siso_almost_equal(1, p[2, 0])
+        self.siso_almost_equal(0, p[2, 1])
+        self.siso_almost_equal(0, p[3, 0])
+        self.siso_almost_equal(1, p[3, 1])
         # u->z3 should be w3*g
-        self.siso_almost_equal(w3*g[0,0], p[0,2])
-        self.siso_almost_equal(w3*g[0,1], p[0,3])
-        self.siso_almost_equal(w3*g[1,0], p[1,2])
-        self.siso_almost_equal(w3*g[1,1], p[1,3])
+        self.siso_almost_equal(w3 * g[0, 0], p[0, 2])
+        self.siso_almost_equal(w3 * g[0, 1], p[0, 3])
+        self.siso_almost_equal(w3 * g[1, 0], p[1, 2])
+        self.siso_almost_equal(w3 * g[1, 1], p[1, 3])
         # # u->v should be -g
-        self.siso_almost_equal(-g[0,0], p[2,2])
-        self.siso_almost_equal(-g[0,1], p[2,3])
-        self.siso_almost_equal(-g[1,0], p[3,2])
-        self.siso_almost_equal(-g[1,1], p[3,3])
-
+        self.siso_almost_equal(-g[0, 0], p[2, 2])
+        self.siso_almost_equal(-g[0, 1], p[2, 3])
+        self.siso_almost_equal(-g[1, 0], p[3, 2])
+        self.siso_almost_equal(-g[1, 1], p[3, 3])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testMimoW123(self):
-        "MIMO plant with all weights"
+    def test_augw_mimo_w123(self):
+        """MIMO plant with all weights"""
         from control import augw, ss, append
-        g = ss([[-1.,-2],[-3,-4]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]],
-               [[1.,0.],[0.,1.]])
+        g = ss([[-1., -2], [-3, -4]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]],
+               [[1., 0.], [0., 1.]])
         # this should be expaned to w1*I
-        w1 = ss([-2.],[2.],[1.],[2.])
+        w1 = ss([-2.], [2.], [1.], [2.])
         # diagonal weighting
-        w2 = append(ss([-3.],[3.],[1.],[3.]), ss([-4.],[4.],[1.],[4.]))
+        w2 = append(ss([-3.], [3.], [1.], [3.]), ss([-4.], [4.], [1.], [4.]))
         # full weighting
-        w3 = ss([[-4.,-5],[-6,-7]],
-               [[2.,3.],[5.,7.]],
-               [[11.,13.],[17.,19.]],
-               [[23.,29.],[31.,37.]])
-        p = augw(g,w1,w2,w3)
-        self.assertEqual(8,p.outputs)
-        self.assertEqual(4,p.inputs)
+        w3 = ss([[-4., -5], [-6, -7]],
+                [[2., 3.], [5., 7.]],
+                [[11., 13.], [17., 19.]],
+                [[23., 29.], [31., 37.]])
+        p = augw(g, w1, w2, w3)
+        self.assertEqual(8, p.outputs)
+        self.assertEqual(4, p.inputs)
         # w->z1 should be w1
-        self.siso_almost_equal(w1, p[0,0])
-        self.siso_almost_equal(0,  p[0,1])
-        self.siso_almost_equal(0,  p[1,0])
-        self.siso_almost_equal(w1, p[1,1])
+        self.siso_almost_equal(w1, p[0, 0])
+        self.siso_almost_equal(0, p[0, 1])
+        self.siso_almost_equal(0, p[1, 0])
+        self.siso_almost_equal(w1, p[1, 1])
         # w->z2 should be 0
-        self.siso_almost_equal(0, p[2,0])
-        self.siso_almost_equal(0, p[2,1])
-        self.siso_almost_equal(0, p[3,0])
-        self.siso_almost_equal(0, p[3,1])
+        self.siso_almost_equal(0, p[2, 0])
+        self.siso_almost_equal(0, p[2, 1])
+        self.siso_almost_equal(0, p[3, 0])
+        self.siso_almost_equal(0, p[3, 1])
         # w->z3 should be 0
-        self.siso_almost_equal(0, p[4,0])
-        self.siso_almost_equal(0, p[4,1])
-        self.siso_almost_equal(0, p[5,0])
-        self.siso_almost_equal(0, p[5,1])
+        self.siso_almost_equal(0, p[4, 0])
+        self.siso_almost_equal(0, p[4, 1])
+        self.siso_almost_equal(0, p[5, 0])
+        self.siso_almost_equal(0, p[5, 1])
         # w->v should be I
-        self.siso_almost_equal(1, p[6,0])
-        self.siso_almost_equal(0, p[6,1])
-        self.siso_almost_equal(0, p[7,0])
-        self.siso_almost_equal(1, p[7,1])
+        self.siso_almost_equal(1, p[6, 0])
+        self.siso_almost_equal(0, p[6, 1])
+        self.siso_almost_equal(0, p[7, 0])
+        self.siso_almost_equal(1, p[7, 1])
 
         # u->z1 should be -w1*g
-        self.siso_almost_equal(-w1*g[0,0], p[0,2])
-        self.siso_almost_equal(-w1*g[0,1], p[0,3])
-        self.siso_almost_equal(-w1*g[1,0], p[1,2])
-        self.siso_almost_equal(-w1*g[1,1], p[1,3])
+        self.siso_almost_equal(-w1 * g[0, 0], p[0, 2])
+        self.siso_almost_equal(-w1 * g[0, 1], p[0, 3])
+        self.siso_almost_equal(-w1 * g[1, 0], p[1, 2])
+        self.siso_almost_equal(-w1 * g[1, 1], p[1, 3])
         # u->z2 should be w2
-        self.siso_almost_equal(w2[0,0], p[2,2])
-        self.siso_almost_equal(w2[0,1], p[2,3])
-        self.siso_almost_equal(w2[1,0], p[3,2])
-        self.siso_almost_equal(w2[1,1], p[3,3])
+        self.siso_almost_equal(w2[0, 0], p[2, 2])
+        self.siso_almost_equal(w2[0, 1], p[2, 3])
+        self.siso_almost_equal(w2[1, 0], p[3, 2])
+        self.siso_almost_equal(w2[1, 1], p[3, 3])
         # u->z3 should be w3*g
-        w3g = w3*g;
-        self.siso_almost_equal(w3g[0,0], p[4,2])
-        self.siso_almost_equal(w3g[0,1], p[4,3])
-        self.siso_almost_equal(w3g[1,0], p[5,2])
-        self.siso_almost_equal(w3g[1,1], p[5,3])
+        w3g = w3 * g;
+        self.siso_almost_equal(w3g[0, 0], p[4, 2])
+        self.siso_almost_equal(w3g[0, 1], p[4, 3])
+        self.siso_almost_equal(w3g[1, 0], p[5, 2])
+        self.siso_almost_equal(w3g[1, 1], p[5, 3])
         # u->v should be -g
-        self.siso_almost_equal(-g[0,0], p[6,2])
-        self.siso_almost_equal(-g[0,1], p[6,3])
-        self.siso_almost_equal(-g[1,0], p[7,2])
-        self.siso_almost_equal(-g[1,1], p[7,3])
-
+        self.siso_almost_equal(-g[0, 0], p[6, 2])
+        self.siso_almost_equal(-g[0, 1], p[6, 3])
+        self.siso_almost_equal(-g[1, 0], p[7, 2])
+        self.siso_almost_equal(-g[1, 1], p[7, 3])
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testErrors(self):
-        "Error cases handled"
-        from control import augw,ss
+    def test_augw_errors(self):
+        """Error cases handled"""
+        from control import augw, ss
         # no weights
-        g1by1 = ss(-1,1,1,0)
-        g2by2 = ss(-np.eye(2),np.eye(2),np.eye(2),np.zeros((2,2)))
-        self.assertRaises(ValueError,augw,g1by1)
+        g1by1 = ss(-1, 1, 1, 0)
+        g2by2 = ss(-np.eye(2), np.eye(2), np.eye(2), np.zeros((2, 2)))
+        self.assertRaises(ValueError, augw, g1by1)
         # mismatched size of weight and plant
-        self.assertRaises(ValueError,augw,g1by1,w1=g2by2)
-        self.assertRaises(ValueError,augw,g1by1,w2=g2by2)
-        self.assertRaises(ValueError,augw,g1by1,w3=g2by2)
+        self.assertRaises(ValueError, augw, g1by1, w1=g2by2)
+        self.assertRaises(ValueError, augw, g1by1, w2=g2by2)
+        self.assertRaises(ValueError, augw, g1by1, w3=g2by2)
+
 
 class TestMixsyn(unittest.TestCase):
-    "Test control.robust.mixsyn"
+    """Test control.robust.mixsyn"""
+
     # it's a relatively simple wrapper; compare results with augw, hinfsyn
     @unittest.skipIf(not slycot_check(), "slycot not installed")
-    def testSiso(self):
-        "mixsyn with SISO system"
+    def test_mixsyn_siso(self):
+        """mixsyn with SISO system"""
         from control import tf, augw, hinfsyn, mixsyn
         from control import ss
         # Skogestad+Postlethwaite, Multivariable Feedback Control, 1st Ed., Example 2.11
         s = tf([1, 0], 1)
         # plant
-        g = 200/(10*s+1)/(0.05*s+1)**2
+        g = 200 / (10 * s + 1) / (0.05 * s + 1) ** 2
         # sensitivity weighting
         M = 1.5
         wb = 10
         A = 1e-4
-        w1 = (s/M+wb)/(s+wb*A)
+        w1 = (s / M + wb) / (s + wb * A)
         # KS weighting
         w2 = tf(1, 1)
 
@@ -348,7 +344,7 @@ def testSiso(self):
         kref, clref, gam, rcond = hinfsyn(p, 1, 1)
         ktest, cltest, info = mixsyn(g, w1, w2)
         # check similar to S+P's example
-        np.testing.assert_allclose(gam, 1.37, atol = 1e-2)
+        np.testing.assert_allclose(gam, 1.37, atol=1e-2)
 
         # mixsyn is a convenience wrapper around augw and hinfsyn, so
         # results will be exactly the same.  Given than, use the lazy
@@ -368,5 +364,29 @@ def testSiso(self):
         np.testing.assert_allclose(rcond, info[1])
 
 
+class TestMakeweight(unittest.TestCase):
+    def test_makeweight_invalid_inputs(self):
+        self.assertRaises(ControlArgument, control.makeweight, None, 1, 1)
+        self.assertRaises(ControlArgument, control.makeweight, 1, [1], 1)
+        self.assertRaises(ControlArgument, control.makeweight, 1, 1, np.array(1))
+
+    def test_makeweight_out_of_bound_inputs(self):
+        self.assertRaises(ValueError, control.makeweight, -1, 1, 1)
+        self.assertRaises(ValueError, control.makeweight, 1, 0, 1)
+        self.assertRaises(ValueError, control.makeweight, 1, 1, -1)
+        self.assertRaises(ValueError, control.makeweight, 1, 1, 1, 0)
+        self.assertRaises(ControlNotImplemented, control.makeweight, 1, 1, 1, 1, 0)
+
+    def test_makeweight_first_order(self):
+        W = control.makeweight(100, .4, .10)
+        np.testing.assert_almost_equal(abs(control.evalfr(W, .4j * 1e-6)), 100)
+        np.testing.assert_almost_equal(abs(control.evalfr(W, .4j*1e6)), .10)
+
+    def test_makeweight_second_order(self):
+        W = control.makeweight(100, .4, .10, order=2)
+        np.testing.assert_almost_equal(abs(control.evalfr(W, .4j*1e-6)), 100)
+        np.testing.assert_almost_equal(abs(control.evalfr(W, .4j*1e6)), .10)
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/doc/control.rst b/doc/control.rst
index 1d0b14644..d794f5e96 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -143,6 +143,7 @@ Utility functions and conversions
     issiso
     issys
     mag2db
+    makeweight
     observable_form
     pade
     reachable_form
diff --git a/examples/robust_mimo.py b/examples/robust_mimo.py
index 15abc9263..370ac4949 100644
--- a/examples/robust_mimo.py
+++ b/examples/robust_mimo.py
@@ -1,7 +1,8 @@
-"""Demonstrate mixed-sensitivity H-infinity design for a MIMO plant.
+"""robust_mimo.py
 
-Based on Example 3.8 from Multivariable Feedback Control, Skogestad
-and Postlethwaite, 1st Edition.
+Demonstrate mixed-sensitivity H-infinity design for a MIMO plant.
+
+Based on Example 3.8 from Multivariable Feedback Control, Skogestad and Postlethwaite, 1st Edition.
 """
 
 import os
@@ -9,51 +10,41 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-from control import tf, ss, mixsyn, feedback, step_response
-
-def weighting(wb,m,a):
-    """weighting(wb,m,a) -> wf
-    wb - design frequency (where |wf| is approximately 1)
-    m - high frequency gain of 1/wf; should be > 1
-    a - low frequency gain of 1/wf; should be < 1
-    wf - SISO LTI object
-    """
-    s = tf([1,0],[1])
-    return (s/m+wb)/(s+wb*a)
+from control import tf, ss, mixsyn, feedback, step_response, makeweight
 
 
 def plant():
     """plant() -> g
     g - LTI object; 2x2 plant with a RHP zero, at s=0.5.
     """
-    den = [0.2,1.2,1]
-    gtf=tf([[[1],[1]],
-            [[2,1],[2]]],
-           [[den,den],
-            [den,den]])
+    den = [0.2, 1.2, 1]
+    gtf = tf([[[1], [1]],
+              [[2, 1], [2]]],
+             [[den, den],
+              [den, den]])
     return ss(gtf)
 
 
 # as of this writing (2017-07-01), python-control doesn't have an
 # equivalent to Matlab's sigma function, so use a trivial stand-in.
-def triv_sigma(g,w):
+def triv_sigma(g, w):
     """triv_sigma(g,w) -> s
     g - LTI object, order n
     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)
-    sv = np.linalg.svd(sjw,compute_uv=False)
+    m, p, _ = g.freqresp(w)
+    sjw = (m * np.exp(1j * p * np.pi / 180)).transpose(2, 0, 1)
+    sv = np.linalg.svd(sjw, compute_uv=False)
     return sv
 
 
 def analysis():
     """Plot open-loop responses for various inputs"""
-    g=plant()
+    g = plant()
 
-    t = np.linspace(0,10,101)
-    _, yu1 = step_response(g,t,input=0)
-    _, yu2 = step_response(g,t,input=1)
+    t = np.linspace(0, 10, 101)
+    _, yu1 = step_response(g, t, input=0)
+    _, yu2 = step_response(g, t, input=1)
 
     yu1 = yu1
     yu2 = yu2
@@ -63,35 +54,35 @@ def analysis():
     yuz = yu1 - yu2
 
     plt.figure(1)
-    plt.subplot(1,3,1)
-    plt.plot(t,yu1[0],label='$y_1$')
-    plt.plot(t,yu1[1],label='$y_2$')
+    plt.subplot(1, 3, 1)
+    plt.plot(t, yu1[0], label='$y_1$')
+    plt.plot(t, yu1[1], label='$y_2$')
     plt.xlabel('time')
     plt.ylabel('output')
-    plt.ylim([-1.1,2.1])
+    plt.ylim([-1.1, 2.1])
     plt.legend()
     plt.title('o/l response\nto input [1,0]')
 
-    plt.subplot(1,3,2)
-    plt.plot(t,yu2[0],label='$y_1$')
-    plt.plot(t,yu2[1],label='$y_2$')
+    plt.subplot(1, 3, 2)
+    plt.plot(t, yu2[0], label='$y_1$')
+    plt.plot(t, yu2[1], label='$y_2$')
     plt.xlabel('time')
     plt.ylabel('output')
-    plt.ylim([-1.1,2.1])
+    plt.ylim([-1.1, 2.1])
     plt.legend()
     plt.title('o/l response\nto input [0,1]')
 
-    plt.subplot(1,3,3)
-    plt.plot(t,yuz[0],label='$y_1$')
-    plt.plot(t,yuz[1],label='$y_2$')
+    plt.subplot(1, 3, 3)
+    plt.plot(t, yuz[0], label='$y_1$')
+    plt.plot(t, yuz[1], label='$y_2$')
     plt.xlabel('time')
     plt.ylabel('output')
-    plt.ylim([-1.1,2.1])
+    plt.ylim([-1.1, 2.1])
     plt.legend()
     plt.title('o/l response\nto input [1,-1]')
 
 
-def synth(wb1,wb2):
+def synth(wb1, wb2):
     """synth(wb1,wb2) -> k,gamma
     wb1: S weighting frequency
     wb2: KS weighting frequency
@@ -100,30 +91,30 @@ def synth(wb1,wb2):
     with loop closed through design
     """
     g = plant()
-    wu = ss([],[],[],np.eye(2))
-    wp1 = ss(weighting(wb=wb1, m=1.5, a=1e-4))
-    wp2 = ss(weighting(wb=wb2, m=1.5, a=1e-4))
+    wu = ss([], [], [], np.eye(2))
+    wp1 = ss(1/makeweight(wc=wb1, hfgain=1.5, dcgain=1e-4))
+    wp2 = ss(1/makeweight(wc=wb2, hfgain=1.5, dcgain=1e-4))
     wp = wp1.append(wp2)
-    k,_,info = mixsyn(g,wp,wu)
+    k, _, info = mixsyn(g, wp, wu)
     return k, info[0]
 
 
-def step_opposite(g,t):
+def step_opposite(g, t):
     """reponse to step of [-1,1]"""
-    _, yu1 = step_response(g,t,input=0)
-    _, yu2 = step_response(g,t,input=1)
+    _, yu1 = step_response(g, t, input=0)
+    _, yu2 = step_response(g, t, input=1)
     return yu1 - yu2
 
 
 def design():
     """Show results of designs"""
     # equal weighting on each output
-    k1, gam1 = synth(0.25,0.25)
+    k1, gam1 = synth(0.25, 0.25)
     # increase "bandwidth" of output 2 by moving crossover weighting frequency 100 times higher
-    k2, gam2 = synth(0.25,25)
+    k2, gam2 = synth(0.25, 25)
     # now weight output 1 more heavily
     # won't plot this one, just want gamma
-    _, gam3 = synth(25,0.25)
+    _, gam3 = synth(25, 0.25)
 
     print('design 1 gamma {:.3g} (Skogestad: 2.80)'.format(gam1))
     print('design 2 gamma {:.3g} (Skogestad: 2.92)'.format(gam2))
@@ -131,42 +122,42 @@ def design():
 
     # do the designs
     g = plant()
-    w = np.logspace(-2,2,101)
-    I = ss([],[],[],np.eye(2))
-    s1 = I.feedback(g*k1)
-    s2 = I.feedback(g*k2)
+    w = np.logspace(-2, 2, 101)
+    I = ss([], [], [], np.eye(2))
+    s1 = I.feedback(g * k1)
+    s2 = I.feedback(g * k2)
 
     # frequency response
-    sv1 = triv_sigma(s1,w)
-    sv2 = triv_sigma(s2,w)
-    
+    sv1 = triv_sigma(s1, w)
+    sv2 = triv_sigma(s2, w)
+
     plt.figure(2)
 
-    plt.subplot(1,2,1)
-    plt.semilogx(w, 20*np.log10(sv1[:,0]), label=r'$\sigma_1(S_1)$')
-    plt.semilogx(w, 20*np.log10(sv1[:,1]), label=r'$\sigma_2(S_1)$')
-    plt.semilogx(w, 20*np.log10(sv2[:,0]), label=r'$\sigma_1(S_2)$')
-    plt.semilogx(w, 20*np.log10(sv2[:,1]), label=r'$\sigma_2(S_2)$')
-    plt.ylim([-60,10])
+    plt.subplot(1, 2, 1)
+    plt.semilogx(w, 20 * np.log10(sv1[:, 0]), label=r'$\sigma_1(S_1)$')
+    plt.semilogx(w, 20 * np.log10(sv1[:, 1]), label=r'$\sigma_2(S_1)$')
+    plt.semilogx(w, 20 * np.log10(sv2[:, 0]), label=r'$\sigma_1(S_2)$')
+    plt.semilogx(w, 20 * np.log10(sv2[:, 1]), label=r'$\sigma_2(S_2)$')
+    plt.ylim([-60, 10])
     plt.ylabel('magnitude [dB]')
-    plt.xlim([1e-2,1e2])
+    plt.xlim([1e-2, 1e2])
     plt.xlabel('freq [rad/s]')
     plt.legend()
     plt.title('Singular values of S')
 
     # time response
- 
+
     # in design 1, both outputs have an inverse initial response; in
     # design 2, output 2 does not, and is very fast, while output 1
     # has a larger initial inverse response than in design 1
-    time = np.linspace(0,10,301)
-    t1 = (g*k1).feedback(I)
-    t2 = (g*k2).feedback(I)
-    
-    y1 = step_opposite(t1,time)
-    y2 = step_opposite(t2,time)
-
-    plt.subplot(1,2,2)
+    time = np.linspace(0, 10, 301)
+    t1 = (g * k1).feedback(I)
+    t2 = (g * k2).feedback(I)
+
+    y1 = step_opposite(t1, time)
+    y2 = step_opposite(t2, time)
+
+    plt.subplot(1, 2, 2)
     plt.plot(time, y1[0], label='des. 1 $y_1(t))$')
     plt.plot(time, y1[1], label='des. 1 $y_2(t))$')
     plt.plot(time, y2[0], label='des. 2 $y_1(t))$')
@@ -177,8 +168,8 @@ def design():
     plt.title('c/l response to reference [1,-1]')
 
 
-analysis()
-design()
-if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
-    plt.show()
-
+if __name__ == "__main__":
+    analysis()
+    design()
+    if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+        plt.show()
diff --git a/examples/robust_siso.py b/examples/robust_siso.py
index 6bc4ab327..013ea821d 100644
--- a/examples/robust_siso.py
+++ b/examples/robust_siso.py
@@ -1,4 +1,6 @@
-"""Demonstrate mixed-sensitivity H-infinity design for a SISO plant.
+"""robust_siso.py
+
+Demonstrate mixed-sensitivity H-infinity design for a SISO plant.
 
 Based on Example 2.11 from Multivariable Feedback Control, Skogestad
 and Postlethwaite, 1st Edition.
@@ -13,16 +15,16 @@
 
 s = tf([1, 0], 1)
 # the plant
-g = 200/(10*s+1)/(0.05*s+1)**2
+g = 200 / (10 * s + 1) / (0.05 * s + 1) ** 2
 # disturbance plant
-gd = 100/(10*s+1)
+gd = 100 / (10 * s + 1)
 
 # first design
 # sensitivity weighting
 M = 1.5
 wb = 10
 A = 1e-4
-ws1 = (s/M+wb)/(s+wb*A)
+ws1 = (s / M + wb) / (s + wb * A)
 # KS weighting
 wu = tf(1, 1)
 
@@ -30,72 +32,72 @@
 
 # sensitivity (S) and complementary sensitivity (T) functions for
 # design 1
-s1 = feedback(1,g*k1)
-t1 = feedback(g*k1,1)
+s1 = feedback(1, g * k1)
+t1 = feedback(g * k1, 1)
 
 # second design
 # this weighting differs from the text, where A**0.5 is used; if you use that,
 # the frequency response doesn't match the figure.  The time responses
 # are similar, though.
-ws2 = (s/M**0.5+wb)**2/(s+wb*A)**2
+ws2 = (s / M ** 0.5 + wb) ** 2 / (s + wb * A) ** 2
 # the KS weighting is the same as for the first design
 
 k2, cl2, info2 = mixsyn(g, ws2, wu)
 
 # S and T for design 2
-s2 = feedback(1,g*k2)
-t2 = feedback(g*k2,1)
+s2 = feedback(1, g * k2)
+t2 = feedback(g * k2, 1)
 
 # frequency response
-omega = np.logspace(-2,2,101)
-ws1mag,_,_ = ws1.freqresp(omega)
-s1mag,_,_ = s1.freqresp(omega)
-ws2mag,_,_ = ws2.freqresp(omega)
-s2mag,_,_ = s2.freqresp(omega)
+omega = np.logspace(-2, 2, 101)
+ws1mag, _, _ = ws1.freqresp(omega)
+s1mag, _, _ = s1.freqresp(omega)
+ws2mag, _, _ = ws2.freqresp(omega)
+s2mag, _, _ = s2.freqresp(omega)
 
 plt.figure(1)
 # text uses log-scaled absolute, but dB are probably more familiar to most control engineers
-plt.semilogx(omega,20*np.log10(s1mag.flat),label='$S_1$')
-plt.semilogx(omega,20*np.log10(s2mag.flat),label='$S_2$')
+plt.semilogx(omega, 20 * np.log10(s1mag.flat), label='$S_1$')
+plt.semilogx(omega, 20 * np.log10(s2mag.flat), label='$S_2$')
 # -1 in logspace is inverse
-plt.semilogx(omega,-20*np.log10(ws1mag.flat),label='$1/w_{P1}$')
-plt.semilogx(omega,-20*np.log10(ws2mag.flat),label='$1/w_{P2}$')
+plt.semilogx(omega, -20 * np.log10(ws1mag.flat), label='$1/w_{P1}$')
+plt.semilogx(omega, -20 * np.log10(ws2mag.flat), label='$1/w_{P2}$')
 
-plt.ylim([-80,10])
-plt.xlim([1e-2,1e2])
+plt.ylim([-80, 10])
+plt.xlim([1e-2, 1e2])
 plt.xlabel('freq [rad/s]')
 plt.ylabel('mag [dB]')
 plt.legend()
 plt.title('Sensitivity and sensitivity weighting frequency responses')
 
 # time response
-time = np.linspace(0,3,201)
-_,y1 = step_response(t1, time)
-_,y2 = step_response(t2, time)
+time = np.linspace(0, 3, 201)
+_, y1 = step_response(t1, time)
+_, y2 = step_response(t2, time)
 
 # gd injects into the output (that is, g and gd are summed), and the
 # closed loop mapping from output disturbance->output is S.
-_,y1d = step_response(s1*gd, time)
-_,y2d = step_response(s2*gd, time)
+_, y1d = step_response(s1 * gd, time)
+_, y2d = step_response(s2 * gd, time)
 
 plt.figure(2)
-plt.subplot(1,2,1)
-plt.plot(time,y1,label='$y_1(t)$')
-plt.plot(time,y2,label='$y_2(t)$')
+plt.subplot(1, 2, 1)
+plt.plot(time, y1, label='$y_1(t)$')
+plt.plot(time, y2, label='$y_2(t)$')
 
-plt.ylim([-0.1,1.5])
-plt.xlim([0,3])
+plt.ylim([-0.1, 1.5])
+plt.xlim([0, 3])
 plt.xlabel('time [s]')
 plt.ylabel('signal [1]')
 plt.legend()
 plt.title('Tracking response')
 
-plt.subplot(1,2,2)
-plt.plot(time,y1d,label='$y_1(t)$')
-plt.plot(time,y2d,label='$y_2(t)$')
+plt.subplot(1, 2, 2)
+plt.plot(time, y1d, label='$y_1(t)$')
+plt.plot(time, y2d, label='$y_2(t)$')
 
-plt.ylim([-0.1,1.5])
-plt.xlim([0,3])
+plt.ylim([-0.1, 1.5])
+plt.xlim([0, 3])
 plt.xlabel('time [s]')
 plt.ylabel('signal [1]')
 plt.legend()