1 算法介紹

說明：1.1節(jié)主要是概括和幫助理解考慮影響因素的BP神經(jīng)網(wǎng)絡(luò)算法原理，即常規(guī)的BP模型訓(xùn)練原理講解（可根據(jù)自身掌握的知識是否跳過）。1.2節(jié)開始講基于歷史值影響的BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型。

使用BP神經(jīng)網(wǎng)絡(luò)進行預(yù)測時，從考慮的輸入指標角度，主要有兩類模型：

1.1 受相關(guān)指標影響的BP神經(jīng)網(wǎng)絡(luò)算法原理

如圖一所示，使用MATLAB的newff函數(shù)訓(xùn)練BP時，可以看到大部分情況是三層的神經(jīng)網(wǎng)絡(luò)（即輸入層，隱含層，輸出層）。這里幫助理解下神經(jīng)網(wǎng)絡(luò)原理：
1）輸入層：相當(dāng)于人的五官，五官獲取外部信息，對應(yīng)神經(jīng)網(wǎng)絡(luò)模型input端口接收輸入數(shù)據(jù)的過程。
2）隱含層：對應(yīng)人的大腦，大腦對五官傳遞來的數(shù)據(jù)進行分析和思考，神經(jīng)網(wǎng)絡(luò)的隱含層hidden Layer對輸入層傳來的數(shù)據(jù)x進行映射，簡單理解為一個公式hiddenLayer_output=F(w*x+b)。其中，w、b叫做權(quán)重、閾值參數(shù)，F(xiàn)()為映射規(guī)則，也叫激活函數(shù)，hiddenLayer_output是隱含層對于傳來的數(shù)據(jù)映射的輸出值。換句話說，隱含層對于輸入的影響因素數(shù)據(jù)x進行了映射，產(chǎn)生了映射值。
3）輸出層：可以對應(yīng)為人的四肢，大腦對五官傳來的信息經(jīng)過思考（隱含層映射）之后，再控制四肢執(zhí)行動作（向外部作出響應(yīng)）。類似地，BP神經(jīng)網(wǎng)絡(luò)的輸出層對hiddenLayer_output再次進行映射，outputLayer_output=w *hiddenLayer_output+b。其中，w、b為權(quán)重、閾值參數(shù)，outputLayer_output是神經(jīng)網(wǎng)絡(luò)輸出層的輸出值（也叫仿真值、預(yù)測值）（理解為，人腦對外的執(zhí)行動作，比如嬰兒拍打桌子）。
4）梯度下降算法：通過計算outputLayer_output和神經(jīng)網(wǎng)絡(luò)模型傳入的y值之間的偏差，使用算法來相應(yīng)調(diào)整權(quán)重和閾值等參數(shù)。這個過程，可以理解為嬰兒拍打桌子，打偏了，根據(jù)偏離的距離遠近，來調(diào)整身體使得再次揮動的胳膊不斷靠近桌子，最終打中。

再舉個例子來加深理解：

圖一所示BP神經(jīng)網(wǎng)絡(luò)，具備輸入層、隱含層和輸出層。BP是如何通過這三層結(jié)構(gòu)來實現(xiàn)輸出層的輸出值outputLayer_output，不斷逼近給定的y值，從而訓(xùn)練得到一個精準的模型的呢？

從圖中串起來的端口，可以想到一個過程：坐地鐵，將圖一想象為一條地鐵線路。王某某坐地鐵回家的一天：在input起點站上車，中途經(jīng)過了很多站（hiddenLayer），然后發(fā)現(xiàn)坐過頭了（outputLayer對應(yīng)現(xiàn)在的位置），那么王某某將會根據(jù)現(xiàn)在的位置離家（目標Target）的距離（誤差Error），返回到中途的地鐵站（hiddenLayer）重新坐地鐵（誤差反向傳遞，使用梯度下降算法更新w和b），如果王某某又一次發(fā)生失誤，那么將再次進行這個調(diào)整的過程。

從在嬰兒拍打桌子和王某某坐地鐵的例子中，思考問題：BP的完整訓(xùn)練，需要先傳入數(shù)據(jù)給input，再經(jīng)過隱含層的映射，輸出層得到BP仿真值，根據(jù)仿真值與目標值的誤差，來調(diào)整參數(shù)，使得仿真值不斷逼近目標值。比如（1）嬰兒受到了外界的干擾因素（x），從而作出反應(yīng)拍桌（predict），大腦不斷的調(diào)整胳膊位置，控制四肢拍準（y、Target）。（2）王某某上車點（x），過站點（predict），不斷返回中途站來調(diào)整位置，到家（y、Target）。

