admin管理员组

文章数量:1130349

  • 状态空间模型
    • 向量和矩阵含义
    • 状态空间模型求解
      • 解析解
    • 状态空间到传递函数
    • 状态空间中的极点
    • 示例
      • 一个简单的例子
      • 从系统图(电气系统)开发状态空间模型
      • 二阶质量-弹簧系统(second-order mass-spring system)
  • 深入理解系统的极点-零点
    • 极点(Pole)的物理/数学含义
    • 零点(Zero)的物理/数学含义
      • SISO 情况下的“零点”
      • MIMO/状态空间下的“传输零点”
      • 举个例子帮助理解
  • 深入理解控制系统稳定性、可控性和可观性
    • 稳定性
      • 系统稳定性分类
        • 稳定系统
        • 边界稳定系统
        • 不稳定系统
    • 可控性与可观性
      • 可控性(Controllability)
      • 可观测性(Observability)
    • 简单示例:二阶质量-弹簧系统
    • 深入理解可控性和可观性
  • 状态空间方程 —— 极点配置
    • 从 PID 到极点配置
    • 举例说明特征值如何影响系统
    • 极点配置如何移动
    • 二阶系统示例
  • 对状态方程进行设计前、仿真前的体检
    • 控制系统设计前的状态方程体检
      • 1. 稳定性(Stability)
      • 2. 可控性(Controllability)
      • 3. 可观测性(Observability)
      • 4. 最小化/约简(Minimal Realization)
      • 5. 传递函数的适定性(Properness & Direct‐Feedthrough)
    • 仿真数值求解前的状态方程体检
      • 刚性(Stiffness)
  • MATLAB 绘制极点-零点分布
    • ✨ 为什么 可以 模型降阶?
    • 根据极点-零点分布进行降阶

在数学表达中,向量 通常用 粗体小写字母 表示,矩阵 通常用 粗体大写字母 表示。

状态空间模型

状态空间模型(State‐Space Model)将一个 线性时不变系统 (LTI, Linear Time-Invariant)的“内部状态”、“外部输入”与“观测输出”用向量和矩阵统一刻画:

x ˙ ( t ) = A   x ( t ) + B   u ( t ) , y ( t ) = C   x ( t ) + D   u ( t ) \dot{\boldsymbol{x}}(t) = \boldsymbol{A}\,\boldsymbol{x}(t) + \boldsymbol{B}\,\boldsymbol{u}(t), \\[10pt] \boldsymbol{y}(t) = \boldsymbol{C}\,\boldsymbol{x}(t) + \boldsymbol{D}\,\boldsymbol{u}(t) x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t)+Du(t)

向量和矩阵含义

  • 状态向量 x ( t ) ∈ R n \boldsymbol{x}(t)\in\mathbb{R}^n x(t)Rn:是一个 n × 1 n \times 1 n×1列向量,表示 系统的内部状态 n n n 是系统的状态维数。
    • 示例:如果是一个 2D 机械系统,可能有 x = [ x 1 x 2 ] = [ 位置 速度 ] \displaystyle \boldsymbol{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} \text{位置} \\ \text{速度} \end{bmatrix} x=[x1x2]=[位置速度]
    • 把系统此刻所有“内部变量”打包,就像把弹簧—质量系统的位置和速度一起装进一个向量里,它包含了系统当前的完整信息。
  • 输入向量 u ( t ) ∈ R m \boldsymbol{u}(t)\in\mathbb{R}^m u(t)Rm:是一个 m × 1 m \times 1 m×1列向量,表示系统的 控制输入外部输入 m m m 是输入的数量。
    • 示例:在机械系统里是推力/扭矩,在电路里是电压/电流。
  • 输出向量 y ( t ) ∈ R p \boldsymbol{y}(t)\in\mathbb{R}^p y(t)Rp:是一个 p × 1 p \times 1 p×1列向量,表示系统的 可观测输出 p p p 是输出的数量。
    • 这是能观测到、关心的量,比如只想测位置,就由 C C C 从完整状态中“筛”出位置分量。
    • 示例:在温度控制系统中, y \boldsymbol{y} y 可能就是测得的温度。

系统矩阵的含义

  • 状态矩阵 A ∈ R n × n A\in\mathbb{R}^{n\times n} ARn×n:描述 系统状态之间的相互关系,决定了系统的固有动态特性(例如,极点位置、稳定性、振荡频率等)

    • 描述系统“自身演化”规律——若无任何外部输入,状态如何因内部动力(弹簧、阻尼、惯性等)自然变化
  • 输入矩阵 B ∈ R n × m B\in\mathbb{R}^{n\times m} BRn×m:描述了 不同输入 u \boldsymbol{u} u 对不同状态 x \boldsymbol{x} x 的影响权重

    • 告诉我们 每一路外部输入 u u u 如何“分配”到各个状态通道上,产生加速或力矩。
  • 输出矩阵 C ∈ R p × n C\in\mathbb{R}^{p\times n} CRp×n:描述了 状态向量 x \boldsymbol{x} x 与输出 y \boldsymbol{y} y 的关系,将完整内部状态投影到我们关心的输出上。

    • C \boldsymbol{C} C 矩阵的每一行对应一个输出变量与所有状态变量之间的关系。
    • i i i 行即对应第 i i i 个输出如何由所有状态分量线性组合而来
  • 直接转换(或馈通)矩阵 D ∈ R p × m \boldsymbol{D}\in\mathbb{R}^{p\times m} DRp×m 描述了 输入 u ( t ) \boldsymbol{u}(t) u(t)输出 y ( t ) \boldsymbol{y}(t) y(t)直接传递(feedthrough)或 立即响应,即在同一时刻没有经过“状态”延迟就能影响输出。

    • 严格本性(Strictly Proper)系统 D = 0 \boldsymbol{D} = \boldsymbol{0} D=0,系统的输出必须“经过”状态动态才能对输入做出反应,传递函数 G ( s ) = C ( s I − A ) − 1 B \displaystyle G(s)=C(sI-A)^{-1}B G(s)=C(sIA)1B 在高频衰减为零。
    • 本性(Proper)系统 D ≠ 0 \boldsymbol{D}\neq \boldsymbol{0} D=0,传递函数 G ( s ) = C ( s I − A ) − 1 B + D \displaystyle G(s)=C(sI-A)^{-1}B + D G(s)=C(sIA)1B+D 在高频极限时收敛到 D D D
    • 校正/补偿在某些控制器结构(如前馈补偿)中,设计一个合适的 D D D 项可以预先补偿已知的输入—输出静态增益,改善稳态性能。

状态空间模型 天然支持 MIMO直接处理多输入多输出,不需像传递函数那样一一写多条关系。

x ˙ = A x + B u \dot x=A x + B u x˙=Ax+Bu 明确指出“当前状态”与“当前输入”如何决定“当前变化率”,可方便处理非零初始条件。

状态空间模型求解

  • 解析解(Analytical Solution)指的是 通过代数运算、积分、微分等数学公式推导出的显式表达式,结果通常是一个函数形式,写成 关于时间 t t t 的表达式

    • 明确的闭式解表达式,适用于纸面推导、公式分析;
    • 通常 无法用于复杂非线性或高维系统
  • 数值解(Numerical Solution)是通过 计算机程序/数值算法(如欧拉法、龙格库塔法、有限差分等),在 离散时间点上近似计算 出的结果。

    • 存在误差,精度依赖于时间步长和算法阶数;
    • 可以应用于几乎任意微分方程(尤其是非线性或高维系统),实际工程中使用最多。
解法类型状态解 x ( t ) \boldsymbol{x}(t) x(t) u ( t ) = u 0 \boldsymbol{u}(t) = \boldsymbol{u}_0 u(t)=u0 为常值向量时输出解 y ( t ) \boldsymbol{y}(t) y(t)
解析解 x ( t ) = e A t x 0 + ∫ 0 t e A ( t − τ ) B u ( τ ) d τ \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \int_0^t e^{\boldsymbol{A}(t - \tau)} \boldsymbol{B} \boldsymbol{u}(\tau) d\tau x(t)=eAtx0+0teA(tτ)Bu(τ)dτ x ( t ) = e A t x 0 + A − 1 ( e A t − I ) B u 0 \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \boldsymbol{A}^{-1}(e^{\boldsymbol{A}t} - \boldsymbol{I})\boldsymbol{B}\boldsymbol{u}_0 x(t)=eAtx0+A1(eAtI)Bu0 y ( t ) = C x ( t ) \boldsymbol{y}(t) = \boldsymbol{C} \boldsymbol{x}(t) y(t)=Cx(t)
数值解 x ( t i ) \boldsymbol{x}(t_i) x(ti) 为一组离散时间点的近似值(如通过 Runge-Kutta 计算) x ( t i + 1 ) ≈ x ( t i ) + h ⋅ ( A x ( t i ) + B u 0 ) \boldsymbol{x}(t_{i+1}) \approx \boldsymbol{x}(t_i) + h \cdot \left( \boldsymbol{A}\boldsymbol{x}(t_i) + \boldsymbol{B} \boldsymbol{u}_0 \right) x(ti+1)x(ti)+h(Ax(ti)+Bu0) y ( t i ) = C x ( t i ) \boldsymbol{y}(t_i) = \boldsymbol{C} \boldsymbol{x}(t_i) y(ti)=Cx(ti)

解析解

给定方程: x ˙ = A x + B u \dot{\boldsymbol{x}} = \boldsymbol{A} \boldsymbol{x} + \boldsymbol{B} \boldsymbol{u} x˙=Ax+Bu,假设初始条件为 x ( 0 ) = x 0 \boldsymbol{x}(0) = \boldsymbol{x}_0 x(0)=x0解析解公式
x ( t ) = e A t x 0 + ∫ 0 t e A ( t − τ ) B u ( τ )   d τ \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \int_0^t e^{\boldsymbol{A}(t - \tau)} \boldsymbol{B} \boldsymbol{u}(\tau) \, d\tau x(t)=eAtx0+0teA(tτ)Bu(τ)dτ
这是利用矩阵指数法推导出来的解析解,前项是自由响应,后项是受控响应(强迫响应)

  • 自然演化 e A t x 0 e^{\boldsymbol{A}t} \boldsymbol{x}_0 eAtx0:若从 t = 0 t=0 t=0 起无任何输入,系统怎样“凭借自身动力”演化到 t t t
  • 卷积累积 ∫ 0 t e A ( t − τ ) B u ( τ )   d τ \int_0^t e^{\boldsymbol{A}(t - \tau)} \boldsymbol{B} \boldsymbol{u}(\tau) \, d\tau 0teA(tτ)Bu(τ)dτ:把每一刻的输入 u ( τ ) \boldsymbol{u}(\tau) u(τ) 按“系统记忆” e A ( t − τ ) e^{\boldsymbol{A}(t - \tau)} eA(tτ) 加权累积,得出总效果。

其中,矩阵指数定义为
e A t    =    ∑ k = 0 ∞ ( A   t ) k k ! . e^{\boldsymbol{A}t}\;=\;\sum_{k=0}^\infty \frac{(\boldsymbol{A}\,t)^k}{k!}. eAt=k=0k!(At)k.

如果 u ( t ) \boldsymbol{u}(t) u(t) 是常值输入(即 u ( t ) = u \boldsymbol{u}(t) = \boldsymbol{u} u(t)=u 是一个常向量),那么数值解为
x ( t ) = e A t x 0 + ( ∫ 0 t e A ( t − τ ) d τ ) B u = e A t x 0 + A − 1 ( e A t − I ) B u \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \left( \int_0^t e^{\boldsymbol{A}(t - \tau)} d\tau \right) \boldsymbol{B} \boldsymbol{u} = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \boldsymbol{A}^{-1} \left( e^{\boldsymbol{A}t} - \boldsymbol{I} \right) \boldsymbol{B} \boldsymbol{u} x(t)=eAtx0+(0teA(tτ)dτ)Bu=eAtx0+A1(eAtI)Bu

注意事项:如果 A A A 不可逆(奇异矩阵),不能直接用 A − 1 A^{-1} A1。这时要用其他方式近似积分,比如数值方法(Runge-Kutta)

推导过程

  • 齐次解(零输入响应):当 u = 0 u=0 u=0 时,系统变为齐次
    x ˙ h ( t ) = A x h ( t ) \dot{x}_h(t) = A x_h(t) x˙h(t)=Axh(t) 其解为: x h ( t ) = e A t x ( 0 ) x_h(t) = e^{A t} x(0) xh(t)=eAtx(0)
  • 特解(零状态响应):由于 u u u 为常量,特解为,
    x p ( t ) = ∫ 0 t e A ( t − τ ) B u   d τ x_p(t) = \int_0^t e^{A (t - \tau)} B u \, d\tau xp(t)=0teA(tτ)Budτ 因为 A A A 是常矩阵,且 u u u 为常量,可得: x p ( t ) = ( ∫ 0 t e A ( t − τ )   d τ ) B u = ( ∫ 0 t e A s   d s ) B u x_p(t) = \left( \int_0^t e^{A (t - \tau)} \, d\tau \right) B u = \left( \int_0^t e^{A s} \, ds \right) B u xp(t)=(0teA(tτ)dτ)Bu=(0teAsds)Bu
    s = t − τ s = t - \tau s=tτ,则: ∫ 0 t e A s   d s = A − 1 ( e A t − I ) \int_0^t e^{A s} \, ds = A^{-1} (e^{A t} - I) 0teAsds=A1(eAtI) 因此,特解为: x p ( t ) = A − 1 ( e A t − I ) B u x_p(t) = A^{-1} (e^{A t} - I) B u xp(t)=A1(eAtI)Bu
  • 将齐次解与特解相加,得到总解
    x ( t ) = x h ( t ) + x p ( t ) = e A t x ( 0 ) + A − 1 ( e A t − I ) B u x(t) = x_h(t) + x_p(t) = e^{A t} x(0) + A^{-1} (e^{A t} - I) B u x(t)=xh(t)+xp(t)=eAtx(0)+A1(eAtI)Bu

状态空间到传递函数

我们从线性时不变系统(LTI)的状态空间模型开始:
x ˙ ( t ) = A x ( t ) + B u ( t ) , y ( t ) = C x ( t ) \dot{\boldsymbol{x}}(t) = \boldsymbol{A}\boldsymbol{x}(t) + \boldsymbol{B}\boldsymbol{u}(t), \quad \boldsymbol{y}(t) = \boldsymbol{C}\boldsymbol{x}(t) x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t)

零初始条件 x ( 0 ) = 0 \boldsymbol{x}(0) = \boldsymbol{0} x(0)=0,对上式两边进行 拉普拉斯变换 L { ⋅ } \mathcal{L}\{\cdot\} L{},就得到了:
s X ( s ) = A X ( s ) + B U ( s ) ⇒ ( s I − A ) X ( s ) = B U ( s ) s\boldsymbol{X}(s) = \boldsymbol{A} \boldsymbol{X}(s) + \boldsymbol{B} \boldsymbol{U}(s) \quad \Rightarrow \quad (s\boldsymbol{I} - \boldsymbol{A}) \boldsymbol{X}(s) = \boldsymbol{B} \boldsymbol{U}(s) sX(s)=AX(s)+BU(s)(sIA)X(s)=BU(s)

进一步得出:
X ( s ) = ( s I − A ) − 1 B U ( s ) \boldsymbol{X}(s) = (s\boldsymbol{I} - \boldsymbol{A})^{-1} \boldsymbol{B} \boldsymbol{U}(s) X(s)=(sIA)1BU(s)

再代入输出方程:
Y ( s ) = C X ( s ) = C ( s I − A ) − 1 B U ( s ) \boldsymbol{Y}(s) = \boldsymbol{C} \boldsymbol{X}(s) = \boldsymbol{C} (s\boldsymbol{I} - \boldsymbol{A})^{-1} \boldsymbol{B} \boldsymbol{U}(s) Y(s)=CX(s)=C(sIA)1BU(s)

我们将系统的 输入 U ( s ) \boldsymbol{U}(s) U(s)输出 Y ( s ) \boldsymbol{Y}(s) Y(s) 的关系写成如下形式:
G ( s ) = Y ( s ) U ( s ) = C ( s I − A ) − 1 B \boldsymbol{G}(s) = \frac{\boldsymbol{Y}(s)}{\boldsymbol{U}(s)} = \boldsymbol{C}(s\boldsymbol{I} - \boldsymbol{A})^{-1} \boldsymbol{B} G(s)=U(s)Y(s)=C(sIA)1B

这个 G ( s ) \boldsymbol{G}(s) G(s) 就是系统的 传递函数矩阵(Transfer Function Matrix),它描述了每个输入如何影响每个输出 —— 是一个 从频域角度刻画系统行为的工具

传递函数中的项 ( s I − A ) − 1 (s\boldsymbol{I} - \boldsymbol{A})^{-1} (sIA)1 是一个 矩阵的逆,它的存在依赖于:
det ⁡ ( s I − A ) ≠ 0 \det(s\boldsymbol{I} - \boldsymbol{A}) \neq 0 det(sIA)=0

因此,若 det ⁡ ( s I − A ) = 0 \det(s\boldsymbol{I} - \boldsymbol{A}) = 0 det(sIA)=0 就表示系统传递函数 G ( s ) \boldsymbol{G}(s) G(s) 的某些分量出现了 极点(即系统响应趋于无穷大)。这时对应的 s s s 值是 系统矩阵 A \boldsymbol{A} A 的特征值

状态空间中的极点

在状态空间模型中:

x ˙ = A x + B u , y = C x + D u \dot{\boldsymbol{x}} = \boldsymbol{A} \boldsymbol{x} + \boldsymbol{B} \boldsymbol{u},\quad \boldsymbol{y} = \boldsymbol{C} \boldsymbol{x} + \boldsymbol{D} \boldsymbol{u} x˙=Ax+Bu,y=Cx+Du

系统的极点(Poles),就是矩阵 A \boldsymbol{A} A特征值(Eigenvalues)。也就是满足:

det ⁡ ( s I − A ) = 0 \det(s\boldsymbol{I} - \boldsymbol{A}) = 0 det(sIA)=0

的那些 s s s 值。

极点意味着什么? 把极点看成一个物理系统“天生的频率特性”。比如一个弹簧-质量系统的震动频率,就由它的质量和刚度决定,对应的就是极点。

  1. 系统的固有动态行为
    极点(Poles)系统的本征频率或固有动态特性,反映了系统在 无输入时(自由响应) 如何自然地演化
    比如你把弹簧系统拨一下,然后放手,它怎么震动、震荡多久、会不会越来越剧烈——全由极点决定。
  2. 稳定性
    • 如果 所有极点实部 < 0:系统 稳定(响应会自然衰减)
    • 如果 存在 极点 实部 > 0:系统不稳定(响应会发散)
    • 如果极点 实部 = 0 且为简单极点边界稳定(比如无阻尼震荡)。
    • 极点有 虚部系统振荡
  3. 响应速度
    • 极点的 实部越小(越负),响应越慢
    • 实部越大(更靠近右侧),系统响应越快或越激烈

一个二阶系统(弹簧-质量):
x ¨ + 2 ζ ω n x ˙ + ω n 2 x = u ( t ) \ddot{x} + 2\zeta\omega_n \dot{x} + \omega_n^2 x = u(t) x¨+2ζωnx˙+ωn2x=u(t)
转换为状态空间后,它的极点为:
s = − ζ ω n ± j ω n 1 − ζ 2 s = -\zeta \omega_n \pm j \omega_n \sqrt{1 - \zeta^2} s=ζωn±jωn1ζ2
在 s 平面上(复平面),这个复数对极点的实部和虚部决定:

  • 衰减速度(实部)
  • 震荡频率(虚部)

示例

一个简单的例子

参考:State Space Representations of Linear Physical Systems

考虑一个由单个四阶微分方程表示的四阶系统,该方程具有输入 x x x 和输出 z z z


我们可以定义 4 个新变量, q 1 q_1 q1 q 4 q_4 q4


但是,


现在,我们可以将 4 阶微分方程重写为 4 个一阶方程,


