Python scipy.signal.cont2discrete() Examples

def test_backward_diff(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5

        ad_truth = 2.0 * np.eye(2)
        bd_truth = 0.5 * np.ones((2, 1))
        cd_truth = np.array([[1.5, 2.0],
                             [2.0, 2.0],
                             [2.0, 0.5]])
        dd_truth = np.array([[0.875],
                             [1.0],
                             [0.295]])

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='backward_diff')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd) 

def cont2disc(self, dt=None):
        """Convert continuous-time SS model into """

        assert self.discr_method is not 'newmark', \
            'For Newmark-beta discretisation, use assemble method directly.'

        if dt is not None:
            self.dt = dt
        else:
            assert self.dt is not None, \
                'Provide time-step for convertion to discrete-time'

        SScont = self.SScont
        tpl = scsig.cont2discrete(
            (SScont.A, SScont.B, SScont.C, SScont.D),
            dt=self.dt, method=self.discr_method)
        self.SSdisc = libss.ss(*tpl[:-1], dt=tpl[-1])
        self.dlti = True 

def linearize(self):
        """Return the discretized, scaled, linearized system.

        Returns
        -------
        Ad : ndarray
            The discrete-time state matrix.

        """
        A = np.array([[0, -1], [1, -1]], dtype=config.np_dtype)
        if self.normalization is not None:
            Tx = np.diag(self.normalization)
            Tx_inv = np.diag(self.inv_norm)
            A = np.linalg.multi_dot((Tx_inv, A, Tx))
        B = np.zeros([2, 1])
        Ad, _, _, _, _ = signal.cont2discrete((A, B, 0, 0), self.dt, method='zoh')
        return Ad 

def test_multioutput(self):
        ts = 0.01  # time step

        tf = ([[1, -3], [1, 5]], [1, 1])
        num, den, dt = c2d(tf, ts)

        tf1 = (tf[0][0], tf[1])
        num1, den1, dt1 = c2d(tf1, ts)

        tf2 = (tf[0][1], tf[1])
        num2, den2, dt2 = c2d(tf2, ts)

        # Sanity checks
        assert_equal(dt, dt1)
        assert_equal(dt, dt2)

        # Check that we get the same results
        assert_allclose(num, np.vstack((num1, num2)), rtol=1e-13)

        # Single input, so the denominator should
        # not be multidimensional like the numerator
        assert_allclose(den, den1, rtol=1e-13)
        assert_allclose(den, den2, rtol=1e-13) 

def test_zerospolesgain(self):
        zeros_c = np.array([0.5, -0.5])
        poles_c = np.array([1.j / np.sqrt(2), -1.j / np.sqrt(2)])
        k_c = 1.0

        zeros_d = [1.23371727305860, 0.735356894461267]
        polls_d = [0.938148335039729 + 0.346233593780536j,
                   0.938148335039729 - 0.346233593780536j]
        k_d = 1.0

        dt_requested = 0.5

        zeros, poles, k, dt = c2d((zeros_c, poles_c, k_c), dt_requested,
                                  method='zoh')

        assert_array_almost_equal(zeros_d, zeros)
        assert_array_almost_equal(polls_d, poles)
        assert_almost_equal(k_d, k)
        assert_almost_equal(dt_requested, dt) 

def test_euler(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5

        ad_truth = 1.5 * np.eye(2)
        bd_truth = 0.25 * np.ones((2, 1))
        cd_truth = np.array([[0.75, 1.0],
                             [1.0, 1.0],
                             [1.0, 0.25]])
        dd_truth = dc

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='euler')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd)
        assert_almost_equal(dt_requested, dt) 

def test_gbt(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5
        alpha = 1.0 / 3.0

        ad_truth = 1.6 * np.eye(2)
        bd_truth = 0.3 * np.ones((2, 1))
        cd_truth = np.array([[0.9, 1.2],
                             [1.2, 1.2],
                             [1.2, 0.3]])
        dd_truth = np.array([[0.175],
                             [0.2],
                             [-0.205]])

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='gbt', alpha=alpha)

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd) 

def test_zoh(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        ad_truth = 1.648721270700128 * np.eye(2)
        bd_truth = 0.324360635350064 * np.ones((2, 1))
        # c and d in discrete should be equal to their continuous counterparts
        dt_requested = 0.5

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested, method='zoh')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cc, cd)
        assert_array_almost_equal(dc, dd)
        assert_almost_equal(dt_requested, dt) 

