简介这是一份面向自动化控制与人工智能交叉领域的二阶多智能体协同控制 MATLAB 源码包聚焦多智能体系统建模、协同控制策略设计与稳定性分析适合高校师生、科研人员及控制工程师学习仿真。系统采用包含位置与速度的二阶动态模型每名智能体不仅感知自身状态还考虑相邻个体影响以通信协作达成同步、编队或路径规划目标。压缩包共 36 个文件压缩后约 232KB含 24 个 mat 数据文件、4 个 m 脚本、1 个 slx/1 个 slxc Simulink 模型及 xml 配置等其中 m 脚本用于协同算法计算、性能数据提取与轨迹可视化slx 模型可直接构建动态系统仿真环境支持交互式参数调整。运行后能够观察多智能体的位置与速度协同演化验证控制器效果分析通信协议与不确定性对系统稳定性和鲁棒性的影响。该资源已有 1800 人浏览/学习适合需要借助 MATLAB/Simulink 深入理解二阶协同控制、开展算法实验与二次开发的读者。1. 二阶多智能体协同在解决什么问题从“追上”到“达成一致”做多智能体协同控制最常遇到的第一道坎不是算法而是“我该用几阶模型”这个看起来不起眼的问题。如果你只需要控制一群点去同一个目标一阶模型就够了但真实机器人、无人机编队里每个智能体都有惯性和速度状态光控制位置不管速度车会冲过头、无人机编队会震荡。二阶多智能体协同控制要解决的就是这个问题在既有位置又有速度的状态下让所有智能体最终位置一致、速度一致同时保证整个过程不发散。它的典型入口是一致性协议加上Laplacian矩阵描述通信拓扑MATLAB是最容易把这些东西跑通验证的环境。这篇文章写给正在搭仿真、准备做毕设或预研实验的人目标是让你看完整套代码能自己在MATLAB里复现并调出稳定的编队轨迹。2. 模型与一致性协议为什么位置误差反馈要配速度阻尼2.1 二阶积分器模型与“速度不可忽略”的物理意义二阶多智能体系统里面每一个智能体 i 的标准模型是x_i_ddot u_i其中 x_i 是位置u_i 是控制输入。这个模型本质上是一个质点带惯性比一阶积分器模型 x_i_dot u_i 多了一层积分关系。你输入的是加速度而不是速度所以控制量作用之后要先变成速度变化再反映到位置上。这个“延迟一拍”的惯性就是真实运动系统的核心特征。如果你用一阶模型去设计编队仿真里一切正常一上真车就发现车在目标点附近来回冲PID参数怎么都压不住。原因就是你把“速度能瞬间改变”当成了前提。二阶模型里速度是被积分出来的状态控制协议必须同时处理位置误差和速度误差这就是一致性协议里“位置误差项 速度阻尼项”的物理来源。在实际工程里这种“欠阻尼”现象在四旋翼编队里尤其明显位置误差大时给大加速度速度一下冲过头然后又反方向拉。所以二阶系统的协议里速度阻尼项不是可加可不加的“设计选项”而是保证稳定性的必要项。2.2 一致性协议的两种写法连续时间与离散时间二阶多智能体一致性协议最常见的连续时间形式是u_i -k1 * sum_{j in N_i} a_ij * (x_i - x_j) - k2 * v_i其中 k1 是位置耦合增益k2 是速度阻尼增益a_ij 是通信拓扑的邻接矩阵元素N_i 是智能体 i 的邻居集合。这个协议的含义很直白每个智能体看着邻居的位置差往共识方向加速同时用自己的速度做阻尼防止冲过头。我一般会把 u_i 拆成两项来看前一项是“协同项”负责把大家拉在一起后一项是“自阻尼项”负责消耗掉系统的动能。只留协同项系统会永远震荡只留阻尼项各智能体各停各的永远不一致。两者缺一不可。离散时间版本则是一个差分方程x_i[k1] x_i[k] h * v_i[k] v_i[k1] v_i[k] h * (-k1 * sum a_ij*(x_i[k] - x_j[k]) - k2 * v_i[k])h 是仿真步长。这个形式更适合写成MATLAB循环也更容易理解状态是怎么一步步更新的。连续协议用 ode45 解离散协议用矩阵迭代解两者在步长足够小时结果一致但实现逻辑完全不同后面我会把两条路都给你。3. 用 MATLAB 从零跑通最小仿真主脚本 邻居关系矩阵3.1 最小闭环ode45 拉普拉斯误差项先给一个可以直接跑通的最小MATLAB脚本。这个脚本里我用四个智能体组成一条链式拓扑1连接22连接33连接4每个智能体只知道自己直接邻居的位置不知道全局目标。注意正是因为每个智能体只知道局部信息这个仿真才有“多智能体协同”的意义如果谁都知道全局目标那就退化成多个独立的跟踪控制问题了。% main_second_order_consensus.m % 二阶多智能体一致性最小仿真R2020a及以上均可用 clear; clc; close all; % 系统参数 n 4; % 智能体数量 k1 1.5; % 位置耦合增益 k2 2.0; % 速度阻尼增益 % 链式拓扑的邻接矩阵(无权无向图) % 1-2-3-4 一条链 A zeros(n, n); A(1,2) 1; A(2,1) 1; A(2,3) 1; A(3,2) 1; A(3,4) 1; A(4,3) 1; % 拉普拉斯矩阵 L D - AD是度矩阵 D diag(sum(A, 2)); L D - A; % 初始状态: x [x1 v1 x2 v2 x3 v3 x4 v4] x0 [0; 0; 1; 0; 3; 0; 5; 0]; % 时间跨度 tspan [0 15]; % 调用ode45积分状态方程写在一个嵌套函数里 [t, X] ode45((t, X) second_order_consensus_dyn(t, X, L, k1, k2), ... tspan, x0); % 提取位置分量(奇数列是位置偶数列是速度) X_pos X(:, 1:2:end); X_vel X(:, 2:2:end); % 绘图 figure; subplot(2,1,1); plot(t, X_pos, LineWidth, 1.5); xlabel(时间 (s)); ylabel(位置); title(四个智能体的位置轨迹); legend(agent1,agent2,agent3,agent4); grid on; subplot(2,1,2); plot(t, X_vel, LineWidth, 1.5); xlabel(时间 (s)); ylabel(速度); title(四个智能体的速度轨迹); legend(agent1,agent2,agent3,agent4); grid on; % 连续时间二阶一致性状态方程 function dXdt second_order_consensus_dyn(~, X, L, k1, k2) n length(L); x_pos X(1:2:end); x_vel X(2:2:end); % 协同项: -L*x_pos 给出每个智能体与其邻居的位置差之和 % 阻尼项: -k2*x_vel 是速度反馈 u -k1 * (L * x_pos) - k2 * x_vel; dXdt zeros(2*n, 1); dXdt(1:2:end) x_vel; % 位置导数等于速度 dXdt(2:2:end) u; % 速度导数等于控制输入 end这段代码里最关键的是u -k1 * (L * x_pos) - k2 * x_vel这一行。L * x_pos算的是每个智能体和所有邻居的位置差的总和比如智能体2和邻居1、3的位置差是(x2 - x1) (x2 - x3)拉普拉斯矩阵正好把这组运算表达成一次矩阵乘法。这就是为什么多智能体控制代码里几乎处处都是Laplacian——它把“邻居关系”和“位置差求和”这两个语义压缩在了一个矩阵里。3.2 输出与验证怎样判断“协同”已经达成代码跑完之后不要只看曲线“好像合到一起了”我一般会用三个量化指标验证第一看所有位置曲线的末端是否收敛到同一个常数。注意这个常数不一定是初始位置的平均值因为二阶系统有阻尼最后停在哪个共同位置取决于初值和增益不一定等于平均值这是和单积分系统的一个明显区别。第二看速度曲线是否全部归零。二阶一致性的收敛标准是位置相等且速度为零缺一不可。如果速度归零但位置不等说明通信拓扑不连通某些智能体成了孤立点。第三算一下稳态误差也就是最后一刻位置的最大差值final_pos X_pos(end, :); err max(final_pos) - min(final_pos); fprintf(最终位置最大差: %.4f\n, err);如果 err 在 1e-4 量级以下系统就算收敛了。如果只到 1e-1 量级就再也压不下去多半是数值积分精度不够把odeset(RelTol, 1e-8, AbsTol, 1e-8)加上再试。3.3 离散化一步结算工程上更常用的矩阵半离散ode45 适合快速验证算法但如果你要跑大批量参数扫参或者后面要接真实控制器的周期采样逻辑用离散递推更顺手。离散化之后每一步状态更新就是一次纯矩阵运算不需要数值积分器介入。% discrete_second_order_consensus.m % 离散时间二阶一致性仿真与上面同参数 clear; clc; close all; n 4; k1 1.5; k2 2.0; h 0.01; % 采样步长 10ms T 15; % 仿真时长 15s steps T / h; % 同样的链式拓扑 A zeros(n, n); A(1,2) 1; A(2,1) 1; A(2,3) 1; A(3,2) 1; A(3,4) 1; A(4,3) 1; L diag(sum(A,2)) - A; % 构建一步递推矩阵 % 状态向量 [x; v] % x[k1] x[k] h*v[k] % v[k1] v[k] - h*k1*L*x[k] - h*k2*v[k] Phi [eye(n), h*eye(n); -h*k1*L, (1-h*k2)*eye(n)]; x [0; 1; 3; 5]; % 初始位置 v zeros(n, 1); % 初始速度为零 X_log zeros(steps1, n); V_log zeros(steps1, n); X_log(1,:) x; V_log(1,:) v; for k 1:steps state Phi * [x; v]; x state(1:n); v state(n1:end); X_log(k1,:) x; V_log(k1,:) v; end t_axis (0:steps) * h; figure; subplot(2,1,1); plot(t_axis, X_log, LineWidth, 1.5); xlabel(时间 (s)); ylabel(位置); title(离散时间位置轨迹); grid on; legend(agent1,agent2,agent3,agent4); subplot(2,1,2); plot(t_axis, V_log, LineWidth, 1.5); xlabel(时间 (s)); ylabel(速度); title(离散时间速度轨迹); grid on; legend(agent1,agent2,agent3,agent4);离散递推的好处是你能直接检查Phi矩阵的特征值。系统稳定的充分且必要条件是Phi的所有特征值落在复平面单位圆内。我在调参时第一件事就是max(abs(eig(Phi)))是否小于 1这一步能瞬间排除掉一大半“参数不对导致发散”的问题不用跑仿真浪费时间。连续时间版本也有同样的判据对系统矩阵求特征根实部必须全为负。两种判据对应两个世界但调参逻辑是相通的——先说结论后面第4章详细展开。4. 增益怎么调k1/k2 比值、Laplacian 代数连通度与稳定域4.1 特征根判据为什么 k1 不能太小二阶一致性系统的核心稳定条件可以从一个简化模型推出来。假设所有智能体通过某种拓扑通信拉普拉斯矩阵 L 的特征值从小到大排列为 lambda_1 0对应一致状态lambda_2 到 lambda_n 是正特征值。协议 u -k1Lx - k2*v 代入模型后系统的模态方程可以等价成s^2 k2s k1lambda_i 0其中 lambda_i 是拉普拉斯矩阵的第 i 个非零特征值。这个方程就是一个标准的二阶系统特征方程稳定性边界是 k1*lambda_i 必须全为正并且阻尼项 k2 的存在使得系统不会等幅振荡。实际调参时你更关心的不是“能不能稳定”而是“是否会超调震荡”。二阶系统的阻尼比给出直观结论如果 k1*lambda_i 相对于 k2^2 偏大阻尼比就小系统会震好几轮才收敛。常见做法是让 k2 和 k1 * lambda_2 保持一定关系lambda_2 是拉普拉斯矩阵的最小非零特征值也就是代数连通度它刻画了整个网络的连通强度。链式拓扑的 lambda_2 很小所以需要相对大的 k1 或者相对小的 k2 来匹配全连接拓扑的 lambda_2 很大同样的 k1 就可能过阻尼甚至震荡。MATLAB里求 lambda_2 只需要一行lambda_2 sort(eig(L), ascend); fprintf(lambda_2 %.4f\n, lambda_2(2)); % 第一个特征值恒为0记住这个数它就是你整定增益的锚点。拓扑一变lambda_2 就变你的 k1/k2 必须跟着变。这也是为什么“换个拓扑结果完全不对”的根本原因——不是代码错了是增益没跟着新拓扑的连通度重新标定。4.2 一张参数表不同拓扑、不同初速下的推荐范围下面这张表是我在调试中最常用的起始参数。注意它不是万能配方而是“从这些值起步再微调”的参考线。表格里的 k1、k2 适用于状态量归一化到 0-10 范围内的仿真如果你的位置量纲到几百k1 要相应调大。拓扑类型智能体数 nlambda_2 近似值k1 推荐起始k2 推荐起始现象与调整方向链式 1-2-3-440.5861.0~1.51.8~2.2收敛慢就加大 k1震荡就加 k2环形每节点连前后邻居5约1.381.5~2.02.2~2.8连通度更高k1 可略加大但别超 3全连接44.02.0~3.03.0~4.0lambda_2 大k1 偏大会剧烈震荡务必先看特征根星形中心连所有节点51.01.2~1.82.0~2.5中心节点负载最重可单独观察它的轨迹初速度不为零的情况需要特别注意。如果每个智能体初始速度不同系统要先消除速度差这个阶段位置差会先变大再收敛属于正常过程。但如果初始速度方向杂乱且幅值大建议把 k2 提高到 k1 的 1.5 倍以上先用阻尼稳住速度再慢慢靠位置耦合拉齐否则会在前 1 秒内飞出很远。4.3 验证方法特征根、超调量与收敛误差三个指标调完参数不要只看曲线“感觉还行”我通常在MATLAB里做三件事量化验证。第一看特征根或转移矩阵谱半径。连续系统用eig(A_sys)看实部离散系统用max(abs(eig(Phi)))看模长。这个指标是稳定性的硬边界任何一组参数如果不满足直接否掉。第二看超调量。取速度曲线最大值与稳态值0的相对偏差记录前 5 秒内速度峰值有没有超过初始速度的最大值。如果超调超过初始极值的 20%说明阻尼偏弱需要加 k2。第三看收敛时间。以“位置最大差首次进入 0.01 以内”为准pos_range max(X_log, [], 2) - min(X_log, [], 2); converged_idx find(pos_range 0.01, 1, first); if ~isempty(converged_idx) fprintf(收敛时间: %.2f s\n, t_axis(converged_idx)); else fprintf(未在仿真时间内收敛\n); end收敛时间这个指标对后续设计非常有用你想让编队在 5 秒内成形收敛时间就必须压到 5 秒以内。单纯看最终稳态姿态度你没法判断系统是否满足时间约束。5. 二阶协同在 MATLAB 里反复翻车的 5 个点现象、原因与解决5.1 现象位置曲线发散数值跑到 1e10 量级原因拉普拉斯矩阵拼错了。最常见的是把度矩阵 D 写成了每个节点与其邻居的权重之和的转置导致L D - A不对称出现了一个正实部特征根。链式拓扑里 D(1,1) 应该是 1不是 2因为节点1只有一个邻居。解决先用eig(L)检查所有非零特征值是否全为正实数值第一个特征值必须等于 0。如果不满足逐行检查 A 矩阵的对称性。5.2 现象曲线振荡但收敛后最终结果正确原因k1 相对 k2 偏大阻尼比不足。系统处在欠阻尼状态收敛路径很长。这不是代码错误是参数没匹配拓扑的代数连通度。解决按第 4.1 节的方法先算 lambda_2然后把 k1 降到接近 k2^2 / lambda_2 的量级。比如 lambda_2 是 0.586k2 是 1.8k1 起始大约在 1.0~1.5不要一上来就用 3。5.3 现象四个智能体中有一个“掉队”最后停在错误位置原因通信拓扑不连通。比如邻接矩阵只把 1 连到 2、2 连到 3却忘了连 4那智能体 4 就是孤立节点没有误差项驱动它移动它会永远停在初始位置。解决检查每个节点的度是否至少为 1。MATLAB 里一行sum(A, 2)就能看到每个节点的邻居数任何一行出现 0 就说明拓扑断了。5.4 现象ode45 和离散递推的结果明显不一致轨迹形状完全不同原因步长 h 取得太大或者连续系统里存在高频模态。二阶系统的模态频率取决于sqrt(k1 * lambda_i)如果 h 大于振荡周期的 1/10离散化误差就会显著扭曲轨迹。解决把 h 降到 0.01验证结果稳定后再逐步加大。连续时间用odeset(RelTol, 1e-8)提高积分精度。5.5 现象初速度不为零时前几步位置冲出初始范围好几倍看起来像“爆炸”原因这是二阶系统惯性的正常表现。位置变化率等于当前速度如果初始速度是 5第一个步长 0.01 秒位置就会跳 0.05连续积累看起来就像“先冲出去”。但如果你用一阶系统的思维去理解就会误判为代码错误。解决先看速度曲线有没有收敛。只要速度先冲到负值再回调说明阻尼在起作用系统是中性的超调阶段。如果速度一路朝一个方向跑停不下来才是真发散回到第 5.1 条查特征根。6. 进阶给协议加通信时延用单步缓存验证鲁棒性真实系统里通信是不可能零延时的WiFi、串口、内部总线都有延迟。给上面的离散模型加一个固定的单步时延是验证协议鲁棒性最轻量的方法。实现上只需要把“邻居位置”从当前步换成上一步就行代码改动极小% delayed_consensus.m % 固定单步通信时延的离散一致性仿真 clear; clc; close all; n 4; k1 1.0; k2 2.0; % 有时延时增益要适当调低 h 0.01; T 20; steps T / h; A zeros(n, n); A(1,2) 1; A(2,1) 1; A(2,3) 1; A(3,2) 1; A(3,4) 1; A(4,3) 1; L diag(sum(A,2)) - A; x [0; 1; 3; 5]; v zeros(n, 1); x_delay x; % 缓存一步之前的位置 X_log zeros(steps1, n); V_log zeros(steps1, n); X_log(1,:) x; V_log(1,:) v; for k 1:steps % 用 x_delay 参与协同项而不是用 x u -k1 * (L * x_delay) - k2 * v; x_new x h * v; v_new v h * u; x_delay x; % 当前步变成下一步的滞后数据 x x_new; v v_new; X_log(k1,:) x; V_log(k1,:) v; end这段代码的用意是模拟“智能体看到的邻居位置总是早了一个采样周期”的通信延迟。你把这段代码跟第 3.3 节的无时延版本对比会发现同样参数下有时延的系统收敛速度会明显变慢甚至出现小幅震荡。这就是时延对稳定域的压缩效应。我的经验是单步时延引入后k1 要下调 20%~30%k2 基本不动才能恢复到接近原来的动态品质。如果加了时延之后系统开始震荡优先检查的依然是特征根只不过这时特征多项式里多了 z^{-1} 项判断标准的说法变了实际调试手段没变画出轨迹、试调 k1、观察前几秒的超调幅度。这三个动作比你翻多少理论推导都更直接。我自己的习惯是每次把仿真参数调稳定之后都会刻意加上一个单步时延再跑一遍这个动作已经帮我拦下了好几次“仿真能过、实车就废”的尴尬。多花两分钟加一个缓存数组比后面在现场查半天黑匣子值太多。希望这套从模型、代码到调试的路径能帮你在二阶多智能体协同这个方向上少走一段弯路。本文还有配套的精品资源点击获取