【MATLAB】應用基因演算法(三):高效率的送貨員

  • 2010
  • 0
  • GA
  • 2018-12-26

最常被提起的最佳化,莫過於旅行者問題,也就是「最佳路徑安排法」。
若把旅行者問題再進一步延伸,就會進化為生產企劃排程問題。

前者是物流業的最愛,而後者則非製造業莫屬了。
身為工程師,怎能不了解其應用呢?

來吧,千篇一律的開場白,當然得先介紹一下夏恩的工作環境:

作業系統:Microsoft Windows 10 專業版
軟體版本:
MATLAB                                  Version 9.5      (R2018b)
Optimization Toolbox               Version 8.2      (R2018b)
Global Optimization Toolbox   Version 4.0      (R2018b)

在本章節會大量使用自定義函數,需要完成的函數有五個,分別為:
1. 初始化函數、2. 交配函數、3. 突變函數、4. 適應函數、5. 繪圖函數。

這麼多東西,若您不想詳看細節,可以直接拉到最後面。
夏恩把完整的測試程式放在那邊囉~

1. 座標化

想要安排最佳路徑,首先得把每個地點座標化才行。
例如大台北地區的商家位置,道路配置等,最好是可以建立一個完整的表格。

ex:
起    訖    路程        時間      其他
A      B     10 km     5 min    常有測速照相
B      A     10 km   15 min    有惱人的單行道
......

這樣可以提高演算法的擴充性,還有多目標運算的能力,計算時直接查表,節省時間。
不過現在只是要做個範例,夏恩就直接用程式隨機生成點位吧~

% 隨機給定 N 個座標,限定在 100~900 之間
N = 40;
Px = randi([100, 900], 1, N);
Py = randi([100, 900], 1, N);

% 看一下這些座標在哪兒
plot(Px, Py, 'o', 'MarkerFaceColor', 'g', 'MarkerSize', 8)
axis([0 1000 0 1000])

uiwait(msgbox('繼續執行'))
close(gcf)

2. 自定義初始化函數

什麼?連初始化函數都要自己寫?夏恩你在跟我開玩笑嗎?

沒開玩笑,旅行者的問題幾乎要從頭到尾打造專屬的基因演算法。
程式寫久了,跳脫預設值,大幅改動設定是家常便飯,這裡只是剛入門而已。

看到這邊,有沒有覺得熱血沸騰,興奮無比?!

好吧,沒有也沒關係,反正就開始吧~

如同先前提到過的自定義繪圖函數,舉凡任何自定義函數,都需要遵守 MATLAB 的規範。
也就是說,這些函數都有固定的輸入參數格式,與輸出格式。

初始化函數的寫法如下:

Population = createFcn(NVARS, FitnessFcn, options)

輸入參數分別為:

NVARS:       變數數量
FitnessFcn:適應函數
options:      參數設定

接著回到旅行者的問題,我們希望初始化的是「前往各地的先後順序」,
假設有 4 個地點,那麼我們要產生的參數是:1234, 4321, 1243,... 共 4!=24 種排法。

在 MATLAB 中,要隨機產生不重複排列的函數是:randperm,下圖為呼叫結果:

根據上面敘述,夏恩提供初始化函數的範例,如下:

% 自定義初始化函數
function Population = createFcn(NVARS, FitnessFcn, options)
    % The arguments to the function are 
    %   NVARS:      Number of variables 
    %   FITNESSFCN: Fitness function 
    %   OPTIONS:    Options structure used by the GA

    % 母群集大小
    pop = options.PopulationSize;
    
    % 隨機陣列
    Population = arrayfun(@(x) {randperm(x)}, ones(pop, 1) .* NVARS);

    % 上一行的簡單易懂寫法
    % Population= cell(pop, 1);
    % for i = 1:pop
    %    Population{i} = randperm(NVARS); 
    % end
end

3. 自定義交配函數

在 GA 中,較常見交配算法是隨機挑選兩條染色體,以某個點位交換基因片段。

在這裡千萬別有這種想法!

回顧一下剛才提到的,每一條染色體都代表了「造訪各地的順序」。
例如 A 染色體是 1 2 3 4;B 染色體是 4 3 2 1,若兩者以中間點進行單點交配,會有什麼結果?

A 染色體變成 1 2 2 1;B 染色體是 4 3 3 4 ...  ( 送貨員表示:_____  )

也就是說染色體互相獨立,不互相影響。
這樣一來,程式的寫法就很有限了,因為我們只能針對同一條染色體做變化。

大概的方法有:1. 隨機交換兩點;2. 隨機交換片段; 3. 隨機打亂 ... 之類的。

