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)