粒子群優(yōu)化算法(MOPSO)原理與PyTorch實(shí)現(xiàn))
1. 項(xiàng)目概述當(dāng)粒子群遇上多目標(biāo)優(yōu)化在工程優(yōu)化和機(jī)器學(xué)習(xí)領(lǐng)域我們常常需要同時(shí)優(yōu)化多個(gè)相互沖突的目標(biāo)函數(shù)。比如設(shè)計(jì)一款電動(dòng)汽車時(shí)既要最大化續(xù)航里程又要最小化電池成本這兩個(gè)目標(biāo)往往難以同時(shí)滿足。傳統(tǒng)粒子群算法(PSO)擅長處理單目標(biāo)優(yōu)化但面對這類多目標(biāo)問題就顯得力不從心。這正是多目標(biāo)粒子群算法(MOPSO)大顯身手的地方。MOPSO通過維護(hù)一個(gè)外部存檔來保存找到的非支配解即Pareto最優(yōu)解并采用特殊的選擇機(jī)制來引導(dǎo)粒子群向真實(shí)的Pareto前沿收斂。與NSGA-II等遺傳算法相比MOPSO具有收斂速度快、參數(shù)調(diào)節(jié)簡單等優(yōu)勢。而使用PyTorch實(shí)現(xiàn)MOPSO則能充分發(fā)揮GPU并行計(jì)算的優(yōu)勢大幅提升算法運(yùn)行效率。2. 核心算法原理拆解2.1 多目標(biāo)優(yōu)化問題數(shù)學(xué)表述一個(gè)典型的多目標(biāo)優(yōu)化問題可以表示為最小化 F(x) [f?(x), f?(x), ..., f?(x)] 滿足 g?(x) ≤ 0, i1,2,...,m h?(x) 0, j1,2,...,p其中x∈??是決策變量F:??→??是目標(biāo)函數(shù)向量。Pareto最優(yōu)解的定義是在可行解空間中如果不存在其他解在所有目標(biāo)上都不差于它且至少在一個(gè)目標(biāo)上嚴(yán)格優(yōu)于它那么這個(gè)解就是Pareto最優(yōu)解。2.2 標(biāo)準(zhǔn)粒子群算法的局限性標(biāo)準(zhǔn)PSO的粒子更新公式為v? w·v? c?·r?·(pbest? - x?) c?·r?·(gbest - x?) x? x? v?其中w是慣性權(quán)重c?、c?是學(xué)習(xí)因子r?、r?是隨機(jī)數(shù)。問題在于多目標(biāo)情況下無法直接確定全局最優(yōu)解gbest因?yàn)榭赡艽嬖诙鄠€(gè)非支配解。2.3 MOPSO的核心改進(jìn)MOPSO引入了三個(gè)關(guān)鍵機(jī)制外部存檔存儲(chǔ)找到的非支配解采用自適應(yīng)網(wǎng)格法維護(hù)存檔多樣性領(lǐng)導(dǎo)者選擇從存檔中隨機(jī)選擇引導(dǎo)粒子(gbest)變異算子防止早熟收斂增強(qiáng)探索能力3. PyTorch實(shí)現(xiàn)詳解3.1 環(huán)境配置與依賴安裝推薦使用conda創(chuàng)建虛擬環(huán)境conda create -n mopso python3.9 conda activate mopso conda install pytorch torchvision torchaudio cudatoolkit11.3 -c pytorch pip install matplotlib numpy pandas3.2 算法核心類設(shè)計(jì)import torch import numpy as np from typing import List, Tuple class MOPSO: def __init__(self, obj_funcs: List[callable], bounds: torch.Tensor, n_particles: int 100, n_iter: int 200, inertia: float 0.7, personal_weight: float 1.5, global_weight: float 1.5): obj_funcs: 目標(biāo)函數(shù)列表 bounds: 決策變量邊界 [n_dim, 2] self.obj_funcs obj_funcs self.bounds bounds.to(device) self.n_particles n_particles self.n_iter n_iter self.w inertia self.c1 personal_weight self.c2 global_weight self.device torch.device(cuda if torch.cuda.is_available() else cpu) # 初始化粒子位置和速度 self.n_dim bounds.shape[0] self.particles torch.rand((n_particles, self.n_dim), deviceself.device) * (bounds[:,1]-bounds[:,0]) bounds[:,0] self.velocities torch.zeros_like(self.particles) # 初始化個(gè)體最優(yōu) self.pbest self.particles.clone() self.pbest_obj self.evaluate(self.particles) # 初始化外部存檔 self.archive [] self.update_archive(self.particles, self.pbest_obj) def evaluate(self, x: torch.Tensor) - torch.Tensor: 評估粒子在多目標(biāo)上的表現(xiàn) return torch.stack([f(x) for f in self.obj_funcs], dim1) def update_archive(self, particles: torch.Tensor, objectives: torch.Tensor): 更新外部存檔 # 非支配排序?qū)崿F(xiàn)... pass def select_leader(self) - torch.Tensor: 從存檔中選擇引導(dǎo)粒子 # 基于擁擠距離的選擇... pass def run(self): 主優(yōu)化循環(huán) for _ in range(self.n_iter): # 選擇領(lǐng)導(dǎo)者 leaders self.select_leader() # 更新速度和位置 r1 torch.rand_like(self.particles) r2 torch.rand_like(self.particles) self.velocities (self.w * self.velocities self.c1 * r1 * (self.pbest - self.particles) self.c2 * r2 * (leaders - self.particles)) self.particles torch.clamp(self.particles self.velocities, self.bounds[:,0], self.bounds[:,1]) # 評估新位置 current_obj self.evaluate(self.particles) # 更新個(gè)體最優(yōu) mask self.is_dominated(current_obj, self.pbest_obj) self.pbest[mask] self.particles[mask] self.pbest_obj[mask] current_obj[mask] # 更新存檔 self.update_archive(self.particles, current_obj) # 應(yīng)用變異算子 self.apply_mutation()3.3 關(guān)鍵技術(shù)實(shí)現(xiàn)細(xì)節(jié)3.3.1 非支配排序?qū)崿F(xiàn)def is_dominated(self, a: torch.Tensor, b: torch.Tensor) - torch.Tensor: 判斷a是否支配b (a不差于b且至少一個(gè)目標(biāo)更好) not_worse (a b).all(dim1) better (a b).any(dim1) return not_worse better def fast_non_dominated_sort(self, objectives: torch.Tensor) - List[torch.Tensor]: 快速非支配排序 fronts [] remaining torch.arange(objectives.shape[0]) while len(remaining) 0: front [] for i in remaining: dominated False for j in remaining: if self.is_dominated(objectives[j], objectives[i]): dominated True break if not dominated: front.append(i) front torch.tensor(front, deviceself.device) fronts.append(front) remaining remaining[~torch.isin(remaining, front)] return fronts3.3.2 自適應(yīng)網(wǎng)格維護(hù)存檔def update_archive(self, particles: torch.Tensor, objectives: torch.Tensor): 使用自適應(yīng)網(wǎng)格法維護(hù)存檔 # 合并新解和現(xiàn)有存檔 all_solutions torch.cat([objectives, torch.stack([a[1] for a in self.archive]) if self.archive else torch.empty((0, len(self.obj_funcs)), deviceself.device)]) # 非支配排序 fronts self.fast_non_dominated_sort(all_solutions) new_archive [] for front in fronts: if len(new_archive) len(front) self.max_archive_size: # 使用擁擠距離選擇最具代表性的解 selected self.select_by_crowding(all_solutions[front]) new_archive.extend(selected) break else: new_archive.extend(front.tolist()) # 更新存檔 self.archive [(particles[i], all_solutions[i]) for i in new_archive]4. 實(shí)戰(zhàn)案例電機(jī)設(shè)計(jì)優(yōu)化4.1 問題描述我們以永磁同步電機(jī)設(shè)計(jì)為例優(yōu)化三個(gè)目標(biāo)最大化效率 η最小化成本 Cost最小化重量 Weight決策變量包括定子外徑 Dso ∈ [100, 200] mm定子內(nèi)徑 Dsi ∈ [50, 150] mm氣隙長度 g ∈ [0.5, 2] mm永磁體厚度 hm ∈ [3, 10] mm4.2 目標(biāo)函數(shù)實(shí)現(xiàn)def efficiency(x: torch.Tensor) - torch.Tensor: 計(jì)算電機(jī)效率 # 簡化的效率計(jì)算公式 Dso, Dsi, g, hm x.T return 0.9 - 0.001*g - 0.0005*hm 0.00001*(Dso-Dsi) def cost(x: torch.Tensor) - torch.Tensor: 計(jì)算電機(jī)成本 Dso, Dsi, g, hm x.T material_cost 0.1*(Dso**2 - Dsi**2) 5*hm return material_cost 50 # 固定加工成本 def weight(x: torch.Tensor) - torch.Tensor: 計(jì)算電機(jī)重量 Dso, Dsi, g, hm x.T return 0.01*(Dso**2 - Dsi**2) 0.05*hm4.3 優(yōu)化過程可視化import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D def plot_pareto_front(archive): 繪制三維Pareto前沿 objs torch.stack([a[1] for a in archive]).cpu().numpy() fig plt.figure(figsize(10,8)) ax fig.add_subplot(111, projection3d) ax.scatter(objs[:,0], objs[:,1], objs[:,2], cr, markero) ax.set_xlabel(Efficiency) ax.set_ylabel(Cost) ax.set_zlabel(Weight) plt.title(Pareto Front for Motor Design) plt.show()5. 性能優(yōu)化技巧5.1 GPU加速策略批量評估將所有粒子的目標(biāo)函數(shù)評估合并為一次矩陣運(yùn)算def evaluate_batch(self, x: torch.Tensor) - torch.Tensor: 批量評估所有粒子 # x shape: [n_particles, n_dim] efficiency 0.9 - 0.001*x[:,2] - 0.0005*x[:,3] 0.00001*(x[:,0]-x[:,1]) cost 0.1*(x[:,0]**2 - x[:,1]**2) 5*x[:,3] 50 weight 0.01*(x[:,0]**2 - x[:,1]**2) 0.05*x[:,3] return torch.stack([efficiency, cost, weight], dim1)內(nèi)存優(yōu)化使用原地操作減少內(nèi)存分配self.velocities.mul_(self.w).add_( self.c1 * r1 * (self.pbest - self.particles), alpha1 ).add_( self.c2 * r2 * (leaders - self.particles), alpha1 )5.2 參數(shù)調(diào)優(yōu)指南參數(shù)推薦范圍影響調(diào)整策略粒子數(shù)50-500探索能力問題維度越高需要越多粒子慣性權(quán)重w0.4-0.9平衡探索與開發(fā)線性遞減效果更好學(xué)習(xí)因子c1,c21.5-2.5收斂速度c1略大于c2增強(qiáng)多樣性存檔大小50-200解的質(zhì)量越大Pareto前沿越完整變異概率0.05-0.2避免早熟初期可設(shè)較高值6. 常見問題與解決方案6.1 收斂問題排查問題1算法過早收斂到局部Pareto前沿檢查存檔中解的分布是否過于集中解決增加變異概率、減小慣性權(quán)重、使用動(dòng)態(tài)參數(shù)調(diào)整問題2Pareto前沿不完整檢查存檔大小是否足夠解決增大存檔容量、增加粒子數(shù)量、延長迭代次數(shù)6.2 數(shù)值穩(wěn)定性問題問題目標(biāo)函數(shù)量綱差異導(dǎo)致某些目標(biāo)被忽略# 標(biāo)準(zhǔn)化處理 def normalize_objectives(self, objs: torch.Tensor) - torch.Tensor: 將各目標(biāo)函數(shù)值歸一化到[0,1]范圍 min_vals objs.min(dim0)[0] max_vals objs.max(dim0)[0] return (objs - min_vals) / (max_vals - min_vals 1e-8)6.3 并行計(jì)算優(yōu)化使用PyTorch的分布式包實(shí)現(xiàn)多GPU訓(xùn)練import torch.distributed as dist def init_process(rank, world_size): dist.init_process_group(gloo, rankrank, world_sizeworld_size) # 分割粒子到不同GPU particles_per_gpu self.n_particles // world_size start rank * particles_per_gpu end (rank 1) * particles_per_gpu if rank ! world_size-1 else self.n_particles # 各GPU處理自己的粒子 local_particles self.particles[start:end] local_velocities self.velocities[start:end] # 處理完成后同步結(jié)果 dist.all_gather(all_particles, local_particles) 多目標(biāo)粒子群優(yōu)化算法(MOPSO)原理與PyTorch實(shí)現(xiàn)
1. 項(xiàng)目概述當(dāng)粒子群遇上多目標(biāo)優(yōu)化在工程優(yōu)化和機(jī)器學(xué)習(xí)領(lǐng)域我們常常需要同時(shí)優(yōu)化多個(gè)相互沖突的目標(biāo)函數(shù)。比如設(shè)計(jì)一款電動(dòng)汽車時(shí)既要最大化續(xù)航里程又要最小化電池成本這兩個(gè)目標(biāo)往往難以同時(shí)滿足。傳統(tǒng)粒子群算法(PSO)擅長處理單目標(biāo)優(yōu)化但面對這類多目標(biāo)問題就顯得力不從心。這正是多目標(biāo)粒子群算法(MOPSO)大顯身手的地方。MOPSO通過維護(hù)一個(gè)外部存檔來保存找到的非支配解即Pareto最優(yōu)解并采用特殊的選擇機(jī)制來引導(dǎo)粒子群向真實(shí)的Pareto前沿收斂。與NSGA-II等遺傳算法相比MOPSO具有收斂速度快、參數(shù)調(diào)節(jié)簡單等優(yōu)勢。而使用PyTorch實(shí)現(xiàn)MOPSO則能充分發(fā)揮GPU并行計(jì)算的優(yōu)勢大幅提升算法運(yùn)行效率。2. 核心算法原理拆解2.1 多目標(biāo)優(yōu)化問題數(shù)學(xué)表述一個(gè)典型的多目標(biāo)優(yōu)化問題可以表示為最小化 F(x) [f?(x), f?(x), ..., f?(x)] 滿足 g?(x) ≤ 0, i1,2,...,m h?(x) 0, j1,2,...,p其中x∈??是決策變量F:??→??是目標(biāo)函數(shù)向量。Pareto最優(yōu)解的定義是在可行解空間中如果不存在其他解在所有目標(biāo)上都不差于它且至少在一個(gè)目標(biāo)上嚴(yán)格優(yōu)于它那么這個(gè)解就是Pareto最優(yōu)解。2.2 標(biāo)準(zhǔn)粒子群算法的局限性標(biāo)準(zhǔn)PSO的粒子更新公式為v? w·v? c?·r?·(pbest? - x?) c?·r?·(gbest - x?) x? x? v?其中w是慣性權(quán)重c?、c?是學(xué)習(xí)因子r?、r?是隨機(jī)數(shù)。問題在于多目標(biāo)情況下無法直接確定全局最優(yōu)解gbest因?yàn)榭赡艽嬖诙鄠€(gè)非支配解。2.3 MOPSO的核心改進(jìn)MOPSO引入了三個(gè)關(guān)鍵機(jī)制外部存檔存儲(chǔ)找到的非支配解采用自適應(yīng)網(wǎng)格法維護(hù)存檔多樣性領(lǐng)導(dǎo)者選擇從存檔中隨機(jī)選擇引導(dǎo)粒子(gbest)變異算子防止早熟收斂增強(qiáng)探索能力3. PyTorch實(shí)現(xiàn)詳解3.1 環(huán)境配置與依賴安裝推薦使用conda創(chuàng)建虛擬環(huán)境conda create -n mopso python3.9 conda activate mopso conda install pytorch torchvision torchaudio cudatoolkit11.3 -c pytorch pip install matplotlib numpy pandas3.2 算法核心類設(shè)計(jì)import torch import numpy as np from typing import List, Tuple class MOPSO: def __init__(self, obj_funcs: List[callable], bounds: torch.Tensor, n_particles: int 100, n_iter: int 200, inertia: float 0.7, personal_weight: float 1.5, global_weight: float 1.5): obj_funcs: 目標(biāo)函數(shù)列表 bounds: 決策變量邊界 [n_dim, 2] self.obj_funcs obj_funcs self.bounds bounds.to(device) self.n_particles n_particles self.n_iter n_iter self.w inertia self.c1 personal_weight self.c2 global_weight self.device torch.device(cuda if torch.cuda.is_available() else cpu) # 初始化粒子位置和速度 self.n_dim bounds.shape[0] self.particles torch.rand((n_particles, self.n_dim), deviceself.device) * (bounds[:,1]-bounds[:,0]) bounds[:,0] self.velocities torch.zeros_like(self.particles) # 初始化個(gè)體最優(yōu) self.pbest self.particles.clone() self.pbest_obj self.evaluate(self.particles) # 初始化外部存檔 self.archive [] self.update_archive(self.particles, self.pbest_obj) def evaluate(self, x: torch.Tensor) - torch.Tensor: 評估粒子在多目標(biāo)上的表現(xiàn) return torch.stack([f(x) for f in self.obj_funcs], dim1) def update_archive(self, particles: torch.Tensor, objectives: torch.Tensor): 更新外部存檔 # 非支配排序?qū)崿F(xiàn)... pass def select_leader(self) - torch.Tensor: 從存檔中選擇引導(dǎo)粒子 # 基于擁擠距離的選擇... pass def run(self): 主優(yōu)化循環(huán) for _ in range(self.n_iter): # 選擇領(lǐng)導(dǎo)者 leaders self.select_leader() # 更新速度和位置 r1 torch.rand_like(self.particles) r2 torch.rand_like(self.particles) self.velocities (self.w * self.velocities self.c1 * r1 * (self.pbest - self.particles) self.c2 * r2 * (leaders - self.particles)) self.particles torch.clamp(self.particles self.velocities, self.bounds[:,0], self.bounds[:,1]) # 評估新位置 current_obj self.evaluate(self.particles) # 更新個(gè)體最優(yōu) mask self.is_dominated(current_obj, self.pbest_obj) self.pbest[mask] self.particles[mask] self.pbest_obj[mask] current_obj[mask] # 更新存檔 self.update_archive(self.particles, current_obj) # 應(yīng)用變異算子 self.apply_mutation()3.3 關(guān)鍵技術(shù)實(shí)現(xiàn)細(xì)節(jié)3.3.1 非支配排序?qū)崿F(xiàn)def is_dominated(self, a: torch.Tensor, b: torch.Tensor) - torch.Tensor: 判斷a是否支配b (a不差于b且至少一個(gè)目標(biāo)更好) not_worse (a b).all(dim1) better (a b).any(dim1) return not_worse better def fast_non_dominated_sort(self, objectives: torch.Tensor) - List[torch.Tensor]: 快速非支配排序 fronts [] remaining torch.arange(objectives.shape[0]) while len(remaining) 0: front [] for i in remaining: dominated False for j in remaining: if self.is_dominated(objectives[j], objectives[i]): dominated True break if not dominated: front.append(i) front torch.tensor(front, deviceself.device) fronts.append(front) remaining remaining[~torch.isin(remaining, front)] return fronts3.3.2 自適應(yīng)網(wǎng)格維護(hù)存檔def update_archive(self, particles: torch.Tensor, objectives: torch.Tensor): 使用自適應(yīng)網(wǎng)格法維護(hù)存檔 # 合并新解和現(xiàn)有存檔 all_solutions torch.cat([objectives, torch.stack([a[1] for a in self.archive]) if self.archive else torch.empty((0, len(self.obj_funcs)), deviceself.device)]) # 非支配排序 fronts self.fast_non_dominated_sort(all_solutions) new_archive [] for front in fronts: if len(new_archive) len(front) self.max_archive_size: # 使用擁擠距離選擇最具代表性的解 selected self.select_by_crowding(all_solutions[front]) new_archive.extend(selected) break else: new_archive.extend(front.tolist()) # 更新存檔 self.archive [(particles[i], all_solutions[i]) for i in new_archive]4. 實(shí)戰(zhàn)案例電機(jī)設(shè)計(jì)優(yōu)化4.1 問題描述我們以永磁同步電機(jī)設(shè)計(jì)為例優(yōu)化三個(gè)目標(biāo)最大化效率 η最小化成本 Cost最小化重量 Weight決策變量包括定子外徑 Dso ∈ [100, 200] mm定子內(nèi)徑 Dsi ∈ [50, 150] mm氣隙長度 g ∈ [0.5, 2] mm永磁體厚度 hm ∈ [3, 10] mm4.2 目標(biāo)函數(shù)實(shí)現(xiàn)def efficiency(x: torch.Tensor) - torch.Tensor: 計(jì)算電機(jī)效率 # 簡化的效率計(jì)算公式 Dso, Dsi, g, hm x.T return 0.9 - 0.001*g - 0.0005*hm 0.00001*(Dso-Dsi) def cost(x: torch.Tensor) - torch.Tensor: 計(jì)算電機(jī)成本 Dso, Dsi, g, hm x.T material_cost 0.1*(Dso**2 - Dsi**2) 5*hm return material_cost 50 # 固定加工成本 def weight(x: torch.Tensor) - torch.Tensor: 計(jì)算電機(jī)重量 Dso, Dsi, g, hm x.T return 0.01*(Dso**2 - Dsi**2) 0.05*hm4.3 優(yōu)化過程可視化import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D def plot_pareto_front(archive): 繪制三維Pareto前沿 objs torch.stack([a[1] for a in archive]).cpu().numpy() fig plt.figure(figsize(10,8)) ax fig.add_subplot(111, projection3d) ax.scatter(objs[:,0], objs[:,1], objs[:,2], cr, markero) ax.set_xlabel(Efficiency) ax.set_ylabel(Cost) ax.set_zlabel(Weight) plt.title(Pareto Front for Motor Design) plt.show()5. 性能優(yōu)化技巧5.1 GPU加速策略批量評估將所有粒子的目標(biāo)函數(shù)評估合并為一次矩陣運(yùn)算def evaluate_batch(self, x: torch.Tensor) - torch.Tensor: 批量評估所有粒子 # x shape: [n_particles, n_dim] efficiency 0.9 - 0.001*x[:,2] - 0.0005*x[:,3] 0.00001*(x[:,0]-x[:,1]) cost 0.1*(x[:,0]**2 - x[:,1]**2) 5*x[:,3] 50 weight 0.01*(x[:,0]**2 - x[:,1]**2) 0.05*x[:,3] return torch.stack([efficiency, cost, weight], dim1)內(nèi)存優(yōu)化使用原地操作減少內(nèi)存分配self.velocities.mul_(self.w).add_( self.c1 * r1 * (self.pbest - self.particles), alpha1 ).add_( self.c2 * r2 * (leaders - self.particles), alpha1 )5.2 參數(shù)調(diào)優(yōu)指南參數(shù)推薦范圍影響調(diào)整策略粒子數(shù)50-500探索能力問題維度越高需要越多粒子慣性權(quán)重w0.4-0.9平衡探索與開發(fā)線性遞減效果更好學(xué)習(xí)因子c1,c21.5-2.5收斂速度c1略大于c2增強(qiáng)多樣性存檔大小50-200解的質(zhì)量越大Pareto前沿越完整變異概率0.05-0.2避免早熟初期可設(shè)較高值6. 常見問題與解決方案6.1 收斂問題排查問題1算法過早收斂到局部Pareto前沿檢查存檔中解的分布是否過于集中解決增加變異概率、減小慣性權(quán)重、使用動(dòng)態(tài)參數(shù)調(diào)整問題2Pareto前沿不完整檢查存檔大小是否足夠解決增大存檔容量、增加粒子數(shù)量、延長迭代次數(shù)6.2 數(shù)值穩(wěn)定性問題問題目標(biāo)函數(shù)量綱差異導(dǎo)致某些目標(biāo)被忽略# 標(biāo)準(zhǔn)化處理 def normalize_objectives(self, objs: torch.Tensor) - torch.Tensor: 將各目標(biāo)函數(shù)值歸一化到[0,1]范圍 min_vals objs.min(dim0)[0] max_vals objs.max(dim0)[0] return (objs - min_vals) / (max_vals - min_vals 1e-8)6.3 并行計(jì)算優(yōu)化使用PyTorch的分布式包實(shí)現(xiàn)多GPU訓(xùn)練import torch.distributed as dist def init_process(rank, world_size): dist.init_process_group(gloo, rankrank, world_sizeworld_size) # 分割粒子到不同GPU particles_per_gpu self.n_particles // world_size start rank * particles_per_gpu end (rank 1) * particles_per_gpu if rank ! world_size-1 else self.n_particles # 各GPU處理自己的粒子 local_particles self.particles[start:end] local_velocities self.velocities[start:end] # 處理完成后同步結(jié)果 dist.all_gather(all_particles, local_particles)