简介面向多智能体网络一致性研究的MATLAB代码包适用于学习分布式控制、一致性协议及群体行为仿真的高校学生、科研人员和工程开发者也可用于课程设计、毕业设计或初步科研探索。代码围绕状态更新、权重设计、网络切换等核心问题展开包含智能体动力学模型、典型一致性算法与结果可视化脚本可帮助验证平均共识、拉普拉斯共识等经典协议在不同初始条件和拓扑下的收敛表现为分析稳态误差与收敛速度提供直接依据。压缩包共五个文件以四个MATLAB脚本为主采用清晰的模块化写法分别负责模型定义、协议实现与仿真分析另含一个Simulink模型支持图形化搭建并拖拽调节参数直观观察多智能体动态演化过程。资源整体仅17KB轻量便携已有213人学习。通过学习这套代码读者能快速掌握一致性算法的实现思路同时获得可扩展的代码框架便于引入异构个体、动态拓扑等更贴近实际的研究因素是一份颇具性价比的学习材料。1. 多智能体网络一致性这套MATLAB代码包能帮你省掉两周搭框架的时间拿到压缩包的时候我先翻了一遍文件列表——square_3D.m、switching.m、weight.m、model.m、untitled.mdl五个文件正好覆盖了多智能体网络一致性研究里最常见的主线静态拓扑下的共识协议、动态切换网络下的收敛、以及Simulink可视化建模。多智能体网络一致性说白了就是让一群各自为政的智能体通过局部信息交互把状态对齐这套代码解决的是从零搭仿真环境这个最耗时的环节。适合正在做分布式控制、编队、共识算法方向的研究生和工程师尤其是那些对着Laplacian矩阵和一致性公式不知道怎么落成可运行代码的人。2. 拆包看结构五个文件不是散的是一条完整的仿真链路2.1 从文件命名反推代码模块model.m与weight.m的职责边界一个多智能体一致性仿真工程无论代码怎么组织都绕不开三个层次模型层负责描述每个智能体怎么运动协议层决定状态怎么更新配置层管理网络拓扑和仿真参数。对照这份资源里的文件model.m属于模型层weight.m和switching.m跨在协议层与配置层之间untitled.mdl和square_3D.m则承担入口角色——前者是Simulink建模入口后者是把前四个模块串成完整仿真的主脚本。我拆这类代码包的习惯是先不看实现细节而是把文件之间的数据依赖画出来。weight.m输出的是邻接矩阵或图拉普拉斯矩阵这个矩阵要喂给model.mswitching.m做的事情是按时间改变这个矩阵square_3D.m在最外层做循环调度把状态向量、矩阵、时间轴搅在一起untitled.mdl则用Simulink的图形化接线方式做了同样的事。搞清楚这条链路你再打开任何一个文件都不会迷路。这里值得先强调一个非常基础的认知一致性协议里最核心的矩阵是图拉普拉斯矩阵L不是邻接矩阵A。A只描述了谁和谁有通信边L才真正进入控制律。很多初次上手的人把A当成控制矩阵直接用结果收敛方向完全不对。这份代码里既然单独设了weight.m说明作者刻意把图构造和动力学解耦了这个设计习惯值得沿用。2.2 model.m一阶积分器模型怎么落到MATLAB函数model.m最可能实现的是连续时间一阶积分器模型这也是多智能体一致性研究里最经典的切入模型。每个智能体的状态变化率等于控制输入而控制输入取邻居状态差的加权和。用公式写就是 dx_i/dt -Σ l_ij * x_j矩阵形式就是 dx -Lx。这段代码的核心逻辑是function dx model(t, x, L) % 智能体动力学一阶积分器 一致性协议 % t : 时间变量保留以兼容 ode45 的函数签名 % x : n 维列向量存放 n 个智能体的标量状态 % L : n x n 图拉普拉斯矩阵由 weight.m 生成 dx -L * x; end这个函数的写法有几个细节值得说。第一t参数虽然函数体里没用到但必须保留因为后面用ode45调用时MATLAB要求被积函数签名必须是 f(t, x) 形式去掉t会直接报错。第二dx -L * x 这行代码里的负号来自Laplacian算子的谱性质——L的特征值非负且有一个零特征值对应全连通时的一致状态负号保证了非零特征值对应的模态随时间指数衰减这正是收敛的来源。实际仿真调用时整个流程是把weight.m生成的L矩阵传进来初始状态随机给一组值用ode45跑一段时间最后看所有智能体的状态曲线是否汇聚到同一个值。常见做法是这样串起来n 5; % 智能体数量 x0 randn(n, 1); % 随机初值模拟任意初始状态 [A, L] weight(n, ring); % 环形拓扑 [t, x] ode45((t, x) model(t, x, L), [0 20], x0); plot(t, x, LineWidth, 1.5); xlabel(时间/s); ylabel(智能体状态); title(环形拓扑下的一致性收敛过程);到这里你能直观看到一个现象不管初值多乱五条曲线最终会挤到同一水平线上。这个水平线的位置取决于L零特征值对应的左特征向量——在无向连通图下最终值就是所有初值的算术平均。很多论文里说的平均一致性指的就是这件事。2.3 weight.m邻接矩阵与图拉普拉斯矩阵的工程化写法weight.m承担的工作是根据拓扑类型生成矩阵这是整个仿真里最容易出隐蔽错误的地方。常见的拓扑无非三类环形、全连接、随机图。实现方式如下function [A, L] weight(n, topology) % n 智能体数量 % topology 拓扑类型ring / complete / random switch topology case ring A zeros(n); for i 1:n j mod(i, n) 1; % 尾节点接回头节点 A(i, j) 1; A(j, i) 1; end case complete A ones(n) - eye(n); % 全连接除自环外全为1 case random rng(42); % 固定随机种子保证结果可复现 A rand(n) 0.6; A (A A) 0; % 对称化保证网络是无向图 A A - diag(diag(A)); % 去掉自环 otherwise error(weight.m: 未知拓扑类型 %s, topology); end D diag(sum(A, 2)); L D - A; % 组合拉普拉斯矩阵 L D - A end里面几个参数需要解释。rng(42) 是我特别建议保留的一行随机图如果不固定随机种子每次运行生成的拓扑都不一样那么你比较不同算法性能时就没法排除网络结构这个变量结果全成了玄学。rand(n) 0.6 里的0.6是连边概率调大得到稠密图调小得到稀疏图但要注意概率低于0.2左右时图很容易不连通一致性收敛就无从谈起。A (A A) 0 这种对称化写法比用 upper transpose 更简洁生成的矩阵天然满足无向图要求。还有一个隐藏细节D diag(sum(A, 2)) 是对每行求和由于A是对称的行和等于列和直接这样写没问题。如果你改成列和 diag(sum(A, 1))在无向图下结果一样但万一A不小心构造得不对称行和与列和就会不同L矩阵也不再是标准的组合拉普拉斯。所以我在代码注释里强调组合拉普拉斯矩阵是有原因的——还有一类叫规范化拉普拉斯矩阵 L_norm D^(-1/2) * L * D^(-1/2)两者的特征值分布不同收敛速度分析与代码写法都不同。这份资源如果以后你想做更深入的理论分析可能会需要补一个函数。3. 一致性协议落地从Laplacian共识到控制律的参数配置3.1 连续时间一致性协议ode45求解与状态更新写法连续时间一致性协议的标准形式是状态导数等于负的Laplacian乘以状态上一章的model.m已经把这个协议写成了MATLAB函数。不过实际仿真时很多人会纠结一个问题ode45的默认精度够不够收敛曲线看起来挺光滑但数值误差会导致最终值偏离理论平均值吗答案是不会太明显但需要注意。ode45用的是自适应步长的Runge-Kutta方法相对误差默认1e-3绝对误差默认1e-6。对于一致性这类线性系统这个精度完全够用。真正需要担心的是时间窗口设置——如果仿真时长太短曲线还没收完你会误判算法不收敛如果太长白白消耗算力。我的做法是先用L的第二小特征值来估算时间常数。对连续时间一致性系统收敛时间常数与第二小特征值 λ₂ 直接相关理论收敛速度由 λ₂ 决定工程估算时取 T_end 10 / λ₂ 左右就能看到明显的收敛尾部。举一个具体配置例子[A, L] weight(6, complete); lambda2 eigs(L, 6, smallestreal); % 取最小的非零特征值 T_end 10 / abs(lambda2(1));这里 eigs 的第三个参数 smallestreal 用来取最小实特征值但要注意 eigs 返回的是按模或按实部排序的多个特征值最小的那个是零特征值你需要取倒数第二个或者直接对全谱做一次 sort。这属于高频翻车点后面避坑章节会细说。除此之外如果系统里有通信时延连续时间模型就不能再用简单的 dx -Lx而要变成带时滞的微分方程MATLAB里要用 dde23 而不是 ode45 来求解。这份资源里没有涉及时延但如果你把代码往实际系统方向扩展迟早会遇到这一步。3.2 离散时间迭代算法步长与收敛性的关系很多实际系统——尤其是数字控制器实现的分布式算法——用的是离散时间迭代而不是连续微分方程。离散一致性算法的更新公式为 x(k1) W * x(k)其中 W 必须是一个双随机矩阵即行和、列和都为1。由 Laplacian 构造 W 的常见做法是 W I - εL其中 ε 是步长。下面是一个可以直接跑通的离散迭代实现function x_hist discrete_consensus(x0, L, epsilon, steps) % x0 n 维初值列向量 % L n x n 图拉普拉斯矩阵 % epsilon 步长必须小于 2/max(eig(L)) % steps 迭代步数 n length(x0); W eye(n) - epsilon * L; % 双随机权重矩阵 x_hist zeros(n, steps 1); x_hist(:, 1) x0; for k 1:steps x_hist(:, k 1) W * x_hist(:, k); end end步长 ε 的选择是有明确理论边界的。离散系统稳定要求 W 的谱半径小于等于1这就推出 ε 必须小于 2/λ_max(L)。实际工程里ε 取 1/λ_max 附近收敛最快取太小收敛很慢取太大会直接发散成一堆振荡值。我的习惯是先算一次特征值再定步长[~, L] weight(10, random); lambda_max eigs(L, 1, largestreal); % 最大特征值 epsilon 0.9 / lambda_max; % 留一点稳定裕度 x_hist discrete_consensus(randn(10, 1), L, epsilon, 500);这也是为什么很多论文里步长喜欢写成 0.9/λ_max 或者 1/(Δ λ_max) 的形式——不是硬凑的是为了在收敛速度和稳定域之间留出安全距离。如果你把 epsilon 设成 1.0/λ_max大概率也能收敛但遇到拓扑切换的瞬间网络拉普拉斯特征值突变就可能翻车。3.3 参数配置清单耦合强度、初值范围与仿真时长仿真参数的配置直接决定了你能不能看到理想的一致性收敛图像。这里给出一份我实际调仿真时用的参数清单覆盖了最常见的几个旋钮参数建议范围影响典型错误智能体数量 n3 ~ 30数量越多图越稀疏时收敛越慢数量太少无法体现网络特性初值范围[-10, 10] 均匀随机不影响收敛性影响曲线观感初值差距太小收敛过程看不出来耦合强度连续系统隐含在L中可整体乘以增益 kk 越大收敛越快但受数值稳定性限制忘记调整 k以为L不能缩放步长 ε离散系统0.1/λ_max ~ 0.9/λ_max决定离散收敛速度与稳定性直接用 0.01 固定值脱离特征值仿真时长 T5/λ₂ ~ 15/λ₂太短截断收敛过程太长空跑对所有拓扑用同一时长随机种子固定一个值保证实验可复现每次运行生成不同拓扑耦合强度这一个参数值得单独说。连续时间一致性协议 dx -Lx 里没有显式的增益系数但你可以把协议改写成 dx -k Lxk 就是全局耦合增益。k 越大状态之间的拉扯力越强收敛越快。但注意ode45的积分器对刚性问题有天然短板k 超过一定阈值后数值误差会急剧增大曲线出现高频抖动。k 取多少需要按 λ₂ 和 λ_max 的比值来定我一般从 k1 开始观察收敛时间再倍增试。4. Simulink与切换拓扑动态网络的仿真细节4.1 untitled.mdlSimulink模型的结构与信号流向untitled.mdl是这份资源里的Simulink模型文件打开后你能看到一个典型的动态系统框图。这类模型最常见的搭建方式是用State-Space模块承载一致性系统的状态方程用Gain模块设置耦合增益用Sum模块做邻居状态差的叠加最后把状态信号引到Scope里观察收敛曲线。模型里的数据流是Constant模块或From Workspace模块提供初值状态空间模块输出状态反馈回路里接一个矩阵增益即 -L形成闭环。拿到手先别急着点Run先做三件事查看模型里有没有配置仿真时长、有没有把L矩阵正确写到增益参数里、Scope模块有没有接入状态信号。很多从别人电脑拷来的Simulink模型会带着调试时的残留配置仿真时长停在0.1s上一运行Scope里只有一条横线误以为代码坏了实际上只是时长没改。Simulink模型与脚本代码的差异在于模型里的矩阵增益是静态的——如果你直接在Gain模块里写死一个L矩阵那拓扑切换就无从谈起。所以切换到动态拓扑时要么用S-Function在每一个仿真步更新L要么干脆用MATLAB脚本循环驱动后者在调试时更直观。这就是为什么资源里同时存在switching.m和untitled.mdl——两者对应的是不同的使用场景。4.2 switching.m周期切换拓扑的实现与索引保护切换拓扑的一致性研究要回答的核心问题是如果智能体之间的通信结构随时间变化一致性还能不能达成答案与切换频率和网络并集连通性有关。switching.m里最标准的做法是预设一组拓扑按固定时间间隔在它们之间切换。关键代码如下function L switching(t, L_set, interval, mode) % t 当前仿真时间 % L_set n x n x k 三维数组每一页是一个拓扑的拉普拉斯矩阵 % interval 拓扑切换间隔秒 % mode cycle 循环切换 / hold 越界后保持不变 k size(L_set, 3); idx floor(t / interval) 1; if strcmp(mode, cycle) idx mod(idx - 1, k) 1; % 循环取模防止越界 else idx min(idx, k); % 停在最后一个拓扑 end L L_set(:, :, idx); end这里 idx 的计算是最容易写错的一行。floor(t/interval)1 在 t0 时返回1对应第一个拓扑tinterval 时返回2切到第二个拓扑。如果用 ceil 而不是 floort0 时就会直接跳到索引1但 tinterval 时会提前切换时序逻辑会乱掉。mod(idx-1, k)1 这种写法是为了让索引在 k 页之间循环——如果不加这层保护仿真时间一旦超过 interval*k索引就越界了MATLAB直接报 Index exceeds array bounds。使用上要把 L_set 预先构造成三维数组。生成方式一般是循环调用 weight.m 多次把每一次的 L 存进一页L_set zeros(n, n, 3); for i 1:3 [~, L] weight(n, random); L_set(:, :, i) L; end然后在整个仿真循环的每一步调用 switching 函数取当前的 L。注意在切换瞬间L 矩阵会突变如果用的是 ode45 连续求解这种突变的矩阵会让积分器在同一时刻反复调整步长仿真变慢是正常现象不要误判为死机。4.3 脚本驱动Simulink仿真的常见联动方式脚本和Simulink联动的标准做法是用 sim 命令把Simulink模型当作一个黑匣子来调用。这种模式的好处是你可以用脚本循环跑几十组参数Simulink负责算脚本负责收集结果。以 untitled.mdl 为例假设模型里已经有 To Workspace 模块把状态写入变量 tout 和 xout脚本侧这样写% 设置Simulink模型的仿真参数 set_param(untitled, StopTime, 20); % 把L矩阵传给模型里的增益模块 L_ring weight(5, ring); set_param(untitled/Gain, Gain, mat2str(L_ring)); % 运行仿真并取回数据 simOut sim(untitled); x simOut.get(xout);set_param 属于Simulink的API第一个参数是模块路径第二个是属性名。mat2str 把矩阵转成字符串因为Gain模块的增益参数是字符串格式而非数值矩阵格式直接传数组会报错。这个坑很隐蔽很多人卡在这儿。如果你的模型里Gain模块显示为矩阵形式直接用 set_param(.../Gain, Gain, mat2str(L)) 即可。使用 sim 命令驱动模型的另一层作用是可以批量做蒙特卡洛实验——不同初值、不同拓扑、不同增益组合各跑一遍把结果存到变量里统一分析。这比在Simulink界面里手工点Run改参数高效得多是后期调参的主力方式。5. 一致性仿真常见问题排查参数翻车与数值陷阱5.1 ode45报错Inputs must be a scalar and a square matrix现象运行主脚本时MATLAB直接报错提示矩阵维度不对或函数输入类型错误。原因大多数情况是 weight.m 返回的 L 矩阵含NaN或Inf或者是 model.m 里的变量名遮蔽了函数名导致 ode45 内部把函数句柄解析成了变量。还有一种是调 ode45 时参数顺序写反了把 (t,x) model(t,x,L) 写成了 (t,x) model(x,t,L)x 维度不对就直接报错。解决第一步在调用 ode45 之前打印 L 的尺寸和内容确认 size(L) 是 n×n第二步在 model.m 第一行加断点看进入函数的 x 是否真的是 n×1 列向量第三步检查函数文件名是否与函数名一致MATLAB 对文件名大小写敏感model.m 里写 function dxmodel(...) 就要求文件必须叫 model.m如果系统里同时有同名脚本变量model也会干扰。5.2 一致性收敛方向反了Laplacian符号的坑现象智能体状态不是互相靠近而是越拉越远曲线呈发散状或全部反向增长。原因控制律符号写反。一致性协议要求 dx -Lx如果你误写成 dx Lx因为 L 有零特征值但其它特征值为正正号会导致非零模态指数增长系统直接发散。另一种可能是你用了 A 替换了 L但邻接矩阵本身没有零均值性质用它做控制矩阵即使符号对了也无法保证收敛到一致值。解决回到 weight.m 的输出确认返回的是 L D - A。如果不确定直接把 L 的特征值打出来eig(L)检查是否有且仅有一个零特征值、其余全部大于0。如果最小特征值不是0比如是负值说明 Laplacian 构造有问题大概率是邻居集合 A 的第 i 行没对称上。5.3 切换拓扑索引越界仿真超时与循环取模现象仿真理前几分钟都正常跑到某个时间点突然报 Index exceeds array bounds错误指向 switching.m 里的 L_set(:, :, idx)。原因仿真的 StopTime 超过了 interval × size(L_set,3)。switching.m 里如果用了 min(idx, k) 的 hold 模式不会越界但如果用了不加保护的写法或者 mod 写成了 idx-k而不是 mod(idx-1,k)1就会在周期边界上翻车。解决两种对策一是仿真时长严格限制在 interval*k 以内二是改为循环模式用 mod(idx-1, k)1 确保索引永远落在1~k之间。我建议直接用 cycle 模式这样无论仿真时长多少都不会越界且符合动态网络周期性切换的研究前提。5.4 离散迭代步长过大导致发散稳定域边界现象离散一致性迭代进行到第几十步时状态值从几百跳变成几万呈爆炸式增长曲线根本没法看。原因步长 ε 超过 2/λ_max(L) 的稳定域。离散权重矩阵 W I - εL 的谱半径超过 1每次迭代都会把状态向边界外推。这个问题在稠密图中非常容易触发因为全连接图的 λ_max nε 如果随手设成 0.1n30 时 ελ_max3早就不稳定了。解决写一个通用函数先求 λ_max 再算 εepsilon 0.9 / eigs(L, 1, largestreal)。如果不想每次手动算可以在 discrete_consensus 内部加一步检测if epsilon 2/max(eig(L))直接 error 提示降低步长。这一步只要几行代码能省掉无数次迷茫盯曲线的时间。5.5 Simulink里Scope画面空白信号路由与仿真时长现象点击Run之后Scope窗口一片空白或者只有一条直线看起来像系统根本没动。原因可能有三层。一是仿真时长没改还停在调试时的 0.01s 级别二是Scope没有正确连接状态输出线信号根本没进来三是模型里的矩阵增益填的不是 L而是全零矩阵。最隐蔽的是第三种因为 G0 的系统状态不演化Scope显示的就是一条平线你甚至会开始怀疑数学是不是错了。解决先检查 StopTime设为 20 以上再检查接线从输出端口到Scope之间是否经过了一个你忘了连的Demux最后还是看增益值双击Gain模块确认里面的矩阵不等于0。三者都查完Scope通常就活了。6. 把仿真从2D推到3Dsquare_3D.m的扩展与收敛验证square_3D.m是这份资源里最有意思的一个文件——它把一维状态扩展到三维坐标从标量一致性升级到三维编队一致性。扩展的核心思路是把每个智能体的状态从标量 x 变成 3×1 向量n 个智能体的状态矩阵变成 3×n。控制协议不变仍然是负Laplacian乘以状态矩阵function square_3D() % 三维空间编队一致性四个智能体从随机位置收敛到一致 n 4; pos [0 0 0; 1 0 0; 1 1 0; 0 1 0]; % 期望的方形队形 x0 pos 2 * randn(3, n); % 给初值加扰动 [~, L] weight(n, ring); dt 0.01; T 15; steps T / dt; x_hist zeros(3, n, steps 1); x_hist(:, :, 1) x0; for k 1:steps x_hist(:, :, k 1) x_hist(:, :, k) - dt * (L * x_hist(:, :, k)); end % 画轨迹 for i 1:n plot3(squeeze(x_hist(1,i,:)), squeeze(x_hist(2,i,:)), ... squeeze(x_hist(3,i,:)), LineWidth, 1.5); hold on; end end注意 L * x_hist(:, :, k) 这一行L是 n×n状态矩阵是 3×n矩阵乘法的维度不允许直接相乘所以必须转置让维度对齐算完再转回来。这是3D实现里最容易写错的地方一写错就是内部矩阵尺寸错误。仿真跑完后验证一致性收敛我习惯看两条曲线一条是所有智能体状态与平均值的最大偏差随时间的变化另一条是状态方差的对数曲线。如果偏差单调下降且对数曲线近似一条直线说明收敛速度确实由 λ₂ 控制如果曲线出现弯折说明拓扑切换或增益设置引入了新的动态。从那以后我每次跑多智能体仿真都强制先过一遍三查查 L 的特征值是否合法、查步长是否在稳定域内、查矩阵维度是否对齐然后才丢进仿真器——这三步几乎堵住了我过去百分之八十的翻车。希望帮到你。本文还有配套的精品资源点击获取
