常微分方程（十六）控制理论基础
Chen Kai BOSS
  2019-06-23 10:30:00 2019-06-23 10:30      6.8k 字  35 分钟  

当你驾驶汽车时，你会不断调整方向盘使车保持在车道中央；空调根据温度反馈自动调节制冷功率；火箭通过精确控制推力保持轨道。这些看似不同的系统有一个共同的数学基础——控制理论。本章我们将探索如何用微分方程描述和设计控制系统，从经典的 PID 控制器到现代状态空间方法，看看数学如何帮助我们驯服复杂的动态系统。

控制理论的基本概念

开环与闭环控制

开环控制：控制信号不依赖于系统输出。 - 例子：微波炉设定时间后加热 - 优点：简单 - 缺点：无法处理干扰

闭环控制（反馈控制）：控制信号根据系统输出进行调整。 - 例子：空调根据室温调节 - 优点：能抵抗干扰，自动修正误差 - 缺点：可能不稳定

控制系统的基本框图

1
2
3
4
5
6
        +-----+     e(t)    +------------+    u(t)    +---------+    y(t)
r(t) -->| + |--(-)-------->| Controller |---------->|  Plant  |-------+-->
        +-----+             +------------+           +---------+       |
           ^                                                           |
           |                     Feedback                              |
           +-----------------------------------------------------------+
  • ：参考输入（期望值）
  • ：系统输出
  • ：误差
  • ：控制输入
  • Plant：被控对象

微分方程描述

大多数物理系统可以用常微分方程描述。例如，一个质量-弹簧-阻尼系统：

$$

m + c + kx = u(t) $$

其中 是外加力（控制输入）， 是位移（输出）。

传递函数与拉普拉斯变换

拉普拉斯变换回顾

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

关键性质： - Extra close brace or missing open brace\mathcal{L}\{\dot{f}} = sF(s) - f(0) - Extra close brace or missing open brace\mathcal{L}\{\ddot{f}} = s^2F(s) - sf(0) - \dot{f}(0) 这使得微分方程变成代数方程！

传递函数

对于线性时不变系统，传递函数定义为（零初始条件）：

$$

G(s) = $$

例子：质量-弹簧-阻尼系统

$$

ms^2X(s) + csX(s) + kX(s) = U(s)

G(s) = = $$

一阶系统

传递函数： $$

G(s) = $$ - ：直流增益 - ：时间常数

阶跃响应： $$

y(t) = K(1 - e^{-t/}) $$ 

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

# 一阶系统
K = 2.0   # 直流增益
tau = 1.0  # 时间常数

# 方法 1：传递函数
num = [K]
den = [tau, 1]
system = lti(num, den)
t, y_tf = step(system)

# 方法 2：微分方程
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))

# 解析解
y_analytical = K * (1 - np.exp(-t_ode/tau))

# 可视化
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)

# 不同时间常数的比较
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()

二阶系统

标准形式：

传递函数： $$

G(s) = $$

参数意义： - ：自然频率 - ：阻尼比 - ：直流增益

极点： $$

s = -_n _n $$

阻尼比的影响： - ：无阻尼（持续振荡） - ：欠阻尼（衰减振荡） - ：临界阻尼 - ：过阻尼

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)

# 不同阻尼比
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)

# 极点位置
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 控制器

基本原理

PID 控制器是工业中最常用的控制器，由三个部分组成：

$$

u(t) = K_p e(t) + K_i _0^t e()d+ K_d $$ - P（比例）：与误差成正比，快速响应但有稳态误差 - I（积分）：消除稳态误差，但可能导致振荡 - D（微分）：预测误差趋势，增加稳定性，但对噪声敏感

传递函数形式

$$

G_c(s) = K_p + + K_d s = $$

各部分的作用

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 控制器仿真"""
    
    # 被控对象：二阶系统
    def plant(y, t, u):
        # y = [x, x_dot]
        x, x_dot = y
        m, c, k = 1.0, 0.5, 1.0  # 质量、阻尼、弹簧常数
        x_ddot = (u - c*x_dot - k*x) / m
        return [x_dot, x_ddot]
    
    # PID 控制器
    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
    
    # 模拟参数
    dt = 0.01
    t_final = 20
    t = np.arange(0, t_final, dt)
    setpoint = 1.0  # 目标位置
    
    # 不同控制器配置
    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]  # 初始状态
        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)
            
            # 积分一步
            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 参数整定

Ziegler-Nichols 方法：

  1. 设，逐渐增大 直到系统恰好持续振荡
  2. 记录此时的（临界增益）和振荡周期3. 根据下表设定 PID 参数：
控制器类型
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 整定演示"""
    
    # 被控对象
    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]
    
    # 找临界增益
    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]  # 初始扰动
            y_history = []
            
            for _ in range(len(t)):
                error = 0 - y[0]  # 目标是 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)
            
            # 检测是否振荡
            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:  # 仍有振幅
                    oscillation_detected.append((Kp, amplitude_variation))
        
        return oscillation_detected
    
    # 演示不同 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()

