简介面向多智能体系统初学者与MATLAB仿真实践者的示例代码包围绕MAS中代理建模、环境感知、消息交互、协调决策与仿真验证等核心环节展开适合希望用脚本方式理解多代理协作机制、搭建小型实验原型的读者。压缩包共4个文件以3个m脚本和1个txt许可说明为主体量约4KB轻量便于查阅与二次修改脚本承担代理状态更新、位置计算及服务端协调等示例功能许可文件用于说明使用授权。已有1806人学习下载可作为课程演示、课题预研或自学练手的参考。内容侧重多智能体服务器如何注册管理代理、分配任务、维护通信通道并定义协同规则同时保留监控评估的扩展空间读者可据此观察服务端创建与跟踪代理实例、按能力派发任务的过程并把代理、环境、交互、协调与决策等概念落到可运行的MATLAB代码骨架上。1. 多智能体 MATLAB 代码拆包先跑通 5 个智能体的编队收敛再谈改协议拿到一份多智能体 MATLAB 代码多数人的第一反应是直接 run 主脚本看那几条曲线收不收敛——收敛就收工散了就怀疑算法写错。我的经验恰好相反MAS 仿真翻车八成不是协议本身的问题而是通信拓扑、采样步长、初始状态这三处没对齐。这套代码真正的价值是把一致性协议、编队偏移、分布式控制器和仿真外壳拆成了互相独立、可单独替换的模块改一处不牵连其余。它适合正在做编队、蜂拥、分布式优化验证的人也适合想把 Simulink 里那堆手工连线换成脚本化流程、方便批量跑参数的场景。2. 一致性协议落到矩阵上拉普拉斯矩阵的三种写法与增益边界多智能体系统在代码层面的全部秘密几乎都压缩在一个矩阵里拉普拉斯矩阵 L。协议写成微分方程是 ẋ -Lx写成代码就是一次矩阵乘法。所以这一章不急着调控制器先把 L 造对——L 错了后面所有参数整定都是白费功夫。2.1 邻接矩阵的三种建法邻接矩阵 A 描述“谁能收到谁的信息”是整个 MAS 的地基。小规模固定拓扑手写边表最直观中等规模按几何距离生成通信半径内的智能体连边大规模或有现成拓扑数据时从稀疏矩阵导入。三种写法各有各的翻车点手写最容易犯的错是只赋值单向。% 建法一手写边表适合 N 10 的固定拓扑 N 5; edges [1 2; 2 3; 3 4; 4 5; 5 1; 1 3]; % 每行一条无向边 A zeros(N); for k 1:size(edges, 1) i edges(k, 1); j edges(k, 2); A(i, j) 1; A(j, i) 1; % 无向图必须成对赋值漏一行就变成有向图 end % 建法二按通信半径生成适合随机撒点场景 function A radius_graph(P, R) N size(P, 1); D squareform(pdist(P)); % 两两欧氏距离 A double(D R D 0); % 排除自身D0 是关键 end % 建法三从稀疏表恢复适合大规模或外部拓扑 % A sparse(rowIdx, colIdx, 1, N, N); A full(A | A);第一种写法里那句成对赋值不是啰嗦是我踩过的坑早期图省事只写 A(i,j)1结果 L 不对称特征值出现负实部曲线前半段看着挺正常后半段直接飞出去。第二种写法的D 0同理——不加这一条对角线会被判成“距离 0 以内”自我连边L 的对角元凭空多出 N收敛速度全乱。第三种写法补一句A | A做对称化防止上游数据只给了半边。提示建完 A 之后立刻断言isequal(A, A)。这一步成本几乎为零能挡掉后面一大半莫名其妙的发散。2.2 从拉普拉斯矩阵到状态方程有了 A拉普拉斯矩阵一行就能算出来L D - AD 是度矩阵。这个式子的物理含义是“我拿自己和所有邻居的差做加权和”所以 L 的每一行和恒为零全 1 向量永远是它的零特征值对应向量——这正是所有智能体最终停在同一值的数学来源。收敛速度由第二小特征值 λ₂ 决定代数连通度越高收敛越快。D diag(sum(A, 2)); % 度矩阵按行求和 L D - A; % 拉普拉斯矩阵 % 三条自检建议写进每个脚本开头 assert(isequal(A, A), 邻接矩阵不对称); assert(all(abs(sum(L, 2)) 1e-12), 拉普拉斯行和不为零); lam eig(L); lam2 sort(real(lam)); fprintf(代数连通度 lambda_2 %.4f\n, lam2(min(2, numel(lam2))));三条自检对应三类错误不对称会让 L 出现复特征值行和不为零说明 D 求错了维度λ₂ 接近零说明图不连通——不连通时一致性协议不是发散而是分裂成几个各自收敛的簇这种“部分收敛”最迷惑人曲线看着像要收敛又永远差一点。连续时间协议 ẋ -Lx 在代码里通常用前向欧拉离散x(k1) (I - h·L)x(k)稳定条件是 h 2/λ_max(L)。这个边界不是理论摆设我见过 h 设成 0.5、λ_max 大约 4 的情形正好压在边上前几步正常二十步后开始振荡放大。2.3 五智能体最小可跑脚本把上面的东西拼起来就是一段不到 20 行、但能完整验证协议的最小脚本。跑它的意义在于在引入编队、饱和、时延之前先确认你的拓扑和步长是干净的。N 5; A zeros(N); edges [1 2; 2 3; 3 4; 4 5; 5 1; 1 3]; for k 1:size(edges, 1) A(edges(k,1), edges(k,2)) 1; A(edges(k,2), edges(k,1)) 1; end L diag(sum(A,2)) - A; x0 [1 3 2 5 4]; % 初始状态故意打散 h 0.05; % 采样步长 T 20; % 仿真时长 K round(T / h); X zeros(N, K1); X(:,1) x0; for k 1:K X(:,k1) (eye(N) - h*L) * X(:,k); end xf mean(x0); % 无向连通图的收敛值即初始均值 err max(abs(X(:,end) - xf)); fprintf(末态最大偏差 %.3e理论 lambda_2 %.4f\n, err, ... sort(eig(L))(2));逻辑是一条直线加权平均 → 矩阵乘 → 迭代。参数上有三处会直接影响结果。h必须小于 2/λ_max实用做法是取 0.1/λ_max宁可多跑几步也别贴边界。T要留足经验值是 10/λ₂ 以上λ₂ 越小需要的时间越长图越“细长”越慢。x0最好手动打散别用rand直接生成——调试阶段你需要能心算出期望均值随机初始值会让你失去这个参照。收敛判据也别只看最后一步。我在脚本里习惯再算一条errHist max(abs(X - mean(x0)), [], 1)把误差历史画出来。真收敛是指数下降的直线半对数坐标下如果误差在某个量级上压不下去多半是 λ₂ 太小或步长偏大导致的残余振荡光看末态一个数看不出来。3. 编队控制的可调旋钮偏移量、饱和限幅与拓扑切换一致性只要求大家“停在一起”编队要求“停成某个形状”。从一致性到编队代码上只多了一个偏移量矩阵但工程上的复杂度翻了一倍——因为偏移量引入了几何约束几何约束又逼着你处理饱和、避碰和拓扑变化。3.1 把方针翻译成偏移量矩阵编队协议的标准写法是 u_i -Σ a_ij[(x_i - p_i) - (x_j - p_j)]其中 p_i 是第 i 个智能体的期望相对位置。写成矩阵就是 u -L(X - P)。这里最容易搞混的是 P 的语义它给的是相对偏移形心在原点还是绝对目标位置代码里必须统一否则整个队形会整体漂移到画面外。% 二维平面上的五个位置pentagon 队形形心在原点 P [ 0.0 1.0; 0.95 0.31; 0.59 -0.81; -0.59 -0.81; -0.95 0.31]; % N×2第 i 行是第 i 个智能体的期望偏移 function Xn formation_step(X, P, L, h) % X: N×2 当前状态P: N×2 期望偏移L: N×N 拉普拉斯 E X - P; % 先减偏移再走一致性 Xn X - h * (L * E); endE X - P这一行的顺序不能反。写成L*X - P是另一个东西那是带绝对参考的跟踪写反了队形会一边收敛一边整体平移看起来像“慢悠悠地飘走”。另外 P 的行数必须等于 N我在批量换 N 的脚本里吃过亏P 是硬编码的 5 行N 改成 7 之后 MATLAB 不报错而是自动做隐式扩展队形直接算成一团乱麻——这种不报错的 bug 最难查。注意换 N 的时候A、L、P、x0 四处的维度必须一起改。建议在函数入口加assert(size(P,1) size(L,1))一秒定位。3.2 饱和限幅与避碰项两个位置别写反真实执行器有速度上限仿真里不加饱和编队初期会出现几十米每秒的瞬时速度看着收敛很漂亮换成实机直接失控。限幅的写法有讲究要等比缩放不能逐分量裁剪否则方向会变队形走形。function Xn step_sat(X, P, L, h, vmax) u -L * (X - P); % 原始速度指令N×2 n sqrt(sum(u.^2, 2)); % 每个智能体的指令模长 scale min(1, vmax ./ max(n, 1e-9)); % 超过上限才缩放 u u .* scale; % 整体缩放保持方向 Xn X h * u; endmax(n, 1e-9)是防零除的后悔药初始位置如果和期望位置完全重合指令为零0/0 会给你 NaN之后整条曲线全是 NaN 且不报错。scale用min(1, ...)保证只往下压不放上去这一点在调试时很关键——如果发现某个智能体始终贴着 vmax 跑说明它被拓扑甩开了问题在图不在限幅。避碰项和限幅是两个独立的东西别塞进同一个 if 里。常见做法是在 u 后面叠加一个斥力项距离小于安全半径 d_safe 时沿相对位置方向加一个反比项。斥力增益不能给大否则会和一致性项打架出现两个智能体高频抖动、位置几乎不动的现象——这种抖动在曲线上看是毛刺在相图上看是一个极限环。3.3 采样步长与拓扑切换的稳定边界固定拓扑下 h 2/λ_max 就够了一旦拓扑会切换通信中断、邻居变化稳定条件变成对所有可能拓扑的并集都成立也就是 h 2/max(λ_max over all modes)。这条规则在实际代码里经常被忽略因为单跑某个拓扑是稳的一切换就炸。% 拓扑集合正常、断边、链路退化三种模式 Amodes {A1, A2, A3}; hmax inf; for m 1:numel(Amodes) Lm diag(sum(Amodes{m},2)) - Amodes{m}; hmax min(hmax, 2 / max(real(eig(Lm)))); end fprintf(拓扑切换下的步长上限 %.4f建议取 %.4f\n, hmax, 0.1*hmax);这段循环算的是“最保守的公共步长”。实用取值我一般按 0.1 倍上限代价是仿真步数变多但换来的是切换瞬间不会出现数值爆点。如果确实需要大 step 提速那就得改成严格离散协议用行随机矩阵而不是 I - hL牺牲一点物理直觉换取无条件稳定。切换逻辑本身还有一个隐藏坑切换时刻落在两个采样点中间时状态更新用的是旧拓扑还是新拓扑仿真里通常按采样点对齐切换只在 k 为整数时生效。如果你的切换由事件触发比如距离阈值一定要把触发时刻量化到采样网格否则同一份代码跑两次结果不同排查起来毫无头绪。4. 仿真工程化参数扫描、数据归档与批量绘图单次仿真跑通只是及格线。真要做对比、写报告、给导师或同事看结论你需要的是“一次配置、批量跑完、结果能追溯”。这一章把脚本改造成可复用的三层结构仿真函数、扫描脚本、绘图脚本。4.1 把单次仿真封成一个函数函数化的第一原则是输入只有参数输出只有数据中间不画图、不打印、不依赖工作区变量。这样它才能被parfor调用也才能被后来的自己信任。function out sim_consensus(cfg) % cfg 字段N, edges, x0, h, T, P, vmax N cfg.N; A zeros(N); for k 1:size(cfg.edges, 1) i cfg.edges(k,1); j cfg.edges(k,2); A(i,j) 1; A(j,i) 1; end L diag(sum(A,2)) - A; K round(cfg.T / cfg.h); X zeros(N, 2, K1); X(:,:,1) cfg.x0; for k 1:K X(:,:,k1) step_sat(squeeze(X(:,:,k)), cfg.P, L, cfg.h, cfg.vmax); end out.X X; out.L L; out.lam2 sort(eig(L)); out.cfg cfg; % 把配置原样存回去方便追溯 out.err squeeze(max(vecnorm(X - cfg.P, 2, 2), [], 1)); end三维数组X(N,2,K1)是这里唯一需要适应的地方第一维智能体、第二维状态分量、第三维时间。好处是squeeze之后每一步都是标准 N×2 矩阵直接喂给step_sat坏处是 N 大、K 大时内存吃得凶N50、K4000 大约 3 MB还能接受再大就应该改成只存抽样点。out.cfg cfg这一句别省——三个月后你看到一份.mat唯一能回答“这组参数是怎么来的”的就是它。4.2 参数扫描与结果归档扫描的写法取决于你要扫几维。单参数比如增益或步长用一层循环双参数步长 × 通信半径用两层结果存成矩阵超过三维就别硬扫了改随机采样。归档用带时间戳的文件名别用result.mat这种名字互相覆盖。hList [0.02 0.05 0.10 0.20 0.30]; res cell(numel(hList), 1); for ii 1:numel(hList) cfg baseCfg; % 基准配置只改一个字段 cfg.h hList(ii); res{ii} sim_consensus(cfg); end stamp datestr(now, yyyymmdd_HHMMSS); save(fullfile(runs, [sweep_h_ stamp .mat]), res, hList);baseCfg这个模式值得单独说一句先在脚本顶部构造一份完整配置循环里只覆写变化的字段。这样新增参数时只需改一处不会出现“某组参数忘了设、用了默认值”的隐蔽错误。我见过最典型的一次是 vmax 忘了传结构体里没有这个字段函数在cfg.vmax处直接报错——这还算好的报错至少让你知道如果代码里写了getfield(cfg,vmax,10)给了默认值那就变成静默跑完、结论全错。归档目录我会在脚本开头if ~exist(runs,dir), mkdir(runs); end别指望手动建。4.3 批量绘图一眼看出哪组参数崩了绘图的目的不是好看是让“哪组参数崩了”这件事在 0.5 秒内可判断。做法是把所有组的误差曲线画在同一张半对数图上收敛的是直线下降振荡的是锯齿发散的直接冲顶。figure(Position, [100 100 720 480]); hold on; for ii 1:numel(res) semilogy(res{ii}.err, LineWidth, 1.2, ... DisplayName, sprintf(h %.2f, hList(ii))); end set(gca, YScale, log); grid on; xlabel(迭代步 k); ylabel(编队误差); legend(Location, southwest);YScale设成 log 是这一步的核心。线性坐标下误差从 5 降到 0.01 和从 5 降到 0.5 看起来差不多看不出收敛速率差异对数坐标下斜率直接就是收敛速度三组参数谁快谁慢一目了然。另外DisplayName一定要写不然 legend 出来的是一串data1、data2回头对不上。如果某条曲线在中段出现平台然后再次下降那通常是饱和起了作用前期被 vmax 卡住误差下降变慢进入线性区后恢复指数收敛。这个“拐点”位置本身是有信息的它告诉你限幅是不是给紧了。5. 排查收敛曲线不对时先查这五件事排查顺序很重要。我的习惯是先看矩阵再看数值最后才怀疑协议——因为协议是最不容易写错的部分而矩阵和步长是最容易在复制粘贴中出问题的部分。5.1 矩阵与拓扑层面的三个坑现象一曲线看着在收敛但每个智能体停在不同的值上。原因是图不连通或者本该无向的图被写成了有向且非强连通。无向连通图的一致性值是全体初始状态的算术平均一旦分裂成两个连通分量就变成两个各自的平均值。解决办法是跑一遍连通性检查numel(unique(conncomp(graph(A)))) 1不连通就回到边表找漏掉的边。现象二前期正常十几步之后突然发散且越跑越陡。步长超过稳定边界。查法很直接h和2/max(real(eig(L)))比一下超了就减半。我踩过一次更隐蔽的版本——代码里 H 是大写的局部变量、h 是小写的全局参数循环里用的是 H检查时改的是 h改了半天没反应血泪经验是变量名别玩这种字母大小写的玄学。现象三改了 N 之后报维度错误或者不报错但队形完全乱掉。硬编码索引没跟着改。典型位置是 A 的边表、P 的行数、x0 的长度以及eye(N)里那个 N。解决方式是全脚本搜一遍数字常量把与规模相关的都换成size(...)或者从 cfg 里取。5.2 数值与工程层面的两个坑现象四曲线里出现 NaN从某一步开始整条线消失。零除或溢出。最可能是初始位置重合导致距离为零或者限幅里vmax / n的 n 是 0。解决办法是给所有分母加max(x, eps)同时在循环里加一句assert(all(isfinite(u(:))), step %d 出现 NaN, k)——早报错比跑完再看图省时间得多。现象五脚本跑和 Simulink 跑结果对不上形状相似但数值差一截。求解器步长和采样步长不一致。Simulink 默认变步长求解器在饱和切换点会自动缩步脚本里是固定步长两者必然有差异。办法是把 Simulink 里的求解器改成固定步长步长填成和脚本一致再对比。如果还差就查模块里的采样时间有没有设成继承-1继承来的采样时间经常和你想的不是一个值。6. 进阶替换点把周期通信换成事件触发附一份验收清单周期通信的假设是“每个采样点所有智能体都要和邻居交换一次信息”在仿真里不要钱在实机上就是通信带宽和功耗。事件触发一致性把这条假设松掉只在状态偏差超过阈值时才通信其余时刻用保持器维持上一次的值。代码上的改动量小得出奇大约十几行但收益在规模大、通信受限的场景里非常明显。6.1 事件触发一致性的替换点核心是把“当前状态”换成“上次广播的状态”再维护一个触发误差。触发条件通常是 ‖e_i‖ σ‖z_i‖e_i 是当前状态与上次广播状态之差z_i 是邻居状态的加权和σ 是触发阈值越小越接近周期通信。function [X, nTrig] et_consensus(X0, L, h, K, sigma) N size(X0, 1); X zeros([size(X0), K1]); X(:,:,1) X0; Xhat X0; % 上次广播出去的状态 nTrig zeros(N, 1); % 每个智能体的触发次数 for k 1:K Xc X(:,:,k); z -L * Xhat; % 用广播状态算控制量 e Xc - Xhat; % 触发误差 % 逐智能体判断是否触发触发者更新广播值 need vecnorm(e, 2, 2) sigma * vecnorm(z, 2, 2); Xhat(need, :) Xc(need, :); nTrig nTrig need; X(:,:,k1) Xc h * z; end endz -L * Xhat用的是Xhat而不是Xc这是整个替换的关键——如果这里写成Xc那事件触发形同虚设你还是每步都在用最新状态。need是一个 N×1 的逻辑向量用逻辑索引整行更新比 for 循环快也比循环不容易写错下标。nTrig返回每个智能体的触发次数这个数字很能说明问题如果某个智能体的触发次数接近 K说明它一直在被邻居甩开问题在拓扑不在阈值。调参上有两条经验。sigma从 0.1 开始试太大会出现“长时间不通信、队形明显滞后”的现象太小则退化成周期通信、省不下什么。另外事件触发会引入额外的收敛时间同样 T 下误差可能比周期通信大一个量级比较两者时要在相同触发次数下比而不是相同时间内比——这一点在写结论时经常被忽略。6.2 一套我自己在用的验收清单代码改完之后跑一遍下面这几项比反复看曲线靠谱。这套清单是从多次“明明仿真收敛、换组参数就崩”的经历里攒出来的。检查项判据不通过时先看哪邻接矩阵对称性isequal(A, A)为真边表是否只写了单向拉普拉斯行和每行绝对值 1e-12度矩阵是否按行求和代数连通度 λ₂明显大于 0 1e-3图是否连通、有无孤立点步长裕度h 0.1 × 2/λ_max拓扑切换时取所有模式的公共上限有限性全程无 NaN / Inf分母是否加了 eps收敛值无向连通时等于初始均值是否有向图但按无向算饱和触顶率大部分时间不贴 vmaxvmax 给紧了或拓扑把某个体甩开触发占比事件触发远小于 1且各体差异不大sigma 是否过小、某体是否孤立表格里最容易被跳过的是“饱和触顶率”但它往往是最早暴露拓扑问题的指标。我一般会在仿真函数里顺手记一条satRatio mean(n vmax)跑完看这个数超过 0.5 就说明整段时间都在限幅区控制器根本没按线性协议工作这时候去调增益全是白费。还有一处建议把 λ₂、h 上限、饱和触顶率、末态误差这四个数打在一行日志里输出。曲线可以骗人尤其是只看末态的时候这四个标量不太会。我自从被一次“图不连通但曲线看起来收敛”坑掉整整两天之后每次改完拓扑做的第一件事都是先跑那三条断言再画图。断言不过就不看曲线因为那时候的曲线只会给你错误的信息。希望帮到你。本文还有配套的精品资源点击获取 多智能体MATLAB代码拆解：一致性协议、编队控制与事件触发
简介面向多智能体系统初学者与MATLAB仿真实践者的示例代码包围绕MAS中代理建模、环境感知、消息交互、协调决策与仿真验证等核心环节展开适合希望用脚本方式理解多代理协作机制、搭建小型实验原型的读者。压缩包共4个文件以3个m脚本和1个txt许可说明为主体量约4KB轻量便于查阅与二次修改脚本承担代理状态更新、位置计算及服务端协调等示例功能许可文件用于说明使用授权。已有1806人学习下载可作为课程演示、课题预研或自学练手的参考。内容侧重多智能体服务器如何注册管理代理、分配任务、维护通信通道并定义协同规则同时保留监控评估的扩展空间读者可据此观察服务端创建与跟踪代理实例、按能力派发任务的过程并把代理、环境、交互、协调与决策等概念落到可运行的MATLAB代码骨架上。1. 多智能体 MATLAB 代码拆包先跑通 5 个智能体的编队收敛再谈改协议拿到一份多智能体 MATLAB 代码多数人的第一反应是直接 run 主脚本看那几条曲线收不收敛——收敛就收工散了就怀疑算法写错。我的经验恰好相反MAS 仿真翻车八成不是协议本身的问题而是通信拓扑、采样步长、初始状态这三处没对齐。这套代码真正的价值是把一致性协议、编队偏移、分布式控制器和仿真外壳拆成了互相独立、可单独替换的模块改一处不牵连其余。它适合正在做编队、蜂拥、分布式优化验证的人也适合想把 Simulink 里那堆手工连线换成脚本化流程、方便批量跑参数的场景。2. 一致性协议落到矩阵上拉普拉斯矩阵的三种写法与增益边界多智能体系统在代码层面的全部秘密几乎都压缩在一个矩阵里拉普拉斯矩阵 L。协议写成微分方程是 ẋ -Lx写成代码就是一次矩阵乘法。所以这一章不急着调控制器先把 L 造对——L 错了后面所有参数整定都是白费功夫。2.1 邻接矩阵的三种建法邻接矩阵 A 描述“谁能收到谁的信息”是整个 MAS 的地基。小规模固定拓扑手写边表最直观中等规模按几何距离生成通信半径内的智能体连边大规模或有现成拓扑数据时从稀疏矩阵导入。三种写法各有各的翻车点手写最容易犯的错是只赋值单向。% 建法一手写边表适合 N 10 的固定拓扑 N 5; edges [1 2; 2 3; 3 4; 4 5; 5 1; 1 3]; % 每行一条无向边 A zeros(N); for k 1:size(edges, 1) i edges(k, 1); j edges(k, 2); A(i, j) 1; A(j, i) 1; % 无向图必须成对赋值漏一行就变成有向图 end % 建法二按通信半径生成适合随机撒点场景 function A radius_graph(P, R) N size(P, 1); D squareform(pdist(P)); % 两两欧氏距离 A double(D R D 0); % 排除自身D0 是关键 end % 建法三从稀疏表恢复适合大规模或外部拓扑 % A sparse(rowIdx, colIdx, 1, N, N); A full(A | A);第一种写法里那句成对赋值不是啰嗦是我踩过的坑早期图省事只写 A(i,j)1结果 L 不对称特征值出现负实部曲线前半段看着挺正常后半段直接飞出去。第二种写法的D 0同理——不加这一条对角线会被判成“距离 0 以内”自我连边L 的对角元凭空多出 N收敛速度全乱。第三种写法补一句A | A做对称化防止上游数据只给了半边。提示建完 A 之后立刻断言isequal(A, A)。这一步成本几乎为零能挡掉后面一大半莫名其妙的发散。2.2 从拉普拉斯矩阵到状态方程有了 A拉普拉斯矩阵一行就能算出来L D - AD 是度矩阵。这个式子的物理含义是“我拿自己和所有邻居的差做加权和”所以 L 的每一行和恒为零全 1 向量永远是它的零特征值对应向量——这正是所有智能体最终停在同一值的数学来源。收敛速度由第二小特征值 λ₂ 决定代数连通度越高收敛越快。D diag(sum(A, 2)); % 度矩阵按行求和 L D - A; % 拉普拉斯矩阵 % 三条自检建议写进每个脚本开头 assert(isequal(A, A), 邻接矩阵不对称); assert(all(abs(sum(L, 2)) 1e-12), 拉普拉斯行和不为零); lam eig(L); lam2 sort(real(lam)); fprintf(代数连通度 lambda_2 %.4f\n, lam2(min(2, numel(lam2))));三条自检对应三类错误不对称会让 L 出现复特征值行和不为零说明 D 求错了维度λ₂ 接近零说明图不连通——不连通时一致性协议不是发散而是分裂成几个各自收敛的簇这种“部分收敛”最迷惑人曲线看着像要收敛又永远差一点。连续时间协议 ẋ -Lx 在代码里通常用前向欧拉离散x(k1) (I - h·L)x(k)稳定条件是 h 2/λ_max(L)。这个边界不是理论摆设我见过 h 设成 0.5、λ_max 大约 4 的情形正好压在边上前几步正常二十步后开始振荡放大。2.3 五智能体最小可跑脚本把上面的东西拼起来就是一段不到 20 行、但能完整验证协议的最小脚本。跑它的意义在于在引入编队、饱和、时延之前先确认你的拓扑和步长是干净的。N 5; A zeros(N); edges [1 2; 2 3; 3 4; 4 5; 5 1; 1 3]; for k 1:size(edges, 1) A(edges(k,1), edges(k,2)) 1; A(edges(k,2), edges(k,1)) 1; end L diag(sum(A,2)) - A; x0 [1 3 2 5 4]; % 初始状态故意打散 h 0.05; % 采样步长 T 20; % 仿真时长 K round(T / h); X zeros(N, K1); X(:,1) x0; for k 1:K X(:,k1) (eye(N) - h*L) * X(:,k); end xf mean(x0); % 无向连通图的收敛值即初始均值 err max(abs(X(:,end) - xf)); fprintf(末态最大偏差 %.3e理论 lambda_2 %.4f\n, err, ... sort(eig(L))(2));逻辑是一条直线加权平均 → 矩阵乘 → 迭代。参数上有三处会直接影响结果。h必须小于 2/λ_max实用做法是取 0.1/λ_max宁可多跑几步也别贴边界。T要留足经验值是 10/λ₂ 以上λ₂ 越小需要的时间越长图越“细长”越慢。x0最好手动打散别用rand直接生成——调试阶段你需要能心算出期望均值随机初始值会让你失去这个参照。收敛判据也别只看最后一步。我在脚本里习惯再算一条errHist max(abs(X - mean(x0)), [], 1)把误差历史画出来。真收敛是指数下降的直线半对数坐标下如果误差在某个量级上压不下去多半是 λ₂ 太小或步长偏大导致的残余振荡光看末态一个数看不出来。3. 编队控制的可调旋钮偏移量、饱和限幅与拓扑切换一致性只要求大家“停在一起”编队要求“停成某个形状”。从一致性到编队代码上只多了一个偏移量矩阵但工程上的复杂度翻了一倍——因为偏移量引入了几何约束几何约束又逼着你处理饱和、避碰和拓扑变化。3.1 把方针翻译成偏移量矩阵编队协议的标准写法是 u_i -Σ a_ij[(x_i - p_i) - (x_j - p_j)]其中 p_i 是第 i 个智能体的期望相对位置。写成矩阵就是 u -L(X - P)。这里最容易搞混的是 P 的语义它给的是相对偏移形心在原点还是绝对目标位置代码里必须统一否则整个队形会整体漂移到画面外。% 二维平面上的五个位置pentagon 队形形心在原点 P [ 0.0 1.0; 0.95 0.31; 0.59 -0.81; -0.59 -0.81; -0.95 0.31]; % N×2第 i 行是第 i 个智能体的期望偏移 function Xn formation_step(X, P, L, h) % X: N×2 当前状态P: N×2 期望偏移L: N×N 拉普拉斯 E X - P; % 先减偏移再走一致性 Xn X - h * (L * E); endE X - P这一行的顺序不能反。写成L*X - P是另一个东西那是带绝对参考的跟踪写反了队形会一边收敛一边整体平移看起来像“慢悠悠地飘走”。另外 P 的行数必须等于 N我在批量换 N 的脚本里吃过亏P 是硬编码的 5 行N 改成 7 之后 MATLAB 不报错而是自动做隐式扩展队形直接算成一团乱麻——这种不报错的 bug 最难查。注意换 N 的时候A、L、P、x0 四处的维度必须一起改。建议在函数入口加assert(size(P,1) size(L,1))一秒定位。3.2 饱和限幅与避碰项两个位置别写反真实执行器有速度上限仿真里不加饱和编队初期会出现几十米每秒的瞬时速度看着收敛很漂亮换成实机直接失控。限幅的写法有讲究要等比缩放不能逐分量裁剪否则方向会变队形走形。function Xn step_sat(X, P, L, h, vmax) u -L * (X - P); % 原始速度指令N×2 n sqrt(sum(u.^2, 2)); % 每个智能体的指令模长 scale min(1, vmax ./ max(n, 1e-9)); % 超过上限才缩放 u u .* scale; % 整体缩放保持方向 Xn X h * u; endmax(n, 1e-9)是防零除的后悔药初始位置如果和期望位置完全重合指令为零0/0 会给你 NaN之后整条曲线全是 NaN 且不报错。scale用min(1, ...)保证只往下压不放上去这一点在调试时很关键——如果发现某个智能体始终贴着 vmax 跑说明它被拓扑甩开了问题在图不在限幅。避碰项和限幅是两个独立的东西别塞进同一个 if 里。常见做法是在 u 后面叠加一个斥力项距离小于安全半径 d_safe 时沿相对位置方向加一个反比项。斥力增益不能给大否则会和一致性项打架出现两个智能体高频抖动、位置几乎不动的现象——这种抖动在曲线上看是毛刺在相图上看是一个极限环。3.3 采样步长与拓扑切换的稳定边界固定拓扑下 h 2/λ_max 就够了一旦拓扑会切换通信中断、邻居变化稳定条件变成对所有可能拓扑的并集都成立也就是 h 2/max(λ_max over all modes)。这条规则在实际代码里经常被忽略因为单跑某个拓扑是稳的一切换就炸。% 拓扑集合正常、断边、链路退化三种模式 Amodes {A1, A2, A3}; hmax inf; for m 1:numel(Amodes) Lm diag(sum(Amodes{m},2)) - Amodes{m}; hmax min(hmax, 2 / max(real(eig(Lm)))); end fprintf(拓扑切换下的步长上限 %.4f建议取 %.4f\n, hmax, 0.1*hmax);这段循环算的是“最保守的公共步长”。实用取值我一般按 0.1 倍上限代价是仿真步数变多但换来的是切换瞬间不会出现数值爆点。如果确实需要大 step 提速那就得改成严格离散协议用行随机矩阵而不是 I - hL牺牲一点物理直觉换取无条件稳定。切换逻辑本身还有一个隐藏坑切换时刻落在两个采样点中间时状态更新用的是旧拓扑还是新拓扑仿真里通常按采样点对齐切换只在 k 为整数时生效。如果你的切换由事件触发比如距离阈值一定要把触发时刻量化到采样网格否则同一份代码跑两次结果不同排查起来毫无头绪。4. 仿真工程化参数扫描、数据归档与批量绘图单次仿真跑通只是及格线。真要做对比、写报告、给导师或同事看结论你需要的是“一次配置、批量跑完、结果能追溯”。这一章把脚本改造成可复用的三层结构仿真函数、扫描脚本、绘图脚本。4.1 把单次仿真封成一个函数函数化的第一原则是输入只有参数输出只有数据中间不画图、不打印、不依赖工作区变量。这样它才能被parfor调用也才能被后来的自己信任。function out sim_consensus(cfg) % cfg 字段N, edges, x0, h, T, P, vmax N cfg.N; A zeros(N); for k 1:size(cfg.edges, 1) i cfg.edges(k,1); j cfg.edges(k,2); A(i,j) 1; A(j,i) 1; end L diag(sum(A,2)) - A; K round(cfg.T / cfg.h); X zeros(N, 2, K1); X(:,:,1) cfg.x0; for k 1:K X(:,:,k1) step_sat(squeeze(X(:,:,k)), cfg.P, L, cfg.h, cfg.vmax); end out.X X; out.L L; out.lam2 sort(eig(L)); out.cfg cfg; % 把配置原样存回去方便追溯 out.err squeeze(max(vecnorm(X - cfg.P, 2, 2), [], 1)); end三维数组X(N,2,K1)是这里唯一需要适应的地方第一维智能体、第二维状态分量、第三维时间。好处是squeeze之后每一步都是标准 N×2 矩阵直接喂给step_sat坏处是 N 大、K 大时内存吃得凶N50、K4000 大约 3 MB还能接受再大就应该改成只存抽样点。out.cfg cfg这一句别省——三个月后你看到一份.mat唯一能回答“这组参数是怎么来的”的就是它。4.2 参数扫描与结果归档扫描的写法取决于你要扫几维。单参数比如增益或步长用一层循环双参数步长 × 通信半径用两层结果存成矩阵超过三维就别硬扫了改随机采样。归档用带时间戳的文件名别用result.mat这种名字互相覆盖。hList [0.02 0.05 0.10 0.20 0.30]; res cell(numel(hList), 1); for ii 1:numel(hList) cfg baseCfg; % 基准配置只改一个字段 cfg.h hList(ii); res{ii} sim_consensus(cfg); end stamp datestr(now, yyyymmdd_HHMMSS); save(fullfile(runs, [sweep_h_ stamp .mat]), res, hList);baseCfg这个模式值得单独说一句先在脚本顶部构造一份完整配置循环里只覆写变化的字段。这样新增参数时只需改一处不会出现“某组参数忘了设、用了默认值”的隐蔽错误。我见过最典型的一次是 vmax 忘了传结构体里没有这个字段函数在cfg.vmax处直接报错——这还算好的报错至少让你知道如果代码里写了getfield(cfg,vmax,10)给了默认值那就变成静默跑完、结论全错。归档目录我会在脚本开头if ~exist(runs,dir), mkdir(runs); end别指望手动建。4.3 批量绘图一眼看出哪组参数崩了绘图的目的不是好看是让“哪组参数崩了”这件事在 0.5 秒内可判断。做法是把所有组的误差曲线画在同一张半对数图上收敛的是直线下降振荡的是锯齿发散的直接冲顶。figure(Position, [100 100 720 480]); hold on; for ii 1:numel(res) semilogy(res{ii}.err, LineWidth, 1.2, ... DisplayName, sprintf(h %.2f, hList(ii))); end set(gca, YScale, log); grid on; xlabel(迭代步 k); ylabel(编队误差); legend(Location, southwest);YScale设成 log 是这一步的核心。线性坐标下误差从 5 降到 0.01 和从 5 降到 0.5 看起来差不多看不出收敛速率差异对数坐标下斜率直接就是收敛速度三组参数谁快谁慢一目了然。另外DisplayName一定要写不然 legend 出来的是一串data1、data2回头对不上。如果某条曲线在中段出现平台然后再次下降那通常是饱和起了作用前期被 vmax 卡住误差下降变慢进入线性区后恢复指数收敛。这个“拐点”位置本身是有信息的它告诉你限幅是不是给紧了。5. 排查收敛曲线不对时先查这五件事排查顺序很重要。我的习惯是先看矩阵再看数值最后才怀疑协议——因为协议是最不容易写错的部分而矩阵和步长是最容易在复制粘贴中出问题的部分。5.1 矩阵与拓扑层面的三个坑现象一曲线看着在收敛但每个智能体停在不同的值上。原因是图不连通或者本该无向的图被写成了有向且非强连通。无向连通图的一致性值是全体初始状态的算术平均一旦分裂成两个连通分量就变成两个各自的平均值。解决办法是跑一遍连通性检查numel(unique(conncomp(graph(A)))) 1不连通就回到边表找漏掉的边。现象二前期正常十几步之后突然发散且越跑越陡。步长超过稳定边界。查法很直接h和2/max(real(eig(L)))比一下超了就减半。我踩过一次更隐蔽的版本——代码里 H 是大写的局部变量、h 是小写的全局参数循环里用的是 H检查时改的是 h改了半天没反应血泪经验是变量名别玩这种字母大小写的玄学。现象三改了 N 之后报维度错误或者不报错但队形完全乱掉。硬编码索引没跟着改。典型位置是 A 的边表、P 的行数、x0 的长度以及eye(N)里那个 N。解决方式是全脚本搜一遍数字常量把与规模相关的都换成size(...)或者从 cfg 里取。5.2 数值与工程层面的两个坑现象四曲线里出现 NaN从某一步开始整条线消失。零除或溢出。最可能是初始位置重合导致距离为零或者限幅里vmax / n的 n 是 0。解决办法是给所有分母加max(x, eps)同时在循环里加一句assert(all(isfinite(u(:))), step %d 出现 NaN, k)——早报错比跑完再看图省时间得多。现象五脚本跑和 Simulink 跑结果对不上形状相似但数值差一截。求解器步长和采样步长不一致。Simulink 默认变步长求解器在饱和切换点会自动缩步脚本里是固定步长两者必然有差异。办法是把 Simulink 里的求解器改成固定步长步长填成和脚本一致再对比。如果还差就查模块里的采样时间有没有设成继承-1继承来的采样时间经常和你想的不是一个值。6. 进阶替换点把周期通信换成事件触发附一份验收清单周期通信的假设是“每个采样点所有智能体都要和邻居交换一次信息”在仿真里不要钱在实机上就是通信带宽和功耗。事件触发一致性把这条假设松掉只在状态偏差超过阈值时才通信其余时刻用保持器维持上一次的值。代码上的改动量小得出奇大约十几行但收益在规模大、通信受限的场景里非常明显。6.1 事件触发一致性的替换点核心是把“当前状态”换成“上次广播的状态”再维护一个触发误差。触发条件通常是 ‖e_i‖ σ‖z_i‖e_i 是当前状态与上次广播状态之差z_i 是邻居状态的加权和σ 是触发阈值越小越接近周期通信。function [X, nTrig] et_consensus(X0, L, h, K, sigma) N size(X0, 1); X zeros([size(X0), K1]); X(:,:,1) X0; Xhat X0; % 上次广播出去的状态 nTrig zeros(N, 1); % 每个智能体的触发次数 for k 1:K Xc X(:,:,k); z -L * Xhat; % 用广播状态算控制量 e Xc - Xhat; % 触发误差 % 逐智能体判断是否触发触发者更新广播值 need vecnorm(e, 2, 2) sigma * vecnorm(z, 2, 2); Xhat(need, :) Xc(need, :); nTrig nTrig need; X(:,:,k1) Xc h * z; end endz -L * Xhat用的是Xhat而不是Xc这是整个替换的关键——如果这里写成Xc那事件触发形同虚设你还是每步都在用最新状态。need是一个 N×1 的逻辑向量用逻辑索引整行更新比 for 循环快也比循环不容易写错下标。nTrig返回每个智能体的触发次数这个数字很能说明问题如果某个智能体的触发次数接近 K说明它一直在被邻居甩开问题在拓扑不在阈值。调参上有两条经验。sigma从 0.1 开始试太大会出现“长时间不通信、队形明显滞后”的现象太小则退化成周期通信、省不下什么。另外事件触发会引入额外的收敛时间同样 T 下误差可能比周期通信大一个量级比较两者时要在相同触发次数下比而不是相同时间内比——这一点在写结论时经常被忽略。6.2 一套我自己在用的验收清单代码改完之后跑一遍下面这几项比反复看曲线靠谱。这套清单是从多次“明明仿真收敛、换组参数就崩”的经历里攒出来的。检查项判据不通过时先看哪邻接矩阵对称性isequal(A, A)为真边表是否只写了单向拉普拉斯行和每行绝对值 1e-12度矩阵是否按行求和代数连通度 λ₂明显大于 0 1e-3图是否连通、有无孤立点步长裕度h 0.1 × 2/λ_max拓扑切换时取所有模式的公共上限有限性全程无 NaN / Inf分母是否加了 eps收敛值无向连通时等于初始均值是否有向图但按无向算饱和触顶率大部分时间不贴 vmaxvmax 给紧了或拓扑把某个体甩开触发占比事件触发远小于 1且各体差异不大sigma 是否过小、某体是否孤立表格里最容易被跳过的是“饱和触顶率”但它往往是最早暴露拓扑问题的指标。我一般会在仿真函数里顺手记一条satRatio mean(n vmax)跑完看这个数超过 0.5 就说明整段时间都在限幅区控制器根本没按线性协议工作这时候去调增益全是白费。还有一处建议把 λ₂、h 上限、饱和触顶率、末态误差这四个数打在一行日志里输出。曲线可以骗人尤其是只看末态的时候这四个标量不太会。我自从被一次“图不连通但曲线看起来收敛”坑掉整整两天之后每次改完拓扑做的第一件事都是先跑那三条断言再画图。断言不过就不看曲线因为那时候的曲线只会给你错误的信息。希望帮到你。本文还有配套的精品资源点击获取