布谷鸟算法的Pareto多目标函数优化

Halcom

布谷鸟算法的Pareto多目标函数优化
https://pan.baidu.com/s/1G_XE4C9kUwYJk3le3L9PiQ
https://pan.baidu.com/s/1H1qHD1u6Bb1awo1-6gIh8w
具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1

1. % 布谷鸟的pareto多目标函数优化分析
   
2. % cuckoo_search -- CS算法
   
3. % web('halcom.cn')
   
4. %% 清空环境
   
5. clc % 清屏
   
6. clear all; % 删除workplace变量
   
7. close all; % 关掉显示图形窗口
   
8. tic
   
9. %% 参数初始化
   
10. global nObj nvar
    
11. nObj = 2;   % 2个目标
    
12. nvar = 2;   % 2个未知量
    
13. maxiter = 50;     % 最大迭代次数
    
14. sizepop = 100;     % 种群规模
    
15. popmin1 = -1;  popmax1 = 1; % x1
    
16. popmin2 = -1;  popmax2 = 1; % x2
    
17. pa=0.25;        % 发现新巢的概率Discovery rate of alien eggs/solutions
    
18. repository = 50;    % 非支配解种群数量
    
19. nGrid=7;            % 每个目标函数解集域的扩展网格大小，单周期数
    
20. alpha=0.1;          % 非支配解的解集的扩展因子
    
21. beta=2;             % 非支配解的选择因子
    
22. gamma=2;            % 非支配解的淘汰因子
    
23. % 初始化种群
    
24. x.Pos=[];  % 未知量的解
    
25. x.fitness=[];         % 粒子适应度值
    
26. x.IsDominated=[];        % 是否非支配解，支配解=1, 非支配解=0
    
27. x.GridIndex=[];          % 非支配解的网格化参数值-索引总循环数（并排目标函数的序列化内）
    
28. x.GridSubIndex=[];       % 每一个非支配解的网格化索引值（扩展取值范围内）
    
29. pop=repmat(x,sizepop,1); % 初始化的种群
    
30. clear x
    
31. for i=1:sizepop
    
32.     pop(i).Pos(1) = unifrnd(popmin1,popmax1,1);
    
33.     pop(i).Pos(2) = unifrnd(popmin2,popmax2,1);
    
34.     % 适应度函数
    
35.     pop(i).fitness = fun(pop(i).Pos);
    
36. end
    
37. % 确定种群之间是否存在支配解
    
38. pop=JudgePopDomination(pop);
    
39. % 然后获取非支配解
    
40. rep=pop(~[pop.IsDominated]);
    
41. % 扩展非支配解的适应度函数取值范围
    
42. Grid=CreateGrids(rep,nGrid,alpha);
    
43. % 计算扩展非支配解的GridIndex
    
44. for i=1:numel(rep)
    
45.     rep(i)=FindGridIndex(rep(i),Grid);
    
46. end
    
47. %% CS算法的多目标函数寻优主程序
    
48. for i=1:maxiter
    
49.    
    
50.     % 粒子位置更新
    
51.     new_nest=[];
    
52.     new_nest = get_cuckoos(pop,[popmin1,popmin2],[popmax1,popmax2], rep, beta); % Levy flights 粒子位置更新
    
53.     [pop, new_nest]=get_best_nest(pop, new_nest);
    
54. 
    
55.     % 增加非支配解
    
56.     rep =[rep; pop(~[pop.IsDominated])];
    
57.     % 非支配解是否有支配解
    
58.     rep=JudgePopDomination(rep);
    
59.     % 删除支配解
    
60.     rep=rep(~[rep.IsDominated]);
    
61.    
    
62.     % 寻找新的鸟巢--新的解
    
63.     new_nest=[];
    
64.     new_nest=empty_nests(pop, [popmin1,popmin2],[popmax1,popmax2], pa) ;
    
65.     [pop, new_nest]=get_best_nest(pop, new_nest);
    
66.    
    
67.     % 增加非支配解
    
68.     rep =[rep; pop(~[pop.IsDominated])];
    
69.     % 非支配解是否有支配解
    
70.     rep=JudgePopDomination(rep);
    
71.     % 删除支配解
    
72.     rep=rep(~[rep.IsDominated]);
    
73.     % 扩展非支配解的适应度函数取值范围
    
74.     Grid=CreateGrids(rep,nGrid,alpha);
    
75.     for k=1:numel(rep)
    
76.         rep(k)=FindGridIndex(rep(k),Grid);
    
77.     end
    
78.     % 非支配解是否达到最大的非支配解的种群数
    
79.     if numel(rep)>repository
    
80.         % 将超过的种群删除
    
81.         Extra=numel(rep)-repository;
    
82.         for k=1:Extra
    
83.             rep=DeleteOneRepMemebr(rep,gamma);
    
84.         end
    
85.     end
    
86.      
    
87.     figure(1);
    
88.     Plotfitness(pop,rep);
    
89.     pause(0.01);
    
90. end
    
91. toc
    
92. % % 确定种群之间是否存在支配解
    
93. % pop=JudgePopDomination(pop);
    
94. % % 增加非支配解
    
95. % rep =[rep; pop(~[pop.IsDominated])];
    
