【协同任务】基于matlab遗传算法考虑分配次序的多无人机协同任务分配【含Matlab源码 143期】

一、多无人机协同作业简介

0 引言
多架无人机组成无人机集群可以协同完成任务，是未来无人机的发展方向。组成无人机集群的多架无人机通过机间链路互相通信实现协作，可以迅速准确地执行路径规划、协同侦察、协同感知和协同攻击等复杂任务。
为实现无人机集群协作的诱人前景，国内外都积极开展了相关研究工作。美国方面，美国国防预先研究计划局(DARPA)于2015年推出“小精灵”项目，计划研制具备自组织和智能协同能力的无人机蜂群系统。美国防部战略能力办公室(SCO)2014年启动了“无人机蜂群”项目，旨在通过有人机空射“灰山鹑”微型无人机蜂群执行低空态势感知和干扰任务。美国海军研究局(ONR)于2015年公布了“低成本无人机蜂群”(LOCUST)项目，研发可快速连续发射的无人机蜂群，无人机之间利用近距离射频网络共享态势信息，协同执行掩护、攻击或防御任务。2017年，在 DARPA 会议中心举办“进攻性集群战术”(OFFSET)项目的提案人活动，目标是发展基于游戏的开放架构，为城市作战的无人集群系统生成、评估和集成集群战术。
欧洲方面，2016 年，欧洲防务局启动了“欧洲蜂群”项目，开展了无人机蜂群的自主决策、协同飞行等关键技术研究。2016年，英国国防部发起无人机蜂群竞赛，参赛的多个团队控制无人机蜂群实现了通信中继、协同干扰、目标跟踪定位和区域测绘等任务。2017 年，俄罗斯无线电电子技术集团对外发表研究计划称，在战斗机上装载多架蜂群无人机可实现协同侦察和攻击的新型作战样式。
国内也相继展开相关研究。最近，中国电科(CETC)电子科学研究院发布了陆军协同无人机“蜂群”视频，引起广泛关注。
面对这一重要课题，本文研究总结了无人机协同应用的发展趋势，对其当前研究进展和发展方向进行了探讨，并提出无人机集群任务协作的发展趋势是多智能体协同。

1 无人机协同应用发展趋势
对现有研究工作的分析如图1所示，无人机协同应用大体上可以分为3个阶段，分别是分布式协同，群体智能协同和未来的多智能体协同。

图1 无人机协同技术的发展趋势
无人机集群协同发展的第1个阶段是简单分布式协同。在该阶段，协同任务根据执行条件，预先在简单连接和组合的集群成员之间经过计算处理后分配执行。集群基本没有能力根据环境与目标的变化动态调整任务分配，各无人机分担的任务通常是确定的。

鉴于预分配方式的局限性，受生物集群活动的启发，群体智能被应用于无人机集群，使无人机集群协同发展到第2个阶段——群体智能协同。在该阶段中，各个无人机节点被赋予初级智能，能够进行简单的认知和决策；通过集群个体之间更为紧密的耦合，可以根据执行中的反馈调整优化方式或者优化目标，使整个无人机集群系统有能力构成自组织、高稳定的分布式系统。群体智能协同阶段当前正处于研究和应用迅速发展时期。

随着节点计算能力的进一步提升和人工智能技术的飞速发展，无人机协同即将进入发展的第3个阶段——多智能体协同。在多智能体协同发展阶段，集群中的各个无人机都将是一个独立的综合智能体，具有多维度认知计算和高级智能处理能力，从而实现更高效的自主学习和决策，并在此基础上，完成更复杂、更艰巨的任务。

2 分布式协同
从无人机集群出现开始就被用于解决协同路径规划、协同感知和协同任务规划等分布式协作任务。早期的无人机分布式协同任务通常都是提前进行充分的计算和分配，无人机节点按照既定算法或者方案予以执行。根据计算结果，分布式协同无人机群可组织执行配置好的任务，如图2所示。

