基于梦境优化算法（Dream Optimization Algorithm, DOA）的多个无人机协同路径规划（可以自定义无人机数量及起始点），MATLAB代码_梦 优化算法

一、梦境优化算法（DOA）

算法原理

算法灵感

梦境优化算法（Dream Optimization Algorithm, DOA）是一种新型的元启发式算法（智能优化算法），其灵感来源于人类梦境的启发。在有做梦经历的快速眼动睡眠期间，低频脑电波的功率降低，而高频脑电波的功率增加，这表明在做梦经历期间大脑的神经兴奋更大。

初始化

与其他元启发式算法类似，在初始化阶段，DOA首先在搜索空间内生成一个随机种群作为初始种群，从而开始算法的优化过程。初始总体的计算公式如下：
Xi=Xl+rand×(Xu−Xl)X_i = X_l + rand \\times (X_u - X_l)Xi​=Xl​+rand×(Xu​−Xl​)
式中，$N$为个体数，即种群大小；XiX_iXi​是种群中的第$i$个个体；XlX_lXl​和XuX_uXu​分别表示搜索空间的下边界和上边界；$rand$是一个dim维向量，每个维度是0到1之间的随机数。

勘探阶段

在勘探阶段（迭代计数从0到TdT_dTd​），我们首先根据记忆能力的差异将群体分为5组，每组中的个体数量大致相等。每组中的个体根据其记忆能力的不同，以不同的方式更新其位置。具体来说，每组中的个体根据其记忆能力的强弱，以不同的概率在全局最优解和局部最优解之间进行搜索。这种分组策略有助于在算法的早期阶段扩大搜索范围，增加找到全局最优解的可能性。

开发阶段

在开发阶段（迭代计数从TdT_dTd​到Td+TeT_d + T_eTd​+Te​），算法逐渐减少对全局搜索的依赖，转而更加关注局部搜索。此时，种群中的个体根据其当前的位置和适应度值，以一定的概率在全局最优解和局部最优解之间进行搜索。这种策略有助于在算法的中期阶段平衡全局搜索和局部搜索，进一步提高算法的搜索效率和精度。

更新阶段

在更新阶段（迭代计数从Td+TeT_d + T_eTd​+Te​到TmaxT_{max}Tmax​），算法主要进行局部搜索，以进一步优化解的质量。此时，种群中的个体根据其当前的位置和适应度值，以较高的概率在局部最优解附近进行搜索。这种策略有助于在算法的后期阶段精细调整解，提高解的精度。

算法流程

1. 初始化：在搜索空间内随机生成初始种群。
2. 勘探阶段：
   • 根据记忆能力将群体分为5组。
   • 每组中的个体根据记忆能力的不同，以不同的概率在全局最优解和局部最优解之间进行搜索。
3. 开发阶段：
   • 种群中的个体根据当前的位置和适应度值，以一定的概率在全局最优解和局部最优解之间进行搜索。
4. 更新阶段：
   • 种群中的个体根据当前的位置和适应度值，以较高的概率在局部最优解附近进行搜索。
5. 终止条件：达到最大迭代次数TmaxT_{max}Tmax​或满足其他终止条件。

算法描述

梦境优化算法（DOA）通过模拟人类梦境中的记忆和遗忘过程，结合基本的记忆策略和遗忘补充策略，平衡探索和利用，从而在优化过程中有效地搜索全局最优解。该算法在不同的阶段采用不同的搜索策略，初期扩大搜索范围，中期平衡全局和局部搜索，后期精细调整解，具有较强的全局搜索能力和良好的收敛性能。
参考文献：
[1]Lang Y, Gao Y. Dream Optimization Algorithm (DOA): A novel metaheuristic optimization algorithm inspired by human dreams and its applications to real-world engineering problems[J]. Computer Methods in Applied Mechanics and Engineering, 2025, 436: 117718.
原文链接：https://blog.csdn.net/weixin_46204734/article/details/145667749

2. 无人机路径规划数学模型

2.1 路径最优性

为了提高无人机的操作效率，规划的路径需要在特定的应用标准下达到最优。在我们的研究中，主要关注空中摄影、测绘和表面检查，因此选择最小化路径长度作为优化目标。由于无人机通过地面控制站（GCS）进行控制，飞行路径 XiX_iXi​ 被表示为无人机需要飞越的一系列 $n$ 个航路点的列表。每个航路点对应于搜索地图中的一个路径节点，其坐标为 Pij=(xij,yij,zij)P_{ij} = (x_{ij}, y_{ij}, z_{ij})Pij​=(xij​,yij​,zij​)。通过表示两个节点之间的欧几里得距离为 $| \\overrightarrow{P_{ij}P_{i,j+1}} |，与路径长度相关的成本 F1F_1F1​ 可以计算为：

