【多智能体】二阶多智能体系统的最优时不变分布式编队跟踪Matlab复现

✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

👇 关注我领取海量matlab电子书和数学建模资料

🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。

🔥 内容介绍

本文探讨了最优时不变编队跟踪问题，旨在为具有二阶积分器动力学特性的多智能体系统提供分布式解决方案。在相关文献中，大多数关于多智能体编队跟踪的研究在探究分布式反馈控制律时并未考虑能量问题。为解决这一关键设计要点，我们通过在优化成本中精心选择特定且关键的势函数，对一个涵盖轨迹跟踪、基于距离的编队控制以及输入能量最小化的优化问题进行了形式化处理并给出解决方案，以此做出贡献。为此，我们展示了如何借助基于投影算子的轨迹优化牛顿法（PRONTO）以集中式方式计算逆动力学。更重要的是，我们将这种离线解决方案作为通用参考，来设计一种稳定的在线分布式控制律。最后，给出了涉及三维空间中沿直线路径的立方体编队的数值示例，以验证所提出的控制策略。

⛳️ 运行结果

📣 部分代码

close allclear all, clc​global dt nAg M kr ka QBp QBdp R kF kA Xdes dijs G flag33 flag50 flag90​dt = 0.01;tl = 200;T0 = (0:dt:tl)';T = length(T0);​nAg = 4;I_nAg = eye(nAg);M = 2;I_M = eye(M);​QBp = 10^1;QBdp = 10^1;R = 10^0;kF = 10^0;kA = 5*10^-1;kr = 1;ka = 50;​% initial conditionsvel = 1;x0p = [-2 1 -3 -1 2 -2 0 0];x0p = x0p/(2*norm(x0p));x0p = rand(1,8)*10;x0v = 5*[0 -vel 0 -vel 0 -vel 0 -vel];​% desired trajectory for the centroidxd0p = zeros(1,M);xd0v = 1*[vel 0];Xdesp = repmat(xd0p,length(T0),1) + T0*xd0v;Xdesv = repmat(xd0v,length(T0),1);Xdes = [ Xdesp Xdesv ];​% framework definitiondist_coeff = 1.05;dd = 5;% dijs = [0 dd dd;% dd 0 dd;% dd dd 0];dd1_2 = sqrt(2)*dd;max_dijs = dd1_2;INF = eps+dist_coeff*dd1_2;dijs = [0 dd INF dd;dd 0 dd dd1_2;INF dd 0 dd;dd dd1_2 dd 0];​A = zeros(nAg);for i = 1:nAgfor j = 1:i-1if dijs(i,j) > 0A(i,j) = 1;A(j,i) = 1;endendendG = graph(A);degs = sum(A);D = diag(degs);L = D-A;L_ = D^(-1/2)*L*D^(-1/2);eigsL_ = sort(eig(L_));lambda1 = eigsL_(2);lambdan_1 = eigsL_(end);muL_ = (lambda1+lambdan_1)/2;eta_star = 1-1/muL_;if eta_star < 0eta_star = 0;end% if topology remains constant:F = eye(nAg)*eta_star+(D^-1*A)*(1-eta_star);​​%% -------------------------------------------------------% parametersparams.dt = dt;params.nAg = nAg;params.M = M;params.QBp = QBp;params.QBdp = QBdp;params.R = R;params.kF = kF;params.kA = kA;params.Xdes = Xdes;params.dijs = dijs;params.G = G;params.tl = tl;params.F = kron(F,I_M);params.dist_coeff = dist_coeff;params.max_dist = max_dijs; % max(dijs(:));​% flags to show completionflag33 = 0;flag50 = 0;flag90 = 0;​% dynamics integration[t,x_history] = ode45(@(t,x)distr_dyn(t,x,params),T0',[x0p x0v]');​​%% display final resultsfinal_state = x_history(end,:)';disp(final_state)p1 = x_history(end,1:2)'p2 = x_history(end,3:4)'p3 = x_history(end,5:6)'p4 = x_history(end,7:8)'e12 = p1-p2;e13 = p1-p3;e23 = p2-p3;e14 = p1-p4;e24 = p2-p4;e34 = p3-p4;pB = (p1+p2+p3+p4)/4Ne12 = norm(e12)Ne13 = norm(e13)Ne23 = norm(e23)Ne14 = norm(e14)Ne24 = norm(e24)Ne34 = norm(e34)​​%% fase diagram for positions​​​figuregrid onhold onplot(Xdes(:,1),Xdes(:,2),'k')for i = 1:nAgplot(x_history(:,(i-1)*M+1),x_history(:,(i-1)*M+2),'r')endfor i = 1:nAgplot(x_history(end,(i-1)*M+1),x_history(end,(i-1)*M+2),'b*','linewidth',3)end% inserire qui un confronto tra grafici... prendere i dati altrui potrebbe% essere una buona idea per poi plottare la traiettoria distribuita su% quella centralizzata -> fare un esempio con triangolo equilatero% - traiettoria% - funzionale di costo% - input speso% - ripartizione tra formation vs tracking% * consensus solo per il distribuito% nota: mettere il peso finale per la final boundary condition nel caso% centralizzato!! (peso finale appropriato, no un punto a caso per la% traiettoria desiderata)axis equal​%% evolutions of the distance errors (sigma_dij-s)​figuregrid onhold ondist_errors = zeros(T,nAg*(nAg-1)/2);​k = 0;for i = 2:nAgfor j = 1:i-1dij = dijs(i,j);if dij >= 0for tt = 1:Tp_i_t = x_history(tt,(i-1)*M+1:(i-1)*M+M);p_j_t = x_history(tt,(j-1)*M+1:(j-1)*M+M);sij_t = norm(p_i_t-p_j_t)^2;dist_errors(tt,k+1) = sigma(sij_t,dij,0);endelsedist_errors(:,k+1) = -ones(T,1);endk = k+1;endend​k = 0;for i = 2:nAgfor j = 1:i-1tratto = '-';if dijs(i,j) < 0tratto = '--';endplot(T0,dist_errors(:,k+1),tratto,'linewidth',1.5)k = k+1;endend​title('zero-order consensus')xlabel('$t$','interpreter','latex')ylabel('$\sigma(s_{ij})\qquad$','interpreter','latex')set(gca,'fontsize',25)set(get(gca,'ylabel'),'rotation',0)​​%% evolutions of the first derivatives of the distance errors (sigma_dij-s)​figuregrid onhold onddist_errors = zeros(T,nAg*(nAg-1)/2);​k = 0;for i = 2:nAgfor j = 1:i-1dij = dijs(i,j);if dij >= 0for tt = 1:Tp_i_t = x_history(tt,(i-1)*M+1:(i-1)*M+M);p_j_t = x_history(tt,(j-1)*M+1:(j-1)*M+M);sij_t = norm(p_i_t-p_j_t)^2;ddist_errors(tt,k+1) = sigma(sij_t,dij,1);endelseddist_errors(:,k+1) = zeros(T,1);endk = k+1;endend​k = 0;for i = 2:nAgfor j = 1:i-1tratto = '-';if dijs(i,j) < 0tratto = '--';endplot(T0,ddist_errors(:,k+1),tratto,'linewidth',1.5)k = k+1;endend​title('first-order consensus')xlabel('$t$','interpreter','latex')ylabel('$\frac{\partial\sigma}{\partial s_{ij}}\qquad$','interpreter','latex')set(gca,'fontsize',25)set(get(gca,'ylabel'),'rotation',0)​%% evolutions of the second derivatives of the distance errors (sigma_dij-s)​figuregrid onhold ondddist_errors = zeros(T,nAg*(nAg-1)/2);​k = 0;for i = 2:nAgfor j = 1:i-1dij = dijs(i,j);if dij >= 0for tt = 1:Tp_i_t = x_history(tt,(i-1)*M+1:(i-1)*M+M);p_j_t = x_history(tt,(j-1)*M+1:(j-1)*M+M);sij_t = norm(p_i_t-p_j_t)^2;dddist_errors(tt,k+1) = sigma(sij_t,dij,2);endelsedddist_errors(:,k+1) = zeros(T,1);endk = k+1;endend​k = 0;for i = 2:nAgfor j = 1:i-1tratto = '-';if dijs(i,j) < 0tratto = '--';endplot(T0,dddist_errors(:,k+1),tratto,'linewidth',1.5)k = k+1;endend​title('second-order consensus')xlabel('$t$','interpreter','latex')ylabel('$\frac{\partial^2\sigma}{\partial s_{ij}^2}\qquad$','interpreter','latex')set(gca,'fontsize',25)set(get(gca,'ylabel'),'rotation',0)

🔗 参考文献

🍅往期回顾扫扫下方二维码