【VRP问题】基于matlab模拟退火算法求解带容量的VRP问题（多车型）

原创于 2021-07-09 15:03:24 发布 · 27 阅读

CC 4.0 BY-SA版权

文章标签：

一、简介

模拟退火算法介绍【VRP问题】基于matlab模拟退火算法求解带容量的VRP问题(多车型)_模拟退火算法
3 模拟退火算法的参数
模拟退火是一种优化算法，它本身是不能独立存在的，需要有一个应用场合，其中温度就是模拟退火需要优化的参数，如果它应用到了聚类分析中，那么就是说聚类分析中有某个或者某几个参数需要优化，而这个参数，或者参数集就是温度所代表的。它可以是某项指标，某项关联度，某个距离等等。

二、源代码

%%
%

clc;
clear;
close all;

%% Problem Definition

model=SelectModel();        % Select Model of the Problem

model.eta=0.1;

CostFunction=@(q) MyCost(q,model);       % Cost Function

%% SA Parameters

MaxIt=500;     % Maximum Number of Iterations

MaxIt2=80;      % Maximum Number of Inner Iterations

T0=100;         % Initial Temperature

alpha=0.98;     % Temperature Damping Rate


%% Initialization

% Create Initial Solution
x.Position=CreateRandomSolution(model);
[x.Cost, x.Sol]=CostFunction(x.Position);

% Update Best Solution Ever Found
BestSol=x;

% Array to Hold Best Cost Values
BestCost=zeros(MaxIt,1);

% Set Initial Temperature
T=T0;


%% SA Main Loop

for it=1:MaxIt
    for it2=1:MaxIt2
        
        % Create Neighbor
        xnew.Position=CreateNeighbor(x.Position);
        [xnew.Cost, xnew.Sol]=CostFunction(xnew.Position);
        
        if xnew.Cost<=x.Cost
            % xnew is better, so it is accepted
            x=xnew;
            
        else
            % xnew is not better, so it is accepted conditionally
            delta=xnew.Cost-x.Cost;
            p=exp(-delta/T);
            
            if rand<=p
                x=xnew;
            end
            
        end
        
        % Update Best Solution
        if x.Cost<=BestSol.Cost
            BestSol=x;
        end
        
    end
    
    % Store Best Cost
    BestCost(it)=BestSol.Cost;

    
    % Display Iteration Information
    if BestSol.Sol.IsFeasible
        FLAG=' *';
    else
        FLAG='';
    end
    disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it)) FLAG]);
    
    % Reduce Temperature
    T=alpha*T;
    
    % Plot Solution
    figure(1);
    PlotSolution(BestSol.Sol,model);
    pause(0.01);
    
end
%
 

function model=CreateRandomModel(I,J)

    rmin=10;
    rmax=25;
    r=randi([rmin rmax],1,I);
    
    TotalDemand=sum(r);
    cmean=TotalDemand/J;
    cmin=round(cmean);
    cmax=round(1.25*cmean);
    c=randi([cmin cmax],1,J);
    
    xmin=0;
    xmax=200;
    
    ymin=0;
    ymax=100;
    
    x=randi([xmin xmax],1,I);
    y=randi([ymin ymax],1,I);
    
    alpha_x=0.1;
    xm=(xmin+xmax)/2;
    dx=xmax-xmin;
    x0min=round(xm-alpha_x*dx);
    x0max=round(xm+alpha_x*dx);
    
    alpha_y=0.1;
    ym=(ymin+ymax)/2;
    dy=ymax-ymin;
    y0min=round(ym-alpha_y*dy);
    y0max=round(ym+alpha_y*dy);
    
    x0=randi([x0min x0max]);
    y0=randi([y0min y0max]);
    
    d=zeros(I,I);
    d0=zeros(1,I);
    for i=1:I
        for i2=i+1:I
            d(i,i2)=sqrt((x(i)-x(i2))^2+(y(i)-y(i2))^2);
            d(i2,i)=d(i,i2);
        end
        
        d0(i)=sqrt((x(i)-x0)^2+(y(i)-y0)^2);
    end
    
    eta=0.5;
    
    model.I=I;
    model.J=J;
    model.r=r;
    model.c=c;
    model.xmin=xmin;
    model.xmax=xmax;
    model.ymin=ymin;
    model.ymax=ymax;
    model.x=x;
    model.y=y;
    model.x0=x0;
    model.y0=y0;
    model.d=d;
    model.d0=d0;
    model.eta=eta;

end
%


function qnew=CreateNeighbor(q)

    m=randi([1 3]);
    
    switch m
        case 1
            % Do Swap
            qnew=Swap(q);
            
        case 2
            % Do Reversion
            qnew=Reversion(q);
            
        case 3
            % Do Insertion
            qnew=Insertion(q);
    end

end

function qnew=Swap(q)

    n=numel(q);
    
    i=randsample(n,2);
    i1=i(1);
    i2=i(2);
    
    qnew=q;
    qnew([i1 i2])=q([i2 i1]);
    
end

function qnew=Reversion(q)

    n=numel(q);
    
    i=randsample(n,2);
    i1=min(i(1),i(2));
    i2=max(i(1),i(2));
    
    qnew=q;
    qnew(i1:i2)=q(i2:-1:i1);

end

function qnew=Insertion(q)

    n=numel(q);
    
    i=randsample(n,2);
    i1=i(1);
    i2=i(2);
    
    if i1<i2
        qnew=[q(1:i1-1) q(i1+1:i2) q(i1) q(i2+1:end)];
    else
        qnew=[q(1:i2) q(i1) q(i2+1:i1-1) q(i1+1:end)];
    end

end