def simulate_linear_system(sys, u, t=None, x0=None, per_channel=False):
    '\n    Compute the linear model response to an input array sampled at given time\n    instances.\n\n    Parameters\n    ----------\n    sys : {State, Transfer}\n        The system model to be simulated\n    u : array_like\n        The real-valued input sequence to force the model. 1D arrays for single\n        input models and 2D arrays that has as many columns as the number of\n        inputs are valid inputs.\n    t : array_like, optional\n        The real-valued sequence to be used for the evolution of the system.\n        The values should be equally spaced otherwise an error is raised. For\n        discrete time models increments different than the sampling period also\n        raises an error. On the other hand for discrete models this can be\n        omitted and a time sequence will be generated automatically.\n    x0 : array_like, optional\n        The initial condition array. If omitted an array of zeros is assumed.\n        Note that Transfer models by definition assume zero initial conditions\n        and will raise an error.\n    per_channel : bool, optional\n        If this is set to True and if the system has multiple inputs, the\n        response of each input is returned individually. For example, if a\n        system has 4 inputs and 3 outputs then the response shape becomes\n        (num, p, m) instead of (num, p) where k-th slice (:, :, k) is the\n        response from the k-th input channel. For single input systems, this\n        keyword has no effect.\n\n    Returns\n    -------\n    yout : ndarray\n        The resulting response array. The array is 1D if sys is SISO and\n        has p columns if sys has p outputs.\n    tout : ndarray\n        The time sequence used in the simulation. If the parameter t is not\n        None then a copy of t is given.\n\n    Notes\n    -----\n    For Transfer models, first conversion to a state model is performed and\n    then the resulting model is used for computations.\n\n    '
    _check_for_state_or_transfer(sys)
    if (x0 is not None):
        if sys._isgain:
            raise ValueError("Static system models can't have initial conditions set.")
        if isinstance(sys, Transfer):
            raise ValueError("Transfer models can't have initial conditions set.")
        x0 = np.asarray(x0, dtype=float).squeeze()
        if (x0.ndim > 1):
            raise ValueError('Initial condition can only be a 1D array.')
        else:
            x0 = x0[:, None]
        if (sys.NumberOfStates != x0.size):
            raise ValueError('The initial condition size does not match the number of states of the model.')
    try:
        _check_for_state(sys)
    except ValueError:
        sys = transfer_to_state(sys)
    (n, m) = (sys.NumberOfStates, sys.shape[1])
    is_discrete = (sys.SamplingSet == 'Z')
    u = np.asarray(u, dtype=float).squeeze()
    if (u.ndim == 1):
        u = u[:, None]
    t = _check_u_and_t_for_simulation(m, sys._dt, u, t, is_discrete)
    if sys._isgain:
        if sys._isSISO:
            yout = (u * sys.d.squeeze())
        else:
            if (m == 1):
                per_channel = False
            if per_channel:
                yout = np.einsum('ij,jk->ikj', u, sys.d.T)
            else:
                yout = (u @ sys.d.T)
    else:
        if (not is_discrete):
            sys = discretize(sys, (t[1] - t[0]), method='zoh')
        sample_num = len(u)
        (a, b, c, d) = sys.matrices
        M_u = np.block([b.T, d.T])
        at = a.T
        if (m == 1):
            per_channel = False
        if per_channel:
            xout = np.empty([sample_num, n, m], dtype=float)
            for col in range(m):
                xout[0, :, col] = (0.0 if (x0 is None) else x0.T)
                Bu = (u[:, [col]] @ b.T[[col], :])
                for row in range(1, sample_num):
                    xout[row, :, col] = ((xout[(row - 1), :, col] @ at) + Bu[(row - 1)])
            yout = (np.einsum('ijk,jl->ilk', xout, c.T) + np.einsum('ij,jk->ikj', u, d.T))
        else:
            BDu = (u @ M_u)
            xout = np.empty([sample_num, n], dtype=float)
            xout[0] = (0.0 if (x0 is None) else x0.T)
            for row in range(1, sample_num):
                xout[row] = ((xout[(row - 1)] @ at) + BDu[(row - 1), :n])
            yout = ((xout @ c.T) + BDu[:, n:])
    return (yout, t)
Compute the linear model response to an input array sampled at given time
instances.

Parameters
----------
sys : {State, Transfer}
    The system model to be simulated
u : array_like
    The real-valued input sequence to force the model. 1D arrays for single
    input models and 2D arrays that has as many columns as the number of
    inputs are valid inputs.
t : array_like, optional
    The real-valued sequence to be used for the evolution of the system.
    The values should be equally spaced otherwise an error is raised. For
    discrete time models increments different than the sampling period also
    raises an error. On the other hand for discrete models this can be
    omitted and a time sequence will be generated automatically.
x0 : array_like, optional
    The initial condition array. If omitted an array of zeros is assumed.
    Note that Transfer models by definition assume zero initial conditions
    and will raise an error.
