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

Notes

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