Ordinary Differential Equations (16): Fundamentals of Control Theory
Chen Kai BOSS

When you drive a car, you constantly adjust the steering wheel to keep the vehicle centered in the lane; an air conditioner automatically regulates cooling power based on temperature feedback; a rocket precisely controls thrust to maintain orbit. These seemingly different systems share a common mathematical foundation — control theory. In this chapter, we explore how differential equations describe and design control systems, from classical PID controllers to modern state-space methods, seeing how mathematics helps us tame complex dynamical systems.

Basic Concepts of Control Theory

Open-Loop vs Closed-Loop Control

Open-loop control: The control signal does not depend on system output. - Example: A microwave heats for a set time - Advantage: Simple - Disadvantage: Cannot handle disturbances

Closed-loop control (feedback control): The control signal adjusts based on system output. - Example: Air conditioner adjusts based on room temperature - Advantage: Can resist disturbances, automatically corrects errors - Disadvantage: May become unstable

Basic Block Diagram of Control Systems

1
2
3
4
5
6
        +-----+     e(t)    +------------+    u(t)    +---------+    y(t)
r(t) -->| + |--(-)-------->| Controller |---------->| Plant |-------+-->
+-----+ +------------+ +---------+ |
^ |
| Feedback |
+-----------------------------------------------------------+
  • : Reference input (setpoint) -: System output -: Error -: Control input
  • Plant: The controlled object

Differential Equation Description

Most physical systems can be described by ordinary differential equations. For example, a mass-spring-damper system:whereis the applied force (control input) andis the displacement (output).

Transfer Functions and Laplace Transform

Laplace Transform ReviewExtra close brace or missing open brace\mathcal{L}\{f(t)} = F(s) = \int_0^\infty f(t)e^{-st}dt

Key properties: - -This transforms differential equations into algebraic equations!

Transfer Functions

For linear time-invariant systems, the transfer function is defined as (with zero initial conditions):

Example: Mass-spring-damper system

First-Order SystemsTransfer function:

-: DC gain -: Time constant

Step response:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.signal import lti, step

# First-order system
K = 2.0 # DC gain
tau = 1.0 # Time constant

# Method 1: Transfer function
num = [K]
den = [tau, 1]
system = lti(num, den)
t, y_tf = step(system)

# Method 2: Differential equation
def first_order_system(y, t, K, tau, u):
dydt = (K*u - y) / tau
return dydt

t_ode = np.linspace(0, 8, 200)
y_ode = odeint(first_order_system, 0, t_ode, args=(K, tau, 1))

# Analytical solution
y_analytical = K * (1 - np.exp(-t_ode/tau))

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