def test_zoh(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        ad_truth = 1.648721270700128 * np.eye(2)
        bd_truth = 0.324360635350064 * np.ones((2, 1))
        # c and d in discrete should be equal to their continuous counterparts
        dt_requested = 0.5

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested, method='zoh')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cc, cd)
        assert_array_almost_equal(dc, dd)
        assert_almost_equal(dt_requested, dt) 

def test_zerospolesgain(self):
        zeros_c = np.array([0.5, -0.5])
        poles_c = np.array([1.j / np.sqrt(2), -1.j / np.sqrt(2)])
        k_c = 1.0

        zeros_d = [1.23371727305860, 0.735356894461267]
        polls_d = [0.938148335039729 + 0.346233593780536j,
                   0.938148335039729 - 0.346233593780536j]
        k_d = 1.0

        dt_requested = 0.5

        zeros, poles, k, dt = c2d((zeros_c, poles_c, k_c), dt_requested,
                                  method='zoh')

        assert_array_almost_equal(zeros_d, zeros)
        assert_array_almost_equal(polls_d, poles)
        assert_almost_equal(k_d, k)
        assert_almost_equal(dt_requested, dt) 

def test_backward_diff(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5

        ad_truth = 2.0 * np.eye(2)
        bd_truth = 0.5 * np.ones((2, 1))
        cd_truth = np.array([[1.5, 2.0],
                             [2.0, 2.0],
                             [2.0, 0.5]])
        dd_truth = np.array([[0.875],
                             [1.0],
                             [0.295]])

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='backward_diff')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd) 

def test_gbt(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5
        alpha = 1.0 / 3.0

        ad_truth = 1.6 * np.eye(2)
        bd_truth = 0.3 * np.ones((2, 1))
        cd_truth = np.array([[0.9, 1.2],
                             [1.2, 1.2],
                             [1.2, 0.3]])
        dd_truth = np.array([[0.175],
                             [0.2],
                             [-0.205]])

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='gbt', alpha=alpha)

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd) 

def test_euler(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5

        ad_truth = 1.5 * np.eye(2)
        bd_truth = 0.25 * np.ones((2, 1))
        cd_truth = np.array([[0.75, 1.0],
                             [1.0, 1.0],
                             [1.0, 0.25]])
        dd_truth = dc

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='euler')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd)
        assert_almost_equal(dt_requested, dt) 

def test_discrete_approx(self):
        """
        Test that the solution to the discrete approximation of a continuous
        system actually approximates the solution to the continuous sytem.
        This is an indirect test of the correctness of the implementation
        of cont2discrete.
        """

        def u(t):
            return np.sin(2.5 * t)

        a = np.array([[-0.01]])
        b = np.array([[1.0]])
        c = np.array([[1.0]])
        d = np.array([[0.2]])
        x0 = 1.0

        t = np.linspace(0, 10.0, 101)
        dt = t[1] - t[0]
        u1 = u(t)

        # Use lsim2 to compute the solution to the continuous system.
        t, yout, xout = lsim2((a, b, c, d), T=t, U=u1, X0=x0,
                              rtol=1e-9, atol=1e-11)

        # Convert the continuous system to a discrete approximation.
        dsys = c2d((a, b, c, d), dt, method='bilinear')

        # Use dlsim with the pairwise averaged input to compute the output
        # of the discrete system.
        u2 = 0.5 * (u1[:-1] + u1[1:])
        t2 = t[:-1]
        td2, yd2, xd2 = dlsim(dsys, u=u2.reshape(-1, 1), t=t2, x0=x0)

        # ymid is the average of consecutive terms of the "exact" output
        # computed by lsim2.  This is what the discrete approximation
        # actually approximates.
        ymid = 0.5 * (yout[:-1] + yout[1:])

        assert_allclose(yd2.ravel(), ymid, rtol=1e-4) 

