无人机集群协同对抗

摘           要：

本文针对非线性约束条件下红蓝双方无人机集群协同对抗的最优规划问题，结合贪婪队形、非线性规划、内点法、蒙特卡洛方法和全联立正交配置有限元法，构建了无人机集群协同对抗推演模型。

针对问题一，构建了基于蒙特卡洛法的蓝方突防概率推演模型。我们综合考虑了基于贪婪队形的红方协作拦截策略，和基于内点法的蓝方突防求解模型。首先分析红方在集群约束下的贪婪队形和拦截策略，然后我们基于内点法联合MATLAB 和 AMPL 编程求解蓝方突防的非线性规划解，最后基于蒙特卡洛法对整个蓝方突防过程进行推演，得到突防概率图。通过推演，我们提出了区域P和Q两块突防成功率最高的区域作为所求区域，并对蓝方的最优突防策略进行分析。

针对问题二，构建了基于正则化目标函数的蓝方突防策略推演模型。我们考虑通过改进第一问的模型进行求解，首先针对所求时间最短突防策略，基于正则化思想改进优化的目标函数。然后在固定的通道宽度下对新的约束和目标函数进行攻防策略推演。针对所求带宽下限，对每一轮推演使用的带宽变化量进行渐进松弛，不断逼近所求下限。获得的带宽下限 Mmin为91km，并对不同条件下的蓝方最短时间突防策略进行分析。

针对问题三，构建了基于全联立正交配置有限元法的非线性求解优化模型。综合考虑多重非线性约束条件和红方多波次释放无人机的复杂环境问题，我们通过全联立正交配置有限元法将问题中的所有状态量离散化，提高非线性优化模型的推演效率。然后在固定的通道宽度下对红蓝双方的攻防过程进行推演。针对所求带宽上限，对每一轮推演使用的带宽变化量进行渐进紧缩，不断逼近所求上限。获得的带宽上限 Mmax为39km，并对不同条件下的红方最优拦截策略进行分析。

针对问题四，构建了基于 NLP 的红蓝双方攻防策略推演模型。综合考虑前述的多重非线性约束和红蓝双方无人机集群博弈的复杂环境，我们通过前述建立的红蓝双方攻防策略推演模型，对红蓝无人机集群的协同突防和协同拦截状态进行推演，最后对不同条件下的红方最优拦截策略和蓝方最优突防策略进行分析。

整体来说，红蓝双方的最优攻防策略均以红蓝双方无人机集群的距离为评价标准，针对蓝方突防，尽量选择远离红方无人机的方向突防； 针对红方拦截，尽量选择朝向蓝方无人机的移动方向拦截、

关键词：非线性约束 最优规划 集群协同 蒙特卡洛方法 全联立正交配置有限元法

1. 问题重述

人工智能技术的不断发展催生了无人战车、无人机、无人艇等新型作战力量，同时也推动着作战方式从单体对抗的智能化到群体协同对抗的智能化方向发展。在复杂的战场环境下，通过多架次无人机集群的协同探测、协同进攻、协同拦截和协同躲避，能有效提高任务的完成效果。

本题即研究在一定战场环境约束和飞行状态约束条件下，攻守无人机集群的最优协同突防和协同拦截策略问题。

攻防过程中，红蓝双方无人机需要满足的约束如表1所示。

表 1 红蓝双方集群攻防约束表

如图1所示，红蓝双方在平面矩形区域ABCD 内进行无人机集群的协同对抗，蓝色方为进攻方，进攻出发位置根据需求变化，进攻目标位置为边界CD； 红色方为防守方，从边界CD 出发对蓝方进行拦截。

在满足表1所示约束的情况下，针对以下问题对红蓝双方在攻防区域内的攻防策略进行建模分析：

问题1. 对抗伊始，在约束条件下，针对设定：

(A) 红方 2 个无人机集群分别以边界CD上的G₁和G₂为中心的圆周队形进入攻防区域, 圆周半径为100m, DG₁=20km,G₁G₂=30km,CG₂=20km,,即通道带宽M=70km;

(B) 蓝方无人机可从矩形区域ABCD内任意位置突防。

分析蓝方总能采用合适的策略躲避拦截，成功突防的进入位置，并讨论相应的最优突防策略。

