多智能体系统故障检测：基于中间变量观测器的分布式估计方法

简介面向多智能体系统与故障诊断领域研究者的技术资料聚焦无向拓扑下线性多智能体系统的执行器故障检测问题。内容基于中间变量观测器方法通过构建维数扩展的虚拟系统同时实现状态与故障的联合估计有效规避传统观测器对匹配条件的依赖并利用分布式残差完成自身及邻域智能体的故障判别。理论部分依托Lyapunov稳定性分析给出估计误差有界性证明附带的Python代码则完整覆盖系统建模、观测器增益设计、故障检测逻辑与仿真验证包含虚拟系统构建、中间变量技术、LMI设计等关键实现。包体为单个docx文件压缩包大小51KB将理论推导与可运行代码整合于同一文档便于对照阅读与复现。目前已有56人浏览学习适合具备自动控制理论基础、希望深入理解中间变量观测器设计及多智能体故障检测算法的研究人员或工程师参考使用。通过调整系统参数与通信拓扑可进一步评估算法在不同场景下的鲁棒性与适用性。1. 线性多智能体系统故障检测中间变量观测器到底拆掉了什么做线性多智能体系统故障检测绕不开一个尴尬每个节点只能看到自己的量测可故障经常出在邻居身上或者出在自己身上但被扰动盖住。直接给每个节点配一个状态观测器估计结果往往跟着故障一起漂误报漏报一起上。中间变量观测器把故障从状态方程里显式拆出来变成“先估计中间变量再反推故障”的两步走配合无向拓扑下的邻居信息交换就能做到状态与故障联合估计再落回分布式残差检测。这套方法不需要全局信息适合编队控制、微电网分布式电源、车联网协同这类场景。本文按建模、增益求解、MATLAB 仿真、残差检测和调参避坑的顺序展开代码可以直接改成自己的系统矩阵。2. 无向拓扑下的问题建模拉普拉斯矩阵与中间变量的引入建模不提清楚后面全是白干。这一章先把数学模型立住重点回答三个问题无向图用什么矩阵描述、中间变量为什么能让故障“可见”、联合估计能做的条件是什么。很多新手上来就写 MATLAB结果观测器不收敛回头才发现是图不连通或者增广系统不可观这一步省不得。2.1 无向图与拉普拉斯矩阵连通性为什么是收敛的前提无向拓扑用图 G(V,E) 描述节点集合 V 对应智能体边集合 E 对应通信链路。图的代数性质由拉普拉斯矩阵 L 刻画L D - A_g其中 A_g 是邻接矩阵A_g(i,j)1 表示 i 与 j 能通信无向图保证 A_g 对称D 是度矩阵D(i,i) 等于节点 i 的邻居数量。这里有个需要盯住的细节邻接矩阵与后面系统矩阵 A 撞名我在代码里统一用 A_g 表示图邻接用 A 表示单智能体动态矩阵避免混用。L 的特征值直接决定分布式观测器能不能收敛。无向连通图有一个关键性质0 是 L 的单特征值对应特征向量是全 1 向量其余特征值都大于 0其中第二小特征值 λ2 叫代数连通度。对 4 节点环加对角线的拓扑L 的特征值是 0、2、2、4λ22收敛速度的下限就由它决定。如果图断开0 特征值会出现多个每个连通子图各估计各的故障位置永远对不齐。所以拿到一个拓扑第一件事就是打印eig(L)看正特征值个数少一个都直接换拓扑不值得往后跑。2.2 中间变量观测器的数学变换从“估计故障”到“估计中间变量”先给单节点模型。每个智能体 i 的动态是线性时不变系统带一个故障输入 f_i 和一个扰动 d_i[ \dot{x}_i A x_i B u_i D f_i d_i, \quad y_i C x_i v_i ]x_i 是 n 维状态y_i 是 m 维量测f_i 是故障信号常值或缓变比如执行器卡死、增益突变v_i 是量测噪声。故障通过输入矩阵 D 进入状态方程目标是只靠本地输出 y_i 和邻居的估计值把 x_i 和 f_i 同时估出来。难点在 f_i 没有自己的动态方程。如果直接把 f 当成状态扩进去增广系统是[ \begin{bmatrix}\dot{x}\ \dot{f}\end{bmatrix} \begin{bmatrix}A D\ 0 0\end{bmatrix} \begin{bmatrix}x\ f\end{bmatrix} \begin{bmatrix}B\ 0\end{bmatrix} u ]这个形式看着简洁但增广后的能观性经常不满足尤其是 D 的方向与 C 的观测方向不匹配时故障对输出其实不可见。中间变量观测器的做法是引入一个变换把故障“编码”进一个新状态里[ \tau_i f_i - K x_i ]K 是需要设计的矩阵作用相当于换基底。反解出 f_i τ_i Kx_i代回状态方程得到[ \dot{x}_i (A D K) x_i B u_i D \tau_i d_i ]对 τ_i 求导假设故障是常值或缓慢变化即 ḟ_i ≈ 0[ \dot{\tau}_i -K \left( (A D K) x_i B u_i D \tau_i \right) ]合并成增广形式令 z_i [x_i; τ_i]得到[ \dot{z}_i A_a z_i B_a u_i, \quad y_i C_a z_i v_i ]其中[ A_a \begin{bmatrix} A D K D \ -K(A D K) -K D \end{bmatrix}, \quad B_a \begin{bmatrix} B \ -K B \end{bmatrix}, \quad C_a \begin{bmatrix} C 0 \end{bmatrix} ]这个变换的价值在于原来故障 f 既没有动态方程也不能直接测量现在它被编码进 τ而 τ 作为增广状态有了明确的动态方程观测器只需估计 z最后用 f̂ τ̂ Kx̂ 反算故障。K 成了设计师手里的旋钮可以改变 A_a 的结构把故障信息“搬”到可观测的通道上这是直接增广法没有的余地。2.3 状态与故障联合估计的可观测性条件观测器能工作的前提是 (C_a, A_a) 能观。用能观性矩阵判定最省事rank(obsv(Aa, Ca)) 必须等于增广状态维数 n1。对单故障标量 f增广维数是 n1秩不够就说明故障方向有一部分落在不可观子空间里输出修正拉不回来。K 的作用在这里体现得最明显。直接增广相当于 K0A_a 变成块三角矩阵能观性完全由原系统 (C, A) 和故障方向 D 决定。引入 K 之后A_a 的左上块变成 AD K相当于把故障通道的“影子”注入原动态即使 C 测不到 D 的方向只要 K 把 D 的信息混进可测状态能观性就可能被救回来。所以调 K 不是玄学本质是在调整故障到输出的信息通道。3. 用 MATLAB 实现分布式联合估计增益求解与主循环代码这一章直接给能跑的 MATLAB 代码。系统取 4 个节点、单智能体 4 维状态、2 维量测故障是标量。代码分四步初始化图与系统、求解中间变量增益 K、求解观测器增益 L、跑仿真主循环。每一步之后都解释参数怎么改以及改动会带来什么连锁反应。3.1 分布式观测器方程本地项与邻居共识项每个节点上的观测器方程是同一个模板包含三个部分[ \dot{\hat{z}}_i A_a \hat{z}i B_a u_i L_i (y_i - C_a \hat{z}i) \theta \sum{j \in N_i} a{ij} (\hat{z}_j - \hat{z}_i) ]第一项是本地模型传播告诉观测器系统“正常时”怎么动第二项是输出新息修正用真实量测和估计输出的差把估计拉回真实轨迹第三项是邻居共识项把邻居的估计值与本地的估计值做差乘耦合强度 θ 后作为额外修正。要特别强调节点之间交换的是估计向量 z_j不是量测 y_j也不是真实状态。这个限制是分布式方法和集中式方法的本质区别代码里要时刻检查每个变量到底来自哪里。θ 的取值直接影响收敛行为。太小邻居信息拉不齐各节点估计各的太大共识项变成强反馈估计值可能高频震荡。我的习惯是从 θ10 起步画误差曲线看趋势再调后面第五章专门讲这个坑。3.2 初始化代码图拓扑、系统矩阵与故障注入% 1) 无向连通图4节点“环对角线”保证拉普拉斯只有1个零特征值 N 4; A_g zeros(N); A_g(1,2)1; A_g(2,3)1; A_g(3,4)1; A_g(1,4)1; A_g A_g A_g; L_g diag(sum(A_g,2)) - A_g; fprintf(L 特征值: %g\n, eig(L_g).); % 2) 单智能体开环模型4状态2输出故障通道D在最后一个状态 A [0 1 0 0; -1 -0.6 1 0; 0 0 0 1; 0 0 -2 -1]; B [0; 1; 0; 1]; C [1 0 0 0; 0 1 0 0]; D [0; 0; 0; 1]; % 故障只作用于第4个状态方程 u 0; % 零输入专注看故障检测能力这里选了个“有心机”的组合C 能测到前两个状态D 却作用在第四个状态上直接增广大概率不可观正好演示中间变量观测器的价值。故障注入函数单独定义用一个缓变的斜升信号模拟真实的缓慢失效过程% 故障5秒后开始斜升10秒后稳定在1 fault (t) (t 5) .* min(t - 5, 1);用缓变信号而不是阶跃是为了验证中间变量观测器对“慢变故障”假设的容忍度这也是后面 5.2 节要处理的问题。3.3 中间变量增益 K 与观测器增益 L 的求解% 3) 中间变量增益 K把故障信息混进可测状态方向 K [2 0 0 1]; % 1x4将故障与x1、x4混合 Aa [A D*K, D; -K*(A D*K), -K*D]; Ba [B; -K*B]; Ca [C, zeros(2,1)]; % 能观性检查增广状态共5维秩必须是5 fprintf(obsv rank %d\n, rank(obsv(Aa, Ca))); assert(rank(obsv(Aa, Ca)) 5, 增广系统不可观换K或加传感器); % 4) 观测器增益 L极点配置极点略快于系统动态 L place(Aa., Ca.).;K 是 1×4 行向量因为故障是标量τ 也是标量。K 取 [2 0 0 1] 的含义是让故障与 x1 和 x4 混合x1 恰好是第一路量测能直接看到的状态这样故障信息就通过 K 进入到了可观测通道。如果你把 K 改成 [0 0 0 1]对应直接增广obsv 矩阵可能会降秩这就是 5.1 节要展开的典型案例。place是极点配置函数要求对 (Aa, Ca) 能观所以先转置再求最后再转置回来。如果place报错基本都是能观性检查没通过回头调 K不要硬撑。place 默认选的极点可能落在系统动态附近想更快跟踪可以显式给极点向量比如place(Aa., Ca., [-3 -4 -5 -6 -7]).代价是对噪声更敏感这个权衡要靠仿真效果定。3.4 仿真主循环状态与故障联合估计% 5) 仿真主循环欧拉积分步长0.01 dt 0.01; T 20; nT T / dt; x zeros(4, N, nT); % 真实状态 z zeros(5, N, nT); % 观测器增广状态 [x; tau] x(:, :, 1) 0.2 * ones(4, N); % 真实状态初值偏离0 z(:, :, 1) zeros(5, N); % 观测器从0开始被迫收敛 theta 10; for k 2:nT t (k - 1) * dt; f_now fault(t); for i 1:N % 真实节点各自独立运行只有本地故障 x(:, i, k) x(:, i, k-1) dt * (A * x(:, i, k-1) B * u D * f_now); % 观测器输出修正 邻居共识 innov C * x(:, i, k) - Ca * z(:, i, k-1); sum_n zeros(5, 1); for j 1:N if A_g(i, j) sum_n sum_n (z(:, j, k-1) - z(:, i, k-1)); end end z(:, i, k) z(:, i, k-1) dt * (Aa*z(:, i, k-1) Ba*u L*innov theta*sum_n); end end % 6) 从中间变量反算故障估计f_hat tau_hat K * x_hat f_hat zeros(N, nT); for i 1:N f_hat(i, :) squeeze(z(5, i, :)). K * squeeze(z(1:4, i, :)); end % 画图对比plot(t, f_hat); hold on; plot(t, fault(t), k--);真实状态更新和观测器更新分开写是为了明确信息边界真实节点只用自己的 x 和 f观测器只用自己的 y 和邻居的 z 估计。代码里 z(:,j,k-1) 用的是上一拍邻居估计值叫零阶保持符合分布式系统里通信一拍延迟的物理现实如果你用同一拍 k 的数据相当于同步更新假设实验结果会比真实系统乐观这个细节很容易被忽略。欧拉积分步长取 0.01 对这套极点位置够用如果把系统矩阵改得更快特征值更大步长要相应缩小否则数值发散。跑完看 f_hat 曲线理想情况 5 秒故障出现后估计值应该在一个短暂暂态后跟上真实故障滞后时间就是检测延迟的基础。4. 分布式残差检测设计阈值、滑动窗口与报警逻辑观测器收敛后状态和故障估计都有了但这还不能直接叫“检测”。检测意味着要从残差里找出一个自动判据告诉系统“现在该报故障了”。这一章把残差生成、阈值设计、报警逻辑三块分开做每块都有对应的 MATLAB 代码和调参建议。4.1 残差生成新息残差与邻居加权残差最自然的残差是本地新息也就是观测器输出与真实输出的差[ r_i(t) y_i(t) - C_a \hat{z}_i(t) ]无故障时观测器收敛后 r_i 应该趋近于零剩下的只是量测噪声故障出现后估计误差被故障推偏r_i 就会出现系统性偏差。但这里有个问题故障如果发生在节点 j而 j 的观测器通过邻居项会把偏差传给邻居 i所以单独看一个节点的残差很难区分“我自己坏了”还是“邻居坏了我被带偏了”。我的做法是对本地残差做一个邻居加权形成分布式残差[ \bar{r}i(t) r_i(t) \beta \sum{j \in N_i} a_{ij} (r_i(t) - r_j(t)) ]β 通常取 1。这项的作用是放大“差异化”如果只有节点 i 的残差偏差明显而邻居正常加权项会进一步放大这个差异如果所有节点同时被共同噪声推动比如同一时刻的外部扰动差异项会互相抵消。这样无故障时的公共偏差不容易误报故障时的局部偏差更容易被抓住。4.2 阈值设计卡方门限与 3σ 经验值残差是多维向量直接比大小不靠谱我的做法是先转成标量统计量。假设量测噪声 v_i 服从零均值高斯分布协方差 Σ_v那么无故障时残差近似服从同样分布加权后的统计量[ J_i(t) \bar{r}_i^\top \Sigma_v^{-1} \bar{r}_i ]服从自由度为 m 的卡方分布m 等于 C 的行数也就是量测维数。给定误报率 α阈值可以查卡方表得到。以 2 维量测为例α0.01 时阈值是 9.21α0.001 时是 13.82。这个表可以直接用误报率 α卡方阈值自由度2适用场景0.055.99对延迟敏感可容忍误报0.019.21常规工业检测0.00113.82误报代价高的场合理论卡方阈值有个前提残差是独立同分布的高斯白噪声。实际仿真里残差往往时间相关直接套理论值会偏乐观。我的习惯是先用理论值跑一遍无故障工况统计 J 的实际分布取 99% 分位当阈值相当于用数据给理论值打补丁。如果两者差太多说明残差里还有未建模的动态要先回第三章查观测器收敛质量而不是盲目放宽阈值。4.3 检测判定代码滑动窗口加连续计数% 检测统计量和报警逻辑 Sigma_v 0.01 * eye(2); % 量测噪声协方差 alpha 0.01; thres chi2inv(1 - alpha, size(C, 1)); % 自由度2阈值9.21 win 10; % 滑动窗口长度 cnt zeros(N, 1); alarm zeros(N, nT); for k 1:nT for i 1:N % 取第i节点k时刻的分布式残差算马氏距离类型的统计量 rr squeeze(r_bar(k, i, :)); % 2x1 J rr / Sigma_v * rr; if J thres cnt(i) cnt(i) 1; else cnt(i) max(0, cnt(i) - 1); % 回滞不立刻清零 end alarm(i, k) cnt(i) 3; % 连续3拍超阈才确认 end end计数器写法有三个细节值得说明。第一用cnt 3而不是J thres单拍判断是为了滤掉单点毛刺第二未超阈时cnt减 1 而不是清零这叫回滞让持续故障的报警能稳定保持不会因为一两拍小于阈值就抖动第三win10这个参数对应“连续 10 拍内有 3 拍超阈”才确认实际效果相当于一个低通滤波把检测延迟和误报率绑在一起调。想要更快报警就减小 win 或把连续次数从 3 降到 2代价是误报概率上升想要更稳就反过来。这个权衡没有标准答案只能按你们系统的误报代价来定。5. 避坑与排查无向拓扑联合估计的 5 个典型问题这一章全部来自实际调试过程每条都是“现象 → 原因 → 解决”的结构。能避开的坑提前避开避不开的至少要知道往哪个方向查。5.1 故障估计发散但状态收敛问题多半出在可观测性现象状态估计跟随得很好误差曲线贴着零走但 f_hat 数值越来越大甚至到十几。原因增广系统的能观性秩不足状态方向勉强能观故障方向却落在不可观子空间观测器对 f 没有修正能力只能靠模型积分慢慢漂。解决打印rank(obsv(Aa, Ca))如果小于 n1调 K。一个非常直观的例子把上一章代码里的 K 从[2 0 0 1]改成[0 0 0 1]运行后 f_hat 就会发散。原因是 K 只有 x4 分量而 x4 不在 C 的观测范围内故障信息没有进入可测通道。改回[2 0 0 1]把故障混进 x1立刻收敛。每次换系统矩阵后先跑这一条检查再继续能省一整晚。5.2 缓变故障估计滞后中间变量的慢变假设被打破现象阶跃故障估计收敛很快斜坡故障的估计值总是差一截检测延迟变大残差里像是有一块“甩不掉”的偏差。原因中间变量观测器推导时假设 ḟ ≈ 0斜坡故障的导数是个常数相当于在 τ 的动态方程里注入了一个未建模输入。解决最实用的一招是把故障模型加一阶积分链把 f 的导数也当成状态。具体做法是重新定义中间变量为 τ [f; ḟ] - Kx增广状态变成 [x; f; ḟ]维数加一代价是 LMI 和矩阵维度都变大。对慢变故障也可以先接受滞后在检测逻辑里加一个故障发生时刻补偿但对大多数要求不高的情况滞后几百毫秒在可接受范围。5.3 邻居耦合项让估计震荡θ 不是越大越好现象θ 从 10 调到 50状态估计开始高频抖动甚至出现局部振荡。原因共识项相当于在每个节点上引入一个与邻居估计差成正比的负反馈回路θ 过大时这个反馈环路的增益过高和观测器本身的极点互相作用形成振荡。解决固定 L单独扫 θ从 1 到 50 步长 1画状态估计误差的 RMS选最低点。更系统的方式是让 θ 和拉普拉斯矩阵的第二小特征值挂钩取 θ c / λ2c 从 1 扫到 10。在 4 节点示例里 λ22c20 对应 θ10和代码里的取值一致。如果换了拓扑特征值变了 θ 却不变这就是翻车点。5.4 误报压不住阈值不能拍脑袋现象无故障工况下报警器频繁触发把阈值调大后确实不报了但真正故障时灵敏度也不够了。原因理论卡方阈值假设残差是白噪声但实际上观测器没有完全收敛时残差里的暂态分量会让 J 统计量远高于稳态值这部分“合法”偏差被误判成故障。解决用离线蒙特卡洛标定阈值。跑 10 次无故障仿真每次用不同的随机噪声种子把每次的 J 值全部收集起来取 99% 分位当阈值。这个数据驱动的阈值天然包含了暂态影响比理论卡方值可靠得多。注意 J 值要跨节点收集因为分布式残差在节点间相关只看单节点会低估整体误报率。5.5 仿真里偷用了全局信息分布式实现的自检现象离线仿真效果很好检测延迟 0.2 秒误报率极低换到半实物平台或者和通信模块联调性能直接崩。原因仿真代码里很可能用了全局变量或未来信息。常见偷用包括观测器更新时取了真实状态 x(:,j,k) 当邻居状态、残差计算时用了故障注入信号 f_now、或者用零延迟同步通信代替一拍延迟。解决给代码做一个“信息边界”检查。列出每个变量来源只允许两类数据进入观测器本地量测 y_i 和上一拍邻居估计 z_j。我在实际项目里会把通信封装成一个黑匣子函数只暴露这两个入参谁想拿别的数据都拿不到。这个习惯能提前暴露大部分分布式实现问题。6. 验证与进阶用蒙特卡洛把检测延迟和误报率量化到最后这一步要做的不是看一张仿真图说“看起来收敛了”而是把检测性能量化成三个数检出率、误报率、检测延迟。蒙特卡洛的思路很简单用随机化把偶然性摊平。故障工况跑 50 次每次故障发生时刻在 5 到 10 秒之间随机取噪声种子随机统计报警信号从故障注入到首次确认报警的平均时间就是平均检测延迟无故障工况同样跑 50 次统计报警次数占总运行拍数的比例就是误报率。MATLAB 里把第三章的主循环包成函数输入是系统矩阵、K、θ、阈值输出是报警时间序列外层循环调 50 次即可。这样做的好处是改参数后不需要肉眼判断好坏直接看数字对比。我个人的验证顺序是先跑无故障蒙特卡洛标定误报率确认阈值贴着 1% 走再跑故障蒙特卡洛看检测延迟如果延迟超过系统能容忍的上限优先调 K 和观测器极点而不是动阈值最后对比不同拓扑——环状、星状、完全图——确认拓扑变化没有让性能崩盘。无向连通图里代数连通度越大收敛越快但检测延迟不一定成比例下降因为故障检测主要靠本地残差共识项带来的增益有限。这个结论让我在做方案评审时能顶住“拓扑越密越好”的直觉压力。蒙特卡洛跑完还有一个隐藏收益它能暴露参数敏感度。如果阈值从 9.2 调到 10.5 检出率就崩了说明残差统计量落在阈值边缘附近鲁棒性不足需要回头改善观测器而不是继续调阈值如果阈值从 9 到 15 性能都很平说明系统裕度大可以大胆放宽。我现在每换一套系统矩阵固定动作就是先跑能观性秩、再跑一次 20 轮蒙特卡洛这两个检查不过关后面调任何参数都是白费力气。希望帮到你。本文还有配套的精品资源点击获取