From 747935398e78fa83e17aa7d4085a449dcedea27c Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Wed, 6 Nov 2024 08:23:37 -0800
Subject: [PATCH 1/8] create find_operating_point with updated documentation

---
 control/nlsys.py | 105 ++++++++++++++++++++++++++++++++---------------
 1 file changed, 72 insertions(+), 33 deletions(-)

diff --git a/control/nlsys.py b/control/nlsys.py
index accd24c0f..768dfa00a 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -13,7 +13,7 @@
 :class:`~control.NonlinearIOSystem` class that represents (possibly nonlinear)
 input/output systems.  The :class:`~control.NonlinearIOSystem` class is a
 general class that defines any continuous or discrete time dynamical system.
-Input/output systems can be simulated and also used to compute equilibrium
+Input/output systems can be simulated and also used to compute operating
 points and linearizations.
 
 """
@@ -1663,69 +1663,98 @@ def ivp_rhs(t, x):
         success=soln.success, message=message)
 
 
-def find_eqpt(sys, x0, u0=None, y0=None, t=0, params=None,
-              iu=None, iy=None, ix=None, idx=None, dx0=None,
-              return_y=False, return_result=False):
-    """Find the equilibrium point for an input/output system.
+def find_operating_point(
+        sys, x0, u0=None, y0=None, t=0, params=None,
+        iu=None, iy=None, ix=None, idx=None, dx0=None, root_method=None,
+        root_kwargs=None, return_y=False, return_result=False):
+    """Find an operating point for an input/output system.
 
-    Returns the value of an equilibrium point given the initial state and
-    either input value or desired output value for the equilibrium point.
+    An operating point for a nonlinear system is a state `xop` and input
+    `uop` around which a nonlinear system operates.  This point is most
+    commonly an equilibrium point for the system, but in some cases a
+    non-equilibrium operating point can be used.
+
+    This function attempts to find an operating point given a specification
+    for the desired inputs, outputs, states, or state updates of the system.
+
+    In its simplest form, `find_operating_point` finds an equilibrium point
+    given either the desired input or desired output:
+
+        xeq, ueq = find_operating_point(sys, x0, u0)
+        xeq, ueq = find_operating_point(sys, x0, u0, y0)
+
+    The first form finds an equilibrium point for a given input u0 based on
+    an initial guess x0.  The second form fixes the desired output values
+    and uses x0 and u0 as an initial guess to find the equilibrium point.
+    If no equilibrium point can be found, the function returns the
+    operating point that minimizes the state update (state derivative for
+    continuous time systems, state difference for discrete time systems).
+
+    More complex operating points can be found by specifying which states,
+    inputs, or outputs should be used in computing the operating point, as
+    well as desired values of the states, inputs, outputs, or state
+    updates.
 
     Parameters
     ----------
     sys : NonlinearIOSystem
-        I/O system for which the equilibrium point is sought.
+        I/O system for which the operating point is sought.
     x0 : list of initial state values
-        Initial guess for the value of the state near the equilibrium point.
+        Initial guess for the value of the state near the operating point.
     u0 : list of input values, optional
-        If `y0` is not specified, sets the equilibrium value of the input.  If
-        `y0` is given, provides an initial guess for the value of the input.
-        Can be omitted if the system does not have any inputs.
+        If `y0` is not specified, sets the value of the input.  If `y0` is
+        given, provides an initial guess for the value of the input.  Can
+        be omitted if the system does not have any inputs.
     y0 : list of output values, optional
         If specified, sets the desired values of the outputs at the
-        equilibrium point.
+        operating point.
     t : float, optional
-        Evaluation time, for time-varying systems
+        Evaluation time, for time-varying systems.
     params : dict, optional
         Parameter values for the system.  Passed to the evaluation functions
         for the system as default values, overriding internal defaults.
     iu : list of input indices, optional
         If specified, only the inputs with the given indices will be fixed at
-        the specified values in solving for an equilibrium point.  All other
+        the specified values in solving for an operating point.  All other
         inputs will be varied.  Input indices can be listed in any order.
     iy : list of output indices, optional
         If specified, only the outputs with the given indices will be fixed at
-        the specified values in solving for an equilibrium point.  All other
+        the specified values in solving for an operating point.  All other
         outputs will be varied.  Output indices can be listed in any order.
     ix : list of state indices, optional
         If specified, states with the given indices will be fixed at the
-        specified values in solving for an equilibrium point.  All other
+        specified values in solving for an operating point.  All other
         states will be varied.  State indices can be listed in any order.
     dx0 : list of update values, optional
         If specified, the value of update map must match the listed value
-        instead of the default value of 0.
+        instead of the default value for an equilibrium point.
     idx : list of state indices, optional
         If specified, state updates with the given indices will have their
         update maps fixed at the values given in `dx0`.  All other update
-        values will be ignored in solving for an equilibrium point.  State
+        values will be ignored in solving for an operating point.  State
         indices can be listed in any order.  By default, all updates will be
-        fixed at `dx0` in searching for an equilibrium point.
+        fixed at `dx0` in searching for an operating point.
+    root_method : str, optonal
+        Method to find the operating point.  If specified, this parameter
+        is passed to the :func:`scipy.optimize.root` function.
+    root_kwargs : dict, optional
+        Additional keyword arguments to pass :func:`scipy.optimize.root`.
     return_y : bool, optional
-        If True, return the value of output at the equilibrium point.
+        If True, return the value of output at the operating point.
     return_result : bool, optional
         If True, return the `result` option from the
-        :func:`scipy.optimize.root` function used to compute the equilibrium
-        point.
+        :func:`scipy.optimize.root` function used to compute the
+        operating point.
 
     Returns
     -------
-    xeq : array of states
+    xop : array of states
         Value of the states at the equilibrium point, or `None` if no
         equilibrium point was found and `return_result` was False.
-    ueq : array of input values
+    uop : array of input values
         Value of the inputs at the equilibrium point, or `None` if no
         equilibrium point was found and `return_result` was False.
-    yeq : array of output values, optional
+    yop : array of output values, optional
         If `return_y` is True, returns the value of the outputs at the
         equilibrium point, or `None` if no equilibrium point was found and
         `return_result` was False.
@@ -1744,6 +1773,11 @@ def find_eqpt(sys, x0, u0=None, y0=None, t=0, params=None,
     """
     from scipy.optimize import root
 
+    # Process arguments for the root function
+    root_kwargs = dict() if root_kwargs is None else root_kwargs
+    if root_method:
+        root_kwargs['method'] = root_method
+
     # Figure out the number of states, inputs, and outputs
     x0, nstates = _process_vector_argument(x0, "x0", sys.nstates)
     u0, ninputs = _process_vector_argument(u0, "u0", sys.ninputs)