这可以简洁地写成状态空间格式,


其中,

从系统图(电气系统)开发状态空间模型

要为电气系统开发状态空间系统,尝试选择 电容器的电压和电感器的电流作为状态变量。回忆一下,

如果我们可以写出电感两端的电压方程,那么当我们除以电感时,它就变成了状态方程。

即,如果我们有一个关于 e inductor e_\text{inductor} einductor 的方程,除以 L L L,它就变成了关于 d i inductor d t \frac{\text{d}i_\text{inductor}}{\text dt} dtdiinductor 的方程,这是我们状态变量之一。

同样,如果我们能写出通过电容的电流方程,并除以电容,它就变成了关于 e capacitor e_\text{capacitor} ecapacitor 的状态方程。

这最好通过一个例子来说明。推导出下图所示系统的状态空间模型。输入是 i a i_a ia ,输出是 e 2 e_2 e2

该系统有三个储能元件,因此我们预计有三个状态方程。尝试选择 i 1 i_1 i1 i 2 i_2 i2 e 1 e_1 e1 作为状态变量。现在我们 想要它们的导数方程

  • 电感 L 2 L_2 L2 两端的电压是 e 1 e_1 e1 (这是我们状态变量之一)。

    因此,我们第一个状态变量方程是,
  • 如果我们将电流求和到标记为 n 1 n_1 n1 的节点,我们得到,

    该方程包含我们的输入( i a i_a ia)、两个状态变量( i L 2 i_{L2} iL2 i L 1 i_{L1} iL1)以及通过电容器的电流。因此,我们可以得到第二个状态方程,
  • 我们的第三个,也是最后一个状态方程,是通过将 L 1 L_1 L1 上的电压(即 e 2 e_2 e2 )用我们的其他状态变量表示得到的,
  • 我们还需要一个输出方程:

因此,我们的状态空间表示变为

这种方法并不总是容易得到一组状态方程。在某些情况下,开发传递函数模型并将其转换为状态空间模型更容易。

二阶质量-弹簧系统(second-order mass-spring system)

一个经典的质量-弹簧(可带阻尼)系统的运动微分方程为:

