基于多目标向日葵优化算法（Multi-objective Sunflower Optimization，MOSFO）的移动机器人路径规划研究，MATLAB代码

一、机器人路径规划介绍

移动机器人路径规划是机器人研究的重要分支，是对其进行控制的基础。根据环境信息的已知程度不同，路径规划分为基于环境信息已知的全局路径规划和基于环境信息未知或局部已知的局部路径规划。随着科技的快速发展以及机器人的大量应用，人们对机器人的要求也越来越高，尤其表现在对机器人的智能化方面的要求，而机器人自主路径规划是实现机器人智能化的重要步骤，路径规划是指规划机器人从起点位置出发，无碰撞、安全到达指定目标位置的最优路径。目前，常用的移动机器人全局路径规划方法很多，如栅格法和人工势场法。对于栅格法，当空间增大时，所需存储空间剧增，决策速度下降；而人工势场法容易产生局部最优解问题和死锁现象。随着智能控制技术的发展，出现了如遗传算法、粒子群优化算法、麻雀搜索算法、灰狼优化算法、鲸鱼优化算法等。

参考文献：

[1]史恩秀,陈敏敏,李俊,等.基于蚁群算法的移动机器人全局路径规划方法研究[J].农业机械学报, 2014, 45(6):5.DOI:CNKI:SUN:NYJX.0.2014-06-009.

[2]朱庆保,张玉兰.基于栅格法的机器人路径规划蚁群算法[J].机器人, 2005, 27(2):5.DOI:10.3321/j.issn:1002-0446.2005.02.008.

[3]曹新亮,王智文,冯晶,等.基于改进蚁群算法的机器人全局路径规划研究[J].计算机工程与科学, 2020, 42(3):7.DOI:CNKI:SUN:JSJK.0.2020-03-027.

二、栅格地图环境搭建

首先建立移动机器人工作环境，设移动机器人的工作空间为二维空间（记为RS），工作环境中的障碍物即为机床。在机器人运动过程中，障碍物为静止且大小不发生变化。按栅格法划分RS，移动机器人在栅格间行走。无障碍物的栅格为可行栅格，有障碍物的栅格为不可行栅格。栅格集包含所有栅格。栅格标识有：直角坐标法和序号法。本文采用序号标识法。

在移动机器人工作空间下按从左到右，从上到下的顺序，依次标记为序号1，2，3，···，n，每一个序号代表一个栅格。为了避免移动机器人与障碍物发生碰撞，可以将障碍物膨胀，障碍物在占原有栅格的同时，再占多个栅格，按个栅格算。这种划分方法简单实用，能够满足环境模型与真实情况相符。从而使移动机器人在路径规划时畅通无阻。令S＝｛1，2，3，···，N｝为栅格序号集。根据上述对应关系，可知g（0，0）的序号为1，g（1，0）序号为2，直至g（X，Y）的序号为n。规划起始位置、目标位置均为任意且都属于S（但不在同一栅格内）。

在实际工作环境中，移动机器人工作环境是复杂多变的，且为三维空间。为了便于研究，本文对环境进行简化建模。栅格法是一种常用的环境表示方法，因其简单方便（二维环境），环境建模的复杂性小，因而本文环境建模采用栅格法。在栅格地图中，工作环境被划分为很多栅格，其中包括有障碍物和无障碍的栅格，在仿真程序中用0表示此栅格无障碍物，机器人可以通过此栅格，用1表示栅格有障碍物，机器人无法通过，需选择其他栅格。栅格的尺寸大小可根据工作环境中的障碍物尺寸以及安全距离进行设置。为了实现程序仿真，需要对栅格进行标识，如下图所示，以20x20的栅格环境为例来说明。

在这里插入图片描述

如上图所示，白色栅格表示无障碍物的栅格，黑色栅格则表示有障碍物的栅格，在地图中对每个栅格编号，不同序号的栅格在坐标系中的坐标可用下式来表示：

$\text{mod}(N_i / N) - 0.5$

$\text{ceil}(N_i / N) + 0.5$

其中， $\text{mod}$ 为取余运算， $\text{ceil}$ 表示向后取整， $N_i$ 是对应栅格的标号， $N$ 表示每列的栅格数量，取栅格中心位置作为栅格在坐标系中的坐标。这样机器人全局路径规划的问题就转变成了利用算法在栅格地图上寻找由起始点到目标点的有序的栅格子集，这些栅格子集的中心连线便是算法寻找的路径。

参考文献：

[1]史恩秀,陈敏敏,李俊,等.基于蚁群算法的移动机器人全局路径规划方法研究[J].农业机械学报, 2014, 45(6):5.DOI:CNKI:SUN:NYJX.0.2014-06-009.

[2]曹新亮,王智文,冯晶,等.基于改进蚁群算法的机器人全局路径规划研究[J].计算机工程与科学, 2020, 42(3):7.DOI:CNKI:SUN:JSJK.0.2020-03-027.

三、机器人路径规划多目标模型

3.1 路径成本

