模拟退火是一种元启发式优化方法，适用于解决距离约束车辆路径问题（DVRP），目标是优化车辆路线以最小化总距离，同时确保每辆车不超过最大行驶距离。代码通过读取Excel数据、初始化随机解、迭代优化和可视化结果来实现该算法，研究表明初始解生成和Metropolis准则在算法性能中起关键作用。专注于距离约束而非容量约束，与传统VRP有所不同，可能存在潜在改进空间，如更复杂的降温计划或初始解生成方法。

数据加载与初始化
从Excel文件“data”中读取需求点坐标（City）、需求量（Demand，未使用）和距离限制（Travelcon）。计算距离矩阵（Distance）以确定各点间距离。

%% 加载数据-
City=xlsread('data','City');%需求点经纬度，用于画实际路径的XY坐标
Demand=xlsread('data','Demand');%需求
Travelcon=xlsread('data','Travelcon');%距离限制
Distance=CalculateDist(City);

初始解生成
生成随机TSP路线（TSProute），然后转换为满足距离约束的VRP路线（S1）。通过检查当前行驶距离是否超过Travelcon，决定是否返回配送中心（标记为1）并启动新车辆。

模拟退火循环
-从高初始温度（$T0=1000$）开始，逐步降温（$q=0.9$）至终止温度（$Tend=1e-3$）。每个温度下进行$L=200$次迭代，生成新解（NewSolution），并根据Metropolis准则决定是否接受新解。记录每次迭代的最优路线和目标值（总距离），更新最佳解。

输出与可视化
输出历史最短距离（mindisever）和对应路径（bestroute），清理冗余配送中心标记。绘制优化过程迭代图和最优解的实际路线图。

数学模型

距离约束车辆路径问题（DVRP）可形式化为以下优化问题：

$V$: 所有节点集，包括配送中心（0）和客户点。$C$: 客户点集，C=V∖{0}C = V\setminus \{0\}C=V∖{0}。

$d_{ij}$: 节点$i$到节点$j$的距离。$Q$: 每辆车最大行驶距离（Travelcon）。
$x_{ij}$: 若车辆从节点$i$到节点$j$，则$x_{ij}=1$，否则为0。路线序列$S$表示车辆访问顺序，包含多次返回配送中心的标记。

目标函数：
$$\min \sum _{i\in V}\sum _{j\in V}{d}_{ij}{x}_{ij}$$
约束条件：
每个客户点被访问一次。每辆车的路线从配送中心开始并结束，且总距离不超过$Q$。路线必须连续，满足距离约束。代码中，路线$S1$表示所有车辆的合并序列，包含配送中心（1）作为分隔符，计算总距离为各段路线的距离和。

数据加载与距离计算

代码从Excel文件读取数据，City为需求点和配送中心的坐标，Demand未使用，Travelcon为距离限制。Distance矩阵通过CalculateDist(City)计算，可能使用欧几里得距离公式：
$$d_{ij}=\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}$$
其中$(x_i,y_i)$为节点$i$的坐标。

初始解生成

初始解通过随机生成TSP路线（TSProute = [0, randperm(CityNum)] + 1）开始，转换为VRP路线（S1）：
（1）遍历TSProute的客户顺序，检查当前行驶距离加上从最后一个节点到新客户再返回配送中心的距离是否超过Travelcon。
（2）若超过，则在S1中插入1（返回配送中心），并从配送中心开始新车辆；否则直接添加客户。
示例：若S1 = [1,2,3,1,4,1]，表示车辆1路线为1→2→3→1，车辆2路线为1→4→1。

模拟退火过程

模拟退火从高温度$T0=1000$开始，逐步降温（$T0=q\times T0，q=0.9$），直至$Tend=1e-3$。每个温度下：
进行$L=200$次迭代，生成新解S2（NewSolution函数，未详细提供，可能通过交换或插入操作扰动S1）。使用Metropolis准则决定是否接受S2：若S2总距离小于S1，则接受。否则以概率$e^{-\Delta E/T}$接受，其中$\Delta E$为S2与S1的总距离差，T为当前温度。记录每次迭代的最优解，更新历史最佳解。

表1：参数设置与意义

运行结果
历时 1.016268 秒。
Total Distance = 935.5531 km 
最优路径:
0 -> 18 -> 1 -> 20 -> 17 -> 12 -> 8 -> 0 -> 2 -> 3 -> 4 -> 5 -> 6 -> 7 -> 16 -> 15 -> 11 -> 10 -> 9 -> 0 -> 13 -> 14 -> 19 -> 0
-------------------------------------------------------------
Route of Vehichle No.1: 0 -> 18 -> 1 -> 20 -> 17 -> 12 -> 8 -> 0  
Distance traveled: 329.27 km;  
-------------------------------------------------------------
Route of Vehichle No.2: 0 -> 2 -> 3 -> 4 -> 5 -> 6 -> 7 -> 16 -> 15 -> 11 -> 10 -> 9 -> 0  
Distance traveled: 378.23 km;  
-------------------------------------------------------------
Route of Vehichle No.3: 0 -> 13 -> 14 -> 19 -> 0  
Distance traveled: 228.05 km;  
-------------------------------------------------------------