這邊夏恩使用隨機交換片段來做交配函數,其標準格式如下:

crossChild  = crossoverFcn(parents, options, NVARS, FitnessFcn, thisScore, thisPopulation)

輸入的參數分別為:

parents:             由選擇函數挑選的親代
options:             參數設定

NVARS:              變數數量
FitnessFcn:       適應函數
thisScore:          每條染色體的適應值
thisPopulation: 全部染色體

接著是交配函數的範例,交配的方法不會只有一種,請各位看官自行發揮想像力:

% 自定義交配函數
% 隨機產生 N 個隨機排列,其中 N 為母群集大小
function crossChild  = crossoverFcn(parents, options, NVARS, FitnessFcn, thisScore, thisPopulation)
    % The arguments to the function are 
    %   PARENTS:        Parents chosen by the selection function
    %   OPTIONS:        Options created from OPTIMOPTIONS
    %   NVARS:          Number of variables 
    %   FITNESSFCN:     Fitness function  
    %   THISSCORE:      Vector of scores of the current population 
    %   THISPOPULATION: Matrix of individuals in the current population

    % 由於預設為兩染色體交配,故子群集數量為母群集的一半。
    % 在此仍保持這種寫法,程式才不會出錯。
    kid = length(parents)/2;
    crossChild = cell(kid, 1);
    
    idx = 1;
    for i = 1:kid
        child = thisPopulation{ parents(idx) };
        idx = idx + 2;
        
        % 隨機給定兩點位,並左右翻轉排列順序
        p1 = ceil((length(child) -1) * rand);
        p2 = p1 + ceil((length(child) - p1- 1) * rand);
        
        child(p1:p2) = fliplr(child(p1:p2)); 
        
        % 輸出
        crossChild{i} = child;
    end
end

4.自定義突變函數

突變和交配的內容其實可以交換,反正都是對同一條染色體做事情。
不過要注意的是交配函數每次變動的染色體數量較多,突變較少。
因此,不同的算法搭配可能會使演算法的效率有相當大的差異。

其實把交配和突變函數的內容互換,夏恩也有嘗試過,但發現搜尋效率會變差。
總之,先來看看突變函數的標準格式吧,若您有空也可以自行開發高效率的寫法:

mutChild  = mutateFcn(parents, options, NVARS, FitnessFcn, state, thisScore, thisPopulation)

輸入的參數分別為:

parents:             由選擇函數挑選的親代
options:             參數設定

NVARS:              變數數量
FitnessFcn:       適應函數
state:                  目前函數相關訊息

thisScore:          每條染色體的適應值
thisPopulation: 全部染色體

在突變函數中,夏恩使用點對點交換的寫法,範例程式請參考:

% 自定義突變函數
function mutChild = mutateFcn(parents, options, NVARS, FitnessFcn, state, thisScore, thisPopulation)
    % The arguments to the function are 
    %   PARENTS:        Parents chosen by the selection function
    %   OPTIONS:        Options created from OPTIMOPTIONS
    %   NVARS:          Number of variables 
    %   FITNESSFCN:     Fitness function 
    %   STATE:          State structure used by the GA solver 
    %   THISSCORE:      Vector of scores of the current population 
    %   THISPOPULATION: Matrix of individuals in the current population
    
    % 在 ga 函數中會挑出要突變染色體,作為 parents 輸入,
    mutChild= cell(length(parents), 1);
    
    for i=1:length(parents) 
        parent = thisPopulation{ parents(i) };
        child = parent;
        
        % 隨機選擇染色體中的兩個點位,並交換基因位置
        p = randi([1, NVARS], [1, 2]);
        child(p(1)) = parent(p(2));
        child(p(2)) = parent(p(1));

        mutChild{i} = child;
    end
end

5. 自定義繪圖函數

這個函數常常出現,夏恩就不贅述了。

自定繪圖函數有固定格式,如下:

function state = plotfun(options,state,flag)

三個輸入的參數分別為:

options:最新的參數設定。
state     :最新子代的資訊。 
flag       :目前算法的狀態,例如:'init'、'iter'...等。

這裡仍然需要帶入外部參數,因此使用匿名函數的代入法。
若您還不熟悉使用方法,請參閱上一篇:找尋正確字串

繪圖函數範例程式如下:

% 自定義繪圖函數
function state = plotFcn(options, state, flag, Px, Py)
    % 找尋最佳適應值
    [score, i] = min(state.Score);
    bestgene = state.Population{i, :};
    
    % 製圖
    plot(Px(bestgene), Py(bestgene), '-o', 'MarkerFaceColor', 'g', ...
                                           'LineWidth', 2, ...
                                           'MarkerSize', 8, ...
                                           'Parent', gca);
    axis([0 1000 0 1000])
    title(sprintf('Score = %.2f', score))
