1. 从“单输入单输出”到“多输入多输出”为什么我们需要传递函数矩阵如果你是从经典控制理论比如拉普拉斯变换、伯德图、奈奎斯特判据一路学过来的那么“传递函数”这个概念对你来说应该像呼吸一样自然。一个输入一个输出一个函数关系清晰明了。我们分析系统的稳定性、动态性能都围绕着这个单一的G(s)展开。但现实世界的系统尤其是现代工程系统很少有这么“单纯”的。想想看一台无人机你需要同时控制它的俯仰、横滚、偏航和高度。你的操纵指令四个电机的转速是多个输入它的姿态和位置是多个输出。一个化工反应釜你需要控制温度、压力、pH值和进料流量。加热功率、冷却水阀、酸碱泵和进料泵是你的输入传感器读数是你的输出。一台机械臂每个关节的电机扭矩是输入末端执行器的位置和姿态是输出。这些系统都有一个共同点多输入多输出英文简称MIMO。在经典控制里我们可能会尝试为每一个输出单独设计一个控制器去调节某一个输入这叫做“单回路控制”。但问题很快就来了——这些回路之间是耦合的。你调节无人机第一个电机想让它抬头结果它可能不仅抬头还开始往左偏航。你加大反应釜的加热功率想升温压力可能也跟着飙升。这种耦合意味着你不能把系统简单地拆成几个独立的单变量系统来处理。它们是一个整体牵一发而动全身。这时沿用单个传递函数的思路就捉襟见肘了。我们需要一个数学工具能够同时描述所有输入对所有输出的影响关系。这个工具就是传递函数矩阵。你可以把它想象成一张关系网或者一个表格。假设系统有p个输入q个输出那么这个传递函数矩阵G(s)就是一个q x p的矩阵。矩阵中的每一个元素G_ij(s)都是一个传递函数它单独地描述了第j个输入对第i个输出的影响而此时其他所有输入都为零这是线性系统叠加原理的前提。所以当我们从“现控理论”的视角重新审视系统时传递函数矩阵就是我们连接系统外部描述输入输出关系和内部状态描述状态空间方程的一座关键桥梁也是处理MIMO系统分析与设计问题的起点。它保留了传递函数直观的频率域特性又将系统的多变量特性用矩阵这一强大的数学语言清晰地表达了出来。2. 传递函数矩阵的核心定义与数学表达理解了为什么需要它我们来看看它具体是什么。我们从最通用的线性时不变系统的状态空间描述出发这是现代控制理论的基础模型状态方程 ẋ(t) A x(t) B u(t) 输出方程 y(t) C x(t) D u(t)其中x(t)是n维状态向量。u(t)是p维输入向量。y(t)是q维输出向量。A, B, C, D是维数匹配的常数矩阵。我们对上述方程两边进行拉普拉斯变换并假设初始状态x(0) 0这是为了聚焦于输入输出的传递特性就像经典控制里一样。得到sX(s) A X(s) B U(s) Y(s) C X(s) D U(s)由第一个式子解出X(s)(sI - A) X(s) B U(s) X(s) (sI - A)^{-1} B U(s)这里I是n x n的单位矩阵(sI - A)^{-1}就是预解矩阵它在系统分析中至关重要。将X(s)代入输出方程Y(s) C (sI - A)^{-1} B U(s) D U(s) [C (sI - A)^{-1} B D] U(s)于是我们得到了输入U(s)和输出Y(s)在复频域的关系Y(s) G(s) U(s)其中传递函数矩阵G(s)被定义为G(s) C (sI - A)^{-1} B D这是一个q x p的矩阵。它的每一个元素G_ij(s)都是一个关于复变量s的有理分式函数分子分母都是s的多项式。G_ij(s)表示的是当只有第j个输入U_j(s)作用时所产生第i个输出Y_i(s)的传递函数。一个关键的理解点(sI - A)^{-1}包含了系统的所有动态模态信息极点而矩阵C和B则决定了这些模态如何被输入激发以及如何被输出观测。矩阵D代表了输入到输出的直接馈通在很多物理系统中如电路、机械D常常是零矩阵因为输入通常不会瞬间无延迟地影响输出。注意这里假设了(sI - A)是可逆的这要求s不是系统矩阵A的特征值。系统矩阵A的特征值正是传递函数矩阵G(s)的极点分母为零的点它们决定了系统的稳定性。3. 如何计算与解读传递函数矩阵一个详细的计算实例定义看起来有点抽象我们通过一个具体的、简化的双输入双输出系统来亲手算一遍感受一下这个过程。考虑一个耦合的弹簧-质量块系统或者一个简单的电路网络我们可以抽象出如下状态空间模型为了计算方便数字是假设的设系统方程为A [ -2 1 0; 1 -3 1; 0 1 -1 ]; B [ 1 0; 0 1; 1 0 ]; C [ 1 0 0; 0 1 0 ]; D [ 0 0; 0 0 ];这里状态x是3维输入u是2维输出y是2维。D矩阵为零表示没有直接馈通。我们的目标是求出G(s) C (sI - A)^{-1} B。第一步构造 (sI - A)sI - A [ s2 -1 0; -1 s3 -1; 0 -1 s1 ];第二步求 (sI - A) 的逆矩阵求一个3x3矩阵的逆我们可以用伴随矩阵法(sI - A)^{-1} adj(sI - A) / det(sI - A)。 先求行列式Δ(s) det(sI - A)Δ(s) (s2)*[(s3)(s1) - (-1)(-1)] - (-1)*[(-1)(s1) - (0)(-1)] 0*... (s2)[(s3)(s1) - 1] 1*[-(s1)] (s2)(s^2 4s 3 - 1) - (s1) (s2)(s^2 4s 2) - s - 1 s^3 4s^2 2s 2s^2 8s 4 - s - 1 s^3 6s^2 9s 3所以系统的特征多项式是s^3 6s^2 9s 3它的根就是系统的极点。再求伴随矩阵adj(sI - A)。这需要计算所有代数余子式过程略这是基本功练习。假设我们算得adj(sI - A) [ (s3)(s1)-1 (s1) 1; (s1) (s2)(s1) (s2); 1 (s2) (s2)(s3)-1 ] [ s^24s2 s1 1; s1 s^23s2 s2; 1 s2 s^25s5 ]因此(sI - A)^{-1} (1 / Δ(s)) * [ s^24s2, s1, 1; s1, s^23s2, s2; 1, s2, s^25s5 ]第三步计算 C (sI - A)^{-1} B这是一个连续的矩阵乘法。C是 2x3逆矩阵是 3x3B是 3x2。所以结果G(s)是 2x2。 我们先计算中间结果T(s) (sI - A)^{-1} B。B矩阵有两列我们分别计算。 令B1 [1; 0; 1](B的第一列)B2 [0; 1; 0](B的第二列)。计算T1(s) (sI - A)^{-1} * B1T1(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [1; (1/Δ(s)) * [ (s^24s2)*1 (s1)*0 1*1; s1, s^23s2, s2; 0; (s1)*1 (s^23s2)*0 (s2)*1; 1, s2, s^25s5 ] 1] 1*1 (s2)*0 (s^25s5)*1 ] (1/Δ(s)) * [ s^24s2 1; s1 s2; 1 s^25s5 ] (1/Δ(s)) * [ s^24s3; 2s3; s^25s6 ]同理计算T2(s) (sI - A)^{-1} * B2T2(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [0; (1/Δ(s)) * [ (s^24s2)*0 (s1)*1 1*0; s1, s^23s2, s2; 1; (s1)*0 (s^23s2)*1 (s2)*0; 1, s2, s^25s5 ] 0] 1*0 (s2)*1 (s^25s5)*0 ] (1/Δ(s)) * [ s1; s^23s2; s2 ]所以T(s) [T1(s), T2(s)] (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2; s^25s6, s2 ]现在左乘C矩阵。C矩阵只取前两行因为输出y只对应前两个状态所以G(s) C * T(s)相当于取T(s)的前两行。G(s) (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2 ]其中Δ(s) s^3 6s^2 9s 3。解读G(s) 这个 2x2 的传递函数矩阵G(s)告诉我们G_11(s) (s^24s3) / (s^36s^29s3)输入u1对输出y1的传递函数。G_12(s) (s1) / (s^36s^29s3)输入u2对输出y1的传递函数。G_21(s) (2s3) / (s^36s^29s3)输入u1对输出y2的传递函数。G_22(s) (s^23s2) / (s^36s^29s3)输入u2对输出y2的传递函数。关键观察共同分母所有四个传递函数都有相同的分母多项式Δ(s)。这个多项式正是矩阵A的特征多项式。这意味着整个MIMO系统的极点决定稳定性和基本动态是由矩阵A决定的是所有输入输出通道所共享的。这是MIMO系统与多个独立SISO系统本质不同的地方。耦合体现非对角线元素G_12(s)和G_21(s)不为零。这说明输入u1会影响输出y2输入u2也会影响输出y1。系统是耦合的。相对阶G_11和G_22的分子是二阶分母是三阶相对阶为1。G_12和G_21的分子是一阶相对阶为2。不同通道的动态特性可能不同。实操心得对于维数高于3的系统手工计算逆矩阵会非常繁琐且容易出错。在实际工程或学习中我们强烈依赖数学工具。在MATLAB/Octave中计算传递函数矩阵极其简单。假设你已经定义了A, B, C, D矩阵只需一行命令G ss(A, B, C, D); G_tf tf(G);。ss创建状态空间模型tf将其转换为传递函数矩阵形式。对于我们的例子在MATLAB中验证上述计算是很好的练习。手工推导的目的在于深刻理解其数学本质和耦合关系的来源而不是用于解决大规模问题。4. 传递函数矩阵的性质与在系统分析中的核心作用得到了传递函数矩阵我们能用它做什么它不仅仅是状态空间模型的一种等价表示更是我们分析MIMO系统特性的强大工具。4.1 极点与零点MIMO系统的扩展定义对于SISO系统极点就是传递函数分母多项式的根零点就是分子多项式的根。对于MIMO系统由于G(s)是一个矩阵极点和零点的定义需要扩展。极点如前所述传递函数矩阵G(s)的所有元素的分母多项式在约去公因子后的根就是系统的极点。更本质地说系统的极点就是系统矩阵A的特征值。这一定义从状态空间模型出发更为根本。所有输入输出通道共享同一组极点这决定了系统的稳定性所有极点均具有负实部则渐近稳定和基本响应模态如振荡频率、衰减速度。传输零点这是一个MIMO系统特有的、非常重要的概念。它不是单个元素分子为零的点。传输零点的定义是存在一个非零的复频率z和一个非零的输入向量U_0使得在零初始状态下系统的输出Y(s)恒为零。即G(z) U_0 0。这意味着在频率z处存在某种特定的输入信号组合其效果被系统内部完全“抵消”了无法在输出端被观测到。计算对于D0的系统传输零点z是使得复合矩阵P(s) [sI-A, -B; C, 0]降秩的s值。在MATLAB中可以使用tzero或zero函数直接计算。物理意义传输零点反映了系统输入输出之间的阻塞特性。例如在飞机控制中某些特定的舵面偏转组合可能无法引起飞机姿态的变化在某个频率下这个频率就是传输零点。传输零点会影响系统的可控制性、可观测性以及控制性能的极限例如非最小相位系统的右半平面零点会限制控制带宽。4.2 稳定性分析基于极点的稳定性判据在MIMO系统中依然直接适用线性时不变系统渐近稳定的充要条件是其传递函数矩阵的所有极点即矩阵A的所有特征值都具有负实部。通过传递函数矩阵G(s)我们可以计算其特征多项式即行列式det(sI-A)或G(s)各元素分母的最小公倍式然后应用劳斯判据、赫尔维茨判据或直接求根来判断稳定性。在MATLAB中使用pole(G)或eig(A)直接获取极点。注意对于MIMO系统不能仅仅通过观察G(s)某个对角元比如G_11(s)的稳定性来判断整个系统的稳定性。因为非对角元的耦合可能引入不稳定的隐藏模态与状态空间中的能控性/能观性相关。必须检查系统矩阵A的全部特征值。4.3 频域分析从伯德图到奇异值图这是传递函数矩阵威力巨大的地方。在SISO系统中我们绘制单个传递函数G(jω)的伯德图幅频和相频特性。在MIMO系统中G(jω)是一个复数矩阵。我们如何绘制它的频率响应逐个元素法可以为G(s)的每一个元素G_ij(jω)绘制伯德图。这能让我们看清每一对输入输出通道在频域的特性。这对于理解耦合的强度比如G_21相对于G_22的幅值大小很有帮助。在MATLAB中bode(G)命令会自动为所有通道生成伯德图。奇异值分析法更强大的工具对于MIMO系统更本质的频域分析工具是奇异值。对于每个频率点ω我们计算复数矩阵G(jω)的奇异值。奇异值总是非负的实数。假设G是q x p矩阵那么它在每个频率点有min(p, q)个奇异值记为σ_1(ω) ≥ σ_2(ω) ≥ ... ≥ σ_min(p,q)(ω) ≥ 0。最大奇异值σ_max(ω)可以理解为系统在频率ω处对所有可能输入方向向量的最大增益。最小奇异值σ_min(ω)可以理解为系统在频率ω处对所有可能输入方向的最小增益。绘制σ_max(ω)和σ_min(ω)随频率变化的曲线就得到了MIMO系统的奇异值伯德图。为什么重要在鲁棒控制和系统性能分析中奇异值提供了关键的洞察。例如σ_min在低频段的大小反映了系统抗干扰和解耦的能力σ_max在高频段的衰减速率反映了系统对噪声和未建模动态的鲁棒性。闭环系统的稳定裕度也可以用开环传递函数矩阵的奇异值来评估。在MATLAB中可以使用sigma(G)命令绘制奇异值图。4.4 能控性与能观性与传递函数矩阵的关联能控性和能观性是状态空间模型的核心概念它们与传递函数矩阵有着深刻的联系。能控性系统是否能在有限时间内通过合适的输入u(t)将状态从任意初始点驱动到原点。状态空间判据是能控性矩阵[B, AB, A^2B, ..., A^{n-1}B]满秩。能观性系统是否能在有限时间内通过输出的观测y(t)唯一地确定初始状态x(0)。状态空间判据是能观性矩阵[C; CA; CA^2; ...; CA^{n-1}]满秩。从传递函数矩阵的角度看如果系统是状态空间能控且能观的那么其传递函数矩阵G(s)将是既约的即没有零极点对消。如果发生了零极点对消则意味着系统不是完全能控和/或完全能观的对消掉的模态在输入输出描述中“消失”了但它们可能隐藏在系统内部影响实际动态甚至稳定性。在计算G(s) C(sI-A)^{-1}B D时如果最终得到的各元素传递函数有公因子被约去这些被约去的因子就对应着不能控或不能观的模态。因此在由传递函数矩阵反推或实现状态空间模型时必须非常小心要确保得到的是最小实现即完全能控且完全能观。5. 从理论到实践在MATLAB/SciPy中操作传递函数矩阵理论学习之后我们必须能在工具中实现它。这里以MATLAB及其开源替代品GNU Octave和Python的SciPy库为例展示基本操作。5.1 MATLAB/Octave 环境1. 定义系统与计算传递函数矩阵% 定义状态空间矩阵 (沿用之前的例子) A [-2, 1, 0; 1, -3, 1; 0, 1, -1]; B [1, 0; 0, 1; 1, 0]; C [1, 0, 0; 0, 1, 0]; D zeros(2,2); % 创建状态空间模型对象 sys_ss ss(A, B, C, D); % 转换为传递函数矩阵形式 sys_tf tf(sys_ss); disp(传递函数矩阵 G(s):); disp(sys_tf); % 输出会显示一个2x2的tf数组每个元素都是tf对象。 % 例如从输入 1 到输出 1: s^2 4 s 3 / s^3 6 s^2 9 s 3 % 从输入 2 到输出 1: s 1 / s^3 6 s^2 9 s 3 ... 等等2. 提取特定通道的传递函数% 提取输入1到输出1的SISO传递函数 G11 sys_tf(1, 1); % 提取输入2到输出2的SISO传递函数 G22 sys_tf(2, 2);3. 计算极点和零点% 计算系统极点 (即A的特征值) poles pole(sys_ss); % 或 pole(sys_tf) disp(系统极点); disp(poles); % 计算系统传输零点 zeros tzero(sys_ss); % 或 zero(sys_tf) disp(系统传输零点); disp(zeros);4. 频域分析绘图% 绘制所有通道的伯德图 (4个子图) figure; bode(sys_ss); % 或 bode(sys_tf) grid on; title(各通道伯德图); % 绘制奇异值图 (MIMO系统频域分析利器) figure; sigma(sys_ss); % 或 sigma(sys_tf) grid on; title(系统奇异值图);5. 时域仿真% 定义时间向量 t 0:0.01:10; % 定义输入信号u1是阶跃u2是正弦波 u [ones(size(t)); sin(2*t)]; % 注意维度是 2 x length(t) % 进行仿真 [y, t_out, x] lsim(sys_ss, u, t); % 绘制输出响应 figure; subplot(2,1,1); plot(t_out, y(:,1)); ylabel(y1); grid on; legend(输出1); subplot(2,1,2); plot(t_out, y(:,2)); ylabel(y2); xlabel(时间 (s)); grid on; legend(输出2);5.2 Python (SciPy Matplotlib) 环境Python在科学计算领域应用广泛控制库虽不如MATLAB专业但基础功能完备。import numpy as np import matplotlib.pyplot as plt from scipy import signal # 1. 定义状态空间矩阵 A np.array([[-2, 1, 0], [1, -3, 1], [0, 1, -1]]) B np.array([[1, 0], [0, 1], [1, 0]]) C np.array([[1, 0, 0], [0, 1, 0]]) D np.array([[0, 0], [0, 0]]) # 2. 创建状态空间系统 sys_ss signal.StateSpace(A, B, C, D) # 3. 转换为传递函数形式 (注意scipy.signal 的 tf 表示是针对SISO的) # 对于MIMO我们需要逐个元素计算或使用专门的库如 control需安装pip install control # 这里演示使用 control 库 (如果已安装) try: import control as ct # 创建控制系统库的状态空间对象 sys_ct ct.ss(A, B, C, D) # 计算传递函数矩阵 (返回一个 TransferFunction 对象列表的列表) sys_tf ct.tf(sys_ct) print(传递函数矩阵:) print(sys_tf) except ImportError: print(control 库未安装。对于MIMO传递函数建议安装 control 库。) # 手动计算其中一个通道作为示例 (例如 G11) # G(s) C(sI-A)^-1 B D, 我们计算第一个元素 # 这里省略手动计算逆矩阵的代码较为复杂。 # 作为替代我们可以直接进行时域/频域仿真。 # 4. 计算极点和零点 (使用 scipy.signal 或 control) # 极点就是 A 的特征值 poles np.linalg.eigvals(A) print(\n系统极点 (A的特征值):) print(poles) # 传输零点计算较复杂control库有提供 try: import control as ct zeros ct.zero(sys_ct) print(\n系统传输零点 (来自 control 库):) print(zeros) except: print(\n传输零点计算需要 control 库。) # 5. 时域仿真 (使用 scipy.signal) t np.arange(0, 10, 0.01) # 输入u1 阶跃 u2 正弦波 u np.vstack([np.ones_like(t), np.sin(2*t)]) # 形状 (2, len(t)) # lsim 需要输入 u 的形状为 (len(t), num_inputs) t_out, y, x signal.lsim(sys_ss, Uu.T, Tt) # 注意 u.T 转置 # 绘制结果 fig, (ax1, ax2) plt.subplots(2, 1, figsize(10, 6)) ax1.plot(t_out, y[:, 0]) ax1.set_ylabel(y1) ax1.grid(True) ax1.legend([输出1]) ax2.plot(t_out, y[:, 1]) ax2.set_ylabel(y2) ax2.set_xlabel(时间 [s]) ax2.grid(True) ax2.legend([输出2]) plt.suptitle(系统时域响应) plt.tight_layout() plt.show() # 6. 频域分析 - 伯德图 (scipy.signal 的 bode 是 SISO 的) # 我们可以循环绘制每个通道 try: import control as ct plt.figure() ct.bode_plot(sys_ct) # control库可以绘制MIMO伯德图 plt.suptitle(各通道伯德图 (control库)) plt.tight_layout() plt.show() except: print(绘制MIMO伯德图需要 control 库。)实操心得与避坑指南工具选择对于严肃的控制系统分析与设计MATLAB及其控制系统工具箱仍然是行业标杆文档齐全、函数丰富。Python的control库功能日益完善对于学习、研究和轻量级应用是完全足够的且免费开源。SciPy.signal更侧重于信号处理对MIMO支持较弱。维度匹配在定义矩阵A, B, C, D时务必反复检查维度。(n x n), (n x p), (q x n), (q x p)。一个常见的错误是B或C矩阵的维度定义反了。仿真输入使用lsim或signal.lsim进行仿真时输入矩阵U的维度容易出错。在MATLAB中U的列数等于输入数行数等于时间点数量。在SciPy的signal.lsim中U的形状是(len(t), num_inputs)。务必查阅文档或打印矩阵形状来确认。最小实现如果你是从物理方程推导或系统辨识得到的传递函数矩阵在转换为状态空间模型时使用minreal函数MATLAB或minreal方法control库来获取最小实现消除不能控或不能观的模态避免后续分析与设计出现问题。6. 传递函数矩阵的局限性与状态空间法的优势尽管传递函数矩阵是分析MIMO系统的有力工具但它并非万能也存在固有的局限性。理解这些局限性正是我们深入学习现代控制理论状态空间法的动力。1. 仅适用于线性时不变系统传递函数建立在拉普拉斯变换的基础上其核心是叠加原理和常系数线性微分方程。对于非线性系统、时变系统传递函数矩阵的定义不再成立。而状态空间方程ẋ f(x, u, t), y g(x, u, t)在形式上可以描述更广泛的系统。2. 丢失内部状态信息只反映输入输出关系这是传递函数方法最根本的局限。G(s)描述的是“黑箱”的外部特性。如果系统不是完全能控和完全能观的那么传递函数矩阵就无法反映系统的全部动态。那些不能控或不能观的模态在G(s)中会被对消掉。但从内部看这些隐藏模态可能是不稳定的会导致实际系统出现问题。状态空间法通过直接研究状态向量x(t)能完整揭示系统的内部行为。3. 难以处理非零初始条件传递函数分析通常假设零初始条件。对于非零初始状态的系统响应传递函数方法处理起来比较麻烦。而状态空间方程结合初始状态x(0)可以非常自然地求解全响应。4. 多变量系统分析与综合的复杂性虽然传递函数矩阵将多变量关系封装了起来但在进行控制器设计时如经典的频域设计法直接处理一个矩阵函数仍然非常复杂。例如如何为MIMO系统设计一个PID控制器每个通道单独设计往往会因为耦合而失败。而状态空间法通过状态反馈u -Kx和观测器设计提供了系统化的多变量控制器设计框架如线性二次型调节器LQR、极点配置等概念上更统一。5. 对系统结构的洞察较弱状态空间表示法通过A, B, C, D矩阵清晰地分离了系统动态(A)、输入影响(B)、输出测量(C)和直接传递(D)。这种结构化的表示更容易与物理模型对应也便于进行能控性、能观性等结构性质的分析。那么为什么我们还要学习传递函数矩阵因为它是连接经典控制与现代控制的桥梁。它保留了频域分析的直观性伯德图、奈奎斯特图为理解系统的频率响应、带宽、鲁棒稳定性提供了图形化工具。许多先进的多变量频域设计方法如H∞鲁棒控制也是建立在传递函数矩阵的基础之上。在实际工程中常常是频域方法基于传递函数矩阵和时域方法基于状态空间结合使用取长补短。因此传递函数矩阵是现代控制理论工具箱中不可或缺的一部分。它让我们能够用熟悉的频域语言去理解和分析复杂的多变量系统同时又提醒我们其边界所在引导我们走向更深刻、更强大的状态空间分析领域。掌握了它你就拿到了理解和处理现实世界中耦合、多变量系统的第一把钥匙。
从SISO到MIMO：传递函数矩阵的核心原理、计算与MATLAB实践
1. 从“单输入单输出”到“多输入多输出”为什么我们需要传递函数矩阵如果你是从经典控制理论比如拉普拉斯变换、伯德图、奈奎斯特判据一路学过来的那么“传递函数”这个概念对你来说应该像呼吸一样自然。一个输入一个输出一个函数关系清晰明了。我们分析系统的稳定性、动态性能都围绕着这个单一的G(s)展开。但现实世界的系统尤其是现代工程系统很少有这么“单纯”的。想想看一台无人机你需要同时控制它的俯仰、横滚、偏航和高度。你的操纵指令四个电机的转速是多个输入它的姿态和位置是多个输出。一个化工反应釜你需要控制温度、压力、pH值和进料流量。加热功率、冷却水阀、酸碱泵和进料泵是你的输入传感器读数是你的输出。一台机械臂每个关节的电机扭矩是输入末端执行器的位置和姿态是输出。这些系统都有一个共同点多输入多输出英文简称MIMO。在经典控制里我们可能会尝试为每一个输出单独设计一个控制器去调节某一个输入这叫做“单回路控制”。但问题很快就来了——这些回路之间是耦合的。你调节无人机第一个电机想让它抬头结果它可能不仅抬头还开始往左偏航。你加大反应釜的加热功率想升温压力可能也跟着飙升。这种耦合意味着你不能把系统简单地拆成几个独立的单变量系统来处理。它们是一个整体牵一发而动全身。这时沿用单个传递函数的思路就捉襟见肘了。我们需要一个数学工具能够同时描述所有输入对所有输出的影响关系。这个工具就是传递函数矩阵。你可以把它想象成一张关系网或者一个表格。假设系统有p个输入q个输出那么这个传递函数矩阵G(s)就是一个q x p的矩阵。矩阵中的每一个元素G_ij(s)都是一个传递函数它单独地描述了第j个输入对第i个输出的影响而此时其他所有输入都为零这是线性系统叠加原理的前提。所以当我们从“现控理论”的视角重新审视系统时传递函数矩阵就是我们连接系统外部描述输入输出关系和内部状态描述状态空间方程的一座关键桥梁也是处理MIMO系统分析与设计问题的起点。它保留了传递函数直观的频率域特性又将系统的多变量特性用矩阵这一强大的数学语言清晰地表达了出来。2. 传递函数矩阵的核心定义与数学表达理解了为什么需要它我们来看看它具体是什么。我们从最通用的线性时不变系统的状态空间描述出发这是现代控制理论的基础模型状态方程 ẋ(t) A x(t) B u(t) 输出方程 y(t) C x(t) D u(t)其中x(t)是n维状态向量。u(t)是p维输入向量。y(t)是q维输出向量。A, B, C, D是维数匹配的常数矩阵。我们对上述方程两边进行拉普拉斯变换并假设初始状态x(0) 0这是为了聚焦于输入输出的传递特性就像经典控制里一样。得到sX(s) A X(s) B U(s) Y(s) C X(s) D U(s)由第一个式子解出X(s)(sI - A) X(s) B U(s) X(s) (sI - A)^{-1} B U(s)这里I是n x n的单位矩阵(sI - A)^{-1}就是预解矩阵它在系统分析中至关重要。将X(s)代入输出方程Y(s) C (sI - A)^{-1} B U(s) D U(s) [C (sI - A)^{-1} B D] U(s)于是我们得到了输入U(s)和输出Y(s)在复频域的关系Y(s) G(s) U(s)其中传递函数矩阵G(s)被定义为G(s) C (sI - A)^{-1} B D这是一个q x p的矩阵。它的每一个元素G_ij(s)都是一个关于复变量s的有理分式函数分子分母都是s的多项式。G_ij(s)表示的是当只有第j个输入U_j(s)作用时所产生第i个输出Y_i(s)的传递函数。一个关键的理解点(sI - A)^{-1}包含了系统的所有动态模态信息极点而矩阵C和B则决定了这些模态如何被输入激发以及如何被输出观测。矩阵D代表了输入到输出的直接馈通在很多物理系统中如电路、机械D常常是零矩阵因为输入通常不会瞬间无延迟地影响输出。注意这里假设了(sI - A)是可逆的这要求s不是系统矩阵A的特征值。系统矩阵A的特征值正是传递函数矩阵G(s)的极点分母为零的点它们决定了系统的稳定性。3. 如何计算与解读传递函数矩阵一个详细的计算实例定义看起来有点抽象我们通过一个具体的、简化的双输入双输出系统来亲手算一遍感受一下这个过程。考虑一个耦合的弹簧-质量块系统或者一个简单的电路网络我们可以抽象出如下状态空间模型为了计算方便数字是假设的设系统方程为A [ -2 1 0; 1 -3 1; 0 1 -1 ]; B [ 1 0; 0 1; 1 0 ]; C [ 1 0 0; 0 1 0 ]; D [ 0 0; 0 0 ];这里状态x是3维输入u是2维输出y是2维。D矩阵为零表示没有直接馈通。我们的目标是求出G(s) C (sI - A)^{-1} B。第一步构造 (sI - A)sI - A [ s2 -1 0; -1 s3 -1; 0 -1 s1 ];第二步求 (sI - A) 的逆矩阵求一个3x3矩阵的逆我们可以用伴随矩阵法(sI - A)^{-1} adj(sI - A) / det(sI - A)。 先求行列式Δ(s) det(sI - A)Δ(s) (s2)*[(s3)(s1) - (-1)(-1)] - (-1)*[(-1)(s1) - (0)(-1)] 0*... (s2)[(s3)(s1) - 1] 1*[-(s1)] (s2)(s^2 4s 3 - 1) - (s1) (s2)(s^2 4s 2) - s - 1 s^3 4s^2 2s 2s^2 8s 4 - s - 1 s^3 6s^2 9s 3所以系统的特征多项式是s^3 6s^2 9s 3它的根就是系统的极点。再求伴随矩阵adj(sI - A)。这需要计算所有代数余子式过程略这是基本功练习。假设我们算得adj(sI - A) [ (s3)(s1)-1 (s1) 1; (s1) (s2)(s1) (s2); 1 (s2) (s2)(s3)-1 ] [ s^24s2 s1 1; s1 s^23s2 s2; 1 s2 s^25s5 ]因此(sI - A)^{-1} (1 / Δ(s)) * [ s^24s2, s1, 1; s1, s^23s2, s2; 1, s2, s^25s5 ]第三步计算 C (sI - A)^{-1} B这是一个连续的矩阵乘法。C是 2x3逆矩阵是 3x3B是 3x2。所以结果G(s)是 2x2。 我们先计算中间结果T(s) (sI - A)^{-1} B。B矩阵有两列我们分别计算。 令B1 [1; 0; 1](B的第一列)B2 [0; 1; 0](B的第二列)。计算T1(s) (sI - A)^{-1} * B1T1(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [1; (1/Δ(s)) * [ (s^24s2)*1 (s1)*0 1*1; s1, s^23s2, s2; 0; (s1)*1 (s^23s2)*0 (s2)*1; 1, s2, s^25s5 ] 1] 1*1 (s2)*0 (s^25s5)*1 ] (1/Δ(s)) * [ s^24s2 1; s1 s2; 1 s^25s5 ] (1/Δ(s)) * [ s^24s3; 2s3; s^25s6 ]同理计算T2(s) (sI - A)^{-1} * B2T2(s) (1/Δ(s)) * [ s^24s2, s1, 1; * [0; (1/Δ(s)) * [ (s^24s2)*0 (s1)*1 1*0; s1, s^23s2, s2; 1; (s1)*0 (s^23s2)*1 (s2)*0; 1, s2, s^25s5 ] 0] 1*0 (s2)*1 (s^25s5)*0 ] (1/Δ(s)) * [ s1; s^23s2; s2 ]所以T(s) [T1(s), T2(s)] (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2; s^25s6, s2 ]现在左乘C矩阵。C矩阵只取前两行因为输出y只对应前两个状态所以G(s) C * T(s)相当于取T(s)的前两行。G(s) (1/Δ(s)) * [ s^24s3, s1; 2s3, s^23s2 ]其中Δ(s) s^3 6s^2 9s 3。解读G(s) 这个 2x2 的传递函数矩阵G(s)告诉我们G_11(s) (s^24s3) / (s^36s^29s3)输入u1对输出y1的传递函数。G_12(s) (s1) / (s^36s^29s3)输入u2对输出y1的传递函数。G_21(s) (2s3) / (s^36s^29s3)输入u1对输出y2的传递函数。G_22(s) (s^23s2) / (s^36s^29s3)输入u2对输出y2的传递函数。关键观察共同分母所有四个传递函数都有相同的分母多项式Δ(s)。这个多项式正是矩阵A的特征多项式。这意味着整个MIMO系统的极点决定稳定性和基本动态是由矩阵A决定的是所有输入输出通道所共享的。这是MIMO系统与多个独立SISO系统本质不同的地方。耦合体现非对角线元素G_12(s)和G_21(s)不为零。这说明输入u1会影响输出y2输入u2也会影响输出y1。系统是耦合的。相对阶G_11和G_22的分子是二阶分母是三阶相对阶为1。G_12和G_21的分子是一阶相对阶为2。不同通道的动态特性可能不同。实操心得对于维数高于3的系统手工计算逆矩阵会非常繁琐且容易出错。在实际工程或学习中我们强烈依赖数学工具。在MATLAB/Octave中计算传递函数矩阵极其简单。假设你已经定义了A, B, C, D矩阵只需一行命令G ss(A, B, C, D); G_tf tf(G);。ss创建状态空间模型tf将其转换为传递函数矩阵形式。对于我们的例子在MATLAB中验证上述计算是很好的练习。手工推导的目的在于深刻理解其数学本质和耦合关系的来源而不是用于解决大规模问题。4. 传递函数矩阵的性质与在系统分析中的核心作用得到了传递函数矩阵我们能用它做什么它不仅仅是状态空间模型的一种等价表示更是我们分析MIMO系统特性的强大工具。4.1 极点与零点MIMO系统的扩展定义对于SISO系统极点就是传递函数分母多项式的根零点就是分子多项式的根。对于MIMO系统由于G(s)是一个矩阵极点和零点的定义需要扩展。极点如前所述传递函数矩阵G(s)的所有元素的分母多项式在约去公因子后的根就是系统的极点。更本质地说系统的极点就是系统矩阵A的特征值。这一定义从状态空间模型出发更为根本。所有输入输出通道共享同一组极点这决定了系统的稳定性所有极点均具有负实部则渐近稳定和基本响应模态如振荡频率、衰减速度。传输零点这是一个MIMO系统特有的、非常重要的概念。它不是单个元素分子为零的点。传输零点的定义是存在一个非零的复频率z和一个非零的输入向量U_0使得在零初始状态下系统的输出Y(s)恒为零。即G(z) U_0 0。这意味着在频率z处存在某种特定的输入信号组合其效果被系统内部完全“抵消”了无法在输出端被观测到。计算对于D0的系统传输零点z是使得复合矩阵P(s) [sI-A, -B; C, 0]降秩的s值。在MATLAB中可以使用tzero或zero函数直接计算。物理意义传输零点反映了系统输入输出之间的阻塞特性。例如在飞机控制中某些特定的舵面偏转组合可能无法引起飞机姿态的变化在某个频率下这个频率就是传输零点。传输零点会影响系统的可控制性、可观测性以及控制性能的极限例如非最小相位系统的右半平面零点会限制控制带宽。4.2 稳定性分析基于极点的稳定性判据在MIMO系统中依然直接适用线性时不变系统渐近稳定的充要条件是其传递函数矩阵的所有极点即矩阵A的所有特征值都具有负实部。通过传递函数矩阵G(s)我们可以计算其特征多项式即行列式det(sI-A)或G(s)各元素分母的最小公倍式然后应用劳斯判据、赫尔维茨判据或直接求根来判断稳定性。在MATLAB中使用pole(G)或eig(A)直接获取极点。注意对于MIMO系统不能仅仅通过观察G(s)某个对角元比如G_11(s)的稳定性来判断整个系统的稳定性。因为非对角元的耦合可能引入不稳定的隐藏模态与状态空间中的能控性/能观性相关。必须检查系统矩阵A的全部特征值。4.3 频域分析从伯德图到奇异值图这是传递函数矩阵威力巨大的地方。在SISO系统中我们绘制单个传递函数G(jω)的伯德图幅频和相频特性。在MIMO系统中G(jω)是一个复数矩阵。我们如何绘制它的频率响应逐个元素法可以为G(s)的每一个元素G_ij(jω)绘制伯德图。这能让我们看清每一对输入输出通道在频域的特性。这对于理解耦合的强度比如G_21相对于G_22的幅值大小很有帮助。在MATLAB中bode(G)命令会自动为所有通道生成伯德图。奇异值分析法更强大的工具对于MIMO系统更本质的频域分析工具是奇异值。对于每个频率点ω我们计算复数矩阵G(jω)的奇异值。奇异值总是非负的实数。假设G是q x p矩阵那么它在每个频率点有min(p, q)个奇异值记为σ_1(ω) ≥ σ_2(ω) ≥ ... ≥ σ_min(p,q)(ω) ≥ 0。最大奇异值σ_max(ω)可以理解为系统在频率ω处对所有可能输入方向向量的最大增益。最小奇异值σ_min(ω)可以理解为系统在频率ω处对所有可能输入方向的最小增益。绘制σ_max(ω)和σ_min(ω)随频率变化的曲线就得到了MIMO系统的奇异值伯德图。为什么重要在鲁棒控制和系统性能分析中奇异值提供了关键的洞察。例如σ_min在低频段的大小反映了系统抗干扰和解耦的能力σ_max在高频段的衰减速率反映了系统对噪声和未建模动态的鲁棒性。闭环系统的稳定裕度也可以用开环传递函数矩阵的奇异值来评估。在MATLAB中可以使用sigma(G)命令绘制奇异值图。4.4 能控性与能观性与传递函数矩阵的关联能控性和能观性是状态空间模型的核心概念它们与传递函数矩阵有着深刻的联系。能控性系统是否能在有限时间内通过合适的输入u(t)将状态从任意初始点驱动到原点。状态空间判据是能控性矩阵[B, AB, A^2B, ..., A^{n-1}B]满秩。能观性系统是否能在有限时间内通过输出的观测y(t)唯一地确定初始状态x(0)。状态空间判据是能观性矩阵[C; CA; CA^2; ...; CA^{n-1}]满秩。从传递函数矩阵的角度看如果系统是状态空间能控且能观的那么其传递函数矩阵G(s)将是既约的即没有零极点对消。如果发生了零极点对消则意味着系统不是完全能控和/或完全能观的对消掉的模态在输入输出描述中“消失”了但它们可能隐藏在系统内部影响实际动态甚至稳定性。在计算G(s) C(sI-A)^{-1}B D时如果最终得到的各元素传递函数有公因子被约去这些被约去的因子就对应着不能控或不能观的模态。因此在由传递函数矩阵反推或实现状态空间模型时必须非常小心要确保得到的是最小实现即完全能控且完全能观。5. 从理论到实践在MATLAB/SciPy中操作传递函数矩阵理论学习之后我们必须能在工具中实现它。这里以MATLAB及其开源替代品GNU Octave和Python的SciPy库为例展示基本操作。5.1 MATLAB/Octave 环境1. 定义系统与计算传递函数矩阵% 定义状态空间矩阵 (沿用之前的例子) A [-2, 1, 0; 1, -3, 1; 0, 1, -1]; B [1, 0; 0, 1; 1, 0]; C [1, 0, 0; 0, 1, 0]; D zeros(2,2); % 创建状态空间模型对象 sys_ss ss(A, B, C, D); % 转换为传递函数矩阵形式 sys_tf tf(sys_ss); disp(传递函数矩阵 G(s):); disp(sys_tf); % 输出会显示一个2x2的tf数组每个元素都是tf对象。 % 例如从输入 1 到输出 1: s^2 4 s 3 / s^3 6 s^2 9 s 3 % 从输入 2 到输出 1: s 1 / s^3 6 s^2 9 s 3 ... 等等2. 提取特定通道的传递函数% 提取输入1到输出1的SISO传递函数 G11 sys_tf(1, 1); % 提取输入2到输出2的SISO传递函数 G22 sys_tf(2, 2);3. 计算极点和零点% 计算系统极点 (即A的特征值) poles pole(sys_ss); % 或 pole(sys_tf) disp(系统极点); disp(poles); % 计算系统传输零点 zeros tzero(sys_ss); % 或 zero(sys_tf) disp(系统传输零点); disp(zeros);4. 频域分析绘图% 绘制所有通道的伯德图 (4个子图) figure; bode(sys_ss); % 或 bode(sys_tf) grid on; title(各通道伯德图); % 绘制奇异值图 (MIMO系统频域分析利器) figure; sigma(sys_ss); % 或 sigma(sys_tf) grid on; title(系统奇异值图);5. 时域仿真% 定义时间向量 t 0:0.01:10; % 定义输入信号u1是阶跃u2是正弦波 u [ones(size(t)); sin(2*t)]; % 注意维度是 2 x length(t) % 进行仿真 [y, t_out, x] lsim(sys_ss, u, t); % 绘制输出响应 figure; subplot(2,1,1); plot(t_out, y(:,1)); ylabel(y1); grid on; legend(输出1); subplot(2,1,2); plot(t_out, y(:,2)); ylabel(y2); xlabel(时间 (s)); grid on; legend(输出2);5.2 Python (SciPy Matplotlib) 环境Python在科学计算领域应用广泛控制库虽不如MATLAB专业但基础功能完备。import numpy as np import matplotlib.pyplot as plt from scipy import signal # 1. 定义状态空间矩阵 A np.array([[-2, 1, 0], [1, -3, 1], [0, 1, -1]]) B np.array([[1, 0], [0, 1], [1, 0]]) C np.array([[1, 0, 0], [0, 1, 0]]) D np.array([[0, 0], [0, 0]]) # 2. 创建状态空间系统 sys_ss signal.StateSpace(A, B, C, D) # 3. 转换为传递函数形式 (注意scipy.signal 的 tf 表示是针对SISO的) # 对于MIMO我们需要逐个元素计算或使用专门的库如 control需安装pip install control # 这里演示使用 control 库 (如果已安装) try: import control as ct # 创建控制系统库的状态空间对象 sys_ct ct.ss(A, B, C, D) # 计算传递函数矩阵 (返回一个 TransferFunction 对象列表的列表) sys_tf ct.tf(sys_ct) print(传递函数矩阵:) print(sys_tf) except ImportError: print(control 库未安装。对于MIMO传递函数建议安装 control 库。) # 手动计算其中一个通道作为示例 (例如 G11) # G(s) C(sI-A)^-1 B D, 我们计算第一个元素 # 这里省略手动计算逆矩阵的代码较为复杂。 # 作为替代我们可以直接进行时域/频域仿真。 # 4. 计算极点和零点 (使用 scipy.signal 或 control) # 极点就是 A 的特征值 poles np.linalg.eigvals(A) print(\n系统极点 (A的特征值):) print(poles) # 传输零点计算较复杂control库有提供 try: import control as ct zeros ct.zero(sys_ct) print(\n系统传输零点 (来自 control 库):) print(zeros) except: print(\n传输零点计算需要 control 库。) # 5. 时域仿真 (使用 scipy.signal) t np.arange(0, 10, 0.01) # 输入u1 阶跃 u2 正弦波 u np.vstack([np.ones_like(t), np.sin(2*t)]) # 形状 (2, len(t)) # lsim 需要输入 u 的形状为 (len(t), num_inputs) t_out, y, x signal.lsim(sys_ss, Uu.T, Tt) # 注意 u.T 转置 # 绘制结果 fig, (ax1, ax2) plt.subplots(2, 1, figsize(10, 6)) ax1.plot(t_out, y[:, 0]) ax1.set_ylabel(y1) ax1.grid(True) ax1.legend([输出1]) ax2.plot(t_out, y[:, 1]) ax2.set_ylabel(y2) ax2.set_xlabel(时间 [s]) ax2.grid(True) ax2.legend([输出2]) plt.suptitle(系统时域响应) plt.tight_layout() plt.show() # 6. 频域分析 - 伯德图 (scipy.signal 的 bode 是 SISO 的) # 我们可以循环绘制每个通道 try: import control as ct plt.figure() ct.bode_plot(sys_ct) # control库可以绘制MIMO伯德图 plt.suptitle(各通道伯德图 (control库)) plt.tight_layout() plt.show() except: print(绘制MIMO伯德图需要 control 库。)实操心得与避坑指南工具选择对于严肃的控制系统分析与设计MATLAB及其控制系统工具箱仍然是行业标杆文档齐全、函数丰富。Python的control库功能日益完善对于学习、研究和轻量级应用是完全足够的且免费开源。SciPy.signal更侧重于信号处理对MIMO支持较弱。维度匹配在定义矩阵A, B, C, D时务必反复检查维度。(n x n), (n x p), (q x n), (q x p)。一个常见的错误是B或C矩阵的维度定义反了。仿真输入使用lsim或signal.lsim进行仿真时输入矩阵U的维度容易出错。在MATLAB中U的列数等于输入数行数等于时间点数量。在SciPy的signal.lsim中U的形状是(len(t), num_inputs)。务必查阅文档或打印矩阵形状来确认。最小实现如果你是从物理方程推导或系统辨识得到的传递函数矩阵在转换为状态空间模型时使用minreal函数MATLAB或minreal方法control库来获取最小实现消除不能控或不能观的模态避免后续分析与设计出现问题。6. 传递函数矩阵的局限性与状态空间法的优势尽管传递函数矩阵是分析MIMO系统的有力工具但它并非万能也存在固有的局限性。理解这些局限性正是我们深入学习现代控制理论状态空间法的动力。1. 仅适用于线性时不变系统传递函数建立在拉普拉斯变换的基础上其核心是叠加原理和常系数线性微分方程。对于非线性系统、时变系统传递函数矩阵的定义不再成立。而状态空间方程ẋ f(x, u, t), y g(x, u, t)在形式上可以描述更广泛的系统。2. 丢失内部状态信息只反映输入输出关系这是传递函数方法最根本的局限。G(s)描述的是“黑箱”的外部特性。如果系统不是完全能控和完全能观的那么传递函数矩阵就无法反映系统的全部动态。那些不能控或不能观的模态在G(s)中会被对消掉。但从内部看这些隐藏模态可能是不稳定的会导致实际系统出现问题。状态空间法通过直接研究状态向量x(t)能完整揭示系统的内部行为。3. 难以处理非零初始条件传递函数分析通常假设零初始条件。对于非零初始状态的系统响应传递函数方法处理起来比较麻烦。而状态空间方程结合初始状态x(0)可以非常自然地求解全响应。4. 多变量系统分析与综合的复杂性虽然传递函数矩阵将多变量关系封装了起来但在进行控制器设计时如经典的频域设计法直接处理一个矩阵函数仍然非常复杂。例如如何为MIMO系统设计一个PID控制器每个通道单独设计往往会因为耦合而失败。而状态空间法通过状态反馈u -Kx和观测器设计提供了系统化的多变量控制器设计框架如线性二次型调节器LQR、极点配置等概念上更统一。5. 对系统结构的洞察较弱状态空间表示法通过A, B, C, D矩阵清晰地分离了系统动态(A)、输入影响(B)、输出测量(C)和直接传递(D)。这种结构化的表示更容易与物理模型对应也便于进行能控性、能观性等结构性质的分析。那么为什么我们还要学习传递函数矩阵因为它是连接经典控制与现代控制的桥梁。它保留了频域分析的直观性伯德图、奈奎斯特图为理解系统的频率响应、带宽、鲁棒稳定性提供了图形化工具。许多先进的多变量频域设计方法如H∞鲁棒控制也是建立在传递函数矩阵的基础之上。在实际工程中常常是频域方法基于传递函数矩阵和时域方法基于状态空间结合使用取长补短。因此传递函数矩阵是现代控制理论工具箱中不可或缺的一部分。它让我们能够用熟悉的频域语言去理解和分析复杂的多变量系统同时又提醒我们其边界所在引导我们走向更深刻、更强大的状态空间分析领域。掌握了它你就拿到了理解和处理现实世界中耦合、多变量系统的第一把钥匙。
