Reimplementation of 2D phase plots #980

murrayrm · 2024-03-28T01:50:35Z

This PR reimplements the phase plot functionality of python-control by creating a new function phase_plane_plot and creating a module phaseplot that provides more specialized functionality. The legacy phase_plot function is still available (with a deprecation warning).

From the documentation:

The default method for generating a phase plane plot is to provide a 2D dynamical system along with a range of coordinates and time limit:

sys = ct.nlsys(
    lambda t, x, u, params: np.array([[0, 1], [-1, -1]]) @ x, 
    states=['position', 'velocity'], inputs=0, name='damped oscillator')
axis_limits = [-1, 1, -1, 1]
T = 8
ct.phase_plane_plot(sys, axis_limits, T)

By default, the plot includes streamlines generated from starting points on limits of the plot, with arrows showing the flow of the system, as well as any equilibrium points for the system. A variety of options are available to modify the information that is plotted, including plotting a grid of vectors instead of streamlines and turning on and off various features of the plot.

To illustrate some of these possibilities, consider a phase plane plot for an inverted pendulum system, which is created using a mesh grid:

def invpend_update(t, x, u, params):
    m, l, b, g = params['m'], params['l'], params['b'], params['g']
    return [x[1], -b/m * x[1] + (g * l / m) * np.sin(x[0]) + u[0]/m]
invpend = ct.nlsys(invpend_update, states=2, inputs=1, name='invpend')
    
ct.phase_plane_plot(
    invpend, [-2*pi, 2*pi, -2, 2], 5,
    gridtype='meshgrid', gridspec=[5, 8], arrows=3,
    plot_equilpoints={'gridspec': [12, 9]},
    params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1})
plt.xlabel(r"$\theta$ [rad]")
plt.ylabel(r"$\dot\theta$ [rad/sec]")

This figure shows several features of more complex phase plane plots: multiple equilibrium points are shown, with saddle points showing separatrices, and streamlines generated along a 5x8 mesh of initial conditions. At each mesh point, a streamline is created that goes 5 time units forward and backward in time. A separate grid specification is used to find equilibrium points and separatrices (since the course grid spacing of 5x8 does not find all possible equilibrium points). Together, the multiple features in the phase plane plot give a good global picture of the topological strucrure of solutions of the dynamical system.

Phase plots can be buit up by hand using a variety of helper functions that are part of the control.phaseplot (pp) module:

import control.phaseplot as pp

def oscillator_update(t, x, u, params):
    return [x[1] + x[0] * (1 - x[0]**2 - x[1]**2),
            -x[0] + x[1] * (1 - x[0]**2 - x[1]**2)]
oscillator = ct.nlsys(
    oscillator_update, states=2, inputs=0, name='nonlinear oscillator')

ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9)
pp.streamlines(
    oscillator, np.array([[0, 0]]), 1.5,
    gridtype='circlegrid', gridspec=[0.5, 6], dir='both')
pp.streamlines(
    oscillator, np.array([[1, 0]]), 2*pi, arrows=6, color='b')
plt.gca().set_aspect('equal')

The following helper functions are available:

pp.equilpoints
pp.separatrices
pp.streamlines
pp.vectorfield

The phase_plane_plot function calls these helper functions based on the options it is passed.