图2 分布式协同阶段的任务执行方式
2.1 协同路径规划
协同路径规划中，要求在给定目标或者搜索目标后，无人机节点根据当前的任务状态来确定飞行路径。针对协同搜索和跟踪任务的路线规划，无人机集群可采用最大化目标功能来检测最重要的目标，并在关键时刻对其进行跟踪，从而获得最有价值的信息。而协同搜索的路径规划可以分成无人机工作区间划分和全区域覆盖搜索路径规划2个子问题，将多机协同搜索转化为子区域上的单机搜索，对目标区域快速进行划分并生成飞行路线。基于改进遗传算法的多无人机协同侦察航迹规划算法，可用于解决面向复杂战场环境中高效侦察多种类型目标过程中的路径规划问题，并能够有效地提高航迹规划精度和效率。

2.2 协同感知
协同感知是多无人机集群共同探测感知某一目标区域状态的任务形式。在这类任务中最常见的是协同频谱感知。针对协同频谱感知任务特点，采用最佳融合准则的分布式协同任务执行方案可以优化检测性能，使协作频谱感知总错误率达到最小，并能降低协作感知时间，节省感知过程开销。

2.3 协同任务规划
协同任务规划要求集群系统能够根据目标任务和执行情况，对任务进行系统分配。例如，针对协同打击任务，通过建立无人机毁伤代价指标函数、航程代价指标函数和价值收益指标函数，可以实现多无人机协同打击任务的分配[5];而通过建立多目标优化模型并采用遗传算法，能有效提高任务的完成效率。针对协同搜救任务，使用一种新的基于通信保持的拍卖方法的自适应反馈调节遗传算法，能够改善传统遗传算法存在易陷入局部最优的弱点。

通过以上3种类型协同任务的相关研究工作可以看出，无人机集群的分布式协同方式虽然充分考虑了“分布式”特点，能够根据任务目标和集群特征，设置有效的目标函数和优化方法去寻求最优或者较优的结果，但是其任务执行环境和求解目标需要在任务执行之前进行优化计算，然后分配执行，不能适应实际中动态的任务目标和环境变化，缺乏“智能性”的感知和适应性行为。随着人们对“蜂群”“鸟群”等生物群体智能研究的深入，群体智能协同被进一步引入到无人机的协同中。

3 群体智能协同
“鸟群”“蚁群”等生物群体，虽然其中的个体智能有限，但是群体却展现出高度的自组织性，这一特点与无人机集群自主协同的需求相符，因而群体智能在无人机协同应用领域也得到了广泛研究，使无人机集群协同具有了初步的智能性。具备群体智能的无人机集群系统在任务拆解和执行过程中引入了群体反馈和适变能力，可执行较复杂的动态任务，其过程如图3所示。

图3 群体智能协同阶段的任务执行方式
3.1 协同路径规划
山区路径规划任务受到地形特征的强烈影响，难以预先确定性分配执行，简单的分布式协同无法胜任。例如，为了执行山区应急物资运输的路径规划任务，一种考虑路径安全度的改进蚁群算法可快速收敛并生成更短路径。蚁群算法也被用于无人机协同飞抵空战场的航迹规划任务中;一种基于改进混沌蚁群算法更能克服传统蚁群智能算法中易出现局部极值、收敛效率低的缺陷，提升算法的全局寻优能力和搜索效率。针对协同攻击移动目标场景的航迹规划，另一种改进的蚁群算法建立了结合任务分配的无人机群协同航迹规划模型，可以快速地对地面多个移动目标规划出有效的航迹。

3.2 协同感知
在协同感知任务中，群体智能也得到应用。针对无人机集群通信场景和需求，出现了考虑结合认知无线电技术指导下的智能通信思想和汇聚有限智慧的群体智能理论方法，构建了群体智能协同通信模型和智能协同感知模型。

3.3 协同任务规划
协同作战是协同任务规划中的典型场景，结合群体智能优化算法的优势，基于粒子群-整数编码狼群算法的集群组网任务分配算法适合解决此类协同问题;由于无人机集群协同决策困难，还可以结合狼群算法的认知与协作能力，实现在复杂环境下迅速对目标进行跟踪和包围。这种协同任务，是第一类简单协同所不能胜任的。

