Skip to content

Commit 19edd7b

Browse files
murrayrmmurrayrm
committed
* Added first cut at phase portrait command (PhasePlot) from MATLAB
 code; needs more work to pythonize the call signatures * Added simple unit tests for margin command See ChangeLog for more detailed list of changes
1 parent 396b376 commit 19edd7b

File tree

6 files changed

+518
-0
lines changed

6 files changed

+518
-0
lines changed

‎ChangeLog

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,17 @@
1+
2011-07-24 Richard Murray <murray@malabar.local>
2+
3+
* tests/margin_test.py: added simple unit tests for margin functions
4+
(initial versions just call functions; some comparisons missing)
5+
6+
* examples/README: added missing README file
7+
8+
* examples/phaseplots.py: FBS examples for phaseplot
9+
10+
* tests/phaseplot_test.py: unit tests for phaseplot
11+
12+
* src/phaseplot.py: initial cut at phase portrait function, built
13+
from amphaseplot (Feeback Systems [FBS], Astrom and Murray, 2008)
14+
115
2011-07-15 Richard Murray <murray@malabar.local>
216

317
* tests/matlab_test.py (TestMatlab): added unittest for margin()

‎examples/README

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
This directory contains worked examples using the python-control
2+
library. Each example should work by running 'ipython -pylab' and
3+
then running the given file (make sure to have the python-control
4+
module in your path).

‎examples/phaseplots.py

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,120 @@
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+
7+
import numpy as np
8+
import matplotlib.pyplot as mpl
9+
from control.phaseplot import PhasePlot
10+
from numpy import pi
11+
12+
# Clear out any figures that are present
13+
mpl.close('all')
14+
15+
#
16+
# Inverted pendulum
17+
#
18+
19+
# Define the ODEs for a damped (inverted) pendulum
20+
def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
21+
return (x[1], -b/m*x[1] + (g*l/m) * np.sin(x[0]))
22+
23+
# Set up the figure the way we want it to look
24+
mpl.figure(); mpl.clf();
25+
mpl.axis([-2*pi, 2*pi, -2.1, 2.1]);
26+
mpl.title('Inverted pendlum')
27+
28+
# Outer trajectories
29+
PhasePlot(invpend_ode,
30+
'logtime', (3, 0.7), None,
31+
[ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
32+
[-1, 2.1], [4.2, 2.1], [5, 2.1],
33+
[2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
34+
[1, -2.1], [-4.2, -2.1], [-5, -2.1] ],
35+
np.linspace(0, 40, 200))
36+
37+
# Separatrices
38+
mpl.hold(True);
39+
PhasePlot(invpend_ode, 'auto', 0, None, [[-2.3056, 2.1], [2.3056, -2.1]], 6)
40+
mpl.show();
41+
42+
#
43+
# Systems of ODEs: damped oscillator example (simulation + phase portrait)
44+
#
45+
46+
def oscillator_ode(x, t, m=1., b=1, k=1):
47+
return (x[1], -k/m*x[0] - b/m*x[1])
48+
49+
# Generate a vector plot for the damped oscillator
50+
mpl.figure(); mpl.clf();
51+
PhasePlot(oscillator_ode, [-1, 1, 10], [-1, 1, 10], 0.15);
52+
mpl.hold(True); mpl.plot([0], [0], '.');
53+
# a=gca; set(a,'FontSize',20); set(a,'DataAspectRatio',[1,1,1]);
54+
mpl.xlabel('x1'); mpl.ylabel('x2');
55+
56+
# Generate a phase plot for the damped oscillator
57+
mpl.figure(); mpl.clf();
58+
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
59+
PhasePlot(oscillator_ode,
60+
'timepts', [0.25, 0.8, 2, 3], None, [
61+
[-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
62+
[1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
63+
], np.linspace(0, 8, 80))
64+
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
65+
# set(gca,'DataAspectRatio',[1,1,1]);
66+
mpl.xlabel('x1'); mpl.ylabel('x2');
67+
68+
mpl.show()
69+
70+
#
71+
# Stability definitions
72+
#
73+
# This set of plots illustrates the various types of equilibrium points.
74+
#
75+
76+
# Saddle point vector field
77+
def saddle_ode(x, t):
78+
return (x[0] - 3*x[1], -3*x[0] + x[1]);
79+
80+
# Asy stable
81+
m = 1; b = 1; k = 1; # default values
82+
mpl.figure(); mpl.clf();
83+
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
84+
PhasePlot(oscillator_ode, 'timepts', [0.3, 1, 2, 3], None,
85+
[[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1],
86+
[1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1],
87+
[-1.3,-1]],
88+
np.linspace(0, 10, 100), parms = (m, b, k));
89+
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
90+
# set(gca,'FontSize', 16);
91+
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
92+
93+
# Saddle
94+
mpl.figure(); mpl.clf();
95+
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
96+
PhasePlot(saddle_ode, 'timepts', [0.2, 0.5, 0.8], None,
97+
[ [-1, -1], [1, 1],
98+
[-1, -0.95], [-1, -0.9], [-1, -0.8], [-1, -0.6], [-1, -0.4], [-1, -0.2],
99+
[-0.95, -1], [-0.9, -1], [-0.8, -1], [-0.6, -1], [-0.4, -1], [-0.2, -1],
100+
[1, 0.95], [1, 0.9], [1, 0.8], [1, 0.6], [1, 0.4], [1, 0.2],
101+
[0.95, 1], [0.9, 1], [0.8, 1], [0.6, 1], [0.4, 1], [0.2, 1],
102+
[-0.5, -0.45], [-0.45, -0.5], [0.5, 0.45], [0.45, 0.5],
103+
[-0.04, 0.04], [0.04, -0.04] ], np.linspace(0, 2, 20));
104+
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
105+
# set(gca,'FontSize', 16);
106+
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
107+
108+
# Stable isL
109+
m = 1; b = 0; k = 1; # zero damping
110+
mpl.figure(); mpl.clf();
111+
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
112+
PhasePlot(oscillator_ode, 'timepts',
113+
[pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi], None,
114+
[ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
115+
np.linspace(0, 20, 200), parms = (m, b, k));
116+
mpl.hold(True); mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3);
117+
# set(gca,'FontSize', 16);
118+
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
119+
120+
mpl.show()

0 commit comments

Comments
 (0)