问题 2. 对抗伊始，在约束条件下，针对设定：

(A) 蓝方无人机边界 AB的中心点进入攻防区域；

(B) 红方2个无人机集群分别以G₁和G₂为中心的圆周队形进入攻防区域，圆周半径为100m, G₁和G₂位于边界CD上, 具体位置未知。

分析通道带宽M的下限 Mmin的存在性，满足当M比Mmin大时，蓝方无人机一定能突防成功，并给出该情形下蓝方无人机时间最短的突防策略。

问题3. 在约束条件下，通道带宽M=70km，针对设定：

(A) 红方每架运载机分两波次共发射 10 架无人机，发射的无人机集群初始队形满足表1中的编队约束，运载机与圆周中心的距离为2km，集群队形可调整；

对抗伊始，

(B) 蓝方无人机从边界 AB 中心进入；

(C) 红方两架运载机分别从边界CD上G₁和G₂进入，同时发射第一波次的无人机集群，运载机和无人机集群中心具体位置未知；

(D) 运载机在保证与第一波次集群满足约束条件下，发射第二波次无人机集群。

讨论红方的最优拦截策略，包括：两架运载机两个波次发射的无人机数量、每架运载机第二波次发射的时刻和位置以及第二波次发射的无人机集群的中心位置。

分析通道带宽M的上限 Mmax的存在性，满足当M 小于 Mmax E时，红方无人机一定能拦截成功。

问题4. 对抗伊始，在约束条件下，针对设定：

(A) 蓝方突防集群的3架无人机从边界AB一侧开始突防，通道带宽M=100km；

                α(E)           蓝方无人机突防成功率

算式·的范数

3. 问题一模型的建立与求解

3.1针对问题一的分析

在无人机协同对抗的环境中，无人机单体可以相互通信、控制，避免碰撞的同时还可以构成集群，进行一定的编队组合，从而最大化躲避、拦截等决策的收益。满足表1中编队约束的队形示意图如图2 所示：

问题一固定了整个对抗平面的攻击纵深L，通道宽度M，红方运载机释放无人机的位置点G₁和G₂，要求我们设计模型讨论满足约束时蓝方无人机从对抗平面内的哪些位置进入时，总能成功突防，并且讨论蓝方无人机的最优突防策略。

在这种情景下，蓝方总能找到成功突防的策略进行突防，可以转化为计算在不同位置进入时蓝方突防的成功率α(E)。为了找到蓝方从攻防区域的何处进入才能使蓝方突防成功率α(E)最高，针对多约束条件下的动态博弈策略问题，我们参考“有限元”的思想，构建基于蒙特卡洛的蓝方突防概率模型，解决非线性约束的动态规划问题，利用模拟出的不同位置下蓝方突防成功率对最优突防位置点进行分析。

针对每一个位置点进入后的红蓝攻防动态决策过程，参考以下约束进行非线性动态规划：

 s.t.VP=200m/s,VH=300m/s,VE=250m/s,

 Rp≥300m,RH≥1000m,RE≥500m,

 M=70000m,tE<360s,

i,j=1,2,3,4,5,j≠i,m=2,3,4,I=1,2,..., ki,

 Pkm-Pkm-1=Pkm-Pkm+1,PI-PJ≥30m,

                       |Hk-PI||>100m,Hk-E1>5000m,                         (1)

 if‖Hₖ-Pₖₗ‖>10000m,

 then error,

 if||Pki-Pbj||>200morPki-Pkj1||>200m,

 then error

其中，ki为红方第k个集群内的第i个无人机，kj₁为红方第k个集群内除第i个无人机外的某一个标号为j₁的无人机

蓝方的优化目标函数为：

                         maxdk,P-lE=E1-Pk*-lE,k=1,2                         (2)

红方的优化目标函数为：

                            mindE,P=Pk*-E1,k=1,2                             (3)

根据这些目标和约束，我们可以对动态过程进行模拟，得到使得蓝方突防成功率最高的区域。

3.2攻防区域的数值化