@@ -1769,7 +1803,7 @@ def state_rhs(z): return sys._rhs(t, z, u0) - z
             else:
                 def state_rhs(z): return sys._rhs(t, z, u0)
 
-            result = root(state_rhs, x0)
+            result = root(state_rhs, x0, **root_kwargs)
             z = (result.x, u0, sys._out(t, result.x, u0))
 
         else:
@@ -1786,9 +1820,10 @@ def rootfun(z):
                     return np.concatenate(
                         (sys._rhs(t, x, u), sys._out(t, x, u) - y0), axis=0)
 
-            z0 = np.concatenate((x0, u0), axis=0)   # Put variables together
-            result = root(rootfun, z0)              # Find the eq point
-            x, u = np.split(result.x, [nstates])    # Split result back in two
+            # Find roots with (x, u) as free variables
+            z0 = np.concatenate((x0, u0), axis=0)
+            result = root(rootfun, z0, **root_kwargs)
+            x, u = np.split(result.x, [nstates])
             z = (x, u, sys._out(t, x, u))
 
     else:
@@ -1903,7 +1938,7 @@ def rootfun(z):
         z0 = np.concatenate((x[state_vars], u[input_vars]), axis=0)
 
         # Finally, call the root finding function
-        result = root(rootfun, z0)
+        result = root(rootfun, z0, **root_kwargs)
 
         # Extract out the results and insert into x and u
         x[state_vars] = result.x[:nstate_vars]
@@ -2509,7 +2544,7 @@ def _find_output_or_input_signal(spec):
                             (syslist[isys].name,
                              syslist[isys].input_labels[isig], gain))
                 return signal_list
-                
+
             if isinstance(connection, list):
                 # Passed a list => create input map
                 dprint(f"  detected output list")
@@ -2519,7 +2554,7 @@ def _find_output_or_input_signal(spec):
                 new_outlist.append(signal_list)
             else:
                 new_outlist += _find_output_or_input_signal(connection)
-                
+
     outlist, outputs = new_outlist, new_outputs
     dprint(f"  {outlist=}\n  {outputs=}")
 
@@ -2679,3 +2714,7 @@ def connection_table(sys, show_names=False, column_width=32):
         "an InterconnectedSystem."
 
     sys.connection_table(show_names=show_names, column_width=column_width)
+
+
+# Short versions of function call
+find_eqpt = find_operating_point

From fa162ece64fc79ac6ac8f79ba45e80ff5c8aa719 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Wed, 6 Nov 2024 14:08:37 -0800
Subject: [PATCH 2/8] add OperatingPoint object + unit test

---
 control/nlsys.py            | 65 ++++++++++++++++++++++++++++++++++++-
 control/tests/iosys_test.py | 36 ++++++++++++++++++++
 2 files changed, 100 insertions(+), 1 deletion(-)

diff --git a/control/nlsys.py b/control/nlsys.py
index 768dfa00a..33d9379c7 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -32,7 +32,7 @@
 
 __all__ = ['NonlinearIOSystem', 'InterconnectedSystem', 'nlsys',
            'input_output_response', 'find_eqpt', 'linearize',
-           'interconnect', 'connection_table']
+           'interconnect', 'connection_table', 'find_operating_point']
 
 
 class NonlinearIOSystem(InputOutputSystem):
@@ -1663,6 +1663,60 @@ def ivp_rhs(t, x):
         success=soln.success, message=message)
 
 