状态空间方法

状态空间表示

任何线性时不变系统可以写成：

其中： - ：状态向量 - ：输入向量 - ：输出向量 - ：系统矩阵 - ：输入矩阵 - ：输出矩阵 - ：前馈矩阵

从微分方程到状态空间

例子： 定义状态：

$$

y =

$$

稳定性分析

系统稳定当且仅当 的所有特征值有负实部。

Lyapunov 稳定性：如果存在正定矩阵 使得（ 正定），则系统稳定。

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():
    """状态空间分析"""
    
    # 系统参数
    m, c, k = 1.0, 0.5, 2.0
    
    # 状态空间矩阵
    A = np.array([[0, 1],
                  [-k/m, -c/m]])
    B = np.array([[0], [1/m]])
    C = np.array([[1, 0]])
    D = np.array([[0]])
    
    # 特征值分析
    eigenvalues = np.linalg.eigvals(A)
    print("系统矩阵 A:")
    print(A)
    print(f"\n 特征值: {eigenvalues}")
    print(f"系统稳定: {all(ev.real < 0 for ev in eigenvalues)}")
    
    # Lyapunov 分析
    Q = np.eye(2)
    P = solve_continuous_lyapunov(A.T, -Q)
    print(f"\nLyapunov 矩阵 P:")
    print(P)
    print(f"P 正定: {all(np.linalg.eigvals(P) > 0)}")
    
    # 仿真
    def state_space_ode(x, t, A, B, u):
        return A @ x + B.flatten() * u
    
    t = np.linspace(0, 15, 500)
    x0 = [1, 0]  # 初始条件
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 自由响应
    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)
    
    # 相图
    ax2 = axes[0, 1]
    # 多个初始条件
    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')
    
    # 阶跃响应
    ax3 = axes[1, 0]
    
    # 用阶跃输入仿真
    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  # 阶跃输入
        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)
    
    # 极点位置
    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()

可控性与可观性

可控性

定义：如果可以通过适当的输入 将系统从任意初始状态转移到任意目标状态，则系统是可控的。

可控性矩阵：

可控性条件：系统可控当且仅当。

可观性

定义：如果可以从输出 和输入 确定初始状态，则系统是可观的。

可观性矩阵：

可观性条件：系统可观当且仅当。

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():
    """可控性和可观性分析"""
    
    # 系统 1：完全可控可观
    A1 = np.array([[0, 1], [-2, -3]])
    B1 = np.array([[0], [1]])
    C1 = np.array([[1, 0]])
    
    # 系统 2：不可控
    A2 = np.array([[1, 0], [0, 2]])
    B2 = np.array([[1], [0]])
    C2 = np.array([[1, 1]])
    
    # 系统 3：不可观
    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 = 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 = 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"  可控性矩阵秩: {ctrl_rank}/{n} -> {'可控' if ctrl_rank == n else '不可控'}")
        print(f"  可观性矩阵秩: {obs_rank}/{n} -> {'可观' if obs_rank == n else '不可观'}")
        
        return ctrl_rank == n, obs_rank == n
    
    systems = [
        (A1, B1, C1, "系统 1（完全可控可观）"),
        (A2, B2, C2, "系统 2（不可控）"),
        (A3, B3, C3, "系统 3（不可观）"),
    ]
    
    for A, B, C, name in systems:
        analyze_system(A, B, C, name)

controllability_observability()

状态反馈与极点配置

状态反馈

如果所有状态可测量，可以使用状态反馈：

闭环系统：

极点配置定理

定理：如果 可控，则存在反馈增益 使得 的特征值可以被任意指定。

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():
    """极点配置演示"""
    
    # 不稳定系统
    A = np.array([[0, 1], [2, 0]])  # 有正特征值
    B = np.array([[0], [1]])
    C = np.array([[1, 0]])
    
    # 原始特征值
    orig_poles = np.linalg.eigvals(A)
    print(f"原始极点: {orig_poles}")
    print(f"系统不稳定（有正实部极点）")
    
    # 期望极点
    desired_poles = np.array([-2, -3])
    
    # 计算反馈增益
    result = place_poles(A, B, desired_poles)
    K = result.gain_matrix
    print(f"\n 期望极点: {desired_poles}")
    print(f"反馈增益 K = {K}")
    
    # 验证
    A_cl = A - B @ K
    achieved_poles = np.linalg.eigvals(A_cl)
    print(f"实现的极点: {achieved_poles}")
    
    # 仿真比较
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    t = np.linspace(0, 5, 500)
    x0 = [1, 0]
    
    # 开环响应
    def open_loop(x, t):
        return A @ x
    
    sol_open = odeint(open_loop, x0, t)
    
    # 闭环响应
    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)
    
    # 极点位置
    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 最优控制

