數學建模方法-遺傳演算法(實戰篇part 2)

一、引言

　　在上一個篇章中，我們用遺傳演算法來計算一個一元函數的最大值，但是，有人會講，這樣是不是有些大材小用了，明明我可以用更少的代碼來實現求最大值的功能。確實，將遺傳演算法用在那裡確實大材小用了，但是，博主的目的並不是為求最大值，而是為了給大家展示，遺傳演算法是一種可行的演算法，並且博主編寫的matlab代碼是可行的。現在，我們就要利用上一篇寫的代碼作為基本代碼，來解決本篇的問題。放心，本篇的問題可不是用普通的方法就能就解決的。遺傳演算法的春天終於來了。讓我們開始吧。。

二、問題描述和建模

　　負載平衡調度問題：假設有N個任務，需要負載平衡器分配給M個伺服器節點去處理。每個任務的任務長度、每台伺服器節點（下面簡稱“節點”）的處理速度已知，請給出一種任務分配方式，使得所有任務的總處理時間最短。

　　好，那麼現在題目有了，但是我們要想辦法對題目進行數學建模。顯然，我們需要2個矩陣，一個存放任務（含任務數量、任務長度），一個存放伺服器節點（含節點數量、處理速度）。

（1）任務長度矩陣（簡稱：任務矩陣）

　　我們將所有任務的任務長度用矩陣tasks表示，其中tasks(i)中的i表示任務的編號，tasks(i)表示任務i的任務長度。比如

　　tasks = [2, 4, 6, 8]                

　　表明，任務1的長度是2；任務2的長度是4，以此類推。

（2）節點處理速度矩陣（簡稱：節點矩陣）

　　我們將所有伺服器節點的處理速度用nodes表示，其中nodes(i)中的i表示節點的編號，nodes(i)表示節點j的處理速度。比如

　　nodes = [2, 1]

　　表明，節點1的速度是2；節點2的速度是1，以此類推。

（3）任務處理時間矩陣

　　當任務矩陣task和節點矩陣nodes確定下來後，那麼所有任務分配給所有節點的任務處理時間就可以確定了，我們用矩陣timeMat來存所有節點的任務處理時間，那麼timeMat(i, j)表示任務i分配給節點j處理所需的時間，它通過以下公式計算：

　　timeMat(i, j) = task(i) / nodes(j);

　　有了上面三個矩陣後，我們就算是成功將題目用數學的方式描述出來了。

　　我們用matlab代碼建立一個createData.m檔案，用來對題目進行建模並且儲存資料。如下：

%% I. 清空自變數clear allclc%% II. 題目建模% 參數定義taskNum = 100;% 任務數量taskLengthRange = [10, 100];% 任務長度範圍nodeNum = 10;% 節點數量nodeSpeedRange = [10, 100];% 節點長度範圍tasks = initRandomMat(taskNum, taskLengthRange);% 構建任務矩陣（含任務數量和任務長度）nodes = initRandomMat(nodeNum, nodeSpeedRange);% 構建節點矩陣（含節點數量和節點處理速度）timeMat = zeros(taskNum, nodeNum);% 構建處理時間矩陣for i = 1: taskNum    for j = 1: nodeNum        timeMat(i, j) = tasks(1, i) / nodes(1, j);    endend    %% III. 資料儲存save('data', 'tasks', 'nodes', 'timeMat');

　　其中initRandomMat()函數代碼如下：

% [name] --initRandomMat(初始化隨機矩陣)% [function]--% [input]--length (數組的長度)%--range (數組元素的取值範圍range(1) to range(2))% [output]--randomMatfunction randomMat = initRandomMat(length, range)    randomMat = zeros(1, length);    for i = 1: length        randomMat(1, i) = round(rand() * (range(2) - range(1)) + range(1));    end

　　三、遺傳演算法解決問題

　　還是老樣子，我們先回顧下遺傳演算法的步驟，如下：

　　1.  種群初始化
　　2. 計算每個種群的適應度值
　　3. 選擇(Selection)
　　4. 交叉(Crossover)
　　5. 變異(Mutation)
　　6. 重複2-5步直至達到進化次數

　　正如我們上一篇所說，為了讓我們的代碼有更好的包容性，我們需要先寫幾個函數包，到時再main.m裡面引用，就很方便了。

　　（1）種群初始化函數popInit()，能根據提供的種群大小和染色體長度產生初始的種群。代碼如下：

% [name]    --      popInit(種群初始化函數)% [function]--      構建種群矩陣，其中行數為種群個數，列數為染色體長度（即基因的個數），並隨機分配染色體的樣式% [input]   --      1. popSize（種群大小/個數)%           --      2. cLength（染色體長度）% [output]  --      popMat（種群矩陣） function popMat = popInit(popSize, cLength)    popMat = zeros(popSize, cLength);           % 預分配記憶體    for i = 1: popSize        for j = 1: cLength            popMat(i, j) = round(rand);         % rand產生(0,1)之間的隨機數，round()是四捨五入函數        end    end clear i;clear j;