虽然无人机集群和群体智能的结合可以充分发挥无人机集群优势，增强分布式协同智能性，可在任务执行过程中与环境和任务执行中间过程产生一定交互和反馈，使之具备一定的自适应能力，但是，这种智能仍然是非常有限的，其本质仍然是基于特定计算模式和反馈模式下的分布式优化算法。

4 多智能体协同
随着人工智能技术和节点自身算力的不断增强，未来无人机集群中的个体将具备更强的智能性，能够独立对环境和任务进行感知和评估，实现多个智能体之间的交互和协同，从而具备多智能体协同能力。

近年来，人工智能领域研究取得了突破性进展。其中，深度强化学习在诸多领域得到了成功应用。无线通信网络基于多智能体深度强化学习的资源分配技术也得到了深入研究。多智能体深度强化学习模型早已被用于解决车联网中频谱资源分配问题,这种应用已经与无人机集群系统颇为接近。例如一种基于多智能体深度强化学习的分布式动态功率分配方案。基于多智能体深度强化学习的策略还可被用于二者的结合——用无人机辅助车辆网络进行多维资源管理。

虽然基于强化学习的多智能体通信网络资源分配问题已经得到了广泛研究，但是由于网络特性的不同，传统的针对其他通信网络的研究成果不能直接用于无人机集群网络。因而基于强化学习的多智能体自主协同应用逐渐成为未来无人机多智能体自主协同的一个研究热点。针对无人机群的通信网络资源动态分配问题而提出的多智能体深度强化学习方案也陆续出现,例如，一种基于多智能体深度强化学习的分布式干扰协调策略被用于受到干扰的无人机网络中的文件下载业务。智能体在适用于无人机网络特点的独立强化学习中，其行为策略通常只能根据它们对全局环境的局部个体观察来制定。针对这种局限性，联合采用2种不同规模的智能体可解决智能体之间的通信问题。

无人机集群协同，经常会处理动态高维离散和连续动作状态空间的优化求解问题，近来出现的演员-评论家算法是深度强化学习的一个新兴方向，结合了基于值函数和基于策略函数的深度强化学习两大分支的优势，非常适用于无人机集群的智能协同。利用演员-评论家算法，在无线信道和可再生能源再生率都是随机变化，且环境动态变化条件下，可寻求资源分配最佳策略,如用于解决车联网中复杂动态环境下的资源分配问题。应用设备到设备(D2D)网络的异构蜂窝网络环境下，基于演员-评论家算法的策略可用于智能化节能模式选择和资源分配。

随着节点智能算力的不断增强，无人机集群中的每架无人机可以作为一个具有深度强化学习能力的智能体，而整个集群可通过合作构成多智能体。相邻无人机之间通过通信网络进行信息交换与分发。如图4所示，每架无人机都与局部环境相互作用，根据从周围环境或者同伴无人机得到的信息，针对承载的任务需求，通过深度强化学习，智能地产生动作策略，进行自身资源与行为的分配与调整，进而与环境和同伴互动，并获得个体奖励。

图4 基于多智能体的无人机集群
每架无人机的深度强化学习智能体由2个深度神经网络构成，包括演员网络和评论家网络，如图5所示。

图5 基于演员评论家算法的无人机智能体
演员网络负责输出动作，评论家网络负责评价演员的动作，以获得相互促进的效果。与传统的深度强化学习方法相比，演员-评论家算法同时吸取了基于值函数方法和基于策略函数方法的优点，从价值和策略两方面来训练提升智能体，训练的更快，效果也更好。通过训练和学习，期望智能体的评论家网络可以获得最佳效用评估函数：

智能体从周围环境中观察得到状态信息St ,例如干扰状态、相邻无人机等。演员网络拟合动作策略函数π(St;ωπ),根据状态信息St,动作策略函数输出当前时隙的动作at,也就是资源分配结果，并应用到环境中得到个体的即时奖励rt。奖励由奖励函数计算得到，负责反馈一个即时的、合理的、具有指导意义的奖励值，从而激励智能体向目标更新策略。评论家网络拟合效用评估函数Q(St,at,ω),负责预测和评估在当前状态St采用动作at所能得到的状态动作价值，即长期性回报Rt为：

