【MATLAB】應用基因演算法(一):Peaks 函數求極值

  • 6663
  • 0
  • GA
  • 2019-11-18

每當遇到最佳化的問題,夏恩最常使用的方法就是基因演算法(Genetic Algorithm, GA)。
隨著手邊的專案愈來愈多,用法也愈來愈多元,在此簡單地分享幾個常見的用法。

目次

1. 【MATLAB】應用基因演算法(一):Peaks 函數求極值
2. 【MATLAB】應用基因演算法(二):找尋正確字串
3. 【MATLAB】應用基因演算法(三):高效率的送貨員

1. GA 調用方法

照例,先介紹一下夏恩的工作環境:

作業系統:Microsoft Windows 10 專業版
軟體版本:
MATLAB                                     Version 9.3         (R2017b)
Optimization Toolbox                  Version 8.0         (R2017b)
Global Optimization Toolbox       Version 3.4.3      (R2017b)

夏恩直接省略 GA 的理論介紹了,這些在維基百科或是各大搜尋引擎都能找到一卡車的資料。
但這裡雖然不提 GA 的原理,不表示原理不重要。GA 在實務應用上需要很好的想像力。

演算法是死的,問題是活的...噢不對,應該說,問題都是來找碴的!
若是不了解 GA 的構成法則,會大幅地限制解題的想像力。

首先介紹在 MATLAB 中的 GA 函數調用方法,如下:

x = ga(fitnessfcn,nvars)
x = ga(fitnessfcn,nvars,A,b)
x = ga(fitnessfcn,nvars,A,b,Aeq,beq)
x = ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB)
x = ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB,nonlcon)
x = ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB,nonlcon,options)
x = ga(fitnessfcn,nvars,A,b,[],[],LB,UB,nonlcon,IntCon)
x = ga(fitnessfcn,nvars,A,b,[],[],LB,UB,nonlcon,IntCon,options)

[x,fval] = ga(fitnessfcn,nvars,...)
[x,fval,exitflag] = ga(fitnessfcn,nvars,...)
[x,fval,exitflag,output] = ga(fitnessfcn,nvars,...)
[x,fval,exitflag,output,population] = ga(fitnessfcn,nvars,...)
[x,fval,exitflag,output,population,scores] = ga(fitnessfcn,nvars,...)

以夏恩為例子,最常使用的方法為:

[x, fval] = ga(fitnessfun, nvars, A, B, C, D, LB, UB, [ ], IntCon, options)
其中,限制式為:A*x < B、C*x = D。

相對應的中文翻譯為:

[最佳解, 適應值] = ga(適應函數,變數數量,A,B,C,D,下界,上界,非線性約束式,整數欄位,參數)
其中 A、B 為不等式係數;C、D為等式係數。

到這邊,看完應該還是不知道該怎麼用,對吧?
所以直接來看看第一個範例吧!

2. Peaks 函數求最小值

在 MATLAB 中,有個超好用的函數:Peaks。
這個函數看起來很漂亮,寫起來很複雜:

除了可以拿來測試立體繪圖的功能外,還能拿來算算 GA 最佳化,是個很棒的傢伙。

在指令列輸入 peaks,即可看到下圖:

在開始最佳化計算之前,首先是參數設定:

options = optimoptions( ...
             'ga', ...                                    % 最佳化算法
             'PopulationSize', 60, ...                    % 染色體數量
             'MaxGenerations', 100, ...                   % 最大繁衍代數
             'PlotFcn', {@gaplotbestf, @my_plot}, ...     % 繪圖函數
             'CrossoverFraction', 0.8, ...                % 交配率
             'Display', 'iter');                          % 結果展示方式

% 突變率設定方式
% options = optimoptions('MutationFcn', {@mutationuniform, 0.1})

針對上述的程式碼,夏恩補充幾點:

一、關於突變率的設定,並不存在 'MutationFraction' 這個參數,預設只有 'MutationFcn',若要自行設定突變率,請使用 {@mutationgaussian, scale, shrink},或是更改突變函數,使用 {@mutationuniform, rate}。

※ 2019.07.24 更新 ※
經貼心的讀者提醒,在 R2019a 中,若給定上下邊界的限制,就不能使用 mutationuniform 這個方法!
此時只能使用預設的 mutationgaussian,或自行撰寫突變函數,這部分請小心使用。

二、染色體數量、繁衍代數等,都有預設值,請參考 GA 手冊

三、展示方式有四種,分別為:off | iter | diagnose | final。預設為 final。

3. 適應函數 & 自定義繪圖函數

再來是自定適應函數:

% 自定義適應函數
function Z = FitnessFcn(x)

    X = x(1);
    Y = x(2);

    % 以下為 Peaks 函數
    Z =  3*(1-X).^2.*exp(-(X.^2) - (Y+1).^2) ... 
         - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) ... 
         - 1/3*exp(-(X+1).^2 - Y.^2);