ax1 = axes[0]
ax1.plot(t, y_tf, 'b-', linewidth=2, label='Transfer Function')
ax1.plot(t_ode, y_ode, 'r--', linewidth=2, label='ODE Solution')
ax1.plot(t_ode, y_analytical, 'g:', linewidth=2, label='Analytical')
ax1.axhline(y=K, color='k', linestyle='--', alpha=0.5, label=f'Final value = {K}')
ax1.axhline(y=K*0.632, color='orange', linestyle=':', alpha=0.5)
ax1.axvline(x=tau, color='orange', linestyle=':', alpha=0.5)
ax1.plot(tau, K*0.632, 'ro', markersize=10)
ax1.text(tau+0.1, K*0.632, f'τ = {tau}', fontsize=10)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Output y(t)', fontsize=12)
ax1.set_title('First-Order System Step Response', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Comparison of different time constants
ax2 = axes[1]
tau_values = [0.5, 1.0, 2.0, 4.0]
colors = plt.cm.viridis(np.linspace(0, 0.8, len(tau_values)))

for tau_val, color in zip(tau_values, colors):
y = K * (1 - np.exp(-t_ode/tau_val))
ax2.plot(t_ode, y, color=color, linewidth=2, label=f'τ = {tau_val}')

ax2.axhline(y=K, color='k', linestyle='--', alpha=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Output y(t)', fontsize=12)
ax2.set_title('Effect of Time Constant', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('first_order_system.png', dpi=150, bbox_inches='tight')
plt.show()

Second-Order Systems

Standard form:Transfer function:

Parameter meanings: -: Natural frequency -: Damping ratio -: DC gain

Poles:

Effect of damping ratio: -: Undamped (sustained oscillation) -: Underdamped (decaying oscillation) -: Critically damped -: Overdamped

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
from scipy.signal import lti, step

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

omega_n = 1.0
K = 1.0
t = np.linspace(0, 15, 500)

# Different damping ratios
ax1 = axes[0]
zeta_values = [0.1, 0.3, 0.5, 0.7, 1.0, 2.0]
colors = plt.cm.coolwarm(np.linspace(0, 1, len(zeta_values)))

for zeta, color in zip(zeta_values, colors):
num = [K * omega_n**2]
den = [1, 2*zeta*omega_n, omega_n**2]
system = lti(num, den)
t_out, y = step(system, T=t)
ax1.plot(t_out, y, color=color, linewidth=2, label=f'ζ = {zeta}')

ax1.axhline(y=K, color='k', linestyle='--', alpha=0.5)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Output y(t)', fontsize=12)
ax1.set_title('Second-Order System: Effect of Damping Ratio', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Pole locations
ax2 = axes[1]
theta = np.linspace(0, np.pi, 100)
ax2.plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3)
ax2.plot(np.cos(theta), -np.sin(theta), 'k--', alpha=0.3)

for zeta, color in zip(zeta_values, colors):
if zeta < 1:
real_part = -zeta * omega_n
imag_part = omega_n * np.sqrt(1 - zeta**2)
ax2.plot(real_part, imag_part, 'o', color=color, markersize=10)
ax2.plot(real_part, -imag_part, 'o', color=color, markersize=10)
elif zeta == 1:
ax2.plot(-omega_n, 0, 'o', color=color, markersize=10)
else:
s1 = -zeta*omega_n + omega_n*np.sqrt(zeta**2 - 1)
s2 = -zeta*omega_n - omega_n*np.sqrt(zeta**2 - 1)
ax2.plot(s1, 0, 'o', color=color, markersize=10)
ax2.plot(s2, 0, 'o', color=color, markersize=10)

ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)
ax2.fill_between([-3, 0], [-3, -3], [3, 3], alpha=0.1, color='green', label='Stable region')
ax2.set_xlabel('Real Part', fontsize=12)
ax2.set_ylabel('Imaginary Part', fontsize=12)
ax2.set_title('Pole Locations in s-plane', fontsize=14, fontweight='bold')
ax2.set_xlim(-3, 1)
ax2.set_ylim(-2, 2)
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal')

plt.tight_layout()
plt.savefig('second_order_system.png', dpi=150, bbox_inches='tight')
plt.show()

PID Controllers

Basic Principle

The PID controller is the most widely used controller in industry, consisting of three components: - P (Proportional): Proportional to error, provides fast response but has steady-state error - I (Integral): Eliminates steady-state error, but may cause oscillation - D (Derivative): Predicts error trend, increases stability, but sensitive to noise

Transfer Function Form

Role of Each Component

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def pid_simulation():
"""PID controller simulation"""

# Plant: second-order system
def plant(y, t, u):
# y = [x, x_dot]
x, x_dot = y
m, c, k = 1.0, 0.5, 1.0 # mass, damping, spring constant
x_ddot = (u - c*x_dot - k*x) / m
return [x_dot, x_ddot]

# PID controller
class PIDController:
def __init__(self, Kp, Ki, Kd):
self.Kp = Kp
self.Ki = Ki
self.Kd = Kd
self.integral = 0
self.prev_error = 0

def compute(self, error, dt):
self.integral += error * dt
derivative = (error - self.prev_error) / dt if dt > 0 else 0
self.prev_error = error

u = self.Kp * error + self.Ki * self.integral + self.Kd * derivative
return u

def reset(self):
self.integral = 0
self.prev_error = 0

# Simulation parameters
dt = 0.01
t_final = 20
t = np.arange(0, t_final, dt)
setpoint = 1.0 # Target position

# Different controller configurations
configs = [
('P only', 2.0, 0, 0),
('PI', 2.0, 1.0, 0),
('PD', 2.0, 0, 1.0),
('PID', 2.0, 1.0, 1.0),
]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

for ax, (name, Kp, Ki, Kd) in zip(axes.flat, configs):
pid = PIDController(Kp, Ki, Kd)

y = [0, 0] # Initial state
y_history = [y[0]]
u_history = []

for i in range(1, len(t)):
error = setpoint - y[0]
u = pid.compute(error, dt)
u_history.append(u)

# Integrate one step
from scipy.integrate import odeint
y_new = odeint(plant, y, [0, dt], args=(u,))[-1]
y = list(y_new)
y_history.append(y[0])

ax.plot(t, y_history, 'b-', linewidth=2, label='Output')
ax.axhline(y=setpoint, color='r', linestyle='--', linewidth=2, label='Setpoint')
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Position', fontsize=12)
ax.set_title(f'{name}: Kp={Kp}, Ki={Ki}, Kd={Kd}', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(-0.5, 2)

plt.tight_layout()
plt.savefig('pid_control.png', dpi=150, bbox_inches='tight')
plt.show()

pid_simulation()

PID Parameter Tuning

Ziegler-Nichols Method:

  1. Set, gradually increaseuntil the system just sustains oscillation
  2. Record this(ultimate gain) and oscillation period3. Set PID parameters according to the following table:
Controller Type
P 0 0
PI 0
PID
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def ziegler_nichols_tuning():
"""Ziegler-Nichols tuning demonstration"""

# Plant
def plant_dynamics(state, t, u, m=1.0, c=0.2, k=1.0):
x, v = state
return [v, (u - c*v - k*x)/m]

# Find critical gain
def find_critical_gain():
Kp_test = np.linspace(0.1, 10, 100)
oscillation_detected = []

for Kp in Kp_test:
dt = 0.01
t = np.arange(0, 50, dt)
y = [0.1, 0] # Initial disturbance
y_history = []

for _ in range(len(t)):
error = 0 - y[0] # Target is 0
u = Kp * error
y_new = odeint(plant_dynamics, y, [0, dt], args=(u,))[-1]
y = list(y_new)
y_history.append(y[0])

y_history = np.array(y_history)

# Check for oscillation
if len(y_history) > 100:
late_signal = y_history[-500:]
amplitude_variation = np.max(late_signal) - np.min(late_signal)
if amplitude_variation > 0.01: # Still has amplitude
oscillation_detected.append((Kp, amplitude_variation))

return oscillation_detected

# Demonstrate response for different Kp
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

Kp_values = [1.0, 3.0, 5.0, 8.0]
dt = 0.01
t = np.arange(0, 30, dt)

for ax, Kp in zip(axes.flat, Kp_values):
y = [0.5, 0]
y_history = []

for _ in range(len(t)):
error = 0 - y[0]
u = Kp * error
y_new = odeint(plant_dynamics, y, [0, dt], args=(u,))[-1]
y = list(y_new)
y_history.append(y[0])

ax.plot(t, y_history, 'b-', linewidth=2)
ax.axhline(y=0, color='r', linestyle='--', alpha=0.5)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Position', fontsize=12)
ax.set_title(f'Kp = {Kp}', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)

plt.suptitle('Finding Critical Gain (Ziegler-Nichols)', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('ziegler_nichols.png', dpi=150, bbox_inches='tight')
plt.show()

ziegler_nichols_tuning()

State-Space Methods

State-Space Representation

Any linear time-invariant system can be written as: where: -: State vector -: Input vector -: Output vector -: System matrix -: Input matrix -: Output matrix -: Feedforward matrix

From Differential Equations to State Space

Example:Define state:

Stability Analysis

The system is stable if and only if all eigenvalues ofhave negative real parts.

Lyapunov stability: If there exists a positive definite matrixsuch that(positive definite), then the system is stable.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
import numpy as np
from scipy.linalg import solve_continuous_lyapunov

def state_space_analysis():
"""State space analysis"""

# System parameters
m, c, k = 1.0, 0.5, 2.0

# State-space matrices
A = np.array([[0, 1],
[-k/m, -c/m]])
B = np.array([[0], [1/m]])
C = np.array([[1, 0]])
D = np.array([[0]])

# Eigenvalue analysis
eigenvalues = np.linalg.eigvals(A)
print("System matrix A:")
print(A)
print(f"\nEigenvalues: {eigenvalues}")
print(f"System stable: {all(ev.real < 0 for ev in eigenvalues)}")

# Lyapunov analysis
Q = np.eye(2)
P = solve_continuous_lyapunov(A.T, -Q)
print(f"\nLyapunov matrix P:")
print(P)
print(f"P positive definite: {all(np.linalg.eigvals(P) > 0)}")

# Simulation
def state_space_ode(x, t, A, B, u):
return A @ x + B.flatten() * u

t = np.linspace(0, 15, 500)
x0 = [1, 0] # Initial conditions

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Free response
ax1 = axes[0, 0]
sol = odeint(state_space_ode, x0, t, args=(A, B, 0))
ax1.plot(t, sol[:, 0], 'b-', linewidth=2, label='Position x')
ax1.plot(t, sol[:, 1], 'r-', linewidth=2, label='Velocity v')
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('State', fontsize=12)
ax1.set_title('Free Response (u = 0)', fontsize=12, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# Phase portrait
ax2 = axes[0, 1]
# Multiple initial conditions
for x0_val in [(1, 0), (0, 1), (-1, 0), (0, -1), (0.5, 0.5)]:
sol = odeint(state_space_ode, x0_val, t, args=(A, B, 0))
ax2.plot(sol[:, 0], sol[:, 1], linewidth=1.5, alpha=0.7)
ax2.plot(x0_val[0], x0_val[1], 'ko', markersize=6)

ax2.plot(0, 0, 'r*', markersize=15, label='Equilibrium')
ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('Velocity v', fontsize=12)
ax2.set_title('Phase Portrait', fontsize=12, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal')

# Step response
ax3 = axes[1, 0]

# Simulate with step input
x = np.array([0.0, 0.0])
x_history = [x.copy()]
dt = t[1] - t[0]

for _ in range(len(t)-1):
u = 1.0 # Step input
x_dot = A @ x + B.flatten() * u
x = x + x_dot * dt
x_history.append(x.copy())

x_history = np.array(x_history)
ax3.plot(t, x_history[:, 0], 'b-', linewidth=2)
ax3.axhline(y=1/k, color='r', linestyle='--', alpha=0.5, label=f'Steady state = {1/k:.2f}')
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Position', fontsize=12)
ax3.set_title('Step Response', fontsize=12, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# Pole locations
ax4 = axes[1, 1]
ax4.plot(eigenvalues.real, eigenvalues.imag, 'rx', markersize=15, mew=3)
ax4.axhline(y=0, color='k', linewidth=0.5)
ax4.axvline(x=0, color='k', linewidth=0.5)
ax4.fill_between([-2, 0], [-2, -2], [2, 2], alpha=0.1, color='green', label='Stable region')
ax4.set_xlabel('Real Part', fontsize=12)
ax4.set_ylabel('Imaginary Part', fontsize=12)
ax4.set_title('Pole Locations', fontsize=12, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(-2, 1)
ax4.set_ylim(-2, 2)

plt.tight_layout()
plt.savefig('state_space.png', dpi=150, bbox_inches='tight')
plt.show()

state_space_analysis()

Controllability and Observability

Controllability

Definition: A system is controllable if it can be transferred from any initial state to any target state through appropriate input.

Controllability matrix:

Controllability condition: The system is controllable if and only if.

Observability

Definition: A system is observable if the initial state can be determined from the outputand input.

Observability matrix:

Observability condition: The system is observable if and only if.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def controllability_observability():
"""Controllability and observability analysis"""

# System 1: Fully controllable and observable
A1 = np.array([[0, 1], [-2, -3]])
B1 = np.array([[0], [1]])
C1 = np.array([[1, 0]])

# System 2: Not controllable
A2 = np.array([[1, 0], [0, 2]])
B2 = np.array([[1], [0]])
C2 = np.array([[1, 1]])

# System 3: Not observable
A3 = np.array([[1, 1], [0, 2]])
B3 = np.array([[1], [1]])
C3 = np.array([[1, 0]])

def analyze_system(A, B, C, name):
n = A.shape[0]

# Controllability matrix
controllability = B
Ai_B = B
for i in range(1, n):
Ai_B = A @ Ai_B
controllability = np.hstack([controllability, Ai_B])

ctrl_rank = np.linalg.matrix_rank(controllability)

# Observability matrix
observability = C
C_Ai = C
for i in range(1, n):
C_Ai = C_Ai @ A
observability = np.vstack([observability, C_Ai])

obs_rank = np.linalg.matrix_rank(observability)

print(f"\n{name}:")
print(f" A = {A.tolist()}")
print(f" B = {B.T.tolist()}")
print(f" C = {C.tolist()}")
print(f" Controllability matrix rank: {ctrl_rank}/{n} -> {'Controllable' if ctrl_rank == n else 'Not controllable'}")
print(f" Observability matrix rank: {obs_rank}/{n} -> {'Observable' if obs_rank == n else 'Not observable'}")

return ctrl_rank == n, obs_rank == n

systems = [
(A1, B1, C1, "System 1 (Fully controllable and observable)"),
(A2, B2, C2, "System 2 (Not controllable)"),
(A3, B3, C3, "System 3 (Not observable)"),
]

for A, B, C, name in systems:
analyze_system(A, B, C, name)

controllability_observability()

State Feedback and Pole Placement

State Feedback

If all states are measurable, state feedback can be used:Closed-loop system:

Pole Placement Theorem

Theorem: Ifis controllable, there exists a feedback gainsuch that the eigenvalues ofcan be arbitrarily assigned.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
from scipy.signal import place_poles

def pole_placement_demo():
"""Pole placement demonstration"""

# Unstable system
A = np.array([[0, 1], [2, 0]]) # Has positive eigenvalue
B = np.array([[0], [1]])
C = np.array([[1, 0]])

# Original eigenvalues
orig_poles = np.linalg.eigvals(A)
print(f"Original poles: {orig_poles}")
print(f"System unstable (has positive real part poles)")

# Desired poles
desired_poles = np.array([-2, -3])

# Compute feedback gain
result = place_poles(A, B, desired_poles)
K = result.gain_matrix
print(f"\nDesired poles: {desired_poles}")
print(f"Feedback gain K = {K}")

# Verify
A_cl = A - B @ K
achieved_poles = np.linalg.eigvals(A_cl)
print(f"Achieved poles: {achieved_poles}")

# Simulation comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

t = np.linspace(0, 5, 500)
x0 = [1, 0]

# Open-loop response
def open_loop(x, t):
return A @ x

sol_open = odeint(open_loop, x0, t)

# Closed-loop response
def closed_loop(x, t):
return (A - B @ K) @ x

sol_closed = odeint(closed_loop, x0, t)

ax1 = axes[0]
ax1.plot(t, sol_open[:, 0], 'r-', linewidth=2, label='Open-loop (unstable)')
ax1.plot(t, sol_closed[:, 0], 'b-', linewidth=2, label='Closed-loop (stabilized)')
ax1.axhline(y=0, color='k', linestyle='--', alpha=0.3)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Position x', fontsize=12)
ax1.set_title('Stabilization via State Feedback', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-5, 15)

# Pole locations
ax2 = axes[1]
ax2.plot(orig_poles.real, orig_poles.imag, 'rx', markersize=15, mew=3, label='Open-loop poles')
ax2.plot(achieved_poles.real, achieved_poles.imag, 'bo', markersize=12, label='Closed-loop poles')
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)
ax2.fill_between([-4, 0], [-3, -3], [3, 3], alpha=0.1, color='green', label='Stable region')
ax2.set_xlabel('Real Part', fontsize=12)
ax2.set_ylabel('Imaginary Part', fontsize=12)
ax2.set_title('Pole Locations', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-4, 3)
ax2.set_ylim(-3, 3)

plt.tight_layout()
plt.savefig('pole_placement.png', dpi=150, bbox_inches='tight')
plt.show()

pole_placement_demo()

LQR Optimal Control

Problem Formulation

The Linear Quadratic Regulator (LQR) finds the control law that minimizes the following cost function: -: State weight matrix (positive semi-definite) -: Control weight matrix (positive definite)

Optimal Solution

The optimal control law is state feedback:where, andis the solution to the algebraic Riccati equation:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
from scipy.linalg import solve_continuous_are

def lqr_demo():
"""LQR controller demonstration"""

# Linearized inverted pendulum model
M = 1.0 # Cart mass
m = 0.1 # Pendulum mass
l = 0.5 # Pendulum length
g = 9.81 # Gravitational acceleration

# State: [x, x_dot, theta, theta_dot]
A = np.array([
[0, 1, 0, 0],
[0, 0, -m*g/M, 0],
[0, 0, 0, 1],
[0, 0, (M+m)*g/(M*l), 0]
])

B = np.array([[0], [1/M], [0], [-1/(M*l)]])

# Open-loop eigenvalues
open_poles = np.linalg.eigvals(A)
print(f"Open-loop poles: {open_poles}")
print(f"System unstable: {any(p.real > 0 for p in open_poles)}")

# LQR design
Q = np.diag([10, 1, 100, 10]) # State weights
R = np.array([[1]]) # Control weight

# Solve Riccati equation
P = solve_continuous_are(A, B, Q, R)
K = np.linalg.inv(R) @ B.T @ P

print(f"\nLQR gain K = {K.flatten()}")

# Closed-loop poles
A_cl = A - B @ K
closed_poles = np.linalg.eigvals(A_cl)
print(f"Closed-loop poles: {closed_poles}")

# Simulation
def inverted_pendulum(x, t, K, A, B):
u = -K @ x
return (A @ x + B @ u.T).flatten()

t = np.linspace(0, 10, 1000)
x0 = [0, 0, 0.2, 0] # Initial angle deviation

sol = odeint(inverted_pendulum, x0, t, args=(K, A, B))
u_history = [-K @ sol[i] for i in range(len(t))]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Position
ax1 = axes[0, 0]
ax1.plot(t, sol[:, 0], 'b-', linewidth=2)
ax1.set_xlabel('Time (s)', fontsize=12)
ax1.set_ylabel('Cart Position (m)', fontsize=12)
ax1.set_title('Cart Position', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# Angle
ax2 = axes[0, 1]
ax2.plot(t, sol[:, 2]*180/np.pi, 'r-', linewidth=2)
ax2.set_xlabel('Time (s)', fontsize=12)
ax2.set_ylabel('Pendulum Angle (deg)', fontsize=12)
ax2.set_title('Pendulum Angle', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Control force
ax3 = axes[1, 0]
ax3.plot(t, [u[0][0] for u in u_history], 'g-', linewidth=2)
ax3.set_xlabel('Time (s)', fontsize=12)
ax3.set_ylabel('Control Force (N)', fontsize=12)
ax3.set_title('Control Input', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# Pole comparison
ax4 = axes[1, 1]
ax4.plot(open_poles.real, open_poles.imag, 'rx', markersize=15, mew=3, label='Open-loop')
ax4.plot(closed_poles.real, closed_poles.imag, 'bo', markersize=12, label='Closed-loop (LQR)')
ax4.axhline(y=0, color='k', linewidth=0.5)
ax4.axvline(x=0, color='k', linewidth=0.5)
ax4.fill_between([-10, 0], [-10, -10], [10, 10], alpha=0.1, color='green')
ax4.set_xlabel('Real Part', fontsize=12)
ax4.set_ylabel('Imaginary Part', fontsize=12)
ax4.set_title('Pole Locations', fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(-10, 5)
ax4.set_ylim(-10, 10)

plt.suptitle('LQR Control of Inverted Pendulum', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('lqr_control.png', dpi=150, bbox_inches='tight')
plt.show()

lqr_demo()

Observer Design

Why Do We Need Observers?

In practical systems, not all states can be directly measured. An observer (or state estimator) reconstructs states from measurable inputs and outputs.

Luenberger Observerwhereis the estimated state andis the observer gain.

Estimation error dynamics:By choosing, the error convergence rate can be arbitrarily assigned.

Separation Principle

Theorem: For linear systems, the state feedback controller and observer can be designed independently. The closed-loop poles are the union of controller poles and observer poles.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
def observer_demo():
"""Observer design demonstration"""

# System
A = np.array([[0, 1], [-2, -3]])
B = np.array([[0], [1]])
C = np.array([[1, 0]])

# Controller gain (via pole placement)
desired_ctrl_poles = np.array([-2, -3])
result_ctrl = place_poles(A, B, desired_ctrl_poles)
K = result_ctrl.gain_matrix

# Observer gain (typically 2-5 times faster than controller)
desired_obs_poles = np.array([-10, -12])
# Dual system pole placement
result_obs = place_poles(A.T, C.T, desired_obs_poles)
L = result_obs.gain_matrix.T

print(f"Controller gain K = {K}")
print(f"Observer gain L = {L.flatten()}")

# Full system simulation (true state + estimated state)
def full_system(state, t, A, B, C, K, L, r):
x = state[:2] # True state
x_hat = state[2:] # Estimated state

y = C @ x
y_hat = C @ x_hat

u = -K @ x_hat + r # Feedback using estimated state

# True system
x_dot = A @ x + B @ u.flatten()

# Observer
x_hat_dot = A @ x_hat + B @ u.flatten() + L @ (y - y_hat)

return np.concatenate([x_dot, x_hat_dot])

t = np.linspace(0, 5, 500)
x0 = [1, 0] # True initial state
x_hat0 = [0, 0] # Estimated initial state (unknown true state)
state0 = x0 + x_hat0

sol = odeint(full_system, state0, t, args=(A, B, C, K, L, 0))

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# State 1
ax1 = axes[0, 0]
ax1.plot(t, sol[:, 0], 'b-', linewidth=2, label='True x ₁')
ax1.plot(t, sol[:, 2], 'r--', linewidth=2, label='Estimated x ̂₁')
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('State 1', fontsize=12)
ax1.set_title('State 1: True vs Estimated', fontsize=12, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# State 2
ax2 = axes[0, 1]
ax2.plot(t, sol[:, 1], 'b-', linewidth=2, label='True x ₂')
ax2.plot(t, sol[:, 3], 'r--', linewidth=2, label='Estimated x ̂₂')
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('State 2', fontsize=12)
ax2.set_title('State 2: True vs Estimated', fontsize=12, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# Estimation error
ax3 = axes[1, 0]
e1 = sol[:, 0] - sol[:, 2]
e2 = sol[:, 1] - sol[:, 3]
ax3.plot(t, e1, 'b-', linewidth=2, label='Error in x ₁')
ax3.plot(t, e2, 'r-', linewidth=2, label='Error in x ₂')
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Estimation Error', fontsize=12)
ax3.set_title('Estimation Error (converges to 0)', fontsize=12, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# Pole locations
ax4 = axes[1, 1]
ctrl_poles = np.linalg.eigvals(A - B @ K)
obs_poles = np.linalg.eigvals(A - L @ C)

ax4.plot(ctrl_poles.real, ctrl_poles.imag, 'bo', markersize=12, label='Controller poles')
ax4.plot(obs_poles.real, obs_poles.imag, 'rs', markersize=12, label='Observer poles')
ax4.axhline(y=0, color='k', linewidth=0.5)
ax4.axvline(x=0, color='k', linewidth=0.5)
ax4.set_xlabel('Real Part', fontsize=12)
ax4.set_ylabel('Imaginary Part', fontsize=12)
ax4.set_title('Closed-Loop Poles', fontsize=12, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(-15, 2)
ax4.set_ylim(-5, 5)

plt.suptitle('Observer-Based Control (Separation Principle)', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('observer.png', dpi=150, bbox_inches='tight')
plt.show()

observer_demo()

Frequency Domain Analysis

Bode Plots

Bode plots show the frequency response of a system: - Magnitude plot:(dB) vs(log) - Phase plot:(degrees) vs(log)

Stability Margins

Gain margin: How much gain can be added before the system becomes unstable, measured when the phase is.

Phase margin: How much phase lag can be tolerated before instability, measured when the gain is 0 dB.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
from scipy import signal

def frequency_analysis():
"""Frequency domain analysis demonstration"""

# System transfer function
num = [1]
den = [1, 2, 1] # (s+1)^2
sys = signal.TransferFunction(num, den)

# Bode plot
w, mag, phase = signal.bode(sys, n=500)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Magnitude
ax1 = axes[0, 0]
ax1.semilogx(w, mag, 'b-', linewidth=2)
ax1.axhline(y=-3, color='r', linestyle='--', alpha=0.5, label='-3 dB')
ax1.set_xlabel('Frequency (rad/s)', fontsize=12)
ax1.set_ylabel('Magnitude (dB)', fontsize=12)
ax1.set_title('Bode Plot: Magnitude', fontsize=12, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# Phase
ax2 = axes[0, 1]
ax2.semilogx(w, phase, 'r-', linewidth=2)
ax2.set_xlabel('Frequency (rad/s)', fontsize=12)
ax2.set_ylabel('Phase (degrees)', fontsize=12)
ax2.set_title('Bode Plot: Phase', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Open-loop transfer function (with controller)
Kp = 10
num_ol = [Kp]
den_ol = [1, 2, 1, 0] # Kp / (s(s+1)^2)
sys_ol = signal.TransferFunction(num_ol, den_ol)

w_ol, mag_ol, phase_ol = signal.bode(sys_ol, n=500)

# Open-loop Bode plot
ax3 = axes[1, 0]
ax3.semilogx(w_ol, mag_ol, 'b-', linewidth=2)
ax3.axhline(y=0, color='k', linestyle='--', alpha=0.5, label='0 dB')
ax3.set_xlabel('Frequency (rad/s)', fontsize=12)
ax3.set_ylabel('Magnitude (dB)', fontsize=12)
ax3.set_title('Open-Loop Bode: Magnitude', fontsize=12, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
ax4.semilogx(w_ol, phase_ol, 'r-', linewidth=2)
ax4.axhline(y=-180, color='k', linestyle='--', alpha=0.5, label='-180°')
ax4.set_xlabel('Frequency (rad/s)', fontsize=12)
ax4.set_ylabel('Phase (degrees)', fontsize=12)
ax4.set_title('Open-Loop Bode: Phase', fontsize=12, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('frequency_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

frequency_analysis()

Application Example: Quadrotor Attitude Control

Simplified Model

Consider attitude control of a quadrotor along one axis:whereis the pitch angle andis the torque.

PID Control Design

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def quadrotor_attitude_control():
"""Quadrotor attitude control"""

I = 0.01 # Moment of inertia

def quadrotor(state, t, tau):
theta, omega = state
theta_ddot = tau / I
return [omega, theta_ddot]

# PID controller
class AttitudeController:
def __init__(self, Kp, Ki, Kd):
self.Kp = Kp
self.Ki = Ki
self.Kd = Kd
self.integral = 0
self.prev_error = 0

def compute(self, theta_des, theta, omega, dt):
error = theta_des - theta
self.integral += error * dt
derivative = -omega # Use angular velocity as derivative

tau = self.Kp * error + self.Ki * self.integral + self.Kd * derivative
return tau

# Simulation parameters
dt = 0.001
t_final = 5
t = np.arange(0, t_final, dt)

# Desired trajectory: step changes
theta_des = np.zeros_like(t)
theta_des[t > 1] = 0.5 # 30 degree step
theta_des[t > 3] = 0 # Return to level

# Different controllers
controllers = [
('Under-tuned', 0.1, 0, 0.01),
('Well-tuned', 0.5, 0.1, 0.05),
('Aggressive', 2.0, 0.5, 0.1),
]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

ax1 = axes[0, 0]
ax2 = axes[0, 1]
ax3 = axes[1, 0]

colors = ['blue', 'green', 'red']

for (name, Kp, Ki, Kd), color in zip(controllers, colors):
controller = AttitudeController(Kp, Ki, Kd)

state = [0, 0]
theta_history = [0]
tau_history = []

for i in range(1, len(t)):
tau = controller.compute(theta_des[i], state[0], state[1], dt)
tau = np.clip(tau, -0.1, 0.1) # Torque limit
tau_history.append(tau)

state_new = odeint(quadrotor, state, [0, dt], args=(tau,))[-1]
state = list(state_new)
theta_history.append(state[0])

ax1.plot(t, np.array(theta_history)*180/np.pi, color=color,
linewidth=2, label=name)
ax2.plot(t[1:], tau_history, color=color, linewidth=1.5,
label=name, alpha=0.8)

ax1.plot(t, theta_des*180/np.pi, 'k--', linewidth=2, label='Desired')
ax1.set_xlabel('Time (s)', fontsize=12)
ax1.set_ylabel('Pitch Angle (deg)', fontsize=12)
ax1.set_title('Attitude Response', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

ax2.set_xlabel('Time (s)', fontsize=12)
ax2.set_ylabel('Control Torque (Nm)', fontsize=12)
ax2.set_title('Control Input', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# Noise robustness test
ax3 = axes[1, 0]

controller = AttitudeController(0.5, 0.1, 0.05)
np.random.seed(42)

state = [0, 0]
theta_history = [0]

for i in range(1, len(t)):
# Add measurement noise
theta_measured = state[0] + 0.01 * np.random.randn()
omega_measured = state[1] + 0.1 * np.random.randn()

# Add disturbance torque
disturbance = 0.005 * np.sin(2*np.pi*0.5*t[i])

tau = controller.compute(theta_des[i], theta_measured, omega_measured, dt)
tau = np.clip(tau, -0.1, 0.1)

state_new = odeint(quadrotor, state, [0, dt], args=(tau + disturbance,))[-1]
state = list(state_new)
theta_history.append(state[0])

ax3.plot(t, np.array(theta_history)*180/np.pi, 'b-', linewidth=1.5, label='With noise & disturbance')
ax3.plot(t, theta_des*180/np.pi, 'k--', linewidth=2, label='Desired')
ax3.set_xlabel('Time (s)', fontsize=12)
ax3.set_ylabel('Pitch Angle (deg)', fontsize=12)
ax3.set_title('Robustness Test', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# Quadrotor schematic
ax4 = axes[1, 1]
ax4.set_xlim(-2, 2)
ax4.set_ylim(-1.5, 1.5)

# Draw body
from matplotlib.patches import Rectangle, FancyArrowPatch

body = Rectangle((-0.5, -0.1), 1.0, 0.2, fill=True, color='gray', alpha=0.7)
ax4.add_patch(body)

# Draw propellers
prop_positions = [(-0.8, 0), (0.8, 0)]
for x, y in prop_positions:
ax4.plot([x-0.3, x+0.3], [y+0.15, y+0.15], 'b-', linewidth=3)
ax4.arrow(x, y+0.2, 0, 0.3, head_width=0.1, head_length=0.05, fc='red', ec='red')

ax4.annotate('', xy=(0.5, 0.8), xytext=(0, 0.3),
arrowprops=dict(arrowstyle='->', color='green', lw=2))
ax4.text(0.6, 0.7, 'θ (pitch)', fontsize=12, color='green')

ax4.set_aspect('equal')
ax4.axis('off')
ax4.set_title('Quadrotor Attitude Control', fontsize=12, fontweight='bold')

plt.tight_layout()
plt.savefig('quadrotor_control.png', dpi=150, bbox_inches='tight')
plt.show()

quadrotor_attitude_control()

Summary

In this chapter, we covered the core content of control theory:

  1. Transfer functions: Transform differential equations into algebraic problems
  2. PID control: The most commonly used controller in industry
  3. State-space methods: Foundation of modern control theory
  4. Controllability and observability: Can the system be controlled and observed
  5. Pole placement: Arbitrarily assign closed-loop poles
  6. LQR optimal control: Optimal control minimizing a cost function
  7. Observer design: Reconstruct states from outputs
  8. Frequency domain analysis: Understand systems from frequency response perspective

Control theory elevates differential equations from a descriptive tool to a design tool — we can not only analyze systems but also actively design system behavior. From aircraft autopilots to smartphone image stabilization, control theory applications are ubiquitous.

Exercises

Basic Problems

  1. Find the step response of the first-order systemand plot the curve.

  2. For the second-order system:

    • Calculate the natural frequencyand damping ratio - Determine if the system is underdamped, critically damped, or overdamped
    • Calculate the step response overshoot
  3. Design a PID controller for the following system to achieve: overshoot < 10%, settling time < 2 seconds:4. Convert the following differential equation to state-space form:5. Determine the controllability and observability of the following system:

Advanced Problems

  1. Cascade control: Design a cascade PID controller for a temperature system:

    • Inner loop: Heater power control
    • Outer loop: Temperature control Compare performance of single-loop and cascade control.
  2. System identification: Given step response data, estimate the parametersandof a first-order system.

  3. Prove that for the LQR problem, ifis controllable andis observable, then the closed-loop system is asymptotically stable.

  4. Discretization: Discretize the continuous-time systemto. Discuss the choice of sampling period.

  5. Design a Luenberger observer with error convergence rate 3 times faster than the controller poles.

Programming Problems

  1. Implement a general PID controller class including:

    • Anti-windup
    • Derivative filtering
    • Output limiting
  2. Write a program for automatic Ziegler-Nichols parameter tuning.

  3. Implement a complete inverted pendulum simulator including:

    • Nonlinear dynamics
    • LQR controller
    • State observer
    • Visualization animation
  4. Train a neural network controller using reinforcement learning (e.g., PPO) and compare performance with LQR.

  5. Implement a digital filter design tool supporting Butterworth, Chebyshev, and other types.

Discussion Questions

  1. PID controllers are so simple — why are they still the most commonly used controllers in industry?

  2. Discuss the effect of model uncertainty on controller performance. How can robust controllers be designed?

  3. Compare the advantages and disadvantages of pole placement and LQR. When should each method be used?

  4. Why are observer poles typically designed to be faster than controller poles?

  5. Discuss the relationship between machine learning methods (such as deep reinforcement learning) and traditional control theory, and their respective application scenarios.

References

  1. Ogata, K. (2010). Modern Control Engineering. Prentice Hall.

  2. Franklin, G. F., Powell, J. D., & Emami-Naeini, A. (2015). Feedback Control of Dynamic Systems. Pearson.

  3. Å str ö m, K. J., & Murray, R. M. (2008). Feedback Systems: An Introduction for Scientists and Engineers. Princeton University Press.

  4. Dorf, R. C., & Bishop, R. H. (2017). Modern Control Systems. Pearson.

  5. Skogestad, S., & Postlethwaite, I. (2005). Multivariable Feedback Control: Analysis and Design. Wiley.

  6. Lewis, F. L., Vrabie, D., & Syrmos, V. L. (2012). Optimal Control. Wiley.

  7. Khalil, H. K. (2002). Nonlinear Systems. Prentice Hall.


Previous Chapter: ← Chapter 15: Population Dynamics

Next Chapter: → Chapter 17: Physics and Engineering Applications


This article is Chapter 16 of the "erta Ordinary Differential Equations" series.

  • Post title:Ordinary Differential Equations (16): Fundamentals of Control Theory
  • Post author:Chen Kai
  • Create time:2019-06-23 10:30:00
  • Post link:https://www.chenk.top/ode-chapter-16-control-theory/
  • Copyright Notice:All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments