問題描述

某市引進一架專業(yè)大型無人機用于緊急狀態(tài)下的藥品投遞，每個站點只能投放一次，可選擇指派任意站點的無人機起飛出發(fā)完成投遞任務，但必須在配送完畢后返回原來的站點。站點地理位置坐標（單位為公理）如下圖所示。每個站點及容納的病人數(shù)量見附件.mat數(shù)據(jù)，現(xiàn)要求通過數(shù)學建模，在保證速度和優(yōu)先救治病人數(shù)量多的站點前提下，提供藥品緊急配送策略，具體問題如下：

請制訂無人機的飛行路線，使盡可能多的病人盡早得到救治。

模擬退火算法

模擬退火算法是一種基于隨機搜索的優(yōu)化算法，其靈感來源于鑄煉工具時，先將固體加熱至高溫后緩慢冷卻，就能使固體獲得更好的韌性（可彎折），其原理是退火降溫過程中固體內部的晶粒能夠以較低的能量重新排列，模擬退火算法就是模擬加熱后的退火過程。它的目標是在解空間中尋找全局最優(yōu)解或接近最優(yōu)解的解（不一定能找到全局最優(yōu)，可能是比較接近全局最優(yōu)的局部最優(yōu)）。

在模擬退火算法中，首先隨機生成一個初始解，然后通過一系列的狀態(tài)轉移來搜索更優(yōu)的解。狀態(tài)轉移的過程（即馬爾可夫鏈）中，算法會引入一個控制參數(shù)（溫度），使得算法在搜索過程中能夠以一定的概率接受劣解，從而避免陷入局部最優(yōu)解。隨著搜索的進行，溫度會逐漸降低，接受劣解的概率也會逐漸減小，最終收斂到全局最優(yōu)解或接近最優(yōu)解。

Metropolis準則

模擬退火的主要思想是:在搜索區(qū)間隨機游走（即隨機選擇點)，再利用Metropolis抽樣準則，使隨機游走逐漸收斂于局部最優(yōu)解。而溫度是Metropolis算法中的一個重要控制參數(shù)，可以認為這個參數(shù)的大小控制了隨機過程向局部或全局最優(yōu)解移動的快慢。

Metropolis準則是用于決定在Metropolis算法中是否接受狀態(tài)轉移的準則。根據(jù)Metropolis準則，狀態(tài)轉移被接受的概率為1或 $exp(?\Delta E/T))$，其中ΔE表示當前狀態(tài)與下一個狀態(tài)的能量差，T表示當前的溫度。

具體來說，
如果ΔE小于0，即下一個狀態(tài)的能量比當前狀態(tài)更低，那么狀態(tài)轉移被接受（被接受的概率為1）；
如果ΔE大于0，即下一個狀態(tài)的能量比當前狀態(tài)更高，那么狀態(tài)轉移被接受的概率為$exp(?\Delta E/T)$（即有一定概率放棄狀態(tài)轉移）.

當溫度較高時，接受劣解的概率較高，有助于跳出局部最優(yōu)解；當溫度較低時，接受劣解的概率較低，有助于收斂到全局最優(yōu)解。

Metropolis準則的核心思想是通過接受概率來探索狀態(tài)空間，從而達到對系統(tǒng)的采樣和模擬。這種準則在Metropolis算法中起到了重要的作用，使得算法能夠在搜索過程中兼顧探索和利用，從而更好地找到系統(tǒng)的平衡分布。

算法流程圖：

Demo1：只考慮累計距離，通過模擬退火算法求解最短路徑

先通過Demo1，感受一下模擬退火算法處理旅行商問題的優(yōu)化過程。

matlab代碼：

clc,clear,close all;
%% ------------------------------------------------------------------------
%加載數(shù)據(jù)
load data_all.mat
data = zeros(25,2);
data(:,1) = data_all(:,2);
data(:,2) = data_all(:,3);
%% ------------------------------------------------------------------------
%模擬退火參數(shù)設置
n = size(data,1);  %站點數(shù)量
T = 100*n;         %初始溫度
L = 100;           %馬爾可夫鏈長度
K = 0.99;          %衰減參數(shù)
%% ------------------------------------------------------------------------
%站點坐標結構體
cp = struct([]);
for i = 1:n
    cp(i).x = data(i,1);
    cp(i).y = data(i,2);
end
le = 1;            %迭代次數(shù)
len(le) = funcp(cp,n);
figure(1);
while T > 0.001    %停止迭代溫度
    %多次迭代擾動，溫度降低前多次實驗
    for i = 1:L
        %計算原路線距離
        len1 = funcp(cp,n);
        %產生隨機擾動
        %隨機置換兩個不同站點的坐標
        p1 = floor(1+n*rand());
        p2 = floor(1+n*rand());
        while p1 == p2
            p1 = floor(1+n*rand());
            p2 = floor(1+n*rand());
        end
        tmp_cp = cp;
        tmp = tmp_cp(p1);
        tmp_cp(p1) = tmp_cp(p2);
        tmp_cp(p2) = tmp;
        %% ----------------------------------------------------------------
        %計算新路線總距離
        len2 = funcp(tmp_cp,n);
        %% ----------------------------------------------------------------
        %新老距離的差值，相當于能量
        delta_e = len2-len1;
        %% ----------------------------------------------------------------
        %若新路線好于舊路線，用新路線代替舊路線
        if delta_e < 0
            cp = tmp_cp;
        else
            %依概率選擇是否接收新解
            if exp(-delta_e/T) > rand()
                cp = tmp_cp;
            end
        end
    end
    le = le+1;
    %% --------------------------------------------------------------------
    %計算新路線距離
    len(le) = funcp(cp,n);
    %溫度不斷下降
    T =T*K;
    for i = 1:n-1
        plot([cp(i).x, cp(i+1).x], [cp(i).y, cp(i+1).y], "o-","LineWidth", 1.2);
        hold on;
        grid on;
    end
    plot([cp(n).x, cp(1).x], [cp(n).y, cp(1).y],'ro--',"LineWidth", 1.5);
    title(['優(yōu)化最短距離：', num2str(len(le))]);
    disp(['優(yōu)化最短距離：', num2str(len(le))])
    hold off;
    pause(0.001);
end
figure(2);
plot(len, 'LineWidth', 1.1)
grid on;
xlabel('迭代次數(shù)')
ylabel('目標函數(shù)值')
title('適應度進行曲線')
disp('優(yōu)化結束')
%% ------------------------------------------------------------------------
%計算路線總長度函數(shù)
function len = funcp(cp,n)
len = 0;
for i = 1:n-1
    len = len+sqrt((cp(i).x-cp(i+1).x)^2+(cp(i).y-cp(i+1).y)^2);  %累計歐氏距離
end
len = len+sqrt((cp(n).x-cp(1).x)^2+(cp(n).y-cp(1).y)^2);
end

最優(yōu)解之一：

由于模擬退火算法的優(yōu)化過程具有隨機性，所以得到的解不一定就是全局最優(yōu)，但肯定比瞎猜一個局部最優(yōu)要好?？梢灾貜瓦\行幾次程序，觀察最優(yōu)解的變化，正常應該差距不大。

只考慮距離的情況不用考慮路線順序，因為正反走都是一樣的距離。

適應度進行曲線：

可見適應度隨著退火過程逐漸趨于穩(wěn)定的最優(yōu)解，算法具有收斂性。

Demo2：考慮每個站點的病人數(shù)量，優(yōu)先給病人數(shù)量多的站點配送，同時兼顧累計最短距離

新建目標值變量tar，tar包含了累計距離變量len的大小，且加入了病人數(shù)量這個變量var，并乘以權重10，用1000-var是為了反向考慮優(yōu)化過程中需要使得目標值tar越來越?。?/p>

tar = tar+sqrt((cp(i).x-cp(i+1).x)^2+(cp(i).y-cp(i+1).y)^2)+10*(1000-var(i)); 

matlab代碼：

clc,clear,close all;
%% ------------------------------------------------------------------------
%加載數(shù)據(jù)
load data_all.mat
data = zeros(25,2);
data(:,1) = data_all(:,2);
data(:,2) = data_all(:,3);
var = data_all(:,4);
%% ------------------------------------------------------------------------
%模擬退火參數(shù)設置
n = size(data,1);  %站點數(shù)量
T = 100*n;         %初始溫度
L = 100;           %馬爾可夫鏈長度
K = 0.99;          %衰減參數(shù)
%% ------------------------------------------------------------------------
%站點坐標結構體
cp = struct([]);
for i = 1:n
    cp(i).x = data(i,1);
    cp(i).y = data(i,2);