+class OperatingPoint(object):
+    """A class for representing the operating point of a nonlinear I/O system.
+
+    The ``OperatingPoint`` class stores the operating point of a nonlinear
+    system, which consists of the state and input for a nonlinear system.
+    The main use for this class is as the return object for the
+    :func:`find_operating_point` function.
+
+    Attributes
+    ----------
+    xop : array
+        State vector at the operating point.
+    uop : array
+        Input vector at the operating point.
+    result : :class:`scipy.optimize.OptimizeResult`, optional
+        Result from the :func:`scipy.optimize.root` function, if available.
+
+    """
+    def __init__(
+            self, xop, uop=None, yop=None, result=None,
+            return_y=False, return_result=False):
+        self.xop = xop
+        self.uop = uop
+
+        if yop is None and return_y and not return_result:
+            raise SystemError("return_y specified by no y0 value")
+        self.yop = yop
+        self.return_y = return_y
+
+        if result is None and return_result:
+            raise SystemError("return_result specified by no result value")
+        self.result = result
+        self.return_result = return_result
+
+    # Implement iter to allow assigning to a tuple
+    def __iter__(self):
+        if self.return_y and self.return_result:
+            return iter((self.xop, self.uop, self.yop, self.result))
+        elif self.return_y:
+            return iter((self.xop, self.uop, self.yop))
+        elif self.return_result:
+            return iter((self.xop, self.uop, self.result))
+        else:
+            return iter((self.xop, self.uop))
+
+    # Implement (thin) getitem to allow access via legacy indexing
+    def __getitem__(self, index):
+        return list(self.__iter__())[index]
+
+    # Implement (thin) len to emulate legacy return value
+    def __len__(self):
+        return len(list(self.__iter__()))
+
+
 def find_operating_point(
         sys, x0, u0=None, y0=None, t=0, params=None,
         iu=None, iy=None, ix=None, idx=None, dx0=None, root_method=None,
@@ -1946,6 +2000,15 @@ def rootfun(z):
         z = (x, u, sys._out(t, x, u))
 
     # Return the result based on what the user wants and what we found
+    if return_result or result.success:
+        return OperatingPoint(
+            z[0], z[1], z[2], result, return_y, return_result)
+    else:
+        # Something went wrong, don't return anything
+        return OperatingPoint(
+            None, None, None, result, return_y, return_result)
+
+    # TODO: remove code when ready
     if not return_y:
         z = z[0:2]              # Strip y from result if not desired
     if return_result:
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index dd30ea71e..beca8d445 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -2087,6 +2087,42 @@ def test_find_eqpt(x0, ix, u0, iu, y0, iy, dx0, idx, dt, x_expect, u_expect):
     np.testing.assert_allclose(np.array(ueq), u_expect, atol=1e-6)
 
 
+# Test out new operating point version of find_eqpt
+def test_find_operating_point():
+    dt = 1
+    sys = ct.NonlinearIOSystem(
+        eqpt_rhs, eqpt_out, dt=dt, states=3, inputs=2, outputs=2)
+
+    # Conditions that lead to no exact solution (from previous unit test)
+    x0 = 0;       ix = None
+    u0 = [-1, 0]; iu = None
+    y0 = None;    iy = None
+    dx0 = None;  idx = None
+
+    # Default version: no equilibrium solution => returns None
+    op_point = ct.find_operating_point(
+        sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx)
+    assert op_point.xop is None
+    assert op_point.uop is None
+    assert op_point.result.success is False
+
+    # Change the method to Levenberg-Marquardt (gives nearest point)
+    op_point = ct.find_operating_point(
+        sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx,
+        root_method='lm')
+    assert op_point.xop is not None
+    assert op_point.uop is not None
+    assert op_point.result.success is True
+
+    # Make sure we get a solution if we ask for the result explicitly
+    op_point = ct.find_operating_point(
+        sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx,
+        return_result=True)
+    assert op_point.xop is not None
+    assert op_point.uop is not None
+    assert op_point.result.success is False
+
+
 def test_iosys_sample():
     csys = ct.rss(2, 1, 1)
     dsys = csys.sample(0.1)

From 223acc64f92e7119a587c93f7e62ebc3a949973d Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Wed, 6 Nov 2024 21:31:46 -0800
Subject: [PATCH 3/8] allow OperatingPoint for linearize + code and
 documentation tweaks

---
 control/matlab/__init__.py  |  2 ++
 control/nlsys.py            | 52 ++++++++++++++++++++++---------------
 control/tests/iosys_test.py | 32 ++++++++++++++++++-----
 doc/control.rst             |  2 +-
 doc/iosys.rst               | 18 ++++++-------
 5 files changed, 69 insertions(+), 37 deletions(-)

diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py
index b02d16d53..98e6babc7 100644
--- a/control/matlab/__init__.py
+++ b/control/matlab/__init__.py
@@ -84,10 +84,12 @@
 from ..dtime import c2d
 from ..sisotool import sisotool
 from ..stochsys import lqe, dlqe
+from ..nlsys import find_operating_point
 
 # Functions that are renamed in MATLAB
 pole, zero = poles, zeros
 freqresp = frequency_response
+trim = find_operating_point
 
 # Import functions specific to Matlab compatibility package
 from .timeresp import *
diff --git a/control/nlsys.py b/control/nlsys.py
index 33d9379c7..761525d4d 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -515,7 +515,7 @@ def feedback(self, other=1, sign=-1, params=None):
         # Return the newly created system
         return newsys
 
-    def linearize(self, x0, u0, t=0, params=None, eps=1e-6,
+    def linearize(self, x0, u0=None, t=0, params=None, eps=1e-6,
                   copy_names=False, **kwargs):
         """Linearize an input/output system at a given state and input.
 
@@ -526,6 +526,14 @@ def linearize(self, x0, u0, t=0, params=None, eps=1e-6,
         """
         from .statesp import StateSpace
 
+        # Allow first argument to be an operating point
+        if isinstance(x0, OperatingPoint):
+            if u0 is None:
+                u0 = x0.inputs
+            x0 = x0.states
+        elif u0 is None:
+            u0 = 0
+
         #
         # If the linearization is not defined by the subclass, perform a
         # numerical linearization use the `_rhs()` and `_out()` member
@@ -1673,19 +1681,19 @@ class OperatingPoint(object):
 
     Attributes
     ----------
-    xop : array
+    states : array
         State vector at the operating point.
-    uop : array
+    inputs : array
         Input vector at the operating point.
     result : :class:`scipy.optimize.OptimizeResult`, optional
         Result from the :func:`scipy.optimize.root` function, if available.
 
     """
     def __init__(
-            self, xop, uop=None, yop=None, result=None,
+            self, states, inputs=None, yop=None, result=None,
             return_y=False, return_result=False):
-        self.xop = xop
-        self.uop = uop
+        self.states = states
+        self.inputs = inputs
 
         if yop is None and return_y and not return_result:
             raise SystemError("return_y specified by no y0 value")
@@ -1700,13 +1708,13 @@ def __init__(
     # Implement iter to allow assigning to a tuple
     def __iter__(self):
         if self.return_y and self.return_result:
-            return iter((self.xop, self.uop, self.yop, self.result))
+            return iter((self.states, self.inputs, self.yop, self.result))
         elif self.return_y:
-            return iter((self.xop, self.uop, self.yop))
+            return iter((self.states, self.inputs, self.yop))
         elif self.return_result:
-            return iter((self.xop, self.uop, self.result))
+            return iter((self.states, self.inputs, self.result))
         else:
-            return iter((self.xop, self.uop))
+            return iter((self.states, self.inputs))
 
     # Implement (thin) getitem to allow access via legacy indexing
     def __getitem__(self, index):
@@ -1723,10 +1731,10 @@ def find_operating_point(
         root_kwargs=None, return_y=False, return_result=False):
     """Find an operating point for an input/output system.
 
-    An operating point for a nonlinear system is a state `xop` and input
-    `uop` around which a nonlinear system operates.  This point is most
-    commonly an equilibrium point for the system, but in some cases a
-    non-equilibrium operating point can be used.
+    An operating point for a nonlinear system is a state and input around
+    which a nonlinear system operates.  This point is most commonly an
+    equilibrium point for the system, but in some cases a non-equilibrium
+    operating point can be used.
 
     This function attempts to find an operating point given a specification
     for the desired inputs, outputs, states, or state updates of the system.
@@ -1802,10 +1810,10 @@ def find_operating_point(
 
     Returns
     -------
-    xop : array of states
+    states : array of states
         Value of the states at the equilibrium point, or `None` if no
         equilibrium point was found and `return_result` was False.
-    uop : array of input values
+    inputs : array of input values
         Value of the inputs at the equilibrium point, or `None` if no
         equilibrium point was found and `return_result` was False.
     yop : array of output values, optional
@@ -2034,12 +2042,13 @@ def linearize(sys, xeq, ueq=None, t=0, params=None, **kw):
     ----------
     sys : InputOutputSystem
         The system to be linearized.
-    xeq : array
-        The state at which the linearization will be evaluated (does not need
-        to be an equilibrium state).
-    ueq : array
+    xeq : array or :class:`~control.OperatingPoint`
+        The state or operating point at which the linearization will be
+        evaluated (does not need to be an equilibrium state).
+    ueq : array, optional
         The input at which the linearization will be evaluated (does not need
-        to correspond to an equlibrium state).
+        to correspond to an equlibrium state).  Can be omitted if `xeq` is
+        an :class:`~control.OperatingPoint`.  Defaults to 0.
     t : float, optional
         The time at which the linearization will be computed (for time-varying
         systems).
@@ -2074,6 +2083,7 @@ def linearize(sys, xeq, ueq=None, t=0, params=None, **kw):
         Description of the system outputs.  Same format as `inputs`.
     states : int, list of str, or None, optional
         Description of the system states.  Same format as `inputs`.
+
     """
     if not isinstance(sys, InputOutputSystem):
         raise TypeError("Can only linearize InputOutputSystem types")
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index beca8d445..608fdc58c 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -2102,27 +2102,47 @@ def test_find_operating_point():
     # Default version: no equilibrium solution => returns None
     op_point = ct.find_operating_point(
         sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx)
-    assert op_point.xop is None
-    assert op_point.uop is None
+    assert op_point.states is None
+    assert op_point.inputs is None
     assert op_point.result.success is False
 
     # Change the method to Levenberg-Marquardt (gives nearest point)
     op_point = ct.find_operating_point(
         sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx,
         root_method='lm')
-    assert op_point.xop is not None
-    assert op_point.uop is not None
+    assert op_point.states is not None
+    assert op_point.inputs is not None
     assert op_point.result.success is True
 
     # Make sure we get a solution if we ask for the result explicitly
     op_point = ct.find_operating_point(
         sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx,
         return_result=True)
-    assert op_point.xop is not None
-    assert op_point.uop is not None
+    assert op_point.states is not None
+    assert op_point.inputs is not None
     assert op_point.result.success is False
 
 
+def test_operating_point():
+    dt = 1
+    sys = ct.NonlinearIOSystem(
+        eqpt_rhs, eqpt_out, dt=dt, states=3, inputs=2, outputs=2)
+
+    # Find the operating point near the origin
+    op_point = ct.find_operating_point(sys, 0, 0)
+
+    # Linearize the old fashioned way
+    linsys_orig = ct.linearize(sys, op_point.states, op_point.inputs)
+
+    # Linearize around the operating point
+    linsys_oppt = ct.linearize(sys, op_point)
+
+    np.testing.assert_allclose(linsys_orig.A, linsys_oppt.A)
+    np.testing.assert_allclose(linsys_orig.B, linsys_oppt.B)
+    np.testing.assert_allclose(linsys_orig.C, linsys_oppt.C)
+    np.testing.assert_allclose(linsys_orig.D, linsys_oppt.D)
+
+
 def test_iosys_sample():
     csys = ct.rss(2, 1, 1)
     dsys = csys.sample(0.1)
diff --git a/doc/control.rst b/doc/control.rst
index 1544f93d0..dd418f2af 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -144,7 +144,7 @@ Nonlinear system support
    :toctree: generated/
 
     describing_function
-    find_eqpt
+    find_operating_point
     linearize
     input_output_response
     summing_junction
diff --git a/doc/iosys.rst b/doc/iosys.rst
index 9008a0e56..db504c0ab 100644
--- a/doc/iosys.rst
+++ b/doc/iosys.rst
@@ -16,12 +16,13 @@ function::
   resp = ct.input_output_response(io_sys, T, U, X0, params)
   t, y, x = resp.time, resp.outputs, resp.states
 
-An input/output system can be linearized around an equilibrium point to obtain
-a :class:`~control.StateSpace` linear system.  Use the
-:func:`~control.find_eqpt` function to obtain an equilibrium point and the
-:func:`~control.linearize` function to linearize about that equilibrium point::
+An input/output system can be linearized around an equilibrium point
+to obtain a :class:`~control.StateSpace` linear system.  Use the
+:func:`~control.find_operating_point` function to obtain an
+equilibrium point and the :func:`~control.linearize` function to
+linearize about that equilibrium point::
 
-  xeq, ueq = ct.find_eqpt(io_sys, X0, U0)
+  xeq, ueq = ct.find_operating_point(io_sys, X0, U0)
   ss_sys = ct.linearize(io_sys, xeq, ueq)
 
 Input/output systems are automatically created for state space LTI systems
@@ -123,9 +124,8 @@ system and computing the linearization about that point.
 
 .. code-block:: python
 
-  eqpt = ct.find_eqpt(io_predprey, X0, 0)
-  xeq = eqpt[0]                         # choose the nonzero equilibrium point
-  lin_predprey = ct.linearize(io_predprey, xeq, 0)
+  eqpt = ct.find_operating_point(io_predprey, X0, 0)
+  lin_predprey = ct.linearize(io_predprey, eqpt)
 
 We next compute a controller that stabilizes the equilibrium point using
 eigenvalue placement and computing the feedforward gain using the number of
@@ -548,7 +548,7 @@ Module classes and functions
 .. autosummary::
    :toctree: generated/
 
-   ~control.find_eqpt
+   ~control.find_operating_point
    ~control.interconnect
    ~control.input_output_response
    ~control.linearize

From 5eaee668fe3166318a1907481f95c02afe444e2c Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Thu, 7 Nov 2024 19:26:43 -0800
Subject: [PATCH 4/8] updated OperatingPoint attribute names + docstring edits

---
 control/nlsys.py | 64 ++++++++++++++++++++++++++++--------------------
 1 file changed, 37 insertions(+), 27 deletions(-)

diff --git a/control/nlsys.py b/control/nlsys.py
index 761525d4d..0af79c8d8 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -519,26 +519,25 @@ def linearize(self, x0, u0=None, t=0, params=None, eps=1e-6,
                   copy_names=False, **kwargs):
         """Linearize an input/output system at a given state and input.
 
-        Return the linearization of an input/output system at a given state
-        and input value as a StateSpace system.  See
-        :func:`~control.linearize` for complete documentation.
+        Return the linearization of an input/output system at a given
+        operating point (or state and input value) as a StateSpace system.
+        See :func:`~control.linearize` for complete documentation.
 
         """
+        #
+        # Default method: if the linearization is not defined by the
+        # subclass, perform a numerical linearization use the `_rhs()` and
+        # `_out()` member functions.
+        #
         from .statesp import StateSpace
 
         # Allow first argument to be an operating point
         if isinstance(x0, OperatingPoint):
-            if u0 is None:
-                u0 = x0.inputs
+            u0 = x0.inputs if u0 is None else u0
             x0 = x0.states
         elif u0 is None:
             u0 = 0
 
-        #
-        # If the linearization is not defined by the subclass, perform a
-        # numerical linearization use the `_rhs()` and `_out()` member
-        # functions.
-        #
         # Process nominal states and inputs
         x0, nstates = _process_vector_argument(x0, "x0", self.nstates)
         u0, ninputs = _process_vector_argument(u0, "u0", self.ninputs)
@@ -1675,9 +1674,10 @@ class OperatingPoint(object):
     """A class for representing the operating point of a nonlinear I/O system.
 
     The ``OperatingPoint`` class stores the operating point of a nonlinear
-    system, which consists of the state and input for a nonlinear system.
-    The main use for this class is as the return object for the
-    :func:`find_operating_point` function.
+    system, consisting of the state and input vectors for the system.  The
+    main use for this class is as the return object for the
+    :func:`find_operating_point` function and as an input to the
+    :func:`linearize` function.
 
     Attributes
     ----------
@@ -1690,14 +1690,14 @@ class OperatingPoint(object):
 
     """
     def __init__(
-            self, states, inputs=None, yop=None, result=None,
+            self, states, inputs=None, outputs=None, result=None,
             return_y=False, return_result=False):
         self.states = states
         self.inputs = inputs
 
-        if yop is None and return_y and not return_result:
+        if outputs is None and return_y and not return_result:
             raise SystemError("return_y specified by no y0 value")
-        self.yop = yop
+        self.outputs = outputs
         self.return_y = return_y
 
         if result is None and return_result:
@@ -1708,9 +1708,9 @@ def __init__(
     # Implement iter to allow assigning to a tuple
     def __iter__(self):
         if self.return_y and self.return_result:
-            return iter((self.states, self.inputs, self.yop, self.result))
+            return iter((self.states, self.inputs, self.outputs, self.result))
         elif self.return_y:
-            return iter((self.states, self.inputs, self.yop))
+            return iter((self.states, self.inputs, self.outputs))
         elif self.return_result:
             return iter((self.states, self.inputs, self.result))
         else:
@@ -1816,7 +1816,7 @@ def find_operating_point(
     inputs : array of input values
         Value of the inputs at the equilibrium point, or `None` if no
         equilibrium point was found and `return_result` was False.
-    yop : array of output values, optional
+    outputs : array of output values, optional
         If `return_y` is True, returns the value of the outputs at the
         equilibrium point, or `None` if no equilibrium point was found and
         `return_result` was False.
@@ -1826,11 +1826,21 @@ def find_operating_point(
 
     Notes
     -----
-    For continuous time systems, equilibrium points are defined as points for
-    which the right hand side of the differential equation is zero:
-    :math:`f(t, x_e, u_e) = 0`. For discrete time systems, equilibrium points
-    are defined as points for which the right hand side of the difference
-    equation returns the current state: :math:`f(t, x_e, u_e) = x_e`.
+    For continuous time systems, equilibrium points are defined as points
+    for which the right hand side of the differential equation is zero:
+    :math:`f(t, x_e, u_e) = 0`. For discrete time systems, equilibrium
+    points are defined as points for which the right hand side of the
+    difference equation returns the current state: :math:`f(t, x_e, u_e) =
+    x_e`.
+
+    Operating points are found using the :func:`scipy.optimize.root`
+    function, which will attempt to find states and inputs that satisfy the
+    specified constraints.  If no solution is found and `return_result` is
+    `False`, the returned state and input for the operating point will be
+    `None`.  If `return_result` is `True`, then the return values from
+    :func:`scipy.optimize.root` will be returned (but may not be valid).
+    If `root_method` is set to `lm`, then the least squares solution (in
+    the free variables) will be returned.
 
     """
     from scipy.optimize import root
@@ -2034,9 +2044,9 @@ def rootfun(z):
 def linearize(sys, xeq, ueq=None, t=0, params=None, **kw):
     """Linearize an input/output system at a given state and input.
 
-    This function computes the linearization of an input/output system at a
-    given state and input value and returns a :class:`~control.StateSpace`
-    object.  The evaluation point need not be an equilibrium point.
+    Compute the linearization of an I/O system at an operating point (state
+    and input) and returns a :class:`~control.StateSpace` object.  The
+    operating point need not be an equilibrium point.
 
     Parameters
     ----------

From aed5706b5bf9ae47a2743a1c6ba476c05eb49473 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Thu, 7 Nov 2024 22:06:08 -0800
Subject: [PATCH 5/8] clean up docstrings, examples

---
 control/nlsys.py            | 37 ++++++++++++++++++-------------------
 control/phaseplot.py        |  5 +++--
 control/tests/iosys_test.py |  4 ++++
 doc/iosys.rst               |  1 +
 examples/cruise-control.py  |  4 ++--
 5 files changed, 28 insertions(+), 23 deletions(-)

diff --git a/control/nlsys.py b/control/nlsys.py
index 0af79c8d8..cd5acbf9a 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -32,7 +32,8 @@
 
 __all__ = ['NonlinearIOSystem', 'InterconnectedSystem', 'nlsys',
            'input_output_response', 'find_eqpt', 'linearize',
-           'interconnect', 'connection_table', 'find_operating_point']
+           'interconnect', 'connection_table', 'OperatingPoint',
+           'find_operating_point']
 
 
 class NonlinearIOSystem(InputOutputSystem):
@@ -92,7 +93,7 @@ class NonlinearIOSystem(InputOutputSystem):
         generic name <sys[id]> is generated with a unique integer id.
 
     params : dict, optional
-        Parameter values for the system.  Passed to the evaluation functions
+       Parameter values for the system.  Passed to the evaluation functions
         for the system as default values, overriding internal defaults.
 
     See Also
@@ -1685,6 +1686,8 @@ class OperatingPoint(object):
         State vector at the operating point.
     inputs : array
         Input vector at the operating point.
+    outputs : array, optional
+        Output vector at the operating point.
     result : :class:`scipy.optimize.OptimizeResult`, optional
         Result from the :func:`scipy.optimize.root` function, if available.
 
@@ -1740,7 +1743,7 @@ def find_operating_point(
     for the desired inputs, outputs, states, or state updates of the system.
 
     In its simplest form, `find_operating_point` finds an equilibrium point
-    given either the desired input or desired output:
+    given either the desired input or desired output::
 
         xeq, ueq = find_operating_point(sys, x0, u0)
         xeq, ueq = find_operating_point(sys, x0, u0, y0)
@@ -1810,19 +1813,15 @@ def find_operating_point(
 
     Returns
     -------
-    states : array of states
-        Value of the states at the equilibrium point, or `None` if no
-        equilibrium point was found and `return_result` was False.
-    inputs : array of input values
-        Value of the inputs at the equilibrium point, or `None` if no
-        equilibrium point was found and `return_result` was False.
-    outputs : array of output values, optional
-        If `return_y` is True, returns the value of the outputs at the
-        equilibrium point, or `None` if no equilibrium point was found and
-        `return_result` was False.
-    result : :class:`scipy.optimize.OptimizeResult`, optional
-        If `return_result` is True, returns the `result` from the
-        :func:`scipy.optimize.root` function.
+    op_point : OperatingPoint
+        The solution represented as an `OperatingPoint` object.  The main
+        attributes are `states` and `inputs`, which represent the state
+        and input arrays at the operating point.  See
+        :class:`OperatingPoint` for a description of other attributes.
+
+        If accessed as a tuple, returns `states`, `inputs`, and optionally
+        `outputs` and `result` based on the `return_y` and `return_result`
+        parameters.
 
     Notes
     -----
@@ -1839,7 +1838,7 @@ def find_operating_point(
     `False`, the returned state and input for the operating point will be
     `None`.  If `return_result` is `True`, then the return values from
     :func:`scipy.optimize.root` will be returned (but may not be valid).
-    If `root_method` is set to `lm`, then the least squares solution (in
+    If `root_method` is set to 'lm', then the least squares solution (in
     the free variables) will be returned.
 
     """
@@ -1956,8 +1955,8 @@ def rootfun(z):
         # * output_vars: indices of outputs that must be constrained
         #
         # This index lists can all be precomputed based on the `iu`, `iy`,
-        # `ix`, and `idx` lists that were passed as arguments to `find_eqpt`
-        # and were processed above.
+        # `ix`, and `idx` lists that were passed as arguments to
+        # `find_operating_point` and were processed above.
 
         # Get the states and inputs that were not listed as fixed
         state_vars = (range(nstates) if not len(ix)
diff --git a/control/phaseplot.py b/control/phaseplot.py
index abc050ffe..6f7ed355f 100644
--- a/control/phaseplot.py
+++ b/control/phaseplot.py
@@ -39,7 +39,8 @@
 from .ctrlplot import ControlPlot, _add_arrows_to_line2D, _get_color, \
     _process_ax_keyword, _update_plot_title
 from .exception import ControlNotImplemented
-from .nlsys import NonlinearIOSystem, find_eqpt, input_output_response
+from .nlsys import NonlinearIOSystem, find_operating_point, \
+    input_output_response
 
 __all__ = ['phase_plane_plot', 'phase_plot', 'box_grid']
 
@@ -853,7 +854,7 @@ def _find_equilpts(sys, points, params=None):
     equilpts = []
     for i, x0 in enumerate(points):
         # Look for an equilibrium point near this point
-        xeq, ueq = find_eqpt(sys, x0, 0, params=params)
+        xeq, ueq = find_operating_point(sys, x0, 0, params=params)
 
         if xeq is None:
             continue            # didn't find anything
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 608fdc58c..09f5596ec 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -2142,6 +2142,10 @@ def test_operating_point():
     np.testing.assert_allclose(linsys_orig.C, linsys_oppt.C)
     np.testing.assert_allclose(linsys_orig.D, linsys_oppt.D)
 
+    # Call find_operating point with method and keyword arguments
+    op_point = ct.find_operating_point(
+        sys, 0, 0, root_method='lm', root_kwargs={'tol': 1e-6})
+
 
 def test_iosys_sample():
     csys = ct.rss(2, 1, 1)
diff --git a/doc/iosys.rst b/doc/iosys.rst
index db504c0ab..f2dfbff4d 100644
--- a/doc/iosys.rst
+++ b/doc/iosys.rst
@@ -544,6 +544,7 @@ Module classes and functions
    ~control.InterconnectedSystem
    ~control.LinearICSystem
    ~control.NonlinearIOSystem
+   ~control.OperatingPoint
 
 .. autosummary::
    :toctree: generated/
diff --git a/examples/cruise-control.py b/examples/cruise-control.py
index 7c2e562a1..b8a91eeab 100644
--- a/examples/cruise-control.py
+++ b/examples/cruise-control.py
@@ -165,7 +165,7 @@ def motor_torque(omega, params={}):
 
 for m in (1200, 1600, 2000):
     # Compute the equilibrium state for the system
-    X0, U0 = ct.find_eqpt(
+    X0, U0 = ct.find_operating_point(
         cruise_tf, [0, vref[0]], [vref[0], gear[0], theta0[0]], 
         iu=[1, 2], y0=[vref[0], 0], iy=[0], params={'m': m})
 
@@ -347,7 +347,7 @@ def cruise_plot(sys, t, y, label=None, t_hill=None, vref=20, antiwindup=False,
 
 # Compute the equilibrium throttle setting for the desired speed (solve for x
 # and u given the gear, slope, and desired output velocity)
-X0, U0, Y0 = ct.find_eqpt(
+X0, U0, Y0 = ct.find_operating_point(
     cruise_pi, [vref[0], 0], [vref[0], gear[0], theta0[0]],
     y0=[0, vref[0]], iu=[1, 2], iy=[1], return_y=True)
 

From 1a2770c9c34295c5928d197f80e0fed85495aae5 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Fri, 8 Nov 2024 20:55:17 -0800
Subject: [PATCH 6/8] return_y -> return_outputs + related updates

---
 control/nlsys.py                 | 40 ++++++++++++++++++--------------
 control/tests/docstrings_test.py |  1 +
 control/tests/iosys_test.py      | 40 ++++++++++++++++++++++++++++++--
 control/tests/kwargs_test.py     |  3 +++
 4 files changed, 65 insertions(+), 19 deletions(-)

diff --git a/control/nlsys.py b/control/nlsys.py
index cd5acbf9a..d463ea7fc 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -1694,14 +1694,14 @@ class OperatingPoint(object):
     """
     def __init__(
             self, states, inputs=None, outputs=None, result=None,
-            return_y=False, return_result=False):
+            return_outputs=False, return_result=False):
         self.states = states
         self.inputs = inputs
 
-        if outputs is None and return_y and not return_result:
-            raise SystemError("return_y specified by no y0 value")
+        if outputs is None and return_outputs and not return_result:
+            raise SystemError("return_outputs specified by no y0 value")
         self.outputs = outputs
-        self.return_y = return_y
+        self.return_outputs = return_outputs
 
         if result is None and return_result:
             raise SystemError("return_result specified by no result value")
@@ -1710,9 +1710,9 @@ def __init__(
 
     # Implement iter to allow assigning to a tuple
     def __iter__(self):
-        if self.return_y and self.return_result:
+        if self.return_outputs and self.return_result:
             return iter((self.states, self.inputs, self.outputs, self.result))
-        elif self.return_y:
+        elif self.return_outputs:
             return iter((self.states, self.inputs, self.outputs))
         elif self.return_result:
             return iter((self.states, self.inputs, self.result))
@@ -1729,9 +1729,9 @@ def __len__(self):
 
 
 def find_operating_point(
-        sys, x0, u0=None, y0=None, t=0, params=None,
-        iu=None, iy=None, ix=None, idx=None, dx0=None, root_method=None,
-        root_kwargs=None, return_y=False, return_result=False):
+        sys, x0, u0=None, y0=None, t=0, params=None, iu=None, iy=None,
+        ix=None, idx=None, dx0=None, root_method=None, root_kwargs=None,
+        return_outputs=None, return_result=None, **kwargs):
     """Find an operating point for an input/output system.
 
     An operating point for a nonlinear system is a state and input around
@@ -1804,8 +1804,8 @@ def find_operating_point(
         is passed to the :func:`scipy.optimize.root` function.
     root_kwargs : dict, optional
         Additional keyword arguments to pass :func:`scipy.optimize.root`.
-    return_y : bool, optional
-        If True, return the value of output at the operating point.
+    return_outputs : bool, optional
+        If True, return the value of outputs at the operating point.
     return_result : bool, optional
         If True, return the `result` option from the
         :func:`scipy.optimize.root` function used to compute the
@@ -1820,8 +1820,8 @@ def find_operating_point(
         :class:`OperatingPoint` for a description of other attributes.
 
         If accessed as a tuple, returns `states`, `inputs`, and optionally
-        `outputs` and `result` based on the `return_y` and `return_result`
-        parameters.
+        `outputs` and `result` based on the `return_outputs` and
+        `return_result` parameters.
 
     Notes
     -----
@@ -1844,6 +1844,12 @@ def find_operating_point(
     """
     from scipy.optimize import root
 
+    # Process keyword arguments
+    return_outputs = config._process_legacy_keyword(
+        kwargs, 'return_y', 'return_outputs', return_outputs)
+    if kwargs:
+        raise TypeError("unrecognized keyword(s): " + str(kwargs))
+
     # Process arguments for the root function
     root_kwargs = dict() if root_kwargs is None else root_kwargs
     if root_method:
@@ -2019,14 +2025,14 @@ def rootfun(z):
     # Return the result based on what the user wants and what we found
     if return_result or result.success:
         return OperatingPoint(
-            z[0], z[1], z[2], result, return_y, return_result)
+            z[0], z[1], z[2], result, return_outputs, return_result)
     else:
         # Something went wrong, don't return anything
         return OperatingPoint(
-            None, None, None, result, return_y, return_result)
+            None, None, None, result, return_outputs, return_result)
 
     # TODO: remove code when ready
-    if not return_y:
+    if not return_outputs:
         z = z[0:2]              # Strip y from result if not desired
     if return_result:
         # Return whatever we got, along with the result dictionary
@@ -2036,7 +2042,7 @@ def rootfun(z):
         return z
     else:
         # Something went wrong, don't return anything
-        return (None, None, None) if return_y else (None, None)
+        return (None, None, None) if return_outputs else (None, None)
 
 
 # Linearize an input/output system
diff --git a/control/tests/docstrings_test.py b/control/tests/docstrings_test.py
index 70a0cc3c1..2b1368aeb 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -53,6 +53,7 @@
     control.create_estimator_iosystem: ['state_labels'],    # deprecated
     control.bode_plot: ['sharex', 'sharey', 'margin_info'], # deprecated
     control.eigensys_realization: ['arg'],                  # quasi-positional
+    control.find_operating_point: ['method'],               # internal use
 }
 
 # Decide on the level of verbosity (use -rP when running pytest)
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 09f5596ec..d2570be01 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -10,10 +10,11 @@
 
 import re
 import warnings
-import pytest
+from math import sqrt
 
 import numpy as np
-from math import sqrt
+import pytest
+import scipy
 
 import control as ct
 
@@ -2088,6 +2089,7 @@ def test_find_eqpt(x0, ix, u0, iu, y0, iy, dx0, idx, dt, x_expect, u_expect):
 
 
 # Test out new operating point version of find_eqpt
+# TODO: add return_)y tests
 def test_find_operating_point():
     dt = 1
     sys = ct.NonlinearIOSystem(
@@ -2122,6 +2124,10 @@ def test_find_operating_point():
     assert op_point.inputs is not None
     assert op_point.result.success is False
 
+    # Check to make sure unknown keywords are caught
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        ct.find_operating_point(sys, x0, u0, unknown=None)
+
 
 def test_operating_point():
     dt = 1
@@ -2146,6 +2152,36 @@ def test_operating_point():
     op_point = ct.find_operating_point(
         sys, 0, 0, root_method='lm', root_kwargs={'tol': 1e-6})
 
+    # Make sure we can get back the right arguments in a tuple
+    op_point = ct.find_operating_point(sys, 0, 0, return_outputs=True)
+    assert len(op_point) == 3
+    assert isinstance(op_point[0], np.ndarray)
+    assert isinstance(op_point[1], np.ndarray)
+    assert isinstance(op_point[2], np.ndarray)
+
+    with pytest.warns(FutureWarning, match="return_outputs"):
+        op_point = ct.find_operating_point(sys, 0, 0, return_y=True)
+        assert len(op_point) == 3
+        assert isinstance(op_point[0], np.ndarray)
+        assert isinstance(op_point[1], np.ndarray)
+        assert isinstance(op_point[2], np.ndarray)
+
+    # Make sure we can get back the right arguments in a tuple
+    op_point = ct.find_operating_point(sys, 0, 0, return_result=True)
+    assert len(op_point) == 3
+    assert isinstance(op_point[0], np.ndarray)
+    assert isinstance(op_point[1], np.ndarray)
+    assert isinstance(op_point[2], scipy.optimize.OptimizeResult)
+
+    # Make sure we can get back the right arguments in a tuple
+    op_point = ct.find_operating_point(
+        sys, 0, 0, return_result=True, return_outputs=True)
+    assert len(op_point) == 4
+    assert isinstance(op_point[0], np.ndarray)
+    assert isinstance(op_point[1], np.ndarray)
+    assert isinstance(op_point[2], np.ndarray)
+    assert isinstance(op_point[3], scipy.optimize.OptimizeResult)
+
 
 def test_iosys_sample():
     csys = ct.rss(2, 1, 1)
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
index 020910e73..b98509d65 100644
--- a/control/tests/kwargs_test.py
+++ b/control/tests/kwargs_test.py
@@ -24,6 +24,7 @@
 import control.tests.frd_test as frd_test
 import control.tests.freqplot_test as freqplot_test
 import control.tests.interconnect_test as interconnect_test
+import control.tests.iosys_test as iosys_test
 import control.tests.optimal_test as optimal_test
 import control.tests.statefbk_test as statefbk_test
 import control.tests.stochsys_test as stochsys_test
@@ -252,6 +253,8 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     'dlqr': test_unrecognized_kwargs,
     'drss': test_unrecognized_kwargs,
     'feedback': test_unrecognized_kwargs,
+    'find_eqpt': iosys_test.test_find_operating_point,
+    'find_operating_point': iosys_test.test_find_operating_point,
     'flatsys.flatsys': test_unrecognized_kwargs,
     'frd': frd_test.TestFRD.test_unrecognized_keyword,
     'gangof4': test_matplotlib_kwargs,

From 47e1c0955ca5b400571ac408ff0c5eeb27371b49 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Sun, 10 Nov 2024 19:00:00 -0800
Subject: [PATCH 7/8] fix small typos in documentation (unrelated to main PR)

---
 doc/conventions.rst | 6 +++---
 doc/intro.rst       | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/doc/conventions.rst b/doc/conventions.rst
index 21f3ab82b..fb1f0715f 100644
--- a/doc/conventions.rst
+++ b/doc/conventions.rst
@@ -289,10 +289,10 @@ element of the `control.config.defaults` dictionary::
 
     ct.config.defaults['module.parameter'] = value
 
-The `~control.config.set_defaults` function can also be used to set multiple
-configuration parameters at the same time::
+The :func:`~control.set_defaults` function can also be used to
+set multiple configuration parameters at the same time::
 
-    ct.config.set_defaults('module', param1=val1, param2=val2, ...]
+    ct.set_defaults('module', param1=val1, param2=val2, ...]
 
 Finally, there are also functions available set collections of variables based
 on standard configurations.
diff --git a/doc/intro.rst b/doc/intro.rst
index 2287bbac4..a45bb030e 100644
--- a/doc/intro.rst
+++ b/doc/intro.rst
@@ -56,7 +56,7 @@ they are not already present.
 .. note::
    Mixing packages from conda-forge and the default conda channel
    can sometimes cause problems with dependencies, so it is usually best to
-   instally NumPy, SciPy, and Matplotlib from conda-forge as well.
+   install NumPy, SciPy, and Matplotlib from conda-forge as well.
 
 To install using pip::
 

From e6fdc17258dee3c3c7df0f48ed9d3b6718c42fbb Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Sun, 10 Nov 2024 19:16:22 -0800
Subject: [PATCH 8/8] updated per @slivingston review comments

---
 control/nlsys.py            | 12 ++++++------
 control/tests/iosys_test.py |  3 +--
 examples/cruise-control.py  |  2 +-
 examples/cruise.ipynb       |  2 +-
 4 files changed, 9 insertions(+), 10 deletions(-)

diff --git a/control/nlsys.py b/control/nlsys.py
index d463ea7fc..835c16ef6 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -93,7 +93,7 @@ class NonlinearIOSystem(InputOutputSystem):
         generic name <sys[id]> is generated with a unique integer id.
 
     params : dict, optional
-       Parameter values for the system.  Passed to the evaluation functions
+        Parameter values for the system.  Passed to the evaluation functions
         for the system as default values, overriding internal defaults.
 
     See Also
@@ -1671,7 +1671,7 @@ def ivp_rhs(t, x):
         success=soln.success, message=message)
 
 
-class OperatingPoint(object):
+class OperatingPoint():
     """A class for representing the operating point of a nonlinear I/O system.
 
     The ``OperatingPoint`` class stores the operating point of a nonlinear
@@ -1693,18 +1693,18 @@ class OperatingPoint(object):
 
     """
     def __init__(
-            self, states, inputs=None, outputs=None, result=None,
+            self, states, inputs, outputs=None, result=None,
             return_outputs=False, return_result=False):
         self.states = states
         self.inputs = inputs
 
         if outputs is None and return_outputs and not return_result:
-            raise SystemError("return_outputs specified by no y0 value")
+            raise ValueError("return_outputs specified but no outputs value")
         self.outputs = outputs
         self.return_outputs = return_outputs
 
         if result is None and return_result:
-            raise SystemError("return_result specified by no result value")
+            raise ValueError("return_result specified but no result value")
         self.result = result
         self.return_result = return_result
 
@@ -1799,7 +1799,7 @@ def find_operating_point(
         values will be ignored in solving for an operating point.  State
         indices can be listed in any order.  By default, all updates will be
         fixed at `dx0` in searching for an operating point.
-    root_method : str, optonal
+    root_method : str, optional
         Method to find the operating point.  If specified, this parameter
         is passed to the :func:`scipy.optimize.root` function.
     root_kwargs : dict, optional
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index d2570be01..baaee03f6 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -2089,7 +2089,6 @@ def test_find_eqpt(x0, ix, u0, iu, y0, iy, dx0, idx, dt, x_expect, u_expect):
 
 
 # Test out new operating point version of find_eqpt
-# TODO: add return_)y tests
 def test_find_operating_point():
     dt = 1
     sys = ct.NonlinearIOSystem(
@@ -2148,7 +2147,7 @@ def test_operating_point():
     np.testing.assert_allclose(linsys_orig.C, linsys_oppt.C)
     np.testing.assert_allclose(linsys_orig.D, linsys_oppt.D)
 
-    # Call find_operating point with method and keyword arguments
+    # Call find_operating_point with method and keyword arguments
     op_point = ct.find_operating_point(
         sys, 0, 0, root_method='lm', root_kwargs={'tol': 1e-6})
 
diff --git a/examples/cruise-control.py b/examples/cruise-control.py
index b8a91eeab..074be8fa8 100644
--- a/examples/cruise-control.py
+++ b/examples/cruise-control.py
@@ -349,7 +349,7 @@ def cruise_plot(sys, t, y, label=None, t_hill=None, vref=20, antiwindup=False,
 # and u given the gear, slope, and desired output velocity)
 X0, U0, Y0 = ct.find_operating_point(
     cruise_pi, [vref[0], 0], [vref[0], gear[0], theta0[0]],
-    y0=[0, vref[0]], iu=[1, 2], iy=[1], return_y=True)
+    y0=[0, vref[0]], iu=[1, 2], iy=[1], return_outputs=True)
 
 # Now simulate the effect of a hill at t = 5 seconds
 plt.figure()
diff --git a/examples/cruise.ipynb b/examples/cruise.ipynb
index 4f1c152f9..64eafbaf7 100644
--- a/examples/cruise.ipynb
+++ b/examples/cruise.ipynb
@@ -804,7 +804,7 @@
     "# Compute the equilibrium throttle setting for the desired speed\n",
     "X0, U0, Y0 = ct.find_eqpt(\n",
     "    cruise_pi, [vref[0], 0], [vref[0], gear[0], theta0[0]],\n",
-    "    y0=[0, vref[0]], iu=[1, 2], iy=[1], return_y=True)\n",
+    "    y0=[0, vref[0]], iu=[1, 2], iy=[1], return_outputs=True)\n",
     "\n",
     "# Now simulate the effect of a hill at t = 5 seconds\n",
     "plt.figure()\n",