end

6. 適應函數

這裡大概是本章最簡單的部分了。

適應值在旅行者問題中,就是「總路徑長」。
因此,我們要計算的就是每條染色體所提供的排列方法之路徑長度。

長度公式,還記得怎麼算吧?

D = sqrt( X^2 + Y^2 )

到這邊您可以發現,若能用查表的方法,就可省去這段計算,
當需要最佳化的點位達到某個數量時,節省下來的計算量相當可觀。

但是相對的,要建立表格,可能得付出很高的成本,這裡就見仁見智,端看各自需求。

適應函數的範例程式如下,同樣在這裡也用到匿名函數代參數的方法:

%% 適應函數
function score = FitnessFcn(x, Px, Py)

    % 資料型態轉換
    x = cell2mat(x);
    
    score = 0;
    for i = 1:length(x)
        if i==1
            % 考慮從終點回到起點的距離
            X = Px(x(i)) - Px(x(end));
            Y = Py(x(i)) - Py(x(end));
        else
            % 兩點間座標計算
            X = Px(x(i)) - Px(x(i-1));
            Y = Py(x(i)) - Py(x(i-1));
        end
        
        % 距離累加
        dist = sqrt(X.^2 + Y.^2);
        score = score + dist;
    end
end

7. 主程式,GA 參數設定

在此我們來設定一下剛才提到的匿名函數,以及 GA 的相關參數。

% 參數個數
var_Num = N;

% 繪圖函數
myPlotFcn = @(options, state, flag) plotFcn(options, state, flag, Px, Py);

% 適應函數
myFitnessFcn = @(x) FitnessFcn(x, Px, Py);

options = optimoptions( ...
            'ga', ...                                     % 最佳化算法
            'PopulationType',         'custom', ...       % 自定初始化函數
            'InitialPopulationRange', [1; var_Num], ...   % 初始化範圍
            'PopulationSize',         var_Num*5, ...      % 染色體數量
            'MaxGenerations',         var_Num*20, ...     % 最大繁衍代數
            'MaxStallGenerations',    var_Num*10, ...     % 最大停滯代數
            'CreationFcn',            @createFcn, ...     % 自定義初始化函數
            'CrossoverFcn',           @crossoverFcn, ...  % 自定義交配函數
            'MutationFcn',            @mutateFcn, ...     % 自定義突變函數
            'PlotFcn',                myPlotFcn);         % 自定義繪圖函數

8. Go~

x = ga(myFitnessFcn, var_Num, options);

以下為程式執行結果:

9. 小結

在實際使用上,旅行者問題只是圖論的一小部分,
像是先前提到的生企排程,或是一筆畫問題、七橋問題等,都是相關的變體。
如果要使用窮舉法來處理旅行者問題,那大概只要 20 個座標點就可以癱瘓一般的電腦了。
因此,熟悉 GA 算法,詳細了解內部的參數設定,會有助於我們在未來面對更困難的情境與限制條件。

用高效率的方法,來解決各種棘手的問題,才是我們工程師存在的意義。

參考資料

1. Graph theory
2. Traveling Salesman Problem
3. Genetic Algorithm Options

附錄

%% 路徑最佳化
% 功能:搜尋最短遍歷路徑
%
%  2018.11.02, Shayne.

% 隨機給定 N 個座標,限制範圍在 100~900 之間
N = 40;
Px = randi([100, 900], 1, N);
Py = randi([100, 900], 1, N);

% 看一下這些座標在哪兒
plot(Px, Py, 'o', 'MarkerFaceColor', 'g', 'MarkerSize', 8)
axis([0 1000 0 1000])

uiwait(msgbox('繼續執行'))
close(gcf)

% 參數個數
var_Num = N;

% 繪圖函數
myPlotFcn = @(options, state, flag) plotFcn(options, state, flag, Px, Py);

% 適應函數
myFitnessFcn = @(x) FitnessFcn(x, Px, Py);

options = optimoptions( ...
            'ga', ...                                     % 最佳化算法
            'PopulationType',         'custom', ...       % 自定初始化函數
            'InitialPopulationRange', [1; var_Num], ...   % 初始化範圍
            'PopulationSize',         var_Num*2, ...      % 染色體數量
            'MaxGenerations',         var_Num*20, ...     % 最大繁衍代數
            'MaxStallGenerations',    var_Num*10, ...     % 最大停滯代數
            'CreationFcn',            @createFcn, ...     % 自定義初始化函數
            'CrossoverFcn',           @crossoverFcn, ...  % 自定義交配函數
            'MutationFcn',            @mutateFcn, ...     % 自定義突變函數
            'PlotFcn',                myPlotFcn);         % 自定義繪圖函數