当机器人从起点向目标点移动时，通常选择最短的一条路径。将一条路径上每条线段的长度累加求和即可得到路径总长度。任何两点形成的线段都是根据欧氏距离来计算的，其中， $[P_0, P_1, \ldots, P_i, P_{i+1}, \ldots, P_n, P_{n+1}]$ 代表路径 $P$ ， $S = P_0$ 代表起点， $T = P_{n+1}$ 代表目标点。路径长度目标计算方法如下式：

$\text{Length}(P) = \sum_{i=0}^{n} d(P_i, P_{i+1})$

其中， $P_i = (x_i, y_i)$ 和 $P_{i+1} = (x_{i+1}, y_{i+1})$ 是路径中的两个连续点， $d(P_i, P_{i+1})$ 是路径中的线段距离； $\text{Length}(P)$ 表示由所有线段相加得到的总路径长度； $n$ 表示路径中的点数量。

3.2 平滑成本

路径平滑度表示路径的弯曲程度，换句话说，只有当路径是平滑的，机器人在移动时才会使用较少的能量。为了衡量可行路径的平滑度，使用两个连续路段之间的夹角 $\text{Ang}[P_i, P_{i+1}, P_{i+2}]$ 表示。路径平滑度的计算方法如下式：

$\text{Smoothness}(P) = \sum_{i=0}^{n-2} \left| \theta_i \right|$

其中， $\theta_i$ 是路径上三个相邻点 $P_i, P_{i+1}, P_{i+2}$ 之间的夹角。

3.3 目标函数

移动机器人路径规划的目标函数 $f_1$ 和 $f_2$ 分别是路径成本最小和平滑成本最小：

$f_1 = \min(\text{Length}(P))$

$f_2 = \min(\text{Smoothness}(P))$

参考文献：

[1]于振翱. 面向多目标优化的移动机器人路径规划方法研究[D]. 山东:聊城大学,2023.

[2]杨嘉. 基于改进NSGA-Ⅱ算法的移动机器人路径规划研究[D]. 陕西:长安大学,2021.

四、多目标向日葵优化算法求解移动机器人路径规划

4.1多目标向日葵优化算法

多目标向日葵优化算法（Multi-objective Sunflower Optimization，MOSFO）是一种用于解决多目标优化问题的强大学术研究算法，由 Joao Luiz Junho Pereira 等人于 2023 年提出。该算法是具有超立方和约束的多目标元启发式算法，其灵感源于向日葵围绕太阳的向光性生命周期。
算法原理
MOSFO 通过模拟向日葵的生长和追光机制来寻找最优解集合。在自然界中，向日葵会随着太阳位置的变化而转动自己的花盘以获得最佳的光照条件。同样地，在多目标优化问题中，MOSFO 旨在通过模拟这一过程来寻找问题的最优解集合。尤其适用于需要同时考虑多个目标函数的优化问题，能够为决策者提供一组具有不同目标权衡的解。
算法流程
初始化种群：在搜索空间中随机生成初始种群，种群中的每个个体代表一个可能的解。
适应度评估：计算每个个体针对不同目标函数的适应度值。
更新种群：根据向日葵的向光性行为更新种群的位置，通过选择邻居并朝向“太阳”（即最优解的方向）移动。
保持多样性：采用某种机制保持种群的多样性，例如外部存档或 crowding distance 方法，以探索新的区域。
终止条件：达到最大迭代次数或其他终止条件后，输出帕累托前沿。
优势
强大的搜索能力：MOSFO 在解决含有多个目标的复杂工程问题上表现出色，能在帕累托前沿上找到分布更广泛、质量更高的解。
效率高：相比其他算法，MOSFO 在保持解的质量的同时，能以更快的速度收敛，具有更低的计算成本。
实现简单：编程模型比大多数进化算法更简单，更容易在不同领域应用。
适应性强：对不同目标数量的问题都适用，尤其在目标数量较多时表现更优
参考文献：
[1] Pereira J L J , Gomes G F .Multi‐objective sunflower optimization: A new hypercubic meta‐heuristic for constrained engineering problems[J].Expert Systems, 2023, 40(8).DOI:10.1111/exsy.13331.

4.2 部分代码

close all
clear
clc
dbstop if all error
addpath("./ /")
global G S E
MultiObj= fun_info();%获取无人机模型信息
params.maxgen=100;  % 最大迭代次数
params.Np=50;      % 种群大小
params.Nr=100;      %外部存档大小（不得小于种群大小）
[Xbest,Fbest] = (params,MultiObj);
%% 获取算法得到的所有路径存于Result中
for i=1:size(Xbest,1)
global_best = round(Xbest(i,:));route = [S(1) global_best E(1)];path=generateContinuousRoute(route,G);% path=shortenRoute(path);path=GenerateSmoothPath(path,G);Result(i).BestPosition= path;Result(i).BestFit=Fbest(i,:);
end
%% 获取 路径成本最小 与 平滑成本最小
if size(Fbest,1)>1idx=min(Fbest);
elseidx=Fbest;
end
KK=find(Fbest(:,1)==idx(1));
IDX(1)=KK(1);
KK=find(Fbest(:,2)==idx(2));
IDX(2)=KK(1);
%%  画图pareto前沿图
dbclear all
figure
plot(Fbest(:,1),Fbest(:,2),'r*');
xlabel('路径成本')
ylabel('平滑成本')
legend('')