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.