在這些環(huán)節(jié)中，涉及了影響因素數(shù)據(jù)x，目標值數(shù)據(jù)y（Target）。根據(jù)x，y，使用BP算法來尋求x與y之間存在的規(guī)律，實現(xiàn)由x來映射逼近y，這就是BP神經(jīng)網(wǎng)絡(luò)算法的作用。再多說一句，上述講的過程，都是BP模型訓(xùn)練，那么最終得到的模型雖然訓(xùn)練準確，但是找到的規(guī)律（bp network）是否準確與可靠呢。于是，我們再給x1到訓(xùn)練好的bp network中，得到相應(yīng)的BP輸出值（預(yù)測值）predict1，通過作圖，計算Mse，Mape，R方等指標，來對比predict1和y1的接近程度，就可以知道模型是否預(yù)測準確。這是BP模型的測試過程，即實現(xiàn)對數(shù)據(jù)的預(yù)測，并且對比實際值檢驗預(yù)測是否準確。

圖一 3層BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖

1.2 基于歷史值影響的BP神經(jīng)網(wǎng)絡(luò)

以電力負荷預(yù)測問題為例，進行兩種模型的區(qū)分。在預(yù)測某個時間段內(nèi)的電力負荷時：

一種做法，是考慮?t?時刻的氣候因素指標，比如該時刻的空氣濕度x1，溫度x2，以及節(jié)假日x3等的影響，對?t?時刻的負荷值進行預(yù)測。這是前面1.1所說的模型。

另一種做法，是認為電力負荷值的變化，與時間相關(guān)，比如認為t-1，t-2，t-3時刻的電力負荷值與t時刻的負荷值有關(guān)系，即滿足公式y(tǒng)(t)=F(y(t-1),y(t-2),y(t-3))。采用BP神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練模型時，則輸入到神經(jīng)網(wǎng)絡(luò)的影響因素值為歷史負荷值y(t-1),y(t-2),y(t-3)，特別地，3叫做自回歸階數(shù)或者延遲。給到神經(jīng)網(wǎng)絡(luò)中的目標輸出值為y(t)。

?1.3 帝國殖民競爭算法

?

?

2 部分代碼