F1(X)=∑j=1n−1∥PijPi,j+1→∥F_1(X) = \\sum_{j=1}^{n-1} \\| \\overrightarrow{P_{ij}P_{i,j+1}} \\|F1​(X)=j=1∑n−1​∥Pij​Pi,j+1​​∥

2.2 安全性和可行性约束

除了最优性之外，规划的路径还需要确保无人机的安全操作，引导其避开操作空间中可能出现的威胁，这些威胁通常由障碍物引起。设 $K$ 为所有威胁的集合，每个威胁被假设为一个圆柱体，其投影的中心坐标为 CkC_kCk​，半径为 RkR_kRk​，如下图 所示。

对于给定的路径段 ∥PijPi,j+1→∥\\| \\overrightarrow{P_{ij}P_{i,j+1}} \\|∥Pij​Pi,j+1​​∥，其相关的威胁成本与它到 CkC_kCk​ 的距离 dkd_kdk​ 成正比。考虑到无人机的直径 $D$ 和到碰撞区域的危险距离 $S$，威胁成本 F2F_2F2​ 在障碍物集合 $K$ 上计算如下：

F2(Xi)=∑j=1n−1∑k=1KTk(PijPi,j+1→),F_2(X_i) = \\sum_{j=1}^{n-1} \\sum_{k=1}^K T_k(\\overrightarrow{P_{ij}P_{i,j+1}}),F2​(Xi​)=j=1∑n−1​k=1∑K​Tk​(Pij​Pi,j+1​​),

其中