蓝方在某一确定位置进入攻防区域后，可以视为一种分阶段递进的动态决策过程。首先红方基于探测到的蓝方位置进行拦截决策，在飞向蓝方进行拦截的过程中通过变更队形以达到最大化拦截面积的效果，这可以提高动态对抗中红方拦截的成功率。而蓝方可以通过对红方这一决策的探测进行有针对性的动态避障。我们可以利用数值方法对红蓝双方的动态对抗过程进行规划模拟。建立二维直角坐标系对整个区域进行数值化分析的步骤如下：

Step1. 选取坐标原点，建立二维平面直角坐标系.

Step2. 依据空间尺度关系和运动约束构建目标物的动态模型.

Step3. 对目标物的动态模型进行递进规划，完成数值模拟.

根据上述过程，我们可以对整个攻防区域ABCD 进行数值化。

如图3所示，以点B为原点，以边界BC为x轴，以边界BA为y轴建立平面直角坐标系，可以假设红方运载机 FY01及其发射集群内的无人机以45°的夹角进入对抗区域，运载机FY02及其发射集群内无人机的进入方向与FY01集群的方向关于边界CD的中垂线对称。

红方无人机集群的运动状态模型可以用P表示，红方运载机的运动状态模型可以用H表示，蓝方无人机的运动状态模型可以用E表示：

                                    P=VPRptpT                              (4)

                                    H=VHRHtHT                             (5)

                                    E=VEREtET                              (6)

针对各集群的状态量可以进行递进的数值模拟。

3.3基于贪婪队形的红方协作拦截策略

3.3.1基于队形的红方无人机集群有效拦截面扩展

