A list,
Particle swarm optimization (PSO) is an algorithm based on swarm search to deal with optimization problems in continuous or discrete space. Particle swarm optimization algorithm uses particles constantly moving in the solution space as the optimization group, each particle has two properties of position and velocity (the dimension of position and velocity is the same as the dimension of space), the position of the particle represents a feasible solution, and the velocity represents the difference between the next feasible solution found. Each particle adjusts its speed (according to a certain rule) to find a better solution according to the optimal solution it has already found and the optimal solution the group has currently found. Pso’s basic process initialization constant and each record matrix: CostFunction calculation function CostFunction; Parameter number of cost function nVar; [VarMin,VarMax] Maximum number of iterations iter_max Number of particle populations nPop particle velocity boundary [VelMin,VelMax] Single particle optimal acceleration constant C1 Global optimal acceleration constant C2(generally C1 = C2 ∈[0,4] C1=C2\in[0,4]C1=C2∈[0,4]) inertia factor W Particle position recording matrix PopPosition(nPop row nVar column, PopVelocity(nPop row nVar column, each row is the velocity component of a particle) Single particle optimal position record matrix Pbest(nPop row nVar column, Each row is a particle historically optimal position coordinates) global historically optimal position record matrix Gbest(a row of nVar columns, only record the optimal solution) historically optimal value record matrix (iter_max row column, the optimal solution value of each generation) initializes the particle position and speed: randomly initializes the position of each particle in the particle solution space; The initial velocity of each particle is initialized randomly in the range of each velocity component of the particle. Calculate the cost function value corresponding to the current position of each particle, and compare with the historical optimal solution of the particle, update the Pbest matrix; At the same time, compare with Gbest, update Gbest value; After updating the Pbest matrix and Gbest matrix, press the formula to update the particle speed: V I +1 = WV I + C 1R 1(P I − X I) + C 2 R 2 (G − X I) V_{I +1} = WV_i + C_1R_1(p_i-x_i) + C_2R_2(G-x_i)Vi+1=WVi+C1R1(Pi−Xi)+C2R2(G−Xi) where, W is the inertia factor, representing the ability of the particle to maintain the original speed constant; C1 is the single-particle acceleration constant, which represents the particle’s ability to solve its own historical optimal acceleration. C2 global optimal acceleration constant, which represents the ability of the particle to accelerate to the current global optimal solution; R1 and R2 are two random numbers between 0 and 1. After calculating the velocity, be careful to limit the velocity boundary. Once the velocity is greater than the boundary value, make it equal to the boundary value. After updating the speed, the next particle position should be calculated and saved as the current particle position: X I +1=X I +V I +1 X_{I +1} = X_i + V_{I +1}Xi+1=Xi+Vi+1 Similarly, the boundary should be limited after calculating the position. Loop iteration and output loop iteration calculate cost function and update speed and position, reach the maximum algebra can output results. When using object-oriented language to implement the algorithm, we can define the particle object, which should have attributes: current position, current speed, current generation value, historical optimal position, historical optimal generation value, and should have methods: speed update method, position update method, generation value update method.
Ii. Source code
% tic clear CLC %% r101= importData ('r101.txt');
cap=200; % Maximum vehicle load %% Extracted data information E= R101 (1.5); % Distribution center time window start time L=r101(1.6); % End time of distribution center time window vertexs=r101(:,2:3); % customer=vertexs(2:end,:); Cusnum =size(customer,1); % Number of customers v_num=25; % Max Volume = R101 (2:end,4); % demand a=r101(2:end,5); % customer time window start time [a[I],b[I]] b=r101(2:end,6); % Customer time window end time [a[I],b[I]] s=r101(2:end,7); % customer service time h=pdist(vertexs); dist=squareform(h); % distance matrix %% Particle swarm optimization parameter alpha=10; % Penalty function coefficient belta= for capacity constraint violations100; % Penalty function coefficient NIND= for violating the time window constraint50; % Number of particles MAXGEN=100; % Number of iterations w=1; % inertia factor wdamp=0.99; % inertia factor decay rate C1 =1.5; % Individual learning factor c2=2.0; % Global learning factor XvMin=1; % Xv floor XvMax = v_num; % Xv ceiling XrMin =1; % XrMax = cusnum Xr lower limit; The upper limit of % Xr VvMin = - (v_num- 1); % Vv floor VvMax = v_num- 1; The upper limit of % Vv VrMin = - (cusnum- 1); % Vr floor VrMax = cusnum- 1; Init_vc =init(Cusnum, A, Demands,cap); % to construct initial solution population = initpopCW (init_vc NIND, cusnum, XvMin, XvMax, XrMin, XrMax, VvMin, VvMax, VrMin, VrMax); ObjV=calObj(population,v_num,cap,demands,a,b,L,s,dist,alpha,belta); % Calculate the objective function value of each particle gbest_pos=population{1.1} (1:2, :); % assume that the first particle position is globally optimal gbest_obj=ObjV(1.1); % pbest_pos=cell(NIND,1); % initialize the individual optimal position of each particle pbest_obj=ObjV; % initializes the individually optimal objective function value for each particlefor i=1:NIND
particle=population{i,1}; % th particle position=particle(1:2, :); % position of the ith particle pbest_pos{I,1}=position; % initializes the individual optimum of this particleif pbest_obj(i,1)<gbest_obj % updates the globally optimal particle gbest_obj in the initial population=pbest_obj(i,1); gbest_pos=position; end end globalVC=decode(gbest_pos,v_num); % Initial global optimal distribution scheme % Count the total cost, the number of vehicles used, the total distance traveled, the number of routes violating constraints and the number of customers violating constraints of a distribution scheme [cost,NV,TD,vio_NV,vio_cus]=statistic(globalVC,a,b,s,L,dist,demands,cap,alpha,belta); disp([The total cost of the initial global optimal solution is:,num2str(cost),', the number of vehicles used is:,num2str(NV),...
', the total distance travelled is:,num2str(TD),', the number of paths violating the constraint is:,num2str(vio_NV),', the number of customers violating the restriction is:,num2str(vio_cus)]);
BestCost=zeros(MAXGEN,1); % records the globally optimal objective function value of each iteration %% main loop gen=1; % counterwhileGen <=MAXGEN %% updates the position and speed of each particlefor i=1:NIND
particle=population{i,1}; % th particle position=particle(1:2, :); % position of the ith particle Xv=position(1, :); Xr=position(2, :); velocity=particle(3:4, :); % the velocity of the ith particle Vv=velocity(1, :); Vr=velocity(2, :); %% update velocity=w*velocity+ +c1*rand([2,cusnum]).*(pbest_pos{i,1}-position)...
+c2*rand([2,cusnum]).*(gbest_pos-position); %% Speed out of bounds1,:)=max(velocity(1,:),VvMin);
velocity(1,:)=min(velocity(1,:),VvMax);
velocity(2,:)=max(velocity(2,:),VrMin);
velocity(2,:)=min(velocity(2,:),VrMax); Position =position+velocity; position(1, :) =ceil(position(1:)); IsOutside=(position(1,:)<XvMin | position(1,:)>XvMax | position(2,:)<XrMin | position(2,:)>XrMax); velocity(IsOutside)=-velocity(IsOutside); %% position out of bounds handle position(1,:)=max(position(1,:),XvMin);
position(1,:)=min(position(1,:),XvMax);
position(2,:)=max(position(2,:),XrMin);
position(2,:)=min(position(2,:),XrMax); %% relocate the distribution scheme decoded by the I particle VC=decode(position,v_num); RVC=relocate(VC,cap,demands,a,b,L,s,dist,alpha,belta); position=change(RVC,cusnum); % % calculation of the ith particle cost objective function value = costFuction (RVC, a, b, s, L, dist, demands, cap, alpha, belta); %% update individual optimalif cost<pbest_obj(i,1)
pbest_pos{i,1}=position;
pbest_obj(i,1)=cost; %% updates global optimalityif pbest_obj(i,1)<gbest_obj
gbest_pos=pbest_pos{i,1};
gbest_obj=pbest_obj(i,1); End end %% update particle=[position;velocity]; population{i,1}=particle; End %% record the global optimal solution of each generation and print BestCost(gen)=gbest_obj; globalVC=decode(gbest_pos,v_num); % Initial global optimal distribution scheme % Count the total cost, the number of vehicles used, the total distance traveled, the number of routes violating constraints and the number of customers violating constraints of a distribution scheme [cost,NV,TD,vio_NV,vio_cus]=statistic(globalVC,a,b,s,L,dist,demands,cap,alpha,belta); disp(['the first',num2str(gen),The total cost of the global optimal solution of generation is:,num2str(cost),', the number of vehicles used is:,num2str(NV),...
', the total distance travelled is:,num2str(TD),', the number of paths violating the constraint is:,num2str(vio_NV),', the number of customers violating the restriction is:,num2str(vio_cus)]); w=w*wdamp; %% update counter gen=gen+1; End %% decode global optimal particle into global optimal distribution scheme globalVC=decode(gBEST_POS, V_NUM); % % of global optimal distribution scheme for local search globalVC = relocate_gbest (globalVC, cap, demands, a, b, L, s, dist, alpha, belta); %% Counted the total cost, the number of vehicles used, the total distance traveled, the number of routes violating constraints and the number of customers violating constraints of the global optimal distribution scheme [cost,NV,TD,vio_NV,vio_cus]=statistic(globalVC,a,b,s,L,dist,demands,cap,alpha,belta); disp(['The total cost of the final global optimal solution is:',num2str(cost),', the number of vehicles used is:,num2str(NV),...
', the total distance travelled is:,num2str(TD),', the number of paths violating the constraint is:,num2str(vio_NV),', the number of customers violating the restriction is:,num2str(vio_cus)]); Draw_Best (globalVC, Vertexs); %% Results figure; %plot(BestCost,'LineWidth'.2);
semilogy(BestCost,'LineWidth'.2);
xlabel('Number of iterations');
ylabel('Globally optimal objective function value');
grid on;
toc
Copy the code3. Operation results
Fourth, note
Version: 2014 a