Phase plot examples¶
Code¶
1# phaseplots.py - examples of phase portraits
2# RMM, 24 July 2011
3#
4# This file contains examples of phase portraits pulled from "Feedback
5# Systems" by Astrom and Murray (Princeton University Press, 2008).
6
7import os
8
9import numpy as np
10import matplotlib.pyplot as plt
11from control.phaseplot import phase_plot
12from numpy import pi
13
14# Clear out any figures that are present
15plt.close('all')
16
17#
18# Inverted pendulum
19#
20
21# Define the ODEs for a damped (inverted) pendulum
22def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
23 return x[1], -b/m*x[1] + (g*l/m)*np.sin(x[0])
24
25
26# Set up the figure the way we want it to look
27plt.figure()
28plt.clf()
29plt.axis([-2*pi, 2*pi, -2.1, 2.1])
30plt.title('Inverted pendulum')
31
32# Outer trajectories
33phase_plot(
34 invpend_ode,
35 X0=[[-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
36 [-1, 2.1], [4.2, 2.1], [5, 2.1],
37 [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
38 [1, -2.1], [-4.2, -2.1], [-5, -2.1]],
39 T=np.linspace(0, 40, 200),
40 logtime=(3, 0.7)
41)
42
43# Separatrices
44phase_plot(invpend_ode, X0=[[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0)
45
46#
47# Systems of ODEs: damped oscillator example (simulation + phase portrait)
48#
49
50def oscillator_ode(x, t, m=1., b=1, k=1):
51 return x[1], -k/m*x[0] - b/m*x[1]
52
53
54# Generate a vector plot for the damped oscillator
55plt.figure()
56plt.clf()
57phase_plot(oscillator_ode, [-1, 1, 10], [-1, 1, 10], 0.15)
58#plt.plot([0], [0], '.')
59# a=gca; set(a,'FontSize',20); set(a,'DataAspectRatio',[1,1,1])
60plt.xlabel('$x_1$')
61plt.ylabel('$x_2$')
62plt.title('Damped oscillator, vector field')
63
64# Generate a phase plot for the damped oscillator
65plt.figure()
66plt.clf()
67plt.axis([-1, 1, -1, 1]) # set(gca, 'DataAspectRatio', [1, 1, 1]);
68phase_plot(
69 oscillator_ode,
70 X0=[
71 [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
72 [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
73 ],
74 T=np.linspace(0, 8, 80),
75 timepts=[0.25, 0.8, 2, 3]
76)
77plt.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3)
78# set(gca, 'DataAspectRatio', [1,1,1])
79plt.xlabel('$x_1$')
80plt.ylabel('$x_2$')
81plt.title('Damped oscillator, vector field and stream lines')
82
83#
84# Stability definitions
85#
86# This set of plots illustrates the various types of equilibrium points.
87#
88
89
90def saddle_ode(x, t):
91 """Saddle point vector field"""
92 return x[0] - 3*x[1], -3*x[0] + x[1]
93
94
95# Asy stable
96m = 1
97b = 1
98k = 1 # default values
99plt.figure()
100plt.clf()
101plt.axis([-1, 1, -1, 1]) # set(gca, 'DataAspectRatio', [1 1 1]);
102phase_plot(
103 oscillator_ode,
104 X0=[
105 [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.7, 1], [1, 1], [1.3, 1],
106 [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.7, -1], [-1, -1],
107 [-1.3, -1]
108 ],
109 T=np.linspace(0, 10, 100),
110 timepts=[0.3, 1, 2, 3],
111 parms=(m, b, k)
112)
113plt.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3)
114# plt.set(gca,'FontSize', 16)
115plt.xlabel('$x_1$')
116plt.ylabel('$x_2$')
117plt.title('Asymptotically stable point')
118
119# Saddle
120plt.figure()
121plt.clf()
122plt.axis([-1, 1, -1, 1]) # set(gca, 'DataAspectRatio', [1 1 1])
123phase_plot(
124 saddle_ode,
125 scale=2,
126 timepts=[0.2, 0.5, 0.8],
127 X0=[
128 [-1, -1], [1, 1],
129 [-1, -0.95], [-1, -0.9], [-1, -0.8], [-1, -0.6], [-1, -0.4], [-1, -0.2],
130 [-0.95, -1], [-0.9, -1], [-0.8, -1], [-0.6, -1], [-0.4, -1], [-0.2, -1],
131 [1, 0.95], [1, 0.9], [1, 0.8], [1, 0.6], [1, 0.4], [1, 0.2],
132 [0.95, 1], [0.9, 1], [0.8, 1], [0.6, 1], [0.4, 1], [0.2, 1],
133 [-0.5, -0.45], [-0.45, -0.5], [0.5, 0.45], [0.45, 0.5],
134 [-0.04, 0.04], [0.04, -0.04]
135 ],
136 T=np.linspace(0, 2, 20)
137)
138plt.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3)
139# set(gca,'FontSize', 16)
140plt.xlabel('$x_1$')
141plt.ylabel('$x_2$')
142plt.title('Saddle point')
143
144# Stable isL
145m = 1
146b = 0
147k = 1 # zero damping
148plt.figure()
149plt.clf()
150plt.axis([-1, 1, -1, 1]) # set(gca, 'DataAspectRatio', [1 1 1]);
151phase_plot(
152 oscillator_ode,
153 timepts=[pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6,
154 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi],
155 X0=[[0.2, 0], [0.4, 0], [0.6, 0], [0.8, 0], [1, 0], [1.2, 0], [1.4, 0]],
156 T=np.linspace(0, 20, 200),
157 parms=(m, b, k)
158)
159plt.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3)
160# plt.set(gca,'FontSize', 16)
161plt.xlabel('$x_1$')
162plt.ylabel('$x_2$')
163plt.title('Undamped system\nLyapunov stable, not asympt. stable')
164
165if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
166 plt.show()
Notes¶
1. The environment variable PYCONTROL_TEST_EXAMPLES is used for testing to turn off plotting of the outputs.