A list,
Simulated annealing algorithm
The famous simulated annealing algorithm, which is a kind of approximate optimization problem based on Monte Carlo design method.
A bit of history — skip it if you’re not interested
In 1953, American physicist N.Metropolis and his colleagues published an article on studying complex systems and calculating the energy distribution in them. They used monte Carlo simulation method to calculate the energy distribution of molecules in multi-subsystem. This is the beginning of the problem we’ll be exploring in this article. In fact, one of the terms commonly used in simulated annealing is the Metropolis criterion, which we’ll cover later.
In 1983, IBM physicists S.Kirkpatrick, C. D. Gelatt and M. P. Vecchi published an influential paper in Science: Optimization by Simulated Annealing. When they discuss a spin glass system using the methods of Metropolis et al., Find that the energy of physical systems is quite similar to the cost function of some combinatorial optimization problems (the famous traveling salesman problem TSP is a representative example) : seeking the lowest cost is like seeking the lowest energy. Therefore, they develop a set of algorithms based on Metropolis method and use it to solve combinatorial problems and seek optimal solutions.
At about the same time, the European physicist V.Carny published almost identical results, but they were discovered independently; It was just that Carny was “unlucky” that few people noticed his masterpiece at the time; Perhaps it is fair to say that while Science is widely known around the world, Carny published his work in another specialized academic journal, J.opt.Theory Appl., with a small circulation, which did not attract the attention it deserved.
Inspired by the Monte Carlo simulations used by Metropolis et al., Kirkpatrick coined the term “simulated annealing” because it resembles the process of annealing objects. Finding the optimal solution (maximum value) of the problem is similar to finding the lowest energy of the system. So as the system cools, the energy drops, and in the same sense, the solution to the problem “drops” to the maximum.
A, what is annealing – physical origin
In thermodynamics, annealing refers to the physical phenomenon that the object gradually cools down. The lower the temperature is, the lower the energy state of the object will be. When it is low enough, the liquid begins to condense and crystallize, and in the crystallized state, the energy state of the system is lowest. In slow cooling (i.e., annealing), nature can “find” the lowest energy state: crystallization. However, if the process is quenched too fast, quenching, or quenching, can result in amorphous particles that are not the lowest energy state.
As shown below, first (left) the object is in an amorphous state. We heat the solid to a high enough temperature (middle image) and then allow it to cool, or to annealing (right image). When the solid is heated, the internal particles become disordered and their internal energy increases. When the solid is cooled slowly, the particles become orderly and reach the equilibrium state at every temperature. Finally, they reach the ground state at room temperature and their internal energy decreases to a minimum (at this time, the object is in the form of crystal).
Nature, it seems, knows how to do it slowly: cool things down so that the molecules have enough time to settle down at each temperature, and gradually, eventually, the lowest energy state is the most stable system.
2. Simulate Anneal
If you’re still confused about the physics of annealing, there’s a simpler way to understand it. Imagine that we now have a function like this, and we want to find the (global) optimal solution for the function. If you take the Greedy strategy, you start at point A and continue the heuristic if the function continues to decrease. By the time we reach point B, it’s obvious that the quest is over (because it just keeps getting bigger no matter what direction we go). So we’re going to have to find a local final solution, B.
Simulated annealing is also an Greedy algorithm, but its search process introduces a random element. The simulated annealing algorithm accepts a solution worse than the current one with a certain probability, so it is possible to jump out of the local optimal solution and reach the global optimal solution. In the above figure, for example, the simulated annealing algorithm will accept to move to the right with a certain probability after finding the local optimal solution B. Maybe after a few of these non-locally optimal moves you get to a peak point between B and C, and you get out of the local minimum B.
According to the Metropolis criterion, the probability of a particle reaching equilibrium at temperature T is exp(- δ E/(kT)), where E is the internal energy at temperature T, δ E is its change, and K is Boltzmann constant. The Metropolis criterion is often expressed as
The Metropolis criterion states that at temperature T, the probability of a cooling with an energy difference of dE is P(dE), expressed as P(dE) = exp(dE/(kT)). Where k is a constant, exp represents the natural exponent, and dE<0. So P is positively correlated with T. This formula says: The higher the temperature is, the greater the probability of a first cooling with an energy difference of dE; The lower the temperature, the less likely there is to be a cooling. And since dE is always less than 0 (because annealing is a gradual decrease in temperature), dE/kT is less than 0, so P(dE) is in the range of (0,1). As T goes down, P(dE) goes down.
We regard a move to the worse solution as a temperature jump, and we accept such a move with probability P(dE). In other words, when solid annealing is used to simulate combinatorial optimization problem, internal energy E is simulated as objective function value F and temperature T is evolved into control parameter T, namely, simulated annealing algorithm to solve combinatorial optimization problem is obtained: Starting from the initial solution I and the initial value t of the control parameter, the current solution is repeated with the iteration of “generating a new solution → calculating the difference of the objective function → accepting or discarding”, and the value of T is gradually attenuated. The current solution at the end of the algorithm is the approximate optimal solution, which is a heuristic random search process based on monte Carlo iterative solution. The annealing process is controlled by the Cooling Schedule, including the initial value T of the control parameters and its attenuation factor δ T, the number of iterations L at each value t and the stopping condition S.
To sum it up:
If f(Y(I +1)) <= f(Y(I)) (that is, the better solution is obtained after the move), the move is always accepted; If f(Y(I +1)) > f(Y(I)) (that is, the solution after the move is worse than the current solution), then the move is accepted with a certain probability, and this probability decreases gradually with time (gradually decreases to become stable), which is equivalent to the small wave peak between B and BC in the figure above, The probability of each move to the right (that is, accepting a worse value) is gradually decreasing. If the slope is extremely long, there’s a good chance we won’t get over it. If it’s not too long, it’s probably going to flip over it, depending on how you set the attenuation t.
Here’s an interesting analogy between the common Greedy algorithm and simulated annealing:
Ordinary Greedy algorithm: The rabbit jumps lower than it is now. It found the lowest valley not far away. But this valley is not necessarily the lowest. This is the ordinary Greedy algorithm, which does not guarantee that the local optimum is the global optimum.
Simulated annealing: The rabbit gets drunk. It jumps randomly for a long time. During this time, it may go low, or it may step flat. However, it gradually woke up and jumped in the lowest direction. This is simulated annealing.
Vehicle Routing Problem (VRP)
Vehicle routing problem (VRP) was first proposed by Dantzig and Ramser in 1959. It refers to a certain number of customers with different demand for goods. Distribution center provides goods to customers, a fleet is responsible for distributing goods, and appropriate driving routes are organized. And under certain constraints, to achieve such as the shortest distance, the least cost, the least time and so on.
Vehicle Routing Problem with Time Windows (VRPTW)
Due to the continuous development of VRP problem, considering that demand points have requirements on vehicle arrival Time, the vehicle routing problem with Time window is added into the vehicle routing problem, which becomes VRP with Time Windows (VRPTW). The Vehicle Routing problem with Time Windows (VRPTW) is a time window constraint for the customer to be visited on the VRP. In the VRPTW problem, in addition to the driving cost, the cost function also includes the waiting time caused by arriving at a customer early and the service time required by the customer. In VRPTW, vehicles must not only meet the constraints of VRP problems, but also meet the Time Window constraints of demand points. The Time Window constraints of demand points can be divided into two types, one is Hard Time Window, which requires that vehicles must arrive within the Time Window. Early arrival must wait, while late arrival is rejected. The other is the Soft Time Window. It is not necessary to arrive in the Time Window, but to arrive outside the Time Window must be punished. The biggest difference between the Soft Time Window and the hard Time Window is to replace waiting and rejection with punishment.
Model 1 Problem definition:
The VRPTW is defined by a fleet of vehicles K={1,… A set of customers C={1,… ,n} , and a directed graph G, Typically the fleet is considered to be homogeneous, that is, All vehicles are identical. The graph consists of | | C + 2 are, where The customers are denoted 1, 2,… ; 20 performers by the vertices 0 (” the starting depot “) and N + 1 (” the returning depot “) All vertices, that is, 0,1… , n+1 is denoted N. The set of arcs, A, represents direct connections between the depot and the customers and among the customers. There are no arcs ending at Vertex 0 or originating from vertex n + 1. With each arc (I,j), where I ≠j, we associate a cost cij and a time tij, which may include service time at customer i.
Each vehicle has a capacity Q and each customer i a demand qi . Each customer i has a time window [ai,bi] and a vehicle must arrive at the customer before bi . If it arrives before the time window opens, it has to wait until ai to service the customer. The time windows for both depots are assumed to be identical to [a0,b0] which represents the scheduling horizon. The vehicles may not leave the depot before a0 and must return at the latest at time bn+1 .
It is assumed that Q, ai , bi , qi , cij are non-negative integers and tij are positive integers. Note that this assumption is necessary to develop an algorithm for the shortest path with resource constraints used in the column generation approach presented later. Furthermore it is assumed that the triangle inequality is satisfied for both cij and tij .
The decision variable s ik is defined for each vertex i and each vehicle k and denotes the time vehicle k starts to service customer i. In case vehicle k does not service customer i, s ik has no meaning and consequently it’s value is considered irrelevant. We assume a0 = 0 and therefore s 0k = 0, for all k.
Model 2(refer to 2017 A generalized Formulation for Vehicle Routing Problems) :
The model is a 2-dimensional decision variable
Iii. Source code
Tic clear CLC %% tic clear CLC %%'c101.txt');
cap=200; %% Extract data information E=c101(1.5); % Distribution center time window start time L=c101(1.6); % End time of distribution center time window vertexs=c101(:,2:3); % customer=vertexs(2:end,:); Cusnum =size(customer,1); % Number of customers v_num=5; % Max Number of vehicles = C101 (2:end,4); % demand a=c101(2:end,5); % Customer time window start time [a[I],b[I]] b=c101(2:end,6); % Customer time window end time [a[I],b[I]] s=c101(2:end,7); % customer service time h=pdist(vertexs); dist=squareform(h); % distance matrix %% simulated annealing parameter Belta =10; % Penalty function coefficient for capacity constraint violations gama=100; % Penalty function coefficient for violation of time window constraint MaxOutIter=2000; % Maximum number of iterations of the outer loop MaxInIter=300; % Maximum number of iterations of inner loop T0=100; % Initial temperature alpha=0.99; % Cooling factor pSwap=0.2; % Probability of selecting the swap structure pReversion=0.5; % Probability of selecting reverse structure pInsertion=1-pSwap-pReversion; % Select the probability of inserting structure N=cusnum+v_num- 1; % solution length = number of customers + maximum number of vehicles used- 1The initial solution currS= Randperm (N) is constructed randomly. % random structure initial solution [currVC, NV, TD, violate_num violate_cus] = decode (currS cusnum, cap, demands, a, b, L, s, dist); % To decode the initial solution % to find the cost of the initial distribution scheme = the total cost of the vehicle + the sum of the capacity constraints violated by Belta * + the sum of time window constraints violated by GAMA * currCost=costFuction(currVC,a,b,s,L,dist,demands,cap,belta,gama); Sbest=currS; % Initially assign the global optimal solution to the initial solution bestVC=currVC; % The global optimal distribution scheme is initially assigned as initial distribution scheme bestCost=currCost; % Initially assign the total cost of the global optimal solution to the total cost of the initial solution BestCost= Zeros (MaxOutIter,1); % Record the total cost of the global optimal solution for each generation T=T0; % temperature initialization %% simulated annealingfor outIter=1:MaxOutIter
for inIter=1:MaxInIter newS=Neighbor(currS,pSwap,pReversion,pInsertion); % % after neighborhood structure of the new solution for the new decoding = + belta * against the cost of the vehicle capacity constraints, the sum of + of gama * in violation of the sum of the time window constraints newVC = decode (newS, cusnum, cap, demands, a, b, L, s, dist); newCost=costFuction(newVC,a,b,s,L,dist,demands,cap,belta,gama); If the new solution is better than the current solution, update the current solution, and the total cost of the current solutionif newCost<=currCost
currS=newS;
currVC=newVC;
currCost=newCost;
elseIf the new solution is not as good as the current solution, the annealing criterion is adopted to accept the new solution with a certain probability delta=(newcost-currcost)/currCost; % Calculate the percentage difference between the total cost of the new solution and the current solution P=exp(-delta/T); % computes the probability of accepting the new solution % if0~1Is less than P, the new solution is accepted, and the current solution is updated, as well as the total cost of the current solutionifrand<=P currS=newS; currVC=newVC; currCost=newCost; End end % Compares the current solution with the global optimal solution. If the current solution is better, the global optimal solution is updated, as well as the total cost of the global optimal solutionifcurrCost<=bestCost Sbest=currS; bestVC=currVC; bestCost=currCost; End end % Record the total cost of the global optimal solution of each iteration of the outer loop BestCost(outIter)= BestCost; % displays the total cost of the global optimal solution for each iteration of the outer loop disp(['the first',num2str(outIter),Global optimal solution of generation:]) [bestVC,bestNV,bestTD,best_vionum,best_viocus]=decode(Sbest,cusnum,cap,demands,a,b,L,s,dist); % decoding of global optimal solution disp(['Number of Vehicles used:',num2str(bestNV),', total distance travelled by vehicle:,num2str(bestTD),', number of paths violating constraints: ',num2str(best_vionum),', number of customers violating constraints: ',num2str(best_viocus)]);
fprintf('\n'T=alpha*T; End %% Print the total cost variation trend figure of the global optimal solution for each iteration of the outer loop; plot(BestCost,'LineWidth'.1);
title('Trend chart of total cost of global Optimal solution')
xlabel('Number of iterations');
ylabel('Total cost'); Draw_Best (bestVC,vertexs); tocCopy the code3. Operation results
Fourth, note
Version: 2014 a