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%%%%%%%%%%%%%%%%%%