MATLAB 实现基于多目标粒子群算法（MOPSO）的无人机三维路径规划

% 初始化粒子群参数
numParticles = 50; % 粒子数量
dim = 3; % 粒子维度，对应三维空间坐标

% 初始化粒子位置和速度
positions = rand(numParticles, dim) * 100; % 粒子初始位置随机分布
velocities = zeros(numParticles, dim); % 粒子速度初始为 0

% 初始化个体最优和全局最优存储结构
pbest = positions;
pbest_scores = inf(numParticles, 2); % 存储多目标得分，如路径长度和避障风险

% 外部存档用于存储非支配解
archive.positions = [];

% 目标函数定义（示例：路径长度和避障风险）
function scores = objectiveFunction(pos)
    pathLength = sum(sqrt(sum(diff(pos).^2, 2)));
    obstacleRisk = computeObstacleRisk(pos);
    scores = [pathLength, obstacleRisk];
end

% 粒子速度更新参数
w_max = 0.9;
w_min = 0.4;
c1 = 2; c2 = 2;

% 迭代过程主循环
for iter = 1:maxIter
    w = w_max - (w_max - w_min) * (iter / maxIter);
    
    for i = 1:numParticles
        currentPos = positions(i, :);
        currentScore = objectiveFunction(currentPos);
        
        % 更新个体最优
        if dominates(currentScore, pbest_scores(i, :))
            pbest(i, :) = currentPos;
            pbest_scores(i, :) = currentScore;
        elseif ~dominates(pbest_scores(i, :), currentScore)
            % 非支配关系时可采用拥挤度等策略进行选择
        end
        
        % 选择全局最优粒子
        gbest = selectGlobalBest(archive);
        
        % 更新速度
        r1 = rand(1, dim);
        velocities(i, :) = w * velocities(i, :) ... + c1 * r1 .* (pbest(i, :) - currentPos) ... + c2 * rand(1, dim) .* (gbest - currentPos);
        
        % 限制速度范围（避免爆炸）
        maxV = 10;
        velocities(i, :) = min(max(velocities(i, :), -maxV), maxV);
        
        % 更新位置
        positions(i, :) = positions(i, :) + velocities(i, :);
        
        % 限制位置边界（保证在搜索空间内）
        positions(i, :) = min(max(positions(i, :), 0), 100);
    end
    
    % 更新非支配解集（外部存档）
    archive = updateArchive(archive, positions, pbest_scores);
end

% 支配关系判断函数
function flag = dominates(score1, score2)
    % 若 score1 需全部目标不劣且至少一个目标更优
    flag = all(score1 <= score2) && any(score1 < score2);
end

% 选择全局最优函数（基于拥挤度距离或其他策略）
function gbest = selectGlobalBest(archive)
    % 简单示例：随机选择一个存档中的非支配解作为全局最优
    if isempty(archive.positions)
        gbest = []; return;
    end
    idx = randi(size(archive.positions, 1));
    gbest = archive.positions(idx, :);
end

% 更新存档函数（维护非支配解集）
function archive = updateArchive(archive, positions, scores)
    combinedPositions = [archive.positions; positions];
    combinedScores = [archive.scores; scores];
    ndIdx = nondominatedSort(combinedScores);
    archive.positions = combinedPositions(ndIdx, :);
    archive.scores = combinedScores(ndIdx, :);
end

% 非支配排序函数示例
function ndIdx = nondominatedSort(scores)
    num = size(scores, 1);
    dominated = false(num, 1);
    for i = 1:num
        for j = 1:num
            if i ~= j && dominates(scores(j, :), scores(i, :))
                dominated(i) = true;
                break;
            end
        end
    end
    ndIdx = find(~dominated);
end

% 计算粒子与障碍物距离的示例函数
function distances = computeObstacleRisk(pos)
    global obstacles;
    distances = zeros(size(obstacles, 1), 1);
    for k = 1:size(obstacles, 1)
        distances(k) = norm(pos - obstacles(k, :));
    end
    distances = min(distances);
end

以上代码示例展示了多目标粒子群算法的核心组成部分，包括初始化粒子群、目标函数定义、多目标支配关系判断、速度与位置更新、外部存档维护和非支配排序机制。