% 清除環(huán)境變量 close all clc; clear %% Problem Statement ProblemParams.CostFuncName = 'fitcal'; ? ?% You should state the name of your cost function here. ProblemParams.CostFuncExtraParams = []; ProblemParams.NPar = 31; ? ? ? ? ? ? ? ? ? ? ? ? ? % Number of optimization variables of your objective function. "NPar" is the dimention of the optimization problem. ProblemParams.VarMin = -1; ? ? ? ? ? ? ? ? ? ? ? ? % Lower limit of the optimization parameters. You can state the limit in two ways. 1) ? 2) ProblemParams.VarMax = 1; ? ? ? ? ? ? ? ? ? ? ? % Lower limit of the optimization parameters. You can state the limit in two ways. 1) ? 2) % Modifying the size of VarMin and VarMax to have a general form if numel(ProblemParams.VarMin)==1 ? %numel 數(shù)組元素個數(shù)計數(shù) ? ?ProblemParams.VarMin=repmat(ProblemParams.VarMin,1,ProblemParams.NPar); %復(fù)制矩陣，行數(shù)不變，仍然是roblemParams.VarMin，列數(shù)重復(fù)ProblemParams.NPar遍ProblemParams.VarMin ? ?ProblemParams.VarMax=repmat(ProblemParams.VarMax,1,ProblemParams.NPar); end ProblemParams.SearchSpaceSize = ProblemParams.VarMax - ProblemParams.VarMin; %搜索區(qū)間 %% Algorithmic Parameter Setting AlgorithmParams.NumOfCountries = 200; ? ? ? ? ? ? ? % Number of initial countries. AlgorithmParams.NumOfInitialImperialists = 10; ? ? ?% Number of Initial Imperialists. AlgorithmParams.NumOfAllColonies = AlgorithmParams.NumOfCountries - AlgorithmParams.NumOfInitialImperialists; AlgorithmParams.NumOfDecades = 100; ? ? ?%迭代次數(shù) AlgorithmParams.RevolutionRate = 0.3; ? ? ? ? ? ? ? % Revolution is the process in which the socio-political characteristics of a country change suddenly. AlgorithmParams.AssimilationCoefficient = 2; ? ? ? ?% In the original paper assimilation coefficient is shown by "beta". ?每次趨近的系數(shù) AlgorithmParams.AssimilationAngleCoefficient = .5; ?% In the original paper assimilation angle coefficient is shown by "gama". ?夾角度數(shù) AlgorithmParams.Zeta = 1; ? ? ? ? ? ? ? ? ? ? ? ?% Total Cost of Empire = Cost of Imperialist + Zeta * mean(Cost of All Colonies); AlgorithmParams.DampRatio = 0.99; AlgorithmParams.StopIfJustOneEmpire = false; ? ? ? ? % Use "true" to stop the algorithm when just one empire is remaining. Use "false" to continue the algorithm. 停止迭代的標志 AlgorithmParams.UnitingThreshold = 0.0001; ? ? ? ? ? ?% The percent of Search Space Size, which enables the uniting process of two Empires. zarib = 1.05; ? ? ? ? ? ? ? ? ? ? ? % **** Zarib is used to prevent the weakest impire to have a probability equal to zero alpha = 0.03; ? ? ? ? ? ? ? ? ? ? ? ?% **** alpha is a number in the interval of [0 1] but alpha<<1. alpha denotes the importance of mean minimum compare to the global mimimum. %% Display Setting DisplayParams.PlotEmpires = true; ? ?% "true" to plot. "false" to cancel ploting. 殖民者參數(shù) if DisplayParams.PlotEmpires ? ?DisplayParams.EmpiresFigureHandle = figure('Name','Plot of Empires','NumberTitle','off'); ? ?DisplayParams.EmpiresAxisHandle = axes; end DisplayParams.PlotCost = true; ? ?% "true" to plot. "false" 消耗參數(shù) if DisplayParams.PlotCost ? ?DisplayParams.CostFigureHandle = figure('Name','Plot of Minimum and Mean Costs','NumberTitle','off'); ? ?DisplayParams.CostAxisHandle = axes; end ColorMatrix = [1 ? 0 ? 0 ?; 0 1 ? 0 ? ?; 0 ? 0 1 ? ?; 1 ? 1 ? 0 ?; 1 ? 0 1 ? ?; 0 1 ? 1 ? ?; 1 1 1 ? ? ? ; ? ? ? ? ? ? ? 0.5 0.5 0.5; 0 0.5 0.5 ?; 0.5 0 0.5 ?; 0.5 0.5 0 ?; 0.5 0 0 ? ?; 0 0.5 0 ? ?; 0 0 0.5 ? ? ; ? ? ? ? ? ? ? 1 ? 0.5 1 ?; 0.1*[1 1 1]; 0.2*[1 1 1]; 0.3*[1 1 1]; 0.4*[1 1 1]; 0.5*[1 1 1]; 0.6*[1 1 1]]; DisplayParams.ColorMatrix = [ColorMatrix ; sqrt(ColorMatrix)]; %sqrt 平方根，什么用？ DisplayParams.AxisMargin.Min = ProblemParams.VarMin; DisplayParams.AxisMargin.Max = ProblemParams.VarMax; %% for i=1 %% Creation of Initial Empires InitialCountries = GenerateNewCountry(AlgorithmParams.NumOfCountries , ProblemParams);%建立國家 子函數(shù)調(diào)用 一個列向量 (AlgorithmParams.NumOfCountriesx1) 值為約束范圍內(nèi)的隨機數(shù) % Calculates the cost of each country. The less the cost is, the more is the power. if isempty(ProblemParams.CostFuncExtraParams) ? ?InitialCost = feval(ProblemParams.CostFuncName,InitialCountries); ? ? else ? ?InitialCost = feval(ProblemParams.CostFuncName,InitialCountries,ProblemParams.CostFuncExtraParams); end [InitialCost,SortInd] = sort(InitialCost); ? ? ? ? ? ? ? ? ? ? ? ? ?% Sort the cost in assending order. The best countries will be in higher places 排序 每行從小到大 InitialCountries = InitialCountries(SortInd,:); ? ? ? ? ? ? ? ? ? ? % Sort the population with respect to their cost. 按照相關(guān)耗費給國家排序 %stop 調(diào)試所用 Empires = CreateInitialEmpires(InitialCountries,InitialCost,AlgorithmParams, ProblemParams);%子函數(shù)調(diào)用 得到帝國種群（殖民國加殖民地） %% Main Loop MinimumCost = repmat(nan,AlgorithmParams.NumOfDecades,1); %把nan復(fù)制了AlgorithmParams.NumOfDecades次，一個超長的列向量 MeanCost = repmat(nan,AlgorithmParams.NumOfDecades,1); if DisplayParams.PlotCost ? ?axes(DisplayParams.CostAxisHandle); ? ?if any(findall(0)==DisplayParams.CostFigureHandle) ? ? ? ?h_MinCostPlot=plot(MinimumCost,'r','LineWidth',1.5,'YDataSource','MinimumCost'); ? ? ? ?hold on; ? ? ? ?h_MeanCostPlot=plot(MeanCost,'k:','LineWidth',1.5,'YDataSource','MeanCost'); ? ? ? ?hold off; ? ? ? ?pause(0.05); ? ?end end for Decade = 1:AlgorithmParams.NumOfDecades ? ?AlgorithmParams.RevolutionRate = AlgorithmParams.DampRatio * AlgorithmParams.RevolutionRate;%進化率=0.99*原進化率 ? ?Remained = AlgorithmParams.NumOfDecades - Decade ? ?for ii = 1:numel(Empires) ? ? ? ?%% Assimilation同化; ?Movement of Colonies Toward Imperialists (Assimilation Policy) ? ? ? ?Empires(ii) = AssimilateColonies(Empires(ii),AlgorithmParams,ProblemParams);%子函數(shù)調(diào)用 ? ? ? ?%% Revolution; ?A Sudden Change in the Socio-Political Characteristics ? ?有部分用重新生成的國家替代 ? ? ? ?Empires(ii) = RevolveColonies(Empires(ii),AlgorithmParams,ProblemParams);%子函數(shù)調(diào)用 ? ? ? ? ? ? ? ?%% New Cost Evaluation ? ? ? ?if isempty(ProblemParams.CostFuncExtraParams) ? ? ? ? ? ?Empires(ii).ColoniesCost = feval(ProblemParams.CostFuncName,Empires(ii).ColoniesPosition); ? ? ? ?else ? ? ? ? ? ?Empires(ii).ColoniesCost = feval(ProblemParams.CostFuncName,Empires(ii).ColoniesPosition,ProblemParams.CostFuncExtraParams); ? ? ? ?end ? ? ? ? ? ? ? ? ? ?%% Empire Possession ?(****** Power Possession, Empire Possession) ? ? ? ?Empires(ii) = PossesEmpire(Empires(ii));%子函數(shù)調(diào)用 ? ? ? ? ? ? ? ?%% Computation of Total Cost for Empires ? ? ? ?Empires(ii).TotalCost = Empires(ii).ImperialistCost + AlgorithmParams.Zeta * mean(Empires(ii).ColoniesCost); ? ? ? ? ? ?end ? ?%% Uniting Similiar Empires ? ?Empires = UniteSimilarEmpires(Empires,AlgorithmParams,ProblemParams);%子函數(shù)調(diào)用 ? ?%% Imperialistic Competition ? ?Empires = ImperialisticCompetition(Empires);%子函數(shù)調(diào)用 ? ? ? ?if numel(Empires) == 1 && AlgorithmParams.StopIfJustOneEmpire%如果只剩下一個殖民者 ? ? ? ?break ? ?end ? ?%% Displaying the Results ? ?DisplayEmpires(Empires,AlgorithmParams,ProblemParams,DisplayParams);%子函數(shù)調(diào)用 ? ? ? ?ImerialistCosts = [Empires.ImperialistCost]; ? [MinimumCost(Decade),index] = min(ImerialistCosts);% 調(diào)試 MinimumCost(Decade) = min(ImerialistCosts); ? ?MeanCost(Decade) = mean(ImerialistCosts); ? ? %% 導(dǎo)出 if Decade == AlgorithmParams.NumOfDecades ? positions(i,:)=[Empires(index).ImperialistPosition ]; ? ? results(i,:)=[MinimumCost(Decade)]; ? ? end ? ?if DisplayParams.PlotCost ? ? ? ?refreshdata(h_MinCostPlot); ? ? % ? refreshdata(h_MeanCostPlot); ? ? ? ?drawnow; grid on;hold on; ? ? ? ?xlabel('迭代代數(shù)/次');ylabel('目標函數(shù)值'); ? ? ? ?pause(0.01); ? ?end % ? ? if DisplayParams.PlotCost % ? ? ? ? refreshdata(h_MinCostPlot);%刷新圖片的數(shù)據(jù) % ? ? ? ? refreshdata(h_MeanCostPlot); % ? ? ? ? drawnow;grid on;hold on; % ? ? ? ? xlabel('迭代代數(shù)/次');ylabel('目標函數(shù)值'); % ? ? ? ? pause(0.01); % ? ? end ? ? end % End of Algorithm %% end MinimumCost(end) save Cost MinimumCost save pos positions

3 仿真結(jié)果

4 參考文獻

[1]張鐿議,焦健,汪可,鄭含博,房加珂,周浩.基于帝國殖民競爭算法優(yōu)化支持向量機的電力變壓器故障診斷模型[J].電力自動化設(shè)備,2018,38(01):99-104.

?