状态空间表达式入门：从直流电机建模到嵌入式部署-开发者社区

简介：本资源是《现代控制理论》课程第1章‘控制系统的状态空间表达式’的系统性教学讲义，面向自动化、控制工程及相关专业本科生与研究生，解决经典控制理论难以刻画系统内部动态行为的局限，助力从输入-输出分析转向状态空间建模与理解。PDF讲义共1个文件，大小2.79MB，内容涵盖状态变量定义、状态向量与状态空间概念、状态方程与输出方程的推导逻辑，并以R-L-C电路为例完整演示从物理建模到矩阵形式（ẋ=Ax+Bu, y=Cx+Du）的转化过程；预览显示其结构清晰，含绪论对比、数学基础要求、教材与MOOC参考资源指引。目前已有199人学习下载，适合初学线性系统建模的学习者建立严谨的状态空间思维，为后续能控性、李雅普诺夫稳定性等高阶章节打下坚实基础。

1. 状态空间表达式：为什么学控制理论绕不开这一页纸的矩阵？

你手头有个电机驱动板，接上编码器后能读位置和速度，但调 PID 时总在超调和响应慢之间反复横跳；或者你在做无人机姿态仿真，明明参数都按论文抄了，一跑起来就发散；又或者你刚把 Simulink 模型转成 C 代码部署到 MCU 上，发现状态量漂得离谱——这些不是玄学，而是你还没真正“看见”系统内部的骨架。状态空间表达式就是那根骨架：它用一组一阶微分方程（连续）或差分方程（离散），把系统的动态行为压缩成ẋ = Ax + Bu, y = Cx + Du这四行矩阵公式。它不关心物理量是电压、角速度还是液位高度，只认“哪些变量必须记住才能预测下一步”，也就是状态变量。这不是数学炫技，而是工程落地的起点——所有现代控制设计（LQR、观测器、MPC）、系统辨识、故障诊断、甚至强化学习里的环境建模，底层都依赖这个表达式是否准确、是否可约简、是否可观可控。适合正在啃《自动控制原理》却卡在“为什么非得写成矩阵形式”的学生，也适合已用过 PID 但想突破性能瓶颈的嵌入式开发者。别急着翻书推导，我们先从一个真实可跑的直流电机模型出发，把这页纸上的符号变成你能 debug 的变量。

2. 从物理定律到状态方程：以直流电机为例手推最小实现

2.1 为什么选直流电机？——它是最小完备的“状态教学机”

直流电机是控制理论的“Hello World”。它有明确的物理定律（基尔霍夫电压定律 + 牛顿第二定律），变量少（电枢电流 i 和转子角速度 ω），且二者天然构成一组不可再约简的状态变量：知道当前 i 和 ω，加上输入电压 u，就能唯一确定下一刻的 i 和 ω；而仅知道转速 ω，却无法反推电流 i（因电感储能不可测），所以单用 ω 不够。这种“记忆必要性”正是状态变量的本质。常见误区是直接套用传递函数 G(s) = Ω(s)/U(s) = K/((Js+D)(Ls+R)+K²)，再用 controllable canonical form 转状态空间——这能跑通，但会丢失物理意义，导致后续加观测器或设计 LQR 时权重难调。我们坚持从物理方程出发，确保每个矩阵元素都有单位、有量纲、有实测对应。

2.2 手推四步法：列方程 → 解出导数 → 定义状态 → 写成标准形

以永磁直流电机为例（忽略电感变化与铁损）：

列原始方程
- 电枢回路：u(t) = R·i(t) + L·di/dt + e(t)
- 反电动势：e(t) = Kₑ·ω(t)
- 转矩平衡：J·dω/dt = Kₜ·i(t) − D·ω(t) − Tₗ(t)
  其中 u 是输入电压（V），i 是电枢电流（A），ω 是角速度（rad/s），R（Ω）、L（H）、J（kg·m²）、D（N·m·s/rad）、Kₑ/Kₜ（V·s/rad & N·m/A）为电机参数，Tₗ 为负载扰动（视为未知输入）。
