Phase plot examples

Code

  1# phase_plane_plots.py - phase portrait examples
  2# RMM, 25 Mar 2024
  3#
  4# This file contains a number of examples of phase plane plots generated
  5# using the phaseplot module.  Most of these figures line up with examples
  6# in FBS2e, with different display options shown as different subplots.
  7
  8import time
  9import warnings
 10from math import pi, sqrt
 11
 12import matplotlib.pyplot as plt
 13import numpy as np
 14
 15import control as ct
 16import control.phaseplot as pp
 17
 18#
 19# Example 1: Dampled oscillator systems
 20#
 21
 22# Oscillator parameters
 23damposc_params = {'m': 1, 'b': 1, 'k': 1}
 24
 25# System model (as ODE)
 26def damposc_update(t, x, u, params):
 27    m, b, k = params['m'], params['b'], params['k']
 28    return np.array([x[1], -k/m * x[0] - b/m * x[1]])
 29damposc = ct.nlsys(damposc_update, states=2, inputs=0, params=damposc_params)
 30
 31fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
 32fig.set_tight_layout(True)
 33plt.suptitle("FBS Figure 5.3: damped oscillator")
 34
 35ct.phase_plane_plot(damposc, [-1, 1, -1, 1], 8, ax=ax1)
 36ax1.set_title("boxgrid [-1, 1, -1, 1], 8")
 37
 38ct.phase_plane_plot(damposc, [-1, 1, -1, 1], ax=ax2, gridtype='meshgrid')
 39ax2.set_title("meshgrid [-1, 1, -1, 1]")
 40
 41ct.phase_plane_plot(
 42    damposc, [-1, 1, -1, 1], 4, ax=ax3, gridtype='circlegrid', dir='both')
 43ax3.set_title("circlegrid [0, 0, 1], 4, both")
 44
 45ct.phase_plane_plot(
 46    damposc, [-1, 1, -1, 1], ax=ax4, gridtype='circlegrid',
 47    dir='reverse', gridspec=[0.1, 12], timedata=5)
 48ax4.set_title("circlegrid [0, 0, 0.1], reverse")
 49
 50#
 51# Example 2: Inverted pendulum
 52#
 53
 54def invpend_update(t, x, u, params):
 55    m, l, b, g = params['m'], params['l'], params['b'], params['g']
 56    return [x[1], -b/m * x[1] + (g * l / m) * np.sin(x[0])]
 57invpend = ct.nlsys(
 58    invpend_update, states=2, inputs=0,
 59    params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1})
 60
 61fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
 62fig.set_tight_layout(True)
 63plt.suptitle("FBS Figure 5.4: inverted pendulum")
 64
 65ct.phase_plane_plot(
 66    invpend, [-2*pi, 2*pi, -2, 2], 5, ax=ax1)
 67ax1.set_title("default, 5")
 68
 69ct.phase_plane_plot(
 70    invpend, [-2*pi, 2*pi, -2, 2], gridtype='meshgrid', ax=ax2)
 71ax2.set_title("meshgrid")
 72
 73ct.phase_plane_plot(
 74    invpend, [-2*pi, 2*pi, -2, 2], 1, gridtype='meshgrid',
 75    gridspec=[12, 9], ax=ax3, arrows=1)
 76ax3.set_title("denser grid")
 77
 78ct.phase_plane_plot(
 79    invpend, [-2*pi, 2*pi, -2, 2], 4, gridspec=[6, 6],
 80    plot_separatrices={'timedata': 20, 'arrows': 4}, ax=ax4)
 81ax4.set_title("custom")
 82
 83#
 84# Example 3: Limit cycle (nonlinear oscillator)
 85#
 86
 87def oscillator_update(t, x, u, params):
 88    return [
 89        x[1] + x[0] * (1 - x[0]**2 - x[1]**2),
 90        -x[0] + x[1] * (1 - x[0]**2 - x[1]**2)
 91    ]
 92oscillator = ct.nlsys(oscillator_update, states=2, inputs=0)
 93
 94fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
 95fig.set_tight_layout(True)
 96plt.suptitle("FBS Figure 5.5: Nonlinear oscillator")
 97
 98ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 3, ax=ax1)
 99ax1.set_title("default, 3")