def linearize(self):
        """Return the discretized, scaled, linearized system.

        Returns
        -------
        Ad : ndarray
            The discrete-time state matrix.
        Bd : ndarray
            The discrete-time action matrix.

        """
        m = self.pendulum_mass
        M = self.cart_mass
        L = self.length
        b = self.rot_friction
        g = self.gravity

        A = np.array([[0, 0,                     1, 0                            ],
                      [0, 0,                     0, 1                            ],
                      [0, g * m / M,             0, -b / (M * L)                 ],
                      [0, g * (m + M) / (L * M), 0, -b * (m + M) / (m * M * L**2)]],
                     dtype=config.np_dtype)

        B = np.array([0, 0, 1 / M, 1 / (M * L)]).reshape((-1, self.action_dim))

        if self.normalization is not None:
            Tx, Tu = map(np.diag, self.normalization)
            Tx_inv, Tu_inv = map(np.diag, self.inv_norm)
            A = np.linalg.multi_dot((Tx_inv, A, Tx))
            B = np.linalg.multi_dot((Tx_inv, B, Tu))

        Ad, Bd, _, _, _ = signal.cont2discrete((A, B, 0, 0), self.dt,
                                               method='zoh')
        return Ad, Bd 

def linearize_discretize(self, x_center=None, u_center=None, normalize=True):
        """ Discretize and linearize the system around an equilibrium point

        Parameters
        ----------
        x_center: 2x0 array[float], optional
            The linearization center of the state.
            Default: the origin
        u_center: 1x0 array[float], optional
            The linearization center of the action
            Default: zero
        """
        if x_center is None:
            x_center = self.p_origin
        else:
            raise NotImplementedError(
                "For now we only allow linearization at the origin!")

        if u_center is None:
            u_center = np.zeros((self.n_s,))
        else:
            raise NotImplementedError(
                "For now we only allow linearization at the origin!")

        jac_ct = self._jac_dynamics()

        A_ct = jac_ct[:, :self.n_s]
        B_ct = jac_ct[:, self.n_s:]

        if normalize:
            m_x = np.diag(self.norm[0])
            m_u = np.diag(self.norm[1])
            m_x_inv = np.diag(self.inv_norm[0])
            m_u_inv = np.diag(self.inv_norm[1])
            A_ct = np.linalg.multi_dot((m_x_inv, A_ct, m_x))
            B_ct = np.linalg.multi_dot((m_x_inv, B_ct, m_u))

        ct_input = (A_ct, B_ct, np.eye(self.n_s), np.zeros((self.n_s, self.n_u)))
        A, B, _, _, _ = cont2discrete(ct_input, self.dt)

        return A, B 

def test_simo_tf(self):
        # See gh-5753
        tf = ([[1, 0], [1, 1]], [1, 1])
        num, den, dt = c2d(tf, 0.01)

        assert_equal(dt, 0.01)  # sanity check
        assert_allclose(den, [1, -0.990404983], rtol=1e-3)
        assert_allclose(num, [[1, -1], [1, -0.99004983]], rtol=1e-3) 

def test_discrete_approx(self):
        """
        Test that the solution to the discrete approximation of a continuous
        system actually approximates the solution to the continuous system.
        This is an indirect test of the correctness of the implementation
        of cont2discrete.
        """

        def u(t):
            return np.sin(2.5 * t)

        a = np.array([[-0.01]])
        b = np.array([[1.0]])
        c = np.array([[1.0]])
        d = np.array([[0.2]])
        x0 = 1.0

        t = np.linspace(0, 10.0, 101)
        dt = t[1] - t[0]
        u1 = u(t)

        # Use lsim2 to compute the solution to the continuous system.
        t, yout, xout = lsim2((a, b, c, d), T=t, U=u1, X0=x0,
                              rtol=1e-9, atol=1e-11)

        # Convert the continuous system to a discrete approximation.
        dsys = c2d((a, b, c, d), dt, method='bilinear')

        # Use dlsim with the pairwise averaged input to compute the output
        # of the discrete system.
        u2 = 0.5 * (u1[:-1] + u1[1:])
        t2 = t[:-1]
        td2, yd2, xd2 = dlsim(dsys, u=u2.reshape(-1, 1), t=t2, x0=x0)

        # ymid is the average of consecutive terms of the "exact" output
        # computed by lsim2.  This is what the discrete approximation
        # actually approximates.
        ymid = 0.5 * (yout[:-1] + yout[1:])

        assert_allclose(yd2.ravel(), ymid, rtol=1e-4) 