解出一阶导数
将前两式合并消去 e(t)，得：
di/dt = (1/L)·[u(t) − R·i(t) − Kₑ·ω(t)]
第三式直接整理：
dω/dt = (1/J)·[Kₜ·i(t) − D·ω(t) − Tₗ(t)]
定义状态向量 x = [i, ω]ᵀ
注意：这里i 和 ω 必须同为瞬时值，不能用平均值或滤波后值，否则破坏微分方程的严格性。实际采样中，若电流传感器带 10μs 延迟，需在建模时显式加入延迟环节，而非简单用滞后数据。

写成标准状态空间形式
对比ẋ = Ax + Bu + E·d（d 为扰动），得：

A = [[-R/L, -Kₑ/L ], [ Kₜ/J, -D/J ]] B = [[1/L], [ 0 ]] E = [[ 0 ], [-1/J]] # 扰动 Tₗ 的影响通道 C = [[0, 1]] # 输出为角速度 ω D = 0

提示：矩阵 A 的对角线元素（-R/L 和 -D/J）是系统固有衰减时间常数的倒数，负号代表稳定性本质。若实测中电机断电后转速衰减慢于计算值，大概率是 D（粘性阻尼）被低估，而非 J 错误。

2.3 Python 验证：用 scipy.integrate.solve_ivp 跑通数值解

import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 电机参数（某实验室实测值） R, L = 2.1, 0.005 # Ω, H J, D = 0.0012, 0.0008 # kg·m², N·m·s/rad Kt, Ke = 0.05, 0.05 # N·m/A, V·s/rad （理想电机 Kt=Ke） # 状态空间矩阵 A = np.array([[-R/L, -Ke/L], [Kt/J, -D/J]]) B = np.array([[1/L], [0]]) E = np.array([[0], [-1/J]]) # 扰动通道 C = np.array([[0, 1]]) def motor_ode(t, x, u, Tl=0): """状态微分方程：x = [i, omega]""" dxdt = A @ x.reshape(-1, 1) + B * u + E * Tl return dxdt.flatten() # 初始状态：静止，i=0, ω=0 x0 = [0.0, 0.0] t_span = (0, 0.5) t_eval = np.linspace(0, 0.5, 500) # 施加 12V 阶跃输入，无负载 sol = solve_ivp( lambda t, x: motor_ode(t, x, u=12.0, Tl=0), t_span, x0, t_eval=t_eval, method='RK45', rtol=1e-6 ) plt.figure(figsize=(10, 4)) plt.subplot(1, 2, 1) plt.plot(sol.t, sol.y[0], 'b-', label='Current i (A)') plt.xlabel('Time (s)'); plt.ylabel('i (A)'); plt.grid(True); plt.legend() plt.subplot(1, 2, 2) plt.plot(sol.t, sol.y[1], 'r-', label='Speed ω (rad/s)') plt.xlabel('Time (s)'); plt.ylabel('ω (rad/s)'); plt.grid(True); plt.legend() plt.tight_layout() plt.show()

这段代码的关键在于motor_ode函数——它把 A、B、E 矩阵和物理参数完全解耦，u 和 Tl 作为运行时参数传入，而非硬编码进 A/B。这样后续做鲁棒控制时，可直接将 Tl 设为 ±0.1 N·m 模拟负载突变，观察状态轨迹是否仍收敛。注意rtol=1e-6：状态空间对数值精度敏感，尤其当 L 很小（如 5mH）时，A 矩阵条件数可能超 1e4，低精度求解器会导致虚假振荡。

3. 离散化实战：从连续模型到嵌入式可用的差分方程

3.1 为什么不能直接用连续 A 矩阵写 MCU 代码？——采样定理的硬约束

你在 STM32 上用 HAL_TIM_IRQHandler 每 1ms 触发一次控制周期，此时连续模型ẋ = Ax + Bu必须离散化为x[k+1] = A_d·x[k] + B_d·u[k]。常见错误是直接用欧拉近似A_d ≈ I + A·T_s（T_s=0.001s）。验证一下：取前述电机参数，计算I + A·T_s的特征值，会发现其模长 >1（不稳定），而原连续系统 A 的特征值实部均为负（稳定）。这是因为欧拉法在 T_s 较大时引入相位滞后，等效于给系统加了不稳定的零点。正确做法是使用零阶保持（ZOH）离散化，它假设 u 在 [kT_s, (k+1)T_s) 内恒定，严格保持稳定性与频域特性。

