From b65e71b0de055005e547f2632058ab28f4cde3e7 Mon Sep 17 00:00:00 2001 From: Joey Gould Date: Mon, 29 Jun 2026 17:19:08 -0700 Subject: [PATCH] added examples from 2nd edition --- examples/bicycle.py | 169 +++++++++++++++++ examples/congctrl.py | 39 ++++ examples/cruise.py | 124 ++++++++++++ examples/example-3.15-queuing_systems.py | 84 +++++++++ examples/example-3.17-consensus.py | 61 ++++++ examples/example-3.18-repressilator.py | 126 +++++++++++++ examples/example-3.19-fitzhugh_nagumo.py | 75 ++++++++ examples/example-3.4-predator_prey.py | 61 ++++++ examples/example-5.13-genetic_switch.py | 139 ++++++++++++++ examples/example-5.14-invpend_stabilized.py | 96 ++++++++++ examples/example-5.16-predprey_bif.py | 144 ++++++++++++++ examples/example-5.17-bicycle_stability.py | 153 +++++++++++++++ examples/example-5.18-noise_cancel.py | 100 ++++++++++ examples/example-5.8-tanker_stabilty.py | 106 +++++++++++ examples/example-5.9-limit_cycle_stability.py | 79 ++++++++ examples/example-6.10-compartment_response.py | 85 +++++++++ examples/example-6.8-opamp_bandpass.py | 53 ++++++ examples/example-6.9-afm_freqresp.py | 113 +++++++++++ examples/example-7.4-steering_place.py | 122 ++++++++++++ examples/example-8.10-steering_gainsched.py | 178 ++++++++++++++++++ examples/example-template.py | 34 ++++ examples/fbs.py | 29 +++ examples/figure-1.11-cruise_robustness.py | 58 ++++++ examples/figure-1.18-airfuel_selectors.py | 97 ++++++++++ examples/figure-2.11-2dof_stepresp.py | 60 ++++++ examples/figure-2.12,14-static_nlsys.py | 136 +++++++++++++ examples/figure-2.19-posfbk_saturation.py | 115 +++++++++++ examples/figure-2.8-PI_step_responses.py | 88 +++++++++ examples/figure-2.9-secord_stepresp.py | 93 +++++++++ .../figure-3.11-spring_mass_simulation.py | 86 +++++++++ examples/figure-3.12-frequency_response.py | 106 +++++++++++ examples/figure-3.2-state_model.py | 99 ++++++++++ examples/figure-3.4-io_response.py | 65 +++++++ examples/figure-4.12-congctrl_eqplot.py | 36 ++++ examples/figure-4.13-congctrl_tcpsim.py | 64 +++++++ examples/figure-4.20-predprey_ctstime.py | 84 +++++++++ examples/figure-5-6-lyapunov_stability.py | 30 +++ examples/figure-5.1-damposc_response.py | 31 +++ examples/figure-5.10-congctrl_dynamics.py | 47 +++++ examples/figure-5.11-invpend_linearized.py | 30 +++ examples/figure-5.3-phase_portraits.py | 33 ++++ examples/figure-5.4-invpend_phaseplot.py | 21 +++ examples/figure-5.5-limit_cycle.py | 38 ++++ examples/figure-5.7-stable_eqpt.py | 39 ++++ examples/figure-5.8-asystable_eqpt.py | 28 +++ examples/figure-5.9-unstable_eqpt.py | 29 +++ examples/figure-6.1-superposition.py | 49 +++++ examples/figure-6.14-cruise_linearized.py | 117 ++++++++++++ examples/figure-6.5-modes.py | 61 ++++++ examples/figure-7.7-predprey_place.py | 61 ++++++ examples/predprey.py | 30 +++ examples/springmass.py | 11 ++ examples/steering.py | 68 +++++++ examples/template.py | 17 ++ 54 files changed, 4097 insertions(+) create mode 100644 examples/bicycle.py create mode 100644 examples/congctrl.py create mode 100644 examples/cruise.py create mode 100644 examples/example-3.15-queuing_systems.py create mode 100644 examples/example-3.17-consensus.py create mode 100644 examples/example-3.18-repressilator.py create mode 100644 examples/example-3.19-fitzhugh_nagumo.py create mode 100644 examples/example-3.4-predator_prey.py create mode 100644 examples/example-5.13-genetic_switch.py create mode 100644 examples/example-5.14-invpend_stabilized.py create mode 100644 examples/example-5.16-predprey_bif.py create mode 100644 examples/example-5.17-bicycle_stability.py create mode 100644 examples/example-5.18-noise_cancel.py create mode 100644 examples/example-5.8-tanker_stabilty.py create mode 100644 examples/example-5.9-limit_cycle_stability.py create mode 100644 examples/example-6.10-compartment_response.py create mode 100644 examples/example-6.8-opamp_bandpass.py create mode 100644 examples/example-6.9-afm_freqresp.py create mode 100644 examples/example-7.4-steering_place.py create mode 100644 examples/example-8.10-steering_gainsched.py create mode 100644 examples/example-template.py create mode 100644 examples/fbs.py create mode 100644 examples/figure-1.11-cruise_robustness.py create mode 100644 examples/figure-1.18-airfuel_selectors.py create mode 100644 examples/figure-2.11-2dof_stepresp.py create mode 100644 examples/figure-2.12,14-static_nlsys.py create mode 100644 examples/figure-2.19-posfbk_saturation.py create mode 100644 examples/figure-2.8-PI_step_responses.py create mode 100644 examples/figure-2.9-secord_stepresp.py create mode 100644 examples/figure-3.11-spring_mass_simulation.py create mode 100644 examples/figure-3.12-frequency_response.py create mode 100644 examples/figure-3.2-state_model.py create mode 100644 examples/figure-3.4-io_response.py create mode 100644 examples/figure-4.12-congctrl_eqplot.py create mode 100644 examples/figure-4.13-congctrl_tcpsim.py create mode 100644 examples/figure-4.20-predprey_ctstime.py create mode 100644 examples/figure-5-6-lyapunov_stability.py create mode 100644 examples/figure-5.1-damposc_response.py create mode 100644 examples/figure-5.10-congctrl_dynamics.py create mode 100644 examples/figure-5.11-invpend_linearized.py create mode 100644 examples/figure-5.3-phase_portraits.py create mode 100644 examples/figure-5.4-invpend_phaseplot.py create mode 100644 examples/figure-5.5-limit_cycle.py create mode 100644 examples/figure-5.7-stable_eqpt.py create mode 100644 examples/figure-5.8-asystable_eqpt.py create mode 100644 examples/figure-5.9-unstable_eqpt.py create mode 100644 examples/figure-6.1-superposition.py create mode 100644 examples/figure-6.14-cruise_linearized.py create mode 100644 examples/figure-6.5-modes.py create mode 100644 examples/figure-7.7-predprey_place.py create mode 100644 examples/predprey.py create mode 100644 examples/springmass.py create mode 100644 examples/steering.py create mode 100644 examples/template.py diff --git a/examples/bicycle.py b/examples/bicycle.py new file mode 100644 index 000000000..ce156104a --- /dev/null +++ b/examples/bicycle.py @@ -0,0 +1,169 @@ +# bicycle.py - bicycle dynamics +# RMM, 24 Nov 2024 (from KJA) +# +# Bicyle dynamics +# +# This model describes the dynamics of a bicycle with the feature that one +# of its key properties is due to a feedback mechanism that is created by +# the design of the front fork. This model is described in more detail in +# Section 4.2 of FBS2e. +# +# MATLAB header (bicycle_stabplot.m) +# % Linearized 4th order model and analysis of eigenvalues +# % Equations based on Schwab et al 2004 +# % Run bicycleparameters first +# % kja 040611 +# % Parameters of a bicycle model +# % kja 040613 +# % Basic data is given by 26 parameters + +import numpy as np +import control as ct +from math import pi + +import numpy as np + +# Acceleration of gravity [m/s^2] +g = 9.81 + +# Wheel base [m] +b = 1.00 + +# Trail [m] +c = 0.08 + +# Wheel radii +Rrw = 0.35 +Rfw = 0.35 + +# Head angle [radians] +lambda_angle = np.pi * 70 / 180 + +# Rear frame mass [kg], center of mass [m], and inertia tensor [kgm^2] +mrf = 87 +xrf = 0.491586 +zrf = 1.028138 +Jxxrf = 3.283666 +Jxzrf = 0.602765 +Jyyrf = 3.8795952 +Jzzrf = 0.565929 + +# Front frame mass [kg], center of mass [m], and inertia tensor [kgm^2] +mff = 2 +xff = 0.866 +zff = 0.676 +Jxxff = 0.08 +Jxzff = -0.02 +Jyyff = 0.07 +Jzzff = 0.02 + +# Rear wheel mass [kg], center of mass [m], and inertia tensor [kgm^2] +mrw = 1.5 +Jxxrw = 0.07 +Jyyrw = 0.14 + +# Front wheel mass [kg], center of mass [m], and inertia tensor [kgm^2] +mfw = 1.5 +Jxxfw = 0.07 +Jyyfw = 0.14 + +# Auxiliary variables +xrw = 0 +zrw = Rrw +xfw = b +zfw = Rfw +Jzzrw = Jxxrw +Jzzfw = Jxxfw + +# Total mass +mt = mrf + mrw + mff + mfw + +# Center of mass +xt = (mrf * xrf + mrw * xrw + mff * xff + mfw * xfw) / mt +zt = (mrf * zrf + mrw * zrw + mff * zff + mfw * zfw) / mt + +# Inertia tensor components +Jxxt = ( + Jxxrf + mrf * zrf**2 + + Jxxrw + mrw * zrw**2 + + Jxxff + mff * zff**2 + + Jxxfw + mfw * zfw**2 +) +Jxzt = ( + Jxzrf + mrf * xrf * zrf + + mrw * xrw * zrw + + Jxzff + mff * xff * zff + + mfw * xfw * zfw +) +Jzzt = ( + Jzzrf + mrf * xrf**2 + + Jzzrw + mrw * xrw**2 + + Jzzff + mff * xff**2 + + Jzzfw + mfw * xfw**2 +) + +# Front frame parameters +mf = mff + mfw +xf = (mff * xff + mfw * xfw) / mf +zf = (mff * zff + mfw * zfw) / mf + +Jxxf = ( + Jxxff + mff * (zff - zf)**2 + + Jxxfw + mfw * (zfw - zf)**2 +) +Jxzf = ( + Jxzff + mff * (xff - xf) * (zff - zf) + + mfw * (xfw - xf) * (zfw - zf) +) +Jzzf = ( + Jzzff + mff * (xff - xf)**2 + + Jzzfw + mfw * (xfw - xf)**2 +) + +# Auxiliary variables +d = (xf - b - c) * np.sin(lambda_angle) + zf * np.cos(lambda_angle) +Fll = ( + mf * d**2 + + Jxxf * np.cos(lambda_angle)**2 + + 2 * Jxzf * np.sin(lambda_angle) * np.cos(lambda_angle) + + Jzzf * np.sin(lambda_angle)**2 +) +Flx = mf * d * zf + Jxxf * np.cos(lambda_angle) + Jxzf * np.sin(lambda_angle) +Flz = mf * d * xf + Jxzf * np.cos(lambda_angle) + Jzzf * np.sin(lambda_angle) +gamma = c * np.sin(lambda_angle) / b +Sr = Jyyrw / Rrw +Sf = Jyyfw / Rfw +St = Sr + Sf +Su = mf * d + gamma * mt * xt + +# Matrices for the linearized fourth-order model +M = np.array([ + [Jxxt, -Flx - gamma * Jxzt], + [-Flx - gamma * Jxzt, Fll + 2 * gamma * Flz + gamma**2 * Jzzt] +]) + +K0 = np.array([ + [-mt * g * zt, g * Su], + [g * Su, -g * Su * np.cos(lambda_angle)] +]) + +K2 = np.array([ + [0, -(St + mt * zt) * np.sin(lambda_angle) / b], + [0, (Su + Sf * np.cos(lambda_angle)) * np.sin(lambda_angle) / b] +]) + +c12 = gamma * St + Sf * np.sin(lambda_angle) \ + + Jxzt * np.sin(lambda_angle) / b + gamma * mt * zt +c22 = Flz * np.sin(lambda_angle) / b \ + + gamma * (Su + Jzzt * np.sin(lambda_angle) / b) + +C = np.array([ + [0, -c12], + [gamma * St + Sf * np.sin(lambda_angle), c22] +]) + +def whipple_A(v0): + return np.block([ + [np.zeros((2, 2)), np.eye(2)], + [-np.linalg.inv(M) @ (K0 + K2 * v0**2), -np.linalg.inv(M) @ C * v0] + ]) diff --git a/examples/congctrl.py b/examples/congctrl.py new file mode 100644 index 000000000..790cd37a1 --- /dev/null +++ b/examples/congctrl.py @@ -0,0 +1,39 @@ +# congctrl.py - Congestion control dynamics and control +# RMM, 11 Jul 2006 (from MATLAB) +# +# Congestion control dynamics +# +# This model implements the congestion control dynamics described in +# the text. We assume that we have N identical sources and 1 router. +# +# To allow the model to be used in a variety of ways, we +# allow the number of states and the number of sources to be set +# independently. The length of the state vector, M+1, is used to determine +# the number of simulated sources, with each source representing N/M +# sources. + +import numpy as np +import control as ct + +# Congestion control dynamics +def _congctrl_update(t, x, u, params): + # Number of sources per state of the simulation + M = x.size - 1 + + # Remaining parameters + N = params.get('N', M) # number of sources + rho = params.get('rho', 2e-4) # RED parameter = pbar / (bupper-blower) + c = params.get('c', 10) # link capacity (Mp/ms) + + # Compute the derivative (last state = bdot) + return np.append( + c / x[M] - (rho * c) * (1 + (x[:-1]**2) / 2), + N/M * np.sum(x[:-1]) * c / x[M] - c) + + +# Function to define an I/O system +def create_iosystem(M, N=60, rho=2e-4, c=10): + #! TODO: Check to make sure M and N are compatible + return ct.nlsys( + _congctrl_update, None, states=M+1, + params={'N': N, 'rho': rho, 'c': c}) diff --git a/examples/cruise.py b/examples/cruise.py new file mode 100644 index 000000000..ba63383db --- /dev/null +++ b/examples/cruise.py @@ -0,0 +1,124 @@ +# cruise.py - Cruise control dynamics and control +# RMM, 20 Jun 2021 + +import numpy as np +import matplotlib.pyplot as plt +from math import pi +import control as ct + +# +# Vehicle model +# +# The dynamics for this system are implemented using the I/O systems module +# in python-control. This model is described in detail in Section 4.1 of +# FBS2e. +# + +# Engine model +def motor_torque(omega, params={}): + # Set up the system parameters + Tm = params.get('Tm', 190.) # engine torque constant + omega_m = params.get('omega_m', 420.) # peak engine angular speed + beta = params.get('beta', 0.4) # peak engine rolloff + + return np.clip(Tm * (1 - beta * (omega/omega_m - 1)**2), 0, None) + + +# Vehicle dynamics +def vehicle_update(t, x, u, params={}): + """Vehicle dynamics for cruise control system. + + Parameters + ---------- + x : array + System state: car velocity in m/s + u : array + System input: [throttle, gear, road_slope], where throttle is + a float between 0 and 1, gear is an integer between 1 and 5, + and road_slope is in rad. + + Returns + ------- + float + Vehicle acceleration + + """ + from math import copysign, sin + sign = lambda x: copysign(1, x) # define the sign() function + + # Set up the system parameters + m = params.get('m', 1600.) # vehicle mass, kg + g = params.get('g', 9.8) # gravitational constant, m/s^2 + Cr = params.get('Cr', 0.01) # coefficient of rolling friction + Cd = params.get('Cd', 0.32) # drag coefficient + rho = params.get('rho', 1.3) # density of air, kg/m^3 + A = params.get('A', 2.4) # car area, m^2 + alpha = params.get( + 'alpha', [40, 25, 16, 12, 10]) # gear ratio / wheel radius + + # Define variables for vehicle state and inputs + v = x[0] # vehicle velocity + throttle = np.clip(u[0], 0, 1) # vehicle throttle + gear = u[1] # vehicle gear + theta = u[2] # road slope + + # Force generated by the engine + + omega = alpha[int(gear)-1] * v # engine angular speed + F = alpha[int(gear)-1] * motor_torque(omega, params) * throttle + + # Disturbance forces + # + # The disturbance force Fd has three major components: Fg, the forces due + # to gravity; Fr, the forces due to rolling friction; and Fa, the + # aerodynamic drag. + + # Letting the slope of the road be \theta (theta), gravity gives the + # force Fg = m g sin \theta. + + Fg = m * g * sin(theta) + + # A simple model of rolling friction is Fr = m g Cr sgn(v), where Cr is + # the coefficient of rolling friction and sgn(v) is the sign of v (±1) or + # zero if v = 0. + + Fr = m * g * Cr * sign(v) + + # The aerodynamic drag is proportional to the square of the speed: Fa = + # 1/2 \rho Cd A |v| v, where \rho is the density of air, Cd is the + # shape-dependent aerodynamic drag coefficient, and A is the frontal area + # of the car. + + Fa = 1/2 * rho * Cd * A * abs(v) * v + + # Final acceleration on the car + Fd = Fg + Fr + Fa + dv = (F - Fd) / m + + return dv + +# Vehicle input/output model +vehicle_dynamics = ct.nlsys( + vehicle_update, None, name='vehicle', + inputs = ['u', 'gear', 'theta'], outputs = ['v'], states=['v']) + +# +# PI controller +# +# The controller for the system is designed in Example XX, but is included +# here since it is needed for some earlier examples. +# + +# Nominal controller design for remaining analyses +# Construct a PI controller with rolloff, as a transfer function +Kp = 0.5 # proportional gain +Ki = 0.1 # integral gain +PI_control = ct.tf( + [Kp, Ki], [1, 0.01*Ki/Kp], name='control', inputs='u', outputs='y') + +cruise_PI = ct.interconnect( + (vehicle_dynamics, PI_control), name='cruise', + connections = [['control.u', '-vehicle.v'], ['vehicle.u', 'control.y']], + inplist = ['control.u', 'vehicle.gear', 'vehicle.theta'], + inputs = ['vref', 'gear', 'theta'], + outlist = ['vehicle.v', 'vehicle.u'], outputs = ['v', 'u']) diff --git a/examples/example-3.15-queuing_systems.py b/examples/example-3.15-queuing_systems.py new file mode 100644 index 000000000..93c9e18fa --- /dev/null +++ b/examples/example-3.15-queuing_systems.py @@ -0,0 +1,84 @@ +# example-3.15-queuing_systems.py - Queuing system modeling +# RMM, 29 Aug 2021 +# +# Figure 3.22: Queuing dynamics. (a) The steady-state queue length as a +# function of $\lambda/\mu_{max}$. (b) The behavior of the queue length when +# there is a temporary overload in the system. The solid line shows a +# realization of an event-based simulation, and the dashed line shows the +# behavior of the flow model (3.33). The maximum service rate is $\mu_{max} +# = 1$, and the arrival rate starts at $\lambda = 0.5$. The arrival rate is +# increased to $\lambda = 4$ at time 20, and it returns to $\lambda =0.5$ at +# time 25. +# + +import control as ct +import numpy as np +import matplotlib.pyplot as plt + +# Queing parameters + +# Queuing system model (KJA, 2006) +def queuing_model(t, x, u, params={}): + # Define default parameters + mu = params.get('mu', 1) + + # Get the current load + lambda_ = u + + # Return the change in queue size + return np.array(lambda_ - mu * x[0] / (1 + x[0])) + +# Create I/O system representation +queuing_sys = ct.nlsys( + updfcn=queuing_model, inputs=1, outputs=1, states=1) + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(3, 2) + +# +# (a) The steady-state queue length as a function of $\lambda/\mu_{max}$. +# + +fig.add_subplot(gs[0, 0]) # first row, first column + +# Steady state queue length +x = np.linspace(0.01, 0.99, 100) +plt.plot(x, x / (1 - x), 'b-') + +# Label the plot +plt.xlabel(r"Service rate excess $\lambda/\mu_{max}$") +plt.ylabel(r"Queue length $x_{e}$") +plt.title("Steady-state queue length") + +# +# (b) The behavior of the queue length when there is a temporary overload +# in the system. The solid line shows a realization of an event-based +# simulation, and the dashed line shows the behavior of the flow model +# (3.33). The maximum service rate is $\mu_{max} = 1$, and the arrival +# rate starts at $\lambda = 0.5$. The arrival rate is increased to $\lambda +# = 4$ at time 20, and it returns to $\lambda =0.5$ at time 25. +# + +fig.add_subplot(gs[0, 1]) # first row, first column + +# Construct the loading condition +t = np.linspace(0, 80, 100) + +u =np.ones_like(t) * 0.5 +u[t <= 25] = 4 +u[t < 20] = 0.5 + +# Simulate the system dynamics +response = ct.input_output_response(queuing_sys, t, u) + +# Plot the results +plt.plot(response.time, response.outputs, 'b-') + +# Label the plot +plt.xlabel("Time $t$ [s]") +plt.ylabel(r"Queue length $x_{e}$") +plt.title("Overload condition") + +# Save the figure +plt.savefig("figure-3.22-queuing_dynamics.png", bbox_inches='tight') diff --git a/examples/example-3.17-consensus.py b/examples/example-3.17-consensus.py new file mode 100644 index 000000000..8cfe87edc --- /dev/null +++ b/examples/example-3.17-consensus.py @@ -0,0 +1,61 @@ +# example-3.17-consensus.py - consensus protocols +# RMM, 30 Aug 2021 + +# Figure 3.24: Consensus protocols for sensor networks. (a) A simple sensor +# net- work with five nodes. In this network, node 1 communicates with node +# 2 and node 2 communicates with nodes 1, 3, 4, 5, etc. (b) A simulation +# demonstrating the convergence of the consensus protocol (3.35) to the +# average value of the initial conditions. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt + +# +# System dynamcis +# + +# Construct the Laplacian corresponding to our network +L = np.array([ + [ 1, -1, 0, 0, 0], + [-1, 4, -1, -1, -1], + [ 0, -1, 2, -1, 0], + [ 0, -1, -1, 2, 0], + [ 0, -1, 0, 0, 1] +]) + +# Now generate the discrete time dynamics matrix +gamma = 0.1 +A = np.eye(5) - gamma * L + +# Initial set of measurements for the system +x0 = [10, 15, 25, 35, 40] + +# Create a discrete time system +sys = ct.ss(A, np.zeros(5), np.zeros(5), 0, dt=True) + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(3, 2) + +# +# (b) A simulation demonstrating the convergence of the consensus protocol +# (3.35) to the average value of the initial conditions. +# + +fig.add_subplot(gs[0, 1]) # first row, first column + +# Simulate the system +response = ct.initial_response(sys, 40, x0) + +# Plot th3 results +for i in range(response.nstates): + plt.plot(response.time, response.states[i], 'b-') + +# Label the figure +plt.xlabel("Iteration") +plt.ylabel("Agent states $x_i$") +plt.title("Consensus convergence") + +# Save the figure +plt.savefig("figure-3.24-consensus_dynamics.png", bbox_inches='tight') diff --git a/examples/example-3.18-repressilator.py b/examples/example-3.18-repressilator.py new file mode 100644 index 000000000..a4993b201 --- /dev/null +++ b/examples/example-3.18-repressilator.py @@ -0,0 +1,126 @@ +# example-3.18-repressilator.py - Transcriptional regulation +# RMM, 29 Aug 2021 +# +# Figure 3.26: The repressilator genetic regulatory network. (a) A schematic +# diagram of the repressilator, showing the layout of the genes in the +# plasmid that holds the circuit as well as the circuit diagram +# (center). (b) A simulation of a simple model for the repressilator, +# showing the oscillation of the individual protein concentrations. +# + +import control as ct +import numpy as np +import matplotlib.pyplot as plt + +# +# Repressilator dynamics +# +# This function implements the basic model of the repressilator All +# parameter values were taken from Nature. 2000 Jan 20; 403(6767):335-8. +# +# This model was developed by members of the 2003 Synthetic Biology Class +# on Engineered Blinkers. +# + +# Dynamics for the repressilator +def repressilator(t, x, u, params): + # store the state variables under more meaningful names + mRNA_cI = x[0] + mRNA_lacI = x[1] + mRNA_tetR = x[2] + protein_cI = x[3] + protein_lacI = x[4] + protein_tetR = x[5] + + # + # set the parameter values + # + + # set the max transcription rate in transcripts per second + k_transcription_cI = params.get('k_transcription_cI', 0.5) + k_transcription_lacI = params.get('k_transcription_lacI', 0.5) + k_transcription_tetR = params.get('k_transcription_tetR', 0.5) + + # set the leakage transcription rate (ie transcription rate if + # promoter region bound by repressor) in transcripts per second + k_transcription_leakage = params.get('k_transcription_leakage', 5e-4) + + # Set the mRNA and protein degradation rates (per second) + mRNA_half_life = params.get('mRNA_half_life', 120) # in seconds + k_mRNA_degradation = np.log(2)/mRNA_half_life + protein_half_life = params.get('protein_half_life', 600) # in seconds + k_protein_degradation = np.log(2)/protein_half_life + + # proteins per transcript lifespan + translation_efficiency = params.get('translation_efficiency', 20) + average_mRNA_lifespan = 1/k_mRNA_degradation + + # proteins per transcript per sec + k_translation = translation_efficiency/average_mRNA_lifespan + + # set the Hill coefficients of the repressors + n_tetR = params.get('n_tetR', 2) + n_cI = params.get('n_cI', 2) + n_lacI = params.get('n_lacI', 2) + + # Set the dissociation constant for the repressors to their target promoters + # in per molecule per second + KM_tetR = params.get('KM_tetR', 40) + KM_cI = params.get('KM_cI', 40) + KM_lacI = params.get('KM_lacI', 40) + + # the differential equations governing the state variables: + # mRNA concentration = transcription given repressor concentration - + # mRNA degradation + transcription leakage + dxdt = np.empty(6) + dxdt[0] = k_transcription_cI/(1 + (protein_tetR / KM_tetR) ** n_tetR) - \ + k_mRNA_degradation * mRNA_cI + k_transcription_leakage + dxdt[1] = k_transcription_lacI/(1 + (protein_cI / KM_cI)**n_cI) - \ + k_mRNA_degradation * mRNA_lacI + k_transcription_leakage + dxdt[2] = k_transcription_tetR/(1 + (protein_lacI / KM_lacI) ** n_lacI) - \ + k_mRNA_degradation * mRNA_tetR + k_transcription_leakage + + # protein concentration = translation - protein degradation + dxdt[3] = k_translation*mRNA_cI - k_protein_degradation*protein_cI + dxdt[4] = k_translation*mRNA_lacI - k_protein_degradation*protein_lacI + dxdt[5] = k_translation*mRNA_tetR - k_protein_degradation*protein_tetR + + return dxdt + +# Define the system as an I/O system +sys = ct.nlsys( + updfcn=repressilator, outfcn=lambda t, x, u, params: x[3:], + states=6, inputs=0, outputs=3) + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) + +# +# (b) A simulation of a simple model for the repressilator, showing the +# oscillation of the individual protein concentrations. +# + +fig.add_subplot(gs[0, 1]) # first row, second column + +# Initial conditions and time +t = np.linspace(0, 20000, 1000) +x0 = [1, 0, 0, 200, 0, 0] + +# Integrate the differential equation +response = ct.input_output_response(sys, t, 0, x0) + +# Plot the results (protein concentrations) +plt.plot(response.time/60, response.outputs[0], '-') +plt.plot(response.time/60, response.outputs[1], '--') +plt.plot(response.time/60, response.outputs[2], '-.') + +plt.axis([0, 300, 0, 5000]) +plt.legend(("cI", "lacI", "tetR"), loc='upper right') + +plt.xlabel("Time [min]") # Axis labels +plt.ylabel("Proteins per cell") +plt.title("Repressilator simulation") # Plot title + +# Save the figure +plt.savefig("figure-3.26-repressilator_dynamics.png", bbox_inches='tight') diff --git a/examples/example-3.19-fitzhugh_nagumo.py b/examples/example-3.19-fitzhugh_nagumo.py new file mode 100644 index 000000000..a674d78e3 --- /dev/null +++ b/examples/example-3.19-fitzhugh_nagumo.py @@ -0,0 +1,75 @@ +# example-3.19-fitzhugh_nagumo.py - nerve cell dynamics +# RMM, 30 Aug 2021 +# +# Figure 3.28: Response of a neuron to a current input. The current +# input is shown in (a) and the neuron voltage V in (b). The +# simulation was done using the FitzHugh–Nagumo model (Exercise 3.11). + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +ct.use_fbs_defaults() + +# +# System dynamics +# + +# FitzHugh-Nagumo dynamics (from KJA) +def fitzhugh_nagumo_dynamics(t, x, u, params): + dx = np.zeros(3) + + # Get the system state + V = x[0] + R = x[1] + + # Compute the dim derivative + dx[0] = 10 * (V - (V**3) / 3 - R + u[0]) + dx[1] = 0.8 * (-R + 1.25 * V + 1.5) + dx[2] = 1 + + return dx + +# Set up an input/output system +sys = ct.nlsys( + updfcn=fitzhugh_nagumo_dynamics, states=3, inputs=1, outputs=3) + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(3, 2) + +# +# (a) Current input +# + +fig.add_subplot(gs[0, 0]) # first row, first column + +# Set up input stimulation +t = np.linspace(0, 50, 500) +u = np.zeros_like(t) +u[t >= 5] = 1.5 # start of short input pulse +u[t >= 6] = 0 # end of short input pulse +u[t >= 30] = 1.5 # longer input pulse + +# Initial state +x0 = [-1.5, -3/8, 0] + +response = ct.input_output_response(sys, t, u, x0) + +plt.plot(response.time, response.inputs[0]) +plt.xlabel("Time $t$ [ms]") +plt.ylabel("Current $I$ [mA]") +plt.title("Input stimulation") + +# +# (b) Neuron response +# + +fig.add_subplot(gs[0, 1]) # first row, second column + +plt.plot(response.time, response.states[0]) +plt.xlabel("Time $t$ [ms]") +plt.ylabel("Voltage $V$ [mV]") +plt.title("Neuron response") + +# Save the figure +plt.savefig("figure-3.19-fitzhugh_nagumo.png", bbox_inches='tight') diff --git a/examples/example-3.4-predator_prey.py b/examples/example-3.4-predator_prey.py new file mode 100644 index 000000000..9bc004154 --- /dev/null +++ b/examples/example-3.4-predator_prey.py @@ -0,0 +1,61 @@ +# example-3.4-predator_prey.py - discrete-time simulation of predator–prey model +# RMM, 15 May 2019 +# +# Figure 3.8: Discrete-time simulation of the predator–prey model +# (3.13). Using the parameters a = c = 0.014, bh(u) = 0.6, and dl = +# 0.7 in equation (3.13), the period and magnitude of the lynx and +# hare population cycles approximately match the data in Figure 3.7. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt + +# Define the dynamics of the predator prey system +def predprey(t, x, u, params): + # Parameter setup + a = params.get('a', 0.014) / 365. + bh0 = params.get('bh0', 0.6) / 365. + c = params.get('c', 0.014) / 365. + dl = params.get('dl', 0.7) / 365. + + # Map the states into local variable names + H = x[0] + L = x[1] + + # Compute the input + bhu = bh0 + u[0] + + # Compute the discrete updates + dH = H + bhu * H - a * L * H + dL = L + c * L * H - dl * L + + return [dH, dL] + +# Create a nonlinear I/O system (dt = days) +io_predprey = ct.nlsys( + predprey, None, inputs=['u'], outputs=['H', 'L'], + states=['H', 'L'], name='predprey', dt=1/365) + +X0 = [10, 10] # Initial H, L +T = np.linspace(1845, 1935, 90*365 + 1) # 90 years + +# Simulate the system +response = ct.input_output_response(io_predprey, T, 0, X0) +t, y = response.time, response.outputs + +# Downsample the responses to individual years +yrs = t[::365] +pop = y[:, ::365] + +# Plot the response +plt.subplot(2, 1, 1) # Set the aspect ration to match the text + +plt.plot(yrs, pop[0], 'b-o', markersize=3) +plt.plot(yrs, pop[1], 'r--o', markersize=3) +plt.legend(["Hare", "Lynx"]) + +plt.xlabel("Year") +plt.ylabel("Population") + +# Save the figure +plt.savefig("figure-3.8-predator_prey.png", bbox_inches='tight') diff --git a/examples/example-5.13-genetic_switch.py b/examples/example-5.13-genetic_switch.py new file mode 100644 index 000000000..ff97b652f --- /dev/null +++ b/examples/example-5.13-genetic_switch.py @@ -0,0 +1,139 @@ +# example-5.13-genetic_switch.py - Genetic switch dynamics +# RMM, 17 Nov 2024 +# +# Figure 5.15: Stability of a genetic switch. The circuit diagram in (a) +# represents two proteins that are each repressing the production of the +# other. The inputs u1 and u2 interfere with this repression, allowing the +# circuit dynamics to be modified. The equilibrium points for this circuit +# can be determined by the intersection of the two curves shown in (b). +# +# Figure 5.16: Dynamics of a genetic switch. The phase portrait on the left +# shows that the switch has three equilibrium points, corresponding to +# protein A having a concentration greater than, equal to, or less than +# protein B. The equilibrium point with equal protein concentrations is +# unstable, but the other equilibrium points are stable. The simulation on +# the right shows the time response of the system starting from two +# different initial conditions. The initial portion of the curve +# corresponds to initial concentrations z(0) = (1, 5) and converges to the +# equilibrium point where z1e < z2e. At time t = 10, the concentrations are +# perturbed by +2 in z1 and −2 in z2, moving the state into the region of +# the state space whose solutions converge to the equilibrium point where +# z2e < z1e. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp +from scipy.optimize import fsolve +ct.use_fbs_defaults() + +# +# System dynamics +# + +# Switch parameters +genswitch_params = {'mu': 4, 'n': 2} + +def genswitch_dynamics(t, x, u, params): + mu = params.get('mu') + n = params.get('n') + z1, z2 = x + dz1 = mu / (1 + z2**n) - z1 + dz2 = mu / (1 + z1**n) - z2 + return [dz1, dz2] + +genswitch_model = ct.nlsys( + genswitch_dynamics, None, states=2, inputs=None, params=genswitch_params) + +# +# 5.15 (b) Equilibrium points +# + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) +ax = fig.add_subplot(gs[0, 1]) + +# Generate nullcline +mu, n = genswitch_params['mu'], genswitch_params['n'] +u = np.linspace(0, 5, 100) +f = mu / (1 + u**n) + +# Find equilibrium points +def eq_func(z): + return mu / (1 + z**2) - z + +eqpts = [fsolve(eq_func, 2)[0], fsolve(eq_func, 0)[0], fsolve(eq_func, -2)[0]] + +ax.plot(u, f, 'b-', label='$z_1, f(z_1)$') +ax.plot(f, u, 'r--', label='$z_2, f(z_2)$') +ax.plot([0, 3], [0, 3], 'k-', linewidth=0.5) +ax.scatter(eqpts, eqpts[::-1], color='k') +ax.axis([0, 5, 0, 5]) +ax.set_xlabel('$z_1, f(z_2)$') +ax.set_ylabel('$z_2, f(z_1)$') +ax.legend() +ax.set_title('(b) Equilibrium Curves') + +# Save the first figure +plt.savefig("figure-5.15-genswitch_nullclines.png", bbox_inches='tight') + +# +# 5.16 (a) Phase portrait +# + +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) +ax = fig.add_subplot(gs[0, 0]) # first row, first column + +ct.phase_plane_plot( + genswitch_model, [0, 5, 0, 5], 10, gridspec=[7, 7], ax=ax, + plot_separatrices=False) + +ax.axis('scaled') +ax.axis([0, 5, 0, 5]) +ax.set_xticks([0, 1, 2, 3, 4, 5]) # defaults are different than y axis (?) +ax.set_xlabel("Protein A [scaled]") +ax.set_ylabel("Protein B [scaled]") +ax.set_title('(a) Phase portrait') + +# +# 5.16 (b) Time traces +# + +ax = fig.add_subplot(gs[0, 1]) # first row, second column + +# Solve the ODE for the first segment +sol1 = ct.input_output_response( + genswitch_model, np.linspace(0, 10, 1000), 0, [1, 5]) + +# Solve the ODE for the second segment +sol2 = ct.input_output_response( + genswitch_model, np.linspace(11, 25, 1000), 0, + [sol1.states[0, -1] + 2, sol1.states[1, -1] - 2]) + +# Second plot: Time traces +ax.plot(sol1.time, sol1.outputs[0], 'b-', label='$z_1$') +ax.plot(sol1.time, sol1.outputs[1], 'r--', label='$z_2$') +ax.plot(sol2.time, sol2.outputs[0], 'b-') +ax.plot(sol2.time, sol2.outputs[1], 'r--') +ax.plot( + [sol1.time[-1], sol2.time[0]], [sol1.outputs[0, -1], sol2.outputs[0, -0]], + 'k:') +ax.plot( + [sol1.time[-1], sol2.time[0]], [sol1.outputs[1, -1], sol2.outputs[1, -0]], + 'k:') +ax.scatter( + [sol1.time[-1], sol2.time[0]], [sol1.outputs[0, -1], sol2.outputs[0, 0]], + edgecolors='b', facecolors='none', marker='o') +ax.scatter( + [sol1.time[-1], sol2.time[0]], [sol1.outputs[1, -1], sol2.outputs[1, 0]], + edgecolors='r', facecolors='none', marker='o') +ax.axis([0, 25, 0, 5]) +ax.set_xlabel('Time [scaled]') +ax.set_ylabel('Protein concentrations [scaled]') +ax.legend() +ax.set_title('(b) Simulation time traces') + +# Save the second figure +plt.savefig("figure-5.16-genswitch_dynamics.png", bbox_inches='tight') diff --git a/examples/example-5.14-invpend_stabilized.py b/examples/example-5.14-invpend_stabilized.py new file mode 100644 index 000000000..2ebcf1b11 --- /dev/null +++ b/examples/example-5.14-invpend_stabilized.py @@ -0,0 +1,96 @@ +# example-5.14-invpend_stabilized.py - stabilized inverted pendulum +# RMM, 17 Nov 2024 + +# Figure 5.17: Stabilized inverted pendulum. A control law applies a +# force u at the bottom of the pendulum to stabilize the inverted +# position (a). The phase portrait (b) shows that the equilibrium +# point corresponding to the vertical position is stabilized. The +# shaded region indicates the set of initial conditions that converge +# to the origin. The ellipse corresponds to a level set of a Lyapunov +# function V (x) for which V (x) > 0 and V ̇ (x) < 0 for all points +# inside the ellipse. This can be used as an estimate of the region of +# attraction of the equilibrium point. The actual dynamics of the +# system evolve on a manifold (c). + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +from math import pi, sin, cos, sqrt +ct.use_fbs_defaults() + +# +# System dynamics +# + +# Stablized (and normalized) inverted pendulum dynamics +def balpend_update(t, x, v, params): + a = params.get('a', 2) + u = -2 * a * sin(x[0]) - x[1] * cos(x[0]) + return [x[1], sin(x[0]) + u * cos(x[0])] +balpend = ct.nlsys( + balpend_update, states=2, inputs=0, name='inverted pendulum') + + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 4) +ax = fig.add_subplot(gs[0, 1:3]) + +# Generate the phase plot (easy part) +cplt = ct.phase_plane_plot( + balpend, [-2*pi, 2*pi, -4, 4], 4.5, gridspec=[10, 10], + plot_separatrices={'timedata': np.linspace(0, 20, 500), 'arrows': 1}, + ax=ax) + +# +# Next we need to shade the region of attraction containing the origin. We +# do this by extracting out the bounding separatrices and then filling in +# the region between them, using some messy NumPy/Matplotlib commands... +# + +# Get the separatrices on either side of the origin +for line in cplt.lines[0]: + xdata, ydata = line.get_data() + if line.get_linestyle() != '--' or abs(ydata[-1]) > 0.1: + continue # Not the line we are looking for + + # Extract out the separatrices for the middle region + if xdata[-1] < 0 and xdata[-1] > -2 and ydata[-1] > 0: + ul = (xdata[::-1], ydata[::-1]) # Sort along increasing y + elif xdata[-1] < 0 and xdata[-1] > -2 and ydata[-1] < 0: + ll = (xdata, ydata) + elif xdata[-1] > 0 and xdata[-1] < 2 and ydata[-1] > 0: + ur = (xdata[::-1], ydata[::-1]) # Sort along increasing y + elif xdata[-1] > 0 and xdata[-1] < 2 and ydata[-1] < 0: + lr = (xdata, ydata) + +# Shade the region between the separatrices (including the middle) +ax.fill_betweenx(ul[1], ul[0], np.interp(ul[1], ur[1], ur[0]), color='0.95') +ax.fill_betweenx(lr[1], np.interp(lr[1], ll[1], ll[0]), lr[0], color='0.95') +ax.fill_betweenx( + [ul[1][0], ll[1][-1]], [ul[0][0], ll[0][-1]], [ur[0][0], lr[0][-1]], + color='0.95') + +# Add the Lypaunov level set +A = balpend.linearize(0, 0).A # Linearized dynamics matrix +P = ct.lyap(A.T, np.eye(2)) # Solve Lyapunov equation for P + +# Figure out the states along the level set +rho = 1.2 # value of the level set +xval, yval = [], [] +for theta in np.linspace(0, 2*pi, 100): + # Find the length of state vector at this angle that gives V(x) = 1 + r = 1 / sqrt(np.array([cos(theta), sin(theta)]).T @ P @ + np.array([cos(theta), sin(theta)])) + xval.append(rho * r * cos(theta)) + yval.append(rho * r * sin(theta)) +ax.plot(xval, yval, 'r-', linewidth=2) + +# Label the plot +ax.set_xticks([-2*pi, -pi, 0, pi, 2*pi]) +ax.set_xticklabels([r'$-2\pi$', r'$-pi$', r'$0$', r'$\pi$', r'$2\pi$']) +ax.set_xlabel('$x_1$') +ax.set_ylabel('$x_2$', rotation=0) +ax.set_title("(b) Phase portrait") + +plt.savefig('figure-5.17-balpend_phaseplot.png', bbox_inches='tight') diff --git a/examples/example-5.16-predprey_bif.py b/examples/example-5.16-predprey_bif.py new file mode 100644 index 000000000..6a18a6b60 --- /dev/null +++ b/examples/example-5.16-predprey_bif.py @@ -0,0 +1,144 @@ +# example-5.16-predprey_bif.py - predator-prey stability analysis +# RMM, 18 Nov 2024 +# +# Figure 5.18: Bifurcation analysis of the predator–prey system. (a) +# Parametric stability diagram showing the regions in parameter space for +# which the system is stable. (b) Bifurcation diagram showing the location +# and stability of the equilib- rium point as a function of a. The solid +# line represents a stable equilibrium point, and the dashed line +# represents an unstable equilibrium point. The dash-dotted lines indicate +# the upper and lower bounds for the limit cycle at that parameter value +# (computed via simulation). The nominal values of the parameters in the +# model are a = 3.2, b = 0.6, c = 50, d = 0.56, k = 125, and r = 1.6. + +import control as ct +import numpy as np +import scipy +import matplotlib.pyplot as plt +ct.use_fbs_defaults() + +# +# System dynamics +# + +from predprey import predprey, predprey_params + +# Create a function to compute the real part of the largest eigenvalue +def maxeig(a, c): + # Initialize parameter values + r, d, b, k = map(predprey_params.get, ['r', 'd', 'b', 'k']) + params = {'a': a, 'c': c} + + # Equilibrium point from equations (4.33) and (4.34) + xeq = [(c*d) / (a*b - d), (b*c*r)*(a*b*k - c*d - d*k)/(k * (a*b - d)**2)] + + # Linearization + A = ct.linearize(predprey, xeq, 0, params=params).A + + return np.max(np.linalg.eig(A).eigenvalues.real) + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) + +# +# (a) Stability diagram +# +# To find the boundaries of stability in terms of the parameters $a$ and +# $c$, we scan over the $a$ parameter and find the values of $c$ that cause +# the real part of one of the eigenvalues to be zero. +# + +ax = fig.add_subplot(gs[0, 0]) # first row, first column +ax.set_title("(a) Stability diagram") + +avals = np.linspace(1.3, 4, 50) +cvals_lower, cvals_upper = [], [] +last_lower, last_upper = 5, 100 +for a in avals: + sol1 = scipy.optimize.root( + lambda c: maxeig(a, c), last_lower, method='broyden1') + if sol1.success: + last_lower = sol1.x.item() + cvals_lower.append(last_lower) + else: + print(sol1.message) + cvals_lower.append(np.nan) + + sol2 = scipy.optimize.root(lambda c: maxeig(a, c), last_upper) + if sol2.success: + last_upper = sol2.x.item() + cvals_upper.append(last_upper) + else: + cvals_upper.append(np.nan) + +ax.plot(avals, cvals_lower, 'k', linewidth=0.5) +ax.plot(avals, cvals_upper, 'k', linewidth=0.5) +ax.fill_between(avals, cvals_lower, cvals_upper, color='0.9') + +ax.set_xlabel("$a$") +ax.set_ylabel("$c$", rotation=0) +ax.text(1.4, 160, "Unstable") +ax.text(2.2, 100, "Stable") +ax.text(3, 25, "Unstable") +ax.axis('tight') +ax.axis([1.35, 4, 0, 200]) + +# +# (b) Bifurcation diagram +# + +ax = fig.add_subplot(gs[0, 1]) # first row, second column +ax.set_title("(b) Bifurcation diagram") + +# Create lists to hold the values of the different branches on the plot +stable_H, unstable_H = [], [] # Equilibrium point +lower_H, upper_H = [], [] # Limit cycle bounds + +# Set the values of 'a' to be denser near the bifurcation point +avals = np.hstack( + [np.linspace(1.35, 2, 10), np.linspace(2, 4, 100), np.linspace(4, 8, 20)]) + +# Set up the remaining parameters for the simulation +timepts = np.linspace(0, 300, 5000) +params = predprey_params + +# Compute the branches of the bifurcation diagram +for a in avals: + # Set the parameter values + params['a'] = a + + # Equilibrium point from equations (4.33) and (4.34) + r, d, b, k, c = map(predprey_params.get, ['r', 'd', 'b', 'k', 'c']) + xeq = [(c*d) / (a*b - d), (b*c*r)*(a*b*k - c*d - d*k)/(k * (a*b - d)**2)] + + # Check stability + if maxeig(a, params['c']) < 0: + # Stable branch + stable_H.append(xeq[0]) + [vlist.append(np.nan) for vlist in [unstable_H, lower_H, upper_H]] + else: + # Unstable branch + unstable_H.append(xeq[0]) + stable_H.append(np.nan) + + # Run a simulation to figure out size of the limit cycle + resp = ct.input_output_response( + predprey, timepts, X0=np.array(xeq) + 0.1, params=params) + lower_H.append(np.min(resp.outputs[0, -500:])) + upper_H.append(np.max(resp.outputs[0, -500:])) + +# Plot the different branches +ax.plot(avals, stable_H, 'b-') +ax.plot(avals, unstable_H, 'r--') +ax.plot(avals, lower_H, 'k-.') +ax.plot(avals, upper_H, 'k-.') + +# Label the plot +ax.set_xlabel("$a$") +ax.set_ylabel("$H$", rotation=0) +ax.axis('tight') +ax.axis([1.35, 8, 0, 150]) + +# Save the figure +plt.savefig("figure-5.18-predprey_bif.png", bbox_inches='tight') diff --git a/examples/example-5.17-bicycle_stability.py b/examples/example-5.17-bicycle_stability.py new file mode 100644 index 000000000..540d39dd2 --- /dev/null +++ b/examples/example-5.17-bicycle_stability.py @@ -0,0 +1,153 @@ +# example-5.17-bicycle_stability.py - Root locus diagram for a bicycle model +# RMM, 24 Nov 2024 +# +# Figure 5.19: Stability plots for a bicycle moving at constant +# velocity. The plot in (a) shows the real part of the system eigenvalues +# as a function of the bicycle velocity v0. The system is stable when all +# eigenvalues have negative real part (shaded region). The plot in (b) +# shows the locus of eigenvalues on the complex plane as the velocity v is +# varied and gives a different view of the stability of the system. This +# type of plot is called a root locus diagram. +# +# Notes: +# +# 1. The line styles used in this plot are slightly different than in the +# book. Solid lines are used for real-valued eigenvalues and dashed +# lines are used for the real part of complex-valued eigenvalues. +# +# 2. This code relies on features on python-control-0.10.2, which is +# currently under development. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +from math import isclose +ct.use_fbs_defaults() + +# +# System dynamics +# + +from bicycle import whipple_A + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) + +# +# (a) Stability diagram +# + +ax = fig.add_subplot(gs[0, 0]) # first row, first column +ax.set_title("(a) Stability diagram") + +# Compute the eigenvalues as a function of velocity +v0_vals = np.linspace(-15, 15, 500) +eig_vals = [] +for v0 in v0_vals: + A = whipple_A(v0) + eig_vals.append(np.sort(np.linalg.eig(A).eigenvalues)) + +# Initialize lists to categorize eigenvalues +eigs_real_stable = [] +eigs_complex_stable = [] +eigs_real_unstable = [] +eigs_complex_unstable = [] + +# Keep track of region in which all eigenvalues are stable +stable_beg = stable_end = None + +# Process each set of eigenvalues +for i, eig_set in enumerate(eig_vals): + # Create arrays filled with NaN for each category + real_stable = np.full(eig_set.shape, np.nan) + complex_stable = np.full(eig_set.shape, np.nan) + real_unstable = np.full(eig_set.shape, np.nan) + complex_unstable = np.full(eig_set.shape, np.nan) + + # Classify eigenvalues + for j, eig in enumerate(eig_set): + if isclose(eig.imag, 0): # Real eigenvalue + if eig.real < 0: + real_stable[j] = eig.real + else: + real_unstable[j] = eig.real + else: # Complex eigenvalue + if eig.real < 0: + complex_stable[j] = eig.real + else: + complex_unstable[j] = eig.real + + # Append categorized arrays to respective lists + eigs_real_stable.append(real_stable) + eigs_complex_stable.append(complex_stable) + eigs_real_unstable.append(real_unstable) + eigs_complex_unstable.append(complex_unstable) + + # Look for regions where everything is stable + if stable_beg is None and all(eig_set.real < 0): + stable_beg = i + elif stable_beg and stable_end is None and any(eig_set.real > 0): + stable_end = i + +# Plot the stability diagram +ax.plot(v0_vals, eigs_real_stable, 'b-') +ax.plot(v0_vals, eigs_real_unstable, 'r-') +ax.plot(v0_vals, eigs_complex_stable, 'b--') +ax.plot(v0_vals, eigs_complex_unstable, 'r--') + +# Add in the coordinate axes +ax.axhline(color='k', linewidth=0.5) +ax.axvline(color='k', linewidth=0.5) + +# Label and shade stable and unstable regions +ax.text(-12, 8, "Unstable") +ax.fill_betweenx( + [-15, 15], [v0_vals[stable_beg], v0_vals[stable_beg]], + [v0_vals[stable_end], v0_vals[stable_end]], color='0.9') +ax.text(7.2, 6, "Stable", rotation=90) +ax.text(11.7, 5, "Unstable", rotation=90) + +# Label the axes +ax.set_xlabel(r"Velocity $v_0$ [m/s]") +ax.set_ylabel(r"$\text{Re}\,\lambda$") +ax.axis('scaled') +ax.axis([-15, 15, -15, 15]) + +# +# (b) Root locus diagram +# + +ax = fig.add_subplot(gs[0, 1]) # first row, second column +ax.set_title("(b) Root locus diagram") + +# Generate the root locus diagram via the root_locus_plot functionality +pos_idx = np.argmax(v0_vals >= 0) +poles = eig_vals[pos_idx] +loci = np.array(eig_vals[pos_idx:]) +rl_map = ct.PoleZeroData(poles, [], v0_vals[pos_idx:], loci) +rl_map.plot(ax=ax) + +# Add in the coordinate axes +ax.axhline(color='k', linewidth=0.5) +ax.axvline(color='k', linewidth=0.5) + +# Label the real axes of the plot +ax.text(-12.5, -2, r"$\leftarrow v_0$") +ax.text(-3.5, -2, r"$v_0 \rightarrow$") + +# Label the crossover points +xo_idx = np.argmax(rl_map.loci[:, 3].real < 0) +ax.plot(0, rl_map.loci[xo_idx, 2].imag, 'bo', markersize=3) +ax.plot(0, rl_map.loci[xo_idx, 3].imag, 'bo', markersize=3) +ax.text(1, rl_map.loci[xo_idx, 2].imag, r"$v_0 = 6.1$") +ax.text(1, rl_map.loci[xo_idx, 3].imag, r"$v_0 = 6.1$") + +# Label the axes +ax.set_xlabel(r"$\text{Re}\,\lambda$") +ax.set_ylabel(r"$\text{Im}\,\lambda$") +ax.set_box_aspect(1) +ax.axis([-15, 15, -10, 10]) + +# Save the figure +plt.savefig("figure-5.19-bicycle_stability.png", bbox_inches='tight') diff --git a/examples/example-5.18-noise_cancel.py b/examples/example-5.18-noise_cancel.py new file mode 100644 index 000000000..2f25fd6cd --- /dev/null +++ b/examples/example-5.18-noise_cancel.py @@ -0,0 +1,100 @@ +# example-5.18-noise_cancel.py - Noise cancellation +# RMM, 24 Nov 2024 +# +# Figure 5.21: Simulation of noise cancellation. The upper left figure +# shows the headphone signal without noise cancellation, and the lower left +# figure shows the signal with noise cancellation. The right figures show +# the parameters a and b of the filter. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +from math import pi +ct.use_fbs_defaults() + +# +# System dynamics +# + +# Headphone dynamics +headphone_params = {'a0': -0.75, 'b0': 0.9} +def headphone_update(t, z, n, params): + return params['a0'] * z[0] + params['b0'] * n[0] +headphone = ct.nlsys( + headphone_update, inputs='n', states='z', params=headphone_params, + name='headphone') + +# Filter dynamics +def filter_update(t, w, u, params): + n, a, b = u + return a * w + b * n +filter = ct.nlsys( + filter_update, inputs=['n', 'a', 'b'], states='w', name='filter') + +# Controller dynamics +control_params = {'alpha': 1} +def control_update(t, x, u, params): + n, e, w = u + a, b = x + return [ + params['alpha'] * w * e, + params['alpha'] * n * e + ] +control = ct.nlsys( + control_update, inputs=['n', 'e', 'w'], states=['a', 'b'], name='control', + params=control_params) + +# Create summing junction to add all of the signal together +summer = ct.summing_junction(inputs=['z', 'S', '-w'], outputs='e') + +# Interconnected system +sys = ct.interconnect( + [headphone, filter, control, summer], name='noise_cancel', + inputs=['S', 'n'], outputs=['e', 'a', 'b']) + +# Create the signal and noise +timepts = np.linspace(0, 200, 2000) +signal = np.sin(0.1 * 2 * pi * timepts) # sinewave with frequency 0.1 Hz +noise = ct.white_noise(timepts, 5) # white noise with covariance 5 + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(3, 2) + +# No noise cancellation +resp_off = ct.input_output_response( + sys, timepts, [signal, noise], params={'alpha': 0}) + +ax = fig.add_subplot(gs[0, 0]) +ax.plot(resp_off.time, resp_off.outputs[0]) +ax.axis('tight') +ax.axis([0, 200, -5, 5]) +ax.set_ylabel("No cancellation") + +resp_on = ct.input_output_response( + sys, timepts, [signal, noise], params={'alpha': 1e-2}) + +ax = fig.add_subplot(gs[1, 0]) +ax.plot(resp_on.time, resp_on.outputs[0], label='e') +# ax.plot(resp_on.time, signal, label='S') +ax.axis('tight') +ax.axis([0, 200, -5, 5]) +ax.set_ylabel("Cancellation") +ax.set_xlabel("Time $t$ [s]") +# ax.legend() + +ax = fig.add_subplot(gs[0, 1]) +ax.plot(resp_on.time, resp_on.outputs[1]) +ax.axis('tight') +ax.axis([0, 200, -1.1, 0]) +ax.set_ylabel("$a$", rotation=0) + +ax = fig.add_subplot(gs[1, 1]) +ax.plot(resp_on.time, resp_on.outputs[2]) +ax.axis('tight') +ax.axis([0, 200, 0, 1.1]) +ax.set_ylabel("$b$", rotation=0) +ax.set_xlabel("Time $t$ [s]") + +# Save the figure +plt.savefig("figure-5.21-noise_cancel.png", bbox_inches='tight') diff --git a/examples/example-5.8-tanker_stabilty.py b/examples/example-5.8-tanker_stabilty.py new file mode 100644 index 000000000..1990b72d4 --- /dev/null +++ b/examples/example-5.8-tanker_stabilty.py @@ -0,0 +1,106 @@ +# example-5.8-tanker_stability.py - Stability of a tanker +# RMM, 16 Nov 2024 +# +# Figure 3.15: Stability analysis for a tanker. The rudder characteristics +# are shown in (a), where the equilibrium points are marked by circles, and +# the tanker trajec- tories are shown in (b). + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +ct.use_fbs_defaults() + +# +# System dynamics +# + +def tanker_dynamics(t, x, u, params): + # Parameter values + a1 = params.get('a1', -0.6) + a2 = params.get('a2', -0.3) + a3 = params.get('a3', -5) + a4 = params.get('a4', -2) + alpha = params.get('alpha', -2) + b1 = params.get('b1', 0.1) + b2 = params.get('b2', -0.8) + + v, r = x[0], x[1] # velocity and turning rate + delta = u[0] # rudder angle + + return [ + a1 * v + a2 * r + alpha * v * abs(v) + b1 * delta, + a3 * v + a4 * r + b2 * delta + ] +tanker_model = ct.nlsys( + tanker_dynamics, None, inputs='delta', states=['v', 'r']) + + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) + +# +# (a) Rudder curve +# +ax = fig.add_subplot(gs[0, 0]) # first row, first column + +# Compute the input for each turning rate +rvec = np.linspace(-0.4, 0.4, 50) +delta_list = [] +for r in rvec: + # Solve for the different equilibrium solutions + eqpt = ct.find_operating_point( + tanker_model, [0, 0], 0, + y0=[0, r], iy=[1], # Look for the desired turning rate + ) + delta_list.append(eqpt.inputs[0]) +dvec = np.array(delta_list) +ax.plot(dvec, rvec) + +# Add the equilibrium points at zero +for x0 in [[0.1, 0.1], [0, 0], [-0.1, -0.1]]: + eqpt = ct.find_operating_point(tanker_model, x0, 0) + ax.scatter( + eqpt.inputs[0], eqpt.outputs[1], facecolors='none', edgecolors='b') + +# Add labels and axis lines +ax.set_title("(a) Rudder curve") +ax.set_xlabel(r"Rudder angle $\delta$") +ax.set_ylabel(r"Noramlized turning rate $r$") +ax.plot([-0.1, 0.1], [0, 0], 'k', linewidth=0.5) +ax.plot([0, 0], [-0.4, 0.4], 'k', linewidth=0.5) +ax.axis([-0.1, 0.1, -0.4, 0.4]) + +# +# (b) Tanker trajectories +# +from math import sin, cos +ax = fig.add_subplot(gs[0, 1]) # first row, second column + +# Create a full tanker model, including position and orientation +def full_tanker_dynamics(t, x, u, params): + vdot, rdot = tanker_dynamics(t, x[3:], u, params) + theta, v, r = x[2], x[3], x[4] + return [ + cos(theta) + v * sin(theta), -sin(theta) + v * cos(theta), + r, vdot, rdot] +full_tanker_model = ct.nlsys( + full_tanker_dynamics, None, inputs='delta', + states=['x', 'y', 'theta', 'v', 'r']) + +# Create simulations and plot them +timepts = np.linspace(0, 100, 100) +for r0, linestyle in zip([0.1, 0, -0.1], ['b-', 'b--', 'b-']): + response = ct.input_output_response( + full_tanker_model, timepts, 0, [0, 0, 0, 0, r0]) + ax.plot(response.outputs[0], response.outputs[1], linestyle) + +# Add labels and axis lines +ax.set_title("(b) Tanker trajectories") +ax.set_xlabel("$x$") +ax.set_ylabel("$y$") +ax.axis('scaled') +ax.axis([0, 40, -20, 20]) + +# Save the figure +plt.savefig("figure-5.12-tanker_stability.png", bbox_inches='tight') diff --git a/examples/example-5.9-limit_cycle_stability.py b/examples/example-5.9-limit_cycle_stability.py new file mode 100644 index 000000000..45d50b1e3 --- /dev/null +++ b/examples/example-5.9-limit_cycle_stability.py @@ -0,0 +1,79 @@ +# example-5.9-limit_cycle_stability.py - Solution curves for stable limit cycle +# RMM, 16 Nov 2024 (with initial converstion from MATLAB by ChatGPT) +# +# Figure 5.13: Solution curves for a stable limit cycle. The phase portrait +# on the left shows that the trajectory for the system rapidly converges to +# the stable limit cycle. The starting points for the trajectories are +# marked by circles in the phase portrait. The time domain plots on the +# right show that the states do not converge to the solution but instead +# maintain a constant phase error. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp +ct.use_fbs_defaults() + +# +# System dynamics +# + +def limcyc(t, x): + E = x[0]**2 + x[1]**2 + xdot1 = x[1] + x[0] * (1 - E) + xdot2 = -x[0] + x[1] * (1 - E) + return [xdot1, xdot2] + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(8, 10) + +# Simulation settings +ts = 20 +tspan = (0, ts) +x01 = [0, 2] +x02 = [1/np.sqrt(2), 1/np.sqrt(2)] + +# +# (a) Phase plane plot +# + +ax = fig.add_subplot(gs[1:5, :4]) # left plot + +# Solve the differential equations +sol1 = solve_ivp(limcyc, tspan, x01, t_eval=np.linspace(*tspan, 1000)) +sol2 = solve_ivp(limcyc, tspan, x02, t_eval=np.linspace(*tspan, 1000)) + +# Plot phase plane +ax.plot(sol1.y[0], sol1.y[1], 'r--') +ax.plot(sol2.y[0], sol2.y[1], 'b-') +ax.scatter(0, 2, facecolors='none', edgecolors='r') +ax.scatter(1/np.sqrt(2), 1/np.sqrt(2), facecolors='none', edgecolors='b') + +ax.axhline(0, color='k') +ax.axvline(0, color='k') +ax.axis('square') +ax.axis([-1.2, 2.2, -1.2, 2.2]) +ax.set_xlabel('$x_1$') +ax.set_ylabel('$x_2$', rotation=0) + +# +# (b) Time traces +# + +ax = fig.add_subplot(gs[1:3, 5:]) # upper right plot +ax.plot(sol1.t, sol1.y[0], 'r--') +ax.plot(sol2.t, sol2.y[0], 'b-') +ax.set_xlim([0, 20]) +ax.set_ylim([-1.2, 2]) +ax.set_ylabel('$x_1$', rotation=0) + +ax = fig.add_subplot(gs[3:5, 5:]) # upper right plot +ax.plot(sol1.t, sol1.y[1], 'r--') +ax.plot(sol2.t, sol2.y[1], 'b-') +ax.set_xlim([0, 20]) +ax.set_ylim([-1.2, 2]) +ax.set_ylabel('$x_2$', rotation=0) +ax.set_xlabel('Time $t$') + +plt.savefig('figure-5.13-limit_cycle_stability.png', bbox_inches='tight') diff --git a/examples/example-6.10-compartment_response.py b/examples/example-6.10-compartment_response.py new file mode 100644 index 000000000..1e84ac39f --- /dev/null +++ b/examples/example-6.10-compartment_response.py @@ -0,0 +1,85 @@ +# example-6.7-compartment_response.py - Compartment model response +# RMM, 24 Nov 2024 +# +# Figure 6.10: Response of a compartment model to a constant drug +# infusion. A simple diagram of the system is shown in (a). The step +# response (b) shows the rate of concentration buildup in compartment 2. In +# (c) a pulse of initial concentration is used to speed up the response. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +ct.use_fbs_defaults() + +# +# System dynamics +# + +# Parameter settings for the model +k0 = 0.1 +k1 = 0.1 +k2 = 0.5 +b0 = 1.5 + +# Compartment model definition +Ada = np.array([[-k0 - k1, k1], [k2, -k2]]) +Bda = np.array([[b0], [0]]) +Cda = np.array([[0, 1]]) +compartment = ct.ss(Ada, Bda, Cda, 0); + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(4, 3) + +# +# (a) Step input +# + +timepts = np.linspace(0, 50) +input = np.ones(timepts.size) * 0.1 +response = ct.forced_response(compartment, timepts, input) + +ax = fig.add_subplot(gs[0, 1]) +ax.set_title("(b) Step input") + +ax.plot(response.time, response.outputs) +ax.axhline(response.outputs[-1], color='k', linewidth=0.5) +ax.set_ylabel("Concentration $C_2$") +ax.axis('tight') +ax.axis([0, 50, 0, 2]) + +ax = fig.add_subplot(gs[1, 1]) +ax.plot(response.time, response.inputs) +ax.set_xlabel("Time $t$ [min]") +ax.set_ylabel("Input dosage") +ax.axis('tight') +ax.axis([0, 50, 0, 0.4]) + +# +# (b) Pulse input +# + +timepts = np.linspace(0, 50, 200) +input = np.ones(timepts.size) * 0.1 +input[:20] = 0.3 # Increase value for first 5 seconds +response = ct.forced_response(compartment, timepts, input) + +ax = fig.add_subplot(gs[0, 2]) +ax.set_title("(c) Pulse input") + +ax.plot(response.time, response.outputs) +ax.axhline(response.outputs[-1], color='k', linewidth=0.5) +ax.set_ylabel("Concentration $C_2$") +ax.axis('tight') +ax.axis([0, 50, 0, 2]) + +ax = fig.add_subplot(gs[1, 2]) +ax.plot(response.time, response.inputs) +ax.set_xlabel("Time $t$ [min]") +ax.set_ylabel("Input dosage") +ax.axis('tight') +ax.axis([0, 50, 0, 0.4]) + +# Save the figure +fig.align_ylabels() +plt.savefig("figure-6.10-compartment_response.png", bbox_inches='tight') diff --git a/examples/example-6.8-opamp_bandpass.py b/examples/example-6.8-opamp_bandpass.py new file mode 100644 index 000000000..1c6da03c9 --- /dev/null +++ b/examples/example-6.8-opamp_bandpass.py @@ -0,0 +1,53 @@ +# example-6.8-opamp_bandpass.py - +# RMM, 28 Nov 2028 +# +# Figure 6.12: Active band-pass filter. The circuit diagram (a) shows an op +# amp with two RC filters arranged to provide a band-pass filter. The plot +# in (b) shows the gain and phase of the filter as a function of frequency. +# Note that the phase starts at -90 deg due to the negative gain of the +# operational amplifier. + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +ct.use_fbs_defaults() + +# +# System dynamics +# + +# Filter parameters +R1 = 100 +R2 = 5000 +C1 = 100e-6 +C2 = 100e-6 + +# State-space form of the solution +A = np.array([[-1 / (R1 * C1), 0], + [1 / (R1 * C2), -1 / (R2 * C2)]]) +B = np.array([[1 / (R1 * C1)], + [-1 / (R1 * C2)]]) +C = np.array([[0, 1]]) +D = np.array([[0]]) + +# Create the state-space system +sys = ct.ss(A, B, C, D) + +# +# (b) Frequency response +# + +# Generate the Bode plot data +omega = np.logspace(-1, 3, 100) # Frequency range from 0.1 to 1000 rad/s +cplt = ct.frequency_response(sys, omega).plot( + initial_phase=-90, title="(b) Frequency response") + +# Plot the unit gain line +cplt.axes[0, 0].axhline(1, color='k', linestyle='-', linewidth=0.5) + +# Print system poles for bandwidth computation +print("System poles (for computing bandwidth):") +print(sys.poles()) + +# Save the figure +plt.savefig("figure-6.8-opamp_bandpass.png", bbox_inches='tight') diff --git a/examples/example-6.9-afm_freqresp.py b/examples/example-6.9-afm_freqresp.py new file mode 100644 index 000000000..253d50afb --- /dev/null +++ b/examples/example-6.9-afm_freqresp.py @@ -0,0 +1,113 @@ +# example-6.9-afm_freqresp.py - +# RMM, 28 Nov 2024 +# +# Figure 6.13: AFM frequency response. (a) A block diagram for the vertical +# dynamics of an atomic force microscope in contact mode. The plot in (b) +# shows the gain and phase for the piezo stack. The response contains two +# frequency peaks at resonances of the system, along with an antiresonance +# at ω = 268 krad/s. The combination of a resonant peak followed by an +# antiresonance is common for systems with multiple lightly damped +# modes. The dashed horizontal line represents the gain equal to the zero +# frequency gain divided by sqrt(2). + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +from math import sqrt +ct.use_fbs_defaults() + +# +# System dynamics +# + +# System parameters +m1 = 0.15e-3 +m2 = 1e-3 +f1 = 40.9e3 +f2 = 41.6e3 +f3 = 120e3 + +# Derived quantities +w1 = 2 * np.pi * f1 +w2 = 2 * np.pi * f2 +w3 = 2 * np.pi * f3 +z1 = 0.1 +z3 = 0.1 +k2 = w2**2 * m2 +c2 = m2 * z1 * w1 + +# System matrices +A = np.array([[0, 1, 0, 0], + [-k2 / (m1 + m2), -c2 / (m1 + m2), 1 / m2, 0], + [0, 0, 0, 1], + [0, 0, -w3**2, -2 * z3 * w3]]) +B = np.array([[0], [0], [0], [w3**2]]) +C = (m2 / (m1 + m2)) * np.array([m1 * k2 / (m1 + m2), m1 * c2 / (m1 + m2), 1, 0]) +D = np.array([[0]]) + +# Create the state-space system +sys = ct.ss(A, B, C, D) + +# +# (b) Frequency response +# + +# Generate the frequency response +omega = np.logspace(4, 7, 10000) # Frequency range from 10^4 to 10^7 rad/s +freqresp = ct.frequency_response(sys, omega) +mag, phase = freqresp.magnitude, freqresp.phase + +# Create the Bode plot +cplt = freqresp.plot() +mag_ax, phase_ax = cplt.axes[:, 0] +cplt.set_plot_title("(b) Frequency response") + +# Locate peaks and valleys +# Find the first peak +max_mag = 0 +omega1 = None +for i, m in enumerate(mag): + if m > max_mag: + max_mag = m + omega1 = omega[i] + elif m < max_mag: + break + +# Find the first valley +min_mag = max_mag +omega2 = None +for i, m in enumerate(mag): + if m < min_mag: + min_mag = m + omega2 = omega[i] + elif m > min_mag: + break + +# Find the second peak (must be higher than first) +omega3 = None +for i, m in enumerate(mag): + if m > max_mag: + max_mag = m + omega3 = omega[i] + elif m < max_mag and omega3 is not None: + break + +# Print peaks and valley frequencies +print(f"Peaks at {omega1:.2e}, {omega3:.2e}; valley at {omega2:.2e}") + +# Add lines to mark the frequencies +mag_ax.axhline(1 / sqrt(2), color='r', linestyle='--', linewidth=0.5) + +mag_ax.axvline(omega1, color='k', linestyle='--', linewidth=0.5) +mag_ax.text(omega1 * 1.1, 2, r"$M_{r1}$") +mag_ax.text(omega1 * 1.1, 0.07, r"$\omega = \omega_{r1}$") + +mag_ax.axvline(omega3, color='k', linestyle='--', linewidth=0.5) +mag_ax.text(omega3 * 1.2, 3, r"$M_{r2}$") +mag_ax.text(omega3 * 1.1, 0.07, r"$\omega = \omega_{r2}$") + +phase_ax.axvline(omega1, color='k', linestyle='--', linewidth=0.5) +phase_ax.axvline(omega3, color='k', linestyle='--', linewidth=0.5) + +# Save the figure +plt.savefig("figure-6.9-afm_freqresp.png", bbox_inches='tight') diff --git a/examples/example-7.4-steering_place.py b/examples/example-7.4-steering_place.py new file mode 100644 index 000000000..5f789dd26 --- /dev/null +++ b/examples/example-7.4-steering_place.py @@ -0,0 +1,122 @@ +# example-7.4-steering_place.py - +# RMM, 28 Nov 2024 +# +# Figure 7.6: State feedback control of a steering system. Unit step +# responses (from zero initial condition) obtained with controllers +# designed with zeta_c = 0.7 and omega_c = 0.5, 0.7, and 1 [rad/s] are +# shown in (a). The dashed lines indicate ±5% deviations from the +# setpoint. Notice that response speed increases with increasing omega_c, +# but that large omega_c also give large initial control actions. Unit step +# responses obtained with a controller designed with omega_c = 0.7 and +# zeta_c = 0.5, 0.7, and 1 are shown 2in (b). + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +ct.use_fbs_defaults() + +# +# System dynamics +# + +# Get the normalized linear dynamics +from steering import linearize_lateral +sys = linearize_lateral(normalize=True, output_full_state=True) + +# Function to place the poles at desired values +def steering_place(sys, omega, zeta): + # Get the pole locations based on omega and zeta + desired_poly = np.polynomial.Polynomial([omega**2, 2 * zeta * omega, 1]) + return ct.place(sys.A, sys.B, desired_poly.roots()) + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(3, 14) # allow some space in the middle +left = slice(0, 6) +right = slice(7, 13) + +# +# (a) Unit step response for varying omega_c +# + +ax_pos = fig.add_subplot(gs[0, left]) +ax_delta = fig.add_subplot(gs[1, left]) +ax_pos.set_title( + r"(a) Unit step response for varying $\omega_c$", size='medium') + +timepts = np.linspace(0, 20) +zeta_c = 0.7 +for omega_c in [0.5, 0.7, 1]: + # Compute the gains for the controller + K = steering_place(sys, omega_c, zeta_c) + kf = omega_c**2 + + # Compute the closed loop system + ctrl, clsys = ct.create_statefbk_iosystem( + sys, K, kf, feedfwd_pattern='refgain') + + # Simulate the closed loop dynamics + response = ct.forced_response(clsys, timepts, 1, X0=0) + + ax_pos.plot(response.time, response.outputs['y'], 'b') + ax_delta.plot(response.time, response.outputs['delta'], 'b') + +# Label the plot +ax_pos.set_ylabel(r"Lateral position $y/b$") +ax_delta.set_xlabel(r"Normalized time $v_0 t/b$") +ax_delta.set_ylabel(r"Steering angle $\delta$ [rad]") + +ax_pos.axhline(0.95, color='k', linestyle='--', linewidth=0.5) +ax_pos.axhline(1.05, color='k', linestyle='--', linewidth=0.5) +ax_pos.annotate( + "", xy=(3.5, 0.3), xytext=[0, 0.8], arrowprops={'arrowstyle': '<-'}) +ax_pos.text(3.6, 0.25, r"$\omega_c$") + +ax_delta.annotate( + "", xy=(4, 0.1), xytext=(0, -0.2), arrowprops={'arrowstyle': '<-'}) +ax_delta.text(3.5, 0.15, r"$\omega_c$") + +# +# (b) Unit step response for varying zeta_c +# + +ax_pos = fig.add_subplot(gs[0, right], sharey=ax_pos) +ax_delta = fig.add_subplot(gs[1, right], sharey=ax_delta) +ax_pos.set_title( + r"(b) Unit step response for varying $\zeta_c$", size='medium') + +timepts = np.linspace(0, 20) +omega_c = 0.7 +for zeta_c in [0.5, 0.7, 1]: + # Compute the gains for the controller + K = steering_place(sys, omega_c, zeta_c) + kf = omega_c**2 + + # Compute the closed loop system + ctrl, clsys = ct.create_statefbk_iosystem( + sys, K, kf, feedfwd_pattern='refgain') + + # Simulate the closed loop dynamics + response = ct.forced_response(clsys, timepts, 1, X0=0) + + ax_pos.plot(response.time, response.outputs['y'], 'b') + ax_delta.plot(response.time, response.outputs['delta'], 'b') + +# Label the plot +ax_pos.set_ylabel(r"Lateral position $y/b$") +ax_delta.set_xlabel(r"Normalized time $v_0 t/b$") +ax_delta.set_ylabel(r"Steering angle $\delta$ [rad]") + +ax_pos.axhline(0.95, color='k', linestyle='--', linewidth=0.5) +ax_pos.axhline(1.05, color='k', linestyle='--', linewidth=0.5) +ax_pos.annotate( + "", xy=(1.7, 1.15), xytext=(5.2, 0.6), arrowprops={'arrowstyle': '<-'}) +ax_pos.text(5.5, 0.5, r"$\zeta_c$") + +ax_delta.annotate( + "", xy=(5.5, -0.2), xytext=(3.5, 0.15), arrowprops={'arrowstyle': '<-'}) +ax_delta.text(3, 0.2, r"$\zeta_c$") + +# Save the figure +fig.align_ylabels() +plt.savefig("figure-7.6-steering_place.png", bbox_inches='tight') diff --git a/examples/example-8.10-steering_gainsched.py b/examples/example-8.10-steering_gainsched.py new file mode 100644 index 000000000..8c151333f --- /dev/null +++ b/examples/example-8.10-steering_gainsched.py @@ -0,0 +1,178 @@ +# example-8.10-steering_gainsched.py - gain scheduling for vehicle steering +# RMM, 8 May 2019 + +import numpy as np +import control as ct +from cmath import sqrt +import matplotlib.pyplot as plt +import fbs # FBS plotting customizations + +# +# Vehicle steering dynamics +# +# The vehicle dynamics are given by a simple bicycle model. We take the state +# of the system as (x, y, theta) where (x, y) is the position of the vehicle +# in the plane and theta is the angle of the vehicle with respect to +# horizontal. The vehicle input is given by (v, phi) where v is the forward +# velocity of the vehicle and phi is the angle of the steering wheel. The +# model includes saturation of the vehicle steering angle. +# +# System state: x, y, theta +# System input: v, phi +# System output: x, y +# System parameters: wheelbase, maxsteer +# +def vehicle_update(t, x, u, params): + # Get the parameters for the model + l = params.get('wheelbase', 3.) # vehicle wheelbase + phimax = params.get('maxsteer', 0.5) # max steering angle (rad) + + # Saturate the steering input + phi = np.clip(u[1], -phimax, phimax) + + # Return the derivative of the state + return np.array([ + np.cos(x[2]) * u[0], # xdot = cos(theta) v + np.sin(x[2]) * u[0], # ydot = sin(theta) v + (u[0] / l) * np.tan(phi) # thdot = v/l tan(phi) + ]) + +def vehicle_output(t, x, u, params): + return x # return x, y, theta (full state) + +# Define the vehicle steering dynamics as an input/output system +vehicle = ct.nlsys( + vehicle_update, vehicle_output, states=3, name='vehicle', + inputs=('v', 'phi'), + outputs=('x', 'y', 'theta')) + +# +# Gain scheduled controller +# +# For this system we use a simple schedule on the forward vehicle velocity and +# place the poles of the system at fixed values. The controller takes the +# current and desired vehicle position and orientation plus the velocity +# velocity as inputs, and returns the velocity and steering commands. +# +# System state: none +# System input: x, y, theta, xd, yd, thetad, vd, phid +# System output: v, phi +# System parameters: longpole, latomega_c, latzeta_c +# +def control_output(t, x, u, params): + # Get the controller parameters + longpole = params.get('longpole', -2.) + latomega_c = params.get('latomega_c', 2) + latzeta_c = params.get('latzeta_c', 0.5) + l = params.get('wheelbase', 3) + + # Extract the system inputs and compute the errors + x, y, theta, xd, yd, thetad, vd, phid = u + ex, ey, etheta = x - xd, y - yd, theta - thetad + + # Determine the controller gains + lambda1 = -longpole + a1 = 2 * latzeta_c * latomega_c + a2 = latomega_c**2 + + # Compute and return the control law + v = -lambda1 * ex # leave off feedforward to generate transient + if vd != 0: + phi = phid - ((a2 * l) / vd**2) * ey - ((a1 * l) / vd) * etheta + else: + # We aren't moving, so don't turn the steering wheel + phi = phid + + return np.array([v, phi]) + +# Define the controller as an input/output system +controller = ct.nlsys( + None, control_output, name='controller', # static system + inputs=('x', 'y', 'theta', 'xd', 'yd', 'thetad', # system inputs + 'vd', 'phid'), + outputs=('v', 'phi') # system outputs +) + +# +# Reference trajectory subsystem +# +# The reference trajectory block generates a simple trajectory for the system +# given the desired speed (vref) and lateral position (yref). The trajectory +# consists of a straight line of the form (vref * t, yref, 0) with nominal +# input (vref, 0). +# +# System state: none +# System input: vref, yref +# System output: xd, yd, thetad, vd, phid +# System parameters: none +# +def trajgen_output(t, x, u, params): + vref, yref = u + return np.array([vref * t, yref, 0, vref, 0]) + +# Define the trajectory generator as an input/output system +trajgen = ct.nlsys( + None, trajgen_output, name='trajgen', + inputs=('vref', 'yref'), + outputs=('xd', 'yd', 'thetad', 'vd', 'phid')) + +# +# System construction +# +# The input to the full closed loop system is the desired lateral position and +# the desired forward velocity. The output for the system is taken as the +# full vehicle state plus the velocity of the vehicle. +# +# We construct the system using the interconnect function and using signal +# labels to keep track of everything. + +steering = ct.interconnect( + # List of subsystems + (trajgen, controller, vehicle), name='steering', + + # System inputs + inplist=['trajgen.vref', 'trajgen.yref'], + inputs=['yref', 'vref'], + + # System outputs + outlist=['vehicle.x', 'vehicle.y', 'vehicle.theta', 'controller.v', + 'controller.phi'], + outputs=['x', 'y', 'theta', 'v', 'phi'] +) + +# Set up the simulation conditions +yref = 1 +T = np.linspace(0, 5, 100) + +# Set up a figure for plotting the results +fbs.figure('mlh') + +# Plot the reference trajectory for the y position +plt.plot([0, 5], [yref, yref], 'k-', linewidth=0.6) + +# Find the signals we want to plot +y_index = steering.find_output('y') +v_index = steering.find_output('v') + +# Do an iteration through different speeds +for vref in [8, 10, 12]: + # Simulate the closed loop controller response + tout, yout = ct.input_output_response( + steering, T, [vref * np.ones(len(T)), yref * np.ones(len(T))]) + + # Plot the reference speed + plt.plot([0, 5], [vref, vref], 'k-', linewidth=0.6) + + # Plot the system output + y_line, = plt.plot(tout, yout[y_index, :], 'r-') # lateral position + v_line, = plt.plot(tout, yout[v_index, :], 'b--') # vehicle velocity + +# Add axis labels +plt.xlabel('Time [s]') +plt.ylabel(r'$\dot x$ [m/s], $y$ [m]') +plt.legend((v_line, y_line), (r'$\dot x$', '$y$'), + loc='center right', frameon=False) + +# Save the figure +fbs.savefig('figure-8.13-steering_gainsched.png') # PNG for web +fbs.savefig('steering-gainsched.eps') # EPS for book diff --git a/examples/example-template.py b/examples/example-template.py new file mode 100644 index 000000000..6750effc5 --- /dev/null +++ b/examples/example-template.py @@ -0,0 +1,34 @@ +# example-N.mm-short_title.py - +# RMM, dd MMM yyyy +# +# Figure 3.15: + +import control as ct +import numpy as np +import matplotlib.pyplot as plt +ct.use_fbs_defaults() + +# +# System dynamics +# + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) + +# +# (a) Description +# + +ax = fig.add_subplot(gs[0, 0]) # first row, first column +ax.set_title("(a) Description") + +# +# (b) Description +# + +ax = fig.add_subplot(gs[0, 1]) # first row, second column +ax.set_title("(b) Description") + +# Save the figure +plt.savefig("figure-N.mm-short_title.png", bbox_inches='tight') diff --git a/examples/fbs.py b/examples/fbs.py new file mode 100644 index 000000000..356c1eed2 --- /dev/null +++ b/examples/fbs.py @@ -0,0 +1,29 @@ +# fbs.py - FBS customization +# RMM, 8 Oct 2021 + +import matplotlib.pyplot as plt + +# Set the fonts to match the main text +plt.rc('font', family='Times New Roman', weight='normal', size=11) +plt.rcParams['mathtext.fontset'] = 'cm' + + +# Create a new figure of a specified size +def figure(size='mlh'): + if size == 'mlf' or size == '111': + plt.figure(figsize=[3.4, 5.1]) + elif size == 'mlh' or size == '221': + plt.figure(figsize=[3.4, 2.55]) + elif size == '211': + plt.figure(figsize=[6.8, 2.55]) + elif size == 'mlt' or size == '321': + plt.figure(figsize=[3.4, 1.7]) + else: + raise ValueError("unknown figure size") + return plt.gca() + + +# Print a figure +def savefig(name, pad=0.1, **kwargs): + plt.tight_layout(pad=pad) # clean up plots + plt.savefig(name) # save to file diff --git a/examples/figure-1.11-cruise_robustness.py b/examples/figure-1.11-cruise_robustness.py new file mode 100644 index 000000000..dc6bdaf3b --- /dev/null +++ b/examples/figure-1.11-cruise_robustness.py @@ -0,0 +1,58 @@ +# figure-1.11-cruise_robustness.py - Cruise control w/ range of vehicle masses +# RMM, 20 Jun 2021 +# +# This script generates the response of the system to a 4 deg hill with +# different vehicles masses. The figure shows how the velocity changes when +# the car travels on a horizontal road and the slope of the road changes to +# a constant uphill slope. The three different curves correspond to +# differing masses of the vehicle, between 1200 and 2000 kg, demonstrating +# that feedback can indeed compensate for the changing slope and that the +# closed loop system is robust to a large change in the vehicle +# characteristics. + +# Package import +import numpy as np +import matplotlib.pyplot as plt +import control as ct +import cruise # vehicle dynamics, PI controller + +# Define the time and input vectors +T = np.linspace(0, 25, 101) +vref = 20 * np.ones(T.shape) +gear = 4 * np.ones(T.shape) +theta0 = np.zeros(T.shape) + +# Now simulate the effect of a hill at t = 5 seconds +theta_hill = np.array([ + 0 if t <= 5 else + 4./180. * np.pi * (t-5) if t <= 6 else + 4./180. * np.pi for t in T]) + +# Create the plot and add a line at the reference speed +plt.subplot(2, 1, 1) +plt.axis([0, T[-1], 18.75, 20.25]) +plt.plot([T[0], T[-1]], [vref[0], vref[-1]], 'k-') # reference velocity +plt.plot([5, 5], [18.75, 20.25], 'k:') # disturbance start + +masses = [1200, 1600, 2000] +for i, m in enumerate(masses): + # Compute the equilibrium state for the system + X0, U0 = ct.find_eqpt( + cruise.cruise_PI, + [vref[0], 0], + [vref[0], gear[0], theta0[0]], + iu=[1, 2], y0=[vref[0], 0], iy=[0], params={'m': m}) + + # Simulate the effect of a hill + t, y = ct.input_output_response( + cruise.cruise_PI, T, + [vref, gear, theta_hill], + X0, params={'m':m}) + + # Plot the response for this mass + plt.plot(t, y[cruise.cruise_PI.find_output('v')], label='m = %d' % m) + +# Add labels and legend to the plot +plt.xlabel('Time [s]') +plt.ylabel('Velocity [m/s]') +plt.legend() diff --git a/examples/figure-1.18-airfuel_selectors.py b/examples/figure-1.18-airfuel_selectors.py new file mode 100644 index 000000000..b753b1c8f --- /dev/null +++ b/examples/figure-1.18-airfuel_selectors.py @@ -0,0 +1,97 @@ +# figure-1.18-airfuel_selectors.py - Air-fuel control example +# RMM, 20 Jun 2021 +# +# Air–fuel controller based on selectors. The right figure shows a +# simulation where the power reference r is changed stepwise at t = 1 and t +# = 15. Notice that the normalized air flow is larger than the normalized +# fuel flow both for increasing and decreasing reference steps. + +# Package import +import numpy as np +import matplotlib.pyplot as plt +import control as ct +import cruise + +# +# Air and fuel (oil) dynamics and controllers +# +# These dynamics come from Karl Astrom and are embedded in a SIMULINK +# diagram used for the initial part of the book. The basic structure for +# both the air and fuel controllers is a PI controller with output feedback +# for the proportional term and integral feedback on the error. This cuts +# the feedthrough term for the proportional feedback and gives a smoother +# step response (see Figure 11.1b for the basic structure). +# + +# Min selector for oil PI controller input +min_block = ct.nlsys( + updfcn=None, outfcn=lambda t, x, u, params: min(u), + name='min', inputs=['u1', 'u2'], outputs='y') + +# Max selector for air PI controller input +max_block = ct.nlsys( + updfcn=None, outfcn=lambda t, x, u, params: max(u), + name='max', inputs=['u1', 'u2'], outputs='y') + +# Oil and air flow dynamics (from KJA SIMULINK diagram) +Po = ct.tf([1], [1, 1]) +Pa = ct.tf([4], [1, 4]) + +# PI controller for oil flow +kpo = 2; kio = 4 +Cio = ct.tf([kio], [1, 0]) +Cpo = kpo +oil_block = ct.tf( + Po * Cio / (1 + Po * (Cio + Cpo)), + name="oil", inputs='r', outputs='y') + +# PI controller for air flow +kpa = 1; kia = 1 +Cia = ct.tf([kia], [1, 0]) +Cpa = kpa +air_block = ct.tf( + Pa * Cia / (1 + Pa * (Cia + Cpa)), + name="air", inputs='r', outputs='y') + +# +# Air-fuel selector dynamics +# +# The selector dynamics are based on the diagram Figure 1.18a, where we have +# already pre-computing the transfer function around the process/controller +# pairs (so the air and oil blocks have input 'R' and output 'Y' from the +# diagram). We use the interconnect function along with named signals to +# set everything up. +# + +airfuel = ct.interconnect( + [min_block, max_block, oil_block, air_block], + connections = ( + ['oil.r', 'min.y'], + ['air.r', 'max.y'], + ['min.u2', 'air.y'], + ['max.u1', 'oil.y']), + inplist = [['min.u1', 'max.u2']], inputs='ref', + outlist = ['air.y', 'oil.y'], outputs=['air', 'oil']) + +# +# Input/output response +# +# Finally, we simulate the dynamics with an input singla as showin in Figure +# 1.18b, consisting of a step increase from 0 to 1 at time t = 1 sec and +# then a decrease from 1 to 0.5 at time t = 15 sec. +# + +T = np.linspace(0, 30, 101) +ref = np.array([ + 0 if t <= 1 else + 1 if t <= 15 else + 0.5 for t in T]) +t, y = ct.input_output_response(airfuel, T, ref) + +# Plot the results +plt.subplot(2, 2, 1) +plt.plot(t, ref, t, y[0], t, y[1]) + +plt.legend(['ref', 'air', 'fuel']) +plt.xlabel('Time $t$ [sec]') +plt.ylabel('Normalized signals') diff --git a/examples/figure-2.11-2dof_stepresp.py b/examples/figure-2.11-2dof_stepresp.py new file mode 100644 index 000000000..b7c0346d8 --- /dev/null +++ b/examples/figure-2.11-2dof_stepresp.py @@ -0,0 +1,60 @@ +# figure-2.11-2dof_stepresp.py - step responses for two DOF system +# RMM, 21 Jun 2021 +# +# Response to a step change in the reference signal for a system with a PI +# controller having two degrees of freedom. The process transfer function is +# P(s) = 1/s and the controller gains are kp = 1.414, ki = 1, and β = 0, +# 0.5, and 1. +# + +import numpy as np +import matplotlib as mpl +import matplotlib.pyplot as plt +import control as ct + +s = ct.TransferFunction.s # define the differentiation operator +mpl.rcParams['text.usetex'] = True # use LaTeX for formatting legend strings + +# Process model: integrator dynamics +P = ct.tf([1], [1, 0]) + +# Set the simulation time vector +time = np.linspace(0, 10, 100) + +# +# Beta sweep +# + +# Choose gains to use +beta_list = [0, 0.5, 1] +kp = 1.414 +ki = 1 + +for beta in beta_list: + C1 = beta * kp + ki / s + C2 = (1 - beta) * kp + + Gyr = P * C1 / (1 + P * (C1 + C2)) + Gur = C1 / (1 + P * (C1 + C2)) + + t, y = ct.step_response(Gyr, time) + t, u = ct.step_response(Gur, time) + + if 'w_y_ax' not in locals(): + w_y_ax = plt.subplot(3, 2, 1) + plt.plot([time[0], time[-1]], [1, 1], 'k-') + plt.xlabel('Time $t$') + plt.ylabel('Output $y$') + w_y_ax.plot(t, y) + + if 'w_u_ax' not in locals(): + w_u_ax = plt.subplot(3, 2, 2) + plt.xlabel('Time $t$') + plt.ylabel('Input $u$') + w_u_ax.plot(t, u, label=r"$\beta = %g$" % beta) + +# Label the omega sweep curves +w_u_ax.legend(loc="upper right") + +# Overalll figure formating +plt.tight_layout() diff --git a/examples/figure-2.12,14-static_nlsys.py b/examples/figure-2.12,14-static_nlsys.py new file mode 100644 index 000000000..ed26442ef --- /dev/null +++ b/examples/figure-2.12,14-static_nlsys.py @@ -0,0 +1,136 @@ +# figure-2.12,14-static_nlsys.py - static nonlinear feedback system +# RMM, 21 Jun 2021 +# +# Figure 2.12: Responses of a static nonlinear system. The left figure shows +# the in- put/output relations of the open loop systems and the right figure +# shows responses to the input signal (2.38). The ideal response is shown +# with solid bold lines. The nominal response of the nonlinear system is +# shown using dashed bold lines and the responses for different parameter +# values are shown using thin lines. Notice the large variability in the +# responses. +# +# Figure 2.14: Responses of the systems with integral feedback (ki = +# 1000). The left figure shows the input/output relationships for the closed +# loop systems, and the center figure shows responses to the input signal +# (2.38) (compare to the corresponding responses in Figure 2.12. The right +# figure shows the individual errors (solid lines) and the approximate error +# given by equation (2.42) (dashed line). +# +# Intial code contributed by Adam Matic, 26 May 2021. +# + +import numpy as np +import matplotlib.pyplot as plt +import control as ct + +# Static nonlinearity +def F(u, alpha, beta): + return alpha * (u + beta * (u ** 3)) + +# Reference signal +t = np.linspace(0, 6, 300) +r = np.sin(t) + np.sin(np.pi * t) + np.sin((np.pi**2) * t) + +# +# Open loop response +# + +# Set up the plots for Figure 2.12 +fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(6, 3)) + +# Generate the input/output curves and system responses +for a in [0.1, 0.2, 0.5]: + for b in [0, 0.5, 1, 2]: + y = F(r, a, b) + ax1.plot(r, y, 'r', linewidth=0.5); + ax2.plot(t, y, 'r', linewidth=0.5) + +# Generate the nominal response +y = F(r, 0.2, 1) +ax1.plot(r, y, 'b--', linewidth=1.5); +ax2.plot(t, y, 'b--', linewidth=1.5) + +# Left plot labels +ax1.set_title("I/O relationships") +ax1.set_xlabel("Input $u$") +ax1.set_ylabel("Output $y$") + +# Draw reference line, set axis limits +ax1.plot(y, y, 'k-', linewidth=1.5) +ax1.set_ylim(-3, 3) +ax1.set_xlim(-2.5, 2.5) + +# Right plot labels +ax2.set_title("Output signals") +ax2.set_xlabel("Time $t$") +ax2.set_ylabel("Output $y$") + +# Draw reference line, set axis limits +ax2.plot(t, r, 'k-', linewidth=1.5) +ax2.set_ylim(-1, 5) +ax2.set_xlim(0, 2) + +plt.tight_layout() + +# +# Closed loop response +# + +# Create an I/O system representing the static nonlinearity +P = ct.nlsys( + updfcn=None, + outfcn=lambda t, x, u, params: F(u, params['a'], params['b']), + inputs=['u'], outputs=['y'], name='P') + +# Integral controller +ki = 1000 +C = ct.tf([ki], [1, 0]) + +# Closed loop system +sys = ct.feedback(P * C, 1) + +# Set up the plots for Figure 2.12 +fig, [ax1, ax2, ax3] = plt.subplots(1, 3, figsize=(6, 3)) + +# Generate the input/output curves and system responses +for a in [0.1, 0.2, 0.5]: + for b in [0, 0.5, 1, 2]: + # Simulate the system dynamics + t, y = ct.input_output_response(sys, t, r, params={'a':a, 'b':b}) + + ax1.plot(r, y, 'r', linewidth=0.5); + ax2.plot(t, y, 'r', linewidth=0.5) + ax3.plot(t, r-y, 'r', linewidth=0.5) + +# Left plot labels +ax1.set_title("I/O relationships") +ax1.set_xlabel("Input $u$") +ax1.set_ylabel("Output $y$") + +# Draw reference line, set axis limits +ax1.plot(y, y, 'k-', linewidth=1.5) +ax1.set_ylim(-3, 3) +ax1.set_xlim(-2.5, 2.5) + +# Center plot labels +ax2.set_title("Output signals") +ax2.set_xlabel("Time $t$") +ax2.set_ylabel("Output $y$") + +# Draw reference line, set axis limits +ax2.plot(t, r, 'k-', linewidth=1) +ax2.set_ylim(-1, 5) +ax2.set_xlim(0, 2) + +# Right plot labels +ax3.set_title("Error") +ax3.set_xlabel("Time $t$") +ax3.set_ylabel("Error $e$") + +# Draw bounding line, set axis limits +rdot = np.diff(r)/(t[1] - t[0]) # Approximation of derivative +bmin = 0.1 # See FBS2e, below equation (2.40) +ax3.plot(t[:-1], rdot/(bmin * ki), 'b--', linewidth=1.5) +ax3.set_xlim(0, 2) + +plt.tight_layout() diff --git a/examples/figure-2.19-posfbk_saturation.py b/examples/figure-2.19-posfbk_saturation.py new file mode 100644 index 000000000..d96f1500a --- /dev/null +++ b/examples/figure-2.19-posfbk_saturation.py @@ -0,0 +1,115 @@ +# figure-2.19-posfbk_saturation.py - positive feedback with saturation +# RMM, 22 Jun 2021 +# +# Figure 2.19: System with positive feedback and saturation. (a) For a +# fixed reference value r, the intersections with the curve r = G(y) +# corresponds to equilibrium points for the system. Equilibrium points +# at selected values of r are shown by circles (note that for some +# reference values there are multiple equilibrium points). Arrows +# indicate the sign of the derivative of y away from the equilibrium +# points, with the solid portions of r = G(y) representing stable +# equilibrium points and dashed portions representing unstable +# equilibrium points. (b) The hysteretic input/output map given by y = +# G+(r), showing that some values of r have single equilibrium points +# while others have two possible (stable) steady-state output val- +# ues. (c) Simulation of the system dynamics showing the reference r +# (dashed curve) and the output y (solid curve). + +import numpy as np +import matplotlib.pyplot as plt +import control as ct + +# Nonlinearity: forward and reverse +def F(x): + return x / (1 + np.abs(x)) + +def Finv(y): + return y / (1 - np.abs(y)) + +# Equilibrium calculation +a = 1; b = 4; # parameters for the system dynamics +def G(y, a=a, b=b): + return a * Finv(y) / b - y + +# +# Stable and unstable equilibrium points +# +plt.subplot(2, 2, 1) +plt.title('Stable and unstable eq points') +plt.xlabel('$y$') +plt.ylabel('$r = G(y)$') + +# Define the stable and unstable branches +ymax = -1 + np.sqrt(a/b) # First maximum (negative value) +y_stable = np.linspace(-ymax, 0.85, 100) +y_unstable = np.linspace(ymax, -ymax, 100) + +# Plot the locus of equilibrium piots +plt.plot(-y_stable, G(-y_stable), 'b-') +plt.plot(y_unstable, G(y_unstable), 'b--') +plt.plot(y_stable, G(y_stable), 'b-') + +# Plot individual equlibrium points +for r in [0, G(ymax), 0.5]: + # Find the values intersecting this value of r + y_left_arg = np.argwhere(G(-y_stable) >= r) + if y_left_arg.size > 0: + y = -y_stable[y_left_arg[-1].item()] + plt.plot(y, r, 'o', fillstyle='none') + + y_center_arg = np.argwhere(G(y_unstable) >= r) + if y_center_arg.size > 0: + y = y_unstable[y_center_arg[-1].item()] + plt.plot(y, r, 'o', fillstyle='none') + + y_right_arg = np.argwhere(G(y_stable) <= r) + if y_right_arg.size > 0: + y = y_stable[y_right_arg[-1].item()] + plt.plot(y, r, 'o', fillstyle='none') + +# +# Hysteretic input/output map y=G+(r) +# +plt.subplot(2, 2, 2) +plt.title('Hysteretic input/output map') +plt.xlabel('$r$') +plt.ylabel(r'$y = G^\dagger(r)$') + +# Plot y versus r (multi-valued) +plt.plot(G(y_stable), y_stable, 'b-') # Upper branch +plt.plot(G(-y_stable), -y_stable, 'b-') # Lower branch + +# Transition lines (look for intersection on opposite branch) +plt.plot( + [G(y_stable[0]), G(y_stable[0])], + [y_stable[0], + -y_stable[np.argwhere(G(-y_stable) > G(y_stable[0]))[-1].item()]], + 'b--') +plt.plot( + [G(-y_stable[0]), G(-y_stable[0])], + [-y_stable[0], + y_stable[np.argwhere(G(y_stable) < G(-y_stable[0]))[-1].item()]], + 'b--') + +# +# Input/output behavior +# +plt.subplot(2, 1, 2) +plt.title('Input/output behavior') + +# Closed loop dynamics +linsys = ct.tf([b], [1, a]) +nonlin = ct.nlsys( + updfcn=None, outfcn=lambda t, x, u, params: F(u), + inputs=1, outputs=1) +posfbk = ct.feedback(nonlin * linsys, 1, 1) + +# Simulate and plot +t = np.linspace(0, 100, 100) +r = 4 * np.sin(0.2 * t) + np.sin(0.1 * t) +t, y = ct.input_output_response(posfbk, t, r) +plt.plot(t, r, 'b--') +plt.plot(t, y, 'r') +plt.plot(t, np.zeros_like(t), 'k') + +plt.tight_layout() diff --git a/examples/figure-2.8-PI_step_responses.py b/examples/figure-2.8-PI_step_responses.py new file mode 100644 index 000000000..9a2f27f65 --- /dev/null +++ b/examples/figure-2.8-PI_step_responses.py @@ -0,0 +1,88 @@ +# figure-2.8-PI_step_reesponses.py - step responses for P/PI controllers +# RMM, 21 Jun 2021 +# +# Step responses for a first-order, closed loop system with proportional +# control and PI control. The process transfer function is P = 2/(s + 1). +# The controller gains for proportional control are k_p = 0, 0.5, 1, and +# 2. The PI controller is designed using equation (2.28) with zeta_c = 0.707 +# and omega_c = 0.707, 1, and 2, which gives the controller parameters k_p = +# 0, 0.207, and 0.914 and k_i = 0.25, 0.50, and 2. +# + +import numpy as np +import matplotlib.pyplot as plt +import control as ct + +# Process model +b = 2; a = 1 +P = ct.tf([b], [1, a]) + +# Set the simulation time vector +time = np.linspace(0, 8, 100) + +# +# Proportional control +# + +# Choose gains to use +kp_gains = [0, 0.5, 1, 2] + +for kp in kp_gains: + Gyv = ct.tf([b], [1, a + b*kp]) + Guv = ct.tf([-b*kp], [1, a + b*kp], dt=0) # force kp=0 to be cts time + + t, y = ct.step_response(Gyv, time) + t, u = ct.step_response(Guv, time) + + if 'p_y_ax' not in locals(): + p_y_ax = plt.subplot(3, 2, 1) + plt.ylabel('Output $y$') + plt.title('Proportional control') + p_y_ax.plot(t, y) + + if 'p_u_ax' not in locals(): + p_u_ax = plt.subplot(3, 2, 3) + plt.ylabel('Input $u$') + plt.xlabel('Normalized time $at$') + p_u_ax.plot(t, u, label="kp = %0.3g" % kp) + +# Label proportional control curves +p_u_ax.legend() + +# +# PI control +# + +# Figure out frequency of critical damping +zeta = 0.707 +wc = a / 2 / zeta + +# Plot results for different resonate frequencies +wc_list = [wc, 1, 2] +for wc in wc_list: + kp = (2 * zeta * wc - a) / b + ki = wc**2 / b + + Gyv = ct.tf([b, 0], [1, a + b*kp, b*ki]) + Guv = -ct.tf([b*kp, b*ki], [1, a + b*kp, b*ki], dt=0) + + t, y = ct.step_response(Gyv, time) + t, u = ct.step_response(Guv, time) + + if 'pi_y_ax' not in locals(): + pi_y_ax = plt.subplot(3, 2, 2) + plt.ylabel('Output $y$') + plt.title('Proportional-integral control') + pi_y_ax.plot(t, y) + + if 'pi_u_ax' not in locals(): + pi_u_ax = plt.subplot(3, 2, 4) + plt.ylabel('Input $u$') + plt.xlabel('Normalized time $at$') + pi_u_ax.plot(t, u, label="wc = %0.3g" % wc) + +# Label PI curves +pi_u_ax.legend() + +# Overall figure labeling +plt.tight_layout() diff --git a/examples/figure-2.9-secord_stepresp.py b/examples/figure-2.9-secord_stepresp.py new file mode 100644 index 000000000..0221fac1b --- /dev/null +++ b/examples/figure-2.9-secord_stepresp.py @@ -0,0 +1,93 @@ +# figure-2.9-secord_stepresp.py - step responses for second order systems +# RMM, 21 Jun 2021 +# +# Responses to a unit step change in the reference signal for different +# values of the design parameters \omega_c and \zeta_c. The left column +# shows responses for fixed \zeta_c = 0.707 and \omega_c = 1, 2, and 5. The +# right figure column responses for \omega_c = 2 and \zeta_c = 0.5, 0.707, +# and 1. The process parameters are a = b = 1. The initial value of the +# control signal is kp. +# + +import numpy as np +import matplotlib.pyplot as plt +import control as ct + +# Process model +b = 1; a = 1 +P = ct.tf([b], [1, a]) + +# Set the simulation time vector +time = np.linspace(0, 6, 100) + +# +# Omega sweep +# + +# Choose gains to use +wc_list = [1, 2, 5] +zc = 0.707 + +for wc in wc_list: + kp = (2 * zc * wc - a) / b + ki = wc**2 + C = ct.tf([kp, ki], [1, 0]) + + Gyr = P*C / (1 + P*C) + Gur = C / (1 + P*C) + + t, y = ct.step_response(Gyr, time) + t, u = ct.step_response(Gur, time) + + if 'w_y_ax' not in locals(): + w_y_ax = plt.subplot(3, 2, 1) + plt.ylabel('Output $y$') + plt.title(r"Sweep $\omega_c$, $\zeta_c = %g$" % zc) + w_y_ax.plot(t, y) + + if 'w_u_ax' not in locals(): + w_u_ax = plt.subplot(3, 2, 3) + plt.ylabel('Input $u$') + plt.xlabel(r'Normalized time $\omega_c t$') + w_u_ax.plot(t, u, label=r"$\omega_c = %g$" % wc) + +# Label the omega sweep curves +w_u_ax.legend(loc="upper right") + +# +# Zeta sweep +# + +# Figure out frequency of critical damping +wc = 2 +zc_list = [0.5, 0.707, 1] + +# Plot results for different resonate frequencies +for zc in zc_list: + kp = (2 * zc * wc - a) / b + ki = wc**2 + C = ct.tf([kp, ki], [1, 0]) + + Gyr = P*C / (1 + P*C) + Gur = C / (1 + P*C) + + t, y = ct.step_response(Gyr, time) + t, u = ct.step_response(Gur, time) + + if 'z_y_ax' not in locals(): + z_y_ax = plt.subplot(3, 2, 2) + plt.ylabel('Output $y$') + plt.title(r"Sweep $\zeta_c$, $\omega_c = %g$" % wc) + z_y_ax.plot(t, y) + + if 'z_u_ax' not in locals(): + z_u_ax = plt.subplot(3, 2, 4) + plt.ylabel('Input $u$') + plt.xlabel(r'Normalized time $\omega_c t$') + z_u_ax.plot(t, u, label=r"$\zeta_c = %g$" % zc) + +# Label the zeta sweep curves +z_u_ax.legend(loc="upper right") + +# Overalll figure labeling +plt.tight_layout() diff --git a/examples/figure-3.11-spring_mass_simulation.py b/examples/figure-3.11-spring_mass_simulation.py new file mode 100644 index 000000000..5d6fe3a29 --- /dev/null +++ b/examples/figure-3.11-spring_mass_simulation.py @@ -0,0 +1,86 @@ +# figure-3.11-spring_mass_simulation.py - forced spring–mass simulation approx +# RMM, 28 Aug 2021 +# +# Figure 3.11: Simulation of the forced spring–mass system with +# different simulation time constants. The solid line represents the +# analytical solution. The dashed lines represent the approximate +# solution via the method of Euler integration, using decreasing step +# sizes. +# + +import numpy as np +import matplotlib.pyplot as plt +import control as ct +ct.use_fbs_defaults() + +# Parameters defining the system +m = 250 # system mass +k = 40 # spring constant +b = 60 # damping constant + +# System matrices +A = [[0, 1], [-k/m, -b/m]] +B = [0, 1/m] +C = [1, 0] +sys = ct.ss(A, B, C, 0) + +# +# Discrete time simulation +# +# This section explores what happens when we discretize the ODE +# and convert it to a discrete time simulation. +# +Af = 20 # forcing amplitude +omega = 0.5 # forcing frequency + +# Sinusoidal forcing function +t = np.linspace(0, 100, 1000) +u = Af * np.sin(omega * t) + +# Simulate the system using standard MATLAB routines +response = ct.forced_response(sys, t, u) +ts, ys = response.time, response.outputs + +# +# Now generate some simulations manually +# + +# Time increments for discrete approximations +hvec = [1, 0.5, 0.1] # h must be a multiple of 0.1 +max_len = int(t[-1] / min(hvec)) # maximum number of time steps + +# Create arrays for storing results +td = np.zeros((len(hvec), max_len)) # discrete time instants +yd = np.zeros((len(hvec), max_len)) # output at discrete time instants + +# Discrete time simulations +maxi = [] # list to store maximum index +for iter, h in enumerate(hvec): + maxi.append(round(t[-1] / h)) # save maximum index for this h + x = np.zeros((2, maxi[-1] + 1)) # create an array to store the state + + # Compute the discrete time Euler approximation of the dynamics + for i in range(maxi[-1]): + offset = int(h/0.1 * i) # input offset + x[:, i+1] = x[:, i] + h * (sys.A @ x[:, i] + (sys.B * u[offset])[:, 0]) + td[iter, i] = (i-1) * h + yd[iter, i] = (sys.C @ x[:, i]).item() + +# Plot the results +plt.subplot(2, 1, 1) + +simh = plt.plot( + td[0, 0:maxi[0]], yd[0, 0:maxi[0]], 'g+--', + td[1, 0:maxi[1]], yd[1, 0:maxi[1]], 'ro--', + td[2, 0:maxi[2]], yd[2, 0:maxi[2]], 'b--', + markersize=4, linewidth=1 +) +analh = plt.plot(ts, ys, 'k-', linewidth=1) + +plt.xlabel("Time [s]") +plt.ylabel("Position $q$ [m]") +plt.axis([0, 50, -2, 2]) +plt.legend(["$h = %g$" % h for h in hvec] + ["analytical"], loc="lower left") + +# Save the figure +plt.savefig("figure-3.11-spring_mass_simulation.png", bbox_inches='tight') diff --git a/examples/figure-3.12-frequency_response.py b/examples/figure-3.12-frequency_response.py new file mode 100644 index 000000000..6707d80df --- /dev/null +++ b/examples/figure-3.12-frequency_response.py @@ -0,0 +1,106 @@ +# figure-3.12-frequency_response.py - frequency response computed by simulation +# +# Figure 3.12: A frequency response (gain only) computed by measuring +# the response of individual sinusoids. The figure on the left shows +# the response of the system as a function of time to a number of +# different unit magnitude inputs (at different frequencies). The +# figure on the right shows this same data in a different way, with +# the magnitude of the response plotted as a function of the input +# frequency. The filled circles correspond to the particular +# frequencies shown in the time responses. +# + +import numpy as np +import matplotlib.pyplot as plt +import control as ct +ct.use_fbs_defaults() + +# System definition - third order, state space system +A = [[-0.2, 2, 0], [-0.5, -0.2, 4], [0, 0, -10]] +B = [0, 0, 1] +C = [2.6, 0, 0] +sys = ct.ss(A, B, C, 0) * 1.4 # state space object (with tweaked scale) + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(2, 2) + +# +# (a) The response of the system as a function of time to a number of +# different unit magnitude inputs (at different frequencies). +# + +fig.add_subplot(gs[0, 0]) # first row, first column + +# List of frequencies for the time simulations (and frequency response points) +omega_time = [0.1, 0.4, 1, 3] +mag_time = [] # list to store magnitude of responses + +# Manual computation of the frequency response +for omega in omega_time: + # Compute out the time vector and inputs + t = np.linspace(0, 50, 1000) + u = np.sin(omega * t) + + # Simulate the system + response = ct.forced_response(sys, t, u) + + # Plot the output + plt.plot(response.time, response.outputs, 'b-') + + # Compute the magnitude of the response (avoiding initial transient) + mag_time.append(max(response.outputs[500:])) + +# Add grid lines +plt.xticks([0, 10, 20, 30, 40, 50]) +plt.grid(which='major') + +# Label the plot +plt.xlabel("Time [s]") +plt.ylabel("Output, $y$") +plt.title("Time domain simulations") + +# +# (b) The same data in a different way, with the magnitude of the +# response plotted as a function of the input frequency. The filled +# circles correspond to the particular frequencies shown in the time +# responses. +# + +fig.add_subplot(gs[0, 1]) # first row, second column + +# List of frequencies to compute the frequency response +omega_freq = [0.1, 0.2, 0.4, 1, 1.6, 3, 8, 10] +mark_index = [1, 4, 6] # frequencies to mark on the plot +mag_freq = [] # list to store magnitude of responses + +# Manual computation of the frequency response +for omega in omega_freq: + # Compute out the time vector and inputs + t = np.linspace(0, 50, 1000) + u = np.sin(omega * t) + + # Simulate the system + response = ct.forced_response(sys, t, u) + + # Compute the magnitude of the response (avoiding initial transient) + mag_freq.append(max(response.outputs[500:])) + +# Figure out which frequency points to mark +omega_mark = np.array(omega_freq)[mark_index] +mag_mark = np.array(mag_freq)[mark_index] + +# Plot the results +plt.loglog(omega_freq, mag_freq, 'b-') +plt.loglog(omega_mark, mag_mark, 'bo', markerfacecolor='none') +plt.loglog(omega_time, mag_time, 'bo') + +plt.grid(axis='x', which='minor') +plt.grid(axis='y', which='major') + +# Label the plot +plt.xlabel("Frequency [rad/s]") +plt.ylabel("Gain (log scale)") +plt.title("Frequency response") + +plt.savefig("figure-3.12-frequency_response.png", bbox_inches='tight') diff --git a/examples/figure-3.2-state_model.py b/examples/figure-3.2-state_model.py new file mode 100644 index 000000000..783319ad5 --- /dev/null +++ b/examples/figure-3.2-state_model.py @@ -0,0 +1,99 @@ +# figure-3.2-state_mode.py - illustration of a state model +# RMM, 2 Jul 2021 +# +# Figure 3.2: Illustration of a state model. A state model gives the rate of +# change of the state as a function of the state. The plot on the left shows +# the evolution of the state as a function of time. The plot on the right, +# called a phase portrait, shows the evolution of the states relative to +# each other, with the velocity of the state denoted by arrows. +# + +import numpy as np +import scipy as sp +import matplotlib.pyplot as plt +import control as ct + +# +# Spring mass system with nonlinear dampling +# +# This function gives the dynamics for a dampled oscillator with nonlinear +# damping. The states of the system are +# +# x[0] position +# x[1] velocity +# +# The nonlinear damping is implemented as a change in the linear damping +# coefficient at a small velocity. This is intended to roughly correspond +# to some sort of stiction (and give an interesting phase portrait). The +# default parameters for the system are given by +# +# m = 1 mass, kg +# k = 1 spring constant, N/m +# b1 = 1 damping constant near origin, N-sec/m +# b2 = 0.01 damping constant away from origin, N-sec/m +# dth = 0.5 threshold for switching between damping +# +# This corresponds to a fairly lightly damped oscillator away from the origin. + +def _nlspringmass(t, x, u, params): + m = params.get('m', 1) + k = params.get('k', 1) + b1 = params.get('b1', 2) + b2 = params.get('b2', 0.01) + dth = params.get('dth', 0.2) + + # Compute the friction force + if abs(x[1]) < dth: + Fb = b1 * x[1]; + elif x[1] < 0: + Fb = -b1 * dth \ + + b2 * (x[1] + dth); + else: + Fb = b1 * dth \ + + b2 * (x[1] - dth); + + # Return the time derivative of the state + return np.array([x[1], -k/m * x[0] - Fb/m]) +nlspringmass = ct.nlsys(_nlspringmass, None, states=2, inputs=0, outputs=2) + +# +# (a) Simulation of the nonlinear spring mass system +# +plt.subplot(2, 2, 1) + +t = np.linspace(0, 16, 100) +resp = ct.input_output_response(nlspringmass, t, 0, [2, 0]) +y = resp.outputs + +plt.plot(t, y[0], '-', t, y[1], '--') +plt.xlabel('Time $t$ [s]') +plt.ylabel(r'Position $q$ [m], velocity $\dot q$̇ [m/s]') +plt.title('Time plot') +plt.legend(['Position $q$', 'Velocity $v$']) + +# +# (b) Generate a phase plot for the damped oscillator +# +ax = plt.subplot(2, 2, 2) + +cplt = ct.phase_plane_plot( + nlspringmass, # dynamics + [-1, 1, -1, 1], # bounds of the plot + gridspec=[8, 8], # number of points for vectorfield + plot_vectorfield=True, # plot vectorfield + plot_streamlines=False, # plot streamlines separately + plot_separatrices=False, # leave off separatrices + ax=ax +) +ct.phaseplot.streamlines( # Plot streamlines from selected points + nlspringmass, + np.array([[-1, 0.4], [0.1, 1], [1, -0.4], [-0.1, -1]]), + 10, ax=ax +) + +plt.xlabel('Position $q$ [m]') +plt.ylabel(r'Velocity $\dot q$ [m/s]') +plt.title('Phase portrait') +plt.axis([-1, 1, -1, 1]) + +plt.tight_layout() diff --git a/examples/figure-3.4-io_response.py b/examples/figure-3.4-io_response.py new file mode 100644 index 000000000..b39a9c969 --- /dev/null +++ b/examples/figure-3.4-io_response.py @@ -0,0 +1,65 @@ +# figure-3-1-io_response.py - input/output response of a linear system +# RMM, 28 Aug 2021 +# +# Figure 3.4: Input/output response of a linear system. The step +# response (a) shows the output of the system due to an input that +# changes from 0 to 1 at time t = 5 s. The frequency response (b) +# shows the amplitude gain and phase change due to a sinusoidal input +# at different frequencies. +# + +import numpy as np +import scipy as sp +import matplotlib.pyplot as plt +import control as ct +ct.use_fbs_defaults() # Use settings to match FBS + +# System definition - third order, state space system +A = [[-0.2, 2, 0], [-0.5, -0.2, 4], [0, 0, -10]] +B = [0, 0, 1] +C = [2.6, 0, 0] +sys = ct.ss(A, B, C, 0) # state space object + +# Set up the plotting grid to match the layout in the book +fig = plt.figure(constrained_layout=True) +gs = fig.add_gridspec(4, 2) + +# +# (a) Step response showing the output of the system due to an input +# that changes from 0 to 1 at time t = 5 s +# + +fig.add_subplot(gs[0:2, 0]) # first column + +# Create an input signal that is zero until time t = 5 +t = np.linspace(0, 30, 100) +u = np.ones_like(t) +u[t < 5] = 0 + +# Compute the response +response = ct.forced_response(sys, t, u) +y = response.outputs + +# Plot the response +plt.plot(t, u, 'b--', label="Input") +plt.plot(t, y, 'r-', label="Output") +plt.xlabel("Time (sec)") +plt.ylabel("Input, output") +plt.title("Step response") +plt.legend() + +# +# (b) Frequency` response showing the amplitude gain and phase change +# due to a sinusoidal input at different frequencies +# + +# Set up the axes for plotting (labels are recognized by bode_plot()) +mag = fig.add_subplot(gs[0, 1], label='control-bode-magnitude') +phase = fig.add_subplot(gs[1, 1], label='control-bode-phase') + +# Generate the Bode plot +ct.bode_plot(sys) + +# Adjust the appearance to match the book +mag.xaxis.set_ticklabels([]) +mag.set_title("Frequency response") diff --git a/examples/figure-4.12-congctrl_eqplot.py b/examples/figure-4.12-congctrl_eqplot.py new file mode 100644 index 000000000..bd034dabc --- /dev/null +++ b/examples/figure-4.12-congctrl_eqplot.py @@ -0,0 +1,36 @@ +# congctrl_eqplot.py - congestion control equilibrium point plot +# RMM, 29 Jun 2007 (converted from MATLAB) +# +# The equilibrium buffer size be for a set of N identical computers sending +# packets through a single router with drop probability ρb. +# + +import matplotlib.pyplot as plt +import numpy as np +import scipy.optimize +import fbs # FBS plotting customizations + +# Range of values to plot (\alpha = 1/(2\rho^2 N^2) +alpha_vals = np.logspace(-2, 4) + +# Solve for the equilibrium value of \rho b_e +bratio_vals = [] +for alpha in alpha_vals: + # Define a function for the equilibrium point (equation (4.22)) + def equilibrium(bratio): + return alpha * bratio**3 + bratio - 1 + + bratio = scipy.optimize.fsolve(equilibrium, 0) + bratio_vals.append(bratio) + +# Set up a figure for plotting the results +fbs.figure('mlh') + +# Plot the equilibrium buffer length +plt.semilogx(alpha_vals, bratio_vals) +plt.xlabel(r"$1/(2 \rho^2 N^2)$ (log scale)") +plt.ylabel(r"$\rho b_{e}$") +plt.title("Operating point") + +# Save the figure +fbs.savefig('figure-4.12-congctrl_eqplot.png') # PNG for web diff --git a/examples/figure-4.13-congctrl_tcpsim.py b/examples/figure-4.13-congctrl_tcpsim.py new file mode 100644 index 000000000..a5a0e15a7 --- /dev/null +++ b/examples/figure-4.13-congctrl_tcpsim.py @@ -0,0 +1,64 @@ +# tcpsim.m - congestion control simulation +# RMM, 9 Sep 2006 (from MATLAB) + +import matplotlib.pyplot as plt +import numpy as np +import random +import control as ct +import congctrl +import fbs # FBS plotting customizations + +# Create an I/O system for simulation +M = 6 # simulate 6 aggregated sources +N = 60 # modeling 60 independent sources +congctrl_sys = congctrl.create_iosystem(M, N=N) + +# Find the equilibrium point +xeq, ueq = ct.find_eqpt(congctrl_sys, np.ones(M+1), 0) +weq, beq = xeq[0], xeq[-1] # first N states are identical + +# Compute a perturbation for the initial condition +random.seed(6) # for repeatability +w0 = np.array([ + weq + 2 * (random.random() - 0.5) * weq for i in range(M)]) +b0 = beq/2 + +# Run a simulation +tvec = np.linspace(0, 500, 100) +resp = ct.input_output_response( + congctrl_sys, tvec, U=0, X0=[w0, b0]) + +# Plot the results +# Set up a figure for plotting the results +fbs.figure('mlh') + +for i in range(M): + plt.plot(resp.time, resp.states[i], 'k') +plt.plot(resp.time, resp.states[-1] / 20, 'b') + +# Now change the number of sources and rerun from the old initial condition +M = 4 # simulate 4 aggregated sources +N = 40 # modeling 40 independent sources +congctrl_pert = congctrl.create_iosystem(M, N=N) + +# Run a simulation starting from where we left off +tvec = np.linspace(500, 1000, 100) +pert = ct.input_output_response( + congctrl_pert, tvec, U=0, X0=[resp.states[0:4, -1], resp.states[-1, -1]]) + +# Plot the results +for i in range(M): + plt.plot(pert.time, pert.states[i], 'k') +plt.plot(pert.time, pert.states[-1] / 20, 'b') + +# Label the plots and make them pretty +plt.axis([0, 1000, 0, 20]) +plt.xlabel("Time $t$ [ms]") +plt.ylabel("States $w_i$ [pkts/ms], $b$ [pkts]") +plt.title("Time response") +plt.text(250, 3, "$w_1$ - $w_{60}$") +plt.text(700, 4, "$w_1$ - $w_{40}$") +plt.text(700, 10, "$b$") + +# Save the figure +fbs.savefig('figure-4.13-congctrl_tcpsim.png') # PNG for web diff --git a/examples/figure-4.20-predprey_ctstime.py b/examples/figure-4.20-predprey_ctstime.py new file mode 100644 index 000000000..8e74da351 --- /dev/null +++ b/examples/figure-4.20-predprey_ctstime.py @@ -0,0 +1,84 @@ +# predprey_ctstime.py - Predator-prey model in continuous time +# RMM, 28 May 2023 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# Define the dynamis for the predator-prey system (no input) +def predprey_update(t, x, u, params={}): + """Predator prey dynamics""" + r = params.get('r', 1.6) + d = params.get('d', 0.56) + b = params.get('b', 0.6) + k = params.get('k', 125) + a = params.get('a', 3.2) + c = params.get('c', 50) + + # Dynamics for the system + dx0 = r * x[0] * (1 - x[0]/k) - a * x[1] * x[0]/(c + x[0]) + dx1 = b * a * x[1] * x[0] / (c + x[0]) - d * x[1] + + return np.array([dx0, dx1]) + +# Create a nonlinear I/O system +predprey_sys = ct.nlsys(predprey_update, states=2) + +# Simulate a trajectory leading to a limit cycle +timepts = np.linspace(0, 70, 500) +sim = ct.input_output_response(predprey_sys, timepts, 0, [25, 20]) + +# Plot the results +fbs.figure('mlh') # FBS conventions +plt.plot(sim.time, sim.states[0], 'b-', label="Hare") +plt.plot(sim.time, sim.states[1], 'r--', label="Lynx") +plt.legend() +plt.xlabel("Time $t$ [years]") +plt.ylabel("Population") +plt.title("Time response") + +# Save the figure +fbs.savefig('figure-4.20-predprey_ctstime-sim.png') # PNG for web + +# Generate a phase portrait +fbs.figure('mlh') +ct.phaseplot.equilpoints(predprey_sys, [-5, 126, -5, 100]) +ct.phaseplot.streamlines( + predprey_sys, np.array([ + [0, 100], [1, 0], + ]), 10, color='b') +ct.phaseplot.streamlines( + predprey_sys, np.array([[124, 1]]), np.linspace(0, 10, 500), color='b') +ct.phaseplot.streamlines( + predprey_sys, np.array([[125, 25], [125, 50], [125, 75]]), 3, color='b') +ct.phaseplot.streamlines(predprey_sys, np.array([2, 8]), 6, color='b') +ct.phaseplot.streamlines( + predprey_sys, np.array([[20, 30]]), np.linspace(0, 65, 500), + gridtype='circlegrid', gridspec=[2, 1], arrows=10, color='r') +ct.phaseplot.vectorfield(predprey_sys, [5, 125, 5, 100], gridspec=[20, 20]) + +# Add the limit cycle +resp1 = ct.initial_response(predprey_sys, np.linspace(0, 100), [20, 75]) +resp2 = ct.initial_response( + predprey_sys, np.linspace(0, 20, 500), resp1.states[:, -1]) +plt.plot(resp2.states[0], resp2.states[1], color='k') + +# Legacy code +# def pp_ode(x, t): +# return predprey_update(t, x, 0, {}) +# ct.phase_plot(pp_ode, [0, 60, 7], [0, 50, 6]) +# ct.phase_plot(pp_ode, [0, 60, 7], [60, 100, 4]) +# ct.phase_plot(pp_ode, [70, 120, 6], [0, 50, 6]) +# ct.phase_plot(pp_ode, [70, 120, 6], [60, 100, 4]) + +# # Plot the limit cycle +# ct.phase_plot(pp_ode, X0=sim.states[:, -1:].T, T=20) # limit cycle +# ct.phase_plot(pp_ode, X0=[[120, 32], [120, 60]], T=20) # outside trajectories +# ct.phase_plot(pp_ode, X0=[[19, 30]], T=75) # inside trajectories + +# Label the plot +plt.xlabel("Hares") +plt.ylabel("Lynxes") +plt.title("Phase portrait") +fbs.savefig('figure-4.20-predprey_ctstime-pp.png') # PNG for web diff --git a/examples/figure-5-6-lyapunov_stability.py b/examples/figure-5-6-lyapunov_stability.py new file mode 100644 index 000000000..f81f97a0d --- /dev/null +++ b/examples/figure-5-6-lyapunov_stability.py @@ -0,0 +1,30 @@ +# lyapunov_stability.py - illustration of Lyapunov stability +# RMM, 6 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +t = np.linspace(0, 6) +x0 = np.sin(t) + 0.8 * np.cos(2.2 * t) + 2.5 + +# Plot the centerline and bounds +fbs.figure('211') +plt.plot(t, x0, 'k') +plt.plot(t, x0 + 0.5, 'r', t, x0 - 0.5, 'r') + +# Plot the signal +x = x0 - 0.4 * np.sin(t - 0.3) +plt.plot(t, x, 'b--') + +# Label the axes +plt.xlabel("Time $t$") +plt.ylabel("State $x$") + +# Add some arrows and label the range of stability +plt.arrow(1.5, 1.3, 0, 0.8, width=0.01, head_width=0.05) +plt.arrow(1.5, 4.15, 0, -0.8, width=0.01, head_width=0.05) +plt.text(1.6, 1.4, "$\\epsilon$") + +fbs.savefig('figure-5.6-lyapunov_stability.png') diff --git a/examples/figure-5.1-damposc_response.py b/examples/figure-5.1-damposc_response.py new file mode 100644 index 000000000..909dec7d2 --- /dev/null +++ b/examples/figure-5.1-damposc_response.py @@ -0,0 +1,31 @@ +# damposc_response.py - dampled oscillator response +# RMM, 28 May 2023 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# Oscillator parameters +damposc_params = {'m': 1, 'b': 0.2, 'k': 1} + +# System model (as ODE) +def damposc_update(t, x, u, params): + m, b, k = params['m'], params['b'], params['k'] + return np.array([x[1], -k/m * x[0] - b/m * x[1]]) +damposc = ct.nlsys(damposc_update, params=damposc_params) + +# Simulate the response +tvec = np.linspace(0, 20, 100) +response = ct.input_output_response(damposc, tvec, 0, X0=[1, 0]) + +# Plot the states +fbs.figure('211') +plt.plot(response.time, response.states[0], 'b-') +plt.plot(response.time, response.states[1], 'b--') +plt.plot([response.time[0], response.time[-1]], [0, 0], 'k-', linewidth=0.75) +plt.xlabel('Time $t$ [s]') +plt.ylabel('States $x_1$, $x_2$') +plt.title( + "Response of the damped oscillator to the initial condition x0 = (1, 0)") +fbs.savefig('figure-5.1-damposc_response-time.png') diff --git a/examples/figure-5.10-congctrl_dynamics.py b/examples/figure-5.10-congctrl_dynamics.py new file mode 100644 index 000000000..9d6e7673c --- /dev/null +++ b/examples/figure-5.10-congctrl_dynamics.py @@ -0,0 +1,47 @@ +# congctrl-dynamics.py - phase plots for congestion control dynamis +# RMM, 7 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# Define the system dynamics +def _congctrl_update(t, x, u, params): + # Number of sources per state of the simulation + M = x.size - 1 # general case + assert M == 1 # make sure nothing funny happens here + + # Remaining parameters + N = params.get('N', M) # number of sources + rho = params.get('rho', 2e-4) # RED parameter = pbar / (bupper-blower) + c = params.get('c', 10) # link capacity (Mp/ms) + + # Compute the derivative (last state = bdot) + return np.append( + c / x[M] - (rho * c) * (1 + (x[:-1]**2) / 2), + N/M * np.sum(x[:-1]) * c / x[M] - c) + +congctrl = ct.nlsys( + _congctrl_update, states=2, inputs=0, + params={'N': 60, 'rho': 2e-4, 'c': 10}) + +fbs.figure() +ct.phase_plane_plot(congctrl, [0, 10, 10, 500], 100) +plt.axis([0, 10, 0, 500]) +plt.suptitle("") +plt.title("$\\rho = 2 \\times 10^{-4}$, $c = 10$ pkts/msec") +plt.xlabel("Window size, $w$ [pkts]") +plt.ylabel("Buffer size, $b$ [pkts]") +fbs.savefig('figure-5.10-congctrl_dynamics-pp1.png') + +fbs.figure() +ct.phase_plane_plot( + congctrl, [0, 10, 10, 500], 100, + params={'rho': 4e-4, 'c': 20}) +plt.axis([0, 10, 0, 500]) +plt.suptitle("") +plt.title("$\\rho = 4 \\times 10^{-4}$, $c = 20$ pkts/msec") +plt.xlabel("Window size, $w$ [pkts]") +plt.ylabel("Buffer size, $b$ [pkts]") +fbs.savefig('figure-5.10-congctrl_dynamics-pp2.png') diff --git a/examples/figure-5.11-invpend_linearized.py b/examples/figure-5.11-invpend_linearized.py new file mode 100644 index 000000000..d87c694bd --- /dev/null +++ b/examples/figure-5.11-invpend_linearized.py @@ -0,0 +1,30 @@ +# invepend_linearized.py - nonlinear vs linear inverted pendulum dynamics +# RMM, 7 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +from math import pi +import control as ct +import fbs # FBS plotting customizations + +def invpend_update(t, x, u, params): + m, l, b, g = params['m'], params['l'], params['b'], params['g'] + return [x[1], -b/m * x[1] + (g * l / m) * np.sin(x[0])] +invpend = ct.nlsys( + invpend_update, states=2, inputs=0, name='inverted pendulum', + params={'m': 1, 'l': 1, 'b': 0.5, 'g': 1}) + +fbs.figure() +ct.phase_plane_plot( + invpend, [0, 2*pi, -2, 2], 6, gridspec=[6, 5], + plot_separatrices={'timedata': 20, 'arrows': 4}) +fbs.savefig('figure-5.10-invpend_linearized-nl.png') + +# Create a linearized model +linsys = invpend.linearize([pi, 0], 0) + +fbs.figure() +ct.phase_plane_plot( + linsys, [-pi, pi, -2, 2], 10, gridspec=[5, 2], + plot_separatrices={'timedata': 20, 'arrows': 4}) +fbs.savefig('figure-5.10-invpend_linearized-ln.png') diff --git a/examples/figure-5.3-phase_portraits.py b/examples/figure-5.3-phase_portraits.py new file mode 100644 index 000000000..13ead3445 --- /dev/null +++ b/examples/figure-5.3-phase_portraits.py @@ -0,0 +1,33 @@ +# phase_portraits.py - phase portrait examples +# RMM, 6 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# Oscillator parameters +damposc_params = {'m': 1, 'b': 1, 'k': 1} + +# System model (as ODE) +def damposc_update(t, x, u, params): + m, b, k = params['m'], params['b'], params['k'] + return np.array([x[1], -k/m * x[0] - b/m * x[1]]) +damposc = ct.nlsys(damposc_update, states=2, params=damposc_params) + +# Set the limits for the plot +limits = [-1, 1, -1, 1] + +# Vector field +fbs.figure('mlh') +ct.phase_plane_plot( + damposc, limits, plot_streamlines=False, + plot_vectorfield=True, gridspec=[15, 12]) +plt.suptitle(""); plt.title("Vector field") +fbs.savefig('figure-5.3-phase_portraits-vf.png') + +# Streamlines +fbs.figure('mlh') +ct.phase_plane_plot(damposc, limits, 8) +plt.suptitle(""); plt.title("Phase portrait") +fbs.savefig('figure-5.3-phase_portraits-sl.png') diff --git a/examples/figure-5.4-invpend_phaseplot.py b/examples/figure-5.4-invpend_phaseplot.py new file mode 100644 index 000000000..dcc7dcd99 --- /dev/null +++ b/examples/figure-5.4-invpend_phaseplot.py @@ -0,0 +1,21 @@ +# invepend_phaseplot.py - inverted pendulum phase plots +# RMM, 6 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +from math import pi +import control as ct +import fbs # FBS plotting customizations + +def invpend_update(t, x, u, params): + m, l, b, g = params['m'], params['l'], params['b'], params['g'] + return [x[1], -b/m * x[1] + (g * l / m) * np.sin(x[0])] +invpend = ct.nlsys( + invpend_update, states=2, inputs=0, name='inverted pendulum', + params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1}) + +fbs.figure() +ct.phase_plane_plot( + invpend, [-2*pi, 2*pi, -2, 2], 4, gridspec=[6, 6], + plot_separatrices={'timedata': 20, 'arrows': 4}) +fbs.savefig('figure-5.4-invpend_phaseplot.png') diff --git a/examples/figure-5.5-limit_cycle.py b/examples/figure-5.5-limit_cycle.py new file mode 100644 index 000000000..d2b29ac04 --- /dev/null +++ b/examples/figure-5.5-limit_cycle.py @@ -0,0 +1,38 @@ +# limit_cycle.py - nonlinear oscillator (limit cycle) phase plot +# RMM, 6 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +from math import pi +import control as ct +import control.phaseplot as pp +import fbs # FBS plotting customizations + +def oscillator_update(t, x, u, params): + return [ + x[1] + x[0] * (1 - x[0]**2 - x[1]**2), + -x[0] + x[1] * (1 - x[0]**2 - x[1]**2) + ] +oscillator = ct.nlsys(oscillator_update, states=2, inputs=0, name='oscillator') + +fbs.figure() +ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9) +pp.streamlines( + oscillator, np.array([[0, 0]]), 1.5, + gridtype='circlegrid', gridspec=[0.5, 6], dir='both') +pp.streamlines( + oscillator, np.array([[1, 0]]), 2*pi, arrows=6, color='b') +plt.gca().set_aspect('equal') +plt.suptitle("") +fbs.savefig('figure-5.5-limit_cycle-pp.png') + +fbs.figure() +plt.axis([0, 30, -2, 2]) +timepts = np.linspace(0, 30) +response = ct.input_output_response(oscillator, timepts, 0, [0.1, 1]) +plt.plot(response.time, response.outputs[0], 'b', label="$x_1$") +plt.plot(response.time, response.outputs[1], 'r--', label="$x_2$") +plt.xlabel("Time $t$") +plt.ylabel("$x_1, x_2$") +plt.legend(loc='upper right', ncols=2, frameon=False) +fbs.savefig('figure-5.5-limit_cycle-time.png') diff --git a/examples/figure-5.7-stable_eqpt.py b/examples/figure-5.7-stable_eqpt.py new file mode 100644 index 000000000..4eefe7e32 --- /dev/null +++ b/examples/figure-5.7-stable_eqpt.py @@ -0,0 +1,39 @@ +# stable_eqpt.py - plots for stable equlibrium point +# RMM, 6 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +from math import pi +import control as ct +import control.phaseplot as pp +import fbs # FBS plotting customizations + +m, b, k = 1, 0, 2 +linsys = ct.ss([[0, 1], [-k/m, -b/m]], [[0], [1]], np.eye(2), 0) + +# Draw the phase portrait +fbs.figure() +ct.phase_plane_plot(linsys, [-1, 1, -1, 1], 1, plot_streamlines=False) +pp.streamlines( + linsys, np.array([[0.2, 0], [0.4, 0], [0.6, 0], [0.8, 0], [1, 0]]), + 4.5, arrows=6) +plt.gca().set_aspect('equal') +plt.suptitle("") + +# Add some level sets +theta = np.linspace(0, 2*pi) +plt.plot(0.2 * np.sin(theta), 0.2 * np.cos(theta), 'r--') +plt.plot(0.3 * np.sin(theta), 0.3 * np.cos(theta), 'r--') + +fbs.savefig('figure-5.7-stable_eqpt-pp.png') + +fbs.figure('321') +plt.axis([0, 10, -2.5, 2.5]) +timepts = np.linspace(0, 10) +response = ct.input_output_response(linsys, timepts, 0, [1, 0]) +plt.plot(response.time, response.outputs[0], 'b', label="$x_1$") +plt.plot(response.time, response.outputs[1], 'r--', label="$x_2$") +plt.xlabel("Time $t$") +plt.ylabel("$x_1, x_2$") +plt.legend(loc='upper right', ncols=2, frameon=False) +fbs.savefig('figure-5.7-stable_eqpt-time.png') diff --git a/examples/figure-5.8-asystable_eqpt.py b/examples/figure-5.8-asystable_eqpt.py new file mode 100644 index 000000000..142c2d399 --- /dev/null +++ b/examples/figure-5.8-asystable_eqpt.py @@ -0,0 +1,28 @@ +# asystable_eqpt.py - plots for stable equlibrium point +# RMM, 7 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +m, b, k = 1, 1, 1 +linsys = ct.ss([[0, 1], [-k/m, -b/m]], [[0], [1]], np.eye(2), 0) + +# Draw the phase portrait +fbs.figure() +ct.phase_plane_plot(linsys, [-1, 1, -1, 1], 5) +plt.gca().set_aspect('equal') +plt.suptitle("") +fbs.savefig('figure-5.8-asystable_eqpt-pp.png') + +fbs.figure('321') +plt.axis([0, 10, -0.6, 1]) +timepts = np.linspace(0, 10) +response = ct.input_output_response(linsys, timepts, 0, [1, 0]) +plt.plot(response.time, response.outputs[0], 'b', label="$x_1$") +plt.plot(response.time, response.outputs[1], 'r--', label="$x_2$") +plt.xlabel("Time $t$") +plt.ylabel("$x_1, x_2$") +plt.legend(loc='upper right', ncols=2, frameon=False) +fbs.savefig('figure-5.8-asystable_eqpt-time.png') diff --git a/examples/figure-5.9-unstable_eqpt.py b/examples/figure-5.9-unstable_eqpt.py new file mode 100644 index 000000000..962edd0b4 --- /dev/null +++ b/examples/figure-5.9-unstable_eqpt.py @@ -0,0 +1,29 @@ +# unstable_eqpt.py - plots for stable equlibrium point +# RMM, 7 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +saddle = ct.ss([[1, -3], [-3, 1]], [[0], [1]], np.eye(2), 0) + +# Draw the phase portrait +fbs.figure() +ct.phase_plane_plot( + saddle, [-1, 1, -1, 1], 0.4, + gridtype='meshgrid', gridspec=[6, 6]) +plt.gca().set_aspect('equal') +plt.suptitle("") +fbs.savefig('figure-5.9-unstable_eqpt-pp.png') + +fbs.figure('321') +plt.axis([0, 3, -100, 100]) +timepts = np.linspace(0, 3) +response = ct.input_output_response(saddle, timepts, 0, [1, 0]) +plt.plot(response.time, response.outputs[0], 'b', label="$x_1$") +plt.plot(response.time, response.outputs[1], 'r--', label="$x_2$") +plt.xlabel("Time $t$") +plt.ylabel("$x_1, x_2$") +plt.legend(loc='upper right', ncols=1, frameon=False) +fbs.savefig('figure-5.9-unstable_eqpt-time.png') diff --git a/examples/figure-6.1-superposition.py b/examples/figure-6.1-superposition.py new file mode 100644 index 000000000..363fed1b0 --- /dev/null +++ b/examples/figure-6.1-superposition.py @@ -0,0 +1,49 @@ +# superposition.py - superposition of homogeneous and particular solutions +# RMM, 19 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# Spring mass system +from springmass import springmass # use spring mass dynamics +sys = springmass / springmass(0).real # normalize the response to 1 +X0 = [2, -1] # initial condition + +# Create input vectors +tvec = np.linspace(0, 60, 100) +u1 = 0 * tvec +u2 = np.hstack([tvec[0:50]/tvec[50], 1 - tvec[0:50]/tvec[50]]) + +# Run simulations for the different cases +homogeneous = ct.forced_response(sys, tvec, u1, X0=X0) +particular = ct.forced_response(sys, tvec, u2) +complete = ct.forced_response(sys, tvec, u1 + u2, X0=X0) + +# Plot results +fig, axs = plt.subplots(3, 3, figsize=[8, 4], layout='tight') +for i, resp in enumerate([homogeneous, particular, complete]): + axs[i, 0].plot(resp.time, resp.inputs) + axs[0, 0].set_title("Input $u$") + axs[i, 0].set_ylim(-2, 2) + + axs[i, 1].plot( + resp.time, resp.states[0], 'b', + resp.time, resp.states[1], 'r--') + axs[0, 1].set_title("States $x_1$, $x_2$") + axs[i, 1].set_ylim(-2, 2) + + axs[i, 2].plot(resp.time, resp.outputs) + axs[0, 2].set_title("Output $y$") + axs[i, 2].set_ylim(-2, 2) + +# Label the plots +axs[0, 0].set_ylabel("Homogeneous") +axs[1, 0].set_ylabel("Particular") +axs[2, 0].set_ylabel("Complete") +for i in range(3): + axs[2, i].set_xlabel("Time $t$ [s]") + +# Save the figure +fbs.savefig('figure-6.1-superposition.png') diff --git a/examples/figure-6.14-cruise_linearized.py b/examples/figure-6.14-cruise_linearized.py new file mode 100644 index 000000000..702de0895 --- /dev/null +++ b/examples/figure-6.14-cruise_linearized.py @@ -0,0 +1,117 @@ +# cruise_linearized.py - linear versus nonlinear response, cruise w/ PI +# RMM, 20 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# System definition +from cruise import vehicle_dynamics as vehicle + +# Figure out the equilibrium point for the system at 20 m/s +xe, ue = ct.find_eqpt(vehicle, 20, u0=[0, 4, 0], iu=[1, 2], y0=20, iy=[0]) + +# Linearized dynamics +vehicle_lin = vehicle.linearize(xe, ue) + +# Controller: PI + antiwindup +ctrl_params = {'kp': 0.5, 'ki': 0.1, 'kaw': 2} + +def ctrl_update(t, x, u, params): + e = u[1] - u[0] # v - vref + v_nom = -params['kp'] * e + x[0] # nominal control input (PI) + v_sat = np.clip(v_nom, 0, 1) # clipped control input + return -params['ki'] * e + params['kaw'] * (v_sat - v_nom) + +def ctrl_output(t, x, u, params): + e = u[1] - u[0] # v - vref + v_nom = -params['kp'] * e + x[0] # nominal control input (PI) + v_sat = np.clip(v_nom, 0, 1) # clipped control input + return v_sat + +ctrl = ct.nlsys( + ctrl_update, ctrl_output, states=1, name='ctrl', + inputs=['vref', 'v'], outputs='u', params=ctrl_params) + +# Figure out the equilibrium point for the system at 20 m/s +xe, ue = ct.find_eqpt(vehicle, 20, u0=[0, 4, 0], iu=[1, 2], y0=20, iy=[0]) + +# Linearized dynamics +vehicle_lin = vehicle.linearize( + xe, ue, inputs=vehicle.input_labels, outputs=vehicle.output_labels) + +# Controller: PI + antiwindup +ctrl_params = {'kp': 0.5, 'ki': 0.1, 'kaw': 2} + +def ctrl_update(t, x, u, params): + e = u[1] - u[0] # v - vref + v_nom = -params['kp'] * e + x[0] # nominal control input (PI) + v_sat = np.clip(v_nom, 0, 1) # clipped control input + return -params['ki'] * e + params['kaw'] * (v_sat - v_nom) + +def ctrl_output(t, x, u, params): + e = u[1] - u[0] # v - vref + v_nom = -params['kp'] * e + x[0] # nominal control input (PI) + v_sat = np.clip(v_nom, 0, 1) # clipped control input + return v_sat + +ctrl = ct.nlsys( + ctrl_update, ctrl_output, states=1, name='ctrl', + inputs=['vref', 'v'], outputs='u', params=ctrl_params) + +# Full system (linear and nonlinear) +nlsys = ct.interconnect( + [vehicle, ctrl], inputs=['vref', 'gear', 'theta'], outputs=['v', 'u']) + +lnsys = ct.interconnect( + [vehicle_lin, ctrl], inputs=['vref', 'gear', 'theta'], outputs=['v', 'u']) + +# Compute system response: flat then a hill +T1 = np.linspace(0, 5, 40) # Flat section +T2 = np.linspace(5, 30) # Hill section + +# Nonlinear response +nl_resp1 = ct.input_output_response(nlsys, T1, [20, 4, 0], X0=[20, ue[0]]) +nl_resp2a = ct.input_output_response( + nlsys, T2, [20, 4, 0.07], X0=nl_resp1.states[:, -1]) +nl_resp2b = ct.input_output_response( + nlsys, T2, [20, 4, 0.105], X0=nl_resp1.states[:, -1]) + +# Linear response +ln_resp1 = ct.input_output_response(lnsys, T1, [20, 4, 0], X0=[20, ue[0]]) +ln_resp2a = ct.input_output_response( + lnsys, T2, [20, 4, 0.07], X0=ln_resp1.states[:, -1]) +ln_resp2b = ct.input_output_response( + lnsys, T2, [20, 4, 0.105], X0=ln_resp1.states[:, -1]) + +# Plot the velocity response +fig, axs = plt.subplots(2, 1, figsize=[3.4, 3.4], sharex=True) + +axs[0].plot(nl_resp1.time, nl_resp1.outputs[0], 'b') +axs[0].plot(nl_resp2a.time, nl_resp2a.outputs[0], 'b') +axs[0].plot(nl_resp2b.time, nl_resp2b.outputs[0], 'b') + +axs[0].plot(ln_resp2a.time, ln_resp2a.outputs[0], 'r--') +axs[0].plot(nl_resp2b.time, ln_resp2b.outputs[0], '--') + +axs[0].axis([0, 30, 18.5, 20.5]) +axs[0].set_yticks([19, 20]) +axs[0].set_ylabel("Velocity $v$ [m/s]") + +# Plot the throttle command +axs[1].plot(nl_resp1.time, nl_resp1.outputs[1], 'b') +axs[1].plot(nl_resp2a.time, nl_resp2a.outputs[1], 'b') +axs[1].plot(nl_resp2b.time, nl_resp2b.outputs[1], 'b') + +axs[1].plot(ln_resp2a.time, ln_resp2a.outputs[1], 'r--') +axs[1].plot(nl_resp2b.time, ln_resp2b.outputs[1], '--') + +axs[1].set_title(" ") # Hack to adjust spacing between plots +axs[1].set_ylim([0, 1.25]) +axs[1].set_yticks([0, 0.5, 1]) +axs[1].set_yticklabels(["0", "0.5", "1"]) +axs[1].set_ylabel("Throttle $u$") +axs[1].set_xlabel("Time $t$ [s]") + +fbs.savefig('figure-6.14-cruise_linearized.png') diff --git a/examples/figure-6.5-modes.py b/examples/figure-6.5-modes.py new file mode 100644 index 000000000..82bc0e1c6 --- /dev/null +++ b/examples/figure-6.5-modes.py @@ -0,0 +1,61 @@ +# modes.py - illustration of modes for a second order system +# RMM, 19 Apr 2024 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# System definition +k0, k1, k2, b0 = 0.1, 0.1, 0.5, 1.5 +A = [[-k0 - k1, k1], [k2, -k2]] +B = [[b0], [0]] +C = [[1, 0], [0, 1]] +sys = ct.ss(A, B, C, 0) + +# Generate a phase plot for the system +fig, ax = plt.subplots(1, 1, figsize=[3.4, 3.4]) +ct.phase_plane_plot(sys, [-1, 1, -1, 1], 5, gridspec=[7, 4]) +ax.set_aspect('equal') + +# Label the figure +plt.suptitle("") +plt.text(-1, -0.8, "Slow") +plt.text(-0.2, 0.8, "Fast") +fbs.savefig('figure-6.5-modes-pp.png') + +# Time domain simulations +evals, evecs = np.linalg.eig(A) +X0s, X0f = evecs[:, 0], evecs[:, 1] # Fast and slow modes +tvec = np.linspace(0, 50, endpoint=True) + +# Slow mode +ax = fbs.figure('321') +resp_slow = ct.initial_response(sys, tvec, X0=X0s) + +ax.set_xlim([0, 50]) +ax.set_ylim([0, 1]) +ax.plot(tvec, resp_slow.states[0], 'b', label="$x_1$") +ax.plot(tvec, resp_slow.states[1], 'r--', label="$x_2$") +ax.text(10, 0.75, "Slow mode") +ax.set_xlabel("Time $t$ [s]") +ax.set_ylabel("$x_1, x_2$") +ax.legend(frameon=False) + +plt.tight_layout() +fbs.savefig('figure-6.5-modes-slow.png') + +# Fast mode +ax = fbs.figure('321') +resp_fast = ct.initial_response(sys, tvec, X0=X0f) + +ax.set_xlim([0, 50]) +ax.set_ylim([-0.25, 1]) +ax.plot(tvec, resp_fast.states[0], 'b', label="$x_1$") +ax.plot(tvec, resp_fast.states[1], 'r--', label="$x_2$") +ax.text(10, 0.75, "Fast mode") +ax.set_xlabel("Time $t$ [s]") +ax.set_ylabel("$x_1, x_2$") + +plt.tight_layout() +fbs.savefig('figure-6.5-modes-fast.png') diff --git a/examples/figure-7.7-predprey_place.py b/examples/figure-7.7-predprey_place.py new file mode 100644 index 000000000..4923d561b --- /dev/null +++ b/examples/figure-7.7-predprey_place.py @@ -0,0 +1,61 @@ +# predprey_place.py - stabilization via state space feedback +# RMM, 18 Jan 2026 + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# System definition +from predprey import predprey + +# Find the equilibrium point and linearize +xe, ue = ct.find_eqpt(predprey, [20, 30], 0) +sys = predprey.linearize(xe, ue) + +# Eigenvalue placement +K = ct.place(sys.A, sys.B, [-0.1, -0.2]) +kf = -1 / (sys.C[1] @ np.linalg.inv(sys.A - sys.B @ K) @ sys.B) + +ctrl = ct.nlsys( + None, lambda t, x, u, params: -K @ (u[0:2] - xe) + kf * (u[2] - xe[1]), + inputs=['H', 'L', 'r'], outputs=['u'], +) +clsys = ct.interconnect( + [predprey, ctrl], inputs=['r'], outputs=['H', 'L', 'u'], + name='predprey_clsys' +) + +# Compute initial condition response +T = np.linspace(0, 100, 1000) +response = ct.input_output_response(clsys, T, 30, [20, 15]) + +# Plot results +fbs.figure('mlh') # Create figure using FBS defaults +plt.plot(response.time, response.outputs[0], 'b-', label="Hare") +plt.plot(response.time, response.outputs[1], 'r--', label="Lynx") + +plt.xlabel("Time [years]") +plt.ylabel("Population") +plt.legend(frameon=False) + +plt.tight_layout() +fbs.savefig('figure-7.7-predprey_place-time.png') + +# Generate a phase portrait +fbs.figure('mlh') # Create figure using FBS defaults + +# Create a function for the phase portrait that includes the reference input +ppfcn = lambda t, x: clsys.dynamics(t, x, [30], {}) + +ct.phase_plane_plot( + ppfcn, [1, 100, 1, 100], 30, plot_separatrices=False) +# Turn off vector field lines since phase_plan_plot has arrows +# ct.phaseplot.vectorfield(ppfcn, [0, 100, 0, 100], gridspec=[20, 20]) + +plt.xlabel("Hares") +plt.ylabel("Lynxes") +plt.suptitle("") + +plt.tight_layout() +fbs.savefig('figure-7.7-predprey_place-pp.png') diff --git a/examples/predprey.py b/examples/predprey.py new file mode 100644 index 000000000..bf92cb197 --- /dev/null +++ b/examples/predprey.py @@ -0,0 +1,30 @@ +# predprey.py - predator-prey dynamics +# RMM, 21 Apr 2024 +# +# Predator-prey dynamics +# +# This model implements the dynamics of a predator-prey system, as +# described in Section 2.7 of FBS. + +import numpy as np +import control as ct + +# Define the dynamis for the predator-prey system (no input) +predprey_params = {'r': 1.6, 'd': 0.56, 'b': 0.6, 'k': 125, 'a': 3.2, 'c': 50} +def predprey_update(t, x, u, params): + """Predator prey dynamics""" + r, d, b, k, a, c = map(params.get, ['r', 'd', 'b', 'k', 'a', 'c']) + u = np.atleast_1d(u) # Fix python-control bug + u = np.clip(u, 0, 4*r) # constraints used in FBS 2e + # u = np.clip(u, -r, 3*r) # constrain the input to keep r_eff >= 0 + + # Dynamics for the system + dH = (r + u[0]) * x[0] * (1 - x[0]/k) - a * x[1] * x[0]/(c + x[0]) + dL = b * a * x[1] * x[0] / (c + x[0]) - d * x[1] + + return np.array([dH, dL]) + +# Create a nonlinear I/O system +predprey = ct.nlsys( + predprey_update, name='predprey', params=predprey_params, + states=['H', 'L'], inputs='u', outputs=['H', 'L']) diff --git a/examples/springmass.py b/examples/springmass.py new file mode 100644 index 000000000..44aea7251 --- /dev/null +++ b/examples/springmass.py @@ -0,0 +1,11 @@ +# springmass.py - Spring mass dynamics +# RMM, 19 Apr 2024 + +import control as ct + +m, k, b = 250, 40, 60 +A = [[0, 1], [-k/m, -b/m]] +B = [[0], [1/m]] +C = [[1, 0]] + +springmass = ct.ss(A, B, C, 0) diff --git a/examples/steering.py b/examples/steering.py new file mode 100644 index 000000000..7c60aecab --- /dev/null +++ b/examples/steering.py @@ -0,0 +1,68 @@ +import numpy as np +import control as ct + +# Define the nonlinear dynamics for the vehicle steering system +def steering_update(t, x, u, params): + """ + Nonlinear dynamics of the vehicle steering system. + + Parameters: + t : float + Current time (not used in the dynamics, but included for compatibility). + x : array_like + Current state vector: [x1, x2, x3]. + u : array_like + Input vector: [velocity, steering angle]. + params : dict + Dictionary of parameters (a, b, maxsteer). + + Returns: + dx : array_like + State derivatives: [dx1, dx2, dx3]. + """ + a = params.get('a', 1.5) # Offset to vehicle reference point [m] + b = params.get('b', 3.0) # Vehicle wheelbase [m] + maxsteer = params.get('maxsteer', 0.5) # Max steering angle [rad] + + # Saturate the steering input + delta = np.clip(u[1], -maxsteer, maxsteer) + + # System dynamics + alpha = np.arctan2(a * np.tan(delta), b) + dx1 = u[0] * np.cos(x[2] + alpha) # x velocity + dx2 = u[0] * np.sin(x[2] + alpha) # y velocity + dx3 = (u[0] / b) * np.tan(delta) # Angular velocity + return [dx1, dx2, dx3] + +# Create a nonlinear input/output system +steering = ct.nlsys( + steering_update, # Update function for system dynamics + inputs=['v', 'delta'], # Inputs: velocity and steering angle + outputs=None, # Outputs are the same as the states + states=['x', 'y', 'theta'], # States: x, y, and theta (angle) + name='steering', + params={'a': 1.5, 'b': 3.0, 'maxsteer': 0.5} # Default parameters +) + +# Generate the linearization at a given velocity +def linearize_lateral(v0=10, normalize=False, output_full_state=False): + # Compute the linearization at the given velocity + linsys = ct.linearize(steering, 0, [v0, 0]) + + # Extract out the lateral dynamics + latsys = ct.model_reduction( + linsys, elim_states=[0], keep_inputs=[1], keep_outputs=[1], + method='truncate', warn_unstable=False) + + # Normalize coordinates if desired + if normalize: + b = steering.params['b'] + latsys = ct.similarity_transform( + latsys, np.array([[1/b, 0], [0, 1]]), timescale=v0/b) + + C = np.eye(2) if output_full_state else np.array([[1, 0]]) + + # Normalized system with (normalized) lateral offset as output + return ct.ss( + latsys.A, latsys.B, C, 0, inputs='delta', + outputs=['y', 'theta'] if output_full_state else 'y') diff --git a/examples/template.py b/examples/template.py new file mode 100644 index 000000000..9dc112355 --- /dev/null +++ b/examples/template.py @@ -0,0 +1,17 @@ +# figure_name.py - short description +# RMM, date + +import matplotlib.pyplot as plt +import numpy as np +import control as ct +import fbs # FBS plotting customizations + +# System definition + +# Compute system response + +# Plot results +fbs.figure('mlh') # Create figure using FBS defaults + +plt.tight_layout() +fbs.savefig('figure-N.m-figure_name-panel.png')