100ax1.set_aspect('equal')
101
102try:
103    ct.phase_plane_plot(
104        oscillator, [-1.5, 1.5, -1.5, 1.5], 1, gridtype='meshgrid',
105        dir='forward', ax=ax2)
106except RuntimeError as inst:
107    axs[0,1].text(0, 0, "Runtime Error")
108    warnings.warn(inst.__str__())
109ax2.set_title("meshgrid, forward, 0.5")
110ax2.set_aspect('equal')
111
112ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], ax=ax3)
113pp.streamlines(
114    oscillator, [-0.5, 0.5, -0.5, 0.5], dir='both', ax=ax3)
115ax3.set_title("outer + inner")
116ax3.set_aspect('equal')
117
118ct.phase_plane_plot(
119    oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9, ax=ax4)
120pp.streamlines(
121    oscillator, np.array([[0, 0]]), 1.5,
122    gridtype='circlegrid', gridspec=[0.5, 6], dir='both', ax=ax4)
123pp.streamlines(
124    oscillator, np.array([[1, 0]]), 2*pi, arrows=6, ax=ax4, color='b')
125ax4.set_title("custom")
126ax4.set_aspect('equal')
127
128#
129# Example 4: Simple saddle
130#
131
132def saddle_update(t, x, u, params):
133    return [x[0] - 3*x[1], -3*x[0] + x[1]]
134saddle = ct.nlsys(saddle_update, states=2, inputs=0)
135
136fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
137fig.set_tight_layout(True)
138plt.suptitle("FBS Figure 5.9: Saddle")
139
140ct.phase_plane_plot(saddle, [-1, 1, -1, 1], ax=ax1)
141ax1.set_title("default")
142
143ct.phase_plane_plot(
144    saddle, [-1, 1, -1, 1], 0.5, gridtype='meshgrid', ax=ax2)
145ax2.set_title("meshgrid")
146
147ct.phase_plane_plot(
148    saddle, [-1, 1, -1, 1], gridspec=[16, 12], ax=ax3, 
149    plot_vectorfield=True, plot_streamlines=False, plot_separatrices=False)
150ax3.set_title("vectorfield")
151
152ct.phase_plane_plot(
153    saddle, [-1, 1, -1, 1], 0.3,
154    gridtype='meshgrid', gridspec=[5, 7], ax=ax4)
155ax3.set_title("custom")
156
157#
158# Example 5: Internet congestion control
159#
160
161def _congctrl_update(t, x, u, params):
162    # Number of sources per state of the simulation
163    M = x.size - 1                      # general case
164    assert M == 1                       # make sure nothing funny happens here
165
166    # Remaining parameters
167    N = params.get('N', M)              # number of sources
168    rho = params.get('rho', 2e-4)       # RED parameter = pbar / (bupper-blower)
169    c = params.get('c', 10)             # link capacity (Mp/ms)
170
171    # Compute the derivative (last state = bdot)
172    return np.append(
173        c / x[M] - (rho * c) * (1 + (x[:-1]**2) / 2),
174        N/M * np.sum(x[:-1]) * c / x[M] - c)
175congctrl = ct.nlsys(
176    _congctrl_update, states=2, inputs=0,
177    params={'N': 60, 'rho': 2e-4, 'c': 10})
178
179fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
180fig.set_tight_layout(True)
181plt.suptitle("FBS Figure 5.10: Congestion control")
182
183try:
184    ct.phase_plane_plot(
185        congctrl, [0, 10, 100, 500], 120, ax=ax1)
186except RuntimeError as inst:
187    ax1.text(5, 250, "Runtime Error")
188    warnings.warn(inst.__str__())
189ax1.set_title("default, T=120")
190
191try:
192    ct.phase_plane_plot(
193        congctrl, [0, 10, 100, 500], 120,
194        params={'rho': 4e-4, 'c': 20}, ax=ax2)
195except RuntimeError as inst:
196    ax2.text(5, 250, "Runtime Error")
197    warnings.warn(inst.__str__())
198ax2.set_title("updated param")
199
200ct.phase_plane_plot(
201    congctrl, [0, 10, 100, 500], ax=ax3,
202    plot_vectorfield=True, plot_streamlines=False)
203ax3.set_title("vector field")
204
205ct.phase_plane_plot(
206    congctrl, [2, 6, 200, 300], 100,
207    params={'rho': 4e-4, 'c': 20},
208    ax=ax4, plot_vectorfield={'gridspec': [12, 9]})
209ax4.set_title("vector field + streamlines")
210
211#
212# End of examples
213#
214
215plt.show(block=False)

Notes

1. The environment variable PYCONTROL_TEST_EXAMPLES is used for testing to turn off plotting of the outputs.