per_channel : bool, optional
    If this is set to True and if the system has multiple inputs, the
    response of each input is returned individually. For example, if a
    system has 4 inputs and 3 outputs then the response shape becomes
    (num, p, m) instead of (num, p) where k-th slice (:, :, k) is the
    response from the k-th input channel. For single input systems, this
    keyword has no effect.

Returns
-------
yout : ndarray
    The resulting response array. The array is 1D if sys is SISO and
    has p columns if sys has p outputs.
tout : ndarray
    The time sequence used in the simulation. If the parameter t is not
    None then a copy of t is given.

Notes
-----
For Transfer models, first conversion to a state model is performed and
then the resulting model is used for computations.
def simulate_linear_system(sys, u, t=None, x0=None, per_channel=False):
    '\n    Compute the linear model response to an input array sampled at given time\n    instances.\n\n    Parameters\n    ----------\n    sys : {State, Transfer}\n        The system model to be simulated\n    u : array_like\n        The real-valued input sequence to force the model. 1D arrays for single\n        input models and 2D arrays that has as many columns as the number of\n        inputs are valid inputs.\n    t : array_like, optional\n        The real-valued sequence to be used for the evolution of the system.\n        The values should be equally spaced otherwise an error is raised. For\n        discrete time models increments different than the sampling period also\n        raises an error. On the other hand for discrete models this can be\n        omitted and a time sequence will be generated automatically.\n    x0 : array_like, optional\n        The initial condition array. If omitted an array of zeros is assumed.\n        Note that Transfer models by definition assume zero initial conditions\n        and will raise an error.\n    per_channel : bool, optional\n        If this is set to True and if the system has multiple inputs, the\n        response of each input is returned individually. For example, if a\n        system has 4 inputs and 3 outputs then the response shape becomes\n        (num, p, m) instead of (num, p) where k-th slice (:, :, k) is the\n        response from the k-th input channel. For single input systems, this\n        keyword has no effect.\n\n    Returns\n    -------\n    yout : ndarray\n        The resulting response array. The array is 1D if sys is SISO and\n        has p columns if sys has p outputs.\n    tout : ndarray\n        The time sequence used in the simulation. If the parameter t is not\n        None then a copy of t is given.\n\n    Notes\n    -----\n    For Transfer models, first conversion to a state model is performed and\n    then the resulting model is used for computations.\n\n    '
    _check_for_state_or_transfer(sys)
    if (x0 is not None):
        if sys._isgain:
            raise ValueError("Static system models can't have initial conditions set.")
        if isinstance(sys, Transfer):
            raise ValueError("Transfer models can't have initial conditions set.")
        x0 = np.asarray(x0, dtype=float).squeeze()
        if (x0.ndim > 1):
            raise ValueError('Initial condition can only be a 1D array.')
        else:
            x0 = x0[:, None]
        if (sys.NumberOfStates != x0.size):
            raise ValueError('The initial condition size does not match the number of states of the model.')
    try:
        _check_for_state(sys)
    except ValueError:
        sys = transfer_to_state(sys)
    (n, m) = (sys.NumberOfStates, sys.shape[1])
    is_discrete = (sys.SamplingSet == 'Z')
    u = np.asarray(u, dtype=float).squeeze()
    if (u.ndim == 1):
        u = u[:, None]
    t = _check_u_and_t_for_simulation(m, sys._dt, u, t, is_discrete)
    if sys._isgain:
        if sys._isSISO:
            yout = (u * sys.d.squeeze())
        else:
            if (m == 1):
                per_channel = False
            if per_channel:
                yout = np.einsum('ij,jk->ikj', u, sys.d.T)
            else:
                yout = (u @ sys.d.T)
    else:
        if (not is_discrete):
            sys = discretize(sys, (t[1] - t[0]), method='zoh')
        sample_num = len(u)
        (a, b, c, d) = sys.matrices
        M_u = np.block([b.T, d.T])
        at = a.T
        if (m == 1):
            per_channel = False
        if per_channel:
            xout = np.empty([sample_num, n, m], dtype=float)
            for col in range(m):
                xout[0, :, col] = (0.0 if (x0 is None) else x0.T)
                Bu = (u[:, [col]] @ b.T[[col], :])
                for row in range(1, sample_num):
                    xout[row, :, col] = ((xout[(row - 1), :, col] @ at) + Bu[(row - 1)])
            yout = (np.einsum('ijk,jl->ilk', xout, c.T) + np.einsum('ij,jk->ikj', u, d.T))
        else:
            BDu = (u @ M_u)
            xout = np.empty([sample_num, n], dtype=float)
            xout[0] = (0.0 if (x0 is None) else x0.T)
            for row in range(1, sample_num):
                xout[row] = ((xout[(row - 1)] @ at) + BDu[(row - 1), :n])
            yout = ((xout @ c.T) + BDu[:, n:])
    return (yout, t)