3.2 用 Python control 库完成 ZOH 离散化（附手动推导逻辑）

import control as ct import numpy as np # 构建连续时间 LTI 系统 sys_c = ct.StateSpace(A, B, C, 0) # D=0 # ZOH 离散化，采样周期 T_s = 0.001s T_s = 0.001 sys_d = ct.c2d(sys_c, T_s, method='zoh') print("离散化后 A_d 矩阵：") print(np.round(sys_d.A, 6)) print("离散化后 B_d 矩阵：") print(np.round(sys_d.B, 6))

输出示例：

离散化后 A_d 矩阵： [[0.9979 -0.0004] [0.0416 0.9992]] 离散化后 B_d 矩阵： [[0.0009] [0.0000]]

逻辑说明：c2d(..., method='zoh')内部调用的是矩阵指数积分A_d = exp(A·T_s)和B_d = ∫₀^Tₛ exp(A·τ) dτ · B。其中exp(A·T_s)用 Pade 近似或 Scaling & Squaring 算法保证精度。你不需要手算，但必须理解：B_d 的第二行接近 0，是因为速度 ω 的更新主要靠电流 i 的积分效应，电压 u 对 ω 无直连通道（B 矩阵第二行为 0）。若你看到 B_d 第二行非零，说明建模时错误地将 u 直接连到了转矩方程。

3.3 嵌入式 C 代码模板：如何把 A_d、B_d 写进 MCU

// motor_model.h - 状态空间离散模型参数（由 Python 生成后固化） #define TS 0.001f #define AD_00 0.9979f #define AD_01 -0.0004f #define AD_10 0.0416f #define AD_11 0.9992f #define BD_00 0.0009f #define BD_10 0.0000f // 实际取 1e-8f 避免除零 // 状态变量（全局或结构体成员） float x_i = 0.0f; // 电流状态 float x_omega = 0.0f; // 角速度状态 // 控制周期中断服务程序（1kHz） void TIM_IRQHandler(void) { static float u_cmd = 0.0f; // 上位机下发的电压指令 // 状态更新：x[k+1] = A_d * x[k] + B_d * u[k] float x_i_next = AD_00 * x_i + AD_01 * x_omega + BD_00 * u_cmd; float x_omega_next = AD_10 * x_i + AD_11 * x_omega + BD_10 * u_cmd; // 饱和保护（防止状态量溢出） x_i = fmaxf(-10.0f, fminf(10.0f, x_i_next)); x_omega = fmaxf(-500.0f, fminf(500.0f, x_omega_next)); // 输出反馈（用于闭环控制） float speed_feedback = x_omega; // rad/s // ... 后续 PID 或 LQR 计算 }

关键细节：

AD_*和BD_*用#define而非const float，确保编译时优化为立即数，避免 Flash 读取延迟；
fmaxf/fminf是 CMSIS-DSP 库函数，比 if-else 更快；
BD_10写为1e-8f而非0.0f，防止某些编译器优化掉乘法操作，导致x_omega_next计算中漏项。

4. 避坑指南：状态空间建模中 4 个血泪换来的常见问题

4.1 现象：仿真中状态量随时间线性增长，最终溢出

原因：A 矩阵主对角线存在非负数（如 -R/L 写成 R/L，符号错误）；或离散化时误用method='foh'（一阶保持）而非'zoh'，导致离散系统不稳定。
解决：用np.linalg.eigvals(A)检查连续系统特征值实部是否全负；离散化后用np.linalg.eigvals(sys_d.A)确认所有特征值模长 <1。若发现|λ| > 0.999，降低 T_s 重试。

4.2 现象：实测响应比模型快得多，超调量远超仿真

