%Generic Algorithm for function f(x1,x2) optimum
clear all;
close all;

%Parameters
Size=50;     %群体大小
G=500;       %终止进化代数

E=round(rand(Size,28));    %Initial Code 28位 11的横坐标 11的纵坐标 6的回收单价

%% 下面这一小段不在任何函数里面，是公开的
AA1= xlsread('RC207.xlsx');
AA1(: , 1 ) = [];                %去掉第一列
% AA1( 1, : )=[];          %去掉第一行
model.CustX = AA1(:, 1)  ;        %顾客X坐标
model.CustY = AA1(:, 2) ;         %顾客Y坐标
model.weight= AA1(:,3);           %顾客持有纸壳重量
model.Custprice = AA1( :,4) ;       %顾客纸壳价钱心理预期
model.numCustomer  = numel( model.CustX  ); %顾客数量
model.numCenter=1;%回收中心数量（暂定1）
%%
%Main Program
for k=1:1:G
    
    time(k)=k; %记录代数，用于输出
    particle=[];
    for s=1:1:Size
%         m=E(s,:);
%         y1=0;y2=0;
        
        %Uncoding 解码
%         index=1;
%         for i=1:6:6*2*CodeL
%             number=m(i)*2^5+m(i+1)*2^4+m(i+2)*2^3+m(i+3)*2^2+m(i+4)*2^1+m(i+5)*2^0;
            
          for i=1:28
          particle(s).position(i)=E(s,i);
          end
%             index=index+1;
%         end
        [particle(s).profit,particle(s).Solution]=Model(particle(s).position,model);
        F(s)=particle(s).profit;

    end
    


    %****** Step 1 : Evaluate BestJ ******
    BestJ(k)=max(F);    
    
    fi=F;                          %Fitness Function  F是[1*12]的适应度数组
    [Oderfi,Indexfi]=sort(fi);     %Arranging fi small to bigger [Oderfi,Indexfi]前面是size长的小到大排序，后面是size长的其原本序号，也就是在E中第几（s）行
    Bestfi=Oderfi(Size);           %Let Bestfi=max(fi)  最大值
    BestS=E(Indexfi(Size),:);      %Let BestS=E(m), m is the Indexfi belong to max(fi)  种群E里面的最大值的粒子码
    Bestparticle=particle(Indexfi(Size)); 
    BestSolution=particle(Indexfi(Size)).Solution;
    bfi(k)=Bestfi;               %本代最大适应度值
    %先将上一代传来的种群E进行适应度计算、记载，在进行下面的各项操作
    
    %****** Step 2 : Select and Reproduct Operation******
       fi_sum=sum(fi); %求种群适应度之和
       fi_Size=(Oderfi/fi_sum)*Size;  
       
       fi_S=floor(fi_Size);        %Selecting Bigger fi value 排好的Oderfi按照floor筛选变成一串00001111，0被剔除，1被保留
       
       kk=1;
       for i=1:1:Size %对于排序后的每一个粒子，其fi_S要么为0要么为1 ，下一排代码fi_S=0就是不操作，1就是操作
          for j=1:1:fi_S(i)        %Select and Reproduce 
           TempE(kk,:)=E(Indexfi(i),:);  %TempE就是把E从小到大排了并取后几个保存，Indexfi(i)代表着初始粒子在粒子群中的位置
             kk=kk+1;              %kk is used to reproduce
          end
       end
       for i=kk:Size
           TempE(i,:)=E(i,:);
       end
           
       
    %************ Step 3 : Crossover Operation ************
    pc=0.60;
    n=ceil(20*rand); %n是1-20随机整数
    for i=1:2:(Size-1) %选中所有单数排，用于和后一排交叉
        temp=rand;
        if pc>temp                  %Crossover Condition
        for j=n:1:20
            TempE(i,j)=TempE(i+1,j);
            TempE(i+1,j)=TempE(i,j);
        end
        end
    end
    TempE(Size,:)=BestS; %保证适应度最大的那个粒子不被改变
    E=TempE;
       
    %************ Step 4: Mutation Operation **************
    %pm=0.001;
    %pm=0.001-[1:1:Size]*(0.001)/Size; %Bigger fi, smaller Pm
    %pm=0.0;    %No mutation
    pm=0.1;     %Big mutation
    
       for i=1:1:Size
          for j=1:1:28
             temp=rand;
             if pm>temp                %Mutation Condition
                if TempE(i,j)==0
                   TempE(i,j)=1;
                else
                   TempE(i,j)=0;
                end
            end
          end
       end
       
    %Guarantee TempPop(30,:) is the code belong to the best individual(max(fi))
    TempE(Size,:)=BestS; %保证适应度最大的那个粒子不被改变
    E=TempE;
   
    disp(['Iteration ' num2str(k) ': BestProfit = ' num2str(Bestfi)]);