def test_gbt_with_sio_tf_and_zpk(self):
        """Test method='gbt' with alpha=0.25 for tf and zpk cases."""
        # State space coefficients for the continuous SIO system.
        A = -1.0
        B = 1.0
        C = 1.0
        D = 0.5

        # The continuous transfer function coefficients.
        cnum, cden = ss2tf(A, B, C, D)

        # Continuous zpk representation
        cz, cp, ck = ss2zpk(A, B, C, D)

        h = 1.0
        alpha = 0.25

        # Explicit formulas, in the scalar case.
        Ad = (1 + (1 - alpha) * h * A) / (1 - alpha * h * A)
        Bd = h * B / (1 - alpha * h * A)
        Cd = C / (1 - alpha * h * A)
        Dd = D + alpha * C * Bd

        # Convert the explicit solution to tf
        dnum, dden = ss2tf(Ad, Bd, Cd, Dd)

        # Compute the discrete tf using cont2discrete.
        c2dnum, c2dden, dt = c2d((cnum, cden), h, method='gbt', alpha=alpha)

        assert_allclose(dnum, c2dnum)
        assert_allclose(dden, c2dden)

        # Convert explicit solution to zpk.
        dz, dp, dk = ss2zpk(Ad, Bd, Cd, Dd)

        # Compute the discrete zpk using cont2discrete.
        c2dz, c2dp, c2dk, dt = c2d((cz, cp, ck), h, method='gbt', alpha=alpha)

        assert_allclose(dz, c2dz)
        assert_allclose(dp, c2dp)
        assert_allclose(dk, c2dk) 

def test_gbt_with_sio_tf_and_zpk(self):
        """Test method='gbt' with alpha=0.25 for tf and zpk cases."""
        # State space coefficients for the continuous SIO system.
        A = -1.0
        B = 1.0
        C = 1.0
        D = 0.5

        # The continuous transfer function coefficients.
        cnum, cden = ss2tf(A, B, C, D)

        # Continuous zpk representation
        cz, cp, ck = ss2zpk(A, B, C, D)

        h = 1.0
        alpha = 0.25

        # Explicit formulas, in the scalar case.
        Ad = (1 + (1 - alpha) * h * A) / (1 - alpha * h * A)
        Bd = h * B / (1 - alpha * h * A)
        Cd = C / (1 - alpha * h * A)
        Dd = D + alpha * C * Bd

        # Convert the explicit solution to tf
        dnum, dden = ss2tf(Ad, Bd, Cd, Dd)

        # Compute the discrete tf using cont2discrete.
        c2dnum, c2dden, dt = c2d((cnum, cden), h, method='gbt', alpha=alpha)

        assert_allclose(dnum, c2dnum)
        assert_allclose(dden, c2dden)

        # Convert explicit solution to zpk.
        dz, dp, dk = ss2zpk(Ad, Bd, Cd, Dd)

        # Compute the discrete zpk using cont2discrete.
        c2dz, c2dp, c2dk, dt = c2d((cz, cp, ck), h, method='gbt', alpha=alpha)

        assert_allclose(dz, c2dz)
        assert_allclose(dp, c2dp)
        assert_allclose(dk, c2dk) 

def test_bilinear(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5

        ad_truth = (5.0 / 3.0) * np.eye(2)
        bd_truth = (1.0 / 3.0) * np.ones((2, 1))
        cd_truth = np.array([[1.0, 4.0 / 3.0],
                             [4.0 / 3.0, 4.0 / 3.0],
                             [4.0 / 3.0, 1.0 / 3.0]])
        dd_truth = np.array([[0.291666666666667],
                             [1.0 / 3.0],
                             [-0.121666666666667]])

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='bilinear')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd)
        assert_almost_equal(dt_requested, dt)

        # Same continuous system again, but change sampling rate

        ad_truth = 1.4 * np.eye(2)
        bd_truth = 0.2 * np.ones((2, 1))
        cd_truth = np.array([[0.9, 1.2], [1.2, 1.2], [1.2, 0.3]])
        dd_truth = np.array([[0.175], [0.2], [-0.205]])

        dt_requested = 1.0 / 3.0

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='bilinear')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd)
        assert_almost_equal(dt_requested, dt) 