m   q ¨ ( t ) + d   q ˙ ( t ) + k   q ( t ) = u ( t ) m\,\ddot{q}(t) + d\,\dot{q}(t) + k\,q(t) = u(t) mq¨(t)+dq˙(t)+kq(t)=u(t)

  • q ( t ) q(t) q(t):质量块的位置(系统的输出)
  • m m m:质量块的质量
  • d d d:阻尼系数(无阻尼时 d = 0 d = 0 d=0
  • k k k:弹簧刚度
  • u ( t ) u(t) u(t):系统所受的外力(输入)

我们定义状态变量为:

x = [ x 1 x 2 ] = [ q q ˙ ] \boldsymbol{x} = \begin{bmatrix} x_1 \\ x_2\\ \end{bmatrix} = \begin{bmatrix} q \\ \dot{q}\\ \end{bmatrix} x=[x1x2]=[qq˙]

根据微分方程可得:

x ˙ 1 = x 2 , x ˙ 2 = − k m x 1 − d m x 2 + 1 m u \dot{x}_1 = x_2,\quad\\[10pt] \dot{x}_2 = -\frac{k}{m}x_1 - \frac{d}{m}x_2 + \frac{1}{m}u x˙1=x2,x˙2=mkx1mdx2+m1u

用矩阵形式写出状态空间模型为:

x ˙ = [ 0 1 − k m − d m ] ⏟ A   x + [ 0 1 m ] ⏟ B   u y = [ 1 0 ] ⏟ C   x + [ 0 ] ⏟ D   u \dot{\boldsymbol{x}} = \underbrace{\begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & -\frac{d}{m}\\ \end{bmatrix}}_{\boldsymbol{A}}\,\boldsymbol{x} + \underbrace{\begin{bmatrix} 0 \\ \frac{1}{m}\\ \end{bmatrix}}_{\boldsymbol{B}}\,u \\[10pt] y = \underbrace{\begin{bmatrix} 1 & 0\\ \end{bmatrix}}_{\boldsymbol{C}}\,\boldsymbol{x} + \underbrace{\begin{bmatrix} 0 \end{bmatrix}}_{\boldsymbol{D}}\,u x˙=A [0mk1md]x+B [0m1]uy=C [10]x+D [0]u

  • A \boldsymbol{A} A:描述状态之间的动态关系
  • B \boldsymbol{B} B:描述输入对状态变化的影响
  • C \boldsymbol{C} C:提取输出(我们这里只观测位移)
  • D = 0 \boldsymbol{D} = 0 D=0:表示没有输入直接传递到输出(无“直通项”)

无阻尼时( d = 0 d = 0 d=0,矩阵 A \boldsymbol{A} A 简化为:

A = [ 0 1 − k m 0 ] \boldsymbol{A} = \begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & 0\\ \end{bmatrix} A=[0mk10]

深入理解系统的极点-零点

  • 极点 告诉你系统自己会怎么“颤动/衰减”——对应齐次响应,是所有模态的汇总;
  • 零点 告诉你在哪些频率上“输入到输出的链路被一脚踹断”,系统对应哪些激励会被完全消除;
  • 极点–零点相同 说明 “内部有这个模态,但外部完全看不到”,正是降阶的目标

极点(Pole)的物理/数学含义

  1. 自然模态(Natural Modes)
    • 极点 s i s_i si 就是矩阵 A A A 的特征值每一个特征值对应一个“模态”——系统在无外力(  ⁣ u = 0 \!u=0 u=0)下,会 沿着对应的特征向量方向,以 e s i t e^{s_i t} esit 的形式 自然衰减或振荡
    • 如果 s i s_i si 是实数,模式是单指数衰减 e s i t e^{s_i t} esit;如果是复数对 α ± j β \alpha\pm j\beta α±jβ,是衰减振荡 e α t cos ⁡ ( β t ) e^{\alpha t}\cos(\beta t) eαtcos(βt)
  2. 时域含义
    • ℜ ( s i ) < 0 \Re(s_i)<0 (si)<0 → 模式衰减,且 ∣ ℜ ( s i ) ∣ |\Re(s_i)| ∣ℜ(si) 越大,衰得越快,系统稳定;
    • ℜ ( s i ) > 0 \Re(s_i)>0 (si)>0模式发散,系统不稳定
    • ℜ ( s i ) = 0 \Re(s_i)=0 (si)=0 ℑ ( s i ) ≠ 0 \Im(s_i)\neq0 (si)=0 → 永久振荡(边界稳定)。
  3. 频域含义
    • 在频率响应里,离开输入的频率越靠近某个复数极点的虚部,系统对该频率的放大倍数越大(共振/峰值)。

零点(Zero)的物理/数学含义

  1. 传输零点(Transmission Zeros)
    • 零点 z j z_j zj 是这样一类频率:即使在输入端 U ( s ) U(s) U(s) 做脉冲/正弦激励,输出 Y ( s ) Y(s) Y(s) 却在此 s = z j s=z_j s=zj 频率处被完全“抑制”(幅度为零)
    • 数学上,零点是使得传递函数矩阵 G ( s ) G(s) G(s) 在某个方向上降秩(或者 分子多项式为零)的 s s s 值。
  2. 输入–输出影响
    • 有零点的模式往往是那种“输入作用到这个模态上,会被系统某种耦合机制立刻在输出端给抵消掉”,即 在传递路径上打了一个“死结”

SISO 情况下的“零点”

对单输入单输出(SISO)系统,传递函数是一个分式

G ( s ) = N ( s ) D ( s ) . G(s) = \frac{N(s)}{D(s)}. G(s)=D(s)N(s).

  • 极点:使 D ( s ) = 0 D(s)=0 D(s)=0 的那些 s s s 值,它们决定 自然模态
  • 零点:使 N ( s ) = 0 N(s)=0 N(s)=0 的那些 s s s 值,它们决定 在对应频率上,输出恰好被“抑制”到零

物理直观:当你给系统一个频率为 s = j ω s=j\omega s= 的正弦输入时,输出幅度 ∣ G ( j ω ) ∣ \bigl|G(j\omega)\bigr| G() 在零点频率 ω 0 \omega_0 ω0 处正好为零——就像你对着一个带有“峰谷滤波器”说话,那些“谷”频率会被完全消除,什么声音都传不下去。

从传递函数出发(SISO) 推导,
Y ( s ) = G ( s )   U ( s ) , U ( s ) ≠ 0. Y(s) = G(s)\,U(s),\quad U(s)\neq 0. Y(s)=G(s)U(s),U(s)=0.
G ( s 0 ) = 0 G(s_0)=0 G(s0)=0,则对任意非零 U ( s ) U(s) U(s),都有 Y ( s 0 ) = 0 Y(s_0)=0 Y(s0)=0。这就是“在 s 0 s_0 s0 频率上零增益”。

MIMO/状态空间下的“传输零点”

对于多输入多输出或者状态空间模型,不能简单地看“分子多项式”。我们用下面的 Rosenbrock 阶跃矩阵 定义零点:

[ s I − A − B C D ] ⏟ P ( s ) \underbrace{\begin{bmatrix} sI - A & -B\\[6pt] C & D \end{bmatrix}}_{\displaystyle P(s)} P(s) [sIACBD]

传输零点 s = z s=z s=z 的定义:存在非零向量 ( x 0 u 0 ) ≠ 0 \begin{pmatrix}x_0\\u_0\end{pmatrix}\neq0 (x0u0)=0,使得
P ( z )   ( x 0 u 0 ) = [ z I − A − B C D ] ( x 0 u 0 ) = 0. P(z)\,\begin{pmatrix} x_0 \\[3pt] u_0 \end{pmatrix} = \begin{bmatrix} zI - A & -B\\[3pt] C & D \end{bmatrix} \begin{pmatrix}x_0\\u_0\end{pmatrix} = 0. P(z)(x0u0)=[zIACBD](x0u0)=0.
这意味着:
{ ( z I − A )   x 0    =    B   u 0 , C   x 0 + D   u 0    =    0. \begin{cases} (zI - A)\,x_0 \;=\; B\,u_0,\\[3pt] C\,x_0 + D\,u_0 \;=\; 0. \end{cases} {(zIA)x0=Bu0,Cx0+Du0=0.

  • 第一行    ( z I − A ) x 0 = B   u 0 \;(zI - A)x_0 = B\,u_0 (zIA)x0=Bu0:表明当输入用复杂指数 u ( t ) = u 0 e z t u(t)=u_0 e^{z t} u(t)=u0ezt 激励时,可在内部激发出一个状态模态 x ( t ) = x 0 e z t x(t)=x_0 e^{z t} x(t)=x0ezt
  • 第二行    C x 0 + D u 0 = 0 \;C x_0 + D u_0 = 0 Cx0+Du0=0:表明对于这对 ( x 0 , u 0 ) (x_0,u_0) (x0,u0),输出 y ( t ) = C x ( t ) + D u ( t ) y(t)=C x(t)+D u(t) y(t)=Cx(t)+Du(t) 恰好恒为 0,输入虽作用于状态,但从状态流到输出的通路被“死结”拦住

因此, P ( s ) P(s) P(s) 的行列式
det ⁡ P ( s ) = 0 \det P(s) = 0 detP(s)=0
所求解的 s s s 就是一组传输零点。

从状态空间刻画(MIMO) 推导,令初始条件为 0,令输入 u ( t ) = u 0 e s t u(t)=u_0 e^{s t} u(t)=u0est。假设存在解
x ( t )    =    x 0 e s t , y ( t ) = C x ( t ) + D u ( t ) . x(t) \;=\; x_0 e^{s t},\quad y(t) = Cx(t)+Du(t). x(t)=x0est,y(t)=Cx(t)+Du(t).
带入微分方程:
s   x 0 e s t = A   x 0 e s t + B   u 0 e s t    ⟹    ( s I − A ) x 0 = B   u 0 . s\,x_0 e^{st} = A\,x_0 e^{st} + B\,u_0 e^{st} \;\Longrightarrow\; (sI - A)x_0 = B\,u_0. sx0est=Ax0est+Bu0est(sIA)x0=Bu0.
同时,
y ( t ) = ( C x 0 + D u 0 ) e s t . y(t)=\bigl(Cx_0 + D u_0\bigr)e^{st}. y(t)=(Cx0+Du0)est.
若想让 y ( t ) ≡ 0 y(t)\equiv0 y(t)0,必须有
C x 0 + D u 0 = 0. C x_0 + D u_0 = 0. Cx0+Du0=0.
非零解 ( x 0 , u 0 ) (x_0,u_0) (x0,u0) 存在的条件,正是 det ⁡ P ( s ) = 0 \det P(s)=0 detP(s)=0,即传输零点

举个例子帮助理解

下面举个具体的例子帮助理解,参考 How Transfer Function Zeros Affect Transient Response


这些脉冲响应来自拥有 相同极点 的系统,那么它们究竟如何产生如此巨大的差异呢?

首先注意到,它们的 传递函数分子(零点位置)是不同的。我们把这两组传递函数 泛化为二阶系统

左侧系统在 s = − z s=-z s=z 处有一个零点;右侧系统没有任何零点。

问题就转化为:一个零点的存在,为什么会导致如此不同的响应?

首先看 无零点系统 H 1 ( s ) H_1(s) H1(s) 的脉冲响应

其中极点为复共轭对(阻尼比 ζ < 1 \zeta<1 ζ<1),可分解并做部分分式展开,取逆拉普拉斯变换后得到时域脉冲响应并简化,


注意, y 1 ( t ) y_1(t) y1(t) 是指数和正弦波的乘积,指数项(极点的实部) 控制衰减或增长速度,正弦项(极点的虚部) 控制振荡频率

现在已经确定了无零脉冲响应,现在来确定 − z -z z 处有一个零点系统 H 2 ( s ) H_2(s) H2(s) 的脉冲响应

同样对分母分解、做部分分式展开,并取逆拉普拉斯变换,化简后可得,

事情开始变得有趣了,

  • 看一下第二项 z − a b \frac{z-a}{b} bza随着 z z z 接近 a a a,频率响应 y 1 y_1 y1 会缩放,特别是当 z z z 接近于 p p p 的实部 a a a 时,极点位于 − p -p p ( z − a ) (z-a) (za) 项会消失,该 模态与 y 1 y_1 y1 完全抵消
  • 第二个观察结果是,零点不会影响衰减或放大的速率(即不会改变极点的实部)

根据这些观察结果,现在要问的是,零点位置如何影响脉冲响应

叠加各个零点位置的响应

将系统极点固定在复平面上的某对位置,依次 把零点从右半平面向左半平面移动,可以同时绘制其脉冲响应,发现:

  • 零点最初位于右半平面,这会导致瞬态响应下降;
  • 当零点向左移动到左平面时,瞬态响应上升。


对两种情况的脉冲响应进行扩展,现在将单零点响应扩展为三个项,看看单零点响应中有什么?

原来的无零点项按 z z z 缩放。再看当零点从右半平面移动到左半平面时,这个项如何影响瞬态响应,

零点的位置 s s s 与第三项中包含的 z z z 值相反,

  • 零点位于左半平面时 z > 0 z > 0 z>0,因此添加了第三项,这就是为什么 z z z 变大时,瞬态响应会增大(过冲)


因此在单零瞬态响应中,我们有按 z z z 缩放的原始无零响应,以及 由于零点而产生的额外瞬态响应

这个额外的瞬态响应是什么?

对每个项进行拉普拉斯变换,我们发现瞬态是无零响应的导数,

注意传递函数 H 2 H_2 H2 现在表示为 H 1 H_1 H1 的函数,零点的影响是缩放无零响应,并将无零响应的导数添加到它,

现在考虑 H 2 H_2 H2 行为如何随着 z z z 的变化而变化,

  • 随着 z z z 的幅度变大,缩放的无零响应占主导地位;
  • 随着 z z z 幅度变小,无零响应的导数占主导地位。

在这两个极端之间, H 2 H_2 H2 H 1 H_1 H1 及其导数的组合。因此,在 y 2 y_2 y2 的时间域中,我们查看的是 y 1 y_1 y1 及其导数的线性组合。

H 2 H_2 H2 中分解出 H 1 H_1 H1,并注意 H 2 H_2 H2 是两个传递函数的乘积。

  • H 1 H_1 H1 将输入 U U U 转换为其输出 Y 1 Y_1 Y1
  • s = − z s=-z s=z 处的零点将 H 1 H_1 H1 的输出通过传递函数 s + z s+z s+z 进行变换,乘积是 H 2 H_2 H2,即单零传递函数,或者说是 H 2 H_2 H2 将输入信号 U U U 转换为输出 Y 2 Y_2 Y2

现在,我们便从 频域(零点的代数作用)时域(叠加导数项),再回到频域看其物理意义。这使我们能够从多个角度理解零点对瞬态响应的影响。

在该例子中, H 1 H_1 H1 H 2 H_2 H2 都是二阶系统,也就是说 每个系统都有两个极点(成对的复共轭极点):

  • 极点:一对位于 s = − a ± b j s=-a\pm bj s=a±bj(阻尼比 ζ < 1 \zeta<1 ζ<1)的复共轭极点。
  • 零点
    • 对于 H 1 ( s ) H_1(s) H1(s)(无零点)来说,分子是常数,只有这两个极点。
    • 对于 H 2 ( s ) H_2(s) H2(s)(带零点)来说,在分子又多了一个 ( s + z ) (s+z) (s+z) 项,因此它有一个额外的零点(位于 s = − z s=-z s=z),但极点依旧只有那对 − a ± b j -a\pm bj a±bj

当把这个 零点移动到恰好和极点位置重合 z = a z=a z=a无虚部 b b b 时),那对极点就会 在传递函数里和零点 完全相消——那一对极点在输入—输出上被抹去,系统表现上像降到了一阶或零阶,从输入到输出的路径上,这个模态“隐形”了。

在滤波器设计中,常用 严格本性(无零点)函数来抑制噪声;而一个零点的加入,实际上是对无零点响应的“反向滤波”,它通过缩放和叠加导数可能在信号上引入更大的瞬态振荡。

总结

  • 极点 决定系统的 衰减/振荡 速率(自然模态);
  • 零点 不会改变衰减率,而是表现为振荡,会通过 叠加导数项 引入 超调/下冲,实际上正在转换无零响应;
  • 当零点与极点相撞时,会 相互抵消,该模态在输入–输出通路上完全“隐形”
  • 正是零点位置的变化,让两个拥有 相同极点 的系统展现出 截然不同 的瞬态响应。

深入理解控制系统稳定性、可控性和可观性

参考:What is stability of control system ?

稳定性

稳定性指的是 闭环系统维持平衡状态或在被施加扰动或变化后返回该状态的能力

  1. 内部(Lyapunov)稳定性:考察系统状态 x ˙ = A x \dot x = Ax x˙=Ax没有外部输入下的固有动态

    • 渐近稳定(Asymptotic Stability):所有解 x ( t ) x(t) x(t) t → ∞ t\to\infty t 都趋于零。
    • 李雅普诺夫稳定(Lyapunov Stability):初始扰动足够小,状态始终保持在小范围内(但不一定趋于零)。
    • 不稳定:存在初始条件使得 ∥ x ( t ) ∥ → ∞ \|x(t)\|\to\infty x(t)
  2. BIBO 稳定性(Bounded-Input Bounded-Output Stability)
    有界输入和有界输出(BIBO)稳定性,通常称为输入输出稳定性,评估系统对有界(有限)输入的响应。如果一个系统    y ( t ) = G ( s ) u ( t ) \;y(t)=G(s)u(t) y(t)=G(s)u(t) 对于任何有界输入信号,其输出都保持有界,即不会无限增长,则该系统被认为是 BIBO 稳定的。

注意内部稳定 ⇔ BIBO 稳定(当系统是可控可观且最小实现时) ,否则两者可能不同。

  • 若系统包含不可控不可观测极点,这些极点在输入—输出传递函数中会被零点抵消,看不到(它们属于内部模态,不影响 BIBO)。
  • 最小实现(Minimal Realization)通过 minreal 或可控/可观子空间剔除这些“被抵消”的模态,使得传递函数中的极点集真正反映系统动态。

一个稳定的控制系统是指 输出保持有界,并且 即使受到外部力量或输入变化的干扰,也不会表现出不受控制的或振荡的行为

相位裕度增益裕度 等数学表达式用于量化 负反馈系统的稳定性

系统稳定性分类

判定条件(连续时域)根据
渐近稳定 ℜ ( s i ) < 0 \Re(s_i)<0 (si)<0极点位置
边界稳定 ℜ ( s i ) ≤ 0 \Re(s_i)\le0 (si)0,无重根极点位置
不稳定存在 ℜ ( s i ) > 0 \Re(s_i)>0 (si)>0 或重根极点位置
零点任意位置不影响稳定性
非最小相位零点在右半平 ℜ ( z j ) > 0 \Re(z_j)>0 (zj)>0零点位置影响瞬态
  • 极点 决定指数衰减/增长振荡,是真正的稳定性指标
  • 零点 决定输入—输出路径上的局部“抑制”相位/幅度畸变,但不改变是否收敛。

举个简单的例子,假设我们有一个传递函数 G ( s ) G(s) G(s)
G ( s ) = 1 s + a G(s)=\frac{1}{s+a} G(s)=s+a1
所以 − a -a a 就是我们的极点,为了理解这在时间域中意味着什么,取系统的逆拉普拉斯变换 L − 1 \mathcal L^{-1} L1,得到系统:
f ( t ) = e − a t u ( t ) f(t) = e^{-at}u(t) f(t)=eatu(t)
假设初始条件为 0,可以将 u u u 设为 1,我们的系统就变成了
f ( t ) = e − a t f(t) = e^{-at} f(t)=eat

稳定系统

如果 系统的传递函数的极点位于 s 平面的左侧,则认为 系统是稳定的。但随着极点靠近原点,系统的稳定性会降低。

脉冲响应显示,

  • 位于实轴上的极点具有过阻尼响应
  • 对于第二和第三象限的极点,瞬态响应是正弦波但指数衰减(欠阻尼响应)

接近原点的极被称为主导极(Dominant Poles)。因此,如果一个稳定系统的极位于 p 1 p_1 p1 p 2 p_2 p2(如上图),则 p 1 p_1 p1 被认为是该系统的主导极。 p 2 p_2 p2 被认为是该系统的非主导极(Non-dominant Poles)。在瞬态响应中,由于 p 2 p_2 p2 的贡献衰减得更快,响应将由主导极 p 1 p_1 p1 主导

边界稳定系统

一个系统如果其 极点位于虚轴上且不重复,则称为边界稳定

时域脉冲响应显示在恒定幅度下的振荡。初始条件(输入幅度)定义了幅度

不稳定系统

如果一个系统 在 s 平面的右侧有任何极点,则称为 不稳定。即使只有一个极点在右侧,也可能使系统不稳定。

输出随时间增长,振荡幅度随时间指数增长。在 实际系统中,振幅会上升到电源电压

可控性与可观性

理解并检验可控性与可观测性,是状态空间设计的第一步,也是保证反馈与估计器有效性的核心前提。

  1. 状态反馈:只有在可控系统上,才能通过反馈 u = − K x u=-Kx u=Kx 任意配置闭环极点。
  2. 状态估计:只有在可观测系统上,才能用观测器(如卡尔曼滤波、Luenberger 观测器)准确重建状态。
  3. 最小实现:不可控或不可观测的模态在输入—输出上无贡献,可用 minreal 删除,得到最小等价系统

考虑 LTI 系统的状态空间描述:

{ x ˙ ( t ) = A   x ( t ) + B   u ( t ) , y ( t ) = C   x ( t ) + D   u ( t ) , \begin{cases} \dot x(t) = A\,x(t) + B\,u(t),\\[5pt] y(t) = C\,x(t) + D\,u(t),\\ \end{cases} x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t)+Du(t),

其中 x ∈ R n x\in\mathbb R^n xRn u ∈ R m u\in\mathbb R^m uRm y ∈ R p y\in\mathbb R^p yRp

可控性(Controllability)

可控性(Controllability) 系统从 任意初始状态 x ( 0 ) = x 0 x(0)=x_0 x(0)=x0有限时间 T T T 内,通过恰当选择输入 u ( t ) u(t) u(t),能否将状态驱动到 任意目标状态 x ( T ) = x f x(T)=x_f x(T)=xf

判断条件:Kalman 可控性矩阵法则
C = [   B      A B      A 2 B    …    A n − 1 B   ] ∈ R n × n m \mathcal{C} = [\,B\;\;A B\;\;A^2B\;\dots\;A^{n-1}B\,]\in\mathbb R^{n\times nm} C=[BABA2BAn1B]Rn×nm

判据:若 r a n k ( C ) = n \mathrm{rank}(\mathcal C)=n rank(C)=n,则系统完全可控;否则存在不可控子空间。

定义可控 Gramian 判据:

W c ( T ) = ∫ 0 T e A τ B B ⊤ e A ⊤ τ   d τ . W_c(T) = \int_{0}^{T} e^{A\tau} B B^\top e^{A^\top \tau}\,\mathrm d\tau. Wc(T)=0TeAτBBeAτdτ.

  • 若对某 T > 0 T>0 T>0 W c ( T ) W_c(T) Wc(T) 是可逆的(正定矩阵),系统可控。
  • 物理/几何意义
    • 每个列向量 A k B A^k B AkB 表示“从输入作用一次、两次……到第 k + 1 k+1 k+1 次,对状态的影响方向”。
    • C \mathcal C C 的列空间若能张满整个状态空间,就说明输入足够“驱动”到任何方向。

可观测性(Observability)

可观性(Observability) 能否通过观测输出 y ( t ) y(t) y(t) 在有限时间内唯一重构系统初始状态 x ( 0 ) = x 0 x(0)=x_0 x(0)=x0

判断条件:Kalman 可观测性矩阵法则
O = [ C C   A C   A 2 ⋮ C   A n − 1 ] ∈ R p n × n \mathcal{O} = \begin{bmatrix} C\\ C\,A\\ C\,A^2\\ \vdots\\ C\,A^{n-1}\\ \end{bmatrix} \in\mathbb R^{pn\times n} O= CCACA2CAn1 Rpn×n

判据:若 r a n k ( O ) = n \mathrm{rank}(\mathcal O)=n rank(O)=n,则系统完全可观测;否则存在不可观测子空间。

定义可观测 Gramian 判据:
W o ( T ) = ∫ 0 T e A ⊤ τ C ⊤ C e A   τ   d τ . W_o(T) = \int_{0}^{T} e^{A^\top \tau} C^\top C e^{A\,\tau}\,\mathrm d\tau. Wo(T)=0TeAτCCeAτdτ.

  • 若对某 T > 0 T>0 T>0 W o ( T ) W_o(T) Wo(T) 是可逆的,系统可观测。
  • 物理/几何意义
    • O \mathcal O O 的行向量 C A k C A^k CAk 表示“状态向第 k k k 次演化后,对输出的映射”。
    • 若行空间张满状态空间,说明每个状态分量都能在某个时刻对输出产生可区分影响。

简单示例:二阶质量-弹簧系统

状态空间(无阻尼):

A = [ 0 1 − k / m 0 ] , B = [ 0 1 / m ] , C = [ 1 0 ] A=\begin{bmatrix}0&1\\-k/m&0\end{bmatrix},\quad\\[10pt] B=\begin{bmatrix}0\\1/m\end{bmatrix},\quad\\[10pt] C=\begin{bmatrix}1&0\end{bmatrix} A=[0k/m10],B=[01/m],C=[10]

  • 可控性 C = [ B ,    A B ] = [   [ 0 ; 1 / m ] ,    [ − k / m ; 0 ]   ] \mathcal C=[B,\;AB] = \bigl[\,[0;1/m],\;[-k/m;0]\,\bigr] C=[B,AB]=[[0;1/m],[k/m;0]] rank ⁡ ( C ) = 2 \operatorname{rank}(\mathcal C)=2 rank(C)=2,完全可控。
  • 可观测性 O = [ C C A ] = [   [ 1 , 0 ] ;    [ 0 , 1 ] ] \mathcal O=\begin{bmatrix}C\\CA\end{bmatrix} = \bigl[\,[1,0];\;[0,1]\bigr] O=[CCA]=[[1,0];[0,1]] rank ⁡ ( O ) = 2 \operatorname{rank}(\mathcal O)=2 rank(O)=2,完全可观测。

在数学上, ( A , B ) (A,B) (A,B) 的可控性等价于 ( A ⊤ , C ⊤ ) (A^\top, C^\top) (A,C) 的可观测性。这意味着可控性的分析工具几乎可以“对转”直接用于可观测性。

深入理解可控性和可观性

参考:A Conceptual Approach to Controllability and Observability

首先要了解你想要控制的系统,然后

  • 可以使用 执行器(Actuators) 来控制系统,执行器会传递控制力和能量。
  • 可以使用 传感器(Sensors) 测量系统的特定状态,例如电压、位置和温度,这些传感器会产生系统的输出。
  • 控制器(Controller) 设计归结为 如何利用传感器数据以及参考信号或某种命令来生成正确的执行器命令


如果你的系统没有合适的执行器来影响系统的正确部分,或者你没有安装合适的传感器来控制,那么你的控制器设计注定会失败。如果没有这两个参数,你就无法充分地影响系统。系统将无法控制,无法被观察。

因此,可控性和可观测性是系统如何与执行器和传感器协同工作的条件,它与特定的控制技术(如 PID 或极点配置)无关。

  • 可控性意味着存在控制信号,允许系统在有限的时间内 达到任何状态,这也称为可达性。所以可控性并不意味着必须维持该状态,而是指 可以达到该状态
  • 可观测性意味着所有状态(关心的状态)都可以从系统的输出中得知。也就是说,我们可以通过使用适当的传感器选择和传感器位置来了解系统中每个状态的值。

现在,让我们抛开约束的现实问题,看看汽车的高维状态空间。我们关心的状态是 XY 位置、XY 速度、偏航角和偏航率。

假设,

  • 传感器有车速表来测量速度,我们的眼睛用来评估其他所有东西。
  • 执行器是方向盘、油门和刹车踏板。

通过这套设置,我们 能够以给定的速度和偏航角控制汽车到特定位置

先来看一下 可观测性的概念,并从一个极端的例子开始,以帮助说明其重要性。

我们闭上眼睛,移除所有传感器信息,这会导致 C C C 矩阵归零。

基本上,我们无法知道任何状态。你不知道自己在哪里,行驶速度是多少,或者你面对的方向。在这个例子中,我们假设你感觉不到加速,也听不到引擎的声音,或者其他任何东西。

现在你仍然可以控制汽车,这并没有受到影响,因为你可以转动方向盘并踩下踏板。但是,由于 缺乏可观测性,你不能安全地驾驶它到达目的地。然而,根据我们之前的定义,汽车是可控的

相反的情况,我们移除方向盘和踏板,这会导致 B B B 矩阵归零。

想象一下,如果你撞到一大片冰面,汽车打滑,会发生什么情况。 并且不受转向或制动的影响,即使您无法影响汽车的状态,该 系统仍然是可观察的,因为您可以查看速度表,并且可以通过观察外部来跟踪您的位置。

可以看到拥有一个可观察和控制的系统是多么必要,它们协同工作,即使是 部分无法控制或无法观察 的系统仍然很麻烦,想象一下尝试在方向盘拆下的情况下驾驶汽车,即使仍然可以通过踏板改变其速度,但无法改变偏航状态。

如果任何一个临界状态无法控制,那么整个系统都被视为无法控制,如果任何一个临界状态无法观察,那么整个系统都被视为无法观察。

状态空间方程 —— 极点配置

极点配置(Pole Placement)完全状态反馈 是控制理论中的一种设计方法,用于通过 适当选择控制器增益,使得闭环系统的极点(也即系统特征方程的根)位于预定的位置,从而达到所需的动态性能和稳定性要求。极点的位置直接影响系统的响应特性,如稳定性、速度、阻尼等。

现在,极点配置本身在工业中并没有得到广泛的应用。你可能会更常使用其他方法,例如 LQR 或 H ∞ H_\infin H,但是极点配置值得花一些时间,因为它能让你更好地理解 使用状态空间方程进行反馈控制 的一般方法,这是进入其他方法的垫脚石。

参考:What is Pole Placement (Full State Feedback)

从 PID 到极点配置

首先,我们有一个具有输入 u u u 和输出 y y y 的装置(Plant),目标是 开发一个反馈控制系统将输出驱动到某个期望值

你可能熟悉的一种方法是将输出与参考信号进行比较以获得控制误差,然后你可以开发一个控制器,该控制器使用该误差来生成输入到装置的信号,目标是将误差驱动到零

但对于极点配置,我们将以不同的方式解决这个问题,我们 将反馈 状态向量中每个状态变量的值,而不是反馈输出 y y y

现在,我们要 假设我们知道每个状态的值,即使它不一定是输出 y y y 的一部分,然后 取状态向量并将其乘以一个由一堆不同增益值组成的矩阵 K K K,并 将结果从缩放的参考信号 K r K_r Kr 中减去,并将此结果直接输入到我们的装置中作为输入。

这里没有像在反馈结构框图中标记为控制器的块,图中绿圈整个部分都是控制器。

通过极点配置 可以 计算适当的增益矩阵 K K K确保系统稳定性,参考上的缩放项用于确保稳态误差性能在可接受范围内。

状态方程中的第一项 A x Ax Ax 捕获线性系统的动态特性,第二项是 系统如何响应输入,但系统中的能量如何存储和移动是由 A x Ax Ax 项捕获的,所以你可能会认为在控制器设计方面, A A A 矩阵有一些特殊之处。

是的!反馈控制器都必须 修改 A A A 矩阵才能改变系统的动态特性,尤其是在稳定性方面


A A A 矩阵的特征值是系统的极点,极点的位置决定了线性系统的稳定性,这是极点配置的关键,通过移动极点或闭环 A A A 矩阵的特征值来生成所需的闭环稳定性

举例说明特征值如何影响系统

看下面这个例子, A A A 矩阵转换为使用另一组状态变量来描述系统的矩阵,此转换是使用 变换矩阵 完成的,其列是 A A A 矩阵的特征向量


转换后我们得到的是一个修改后的 A A A 矩阵,由对角线上的复数特征值和其他所有位置的零组成,现在 绿色箭头所指的这两个模型代表完全相同的系统,它们具有 相同的特征值和相同的动态,只是第二个模型使用一组彼此独立变化的状态变量来描述

A A A 矩阵以对角形式写出时,很容易看出我们剩下的是 一组一阶微分方程,其中每个状态的导数仅受该状态的影响,而不受其他任何影响

很酷的是,像这样的微分方程的解的形式是 Z = C e λ T Z=C e^{\lambda T} Z=CeλT,其中 λ \lambda λ 是给定状态变量的特征值。

看几个不同的 λ \lambda λ 值,这样你就可以直观地看到 能量是如何根据复平面内特征值的位置变化的

  • 如果 λ \lambda λ 是负实数,那么这是一个稳定的特征值,因为解是 e λ t → 0 e^{\lambda t} \rightarrow 0 eλt0 t → ∞ t \rightarrow \infin t,任何初始能量都会随时间消散;
  • 但如果 λ \lambda λ 是正数,那么它是不稳定的,因为能量只会随时间增长;
  • 如果有一对 虚数特征值,那么这种模式下的能量将 振荡,因为 e e e 的虚数部分会产生 s i n sin sin c o s cos cos
  • 而实数和虚数的任意组合都会产生振荡和指数能量耗散。

极点配置如何移动

现在我们可以陈述我们要解决的问题了。如果对象的特征值位于复平面中不理想的位置,那么我们可以使用极点配置将它们移动到其他地方

当然,如果它们位于右半平面,这是理想的,因为它们不稳定,但非理想位置也可能意味着存在振荡,你想消除它,或者只是加速或减缓特定模式下的能量耗散。

我们现在可以讨论极点配置如何移动特征值。记住我们一开始画的控制器结构,这会导致输入 u = r K r − K x u= rK_r - K x u=rKrKx,其中 r K r r K_r rKr 是缩放参考值, K x Kx Kx 是状态向量,将其反馈 x x x 乘以增益矩阵 K K K

现在,如果 将此控制输入代入状态方程,将 闭合循环,得到以下状态方程,注意 A A A − B k -Bk Bk 都作用于状态向量,因此我们可以将它们组合起来,得到一个修改后的 A A A 矩阵,这是 闭环 A A A 矩阵,可以通过 选择适当的 K K K 来移动特征值

二阶系统示例

对于简单的系统,这很容易手动完成,让我们尝试一个具有单个输入的二阶系统的示例,

我们可以通过将 det ( A − λ I ) = 0 \text{det} (A - \lambda I) = 0 det(AλI)=0 来找到特征值,然后求解 λ \lambda λ,它们分别位于 − 2 -2 2 + 1 +1 +1 处,由于存在正的实特征值,其中一个模式将爆炸到无穷大,因此整个 系统是不稳定的

让我们使用极点配置设计一个反馈控制器,通过将不稳定的极点移动到左半平面来稳定该系统,我们的闭环 A A A 矩阵是 A − B k A-Bk ABk,增益矩阵 K K K 1 × 2 1\times 2 1×2,因为有一个输出 两个状态。

对得到的新 A − B k A-Bk ABk 进行,像以前一样求解特征值,得到这个特征方程,它是两个增益值的函数。现在 假设我们 希望闭环极点位于 − 1 -1 1 − 2 -2 2,这样特征方程为 λ 2 + 3 λ + 2 = 0 \lambda^2 + 3 \lambda + 2 = 0 λ2+3λ+2=0,此时很容易找到合适的 k 1 k_1 k1 k 2 k_2 k2,使这两个方程相等,我们 只需将系数设置为相等并求解,得到 k 1 = 2 k_1 = 2 k1=2 k 2 = 1 k_2 = 1 k2=1

就是这样,如果我们将这两个增益放在该系统的状态反馈路径中,它将在特征值位于 − 1 -1 1 − 2 -2 2 时保持稳定。

对于具有两个以上状态的系统,所涉及的数学运算开始变得难以理解。想法一样的,只是求解行列式变得不切实际,但在 MATLAB 几乎只用一个命令。

% Define state matrices
A = [0 1; 2 -1];
B = [1; 0];
C = [1 0];
D = 0;

% Create state space object
sys = ss(A, B, C, D);

% Check open loop eigenvalues
E = eig(A);

% Desired closed loop eigenvalues
P = [-2 -1];

% Solve for K using pole placement
K = place(A, B, P);

% Check for closed loop eigenvalues
Acl = A - B*K;
Ecl = eig(Acl);

% Create closed loop system
syscl = ss(Acl, B, C, D);

% Check step response
step(sys)

% Create colsed loop system
syscl = ss(Acl, B, C, D);

% Check step response
step(syscl)

可以预见的是,开环系统的阶跃是不稳定的

闭环系统的阶跃响应看起来要好得多,但它并不完美,而不是像我们预期的那样上升到 1,预计稳态输出仅为 0.5,这最终就是缩放项发挥作用的地方。

对状态方程进行设计前、仿真前的体检

控制系统设计前的状态方程体检

在拿到一组状态空间形式
{ x ˙ = A x + B u , y = C x + D u , \begin{cases} \dot x = A x + B u,\\ y = C x + D u, \end{cases} {x˙=Ax+Bu,y=Cx+Du,
之后,通常要从以下几方面去“体检”系统,才能保证后续的设计(比如状态反馈、观测器设计、滤波器设计等)有的放矢。

检查项判据结果不合格时的“补救”
稳定性 ℜ ( λ i ) < 0 \Re(\lambda_i)<0 (λi)<0设计状态反馈或动态补偿
可控性 rank ⁡ ( C ) = n \operatorname{rank}(\mathcal{C})=n rank(C)=n增执行器/降阶建模
可观测性 rank ⁡ ( O ) = n \operatorname{rank}(\mathcal{O})=n rank(O)=n增传感器/降阶建模
最小实现同时可控可观Kalman 分解去冗余
适定性 & 直通 D D D 是否零,因果性检查置零 D D D 或加滤波

1. 稳定性(Stability)

  • 检查项:矩阵 A A A 的特征值 { λ i } \{\lambda_i\} {λi}
  • 原理:连续系统内在稳定要求所有 ℜ ( λ i ) < 0 \Re(\lambda_i)<0 (λi)<0(特征值的实部);若有 ℜ ( λ i ) > 0 \Re(\lambda_i)>0 (λi)>0,系统本身就不稳定;若有 ℜ ( λ i ) = 0 \Re(\lambda_i)=0 (λi)=0,则可能边界稳定或振荡。
  • 调整方法
    • 若不稳定,设计 状态反馈 u = − K x u=-Kx u=Kx,使闭环矩阵 A c l = A − B K A_{\rm cl}=A-BK Acl=ABK 的特征值都在左半平面(如极点配置法、LQR)。
    • 若有边界或振荡模态,可适当加阻尼或配合动态补偿器。

举个例子,假设状态方程 A A A所有特征值 λ i λ_i λi 的实部均为负(最小约 − 0.000002326 −0.000002326 0.000002326,最大约 − 9.41 × 10 4 −9.41×10^4 9.41×104),则 系统渐近稳定
原理:连续时间系统 x ˙ = A x \dot x = A x x˙=Ax 的解为 x ( t ) = e A t   x ( 0 )   x(t) = e^{A t}\,x(0)\, x(t)=eAtx(0),当且仅当 ∀ i , ℜ ( λ i ) < 0 ∀i, \Re( λ_i)<0 i,(λi)<0 时, e λ i t → 0 e^{\lambda_i t}\to0 eλit0 t → + ∞ t→+∞ t+),系统才是渐近稳定的。
公式:特征方程 det ⁡ ( λ I − A ) = 0 \det(\lambda I - A)=0 det(λIA)=0,解得 λ i λ_i λi,即系统模式指数。

2. 可控性(Controllability)

  • 检查项:构造可控性矩阵
    C = [ B ,    A B ,    A 2 B ,    … ,    A n − 1 B ] , \mathcal{C} = \bigl[B,\;A B,\;A^2 B,\;\dots,\;A^{n-1}B\bigr], C=[B,AB,A2B,,An1B],
    并检验 rank ⁡ ( C ) = ? n \operatorname{rank}(\mathcal{C})\stackrel{?}{=}n rank(C)=?n
  • 原理:若 rank ⁡ ( C ) < n \operatorname{rank}(\mathcal{C})<n rank(C)<n,说明存在状态方向无法被输入 u u u 影响,这些“死角”无法通过任何输入到达。
  • 调整方法
    • 硬件层面:增设或重新布置执行器,使得原来不可控的模态得以激励。
    • 软件层面:对可控子空间进行降阶,只保留可控部分作为等价最小模型;对不可控部分只能当做扰动或噪声处理。

3. 可观测性(Observability)

  • 检查项:构造可观测性矩阵
    O = [ C C A C A 2 ⋮ C A n − 1 ] , \mathcal{O} = \begin{bmatrix}C\\C A\\C A^2\\\vdots\\C A^{n-1}\end{bmatrix}, O= CCACA2CAn1 ,
    并检验 rank ⁡ ( O ) = ? n \operatorname{rank}(\mathcal{O})\stackrel{?}{=}n rank(O)=?n
  • 原理:若 rank ⁡ ( O ) < n \operatorname{rank}(\mathcal{O})<n rank(O)<n,说明存在状态方向即便变化也无法从输出 y y y 中分辨。
  • 调整方法
    • 硬件层面:增设或调整传感器,使得原来不可观测的模态被量测。
    • 软件层面:对可观测子空间降阶处理,只建模可观测部分;对不可观测部分用作未建模动态。

4. 最小化/约简(Minimal Realization)

  • 检查项:可控性 & 可观测性都要满秩;若有任一不满秩,系统不是 最小实现
  • 原理:最小实现去掉了所有不可控或不可观测模态,模型既不冗余,也不过度简化;可提高数值稳健性和计算效率。
  • 调整方法
    • 利用 Kalman 分解Gramian 分解 将系统分解成可控可观、可控不可观、不可控可观、不可控不可观四个子系统,丢弃不可控或不可观子系统。

5. 传递函数的适定性(Properness & Direct‐Feedthrough)

  • 检查项:矩阵 D D D 是否为零?系统阶次是否满足因果性(严格、准)。
  • 原理
    • D ≠ 0 D\ne0 D=0,意味着有“直通”通道,输入 u u u 能瞬时影响输出 y y y
    • 对于物理系统,过高的直通增益可能导致测量噪声放大或数值震荡。
  • 调整方法
    • 若不希望直通,令 D = 0 D=0 D=0
    • 或在设计控制器时,加低通滤波器以抑制高频噪声。

掌握了以上原则和方法,就能针对任意给定的 A , B , C , D A,B,C,D A,B,C,D 矩阵,先“体检”再“对症下药”,确保后续控制器/观测器设计既可行又高效。

仿真数值求解前的状态方程体检

体检项检查方法结果不合格时的调整思路
刚性(Stiffness) τ max ⁡ τ min ⁡ \frac{\tau_{\max}}{\tau_{\min}} τminτmax、谱半径比 max ⁡ ∣ ℜ ( λ ) ∣ min ⁡ ∣ ℜ ( λ ) ∣ \displaystyle\frac{\max | \Re(\lambda) |}{\min | \Re(\lambda)|} min∣ℜ(λ)max∣ℜ(λ)、Gramian 条件数、仿真对比多时尺度分离、模型降阶、数值正则化、选用刚性求解器

刚性(Stiffness)

除了检查极点实部以外,“刚性”(stiffness)确实也是“体检”时一个值得关注的性质,尤其在以下两种情境中尤为重要:

  1. 数值仿真与离散化:当系统的模态(极点)时间常数相差悬殊时,常规的数值积分器(如显式欧拉、Runge–Kutta)往往要求 极小的步长 才能保持稳定 和精度。
  2. 分层/多时尺度控制:当某些模态非常快(刚性模态),而另一些模态又很慢,如果不加以分离、近似或降阶,设计出的控制器可能在工程上既难以实现又成本高昂。

刚性的定义与量化:

  1. 最大时间常数和最小非零时间常数之比

    • 定义每个模态的时间常数为:
      τ i = 1 ∣ ℜ ( λ i ) ∣ , 仅当  ℜ ( λ i ) ≠ 0 \tau_i = \frac{1}{|\Re(\lambda_i)|}, \quad \text{仅当 } \Re(\lambda_i) \neq 0 τi=∣ℜ(λi)1,仅当 (λi)=0
    • 刚性比(stiffness ratio): R = τ max ⁡ τ min ⁡ R = \frac{\tau_{\max}}{\tau_{\min}} R=τminτmax
    • R ≫ 1 R \gg 1 R1(比如 R > 10 3 R > 10^3 R>103 10 5 10^5 105),系统被认为是刚性的(stiff system),可能导致显式数值解发散。

    举个例子,假设系统所有特征值的实部,最小约 − 0.000002326 −0.000002326 0.000002326,最大约 − 9.41 × 10 4 −9.41×10^4 9.41×104,则根据极大/极小特征值可以估算,

    • 最快模态约 τ min ⁡ ≈ 1 / 9.4 × 10 4 ≈ 10 − 5   s \tau_{\min}\approx 1/9.4\times10^4 \approx10^{-5}\,\mathrm{s} τmin1/9.4×104105s
    • 最慢模态约 τ max ⁡ ≈ 1 / 2.3 × 10 − 6 ≈ 4 × 10 5   s \tau_{\max}\approx 1/2.3\times10^{-6}\approx4\times10^5\,\mathrm{s} τmax1/2.3×1064×105s(数天量级)
    • 刚性指标(stiffness ratio) τ max ⁡ τ min ⁡ ∼ 10 10 ≫ 1 , \frac{\tau_{\max}}{\tau_{\min}}\sim10^{10}\gg1, τminτmax10101,
      则为高刚性系统。
  2. 频谱比(Spectral Radius Ratio)
    Stiffness Ratio = max ⁡ i ∣ ℜ ( λ i ) ∣ min ⁡ j   :   ℜ ( λ j ) ≠ 0 ∣ ℜ ( λ j ) ∣ . \text{Stiffness Ratio} = \frac{\max_i \bigl|\Re(\lambda_i)\bigr|}{\min_{j\,:\,\Re(\lambda_j)\neq0}\bigl|\Re(\lambda_j)\bigr|}. Stiffness Ratio=minj:(λj)=0 (λj) maxi (λi) .
    如果该比值远大于 10 2 10^2 102 10 3 10^3 103,就常被认为是刚性系统。

  3. 条件数(Condition Number)

    • 状态转移矩阵 e A t e^{At} eAt 的条件数高,或者

    • Gramian 矩阵(可控 Gramian W c W_c Wc、可观测 Gramian W o W_o Wo)的条件数高,
      都会导致求解或降阶过程中的数值不稳定。

      检查 Gramian 条件数
      W c = ∫ 0 T e A t B B T e A T t   d t , W o = ∫ 0 T e A T t C T C e A t   d t . W_c = \int_0^{T} e^{A t} B B^T e^{A^T t} \,\mathrm{d}t,\quad W_o = \int_0^{T} e^{A^T t} C^T C e^{A t} \,\mathrm{d}t. Wc=0TeAtBBTeATtdt,Wo=0TeATtCTCeAtdt.
      计算 κ ( W c ) \kappa(W_c) κ(Wc) κ ( W o ) \kappa(W_o) κ(Wo)

  4. 也可以从 数值仿真试验 结果来减压系统的刚性,

    • 用不同步长的显式积分器(如 RK4)模拟系统,观察是否存在“爆渐”(数值不稳定)或“过度阻尼”现象;
    • 与刚性专用积分器(如隐式 Euler、Gear’s 方法)对比,若差异显著,则说明刚性问题突出。

刚性不合格时的“对症下药”

  1. 多时尺度分离
    • 奇异摄动方法将系统拆成 “快子系统”“慢子系统”,分别降阶或分别设计控制器,再用级联/分层方式集成。
    • 模态分解:对极值很大的快模态进行定性建模或近似“瞬态消隐”,只对慢模态进行精细设计。
  2. 模型降阶
    • 平衡截断(Balanced Truncation):在保留主要能量传递通道的同时,丢弃对输入—输出映射影响非常小的高频刚性模态
    • 投影法:如 Krylov 子空间法,针对高频模态快速收敛的特点进行降阶。
  3. 数值处理
    • 状态变量尺度变换:对状态量做合理归一化或尺度分割,减少数值差异;
    • 雅可比正则化:在离散化或仿真时,对 A A A 做小量级正则化(如 A + ϵ I A+\epsilon I A+ϵI)以改善条件数。
  4. 控制器/观测器设计
    • 含刚性补偿器:设计时考虑频率加权,避免对极高频模态过度放大;
    • 选用刚性友好算法:如 LQR 中的 Ricatti 方程求解时,选用隐式方法及高精度求解器。
  5. 仿真工具选型
    • 明确使用“刚性求解器”(Stiff Solver),例如 MATLAB 中的 ode15sode23s,或仿真平台相应的隐式积分器。

将“刚性”纳入常规体检清单后,可在早期就识别潜在的数值稳定性与实现难题,从而在模型简化、控制器设计、仿真与在线实现等各个环节做到有的放矢。这样,既能保证系统性能,又可提高仿真与在线运行的鲁棒性和效率。

MATLAB 绘制极点-零点分布

%% 1. 构造系统并算极零点
sys = ss(A,B,C,D);  % 自定义的 A,B,C,D

%% 2. 绘制极点-零点图
figure
hold on;
pzmap(sys)
grid on

%% 3. 叠加阻尼比和自然频率线
sgrid

底图并不是一个“正方形”网格,而是典型的 s‑平面(damping–frequency)网格:

  • 圆弧(constant natural frequency)
  • 放射状直线(constant damping ratio)
  • x 轴、y 轴刻度不等距(非等高宽)

MATLAB 自带了 Control Toolbox 里的 sgrid 函数,专门用来往 s‑平面上画这两种网格。

从极点(Poles)和传输零点(Transmission Zeros)结果来看,主要暴露出以下几个问题:

  1. 高阶系统中存在大量“重合”极点–零点
    • 除了一个 幅值约为 − 9.4108 × 10 4 -9.4108\times10^4 9.4108×104 的快模,几乎所有极点都落在 − 0.0 … ∼ − 1.0 -0.0…\sim-1.0 0.01.0 之间,而且传输零点列表里几乎“抄”了一份极点(同样的值、同样的幅度)——这说明在这些频率点上存在 成对的极点–零点
    • 这正是典型的 “极点–零点相互抵消”(pole–zero cancellation):系统在这些模态上既有极点又有相同位置的零点,结果在开环或闭环响应中,这部分动态都会被“抹去”,对外表现得像不存在这些模态
  2. 模型“非最简”/不必要的高阶
    • 当从某个更复杂的结构(比如带有很多中间变量、滤波器、观测器等)直接写出 A , B , C , D A,B,C,D A,B,C,D 后,很容易出现 非最小实现(non‑minimal realization):可控但不可观测或可观测但不可控的状态分量,被保留在模型里,却对输入—输出响应没有实质贡献
  3. 大量频率接近零的模式
    • 大量接近 0 (如 − 0.0000 -0.0000 0.0000)的极点和零点,往往是数值精度或模型降阶/近似造成的“伪模式”——物理上并不存在这么多几乎为零的慢模态(“伪零点”)。

当看到:

极点: -1.8423, -1.7864, …  
零点: -1.8423, -1.7864,

一一对应,说明:

  1. 这些 极点–零点对 在传递函数里互相抵消,对输入–输出动态没有贡献
  2. 换言之,这部分“内部模态”虽然存在于 A A A 矩阵,却对任何外部输入 u u u 并不会在输出 y y y 上留下痕迹——要么是 不可观测,要么是 不可控,或者两者兼有。

真正对系统输入—输出贡献 的是 被保留下来的“非抵消”极点/零点,对这些模式之外的部分,可以 用 minreal 或降阶方法清理掉。这样既让模型更简洁,也能避免后续控制器设计中出现不必要的麻烦。

  1. 最小化模型: 在 MATLAB 中对状态空间对象调用

    sys_min = minreal(sys);
    

    它会自动做可控性/可观测性分析,去掉所有“成对抵消”的极点–零点,得到最简的等价系统。

  2. 数值条件检查:

    • ctrb(A,B)obsv(A,C) 查看可控/可观矩阵秩,确认是否存在不可控或不可观测子空间。
    rank(ctrb(A,B)), rank(obsv(A,C))
    
    • 若秩小于状态维数,说明有冗余(不可控或不可观),需要先做可控/可观子空间分解。
  3. 平衡降阶(可选) 如果想在 保留主要特性下更进一步降阶,可用

    sys_bal = balreal(sys);
    sys_red = modred(sys_bal, indices_to_remove);
    

    或直接

    sys_red = balred(sys, r);
    

    用 Hankel 奇异值自动去掉贡献最小的模式。

  4. 避免数值爆炸/精度问题:

    • 如果原始 A A A 矩阵元素跨度极大(如首个极点是 10 4 10^4 104,其余都是 10 − 1 10^{-1} 101 量级),建议先 缩放状态 或做 坐标变换 让系统矩阵条件数改善。
    • 保持同一数量级,有助于数值算法(计算特征值、零点)更稳健。

注意,上图是 SISO 系统。
如果是 MIMO 系统,可能会像下图一样只有极点而没有零点,

  • 如果是 MIMO 且模型里 D = 0 D = 0 D=0,说明所有的传递函数都是“真分式”(分子次数 < 分母次数)。此时系统的相对阶数 > 0,多出来的阶数会在 ∞ \infin 处“贡献”零点,不会在复平面上出现任何有限零点
  • MATLAB 的 tzero(sys) 针对的是 整块 MIMO 系统的 传输(invariant)零点,它反映了系统输入到输出的“整体”可控可观特性。即便每个通道(SISO)都有零点,只要它们在组合里没有成为“传输零点”,tzero 也不会把它们列出来

✨ 为什么 可以 模型降阶?

含义
极点 = 模态存在系统里有这个模式
零点 = 模态抵消输入–输出通路上不允许这个模式传播
极点 = 零点模态虽然存在,但完全被抑制,对输入输出行为 没有任何影响
极点 = 零点 为什么无用?不影响响应、不影响控制器、不改变系统性能,删掉更高效
  • 移除“丝毫不影响”输入–输出性能的那些模式,让模型变得更小、更高效
  • 降阶后剩下的极点,才是真正决定系统对外(关心的那个输入–输出通道)动力学性能的 “有效模式”;
  • 去除无用模式还有利于控制器设计(避免去补偿那些根本检测不到或控制不到的动态)。

举一个简单的例子,考虑一个 SISO 系统,传递函数为:

G ( s ) = ( s + 3 ) ( s + 1 ) ( s + 3 ) G(s) = \frac{(s + 3)}{(s + 1)(s + 3)} G(s)=(s+1)(s+3)(s+3)
可以对上面的式子做约分,化简后得到的系统:
G 简化 ( s ) = 1 s + 1 G_{\text{简化}}(s) = \frac{1}{s + 1} G简化(s)=s+11

  • 原本极点有 { − 1 , − 3 } \{-1,-3\} {1,3},零点有 { − 3 } \{-3\} {3}
  • 因为 − 3 -3 3 频率上“零点–极点成对抵消”,最终等价于一个一阶系统 1 / ( s + 1 ) 1/(s+1) 1/(s+1)

为什么这个抵消的模态没用了?

可以从时域和频域两个角度来理解。

  • 时域角度理解
    系统本来包含模态:

    • e − t e^{-t} et(极点 -1)→ 留下来了
    • e − 3 t e^{-3t} e3t(极点 -3)→ 被零点抵消了!

    也就是说:系统本来有两个动态过程在响应输入时发生,一个快的、一个慢的。但是由于 s = − 3 s = -3 s=3 同时是零点,所以该模态虽然 系统内部存在,但它 对输出的影响被输入路径“抵消”掉了
    想象成:这个模态存在于系统内部,但输入 u u u 不会激发它,输出 y y y 也看不到它

  • 频域角度理解
    来看系统的频率响应:

    • 零点是让系统在某个频率 完全无响应
    • 极点是系统在某个频率 响应增强(或延迟)。

    极点和零点相同 时,这两种影响完全抵消掉了。在 s = − 3 s = -3 s=3 附近,刚想要“振起来”,系统立刻把那一部分“消掉”,结果就是没有那部分响应。

在控制系统里,我们只关心 输入 u u u输出 y y y 的通路。如果某个内部状态的动态,输入激不起它,输出也看不到它,那它对控制器设计来说是 完全无用的:无需观测它、无需控制它、不影响性能、不影响稳定性。

所以在建模或简化模型时,完全可以删掉它

根据极点-零点分布进行降阶

给定一个高阶状态空间模型,希望得到一个低阶模型,其动态行为在指定频段内近似等效(比如频率响应、阶跃响应等尽量接近)。

构造降阶模型常用的方法

  • 平衡截断(Balanced Truncation):虽然不是直接基于极点零点,但会保留能量最大的状态,对极点-零点分布也有一定反映。数值稳定。

  • 极点匹配法(Pole-Matching Reduction):直接构造一个只有保留下来的极点的低阶模型。可通过匹配保留的极点和零点重新构建系统。

  • 共振频率保留法(尤其适合振动系统):保留低频处主导极点,丢弃高频极点。

误差评估

  • 频率响应对比(Bode 图)
  • 时域响应对比(阶跃响应、冲激响应)
  • H ∞ \mathcal{H}_\infty H H 2 \mathcal{H}_2 H2范数误差评估

  • 状态空间模型
    • 向量和矩阵含义
    • 状态空间模型求解
      • 解析解
    • 状态空间到传递函数
    • 状态空间中的极点
    • 示例
      • 一个简单的例子
      • 从系统图(电气系统)开发状态空间模型
      • 二阶质量-弹簧系统(second-order mass-spring system)
  • 深入理解系统的极点-零点
    • 极点(Pole)的物理/数学含义
    • 零点(Zero)的物理/数学含义
      • SISO 情况下的“零点”
      • MIMO/状态空间下的“传输零点”
      • 举个例子帮助理解
  • 深入理解控制系统稳定性、可控性和可观性
    • 稳定性
      • 系统稳定性分类
        • 稳定系统
        • 边界稳定系统
        • 不稳定系统
    • 可控性与可观性
      • 可控性(Controllability)
      • 可观测性(Observability)
    • 简单示例:二阶质量-弹簧系统
    • 深入理解可控性和可观性
  • 状态空间方程 —— 极点配置
    • 从 PID 到极点配置
    • 举例说明特征值如何影响系统
    • 极点配置如何移动
    • 二阶系统示例
  • 对状态方程进行设计前、仿真前的体检
    • 控制系统设计前的状态方程体检
      • 1. 稳定性(Stability)
      • 2. 可控性(Controllability)
      • 3. 可观测性(Observability)
      • 4. 最小化/约简(Minimal Realization)
      • 5. 传递函数的适定性(Properness & Direct‐Feedthrough)
    • 仿真数值求解前的状态方程体检
      • 刚性(Stiffness)
  • MATLAB 绘制极点-零点分布
    • ✨ 为什么 可以 模型降阶?
    • 根据极点-零点分布进行降阶

在数学表达中,向量 通常用 粗体小写字母 表示,矩阵 通常用 粗体大写字母 表示。

状态空间模型

状态空间模型(State‐Space Model)将一个 线性时不变系统 (LTI, Linear Time-Invariant)的“内部状态”、“外部输入”与“观测输出”用向量和矩阵统一刻画:

x ˙ ( t ) = A   x ( t ) + B   u ( t ) , y ( t ) = C   x ( t ) + D   u ( t ) \dot{\boldsymbol{x}}(t) = \boldsymbol{A}\,\boldsymbol{x}(t) + \boldsymbol{B}\,\boldsymbol{u}(t), \\[10pt] \boldsymbol{y}(t) = \boldsymbol{C}\,\boldsymbol{x}(t) + \boldsymbol{D}\,\boldsymbol{u}(t) x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t)+Du(t)

向量和矩阵含义

  • 状态向量 x ( t ) ∈ R n \boldsymbol{x}(t)\in\mathbb{R}^n x(t)Rn:是一个 n × 1 n \times 1 n×1列向量,表示 系统的内部状态 n n n 是系统的状态维数。
    • 示例:如果是一个 2D 机械系统,可能有 x = [ x 1 x 2 ] = [ 位置 速度 ] \displaystyle \boldsymbol{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} \text{位置} \\ \text{速度} \end{bmatrix} x=[x1x2]=[位置速度]
    • 把系统此刻所有“内部变量”打包,就像把弹簧—质量系统的位置和速度一起装进一个向量里,它包含了系统当前的完整信息。
  • 输入向量 u ( t ) ∈ R m \boldsymbol{u}(t)\in\mathbb{R}^m u(t)Rm:是一个 m × 1 m \times 1 m×1列向量,表示系统的 控制输入外部输入 m m m 是输入的数量。
    • 示例:在机械系统里是推力/扭矩,在电路里是电压/电流。
  • 输出向量 y ( t ) ∈ R p \boldsymbol{y}(t)\in\mathbb{R}^p y(t)Rp:是一个 p × 1 p \times 1 p×1列向量,表示系统的 可观测输出 p p p 是输出的数量。
    • 这是能观测到、关心的量,比如只想测位置,就由 C C C 从完整状态中“筛”出位置分量。
    • 示例:在温度控制系统中, y \boldsymbol{y} y 可能就是测得的温度。

系统矩阵的含义

  • 状态矩阵 A ∈ R n × n A\in\mathbb{R}^{n\times n} ARn×n:描述 系统状态之间的相互关系,决定了系统的固有动态特性(例如,极点位置、稳定性、振荡频率等)

    • 描述系统“自身演化”规律——若无任何外部输入,状态如何因内部动力(弹簧、阻尼、惯性等)自然变化
  • 输入矩阵 B ∈ R n × m B\in\mathbb{R}^{n\times m} BRn×m:描述了 不同输入 u \boldsymbol{u} u 对不同状态 x \boldsymbol{x} x 的影响权重

    • 告诉我们 每一路外部输入 u u u 如何“分配”到各个状态通道上,产生加速或力矩。
  • 输出矩阵 C ∈ R p × n C\in\mathbb{R}^{p\times n} CRp×n:描述了 状态向量 x \boldsymbol{x} x 与输出 y \boldsymbol{y} y 的关系,将完整内部状态投影到我们关心的输出上。

    • C \boldsymbol{C} C 矩阵的每一行对应一个输出变量与所有状态变量之间的关系。
    • i i i 行即对应第 i i i 个输出如何由所有状态分量线性组合而来
  • 直接转换(或馈通)矩阵 D ∈ R p × m \boldsymbol{D}\in\mathbb{R}^{p\times m} DRp×m 描述了 输入 u ( t ) \boldsymbol{u}(t) u(t)输出 y ( t ) \boldsymbol{y}(t) y(t)直接传递(feedthrough)或 立即响应,即在同一时刻没有经过“状态”延迟就能影响输出。

    • 严格本性(Strictly Proper)系统 D = 0 \boldsymbol{D} = \boldsymbol{0} D=0,系统的输出必须“经过”状态动态才能对输入做出反应,传递函数 G ( s ) = C ( s I − A ) − 1 B \displaystyle G(s)=C(sI-A)^{-1}B G(s)=C(sIA)1B 在高频衰减为零。
    • 本性(Proper)系统 D ≠ 0 \boldsymbol{D}\neq \boldsymbol{0} D=0,传递函数 G ( s ) = C ( s I − A ) − 1 B + D \displaystyle G(s)=C(sI-A)^{-1}B + D G(s)=C(sIA)1B+D 在高频极限时收敛到 D D D
    • 校正/补偿在某些控制器结构(如前馈补偿)中,设计一个合适的 D D D 项可以预先补偿已知的输入—输出静态增益,改善稳态性能。

状态空间模型 天然支持 MIMO直接处理多输入多输出,不需像传递函数那样一一写多条关系。

x ˙ = A x + B u \dot x=A x + B u x˙=Ax+Bu 明确指出“当前状态”与“当前输入”如何决定“当前变化率”,可方便处理非零初始条件。

状态空间模型求解

  • 解析解(Analytical Solution)指的是 通过代数运算、积分、微分等数学公式推导出的显式表达式,结果通常是一个函数形式,写成 关于时间 t t t 的表达式

    • 明确的闭式解表达式,适用于纸面推导、公式分析;
    • 通常 无法用于复杂非线性或高维系统
  • 数值解(Numerical Solution)是通过 计算机程序/数值算法(如欧拉法、龙格库塔法、有限差分等),在 离散时间点上近似计算 出的结果。

    • 存在误差,精度依赖于时间步长和算法阶数;
    • 可以应用于几乎任意微分方程(尤其是非线性或高维系统),实际工程中使用最多。
解法类型状态解 x ( t ) \boldsymbol{x}(t) x(t) u ( t ) = u 0 \boldsymbol{u}(t) = \boldsymbol{u}_0 u(t)=u0 为常值向量时输出解 y ( t ) \boldsymbol{y}(t) y(t)
解析解 x ( t ) = e A t x 0 + ∫ 0 t e A ( t − τ ) B u ( τ ) d τ \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \int_0^t e^{\boldsymbol{A}(t - \tau)} \boldsymbol{B} \boldsymbol{u}(\tau) d\tau x(t)=eAtx0+0teA(tτ)Bu(τ)dτ x ( t ) = e A t x 0 + A − 1 ( e A t − I ) B u 0 \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \boldsymbol{A}^{-1}(e^{\boldsymbol{A}t} - \boldsymbol{I})\boldsymbol{B}\boldsymbol{u}_0 x(t)=eAtx0+A1(eAtI)Bu0 y ( t ) = C x ( t ) \boldsymbol{y}(t) = \boldsymbol{C} \boldsymbol{x}(t) y(t)=Cx(t)
数值解 x ( t i ) \boldsymbol{x}(t_i) x(ti) 为一组离散时间点的近似值(如通过 Runge-Kutta 计算) x ( t i + 1 ) ≈ x ( t i ) + h ⋅ ( A x ( t i ) + B u 0 ) \boldsymbol{x}(t_{i+1}) \approx \boldsymbol{x}(t_i) + h \cdot \left( \boldsymbol{A}\boldsymbol{x}(t_i) + \boldsymbol{B} \boldsymbol{u}_0 \right) x(ti+1)x(ti)+h(Ax(ti)+Bu0) y ( t i ) = C x ( t i ) \boldsymbol{y}(t_i) = \boldsymbol{C} \boldsymbol{x}(t_i) y(ti)=Cx(ti)

解析解

给定方程: x ˙ = A x + B u \dot{\boldsymbol{x}} = \boldsymbol{A} \boldsymbol{x} + \boldsymbol{B} \boldsymbol{u} x˙=Ax+Bu,假设初始条件为 x ( 0 ) = x 0 \boldsymbol{x}(0) = \boldsymbol{x}_0 x(0)=x0解析解公式
x ( t ) = e A t x 0 + ∫ 0 t e A ( t − τ ) B u ( τ )   d τ \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \int_0^t e^{\boldsymbol{A}(t - \tau)} \boldsymbol{B} \boldsymbol{u}(\tau) \, d\tau x(t)=eAtx0+0teA(tτ)Bu(τ)dτ
这是利用矩阵指数法推导出来的解析解,前项是自由响应,后项是受控响应(强迫响应)

  • 自然演化 e A t x 0 e^{\boldsymbol{A}t} \boldsymbol{x}_0 eAtx0:若从 t = 0 t=0 t=0 起无任何输入,系统怎样“凭借自身动力”演化到 t t t
  • 卷积累积 ∫ 0 t e A ( t − τ ) B u ( τ )   d τ \int_0^t e^{\boldsymbol{A}(t - \tau)} \boldsymbol{B} \boldsymbol{u}(\tau) \, d\tau 0teA(tτ)Bu(τ)dτ:把每一刻的输入 u ( τ ) \boldsymbol{u}(\tau) u(τ) 按“系统记忆” e A ( t − τ ) e^{\boldsymbol{A}(t - \tau)} eA(tτ) 加权累积,得出总效果。

其中,矩阵指数定义为
e A t    =    ∑ k = 0 ∞ ( A   t ) k k ! . e^{\boldsymbol{A}t}\;=\;\sum_{k=0}^\infty \frac{(\boldsymbol{A}\,t)^k}{k!}. eAt=k=0k!(At)k.

如果 u ( t ) \boldsymbol{u}(t) u(t) 是常值输入(即 u ( t ) = u \boldsymbol{u}(t) = \boldsymbol{u} u(t)=u 是一个常向量),那么数值解为
x ( t ) = e A t x 0 + ( ∫ 0 t e A ( t − τ ) d τ ) B u = e A t x 0 + A − 1 ( e A t − I ) B u \boldsymbol{x}(t) = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \left( \int_0^t e^{\boldsymbol{A}(t - \tau)} d\tau \right) \boldsymbol{B} \boldsymbol{u} = e^{\boldsymbol{A}t} \boldsymbol{x}_0 + \boldsymbol{A}^{-1} \left( e^{\boldsymbol{A}t} - \boldsymbol{I} \right) \boldsymbol{B} \boldsymbol{u} x(t)=eAtx0+(0teA(tτ)dτ)Bu=eAtx0+A1(eAtI)Bu

注意事项:如果 A A A 不可逆(奇异矩阵),不能直接用 A − 1 A^{-1} A1。这时要用其他方式近似积分,比如数值方法(Runge-Kutta)

推导过程

  • 齐次解(零输入响应):当 u = 0 u=0 u=0 时,系统变为齐次
    x ˙ h ( t ) = A x h ( t ) \dot{x}_h(t) = A x_h(t) x˙h(t)=Axh(t) 其解为: x h ( t ) = e A t x ( 0 ) x_h(t) = e^{A t} x(0) xh(t)=eAtx(0)
  • 特解(零状态响应):由于 u u u 为常量,特解为,
    x p ( t ) = ∫ 0 t e A ( t − τ ) B u   d τ x_p(t) = \int_0^t e^{A (t - \tau)} B u \, d\tau xp(t)=0teA(tτ)Budτ 因为 A A A 是常矩阵,且 u u u 为常量,可得: x p ( t ) = ( ∫ 0 t e A ( t − τ )   d τ ) B u = ( ∫ 0 t e A s   d s ) B u x_p(t) = \left( \int_0^t e^{A (t - \tau)} \, d\tau \right) B u = \left( \int_0^t e^{A s} \, ds \right) B u xp(t)=(0teA(tτ)dτ)Bu=(0teAsds)Bu
    s = t − τ s = t - \tau s=tτ,则: ∫ 0 t e A s   d s = A − 1 ( e A t − I ) \int_0^t e^{A s} \, ds = A^{-1} (e^{A t} - I) 0teAsds=A1(eAtI) 因此,特解为: x p ( t ) = A − 1 ( e A t − I ) B u x_p(t) = A^{-1} (e^{A t} - I) B u xp(t)=A1(eAtI)Bu
  • 将齐次解与特解相加,得到总解
    x ( t ) = x h ( t ) + x p ( t ) = e A t x ( 0 ) + A − 1 ( e A t − I ) B u x(t) = x_h(t) + x_p(t) = e^{A t} x(0) + A^{-1} (e^{A t} - I) B u x(t)=xh(t)+xp(t)=eAtx(0)+A1(eAtI)Bu

状态空间到传递函数

我们从线性时不变系统(LTI)的状态空间模型开始:
x ˙ ( t ) = A x ( t ) + B u ( t ) , y ( t ) = C x ( t ) \dot{\boldsymbol{x}}(t) = \boldsymbol{A}\boldsymbol{x}(t) + \boldsymbol{B}\boldsymbol{u}(t), \quad \boldsymbol{y}(t) = \boldsymbol{C}\boldsymbol{x}(t) x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t)

零初始条件 x ( 0 ) = 0 \boldsymbol{x}(0) = \boldsymbol{0} x(0)=0,对上式两边进行 拉普拉斯变换 L { ⋅ } \mathcal{L}\{\cdot\} L{},就得到了:
s X ( s ) = A X ( s ) + B U ( s ) ⇒ ( s I − A ) X ( s ) = B U ( s ) s\boldsymbol{X}(s) = \boldsymbol{A} \boldsymbol{X}(s) + \boldsymbol{B} \boldsymbol{U}(s) \quad \Rightarrow \quad (s\boldsymbol{I} - \boldsymbol{A}) \boldsymbol{X}(s) = \boldsymbol{B} \boldsymbol{U}(s) sX(s)=AX(s)+BU(s)(sIA)X(s)=BU(s)

进一步得出:
X ( s ) = ( s I − A ) − 1 B U ( s ) \boldsymbol{X}(s) = (s\boldsymbol{I} - \boldsymbol{A})^{-1} \boldsymbol{B} \boldsymbol{U}(s) X(s)=(sIA)1BU(s)

再代入输出方程:
Y ( s ) = C X ( s ) = C ( s I − A ) − 1 B U ( s ) \boldsymbol{Y}(s) = \boldsymbol{C} \boldsymbol{X}(s) = \boldsymbol{C} (s\boldsymbol{I} - \boldsymbol{A})^{-1} \boldsymbol{B} \boldsymbol{U}(s) Y(s)=CX(s)=C(sIA)1BU(s)

我们将系统的 输入 U ( s ) \boldsymbol{U}(s) U(s)输出 Y ( s ) \boldsymbol{Y}(s) Y(s) 的关系写成如下形式:
G ( s ) = Y ( s ) U ( s ) = C ( s I − A ) − 1 B \boldsymbol{G}(s) = \frac{\boldsymbol{Y}(s)}{\boldsymbol{U}(s)} = \boldsymbol{C}(s\boldsymbol{I} - \boldsymbol{A})^{-1} \boldsymbol{B} G(s)=U(s)Y(s)=C(sIA)1B

这个 G ( s ) \boldsymbol{G}(s) G(s) 就是系统的 传递函数矩阵(Transfer Function Matrix),它描述了每个输入如何影响每个输出 —— 是一个 从频域角度刻画系统行为的工具

传递函数中的项 ( s I − A ) − 1 (s\boldsymbol{I} - \boldsymbol{A})^{-1} (sIA)1 是一个 矩阵的逆,它的存在依赖于:
det ⁡ ( s I − A ) ≠ 0 \det(s\boldsymbol{I} - \boldsymbol{A}) \neq 0 det(sIA)=0

因此,若 det ⁡ ( s I − A ) = 0 \det(s\boldsymbol{I} - \boldsymbol{A}) = 0 det(sIA)=0 就表示系统传递函数 G ( s ) \boldsymbol{G}(s) G(s) 的某些分量出现了 极点(即系统响应趋于无穷大)。这时对应的 s s s 值是 系统矩阵 A \boldsymbol{A} A 的特征值

状态空间中的极点

在状态空间模型中:

x ˙ = A x + B u , y = C x + D u \dot{\boldsymbol{x}} = \boldsymbol{A} \boldsymbol{x} + \boldsymbol{B} \boldsymbol{u},\quad \boldsymbol{y} = \boldsymbol{C} \boldsymbol{x} + \boldsymbol{D} \boldsymbol{u} x˙=Ax+Bu,y=Cx+Du

系统的极点(Poles),就是矩阵 A \boldsymbol{A} A特征值(Eigenvalues)。也就是满足:

det ⁡ ( s I − A ) = 0 \det(s\boldsymbol{I} - \boldsymbol{A}) = 0 det(sIA)=0

的那些 s s s 值。

极点意味着什么? 把极点看成一个物理系统“天生的频率特性”。比如一个弹簧-质量系统的震动频率,就由它的质量和刚度决定,对应的就是极点。

  1. 系统的固有动态行为
    极点(Poles)系统的本征频率或固有动态特性,反映了系统在 无输入时(自由响应) 如何自然地演化
    比如你把弹簧系统拨一下,然后放手,它怎么震动、震荡多久、会不会越来越剧烈——全由极点决定。
  2. 稳定性
    • 如果 所有极点实部 < 0:系统 稳定(响应会自然衰减)
    • 如果 存在 极点 实部 > 0:系统不稳定(响应会发散)
    • 如果极点 实部 = 0 且为简单极点边界稳定(比如无阻尼震荡)。
    • 极点有 虚部系统振荡
  3. 响应速度
    • 极点的 实部越小(越负),响应越慢
    • 实部越大(更靠近右侧),系统响应越快或越激烈

一个二阶系统(弹簧-质量):
x ¨ + 2 ζ ω n x ˙ + ω n 2 x = u ( t ) \ddot{x} + 2\zeta\omega_n \dot{x} + \omega_n^2 x = u(t) x¨+2ζωnx˙+ωn2x=u(t)
转换为状态空间后,它的极点为:
s = − ζ ω n ± j ω n 1 − ζ 2 s = -\zeta \omega_n \pm j \omega_n \sqrt{1 - \zeta^2} s=ζωn±jωn1ζ2
在 s 平面上(复平面),这个复数对极点的实部和虚部决定:

  • 衰减速度(实部)
  • 震荡频率(虚部)

示例

一个简单的例子

参考:State Space Representations of Linear Physical Systems

考虑一个由单个四阶微分方程表示的四阶系统,该方程具有输入 x x x 和输出 z z z


我们可以定义 4 个新变量, q 1 q_1 q1 q 4 q_4 q4


但是,


现在,我们可以将 4 阶微分方程重写为 4 个一阶方程,


这可以简洁地写成状态空间格式,


其中,

从系统图(电气系统)开发状态空间模型

要为电气系统开发状态空间系统,尝试选择 电容器的电压和电感器的电流作为状态变量。回忆一下,

如果我们可以写出电感两端的电压方程,那么当我们除以电感时,它就变成了状态方程。

即,如果我们有一个关于 e inductor e_\text{inductor} einductor 的方程,除以 L L L,它就变成了关于 d i inductor d t \frac{\text{d}i_\text{inductor}}{\text dt} dtdiinductor 的方程,这是我们状态变量之一。

同样,如果我们能写出通过电容的电流方程,并除以电容,它就变成了关于 e capacitor e_\text{capacitor} ecapacitor 的状态方程。

这最好通过一个例子来说明。推导出下图所示系统的状态空间模型。输入是 i a i_a ia ,输出是 e 2 e_2 e2

该系统有三个储能元件,因此我们预计有三个状态方程。尝试选择 i 1 i_1 i1 i 2 i_2 i2 e 1 e_1 e1 作为状态变量。现在我们 想要它们的导数方程

  • 电感 L 2 L_2 L2 两端的电压是 e 1 e_1 e1 (这是我们状态变量之一)。

    因此,我们第一个状态变量方程是,
  • 如果我们将电流求和到标记为 n 1 n_1 n1 的节点,我们得到,

    该方程包含我们的输入( i a i_a ia)、两个状态变量( i L 2 i_{L2} iL2 i L 1 i_{L1} iL1)以及通过电容器的电流。因此,我们可以得到第二个状态方程,
  • 我们的第三个,也是最后一个状态方程,是通过将 L 1 L_1 L1 上的电压(即 e 2 e_2 e2 )用我们的其他状态变量表示得到的,
  • 我们还需要一个输出方程:

因此,我们的状态空间表示变为

这种方法并不总是容易得到一组状态方程。在某些情况下,开发传递函数模型并将其转换为状态空间模型更容易。

二阶质量-弹簧系统(second-order mass-spring system)

一个经典的质量-弹簧(可带阻尼)系统的运动微分方程为:

m   q ¨ ( t ) + d   q ˙ ( t ) + k   q ( t ) = u ( t ) m\,\ddot{q}(t) + d\,\dot{q}(t) + k\,q(t) = u(t) mq¨(t)+dq˙(t)+kq(t)=u(t)

  • q ( t ) q(t) q(t):质量块的位置(系统的输出)
  • m m m:质量块的质量
  • d d d:阻尼系数(无阻尼时 d = 0 d = 0 d=0
  • k k k:弹簧刚度
  • u ( t ) u(t) u(t):系统所受的外力(输入)

我们定义状态变量为:

x = [ x 1 x 2 ] = [ q q ˙ ] \boldsymbol{x} = \begin{bmatrix} x_1 \\ x_2\\ \end{bmatrix} = \begin{bmatrix} q \\ \dot{q}\\ \end{bmatrix} x=[x1x2]=[qq˙]

根据微分方程可得:

x ˙ 1 = x 2 , x ˙ 2 = − k m x 1 − d m x 2 + 1 m u \dot{x}_1 = x_2,\quad\\[10pt] \dot{x}_2 = -\frac{k}{m}x_1 - \frac{d}{m}x_2 + \frac{1}{m}u x˙1=x2,x˙2=mkx1mdx2+m1u

用矩阵形式写出状态空间模型为:

x ˙ = [ 0 1 − k m − d m ] ⏟ A   x + [ 0 1 m ] ⏟ B   u y = [ 1 0 ] ⏟ C   x + [ 0 ] ⏟ D   u \dot{\boldsymbol{x}} = \underbrace{\begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & -\frac{d}{m}\\ \end{bmatrix}}_{\boldsymbol{A}}\,\boldsymbol{x} + \underbrace{\begin{bmatrix} 0 \\ \frac{1}{m}\\ \end{bmatrix}}_{\boldsymbol{B}}\,u \\[10pt] y = \underbrace{\begin{bmatrix} 1 & 0\\ \end{bmatrix}}_{\boldsymbol{C}}\,\boldsymbol{x} + \underbrace{\begin{bmatrix} 0 \end{bmatrix}}_{\boldsymbol{D}}\,u x˙=A [0mk1md]x+B [0m1]uy=C [10]x+D [0]u

  • A \boldsymbol{A} A:描述状态之间的动态关系
  • B \boldsymbol{B} B:描述输入对状态变化的影响
  • C \boldsymbol{C} C:提取输出(我们这里只观测位移)
  • D = 0 \boldsymbol{D} = 0 D=0:表示没有输入直接传递到输出(无“直通项”)

无阻尼时( d = 0 d = 0 d=0,矩阵 A \boldsymbol{A} A 简化为:

A = [ 0 1 − k m 0 ] \boldsymbol{A} = \begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & 0\\ \end{bmatrix} A=[0mk10]

深入理解系统的极点-零点

  • 极点 告诉你系统自己会怎么“颤动/衰减”——对应齐次响应,是所有模态的汇总;
  • 零点 告诉你在哪些频率上“输入到输出的链路被一脚踹断”,系统对应哪些激励会被完全消除;
  • 极点–零点相同 说明 “内部有这个模态,但外部完全看不到”,正是降阶的目标

极点(Pole)的物理/数学含义

  1. 自然模态(Natural Modes)
    • 极点 s i s_i si 就是矩阵 A A A 的特征值每一个特征值对应一个“模态”——系统在无外力(  ⁣ u = 0 \!u=0 u=0)下,会 沿着对应的特征向量方向,以 e s i t e^{s_i t} esit 的形式 自然衰减或振荡
    • 如果 s i s_i si 是实数,模式是单指数衰减 e s i t e^{s_i t} esit;如果是复数对 α ± j β \alpha\pm j\beta α±jβ,是衰减振荡 e α t cos ⁡ ( β t ) e^{\alpha t}\cos(\beta t) eαtcos(βt)
  2. 时域含义
    • ℜ ( s i ) < 0 \Re(s_i)<0 (si)<0 → 模式衰减,且 ∣ ℜ ( s i ) ∣ |\Re(s_i)| ∣ℜ(si) 越大,衰得越快,系统稳定;
    • ℜ ( s i ) > 0 \Re(s_i)>0 (si)>0模式发散,系统不稳定
    • ℜ ( s i ) = 0 \Re(s_i)=0 (si)=0 ℑ ( s i ) ≠ 0 \Im(s_i)\neq0 (si)=0 → 永久振荡(边界稳定)。
  3. 频域含义
    • 在频率响应里,离开输入的频率越靠近某个复数极点的虚部,系统对该频率的放大倍数越大(共振/峰值)。

零点(Zero)的物理/数学含义

  1. 传输零点(Transmission Zeros)
    • 零点 z j z_j zj 是这样一类频率:即使在输入端 U ( s ) U(s) U(s) 做脉冲/正弦激励,输出 Y ( s ) Y(s) Y(s) 却在此 s = z j s=z_j s=zj 频率处被完全“抑制”(幅度为零)
    • 数学上,零点是使得传递函数矩阵 G ( s ) G(s) G(s) 在某个方向上降秩(或者 分子多项式为零)的 s s s 值。
  2. 输入–输出影响
    • 有零点的模式往往是那种“输入作用到这个模态上,会被系统某种耦合机制立刻在输出端给抵消掉”,即 在传递路径上打了一个“死结”

SISO 情况下的“零点”

对单输入单输出(SISO)系统,传递函数是一个分式

G ( s ) = N ( s ) D ( s ) . G(s) = \frac{N(s)}{D(s)}. G(s)=D(s)N(s).

  • 极点:使 D ( s ) = 0 D(s)=0 D(s)=0 的那些 s s s 值,它们决定 自然模态
  • 零点:使 N ( s ) = 0 N(s)=0 N(s)=0 的那些 s s s 值,它们决定 在对应频率上,输出恰好被“抑制”到零

物理直观:当你给系统一个频率为 s = j ω s=j\omega s= 的正弦输入时,输出幅度 ∣ G ( j ω ) ∣ \bigl|G(j\omega)\bigr| G() 在零点频率 ω 0 \omega_0 ω0 处正好为零——就像你对着一个带有“峰谷滤波器”说话,那些“谷”频率会被完全消除,什么声音都传不下去。

从传递函数出发(SISO) 推导,
Y ( s ) = G ( s )   U ( s ) , U ( s ) ≠ 0. Y(s) = G(s)\,U(s),\quad U(s)\neq 0. Y(s)=G(s)U(s),U(s)=0.
G ( s 0 ) = 0 G(s_0)=0 G(s0)=0,则对任意非零 U ( s ) U(s) U(s),都有 Y ( s 0 ) = 0 Y(s_0)=0 Y(s0)=0。这就是“在 s 0 s_0 s0 频率上零增益”。

MIMO/状态空间下的“传输零点”

对于多输入多输出或者状态空间模型,不能简单地看“分子多项式”。我们用下面的 Rosenbrock 阶跃矩阵 定义零点:

[ s I − A − B C D ] ⏟ P ( s ) \underbrace{\begin{bmatrix} sI - A & -B\\[6pt] C & D \end{bmatrix}}_{\displaystyle P(s)} P(s) [sIACBD]

传输零点 s = z s=z s=z 的定义:存在非零向量 ( x 0 u 0 ) ≠ 0 \begin{pmatrix}x_0\\u_0\end{pmatrix}\neq0 (x0u0)=0,使得
P ( z )   ( x 0 u 0 ) = [ z I − A − B C D ] ( x 0 u 0 ) = 0. P(z)\,\begin{pmatrix} x_0 \\[3pt] u_0 \end{pmatrix} = \begin{bmatrix} zI - A & -B\\[3pt] C & D \end{bmatrix} \begin{pmatrix}x_0\\u_0\end{pmatrix} = 0. P(z)(x0u0)=[zIACBD](x0u0)=0.
这意味着:
{ ( z I − A )   x 0    =    B   u 0 , C   x 0 + D   u 0    =    0. \begin{cases} (zI - A)\,x_0 \;=\; B\,u_0,\\[3pt] C\,x_0 + D\,u_0 \;=\; 0. \end{cases} {(zIA)x0=Bu0,Cx0+Du0=0.

  • 第一行    ( z I − A ) x 0 = B   u 0 \;(zI - A)x_0 = B\,u_0 (zIA)x0=Bu0:表明当输入用复杂指数 u ( t ) = u 0 e z t u(t)=u_0 e^{z t} u(t)=u0ezt 激励时,可在内部激发出一个状态模态 x ( t ) = x 0 e z t x(t)=x_0 e^{z t} x(t)=x0ezt
  • 第二行    C x 0 + D u 0 = 0 \;C x_0 + D u_0 = 0 Cx0+Du0=0:表明对于这对 ( x 0 , u 0 ) (x_0,u_0) (x0,u0),输出 y ( t ) = C x ( t ) + D u ( t ) y(t)=C x(t)+D u(t) y(t)=Cx(t)+Du(t) 恰好恒为 0,输入虽作用于状态,但从状态流到输出的通路被“死结”拦住

因此, P ( s ) P(s) P(s) 的行列式
det ⁡ P ( s ) = 0 \det P(s) = 0 detP(s)=0
所求解的 s s s 就是一组传输零点。

从状态空间刻画(MIMO) 推导,令初始条件为 0,令输入 u ( t ) = u 0 e s t u(t)=u_0 e^{s t} u(t)=u0est。假设存在解
x ( t )    =    x 0 e s t , y ( t ) = C x ( t ) + D u ( t ) . x(t) \;=\; x_0 e^{s t},\quad y(t) = Cx(t)+Du(t). x(t)=x0est,y(t)=Cx(t)+Du(t).
带入微分方程:
s   x 0 e s t = A   x 0 e s t + B   u 0 e s t    ⟹    ( s I − A ) x 0 = B   u 0 . s\,x_0 e^{st} = A\,x_0 e^{st} + B\,u_0 e^{st} \;\Longrightarrow\; (sI - A)x_0 = B\,u_0. sx0est=Ax0est+Bu0est(sIA)x0=Bu0.
同时,
y ( t ) = ( C x 0 + D u 0 ) e s t . y(t)=\bigl(Cx_0 + D u_0\bigr)e^{st}. y(t)=(Cx0+Du0)est.
若想让 y ( t ) ≡ 0 y(t)\equiv0 y(t)0,必须有
C x 0 + D u 0 = 0. C x_0 + D u_0 = 0. Cx0+Du0=0.
非零解 ( x 0 , u 0 ) (x_0,u_0) (x0,u0) 存在的条件,正是 det ⁡ P ( s ) = 0 \det P(s)=0 detP(s)=0,即传输零点

举个例子帮助理解

下面举个具体的例子帮助理解,参考 How Transfer Function Zeros Affect Transient Response


这些脉冲响应来自拥有 相同极点 的系统,那么它们究竟如何产生如此巨大的差异呢?

首先注意到,它们的 传递函数分子(零点位置)是不同的。我们把这两组传递函数 泛化为二阶系统

左侧系统在 s = − z s=-z s=z 处有一个零点;右侧系统没有任何零点。

问题就转化为:一个零点的存在,为什么会导致如此不同的响应?

首先看 无零点系统 H 1 ( s ) H_1(s) H1(s) 的脉冲响应

其中极点为复共轭对(阻尼比 ζ < 1 \zeta<1 ζ<1),可分解并做部分分式展开,取逆拉普拉斯变换后得到时域脉冲响应并简化,


注意, y 1 ( t ) y_1(t) y1(t) 是指数和正弦波的乘积,指数项(极点的实部) 控制衰减或增长速度,正弦项(极点的虚部) 控制振荡频率

现在已经确定了无零脉冲响应,现在来确定 − z -z z 处有一个零点系统 H 2 ( s ) H_2(s) H2(s) 的脉冲响应

同样对分母分解、做部分分式展开,并取逆拉普拉斯变换,化简后可得,

事情开始变得有趣了,

  • 看一下第二项 z − a b \frac{z-a}{b} bza随着 z z z 接近 a a a,频率响应 y 1 y_1 y1 会缩放,特别是当 z z z 接近于 p p p 的实部 a a a 时,极点位于 − p -p p ( z − a ) (z-a) (za) 项会消失,该 模态与 y 1 y_1 y1 完全抵消
  • 第二个观察结果是,零点不会影响衰减或放大的速率(即不会改变极点的实部)

根据这些观察结果,现在要问的是,零点位置如何影响脉冲响应

叠加各个零点位置的响应

将系统极点固定在复平面上的某对位置,依次 把零点从右半平面向左半平面移动,可以同时绘制其脉冲响应,发现:

  • 零点最初位于右半平面,这会导致瞬态响应下降;
  • 当零点向左移动到左平面时,瞬态响应上升。


对两种情况的脉冲响应进行扩展,现在将单零点响应扩展为三个项,看看单零点响应中有什么?

原来的无零点项按 z z z 缩放。再看当零点从右半平面移动到左半平面时,这个项如何影响瞬态响应,

零点的位置 s s s 与第三项中包含的 z z z 值相反,

  • 零点位于左半平面时 z > 0 z > 0 z>0,因此添加了第三项,这就是为什么 z z z 变大时,瞬态响应会增大(过冲)


因此在单零瞬态响应中,我们有按 z z z 缩放的原始无零响应,以及 由于零点而产生的额外瞬态响应

这个额外的瞬态响应是什么?

对每个项进行拉普拉斯变换,我们发现瞬态是无零响应的导数,

注意传递函数 H 2 H_2 H2 现在表示为 H 1 H_1 H1 的函数,零点的影响是缩放无零响应,并将无零响应的导数添加到它,

现在考虑 H 2 H_2 H2 行为如何随着 z z z 的变化而变化,

  • 随着 z z z 的幅度变大,缩放的无零响应占主导地位;
  • 随着 z z z 幅度变小,无零响应的导数占主导地位。

在这两个极端之间, H 2 H_2 H2 H 1 H_1 H1 及其导数的组合。因此,在 y 2 y_2 y2 的时间域中,我们查看的是 y 1 y_1 y1 及其导数的线性组合。

H 2 H_2 H2 中分解出 H 1 H_1 H1,并注意 H 2 H_2 H2 是两个传递函数的乘积。

  • H 1 H_1 H1 将输入 U U U 转换为其输出 Y 1 Y_1 Y1
  • s = − z s=-z s=z 处的零点将 H 1 H_1 H1 的输出通过传递函数 s + z s+z s+z 进行变换,乘积是 H 2 H_2 H2,即单零传递函数,或者说是 H 2 H_2 H2 将输入信号 U U U 转换为输出 Y 2 Y_2 Y2

现在,我们便从 频域(零点的代数作用)时域(叠加导数项),再回到频域看其物理意义。这使我们能够从多个角度理解零点对瞬态响应的影响。

在该例子中, H 1 H_1 H1 H 2 H_2 H2 都是二阶系统,也就是说 每个系统都有两个极点(成对的复共轭极点):

  • 极点:一对位于 s = − a ± b j s=-a\pm bj s=a±bj(阻尼比 ζ < 1 \zeta<1 ζ<1)的复共轭极点。
  • 零点
    • 对于 H 1 ( s ) H_1(s) H1(s)(无零点)来说,分子是常数,只有这两个极点。
    • 对于 H 2 ( s ) H_2(s) H2(s)(带零点)来说,在分子又多了一个 ( s + z ) (s+z) (s+z) 项,因此它有一个额外的零点(位于 s = − z s=-z s=z),但极点依旧只有那对 − a ± b j -a\pm bj a±bj

当把这个 零点移动到恰好和极点位置重合 z = a z=a z=a无虚部 b b b 时),那对极点就会 在传递函数里和零点 完全相消——那一对极点在输入—输出上被抹去,系统表现上像降到了一阶或零阶,从输入到输出的路径上,这个模态“隐形”了。

在滤波器设计中,常用 严格本性(无零点)函数来抑制噪声;而一个零点的加入,实际上是对无零点响应的“反向滤波”,它通过缩放和叠加导数可能在信号上引入更大的瞬态振荡。

总结

  • 极点 决定系统的 衰减/振荡 速率(自然模态);
  • 零点 不会改变衰减率,而是表现为振荡,会通过 叠加导数项 引入 超调/下冲,实际上正在转换无零响应;
  • 当零点与极点相撞时,会 相互抵消,该模态在输入–输出通路上完全“隐形”
  • 正是零点位置的变化,让两个拥有 相同极点 的系统展现出 截然不同 的瞬态响应。

深入理解控制系统稳定性、可控性和可观性

参考:What is stability of control system ?

稳定性

稳定性指的是 闭环系统维持平衡状态或在被施加扰动或变化后返回该状态的能力

  1. 内部(Lyapunov)稳定性:考察系统状态 x ˙ = A x \dot x = Ax x˙=Ax没有外部输入下的固有动态

    • 渐近稳定(Asymptotic Stability):所有解 x ( t ) x(t) x(t) t → ∞ t\to\infty t 都趋于零。
    • 李雅普诺夫稳定(Lyapunov Stability):初始扰动足够小,状态始终保持在小范围内(但不一定趋于零)。
    • 不稳定:存在初始条件使得 ∥ x ( t ) ∥ → ∞ \|x(t)\|\to\infty x(t)
  2. BIBO 稳定性(Bounded-Input Bounded-Output Stability)
    有界输入和有界输出(BIBO)稳定性,通常称为输入输出稳定性,评估系统对有界(有限)输入的响应。如果一个系统    y ( t ) = G ( s ) u ( t ) \;y(t)=G(s)u(t) y(t)=G(s)u(t) 对于任何有界输入信号,其输出都保持有界,即不会无限增长,则该系统被认为是 BIBO 稳定的。

注意内部稳定 ⇔ BIBO 稳定(当系统是可控可观且最小实现时) ,否则两者可能不同。

  • 若系统包含不可控不可观测极点,这些极点在输入—输出传递函数中会被零点抵消,看不到(它们属于内部模态,不影响 BIBO)。
  • 最小实现(Minimal Realization)通过 minreal 或可控/可观子空间剔除这些“被抵消”的模态,使得传递函数中的极点集真正反映系统动态。

一个稳定的控制系统是指 输出保持有界,并且 即使受到外部力量或输入变化的干扰,也不会表现出不受控制的或振荡的行为

相位裕度增益裕度 等数学表达式用于量化 负反馈系统的稳定性

系统稳定性分类

判定条件(连续时域)根据
渐近稳定 ℜ ( s i ) < 0 \Re(s_i)<0 (si)<0极点位置
边界稳定 ℜ ( s i ) ≤ 0 \Re(s_i)\le0 (si)0,无重根极点位置
不稳定存在 ℜ ( s i ) > 0 \Re(s_i)>0 (si)>0 或重根极点位置
零点任意位置不影响稳定性
非最小相位零点在右半平 ℜ ( z j ) > 0 \Re(z_j)>0 (zj)>0零点位置影响瞬态
  • 极点 决定指数衰减/增长振荡,是真正的稳定性指标
  • 零点 决定输入—输出路径上的局部“抑制”相位/幅度畸变,但不改变是否收敛。

举个简单的例子,假设我们有一个传递函数 G ( s ) G(s) G(s)
G ( s ) = 1 s + a G(s)=\frac{1}{s+a} G(s)=s+a1
所以 − a -a a 就是我们的极点,为了理解这在时间域中意味着什么,取系统的逆拉普拉斯变换 L − 1 \mathcal L^{-1} L1,得到系统:
f ( t ) = e − a t u ( t ) f(t) = e^{-at}u(t) f(t)=eatu(t)
假设初始条件为 0,可以将 u u u 设为 1,我们的系统就变成了
f ( t ) = e − a t f(t) = e^{-at} f(t)=eat

稳定系统

如果 系统的传递函数的极点位于 s 平面的左侧,则认为 系统是稳定的。但随着极点靠近原点,系统的稳定性会降低。

脉冲响应显示,

  • 位于实轴上的极点具有过阻尼响应
  • 对于第二和第三象限的极点,瞬态响应是正弦波但指数衰减(欠阻尼响应)

接近原点的极被称为主导极(Dominant Poles)。因此,如果一个稳定系统的极位于 p 1 p_1 p1 p 2 p_2 p2(如上图),则 p 1 p_1 p1 被认为是该系统的主导极。 p 2 p_2 p2 被认为是该系统的非主导极(Non-dominant Poles)。在瞬态响应中,由于 p 2 p_2 p2 的贡献衰减得更快,响应将由主导极 p 1 p_1 p1 主导

边界稳定系统

一个系统如果其 极点位于虚轴上且不重复,则称为边界稳定

时域脉冲响应显示在恒定幅度下的振荡。初始条件(输入幅度)定义了幅度

不稳定系统

如果一个系统 在 s 平面的右侧有任何极点,则称为 不稳定。即使只有一个极点在右侧,也可能使系统不稳定。

输出随时间增长,振荡幅度随时间指数增长。在 实际系统中,振幅会上升到电源电压

可控性与可观性

理解并检验可控性与可观测性,是状态空间设计的第一步,也是保证反馈与估计器有效性的核心前提。

  1. 状态反馈:只有在可控系统上,才能通过反馈 u = − K x u=-Kx u=Kx 任意配置闭环极点。
  2. 状态估计:只有在可观测系统上,才能用观测器(如卡尔曼滤波、Luenberger 观测器)准确重建状态。
  3. 最小实现:不可控或不可观测的模态在输入—输出上无贡献,可用 minreal 删除,得到最小等价系统

考虑 LTI 系统的状态空间描述:

{ x ˙ ( t ) = A   x ( t ) + B   u ( t ) , y ( t ) = C   x ( t ) + D   u ( t ) , \begin{cases} \dot x(t) = A\,x(t) + B\,u(t),\\[5pt] y(t) = C\,x(t) + D\,u(t),\\ \end{cases} x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t)+Du(t),

其中 x ∈ R n x\in\mathbb R^n xRn u ∈ R m u\in\mathbb R^m uRm y ∈ R p y\in\mathbb R^p yRp

可控性(Controllability)

可控性(Controllability) 系统从 任意初始状态 x ( 0 ) = x 0 x(0)=x_0 x(0)=x0有限时间 T T T 内,通过恰当选择输入 u ( t ) u(t) u(t),能否将状态驱动到 任意目标状态 x ( T ) = x f x(T)=x_f x(T)=xf

判断条件:Kalman 可控性矩阵法则
C = [   B      A B      A 2 B    …    A n − 1 B   ] ∈ R n × n m \mathcal{C} = [\,B\;\;A B\;\;A^2B\;\dots\;A^{n-1}B\,]\in\mathbb R^{n\times nm} C=[BABA2BAn1B]Rn×nm

判据:若 r a n k ( C ) = n \mathrm{rank}(\mathcal C)=n rank(C)=n,则系统完全可控;否则存在不可控子空间。

定义可控 Gramian 判据:

W c ( T ) = ∫ 0 T e A τ B B ⊤ e A ⊤ τ   d τ . W_c(T) = \int_{0}^{T} e^{A\tau} B B^\top e^{A^\top \tau}\,\mathrm d\tau. Wc(T)=0TeAτBBeAτdτ.

  • 若对某 T > 0 T>0 T>0 W c ( T ) W_c(T) Wc(T) 是可逆的(正定矩阵),系统可控。
  • 物理/几何意义
    • 每个列向量 A k B A^k B AkB 表示“从输入作用一次、两次……到第 k + 1 k+1 k+1 次,对状态的影响方向”。
    • C \mathcal C C 的列空间若能张满整个状态空间,就说明输入足够“驱动”到任何方向。

可观测性(Observability)

可观性(Observability) 能否通过观测输出 y ( t ) y(t) y(t) 在有限时间内唯一重构系统初始状态 x ( 0 ) = x 0 x(0)=x_0 x(0)=x0

判断条件:Kalman 可观测性矩阵法则
O = [ C C   A C   A 2 ⋮ C   A n − 1 ] ∈ R p n × n \mathcal{O} = \begin{bmatrix} C\\ C\,A\\ C\,A^2\\ \vdots\\ C\,A^{n-1}\\ \end{bmatrix} \in\mathbb R^{pn\times n} O= CCACA2CAn1 Rpn×n

判据:若 r a n k ( O ) = n \mathrm{rank}(\mathcal O)=n rank(O)=n,则系统完全可观测;否则存在不可观测子空间。

定义可观测 Gramian 判据:
W o ( T ) = ∫ 0 T e A ⊤ τ C ⊤ C e A   τ   d τ . W_o(T) = \int_{0}^{T} e^{A^\top \tau} C^\top C e^{A\,\tau}\,\mathrm d\tau. Wo(T)=0TeAτCCeAτdτ.

  • 若对某 T > 0 T>0 T>0 W o ( T ) W_o(T) Wo(T) 是可逆的,系统可观测。
  • 物理/几何意义
    • O \mathcal O O 的行向量 C A k C A^k CAk 表示“状态向第 k k k 次演化后,对输出的映射”。
    • 若行空间张满状态空间,说明每个状态分量都能在某个时刻对输出产生可区分影响。

简单示例:二阶质量-弹簧系统

状态空间(无阻尼):

A = [ 0 1 − k / m 0 ] , B = [ 0 1 / m ] , C = [ 1 0 ] A=\begin{bmatrix}0&1\\-k/m&0\end{bmatrix},\quad\\[10pt] B=\begin{bmatrix}0\\1/m\end{bmatrix},\quad\\[10pt] C=\begin{bmatrix}1&0\end{bmatrix} A=[0k/m10],B=[01/m],C=[10]

  • 可控性 C = [ B ,    A B ] = [   [ 0 ; 1 / m ] ,    [ − k / m ; 0 ]   ] \mathcal C=[B,\;AB] = \bigl[\,[0;1/m],\;[-k/m;0]\,\bigr] C=[B,AB]=[[0;1/m],[k/m;0]] rank ⁡ ( C ) = 2 \operatorname{rank}(\mathcal C)=2 rank(C)=2,完全可控。
  • 可观测性 O = [ C C A ] = [   [ 1 , 0 ] ;    [ 0 , 1 ] ] \mathcal O=\begin{bmatrix}C\\CA\end{bmatrix} = \bigl[\,[1,0];\;[0,1]\bigr] O=[CCA]=[[1,0];[0,1]] rank ⁡ ( O ) = 2 \operatorname{rank}(\mathcal O)=2 rank(O)=2,完全可观测。

在数学上, ( A , B ) (A,B) (A,B) 的可控性等价于 ( A ⊤ , C ⊤ ) (A^\top, C^\top) (A,C) 的可观测性。这意味着可控性的分析工具几乎可以“对转”直接用于可观测性。

深入理解可控性和可观性

参考:A Conceptual Approach to Controllability and Observability

首先要了解你想要控制的系统,然后

  • 可以使用 执行器(Actuators) 来控制系统,执行器会传递控制力和能量。
  • 可以使用 传感器(Sensors) 测量系统的特定状态,例如电压、位置和温度,这些传感器会产生系统的输出。
  • 控制器(Controller) 设计归结为 如何利用传感器数据以及参考信号或某种命令来生成正确的执行器命令


如果你的系统没有合适的执行器来影响系统的正确部分,或者你没有安装合适的传感器来控制,那么你的控制器设计注定会失败。如果没有这两个参数,你就无法充分地影响系统。系统将无法控制,无法被观察。

因此,可控性和可观测性是系统如何与执行器和传感器协同工作的条件,它与特定的控制技术(如 PID 或极点配置)无关。

  • 可控性意味着存在控制信号,允许系统在有限的时间内 达到任何状态,这也称为可达性。所以可控性并不意味着必须维持该状态,而是指 可以达到该状态
  • 可观测性意味着所有状态(关心的状态)都可以从系统的输出中得知。也就是说,我们可以通过使用适当的传感器选择和传感器位置来了解系统中每个状态的值。

现在,让我们抛开约束的现实问题,看看汽车的高维状态空间。我们关心的状态是 XY 位置、XY 速度、偏航角和偏航率。

假设,

  • 传感器有车速表来测量速度,我们的眼睛用来评估其他所有东西。
  • 执行器是方向盘、油门和刹车踏板。

通过这套设置,我们 能够以给定的速度和偏航角控制汽车到特定位置

先来看一下 可观测性的概念,并从一个极端的例子开始,以帮助说明其重要性。

我们闭上眼睛,移除所有传感器信息,这会导致 C C C 矩阵归零。

基本上,我们无法知道任何状态。你不知道自己在哪里,行驶速度是多少,或者你面对的方向。在这个例子中,我们假设你感觉不到加速,也听不到引擎的声音,或者其他任何东西。

现在你仍然可以控制汽车,这并没有受到影响,因为你可以转动方向盘并踩下踏板。但是,由于 缺乏可观测性,你不能安全地驾驶它到达目的地。然而,根据我们之前的定义,汽车是可控的

相反的情况,我们移除方向盘和踏板,这会导致 B B B 矩阵归零。

想象一下,如果你撞到一大片冰面,汽车打滑,会发生什么情况。 并且不受转向或制动的影响,即使您无法影响汽车的状态,该 系统仍然是可观察的,因为您可以查看速度表,并且可以通过观察外部来跟踪您的位置。

可以看到拥有一个可观察和控制的系统是多么必要,它们协同工作,即使是 部分无法控制或无法观察 的系统仍然很麻烦,想象一下尝试在方向盘拆下的情况下驾驶汽车,即使仍然可以通过踏板改变其速度,但无法改变偏航状态。

如果任何一个临界状态无法控制,那么整个系统都被视为无法控制,如果任何一个临界状态无法观察,那么整个系统都被视为无法观察。

状态空间方程 —— 极点配置

极点配置(Pole Placement)完全状态反馈 是控制理论中的一种设计方法,用于通过 适当选择控制器增益,使得闭环系统的极点(也即系统特征方程的根)位于预定的位置,从而达到所需的动态性能和稳定性要求。极点的位置直接影响系统的响应特性,如稳定性、速度、阻尼等。

现在,极点配置本身在工业中并没有得到广泛的应用。你可能会更常使用其他方法,例如 LQR 或 H ∞ H_\infin H,但是极点配置值得花一些时间,因为它能让你更好地理解 使用状态空间方程进行反馈控制 的一般方法,这是进入其他方法的垫脚石。

参考:What is Pole Placement (Full State Feedback)

从 PID 到极点配置

首先,我们有一个具有输入 u u u 和输出 y y y 的装置(Plant),目标是 开发一个反馈控制系统将输出驱动到某个期望值

你可能熟悉的一种方法是将输出与参考信号进行比较以获得控制误差,然后你可以开发一个控制器,该控制器使用该误差来生成输入到装置的信号,目标是将误差驱动到零

但对于极点配置,我们将以不同的方式解决这个问题,我们 将反馈 状态向量中每个状态变量的值,而不是反馈输出 y y y

现在,我们要 假设我们知道每个状态的值,即使它不一定是输出 y y y 的一部分,然后 取状态向量并将其乘以一个由一堆不同增益值组成的矩阵 K K K,并 将结果从缩放的参考信号 K r K_r Kr 中减去,并将此结果直接输入到我们的装置中作为输入。

这里没有像在反馈结构框图中标记为控制器的块,图中绿圈整个部分都是控制器。

通过极点配置 可以 计算适当的增益矩阵 K K K确保系统稳定性,参考上的缩放项用于确保稳态误差性能在可接受范围内。

状态方程中的第一项 A x Ax Ax 捕获线性系统的动态特性,第二项是 系统如何响应输入,但系统中的能量如何存储和移动是由 A x Ax Ax 项捕获的,所以你可能会认为在控制器设计方面, A A A 矩阵有一些特殊之处。

是的!反馈控制器都必须 修改 A A A 矩阵才能改变系统的动态特性,尤其是在稳定性方面


A A A 矩阵的特征值是系统的极点,极点的位置决定了线性系统的稳定性,这是极点配置的关键,通过移动极点或闭环 A A A 矩阵的特征值来生成所需的闭环稳定性

举例说明特征值如何影响系统

看下面这个例子, A A A 矩阵转换为使用另一组状态变量来描述系统的矩阵,此转换是使用 变换矩阵 完成的,其列是 A A A 矩阵的特征向量


转换后我们得到的是一个修改后的 A A A 矩阵,由对角线上的复数特征值和其他所有位置的零组成,现在 绿色箭头所指的这两个模型代表完全相同的系统,它们具有 相同的特征值和相同的动态,只是第二个模型使用一组彼此独立变化的状态变量来描述

A A A 矩阵以对角形式写出时,很容易看出我们剩下的是 一组一阶微分方程,其中每个状态的导数仅受该状态的影响,而不受其他任何影响

很酷的是,像这样的微分方程的解的形式是 Z = C e λ T Z=C e^{\lambda T} Z=CeλT,其中 λ \lambda λ 是给定状态变量的特征值。

看几个不同的 λ \lambda λ 值,这样你就可以直观地看到 能量是如何根据复平面内特征值的位置变化的

  • 如果 λ \lambda λ 是负实数,那么这是一个稳定的特征值,因为解是 e λ t → 0 e^{\lambda t} \rightarrow 0 eλt0 t → ∞ t \rightarrow \infin t,任何初始能量都会随时间消散;
  • 但如果 λ \lambda λ 是正数,那么它是不稳定的,因为能量只会随时间增长;
  • 如果有一对 虚数特征值,那么这种模式下的能量将 振荡,因为 e e e 的虚数部分会产生 s i n sin sin c o s cos cos
  • 而实数和虚数的任意组合都会产生振荡和指数能量耗散。

极点配置如何移动

现在我们可以陈述我们要解决的问题了。如果对象的特征值位于复平面中不理想的位置,那么我们可以使用极点配置将它们移动到其他地方

当然,如果它们位于右半平面,这是理想的,因为它们不稳定,但非理想位置也可能意味着存在振荡,你想消除它,或者只是加速或减缓特定模式下的能量耗散。

我们现在可以讨论极点配置如何移动特征值。记住我们一开始画的控制器结构,这会导致输入 u = r K r − K x u= rK_r - K x u=rKrKx,其中 r K r r K_r rKr 是缩放参考值, K x Kx Kx 是状态向量,将其反馈 x x x 乘以增益矩阵 K K K

现在,如果 将此控制输入代入状态方程,将 闭合循环,得到以下状态方程,注意 A A A − B k -Bk Bk 都作用于状态向量,因此我们可以将它们组合起来,得到一个修改后的 A A A 矩阵,这是 闭环 A A A 矩阵,可以通过 选择适当的 K K K 来移动特征值

二阶系统示例

对于简单的系统,这很容易手动完成,让我们尝试一个具有单个输入的二阶系统的示例,

我们可以通过将 det ( A − λ I ) = 0 \text{det} (A - \lambda I) = 0 det(AλI)=0 来找到特征值,然后求解 λ \lambda λ,它们分别位于 − 2 -2 2 + 1 +1 +1 处,由于存在正的实特征值,其中一个模式将爆炸到无穷大,因此整个 系统是不稳定的

让我们使用极点配置设计一个反馈控制器,通过将不稳定的极点移动到左半平面来稳定该系统,我们的闭环 A A A 矩阵是 A − B k A-Bk ABk,增益矩阵 K K K 1 × 2 1\times 2 1×2,因为有一个输出 两个状态。

对得到的新 A − B k A-Bk ABk 进行,像以前一样求解特征值,得到这个特征方程,它是两个增益值的函数。现在 假设我们 希望闭环极点位于 − 1 -1 1 − 2 -2 2,这样特征方程为 λ 2 + 3 λ + 2 = 0 \lambda^2 + 3 \lambda + 2 = 0 λ2+3λ+2=0,此时很容易找到合适的 k 1 k_1 k1 k 2 k_2 k2,使这两个方程相等,我们 只需将系数设置为相等并求解,得到 k 1 = 2 k_1 = 2 k1=2 k 2 = 1 k_2 = 1 k2=1

就是这样,如果我们将这两个增益放在该系统的状态反馈路径中,它将在特征值位于 − 1 -1 1 − 2 -2 2 时保持稳定。

对于具有两个以上状态的系统,所涉及的数学运算开始变得难以理解。想法一样的,只是求解行列式变得不切实际,但在 MATLAB 几乎只用一个命令。

% Define state matrices
A = [0 1; 2 -1];
B = [1; 0];
C = [1 0];
D = 0;

% Create state space object
sys = ss(A, B, C, D);

% Check open loop eigenvalues
E = eig(A);

% Desired closed loop eigenvalues
P = [-2 -1];

% Solve for K using pole placement
K = place(A, B, P);

% Check for closed loop eigenvalues
Acl = A - B*K;
Ecl = eig(Acl);

% Create closed loop system
syscl = ss(Acl, B, C, D);

% Check step response
step(sys)

% Create colsed loop system
syscl = ss(Acl, B, C, D);

% Check step response
step(syscl)

可以预见的是,开环系统的阶跃是不稳定的

闭环系统的阶跃响应看起来要好得多,但它并不完美,而不是像我们预期的那样上升到 1,预计稳态输出仅为 0.5,这最终就是缩放项发挥作用的地方。

对状态方程进行设计前、仿真前的体检

控制系统设计前的状态方程体检

在拿到一组状态空间形式
{ x ˙ = A x + B u , y = C x + D u , \begin{cases} \dot x = A x + B u,\\ y = C x + D u, \end{cases} {x˙=Ax+Bu,y=Cx+Du,
之后,通常要从以下几方面去“体检”系统,才能保证后续的设计(比如状态反馈、观测器设计、滤波器设计等)有的放矢。

检查项判据结果不合格时的“补救”
稳定性 ℜ ( λ i ) < 0 \Re(\lambda_i)<0 (λi)<0设计状态反馈或动态补偿
可控性 rank ⁡ ( C ) = n \operatorname{rank}(\mathcal{C})=n rank(C)=n增执行器/降阶建模
可观测性 rank ⁡ ( O ) = n \operatorname{rank}(\mathcal{O})=n rank(O)=n增传感器/降阶建模
最小实现同时可控可观Kalman 分解去冗余
适定性 & 直通 D D D 是否零,因果性检查置零 D D D 或加滤波

1. 稳定性(Stability)

  • 检查项:矩阵 A A A 的特征值 { λ i } \{\lambda_i\} {λi}
  • 原理:连续系统内在稳定要求所有 ℜ ( λ i ) < 0 \Re(\lambda_i)<0 (λi)<0(特征值的实部);若有 ℜ ( λ i ) > 0 \Re(\lambda_i)>0 (λi)>0,系统本身就不稳定;若有 ℜ ( λ i ) = 0 \Re(\lambda_i)=0 (λi)=0,则可能边界稳定或振荡。
  • 调整方法
    • 若不稳定,设计 状态反馈 u = − K x u=-Kx u=Kx,使闭环矩阵 A c l = A − B K A_{\rm cl}=A-BK Acl=ABK 的特征值都在左半平面(如极点配置法、LQR)。
    • 若有边界或振荡模态,可适当加阻尼或配合动态补偿器。

举个例子,假设状态方程 A A A所有特征值 λ i λ_i λi 的实部均为负(最小约 − 0.000002326 −0.000002326 0.000002326,最大约 − 9.41 × 10 4 −9.41×10^4 9.41×104),则 系统渐近稳定
原理:连续时间系统 x ˙ = A x \dot x = A x x˙=Ax 的解为 x ( t ) = e A t   x ( 0 )   x(t) = e^{A t}\,x(0)\, x(t)=eAtx(0),当且仅当 ∀ i , ℜ ( λ i ) < 0 ∀i, \Re( λ_i)<0 i,(λi)<0 时, e λ i t → 0 e^{\lambda_i t}\to0 eλit0 t → + ∞ t→+∞ t+),系统才是渐近稳定的。
公式:特征方程 det ⁡ ( λ I − A ) = 0 \det(\lambda I - A)=0 det(λIA)=0,解得 λ i λ_i λi,即系统模式指数。

2. 可控性(Controllability)

  • 检查项:构造可控性矩阵
    C = [ B ,    A B ,    A 2 B ,    … ,    A n − 1 B ] , \mathcal{C} = \bigl[B,\;A B,\;A^2 B,\;\dots,\;A^{n-1}B\bigr], C=[B,AB,A2B,,An1B],
    并检验 rank ⁡ ( C ) = ? n \operatorname{rank}(\mathcal{C})\stackrel{?}{=}n rank(C)=?n
  • 原理:若 rank ⁡ ( C ) < n \operatorname{rank}(\mathcal{C})<n rank(C)<n,说明存在状态方向无法被输入 u u u 影响,这些“死角”无法通过任何输入到达。
  • 调整方法
    • 硬件层面:增设或重新布置执行器,使得原来不可控的模态得以激励。
    • 软件层面:对可控子空间进行降阶,只保留可控部分作为等价最小模型;对不可控部分只能当做扰动或噪声处理。

3. 可观测性(Observability)

  • 检查项:构造可观测性矩阵
    O = [ C C A C A 2 ⋮ C A n − 1 ] , \mathcal{O} = \begin{bmatrix}C\\C A\\C A^2\\\vdots\\C A^{n-1}\end{bmatrix}, O= CCACA2CAn1 ,
    并检验 rank ⁡ ( O ) = ? n \operatorname{rank}(\mathcal{O})\stackrel{?}{=}n rank(O)=?n
  • 原理:若 rank ⁡ ( O ) < n \operatorname{rank}(\mathcal{O})<n rank(O)<n,说明存在状态方向即便变化也无法从输出 y y y 中分辨。
  • 调整方法
    • 硬件层面:增设或调整传感器,使得原来不可观测的模态被量测。
    • 软件层面:对可观测子空间降阶处理,只建模可观测部分;对不可观测部分用作未建模动态。

4. 最小化/约简(Minimal Realization)

  • 检查项:可控性 & 可观测性都要满秩;若有任一不满秩,系统不是 最小实现
  • 原理:最小实现去掉了所有不可控或不可观测模态,模型既不冗余,也不过度简化;可提高数值稳健性和计算效率。
  • 调整方法
    • 利用 Kalman 分解Gramian 分解 将系统分解成可控可观、可控不可观、不可控可观、不可控不可观四个子系统,丢弃不可控或不可观子系统。

5. 传递函数的适定性(Properness & Direct‐Feedthrough)

  • 检查项:矩阵 D D D 是否为零?系统阶次是否满足因果性(严格、准)。
  • 原理
    • D ≠ 0 D\ne0 D=0,意味着有“直通”通道,输入 u u u 能瞬时影响输出 y y y
    • 对于物理系统,过高的直通增益可能导致测量噪声放大或数值震荡。
  • 调整方法
    • 若不希望直通,令 D = 0 D=0 D=0
    • 或在设计控制器时,加低通滤波器以抑制高频噪声。

掌握了以上原则和方法,就能针对任意给定的 A , B , C , D A,B,C,D A,B,C,D 矩阵,先“体检”再“对症下药”,确保后续控制器/观测器设计既可行又高效。

仿真数值求解前的状态方程体检

体检项检查方法结果不合格时的调整思路
刚性(Stiffness) τ max ⁡ τ min ⁡ \frac{\tau_{\max}}{\tau_{\min}} τminτmax、谱半径比 max ⁡ ∣ ℜ ( λ ) ∣ min ⁡ ∣ ℜ ( λ ) ∣ \displaystyle\frac{\max | \Re(\lambda) |}{\min | \Re(\lambda)|} min∣ℜ(λ)max∣ℜ(λ)、Gramian 条件数、仿真对比多时尺度分离、模型降阶、数值正则化、选用刚性求解器

刚性(Stiffness)

除了检查极点实部以外,“刚性”(stiffness)确实也是“体检”时一个值得关注的性质,尤其在以下两种情境中尤为重要:

  1. 数值仿真与离散化:当系统的模态(极点)时间常数相差悬殊时,常规的数值积分器(如显式欧拉、Runge–Kutta)往往要求 极小的步长 才能保持稳定 和精度。
  2. 分层/多时尺度控制:当某些模态非常快(刚性模态),而另一些模态又很慢,如果不加以分离、近似或降阶,设计出的控制器可能在工程上既难以实现又成本高昂。

刚性的定义与量化:

  1. 最大时间常数和最小非零时间常数之比

    • 定义每个模态的时间常数为:
      τ i = 1 ∣ ℜ ( λ i ) ∣ , 仅当  ℜ ( λ i ) ≠ 0 \tau_i = \frac{1}{|\Re(\lambda_i)|}, \quad \text{仅当 } \Re(\lambda_i) \neq 0 τi=∣ℜ(λi)1,仅当 (λi)=0
    • 刚性比(stiffness ratio): R = τ max ⁡ τ min ⁡ R = \frac{\tau_{\max}}{\tau_{\min}} R=τminτmax
    • R ≫ 1 R \gg 1 R1(比如 R > 10 3 R > 10^3 R>103 10 5 10^5 105),系统被认为是刚性的(stiff system),可能导致显式数值解发散。

    举个例子,假设系统所有特征值的实部,最小约 − 0.000002326 −0.000002326 0.000002326,最大约 − 9.41 × 10 4 −9.41×10^4 9.41×104,则根据极大/极小特征值可以估算,

    • 最快模态约 τ min ⁡ ≈ 1 / 9.4 × 10 4 ≈ 10 − 5   s \tau_{\min}\approx 1/9.4\times10^4 \approx10^{-5}\,\mathrm{s} τmin1/9.4×104105s
    • 最慢模态约 τ max ⁡ ≈ 1 / 2.3 × 10 − 6 ≈ 4 × 10 5   s \tau_{\max}\approx 1/2.3\times10^{-6}\approx4\times10^5\,\mathrm{s} τmax1/2.3×1064×105s(数天量级)
    • 刚性指标(stiffness ratio) τ max ⁡ τ min ⁡ ∼ 10 10 ≫ 1 , \frac{\tau_{\max}}{\tau_{\min}}\sim10^{10}\gg1, τminτmax10101,
      则为高刚性系统。
  2. 频谱比(Spectral Radius Ratio)
    Stiffness Ratio = max ⁡ i ∣ ℜ ( λ i ) ∣ min ⁡ j   :   ℜ ( λ j ) ≠ 0 ∣ ℜ ( λ j ) ∣ . \text{Stiffness Ratio} = \frac{\max_i \bigl|\Re(\lambda_i)\bigr|}{\min_{j\,:\,\Re(\lambda_j)\neq0}\bigl|\Re(\lambda_j)\bigr|}. Stiffness Ratio=minj:(λj)=0 (λj) maxi (λi) .
    如果该比值远大于 10 2 10^2 102 10 3 10^3 103,就常被认为是刚性系统。

  3. 条件数(Condition Number)

    • 状态转移矩阵 e A t e^{At} eAt 的条件数高,或者

    • Gramian 矩阵(可控 Gramian W c W_c Wc、可观测 Gramian W o W_o Wo)的条件数高,
      都会导致求解或降阶过程中的数值不稳定。

      检查 Gramian 条件数
      W c = ∫ 0 T e A t B B T e A T t   d t , W o = ∫ 0 T e A T t C T C e A t   d t . W_c = \int_0^{T} e^{A t} B B^T e^{A^T t} \,\mathrm{d}t,\quad W_o = \int_0^{T} e^{A^T t} C^T C e^{A t} \,\mathrm{d}t. Wc=0TeAtBBTeATtdt,Wo=0TeATtCTCeAtdt.
      计算 κ ( W c ) \kappa(W_c) κ(Wc) κ ( W o ) \kappa(W_o) κ(Wo)

  4. 也可以从 数值仿真试验 结果来减压系统的刚性,

    • 用不同步长的显式积分器(如 RK4)模拟系统,观察是否存在“爆渐”(数值不稳定)或“过度阻尼”现象;
    • 与刚性专用积分器(如隐式 Euler、Gear’s 方法)对比,若差异显著,则说明刚性问题突出。

刚性不合格时的“对症下药”

  1. 多时尺度分离
    • 奇异摄动方法将系统拆成 “快子系统”“慢子系统”,分别降阶或分别设计控制器,再用级联/分层方式集成。
    • 模态分解:对极值很大的快模态进行定性建模或近似“瞬态消隐”,只对慢模态进行精细设计。
  2. 模型降阶
    • 平衡截断(Balanced Truncation):在保留主要能量传递通道的同时,丢弃对输入—输出映射影响非常小的高频刚性模态
    • 投影法:如 Krylov 子空间法,针对高频模态快速收敛的特点进行降阶。
  3. 数值处理
    • 状态变量尺度变换:对状态量做合理归一化或尺度分割,减少数值差异;
    • 雅可比正则化:在离散化或仿真时,对 A A A 做小量级正则化(如 A + ϵ I A+\epsilon I A+ϵI)以改善条件数。
  4. 控制器/观测器设计
    • 含刚性补偿器:设计时考虑频率加权,避免对极高频模态过度放大;
    • 选用刚性友好算法:如 LQR 中的 Ricatti 方程求解时,选用隐式方法及高精度求解器。
  5. 仿真工具选型
    • 明确使用“刚性求解器”(Stiff Solver),例如 MATLAB 中的 ode15sode23s,或仿真平台相应的隐式积分器。

将“刚性”纳入常规体检清单后,可在早期就识别潜在的数值稳定性与实现难题,从而在模型简化、控制器设计、仿真与在线实现等各个环节做到有的放矢。这样,既能保证系统性能,又可提高仿真与在线运行的鲁棒性和效率。

MATLAB 绘制极点-零点分布

%% 1. 构造系统并算极零点
sys = ss(A,B,C,D);  % 自定义的 A,B,C,D

%% 2. 绘制极点-零点图
figure
hold on;
pzmap(sys)
grid on

%% 3. 叠加阻尼比和自然频率线
sgrid

底图并不是一个“正方形”网格,而是典型的 s‑平面(damping–frequency)网格:

  • 圆弧(constant natural frequency)
  • 放射状直线(constant damping ratio)
  • x 轴、y 轴刻度不等距(非等高宽)

MATLAB 自带了 Control Toolbox 里的 sgrid 函数,专门用来往 s‑平面上画这两种网格。

从极点(Poles)和传输零点(Transmission Zeros)结果来看,主要暴露出以下几个问题:

  1. 高阶系统中存在大量“重合”极点–零点
    • 除了一个 幅值约为 − 9.4108 × 10 4 -9.4108\times10^4 9.4108×104 的快模,几乎所有极点都落在 − 0.0 … ∼ − 1.0 -0.0…\sim-1.0 0.01.0 之间,而且传输零点列表里几乎“抄”了一份极点(同样的值、同样的幅度)——这说明在这些频率点上存在 成对的极点–零点
    • 这正是典型的 “极点–零点相互抵消”(pole–zero cancellation):系统在这些模态上既有极点又有相同位置的零点,结果在开环或闭环响应中,这部分动态都会被“抹去”,对外表现得像不存在这些模态
  2. 模型“非最简”/不必要的高阶
    • 当从某个更复杂的结构(比如带有很多中间变量、滤波器、观测器等)直接写出 A , B , C , D A,B,C,D A,B,C,D 后,很容易出现 非最小实现(non‑minimal realization):可控但不可观测或可观测但不可控的状态分量,被保留在模型里,却对输入—输出响应没有实质贡献
  3. 大量频率接近零的模式
    • 大量接近 0 (如 − 0.0000 -0.0000 0.0000)的极点和零点,往往是数值精度或模型降阶/近似造成的“伪模式”——物理上并不存在这么多几乎为零的慢模态(“伪零点”)。

当看到:

极点: -1.8423, -1.7864, …  
零点: -1.8423, -1.7864,

一一对应,说明:

  1. 这些 极点–零点对 在传递函数里互相抵消,对输入–输出动态没有贡献
  2. 换言之,这部分“内部模态”虽然存在于 A A A 矩阵,却对任何外部输入 u u u 并不会在输出 y y y 上留下痕迹——要么是 不可观测,要么是 不可控,或者两者兼有。

真正对系统输入—输出贡献 的是 被保留下来的“非抵消”极点/零点,对这些模式之外的部分,可以 用 minreal 或降阶方法清理掉。这样既让模型更简洁,也能避免后续控制器设计中出现不必要的麻烦。

  1. 最小化模型: 在 MATLAB 中对状态空间对象调用

    sys_min = minreal(sys);
    

    它会自动做可控性/可观测性分析,去掉所有“成对抵消”的极点–零点,得到最简的等价系统。

  2. 数值条件检查:

    • ctrb(A,B)obsv(A,C) 查看可控/可观矩阵秩,确认是否存在不可控或不可观测子空间。
    rank(ctrb(A,B)), rank(obsv(A,C))
    
    • 若秩小于状态维数,说明有冗余(不可控或不可观),需要先做可控/可观子空间分解。
  3. 平衡降阶(可选) 如果想在 保留主要特性下更进一步降阶,可用

    sys_bal = balreal(sys);
    sys_red = modred(sys_bal, indices_to_remove);
    

    或直接

    sys_red = balred(sys, r);
    

    用 Hankel 奇异值自动去掉贡献最小的模式。

  4. 避免数值爆炸/精度问题:

    • 如果原始 A A A 矩阵元素跨度极大(如首个极点是 10 4 10^4 104,其余都是 10 − 1 10^{-1} 101 量级),建议先 缩放状态 或做 坐标变换 让系统矩阵条件数改善。
    • 保持同一数量级,有助于数值算法(计算特征值、零点)更稳健。

注意,上图是 SISO 系统。
如果是 MIMO 系统,可能会像下图一样只有极点而没有零点,

  • 如果是 MIMO 且模型里 D = 0 D = 0 D=0,说明所有的传递函数都是“真分式”(分子次数 < 分母次数)。此时系统的相对阶数 > 0,多出来的阶数会在 ∞ \infin 处“贡献”零点,不会在复平面上出现任何有限零点
  • MATLAB 的 tzero(sys) 针对的是 整块 MIMO 系统的 传输(invariant)零点,它反映了系统输入到输出的“整体”可控可观特性。即便每个通道(SISO)都有零点,只要它们在组合里没有成为“传输零点”,tzero 也不会把它们列出来

✨ 为什么 可以 模型降阶?

含义
极点 = 模态存在系统里有这个模式
零点 = 模态抵消输入–输出通路上不允许这个模式传播
极点 = 零点模态虽然存在,但完全被抑制,对输入输出行为 没有任何影响
极点 = 零点 为什么无用?不影响响应、不影响控制器、不改变系统性能,删掉更高效
  • 移除“丝毫不影响”输入–输出性能的那些模式,让模型变得更小、更高效
  • 降阶后剩下的极点,才是真正决定系统对外(关心的那个输入–输出通道)动力学性能的 “有效模式”;
  • 去除无用模式还有利于控制器设计(避免去补偿那些根本检测不到或控制不到的动态)。

举一个简单的例子,考虑一个 SISO 系统,传递函数为:

G ( s ) = ( s + 3 ) ( s + 1 ) ( s + 3 ) G(s) = \frac{(s + 3)}{(s + 1)(s + 3)} G(s)=(s+1)(s+3)(s+3)
可以对上面的式子做约分,化简后得到的系统:
G 简化 ( s ) = 1 s + 1 G_{\text{简化}}(s) = \frac{1}{s + 1} G简化(s)=s+11

  • 原本极点有 { − 1 , − 3 } \{-1,-3\} {1,3},零点有 { − 3 } \{-3\} {3}
  • 因为 − 3 -3 3 频率上“零点–极点成对抵消”,最终等价于一个一阶系统 1 / ( s + 1 ) 1/(s+1) 1/(s+1)

为什么这个抵消的模态没用了?

可以从时域和频域两个角度来理解。

  • 时域角度理解
    系统本来包含模态:

    • e − t e^{-t} et(极点 -1)→ 留下来了
    • e − 3 t e^{-3t} e3t(极点 -3)→ 被零点抵消了!

    也就是说:系统本来有两个动态过程在响应输入时发生,一个快的、一个慢的。但是由于 s = − 3 s = -3 s=3 同时是零点,所以该模态虽然 系统内部存在,但它 对输出的影响被输入路径“抵消”掉了
    想象成:这个模态存在于系统内部,但输入 u u u 不会激发它,输出 y y y 也看不到它

  • 频域角度理解
    来看系统的频率响应:

    • 零点是让系统在某个频率 完全无响应
    • 极点是系统在某个频率 响应增强(或延迟)。

    极点和零点相同 时,这两种影响完全抵消掉了。在 s = − 3 s = -3 s=3 附近,刚想要“振起来”,系统立刻把那一部分“消掉”,结果就是没有那部分响应。

在控制系统里,我们只关心 输入 u u u输出 y y y 的通路。如果某个内部状态的动态,输入激不起它,输出也看不到它,那它对控制器设计来说是 完全无用的:无需观测它、无需控制它、不影响性能、不影响稳定性。

所以在建模或简化模型时,完全可以删掉它

根据极点-零点分布进行降阶

给定一个高阶状态空间模型,希望得到一个低阶模型,其动态行为在指定频段内近似等效(比如频率响应、阶跃响应等尽量接近)。

构造降阶模型常用的方法

  • 平衡截断(Balanced Truncation):虽然不是直接基于极点零点,但会保留能量最大的状态,对极点-零点分布也有一定反映。数值稳定。

  • 极点匹配法(Pole-Matching Reduction):直接构造一个只有保留下来的极点的低阶模型。可通过匹配保留的极点和零点重新构建系统。

  • 共振频率保留法(尤其适合振动系统):保留低频处主导极点,丢弃高频极点。

误差评估

  • 频率响应对比(Bode 图)
  • 时域响应对比(阶跃响应、冲激响应)
  • H ∞ \mathcal{H}_\infty H H 2 \mathcal{H}_2 H2范数误差评估

本文标签: 可控性零点可观稳定性函数