議、編隊控制與事件觸發(fā))
簡介面向多智能體系統(tǒng)初學者與MATLAB仿真實踐者的示例代碼包圍繞MAS中代理建模、環(huán)境感知、消息交互、協(xié)調決策與仿真驗證等核心環(huán)節(jié)展開適合希望用腳本方式理解多代理協(xié)作機制、搭建小型實驗原型的讀者。壓縮包共4個文件以3個m腳本和1個txt許可說明為主體量約4KB輕量便于查閱與二次修改腳本承擔代理狀態(tài)更新、位置計算及服務端協(xié)調等示例功能許可文件用于說明使用授權。已有1806人學習下載可作為課程演示、課題預研或自學練手的參考。內容側重多智能體服務器如何注冊管理代理、分配任務、維護通信通道并定義協(xié)同規(guī)則同時保留監(jiān)控評估的擴展空間讀者可據此觀察服務端創(chuàng)建與跟蹤代理實例、按能力派發(fā)任務的過程并把代理、環(huán)境、交互、協(xié)調與決策等概念落到可運行的MATLAB代碼骨架上。1. 多智能體 MATLAB 代碼拆包先跑通 5 個智能體的編隊收斂再談改協(xié)議拿到一份多智能體 MATLAB 代碼多數(shù)人的第一反應是直接 run 主腳本看那幾條曲線收不收斂——收斂就收工散了就懷疑算法寫錯。我的經驗恰好相反MAS 仿真翻車八成不是協(xié)議本身的問題而是通信拓撲、采樣步長、初始狀態(tài)這三處沒對齊。這套代碼真正的價值是把一致性協(xié)議、編隊偏移、分布式控制器和仿真外殼拆成了互相獨立、可單獨替換的模塊改一處不牽連其余。它適合正在做編隊、蜂擁、分布式優(yōu)化驗證的人也適合想把 Simulink 里那堆手工連線換成腳本化流程、方便批量跑參數(shù)的場景。2. 一致性協(xié)議落到矩陣上拉普拉斯矩陣的三種寫法與增益邊界多智能體系統(tǒng)在代碼層面的全部秘密幾乎都壓縮在一個矩陣里拉普拉斯矩陣 L。協(xié)議寫成微分方程是 ? -Lx寫成代碼就是一次矩陣乘法。所以這一章不急著調控制器先把 L 造對——L 錯了后面所有參數(shù)整定都是白費功夫。2.1 鄰接矩陣的三種建法鄰接矩陣 A 描述“誰能收到誰的信息”是整個 MAS 的地基。小規(guī)模固定拓撲手寫邊表最直觀中等規(guī)模按幾何距離生成通信半徑內的智能體連邊大規(guī)?；蛴鞋F(xiàn)成拓撲數(shù)據時從稀疏矩陣導入。三種寫法各有各的翻車點手寫最容易犯的錯是只賦值單向。% 建法一手寫邊表適合 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 % 建法三從稀疏表恢復適合大規(guī)模或外部拓撲 % A sparse(rowIdx, colIdx, 1, N, N); A full(A | A);第一種寫法里那句成對賦值不是啰嗦是我踩過的坑早期圖省事只寫 A(i,j)1結果 L 不對稱特征值出現(xiàn)負實部曲線前半段看著挺正常后半段直接飛出去。第二種寫法的D 0同理——不加這一條對角線會被判成“距離 0 以內”自我連邊L 的對角元憑空多出 N收斂速度全亂。第三種寫法補一句A | A做對稱化防止上游數(shù)據只給了半邊。提示建完 A 之后立刻斷言isequal(A, A)。這一步成本幾乎為零能擋掉后面一大半莫名其妙的發(fā)散。2.2 從拉普拉斯矩陣到狀態(tài)方程有了 A拉普拉斯矩陣一行就能算出來L D - AD 是度矩陣。這個式子的物理含義是“我拿自己和所有鄰居的差做加權和”所以 L 的每一行和恒為零全 1 向量永遠是它的零特征值對應向量——這正是所有智能體最終停在同一值的數(shù)學來源。收斂速度由第二小特征值 λ? 決定代數(shù)連通度越高收斂越快。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(代數(shù)連通度 lambda_2 %.4f\n, lam2(min(2, numel(lam2))));三條自檢對應三類錯誤不對稱會讓 L 出現(xiàn)復特征值行和不為零說明 D 求錯了維度λ? 接近零說明圖不連通——不連通時一致性協(xié)議不是發(fā)散而是分裂成幾個各自收斂的簇這種“部分收斂”最迷惑人曲線看著像要收斂又永遠差一點。連續(xù)時間協(xié)議 ? -Lx 在代碼里通常用前向歐拉離散x(k1) (I - h·L)x(k)穩(wěn)定條件是 h 2/λ_max(L)。這個邊界不是理論擺設我見過 h 設成 0.5、λ_max 大約 4 的情形正好壓在邊上前幾步正常二十步后開始振蕩放大。2.3 五智能體最小可跑腳本把上面的東西拼起來就是一段不到 20 行、但能完整驗證協(xié)議的最小腳本。跑它的意義在于在引入編隊、飽和、時延之前先確認你的拓撲和步長是干凈的。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]; % 初始狀態(tài)故意打散 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(末態(tài)最大偏差 %.3e理論 lambda_2 %.4f\n, err, ... sort(eig(L))(2));邏輯是一條直線加權平均 → 矩陣乘 → 迭代。參數(shù)上有三處會直接影響結果。h必須小于 2/λ_max實用做法是取 0.1/λ_max寧可多跑幾步也別貼邊界。T要留足經驗值是 10/λ? 以上λ? 越小需要的時間越長圖越“細長”越慢。x0最好手動打散別用rand直接生成——調試階段你需要能心算出期望均值隨機初始值會讓你失去這個參照。收斂判據也別只看最后一步。我在腳本里習慣再算一條errHist max(abs(X - mean(x0)), [], 1)把誤差歷史畫出來。真收斂是指數(shù)下降的直線半對數(shù)坐標下如果誤差在某個量級上壓不下去多半是 λ? 太小或步長偏大導致的殘余振蕩光看末態(tài)一個數(shù)看不出來。3. 編隊控制的可調旋鈕偏移量、飽和限幅與拓撲切換一致性只要求大家“停在一起”編隊要求“停成某個形狀”。從一致性到編隊代碼上只多了一個偏移量矩陣但工程上的復雜度翻了一倍——因為偏移量引入了幾何約束幾何約束又逼著你處理飽和、避碰和拓撲變化。3.1 把方針翻譯成偏移量矩陣編隊協(xié)議的標準寫法是 u_i -Σ a_ij[(x_i - p_i) - (x_j - p_j)]其中 p_i 是第 i 個智能體的期望相對位置。寫成矩陣就是 u -L(X - P)。這里最容易搞混的是 P 的語義它給的是相對偏移形心在原點還是絕對目標位置代碼里必須統(tǒng)一否則整個隊形會整體漂移到畫面外。% 二維平面上的五個位置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 當前狀態(tài)P: N×2 期望偏移L: N×N 拉普拉斯 E X - P; % 先減偏移再走一致性 Xn X - h * (L * E); endE X - P這一行的順序不能反。寫成L*X - P是另一個東西那是帶絕對參考的跟蹤寫反了隊形會一邊收斂一邊整體平移看起來像“慢悠悠地飄走”。另外 P 的行數(shù)必須等于 N我在批量換 N 的腳本里吃過虧P 是硬編碼的 5 行N 改成 7 之后 MATLAB 不報錯而是自動做隱式擴展隊形直接算成一團亂麻——這種不報錯的 bug 最難查。注意換 N 的時候A、L、P、x0 四處的維度必須一起改。建議在函數(shù)入口加assert(size(P,1) size(L,1))一秒定位。3.2 飽和限幅與避碰項兩個位置別寫反真實執(zhí)行器有速度上限仿真里不加飽和編隊初期會出現(xiàn)幾十米每秒的瞬時速度看著收斂很漂亮換成實機直接失控。限幅的寫法有講究要等比縮放不能逐分量裁剪否則方向會變隊形走形。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, ...)保證只往下壓不放上去這一點在調試時很關鍵——如果發(fā)現(xiàn)某個智能體始終貼著 vmax 跑說明它被拓撲甩開了問題在圖不在限幅。避碰項和限幅是兩個獨立的東西別塞進同一個 if 里。常見做法是在 u 后面疊加一個斥力項距離小于安全半徑 d_safe 時沿相對位置方向加一個反比項。斥力增益不能給大否則會和一致性項打架出現(xiàn)兩個智能體高頻抖動、位置幾乎不動的現(xiàn)象——這種抖動在曲線上看是毛刺在相圖上看是一個極限環(huán)。3.3 采樣步長與拓撲切換的穩(wěn)定邊界固定拓撲下 h 2/λ_max 就夠了一旦拓撲會切換通信中斷、鄰居變化穩(wěn)定條件變成對所有可能拓撲的并集都成立也就是 h 2/max(λ_max over all modes)。這條規(guī)則在實際代碼里經常被忽略因為單跑某個拓撲是穩(wěn)的一切換就炸。% 拓撲集合正常、斷邊、鏈路退化三種模式 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);這段循環(huán)算的是“最保守的公共步長”。實用取值我一般按 0.1 倍上限代價是仿真步數(shù)變多但換來的是切換瞬間不會出現(xiàn)數(shù)值爆點。如果確實需要大 step 提速那就得改成嚴格離散協(xié)議用行隨機矩陣而不是 I - hL犧牲一點物理直覺換取無條件穩(wěn)定。切換邏輯本身還有一個隱藏坑切換時刻落在兩個采樣點中間時狀態(tài)更新用的是舊拓撲還是新拓撲仿真里通常按采樣點對齊切換只在 k 為整數(shù)時生效。如果你的切換由事件觸發(fā)比如距離閾值一定要把觸發(fā)時刻量化到采樣網格否則同一份代碼跑兩次結果不同排查起來毫無頭緒。4. 仿真工程化參數(shù)掃描、數(shù)據歸檔與批量繪圖單次仿真跑通只是及格線。真要做對比、寫報告、給導師或同事看結論你需要的是“一次配置、批量跑完、結果能追溯”。這一章把腳本改造成可復用的三層結構仿真函數(shù)、掃描腳本、繪圖腳本。4.1 把單次仿真封成一個函數(shù)函數(shù)化的第一原則是輸入只有參數(shù)輸出只有數(shù)據中間不畫圖、不打印、不依賴工作區(qū)變量。這樣它才能被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三維數(shù)組X(N,2,K1)是這里唯一需要適應的地方第一維智能體、第二維狀態(tài)分量、第三維時間。好處是squeeze之后每一步都是標準 N×2 矩陣直接喂給step_sat壞處是 N 大、K 大時內存吃得兇N50、K4000 大約 3 MB還能接受再大就應該改成只存抽樣點。out.cfg cfg這一句別省——三個月后你看到一份.mat唯一能回答“這組參數(shù)是怎么來的”的就是它。4.2 參數(shù)掃描與結果歸檔掃描的寫法取決于你要掃幾維。單參數(shù)比如增益或步長用一層循環(huán)雙參數(shù)步長 × 通信半徑用兩層結果存成矩陣超過三維就別硬掃了改隨機采樣。歸檔用帶時間戳的文件名別用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這個模式值得單獨說一句先在腳本頂部構造一份完整配置循環(huán)里只覆寫變化的字段。這樣新增參數(shù)時只需改一處不會出現(xiàn)“某組參數(shù)忘了設、用了默認值”的隱蔽錯誤。我見過最典型的一次是 vmax 忘了傳結構體里沒有這個字段函數(shù)在cfg.vmax處直接報錯——這還算好的報錯至少讓你知道如果代碼里寫了getfield(cfg,vmax,10)給了默認值那就變成靜默跑完、結論全錯。歸檔目錄我會在腳本開頭if ~exist(runs,dir), mkdir(runs); end別指望手動建。4.3 批量繪圖一眼看出哪組參數(shù)崩了繪圖的目的不是好看是讓“哪組參數(shù)崩了”這件事在 0.5 秒內可判斷。做法是把所有組的誤差曲線畫在同一張半對數(shù)圖上收斂的是直線下降振蕩的是鋸齒發(fā)散的直接沖頂。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 看起來差不多看不出收斂速率差異對數(shù)坐標下斜率直接就是收斂速度三組參數(shù)誰快誰慢一目了然。另外DisplayName一定要寫不然 legend 出來的是一串data1、data2回頭對不上。如果某條曲線在中段出現(xiàn)平臺然后再次下降那通常是飽和起了作用前期被 vmax 卡住誤差下降變慢進入線性區(qū)后恢復指數(shù)收斂。這個“拐點”位置本身是有信息的它告訴你限幅是不是給緊了。5. 排查收斂曲線不對時先查這五件事排查順序很重要。我的習慣是先看矩陣再看數(shù)值最后才懷疑協(xié)議——因為協(xié)議是最不容易寫錯的部分而矩陣和步長是最容易在復制粘貼中出問題的部分。5.1 矩陣與拓撲層面的三個坑現(xiàn)象一曲線看著在收斂但每個智能體停在不同的值上。原因是圖不連通或者本該無向的圖被寫成了有向且非強連通。無向連通圖的一致性值是全體初始狀態(tài)的算術平均一旦分裂成兩個連通分量就變成兩個各自的平均值。解決辦法是跑一遍連通性檢查numel(unique(conncomp(graph(A)))) 1不連通就回到邊表找漏掉的邊?，F(xiàn)象二前期正常十幾步之后突然發(fā)散且越跑越陡。步長超過穩(wěn)定邊界。查法很直接h和2/max(real(eig(L)))比一下超了就減半。我踩過一次更隱蔽的版本——代碼里 H 是大寫的局部變量、h 是小寫的全局參數(shù)循環(huán)里用的是 H檢查時改的是 h改了半天沒反應血淚經驗是變量名別玩這種字母大小寫的玄學?，F(xiàn)象三改了 N 之后報維度錯誤或者不報錯但隊形完全亂掉。硬編碼索引沒跟著改。典型位置是 A 的邊表、P 的行數(shù)、x0 的長度以及eye(N)里那個 N。解決方式是全腳本搜一遍數(shù)字常量把與規(guī)模相關的都換成size(...)或者從 cfg 里取。5.2 數(shù)值與工程層面的兩個坑現(xiàn)象四曲線里出現(xiàn) NaN從某一步開始整條線消失。零除或溢出。最可能是初始位置重合導致距離為零或者限幅里vmax / n的 n 是 0。解決辦法是給所有分母加max(x, eps)同時在循環(huán)里加一句assert(all(isfinite(u(:))), step %d 出現(xiàn) NaN, k)——早報錯比跑完再看圖省時間得多?，F(xiàn)象五腳本跑和 Simulink 跑結果對不上形狀相似但數(shù)值差一截。求解器步長和采樣步長不一致。Simulink 默認變步長求解器在飽和切換點會自動縮步腳本里是固定步長兩者必然有差異。辦法是把 Simulink 里的求解器改成固定步長步長填成和腳本一致再對比。如果還差就查模塊里的采樣時間有沒有設成繼承-1繼承來的采樣時間經常和你想的不是一個值。6. 進階替換點把周期通信換成事件觸發(fā)附一份驗收清單周期通信的假設是“每個采樣點所有智能體都要和鄰居交換一次信息”在仿真里不要錢在實機上就是通信帶寬和功耗。事件觸發(fā)一致性把這條假設松掉只在狀態(tài)偏差超過閾值時才通信其余時刻用保持器維持上一次的值。代碼上的改動量小得出奇大約十幾行但收益在規(guī)模大、通信受限的場景里非常明顯。6.1 事件觸發(fā)一致性的替換點核心是把“當前狀態(tài)”換成“上次廣播的狀態(tài)”再維護一個觸發(fā)誤差。觸發(fā)條件通常是 ‖e_i‖ σ‖z_i‖e_i 是當前狀態(tài)與上次廣播狀態(tài)之差z_i 是鄰居狀態(tài)的加權和σ 是觸發(fā)閾值越小越接近周期通信。function [X, nTrig] et_consensus(X0, L, h, K, sigma) N size(X0, 1); X zeros([size(X0), K1]); X(:,:,1) X0; Xhat X0; % 上次廣播出去的狀態(tài) nTrig zeros(N, 1); % 每個智能體的觸發(fā)次數(shù) for k 1:K Xc X(:,:,k); z -L * Xhat; % 用廣播狀態(tài)算控制量 e Xc - Xhat; % 觸發(fā)誤差 % 逐智能體判斷是否觸發(fā)觸發(fā)者更新廣播值 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那事件觸發(fā)形同虛設你還是每步都在用最新狀態(tài)。need是一個 N×1 的邏輯向量用邏輯索引整行更新比 for 循環(huán)快也比循環(huán)不容易寫錯下標。nTrig返回每個智能體的觸發(fā)次數(shù)這個數(shù)字很能說明問題如果某個智能體的觸發(fā)次數(shù)接近 K說明它一直在被鄰居甩開問題在拓撲不在閾值。調參上有兩條經驗。sigma從 0.1 開始試太大會出現(xiàn)“長時間不通信、隊形明顯滯后”的現(xiàn)象太小則退化成周期通信、省不下什么。另外事件觸發(fā)會引入額外的收斂時間同樣 T 下誤差可能比周期通信大一個量級比較兩者時要在相同觸發(fā)次數(shù)下比而不是相同時間內比——這一點在寫結論時經常被忽略。6.2 一套我自己在用的驗收清單代碼改完之后跑一遍下面這幾項比反復看曲線靠譜。這套清單是從多次“明明仿真收斂、換組參數(shù)就崩”的經歷里攢出來的。檢查項判據不通過時先看哪鄰接矩陣對稱性isequal(A, A)為真邊表是否只寫了單向拉普拉斯行和每行絕對值 1e-12度矩陣是否按行求和代數(shù)連通度 λ?明顯大于 0 1e-3圖是否連通、有無孤立點步長裕度h 0.1 × 2/λ_max拓撲切換時取所有模式的公共上限有限性全程無 NaN / Inf分母是否加了 eps收斂值無向連通時等于初始均值是否有向圖但按無向算飽和觸頂率大部分時間不貼 vmaxvmax 給緊了或拓撲把某個體甩開觸發(fā)占比事件觸發(fā)遠小于 1且各體差異不大sigma 是否過小、某體是否孤立表格里最容易被跳過的是“飽和觸頂率”但它往往是最早暴露拓撲問題的指標。我一般會在仿真函數(shù)里順手記一條satRatio mean(n vmax)跑完看這個數(shù)超過 0.5 就說明整段時間都在限幅區(qū)控制器根本沒按線性協(xié)議工作這時候去調增益全是白費。還有一處建議把 λ?、h 上限、飽和觸頂率、末態(tài)誤差這四個數(shù)打在一行日志里輸出。曲線可以騙人尤其是只看末態(tài)的時候這四個標量不太會。我自從被一次“圖不連通但曲線看起來收斂”坑掉整整兩天之后每次改完拓撲做的第一件事都是先跑那三條斷言再畫圖。斷言不過就不看曲線因為那時候的曲線只會給你錯誤的信息。希望幫到你。本文還有配套的精品資源點擊獲取 多智能體MATLAB代碼拆解：一致性協(xié)議、編隊控制與事件觸發(fā)
簡介面向多智能體系統(tǒng)初學者與MATLAB仿真實踐者的示例代碼包圍繞MAS中代理建模、環(huán)境感知、消息交互、協(xié)調決策與仿真驗證等核心環(huán)節(jié)展開適合希望用腳本方式理解多代理協(xié)作機制、搭建小型實驗原型的讀者。壓縮包共4個文件以3個m腳本和1個txt許可說明為主體量約4KB輕量便于查閱與二次修改腳本承擔代理狀態(tài)更新、位置計算及服務端協(xié)調等示例功能許可文件用于說明使用授權。已有1806人學習下載可作為課程演示、課題預研或自學練手的參考。內容側重多智能體服務器如何注冊管理代理、分配任務、維護通信通道并定義協(xié)同規(guī)則同時保留監(jiān)控評估的擴展空間讀者可據此觀察服務端創(chuàng)建與跟蹤代理實例、按能力派發(fā)任務的過程并把代理、環(huán)境、交互、協(xié)調與決策等概念落到可運行的MATLAB代碼骨架上。1. 多智能體 MATLAB 代碼拆包先跑通 5 個智能體的編隊收斂再談改協(xié)議拿到一份多智能體 MATLAB 代碼多數(shù)人的第一反應是直接 run 主腳本看那幾條曲線收不收斂——收斂就收工散了就懷疑算法寫錯。我的經驗恰好相反MAS 仿真翻車八成不是協(xié)議本身的問題而是通信拓撲、采樣步長、初始狀態(tài)這三處沒對齊。這套代碼真正的價值是把一致性協(xié)議、編隊偏移、分布式控制器和仿真外殼拆成了互相獨立、可單獨替換的模塊改一處不牽連其余。它適合正在做編隊、蜂擁、分布式優(yōu)化驗證的人也適合想把 Simulink 里那堆手工連線換成腳本化流程、方便批量跑參數(shù)的場景。2. 一致性協(xié)議落到矩陣上拉普拉斯矩陣的三種寫法與增益邊界多智能體系統(tǒng)在代碼層面的全部秘密幾乎都壓縮在一個矩陣里拉普拉斯矩陣 L。協(xié)議寫成微分方程是 ? -Lx寫成代碼就是一次矩陣乘法。所以這一章不急著調控制器先把 L 造對——L 錯了后面所有參數(shù)整定都是白費功夫。2.1 鄰接矩陣的三種建法鄰接矩陣 A 描述“誰能收到誰的信息”是整個 MAS 的地基。小規(guī)模固定拓撲手寫邊表最直觀中等規(guī)模按幾何距離生成通信半徑內的智能體連邊大規(guī)?；蛴鞋F(xiàn)成拓撲數(shù)據時從稀疏矩陣導入。三種寫法各有各的翻車點手寫最容易犯的錯是只賦值單向。% 建法一手寫邊表適合 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 % 建法三從稀疏表恢復適合大規(guī)模或外部拓撲 % A sparse(rowIdx, colIdx, 1, N, N); A full(A | A);第一種寫法里那句成對賦值不是啰嗦是我踩過的坑早期圖省事只寫 A(i,j)1結果 L 不對稱特征值出現(xiàn)負實部曲線前半段看著挺正常后半段直接飛出去。第二種寫法的D 0同理——不加這一條對角線會被判成“距離 0 以內”自我連邊L 的對角元憑空多出 N收斂速度全亂。第三種寫法補一句A | A做對稱化防止上游數(shù)據只給了半邊。提示建完 A 之后立刻斷言isequal(A, A)。這一步成本幾乎為零能擋掉后面一大半莫名其妙的發(fā)散。2.2 從拉普拉斯矩陣到狀態(tài)方程有了 A拉普拉斯矩陣一行就能算出來L D - AD 是度矩陣。這個式子的物理含義是“我拿自己和所有鄰居的差做加權和”所以 L 的每一行和恒為零全 1 向量永遠是它的零特征值對應向量——這正是所有智能體最終停在同一值的數(shù)學來源。收斂速度由第二小特征值 λ? 決定代數(shù)連通度越高收斂越快。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(代數(shù)連通度 lambda_2 %.4f\n, lam2(min(2, numel(lam2))));三條自檢對應三類錯誤不對稱會讓 L 出現(xiàn)復特征值行和不為零說明 D 求錯了維度λ? 接近零說明圖不連通——不連通時一致性協(xié)議不是發(fā)散而是分裂成幾個各自收斂的簇這種“部分收斂”最迷惑人曲線看著像要收斂又永遠差一點。連續(xù)時間協(xié)議 ? -Lx 在代碼里通常用前向歐拉離散x(k1) (I - h·L)x(k)穩(wěn)定條件是 h 2/λ_max(L)。這個邊界不是理論擺設我見過 h 設成 0.5、λ_max 大約 4 的情形正好壓在邊上前幾步正常二十步后開始振蕩放大。2.3 五智能體最小可跑腳本把上面的東西拼起來就是一段不到 20 行、但能完整驗證協(xié)議的最小腳本。跑它的意義在于在引入編隊、飽和、時延之前先確認你的拓撲和步長是干凈的。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]; % 初始狀態(tài)故意打散 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(末態(tài)最大偏差 %.3e理論 lambda_2 %.4f\n, err, ... sort(eig(L))(2));邏輯是一條直線加權平均 → 矩陣乘 → 迭代。參數(shù)上有三處會直接影響結果。h必須小于 2/λ_max實用做法是取 0.1/λ_max寧可多跑幾步也別貼邊界。T要留足經驗值是 10/λ? 以上λ? 越小需要的時間越長圖越“細長”越慢。x0最好手動打散別用rand直接生成——調試階段你需要能心算出期望均值隨機初始值會讓你失去這個參照。收斂判據也別只看最后一步。我在腳本里習慣再算一條errHist max(abs(X - mean(x0)), [], 1)把誤差歷史畫出來。真收斂是指數(shù)下降的直線半對數(shù)坐標下如果誤差在某個量級上壓不下去多半是 λ? 太小或步長偏大導致的殘余振蕩光看末態(tài)一個數(shù)看不出來。3. 編隊控制的可調旋鈕偏移量、飽和限幅與拓撲切換一致性只要求大家“停在一起”編隊要求“停成某個形狀”。從一致性到編隊代碼上只多了一個偏移量矩陣但工程上的復雜度翻了一倍——因為偏移量引入了幾何約束幾何約束又逼著你處理飽和、避碰和拓撲變化。3.1 把方針翻譯成偏移量矩陣編隊協(xié)議的標準寫法是 u_i -Σ a_ij[(x_i - p_i) - (x_j - p_j)]其中 p_i 是第 i 個智能體的期望相對位置。寫成矩陣就是 u -L(X - P)。這里最容易搞混的是 P 的語義它給的是相對偏移形心在原點還是絕對目標位置代碼里必須統(tǒng)一否則整個隊形會整體漂移到畫面外。% 二維平面上的五個位置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 當前狀態(tài)P: N×2 期望偏移L: N×N 拉普拉斯 E X - P; % 先減偏移再走一致性 Xn X - h * (L * E); endE X - P這一行的順序不能反。寫成L*X - P是另一個東西那是帶絕對參考的跟蹤寫反了隊形會一邊收斂一邊整體平移看起來像“慢悠悠地飄走”。另外 P 的行數(shù)必須等于 N我在批量換 N 的腳本里吃過虧P 是硬編碼的 5 行N 改成 7 之后 MATLAB 不報錯而是自動做隱式擴展隊形直接算成一團亂麻——這種不報錯的 bug 最難查。注意換 N 的時候A、L、P、x0 四處的維度必須一起改。建議在函數(shù)入口加assert(size(P,1) size(L,1))一秒定位。3.2 飽和限幅與避碰項兩個位置別寫反真實執(zhí)行器有速度上限仿真里不加飽和編隊初期會出現(xiàn)幾十米每秒的瞬時速度看著收斂很漂亮換成實機直接失控。限幅的寫法有講究要等比縮放不能逐分量裁剪否則方向會變隊形走形。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, ...)保證只往下壓不放上去這一點在調試時很關鍵——如果發(fā)現(xiàn)某個智能體始終貼著 vmax 跑說明它被拓撲甩開了問題在圖不在限幅。避碰項和限幅是兩個獨立的東西別塞進同一個 if 里。常見做法是在 u 后面疊加一個斥力項距離小于安全半徑 d_safe 時沿相對位置方向加一個反比項。斥力增益不能給大否則會和一致性項打架出現(xiàn)兩個智能體高頻抖動、位置幾乎不動的現(xiàn)象——這種抖動在曲線上看是毛刺在相圖上看是一個極限環(huán)。3.3 采樣步長與拓撲切換的穩(wěn)定邊界固定拓撲下 h 2/λ_max 就夠了一旦拓撲會切換通信中斷、鄰居變化穩(wěn)定條件變成對所有可能拓撲的并集都成立也就是 h 2/max(λ_max over all modes)。這條規(guī)則在實際代碼里經常被忽略因為單跑某個拓撲是穩(wěn)的一切換就炸。% 拓撲集合正常、斷邊、鏈路退化三種模式 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);這段循環(huán)算的是“最保守的公共步長”。實用取值我一般按 0.1 倍上限代價是仿真步數(shù)變多但換來的是切換瞬間不會出現(xiàn)數(shù)值爆點。如果確實需要大 step 提速那就得改成嚴格離散協(xié)議用行隨機矩陣而不是 I - hL犧牲一點物理直覺換取無條件穩(wěn)定。切換邏輯本身還有一個隱藏坑切換時刻落在兩個采樣點中間時狀態(tài)更新用的是舊拓撲還是新拓撲仿真里通常按采樣點對齊切換只在 k 為整數(shù)時生效。如果你的切換由事件觸發(fā)比如距離閾值一定要把觸發(fā)時刻量化到采樣網格否則同一份代碼跑兩次結果不同排查起來毫無頭緒。4. 仿真工程化參數(shù)掃描、數(shù)據歸檔與批量繪圖單次仿真跑通只是及格線。真要做對比、寫報告、給導師或同事看結論你需要的是“一次配置、批量跑完、結果能追溯”。這一章把腳本改造成可復用的三層結構仿真函數(shù)、掃描腳本、繪圖腳本。4.1 把單次仿真封成一個函數(shù)函數(shù)化的第一原則是輸入只有參數(shù)輸出只有數(shù)據中間不畫圖、不打印、不依賴工作區(qū)變量。這樣它才能被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三維數(shù)組X(N,2,K1)是這里唯一需要適應的地方第一維智能體、第二維狀態(tài)分量、第三維時間。好處是squeeze之后每一步都是標準 N×2 矩陣直接喂給step_sat壞處是 N 大、K 大時內存吃得兇N50、K4000 大約 3 MB還能接受再大就應該改成只存抽樣點。out.cfg cfg這一句別省——三個月后你看到一份.mat唯一能回答“這組參數(shù)是怎么來的”的就是它。4.2 參數(shù)掃描與結果歸檔掃描的寫法取決于你要掃幾維。單參數(shù)比如增益或步長用一層循環(huán)雙參數(shù)步長 × 通信半徑用兩層結果存成矩陣超過三維就別硬掃了改隨機采樣。歸檔用帶時間戳的文件名別用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這個模式值得單獨說一句先在腳本頂部構造一份完整配置循環(huán)里只覆寫變化的字段。這樣新增參數(shù)時只需改一處不會出現(xiàn)“某組參數(shù)忘了設、用了默認值”的隱蔽錯誤。我見過最典型的一次是 vmax 忘了傳結構體里沒有這個字段函數(shù)在cfg.vmax處直接報錯——這還算好的報錯至少讓你知道如果代碼里寫了getfield(cfg,vmax,10)給了默認值那就變成靜默跑完、結論全錯。歸檔目錄我會在腳本開頭if ~exist(runs,dir), mkdir(runs); end別指望手動建。4.3 批量繪圖一眼看出哪組參數(shù)崩了繪圖的目的不是好看是讓“哪組參數(shù)崩了”這件事在 0.5 秒內可判斷。做法是把所有組的誤差曲線畫在同一張半對數(shù)圖上收斂的是直線下降振蕩的是鋸齒發(fā)散的直接沖頂。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 看起來差不多看不出收斂速率差異對數(shù)坐標下斜率直接就是收斂速度三組參數(shù)誰快誰慢一目了然。另外DisplayName一定要寫不然 legend 出來的是一串data1、data2回頭對不上。如果某條曲線在中段出現(xiàn)平臺然后再次下降那通常是飽和起了作用前期被 vmax 卡住誤差下降變慢進入線性區(qū)后恢復指數(shù)收斂。這個“拐點”位置本身是有信息的它告訴你限幅是不是給緊了。5. 排查收斂曲線不對時先查這五件事排查順序很重要。我的習慣是先看矩陣再看數(shù)值最后才懷疑協(xié)議——因為協(xié)議是最不容易寫錯的部分而矩陣和步長是最容易在復制粘貼中出問題的部分。5.1 矩陣與拓撲層面的三個坑現(xiàn)象一曲線看著在收斂但每個智能體停在不同的值上。原因是圖不連通或者本該無向的圖被寫成了有向且非強連通。無向連通圖的一致性值是全體初始狀態(tài)的算術平均一旦分裂成兩個連通分量就變成兩個各自的平均值。解決辦法是跑一遍連通性檢查numel(unique(conncomp(graph(A)))) 1不連通就回到邊表找漏掉的邊?，F(xiàn)象二前期正常十幾步之后突然發(fā)散且越跑越陡。步長超過穩(wěn)定邊界。查法很直接h和2/max(real(eig(L)))比一下超了就減半。我踩過一次更隱蔽的版本——代碼里 H 是大寫的局部變量、h 是小寫的全局參數(shù)循環(huán)里用的是 H檢查時改的是 h改了半天沒反應血淚經驗是變量名別玩這種字母大小寫的玄學?，F(xiàn)象三改了 N 之后報維度錯誤或者不報錯但隊形完全亂掉。硬編碼索引沒跟著改。典型位置是 A 的邊表、P 的行數(shù)、x0 的長度以及eye(N)里那個 N。解決方式是全腳本搜一遍數(shù)字常量把與規(guī)模相關的都換成size(...)或者從 cfg 里取。5.2 數(shù)值與工程層面的兩個坑現(xiàn)象四曲線里出現(xiàn) NaN從某一步開始整條線消失。零除或溢出。最可能是初始位置重合導致距離為零或者限幅里vmax / n的 n 是 0。解決辦法是給所有分母加max(x, eps)同時在循環(huán)里加一句assert(all(isfinite(u(:))), step %d 出現(xiàn) NaN, k)——早報錯比跑完再看圖省時間得多?，F(xiàn)象五腳本跑和 Simulink 跑結果對不上形狀相似但數(shù)值差一截。求解器步長和采樣步長不一致。Simulink 默認變步長求解器在飽和切換點會自動縮步腳本里是固定步長兩者必然有差異。辦法是把 Simulink 里的求解器改成固定步長步長填成和腳本一致再對比。如果還差就查模塊里的采樣時間有沒有設成繼承-1繼承來的采樣時間經常和你想的不是一個值。6. 進階替換點把周期通信換成事件觸發(fā)附一份驗收清單周期通信的假設是“每個采樣點所有智能體都要和鄰居交換一次信息”在仿真里不要錢在實機上就是通信帶寬和功耗。事件觸發(fā)一致性把這條假設松掉只在狀態(tài)偏差超過閾值時才通信其余時刻用保持器維持上一次的值。代碼上的改動量小得出奇大約十幾行但收益在規(guī)模大、通信受限的場景里非常明顯。6.1 事件觸發(fā)一致性的替換點核心是把“當前狀態(tài)”換成“上次廣播的狀態(tài)”再維護一個觸發(fā)誤差。觸發(fā)條件通常是 ‖e_i‖ σ‖z_i‖e_i 是當前狀態(tài)與上次廣播狀態(tài)之差z_i 是鄰居狀態(tài)的加權和σ 是觸發(fā)閾值越小越接近周期通信。function [X, nTrig] et_consensus(X0, L, h, K, sigma) N size(X0, 1); X zeros([size(X0), K1]); X(:,:,1) X0; Xhat X0; % 上次廣播出去的狀態(tài) nTrig zeros(N, 1); % 每個智能體的觸發(fā)次數(shù) for k 1:K Xc X(:,:,k); z -L * Xhat; % 用廣播狀態(tài)算控制量 e Xc - Xhat; % 觸發(fā)誤差 % 逐智能體判斷是否觸發(fā)觸發(fā)者更新廣播值 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那事件觸發(fā)形同虛設你還是每步都在用最新狀態(tài)。need是一個 N×1 的邏輯向量用邏輯索引整行更新比 for 循環(huán)快也比循環(huán)不容易寫錯下標。nTrig返回每個智能體的觸發(fā)次數(shù)這個數(shù)字很能說明問題如果某個智能體的觸發(fā)次數(shù)接近 K說明它一直在被鄰居甩開問題在拓撲不在閾值。調參上有兩條經驗。sigma從 0.1 開始試太大會出現(xiàn)“長時間不通信、隊形明顯滯后”的現(xiàn)象太小則退化成周期通信、省不下什么。另外事件觸發(fā)會引入額外的收斂時間同樣 T 下誤差可能比周期通信大一個量級比較兩者時要在相同觸發(fā)次數(shù)下比而不是相同時間內比——這一點在寫結論時經常被忽略。6.2 一套我自己在用的驗收清單代碼改完之后跑一遍下面這幾項比反復看曲線靠譜。這套清單是從多次“明明仿真收斂、換組參數(shù)就崩”的經歷里攢出來的。檢查項判據不通過時先看哪鄰接矩陣對稱性isequal(A, A)為真邊表是否只寫了單向拉普拉斯行和每行絕對值 1e-12度矩陣是否按行求和代數(shù)連通度 λ?明顯大于 0 1e-3圖是否連通、有無孤立點步長裕度h 0.1 × 2/λ_max拓撲切換時取所有模式的公共上限有限性全程無 NaN / Inf分母是否加了 eps收斂值無向連通時等于初始均值是否有向圖但按無向算飽和觸頂率大部分時間不貼 vmaxvmax 給緊了或拓撲把某個體甩開觸發(fā)占比事件觸發(fā)遠小于 1且各體差異不大sigma 是否過小、某體是否孤立表格里最容易被跳過的是“飽和觸頂率”但它往往是最早暴露拓撲問題的指標。我一般會在仿真函數(shù)里順手記一條satRatio mean(n vmax)跑完看這個數(shù)超過 0.5 就說明整段時間都在限幅區(qū)控制器根本沒按線性協(xié)議工作這時候去調增益全是白費。還有一處建議把 λ?、h 上限、飽和觸頂率、末態(tài)誤差這四個數(shù)打在一行日志里輸出。曲線可以騙人尤其是只看末態(tài)的時候這四個標量不太會。我自從被一次“圖不連通但曲線看起來收斂”坑掉整整兩天之后每次改完拓撲做的第一件事都是先跑那三條斷言再畫圖。斷言不過就不看曲線因為那時候的曲線只會給你錯誤的信息。希望幫到你。本文還有配套的精品資源點擊獲取