end
le = 1;            %迭代次數(shù)
[tar(le), len(le)] = funcp(cp,n);
figure(1);
while T > 0.001    %停止迭代溫度
    %多次迭代擾動，溫度降低前多次實驗
    for i = 1:L
        %計算原路線距離
        [tar1, ~] = funcp(cp,n);
        %產生隨機擾動
        %隨機置換兩個不同站點的坐標
        p1 = floor(1+n*rand());
        p2 = floor(1+n*rand());
        while p1 == p2
            p1 = floor(1+n*rand());
            p2 = floor(1+n*rand());
        end
        tmp_cp = cp;
        tmp = tmp_cp(p1);
        tmp_cp(p1) = tmp_cp(p2);
        tmp_cp(p2) = tmp;
        %% ----------------------------------------------------------------
        %計算新路線總距離
        [tar2, ~] = funcp(tmp_cp,n);
        %% ----------------------------------------------------------------
        %新老距離的差值，相當于能量
        delta_e = tar2-tar1;
        %% ----------------------------------------------------------------
        %若新路線好于舊路線，用新路線代替舊路線
        if delta_e < 0
            cp = tmp_cp;
        else
            %依概率選擇是否接收新解
            if exp(-delta_e/T) > rand()
                cp = tmp_cp;
            end
        end
    end
    le = le+1;
    %% --------------------------------------------------------------------
    %計算新路線距離
    [tar(le), len(le)] = funcp(cp,n);
    %溫度不斷下降
    T =T*K;
    for i = 1:n-1
        plot([cp(i).x, cp(i+1).x], [cp(i).y, cp(i+1).y], "o-","LineWidth", 1.2);
        hold on;
        grid on;
    end
    plot([cp(n).x, cp(1).x], [cp(n).y, cp(1).y],'r-.',"LineWidth", 1.5);
    for k = 1:n
        text(data(k,1)+0.5,data(k,2)+0.75,num2str(var(k)),'fontsize',10);
    end
    title(['優(yōu)化最短距離：', num2str(len(le))]);
    disp(['優(yōu)化最短距離：', num2str(len(le))])
    disp(['起點坐標：', num2str([cp(1).x, cp(1).y])])
    disp(['最后一個站點坐標：', num2str([cp(n).x, cp(n).y])])
    hold off;
    pause(0.001);
end
figure(2);
plot(tar, 'LineWidth', 1.1)
grid on;
xlabel('迭代次數(shù)')
ylabel('目標函數(shù)值')
title('適應度進行曲線')
disp('優(yōu)化結束')
%% ------------------------------------------------------------------------
%計算目標函數(shù)和路線總長度函數(shù)
function [tar, len] = funcp(cp,n)
tar = 0;
len = 0;
for i = 1:n-1
    tar = tar+sqrt((cp(i).x-cp(i+1).x)^2+(cp(i).y-cp(i+1).y)^2)+10*(1000-var(i)); 
    len = len+sqrt((cp(i).x-cp(i+1).x)^2+(cp(i).y-cp(i+1).y)^2);
end
tar = tar+sqrt((cp(n).x-cp(1).x)^2+(cp(n).y-cp(1).y)^2)+10*(1000-var(1));
len = len+sqrt((cp(n).x-cp(1).x)^2+(cp(n).y-cp(1).y)^2);
end

最優(yōu)解之一：

查看路線的起點和終點（最后一個配送站點）：

>> disp(['起點坐標：', num2str([cp(1).x, cp(1).y])])
起點坐標：18  28
>> disp(['最后一個站點坐標：', num2str([cp(n).x, cp(n).y])])
最后一個站點坐標：21  30

可在data_all.mat的第一列數(shù)據(jù)中看到各坐標點對應的節(jié)點標號。

適應度進行曲線：

最優(yōu)解搜索過程：

北太天元版本：Baltamatica 北太天元 —— 基于模擬退火算法的旅行商問題文章來源地址http://www.zghlxwxcb.cn/news/detail-517474.html

到了這里，關于Matlab【旅行商問題】—— 基于模擬退火算法的無人機藥品配送路線最優(yōu)化的文章就介紹完了。如果您還想了解更多內容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關文章，希望大家以后多多支持TOY模板網(wǎng)！