1 简介
四旋翼无人机具有成本较低、设备简单、飞行时间灵活等特点，近些年被广泛应用于军事和民用领域，如目标侦察、应急救援、农业植保、无人机灯光表演。随着任务复杂度的增加，单架无人机往往难以满足任务需求，因此无人机集群控制及其应用由此成为目前的研究热点，多无人机集群能够提高执行任务效率和灵活度。无人机队形变换控制方法是实现多无人机编队飞行的前提，集群无人机队形重构问题是我们要考虑的一个重要问题，让每架无人机都能从初始位置无碰撞的到达最终位置，进而保证队形重构过程中代价最小或能耗最优。其中目标分配问题最多利用匈牙利算法进行解决，但是在多无人机轨迹规划时普遍存在计算难度大、规划时间增长、规划效率难以满足实际需求的问题。因此，探索计算简便、效率高的多无人机路径规划算法是目前迫切需要的。
（1）实验描述
十九架无人机组成的编队如图2.1所示，由F队形切换至Z队形，如果忽略本身以及外在约束条件的情况下，将会有多种不同的移动方案。但实际情形下无人机飞行距离越长，耗能也就越多，部分无人机到达目标的距离过长，消耗电力过多，从而提前降落。
图2.1： 编队实景图
使用传统匈牙利算法来解决队形变换的目标分配问题，其编队切换仿真图如下图，其中绿色方块表示无人机在F队形中的位置，红色圆点表示Z队形中无人机集群所要到达的各目标航点，黑色直线表示经分配后无人机由初始位置到达目标位置的直线路径。
图2.2： 匈牙利算法路径图
由上图可以看出，部分无人机到达目的地的直线路径过长，而大部分无人机初始位置与目标位置重合耗能少。移动路径长短不一，造成无人机耗能相差较大。所以我们需要寻求一种移动方式，使各个无人机移动路径长度趋于一致。
（2）实验要求
1、改进分配方案，使无人机在编队切换过程中飞行路径相近，飞行至目标航点时间一致，避免无人机耗能相差较大的问题。
2、集群无人机空中个别无人机能量不足或故障情况需退出，编队无人机数目发生变化，为维持编队队形需重新规划各无人机所在位置。
1 解决方案
集群无人机编队重构中考虑时间的一致性问题，可扩展研究切换过程中各无人机的距离约束，规划编队重构中各无人机路径长度，在同样移动速度下，使集群无人机在避免碰撞的情形下尽量同时到达目标点。通过考虑无人机编队切换的距离大小，若在编队切换中各个无人机移动的距离相近，那么在避免碰撞的情况下到达指定位置，使整个机群飞至目标航点时间趋于一致。
和匈牙利算法类似的是，我们同样使位于变换图形边缘的九个重合目标不发生移动，但其余十个目标均发生移动（匈牙利算法中除边缘九个点还存在两个重合点，详情可见图2）。通过一个10*10的代价矩阵（将移动距离设为代价），利用算法寻求出一个基于移动距离一致性为目标的优化分配方案。
2 部分代码
clear all
%C=[2 15 13 4;10 4 14 15;9 14 16 13;7 8 11 9];
symbol = 'bo'; % 打点颜色符号(b. 蓝点； bo蓝圈)
symbol1 = 'wo'; % 打点颜色符号(w. 白点； wo白圈)
dt = 1; % 采样步距
v = 1; % 速度
max_time = 1000; % 最大运行时间
% 两点之间距离
width = 10;
% 做 8*8 矩阵 19 个点 字符：F
% array_f = [0 1 1 1 1 1 1 0;
% 0 0 1 0 0 0 0 0;
% 0 0 1 0 0 1 0 0;
% 0 0 1 1 1 1 0 0;
% 0 0 0 0 0 1 0 0;
% 0 0 1 0 0 0 0 0;
% 0 0 1 0 0 0 0 0;
% 0 1 1 1 0 0 0 0];
%
% % 做 8*8 矩阵 19 个点 字符：Z
% array_z = [0 1 1 1 1 1 1 0;
% 0 1 0 0 0 1 0 0;
% 0 0 0 0 1 0 0 0;
% 0 0 0 1 0 0 0 0;
% 0 0 0 1 0 0 0 0;
% 0 0 1 0 0 0 0 0;
% 0 0 0 0 0 1 0 0;
% 1 1 1 1 1 1 0 0];
%
% T0=zeros(8,3);
% T1=zeros(8,3);
% m=1;n=1;
% for i=1:8 %% 去除共同点，得到无共同点的F、Z坐标集
% for j=1:8
% if array_f(i,j) ~= array_z(i,j) && array_f(i,j) == 1
% T0(m,:) = [i,j,0];
% m = m+1;
% elseif array_f(i,j) ~= array_z(i,j) && array_z(i,j) == 1
% T1(n,:) = [i,j,0];
% n = n+1;
% elseif array_z(i,j) == array_f(i,j) && array_z(i,j)==1 % 如果为共同点，将矩阵中标记 1变为2
% array_z(i,j) = 2;
% array_f(i,j) = 2;
% end
% end
% end
% array_f
% array_z
% T0
% T1
%------------------------------------------------------------------------
% 测试加中间重合点
array_f =[0 2 2 2 2 2 2 0;
0 0 1 0 0 0 0 0;
0 0 1 0 0 1 0 0;
0 0 1 1 1 1 0 0;
0 0 0 0 0 1 0 0;
0 0 1 0 0 0 0 0;
0 0 1 0 0 0 0 0;
0 2 2 2 0 0 0 0;];
array_z =[0 2 2 2 2 2 2 0;
0 1 0 0 0 1 0 0;
0 0 0 0 1 0 0 0;
0 0 0 1 0 0 0 0;
0 0 0 1 0 0 0 0;
0 0 1 0 0 0 0 0;
0 0 0 0 0 1 0 0;
1 2 2 2 1 1 0 0;];
T0 =[2 3 0;
3 3 0;
3 6 0;
4 3 0;
4 4 0;
4 5 0;
4 6 0;
5 6 0;
6 3 0;
7 3 0;] ;
T1 =[2 2 0;
2 6 0;
3 5 0;
4 4 0;
5 4 0;
6 3 0;
7 6 0;
8 1 0;
8 5 0;
8 6 0;];
%--------------------------------------------------------------------------------------------
%disp(array_z);
t=length(T1);
C=ones(t,t);
for i=1:t
for j=1:t
a=(T0(i,1)-T1(j,1))^2;
b=(T0(i,2)-T1(j,2))^2;
c=(T0(i,3)-T1(j,3))^2;
d=sqrt(a+b+c);
C(j,i)=d;
end
end
c =C;
disp(c);
%----------------------------------------------------------------------------------------
% 进行点匹配
% [Matching,Cost]=Hung4_Average(C);
% for i=1:t
% for j=1:t
% if Matching(i,j)==1
% M(1,i) = C(i,j);
% end
% end
% end
%
% ave_M = mean(M);
%
[M,C_select,cost,Mean_square_erro] = zq_ave(C);
% 将匹配矩阵映射到一维数组中
disp('最优匹配为:');
disp(M);
disp('最优代价为:');
C_select
disp('最优解为:');
cost
disp('最优均方差为:');
Mean_square_erro
%---------------------------------------------------------------------------------------%
% 建模
%场景的范围
xmin = 0;xmax = width * 8 + 20;ymin = 0;ymax = width * 8 + 20;
axis([xmin xmax ymin ymax]); %设定坐标范围
figure(1);
hold on ; %保留绘图内容
% 初始化矩阵
id_sta_addr = zeros(19,2);
id_sto_addr = zeros(19,2);
id_cur_addr = zeros(19,2);
id_tm = zeros(19, 1); %保存移动时间
index = 1;
% 给飞机实时编号，行扫描 安置无人机初始位置 实际上就是按比例显示F
for i = 1: 8
for j = 1: 8
% start 坐标
if array_f(i, j) == 1
% 做矩阵点位置应该实际场景
[id_sta_addr(index, 1), id_sta_addr(index, 2)] = calc(i, j, width);
% 安置飞机
plot(id_sta_addr(index, 1), id_sta_addr(index, 2), 'bo');
index = index + 1;
end
end
end
for i = 1: 8
for j = 1: 8
% start 坐标
if array_f(i, j) == 2
% 做矩阵点位置应该实际场景
[id_sta_addr(index, 1), id_sta_addr(index, 2)] = calc(i,j, width);
% 安置飞机
plot(id_sta_addr(index, 1), id_sta_addr(index, 2), 'ro');
index = index + 1;
end
end
end
pause(2);
%%算法替代部分 这里只是模型验证
% 测试部分--------------------------------------------------
id_sto_addr = id_sta_addr;
% calc(2, 2, width); 中(2,2)表示无人机终止位置坐标，width表示位宽
%id_sto_addr(7, 1)； 中(8,1)表示8 表示 ID， 1 表示 X坐标， 2 表示Y坐标
% 表示 ID 7的点 移动到(2,2) 的位置
for i=1:t
p = M(1,i);
[id_sto_addr(i, 1) , id_sto_addr( i, 2)] =calc(T1(p,1), T1(p,2), width);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算 各个无人机运行时间
for i=1:19
id_tm(i)=sqrt((id_sto_addr(i, 1)-id_sta_addr(i, 1))^2 + (id_sto_addr(i, 2)-id_sta_addr(i, 2))^2) / v; %计算这个 Step 的移动时间
end
%获取最大运行时间
max_tm = max(id_tm);
%max_tm_id = find(max_tm);
%状态 是否已经到达指定位置
id_status = zeros(19,1);
% id scan
id_cur_addr = id_sta_addr;
for t=0:dt:max_time
for index = 1:19
% 单点移动
[id_cur_addr(index, 1), id_cur_addr(index, 2), id_status(index)] = move(index, ...
id_sta_addr, ...
id_sto_addr, ...
id_cur_addr, ...
id_status, ...
t, ...
id_tm, ...
v);
end
pause(0.1);
%判断是否为空矩阵
if isempty(find(id_status == 0))
break;
end
end
fprintf(' 运行时间为：%d s \r\n', t);
3 仿真结果
4 参考文献
[1] 张思宇. 多无人机协同航迹规划及其控制方法研究[D].北京理工大学,2016.
[2] 卢燕梅,宗群,张秀云,鲁瀚辰,张睿隆.集群无人机队形重构及虚拟仿真验证[J].航空学报,2020,41(04):248-259.
[3] 毛琼,李小民,王正军.基于规则的无人机编队队形构建与重构控制方法[J].系统工程与电子技术,2019,41(05):1118-1126.
[4] 赵凌楚.多旋翼无人机编队队形变换算法研究[J].信息通信,2020(03):57-59.
[5] 林倩玉. 多无人机协同编队控制算法研究[D].哈尔滨工业大学,2018.
[6] 吴宇,梁天骄.基于改进一致性算法的无人机编队控制[J/OL].航空学报:1-18[2020-06-26].http://kns.cnki.net/kcms/detail/11.1929.V.20200521.1555.012.html.
[7] 马培蓓,雷明,纪军.基于一致性的多无人机协同编队设计[J].战术导弹技术,2017(02):86-90.
[8] 陈杰敏,吴发林,耿澄浩,徐珊.四旋翼无人机一致性编队飞行控制方法[J].航空兵器,2017(06):25-31.
部分理论引用网络文献，若有侵权联系博主删除。
5 完整MATLAB代码与数据下载地址
见博客第一条
