一、粒子群算法及LSTM簡介
1 粒子群算法簡介
1.1 粒子群算法的概念**
粒子群優化算法(PSO：Particle swarm optimization) 是一種進化計算技術（evolutionary computation）。源于對鳥群捕食的行為研究。粒子群優化算法的基本思想：是通過群體中個體之間的協作和資訊共享來尋找最優解．
PSO的優勢：在于簡單容易實作并且沒有許多參數的調節。目前已被廣泛應用于函數優化、神經網絡訓練、模糊系統控制以及其他遺傳算法的應用領域。
1.2 粒子群算法分析
1.2.1基本思想
粒子群算法通過設計一種無品質的粒子來模拟鳥群中的鳥，粒子僅具有兩個屬性：速度和位置，速度代表移動的快慢，位置代表移動的方向。每個粒子在搜尋空間中單獨的搜尋最優解，并将其記為目前個體極值，并将個體極值與整個粒子群裡的其他粒子共享，找到最優的那個個體極值作為整個粒子群的目前全局最優解，粒子群中的所有粒子根據自己找到的目前個體極值和整個粒子群共享的目前全局最優解來調整自己的速度和位置。下面的動圖很形象地展示了PSO算法的過程：
1.2.2 更新規則
PSO初始化為一群随機粒子(随機解)。然後通過疊代找到最優解。在每一次的疊代中，粒子通過跟蹤兩個“極值”(pbest，gbest)來更新自己。在找到這兩個最優值後，粒子通過下面的公式來更新自己的速度和位置。
公式(1)的第一部分稱為【記憶項】，表示上次速度大小和方向的影響；公式(1)的第二部分稱為【自身認知項】，是從目前點指向粒子自身最好點的一個矢量，表示粒子的動作來源于自己經驗的部分；公式(1)的第三部分稱為【群體認知項】，是一個從目前點指向種群最好點的矢量，反映了粒子間的協同合作和知識共享。粒子就是通過自己的經驗和同伴中最好的經驗來決定下一步的運動。以上面兩個公式為基礎，形成了PSO的标準形式。
公式(2)和 公式(3)被視為标準PSO算法。
1.2.3 PSO算法的流程和僞代碼
2 LSTM簡介
2.1 LSTM控制流程
LSTM的控制流程：是在前向傳播的過程中處理流經細胞的資料，不同之處在于 LSTM 中細胞的結構和運算有所變化。
這一系列運算操作使得 LSTM具有能選擇儲存資訊或遺忘資訊的功能。咋一看這些運算操作時可能有點複雜，但沒關系下面将帶你一步步了解這些運算操作。
2.2 核心概念
LSTM 的核心概念在于細胞狀态以及“門”結構。細胞狀态相當于資訊傳輸的路徑，讓資訊能在序列連中傳遞下去。你可以将其看作網絡的“記憶”。理論上講，細胞狀态能夠将序列處理過程中的相關資訊一直傳遞下去。
是以，即使是較早時間步長的資訊也能攜帶到較後時間步長的細胞中來，這克服了短時記憶的影響。資訊的添加和移除我們通過“門”結構來實作，“門”結構在訓練過程中會去學習該儲存或遺忘哪些資訊。
2.3 Sigmoid
門結構中包含着 sigmoid 激活函數。Sigmoid 激活函數與 tanh 函數類似，不同之處在于 sigmoid 是把值壓縮到 0~1 之間而不是 -1~1 之間。這樣的設定有助于更新或忘記資訊，因為任何數乘以 0 都得 0，這部分資訊就會剔除掉。同樣的，任何數乘以 1 都得到它本身，這部分資訊就會完美地儲存下來。這樣網絡就能了解哪些資料是需要遺忘，哪些資料是需要儲存。
2.4 LSTM門結構
LSTM 有三種類型的門結構：遺忘門、輸入門和輸出門。
2.4.1 遺忘門
遺忘門的功能是決定應丢棄或保留哪些資訊。來自前一個隐藏狀态的資訊和目前輸入的資訊同時傳遞到 sigmoid 函數中去，輸出值介于 0 和 1 之間，越接近 0 意味着越應該丢棄，越接近 1 意味着越應該保留。
2.4.2 輸入門
輸入門用于更新細胞狀态。首先将前一層隐藏狀态的資訊和目前輸入的資訊傳遞到 sigmoid 函數中去。将值調整到 0~1 之間來決定要更新哪些資訊。0 表示不重要，1 表示重要。
其次還要将前一層隐藏狀态的資訊和目前輸入的資訊傳遞到 tanh 函數中去，創造一個新的侯選值向量。最後将 sigmoid 的輸出值與 tanh 的輸出值相乘，sigmoid 的輸出值将決定 tanh 的輸出值中哪些資訊是重要且需要保留下來的。
2.4.3 細胞狀态
下一步，就是計算細胞狀态。首先前一層的細胞狀态與遺忘向量逐點相乘。如果它乘以接近 0 的值，意味着在新的細胞狀态中，這些資訊是需要丢棄掉的。然後再将該值與輸入門的輸出值逐點相加，将神經網絡發現的新資訊更新到細胞狀态中去。至此，就得到了更新後的細胞狀态。
2.4.4 輸出門
輸出門用來确定下一個隐藏狀态的值，隐藏狀态包含了先前輸入的資訊。首先，我們将前一個隐藏狀态和目前輸入傳遞到 sigmoid 函數中，然後将新得到的細胞狀态傳遞給 tanh 函數。
最後将 tanh 的輸出與 sigmoid 的輸出相乘，以确定隐藏狀态應攜帶的資訊。再将隐藏狀态作為目前細胞的輸出，把新的細胞狀态和新的隐藏狀态傳遞到下一個時間步長中去。
讓我們再梳理一下。遺忘門确定前一個步長中哪些相關的資訊需要被保留；輸入門确定目前輸入中哪些資訊是重要的，需要被添加的；輸出門确定下一個隐藏狀态應該是什麼。
二、部分源代碼
%%%% 基于粒子群算法優化lstm預測單序列
clc
clear all
close all
%加載資料，重構為行向量
%加載資料，重構為行向量
data =xlsread('台風日資料2','Sheet1','B2:E481');%把你的負荷資料指派給data變量就可以了。
%%
%序列的前 90% 用于訓練，後 10% 用于測試
numTimeStepsTrain =round(0.8*size(data,1));
dataTrain = data(1:numTimeStepsTrain,:)';
dataTest = data(numTimeStepsTrain+1:end-1,:)';
numTimeStepsTest = size(data,1)-numTimeStepsTrain-1 ;%步數
%% 資料歸一化
[dataTrainStandardized, ps_input] = mapminmax(dataTrain,0,1);
[dataTestStandardized, ps_output] = mapminmax(dataTest,0,1);
%輸入LSTM的時間序列交替一個時間步
XTrain = dataTrainStandardized(2:end,:);
YTrain = dataTrainStandardized(1,:);
[dataTrainStandardized_zong, ps_input_zong] = mapminmax(data',0,1);
XTest_zong = dataTrainStandardized_zong(2:end,end);%測試集輸入
YTest_zong = dataTest(1,1)';%測試集輸出
XTest = dataTestStandardized(2:end,:);%測試集輸入
YTest = dataTest(1,:)';%測試集輸出
%%
%建立LSTM回歸網絡，指定LSTM層的隐含單元個數96*3
%序列預測，是以，輸入一維，輸出一維
numFeatures = 3;%輸入層數
numResponses = 1;%輸出層數
numHiddenUnits = 20*3;
layers = [ ...
sequenceInputLayer(numFeatures)
lstmLayer(numHiddenUnits)
fullyConnectedLayer(numResponses)
regressionLayer];
%% 初始化種群
N = 5; % 初始種群個數
d = 1; % 空間維數
ger =100; % 最大疊代次數
limit = [0.001, 0.01;]; % 設定位置參數限制(矩陣的形式可以多元)
vlimit = [-0.005, 0.005;]; % 設定速度限制
c_1 = 0.8; % 慣性權重
c_2 = 0.5; % 自我學習因子
c_3 = 0.5; % 群體學習因子
for i = 1:d
x(:,i) = limit(i, 1) + (limit(i, 2) - limit(i, 1)) * rand(N, 1);%初始種群的位置
end
v = 0.005*rand(N, d); % 初始種群的速度
xm = x; % 每個個體的曆史最佳位置
ym = zeros(1, d); % 種群的曆史最佳位置
fxm = 1000*ones(N, 1); % 每個個體的曆史最佳适應度
fym = 1000; % 種群曆史最佳适應度
%% 粒子群工作
iter = 1;
times = 1;
record = zeros(ger, 1); % 記錄器
while iter <= ger
iter
for i=1:N
%指定訓練選項，求解器設定為adam， 250 輪訓練。
%梯度門檻值設定為 1。指定初始學習率 0.005，在 125 輪訓練後通過乘以因子 0.2 來降低學習率。
options = trainingOptions('adam', ...
'MaxEpochs',250, ...
'GradientThreshold',1, ...
'InitialLearnRate',x(i,:), ...
'LearnRateSchedule','piecewise', ...
'LearnRateDropPeriod',125, ...
'LearnRateDropFactor',0.1, ...
'Verbose',0);
%訓練LSTM
net = trainNetwork(XTrain,YTrain,layers,options);%訓練網絡
net = resetState(net);
net = predictAndUpdateState(net,XTrain);%
YPred1 = [];
for mm = 1:numTimeStepsTest
[net,YPred1(:,mm)] = predictAndUpdateState(net,XTest(:,mm),'ExecutionEnvironment','cpu');%%預測
end
mint=ps_output.xmin(1);
maxt=ps_output.xmax(1);
YPred=postmnmx(YPred1,mint,maxt);%反歸一化
for mm = 1:numTimeStepsTest
if dataTest(2,mm)>25
YPred(mm)=0;
end
end
%使用先前計算的參數對預測去标準化。
%計算均方根誤差 (RMSE)。
rmse = sqrt(mean((YPred-YTest').^2));%均方差
fx(i) = rmse ; % 個體目前适應度
end
for i = 1:N
if fxm(i) > fx(i)
fxm(i) = fx(i); % 更新個體曆史最佳适應度
xm(i,:) = x(i,:); % 更新個體曆史最佳位置
YPred_best1=YPred;
end
end
if fym > min(fxm)
[fym, nmax] = min(fxm); % 更新群體曆史最佳适應度
ym = xm(nmax, :); % 更新群體曆史最佳位置
YPred_best=YPred_best1;
end
%%
clc
clear all
close all
%加載資料，重構為行向量
data =xlsread('台風日資料2','Sheet1','B2:E481');%把你的負荷資料指派給data變量就可以了。
%%
%序列的前 90% 用于訓練，後 10% 用于測試
numTimeStepsTrain =round(0.8*size(data,1));
dataTrain = data(1:numTimeStepsTrain,:)';
dataTest = data(numTimeStepsTrain+1:end-1,:)';
numTimeStepsTest = size(data,1)-numTimeStepsTrain-1 ;%步數
%% 資料歸一化
[dataTrainStandardized, ps_input] = mapminmax(dataTrain,0,1);
[dataTestStandardized, ps_output] = mapminmax(dataTest,0,1);
%輸入LSTM的時間序列交替一個時間步
XTrain = dataTrainStandardized(2:end,:);
YTrain = dataTrainStandardized(1,:);
[dataTrainStandardized_zong, ps_input_zong] = mapminmax(data',0,1);
XTest_zong = dataTrainStandardized_zong(2:end,end);%測試集輸入
YTest_zong = dataTest(1,1)';%測試集輸出
XTest = dataTestStandardized(2:end,:);%測試集輸入
YTest = dataTest(1,:)';%測試集輸出 三、運作結果
四、matlab版本及參考文獻
1 matlab版本
2014a
2 參考文獻
[1] 包子陽,餘繼周,楊杉.智能優化算法及其MATLAB執行個體（第2版）[M].電子工業出版社，2016.
[2]張岩,吳水根.MATLAB優化算法源代碼[M].清華大學出版社，2017.
[3]周品.MATLAB 神經網絡設計與應用[M].清華大學出版社，2013.
[4]陳明.MATLAB神經網絡原理與執行個體精解[M].清華大學出版社，2013.