式中，γ是折扣因子，γ表示了当前对于未来收益的重视程度，取值在0～1之间，取值0则表示智能体忽视未来收益，只看重当前的收益rt,取值1则表示智能体将未来每个时刻的收益都视为与当前时刻收益一样的重要性。获得最大的长期性回报是智能体的最终目的，这个目标可以根据任务的性质被定义为不同的评判标准。

评论家网络的最佳效用评估函数对应的演员网络动作策略即为最佳动作策略π*。动作策略函数的参数通过策略梯度的方法进行更新，效用评估函数的参数采用最小化损失函数进行更新：

式中，yt为优势函数，用来衡量计算执行动作后的实际效用评估值。通过更新，智能体的演员网络的动作策略输出的动作越来越好，而评论家网络的效用评估也会越来越准确。基于多智能体的无人机群的每架无人机都会朝着收益最大的方向执行动作，从而实现集群的收益最大化。

二、遗传算法简介

首先通过分析UAV分配次序对打击任务总收益的影响, 设计了动态战场环境的更新规则. 将航程代价和任务代价作为惩罚项修正目标函数, 建立了考虑分配次序的UAVs协同目标分配优化模型. 然后针对模型的物理意义改进了遗传算法基因编码方式, 设计了MUCTA遗传算法. 该算法利用状态转移思想, 引进SDR算子获得多种分配次序种群, 同时以单行变异算子修正UAV与目标对应关系, 并采用最优个体法和轮盘赌法筛选子代个体.

1 遗传算法概述
遗传算法（Genetic Algorithm，GA）是进化计算的一部分，是模拟达尔文的遗传选择和自然淘汰的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。该算法简单、通用，鲁棒性强，适于并行处理。

2 遗传算法的特点和应用
遗传算法是一类可用于复杂系统优化的具有鲁棒性的搜索算法，与传统的优化算法相比，具有以下特点：
（1）以决策变量的编码作为运算对象。传统的优化算法往往直接利用决策变量的实际值本身来进行优化计算，但遗传算法是使用决策变量的某种形式的编码作为运算对象。这种对决策变量的编码处理方式，使得我们在优化计算中可借鉴生物学中染色体和基因等概念，可以模仿自然界中生物的遗传和进化激励，也可以很方便地应用遗传操作算子。
（2）直接以适应度作为搜索信息。传统的优化算法不仅需要利用目标函数值，而且搜索过程往往受目标函数的连续性约束，有可能还需要满足“目标函数的导数必须存在”的要求以确定搜索方向。遗传算法仅使用由目标函数值变换来的适应度函数值就可确定进一步的搜索范围，无需目标函数的导数值等其他辅助信息。直接利用目标函数值或个体适应度值也可以将搜索范围集中到适应度较高部分的搜索空间中，从而提高搜索效率。
（3）使用多个点的搜索信息，具有隐含并行性。传统的优化算法往往是从解空间的一个初始点开始最优解的迭代搜索过程。单个点所提供的搜索信息不多，所以搜索效率不高，还有可能陷入局部最优解而停滞；遗传算法从由很多个体组成的初始种群开始最优解的搜索过程，而不是从单个个体开始搜索。对初始群体进行的、选择、交叉、变异等运算，产生出新一代群体，其中包括了许多群体信息。这些信息可以避免搜索一些不必要的点，从而避免陷入局部最优，逐步逼近全局最优解。
（4）使用概率搜索而非确定性规则。传统的优化算法往往使用确定性的搜索方法，一个搜索点到另一个搜索点的转移有确定的转移方向和转移关系，这种确定性可能使得搜索达不到最优店，限制了算法的应用范围。遗传算法是一种自适应搜索技术，其选择、交叉、变异等运算都是以一种概率方式进行的，增加了搜索过程的灵活性，而且能以较大概率收敛于最优解，具有较好的全局优化求解能力。但，交叉概率、变异概率等参数也会影响算法的搜索结果和搜索效率，所以如何选择遗传算法的参数在其应用中是一个比较重要的问题。
综上，由于遗传算法的整体搜索策略和优化搜索方式在计算时不依赖于梯度信息或其他辅助知识，只需要求解影响搜索方向的目标函数和相应的适应度函数，所以遗传算法提供了一种求解复杂系统问题的通用框架。它不依赖于问题的具体领域，对问题的种类有很强的鲁棒性，所以广泛应用于各种领域，包括：函数优化、组合优化生产调度问题、自动控制
、机器人学、图像处理（图像恢复、图像边缘特征提取…)、人工生命、遗传编程、机器学习。