def sample(self, Ts, method='zoh', alpha=None):
        """Convert a continuous time system to discrete time

        Creates a discrete-time system from a continuous-time system by
        sampling.  Multiple methods of conversion are supported.

        Parameters
        ----------
        Ts : float
            Sampling period
        method :  {"gbt", "bilinear", "euler", "backward_diff", "zoh"}
            Which method to use:

            * gbt: generalized bilinear transformation
            * bilinear: Tustin's approximation ("gbt" with alpha=0.5)
            * euler: Euler (or forward differencing) method ("gbt" with
              alpha=0)
            * backward_diff: Backwards differencing ("gbt" with alpha=1.0)
            * zoh: zero-order hold (default)

        alpha : float within [0, 1]
            The generalized bilinear transformation weighting parameter, which
            should only be specified with method="gbt", and is ignored
            otherwise

        Returns
        -------
        sysd : StateSpace
            Discrete time system, with sampling rate Ts

        Notes
        -----
        Uses the command 'cont2discrete' from scipy.signal

        Examples
        --------
        >>> sys = StateSpace(0, 1, 1, 0)
        >>> sysd = sys.sample(0.5, method='bilinear')

        """
        if not self.isctime():
            raise ValueError("System must be continuous time system")

        sys = (self.A, self.B, self.C, self.D)
        Ad, Bd, C, D, dt = cont2discrete(sys, Ts, method, alpha)
        return StateSpace(Ad, Bd, C, D, dt) 

def test_bilinear(self):
        ac = np.eye(2)
        bc = 0.5 * np.ones((2, 1))
        cc = np.array([[0.75, 1.0], [1.0, 1.0], [1.0, 0.25]])
        dc = np.array([[0.0], [0.0], [-0.33]])

        dt_requested = 0.5

        ad_truth = (5.0 / 3.0) * np.eye(2)
        bd_truth = (1.0 / 3.0) * np.ones((2, 1))
        cd_truth = np.array([[1.0, 4.0 / 3.0],
                             [4.0 / 3.0, 4.0 / 3.0],
                             [4.0 / 3.0, 1.0 / 3.0]])
        dd_truth = np.array([[0.291666666666667],
                             [1.0 / 3.0],
                             [-0.121666666666667]])

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='bilinear')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd)
        assert_almost_equal(dt_requested, dt)

        # Same continuous system again, but change sampling rate

        ad_truth = 1.4 * np.eye(2)
        bd_truth = 0.2 * np.ones((2, 1))
        cd_truth = np.array([[0.9, 1.2], [1.2, 1.2], [1.2, 0.3]])
        dd_truth = np.array([[0.175], [0.2], [-0.205]])

        dt_requested = 1.0 / 3.0

        ad, bd, cd, dd, dt = c2d((ac, bc, cc, dc), dt_requested,
                                 method='bilinear')

        assert_array_almost_equal(ad_truth, ad)
        assert_array_almost_equal(bd_truth, bd)
        assert_array_almost_equal(cd_truth, cd)
        assert_array_almost_equal(dd_truth, dd)
        assert_almost_equal(dt_requested, dt) 

def control_systems(request):
    ct_sys, ref = request.param
    Ac, Bc, Cc = ct_sys.data
    Dc = np.zeros((Cc.shape[0], 1))

    Q = np.eye(Ac.shape[0])
    R = np.eye(Bc.shape[1] if len(Bc.shape) > 1 else 1)

    Sc = linalg.solve_continuous_are(Ac, Bc.reshape(-1, 1), Q, R,)
    Kc = linalg.solve(R, Bc.T @ Sc).reshape(1, -1)
    ct_ctr = LTISystem(Kc)

    evals = np.sort(np.abs(
        linalg.eig(Ac, left=False, right=False, check_finite=False)
    ))
    dT = 1/(2*evals[-1])

    Tsim = (8/np.min(evals[~np.isclose(evals, 0)])
            if np.sum(np.isclose(evals[np.nonzero(evals)], 0)) > 0
            else 8
            )

    dt_data = signal.cont2discrete((Ac, Bc.reshape(-1, 1), Cc, Dc), dT)
    Ad, Bd, Cd, Dd = dt_data[:-1]
    Sd = linalg.solve_discrete_are(Ad, Bd.reshape(-1, 1), Q, R,)
    Kd = linalg.solve(Bd.T @ Sd @ Bd + R, Bd.T @ Sd @ Ad)

    dt_sys = LTISystem(Ad, Bd, dt=dT)
    dt_sys.initial_condition = ct_sys.initial_condition
    dt_ctr = LTISystem(Kd, dt=dT)

    yield ct_sys, ct_ctr, dt_sys, dt_ctr, ref, Tsim 

