请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 2293|回复: 59

26粒子群算法的Pareto多目标函数优化

  [复制链接]

570

主题

667

帖子

3394

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3394
发表于 2017-2-7 21:08:31 | 显示全部楼层 |阅读模式
粒子群算法的Pareto多目标函数优化:
百度网盘视频链接:https://pan.baidu.com/s/1miVJFCG密码:5dga

具体链接在halcom.cn论坛,联系人QQ:3283892722
该论坛是一个学习交流平台,我会逐一的和大家分享学习。
欢迎大家录制视频,并提交给我,我来设置视频,你可在论坛进行打赏分享。
视频专用播放器:http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1
  1. % 基于PSO算法的多目标函数寻优
  2. clc,clear,close all
  3. warning off
  4. format longG
  5. % web('halcom.cn')
  6. global nObj
  7. nObj = 2;   % 2个目标
  8. nVar = 2;   % 2个未知量
  9. % PSO 参数
  10. c1 = 1.4995;  
  11. c2 = 1.4995;
  12. Vmin = -1;      % 粒子最大速度
  13. Vmax = 1;       % 粒子最小速度
  14. popmin1 = -1;  popmax1 = 1; % x1
  15. popmin2 = -1;  popmax2 = 1; % x2
  16. nVar=2;         % 未知量个数
  17. maxiter = 50;   % 迭代次数
  18. sizepop = 100;  % 种群数量
  19. repository = 50; % 非支配解种群数量
  20. w = 0.5;         % 初始权重
  21. wdamp=0.99;      % 权重衰减因子
  22. nGrid=7;            % 每个目标函数解集域的扩展网格大小,单周期数
  23. alpha=0.1;          % 非支配解的解集的扩展因子
  24. beta=2;             % 非支配解的选择因子
  25. gamma=2;            % 非支配解的淘汰因子

  26. % 初始化种群
  27. x.Pos=[];  % 未知量的解
  28. x.V=[];    % 粒子速度
  29. x.fitness=[];         % 粒子适应度值
  30. x.Best.Pos=[];        % 最优解
  31. x.Best.fitness=[];    % 最优解对应的适应度值--目标函数值
  32. x.IsDominated=[];        % 是否非支配解,支配解=1, 非支配解=0
  33. x.GridIndex=[];          % 非支配解的网格化参数值-索引总循环数(并排目标函数的序列化内)
  34. x.GridSubIndex=[];       % 每一个非支配解的网格化索引值(扩展取值范围内)
  35. pop=repmat(x,sizepop,1); % 初始化的种群
  36. clear x
  37. for i=1:sizepop
  38.     pop(i).Pos(1) = unifrnd(popmin1,popmax1,1);
  39.     pop(i).Pos(2) = unifrnd(popmin2,popmax2,1);
  40.     pop(i).V = zeros(1,nVar);
  41.     % 适应度函数
  42.     pop(i).fitness = fun(pop(i).Pos);
  43.     % 种群个体最优
  44.     pop(i).Best.Pos=pop(i).Pos;
  45.     % 种群最优适应度值
  46.     pop(i).Best.fitness=pop(i).fitness;
  47. end

  48. % 确定种群之间是否存在支配解
  49. pop=JudgePopDomination(pop);
  50. % 然后获取非支配解
  51. rep=pop(~[pop.IsDominated]);
  52. % 扩展非支配解的适应度函数取值范围
  53. Grid=CreateGrids(rep,nGrid,alpha);
  54. % 计算扩展非支配解的GridIndex
  55. for i=1:numel(rep)
  56.     rep(i)=FindGridIndex(rep(i),Grid);
  57. end

  58. %% AWPSO算法的多目标函数寻优主程序
  59. for i=1:maxiter
  60.     for j=1:sizepop
  61.         % 轮盘赌算子随机选择一个非支配解,扰动解,避免陷入局部最优
  62. % zbest=Selectzbest(rep,beta);
  63.         zbest=rep(1);
  64.         % 粒子速度更新
  65.         pop(j).V = w*pop(j).V +c1*rand(1,nVar).*(pop(j).Best.Pos-pop(j).Pos) ...
  66.             +c2*rand(1,nVar).*(zbest.Pos-pop(j).Pos);
  67.         
  68.         % 粒子位置更新
  69.         pop(j).Pos = pop(j).Pos + 0.5*pop(j).V;
  70.         
  71.         % 粒子取值范围约束
  72.         pop(j).Pos(1) = max(pop(j).Pos(1), popmin1);
  73.         pop(j).Pos(1) = min(pop(j).Pos(1), popmax1);
  74.         pop(j).Pos(2) = max(pop(j).Pos(2), popmin2);
  75.         pop(j).Pos(2) = min(pop(j).Pos(2), popmax2);
  76.         
  77.         % 适应度函数值
  78.         pop(j).fitness = fun(pop(j).Pos);

  79.         % 是否受支配
  80.         if Domination(pop(j),pop(j).Best)
  81.             pop(j).Best.Pos=pop(j).Pos;
  82.             pop(j).Best.fitness=pop(j).fitness;
  83.         else
  84.             if rand<0.5
  85.                 pop(j).Best.Pos=pop(j).Pos;
  86.                 pop(j).Best.fitness=pop(j).fitness;
  87.             end
  88.         end
  89.     end
  90.    
  91.     % 增加非支配解
  92.     rep =[rep; pop(~[pop.IsDominated])];
  93.     % 非支配解是否有支配解
  94.     rep=JudgePopDomination(rep);
  95.     % 删除支配解
  96.     rep=rep(~[rep.IsDominated]);
  97.    
  98.     % 扩展非支配解的适应度函数取值范围
  99.     Grid=CreateGrids(rep,nGrid,alpha);
  100.     for k=1:numel(rep)
  101.         rep(k)=FindGridIndex(rep(k),Grid);
  102.     end
  103.    
  104.     % 非支配解是否达到最大的非支配解的种群数
  105.     if numel(rep)>repository
  106.         % 将超过的种群删除
  107.         Extra=numel(rep)-repository;
  108.         for k=1:Extra
  109.             rep=DeleteOneRepMemebr(rep,gamma);
  110.         end
  111.     end

  112.      w=w*wdamp;  % 惯性权重
  113.      
  114.     figure(1);
  115.     Plotfitness(pop,rep);
  116.     pause(0.01);

  117. end
  118. % % 确定种群之间是否存在支配解
  119. % pop=JudgePopDomination(pop);
  120. % % 增加非支配解
  121. % rep =[rep; pop(~[pop.IsDominated])];
  122. % rep = uniqueRep(rep);
  123. figure(1);
  124. Plotfitness(pop,rep);

  125. % 最优结果
  126. for i=1:size( rep,1 )
  127.     zbest_sol(i,:) = rep(i).Pos;
  128.     if i==1
  129.         fprintf('\n')
  130.         disp('PSO多目标计算的第一组结果如下:')
  131.         disp(['x1 = ',num2str(zbest_sol(i,1))])
  132.         disp(['x2 = ',num2str(zbest_sol(i,2))])
  133.         disp(['最优适应度值为  ',num2str([rep(i).fitness]')])
  134.         fprintf('\n')
  135.     else
  136.         disp('PSO多目标计算结果有多个Pareto解集')
  137.     end
  138. end
复制代码
适应度函数:
  1. function fitness = fun(x)
  2. y1 = (x(1)-0.1)^2+(x(2)-0.9)^2;
  3. y2 = -20*exp(-0.2*sqrt((x(1)^2+x(2)^2)/2))-exp((cos(2*pi*x(1))+cos(2*pi*x(2)))/2)+20+2.71289;
  4. fitness = [y1;y2];
复制代码
适应度区间扩展:
  1. function Grid=CreateGrids(rep,nGrid,alpha)
  2. global nObj
  3. F = [rep.fitness];
  4. Fmin = min(F,[],2);
  5. Fmax = max(F,[],2);

  6. % 适应度值范围扩展
  7. dF=Fmax-Fmin;
  8. Fmin=Fmin-alpha*dF;
  9. Fmax=Fmax+alpha*dF;

  10. for i=1:nObj
  11.     Grid(i).LB = [];
  12.     Grid(i).UB = [];
  13. end

  14. for i=1:nObj
  15.    Fi = linspace( Fmin(i),Fmax(i), nGrid+1 );
  16.    Grid(i).LB = [-inf, Fi];
  17.    Grid(i).UB = [Fi, +inf];
  18. end
复制代码
删除一个Pareto解
  1. function rep=DeleteOneRepMemebr(rep,gamma)

  2.     % 所有的非支配解节点索引数
  3.     GridIndex=[rep.GridIndex];
  4.     % 唯一出现的数值,从小到大排序
  5.     unique_GridIndex=unique(GridIndex);
  6.     % 初始化,唯一出现该数值的次数
  7.     N=zeros(size(unique_GridIndex));
  8.     for k=1:numel(unique_GridIndex)
  9.         N(k)=numel(find(GridIndex==unique_GridIndex(k)));
  10.     end
  11.     % 选择概率
  12.     P=exp(gamma*N);
  13.     P=P/sum(P);
  14.     % 轮赌法随机选择一个节点
  15.     rand_Index=RandWheelSelection(P);
  16.     % 获取对应节点的节点索引数
  17.     rand_GridIndex=unique_GridIndex(rand_Index);
  18.     % 在GridIndex中寻找rand_GridIndex所在的索引值
  19.     index=find(GridIndex==rand_GridIndex);
  20.     % 在【1到numel(index)】之间随机生成一个整数
  21.     integer=randi([1 numel(index)]);
  22.     % zbest
  23.     rep(index(integer)) = [];
复制代码
支配关系:
游客,如果您要查看本帖隐藏内容请回复
区分节点索引:
  1. function rep=FindGridIndex(rep,Grid)
  2. global nObj    % 目标函数个数
  3. nGrid=numel(Grid(1).LB);
  4. rep.GridSubIndex=zeros(1,nObj);
  5. for i=1:nObj
  6.     rep.GridSubIndex(i) = find( rep.fitness(i)<Grid(i).UB,1,'first');
  7. end

  8. % 对GridIndex计算更新--每个目标函数的索引值
  9. rep.GridIndex = rep.GridSubIndex(1);
  10. for i=2:nObj
  11. %     rep.GridIndex = nGrid*( rep.GridIndex-1 ) + rep.GridSubIndex(i);
  12.     rep.GridIndex = rep.GridIndex + rep.GridSubIndex(i);
  13. end
复制代码
判断种群的支配关系:
  1. function pop = JudgePopDomination(pop)
  2. sizepop = numel(pop);  % 种群数量
  3. for i=1:sizepop
  4.     pop(i).IsDominated=0;
  5. end

  6. for i=1:sizepop-1
  7.     for j=i+1:sizepop
  8.         b1 = Domination(pop(i),pop(j));
  9.         b2 = Domination(pop(j),pop(i));
  10.         if b1==1
  11.             pop(j).IsDominated=1;
  12.         end
  13.         if b2==1
  14.             pop(i).IsDominated=1;
  15.         end
  16.     end
  17. end
复制代码
轮赌法:
  1. function rand_Index=RandWheelSelection(P)

  2. K = rand;
  3. Pcum = cumsum(P);

  4. rand_Index = find( K<Pcum,1,'first');
复制代码
画图函数:
  1. function Plotfitness(pop,rep)
  2.     % 删除不合适的解集
  3.     flag=[];
  4.     for i=1:size(pop,1)
  5.         if(pop(i).fitness(1) == 1e6 && pop(i).fitness(2) == 1e6)
  6.             flag = [flag, i];
  7.         end
  8.     end
  9.     pop(flag) = [];
  10.    
  11.     % 支配解
  12.     pop_fitness=[pop.fitness];
  13.     plot(pop_fitness(1,:),pop_fitness(2,:),'ko');
  14.     hold on;
  15.     popBest_fitness=[pop.Best];
  16.     popBest_fitness=[popBest_fitness.fitness];
  17.     plot(popBest_fitness(1,:),popBest_fitness(2,:),'b.');
  18.     % 非支配解
  19.     rep_fitness=[rep.fitness];
  20.     plot(rep_fitness(1,:),rep_fitness(2,:),'r*');
  21.    
  22.     xlabel('1^{st} Objective');
  23.     ylabel('2^{nd} Objective');
  24.    
  25.     grid on;
  26.     hold off;

  27. end
复制代码


回复

使用道具 举报

0

主题

1

帖子

12

积分

新手上路

Rank: 1

积分
12
发表于 2017-12-28 22:52:54 | 显示全部楼层
顶一下,新手刚开始学。。。。。。。
回复 支持 1 反对 0

使用道具 举报

3

主题

28

帖子

93

积分

注册会员

Rank: 2

积分
93
发表于 2017-3-10 20:24:13 | 显示全部楼层
谢谢楼主,求看帖子
回复 支持 反对

使用道具 举报

0

主题

3

帖子

12

积分

新手上路

Rank: 1

积分
12
发表于 2017-3-12 18:05:50 | 显示全部楼层
谢谢分享,看看子程序,继续学习
回复 支持 反对

使用道具 举报

0

主题

2

帖子

10

积分

新手上路

Rank: 1

积分
10
发表于 2017-3-14 10:07:39 | 显示全部楼层
感谢楼主的分享
回复 支持 反对

使用道具 举报

0

主题

3

帖子

10

积分

新手上路

Rank: 1

积分
10
发表于 2017-3-15 19:54:37 | 显示全部楼层
楼主好人。代码学习一下,粒子群算法,赞一个
回复 支持 反对

使用道具 举报

0

主题

3

帖子

10

积分

新手上路

Rank: 1

积分
10
发表于 2017-3-16 19:07:10 | 显示全部楼层
感谢楼主的分享
回复 支持 反对

使用道具 举报

0

主题

18

帖子

58

积分

注册会员

Rank: 2

积分
58
发表于 2017-3-20 15:55:17 | 显示全部楼层

谢谢分享,继续学习
回复 支持 反对

使用道具 举报

0

主题

50

帖子

112

积分

注册会员

Rank: 2

积分
112
发表于 2017-3-24 23:07:04 | 显示全部楼层
粒子群算法多目标函数优化,学习。
回复 支持 反对

使用道具 举报

0

主题

4

帖子

9

积分

新手上路

Rank: 1

积分
9
发表于 2017-3-28 11:04:00 | 显示全部楼层
求看子函数
回复 支持 反对

使用道具 举报

0

主题

11

帖子

34

积分

新手上路

Rank: 1

积分
34
发表于 2017-4-13 21:12:47 | 显示全部楼层
学习一下,主要想要找到约束条件的处理方法
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则


Python|Opencv|MATLAB|Halcom.cn  

GMT+8, 2018-2-22 22:33 , Processed in 0.103769 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表