　　（2）計算每個種群的適應度值函數getfitnessValue()。代碼如下：

% [name] --getfitnessValue(計算種群個體適應度值)% [function]--根據不同的題目要求，設計適應度方程計算的演算法%--在實際的負載平衡調度中，各個節點的任務處理是並行計算的，所以，所有任務的完成時間應該是所有節點任務完成時間的最大值，並非所有任務完成時間的總和。% [input]--timeMat (時間處理矩陣)-- timeMat(i, j): 第i個任務分配到第j個節點，需要的處理時間%--popMat (種群矩陣)-- popMat(i, j): 第i個染色體上, 第j個任務分配到的節點標號%--nodeNumMax (節點的最大標號)% [output]--fitValMat（每個染色體的適應度值）function fitnessValueMatrix = getfitnessValue(popMat, timeMat, nodeNumMax)[popSize, cLength] = size(popMat);% 擷取種群的數目（行）和染色體長度（列）fitnessValueMatrix = zeros(popSize, 1);% 初始化適應度矩陣% 最大的節點標號for i = 1: popSizemaxLength = eps;% 大於0的最小正數for nodeIndex = 1: nodeNumMaxsumLength = 0;                                                            % 清零for j = 1: cLengthif popMat(i, j) == nodeIndexsumLength = sumLength + timeMat(j, nodeIndex);endendif sumLength > maxLengthmaxLength = sumLength;            end endfitnessValueMatrix(i) = maxLength;endclear i;clear j;clear nodeIndex;

　　（3）選擇函數selection()，可以對種群進行選擇，具體演算法和代碼如下：

% [name] --selection(選擇操作)% [function]--採用輪盤賭的一個形式來進行選擇，增加不確定性，這樣種群就不容易趨向局部最優% [input]--1. fitnessValueMatrix (適應度值矩陣)%--2. popMat(未選擇的種群矩陣)%--3. type        --       1: 保留親代最優個體%                                           0: 不保留親代最優個體% [output]--updatedPopMat（經選擇後的種群矩陣）function updatedPopMat = selection(fitnessValueMatrix, popMat, type)[popSize, cLength] = size(popMat);    updatedPopMat = zeros(popSize, cLength);    fitnessValueMatrix = 100 ./ fitnessValueMatrix;% 剔除適應值為負的種群，使其適應值為0for i = 1: popSizeif fitnessValueMatrix(i, 1) < 0fitnessValueMatrix(i, 1) = 0;endend% 輪盤賭演算法P = fitnessValueMatrix / sum(fitnessValueMatrix);Pc = cumsum(P);for i = 1: popSizeindex = find(Pc >= rand);updatedPopMat(i, :) = popMat(index(1), :);end% 是否要保留親本適應度值最高的，若是，則type = 1,否則為0if type[~, bestIndex] = max(fitnessValueMatrix);updatedPopMat(popSize, :) = popMat(bestIndex, :);% 將親代最優染色體放到子代的最後一個個體中endclear i;clear j;

　　（4）交叉函數crossover()，可以對種群進行交叉，交叉的方式又分為單點交叉和多點交叉，根據輸入不同的參數來選擇不同的實現方式，具體演算法和代碼如下：

% [name] --crossover(交叉操作)% [function]--選定交叉點並進行互換% [input]--1. popMat (未交叉的種群矩陣)%--2. type--1: 單點交叉%--2: 多點交叉%--3. crossrate (交叉率) -- 建議值為0.6% [output]--updatedPopMat（經交叉後的種群矩陣）function updatedPopMat = crossover(popMat, type, crossrate)[popSize, cLength] = size(popMat);if type == 1% 單點交叉for i = 1: 2: popSizeif crossrate >= randcrossPosition = round(rand() * (cLength - 2) + 2);% 隨機擷取交叉點，去除0和1% 對 crossPosition及之後的二進位串進行交換for j = crossPosition: cLengthtemp = popMat(i, j);popMat(i, j) = popMat(i + 1, j);popMat(i + 1, j) = temp;endend        end        updatedPopMat = popMat;elseif type == 2% 多點交叉for i = 1: 2: popSizeif crossrate >= randcrossPosition1 = round(rand() * (cLength - 2) + 2);% 第一個交叉點crossPosition2 = round(rand() * (cLength - 2) + 2);% 第二個交叉點first = min(crossPosition1, crossPosition2);last = max(crossPosition1, crossPosition2);for j = first: lasttemp = popMat(i, j);popMat(i, j) = popMat(i + 1, j);popMat(i + 1, j) = temp;endend        end        updatedPopMat = popMat;elseh = errordlg('type的類型只能為1(單點交叉)或者2(多點交叉)', '進行交叉時發生錯誤');    end    clear i;clear j;

　　（5）變異函數mutation()，可以對種群進行變異，具體演算法和代碼如下：

% [name] --mutation(變異操作)% [function]--單點變異：隨機播放變異點，將其變為0或1% [input]--1. popMat (未交叉的種群矩陣)%--2. mutateRate (交叉率) -- 建議值為0.1%--3. randValur (變異後的隨機數範圍0-randValue)% [output]--updatedPopMat（經交叉後的種群矩陣）function updatedPopMat = mutation(popMat, mutateRate, randValue)[popSize, cLength] = size(popMat);for i = 1: popSizeif mutateRate >= randmutatePosition = ceil(rand() * cLength);% 隨機擷取交叉點，去除0% 對mutatePosition點進行變異popMat(i, mutatePosition) = ceil(rand() * randValue);endendupdatedPopMat = popMat;clear i;clear j;

　　（6）畫圖函數plotGraph：用於直觀的顯示在進化的過程中，種群平均總處理時間的變化。代碼如下：

% [name] --plotGraph(畫圖)% [function]--直觀的展示在進化過程中，平均適應度值的趨勢變化% [input]--1. generationTime(進化次數)%2. fitnessAverageValueMatrix(平均適應度值矩陣)% [output]--nonefunction plotGraph(fitnessAverageValueMatrix, generationTime)x = 1: 1: generationTime;y = fitnessAverageValueMatrix;plot(x, y);xlabel('迭代次數');ylabel('平均總處理時間');title('種群平均總處理時間的進化曲線');

　　好了，至此，我們就完成瞭解決本題目需要的函數塊。接下來，我們只需要編寫主函數main.m，針對本章的題目，對其進行求解即可。代碼如下：

%% I. 清空變數clear allclc%% II. 匯入資料和參數初始化%匯入參數資料tasks = load('data', 'tasks');                                                          % 匯入任務矩陣（含任務數量和任務長度）tasks = tasks.tasks;                                                                    % 匯入的是struct格式，因此要再多此一步nodes = load('data', 'nodes');                                                          % 匯入節點矩陣（含節點數量和節點處理速度）nodes = nodes.nodes;timeMat = load('data', 'timeMat');                                                      % 匯入處理時間矩陣timeMat = timeMat.timeMat;taskNum = length(tasks);                                                                % 任務數量nodeNum = length(nodes);                                                                % 節點數量iteratorNum = 1000;                                                                     % 迭代次數popSize = 300;                                                                          % 種群大小（染色體數量）crossRate = 0.6;                                                                        % 交叉機率mutateRate = 0.01;                                                                      % 變異機率 fitnessAverageValueMatrix = zeros(iteratorNum, 1);                                      % 每代平均適應度值fitnessBestValueMatrix = zeros(iteratorNum, 1);                                         % 每代最優適應度值bestIndividual = zeros(iteratorNum, taskNum);                                           % 每代最佳自變數%% III. 初始化種群popMat = popInit(popSize, taskNum, nodeNum);%% IV. 迭代繁衍擷取更好的個體(選擇、交叉、變異)for currentNum = 1: iteratorNum% 求單個染色體的總處理時間和適應度矩陣sumTimeMat = getfitnessValue(popMat, timeMat, nodeNum);% 染色體的總處理時間    % 記錄當前最好的資料    if currentNum == 1        [~, bestIndex] = min(sumTimeMat);        fitnessBestValueMatrix(currentNum) = min(sumTimeMat);        fitnessAverageValueMatrix(currentNum) = mean(sumTimeMat);        bestIndividual(currentNum, :)= popMat(bestIndex, :);% 最佳適應度值所對應的染色體（分配方案）    else        [~, bestIndex] = min(sumTimeMat);        fitnessBestValueMatrix(currentNum) = min(fitnessBestValueMatrix(currentNum - 1), min(sumTimeMat));        fitnessAverageValueMatrix(currentNum) = mean(sumTimeMat);        if fitnessBestValueMatrix(currentNum) == fitnessBestValueMatrix(currentNum - 1)            bestIndividual(currentNum, :) = bestIndividual(currentNum - 1, :);                      elseif fitnessBestValueMatrix(currentNum) == min(sumTimeMat)            bestIndividual(currentNum, :) = popMat(bestIndex, :);                                   else            h = errordlg('請檢查代碼', '意外錯誤');        end     end        if currentNum ~= iteratorNum        % 選擇        popMat = selection(sumTimeMat, popMat, 1);% 保留親代最佳個體        % 交叉        popMat = crossover(popMat, 1, crossRate);% 單點交叉        % 變異        popMat = mutation(popMat, mutateRate, nodeNum);    endend%% V. 畫圖並展示結果 % 畫圖plotGraph(fitnessAverageValueMatrix, iteratorNum);% 展示資料 disp 最短處理時間fitnessBestValueMatrix(currentNum, 1) plan = bestIndividual(currentNum, :);

　　運行上述程式，可以得到：

ans =             9.5278

　　可以得知，經過我們的遺傳演算法，可以設計出一套方案，讓總處理時間大大減小（不能保證是最小，但接近最小）。下面大家可以看看，隨著迭代次數的增加，遺傳演算法對方案的最佳化進而平均總處理時間的變化。

　　可以看到，隨著迭代次數的增加，平均總處理時間從30多降到最後的9點多。