end

這邊要注意的是適應函數的輸入值,不論有幾個變數,通通只能有一個輸入值!!!
從上面的程式中可以看到,若我們給定兩個輸入值,則在適應函數中須使用 x(1),x(2)。
若有 N 個輸入值,則使用 x(1),...,x(n) 來表示不同的輸入值。

可是我就是想要輸入其他東西啊,該怎麼辦呢?
若要輸入其他參數,必須使用 匿名函數 的技巧,這等下一章再聊,這個範例不會用到。

再接著,是繪圖函數,這個大概是本範例中最困難的地方了。
MATLAB 提供很多繪圖的函數供使用者挑選,比較常見的有:

1. @gaplotbestf               :最佳分數及平均分數
2. @gaplotscores            :每代個體的分數
3. @gaplotdistance         :每代個體間的平均距離
4. @gaplotscorediversity:繪製每一代分數的直方圖

此外,本範例所使用的是自定繪圖函數,自定繪圖函數有固定格式,如下:

function state = plotfun(options,state,flag)

三個輸入的參數分別為:

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

從 MATALB 中,可以看到的內含格式如下圖:
左邊的 state 的內容,右邊是 options 的內容。
坦白說,裡面有什麼其實不太重要,有需要再仔細看一下就好。

有了以上的規範,夏恩就依樣畫葫蘆,寫一個自定的繪圖函數,範例如下:

% 自定義繪圖函數,畫出每次迭代的結果。
function state = my_plot(options, state, flag)
    
    % 繪製所有點位的適應值  
    X = state.Population(:, 1);
    Y = state.Population(:, 2);
    Z =  3*(1-X).^2.*exp(-(X.^2) - (Y+1).^2) ... 
       - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) ... 
       - 1/3*exp(-(X+1).^2 - Y.^2);

    % 叫出 peaks 圖形
    peaks;

    % 把所有適應值,繪製到同一張圖面上
    hold on
    plot3(X, Y, Z, 'r.', 'MarkerSize', 10, 'Parent', gca)
    hold off
end

最後,就開始執行 GA 吧!

% 第一項輸入值為適應函數
% 第二項為最佳化的參數個數,在本範例中有X與Y,故輸入2
% 第三~六項為限制式
% 第七、八項為下界與上界
% 第九項為非線性約束式
% 第十項為整數欄位,在本範例中,並不限制一定為整數。
% 第十一項為剛才設定的 GA 參數。
[x, fval] = ga(@FitnessFcn, 2, [], [], [], [], [-3, -3], [3, 3], [], [], options);

以下為執行結果,我們可以清楚地看到紅色點點最後都收斂到最小值的區域。
偶爾出現在奇怪地方的紅點,是突變函數造成的,目的是為了找尋其他區域的最佳解。

好了,搞定一個,那如果要找最大值怎麼辦?
別緊張,把適應函數裡面的輸出值加個負號就行了!

像是: Z = -Z ; 就可以了。

本篇文章的完整程式碼放在附錄,請隨意使用。

下一篇是使用 GA 來找尋正確的字串,夏恩試著用比較簡單的方式來討論重複排列之應用問題。
【MATLAB】應用基因演算法(二):找尋正確字串

附錄

%% Peaks 函數求極值
% 功能:如標題
%
% 2018.03.02, Shayne.

options = optimoptions( ...
             'ga', ...                                    % 最佳化算法
             'PopulationSize', 60, ...                    % 染色體數量
             'MaxGenerations', 100, ...                   % 最大繁衍代數
             'PlotFcn', {@gaplotbestf, @my_plot}, ...     % 繪圖函數
             'CrossoverFraction', 0.8, ...                % 交配率
             'Display', 'iter');                          % 結果展示方式
         
[x, fval] = ga(@FitnessFcn, 2, [], [], [], [], [-3, -3], [3, 3], [], [], options);

% 自定義適應函數
function Z = FitnessFcn(x)

    X = x(1);
    Y = x(2);

    % 以下為 Peaks 函數
    Z =  3*(1-X).^2.*exp(-(X.^2) - (Y+1).^2) ... 
         - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) ... 
         - 1/3*exp(-(X+1).^2 - Y.^2);
end

% 自定義繪圖函數,畫出每次迭代的結果。
function state = my_plot(options, state, flag)
    
    % 繪製所有點位的適應值  
    X = state.Population(:, 1);
    Y = state.Population(:, 2);
    Z =  3*(1-X).^2.*exp(-(X.^2) - (Y+1).^2) ... 
       - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) ... 
       - 1/3*exp(-(X+1).^2 - Y.^2);

    % 叫出 peaks 圖形
    peaks;

    % 把所有適應值,繪製到同一張圖面上
    hold on
    plot3(X, Y, Z, 'r.', 'MarkerSize', 10, 'Parent', gca)
    hold off
end