3 遗传算法的基本流程及实现技术
基本遗传算法（Simple Genetic Algorithms,SGA）只使用选择算子、交叉算子和变异算子这三种遗传算子，进化过程简单，是其他遗传算法的基础。

3.1 遗传算法的基本流程
通过随机方式产生若干由确定长度（长度与待求解问题的精度有关）编码的初始群体；
通过适应度函数对每个个体进行评价，选择适应度值高的个体参与遗传操作，适应度低的个体被淘汰；
经遗传操作（复制、交叉、变异）的个体集合形成新一代种群，直到满足停止准则（进化代数GEN>=?）；
将后代中变现最好的个体作为遗传算法的执行结果。

其中，GEN是当前代数；M是种群规模，i代表种群数量。

3.2 遗传算法的实现技术
基本遗传算法（SGA）由编码、适应度函数、遗传算子（选择、交叉、变异）及运行参数组成。
3.2.1 编码
（1）二进制编码
二进制编码的字符串长度与问题所求解的精度有关。需要保证所求解空间内的每一个个体都可以被编码。
优点：编、解码操作简单，遗传、交叉便于实现
缺点：长度大
（2）其他编码方法
格雷码、浮点数编码、符号编码、多参数编码等
3.2.2 适应度函数
适应度函数要有效反映每一个染色体与问题的最优解染色体之间的差距。
3.2.3选择算子

3.2.4 交叉算子
交叉运算是指对两个相互配对的染色体按某种方式相互交换其部分基因，从而形成两个新的个体；交叉运算是遗传算法区别于其他进化算法的重要特征，是产生新个体的主要方法。在交叉之前需要将群体中的个体进行配对，一般采取随机配对原则。
常用的交叉方式：
单点交叉
双点交叉（多点交叉，交叉点数越多，个体的结构被破坏的可能性越大，一般不采用多点交叉的方式）
均匀交叉
算术交叉
3.2.5 变异算子
遗传算法中的变异运算是指将个体染色体编码串中的某些基因座上的基因值用该基因座的其他等位基因来替换，从而形成一个新的个体。

就遗传算法运算过程中产生新个体的能力方面来说，交叉运算是产生新个体的主要方法，它决定了遗传算法的全局搜索能力；而变异运算只是产生新个体的辅助方法，但也是必不可少的一个运算步骤，它决定了遗传算法的局部搜索能力。交叉算子与变异算子的共同配合完成了其对搜索空间的全局搜索和局部搜索，从而使遗传算法能以良好的搜索性能完成最优化问题的寻优过程。

3.2.6 运行参数

4 遗传算法的基本原理
4.1 模式定理

4.2 积木块假设
具有低阶、定义长度短，且适应度值高于群体平均适应度值的模式称为基因块或积木块。
积木块假设：个体的基因块通过选择、交叉、变异等遗传算子的作用，能够相互拼接在一起，形成适应度更高的个体编码串。
积木块假设说明了用遗传算法求解各类问题的基本思想，即通过积木块直接相互拼接在一起能够产生更好的解。

三、部分源代码

