add warning messages on trajectory errors (+ ability to suppress)

murrayrm · Jun 26, 2024 · Apr 20, 2024 · May 11, 2024 · May 11, 2024 · Jun 5, 2024
commit 1b71fa00ebca08af89f6f96074f89b9b1e4b96b4
diff --git a/control/nlsys.py b/control/nlsys.py
@@ -18,16 +18,17 @@
 
 """
 
-import numpy as np
-import scipy as sp
 import copy
 from warnings import warn
 
+import numpy as np
+import scipy as sp
+
 from . import config
-from .iosys import InputOutputSystem, _process_signal_list, \
-    _process_iosys_keywords, isctime, isdtime, common_timebase, _parse_spec
-from .timeresp import _check_convert_array, _process_time_response, \
-    TimeResponseData
+from .iosys import (InputOutputSystem, _parse_spec, _process_iosys_keywords,
+                    _process_signal_list, common_timebase, isctime, isdtime)
+from .timeresp import (TimeResponseData, _check_convert_array,
+                       _process_time_response)
 
 __all__ = ['NonlinearIOSystem', 'InterconnectedSystem', 'nlsys',
            'input_output_response', 'find_eqpt', 'linearize',
@@ -528,7 +529,6 @@ def linearize(self, x0, u0, t=0, params=None, eps=1e-6,
         # numerical linearization use the `_rhs()` and `_out()` member
         # functions.
         #
-
         # If x0 and u0 are specified as lists, concatenate the elements
         x0 = _concatenate_list_elements(x0, 'x0')
         u0 = _concatenate_list_elements(u0, 'u0')
@@ -1317,7 +1317,7 @@ def nlsys(
 
 
 def input_output_response(
-        sys, T, U=0., X0=0, params=None, ignore_error=False,
+        sys, T, U=0., X0=0, params=None, ignore_errors=False,
         transpose=False, return_x=False, squeeze=None,
         solve_ivp_kwargs=None, t_eval='T', **kwargs):
     """Compute the output response of a system to a given input.
