FBSwiki - User contributions [en]

Figure 7.7: Simulation results for the controlled predator–prey system

2026-01-19T01:51:51Z

Murray:

{{Figure
|Chapter=State Feedback
|Figure number=7
|Sort key=707
|Figure title=Simulation results for the controlled predator–prey system
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/figure-7.7-predprey_place.py
|Requires=predprey.py
}}
[[Image:figure-7.7-predprey_place-time.png]]  
[[Image:figure-7.7-predprey_place-pp.png]]

'''Figure 7.7:''' Simulation results for the controlled predator–prey system. The population of lynxes and hares as a function of time is shown in (a), and a phase portrait for the controlled system is shown in (b). Feedback is used to make the population stable at He = 20.6 and Le = 30.

<nowiki>
# 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')
</nowiki>

File:Figure-7.7-predprey place-pp.png

2026-01-18T16:37:46Z

Murray:

File:Figure-7.7-predprey place-time.png

2026-01-18T16:33:03Z

Murray:

Figure 7.7: Simulation results for the controlled predator–prey system

2026-01-18T16:32:24Z

Murray:

{{Figure
|Chapter=7
|Figure number=7
|Sort key=707
|Figure title=Simulation results for the controlled predator–prey system
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/figure-7.7-predprey_place.py
|Requires=predprey.py
}}
[[Image:figure-7.7-predprey_place-time.png]]  
[[Image:figure-7.7-predprey_place-pp.png]]

'''Figure 7.7:''' Simulation results for the controlled predator–prey system. The population of lynxes and hares as a function of time is shown in (a), and a phase portrait for the controlled system is shown in (b). Feedback is used to make the population stable at He = 20.6 and Le = 30.

<nowiki>
# 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')
</nowiki>

Figure 7.7: Simulation results for the controlled predator–prey system

2026-01-18T16:27:33Z

Murray: Created page with "{{Figure |Chapter=7 |Figure number=7 |Sort key=707 |Figure title=Simulation results for the controlled predator–prey system |GitHub URL=https://github.com/murrayrm/fbs2e-pyt..."

Biomolecular Feedback Systems/index.php/Main Page

2025-07-03T12:04:04Z

Murray: Redirected page to Biomolecular Feedback Systems

#REDIRECT [[Biomolecular Feedback Systems]]

Biomolecular Feedback Systems/

2025-07-03T12:02:27Z

Murray: Redirected page to Biomolecular Feedback Systems

#REDIRECT [[Biomolecular Feedback Systems]]

Biomolecular Feedback Systems

2025-06-24T13:24:51Z

Murray: /* Contents */

{{BFS quick links}}
[[Domitilla Del Vecchio]] and [[User:murray|Richard M. Murray]]
__NOTOC__
This page contains information on the text ''Biomolecular Feedback Systems'' by Domitilla Del Vecchio and Richard M. Murray.

Copyright in this book is held by Princeton University Press, who have kindly agreed to allow us to keep the book available on the web.

This textbook is intended for researchers interested in the application of feedback and control to biomolecular systems. The material has been designed so that it can be used in parallel with ''Feedback Systems'' as part of a course on biomolecular feedback and control systems, or as a standalone reference for readers who have had a basic course in feedback and control theory.

===== News =====
{{#ask: [[Category:Announcements]] [[Book::BFS]]
| sort=Date | order=descening
| format=ul | limit=6
}}

=== Contents ===
{| width=100% border=1
|- valign=top
| width=50% |
* '''{{BFS pdf|Preface and Contents|bfs-frontmatter|14Sep14}}'''
* '''Chapter 1: {{BFS pdf|Introductory Concepts|bfs-intro|14Sep14}}'''
** Systems Biology: Modeling, Analysis and the Role of Feedback
** The Cell as a System
** Control and Dynamical Systems Tools
** Input/Output Modeling
** From Systems to Synthetic Biology
* '''Chapter 2: {{BFS pdf|Core Processes|bfs-coreproc|14Sep14}}'''
** Modeling Chemical Reactions
** Transcription and Translation
** Transcriptional Regulation
** Post-Transcriptional Regulation
** Cellular Subsystems
* '''Chapter 3: {{BFS pdf|Dynamic Behavior|bfs-dynamics|14Sep14}}'''
** Analysis Near Equilibria
** Robustness
** Oscillatory Behavior
** Bifurcations
** Model Reduction Techniques
* '''Chapter 4: {{BFS pdf|Stochastic Modeling and Analysis|bfs-stochastic|14Sep14}}'''
** Stochastic Modeling of Biochemical Systems
** Simulation of Stochastic Systems
** Input/Output Linear Stochastic Systems
| width=50% |
* '''Chapter 5: {{BFS pdf|Biological Circuit Components|bfs-circuits|14Sep14}}'''
** Introduction to Biological Circuit Design
** Negative Autoregulation
** The Toggle Switch
** The Repressilator
** Activator-repressor Clock
** An Incoherent Feedforward Loop (IFFL)
** Bacterial Chemotaxis
* '''Chapter 6: {{BFS pdf|Interconnecting Components|bfs-modules|14Sep14}}'''
** Input/Output Modeling and the Modularity Assumption
** Introduction to Retroactivity
** Retroactivity in Gene Circuits
** Retroactivity in Signaling Systems
** Insulation Devices: Retroactivity Attentuation
** A Case Study on the Use of Insulation Devices
* '''Chapter 7: {{BFS pdf|Design Tradeoffs|bfs-tradeoffs|14Sep14}}'''
** Competition for Shared Cellular Resources
** Stochastic Effects: Design Tradeoffs in Systems with Large Gains
* '''Chapter 8: {{BFS pdf|Resource Competition|bfs-competition|24Jun2025}}'''
** Resource competition in bacterial genetic circuits
** Resource competition in mammalian genetic circuits
** Feedforward control to mitigate the effects of resource sharing
** Feedback control to mitigate the effects of resource sharing
* '''{{BFS pdf|Bibliography and Index|bfs-backmatter|14Sep14}}'''
|}

=== Resources ===
* [[BFS errata|List of errata]]

Figure 6.13: AFM frequency response

2025-02-19T06:09:09Z

Murray:

{{Figure
|Chapter=Linear Systems
|Figure number=6.13
|Sort key=613
|Figure title=AFM frequency response
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-6.9-afm_freqresp.py
}}
[[Image:figure-6.13-afm_freqresp.png]]

'''Figure 6.13:''' AFM frequency response. (a) A block diagram for the vertical dynamics of an atomic force microscope in contact mode [not shown]. 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 <math>\omega = 268</math> 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 <math>\sqrt{2}</math>.

<nowiki>
# example-6.9-afm_freqresp.py - <Short description>
# 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')
</nowiki>

File:Springmass-coupled.png

2024-12-26T17:05:30Z

Murray:

File:Figure-7.6-steering place.png

2024-11-29T01:45:21Z

Murray:

Figure 7.6: State feedback control of a steering system

2024-11-29T01:44:24Z

Murray:

{{Figure
|Chapter=State Feedback
|Figure number=7.6
|Sort key=706
|Figure title=State feedback control of a steering system
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-7.4-steering_place.py
}}
[[Image:figure-7.6-steering_place.png]]

'''Figure 7.6:''' State feedback control of a steering system. Unit step responses (from zero initial condition) obtained with controllers designed with <math>\zeta_c = 0.7</math> and <math>\omega_c = 0.5</math>, <math>0.7</math>, and <math>1</math> [rad/s] are shown in (a). The dashed lines indicate <math>\pm 5</math>% deviations from the setpoint. Notice that response speed increases with increasing <math>\omega_c</math>, but that large <math>\omega_c</math> also give large initial control actions. Unit step responses obtained with a controller designed with <math>\omega_c = 0.7</math> and <math>\zeta_c = 0.5</math>, <math>0.7</math>, and <math>1</math> are shown in (b).

<nowiki>
# example-7.4-steering_place.py - <Short description>
# 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
clsys = ct.feedback(sys, K) * kf

# Simulate the closed loop dynamics
response = ct.forced_response(clsys, timepts, 1, X0=0)

