简介这是一份现代控制理论课程中“控制系统的状态空间表达式”章节的PDF课件内容以状态变量、状态向量与状态空间为核心从R-L-C电路等案例引出状态方程与输出方程的矩阵写法并对经典控制理论同现代控制理论的研究对象、数学工具及方法体系作了对比梳理。资料适合具备线性代数与常微分方程基础、正在学习现代控制理论或准备期末复习的本科学生使用。资源包共1个文件为PDF格式整体大小约2.79MB可直接用于阅读与打印。目前已有199人学习浏览课件页面保留课程要求、参考书目与线上学习资源链接便于延伸自学。通过这份讲义读者可以快速把握状态空间表达式的基本框架与符号约定为后续理解状态方程求解、能控能观性以及李雅普诺夫稳定性等内容奠定基础。1. 状态空间表达式为什么现代控制系统绕不开这个起点刚接触控制的时候大家都习惯用传递函数描述输入输出关系一个方块一个箭头S域里算增益算相位简单直观。可真到做多输入多输出系统、非线性模型或状态反馈控制时传递函数就开始不够用了它只告诉你输出怎么变系统内部哪些变量在动态变化完全是一团黑匣子。状态空间表达式把这个黑匣子拆开用一组一阶微分方程把系统内部的“状态”显式建出来所有影响未来输出的信息都被整理成固定的矩阵形式控制器设计、观测器设计、仿真验证全部围绕这组矩阵展开。这篇笔记就从建模、转换、仿真到落地避坑把状态空间表达式从课本定义变成你电脑上能直接跑通的工具箱。适合正在补现代控制理论的工程师、做嵌入式控制算法或者搞仿真建模的同学你可以先照着后面的操作在自己的模型上跑一遍再回头理解为什么状态变量要这么选。2. 从物理系统建立状态空间表达式状态变量选取与矩阵成形2.1 状态变量怎么选选少了会翻车选多了是浪费状态空间表达式最核心的一步不是写矩阵而是决定哪些量作为状态变量。状态变量的定义是知道这一组变量在当前时刻的值再加上输入信号从现在开始的变化就能唯一确定系统在任意未来时刻的行为。换句话说状态变量包含了系统的全部记忆它是最小的一组内部变量。工程上最常见的翻车方式是只拿输出当状态变量。比如一个质量-弹簧-阻尼系统输出是位移y但如果只把y当作状态变量就没法用一个一阶微分方程描述速度的动态变化。位移的导数vy是另一个独立的记忆单元必须把它也纳入状态向量。反过来如果你把位移和两倍位移同时设为状态这两个变量线性相关A矩阵会奇异仿真直接炸掉或者后续算能控性时发现矩阵秩缺白折腾半天。怎么判断选对了理论上用能控性和能观性检验工程上先数系统里有多少个独立储能元件。机械系统的位移对应势能、速度对应动能电路系统的电容电压和电感电流也是储能元件每个独立储能对应一个状态变量。如果你数了半天发现状态数跟储能元件数对不上大概率是漏了某个动态环节。另外要注意非线性系统在平衡点线性化后状态变量通常是相对于平衡点的增量这时候原始物理量本身可能不是状态它的偏差才是。2.2 一个质量-弹簧-阻尼系统的完整建模过程我们用一个经典二阶机械系统走一遍完整流程。运动方程为m y c y k y f(t)其中y是质量块的位移f(t)是外力m是质量c是阻尼系数k是弹簧刚度。这里输出是位移y输入是外力f系统是二阶的明显需要两个状态变量。定义状态变量x1 y位移 x2 y速度那么直接可以得到一阶微分方程组x1 x2 x2 (f - c x2 - k x1) / m整理成状态空间表达式就是我每次都会先写在草稿纸上的标准矩阵形式。对应代码里就是import numpy as np # 机械系统参数质量、刚度、阻尼 m 2.0 # 质量, kg k 8.0 # 弹簧刚度, N/m c 1.5 # 阻尼系数, N*s/m # 状态向量 x [位移, 速度] A np.array([[0, 1], [-k/m, -c/m]]) B np.array([[0], [1/m]]) C np.array([[1, 0]]) # 输出位移 D np.array([[0]]) print(A , A) print(B , B) print(特征值 , np.linalg.eigvals(A))这段代码的逻辑是A矩阵第一行负责描述x1x2第二行的-k/m和-c/m分别对应弹簧回复力和阻尼对加速度的影响B矩阵只有第二个元素非零表示外力直接作用于加速度通道也就是速度状态的变化率。C矩阵把状态向量映射到输出这里只取位移所以是[1, 0]。D矩阵为零因为外力不会直接出现在输出中。参数说明m、k、c必须统一单位工程上最常见的错误是刚度用N/mm、阻尼用N·s/m混着写最后算出来的特征值完全对不上。你可以把阻尼c设成0再看特征值会发现特征值变成纯虚数系统会持续振荡不衰减这在仿真里非常容易看出来。实际物理系统多少都有阻尼但建模时先跑一遍c0的情况反而能帮你确认矩阵写对了没有。2.3 把微分方程写成矩阵形式标准型和参数对应关系从力学方程到矩阵形式很多新手卡在“为什么A矩阵第二行是[-k/m, -c/m]而不是[-k, -c]”。这是因为状态变量本身就是位移和速度加速度x2 (f - cx2 - kx1)/m把f当成输入u剩下的系数就是-k/m和-c/m。理解这个对应关系比背标准型重要得多因为你面对真实系统时物理参数单位不同代入矩阵的位置也就不同。状态空间表达式的形式不是唯一的。同一个物理系统可以写出无数种A、B、C、D的组合但系统的特征值、传递函数是唯一的。标准型只是把矩阵整理成便于计算的形状最常见的是能控标准型和能观标准型会在第3章讨论。这里提醒一个维度匹配原则A是n×nB是n×mC是p×nD是p×m其中n是状态数m是输入数p是输出数。很多报错都源于维度没对齐B写成了1×n而不是n×1。我在实际建模时会在代码里强制加维度检查把错误拦截在运行之前n_states A.shape[0] n_inputs B.shape[1] n_outputs C.shape[0] assert B.shape (n_states, n_inputs), fB矩阵维度错误: {B.shape} assert C.shape (n_outputs, n_states), fC矩阵维度错误: {C.shape} assert D.shape (n_outputs, n_inputs), fD矩阵维度错误: {D.shape} print(维度检查通过)这段断言代码看起来不起眼但能帮你省掉大量排查矩阵乘法的低级错误。尤其是后面做离散化、状态反馈、观测器设计时任何矩阵维度错位都会在np.dot那一行报出难以理解的错误提前检查是最便宜的后悔药。另外真实工程里很少直接拿到理想的线性模型更多的是非线性微分方程。常见做法是在工作点附近做线性化先写出非线性状态方程然后对状态和输入求偏导数得到雅可比矩阵再代入工作点求出线性化的A和B。小信号模型叠加在工作点之上这就是状态空间表达式在工程中最常见的落地方式。做电机控制、电力电子、机械系统仿真时这套流程比从传递函数出发要通用得多。3. 从传递函数到状态空间三种实现方式与最小实现3.1 能控标准型与能观标准型什么时候用哪个很多控制系统设计是从传递函数开始的毕竟频域辨识工具成熟、实验数据也容易拟合。但后面要设计状态反馈或观测器就必须把传递函数转换成状态空间表达式。这里的核心问题是同一个G(s)可以有很多套状态空间实现你选哪一种直接关系到后续设计是否顺手。能控标准型也叫可控标准型把状态变量定义为输出的逐阶导数组合A矩阵呈现为友矩阵形式B矩阵的最后一个元素为1。这种形式特别适合做状态反馈极点配置因为反馈增益K和特征多项式系数是一一对应的代入公式就能算出结果。能观标准型则是能控标准型的对偶把A矩阵转置C矩阵变成最后一个元素为1适合设计观测器因为观测器增益L也可以通过对偶原理直接套用能控标准型的公式。以三阶系统为例分母是s³ a2 s² a1 s a0分子是b2 s² b1 s b0能控标准型的矩阵如下import numpy as np # 传递函数分母系数从高次到常数已经归一化 a0, a1, a2 2.0, 3.0, 1.0 b0, b1, b2 1.0, 0.5, 0.0 A np.array([[0, 1, 0], [0, 0, 1], [-a0, -a1, -a2]]) B np.array([[0], [0], [1]]) C np.array([[b0, b1, b2]]) D np.array([[0]]) print(能控标准型 A , A) print(B , B)这里A矩阵最后一行是特征多项式的负系数B矩阵末位是1C矩阵直接就是分子系数。如果分母首项系数不为1先整体除以首项系数否则你会发现A矩阵最后一行对不上特征多项式。选择原则很简单只做状态反馈用能控标准型只做观测器设计用能观标准型两个都要做建议先把系统转换到约旦标准型或者直接交给计算机库处理避免自己手算特征向量时出错。3.2 用常见科学计算库一行完成转换手动构造标准型容易但碰到高阶系统或者复杂零点分布时手算就很痛苦。我一般直接用scipy.signal.tf2ss这个函数它返回的是SISO系统的能控标准型实现之后再做后续处理。from scipy.signal import tf2ss # 传递函数 G(s) (s 2) / (s^2 3s 2) num [1, 2] den [1, 3, 2] A, B, C, D tf2ss(num, den) print(A , A) print(B , B) print(C , C) print(D , D)注意scipy.signal.tf2ss返回的矩阵顺序和很多教材上写的能控标准型可能看起来不一样。它把状态向量的顺序做了调整但特征值和输入输出响应完全一致。如果你要跟论文里的公式对应最好先打印出来检查一下A矩阵的样子再决定怎么用。另一个常用工具是control库的ss(tf_obj)但那个库返回的标准型在不同版本上可能会有差异我建议总是先用tf2ss跑一遍再用矩阵特征值验证A的特征值就是传递函数的极点如果对不上说明分子分母输入顺序写错了。3.3 最小实现为什么直接转换会得到冗余状态直接tf2ss得到的矩阵维度等于传递函数分母的阶数但如果分子分母存在公共因子也就是零极点相消这个维度里就混进了不可控或不可观的模态。比如G(s)(s1)/(s²2s1)分子分母都有s1实际系统就是一阶的但直接转换会得到一个二阶状态空间其中一个极点被抵消掉了对应的状态要么控制不到要么观察不到。这个冗余状态在仿真里看不出问题因为输入输出曲线完全正常。但一旦你在这个模型上设计状态反馈你会发现有两个极点可以配置实际上其中一个根本无法移动。更麻烦的是物理系统里那个被抵消的极点可能是不稳定的它被某种近似相消掩盖了最后控制器看着稳定实际系统一上电就发散。解决办法是做最小实现也就是把能控能观的那部分保留下来去掉所有不可控或不可观的子空间。control库有现成的minreal函数底层做的是零极点对消和状态空间化简from control import tf, ss, minreal G tf([1, 1], [1, 2, 1]) # 分子和分母有公共因子 s1 Gm minreal(G) sys_ss ss(Gm) print(最小实现极点数:, Gm.pole().size) print(A , sys_ss.A)这里输出的一阶矩阵就是真正可控可观的完整模型。但注意minreal用的零极点对消是数值方法如果公共因子距离很近但不完全相等比如s0.9999和s1.0001它也会盲目对消导致把原本存在的动态丢掉了。所以我在用minreal之前一定会先看零极点的具体位置人工确认它们确实是近似相消而不是物理上就接近但独立的两个极点。这个判断只有结合系统物理背景才能做完全交给自动化工具是会翻车的。4. 状态空间模型的仿真与求解状态转移矩阵和数值积分4.1 状态转移矩阵与解析解有了A、B、C、D下一步就是求解。线性时不变系统的状态方程xAxBu解析解形式是x(t) e^{A(t-t0)} x(t0) ∫_{t0}^{t} e^{A(t-τ)} B u(τ) dτ其中e^{At}叫状态转移矩阵它把任意时刻的状态映射到未来。这个矩阵的物理意义比公式本身重要它包含了系统在没有外部输入时的全部自由运动信息。如果A的特征值实部都为负状态转移矩阵会随着t增大趋向零系统稳定如果有正实部特征值任何小的初值扰动都会被放大。计算e^{At}的常见方法有三种拉普拉斯变换法适合手算低阶系统特征值特征向量法适合理解结构矩阵指数级数法适合数值计算。工程上最常用的是特征值法如果A可以对角化AV Λ V⁻¹那么e^{At}V e^{Λt} V⁻¹。注意如果A有重特征值或者接近重特征值V会病态这时候数值计算容易不准需要换用Pade近似算法。scipy.linalg.expm内部处理得很好但如果你手写公式就很容易踩到这个坑。4.2 用Python做数值仿真实际工程中解析解只有在线性时不变系统且输入是简单函数时才实用。更通用的做法是直接数值积分。下面这段代码用solve_ivp仿真一个带单位阶跃输入的弹簧阻尼系统并绘制位移响应曲线import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt m, k, c 2.0, 8.0, 1.5 A np.array([[0, 1], [-k/m, -c/m]]) B np.array([[0], [1/m]]) C np.array([[1, 0]]) def state_derivative(t, x): u 1.0 # 单位阶跃输入 return A x B.flatten() * u x0 [0.0, 0.0] t_span [0, 10] t_eval np.linspace(0, 10, 1000) sol solve_ivp(state_derivative, t_span, x0, t_evalt_eval, methodRK45) y C sol.y plt.plot(sol.t, y[0]) plt.xlabel(时间 (s)) plt.ylabel(位移 (m)) plt.grid(True) plt.show()逻辑说明state_derivative函数返回dx/dtsolve_ivp会自动调节步长保证误差。这里用的是RK45即四阶五阶龙格库塔法适合大多数非刚性问题。如果系统特征值差别很大比如一个极点是-1000另一个是-0.01那就是刚性系统RK45会卡在微小步长上仿真慢到怀疑人生这时候把method改成LSODA或者Radau会顺畅很多。参数说明t_eval只是输出点不改变实际积分步长。x0里的初值必须和物理量对应这里位移和速度都从0开始。如果你设置x0[0.1, 0]响应曲线初始点就不是0仔细观察可以从图上确认状态转移矩阵的贡献。另外单位阶跃输入在t0时刻从0跳到1可以不用if语句特殊处理直接在整个时间区间返回常数u即可因为solve_ivp从t0开始积分默认就是阶跃已经生效了。4.3 从仿真曲线判断系统稳定性仿真完成后不能只盯着输出曲线说“看起来稳定”。工程上判断稳定性有三道关第一算A矩阵特征值这是理论依据第二看仿真曲线在初始扰动下是否衰减到零第三换多个初值、多个输入幅值再跑几遍确认不是某个特殊初值下偶然收敛。对线性时不变系统这三道关结果应该完全一致。如果A特征值实部是正的但仿真曲线看起来收敛了几乎可以肯定是代码里A矩阵填错了或者求解器误差大到离谱。如果特征值实部为负但曲线一直小幅振荡那可能是阻尼比很小系统确实稳定但收敛极慢这在实际系统中也很有参考价值。一个常见误用是只用阶跃响应判断稳定性。阶跃输入本身有持续的能量系统可能因为输入的作用看起来输出有界但去掉输入后初值扰动依然发散。尤其对不稳定系统仿真时长不够会误判为稳定。所以我的习惯是先跑零输入响应只看初值x0[1,0]这样的衰减曲线确认A矩阵的特征值是否真的负实部再上阶跃输入分析稳态精度和超调量。零输入响应是状态空间模型独有的诊断手段这是传递函数模型很难直接做到的也是状态空间表达式在工程验证中的一个明显优势。5. 状态空间表达式落地避坑五个最常翻车的细节5.1 现象仿真一开始就输出 NaN做仿真时最常遇到的就是曲线图上第一个点就是NaN后面全空白。最初我以为是求解器坏了换了好几种方法都一样。原因A矩阵本身有奇异点或者初值x0里混入了NaN/Inf还有一种情况是积分步长过大导致数值发散尤其是系统极点远离虚轴时RK45在初始阶段可能直接越过一个很大的变化率。解决先检查初值用np.isfinite(x0)确认没有非有限值。然后打印A矩阵的奇异值如果有接近0的值说明状态变量选取有冗余或模型参数不对。最后把求解器换成LSODA并设置小一点的相对误差比如rtol1e-8。别急着怀疑算法先检查你自己的模型矩阵八成是最基础的地方出了问题。5.2 现象状态变量选错导致系统不可控翻车场景是设计状态反馈控制器时无论怎么配置极点都只能控制住一部分状态剩下的状态在响应曲线里慢慢悠悠自己动根本不听指令。原因状态变量选了线性相关的组合或者漏了关键储能状态使得能控性矩阵Co [B, AB, A²B, ...]秩亏系统存在不可控子空间。解决在设计控制器之前先算能控性矩阵的秩。代码很简单Co np.hstack([B, A B]) # 二阶系统只用到AB高阶继续拼 rank np.linalg.matrix_rank(Co) print(能控性矩阵的秩:, rank, 状态数:, A.shape[0])如果秩小于状态数说明系统不可控那就需要回去重新选择状态变量或者检查是不是把某个物理约束当成了状态。秩刚好等于状态数也不代表万事大吉我还见过数值接近满秩但条件数极大的情况极点配置出来增益巨大系统一上电就饱和。这时候要用svd看奇异值如果最小奇异值比最大奇异值小几个数量级你的系统实际上处于可控性边缘控制器不要指望它能做到任意极点配置。5.3 现象离散化后系统失稳连续域模型仿真稳定结果用零阶保持器离散化后丢到嵌入式平台跑输出直接发散或者高频振荡。原因采样周期选得太大。连续系统的极点离虚轴越近对采样周期的上限要求越苛刻。如果采样周期大于奈奎斯特对应的要求离散化后的极点可能跳出单位圆。解决选采样周期最稳妥的经验是让采样频率至少是系统最高特征频率的10到20倍。具体做法是用scipy.signal.cont2discrete做转换后立刻对比离散和连续极点from scipy.signal import cont2discrete # 连续系统矩阵 system (A, B, C, D, 0) # 最后一个0表示连续系统 dt 0.05 Ad, Bd, Cd, Dd, _ cont2discrete(system, dt, methodzoh) print(连续特征值:, np.linalg.eigvals(A)) print(离散特征值:, np.linalg.eigvals(Ad))如果离散特征值的模长接近1甚至大于1说明dt太大需要减小采样周期。另外零阶保持器和一阶保持器在高增益系统上的行为差异很大不要盲目相信默认选项。实际工程里如果控制器主导极点时间常数是0.1秒采样时间选1毫秒和选10毫秒得到的控制效果差别巨大而这个问题在连续域仿真里完全看不出来只有离散化后才会暴露。5.4 现象最小实现检查时秩判断出错调用minreal或者自己构造能控性矩阵时明明系统是不可控的np.linalg.matrix_rank却说秩等于状态数反过来也可能把可控系统判成不可控。原因matrix_rank的默认容差依赖于最大奇异值是一个相对阈值。如果你的系统矩阵数值尺度差别很大比如刚度是1e6阻尼是1e-3相对阈值会把本应视为零的小奇异值误判成非零。解决改用svd直接看奇异值分布并设置绝对阈值svs np.linalg.svd(Co, compute_uvFalse) tol 1e-10 print(奇异值:, svs) print(有效奇异值数量:, np.sum(svs tol))这里的tol要根据你系统的物理单位来定我一般先用1e-10试如果矩阵条件数实在太大再参考状态变量的数量级调整。还有一种做法是先把系统状态归一化把A矩阵除以参考量让各个状态在同一个量级下再做能控性分析。总之别拿matrix_rank的默认输出当圣旨工程上要结合物理背景解释奇异值什么是0什么只是小不是算法说了算。5.5 现象B矩阵维度不匹配导致矩阵相乘报错最常见的低级但极其耗时的错误np.dot(A, B)报错说维度不一致或者算能控性矩阵时报错。原因B矩阵在定义时写成了[0, 1/m]这样的一维数组而不是[[0], [1/m]]这样的列向量。一维数组的shape是(2,)和A的(2,2)做矩阵乘法时numpy会把一维当行向量或列向量处理结果有时能用有时报错非常迷惑。解决从建模一开始就确保B、C都是二维矩阵写一行reshape或者直接定义成嵌套列表。我更建议在代码里加一个统一的矩阵维度断言前面第2.3节已经写过了。养成打印A.shape、B.shape的习惯出现错误时第一反应不是去查矩阵数学而是看shape80%的报错都能在一分钟内解决。做离散化时cont2discrete要求输入必须是二维矩阵这个坑会在你用了很久一维数组后突然爆发别问我是怎么知道的。6. 把状态空间表达式用起来从状态反馈到观测器设计状态空间表达式的价值不在建模本身而在建模之后那一连串设计工具。到这里你已经有一个能输出状态轨迹的模型接下来最常见的一步是设计状态反馈控制器把闭环系统极点放到期望位置。设计公式很直接令u-Kx则闭环状态方程为x(A-BK)x闭环极点就是A-BK的特征值。只要系统能控就能通过选择K来任意配置闭环极点。实际设计时我一般用scipy.signal.place_poles来算增益矩阵它会自动处理多输入系统的参数自由度避免手算代数黎卡提方程的麻烦。一个完整的验证流程如下import numpy as np from scipy.signal import place_poles m, k, c 2.0, 8.0, 1.5 A np.array([[0, 1], [-k/m, -c/m]]) B np.array([[0], [1/m]]) # 期望闭环极点主导极点阻尼比约0.7自然频率约2 rad/s poles [-1.4 1.4j, -1.4 - 1.4j] result place_poles(A, B, poles) K result.gain_matrix print(状态反馈增益 K , K) print(闭环特征值 , np.linalg.eigvals(A - B K))这里K输出后需要回代验证一下而不是直接相信place_poles的返回。我遇到过因为A矩阵手滑多打了一个负号导致配置出来的极点完全不是预期的情况所以闭环特征值这行打印必须保留。另外极点位置不是越左越好越左意味着增益越大执行器会更容易饱和。工程经验是先把你关心的稳定性裕量折算成极点区域再在这个区域内选一组阻尼比和自然频率都合理的极点。状态反馈的前提是所有状态都可测量但实际系统往往只能测到输出比如只有位移传感器没有速度传感器。这时候需要设计观测器利用模型和输出重建完整状态。观测器增益L可以通过对偶原理计算把能控性问题转成对偶系统A.T和C.T的极点配置问题再把结果转置回来。一个实用的经验法则是观测器极点比控制器极点快3到10倍但也不能太快否则测量噪声会被放大到不可接受。我自己的习惯是每次拿到一个新状态空间模型第一件事就是算能控性和能观性第二件事才算极点配置。这两步加起来不过两分钟却能省下后面至少两天的调参时间。最后再提醒一句状态反馈K输出之后务必把闭环系统重新放回仿真里跑一遍阶跃响应确认超调量、调节时间和执行器饱和都在可接受范围模型里一切正常再烧进硬件。这套流程走顺了状态空间表达式就不再是课本里的一章而是你手上一件顺手的设计工具。希望这一篇能帮你在自己的系统里把它跑通。本文还有配套的精品资源点击获取