Repository: python-control/python-control
Pull request #1074 (2024-12-04 20:28:51)

update model_reduction to allow input/output selection + unstable warning (vs exception)
This PR updates the `model_reduction` function to add some additional functionality.  From the docstring:

> This function produces a reduced-order model of a system by eliminating specified inputs, outputs, and/or states from the original system.  The specific states, inputs, or outputs that are eliminated can be specified by either listing the states, inputs, or outputs to be eliminated or those to be kept.

Specific updates:
* Changed the `ELIM` positional parameter to `elim_states` (positional or keyword).
* Added `keep_states` to allow specifying the list of states to keep instead of those to eliminate.
* Added `elim_inputs`, `elim_outputs`, `keep_inputs`, `keep_outputs` to allow selection of inputs and outputs.
* Specifications of states, inputs, and outputs can be done by a list of integers (original implementation) as well as slices or lists of signal names (new).
* In the case of eliminating states, inputs, or outputs from an unstable system, a warning is issued (versus an exception) and this warning can be suppressed using `warn_unstable=False`.
* Removed the copyright info from the top of `modelsimp.py` since we know just use the `LICENSE` file.
* Add unit tests and documentation updates.
* Cleaned up import order (isort -m2) and removed extraneous spaces.

This PR was motivated by an example from FBS2e (see [steering.py](https://github.com/murrayrm/fbs2e-python/blob/main/steering.py)), where I wanted to extract the dynamics of a subsystem (and the original system had an eigenvalue at the origin).  It is possible to accomplish the same result by using the old implementation (assuming the exception is changed to a warning) and then indexing the system directly to handle the inputs and outputs, but it seemed natural to include in a single function call.



Patch:
diff --git a/control/canonical.py b/control/canonical.py
index a6204432..7be7f88a 100644
--- a/control/canonical.py
+++ b/control/canonical.py
@@ -205,7 +205,7 @@ def similarity_transform(xsys, T, timescale=1, inverse=False):
     timescale : float, optional
         If present, also rescale the time unit to tau = timescale * t
     inverse : bool, optional
-        If True (default), transform so z = T x.  If False, transform
+        If False (default), transform so z = T x.  If True, transform
         so x = T z.
 
     Returns
diff --git a/control/modelsimp.py b/control/modelsimp.py
index 966da1ce..fe519b82 100644
--- a/control/modelsimp.py
+++ b/control/modelsimp.py
@@ -1,4 +1,3 @@
-#! TODO: add module docstring
 # modelsimp.py - tools for model simplification
 #
 # Author: Steve Brunton, Kevin Chen, Lauren Padilla
@@ -6,47 +5,16 @@
 #
 # This file contains routines for obtaining reduced order models
 #
-# Copyright (c) 2010 by California Institute of Technology
-# All rights reserved.
-#
-# Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions
-# are met:
-#
-# 1. Redistributions of source code must retain the above copyright
-#    notice, this list of conditions and the following disclaimer.
-#
-# 2. Redistributions in binary form must reproduce the above copyright
-#    notice, this list of conditions and the following disclaimer in the
-#    documentation and/or other materials provided with the distribution.
-#
-# 3. Neither the name of the California Institute of Technology nor
-#    the names of its contributors may be used to endorse or promote
-#    products derived from this software without specific prior
-#    written permission.
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
-# FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CALTECH
-# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
-# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
-# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
-# SUCH DAMAGE.
-#
-# $Id$
+
+import warnings
 
 # External packages and modules
 import numpy as np
-import warnings
-from .exception import ControlSlycot, ControlArgument, ControlDimension
-from .iosys import isdtime, isctime
-from .statesp import StateSpace
+
+from .exception import ControlArgument, ControlDimension, ControlSlycot
+from .iosys import isctime, isdtime
 from .statefbk import gram
+from .statesp import StateSpace
 from .timeresp import TimeResponseData
 
 __all__ = ['hankel_singular_values', 'balanced_reduction', 'model_reduction',
@@ -108,89 +76,146 @@ def hankel_singular_values(sys):
     return hsv[::-1]
 
 
-def model_reduction(sys, ELIM, method='matchdc'):
-    """Model reduction by state elimination.
-    
-    Model reduction of `sys` by eliminating the states in `ELIM` using a given
-    method.
+def model_reduction(
+        sys, elim_states=None, method='matchdc', elim_inputs=None,
+        elim_outputs=None, keep_states=None, keep_inputs=None,
+        keep_outputs=None, warn_unstable=True):
+    """Model reduction by input, output, or state elimination.
+
+    This function produces a reduced-order model of a system by eliminating
+    specified inputs, outputs, and/or states from the original system.  The
+    specific states, inputs, or outputs that are eliminated can be
+    specified by either listing the states, inputs, or outputs to be
+    eliminated or those to be kept.
+
+    Two methods of state reduction are possible: 'truncate' removes the
+    states marked for elimination, while 'matchdc' replaces the eliminated
+    states with their equilibrium values (thereby keeping the input/output
+    gain unchanged at zero frequency ["DC"]).
 
     Parameters
     ----------
     sys : StateSpace
         Original system to reduce.
-    ELIM : array
-        Vector of states to eliminate.
+    elim_inputs, elim_outputs, elim_states : array of int or str, optional
+        Vector of inputs, outputs, or states to eliminate.  Can be specified
+        either as an offset into the appropriate vector or as a signal name.
+    keep_inputs, keep_outputs, keep_states : array, optional
+        Vector of inputs, outputs, or states to keep.  Can be specified
+        either as an offset into the appropriate vector or as a signal name.
     method : string
-        Method of removing states in `ELIM`: either 'truncate' or
-        'matchdc'.
+        Method of removing states: either 'truncate' or 'matchdc' (default).
+    warn_unstable : bool, option
+        If `False`, don't warn if system is unstable.
 
     Returns
     -------
     rsys : StateSpace
-        A reduced order model.
+        Reduced order model.
 
     Raises
     ------
     ValueError
-        Raised under the following conditions:
+        If `method` is not either 'matchdc' or 'truncate'.
+    NotImplementedError
+        If the 'matchdc' method is used for a discrete time system.
 
-            * if `method` is not either ``'matchdc'`` or ``'truncate'``
-
-            * if eigenvalues of `sys.A` are not all in left half plane
-              (`sys` must be stable)
+    Warns
+    -----
+    UserWarning
+        If eigenvalues of `sys.A` are not all stable.
 
     Examples
     --------
     >>> G = ct.rss(4)
-    >>> Gr = ct.modred(G, [0, 2], method='matchdc')
+    >>> Gr = ct.model_reduction(G, [0, 2], method='matchdc')
     >>> Gr.nstates
     2
 
-    """
+    See Also
+    --------
+    balanced_reduction : Eliminate states using Hankel singular values.
+    minimal_realization : Eliminate unreachable or unobservable states.
 
-    # Check for ss system object, need a utility for this?
+    Notes
+    -----
+    The model_reduction function issues a warning if the system has
+    unstable eigenvalues, since in those situations the stability of the
+    reduced order model may be different than the stability of the full
+    model.  No other checking is done, so users must to be careful not to
+    render a system unobservable or unreachable.
 
-    # TODO: Check for continous or discrete, only continuous supported for now
-    #   if isCont():
-    #       dico = 'C'
-    #   elif isDisc():
-    #       dico = 'D'
-    #   else:
-    if (isctime(sys)):
-        dico = 'C'
-    else:
-        raise NotImplementedError("Function not implemented in discrete time")
+    States, inputs, and outputs can be specified using integer offsets or
+    using signal names.  Slices can also be specified, but must use the
+    Python ``slice()`` function.
+
+    """
+    if not isinstance(sys, StateSpace):
+        raise TypeError("system must be a StateSpace system")
 
     # Check system is stable
-    if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
-        raise ValueError("Oops, the system is unstable!")
-
-    ELIM = np.sort(ELIM)
-    # Create list of elements not to eliminate (NELIM)
-    NELIM = [i for i in range(len(sys.A)) if i not in ELIM]
-    # A1 is a matrix of all columns of sys.A not to eliminate
-    A1 = sys.A[:, NELIM[0]].reshape(-1, 1)
-    for i in NELIM[1:]:
-        A1 = np.hstack((A1, sys.A[:, i].reshape(-1, 1)))
-    A11 = A1[NELIM, :]
-    A21 = A1[ELIM, :]
-    # A2 is a matrix of all columns of sys.A to eliminate
-    A2 = sys.A[:, ELIM[0]].reshape(-1, 1)
-    for i in ELIM[1:]:
-        A2 = np.hstack((A2, sys.A[:, i].reshape(-1, 1)))
-    A12 = A2[NELIM, :]
-    A22 = A2[ELIM, :]
-
-    C1 = sys.C[:, NELIM]
-    C2 = sys.C[:, ELIM]
-    B1 = sys.B[NELIM, :]
-    B2 = sys.B[ELIM, :]
-
-    if method == 'matchdc':
-        # if matchdc, residualize
+    if warn_unstable:
+        if isctime(sys) and np.any(np.linalg.eigvals(sys.A).real >= 0.0) or \
+           isdtime(sys) and np.any(np.abs(np.linalg.eigvals(sys.A)) >= 1):
+            warnings.warn("System is unstable; reduction may be meaningless")
+
+    # Utility function to process keep/elim keywords
+    def _process_elim_or_keep(elim, keep, labels):
+        def _expand_key(key):
+            if key is None:
+                return []
+            elif isinstance(key, str):
+                return labels.index(key)
+            elif isinstance(key, list):
+                return [_expand_key(k) for k in key]
+            elif isinstance(key, slice):
+                return range(len(labels))[key]
+            else:
+                return key
 
+        elim = np.atleast_1d(_expand_key(elim))
+        keep = np.atleast_1d(_expand_key(keep))
+            
+        if len(elim) > 0 and len(keep) > 0:
+            raise ValueError(
+                "can't provide both 'keep' and 'elim' for same variables")
+        elif len(keep) > 0:
+            keep = np.sort(keep).tolist()
+            elim = [i for i in range(len(labels)) if i not in keep]
+        else:
+            elim = [] if elim is None else np.sort(elim).tolist()
+            keep = [i for i in range(len(labels)) if i not in elim]
+        return elim, keep
+
+    # Determine which states to keep
+    elim_states, keep_states = _process_elim_or_keep(
+        elim_states, keep_states, sys.state_labels)
+    elim_inputs, keep_inputs = _process_elim_or_keep(
+        elim_inputs, keep_inputs, sys.input_labels)
+    elim_outputs, keep_outputs = _process_elim_or_keep(
+        elim_outputs, keep_outputs, sys.output_labels)
+
+    # Create submatrix of states we are keeping
+    A11 = sys.A[:, keep_states][keep_states, :]     # states we are keeping
+    A12 = sys.A[:, elim_states][keep_states, :]     # needed for 'matchdc'
+    A21 = sys.A[:, keep_states][elim_states, :]
+    A22 = sys.A[:, elim_states][elim_states, :]
+
+    B1 = sys.B[keep_states, :]
+    B2 = sys.B[elim_states, :]
+
+    C1 = sys.C[:, keep_states]
+    C2 = sys.C[:, elim_states]
+
+    # Figure out the new state space system
+    if method == 'matchdc' and A22.size > 0:
+        if sys.isdtime(strict=True):
+            raise NotImplementedError(
+                "'matchdc' not (yet) supported for discrete time systems")
+        
+        # if matchdc, residualize
         # Check if the matrix A22 is invertible
-        if np.linalg.matrix_rank(A22) != len(ELIM):
+        if np.linalg.matrix_rank(A22) != len(elim_states):
             raise ValueError("Matrix A22 is singular to working precision.")
 
         # Now precompute A22\A21 and A22\B2 (A22I = inv(A22))
@@ -206,22 +231,29 @@ def model_reduction(sys, ELIM, method='matchdc'):
         Br = B1 - A12 @ A22I_B2
         Cr = C1 - C2 @ A22I_A21
         Dr = sys.D - C2 @ A22I_B2
-    elif method == 'truncate':
-        # if truncate, simply discard state x2
+
+    elif method == 'truncate' or A22.size == 0:
+        # Get rid of unwanted states
         Ar = A11
         Br = B1
         Cr = C1
         Dr = sys.D
+
     else:
         raise ValueError("Oops, method is not supported!")
 
+    # Get rid of additional inputs and outputs
+    Br = Br[:, keep_inputs]
+    Cr = Cr[keep_outputs, :]
+    Dr = Dr[keep_outputs, :][:, keep_inputs]
+
     rsys = StateSpace(Ar, Br, Cr, Dr)
     return rsys
 
 
 def balanced_reduction(sys, orders, method='truncate', alpha=None):
     """Balanced reduced order model of sys of a given order.
-    
+
     States are eliminated based on Hankel singular value.
     If sys has unstable modes, they are removed, the
     balanced realization is done on the stable part, then
@@ -276,7 +308,7 @@ def balanced_reduction(sys, orders, method='truncate', alpha=None):
         raise ValueError("supported methods are 'truncate' or 'matchdc'")
     elif method == 'truncate':
         try:
-            from slycot import ab09md, ab09ad
+            from slycot import ab09ad, ab09md
         except ImportError:
             raise ControlSlycot(
                 "can't find slycot subroutine ab09md or ab09ad")
@@ -346,7 +378,7 @@ def balanced_reduction(sys, orders, method='truncate', alpha=None):
 
 def minimal_realization(sys, tol=None, verbose=True):
     """ Eliminate uncontrollable or unobservable states.
-    
+
     Eliminates uncontrollable or unobservable states in state-space
     models or cancelling pole-zero pairs in transfer functions. The
     output sysr has minimal order and the same response
@@ -378,14 +410,14 @@ def _block_hankel(Y, m, n):
     """Create a block Hankel matrix from impulse response."""
     q, p, _ = Y.shape
     YY = Y.transpose(0, 2, 1) # transpose for reshape
-    
+
     H = np.zeros((q*m, p*n))
-    
+
     for r in range(m):
         # shift and add row to Hankel matrix
         new_row = YY[:, r:r+n, :]
         H[q*r:q*(r+1), :] = new_row.reshape((q, p*n))
-            
+
     return H
 
 
@@ -431,7 +463,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
         unit-area impulse response of python-control. Default is True.
     transpose : bool, optional
         Assume that input data is transposed relative to the standard
-        :ref:`time-series-convention`. For TimeResponseData this parameter 
+        :ref:`time-series-convention`. For TimeResponseData this parameter
         is ignored. Default is False.
 
     Returns
@@ -466,7 +498,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
         YY = np.array(arg, ndmin=3)
         if transpose:
             YY = np.transpose(YY)
-    
+
     q, p, l = YY.shape
 
     if m is None:
@@ -476,14 +508,14 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
 
     if m*q < r or n*p < r:
         raise ValueError("Hankel parameters are to small")
-    
+
     if (l-1) < m+n:
         raise ValueError("not enough data for requested number of parameters")
-    
+
     H = _block_hankel(YY[:, :, 1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
     Hf = H[:, :-p] # first p*n columns of H
     Hl = H[:, p:] # last p*n columns of H
-    
+
     U,S,Vh = np.linalg.svd(Hf, True)
     Ur =U[:, 0:r]
     Vhr =Vh[0:r, :]
@@ -500,7 +532,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
 
 def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     """markov(Y, U, [, m])
-    
+
     Calculate the first `m` Markov parameters [D CB CAB ...] from data.
 
     This function computes the Markov parameters for a discrete time
@@ -579,7 +611,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     # Get the system description
     if len(args) < 1:
         raise ControlArgument("not enough input arguments")
-    
+
     if isinstance(args[0], TimeResponseData):
         data = args[0]
         Umat = np.array(data.inputs, ndmin=2)
@@ -639,7 +671,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     # This algorithm sets up the following problem and solves it for
     # the Markov parameters
     #
-    # (l,q)   = (l,p*m) @ (p*m,q) 
+    # (l,q)   = (l,p*m) @ (p*m,q)
     # YY      = UU @ H.T
     #
     # [ y(0)   ]   [ u(0)    0       0                 ] [ D           ]
@@ -675,7 +707,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     # Truncate first t=0 or t=m time steps, transpose the problem for lsq
     YY = Ymat[:, t:].T
     UU = UUT[:, t:].T
-    
+
     # Solve for the Markov parameters from  YY = UU @ H.T
     HT, _, _, _ = np.linalg.lstsq(UU, YY, rcond=None)
     H = HT.T/dt # scaling
diff --git a/control/pzmap.py b/control/pzmap.py
index f1d0ecae..400db280 100644
--- a/control/pzmap.py
+++ b/control/pzmap.py
@@ -79,13 +79,11 @@ class PoleZeroData:
         System name.
     sys : StateSpace or TransferFunction
         System corresponding to the data.
-    sort_loci : bool, optional
-        Set to False to turn off sorting of loci into unique branches.
 
     """
     def __init__(
             self, poles, zeros, gains=None, loci=None, dt=None, sysname=None,
-            sys=None, sort_loci=True):
+            sys=None):
         """Create a pole/zero map object.
 
         Parameters
@@ -104,14 +102,10 @@ class PoleZeroData:
             System corresponding to the data.
 
         """
-        from .rlocus import _RLSortRoots
         self.poles = poles
         self.zeros = zeros
         self.gains = gains
-        if loci is not None and sort_loci:
-            self.loci = _RLSortRoots(loci)
-        else:
-            self.loci = loci
+        self.loci = loci
         self.dt = dt
         self.sysname = sysname
         self.sys = sys
@@ -193,7 +187,7 @@ def pole_zero_plot(
     ----------
     data : List of PoleZeroData objects or LTI systems
         List of pole/zero response data objects generated by pzmap_response()
-        or root_locus_map() that are to be plotted.  If a list of systems
+        or rootlocus_response() that are to be plotted.  If a list of systems
         is given, the poles and zeros of those systems will be plotted.
     grid : bool or str, optional
         If `True` plot omega-damping grid, if `False` show imaginary axis
diff --git a/control/rlocus.py b/control/rlocus.py
index c14d1e31..54c67299 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -88,7 +88,7 @@ def root_locus_map(sysdata, gains=None):
             root_array = _RLSortRoots(root_array)
 
         responses.append(PoleZeroData(
-            sys.poles(), sys.zeros(), kvect, root_array, sort_loci=False,
+            sys.poles(), sys.zeros(), kvect, root_array,
             dt=sys.dt, sysname=sys.name, sys=sys))
 
     if isinstance(sysdata, (list, tuple)):


Test patch:
diff --git a/control/tests/docstrings_test.py b/control/tests/docstrings_test.py
index 27ced105..991ead3e 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -39,7 +39,7 @@ function_docstring_hash = {
     control.ss2tf:                      '48ff25d22d28e7b396e686dd5eb58831',
     control.tf:                         '53a13f4a7f75a31c81800e10c88730ef',
     control.tf2ss:                      '086a3692659b7321c2af126f79f4bc11',
-    control.markov:                     'eda7c4635bbb863ae6659e574285d356',
+    control.markov:                     'a4199c54cb50f07c0163d3790739eafe',
     control.gangof4:                    '0e52eb6cf7ce024f9a41f3ae3ebf04f7',
 }
 
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
index 616ef5f0..dead4eb7 100644
--- a/control/tests/modelsimp_test.py
+++ b/control/tests/modelsimp_test.py
@@ -3,11 +3,15 @@
 RMM, 30 Mar 2011 (based on TestModelSimp from v0.4a)
 """
 
+import math
+import warnings
+
 import numpy as np
 import pytest
 
+import control as ct
 from control import StateSpace, TimeResponseData, c2d, forced_response, \
-    impulse_response, step_response, rss, tf
+    impulse_response, rss, step_response, tf
 from control.exception import ControlArgument, ControlDimension
 from control.modelsimp import balred, eigensys_realization, hsvd, markov, \
     modred
@@ -42,7 +46,7 @@ class TestModelsimp:
                                     inputs=U,
                                     input_labels='u',
                                     )
-        
+
         # setup
         m = 3
         Htrue = np.array([1., 0., 0.])
@@ -101,7 +105,6 @@ class TestModelsimp:
         HT = markov(response, m)
         np.testing.assert_array_almost_equal(HT, np.transpose(Htrue))
         response.transpose=False
-        
 
         # Test example from docstring
         # TODO: There is a problem here, last markov parameter does not fit
@@ -212,7 +215,7 @@ class TestModelsimp:
         Mtrue = np.hstack([Hd.D] + [
             Hd.C @ np.linalg.matrix_power(Hd.A, i) @ Hd.B
             for i in range(m-1)])
-        
+
         Mtrue = np.squeeze(Mtrue)
 
         # Generate input/output data
@@ -238,8 +241,8 @@ class TestModelsimp:
         Mcomp_scaled = markov(response, m, dt=Ts)
 
         np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
-        np.testing.assert_allclose(Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
-        
+        np.testing.assert_allclose(
+            Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
 
     def testERASignature(self):
 
@@ -253,7 +256,7 @@ class TestModelsimp:
         B = np.array([[1.],[1.,]])
         C = np.array([[1., 0.,]])
         D = np.array([[0.,]])
-        
+
         T = np.arange(0,10,1)
         sysd_true = StateSpace(A,B,C,D,True)
         ir_true = impulse_response(sysd_true,T=T)
@@ -262,7 +265,7 @@ class TestModelsimp:
         sysd_est, _  = eigensys_realization(ir_true,r=2)
         ir_est = impulse_response(sysd_est, T=T)
         _, H_est = ir_est
-    
+
         np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
 
         # test ndarray
@@ -270,7 +273,7 @@ class TestModelsimp:
         sysd_est, _  = eigensys_realization(YY_true,r=2)
         ir_est = impulse_response(sysd_est, T=T)
         _, H_est = ir_est
-    
+
         np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
 
         # test mimo
@@ -304,10 +307,10 @@ class TestModelsimp:
         step_true = step_response(sysd_true)
         step_est = step_response(sysd_est)
 
-        np.testing.assert_allclose(step_true.outputs, 
+        np.testing.assert_allclose(step_true.outputs,
                                    step_est.outputs,
                                    rtol=1e-6, atol=1e-8)
-        
+
         # test ndarray
         _, YY_true = ir_true
         sysd_est, _  = eigensys_realization(YY_true,r=4,dt=dt)
@@ -315,7 +318,7 @@ class TestModelsimp:
         step_true = step_response(sysd_true, T=T)
         step_est = step_response(sysd_est, T=T)
 
-        np.testing.assert_allclose(step_true.outputs, 
+        np.testing.assert_allclose(step_true.outputs,
                                    step_est.outputs,
                                    rtol=1e-6, atol=1e-8)
 
@@ -343,7 +346,7 @@ class TestModelsimp:
         np.testing.assert_array_almost_equal(rsys.D, Drtrue, decimal=2)
 
     def testModredUnstable(self):
-        """Check if an error is thrown when an unstable system is given"""
+        """Check if warning is issued when an unstable system is given"""
         A = np.array(
             [[4.5418, 3.3999, 5.0342, 4.3808],
              [0.3890, 0.3599, 0.4195, 0.1760],
@@ -353,7 +356,16 @@ class TestModelsimp:
         C = np.array([[1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0]])
         D = np.array([[0.0, 0.0], [0.0, 0.0]])
         sys = StateSpace(A, B, C, D)
-        np.testing.assert_raises(ValueError, modred, sys, [2, 3])
+
+        # Make sure we get a warning message
+        with pytest.warns(UserWarning, match="System is unstable"):
+            newsys1 = modred(sys, [2, 3])
+
+        # Make sure we can turn the warning off
+        with warnings.catch_warnings():
+            warnings.simplefilter('error')
+            newsys2 = ct.model_reduction(sys, [2, 3], warn_unstable=False)
+            np.testing.assert_equal(newsys1.A, newsys2.A)
 
     def testModredTruncate(self):
         #balanced realization computed in matlab for the transfer function:
@@ -391,7 +403,7 @@ class TestModelsimp:
         B = np.array([[2.], [0.], [0.], [0.]])
         C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
         D = np.array([[0.]])
-        
+
         sys = StateSpace(A, B, C, D)
         orders = 2
         rsys = balred(sys, orders, method='truncate')
@@ -412,7 +424,7 @@ class TestModelsimp:
                 # Apply a similarity transformation
                 Ar, Br, Cr = T @ Ar @ T, T @ Br, Cr @ T
                 break
-            
+
         # Make sure we got the correct answer
         np.testing.assert_array_almost_equal(Ar, Artrue, decimal=2)
         np.testing.assert_array_almost_equal(Br, Brtrue, decimal=4)
@@ -432,12 +444,12 @@ class TestModelsimp:
         B = np.array([[2.], [0.], [0.], [0.]])
         C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
         D = np.array([[0.]])
-        
+
         sys = StateSpace(A, B, C, D)
         orders = 2
         rsys = balred(sys,orders,method='matchdc')
         Ar, Br, Cr, Dr = rsys.A, rsys.B, rsys.C, rsys.D
-        
+
         # Result from MATLAB
         Artrue = np.array(
             [[-4.43094773, -4.55232904],
@@ -445,7 +457,7 @@ class TestModelsimp:
         Brtrue = np.array([[1.36235673], [1.03114388]])
         Crtrue = np.array([[1.36235673, 1.03114388]])
         Drtrue = np.array([[-0.08383902]])
-        
+
         # Look for possible changes in state in slycot
         T1 = np.array([[1, 0], [0, -1]])
         T2 = np.array([[-1, 0], [0, 1]])
@@ -455,9 +467,46 @@ class TestModelsimp:
                 # Apply a similarity transformation
                 Ar, Br, Cr = T @ Ar @ T, T @ Br, Cr @ T
                 break
-            
+
         # Make sure we got the correct answer
         np.testing.assert_array_almost_equal(Ar, Artrue, decimal=2)
         np.testing.assert_array_almost_equal(Br, Brtrue, decimal=4)
         np.testing.assert_array_almost_equal(Cr, Crtrue, decimal=4)
         np.testing.assert_array_almost_equal(Dr, Drtrue, decimal=4)
+
+
+@pytest.mark.parametrize("kwargs, nstates, noutputs, ninputs", [
+    ({'elim_states': [1, 3]}, 3, 3, 3),
+    ({'elim_inputs': [1, 2], 'keep_states': [1, 3]}, 2, 3, 1),
+    ({'elim_outputs': [1, 2], 'keep_inputs': [0, 1],}, 5, 1, 2),
+    ({'keep_states': [2, 0], 'keep_outputs': [0, 1]}, 2, 2, 3),
+    ({'keep_states': slice(0, 4, 2), 'keep_outputs': slice(None, 2)}, 2, 2, 3),
+    ({'keep_states': ['x[0]', 'x[3]'], 'keep_inputs': 'u[0]'}, 2, 3, 1),
+    ({'elim_inputs': [0, 1, 2]}, 5, 3, 0),              # no inputs
+    ({'elim_outputs': [0, 1, 2]}, 5, 0, 3),             # no outputs
+    ({'elim_states': [0, 1, 2, 3, 4]}, 0, 3, 3),        # no states
+    ({'elim_states': [0, 1], 'keep_states': [1, 2]}, None, None, None),
+])
+@pytest.mark.parametrize("method", ['truncate', 'matchdc'])
+def test_model_reduction(method, kwargs, nstates, noutputs, ninputs):
+    sys = ct.rss(5, 3, 3)
+
+    if nstates is None:
+        # Arguments should generate an error
+        with pytest.raises(ValueError, match="can't provide both"):
+            red = ct.model_reduction(sys, **kwargs, method=method)
+        return
+    else:
+        red = ct.model_reduction(sys, **kwargs, method=method)
+
+    assert red.nstates == nstates
+    assert red.ninputs == ninputs
+    assert red.noutputs == noutputs
+
+    if method == 'matchdc':
+        # Define a new system with truncated inputs and outputs
+        # (assumes we always keep the initial inputs and outputs)
+        chk = ct.ss(
+            sys.A, sys.B[:, :ninputs], sys.C[:noutputs, :],
+            sys.D[:noutputs, :][:, :ninputs])
+        np.testing.assert_allclose(red(0), chk(0))