A list,
Simulated annealing algorithm The famous simulated annealing algorithm is a method of approximate solving optimization problems based on Monte Carlo design. A bit of history – if you’re not interested, skip the 1953 paper by THE American physicist N.Metropolis and his colleagues, who studied complex systems and calculated the energy distribution in them. They used the Monte Carlo simulation method to calculate the energy distribution of molecules in a multisubsystem. 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. In thermodynamics, annealing (rolling) refers to the physical phenomenon of the object gradually cooling, the lower the temperature, the energy state of the object will be low; 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.
2.1 Vehicle Routing Problem with Capacity Constraints (CVRP)
The model was difficult to extend to other scenarios of VRP and the execution paths of specific vehicles were not known, so the model continued to be improved.
Ii. Source code
Tic clear CLC %% tic clear CLC %% tic clear CLC %%'rc208.txt');
cap=1000; Vertexs =data(:,2:3); % customer=vertexs(2:end,:); Cusnum =size(customer,1); % Number of customers v_num=25; % Max =data(2:end,4); % demand h=pdist(vertexs); dist=squareform(h); % cost matrix %% simulated annealing parameter Belta =10; % Penalty function coefficient for capacity constraint violation 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, dist); % of initial solution decoding currCost = costFuction (currVC, dist, demands, cap, belta); % The cost of the initial distribution scheme = the total cost of the vehicle +alpha* the sum of the capacity constraints violated 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 to produce new solution newVC = decode (newS, cusnum, cap, demands, dist); % to decrypt data processing newCost = costFuction (newVC, dist, demands, cap, belta); % Find the cost of the initial distribution solution = total cost of vehicle travel +alpha* sum of capacity constraints violated % 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,dist);
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