96. rep = uniqueRep(rep);
    
97. figure(1);
    
98. Plotfitness(pop,rep);
    
99. 
    
100. % 最优结果
     
101. for i=1:size( rep,1 )
     
102.     zbest_sol(i,:) = rep(i).Pos;
     
103.     if i==1
     
104.         fprintf('\n')
     
105.         disp('PSO多目标计算的第一组结果如下：')
     
106.         disp(['x1 = ',num2str(zbest_sol(i,1))])
     
107.         disp(['x2 = ',num2str(zbest_sol(i,2))])
     
108.         disp(['最优适应度值为  ',num2str([rep(i).fitness]')])
     
109.         fprintf('\n')
     
110.     else
     
111.         disp('PSO多目标计算结果有多个Pareto解集')
     
112.     end
     
113. end

复制代码

empty_nests函数：

1. function new_nest=empty_nests(nest, popmin, popmax, pa)
   
2. global nObj nvar
   
3. nest1 = [nest.Pos];
   
4. nest1 = reshape(nest1, [nvar, size(nest,1)])';
   
5. 
   
6. n=size(nest1,1);  % 种群个数
   
7. % 发现新巢的概率Discovery rate of alien eggs/solutions
   
8. K=rand(size(nest1))>pa;
   
9. stepsize = rand*( nest1(randperm(n),:) - nest1(randperm(n),:) );
   
10. new_nest1 = nest1 + stepsize.*K;
    
11. 
    
12. for i=1:n
    
13.     for j=1:size( new_nest1, 2 )
    
14.         if(new_nest1(i, j)>popmax(j))
    
15.             new_nest1(i, j) = popmax(j);
    
16.         end
    
17.         if(new_nest1(i, j)<popmin(j))
    
18.             new_nest1(i, j)=popmin(j);
    
19.         end
    
20.     end
    
21.     new_nest(i).Pos = new_nest1(i,:);
    
22. end

复制代码

get_best_nest函数：

1. function [nest, new_nest]=get_best_nest(nest, new_nest)
   
2. % nest      输入   上一代的解
   
3. % fun：     输入   适应度函数
   
4. % new_nest  输入   新产生的解
   
5. 
   
6. for i=1:size(nest,1)
   
7.     new_nest(i).fitness = fun( new_nest(i).Pos );  % 适应度值
   
8.    
   
9.     % 是否受支配
   
10.     if Domination(new_nest(i), nest(i))
    
11.         nest(i).Pos=new_nest(i).Pos;
    
12.         nest(i).fitness=new_nest(i).fitness;
    
13.     else
    
14.         if rand<0.5
    
15.             nest(i).Pos=new_nest(i).Pos;
    
16.             nest(i).fitness=new_nest(i).fitness;
    
17.         end
    
18.     end
    
19.    
    
20. end
    
21. % % 找当前最好的个体
    
22. % [bestfitness,bestindex]=min(fitness);  % 最优适应度值
    
23. % bestnest = nest(bestindex,:);   % 全局最佳种群

复制代码

get_cuckoos函数：

1. function new_nest = get_cuckoos(pop, popmin, popmax, rep, betaPareto)
   
2. % Levy flights 粒子位置更新
   
3. global nObj nvar
   
4. nest = [pop.Pos];
   
5. nest = reshape(nest, [nvar, size(pop,1)])';
   
6. new_nest = pop;
   
7. pop = [];
   
8. 
   
9. % Levy flights
   
10. n=size(nest,1);   % 种群个数
    
11. % Levy exponent and coefficient
    
12. beta=3/2;
    
13. sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
    
14. 
    
15. % new_nest = nest;
    
16. for i=1:n   % 种群个数
    
17.     pop = nest(i,:);
    
18.     % Levy flights by Mantegna's algorithm
    
19.     u=randn(size(pop))*sigma;
    
20.     v=randn( size(pop) );
    
21.     step=u./abs(v).^(1/beta);
    
22.    
    
23.     % 轮盘赌算子随机选择一个非支配解，扰动解，避免陷入局部最优
    
24.     zbest=Selectzbest(rep, betaPareto);
    
25.     bestnest = zbest.Pos;
    
26.    
    
27.     stepsize=0.01*step.*(pop-bestnest);
    
28.     pop = pop + stepsize.*randn(size(pop));
    
29.     % 防止越界
    
30.     for k=1:length( pop )
    
31.         if pop(k)>popmax(k)
    
32.             pop(k)=popmax(k);
    
33.         end
    
34.         if pop(k)<popmin(k)
    
35.             pop(k)=popmin(k);
    
36.         end
    
37.     end
    
38. %     new_nest(i,:) = pop;  % 更新种群
    
39.     new_nest(i).Pos = pop;   % 更新种群
    
40. end

复制代码

参考：
【1】http://halcom.cn/forum.php?mod=viewthread&tid=357&extra=page%3D1
【2】Solving Multiobjective Optimization Problems Using Artificial Bee Colony Algorithm

George

感谢楼主，学习学习

jackian

非常好的学习资源，感谢分享

jackian

机器码：FF48E712CD53F57C-0-FF48C202D90BF17C

风飞纤手

期待，学习一下~:lol