1, Introduction to particle swarm optimization
Particle swarm optimization algorithm was proposed by Dr. Eberhart and Dr. Kennedy in 1995. It comes from the study of bird predation behavior. Its basic core is to make use of the information sharing of individuals in the group, so as to make the movement of the whole group produce an evolutionary process from disorder to order in the problem-solving space, so as to obtain the optimal solution of the problem. Imagine this scenario: a flock of birds are foraging, and there is a corn field in the distance. All the birds don't know where the corn field is, but they know how far their current position is from the corn field. Then the best strategy to find the corn field, and the simplest and most effective strategy, is to search the surrounding area of the nearest bird group to the corn field.
In PSO, the solution of each optimization problem is a bird in the search space, which is called "particle", and the optimal solution of the problem corresponds to the "corn field" found in the bird swarm. All particles have a position vector (the position of the particle in the solution space) and a velocity vector (which determines the direction and speed of the next flight), and the fitness value of the current position can be calculated according to the objective function, which can be understood as the distance from the "corn field". In each iteration, the examples in the population can learn not only according to their own experience (historical position), but also according to the "experience" of the optimal particles in the population, so as to determine how to adjust and change the flight direction and speed in the next iteration. In this way, the whole population will gradually tend to the optimal solution.
The above explanation may be more abstract. Let's illustrate it with a simple example
There are two people in a lake. They can communicate with each other and detect the lowest point of their position. The initial position is shown in the figure above. Because the right side is deep, the people on the left will move the boat to the right.
Now the left side is deep, so the person on the right will move the boat to the left
Keep repeating the process and the last two boats will meet
A local optimal solution is obtained
Represent each individual as a particle. The position of each individual at a certain time is expressed as x(t), and the direction is expressed as v(t)
p (T) is the optimal solution of individual x at time t, g(t) is the optimal solution of all individuals at time t, v(t) is the direction of individual at time t, and x(t) is the position of individual at time t
The next position is shown in the figure above, which is determined by X, P and G
The particles in the population can find the optimal solution of the problem by constantly learning from the historical information of themselves and the population.
However, in the follow-up research, the table shows that there is a problem in the above original formula: the update of V in the formula is too random, which makes the global optimization ability of the whole PSO algorithm strong, but the local search ability is poor. In fact, we need that PSO has strong global optimization ability in the early stage of algorithm iteration, and the whole population should have stronger local search ability in the later stage of algorithm iteration. Therefore, according to the above disadvantages, shi and Eberhart modified the formula by introducing inertia weight, so as to put forward the inertia weight model of PSO:
The components of each vector are represented as follows
Where W is called the inertia weight of PSO, and its value is between [0,1]. Generally, adaptive value method is adopted, that is, w=0.9 at the beginning, which makes PSO have strong global optimization ability. With the deepening of iteration, the parameter W decreases, so that PSO has strong local optimization ability. When the iteration is over, w=0.1. The parameters c1 and c2 are called learning factors and are generally set to 14961; r1 and r2 are random probability values between [0,1].
The algorithm framework of the whole particle swarm optimization algorithm is as follows:
Step 1 population initialization can be carried out randomly or design a specific initialization method according to the optimized problem, and then calculate the individual fitness value, so as to select the individual local optimal position vector and the global optimal position vector of the population.
Step 2 iteration setting: set the number of iterations and set the current number of iterations to 1
step3 speed update: update the speed vector of each individual
step4 position update: update the position vector of each individual
Step 5 local position and global position vector update: update the local optimal solution of each individual and the global optimal solution of the population
Step 6 termination condition judgment: the maximum number of iterations is reached when judging the number of iterations. If it is satisfied, the global optimal solution is output. Otherwise, continue the iteration and jump to step 3.
For the application of particle swarm optimization algorithm, it is mainly the design of velocity and position vector iterative operators. Whether the iterative operator is effective or not will determine the performance of the whole PSO algorithm, so how to design the iterative operator of PSO is the research focus and difficulty in the application of PSO algorithm.
2, VRP model
(1) Introduction to vehicle path planning
After 60 years of research and development, the research objectives and constraints of vehicle routing problem have changed. It has changed from the initial simple vehicle scheduling problem to a complex system problem. The initial vehicle path planning problem can be described as: there is a starting point and several customer points. Given the geographical location and needs of each point, how to plan the optimal path under the conditions of meeting various constraints, so that it can serve each customer point, and finally return to the starting point. By imposing different constraints and changing the optimization objectives, different kinds of vehicle routing problems can be derived. At the same time, the vehicle routing problem belongs to a typical NP hard problem, and its exact algorithm can solve a small scale, so the heuristic algorithm has become a research hotspot.
(2) Introduction to VRPTW
VRPTW(Vehicle routing problem with time windows) is a vehicle routing problem with time window. It adds time window constraints to each demand point, that is, for each demand point, set the earliest and latest time of service start, and require the vehicle to start serving customers within the time window.
There are two kinds of time window restrictions at the demand point. One is the hard time window, which requires the vehicle to start serving customers within the time window, wait early and reject it if it is late. The other is the soft time window, which does not necessarily start serving customers within the time window, but must be punished if it starts serving outside the time window, Penalty instead of waiting and rejection is the biggest difference between soft time window and hard time window.
The mathematical model of VRPTW is as follows:
(2.2) ensure that each customer is only visited once
(2.3) ensure that the loaded goods do not exceed the capacity
(2.4) (2.5) (2.6) ensure that each vehicle starts from the depot and finally returns to the depot
(2.7) (2.8) ensure that the service starts within the time window
clc %Clear command line window clear %Deletes all variables from the current workspace and frees them from system memory close all %Delete all windows whose handles are not hidden tic % Save current time %% GA-PSO Algorithm solving CDVRP %Input: %City Longitude and latitude of demand point %Distance Distance matrix %Demand Demand at each point %Travelcon Travel constraint %Capacity Vehicle capacity constraint %NIND Population number %MAXGEN Inherited to the third MAXGEN Time substitution program stop %Output: %Gbest Shortest path %GbestDistance Shortest path length %% Load data load('../test_data/City.mat') %The longitude and latitude of the demand point is used to draw the actual path XY coordinate load('../test_data/Distance.mat') %Distance matrix load('../test_data/Demand.mat') %requirement load('../test_data/Travelcon.mat') %Travel constraint load('../test_data/Capacity.mat') %Vehicle capacity constraint %% Initialize problem parameters CityNum=size(City,1)-1; %Number of demand points %% Initialization algorithm parameters NIND=60; %Number of particles MAXGEN=100; %Maximum number of iterations %% 0 matrix initialized for pre allocated memory Population = zeros(NIND,CityNum*2+1); %Pre allocated memory for storing population data PopDistance = zeros(NIND,1); %Pre allocated matrix memory GbestDisByGen = zeros(MAXGEN,1); %Pre allocated matrix memory for i = 1:NIND %% Initialize each particle Population(i,:)=InitPop(CityNum,Distance,Demand,Travelcon,Capacity); %Random generation TSP route %% Calculate path length PopDistance(i) = CalcDis(Population(i,:),Distance,Demand,Travelcon,Capacity); % Calculate path length end %% storage Pbest data Pbest = Population; % initialization Pbest For the current particle collection PbestDistance = PopDistance; % initialization Pbest The objective function value of is the objective function value of the current particle collection %% storage Gbest data [mindis,index] = min(PbestDistance); %get Pbest in Gbest = Pbest(index,:); % initial Gbest particle GbestDistance = mindis; % initial Gbest Objective function value of particle %% Start iteration gen=1; while gen <= MAXGEN %% Update per particle for i=1:NIND %% Particles and Pbest overlapping Population(i,2:end-1)=Crossover(Population(i,2:end-1),Pbest(i,2:end-1)); %overlapping % If the length of the new path becomes shorter, record to Pbest PopDistance(i) = CalcDis(Population(i,:),Distance,Demand,Travelcon,Capacity); %Calculate distance if PopDistance(i) < PbestDistance(i) %If the length of the new path becomes shorter Pbest(i,:)=Population(i,:); %to update Pbest PbestDistance(i)=PopDistance(i); %to update Pbest distance end %% Update iterations gen=gen+1; end %Delete redundant 1 in path for i=1:length(Gbest)-1 if Gbest(i)==Gbest(i+1) Gbest(i)=0; %When the adjacent bits are all 1, the previous one is set to zero end end Gbest(Gbest==0)=; %Delete redundant zero elements Gbest=Gbest-1; % Subtract 1 from each code, which is consistent with the code in the text %% The calculation result data is output to the command line disp('-------------------------------------------------------------') toc %Display run time fprintf('Total Distance = %s km \n',num2str(GbestDistance)) TextOutput(Gbest,Distance,Demand,Capacity) %Show optimal path disp('-------------------------------------------------------------') %% Iterative graph figure plot(GbestDisByGen,'LineWidth',2) %Show the historical change of objective function value xlim([1 gen-1]) %set up x Axis range set(gca, 'LineWidth',1) xlabel('Iterations') ylabel('Min Distance(km)') title('HPSO Process') %% Draw actual route DrawPath(Gbest,City)