Note that unlike other plotting function, phase plane plots do not involve computing a response and then plotting the result via a plot() method. Instead, the plot is generated directly be a call to the phase_plane_plot function (or one of the ct.phaseplot helper functions.

coveralls · 2024-03-28T01:59:38Z

coverage: 94.467% (+0.05%) from 94.422%
when pulling 748e36a on murrayrm:phaseplots_01Jan2024
into 7a70be1 on python-control:main.

slivingston · 2024-03-28T03:02:19Z

control/flatsys/__init__.py

@@ -35,7 +35,7 @@
 # Author: Richard M. Murray
 # Date: 1 Jul 2019

-r"""The :mod:`control.flatsys` package contains a set of classes and functions
+r"""The :mod:`control.flatsys` module contains a set of classes and functions


The original is correct: flatsys is a package according to standard terms: https://docs.python.org/3/glossary.html#term-package

slivingston · 2024-03-28T21:38:01Z

control/phaseplot.py

-import numpy as np
-import matplotlib.pyplot as mpl
+"""The :mod:`control.phaseplot` module contains functions for generating 2D
+phase plots.


Convention is for summary to be 1 line (https://peps.python.org/pep-0257/#multi-line-docstrings).

slivingston · 2024-03-28T21:39:44Z

examples/phase_plane_plots.py

@@ -0,0 +1,215 @@
+# phase_portraits.py - phase portrait examples


Suggested change

# phase_portraits.py - phase portrait examples

# phase_plane_plots.py - phase portrait examples

slivingston · 2024-03-28T21:45:49Z

doc/Makefile

@@ -29,6 +30,9 @@ freqplot-siso_bode-default.png: ../control/tests/freqplot_test.py
 rlocus-siso_ctime-default.png: ../control/tests/rlocus_test.py
 	PYTHONPATH=.. python $<

+phaseplot-dampedosc-default.png: ../control/tests/phaseplot_test.py
+	PYTHONPATH=.. python $<


Should the other images generated by phaseplot_test.py (phaseplot-oscillator-helpers.png, phaseplot-invpend-meshgrid.png) be included here? It does not seem important because the PNG files are being committed, but I note this here in case you wanted to.

Agree that they could be listed, but I didn't list them all out since triggering on the one file will update all files => if the test generation script is updated, all PNG files are updated.

Good point! I will complete review tomorrow.

sawyerbfuller

Fixes some very outdated code. Very nice!

sawyerbfuller · 2024-03-29T23:00:17Z

control/phaseplot.py

+        return new_kwargs
+
+    # Create array for storing outputs
+    out = np.empty(3, dtype=object)


Docstring says this is a list. Seems like a list might be better/simpler here than a numpy array, e.g. out = [[], None, None]

sawyerbfuller · 2024-03-29T23:05:26Z

control/phaseplot.py

+
+    Parameters
+    ----------
+    sys : NonlinearIOSystem or callable(x, t, ...)


What is motivation to maintain the old convention of (x, t, ...) argument order, rather than (t, x, ...) when called with a callable? Backwards compatibility? I am not aware of anywhere else we or scipy adheres to that convention anymore, but I might be missing something.

sawyerbfuller · 2024-03-29T23:10:04Z

control/phaseplot.py

+    vfdata = np.zeros((points.shape[0], 4))
+    for i, x in enumerate(points):
+        vfdata[i, :2] = x
+        vfdata[i, 2:] = sys.dynamics(0, x, 0, params=params)


I think this wouldn't work if sys is just a function - is the callable intended to be a function if it is not a system?

sawyerbfuller · 2024-03-29T23:19:24Z

control/phaseplot.py

+    Parameters
+    ----------
+    sys : NonlinearIOSystem or callable(x, t, ...)
+        Function used to generate phase plane data.


Relatedly, should there be a statement here that only two-state systems can be accepted? (and a test that gives a helpful error?)

murrayrm · 2024-03-30T05:47:11Z

In addition to adressing comments, functionality for handling callables (versus systems) was missing; now added.

sawyerbfuller · 2024-03-30T20:10:53Z

control/phaseplot.py

+    # Create reverse time system, if needed
+    if dir != 'forward':
+        revsys = NonlinearIOSystem(
+            lambda t, x, u, params: -np.array(sys.updfcn(t, x, u, params)),


Suggested change

lambda t, x, u, params: -np.array(sys.updfcn(t, x, u, params)),

lambda t, x, u, params: -np.asarray(sys.updfcn(t, x, u, params)),

This won't re-create the array if it is already one.

relatedly, and this is more of a question: I am confused about the use of each of these methods, all of which do roughly the same thing. Is this correct: _rhs fast and low-level, updfcn just a reference to the function that was passed in, and dynamics: user-friendly function that may not be as fast? Just checking that this is the right one (speed is kind of important in this application)

Correct on the various versions. In general:

_updfcn() should not be called directly (though OK here, because we are creating a new system).

_rhs() should be called inside loops (with a call to _update_params() before the loop starts, which sets the parameter values)

dynamics() is the user-callable function (which calls _update_params and then _rhs)

sawyerbfuller

LGTM

slivingston

Great work! I found a few more misprints, but nothing critical. Locally I merged this into the current tip of main branch and successfully ran the tests and example code, and built the documentation.

slivingston · 2024-03-31T02:20:54Z

examples/phase_plane_plots.py

+# RMM, 25 Mar 2024
+#
+# This file contains a number of examples of phase plane plots generated
+# using the phaseplot module.  Most of these figures lines up with examples


Suggested change

# using the phaseplot module. Most of these figures lines up with examples

# using the phaseplot module. Most of these figures line up with examples

slivingston · 2024-03-31T02:23:31Z

control/phaseplot.py

+phase plots. The base function for creating phase plane portraits is
+:func:`~control.phase_plane_plot`, which generates a phase plane portrait
+for a 2 state I/O system (with no inputs).  In addition, several other
+functions are available to creat ecustomize phase plane plots:


Suggested change

functions are available to creat ecustomize phase plane plots:

functions are available to customize phase plane plots:

slivingston · 2024-03-31T02:24:58Z

control/phaseplot.py

+    Parameters
+    ----------
+    sys : NonlinearIOSystem or callable(t, x, ...)
+        I/O systems or function used to generate phase plane data. If a


Suggested change

I/O systems or function used to generate phase plane data. If a

I/O system or function used to generate phase plane data. If a

slivingston · 2024-03-31T02:57:26Z

doc/plotting.rst

+The :func:`~control.phase_plane_plot` function calls these helper functions
+based on the options it is passed.
+
+Note that unlike other plotting function, phase plane plots do not involve


Suggested change

Note that unlike other plotting function, phase plane plots do not involve

Note that unlike other plotting functions, phase plane plots do not involve

slivingston · 2024-03-31T03:18:27Z

control/phaseplot.py

+    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.streamlines`.


Suggested change

in the dict as keywords to :func:`~control.phaseplot.streamlines`.

in the dict as keywords to :func:`~control.phaseplot.separatrices`.

slivingston · 2024-03-31T04:26:12Z

control/phaseplot.py


-import numpy as np
-import matplotlib.pyplot as mpl
+# Create a system form a callable


Suggested change

# Create a system form a callable

# Create a system from a callable

slivingston · 2024-03-31T04:41:51Z

control/phaseplot.py

+        # if not isinstance(pointdata, (list, tuple)) or len(pointdata) != 4:
+        #     raise ValueError("invalid grid data specification")


Suggested change

# if not isinstance(pointdata, (list, tuple)) or len(pointdata) != 4:

# raise ValueError("invalid grid data specification")

This can be deleted, or the check in the comment can be expanded to handle case of isinstance(pointdata, np.ndarray).

murrayrm · 2024-03-31T05:05:21Z

Thanks for the careful reviews, @sawyerbfuller and @slivingston! Will merge as soon as checks finish up (just in case).

initial reimplementation of phase plots

3ca12ff

This was linked to issues Mar 28, 2024

phase_plot ODE function signature is not compatible with solve_ivp and StateSpace.dynamics #968

Closed

Strange nondeterministic issues when using phase_plot with X, Y and X0 set #970

Closed

murrayrm mentioned this pull request Mar 28, 2024

Version 0.10 changes (tracking issue) #917

Closed

7 tasks

murrayrm added this to the 0.10.0 milestone Mar 28, 2024

slivingston reviewed Mar 28, 2024

View reviewed changes

sawyerbfuller requested changes Mar 29, 2024

View reviewed changes

address @sliivingston, @sawyerbfuller comments

c41d472

sawyerbfuller reviewed Mar 30, 2024

View reviewed changes

sawyerbfuller approved these changes Mar 30, 2024

View reviewed changes

small updates to address @sawyerbfuller comments on updfcn

edb444e

slivingston approved these changes Mar 31, 2024

View reviewed changes

fix typos + other small fixes from @slivingston

748e36a

murrayrm merged commit 8983e39 into python-control:main Mar 31, 2024
22 checks passed

murrayrm deleted the phaseplots_01Jan2024 branch March 31, 2024 05:09

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Reimplementation of 2D phase plots #980

Reimplementation of 2D phase plots #980

murrayrm commented Mar 28, 2024

coveralls commented Mar 28, 2024 •

edited

Loading

slivingston Mar 28, 2024

slivingston Mar 28, 2024

slivingston Mar 28, 2024

slivingston Mar 28, 2024

murrayrm Mar 30, 2024

slivingston Mar 30, 2024

sawyerbfuller left a comment

sawyerbfuller Mar 29, 2024

sawyerbfuller Mar 29, 2024

sawyerbfuller Mar 29, 2024

sawyerbfuller Mar 29, 2024

murrayrm commented Mar 30, 2024

sawyerbfuller Mar 30, 2024

murrayrm Mar 30, 2024

sawyerbfuller left a comment

slivingston left a comment

slivingston Mar 31, 2024

slivingston Mar 31, 2024

slivingston Mar 31, 2024

slivingston Mar 31, 2024

slivingston Mar 31, 2024

slivingston Mar 31, 2024

slivingston Mar 31, 2024

murrayrm commented Mar 31, 2024

		@@ -0,0 +1,215 @@
		# phase_portraits.py - phase portrait examples

	# phase_portraits.py - phase portrait examples
	# phase_plane_plots.py - phase portrait examples

	lambda t, x, u, params: -np.array(sys.updfcn(t, x, u, params)),
	lambda t, x, u, params: -np.asarray(sys.updfcn(t, x, u, params)),

	# using the phaseplot module. Most of these figures lines up with examples
	# using the phaseplot module. Most of these figures line up with examples

	functions are available to creat ecustomize phase plane plots:
	functions are available to customize phase plane plots:

	I/O systems or function used to generate phase plane data. If a
	I/O system or function used to generate phase plane data. If a

	in the dict as keywords to :func:`~control.phaseplot.streamlines`.
	in the dict as keywords to :func:`~control.phaseplot.separatrices`.

	# Create a system form a callable
	# Create a system from a callable

		# if not isinstance(pointdata, (list, tuple)) or len(pointdata) != 4:
		# raise ValueError("invalid grid data specification")

Reimplementation of 2D phase plots #980

Reimplementation of 2D phase plots #980

Conversation

murrayrm commented Mar 28, 2024

coveralls commented Mar 28, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sawyerbfuller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

murrayrm commented Mar 30, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sawyerbfuller left a comment

Choose a reason for hiding this comment

slivingston left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

murrayrm commented Mar 31, 2024

coveralls commented Mar 28, 2024 •

edited

Loading