%本实验目的是验证改进的DE离散目标分配算法的有效性
%实验设置：各种环境，调用目标分配的算法
%目标分配实验的输入:代价矩阵
%利用代价矩阵和映射方法，在离散和连续空间转换
%目标分配实验的输出:有效的基因染色体表示
%比较实验：
%1、可行性分析，三种模型
%2、进化策略对比分析，证明采用双策略的有效性
%3、大数据分析
%4、与其他方法的比较
% 本文只用到了AssignType =2;情形，即UAV数量大于TARGET数量
%使用了改动的遗传算法
%相比较源程序，较大改动部分为：对象执行区（GetFit1和Getfit2两个函数），和差分进化区（全部改动），其他部分有少量改动，我不记得了。很多部分直接删改，可能注释没有对应的上，强烈建议对比源代码学习观看
%本人编程能力较差，编写格式不规范，很多地方注释不足以及随心所欲各种直接删改，不同实验直接改数据，以及出现的重复，各种FOR循环嵌套等等，我自己都觉得乱七八糟。再次强烈建议对比源代码
%建议只使用该代码框架（即阅读源代码即可，本代码可适当看一下就行，估计很费劲），自行删改其中内容，写自己的文章，但注意源代码自行备份，自己的代码尽量多注释以及写实验说明，方便后来的师弟或师妹。%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%清空环境
clear;%定义地形对象
cTerrain = CTerrain3D;
%定义威胁区对象
cTread = CTread;%调整为统一的图形环境
figure('Renderer','opengl');%绘制地形和威胁区cTerrain.draw();
cTread.draw();%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%数据形式说明
%xs 初始点位置集合，三维点集
%yt 目标点位置集合，三维点集
%MxsEna 能力矩阵，集合中的二元组表示最大最小飞行速度
%MxsDIs 最大航程矩阵，集合中的值表示UAV最大可飞的距离
%MytTOrd 任务点间的时序坐标，要求尽量在前面的先执行%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%调用目标分配算法
%AssignType表示不同的模型，1, N=M; 2, N>M; 3, N<M
%Key表示算法执行次数, 0，执行一次，1，执行多次
AssignType =2;
Key =0;%时间测试
Time = cputime;%%    第一种情况的实验
if AssignType == 1%%%%%%%%%%%%%%%%%%%%%%%数据区%%%%%%%%%%%%xs = [%11 2 3;2 17 5; 34 26 7;5 2 2;6 33 3;27 34 6;21 45 5;25 12 3;22 17 11;13 31 12;%45 21 11; 36 22 11;13 21 16;41 25 10;25 25 11;53 12 14;22 12 12;18 32 11; 67 34 12;26 16 15;78 20 15;93 31 12;31 20 13;112 32 15; 134 26 17;45 52 12;36 63 13;67 34 16;21 85 15;32 62 13;];yt = [%40 210 12;170 90 13 ; 80 197 21 ;172 120 13;160 56 13;160 143 21;170 200 21;113 200 12;97 134 16;100 145 11;%81 101 21; 72 152 12; 133 60 11;84 164 21; 150 90 21;146 121 12; 167 54 12; 108 165 12;99 120 21;110 143 11;70 181 12;15 133 13 ; 19 151 12 ;160 192 13;101 120 23;160 113 12;82 101 11;173 140 12;182 65 14; 128 156 12;];  MxsEna=[%0.2 0.3; 0.2 0.4;0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3; 0.2 0.4;% 0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3; 0.2 0.4;0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.3 0.5;0.3 0.6;];MxsDis=[%400 700 650 500 700 900 450 610 400 700%650 500 700 900 450 610 700 900 450 610500 700 300 350 700 900 450 610 450 610];ytW = [%1 3 4 2 1 1 3 2 1 2%3 2 1 3 2 1 2 3 2 1 1 1 1 1 1 1 1 1 1 1];MytTOrd=[3 4; 5 2; 6 8; 7 4;];%引入新的协同矩阵，最大起始时间矩阵Twait = [];Twindow=[];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if Key == 0%%%%%%%%%%%%%%%%%%%%%%算法执行区%%%%%%%%%%%%%%%%%%%调用目标分配算法，并取得解和相关数据,计算算法耗费时间%定义目标分配算法对象DDE1 = DMDEAssignTarget(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd,Twait,Twindow,1);% MC = Cal1.GetMcost();%Cal1 = CallocationTargetOld1(cTerrain);%标签显示文字xlabel('X/km')ylabel('Y/km')zlabel('Z/km')   %DDEAssign1=DDEAssignTarget1(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd);%执行run并返回结果集t1 = cputime;[solU solT solC solF ] = DDE1.run();Time = cputime -t1;%         %结果solUsolTsolCsolFTimeelse if Key == 1for i=1:10DDE1 = DMDEAssignTarget(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd,Twait,Twindow,i);% MC = Cal1.GetMcost();%Cal1 = CallocationTargetOld1(cTerrain);%标签显示文字xlabel('X/km')ylabel('Y/km')zlabel('Z/km')   %DDEAssign1=DDEAssignTarget1(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd);%执行run并返回结果集t1 = cputime;[solU solT solC solF ] = DDE1.run();solu(i,:) = solU;solt(i,:) = solT;solc(i,:) = solC;solf(i,:) = solF;Time(i,:) = cputime -t1;end %i
%         %结果solusoltsolcsolfTimeend % key=1end % key=0else if AssignType == 2%% 第二种实验的情况%%%%%%%%%%%%%%%%%%%%%%%数据区%%%%%%%%%%%%xs = [%11 2 3;2 17 5; 34 26 7;5 2 2;6 33 3;27 34 6;21 45 5;25 12 3;22 17 11;13 31 12;%45 21 11; 36 22 11;13 21 16;41 25 10;25 25 11;53 12 14;22 12 12;18 32 11; 67 34 12;26 16 15;78 40 15;93 31 12;31 25 13;112 32 15; 134 26 17;45 52 12;36 55 13;67 34 16;21 45 15;160 30 13;];yt = [%40 210 12;170 90 13 ; 80 197 21 ;172 120 13;160 56 13;160 143 21;170 200 21;113 200 12;97 134 16;100 145 11;%81 101 21; 72 152 12; 133 60 11;84 164 21; 150 90 21;146 121 12; 167 54 12; 108 165 12;99 120 21;110 143 11;70 151 12;40 155 13 ; 19 151 12 ;150 140 13 ];   MxsEna=[%0.2 0.3; 0.2 0.4;0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3; 0.2 0.4;% 0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3; 0.2 0.4;0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.3 0.5;0.3 0.6;];MxsDis=[%400 700 650 500 700 900 450 610 400 700%650 500 700 900 450 610 700 900 450 610500 700 300 350 700 900 450 610 450 610];ytW = [%1 3 4 2 1 1 3 2 1 2%3 2 1 3 2 1 2 3 2 1 1 1 1 1];%这个数据里都代表的是目标点，是目标点与目标点的关系MytTOrd=[3 2];%引入新的协同矩阵，最大起始时间矩阵Twait = [30 40 20 15 50 70 30 80 90 65];Twindow=[];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if Key == 0%%%%%%%%%%%%%%%%%%%%%%算法执行区%%%%%%%%%%%%%%%%%%%调用目标分配算法，并取得解和相关数据,计算算法耗费时间%定义目标分配算法对象DDE1 = DMDEAssignTarget(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd,Twait,Twindow,1);% MC = Cal1.GetMcost();%Cal1 = CallocationTargetOld1(cTerrain);%标签显示文字xlabel('X/km')ylabel('Y/km')zlabel('Z/km')   %DDEAssign1=DDEAssignTarget1(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd);%执行run并返回结果集[solU solT solC solF   cpop ]  = DDE1.run();Time
cpop;
%         %结果solUsolTsolCsolFTimeelse if Key == 1Time2=0;cpop3=0;for i=1:20DDE1 = DMDEAssignTarget(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd,Twait,Twindow,i);% MC = Cal1.GetMcost();%Cal1 = CallocationTargetOld1(cTerrain);%标签显示文字xlabel('X/km')ylabel('Y/km')zlabel('Z/km')   %DDEAssign1=DDEAssignTarget1(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd);%执行run并返回结果集t1 = cputime;[solU solT solC solF  cpop ] = DDE1.run();solu(i,:) = solU;solt(i,:) = solT;solc(i,:) = solC;solf(i,:) = solF;Time(i,:) = cputime -t1;
%         genX(i,:)=cpop(:,2);cpop2(:,i)=cpop(:,2);Time(i)end %i
%         %结果%             figure(3);
%             grid on;
%             hold on;
%             plot(genX(1:length(genX),1)',genX(1:length(genX),2)','linewidth',2);
%          solu
%          solt
%          solc
%          solfTime=Time2/20;cpop=cpop2/20;end % key=1end % key=0else if AssignType == 3%% 第三种实验的情况%%%%%%%%%%%%%%%%%%%%%%%数据区%%%%%%%%%%%%xs = [%11 2 3;2 17 5; 34 26 7;5 2 2;6 33 3;27 34 6;21 45 5;25 12 3;22 17 11;13 31 12;%45 21 11; 36 22 11;13 21 16;41 25 10;25 25 11;53 12 14;22 12 12;18 32 11; 67 34 12;26 16 15;78 20 15;93 31 12;31 20 13;112 32 15];yt = [%40 210 12;170 90 13 ; 80 197 21 ;172 120 13;160 56 13;160 143 21;170 200 21;113 200 12;97 134 16;100 145 11;%81 101 21; 72 152 12; 133 60 11;84 164 21; 150 90 21;146 121 12; 167 54 12; 108 165 12;99 120 21;110 143 11;70 131 12;15 125 13 ; 19 151 12 ;142 120 13;120 100 23;160 90 12;82 101 11;173 140 12;160 160 14; 45 155 12;];    MxsEna=[%0.2 0.3; 0.2 0.4;0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3; 0.2 0.4;% 0.4 0.75;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.3;0.35 0.45;0.3 0.5;0.3 0.6;0.2 0.5; 0.3 0.4;0.4 0.75;0.3 0.45];MxsDis=[%400 700 650 500 700 900 450 610 400 700%650 500 700 900 450 610 700 900 450 610700 700 900 650];ytW = [%1 3 4 2 1 1 3 2 1 2%3 2 1 3 2 1 2 3 2 1 1 1 1 1 1 1 1 1 1 1];MytTOrd=[3 4; 5 2];%引入新的协同矩阵，最大起始时间矩阵Twait = [30 40 20 15];Twindow=[1300,6000];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if Key == 0%%%%%%%%%%%%%%%%%%%%%%算法执行区%%%%%%%%%%%%%%%%%%%调用目标分配算法，并取得解和相关数据,计算算法耗费时间%定义目标分配算法对象DDE1 = DMDEAssignTarget(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd,Twait,Twindow,1);% MC = Cal1.GetMcost();%Cal1 = CallocationTargetOld1(cTerrain);%标签显示文字xlabel('X/km')ylabel('Y/km')zlabel('Z/km')   %DDEAssign1=DDEAssignTarget1(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd);%执行run并返回结果集t1 = cputime;[solU solT solC solF]  = DDE1.run();Time = cputime -t1;%         %结果solUsolTsolCsolFTimeelse if Key == 1for i=1:20DDE1 = DMDEAssignTarget(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd,Twait,Twindow,i);% MC = Cal1.GetMcost();%Cal1 = CallocationTargetOld1(cTerrain);%标签显示文字xlabel('X/km')ylabel('Y/km')zlabel('Z/km')   %DDEAssign1=DDEAssignTarget1(cTerrain,xs,yt,MxsEna,MxsDis,ytW,MytTOrd);%执行run并返回结果集t1 = cputime;[solU solT solC solF ] = DDE1.run();solu(i,:) = solU;solt(i,:) = solT;solc(i,:) = solC;solf(i,:) = solF;Time(i,:) = cputime -t1;end %i
%         %结果solusoltsolcsolfTimeend % key=1end % key=0end % type =3end % type =2end % type=1

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
[3]巫茜,罗金彪,顾晓群,曾青.基于改进PSO的无人机三维航迹规划优化算法[J].兵器装备工程学报. 2021,42(08)
[4]邓叶,姜香菊.基于改进人工势场法的四旋翼无人机航迹规划算法[J].传感器与微系统. 2021,40(07)
[5]马云红,张恒,齐乐融,贺建良.基于改进A*算法的三维无人机路径规划[J].电光与控制. 2019,26(10)
[6]焦阳.基于改进蚁群算法的无人机三维路径规划研究[J].舰船电子工程. 2019,39(03)