@@ -1393,6 +1393,11 @@ def input_output_response(
         to 'RK45'.
     solve_ivp_kwargs : dict, optional
         Pass additional keywords to :func:`scipy.integrate.solve_ivp`.
+    ignore_errors : bool, optional
+        If ``False`` (default), errors during computation of the trajectory
+        will raise a ``RuntimeError`` exception.  If ``True``, do not raise
+        an exception and instead set ``results.success`` to ``False`` and
+        place an error message in ``results.message``.
 
     Raises
     ------
@@ -1593,8 +1598,12 @@ def ivp_rhs(t, x):
         soln = sp.integrate.solve_ivp(
             ivp_rhs, (T0, Tf), X0, t_eval=t_eval,
             vectorized=False, **solve_ivp_kwargs)
-        if not ignore_error and not soln.success:
-            raise RuntimeError("solve_ivp failed: " + soln.message)
+        if not soln.success:
+            message = "solve_ivp failed: " + soln.message
+            if not ignore_errors:
+                raise RuntimeError(message)
+        else:
+            message = None
 
         # Compute inputs and outputs for each time point
         u = np.zeros((ninputs, len(soln.t)))
@@ -1650,7 +1659,7 @@ def ivp_rhs(t, x):
         u = np.transpose(np.array(u))
 
         # Mark solution as successful
-        soln.success = True     # No way to fail
+        soln.success, message = True, None      # No way to fail
 
     else:                       # Neither ctime or dtime??
         raise TypeError("Can't determine system type")
@@ -1660,7 +1669,8 @@ def ivp_rhs(t, x):
         output_labels=sys.output_labels, input_labels=sys.input_labels,
         state_labels=sys.state_labels, sysname=sys.name,
         title="Input/output response for " + sys.name,
-        transpose=transpose, return_x=return_x, squeeze=squeeze)
+        transpose=transpose, return_x=return_x, squeeze=squeeze,
+        success=soln.success, message=message)
 
 
 def find_eqpt(sys, x0, u0=None, y0=None, t=0, params=None,
@@ -2252,7 +2262,7 @@ def interconnect(
     `outputs`, for more natural naming of SISO systems.
 
     """
-    from .statesp import StateSpace, LinearICSystem, _convert_to_statespace
+    from .statesp import LinearICSystem, StateSpace, _convert_to_statespace
     from .xferfcn import TransferFunction
 
     dt = kwargs.pop('dt', None)         # bypass normal 'dt' processing

diff --git a/control/phaseplot.py b/control/phaseplot.py
@@ -52,7 +52,7 @@
 def phase_plane_plot(
         sys, pointdata=None, timedata=None, gridtype=None, gridspec=None,
         plot_streamlines=True, plot_vectorfield=False, plot_equilpoints=True,
-        plot_separatrices=True, ax=None, **kwargs
+        plot_separatrices=True, ax=None, suppress_warnings=False, **kwargs
 ):
     """Plot phase plane diagram.
 
@@ -88,22 +88,6 @@ def phase_plane_plot(
         Parameters to pass to system. For an I/O system, `params` should be
         a dict of parameters and values. For a callable, `params` should be
         dict with key 'args' and value given by a tuple (passed to callable).
-    plot_streamlines : bool or dict
-        If `True` (default) then plot streamlines based on the pointdata
-        and gridtype.  If set to a dict, pass on the key-value pairs in
-        the dict as keywords to :func:`~control.phaseplot.streamlines`.
-    plot_vectorfield : bool or dict
-        If `True` (default) then plot the vector field based on the pointdata
-        and gridtype.  If set to a dict, pass on the key-value pairs in
-        the dict as keywords to :func:`~control.phaseplot.vectorfield`.
-    plot_equilpoints : bool or dict
-        If `True` (default) then plot equilibrium points based in the phase
-        plot boundary. If set to a dict, pass on the key-value pairs in the
-        dict as keywords to :func:`~control.phaseplot.equilpoints`.
-    plot_separatrices : bool or dict
-        If `True` (default) then plot separatrices starting from each
-        equilibrium point.  If set to a dict, pass on the key-value pairs
-        in the dict as keywords to :func:`~control.phaseplot.separatrices`.
     color : str
         Plot all elements in the given color (use `plot_<fcn>={'color': c}`
         to set the color in one element of the phase plot.
@@ -117,6 +101,27 @@ def phase_plane_plot(
         out[1] = Quiver object (vector field arrows)
         out[2] = list of Line2D objects (equilibrium points)
 
+    Other parameters
+    ----------------
+    plot_streamlines : bool or dict, optional
+        If `True` (default) then plot streamlines based on the pointdata
+        and gridtype.  If set to a dict, pass on the key-value pairs in
+        the dict as keywords to :func:`~control.phaseplot.streamlines`.
+    plot_vectorfield : bool or dict, optional
+        If `True` (default) then plot the vector field based on the pointdata
+        and gridtype.  If set to a dict, pass on the key-value pairs in
+        the dict as keywords to :func:`~control.phaseplot.vectorfield`.
+    plot_equilpoints : bool or dict, optional
+        If `True` (default) then plot equilibrium points based in the phase
+        plot boundary. If set to a dict, pass on the key-value pairs in the
+        dict as keywords to :func:`~control.phaseplot.equilpoints`.
+    plot_separatrices : bool or dict, optional
+        If `True` (default) then plot separatrices starting from each
+        equilibrium point.  If set to a dict, pass on the key-value pairs
+        in the dict as keywords to :func:`~control.phaseplot.separatrices`.
+    suppress_warnings : bool, optional
+        If set to `True`, suppress warning messages in generating trajectories.
+
     """
     # Process arguments
     params = kwargs.get('params', None)
@@ -149,7 +154,8 @@ def _create_kwargs(global_kwargs, local_kwargs, **other_kwargs):
             kwargs, plot_streamlines, gridspec=gridspec, gridtype=gridtype,
             ax=ax)
         out[0] += streamlines(
-            sys, pointdata, timedata, check_kwargs=False, **kwargs_local)
+            sys, pointdata, timedata, check_kwargs=False,
+            suppress_warnings=suppress_warnings, **kwargs_local)
 
         # Get rid of keyword arguments handled by streamlines
         for kw in ['arrows', 'arrow_size', 'arrow_style', 'color',
@@ -203,7 +209,8 @@ def _create_kwargs(global_kwargs, local_kwargs, **other_kwargs):
 
 
 def vectorfield(
-        sys, pointdata, gridspec=None, ax=None, check_kwargs=True, **kwargs):
+        sys, pointdata, gridspec=None, ax=None, suppress_warnings=False,
+        check_kwargs=True, **kwargs):
     """Plot a vector field in the phase plane.
 
     This function plots a vector field for a two-dimensional state
@@ -244,6 +251,11 @@ def vectorfield(
     -------
     out : Quiver
 
+    Other parameters
+    ----------------
+    suppress_warnings : bool, optional
+        If set to `True`, suppress warning messages in generating trajectories.
+
     """
     # Get system parameters
     params = kwargs.pop('params', None)
@@ -283,8 +295,8 @@ def vectorfield(
 
 
 def streamlines(
-        sys, pointdata, timedata=1, gridspec=None, gridtype=None,
-        dir=None, ax=None, check_kwargs=True, **kwargs):
+        sys, pointdata, timedata=1, gridspec=None, gridtype=None, dir=None,
+        ax=None, check_kwargs=True, suppress_warnings=False, **kwargs):
     """Plot stream lines in the phase plane.
 
     This function plots stream lines for a two-dimensional state space
@@ -328,6 +340,11 @@ def streamlines(
     -------
     out : list of Line2D objects
 
+    Other parameters
+    ----------------
+    suppress_warnings : bool, optional
+        If set to `True`, suppress warning messages in generating trajectories.
+
     """
     # Get system parameters
     params = kwargs.pop('params', None)
@@ -373,7 +390,8 @@ def streamlines(
         timepts = _make_timepts(timedata, i)
         traj = _create_trajectory(
             sys, revsys, timepts, X0, params, dir,
-            gridtype=gridtype, gridspec=gridspec, xlim=xlim, ylim=ylim)
+            gridtype=gridtype, gridspec=gridspec, xlim=xlim, ylim=ylim,
+            suppress_warnings=suppress_warnings)
 
         # Plot the trajectory (if there is one)
         if traj.shape[1] > 1:
@@ -465,7 +483,7 @@ def equilpoints(
 
 def separatrices(
         sys, pointdata, timedata=None, gridspec=None, ax=None,
-        check_kwargs=True, **kwargs):
+        check_kwargs=True, suppress_warnings=False, **kwargs):
     """Plot separatrices in the phase plane.
 
     This function plots separatrices for a two-dimensional state space
@@ -509,6 +527,11 @@ def separatrices(
     -------
     out : list of Line2D objects
 
+    Other parameters
+    ----------------
+    suppress_warnings : bool, optional
+        If set to `True`, suppress warning messages in generating trajectories.
+
     """
     # Get system parameters
     params = kwargs.pop('params', None)
@@ -586,13 +609,15 @@ def separatrices(
                 if evals[j].real < 0:
                     traj = _create_trajectory(
                         sys, revsys, timepts, x0, params, 'reverse',
-                        gridtype='boxgrid', xlim=xlim, ylim=ylim)
+                        gridtype='boxgrid', xlim=xlim, ylim=ylim,
+                        suppress_warnings=suppress_warnings)
                     color = stable_color
                     linestyle = '--'
                 elif evals[j].real > 0:
                     traj = _create_trajectory(
                         sys, revsys, timepts, x0, params, 'forward',
-                        gridtype='boxgrid', xlim=xlim, ylim=ylim)
+                        gridtype='boxgrid', xlim=xlim, ylim=ylim,
+                        suppress_warnings=suppress_warnings)
                     color = unstable_color
                     linestyle = '-'
 
@@ -880,17 +905,21 @@ def _get_color(kwargs, ax=None):
 
 
 def _create_trajectory(
-        sys, revsys, timepts, X0, params, dir,
+        sys, revsys, timepts, X0, params, dir, suppress_warnings=False,
         gridtype=None, gridspec=None, xlim=None, ylim=None):
     # Comput ethe forward trajectory
     if dir == 'forward' or dir == 'both':
         fwdresp = input_output_response(
-            sys, timepts, X0=X0, params=params, ignore_error=True)
+            sys, timepts, X0=X0, params=params, ignore_errors=True)
+        if not fwdresp.success and not suppress_warnings:
+            warnings.warn(f"{X0=}, {fwdresp.message}")
 
     # Compute the reverse trajectory
     if dir == 'reverse' or dir == 'both':
         revresp = input_output_response(
-            revsys, timepts, X0=X0, params=params, ignore_error=True)
+            revsys, timepts, X0=X0, params=params, ignore_errors=True)
+        if not revresp.success and not suppress_warnings:
+            warnings.warn(f"{X0=}, {revresp.message}")
 
     # Create the trace to plot
     if dir == 'forward':
@@ -1212,6 +1241,3 @@ def _find(condition):
     Private implementation of deprecated matplotlib.mlab.find
     """
     return np.nonzero(np.ravel(condition))[0]
-
-
-
diff --git a/control/tests/phaseplot_test.py b/control/tests/phaseplot_test.py
@@ -9,14 +9,16 @@
 the figures so that you can check them visually.
 """
 
+import warnings
 
 import matplotlib.pyplot as plt
 import numpy as np
-from numpy import pi
 import pytest
-from control import phase_plot
+from math import pi
+
 import control as ct
 import control.phaseplot as pp
+from control import phase_plot
 
 
 # Legacy tests
@@ -156,7 +158,22 @@ def invpend_ode(t, x, m=0, l=0, b=0, g=0):
         ct.phase_plane_plot(
             invpend_ode, [-5, 5, 2, 2], params={'stuff': (1, 1, 0.2, 1)})
 
-
+    # Warning messages for invalid solutions: nonlinear spring mass system
+    sys = ct.nlsys(
+        lambda t, x, u, params: np.array(
+            [x[1], -0.25 * (x[0] - 0.01 * x[0]**3) - 0.1 * x[1]]),
+        states=2, inputs=0)
+    with pytest.warns(UserWarning, match=r"X0=array\(.*\), solve_ivp failed"):
+        ct.phase_plane_plot(
+            sys, [-12, 12, -10, 10], 15, gridspec=[2, 9],
+            plot_separatrices=False)
+
+    # Turn warnings off
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        ct.phase_plane_plot(
+            sys, [-12, 12, -10, 10], 15, gridspec=[2, 9],
+            plot_separatrices=False, suppress_warnings=True)
 
 
 def test_basic_phase_plots(savefigs=False):