DE（差分進化）優化演算法MATLAB原始碼詳細中文註解

以優化SVR演算法的引數c和g為例，對DE（差分進化）演算法MATLAB原始碼進行了詳細中文註解。
完整程式和示例檔案地址：http://download.csdn.net/detail/u013337691/9671714
百度雲連結: http://pan.baidu.com/s/1dEYAHS9 密碼: 6xw5

function [bestc,bestg,test_pre]=my_DE_SVR(para,input_train,output_train,input_test,output_test)
% 引數向量 parameters [n,N_iteration,beta_min,beta_max,pCR]
% n為種群規模，N_iteration為迭代次數
% beta_min 縮放因子下界 Lower Bound of Scaling Factor
% beta_max=0.8; % 縮放因子上界 Upper Bound of Scaling Factor
% pCR 交叉概率 Crossover Probability
% 要求輸入資料為列向量（矩陣）

%% 資料歸一化
[input_train,rule1]=mapminmax(input_train');
[output_train,rule2]=mapminmax(output_train');
input_test=mapminmax('apply',input_test',rule1);
output_test=mapminmax('apply',output_test',rule2);
input_train=input_train';
output_train=output_train';
input_test=input_test';
output_test=output_test';
%% 利用差分進化（DE）演算法選擇最佳的SVR引數
nPop=para(1); % 種群規模 Population Size
MaxIt=para(2); % 最大迭代次數Maximum Number of Iterations
nVar=2; % 自變數維數，此例需要優化兩個引數c和g Number of Decision Variables
VarSize=[1,nVar]; % 決策變數矩陣大小 Decision Variables Matrix Size
beta_min=para(3); % 縮放因子下界 Lower Bound of Scaling Factor
beta_max=para(4); % 縮放因子上界 Upper Bound of Scaling Factor
pCR=para(5); %  交叉概率 Crossover Probability
lb=[0.01,0.01]; % 引數取值下界
ub=[100,100]; % 引數取值上界
%% 初始化 Initialization
empty_individual.Position=[]; % 種群初始化
empty_individual.Cost=[]; % 種群目標函式值初始化
BestSol.Cost=inf; % 最優值初始化
pop=repmat(empty_individual,nPop,1); % 將儲存種群資訊的結構體擴充套件為結構體矩陣，行數等於種群大小
for i=1:nPop % 遍歷每個個體
    pop(i).Position=init_individual(lb,ub,nVar,1); % 隨機初始化個體    
    pop(i).Cost=fobj(pop(i).Position,input_train,output_train,input_test,output_test); % 計算個體目標函式值
    if pop(i).Cost<BestSol.Cost % 如果個體目標函式值優於當前最優值
        BestSol=pop(i); % 更新最優值
    end    
end
BestCost=zeros(MaxIt,1); % 初始化迭代最優值
%% 主迴圈 DE Main Loop
for it=1:MaxIt
    for i=1:nPop % 遍歷每個個體
        x=pop(i).Position; % 提取個體位置
        % 隨機選擇三個個體以備變異使用
        A=randperm(nPop); % 個體順序重新隨機排列
        A(A==i)=[]; % 當前個體所排位置騰空（產生變異中間體時當前個體不參與）
        a=A(1);
        b=A(2);
        c=A(3);
        % 變異操作 Mutation
        beta=unifrnd(beta_min,beta_max,VarSize); % 隨機產生縮放因子
        y=pop(a).Position+beta.*(pop(b).Position-pop(c).Position); % 產生中間體
        % 防止中間體越界
        y=max(y,lb);
        y=min(y,ub);
        % 交叉操作 Crossover
        z=zeros(size(x)); % 初始化一個新個體
        j0=randi([1,numel(x)]); % 產生一個偽隨機數，即選取待交換維度編號？？？
        for j=1:numel(x) % 遍歷每個維度
            if j==j0 || rand<=pCR % 如果當前維度是待交換維度或者隨機概率小於交叉概率
                z(j)=y(j); % 新個體當前維度值等於中間體對應維度值
            else
                z(j)=x(j); % 新個體當前維度值等於當前個體對應維度值
            end
        end
        NewSol.Position=z; % 交叉操作之後得到新個體
        NewSol.Cost=fobj(NewSol.Position,input_train,output_train,input_test,output_test); % 新個體目標函式值
        if NewSol.Cost<pop(i).Cost % 如果新個體優於當前個體
            pop(i)=NewSol; % 更新當前個體
            if pop(i).Cost<BestSol.Cost % 如果當前個體（更新後的）優於最優個體
               BestSol=pop(i); % 更新最優個體
            end
        end
    end    
    % 儲存當前迭代最優個體函式值 Update Best Cost
    BestCost(it)=BestSol.Cost;  
end
bestc=BestSol.Position(1,1);
bestg=BestSol.Position(1,2);
%% 圖示尋優過程
plot(BestCost);
xlabel('Iteration');
ylabel('Best Val');
grid on;
%% 利用迴歸預測分析最佳的引數進行SVM網路訓練
cmd_cs_svr=['-s 3 -t 2',' -c ',num2str(bestc),' -g ',num2str(bestg)];
model_cs_svr=svmtrain(output_train,input_train,cmd_cs_svr); % SVM模型訓練
%% SVM網路迴歸預測
[output_test_pre,~]=svmpredict(output_test,input_test,model_cs_svr); % SVM模型預測及其精度
test_pre=mapminmax('reverse',output_test_pre',rule2);
test_pre = test_pre';

%% SVR_fitness -- objective function
function fitness=fobj(cv,input_train,output_train,input_test,output_test)
% cv為長度為2的橫向量，即SVR中引數c和v的值

cmd = ['-s 3 -t 2',' -c ',num2str(cv(1)),' -g ',num2str(cv(2))];
model=svmtrain(output_train,input_train,cmd); % SVM模型訓練
[~,fitness]=svmpredict(output_test,input_test,model); % SVM模型預測及其精度
fitness=fitness(2); % 以平均均方誤差MSE作為優化的目標函式值

function x=init_individual(xlb,xub,dim,sizepop)
% 引數初始化函式
% lb：引數下界，行向量
% ub：引數上界，行向量
% dim：引數維度
% sizepop 種群規模
% x：返回sizepop*size(lb,2)的引數矩陣
xRange=repmat((xub-xlb),[sizepop,1]);
xLower=repmat(xlb,[sizepop,1]);
x=rand(sizepop,dim).*xRange+xLower;

clear
clc
close all
load wndspd % 示例資料為風速（時間序列）資料，共144個樣本
%% PSO-SVR
% 訓練/測試資料準備（用前3天預測後一天）,用前100天做訓練資料
input_train(1,:)=wndspd(1:97);
input_train(2,:)=wndspd(2:98);
input_train(3,:)=wndspd(3:99);
output_train=[wndspd(4:100)]';
input_test(1,:)=wndspd(101:end-3);
input_test(2,:)=wndspd(102:end-2);
input_test(3,:)=wndspd(103:end-1);
output_test=[wndspd(104:end)]';
para=[30,200,0.2,0.8,0.2];
[bestc,bestg,test_pre]=my_DE_SVR(para,input_train',output_train',input_test',output_test');
% 預測誤差計算
MSE=mymse(output_test',test_pre)
MAE=mymae(output_test',test_pre)
MAPE=mymape(output_test',test_pre)
FVD=myfvd(output_test',test_pre)
CDFR=mycdfr(output_test',test_pre)
%% 預測結果圖
err_pre=output_test'-test_pre;
figure('Name','測試資料殘差圖')
set(gcf,'unit','centimeters','position',[0.5,5,30,5])
plot(err_pre,'*-');
figure('Name','原始-預測圖')
plot(test_pre,'*r-');hold on;plot(output_test,'bo-');
legend('預測','原始',0)
set(gcf,'unit','centimeters','position',[0.5,13,30,5])
toc

本文參考：http://cn.mathworks.com/matlabcentral/fileexchange/52897-differential-evolution–de-?s_tid=srchtitle

（廣告）歡迎掃描關注微信公眾號：Genlovhyy的資料小站(Gnelovy212)