% 啟動演算法
x = ga(myFitnessFcn, var_Num, options);

% ========================================================================= %
%                             以下為自定義函數區塊
% ========================================================================= %

%% 適應函數
function score = FitnessFcn(x, Px, Py)

    % 資料型態轉換
    x = cell2mat(x);
    
    score = 0;
    for i = 1:length(x)
        if i==1
            % 考慮從終點回到起點的距離
            X = Px(x(i)) - Px(x(end));
            Y = Py(x(i)) - Py(x(end));
        else
            % 兩點間座標計算
            X = Px(x(i)) - Px(x(i-1));
            Y = Py(x(i)) - Py(x(i-1));
        end
        
        % 距離累加
        dist = sqrt(X.^2 + Y.^2);
        score = score + dist;
    end
end

%% 自定義繪圖函數
function state = plotFcn(options, state, flag, Px, Py)
    % 找尋最佳適應值
    [score, i] = min(state.Score);
    bestgene = state.Population{i, :};
    
    % 製圖
    plot(Px(bestgene), Py(bestgene), '-o', 'MarkerFaceColor', 'g', ...
                                           'LineWidth', 2, ...
                                           'MarkerSize', 8, ...
                                           'Parent', gca);
    axis([0 1000 0 1000])
    title(sprintf('Score = %.2f', score))
end

%% 自定義初始化函數
% 隨機產生 N 個隨機排列,其中 N 為母群集大小
function Population = createFcn(NVARS, FitnessFcn, options)
    % The arguments to the function are 
    %   NVARS:      Number of variables 
    %   FITNESSFCN: Fitness function 
    %   OPTIONS:    Options structure used by the GA

    % 母群集大小
    pop = options.PopulationSize;
    
    % 隨機陣列
    Population = arrayfun(@(x) {randperm(x)}, ones(pop, 1) .* NVARS);
    
    % 簡單易懂的寫法
    % Population= cell(pop, 1);
    % for i = 1:pop
    %    Population{i} = randperm(NVARS); 
    % end
end

%% 自定義交配函數
% 隨機產生 N 個隨機排列,其中 N 為母群集大小
function crossChild  = crossoverFcn(parents, options, NVARS, FitnessFcn, thisScore, thisPopulation)
    % The arguments to the function are 
    %   PARENTS:        Parents chosen by the selection function
    %   OPTIONS:        Options created from OPTIMOPTIONS
    %   NVARS:          Number of variables 
    %   FITNESSFCN:     Fitness function 
    %   THISSCORE:      Vector of scores of the current population 
    %   THISPOPULATION: Matrix of individuals in the current population

    % 由於預設為兩染色體交配,故子群集數量為母群集的一半。
    % 在此仍保持這種寫法,程式才不會出錯。
    kid = length(parents)/2;
    crossChild = cell(kid, 1);
    
    idx = 1;
    for i = 1:kid
        child = thisPopulation{ parents(idx) };
        idx = idx + 2;
        
        % 隨機給定兩點位,並左右翻轉排列順序
        p1 = ceil((length(child) -1) * rand);
        p2 = p1 + ceil((length(child) - p1- 1) * rand);
        
        child(p1:p2) = fliplr(child(p1:p2)); 
        
        % 輸出
        crossChild{i} = child;
    end
end

%% 自定義突變函數
function mutChild = mutateFcn(parents, options, NVARS, FitnessFcn, state, thisScore, thisPopulation)
    % The arguments to the function are 
    %   PARENTS:        Parents chosen by the selection function
    %   OPTIONS:        Options created from OPTIMOPTIONS
    %   NVARS:          Number of variables 
    %   FITNESSFCN:     Fitness function 
    %   STATE:          State structure used by the GA solver 
    %   THISSCORE:      Vector of scores of the current population 
    %   THISPOPULATION: Matrix of individuals in the current population
    
    % 在 ga 函數中會挑出要突變染色體,作為 parents 輸入,
    mutChild= cell(length(parents), 1);
    
    for i=1:length(parents) 
        parent = thisPopulation{ parents(i) };
        child = parent;
        
        % 隨機選擇染色體中的兩個點位,並交換基因位置
        p = randi([1, NVARS], [1, 2]);
        child(p(1)) = parent(p(2));
        child(p(2)) = parent(p(1));

        mutChild{i} = child;
    end
end