BugFix: DC gain for discrete-time systems

roryyorke · roryyorke · commit 310f5805f0fd · 2017-01-04T19:42:43.000+02:00
For discrete time systems the DC gain is found at z=1; change dcgain
method in TransferFunction and StateSpace classes for this.

Added tests for static gain, low-pass filter, differencer and
summer for both classes.

If the StateSpace DC gain computation fails due to singularity, return
an array of NaNs of size (output,input), rather than scalar NaN.  Added
separate test for this for continuous- and discrete-time cases.
diff --git a/control/statesp.py b/control/statesp.py
@@ -612,11 +612,14 @@ def dcgain(self):
             at the origin
         """
         try:
-            gain = np.asarray(self.D -
-                              self.C.dot(np.linalg.solve(self.A, self.B)))
+            if self.isctime():
+                gain = np.asarray(self.D -
+                                    self.C.dot(np.linalg.solve(self.A, self.B)))
+            else:
+                gain = self.horner(1)
         except LinAlgError:
-            # zero eigenvalue: singular matrix
-            return np.nan
+            # eigenvalue at DC
+            gain = np.tile(np.nan,(self.outputs,self.inputs))
         return np.squeeze(gain)
 
 
diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
@@ -228,7 +228,8 @@ def testArrayAccessSS(self):
 
         assert sys1.dt == sys1_11.dt
 
-    def test_dcgain(self):
+    def test_dcgain_cont(self):
+        """Test DC gain for continuous-time state-space systems"""
         sys = StateSpace(-2.,6.,5.,0)
         np.testing.assert_equal(sys.dcgain(), 15.)
 
@@ -239,6 +240,43 @@ def test_dcgain(self):
         sys3 = StateSpace(0., 1., 1., 0.)
         np.testing.assert_equal(sys3.dcgain(), np.nan)
 
+    def test_dcgain_discr(self):
+        """Test DC gain for discrete-time state-space systems"""
+        # static gain
+        sys = StateSpace([], [], [], 2, True)
+        np.testing.assert_equal(sys.dcgain(), 2)
+
+        # averaging filter
+        sys = StateSpace(0.5, 0.5, 1, 0, True)
+        np.testing.assert_almost_equal(sys.dcgain(), 1)
+
+        # differencer
+        sys = StateSpace(0, 1, -1, 1, True)
+        np.testing.assert_equal(sys.dcgain(), 0)
+
+        # summer
+        sys = StateSpace(1, 1, 1, 0, True)
+        np.testing.assert_equal(sys.dcgain(), np.nan)
+
+    def test_dcgain_integrator(self):
+        """DC gain when eigenvalue at DC returns appropriately sized array of nan"""
+        # the SISO case is also tested in test_dc_gain_{cont,discr}
+        import itertools
+        # iterate over input and output sizes, and continuous (dt=None) and discrete (dt=True) time
+        for inputs,outputs,dt in itertools.product(range(1,6),range(1,6),[None,True]):
+            states = max(inputs,outputs)
+
+            # a matrix that is singular at DC, and has no "useless" states as in _remove_useless_states
+            a = np.triu(np.tile(2,(states,states)))
+            # eigenvalues all +2, except for ...
+            a[0,0] = 0 if dt is None else 1
+            b = np.eye(max(inputs,states))[:states,:inputs]
+            c = np.eye(max(outputs,states))[:outputs,:states]
+            d = np.zeros((outputs,inputs))
+            sys = StateSpace(a,b,c,d,dt)
+            dc = np.squeeze(np.tile(np.nan,(outputs,inputs)))
+            np.testing.assert_array_equal(dc, sys.dcgain())
+
 
     def test_scalarStaticGain(self):
         """Regression: can we create a scalar static gain?"""
diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
@@ -524,7 +524,8 @@ def testMatrixMult(self):
         np.testing.assert_array_almost_equal(H.num[1][0], H2.num[0][0])
         np.testing.assert_array_almost_equal(H.den[1][0], H2.den[0][0])
 
-    def test_dcgain(self):
+    def test_dcgain_cont(self):
+        """Test DC gain for continuous-time transfer functions"""
         sys = TransferFunction(6, 3)
         np.testing.assert_equal(sys.dcgain(), 2)
 
@@ -540,6 +541,26 @@ def test_dcgain(self):
         expected = [[5, 7, 11], [2, 2, 2]]
         np.testing.assert_array_equal(sys4.dcgain(), expected)
 
+    def test_dcgain_discr(self):
+        """Test DC gain for discrete-time transfer functions"""
+        # static gain
+        sys = TransferFunction(6, 3, True)
+        np.testing.assert_equal(sys.dcgain(), 2)
+
+        # averaging filter
+        sys = TransferFunction(0.5, [1, -0.5], True)
+        np.testing.assert_almost_equal(sys.dcgain(), 1)
+
+        # differencer
+        sys = TransferFunction(1, [1, -1], True)
+        np.testing.assert_equal(sys.dcgain(), np.inf)
+
+        # summer
+        # causes a RuntimeWarning due to the divide by zero
+        sys = TransferFunction([1,-1], [1], True)
+        np.testing.assert_equal(sys.dcgain(), 0)
+
+
 def suite():
     return unittest.TestLoader().loadTestsFromTestCase(TestXferFcn)
 
diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -945,13 +945,23 @@ def sample(self, Ts, method='zoh', alpha=None):
     def dcgain(self):
         """Return the zero-frequency (or DC) gain
 
-        For a transfer function G(s), the DC gain is G(0)
+        For a continous-time transfer function G(s), the DC gain is G(0)
+        For a discrete-time transfer function G(z), the DC gain is G(1)
 
         Returns
         -------
         gain : ndarray
             The zero-frequency gain
         """
+        if self.isctime():
+            return self._dcgain_cont()
+        else:
+            return self(1)
+
+    def _dcgain_cont(self):
+        """_dcgain_cont() -> DC gain as matrix or scalar
+
+        Special cased evaluation at 0 for continuous-time systems"""
         gain = np.empty((self.outputs, self.inputs), dtype=float)
         for i in range(self.outputs):
             for j in range(self.inputs):
