状态空间表达式全解析:从建模到Python仿真验证
简介:这是一份现代控制理论课程第1章“控制系统的状态空间表达式”的配套讲义,面向自动化、控制类专业本科生及考研复习者。内容从“状态”“状态变量”“状态空间”等核心概念切入,阐明一组最少变量在给定初始状态与输入后如何完整决定系统未来行为,进而引出状态方程与输出方程的矩阵化表示。讲义通过R-L-C电路实例,完整演示从电路微分方程到状态空间表达式的整理过程,让读者直观理解状态矩阵、输入矩阵、输出矩阵和直接传输矩阵各自的作用,掌握现代控制理论中最基础的建模方法。同时附带课程要求、主要参考书与章节框架,便于把握学习主线、安排复习节奏。讲义还说明了状态变量的物理意义和选取原则,例题步骤详细,适合课后独立推导。资源共1个PDF文件,大小约为2.79MB,排版清晰,适合打印或配合在线课程使用。已有199人学习,可帮助初学者补强线性代数在控制建模中的实际应用。
1. 状态空间表达式:一套表达式就能描述整个动态系统
某工程师调试一套双电机同步驱动系统,用传递函数矩阵建模,写了整整两页,到设计状态反馈时却发现模型根本没法继续推。换用状态空间表达式,四个矩阵 A、B、C、D 就把系统完整描述出来,仿真、观测器、最优控制全都建立在同一套表达上。这就是为什么控制理论课程把“状态空间表达式”放在第1章:它是后续所有设计方法的入口。这篇笔记针对正在跟这门课、或正在做控制算法落地的工程师,目标是让你看完就能从微分方程写出状态空间表达式,知道参数怎么选、坑在哪里,并能用 Python 快速验证自己写的模型对不对。
2. 状态变量的选取:从微分方程到矩阵微分方程的完整流程
2.1 选状态变量的三条规则:阶数、物理量、线性无关
状态变量的个数是一个硬指标:系统独立储能元件的个数,也就是微分方程的阶数。弹簧-质量系统有两个储能元件,所以需要两个状态变量;RLC 电路有电感和电容两个储能元件,同样需要两个状态变量。如果数错了,后面 A 矩阵的维度一定跟着错,仿真器大概率直接抛异常。
第二条规则是优先选物理上可测量的量,例如位移、速度、电压、电流。原因很直接:状态反馈最终要拿传感器测到的量去做控制,状态变量选成不可测的量,后面还得额外配一个观测器把状态估出来。观测器当然能解决,但那是多出来的设计负担。
第三条规则是状态变量之间必须线性无关。这里说的不是物理上“看起来不同”,而是你选的 n 个变量能张成一个 n 维空间。例如 x1 是位移、x2 是速度,这是标准选法;如果额外把 x3 也选成位移的倍数,那三个变量就只有两个独立维度,构造出的 A 矩阵会奇异,状态转移过程也存在冗余。
实际使用时三条规则会有冲突。优先保证第一条和第三条,第二条在条件允许时尽量满足。如果实在找不到物理可测的状态变量,就把输出及其各阶导数组成一组相变量,这是后面要展开的通用做法。
2.2 从弹簧-质量-阻尼系统推导状态空间表达式的完整算例
拿一个典型二阶系统来做例子。运动微分方程是:
m * y''(t) + c * y'(t) + k * y(t) = f(t)
定义状态 x1 = y(位移),x2 = y'(速度)。直接对时间求导:
x1' = x2 x2' = f/m - c/m * x2 - k/m * x1
写成矩阵形式,A、B、C、D 就齐了:
[x1'] [ 0 1 ] [x1] [ 0 ] [x2'] = [-k/m -c/m] [x2] + [1/m] [f] y = [1 0] [x1] [x2]用代码把矩阵填出来并打印极点是验证模型的第一步:
import numpy as np m = 2.0 # 质量, kg c = 0.5 # 阻尼系数, N·s/m k = 20.0 # 弹簧刚度, N/m A = np.array([[0.0, 1.0], [-k/m, -c/m]]) B = np.array([[0.0], [1.0/m]]) C = np.array([[1.0, 0.0]]) D = np.array([[0.0]]) print("A =", A, sep="\n") print("极点 =", np.linalg.eigvals(A))这段代码里的 m、c、k 就是物理参数,改它们能直接观察对极点位置的影响。A(0,1)=1 表示 x1 的导数等于 x2,这是状态定义决定的,不能随便改;A(1,0)=−k/m 来自弹簧恢复力,A(1,1)=−c/m 来自阻尼力。如果系统单位不统一,例如质量用 kg 而力用 N,那么 k/m 的数值就会差好几个数量级,这是后文要重点说的量纲问题。
2.3 状态变量选取不唯一:换一组状态变量,极点不变
常见的错误是以为状态变量选法是唯一的,所以看到别人写的 A 矩阵和自己不一样就怀疑算错了。实际上只要做非奇异线性变换,系统的外部特性完全不变。
对上面的例子,取一组新状态 z1 = x1,z2 = −2*x1 − x2,也就是做线性变换 z = T x:
T = np.array([[1.0, 0.0], [-2.0, -1.0]]) A_new = T @ A @ np.linalg.inv(T) B_new = T @ B C_new = C @ np.linalg.inv(T) print("新A极点 =", np.linalg.eigvals(A_new)) print("原A极点 =", np.linalg.eigvals(A))A_new 是相似变换后的系统矩阵,特征值不变,说明极点由系统本身决定,跟状态变量选择无关。实际工程中,把状态变量从物理量换成正则坐标或误差量,往往是控制设计的一部分。反馈矩阵 K 的具体数值也会跟着变,但闭环极点位置在设计目标下保持不变。
2.4 多输入多输出时,矩阵维度怎么填不踩空
前面的例子是单输入单输出。如果系统有两个输入 u1、u2,三个测量输出 y1、y2、y3,那 B 矩阵就是 n×2,C 矩阵是 3×n,D 矩阵是 3×2。
这个维度关系其实比矩阵元素值更容易错。最实用的办法是先把端口标在系统方框图上:输入箭头数量决定 B 的列数,输出箭头数量决定 C 的行数。然后对每个输出写一个 y = C x + D u 的等式,逐个核对。写代码时加断言:
assert B.shape == (A.shape[0], u_dim), "B矩阵维度应为 (n, p)" assert C.shape == (y_dim, A.shape[0]), "C矩阵维度应为 (q, n)" assert D.shape == (y_dim, u_dim), "D矩阵维度应为 (q, p)"u_dim 是输入通道数,y_dim 是输出通道数。不要只看矩阵元素值对不对,维度错了仿真器会直接拒绝运行,反而是最容易定位的一类错误。
3. 从传递函数到状态空间:能控标准型、能观标准型与最小实现
很多系统建模的起点是频域测量得到的传递函数,不是微分方程。这时候要用一个公式化的方法把传递函数转成状态空间表达式。最常见的是能控标准型(也叫相变量标准型)和能观标准型。
3.1 能控标准型:把传递函数分母系数直接填进友矩阵
对于 G(s) = (b1·s + b0) / (s² + a1·s + a0),能控标准型的关系是:
A = [[0, 1], [−a0, −a1]] B = [[0], [1]] C = [[b0, b1]] D = [[0]]
为什么这么填?因为取状态 x1 = y,x2 = y'之后,x2 的导数可以写成 −a0·x1 − a1·x2 + u,所以 A 的最后一行是 [−a0, −a1],而 B 只在最后一行为 1。这就是“友矩阵”的结构,很多教材也直接叫它 companion form。
适用场景:做极点配置状态反馈时,这种形式最受欢迎。因为 B 向量非常“干净”,设计反馈增益 K = [k1, k2] 后,闭环系统 A − B·K 的最后一行可以直接写出期望特征多项式的系数。你不需要复杂的坐标变换,直接把期望极点对应的多项式系数填进去就能反解出 K。
3.2 能观标准型:适合观测器设计,与能控标准型对偶
能观标准型的关系是:
A = [[0, −a0], [1, −a1]] B = [[b0], [b1]] C = [[0, 1]] D = [[0]]
它和能控标准型互为对偶:Ao = Acᵀ,Bo = Ccᵀ,Co = Bcᵀ。这个关系不是巧合,它来自线性系统理论里的对偶原理。如果你正在设计状态观测器,用能观标准型可以直接把观测器增益设计问题转换成一个对偶的状态反馈问题。本质上是同一套计算,只是换了视角。
工程上的选择习惯是:做控制器设计用能控标准型,做观测器设计用能观标准型。两种标准型只差一个转置,但直接按目标场景选型能省掉大量中间推导。
3.3 分子分母同阶时,D矩阵不能漏
初学者最容易漏掉的情况:G(s) = (b2·s² + b1·s + b0) / (s² + a1·s + a0),分子和分母同阶。此时不能直接把 b2 作为 C 的分量丢进去,要先做一个长除法:
G(s) = b2 + [(b1 − a1·b2)·s + (b0 − a0·b2)] / (s² + a1·s + a0)
然后对后面这个真分式部分构造能控标准型,直通项 b2 填到 D 矩阵里。工程里做 PID 控制器的状态空间实现时经常会遇到这种情况,因为 PID 的传递函数分子阶数至少是分母阶数减一,微分项一加就分子分母同阶了。D 矩阵在这个场景下代表“输入直接出现在输出端”,是物理上存在的直通路径,不是可以随便忽略的项。
一个小习惯:拿到传递函数后先比较分子分母阶次,如果同阶,先做除法,再构造状态空间表达式。这样 D 矩阵才不会被遗漏,C 矩阵里的系数也才是对的。
3.4 最小实现:先约分再构造,否则白多一阶
如果传递函数分子分母存在公因子,比如 G(s) = (s+1) / (s² + 3s + 2),约分之后其实是 1/(s+2)。直接按系数构造的表达式阶数会是 2,但系统的真正极点只有 −2,另一个状态既不可控也不可观,是一个“隐藏”在表达式里的冗余结构。用这种非最小实现去设计控制器,会出现算出来的反馈增益矩阵实际没有效果的情况。
判断最小实现可以计算可控性矩阵 Qc = [B, A·B, ...] 和可观性矩阵 Qo = [C; C·A; ...] 的秩:
import numpy as np from scipy.signal import tf2ss num = [1, 1] # s + 1 den = [1, 3, 2] # s^2 + 3s + 2 A, B, C, D = tf2ss(num, den) n = A.shape[0] Qc = np.hstack([B, A @ B]) Qo = np.vstack([C, C @ A]) print("可控性秩 =", np.linalg.matrix_rank(Qc)) print("可观性秩 =", np.linalg.matrix_rank(Qo)) print("系统阶数 =", n)直接按未约分传递函数构造时,跑出来阶数是 2,可观性秩往往只有 1,说明存在不可观模态;约分后 num=[1], den=[1,2] 时,阶数和两个秩都是 1。这才是最小实现。我的建模习惯是先做符号约分,再转状态空间,最后跑一遍秩检查,三步做完才敢把这个模型拿去设计控制器。
3.5 反推传递函数:验证 A、B、C、D 是否自洽
模型构造完,可以用公式 G(s) = C(sI − A)⁻¹B + D 反推传递函数,和原始传递函数对比。这一步能同时验证矩阵元素是否填错、D 矩阵是否漏掉、以及状态变量选择是否引入额外模态。
from scipy.signal import ss2tf num_back, den_back = ss2tf(A, B, C, D) print("反推分子 =", num_back[0]) print("反推分母 =", den_back)如果反推出来的分子分母和原始传递函数不一致,问题大概率出在 A 矩阵的符号上。比如把 −a0 填成了 +a0,传递函数会完全变样。ss2tf 返回的 num_back 是一个二维数组,取第一行即可。
4. 状态空间表达式落地必踩的坑:维度、量纲、初值与离散化
状态空间表达式看着只是四个矩阵,但毛刺都在细节里。这一类问题我总结下来集中在四块:维度、量纲、初值、离散化。
4.1 B矩阵列数写错,仿真器直接崩
现象:一个二输入系统,B 矩阵写成了 n×1,仿真一跑就报矩阵维度错误。
原因:没有按输入端口数来确定 B 的列数。系统是单输入还是双输入,B 矩阵列数是硬约束,写错一列都不行。
解决:写模型时先画方框图,数输入箭头。另一个办法是在代码里直接加断言,前面 2.4 节已经写过。这里再强调一遍:维度断言建议放在建模函数的第一行,每次装配模型都会自动检查。
提示:矩阵元素值错了仿真器不一定报错,但维度错了它一定会报错。把断言写进建模流程,等于让编译器帮你检查端口定义。
4.2 量纲不统一造成假发散
现象:理论上稳定的系统,仿真输出却快速发散,或者响应时间比实测快了几个数量级。
原因:A 矩阵里填的是 k/m 和 c/m 的数值,单位制不统一会导致这些比值偏几个数量级。比如质量用克、力用牛顿,算出来的 k/m 数值差了 1000 倍,极点位置就会完全变样。
解决:所有物理量统一用 SI 制,至少在构造矩阵前把所有单位换算好。还有一个排查技巧:把 A 矩阵的每个元素逐一和量纲公式比对。A(0,1) 的量纲是 1/s,A(1,0) 的量纲是 1/s²,A(1,1) 的量纲也是 1/s。如果发现 A(1,0) 的量纲不是 1/s²,那大概率是单位换算没到位。我从这里踩过的坑不算少,现在每次建模都会把 A 矩阵元素连同物理参数一起打印出来,人工扫一眼量纲。
4.3 初值默认是零,实际系统有储能就对不上
现象:仿真阶跃响应形状对,但起始段跟实测曲线差一大截。
原因:仿真工具默认状态初值为零,而实际系统开机时,电容可能已经有电压,机械系统可能已经有初速或初始位移。状态空间表达式本身是带着初始条件 x(t0) 的,但仿真函数不传初值就默认全零。
解决:用 initial_response 或直接在积分器里设置 x0:
import control as ct sys = ct.ss(A, B, C, D) T = np.linspace(0, 10, 1000) t, y = ct.initial_response(sys, T, X0=[0.05, 0.0])这就是为什么我调试状态反馈控制器时,会把阶跃响应和初值响应各跑一次。两个响应不一样,才说明模型包含了完整的动态特性。如果只测阶跃响应,初值问题很容易被漏掉。实测对比时,初值响应常常比阶跃响应更能暴露模型误差。
4.4 公因子未约分,非最小实现带来隐藏模态
现象:tf2ss 得到的系统是 2 阶,但实际物理系统只有 1 个极点;或者状态反馈调参时无论怎么改 K,某个模态始终不受控制。
原因:传递函数分子分母存在公因子,直接构造的表达式包含不可控或不可观模态。
解决:第一步用符号计算对传递函数约分;第二步做可控性和可观性秩检查。如果发现秩亏,不要硬用这个非最小实现做设计,回去先把传递函数约分干净。状态空间下也有自动消去冗余模态的算法,但工程上不要依赖这个,建模阶段把模型整干净才是根治。约分这个动作看起来很基础,但遇到带延迟环节或近似零极点对消的模型时,最容易被忽略。
4.5 离散化时采样周期太大,连续稳定离散发散
现象:连续系统极点都在左半平面,但用零阶保持器离散后,离散极点在单位圆外,仿真发散。
原因:采样周期 Ts 大于系统动态的最快时间常数。连续域离虚轴越远的极点,离散化后对 Ts 越敏感。Ts 太大时,零阶保持器的延迟效应会超过系统的稳定裕度。
解决:先看连续极点的实部,选 Ts 满足 Ts < 1/(10·|Re(λmax)|),也就是让采样频率至少是系统带宽的 10 到 20 倍。然后用 c2d 实测一下离散极点:
from scipy.signal import cont2discrete Ts = 0.05 # 采样周期 50ms Ad, Bd, _, _, _ = cont2discrete((A, B, C, D), Ts, method='zoh') print("离散极点 =", np.linalg.eigvals(Ad))如果离散极点模值大于 1,先把 Ts 调小一个数量级再试。这套经验在做嵌入式控制器参数整定时特别常见,采样周期不是靠拍脑袋定的,而是从连续极点出发反推的。
5. 用仿真把状态空间表达式“焊死”:响应、极点与李雅普诺夫校验
模型写完之后的第一件事不是直接上控制器设计,而是先验证你手里的状态空间表达式能真实复现系统的动态。我现在固定用三步:响应曲线、极点、李雅普诺夫方程。
第一步是阶跃响应。用 python-control 先把连续系统建出来,然后跑一条阶跃响应曲线,跟实验数据或原始微分方程的数值解做对比。如果状态空间表达式是从微分方程推出来的,两者应当完全重合。
import numpy as np import control as ct import matplotlib.pyplot as plt A = np.array([[0, 1], [-10, -1]]) B = np.array([[0], [1]]) C = np.array([[1, 0]]) D = np.array([[0]]) sys = ct.ss(A, B, C, D) t, y = ct.step_response(sys) plt.plot(t, y) plt.xlabel("time (s)") plt.ylabel("output") plt.savefig("step_verify.png")第二步是打印极点,看是否落在左半平面。这个检查只需要一行 np.linalg.eigvals(A),但每次建模都要做。你会发现不同状态变量组合得到的不同状态空间表达式,极点完全一致,这是模型正确的必要条件,也是验证“状态变量选取不唯一”这句话最快的手段。
第三步是我个人的习惯:用李雅普诺夫方程做稳定判据,而不只看极点。原因是极点判据在线性时不变系统里可靠,但遇到参数不确定或非线性系统线性化模型时,直接求解李雅普诺夫方程能顺带给出一个能量函数的轮廓,这对后续做鲁棒性分析很有用。
from scipy.linalg import solve_continuous_lyapunov Q = np.eye(2) P = solve_continuous_lyapunov(A.T, -Q) print("P特征值 =", np.linalg.eigvalsh(P))如果 P 正定,系统在李雅普诺夫意义下稳定。注意这里传的是 A.T,因为该函数解的是 AᵀP + PA = −Q 的形式。不同 scipy 版本对参数的解释略有差异,建议在本地确认矩阵方程后再用到正式项目里。
我自己的习惯是每套状态空间表达式都留一个仿真脚本,里面固定放这三段:阶跃响应、极点打印、李雅普诺夫检查。模型改动任何参数,顺手跑一遍,有问题当场暴露。这个习惯救过我好几次,因为状态空间表达式的矩阵元素太多,肉眼检查总有漏掉的时候,而仿真会诚实地告诉你哪里错了。希望帮到你。
本文还有配套的精品资源,点击获取