问题设定

线性二次调节器（ LQR）寻找使以下代价函数最小的控制律：

$$

J = _0(T Q + ^T R ) dt $$ - ：状态权重矩阵（半正定） - ：控制权重矩阵（正定）

最优解

最优控制律是状态反馈：

其中， 是代数 Riccati 方程的解： $$

A^T P + PA - PBR{-1}BT P + Q = 0 $$ 

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 控制器演示"""
    
    # 倒立摆线性化模型
    M = 1.0   # 小车质量
    m = 0.1   # 摆杆质量
    l = 0.5   # 摆杆长度
    g = 9.81  # 重力加速度
    
    # 状态: [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_poles = np.linalg.eigvals(A)
    print(f"开环极点: {open_poles}")
    print(f"系统不稳定: {any(p.real > 0 for p in open_poles)}")
    
    # LQR 设计
    Q = np.diag([10, 1, 100, 10])  # 状态权重
    R = np.array([[1]])            # 控制权重
    
    # 求解 Riccati 方程
    P = solve_continuous_are(A, B, Q, R)
    K = np.linalg.inv(R) @ B.T @ P
    
    print(f"\nLQR 增益 K = {K.flatten()}")
    
    # 闭环极点
    A_cl = A - B @ K
    closed_poles = np.linalg.eigvals(A_cl)
    print(f"闭环极点: {closed_poles}")
    
    # 仿真
    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]  # 初始角度偏离
    
    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))
    
    # 位置
    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)
    
    # 角度
    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)
    
    # 控制力
    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)
    
    # 极点比较
    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()

观测器设计

为什么需要观测器？

在实际系统中，并非所有状态都可直接测量。观测器（或状态估计器）从可测量的输入和输出重构状态。

Luenberger 观测器

其中 是估计状态， 是观测器增益。

估计误差动态：

通过选择，可以任意配置误差收敛速度。

分离原理

定理：对于线性系统，可以独立设计状态反馈控制器和观测器。闭环极点是控制器极点和观测器极点的并集。

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():
    """观测器设计演示"""
    
    # 系统
    A = np.array([[0, 1], [-2, -3]])
    B = np.array([[0], [1]])
    C = np.array([[1, 0]])
    
    # 控制器增益（通过极点配置）
    desired_ctrl_poles = np.array([-2, -3])
    result_ctrl = place_poles(A, B, desired_ctrl_poles)
    K = result_ctrl.gain_matrix
    
    # 观测器增益（通常比控制器快 2-5 倍）
    desired_obs_poles = np.array([-10, -12])
    # 对偶系统极点配置
    result_obs = place_poles(A.T, C.T, desired_obs_poles)
    L = result_obs.gain_matrix.T
    
    print(f"控制器增益 K = {K}")
    print(f"观测器增益 L = {L.flatten()}")
    
    # 完整系统仿真（真实状态 + 估计状态）
    def full_system(state, t, A, B, C, K, L, r):
        x = state[:2]      # 真实状态
        x_hat = state[2:]  # 估计状态
        
        y = C @ x
        y_hat = C @ x_hat
        
        u = -K @ x_hat + r  # 用估计状态反馈
        
        # 真实系统
        x_dot = A @ x + B @ u.flatten()
        
        # 观测器
        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]          # 真实初始状态
    x_hat0 = [0, 0]      # 估计初始状态（未知真实状态）
    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))
    
    # 状态 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)
    
    # 状态 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)
    
    # 估计误差
    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)
    
    # 极点位置
    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()

频域分析

Bode 图

Bode 图显示系统的频率响应： - 幅值图：（ dB） vs （ log） - 相位图：（度） vs （ log）

稳定裕度

增益裕度：相位为 时，还可以增加多少增益才使系统不稳定。

相位裕度：增益为 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():
    """频域分析演示"""
    
    # 系统传递函数
    num = [1]
    den = [1, 2, 1]  # (s+1)^2
    sys = signal.TransferFunction(num, den)
    
    # Bode 图
    w, mag, phase = signal.bode(sys, n=500)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 幅值
    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)
    
    # 相位
    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)
    
    # 开环传递函数（带控制器）
    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)
    
    # 开环 Bode 图
    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()

应用实例：四旋翼飞行器

简化模型

考虑四旋翼在一个轴上的姿态控制：

$$

I = $$

其中 是俯仰角， 是力矩。

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
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():
    """四旋翼姿态控制"""
    
    I = 0.01  # 转动惯量
    
    def quadrotor(state, t, tau):
        theta, omega = state
        theta_ddot = tau / I
        return [omega, theta_ddot]
    
    # PID 控制器
    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  # 使用角速度作为微分
            
            tau = self.Kp * error + self.Ki * self.integral + self.Kd * derivative
            return tau
    
    # 仿真参数
    dt = 0.001
    t_final = 5
    t = np.arange(0, t_final, dt)
    
    # 期望轨迹：步进变化
    theta_des = np.zeros_like(t)
    theta_des[t > 1] = 0.5  # 30 度步进
    theta_des[t > 3] = 0    # 返回水平
    
    # 不同控制器
    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)  # 力矩限制
            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)
    
    # 噪声鲁棒性测试
    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)):
        # 添加测量噪声
        theta_measured = state[0] + 0.01 * np.random.randn()
        omega_measured = state[1] + 0.1 * np.random.randn()
        
        # 添加干扰力矩
        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)
    
    # 四旋翼示意图
    ax4 = axes[1, 1]
    ax4.set_xlim(-2, 2)
    ax4.set_ylim(-1.5, 1.5)
    
    # 画机身
    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)
    
    # 画螺旋桨
    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()

总结

本章我们学习了控制理论的核心内容：

  1. 传递函数：将微分方程转化为代数问题
  2. PID 控制：工业中最常用的控制器
  3. 状态空间方法：现代控制理论的基础
  4. 可控性与可观性：系统能否被控制和观测
  5. 极点配置：任意指定闭环极点
  6. LQR 最优控制：最小化代价函数的最优控制
  7. 观测器设计：从输出重构状态
  8. 频域分析：从频率响应角度理解系统

控制理论将微分方程从描述工具提升为设计工具——我们不仅能分析系统，还能主动设计系统的行为。从飞机自动驾驶到手机防抖，控制理论的应用无处不在。

练习题

基础题

  1. 求一阶系统 的阶跃响应，并画出曲线。

  2. 对于二阶系统：

    • 计算自然频率 和阻尼比 - 判断系统是欠阻尼、临界阻尼还是过阻尼
    • 计算阶跃响应的超调量

  3. 设计一个 PID 控制器使以下系统的阶跃响应满足：超调量<10%，调节时间<2 秒： $$

G(s) = 将以下微分方程转化为状态空间形式： + 6 + 11 + 6y = u判断以下系统的可控性和可观性：

A = , B = , C =

$$

进阶题

  1. 串级控制：设计一个串级 PID 控制器控制温度系统：

    • 内环：加热器功率控制
    • 外环：温度控制 比较单回路和串级控制的性能。

  2. 系统辨识：给定阶跃响应数据，估计一阶系统的参数 和。

  3. 证明对于 LQR 问题，如果 可控且 可观，则闭环系统渐近稳定。

  4. 离散化：将连续时间系统 离散化为。讨论采样周期的选择。

  5. 设计一个 Luenberger 观测器，使估计误差的收敛速度是控制器极点的 3 倍。

编程题

  1. 实现一个通用的 PID 控制器类，包括：

    • 积分饱和（ anti-windup）
    • 微分滤波
    • 输出限幅

  2. 编写程序自动进行 Ziegler-Nichols 参数整定。

  3. 实现一个完整的倒立摆仿真器，包括：

    • 非线性动力学
    • LQR 控制器
    • 状态观测器
    • 可视化动画

  4. 使用强化学习（如 PPO）训练一个神经网络控制器，与 LQR 比较性能。

  5. 实现一个数字滤波器设计工具，支持 Butterworth 、 Chebyshev 等类型。

思考题

  1. PID 控制器如此简单，为什么仍然是工业中最常用的控制器？

  2. 讨论模型不确定性对控制器性能的影响。如何设计鲁棒控制器？

  3. 比较极点配置和 LQR 的优缺点。什么情况下用哪种方法？

  4. 为什么观测器极点通常设计得比控制器极点快？

  5. 讨论机器学习方法（如深度强化学习）与传统控制理论的关系和各自适用场景。

参考资料

  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.

本文是《常微分方程的世界》系列的第 16 章。

  • 本文标题：常微分方程（十六）控制理论基础
  • 本文作者：Chen Kai
  • 创建时间：2019-06-23 10:30:00
  • 本文链接：https://www.chenk.top/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%EF%BC%88%E5%8D%81%E5%85%AD%EF%BC%89%E6%8E%A7%E5%88%B6%E7%90%86%E8%AE%BA%E5%9F%BA%E7%A1%80/
  • 版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
 评论