%% PSO 粒子群优化算法 clc; clear; close all; close all %% Problem Definition 问题构建 %CostFunction=@(oneparticle)ModelforQ1(oneparticle); % 调用目标函数 oneparticle=[]; nVar =3; % Number of Decision Variables 变量个数，编码长度 VarSize=[1 nVar]; % Size of Decision Variables Matrix 决策变量的维度 VarMin= 0 ; % Lower Bound of Variables 变量 低阶 VarMax=100; % Upper Bound of Variables 变量高阶 rand('seed',sum(clock)); %% Parameters % PSO MaxIt=500; %model.q= Maximum Number of Iterations 算法最大迭代次数 nPop= 50; % Population Size (Swarm Size) 种群数目 wmax=0.9; % Inertia Weight 最大权重 wmin=0.1; Pr=0.9; Pns=0.6; c1=1.49618; % Personal Learning Coefficient 参数C1 c2=1.49618; % Global Learning Coefficient 参数C2 VelMax=0.1*(VarMax-VarMin); % 速度参数 VelMin=-VelMax; %% %% 下面这一小段不在任何函数里面，是公开的 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） %% Initialization 初始化 % empty_particle.Position=[]; % empty_particle.PositionT=[]; % empty_particle.Cost=[]; % empty_particle.Velocity=[]; % empty_particle.sol = [ ]; % empty_particle.Best.Position=[]; % empty_particle.Best.Cost=[]; % empty_particle.Best.sol=[ ]; % particle=repmat(empty_particle,nPop,1); GlobalBest.profit=-inf; GlobalBest.Position=-inf; GlobalBest.sol=-inf; for i=1:nPop % Initialize Position 初始化种群 particle(i).Position(1)=(rand( ) .* ( VarMax - VarMin ) + VarMin) ; particle(i).Position(2)=(rand( ) .* ( VarMax - VarMin ) + VarMin) ; particle(i).Position(3)=(rand() .*10) ; %要根据倒卖价钱定 % Initialize Velocity 初始化速度 particle(i).Velocity= ( rand( VarSize) .* ( VelMax - VelMin ) + VelMin ); % Evaluation 粒子评价 [particle(i).profit,particle(i).Solution] =Model(particle(i).Position,model); %? % Update Personal Best 更新粒子当前最优 particle(i).Best.Position=particle(i).Position; particle(i).Best.profit=particle(i).profit; particle(i).Best.Solution=particle(i).Solution; % Update Global Best 更新全局 最优 if particle(i).Best.profit>GlobalBest.profit GlobalBest.profit=particle(i).Best.profit; GlobalBest.Position=particle(i).Best.Position; GlobalBest.Solution=particle(i).Best.Solution; end warning ('off'); end Bestprofit=zeros(MaxIt,1); t0=cputime; %% PSO Main Loop for it=1:MaxIt w=wmax-(wmax-wmin)/MaxIt*it; for i=1:nPop %Pi(t) % Update Velocity 更新速度 particle(i).Velocity = (w*particle(i).Velocity ... +c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ... +c2*rand(VarSize).*(GlobalBest.Position-particle(i).Position)); % Apply Velocity Limits 速度限制 index = find(particle(i).Velocity(1,:)>VelMax | particle(i).Velocity(1,:)VarMax | particle(i).PositionNORMAL(1,:)particle(i).profit0 particle(i).profit=particle(i).profit; particle(i).Solution=particle(i).Solution; particle(i).Position=particle(i).Position; else particle(i).profit=particle(i).profit0; particle(i).Solution=particle(i).Solution0; particle(i).Position=particle(i).PositionNORMAL; end % if particle(i).Costparticle(i).Best.profit particle(i).Best.Position=particle(i).Position; particle(i).Best.profit=particle(i).profit; particle(i).Best.Solution=particle(i).Solution; end % Update Global Best 更新 种群的全局 最优 if particle(i).Best.profit>GlobalBest.profit GlobalBest=particle(i).Best; end end %% 邻域搜索策略 % for i=1:nPop % b=rand(0,1); % if bVarMax | particle(i).PositionLLL(1,:)VarMax | particle(i).PositionGGG(1,:)=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 totalincentive=totalweight/1000*incentive; 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