%%部分代码
clc %清空命令行窗口
clear %从当前工作区中删除所有变量，并将它们从系统内存中释放
close all %删除其句柄未隐藏的所有图窗
tic % 保存当前时间
%% 模拟退火算法求解DVRP
%输入：
%City           需求点经纬度
%Distance       距离矩阵
%Travelcon      行程约束
%T0             初始温度
%Tend           终止温度
%L              各温度下的迭代次数（链长）
%q              降温速率

%输出：
%bestroute      最短路径
%mindisever     路径长度

%% 加载数据-
City=xlsread('data','City');%需求点经纬度，用于画实际路径的XY坐标
Demand=xlsread('data','Demand');%需求
Travelcon=xlsread('data','Travelcon');%距离限制
Distance=CalculateDist(City);

%% 初始化问题参数
CityNum=size(City,1)-1;    %需求点个数

%% 初始化算法参数
T0=1000;        %初始温度
Tend=1e-3;      %终止温度
L=200;          %各温度下的迭代次数（链长）
q=0.9;          %降温速率
count=0;        %迭代计数

%% 初始解（必须为满足各约束的可行解）
TSProute=[0,randperm(CityNum)]+1; %随机生成一条不包括尾位的TSP路线
S1=ones(1,CityNum*2+1); %初始化VRP路径S1，为其分配内存
    
DisTraveled=0; % 汽车已经行驶距离
k=1;
for j=2:CityNum+1
	k=k+1;%对VRP路径下一点进行赋值
	%若VRProute(k-1)到TSProute(j)超行程约束，则须返回配送中心再行驶至TSProute(j)
	%若VRProute(k-1)到TSProute(j)未超行程约束，但VRProute(k-1)到TSProute(j)到1超行程约束，则须返回配送中心再行驶至TSProute(j)
	if DisTraveled+Distance(S1(k-1),TSProute(j))+Distance(TSProute(j),1) > Travelcon  %这一点配送完成后,再回到配送中心不超行程约束
        S1(k)=1; %经过配送中心
        DisTraveled=Distance(1,TSProute(j)); %距离初始化为配送中心到此点距离
        % 新一辆车再去下一个需求点
        k=k+1;
        S1(k)=TSProute(j);  %TSP路线中此点添加到VRP路线
	else %
        DisTraveled=DisTraveled+Distance(TSProute(j-1),TSProute(j));
        S1(k)=TSProute(j); %TSP路线中此点添加到VRP路线
	end
end

%% 计算迭代的次数Time，即求解 T0 * q^x = Tend
Time=ceil(double(log(Tend/T0)/log(q)));

bestind=zeros(1,CityNum*2+1);      %每代的最优路线矩阵初始化
BestObjByIter=zeros(Time,1);       %每代目标值矩阵初始化

%% 迭代
while T0 > Tend
    count = count+1;     %更新迭代次数
    Population = zeros(L,CityNum*2+1); %为此温度下迭代个体矩阵分配内存
    ObjByIter = zeros(L,1); %为此温度下迭代个体的目标函数值矩阵分配内存
    for k = 1:L
        %% 产生新解
        S2 = NewSolution(S1);
        %% Metropolis法则判断是否接受新解
        [S1,ttlDis] = Metropolis(S1,S2,Distance,Travelcon,T0);  % Metropolis 抽样算法
        ObjByIter(k) = ttlDis;    %此温度下每迭代一次就存储一次目标函数值
        Population(k,:) = S1;          %此温度下每迭代一次就存储一次此代最优个体
    end
    
    %% 记录每次迭代过程的最优路线
    [d0,index] = min(ObjByIter); %取出此温度下所有迭代中最优的一次
    if count == 1 || d0 < BestObjByIter(count-1) %若为第一次迭代或上次迭代比这次更满意
        BestObjByIter(count) = d0;            %如果当前温度下最优路程小于上一路程则记录当前路程
        bestind = Population(index,:);  %记录当前温度的最优路线
    else
        BestObjByIter(count) = BestObjByIter(count-1);  %如果当前温度下最优路程大于上一路程则记录上一路程的目标函数值
    end
    
    T0 = q * T0; %降温
     %% 显示此代信息
    fprintf('Iteration = %d, Min Distance = %.2f km  \n',count,BestObjByIter(count)) %输出当前迭代信息
end

%% 找出历史最短距离和对应路径
mindisever = BestObjByIter(count); %找出历史最优目标函数值
bestroute=bestind; % 取最优个体

%删去路径中多余1
for i=1:length(bestroute)-1
    if bestroute(i)==bestroute(i+1)
        bestroute(i)=0;  %相邻位都为1时前一个置零
    end
end
bestroute(bestroute==0)=[];  %删去多余零元素

bestroute=bestroute-1;  % 编码各减1，与文中的编码一致

%% 计算结果数据输出到命令行
disp('-------------------------------------------------------------')
toc %显示运行时间
fprintf('Total Distance = %s km \n',num2str(mindisever))
TextOutput(Distance,bestroute,Travelcon)  %显示最优路径
disp('-------------------------------------------------------------')

%% 优化过程迭代图
figure
plot(BestObjByIter,'LineWidth',2)
xlim([1 count]) %设置 x 坐标轴范围
set(gca, 'LineWidth',1)
xlabel('Iterations')
ylabel('Min Distance(km)')
title('SA Optimization Process')

%% 绘制最优解的实际路线
DrawPath(bestroute,City)