2020年华为杯数学建模竞赛D题论文和代码

无人机集群协同对抗

摘           要：

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

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

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

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

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

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

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

1. 问题重述

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

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

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

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

对象	参数	约束
对抗区域	攻击纵深BC(AD)	L=50km
	通道宽度AB(DC)	M 为限定约束
红无	运载机 速度VH	VH=300m/s, 方向依转弯半径改变
人机	最小转弯半径RH	RH =1000m
集群	无人机 速度 Vp	Vp=200m/s, 方向依转弯半径改变
	最小转弯半径 Rp	Rp=350m
蓝无人机集群	速度VE	VE=250m/s,方向依转弯半径改变
	最小转弯半径RE	RE=500m
飞行边界	红无人机集群	无边界限制
	蓝无人机集群	不能越过AD 和BC
		(1) 每个集群5架无人机初始位置近似分布在一个圆周上，
		任意相邻2架无人机的间距相同；
	红无人机	(2) 任意2 架无人机的间距大于30m；
	200m。	(3) 每一架无人机与本集群中至少2架无人机的距离不超过
编队	10km;	(1) 与所属无人机集群中至少一架无人机的距离不超过
	红运载机	(2) 与任何一架无人机的距离需大于100m； (3) 与蓝方无人机距离大于5km；
		(4) 每架携带10架无人机， 每波次发射不少于 3架无人机。
	蓝无人机	任意两架无人机的间距大于30m
攻防结果判定	蓝方突防成功	360s内越过边界， 且不被拦截
	红方拦截成功	蓝方突防无人机与红方至少2架无人机的距离均小于300m

如图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；

B) 红方 5 架运载机在约束下从边界CD一侧开始拦截，同时发射第一波次的无人机集群，运载机和无人机集群中心的距离为2km，分两波次发射无人机集群；

(C) 蓝方希望尽可能多的无人机突防成功；

(D) 红方希望尽可能多的拦截蓝方无人机。

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

2. 问题假设及符号说明

2.1 问题假设

1. 假设红、蓝双方的运动只考虑对抗平面上的运动，不考虑俯仰、滚转运动，飞行速度大小不变，方向可根据策略需要改变；

2. 假设运载机在满足约束条件下瞬间布设好无人集群为设定队形，攻防时在约束下转变队形；

3. 假设红、蓝双方无人机的体积和质量分布不影响运动和对抗结果判定，即无人机可视为质点进行分析；

4. 假设红、蓝双方无人机集群在初始进入区域时即探测到对方位置；

5. 假设红、蓝双方攻防结果只有蓝方突防成功和红方拦截成功两种判定情况；

6. 假设红方无人机拦截成功后，被拦截的蓝方无人机退出攻防，实施拦截动作的红方无人机位置不变继续参加攻防。

2.2 符号说明

符号	含义
P₁(E₁)	红(蓝)方第i个无人机
Hi	红方第i个运载机
VP₁(VE₁)	红(蓝)方第i个无人机速度
VHi	红方第i个运载机速度
RP₁(REi)	红(蓝)方第i个无人机转弯半径
RHi	红方第i个运载机转弯半径
tEi	蓝方第i个无人机的突防时间
lEi	蓝方第i个无人机的突防路径长度
SPk	红方第k个集群有效拦截面积
Pk	红方第k个集群队形的尖端点无人机
D	红方集群有效拦截面的最大宽度
M	通道宽度(突防带宽)
L	攻击纵深
N	蒙特卡洛试验次数

 

                α(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)

队形(C)的有效拦截面扩展后的最大宽度 Dc为：

                                   DC=1000m                              (10)

红色无人机群在进行拦截的过程中，保持以“贪婪”队形的最大拦截宽度对蓝方实施拦截，形成基于“贪婪”队形的红方协作拦截策略。

3.4基于内点法的蓝方突防求解模型

针对题目中蓝方无人机突防的约束条件具有较强非线性，甚至体现非凸性及病态性，无人机集群协同合作运动以及不同集群间的对抗规划任务中，碰撞躲避约束条件与拦截条件的非线性耦合进一步增加了求解难度。