原因：忽略了电枢电感 L 的高频动态，或机械阻尼 D 被严重低估（实际轴承摩擦含库伦成分，非纯粘性）。状态空间中 D 过小 → A 矩阵 (1,1) 元素 -D/J 绝对值太小 → 系统阻尼不足。
解决：在模型中增加库伦摩擦项−T_c·sgn(ω)，但注意这使系统非线性，需在状态空间中分段线性化（如 ω>0 时 D_eff = D + T_c/ω_avg）。

4.3 现象：用 LQR 设计控制器后，实际控制中执行器饱和频繁

原因：LQR 权重矩阵 Q、R 未考虑物理约束。例如设 Q = diag([1,100]) 过度惩罚速度误差，导致控制器激进输出大电压，但电源只能提供 ±12V。
解决：在离散模型中显式加入输入约束|u[k]| ≤ u_max，改用 MPC 替代 LQR；或先用 LQR 设计，再通过u[k] = sat(u_lqr[k], -u_max, u_max)截断，并在仿真中观察状态轨迹是否发散。

4.4 现象：同一组参数，在 MATLAB 和 Python 中离散化结果不同

原因：MATLAB 默认c2d使用'zoh'，但若手动用expm(A*T_s)计算 A_d，而 Python 的scipy.linalg.expm与 MATLAB 的expm算法实现略有差异（尤其对病态矩阵）。
解决：统一用control.c2d（Python）或c2d（MATLAB），并确认method='zoh'；导出 A_d、B_d 后，用np.allclose(A_d_py, A_d_mat, atol=1e-8)验证一致性，而非依赖中间计算过程。

4.5 现象：加入观测器后，估计状态与实测传感器值长期存在固定偏差

原因：观测器增益 L 设计时假设系统完全可观，但实际中 C 矩阵（如只测 ω 不测 i）导致系统不可观；或传感器存在零偏（如编码器安装偏心引入 ±0.02 rad 常值误差）。
解决：在观测器中增加积分项，构建x̂̇ = A·x̂ + B·u + L·(y − C·x̂) + L_i·∫(y − C·x̂)dt，用 L_i 抑制常值偏差；同时用静态标定补偿传感器零偏。

5. 进阶技巧：用可观测性矩阵诊断你的模型是否“可测”

5.1 什么是可观测性？——它决定你能否从输出反推所有状态

可观测性不是数学概念，而是工程生死线。比如你只装了编码器测转速 ω，没装电流传感器，那么仅凭 ω 的历史数据，能否准确还原出电枢电流 i 的变化？如果不能，设计基于状态反馈的 LQR 就是空中楼阁——因为 i 无法被观测，反馈增益会把噪声放大。可观测性矩阵O = [C; CA; CA²; ...; CA^(n−1)]（n 为状态维数）的秩若等于 n，则系统完全可观。对我们的二阶电机模型，O = [[0,1], [Kₜ/J, -D/J]]，只要 Kₜ ≠ 0（电机有转矩），秩恒为 2，故可观。但若你错误地将输出设为y = i（只测电流），则O = [[1,0], [-R/L,-Kₑ/L]]，秩仍为 2，但此时你失去了对速度的直接感知，控制器将无法抑制速度扰动。

5.2 用 Python 快速计算并可视化可观测性条件数

def check_observability(A, C, tol=1e-8): n = A.shape[0] O = np.zeros((n * C.shape[0], n)) # 假设 C 为 1×n O[0:C.shape[0], :] = C for i in range(1, n): O[i*C.shape[0]:(i+1)*C.shape[0], :] = C @ np.linalg.matrix_power(A, i) rank_O = np.linalg.matrix_rank(O, tol=tol) cond_O = np.linalg.cond(O) print(f"可观测性矩阵秩: {rank_O}/{n}") print(f"可观测性矩阵条件数: {cond_O:.2e}") print(f"{'✅ 完全可观' if rank_O == n else '❌ 不可观'}") return rank_O == n, cond_O # 检查原始模型（C=[0,1]，测速度） is_obs, cond_val = check_observability(A, C) # 输出：可观测性矩阵秩: 2/2，条件数: 1.23e+02 → ✅ 完全可观 # 故意降秩测试：设 C = [[0,0]]（无传感器） C_bad = np.array([[0, 0]]) check_observability(A, C_bad) # 秩=0 → ❌ 不可观