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.