AMPL(A Mathematical Programming Language) 是一种描述并求解大规模复杂数学问题的建模语言，支持世界上大部分的求解器，主要的特点是对优化问题的数学表达式的简化，使得简明地可读地定义优化问题成为可能。在对蓝方突防过程进行非线性的动态规划时，我们基于AMPL的自动微分功能，实现求解最优突防策略的精确解。

3.4.1 基于非线性规划 NLP 模型

蓝方无人机进入攻防区域后，在约束(1)中的运动约束、边界约束和对抗约束下，不断接近目标边界，优化函数为公式(2)。对于每一个考虑的场景设置相应的性能指标函数J，完整的无人机集群协同对抗运动规划任务可描述为一下最优控制问题[1]：

优化目标：

 

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

3.6.1针对蓝方突防概率结果

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

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

图12中的区域P，Q即为问题需要求解的蓝方进入位置区域。

3.6.2针对蓝方最优突防策略

从推演中，可以看出，红方无人机机群先采取圆形队形向前突进，再利用拦截宽带最大的队形进行拦截。蓝方无人机为避免面对红方两个无人机机群的围堵，尽量不从 G1G2通道突围，而选择向AD(BC)边机动突围，对应的红方无人机机群G1(G2)也朝着AD(BC)边移动围堵，并尽可能接近蓝方无人机。以蓝方无人机与红方无人机机群的尖端点之间的相对距离作为评价标准，红方无人机机群 G1(G2)朝着蓝方无人机方向移动，而相应的另一红方无人机机群 G2(G1)也向着蓝色无人机方向移动，辅助协防，避免蓝方无人机从红方两无人机机群中间位置实现突防，即均朝着使目标函数减小的方向移动，这时，蓝方无人机要朝着避开红方无人机机群 G1(G2)的方向机动，即朝着使目标函数增大的方向移动，实现突防。但由于红方无人机机群中间距离相对较大，蓝方无人机也有从中间通道突防的机会，如图13 所示：

针对此问中，由于蓝方无人机的起始位置未知，我们采用两类蓝方无人机的起始位置进行了策略讨论。

(1) 图12 对应的蓝方无人机位于 G1G2 通道内任意位置起飞的红蓝双方示意图，由于红方无人机机群中间距离相对较大，有从中心突围的可能性，关于策略的表述见表2。

代码

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

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;

85.         if (~ fail   flag)

86.               success   flag( case   id) = 1;

87.              GenerateVideo task04( init  x, init  y, init   theta, init   phy);

88.          else

89.                success   flag( case   id) = 0;

90.          end

91. end

1.  %% Write0bstaclesForReducedNLP. m

2.  function  num   constraints= Write0bstaclesForReducedNLP( x, y, theta, Nv)

3.  global  params  

4.  Nfe= params  . Nfe;

5.  index= ones( Nfe, Nv- 1, Nv- 1).*- 99;

6.  num   conflict   vehicles= ones( Nfe, Nv- 1).*- 99;

7.  num  v2v   constraints= 0;

8.  for  ii= 1: Nfe

9.           for          jj= 1:( Nv- 1)

10.                counter= 0;

11.                for  kk=( jj+ 1): Nv

12.                  if (~ IsConfig1FarFromConfig2([ x( jj, ii), y( jj, ii), theta( jj, ii)],[ x( kk, ii), y( kk, ii), theta( kk, ii)]))

13.                            counter= counter+ 1;

14.                            num  v2v   constraints= num  v2v   constraints+ 1;

15.                            index( ii, jj, counter)= kk;

16.                      end

17.               end

18.               num   conflict   vehicles( ii, jj)= counter;

19.         end

20. end

21.

22. delete( 'IND V2V');

23. fid= fopen( 'IND V2V',  'w');

24. for  ii= 1: Nfe

25.          for  jj= 1:( Nv- 1)

26.                for  kk= 1:( Nv- 1)

27.                      fprintf( fid, '%g  %g  %g  %g\r\n', ii, jj, kk, index( ii, jj, kk));

28.               end

29.          end

30. end

31. fclose( fid);

32.

33. delete( 'NUM IND V2V');