pso_opt.m

基于粒子群优化与灰色系统理论的多目标优化程序

clear;
global  f1 f0 kg kexample nn kdegree nn0;
global chaoV chaoV1; 
global keyH GreyEntropy;
fprintf('\n\t欢迎使用灰色混合离散变量－改进粒子群——混沌复合优化算法软件GS_CH_PSO1.0！\n')
fprintf('\t\t设计：湖南文理学院机械工程系 罗佑新\n')
fprintf('\t\t资助课题：湖南省科技厅计划项目资助\n')
fprintf('\t\t版权所有，非法复制传播必究。2007/08于常德\n')
fprintf(' Gs_CH_PSO1.0采用动态惩罚因子。\n')
fprintf('适应数学模型：    min f(x)   (x为n维，x=[x离散,x连续])\n')
fprintf('                s.t. gi(x)<=0  (i=1,2,……，gk);\n')
fprintf('%                   hj(x)==0  (j=1,2,……，hk).\n') 
fprintf('% 如何使用混沌还要参考其它文献，效率有待提高！.\n') 
fprintf('\n\t程序已经开始运行！\n\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%
m2=30;m=50;c1=2;c2=2;wmax=1.88;wmin=0.8;
vmax=3;kmax=10000;
%%%%%%%%%%%%%%%%%%%%%%%%%%
chaoV=1;%使用混沌chaoV=1;不使用混沌chaoV=0;
chaoV1=1;%emigrant_chao(us,ux,n,k)中变量；
keyH=1;%keyH=1;改进的罚函数（哈尔滨工业大学2006（10）：1716-1718相应动态罚函数）；
       %keyH=0;%未改进的罚函数；
       %注意： 用f(x)单目标计算即kg=3 时,取keyH=0（即用未改进的罚函数就可以了）！！！。
kg=1;% 用关联度计算;kg=3 用f(x)单目标计算。否则另用程序opt_fminimax即kg=2，用fminimax计算
kexample=7;%第一个类例题kexample=1；其余类推。
GreyEntropy=2;%GreyEntropy=1,使用灰熵计算；GreyEntropy=0；不使用灰熵计算；GreyEntropy=2 优度计算。
if chaoV==1,BB='使用混沌变量';else BB='不使用混沌变量';end
if keyH==1,CC=',改进的罚函数（哈尔滨工业大学2006（10）：1716-1718相应动态罚函数）';else, CC=',未改进的罚函数';end
%注意： 用f(x)单目标计算即kg=3 时,取keyH=0（即用未改进的罚函数就可以了）！！！。
m1=1;% 独立运行m1次
kk=1;%kk=1% 哈尔滨工业大学2006（10）：1716-1718相应动态罚函数(好）;
     %kk=2% 采用精确不可微罚函数,动态惩罚因子(引用郭惠昕开发的相应程序修改成gi(x)<=0 )（不好！！！！）!!
kdegree=7;
% kdegree=1 % absolute degree of grey incidence%%%%%%%%%%%%%%%%%
% kdegree=2 % Deng Type degree of grey incidencer%（有待改进）
% kdegree=3 % relative degree of grey incidenc(与其它例题相比例2最快）
% kdegree=4 % compositive  degree of grey incidenc
% kdegree=5 % improved relative degree of grey incidenc
% kdegree=6 %improved absolue degree of grey incidenc
% kdegree=7 % improved  compositive  degree of grey incidenc
% kdegree=8 %B type degree of grey incidenc
% kdegree=9 %C type degree of grey incidenc（效果不佳）
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 nn=1;nn0=2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if kexample==1, %%%%%%%%%%%%%%%%%% the first example(good)%%%%%%%%%%%%%%%%%%%%%%%%
    n=5;ux=[78,33,27,27,27]; us=[102,45,45,45,45];
    flag_x=[0,0,0,0,0];
    dis_date=[0.2,0.25,0.3,0.4,0.5,0.6,0.8,1,1.25,1.5,2,2.5,3,4,5,6,8,10,12,16,20,25,32,40,50];
    float_date=[2,2,2,2,2]; %各连续变量均取2位小数。
 elseif kexample==2    %%%%%%%%%%%%%%%%%%%% the second example (not find)%%%%%%%%%%%%%%%%
   n=6;
   ux=[10,17,0.2,17,10,13]; us=[20,130,0.8,30,15,20];
   flag_x=[1,1,2,0,1,1];
    dis_date=[0.2,0.25,0.3,0.4,0.5,0.6,0.8,1,1.25,1.5,2,2.5,3,4,5,6,8,10,12,16,20,25,32,40,50];
    float_date=[2,2,2,2,2,2]; %各连续变量均取2位小数。
 elseif kexample==3 
    n=7;ux=[2.6,0.3,17,7.3,7.3,2.9,5]; us=[3.6,1.0,28,8.3,8.3,3.9,5.5];
    flag_x=[0,0,0,0,0,0,0];
    dis_date=[0.2,0.25,0.3,0.4,0.5,0.6,0.8,1,1.25,1.5,2,2.5,3,4,5,6,8,10,12,16,20,25,32,40,50];
    float_date=[2,2,2,2,2,2,2]; %各连续变量均取2位小数。
elseif kexample==4
    n=3; ux=[3,18,4]; us=[7,100,10];
    flag_x=[0,1,1];
    dis_date=[0.2,0.25,0.3,0.4,0.5,0.6,0.8,1,1.25,1.5,2,2.5,3,4,5,6,8,10,12,16,20,25,32,40,50];
    float_date=[2,2,2]; %各连续变量均取2位小数。
elseif kexample==5
    n=4; ux=[10,10,0.9,0.9]; us=[80,50,5,5];
    flag_x=[0,0,0,0];
    dis_date=[0.2,0.25,0.3,0.4,0.5,0.6,0.8,1,1.25,1.5,2,2.5,3,4,5,6,8,10,12,16,20,25,32,40,50];
    float_date=[2,2,2,2]; %各连续变量均取2位小数。
elseif kexample==6
    n=5; ux=[0.6,0.6,0.6,0.6,0.6]; us=[2.4,2.4,2.4,2.4,2.4];
    flag_x=[0,0,0,0,0];
    dis_date=[0.2,0.25,0.3,0.4,0.5,0.6,0.8,1,1.25,1.5,2,2.5,3,4,5,6,8,10,12,16,20,25,32,40,50];
    float_date=[2,2,2,2,2];
else %kexample==7
   n=3; ux=[2.5,30,3]; us=[9.5,90,30];
    flag_x=[2,1,2];
    dis_date=[0.2,0.25,0.3,0.4,0.5,0.6,0.8,1,1.25,1.5,2,2.5,3,4,4.5,5,5.5,6,6.5,7,9,10,11,12,14,16,20,25,32,40,50;
              2,3,4,5,6,7,8,8.5,9,9.5,10,11,12,13,13.5,14,14.5,15,16,17,18,19,20,22,25,28,30,32,40,45,50];
    float_date=[2,2,2]; 
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555555
LB=ux;UB=us; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
for run_times=1:m1 % 独立运行m1次
   % 随机生成m个粒子。
   for i=1:m
       for j=1:n
         x(i,j)=LB(j)+(UB(j)-LB(j))*rand;v=0.5*rand(m,n);  
       end
       x(i,:)=discrete(x(i,:),n,flag_x,dis_date,float_date); %离散化处理
   end
   % 计算适应度值
   for i=1:m
     if kk==1
         f(i)= Pso_Opt_fun_1(x(i,:),kmax,LB);
     else
         f(i)= Pso_Opt_fun_2(x(i,:),kmax,kmax);
     end
   end
   % 找出个体极值和全局极值
   pbx=x;pbf=f;[gbf i]=min(pbf);gbx=pbx(i,:);w=wmax;
   for i=1:m
    v(i,:)=w*v(i,:)+c1*rand*(pbx(i,:)-x(i,:))+c2*rand*(gbx-x(i,:));
      for j=1:n
        if v(i,j)>vmax,v(i,j)=vmax;  elseif v(i,j)<-vmax, v(i,j)<-vmax; end
      end
      x(i,:)=x(i,:)+v(i,:);
      for j=1:n
        if x(i,j)<LB(j),x(i,j)=LB(j);elseif x(i,j)>UB(j),x(i,j)=UB(j);end;%%%%%%%%%%%%%%%%%%
      end
      x(i,:)=discrete(x(i,:),n,flag_x,dis_date,float_date); %离散化处理
   end
   % 开始循环迭代
  for k=1:kmax
     w=wmax*exp(-20*(k/kmax)^10);
     %if (k/100-floor(k/100))==0 & run_times==1,fprintf('%5.0f次迭代已完成',k); end
     waitbar(k/kmax);
     %w=wmax-k*(wmax-wmin)/kmax;
     for i=1:m
       if kk==1
         f(i)= Pso_Opt_fun_1(x(i,:),k,gbx);
       else
         f(i)= Pso_Opt_fun_2(x(i,:),k,kmax);
       end
     end
     for i=1:m
        if f(i)<pbf(i)
            pbf(i)=f(i);
            pbx(i,:)=x(i,:);
        end
     end
     [gbf i]=min(pbf);
     gbx=pbx(i,:);
     for i=1:m
         v(i,:)=w*v(i,:)+c1*rand*(pbx(i,:)-x(i,:))+c2*rand*(gbx-x(i,:));
        for j=1:n
            if v(i,j)>vmax
                v(i,j)=vmax;
            elseif v(i,j)<-vmax
                v(i,j)=-vmax;
            end
        end
        x(i,:)=x(i,:)+v(i,:);
      %%%%%%%%%%%%%%%%%%%% 完成混沌变量矩阵。%%%%%%%%%%%%% 
          if k>=m2 
                if (k/m2-floor(k/m2))==0,
                    if chaoV==1
                        x(i,:)=emigrant_chao(UB,LB,n,m);%disp('使用混沌');
                    %else  disp('不使用混沌');
                    end
                end
          end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for j=1:n
            if x(i,j)<LB(j),x(i,j)=LB(j);elseif x(i,j)>UB(j),x(i,j)=UB(j);end;%%%%%%%%%%%%%%%%%%
        end
         x(i,:)=discrete(x(i,:),n,flag_x,dis_date,float_date); %离散化处理
        %%%%%%%%%%%%%%%%%nn0可调%%%%%%%%%%%%
       for i=1:m,FF(i,:)=f1; end
       if nn>nn0
          [Ntemp1,Ntemp2]=size(FF);
         for j=1:Ntemp2
           ftemp1=FF(:,j);ff0=f0(j,1);
           for i=1:Ntemp1
             if ff0>ftemp1(i,1),ff0=ftemp1(i,1); end
           end 
            f0(j,1)=ff0;
         end
       end
     end
  end
  %%%%%%%%%%%%%%%%%%%
  format long 
  H=Pso_Opt_fun_1(gbx,1,gbx);
  F=Pso_Object(gbx); [Gx,Hx]=Pso_restrain(gbx);G=Gx' ;
  %disp(['\n']);
  if m1==1
    disp(['GS_CH_PSO1.0程序第',num2str(run_times),'次运行结果为(',BB,CC,')：']);
    if kg==3,
      disp(['(混合离散变量－改进粒子群——混沌复合优化算法单目标计算)']);
    elseif kg==1
      disp(['(灰色混合离散变量－改进粒子群——混沌复合优化算法多目标计算']);
        if kdegree==1,EE='absolute degree of grey incidence';end
      if kdegree==2,EE='Deng Type degree of grey incidencer,要多运行几次';end
      if kdegree==3,EE='Relative degree of grey incidenc';end
      if kdegree==4,EE='Compositive  degree of grey incidenc';end
      if kdegree==5,EE='Improved relative degree of grey incidenc';end
      if kdegree==6,EE='Improved absolue degree of grey incidenc';end
      if kdegree==7,EE='Improved  compositive  degree of grey incidenc';end
      if kdegree==8,EE='B type degree of grey incidenc,要多运行几次';end
      if kdegree==9,EE='C type degree of grey incidenc,要多运行几次';end
    else
      disp(['kg=2,用fminimax另编程计算']);
    end
    if GreyEntropy==2;
        DD='GreyEntropy=2 优度+灰关联计算';
    elseif GreyEntropy==1
      DD='GreyEntropy=1,使用灰熵灰关联计算';
    else
      DD='GreyEntropy=0；不使用灰熵计算,仅用灰关联'
    end
    if kg==1,disp([DD,EE]);end
  end
  disp(['第',num2str(kexample),'题，第',num2str(run_times),'次计算！']);
  disp(['设计变量上界=[',num2str(UB),']，设计变量下界=[',num2str(LB),']']);
  disp(['设计变量X=[',num2str(gbx),']']);
  disp(['动态罚函数：H=',num2str(H)]);
  disp(['检验：目标函数F=',num2str(F)]);
  disp(['各种目标函数:f1=[',num2str(f1'),']']);
  disp(['各种目标函数最好值:f0=[',num2str(f0'),']'])
  disp(['约束函数：G=[',num2str(G),']']);
  end
%%%%%%%%%%%%%%%%%%%%%2007.8.22%%%%%%%%%%%%%%%%%%