allow OperatingPoint for linearize + code and documentation tweaks

murrayrm · murrayrm · commit 223acc64f92e · 2024-11-06T21:31:46.000-08:00
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