end
 
for i=1:28
          Bestparticle.position(i)=TempE(Size,i);
end
        [Bestparticle.profit,Bestparticle.Solution]=Model(Bestparticle.position,model);
 disp(Bestparticle.position)
% disp('Solution')
% disp( num2str(Bestparticle.Solution.CenterX));
% disp( num2str(Bestparticle.Solution.CenterY));
%            disp( num2str(Bestparticle.Solution.price));
%      disp( num2str(Bestparticle.Solution.collectpercent));
%      disp( num2str(Bestparticle.Solution.collectnumber));
%      disp(num2str(Bestparticle.Solution.dealedpercent));
%      disp( num2str(Bestparticle.Solution.totaldistance));
%      disp( num2str(Bestparticle.Solution.carboncost));
%      disp( num2str(Bestparticle.Solution.recyclingcost));
%      disp( num2str(Bestparticle.Solution.totalweight));
%      disp( num2str(Bestparticle.Solution.totalincentive));
%      disp( num2str(Bestfi));
     disp('Solution')
      disp(['Solution.CenterX'  ' is:' num2str(Bestparticle.Solution.CenterX)]);
    disp(['Solution.CenterY'  ' is:' num2str(Bestparticle.Solution.CenterY)]);
      disp(['Recyling price'  ' is:' num2str(Bestparticle.Solution.price)]);
     disp(['Solution.collectpercent'  ' is:' num2str(Bestparticle.Solution.collectpercent)]);
     disp(['Solution.collectnumber'  ' is:' num2str(Bestparticle.Solution.collectnumber)]);
     disp(['Solution.dealedpercent'  ' is ' num2str(Bestparticle.Solution.dealedpercent)]);
     disp(['Solution.totaldistance'  ' is ' num2str(Bestparticle.Solution.totaldistance)]);
     disp(['Solution.carboncost'  ' is ' num2str(Bestparticle.Solution.carboncost)]);
     disp(['Solution.recyclingcost'  ' is ' num2str(Bestparticle.Solution.recyclingcost)]);
     disp(['Solution.totalweight'  ' is ' num2str(Bestparticle.Solution.totalweight)]);
     disp(['Solution.totalincentive'  ' is ' num2str(Bestparticle.Solution.totalincentive)]);
     disp(['Solution.BestProfit'  ' is ' num2str(Bestfi)]);
      %%
% 目标函数值收敛图
figure('NumberTitle', 'off', 'Name', '算法收敛曲线');
set(gcf,'Color',[1 1 1]);
%plot(BestCost,'LineWidth',2);
plot(bfi,'LineWidth',2);
xlabel('Number of Iterations','fontsize',13 );
ylabel('Value of Fitness','fontsize',13 );
box off ; grid off;

%尝试画二维图
figure('NumberTitle', 'off', 'Name', '回收情况二维图');
Colormattrix=rand(model.numCustomer,3);
for i=1:model.numCustomer
plot(model.CustX(i),model.CustY(i),'LineStyle','none','Marker','o','MarkerSize',8,'MarkerFace',[Colormattrix(i,:)],'LineWidth',2,'MarkerEdge',[rand(1,3)],'LineWidth',1);
hold on;
end
plot(Bestparticle.Solution.CenterX,Bestparticle.Solution.CenterY,'LineStyle','none','Marker','o','MarkerSize',10,'MarkerFace','g','LineWidth',3);
hold on;
for i=1:numel(Bestparticle.Solution.DealedcustX)
    h(i)=plot([Bestparticle.Solution.DealedcustX(i),Bestparticle.Solution.CenterX],[Bestparticle.Solution.DealedcustY(i),Bestparticle.Solution.CenterY]);
    set(h(i),'LineWidth',1.5);
    set(h(i),'Color',[Colormattrix(Bestparticle.Solution.Index(i),:)]);

    hold on;
end
hold off;
xlabel('x','FontSize',15);ylabel('y','FontSize',15)
%%
%以下是建模