初始时刻，红方运载机在释放出无人机集群后，无人机集群首先形成一个半径为 100m的圆形队形，此时红方集群需要实时探测红蓝方之间的距离( dE,P,进行拦截决策。此时红蓝双方的位置示意图如图4所示：

由于红方集群形成了某种拓扑结构，不容易对距离 dE,P进行表示，但双方的相对位置已知，并且对每一个红方无人机来说，其周围半径300m的圆形区域是蓝方无人机的危险区域，所以我们在红方队形拓扑结构的基础上，以每架红无人机为圆心， 以300m为半径，可以得到约束下的红方无人机集群对蓝方无人机的有效拦截面，以初始时刻的圆周队形为例，如图5中所示，这样可以有效表示出红方的拦截效果。

在红方有效拦截面基础上进行扩展，使用合适的外接几何结构可以得到方便描述集群外轮廓的尖端点P²，在所有的尖端点中，我们选取与蓝方无人机方向更接近的尖端点，代入距离dEp作为的近似。这样的选取方法既能保证红方拦截方向的局部最优，又能简化动态规划过程。

此时可以定义距离(dE. p:

                                  dE,P=E1-Pk*||2,k=1,2                         (7)

3.3.2基于贪婪队形的集群协作拦截

按上小节方法选取合适的红蓝方距离dEp后，红方集群将对优化函数(3)进行优化以靠近蓝方无人机，实施拦截。同时，在靠近红方无人机还会不断变更队形来提高拦截概率，主要体现在通过变换队形保持较大的有效拦截面积SPₖ。为了达到最大的拦截概率，我们基于的等边三角形、圆形和正五边形的内外接关系设计了下面几种“贪婪”队形，如图 6所示：

在图6的队形中，P⁺为其中的一个尖端点。红方无人机均以“贪婪条件”，即最大约束条件编组队形，相邻两架无人机的距离d₀都为200m，这样可以在约束条件下尽可能的扩大有效拦截面。以每一种队形都在运动约束下以最宽的拦截距离向蓝方靠近。队形(A)为圆形，有效拦截面扩展后的最大宽度DA为绿色扩展圆的直径，为：

                             DA=50+1055⋅d0≈875.28m                        (8)

同理可以得到，队形(B)的有效拦截面扩展后的最大宽度DB为：

                            DB=32⋅d0+300⋅2≈946.41m                        (9)

3.6蓝方突防概率结果与最优突防策略分析

3.6.1针对蓝方突防概率结果

采用蒙特卡洛随机洒点的方法，结合非线性动态规划得到了蓝方无人机突防情况散点图，如图 10 所示，其中，红色点代表蓝方无人机突防成功的点，蓝点代表蓝方无人机突防失败的点。

根据蓝方无人机突防情况散点图，我们绘制了蓝方无人机从不同位置进入的情况下的三维突防概率图和二维突防概率地形图，如图11，图12所示：

. %%  《无人机集群协同对抗》实现源码

2. %% main 主函数

3. clc; clear; close all;

4. % 基础参数设置 注：为方便计算与作图，本文采用 Km, Km/s,s作为距离，速度，时间的度量单位

5. global params   bandwidth

6. %基于贪婪队形的红方协作拦截策略 队形(A)/(B)/(C)

7. Formation   Red = 0;                % (A):0 / (B):1 / (C):2

8. if Formation   Red == 0

9.        params  . maxwidth = 0.87528;

10. elseif Formation   Red == 1

11.       params  . maxwidth = 0.94641;

12. elseif Formation   Red == 2

13.       params  . maxwidth = 1.0000;

14. end

15. params  . vE  max = 0.25;

16. params  . vP  max = 0.20;

17. params  . vT  max = 0.30;

18. params  . a  max = 0;

19. params  . phy   max = 0;

20. params  . w  max = 0;

21. params  . x  min = 0;

22. params  . x  max = 50;

23. params  . y  min = 0;

24. params  . y  max = 70;

25. params  . x  scale = params  . x  max - params  . x  min;

26. params  . y  scale = params  . y  max - params  . y  min;

27. params  . tf = 1;

28. params  . num   nodes  x = 1500;

29. params  . num   nodes  y = 1500;

30. params  . num   nodes  t = 200;

31. params  . resolution  t = params  . tf / ( params  . num   nodes  t - 1);

32. params  . resolution  x = params  . x  scale / ( params  . num   nodes  x - 1);

33. params  . resolution  y = params  . y  scale / ( params  . num   nodes  y - 1);

34. params  . multiplier   for  H = 2.0;

35. params  . weight   for   time = 5.0;

36. params  . max   cycle = 1000;

37. params  . n  large   value = 10000000;

38. params  . Rcos   ang = params  . radius .* cos( linspace(0, 2 * pi, 30));

39. params  . Rsin   ang = params  . radius .* sin( linspace(0, 2 * pi, 30));

40. params  . Nfe = 100;

41. params  . norm  L  threshold = 1;

42. params  . L  threshold = params  . norm  L  threshold;

43. params  . alpha = 0.5;

44. params  . Iter   max = 10;

45. bandwidth = 100;

46. num   cases = 5000;

47. success   flag = zeros(1, num   cases);

48. cpu   time = zeros(1, num   cases);

49. xv   library = 50* rand(1, num   cases);

50. yv   library = 70* rand(1, num   cases);

51. theta0  library = -1.6+(3.2)* rand(1, num   cases);

52. yvtf   library = abs(50* rand(1, num   cases));

53. thetatf   library = -1.6+(3.2)* rand(1, num   cases);

54.

55. for case   id = 1: num   cases

56.         clc;

57.         str=[' Iterations =' num2str( case   id)];

58.         disp( str);

59.        str1=[' Success flag =' num2str( success   flag)];

60.         disp(str1);

61.         case   str = ['BenchmarkCases\', num2str(1)];

62.         global boundary   configs   obstacles   original   obstacle   layers  

63.         load( case   str);

64.         tic;

65.        WriteBoundaryValuesAndBasicParams( Nv, Nobs);

66.      [ init  x, init  y, init   theta, init  v, init  a, init   phy, init  w] = GenerateInitialGue ss( Nv);

67.         iter = 0; fail   flag = 0;

68.        while (~IsCollisionFree( init  x, init  y, init   theta))

69.              if ( iter > params  . Iter   max)

70.                      fail   flag = 1;

71.                    break;

72.              end

73.             Write0bstaclesForReducedNLP( init  x, init  y, init   theta, Nv);

74.               ! ampl rr. run

75.              load opti   flag. txt

76.              if (~ opti   flag)

77.                    params  . L  threshold = params  . L  threshold * params  . alpha;

78.              else

79.                   [ init  x, init  y, init   theta, init   phy] = UpdateInitialGuess( Nv);

80.                    params  . L  threshold = params  . norm  L  threshold;

81.               end

82.              iter = iter + 1;

83.         end

84.         cpu   time( case   id) = toc;