ax_pos.plot(response.time, response.states[0], 'b')
ax_delta.plot(response.time, (kf - K @ response.states)[0], '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
clsys = ct.feedback(sys, K) * kf

# Simulate the closed loop dynamics
response = ct.forced_response(clsys, timepts, 1, X0=0)

ax_pos.plot(response.time, response.states[0], 'b')
ax_delta.plot(response.time, (kf - K @ response.states)[0], '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.4-steering_place.png", bbox_inches='tight')
</nowiki>

Figure 7.6: State feedback control of a steering system

2024-11-29T01:42:30Z

Murray: Created page with "{{Figure |Chapter=State Feedback |Figure number=7.6 |Sort key=706 |Figure title=State feedback control of a steering system |GitHub URL=https://github.com/murrayrm/fbs2e-pytho..."

{{Figure
|Chapter=State Feedback
|Figure number=7.6
|Sort key=706
|Figure title=State feedback control of a steering system
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-7.4-steering_place.py
}}
[[Image:figure-7.6-steering_place.png]]

'''Figure 7.6:''' State feedback control of a steering system. Unit step responses (from zero initial condition) obtained with controllers designed with <math>\zeta_c = 0.7</math) and <math>\omega_c = 0.5</math>, 0.7, and 1 [rad/s] are shown in (a). The dashed lines indicate <math>\pm 5</math>% deviations from the setpoint. Notice that response speed increases with increasing <math>\omega_c</math>, but that large <math>\omega_c</math> also give large initial control actions. Unit step responses obtained with a controller designed with <math>\omega_c = 0.7</math> and <math>\zeta_c = 0.5</math>, 0.7, and 1 are shown in (b).

<nowiki>
# example-7.4-steering_place.py - <Short description>
# 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
clsys = ct.feedback(sys, K) * kf

# Simulate the closed loop dynamics
response = ct.forced_response(clsys, timepts, 1, X0=0)

ax_pos.plot(response.time, response.states[0], 'b')
ax_delta.plot(response.time, (kf - K @ response.states)[0], '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
clsys = ct.feedback(sys, K) * kf

# Simulate the closed loop dynamics
response = ct.forced_response(clsys, timepts, 1, X0=0)

ax_pos.plot(response.time, response.states[0], 'b')
ax_delta.plot(response.time, (kf - K @ response.states)[0], '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.4-steering_place.png", bbox_inches='tight')
</math>

File:Figure-6.13-afm freqresp.png

2024-11-28T19:22:02Z

Murray:

Figure 6.13: AFM frequency response

2024-11-28T19:21:38Z

Murray: Created page with "{{Figure |Chapter=Linear Systems |Figure number=6.13 |Sort key=613 |Figure title=AFM frequency response |GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-..."

{{Figure
|Chapter=Linear Systems
|Figure number=6.13
|Sort key=613
|Figure title=AFM frequency response
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-6.9-afm_freqresp.py
}}
[[Image:figure-6.13-afm_freqresp.png]]

'''Figure 6.13:''' AFM frequency response. (a) A block diagram for the vertical dynamics of an atomic force microscope in contact mode [not shown]. 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 <math>\omega = 268</math> 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 <math>\sqrt{2}</math>.

<nowiki>
# example-6.9-afm_freqresp.py - <Short description>
# 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[0, 0], freqresp.phase[0, 0]

# 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')
</nowiki>

File:Figure-6.12-opamp bandpass.png

2024-11-28T18:49:39Z

Murray:

Figure 6.12: Active band-pass filter

2024-11-28T18:49:06Z

Murray: Created page with "{{Figure |Chapter=Linear Systems |Figure number=6.12 |Sort key=612 |Figure title=Active band-pass filter |GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example..."

{{Figure
|Chapter=Linear Systems
|Figure number=6.12
|Sort key=612
|Figure title=Active band-pass filter
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-6.8-opamp_bandpass.py
}}
[[Image:figure-6.12-opamp_bandpass.png]]

'''Figure 6.12:''' Active band-pass filter. The circuit diagram (a) [not shown] 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 <math>-90^\circ</math> due to the negative gain of the operational amplifier.

<nowiki>
# example-6.8-opamp_bandpass.py - <Short description>
# 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')
</nowiki>

Errata: sign errors in Example 5.18 (noise cancellation)

2024-11-27T00:44:16Z

Murray:

{{Errata
|Chapter=Dynamic Behavior
|Page number=5-34
|Line=8
|Version=3.1.5
|Date=26 Nov 2024
}}
In Example 5.18 (noise cancellation), there are two sign errors in equation (5.26) that are propagated through the next several lines. The corrected text should read (with changes in red):

Assuming for simplicity that <math>S=0</math>, introduce <math>x_1=e=z-w</math>, <math>x_2 = a - a_0</math>, and <math>x_3 = b - b_0</math>. Then
{| width=100%
|-
|
<center><math>
\frac{dx_1}{dt} = a_0 (z - w) + (a - a_0) w + (b - b_0) n
= a_0 x_1 {\color{red} \boldsymbol{-}} x_2 w {\color{red} \boldsymbol{-}} x_3 n.
</math></center>
| (5.26)
|}
We will achieve noise cancellation if we can find a feedback law for changing the parameters <math>a</math> and <math>b</math> so that the error <math>e</math> goes to zero. To do this we choose
<center><math>
V(x_1,x_2,x_3)=\frac{1}{2} \bigl( \alpha x_1^2+x_2^2+x_3^2 \bigr)
</math></center>
as a candidate Lyapunov function for
equation (5.26). The derivative of <math>V</math> is
<center><math>
\dot V = \alpha x_1 \dot x_1 + x_2 \dot x_2 + x_3 \dot x_3
= \alpha a_0 x_1^2 + x_2 (\dot x_2 {\color{red} \boldsymbol{-}} \alpha w x_1) + x_3 (\dot x_3 {\color{red} \boldsymbol{-}}
\alpha n x_1).
</math></center>
Choosing
{| width=100%
|-
|
<center><math>
\dot a = \dot x_2 = {\color{red}+} \alpha w x_1 = {\color{red}+} \alpha w e,\qquad
\dot b =\dot x_3 = {\color{red}+} \alpha n x_1 = {\color{red}+} \alpha n e,
</math></center>
| (5.27)
|}
we find that <math>\dot V = \alpha a_0 x_1^2 < 0</math>, and it follows that the
quadratic function will decrease as long as <math>e = x_1 = w - z \neq 0</math>.

Errata: sign errors in Example 5.18 (noise cancellation)

2024-11-27T00:43:17Z

Murray:

{{Errata
|Chapter=Dynamic Behavior
|Page number=5-34
|Line=8
|Version=3.1.5
|Date=26 Nov 2024
}}
In Example 5.18 (noise cancellation), there are two sign errors in equation (5.26) that are propagated through the next several lines. The corrected text should read (with changes in red):

Assuming for simplicity that <math>S=0</math>, introduce <math>x_1=e=z-w</math>, <math>x_2 = a - a_0</math>, and <math>x_3 = b - b_0</math>. Then
{| width=100%
|-
|
<center><math>
\frac{dx_1}{dt} = a_0 (z - w) + (a - a_0) w + (b - b_0) n
= a_0 x_1 {\color{red} \boldsymbol{-}} x_2 w {\color{blue}-} x_3 n.
</math></center>
| (5.26)
|}
We will achieve noise cancellation if we can find a feedback law for changing the parameters <math>a</math> and <math>b</math> so that the error <math>e</math> goes to zero. To do this we choose
<center><math>
V(x_1,x_2,x_3)=\frac{1}{2} \bigl( \alpha x_1^2+x_2^2+x_3^2 \bigr)
</math></center>
as a candidate Lyapunov function for
equation (5.26). The derivative of <math>V</math> is
<center><math>
\dot V = \alpha x_1 \dot x_1 + x_2 \dot x_2 + x_3 \dot x_3
= \alpha a_0 x_1^2 + x_2 (\dot x_2 {\color{red}-} \alpha w x_1) + x_3 (\dot x_3 {\color{red}-}
\alpha n x_1).
</math></center>
Choosing
{| width=100%
|-
|
<center><math>
\dot a = \dot x_2 = {\color{red}+} \alpha w x_1 = {\color{red}+} \alpha w e,\qquad
\dot b =\dot x_3 = {\color{red}+} \alpha n x_1 = {\color{red}+} \alpha n e,
</math></center>
| (5.27)
|}
we find that <math>\dot V = \alpha a_0 x_1^2 < 0</math>, and it follows that the
quadratic function will decrease as long as <math>e = x_1 = w - z \neq 0</math>.

Errata: sign errors in Example 5.18 (noise cancellation)

2024-11-27T00:40:38Z

Murray:

{{Errata
|Chapter=Dynamic Behavior
|Page number=5-34
|Line=8
|Version=3.1.5
|Date=26 Nov 2024
}}
In Example 5.18 (noise cancellation), there are two sign errors in equation (5.26) that are propagated through the next several lines. The corrected text should read (with changes in red):

Assuming for simplicity that <math>S=0</math>, introduce <math>x_1=e=z-w</math>, <math>x_2 = a - a_0</math>, and <math>x_3 = b - b_0</math>. Then
{| width=100%
|-
|
<center><math>
\frac{dx_1}{dt} = a_0 (z - w) + (a - a_0) w + (b - b_0) n
= a_0 x_1 {\color{red} \boldsymbol{-}} x_2 w {\color{red}-} x_3 n.
</math></center>
| (5.26)
|}
We will achieve noise cancellation if we can find a feedback law for changing the parameters <math>a</math> and <math>b</math> so that the error <math>e</math> goes to zero. To do this we choose
<center><math>
V(x_1,x_2,x_3)=\frac{1}{2} \bigl( \alpha x_1^2+x_2^2+x_3^2 \bigr)
</math></center>
as a candidate Lyapunov function for
equation (5.26). The derivative of <math>V</math> is
<center><math>
\dot V = \alpha x_1 \dot x_1 + x_2 \dot x_2 + x_3 \dot x_3
= \alpha a_0 x_1^2 + x_2 (\dot x_2 {\color{red}-} \alpha w x_1) + x_3 (\dot x_3 {\color{red}-}
\alpha n x_1).
</math></center>
Choosing
{| width=100%
|-
|
<center><math>
\dot a = \dot x_2 = {\color{red}+} \alpha w x_1 = {\color{red}+} \alpha w e,\qquad
\dot b =\dot x_3 = {\color{red}+} \alpha n x_1 = {\color{red}+} \alpha n e,
</math></center>
| (5.27)
|}
we find that <math>\dot V = \alpha a_0 x_1^2 < 0</math>, and it follows that the
quadratic function will decrease as long as <math>e = x_1 = w - z \neq 0</math>.

Errata: sign errors in Example 5.18 (noise cancellation)

2024-11-27T00:39:11Z

Murray:

{{Errata
|Chapter=Dynamic Behavior
|Page number=5-34
|Line=8
|Version=3.1.5
|Date=26 Nov 2024
}}
In Example 5.18 (noise cancellation), there are two sign errors in equation (5.26) that are propagated through the next several lines. The corrected text should read (with changes in red):

Assuming for simplicity that <math>S=0</math>, introduce <math>x_1=e=z-w</math>, <math>x_2 = a - a_0</math>, and <math>x_3 = b - b_0</math>. Then
{| width=100%
|-
|
<center><math>
\frac{dx_1}{dt} = a_0 (z - w) + (a - a_0) w + (b - b_0) n
= a_0 x_1 {\color{red}-} x_2 w {\color{red}-} x_3 n.
</math></center>
| (5.26)
|}
We will achieve noise cancellation if we can find a feedback law for changing the parameters <math>a</math> and <math>b</math> so that the error <math>e</math> goes to zero. To do this we choose
<center><math>
V(x_1,x_2,x_3)=\frac{1}{2} \bigl( \alpha x_1^2+x_2^2+x_3^2 \bigr)
</math></center>
as a candidate Lyapunov function for
equation (5.26). The derivative of <math>V</math> is
<center><math>
\dot V = \alpha x_1 \dot x_1 + x_2 \dot x_2 + x_3 \dot x_3
= \alpha a_0 x_1^2 + x_2 (\dot x_2 {\color{red}-} \alpha w x_1) + x_3 (\dot x_3 {\color{red}-}
\alpha n x_1).
</math></center>
Choosing
{| width=100%
|-
|
<center><math>
\dot a = \dot x_2 = {\color{red}+} \alpha w x_1 = {\color{red}+} \alpha w e,\qquad
\dot b =\dot x_3 = {\color{red}+} \alpha n x_1 = {\color{red}+} \alpha n e,
</math></center>
| (5.27)
|}
we find that <math>\dot V = \alpha a_0 x_1^2 < 0</math>, and it follows that the
quadratic function will decrease as long as <math>e = x_1 = w - z \neq 0</math>.

Errata: sign errors in Example 5.18 (noise cancellation)

2024-11-27T00:30:05Z

Murray:

{{Errata
|Chapter=Dynamic Behavior
|Page number=5-34
|Line=8
|Version=3.1.5
|Date=26 Nov 2024
}}
In Example 5.18 (noise cancellation), there are two sign errors in equation (5.26) that are propagated through the next several lines. The corrected text should read (with changes in red):

Assuming for simplicity that <math>S=0</math>, introduce <math>x_1=e=z-w</math>, <math>x_2 = a - a_0</math>, and <math>x_3 = b - b_0</math>. Then
{| width=100%
|-
|
<center><math>
\frac{dx_1}{dt} = a_0 (z - w) + (a - a_0) w + (b - b_0) n
= a_0 x_1 - x_2 w - x_3 n.
</math></center>
| (5.26)
|}
We will achieve noise cancellation if we can find a feedback law for changing the parameters <math>a</math> and <math>b</math> so that the error <math>e</math> goes to zero. To do this we choose
<center><math>
V(x_1,x_2,x_3)=\frac{1}{2} \bigl( \alpha x_1^2+x_2^2+x_3^2 \bigr)
</math></center>
as a candidate Lyapunov function for
equation (5.26). The derivative of <math>V</math> is
<center><math>
\dot V = \alpha x_1 \dot x_1 + x_2 \dot x_2 + x_3 \dot x_3
= \alpha a_0 x_1^2 + x_2 (\dot x_2 - \alpha w x_1) + x_3 (\dot x_3 -
\alpha n x_1).
</math></center>
Choosing
{| width=100%
|-
|
<center><math>
\dot a = \dot x_2 = \alpha w x_1 = \alpha w e,\qquad
\dot b =\dot x_3 = \alpha n x_1 = \alpha n e,
</math></center>
| (5.27)
|}
we find that <math>\dot V = \alpha a_0 x_1^2 < 0</math>, and it follows that the
quadratic function will decrease as long as <math>e = x_1 = w - z \neq 0</math>.

Errata: sign errors in Example 5.18 (noise cancellation)

2024-11-27T00:29:24Z

Murray:

{{Errata
|Chapter=Dynamic Behavior
|Page number=5-34
|Line=8
|Version=3.1.5
|Date=26 Nov 2024
}}
In Example 5.18 (noise cancellation), there are two sign errors in equation (5.26) that are propagated through the next several lines. The corrected text should read (with changes in red):

Assuming for simplicity that <math>S=0</math>, introduce <math>x_1=e=z-w</math>, <math>x_2 = a - a_0</math>, and <math>x_3 = b - b_0</math>. Then
{| width=100%
|-
|
<center><math>
\frac{dx_1}{dt} = a_0 (z - w) + (a - a_0) w + (b - b_0) n
= a_0 x_1 - x_2 w - x_3 n.
</math></center>
| (5.26)
|}
We will achieve noise cancellation if we can find a feedback law for changing the parameters <math>a</math> and <math>b</math> so that the error <math>e</math> goes to zero. To do this we choose
<center><math>
V(x_1,x_2,x_3)=\frac{1}{2} \bigl( \alpha x_1^2+x_2^2+x_3^2 \bigr)
</math></center>
as a candidate Lyapunov function for
equation (5.26). The derivative of <math>V</math> is
<center><math>
\dot V = \alpha x_1 \dot x_1 + x_2 \dot x_2 + x_3 \dot x_3
= \alpha a_0 x_1^2 + x_2 (\dot x_2 - \alpha w x_1) + x_3 (\dot x_3 -
\alpha n x_1).
</math></center>
Choosing
<center><math>
\dot a = \dot x_2 = \alpha w x_1 = \alpha w e,\qquad
\dot b =\dot x_3 = \alpha n x_1 = \alpha n e,
</math></center>
we find that <math>\dot V = \alpha a_0 x_1^2 < 0</math>, and it follows that the
quadratic function will decrease as long as <math>e = x_1 = w - z \neq 0</math>.

Errata: sign errors in Example 5.18 (noise cancellation)

2024-11-27T00:27:57Z

Murray: Created page with "{{Errata |Chapter=Dynamic Behavior |Page number=5-34 |Line=8 |Version=3.1.5 |Date=26 Nov 2024 }} In Example 5.18 (noise cancellation), there are two sign errors in equation (5..."

{{Errata
|Chapter=Dynamic Behavior
|Page number=5-34
|Line=8
|Version=3.1.5
|Date=26 Nov 2024
}}
In Example 5.18 (noise cancellation), there are two sign errors in equation (5.26) that are propagated through the next several lines. The corrected text should read (with changes in red):

Assuming for simplicity that <math>S=0</math>, introduce <math>x_1=e=z-w</math>, <math>x_2 = a - a_0</math>, and <math>x_3 = b - b_0</math>. Then
<center><math>
\hfill \frac{dx_1}{dt} = a_0 (z - w) + (a - a_0) w + (b - b_0) n
= a_0 x_1 - x_2 w - x_3 n. \hfill (5.26)
</math></center>
We will achieve noise cancellation if we can find a feedback law for changing the parameters <math>a</math> and <math>b</math> so that the error <math>e</math> goes to zero. To do this we choose
<center><math>
V(x_1,x_2,x_3)=\frac{1}{2} \bigl( \alpha x_1^2+x_2^2+x_3^2 \bigr)
</math></center>
as a candidate Lyapunov function for
equation (5.26). The derivative of <math>V</math> is
<center><math>
\dot V = \alpha x_1 \dot x_1 + x_2 \dot x_2 + x_3 \dot x_3
= \alpha a_0 x_1^2 + x_2 (\dot x_2 - \alpha w x_1) + x_3 (\dot x_3 -
\alpha n x_1).
</math></center>
Choosing
<center><math>
\dot a = \dot x_2 = \alpha w x_1 = \alpha w e,\qquad
\dot b =\dot x_3 = \alpha n x_1 = \alpha n e,
\label{ex:dynamics:noisecancu}
</math></center>
we find that <math>\dot V = \alpha a_0 x_1^2 < 0</math>, and it follows that the
quadratic function will decrease as long as <math>e = x_1 = w - z \neq 0</math>.

File:Figure-6.10-compartment response.png

2024-11-25T06:29:37Z

Murray:

Figure 6.10: Response of a compartment model to a constant drug infusion

2024-11-25T06:29:20Z

Murray: Created page with "{{Figure |Chapter=Linear Systems |Figure number=6.10 |Sort key=610 |Figure title=Response of a compartment model to a constant drug infusion |GitHub URL=https://github.com/mur..."

{{Figure
|Chapter=Linear Systems
|Figure number=6.10
|Sort key=610
|Figure title=Response of a compartment model to a constant drug infusion
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-6.10-compartment_response.py
}}
[[Image:figure-6.10-compartment_response.png]]

'''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.

<nowiki>
# 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')
</nowiki>

Spring-mass system

2024-11-25T05:59:57Z

Murray: Created page with "{| class="wikitable" |- ! GitHub URL | https://github.com/murrayrm/fbs2e-python/blob/main/springmass.py |- ! Requires | python-control |} This pa..."

{| class="wikitable"
|-
! GitHub URL
| https://github.com/murrayrm/fbs2e-python/blob/main/springmass.py
|-
! Requires
| [[https:python-control.org|python-control]]
|}

This page documents the spring-mass system that is used as a running example throughout the text. A description of the dynamics of this system are presented in {{chapter link|System Modeling}}, Example 3.1. This page contains a description of the system, including the models and commands used to generate some of the plots in the text.

Note: the Python code on this page makes use of the [https://python-control.org Python Control Systems Library], an open source software toolbox for control systems analysis.

Figures making use of this model:
{{#ask: [[Category:Figures]] [[Requires::springmass.py]] | format=ul | sort=Sort key}}

Springmass.py

2024-11-25T05:58:19Z

Murray: Redirected page to Spring-mass system

#REDIRECT [[Spring-mass system]]

Figure 6.5: Modes for a second-order system with real eigenvalues

2024-11-25T05:55:44Z

Murray:

{{Figure
|Chapter=Linear Systems
|Figure number=6.5
|Sort key=605
|Figure title=The notion of modes for a second-order system with real eigenvalues.
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/figure-6.5-modes.py
|Requires=springmass.py
}}
{| border=0
|-
| rowspan=2 | [[Image:figure-6.5-springmass_modes-pp.png]]
| [[Image:figure-6.5-springmass_modes-slow.png]]
|-
| [[Image:figure-6.5-springmass_modes-fast.png]]
|}

'''Figure 6.5''': The notion of modes for a second-order system with real eigenvalues. The left figure shows the phase portrait and the modes corresponding to solutions that start on the eigenvectors (bold lines). The corresponding time functions are shown on the right.

<nowiki>
# 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')
</nowiki>

Figure 6.5: Modes for a second-order system with real eigenvalues

2024-11-25T05:55:00Z

Murray:

{{Figure
|Chapter=Linear Systems
|Figure number=6.5
|Sort key=605
|Figure title=The notion of modes for a second-order system with real eigenvalues.
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/figure-6.5-superposition.py
|Requires=springmass.py
}}
{| border=0
|-
| rowspan=2 | [[Image:figure-6.5-springmass_modes-pp.png]]
| [[Image:figure-6.5-springmass_modes-slow.png]]
|-
| [[Image:figure-6.5-springmass_modes-fast.png]]
|}

'''Figure 6.5''': The notion of modes for a second-order system with real eigenvalues. The left figure shows the phase portrait and the modes corresponding to solutions that start on the eigenvectors (bold lines). The corresponding time functions are shown on the right.

<nowiki>
# 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')
</nowiki>

Figure 6.5: Modes for a second-order system with real eigenvalues

2024-11-25T05:54:29Z

Murray:

{{Figure
|Chapter=Linear Systems
|Figure number=6.5
|Sort key=605
|Figure title=The notion of modes for a second-order system with real eigenvalues.
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/figure-6.5-springmass_modes.py
|Requires=springmass.py
}}
{| border=0
|-
| rowspan=2 | [[Image:figure-6.5-springmass_modes-pp.png]]
| [[Image:figure-6.5-springmass_modes-slow.png]]
|-
| [[Image:figure-6.5-springmass_modes-fast.png]]
|}

'''Figure 6.5''': The notion of modes for a second-order system with real eigenvalues. The left figure shows the phase portrait and the modes corresponding to solutions that start on the eigenvectors (bold lines). The corresponding time functions are shown on the right.

<nowiki>
# 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')
</nowiki>

File:Figure-5.21-noise cancel.png

2024-11-25T00:48:45Z

Murray:

Figure 5.21: Simulation of noise cancellation

2024-11-25T00:48:30Z

Murray: Created page with "{{Figure |Chapter=Dynamic Behavior |Figure number=5.21 |Sort key=521 |Figure title=Simulation of noise cancellation |GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/m..."

{{Figure
|Chapter=Dynamic Behavior
|Figure number=5.21
|Sort key=521
|Figure title=Simulation of noise cancellation
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-5.18-noise_cancel.py
}}
[[Image:figure-5.21-noise_cancel.png]]

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 <math>a</math> and <math>b</math> of the filter.

<nowiki>
# 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')
</nowiki>

Bicycle dynamics

2024-11-24T20:45:11Z

Murray: Created page with "{{righttoc}} {| class="wikitable" |- ! GitHub URL | https://github.com/murrayrm/fbs2e-python/blob/main/bicycle.py |- ! Requires | python-control..."

{{righttoc}}

{| class="wikitable"
|-
! GitHub URL
| https://github.com/murrayrm/fbs2e-python/blob/main/bicycle.py
|-
! Requires
| [[https:python-control.org|python-control]]
|}

This page documents the bicycle dynamics model that is used as a running example throughout the text. A detailed description of the dynamics of this system are presented in {{chapter link|Examples}}. This page contains a description of the system, including the models and commands used to generate some of the plots in the text.

Note: the Python code on this page makes use of the [https://python-control.org Python Control Systems Library], an open source software toolbox for control systems analysis.

Figures making use of this model:
{{#ask: [[Category:Figures]] [[Requires::bicycle.py]] | format=ul | sort=Sort key}}

Examples

2024-11-24T20:44:18Z

Murray:

{{Chapter
|Chapter number=4
|Short name=examples
|Previous chapter=System Modeling
|Next chapter=Dynamic Behavior
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/index.php?title=Examples
|Chapter summary=In this chapter we present a collection of examples spanning many different fields of science and engineering. These examples are used throughout the text and in exercises to illustrate different concepts. First-time readers may wish to focus on only a few examples with which they have had the most prior experience or insight to understand the concepts of state, input, output, and dynamics in a familiar setting.
|Chapter contents=# [[Cruise control|Cruise Control]]
# [[Bicycle dynamics|Bicycle Dynamics]]
# Operational Amplifier Circuits
# Computing Systems and Networks
#* Web Server Control
#* Congestion Control
# Atomic Force Microscopy
# Drug Administration
#* Compartment Models
#* Insulin--Glucose Dynamics
# Population Dynamics
#* Logistic Growth Model
#* [[Predator-prey dynamics|Predator-Prey Models]]
:: Exercises
}}

{{Chapter footer
| First edition URL=https://www.cds.caltech.edu/~murray/amwiki/Examples.html
}}

Bicycle.py

2024-11-24T20:43:50Z

Murray: Redirected page to Bicycle dynamics

#REDIRECT [[Bicycle dynamics]]

File:Figure-5.19-bicycle stability.png

2024-11-24T20:43:13Z

Murray:

Figure 5.19: Stability plots for a bicycle moving at constant velocity

2024-11-24T20:42:56Z

Murray:

{{Figure
|Chapter=Dynamic Behavior
|Figure number=5.19
|Sort key=519
|Figure title=Stability plots for a bicycle moving at constant velocity
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-5.17-bicycle_stability.py
|Requires=bicycle.py
}}
[[Image:figure-5.19-bicycle_stability.png]]

'''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 <math>v_0</math>. 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 <math>v</math> is varied and gives a different view of the stability of the system. This type of plot is called a ''root locus diagram''.

Note: 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.

<nowiki>
# 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')
</nowiki>

Figure 5.19: Stability plots for a bicycle moving at constant velocity

2024-11-24T20:42:40Z

Murray: Created page with "{{Figure |Chapter=Dynamic Behavior |Figure number=5.19 |Sort key=519 |Figure title=Stability plots for a bicycle moving at constant velocity |GitHub URL=https://github.com/mur..."

{{Figure
|Chapter=Dynamic Behavior
|Figure number=5.19
|Sort key=519
|Figure title=Stability plots for a bicycle moving at constant velocity
|GitHub URL=https://github.com/murrayrm/fbs2e-python/blob/main/example-5.17-bicycle_stability.py
|Requires=bicycle.py
}}
[[Image:figure-5.18-predprey_bif.png]]

'''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 <math>v_0</math>. 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 <math>v</math> is varied and gives a different view of the stability of the system. This type of plot is called a ''root locus diagram''.

Note: 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.

<nowiki>
# 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')
</nowiki>

Dynamic Behavior

2024-11-24T16:53:35Z

Murray: /* Chapter Summary */

{{Chapter
|Chapter number=5
|Short name=dynamics
|Previous chapter=Examples
|Next chapter=Linear Systems
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/index.php?title=Dynamic_Behavior
|Chapter summary=In this chapter we give a broad discussion of the behavior of dynamical systems, focused on systems modeled by nonlinear differential equations. This allows us to discuss equilibrium points, stability, limit cycles and other key concepts of dynamical systems. We also introduce some methods for analyzing global behavior of solutions.
|Chapter contents=# Solving Differential Equations
# Qualitative Analysis
# Stability
# Lyapunov Stability Analysis
# Parametric and Nonlocal Behavior
# Further Reading
:: Exercises
}}
== Chapter Summary ==
<ol>
<li> We say that <math>x(t)</math> is a solution of a differential equation on the time interval <math>t_0</math> to <math>t_\text{f}</math> with initial value <math>x_0</math> if it satisfies
<center><math>
x(t_0) = x_0 \quad\text{and}\quad \dot x(t) = F(x(t)) \quad\text{for all}\quad t_0 \leq t \leq t_\text{f}.
</math></center>
We will usually assume <math>t_0 = 0</math>. For most differential equations we will encounter, there is a unique solution for a given initial condition. Numerical tools such as MATLAB and Mathematica can be used to obtain numerical solutions for <math>x(t)</math> given the function <math>F(x)</math>.</li>

<li> An ''equilibrium point'' for a dynamical system represents a point <math>x_{e}</math> such that if <math>x(0) = x_\text{e}</math> then <math>x(t) = x_\text{e}</math> for all <math>t</math>. Equilibrium points represent stationary conditions for the dynamics of a system. A ''limit cycle'' for a dynamical system is a solution <math>x(t)</math> which is periodic with some period <math>T</math>, so that <math>x(t + T) = x(t)</math> for all <math>t</math>.</li>

<li>An equilibrium point is (locally) ''stable'' if initial conditions that start near an equilibrium point stay near that equilibrium point. A equilibrium point is (locally) ''asymptotically stable'' if it is stable and, in addition, the state of the system converges to the equilibrium point as time increases. An equilibrium point is ''unstable'' if it is not stable. Similar definitions can be used to define the stability of a limit cycle.</li>

<li> Phase portraits provide a convenient way to understand the behavior of 2-dimensional dynamical systems. A phase portrait is a graphical representation of the dynamics obtained by plotting the state <math>x(t) = (x_1(t), x_2(t))</math> in the plane. This portrait is often augmented by plotting an arrow in the plane corresponding to <math>F(x)</math>, which shows the rate of change of the state. The following diagrams illustrate the basic features of a dynamical systems:
<table border=0>
<tr>
<td align=center width=33%> [[Image:doscpp.png|180px]]</td>
<td align=center width=33%> [[Image:oscpp.png|180px]]</td>
<td align=center width=33%> [[Image:stablepp.png|180px]]</td>
</tr><tr>
<td align=center>An asymptotically stable equilibrium point at <math>x = (0, 0)</math>.</td>
<td align=center>A limit cycle of radius one, with an unstable equilibrium point at <math>x = (0, 0)</math>.</td>
<td align=center>A stable equlibirum point at <math>x = (0, 0)</math> (nearby initial conditions stay nearby).</td>
</tr></table>
</li>

<li>A linear system
<center><math>
\dfrac{dx}{dt} = A x
</math></center>
is asymptotically stable if and only if all eigenvalues of <math>A</math> all have strictly negative real part and is unstable if any eigenvalue of <math>A</math> has strictly positive real part. A nonlinear system can be approximated by a linear system around an equilibrium point by using the relationship
<center><math>
\dot x = F(x_\text{e}) + \left.\frac{\partial F}{\partial x}\right|_{x_\text{e}} (x - x_\text{e}) + \text{higher order terms in}\, (x-x_\text{e}).
</math></center>
Since <math>F(x_e) = 0</math>, we can approximate the system by choosing a new
state variable <math>z = x - x_\text{e}</math> and writing the dynamics as <math>\dot z = A z</math>. The stability of the nonlinear system can be determined in a local neighborhood of the equilibrium point through its linearization.
</li>

<li>A ''Lyapunov function'' is an energy-like function <math>V:\mathbb{R}^n \to \mathbb{R}</math> that can be used to reason about the stability of an equilibrium point. We define the derivative of <math>V</math> along the trajectory of the system as
<center><math>
\dot V(x) = \dfrac{\partial V}{\partial x} \dot x = \dfrac{\partial V}{\partial x} F(x)
</math></center>
Assuming <math>x_\text{e} = 0</math> and <math>V(0) = 0</math>, the following conditions hold:
<center>
{| border=1 style="border-spacing: 2px;"
|-
| Condition on <math>V</math> || Condition on <math>\dot V</math> || Stability
|-
| align=center | <math> V(x) > 0, x \neq 0</math>
| align=center | <math>\dot V(x) \leq 0</math> for all <math>x</math>
| align=left | <math>x_e</math> stable
|-
| align=center | <math>V(x) > 0, x \neq 0 </math>
| align=center | <math>\dot V(x) < 0, x \neq 0</math>
| align=left | <math>x_e</math> asymptotically stable
|}
</center>
Stability of limit cycles can also be studied using Lyapunov functions.
</li>

<li>The ''Krasovskii-LaSalle Principle'' allows one to reason about asymptotic stability even if the time derivative of <math>V</math> is only negative semi-definite (<math>\leq 0</math> rather than <math>< 0</math>). Let <math>V: \mathbb{R}^n \to \mathbb{R}</math> be a ''positive definite function'', <math>V(x) > 0</math> for all <math> x \neq 0</math> and <math>V(0) = 0</math>, such that
<center>
<math>\dot V(x) \leq 0</math> on the compact set <math>\Omega_c = \{x \in \mathbb{R}^n:V(x) \leq c\}</math>.
</center>
Then as <math>t \to \infty</math>, the trajectory of the system will converge to the largest invariant set inside
<center><math>S = \{x \in \Omega_c:\dot V(x) = 0\}</math>.</center>
In particular, if <math>S</math> contains no invariant sets other than <math>x = 0</math>, then 0 is asymptotically stable.
</li>

<li>The ''global behavior'' of a nonlinear system refers to dynamics of the system far away from equilibrium points. The ''region of attraction'' of an asymptotically stable equilirium point refers to the set of all initial conditions that converge to that equilibrium point. An equilibrium point is said to be ''globally asymptotically stable'' if all initial conditions converge to that equilibrium point. Global stability can be checked by finding a Lyapunov function that is globally positive definition with time derivative globally negative definite.</li>
</ol>

{{Chapter footer
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/Dynamic_Behavior.html
}}

Bibliography

2024-11-24T16:44:50Z

Murray:

{{Chapter header
|Chapter number=
|Chapter title={{PAGENAME}}
|Previous chapter=Architecture and System Design
|Next chapter=FBS
}}

The bibliography contains information on all citations in the book. Note that the bibliography uses an alphanumberic tag format, so that the bibliography can be updated as the text is written.

{{FBS pdf|PDF|fbs-backmatter}} ({{FBS date}})

== Online References ==

The following notes, reports and books are available online.


'''[Bec05]''' [http://dx.doi.org/10.1103/RevModPhys.77.783 "Feedback for physicists: A tutorial essay on control"], ''Review of Modern Physics'', 77:783, 2005 ([http://www.sfu.ca/chaos/Publications/papers/rmp_reprint05.pdf reprint]).



'''[BFS]''' ''[[Biomolecular Feedback Systems]]'', D. Del Vecchio and R. M. Murray, 2014.



'''[DFT92]''' [https://www.control.utoronto.ca/people/profs/francis/dft.pdf ''Feedback Control Theory''], J. C. Doyle, B. Francis, A. Tannenbaum, MacMillan, 1992.



'''[Fra15]''' [https://web.archive.org/web/20220423093734/http://www.scg.utoronto.ca/~francis/main.pdf ''Classical Control]'', B. A. Francis, University of Toronto, 2015.



'''[Lew03]''' [https://mast.queensu.ca/~andrew/teaching/pdf/332-notes.pdf ''A Mathematical Approach to Classical Control''], A. D. Lewis, 2003.



'''[LST]''' ''[[Supplement: Linear Systems Theory|Linear Systems Theory]]'' (supplemental notes), R. M. Murray, 2020.



'''[OBC]''' [[Supplement: Optimization-Based Control|''Optimization-Based Control'']] (supplemental notes), R. M. Murray, 2010.



'''[Mur03]''' [http://www.cds.caltech.edu/~murray/cdspanel ''Control in an Information Rich World''], R. M. Murray (ed), SIAM, 2003.



'''[Son98]''' [http://www.sontaglab.org/FTPDIR/sontag_mathematical_control_theory_springer98.pdf ''Mathematical Control Theory''], E. D. Sontag, Springer, 1998.

Bibliography

2024-11-24T16:44:31Z

Murray:

{{Chapter header
|Chapter number=
|Chapter title={{PAGENAME}}
|Previous chapter=Architecture and System Design
|Next chapter=FBS
}}

The bibliography contains information on all citations in the book. Note that the bibliography uses an alphanumberic tag format, so that the bibliography can be updated as the text is written.

{{FBS pdf|PDF|fbs-backmatter}} ({{FBS date}})

== Online References ==

The following notes, reports and books are available online.


'''[Bec05]''' [http://dx.doi.org/10.1103/RevModPhys.77.783 "Feedback for physicists: A tutorial essay on control"], ''Review of Modern Physics'', 77:783, 2005 ([http://www.sfu.ca/chaos/Publications/papers/rmp_reprint05.pdf reprint]).



'''[BFS]''' ''[[Biomolecular Feedback Systems]]'', D. Del Vecchio and R. M. Murray, 2014.



'''[DFT92]''' [https://www.control.utoronto.ca/people/profs/francis/dft.pdf ''Feedback Control Theory''], J. C. Doyle, B. Francis, A. Tannenbaum, MacMillan, 1992.



'''[Fra15]''' [http:web.archive.org/web/20220423093734/http://www.scg.utoronto.ca/~francis/main.pdf ''Classical Control]'', B. A. Francis, University of Toronto, 2015.



'''[Lew03]''' [https://mast.queensu.ca/~andrew/teaching/pdf/332-notes.pdf ''A Mathematical Approach to Classical Control''], A. D. Lewis, 2003.



'''[LST]''' ''[[Supplement: Linear Systems Theory|Linear Systems Theory]]'' (supplemental notes), R. M. Murray, 2020.



'''[OBC]''' [[Supplement: Optimization-Based Control|''Optimization-Based Control'']] (supplemental notes), R. M. Murray, 2010.



'''[Mur03]''' [http://www.cds.caltech.edu/~murray/cdspanel ''Control in an Information Rich World''], R. M. Murray (ed), SIAM, 2003.



'''[Son98]''' [http://www.sontaglab.org/FTPDIR/sontag_mathematical_control_theory_springer98.pdf ''Mathematical Control Theory''], E. D. Sontag, Springer, 1998.

Predator-prey dynamics

2024-11-24T16:40:07Z

Murray:

{| class="wikitable"
|-
! GitHub URL
| https://github.com/murrayrm/fbs2e-python/blob/main/predprey.py
|-
! Requires
| [[https:python-control.org|python-control]]
|}

This page documents the predator-prey model that is used as a running example throughout the text. A detailed description of the dynamics of this system are presented in {{chapter link|Examples}}. This page contains a description of the system, including the models and commands used to generate some of the plots in the text.

Note: the Python code on this page makes use of the [https://python-control.org Python Control Systems Library], an open source software toolbox for control systems analysis.

Figures making use of this model:
{{#ask: [[Category:Figures]] [[Requires::predprey.py]] | format=ul | sort=Sort key}}

Cruise control

2024-11-24T16:39:50Z

Murray:

{{righttoc}}

{| class="wikitable"
|-
! GitHub URL
| https://github.com/murrayrm/fbs2e-python/blob/main/cruise.py
|-
! Requires
| [[https:python-control.org|python-control]]
|}

This page documents the cruise control system that is used as a running example throughout the text. A detailed description of the dynamics of this system are presented in {{chapter link|Examples}}. This page contains a description of the system, including the models and commands used to generate some of the plots in the text.

Note: the Python code on this page makes use of the [https://python-control.org Python Control Systems Library], an open source software toolbox for control systems analysis.

Figures making use of this model:
{{#ask: [[Category:Figures]] [[Requires::cruise.py]] | format=ul | sort=Sort key}}

== Introduction ==

[[Image:cruise-block.png|right|300px|border|A feedback system for controlling the velocity of a vehicle. In the block diagram, the velocity of the vehicle is measured and compared to the desired velocity within the “Compute” block. Based on the difference in the actual and desired velocities, the throttle (or brake) is used to modify the force applied to the vehicle by the engine, drivetrain, and wheels.]]
Cruise control is the term used to describe a control system that regulates the speed of an automobile. Cruise control was commercially introduced in 1958 as an option on the Chrysler Imperial. The basic operation of a cruise controller is to sense the speed of the vehicle, compare this speed to a desired reference, and then accelerate or decelerate the car as required. The figure to the right shows a block diagram of this feedback system.

[[Image:cruise-speedresp.png|right|300px|border|Response to a disturbance. 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.
]]
A simple control algorithm for controlling the speed is to use a "proportional plus integral" feedback. In this algorithm, we choose the amount of gas flowing to the engine based on both the error between the current and desired speed, and the integral of that error. The plot on the right shows the results of this feedback for a step change in the desired speed and a variety of different masses for the car (which might result from having a different number of passengers or towing a trailer). Notice that independent of the mass (which varies by 25% of the total weight of the car), the steady state speed of the vehicle always approaches the desired speed and achieves that speed within approximately 10-15 seconds. Thus the performance of the system is robust with respect to this uncertainty.

== Dynamic model ==

To develop a mathematical model we start with a force balance for the car body. Let <math>v</math> be the speed of the car, <math>m</math> the total mass (including passengers), <math>F</math> the force generated by the contact of the wheels with the road, and <math>F_\text{d}</math> the disturbance force due to gravity, friction and aerodynamic drag. The equation of motion of the car is simply
<center><math>
m \frac{dv}{dt} = F - F_\text{d}.
</math></center>

The force <math>F</math> is generated by the engine, whose torque is proportional to the rate of fuel injection, which is itself proportional to a control signal <math>0 \leq u \leq 1</math> that controls the throttle position. The torque also depends on engine speed <math>\omega</math>. A simple representation of the torque at full throttle is given by the torque curve
<center><math>
T(\omega) = T_\text{m} \left(1 - \beta \biggl(\frac{\omega}{\omega_\text{m}} - 1\biggr)^2
\right),
</math></center>
where the maximum torque <math>T_\text{m}</math> is obtained at engine speed <math>\omega_\text{m}</math>. Typical parameters are <math>T_m = 190</math> Nm, <math>\omega_m</math> = 420 rad/s (about 4000 RPM) and <math>\beta = 0.4</math>.

Let <math>n</math> be the gear ratio and <math>r</math> the wheel radius. The engine speed is related to the velocity through the expression
<center><math>
\omega = \frac{n}{r} v =: \alpha_n v,
</math></center>
and the driving force can be written as
<center><math>
F = \frac{nu}{r} T(\omega) = \alpha_n u T(\alpha_n v).
</math></center>
[[Image:cruise-gearcurves.png|right|400px|border|Torque curves for typical car engine. The graph on the left shows the torque generated by the engine as a function of the angular velocity of the engine, while the curve on the right shows torque as a function of car speed for different gears.]]
Typical values of <math>\alpha_n</math> for gears 1 through 5 are <math>\alpha_1 = 40</math>, <math>\alpha_2 = 25</math>, <math>\alpha_3 = 16</math>, <math>\alpha_4 = 12</math> and <math>\alpha_5 = 10</math>. The inverse of <math>\alpha_n</math> has a physical interpretation as the effective wheel radius. The figure to the right shows the torque as a function of vehicle speed. The figure shows that the effect of the gear is to "flatten" the torque curve so that an almost full torque can be obtained almost over the whole speed range.

The disturbance force <math>F_\text{d}</math> has three major components: <math>F_\text{g}</math>, the forces due to gravity; <math>F_\text{r}</math>, the forces due to rolling friction; and <math>F_\text{a}</math>, the aerodynamic drag. Letting the slope of the
road be <math>\theta</math>, gravity gives the force <math>F_\text{g} = m g \sin\theta</math>, where <math>g = 9.8\,
\text{m}/\text{s}^2</math> is the gravitational constant. A simple model of rolling friction is
<center><math>
F_\text{r} = m g C_\text{r}\, \text{sgn}(v),
</math></center>
where <math>C_\text{r}</math> is the coefficient of rolling friction and sgn(<math>v</math>) is the sign of <math>v</math> or zero if <math>v = 0</math>. A typical value for the coefficient of rolling friction is <math>C_\text{r} = 0.01</math>. Finally, the aerodynamic drag is proportional to the square of the speed:
<center><math>
F_\text{a} = \frac{1}{2} \rho C_\text{d} A v^2,
</math></center>
where <math>\rho</math> is the density of air, <math>C_d</math> is the shape-dependent aerodynamic drag coefficient and <math>A</math> is the frontal area of the car. Typical parameters are <math>\rho = </math> 1.3 k/m<math>{}^3</math>, <math>C_\text{d} = 0.32</math> and <math>A =</math> 2.4 m<math>{}^2</math>.

=== Python code ===

The model for the system above can be built using the [[https:python-control.org|Python Control Toolbox]]. The code blocks in this section can be used to generate the plots above. The vehicle dynamics and PI controller are defined in [[https:github.com/murrayrm/fbs2e-python/blob/main/cruise.py|cruise.py]] (via GitHub).

{{Code block|Package initialization}}
{{Code start}}
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from math import pi
import control as ct
{{Code end}}

{{Code block|Vehicle model}}
{{Code start}}
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
{{Code end}}

{{Code block|Engine model}}
{{Code start}}
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)
{{Code end}}

{{Code block|Input/output model for the vehicle system}}
{{Code start}}
vehicle = ct.NonlinearIOSystem(
vehicle_update, None, name='vehicle',
inputs = ('u', 'gear', 'theta'), outputs = ('v'), states=('v'))
{{Code end}}

{{Code block|Input/output torque curves (plot)}}
{{Code start}}
# Figure 4.2a - single torque curve as function of omega
omega_range = np.linspace(0, 700, 701)
plt.subplot(2, 2, 1)
plt.plot(omega_range, [motor_torque(w) for w in omega_range])
plt.xlabel('Angular velocity $\omega$ [rad/s]')
plt.ylabel('Torque $T$ [Nm]')
plt.grid(True, linestyle='dotted')

# Figure 4.2b - torque curves in different gears, as function of velocity
plt.subplot(2, 2, 2)
v_range = np.linspace(0, 70, 71)
alpha = [40, 25, 16, 12, 10]
for gear in range(5):
omega_range = alpha[gear] * v_range
plt.plot(v_range, [motor_torque(w) for w in omega_range],
color='blue', linestyle='solid')

# Set up the axes and style
plt.axis([0, 70, 100, 200])
plt.grid(True, linestyle='dotted')

# Add labels
plt.text(11.5, 120, '$n$=1')
plt.text(24, 120, '$n$=2')
plt.text(42.5, 120, '$n$=3')
plt.text(58.5, 120, '$n$=4')
plt.text(58.5, 185, '$n$=5')
plt.xlabel('Velocity $v$ [m/s]')
plt.ylabel('Torque $T$ [Nm]')

plt.suptitle('Torque curves for typical car engine');
plt.tight_layout()
plt.savefig("cruise-gearcurves.png", bbox_inches='tight')
{{Code end}}

{{Code block|PI controller}}
{{Code start}}
# Construct a PI controller with rolloff, as a transfer function
Kp = 0.5 # proportional gain
Ki = 0.1 # integral gain
control_pi = ct.tf2io(
ct.TransferFunction([Kp, Ki], [1, 0.01*Ki/Kp]),
name='control', inputs='u', outputs='y')

cruise_pi = ct.InterconnectedSystem(
(vehicle, control_pi), 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'))
{{Code end}}

{{Code block|Plotting function}}
{{Code start}}
# Define a generator for creating a "standard" cruise control plot
def cruise_plot(sys, t, y, t_hill=5, vref=20, antiwindup=False, linetype='b-',
subplots=[None, None]):
# Figure out the plot bounds and indices
v_min = vref-1.2; v_max = vref+0.5; v_ind = sys.find_output('v')
u_min = 0; u_max = 2 if antiwindup else 1; u_ind = sys.find_output('u')

# Make sure the upper and lower bounds on v are OK
while max(y[v_ind]) > v_max: v_max += 1
while min(y[v_ind]) < v_min: v_min -= 1

# Create arrays for return values
subplot_axes = subplots.copy()

# Velocity profile
if subplot_axes[0] is None:
subplot_axes[0] = plt.subplot(2, 1, 1)
else:
plt.sca(subplots[0])
plt.plot(t, y[v_ind], linetype)
plt.plot(t, vref*np.ones(t.shape), 'k-')
plt.plot([t_hill, t_hill], [v_min, v_max], 'k--')
plt.axis([0, t[-1], v_min, v_max])
plt.xlabel('Time $t$ [s]')
plt.ylabel('Velocity $v$ [m/s]')

# Commanded input profile
if subplot_axes[1] is None:
subplot_axes[1] = plt.subplot(2, 1, 2)
else:
plt.sca(subplots[1])
plt.plot(t, y[u_ind], 'r--' if antiwindup else linetype)
plt.plot([t_hill, t_hill], [u_min, u_max], 'k--')
plt.axis([0, t[-1], u_min, u_max])
plt.xlabel('Time $t$ [s]')
plt.ylabel('Throttle $u$')

# Applied input profile
if antiwindup:
plt.plot(t, np.clip(y[u_ind], 0, 1), linetype)
plt.legend(['Commanded', 'Applied'], frameon=False)

return subplot_axes
{{Code end}}

{{Code block|Trajectory description}}
{{Code start}}
# 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)

# Effect of a hill at t = 5 seconds
theta_hill = np.array([
0 if t <= 5 else
4./180. * pi * (t-5) if t <= 6 else
4./180. * pi for t in T])
{{Code end}}

{{Code block|Simulated responses (plot)}}
{{Code start}}
# Simulate and plot
plt.figure()
plt.suptitle('Response to change in road slope')

subplots = [None, None]
linecolor = ['red', 'blue', 'green']
handles = []
for i, m in enumerate([1200, 1600, 2000]):
# Compute the equilibrium state for the system
X0, U0 = ct.find_eqpt(
cruise_pi, [vref[0], 0], [vref[0], gear[0], theta0[0]],
iu=[1, 2], y0=[vref[0], 0], iy=[0], params={'m':m})

t, y = ct.input_output_response(
cruise_pi, T, [vref, gear, theta_hill], X0, params={'m':m})

subplots = cruise_plot(cruise_pi, t, y, t_hill=5, subplots=subplots,
linetype=linecolor[i][0] + '-')
handles.append(mpl.lines.Line2D([], [], color=linecolor[i],
linestyle='-', label="m = %d" % m))

# Add labels to the plots
plt.sca(subplots[0])
plt.ylabel('Speed [m/s]')
plt.legend(handles=handles, frameon=False, loc='lower right');

plt.sca(subplots[1])
plt.ylabel('Throttle')
plt.xlabel('Time [s]');
plt.savefig("cruise-speedresp.png", bbox_inches='tight')
{{Code end}}

=== Linearized Dynamics ===

To explore the behavior of the cruise control system near the equilibrium point we
will linearize the system. A Taylor series expansion of
the dynamics around the equilibrium point gives
<center><math>
\frac{d(v-v_\text{e})}{dt} = -a (v - v_\text{e})
- b_\text{g} (\theta - \theta_\text{e}) + b (u - u_\text{e})
+ \text{higher-order terms,}
</math></center>
where
<center><math>
a = \frac{\rho C_\text{d} A |v_\text{e}|
- u_\text{e} \alpha_n^2 T'(\alpha_n v_\text{e})}{m}, \qquad
b_\text{g} = g \cos{\theta_\text{e}}, \qquad
b = \frac{\alpha_n T(\alpha_n v_\text{e})}{m}.
</math></center>
[[Image:cruise-linear_vs_nonlinear.png|right|border|400px|Simulated response of a vehicle with PI cruise control as it climbs a hill with a slope of 4 degrees (smaller velocity deviation/throttle) and a slope of 6 degrees (larger velocity deviation/throttle). The solid line is the simulation based on a nonlinear model, and the dashed line shows the corresponding simulation using a linear model. The controller gains are kp = 0.5 and ki = 0.1 and include anti-windup compensation (described in more detail in below).]]
Notice that the term corresponding to rolling friction disappears if <math>v > 0</math>. For a car in fourth gear with <math>v_\text{e} = 20</math> m/s, <math>\theta_\text{e} = 0</math>, and the numerical values from above, the equilibrium value for the throttle is <math>u_\text{e}=0.1687</math> and the parameters are <math>a = 0.01</math>, <math>b = 1.32</math>, and <math>b_\text{g} = 9.8</math>. This linear model describes how small perturbations in the velocity about the nominal speed evolve in time.

We apply the PI controller above to the case where the car is running with constant speed on a horizontal road and the system has stabilized so that the vehicle speed and the controller output are constant. The figure to the right shows what happens when the car encounters a hill with a slope of 4 degrees and a hill with a slope of 6 degrees at time <math>t =\,</math> 5 seconds. The results for the nonlinear model are dashed curves and those for the linear model are solid curves. The differences between the curves are very small (especially for <math>\theta =\,</math> 4 degrees), and control design based on the linearized model is thus validated. 

=== Python code ===

{{Code block|Linearized dynamics}}
{{Code start}}
# Find the equilibrium point for the system
Xeq, Ueq = ct.find_eqpt(
vehicle, [vref[0]], [0, gear[0], theta0[0]], y0=[vref[0]], iu=[1, 2])
print("Xeq = ", Xeq, "\nUeq = ", Ueq)

# Compute the linearized system at the eq pt
vehlin = ct.linearize(vehicle, Xeq, [Ueq[0], gear[0], 0], name='vehlin', copy=True)
{{Code end}}

{{Code block|PI controller with anti-windup compensation}}
{{Code start}}
def pi_update(t, x, u, params={}):
# Get the controller parameters that we need
ki = params.get('ki', 0.1)
kaw = params.get('kaw', 2) # anti-windup gain

# Assign variables for inputs and states (for readability)
v = u[0] # current velocity
vref = u[1] # reference velocity
z = x[0] # integrated error

# Compute the nominal controller output (needed for anti-windup)
u_a = pi_output(t, x, u, params)

# Compute anti-windup compensation (scale by ki to account for structure)
u_aw = kaw/ki * (np.clip(u_a, 0, 1) - u_a) if ki != 0 else 0

# State is the integrated error, minus anti-windup compensation
return (vref - v) + u_aw

def pi_output(t, x, u, params={}):
# Get the controller parameters that we need
kp = params.get('kp', 0.5)
ki = params.get('ki', 0.1)

# Assign variables for inputs and states (for readability)
v = u[0] # current velocity
vref = u[1] # reference velocity
z = x[0] # integrated error

# PI controller
return kp * (vref - v) + ki * z

control_pi_aw = ct.NonlinearIOSystem(
pi_update, pi_output, name='control',
inputs = ['v', 'vref'], outputs = ['u'], states = ['z'],
params = {'kp':0.5, 'ki':0.1})

# Create a closed loop controller for the linear system
cruise_lin = ct.InterconnectedSystem(
(vehlin, control_pi_aw), name='cruise_lin',
connections=(
('vehlin.u', 'control.u'),
('control.v', 'vehlin.v')),
inplist=('control.vref', 'vehlin.gear', 'vehlin.theta'),
outlist=('control.u', 'vehlin.v'), outputs=['u', 'v'])

# Create a closed loop controller for the nonlinear system
cruise_nonlin = ct.InterconnectedSystem(
(vehicle, control_pi_aw), name='cruise_nonlin',
connections=(
('vehicle.u', 'control.u'),
('control.v', 'vehicle.v')),
inplist=('control.vref', 'vehicle.gear', 'vehicle.theta'),
outlist=('control.u', 'vehicle.v'), outputs=['u', 'v'])
{{Code end}}

{{Code block|Simulations}}
{{Code start}}
# Linear response
X0_lin, U0 = ct.find_eqpt(
cruise_lin, [vref[0], 0], [vref[0], gear[0], theta0[0]],
y0=[0, vref[0]], iu=[1, 2], iy=[1])

t, y = ct.input_output_response(cruise_lin, T, [vref, gear, theta_hill], X0_lin)
subplots = cruise_plot(cruise_lin, t, y, t_hill=5, linetype='r-')

# Nonlinear response
X0_nonlin, U0 = ct.find_eqpt(
cruise_nonlin, [vref[0], 0], [vref[0], gear[0], theta0[0]],
y0=[0, vref[0]], iu=[1, 2], iy=[1])

t, y = ct.input_output_response(cruise_nonlin, T, [vref, gear, theta_hill], X0_nonlin)
subplots = cruise_plot(cruise_nonlin, t, y, t_hill=5, subplots=subplots, linetype='b--')

# Add a legend to identify linear vs nonlinear
plt.legend(['linear', 'nonlinear'], frameon=False, loc='lower right')

# Add two more simulations for 6 degree hill instead of 4 degree
t, y = ct.input_output_response(cruise_lin, T, [vref, gear, theta_hill*1.5], X0_lin)
subplots = cruise_plot(cruise_lin, t, y, t_hill=5, subplots=subplots, linetype='r-')

t, y = ct.input_output_response(cruise_nonlin, T, [vref, gear, theta_hill*1.5], X0_nonlin)
subplots = cruise_plot(cruise_nonlin, t, y, t_hill=5, subplots=subplots, linetype='b--')

# Add titles and legends and save the figure
plt.suptitle('Linear versus nonlinear model response')
plt.sca(subplots[1]); plt.axis([0, 25, 0, 1.2])

plt.savefig("cruise-linear_vs_nonlinear.png", bbox_inches='tight')
{{Code end}}

== State Space Control ==

The linearized dynamics of
the process around an equilibrium point <math>v_\text{e}</math>, <math>u_\text{e}</math> are given by
<center><math>
\begin{aligned}
\frac{dx}{dt} &= -a x - b_\text{g} \theta + b w, \\
y &= v = x + v_\text{e},
\end{aligned}
</math></center>
where <math>x = v - v_\text{e}</math>, <math> = u - u_\text{e}</math>, <math>m</math> is the mass of
the car, and <math>\theta</math> is the angle of the road. The constants <math>a</math>,
<math>b</math>, and <math>b_\text{g}</math> depend on the properties of the car and are
given above.

If we augment the system with an integrator, the system dynamics become
<center><math>
\begin{aligned}
\frac{dx}{dt} &= -a x - b_\text{g} \theta + b w, \\
\frac{dz}{dt} &= y - v_\text{r} = v_\text{e} + x - v_\text{r},
\end{aligned}
</math></center>
or, in state space form,
<center><math>
\frac{d}{dt} \begin{bmatrix} x \\ z \end{bmatrix} = \begin{bmatrix}
-a & 0 \\
1 & 0
\end{bmatrix} \begin{bmatrix} x \\ z \end{bmatrix} +
\begin{bmatrix} b \\ 0 \end{bmatrix} w +
\begin{bmatrix} -b_\text{g} \\ 0 \end{bmatrix} \theta +
\begin{bmatrix} 0 \\ v_\text{e} - v_\text{r} \end{bmatrix}.
</math></center>
Note that when the system is at equilibrium, we have that <math>\dot z = 0</math>,
which implies that the vehicle speed <math>v = v_\text{e} + x</math> should be
equal to the desired reference speed <math>v_\text{r}</math>. Our controller will be of
the form
<center><math>
\begin{aligned}
\frac{dz}{dt} &= y - v_\text{r}, \\
w &= -k_\text{p} x - k_\text{i} z + k_\text{f} v_\text{r},
\end{aligned}
</math></center>
and the gains <math>k_\text{p}</math>, <math>k_\text{i}</math>, and <math>k_\text{f}</math> will be chosen to
stabilize the system and provide the correct input for the reference speed.

Assume that we wish to design the closed loop system to have the
characteristic polynomial
<center><math>
\lambda(s) = s^2 + a_1 s + a_2.
</math></center>
Setting the disturbance <math>\theta = 0</math>,
the characteristic polynomial of the closed loop system is given by
<center><math>
\det\bigl(sI - (A - BK)\bigr) = s^2 + (b k_\text{p} + a) s + b k_\text{i},
</math></center>
and hence we set
<center><math>
k_\text{p} = \frac{a_1 - a}{b}, \qquad k_\text{i} = \frac{a_2}{b}, \qquad
k_\text{f} = {-1}/\bigl(C (A-BK)^{-1} B\bigr) = \frac{a_1}{b}.
</math></center>
[[Image:cruise-statefbk.png|right|border|400px|Velocity and throttle for a car with cruise control based on state feedback (dashed) and state feedback with integral action (solid). The controller with integral action is able to adjust the throttle to compensate for the effect of the hill and maintain the speed at the reference value of vr = 20 m/s. The controller gains are kp = 0.5 and ki = 0.1.]]
The resulting controller stabilizes the system and hence brings
<math>\dot z = y - v_\text{r}</math> to zero, resulting in perfect tracking.
Notice that even if we have a small error in the values of the
parameters defining the system, as long as the closed loop eigenvalues
are still stable, then the tracking error will approach zero. Thus
the exact calibration required in
our previous approach (using <math>k_\text{f}</math>) is not needed here.
Indeed, we can even choose <math>k_\text{f} = 0</math> and let the feedback
controller do all of the work. However, <math>k_\text{f}</math> does influence
the transient response to reference signals and setting it properly
will generally give a more favorable response.

Integral feedback can also be used to compensate for constant
disturbances. The figure to the right shows the results of
a simulation in which the car encounters a hill with angle
<math>\theta = 4</math> degrees at <math>t = </math> 5 seconds. The
steady-state values of the throttle for a state feedback controller
and a controller with integral action are very close, but the
corresponding values of the car velocity are quite different. The
reason for this is that the zero frequency gain from throttle to
velocity is <math>-b/a=130</math> is high.
The stability of the system is not affected by this external
disturbance, and so we once again see that the car's velocity converges
to the reference speed. This ability to handle constant disturbances
is a general property of controllers with integral feedback.

=== Python code ===

{{Code block|State feedback controller}}
{{Code start}}
def sf_update(t, z, u, params={}):
y, r = u[1], u[2]
return y - r

def sf_output(t, z, u, params={}):
# Get the controller parameters that we need
K = params.get('K', 0)
ki = params.get('ki', 0)
kf = params.get('kf', 0)
xd = params.get('xd', 0)
yd = params.get('yd', 0)
ud = params.get('ud', 0)

# Get the system state and reference input
x, y, r = u[0], u[1], u[2]

return ud - K * (x - xd) - ki * z + kf * (r - yd)

# Create the input/output system for the controller
control_sf = ct.NonlinearIOSystem(
sf_update, sf_output, name='control',
inputs=('x', 'y', 'r'),
outputs=('u'),
states=('z'))

# Create the closed loop system for the state space controller
cruise_sf = ct.InterconnectedSystem(
(vehicle, control_sf), name='cruise',
connections=(
('vehicle.u', 'control.u'),
('control.x', 'vehicle.v'),
('control.y', 'vehicle.v')),
inplist=('control.r', 'vehicle.gear', 'vehicle.theta'),
outlist=('control.u', 'vehicle.v'), outputs=['u', 'v'])

# Construct the gain matrices for the system
A, B, C = vehlin.A, vehlin.B[0, 0], vehlin.C
K = 0.5
kf = -1 / (C * np.linalg.inv(A - B * K) * B)

# Compute the steady state velocity and throttle setting
xd = Xeq[0]
ud = Ueq[0]
yd = vref[-1]
{{Code end}}

{{Code block|Simulation (plot)}}
{{Code start}}
# Response of the system with no integral feedback term
t, y_sfb = ct.input_output_response(
cruise_sf, T, [vref, gear, theta_hill], [Xeq[0], 0],
params={'K':K, 'ki':0.0, 'kf':kf, 'xd':xd, 'ud':ud, 'yd':yd})
subplots = cruise_plot(cruise_sf, t, y_sfb, t_hill=5, linetype='b--')

# Response of the system with state feedback + integral action
t, y_sfb_int = ct.input_output_response(
cruise_sf, T, [vref, gear, theta_hill], [Xeq[0], 0],
params={'K':K, 'ki':0.1, 'kf':kf, 'xd':xd, 'ud':ud, 'yd':yd})
cruise_plot(cruise_sf, t, y_sfb_int, t_hill=5, linetype='b-', subplots=subplots)

# Add title and legend
plt.suptitle('Cruise control with state feedback, integral action')
p_line = mpl.lines.Line2D([], [], color='blue', linestyle='--', label='State feedback')
pi_line = mpl.lines.Line2D([], [], color='blue', linestyle='-', label='w/ integral action')
plt.legend(handles=[p_line, pi_line], frameon=False, loc='lower right')
plt.savefig("cruise-statefbk.png", bbox_inches='tight')
{{Code end}}

== PID Control ==

The dynamics of the system are given by a first-order system with the transfer function
<center><math>
P(s)=\frac{b}{s+a}.
</math></center>
With a PI controller the closed loop system has the characteristic polynomial
<center><math>
s (s + a) + b k_\text{p} s + b k_\text{i} = s^2 + (a + b k_\text{p}) s + b k_\text{i}.
</math></center>
The closed loop poles can thus be assigned arbitrary values by proper
choice of the controller gains <math>k_\text{p}</math> and <math>k_\text{i}</math>.
Requiring that the closed loop system have the characteristic polynomial
<center><math>
p(s) = s^2+a_1s+a_2,
</math></center>
we find that the controller parameters are
<center><math>
k_\text{p}=\frac{a_1-a}{b},\qquad k_\text{i}=\frac{a_2}{b}.
</math></center>
[[Image:cruise-pid-sweep-zeta.png|right|border|400px|Cruise control using PI feedback. The step responses for the error and input illustrate the effect of parameters zeta_c and omega_c on the response of a car with cruise control. The slope of the road changes linearly from 0 degrees to 4 degrees between t = 5 and 6 seconds. Responses for omega_c = 0.5 and zeta_c = 0.5, 1, and 2. Choosing
zeta_c >= 1 gives no overshoot in the velocity v.]]
If we require a response of the closed loop system that is slower than
that of the open loop system, a reasonable choice is <math>a_1 = a + \alpha</math>
and <math>a_2 = \alpha a</math>, where <math>\alpha < a</math> determines the closed loop
response. If a response faster than that of the open loop system is
required, a possible choice is <math>a_1 = 2 \zeta_\text{c} \omega_\text{c}</math> and
<math>a_2 = \omega_\text{c}^2</math>,
where <math>\omega_\text{c}</math> and <math>\zeta_\text{c}</math> are the undamped natural
frequency and damping ratio of the dominant mode.

To design a PI controller we choose <math>\zeta_\text{c} = 1</math> to obtain a response without overshoot. The choice of <math>\omega_\text{c}</math> is a compromise between response speed and control actions: a large value gives a fast response, but it requires fast control action. The largest velocity error decreases with increasing <math>\omega_\text{c}</math>, but the control signal also changes more rapidly. In the simple model it was assumed that the force responds instantaneously to throttle commands. For rapid changes there may be additional dynamics that have to be accounted for. There are also physical limits to the rate of change of the force, which also restricts the admissible value of <math>\omega_\text{c}</math>. A reasonable choice of <math>\omega_\text{c}</math> is in the range 0.5--1.0.

=== Python code ===

{{Code block|Get transfer function parameters}}
{{Code start}}
# Get the transfer function from throttle input + hill to vehicle speed
P = ct.ss2tf(vehlin[0, 0])

# Construction a controller that cancels the pole
a = -P.pole()[0]
b = np.real(P(0)) * a
{{Code end}}

{{Code block|Fix <math>\omega_\text{c}</math> and vary <math>\zeta_\text{c} (plot)</math>}}
{{Code start}}
# Fix \omega_0 and vary \zeta
w0 = 0.5
subplots = [None, None]
for zeta in [0.5, 1, 2]:
# Create the controller transfer function (as an I/O system)
kp = (2*zeta*w0 - a)/b
ki = w0**2 / b
control_tf = ct.tf2io(
ct.TransferFunction([kp, ki], [1, 0.01*ki/kp]),
name='control', inputs='u', outputs='y')

# Construct the closed loop system by interconnecting process and controller
cruise_tf = ct.InterconnectedSystem(
(vehicle, control_tf), 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'))

# Plot the velocity response
X0, U0 = ct.find_eqpt(
cruise_tf, [vref[0], 0], [vref[0], gear[0], theta_hill[0]],
iu=[1, 2], y0=[vref[0], 0], iy=[0])

t, y = ct.input_output_response(cruise_tf, T, [vref, gear, theta_hill], X0)
subplots = cruise_plot(cruise_tf, t, y, t_hill=5, subplots=subplots)

plt.suptitle('PID controller with varying zeta')
plt.savefig("cruise-pid-sweep-zeta.png", bbox_inches='tight')
{{Code end}}

{{Code block|Fix vary <math>\zeta_\text{c}</math> and <math>\omega_\text{c} (plot)</math>}}
{{Code start}}
zeta = 1
subplots = [None, None]
for w0 in [0.2, 0.5, 1]:
# Create the controller transfer function (as an I/O system)
kp = (2*zeta*w0 - a)/b
ki = w0**2 / b
control_tf = ct.tf2io(
ct.TransferFunction([kp, ki], [1, 0.01*ki/kp]),
name='control', inputs='u', outputs='y')

# Construct the closed loop system by interconnecting process and controller
cruise_tf = ct.InterconnectedSystem(
(vehicle, control_tf), 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'))

# Plot the velocity response
X0, U0 = ct.find_eqpt(
cruise_tf, [vref[0], 0], [vref[0], gear[0], theta_hill[0]],
iu=[1, 2], y0=[vref[0], 0], iy=[0])

t, y = ct.input_output_response(cruise_tf, T, [vref, gear, theta_hill], X0)
subplots = cruise_plot(cruise_tf, t, y, t_hill=5, subplots=subplots)

plt.suptitle('PID controller with varying omega')
plt.savefig("cruise-pid-sweep-omega.png", bbox_inches='tight')
{{Code end}}

=== Anti-windup compensation ===

The windup effect is illustrated in the left figure below, which shows what happens when a car encounters a hill that is so steep (6 degrees) that the throttle saturates when the cruise controller attempts to maintain speed. When encountering the slope at time <math>t=5</math>, the velocity decreases and the throttle increases to generate more torque. However, the torque required is so large that the throttle saturates. The error decreases slowly because the torque generated by the engine is just a little larger than the torque required to compensate for gravity. The error is large and the integral continues to build up until the error reaches zero at time 25, but the controller output is still larger than the saturation limit and the actuator remains saturated. The integral term starts to decrease and the velocity settles to the desired value at time <math>t = 40</math>. Also notice the large overshoot.

The right figure below shows what happens when a controller with anti-windup is applied to the system. Because of the feedback from the actuator model, the output of the integrator is quickly reset to a value such that the controller output is at the saturation limit. The behavior is drastically different from that on the left and the large overshoot is avoided. The tracking gain used in the simulation is <math>k_\text{aw} = 2</math> which is an order of magnitude larger than the integral gain <math>k_\text{i} = 0.2</math>.

{| align=center
|-
| width=45% |
[[Image:cruise-windup.png|border|400px|Simulation of PI cruise control with windup. The figure shows the speed v and the throttle u for a car that encounters a slope that is so steep that the throttle saturates. The controller output is a dashed line. The controller parameters are kp = 0.5 and ki = 0.1.]]
| width=10% |
| width=45% |
[[Image:cruise-antiwindup.png|border|400px|Simulation of PI cruise control with anti-windup. The figure shows the speed v and the throttle u for a car that encounters a slope that is so steep that the throttle saturates. The controller output is a dashed line. The controller parameters are kp = 0.5, ki = 0.1. and kaw =2. The anti-windup compensator eliminates the overshoot by preventing the error from building up in the integral term of the controller.]]
|}

=== Python code ===

{{code block|Redefine input with longer duration}}
{{Code start}}
T = np.linspace(0, 50, 101)
vref = 20 * np.ones(T.shape)
theta_hill = [
0 if t <= 5 else
6./180. * pi * (t-5) if t <= 6 else
6./180. * pi for t in T]

# Compute the equilibrium throttle setting for the desired speed
X0, U0 = ct.find_eqpt(
cruise_nonlin, [vref[0], 0], [vref[0], gear[0], theta0[0]],
y0=[0, vref[0]], iu=[1, 2], iy=[1])
{{Code end}}

{{Code block|Effects of windup (plot)}}
{{Code start}}
t, y = ct.input_output_response(
cruise_nonlin, T, [vref, gear, theta_hill], X0,
params={'kaw':0})
cruise_plot(cruise_nonlin, t, y, antiwindup=True);

plt.suptitle('Cruise control with integrator windup')
plt.savefig("cruise-windup.png", bbox_inches='tight')
{{Code end}}

{{Code block|Anti-windup compensation (plot)}}
{{Code start}}
t, y = ct.input_output_response(
cruise_nonlin, T, [vref, gear, theta_hill], X0,
params={'kaw':2})
cruise_plot(cruise_nonlin, t, y, antiwindup=True);

plt.suptitle('Cruise control with anti-windup')
plt.savefig("cruise-antiwindup.png", bbox_inches='tight')
{{Code end}}

== Further Reading ==
* How Stuff Works: [http://auto.howstuffworks.com/cruise-control.htm cruise control]
* Wikipedia: [[Wikipedia:Cruise control|cruise control]]

[[Category:Running examples]]

Cruise control

2024-11-24T16:39:36Z

Murray:

Frequency Domain Design

2024-11-24T16:37:12Z

Murray: Undo revision 1324 by Murray (talk)

{{Chapter
|Chapter number=12
|Short name=loopsyn
|Previous chapter=PID Control
|Next chapter=Robust Performance
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/index.php?title=Frequency_Domain_Design#Frequently_Asked_Questions
|Chapter summary=In this chapter we continue to explore the use of frequency domain techniques with a focus on the design of feedback systems. We begin with a more thorough description of the performance specifications for control systems and then introduce the concept of “loop shaping” as a mechanism for designing controllers in the frequency domain. Additional techniques discussed in this chapter include feedforward compensation, the root locus method, and nested controller design.
|Chapter contents=# Sensitivity Functions
# Performance Specifications
#* Response to Reference Signals
#* Response to Load Disturbances and Measurement Noise
#* Measuring Specifications
# Feedback Design via Loop Shaping
#* Design Considerations
#* Lead and Lag Compensation
# Feedforward Design
#* Combining Feedforward and Feedback
#* Difficulties with Feedforward
#* Approximate Inverses
# The Root Locus Method
# Design Example
# Further Reading
:: Exercises
}}

PID Control

2024-11-24T16:36:59Z

Murray: Undo revision 1323 by Murray (talk)

{{Chapter
|Chapter number=11
|Short name=pid
|Previous chapter=Frequency Domain Analysis
|Next chapter=Frequency Domain Design
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/index.php?title=PID_Control#Frequently_Asked_Questions
|Chapter summary=Proportional-integral-derivative (PID) control is by far the most common way of using feedback in engineering systems. In this chapter we present the basic properties of PID control and the methods for choosing the parameters of the controllers. We also analyze the effects of actuator saturation, an important feature of many feedback systems, and describe methods for compensating for it. Finally, we discuss the implementation of PID controllers as an example of how to implement feedback control systems using analog or digital computation.
|Chapter contents=# Basic Control Functions
# Simple Controllers for Complex Systems
# PID Tuning
#* Ziegler--Nichols' Tuning
#* Tuning Based on the FOTD Model
#* Relay Feedback
# Integral Windup
#* Avoiding Windup
#* Manual Control and Tracking
#* Anti-Windup for General Controllers
# Implementation
#* Filtering the Derivative
#* Setpoint Weighting
#* Implementation Based on Operational Amplifiers
#* Computer Implementation
# Further Reading
:: Exercises
}}

Frequency Domain Analysis

2024-11-24T16:36:49Z

Murray: Undo revision 1322 by Murray (talk)

{{Chapter
|Chapter number=10
|Short name=loopanal
|Previous chapter=Transfer Functions
|Next chapter=PID Control
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/index.php?title=Frequency_Domain_Analysis#Frequently_Asked_Questions
|Chapter summary=In this chapter we study how the stability and robustness of closed loop systems can be determined by investigating how sinusoidal signals of different frequencies propagate around the feedback loop. This technique allows us to reason about the closed loop behavior of a system through the frequency domain properties of the open loop transfer function. The Nyquist stability theorem is a key result that provides a way to analyze stability and introduce measures of degrees of stability.
|Chapter contents=# The Loop Transfer Function
# The Nyquist Criterion
#* The Nyquist Plot
#* The General Nyquist Criterion
#* Conditional Stability
# Stability Margins
# Bode's Relations and Minimum Phase Systems
# Generalized Notions of Gain and Phase
#* System Gain and Passivity
#* Extensions of the Nyquist Criterion
#* Describing Functions
# Further Reading
:: Exercises
}}

Transfer Functions

2024-11-24T16:36:39Z

Murray: Undo revision 1321 by Murray (talk)

{{Chapter
|Chapter number=9
|Short name=xferfcns
|Previous chapter=Output Feedback
|Next chapter=Frequency Domain Analysis
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/index.php?title=Transfer_Functions#Frequently_Asked_Questions
|Chapter summary=This chapter introduces the concept of the transfer function, which is a com- pact description of the input/output relation for a linear time-invariant system. We show how to obtain transfer functions analytically and experimentally. Combining transfer functions with block diagrams gives a powerful algebraic method to analyze linear systems with many blocks. The transfer function allows new interpretations of system dynamics. We also introduce the Bode plot, a powerful graphical rep- resentation of the transfer function that was introduced by Bode to analyze and design feedback amplifiers.
|Chapter contents=# The Loop Transfer Function
# The Nyquist Criterion
#* The Nyquist Plot
#* The General Nyquist Criterion
#* Conditional Stability
# Stability Margins
# Bode's Relations and Minimum Phase Systems
# Generalized Notions of Gain and Phase
#* System Gain and Passivity
#* Extensions of the Nyquist Criterion
#* Describing Functions
# Further Reading
:: Exercises
}}

Frequency Domain Design

2024-11-24T16:31:05Z

Murray:

{{Chapter
|Chapter number=12
|Short name=loopsyn
|Previous chapter=PID Control
|Next chapter=Robust Performance
|First edition URL=https://www.cds.caltech.edu/~murray/amwiki/index.php?title=Frequency_Domain_Design#Frequently_Asked_Questions
|Chapter summary=In this chapter we continue to explore the use of frequency domain techniques with a focus on the design of feedback systems. We begin with a more thorough description of the performance specifications for control systems and then introduce the concept of “loop shaping” as a mechanism for designing controllers in the frequency domain. Additional techniques discussed in this chapter include feedforward compensation, the root locus method, and nested controller design.
|Chapter contents=# Sensitivity Functions
# Performance Specifications
# Feedback Design via Loop Shaping
# Feedforward Design
# The Root Locus Method
# Design Example
# Further Reading
:: Exercises
}}