Tk(PijPi,j+1→)={0,if dk>S+D+Rk(S+D+Rk)−dk,if D+Rk S + D + R_k \\\\(S + D + R_k) - d_k, & \\text{if } D + R_k < d_k \\leq S + D + R_k \\\\\\infty, & \\text{if } d_k \\leq D + R_k\\end{cases}Tk​(Pij​Pi,j+1​​)=⎩⎨⎧​0,(S+D+Rk​)−dk​,∞,​if dk​>S+D+Rk​if D+Rk​<dk​≤S+D+Rk​if dk​≤D+Rk​​

在操作过程中，飞行高度通常被限制在给定的最小和最大高度之间，例如在调查和搜索应用中，需要相机以特定的分辨率和视场收集视觉数据，从而限制飞行高度。设最小和最大高度分别为 hminh_{\\text{min}}hmin​ 和 hmaxh_{\\text{max}}hmax​。与航路点 PijP_{ij}Pij​ 相关的高度成本计算为：

Hij={∣hij−hmax+hmin2∣,if hmin≤hij≤hmax∞,otherwiseH_{ij} =\\begin{cases}|h_{ij} - \\frac{h_{\\text{max}} + h_{\\text{min}}}{2}|, & \\text{if } h_{\\text{min}} \\leq h_{ij} \\leq h_{\\text{max}} \\\\\\infty, & \\text{otherwise}\\end{cases}Hij​={∣hij​−2hmax​+hmin​​∣,∞,​if hmin​≤hij​≤hmax​otherwise​

其中 hijh_{ij}hij​ 表示相对于地面的飞行高度，如下图所示。

可以看出，HijH_{ij}Hij​ 保持平均高度并惩罚超出范围的值。对所有航路点求和得到高度成本：

F3(X)=∑j=1nHijF_3(X) = \\sum_{j=1}^n H_{ij}F3​(X)=j=1∑n​Hij​

平滑成本评估转弯率和爬升率，这对于生成可行路径至关重要。如下图 所示。

转弯角 ϕij\\phi_{ij}ϕij​ 是两个连续路径段 Pij′Pi,j+1′→\\overrightarrow{P\'_{ij}P\'_{i,j+1}}Pij′​Pi,j+1′​​ 和 Pi,j+1′Pi,j+2′→\\overrightarrow{P\'_{i,j+1}P\'_{i,j+2}}Pi,j+1′​Pi,j+2′​​ 在水平面 Oxy 上的投影之间的角度。设 k→\\overrightarrow{k}k 是 z 轴方向的单位向量，投影向量可以计算为：

Pij′Pi,j+1′→=k→×(PijPi,j+1→×k→)\\overrightarrow{P\'_{ij}P\'_{i,j+1}} = \\overrightarrow{k} \\times (\\overrightarrow{P_{ij}P_{i,j+1}} \\times \\overrightarrow{k})Pij′​Pi,j+1′​​=k×(Pij​Pi,j+1​​×k)

因此，转弯角计算为：

ϕij=arctan⁡(∥Pij′Pi,j+1′→×Pi,j+1′Pi,j+2′→∥PijPi,j+1′→⋅Pi,j+1′Pi,j+2′→)\\phi_{ij} = \\arctan\\left( \\frac{\\| \\overrightarrow{P\'_{ij}P\'_{i,j+1}} \\times \\overrightarrow{P\'_{i,j+1}P\'_{i,j+2}} \\|}{\\overrightarrow{P_{ij}P\'_{i,j+1}} \\cdot \\overrightarrow{P\'_{i,j+1}P\'_{i,j+2}}} \\right)ϕij​=arctan​Pij​Pi,j+1′​​⋅Pi,j+1′​Pi,j+2′​​∥Pij′​Pi,j+1′​​×Pi,j+1′​Pi,j+2′​​∥​​

爬升角 ψij\\psi_{ij}ψij​ 是路径段 PijPi,j+1→\\overrightarrow{P_{ij}P_{i,j+1}}Pij​Pi,j+1​​ 与其在水平面上的投影 Pij′Pi,j+1′→\\overrightarrow{P\'_{ij}P\'_{i,j+1}}Pij′​Pi,j+1′​​ 之间的角度，由下式给出：

ψij=arctan⁡(zi,j+1−zij∥Pij′Pi,j+1′→∥)\\psi_{ij} = \\arctan\\left( \\frac{z_{i,j+1} - z_{ij}}{\\| \\overrightarrow{P\'_{ij}P\'_{i,j+1}} \\|} \\right)ψij​=arctan​∥Pij′​Pi,j+1′​​∥zi,j+1​−zij​​​

然后，平滑成本计算为：

F4(X)=a1∑j=1n−2ϕij+a2∑j=1n−1∣ψij−ψj−1∣F_4(X) = a_1 \\sum_{j=1}^{n-2} \\phi_{ij} + a_2 \\sum_{j=1}^{n-1} |\\psi_{ij} - \\psi_{j-1}|F4​(X)=a1​j=1∑n−2​ϕij​+a2​j=1∑n−1​∣ψij​−ψj−1​∣

其中 a1a_1a1​ 和 a2a_2a2​ 分别是转弯角和爬升角的惩罚系数。

2.3 总体成本函数

2.3.1 单个无人成本计算

考虑到路径 $X$ 的最优性、安全性和可行性约束，第$i$ 个无人机总体成本函数可以定义为以下形式：

fi(X)=∑k=14bkFk(Xi)f_i(X) = \\sum_{k=1}^4 b_k F_k(X_i)fi​(X)=k=1∑4​bk​Fk​(Xi​)

其中 bkb_kbk​ 是权重系数，F1(Xi)F_1(X_i)F1​(Xi​) 到 F4(Xi)F_4(X_i)F4​(Xi​) 分别是路径长度、威胁、平滑度和飞行高度相关的成本。决策变量是 $X$，包括 $n$ 个航路点 Pij=(xij,yij,zij)P_{ij} = (x_{ij}, y_{ij}, z_{ij})Pij​=(xij​,yij​,zij​) 的列表，使得 Pij∈OP_{ij} \\in OPij​∈O，其中 $O$ 是无人机的操作空间。根据这些定义，成本函数 $F$ 是完全确定的，可以作为路径规划过程的输入。

2.3.2 多无人机总成本计算

若共有$m$ 个无人机，其总成本为单个无人机成本和，计算公式如下：
fitness(X)=∑i=1mfi(X)fitness(X) = \\sum_{i=1}^mf_i(X)fitness(X)=i=1∑m​fi​(X)
参考文献：
[1] Phung M D , Ha Q P .Safety-enhanced UAV path planning with spherical vector-based particle swarm optimization[J].Applied Soft Computing, 2021(2):107376.DOI:10.1016/j.asoc.2021.107376.

三、部分代码及结果

close allclearclc% dbstop if all errorpop=100;%种群大小（可以修改）maxgen=200;%最大迭代（可以修改）%% 模型建立model=Create_Model();UAVnum=4;%无人机数量(可以修改) 必须与无人机的起始点保持一致%% 初始化每个无人机的模型for i=1:UAVnum ModelUAV(i).model=model;end%% 第一个无人机 起始点start_location = [120;200;100];end_location = [800;800;150];ModelUAV(1).model.start=start_location;ModelUAV(1).model.end=end_location;%% 第二个无人机 起始点start_location = [400;100;100];end_location = [900;600;150];ModelUAV(2).model.start=start_location;ModelUAV(2).model.end=end_location;%% 第三个无人机 起始点start_location = [200;150;150];end_location =[850;750;150];ModelUAV(3).model.start=start_location;ModelUAV(3).model.end=end_location;%% 第四个无人机 起始点start_location = [100;100;150];end_location = [800;730;150];ModelUAV(4).model.start=start_location;ModelUAV(4).model.end=end_location;%% 第5个无人机 起始点% start_location = [500;100;130];% end_location = [850;650;150];% ModelUAV(5).model.start=start_location;% ModelUAV(5).model.end=end_location;% %% 第6个无人机 起始点% start_location = [100;100;150];% end_location = [800;800;150];% ModelUAV(6).model.start=start_location;% ModelUAV(6).model.end=end_location;

五个无人机：

四、完整MATLAB代码见下方名片