def test_mixed_dts(simulation_results):
    """
    A DT equivalent that is sampled twice as fast should match original DT
    equivalent exactly at t= k*dT under the same inputs.
    """
    results, ct_sys, ct_ctr, dt_sys, dt_ctr, ref, Tsim, tspan, intname = \
        simulation_results
    Ac, Bc, Cc = ct_sys.data
    Dc = np.zeros((Cc.shape[0], 1))

    dt_unique_t, dt_unique_t_idx = np.unique(
        results[0].t, return_index=True
    )
    discrete_sel = dt_unique_t_idx[
        (dt_unique_t < (Tsim*7/8)) & (dt_unique_t % dt_sys.dt == 0)
    ]

    scale_factor = 1/2
    Ad, Bd, Cd, Dd, dT = signal.cont2discrete(
        (Ac, Bc.reshape(-1, 1), Cc, Dc),
        dt_sys.dt*scale_factor
    )
    dt_sys2 = LTISystem(Ad, Bd, dt=dT)
    dt_sys2.initial_condition = ct_sys.initial_condition

    intopts = block_diagram.DEFAULT_INTEGRATOR_OPTIONS.copy()
    intopts['name'] = intname

    bd = BlockDiagram(dt_sys2, ref, dt_ctr)
    bd.connect(dt_sys2, ref)
    bd.connect(ref, dt_ctr)
    bd.connect(dt_ctr, dt_sys2)
    res = bd.simulate(tspan, integrator_options=intopts)

    mixed_t_discrete_t_equal_idx = np.where(
        np.equal(*np.meshgrid(res.t, results[0].t[discrete_sel]))
    )[1]

    mixed_unique_t, mixed_unique_t_idx_rev = np.unique(
        res.t[mixed_t_discrete_t_equal_idx][::-1], return_index=True
    )
    mixed_unique_t_idx = (len(mixed_t_discrete_t_equal_idx)
                          - mixed_unique_t_idx_rev - 1)
    mixed_sel = mixed_t_discrete_t_equal_idx[mixed_unique_t_idx]

    npt.assert_allclose(
        res.x[mixed_sel, :], results[0].x[discrete_sel, :],
        atol=TEST_ATOL, rtol=TEST_RTOL
    ) 

def sample(self, Ts, method='zoh', alpha=None):
        """Convert a continuous-time system to discrete time

        Creates a discrete-time system from a continuous-time system by
        sampling.  Multiple methods of conversion are supported.

        Parameters
        ----------
        Ts : float
            Sampling period
        method : {"gbt", "bilinear", "euler", "backward_diff",
                  "zoh", "matched"}
            Method to use for sampling:

            * gbt: generalized bilinear transformation
            * bilinear: Tustin's approximation ("gbt" with alpha=0.5)
            * euler: Euler (or forward difference) method ("gbt" with alpha=0)
            * backward_diff: Backwards difference ("gbt" with alpha=1.0)
            * zoh: zero-order hold (default)

        alpha : float within [0, 1]
            The generalized bilinear transformation weighting parameter, which
            should only be specified with method="gbt", and is ignored
            otherwise.

        Returns
        -------
        sysd : StateSpace system
            Discrete time system, with sampling rate Ts

        Notes
        -----
        1. Available only for SISO systems

        2. Uses the command `cont2discrete` from `scipy.signal`

        Examples
        --------
        >>> sys = TransferFunction(1, [1,1])
        >>> sysd = sys.sample(0.5, method='bilinear')

        """
        if not self.isctime():
            raise ValueError("System must be continuous time system")
        if not self.issiso():
            raise NotImplementedError("MIMO implementation not available")
        if method == "matched":
            return _c2d_matched(self, Ts)
        sys = (self.num[0][0], self.den[0][0])
        numd, dend, dt = cont2discrete(sys, Ts, method, alpha)
        return TransferFunction(numd[0, :], dend, dt)