function [profit,Solution]=Model(oneparticle,model)
% AA1= xlsread('customer.xlsx');
% AA1(: , 1 ) = [];                %去掉第一列
% 
% model.CustX = AA1(:, 1)  ;        %顾客X坐标
% model.CustY = AA1(:, 2) ;         %顾客Y坐标
% model.weight= AA1(:,3);           %顾客持有纸壳重量
% model.Custprice = AA1( :,4) ;       %顾客纸壳价钱心理预期
% model.numCustomer  = numel( model.CustX  ); %顾客数量
% model.numCenter=1;%回收中心数量（暂定1）
%回收中心坐标解码，回收价钱解码
model.CenterX=(oneparticle(1)*2^10+oneparticle(2)*2^9+oneparticle(3)*2^8+oneparticle(4)*2^7+oneparticle(5)*2^6+oneparticle(6)*2^5+oneparticle(7)*2^4+oneparticle(8)*2^3+oneparticle(9)*2^2+oneparticle(10)*2^1+oneparticle(11)*2^0)/10;
model.CenterY=(oneparticle(12)*2^10+oneparticle(13)*2^9+oneparticle(14)*2^8+oneparticle(15)*2^7+oneparticle(16)*2^6+oneparticle(17)*2^5+oneparticle(18)*2^4+oneparticle(19)*2^3+oneparticle(20)*2^2+oneparticle(21)*2^1+oneparticle(22)*2^0)/10;

% model.CenterX=oneparticle(1);
% model.CenterY=oneparticle(2);
model.price=(oneparticle(23)*2^5+oneparticle(24)*2^4+oneparticle(25)*2^3+oneparticle(26)*2^2+oneparticle(27)*2^1+oneparticle(28)*2^0)/10;
%顾客位置到收集点距离
model.Center2CustD=  zeros( model.numCustomer,model.numCenter )  ;
for i=1:  model.numCustomer  
    Q1= [  model.CenterX(1) model.CenterY(1)  ];
    Q2= [  model.CustX(i)   model.CustY(i) ];
    model.Center2CustD(i,1)=norm(     Q1-    Q2   );
    model.Center2CustD(1,i)=model.Center2CustD(i,1);
end

%参数设定
recyclingprice=2.0;
%  price=1.1;
carbonprice=0.1;
incentive=recyclingprice*0;
%初始化结算变量
totalweight=0; %总收集重量
totaldistance=0;%总收集距离
recyclingcost=0; %收掉客户手上纸壳花的钱
totalincentive=0;%总激励
%判断每个客户是否卖出纸壳
t=1;
DealedcustX=[];
DealedcustY=[];
Index=[];
for i=1:model.numCustomer
    if model.price>=model.Custprice(i)
        totalweight=totalweight+model.weight(i);
        totaldistance=totaldistance+model.Center2CustD(i,1)*2;
        recyclingcost=recyclingcost+model.weight(i)*model.price;
        totalincentive=totalincentive+incentive;
        DealedcustX(t)=model.CustX(i);
        DealedcustY(t)=model.CustY(i);
        Index(t)=i; %成交客户t是原来顾客矩阵的第i个
        t=t+1;
    end
%     if price>=model.Custprice(i)
%         totalweight=totalweight+model.weight(i);
%         totaldistance=totaldistance+model.Center2CustD(i,1)*2;
%         recyclingcost=recyclingcost+model.weight(i)*price;
%         DealedcustX(t)=model.CustX(i);
%         DealedcustY(t)=model.CustY(i);
%         t=t+1;
%     end
end
profit=totalweight*recyclingprice+totalincentive-totaldistance*carbonprice-recyclingcost;

Solution.price=model.price;%回收价格
Solution.collectpercent=totalweight/sum(model.weight);%回收重量占总比
Solution.collectnumber=numel(DealedcustX);%回收客户数目
Solution.dealedpercent=numel(DealedcustX)/numel( model.CustX  );%回收客户占总比
Solution.totaldistance=totaldistance;%总距离
Solution.carboncost=totaldistance*carbonprice;%与路程有关的花费
Solution.recyclingcost=recyclingcost;%回收花费
Solution.totalweight=totalweight;%总回收重量
Solution.totalincentive=totalincentive;%总激励
Solution.CenterX=model.CenterX;
Solution.CenterY=model.CenterY; 
Solution.DealedcustX=DealedcustX;%回收客户横坐标
Solution.DealedcustY=DealedcustY;
Solution.Index=Index;%记录成交的第几个客户是原来顾客矩阵的第几个
end