【风电功率预测】基于matlab帝国殖民竞争算法优化BP神经网络风电功率预测【含Matlab源

⛄一、获取代码方式

获取代码方式1：

完整代码已上传我的资源：【风电功率预测】基于matlab帝国殖民竞争算法优化BP神经网络风电功率预测【含Matlab源码 1314期】

⛄二、帝国殖民竞争算法简介

帝国竞争算法是基于人类社会政治进化现象提出的智能优化算法，算法中的每个国家均表示为一个可能的解空间，在模拟社会政治进化的过程中逐渐优化国家的解空间，找到最优解。算法主要包括帝国初始化、帝国同化、帝国革命、殖民地竞争四个步骤[18]:

a）帝国初始化。通过随机生成的方式产生多个国家，根据国家的强弱分为殖民国家与殖民地，各殖民国家及其所属殖民地构成一个帝国。

b）帝国同化。将每个帝国殖民国家的解信息传递给其殖民地，殖民地逐步靠近殖民国家，帝国整体解质量得到逐步优化。

c）帝国革命。该阶段模拟历史进程中的革命行为，一定几率下殖民地会主动优化，甚至超越其所属殖民国家成为帝国的殖民国家。

d）殖民地竞争。该阶段是殖民地再分配的过程，帝国之间通过竞争机制争夺殖民地。理想情况下只会存在一个帝国，即只有一个殖民国家，其他国家均沦为殖民地，此时算法停止，殖民国家即代表算法优化的最优解；若仍存在多个国家则转到帝国同化阶段，依次循环更新。

⛄三、BP神经网络简介

1 BP神经网络概述

BP（Back Propagation）神经网络是1986年由Rumelhart和McCelland为首的科研小组提出，参见他们发表在Nature上的论文 Learning representations by back-propagating errors 。

BP神经网络是一种按误差逆传播算法训练的多层前馈网络，是目前应用最广泛的神经网络模型之一。BP网络能学习和存贮大量的输入-输出模式映射关系，而无需事前揭示描述这种映射关系的数学方程。它的学习规则是使用最速下降法，通过反向传播来不断调整网络的权值和阈值，使网络的误差平方和最小。

2 BP算法的基本思想

上一次我们说到，多层感知器在如何获取隐层的权值的问题上遇到了瓶颈。既然我们无法直接得到隐层的权值，能否先通过输出层得到输出结果和期望输出的误差来间接调整隐层的权值呢？BP算法就是采用这样的思想设计出来的算法，它的基本思想是,学习过程由信号的正向传播与误差的反向传播两个过程组成。

正向传播时，输入样本从输入层传入,经各隐层逐层处理后,传向输出层。若输出层的实际输出与期望的输出(教师信号)不符,则转入误差的反向传播阶段。

反向传播时，将输出以某种形式通过隐层向输入层逐层反传,并将误差分摊给各层的所有单元,从而获得各层单元的误差信号,此误差信号即作为修正各单元权值的依据。这两个过程的具体流程会在后文介绍。

BP算法的信号流向图如下图所示

3 BP网络特性分析——BP三要素

我们分析一个ANN时，通常都是从它的三要素入手，即

1)网络拓扑结构；

2)传递函数；

3)学习算法。

每一个要素的特性加起来就决定了这个ANN的功能特性。所以，我们也从这三要素入手对BP网络的研究。

3.1 BP网络的拓扑结构

上一次已经说了，BP网络实际上就是多层感知器，因此它的拓扑结构和多层感知器的拓扑结构相同。由于单隐层（三层）感知器已经能够解决简单的非线性问题，因此应用最为普遍。三层感知器的拓扑结构如下图所示。

一个最简单的三层BP：

3.2 BP网络的传递函数

BP网络采用的传递函数是非线性变换函数——Sigmoid函数（又称S函数）。其特点是函数本身及其导数都是连续的，因而在处理上十分方便。为什么要选择这个函数，等下在介绍BP网络的学习算法的时候会进行进一步的介绍。

单极性S型函数曲线如下图所示。

双极性S型函数曲线如下图所示。

3.3 BP网络的学习算法

BP网络的学习算法就是BP算法，又叫 δ 算法（在ANN的学习过程中我们会发现不少具有多个名称的术语），以三层感知器为例，当网络输出与期望输出不等时，存在输出误差 E ，定义如下

下面我们会介绍BP网络的学习训练的具体过程。

4 BP网络的训练分解

训练一个BP神经网络，实际上就是调整网络的权重和偏置这两个参数，BP神经网络的训练过程分两部分：

前向传输，逐层波浪式的传递输出值；

逆向反馈，反向逐层调整权重和偏置；

我们先来看前向传输。

前向传输（Feed-Forward前向反馈）

在训练网络之前，我们需要随机初始化权重和偏置，对每一个权重取[ − 1 , 1 ] [-1,1][−1,1]的一个随机实数，每一个偏置取[ 0 , 1 ] [0,1][0,1]的一个随机实数，之后就开始进行前向传输。

神经网络的训练是由多趟迭代完成的，每一趟迭代都使用训练集的所有记录，而每一次训练网络只使用一条记录，抽象的描述如下：

while 终止条件未满足：for record:dataset:trainModel(record)

4.1 逆向反馈（Backpropagation）

4.2 训练终止条件

每一轮训练都使用数据集的所有记录，但什么时候停止，停止条件有下面两种：

设置最大迭代次数，比如使用数据集迭代100次后停止训练

计算训练集在网络上的预测准确率，达到一定门限值后停止训练

5 BP网络运行的具体流程

5.1 网络结构

输入层有n nn个神经元，隐含层有p pp个神经元，输出层有q qq个神经元。

5.2 变量定义

第九步：判断模型合理性

判断网络误差是否满足要求。

当误差达到预设精度或者学习次数大于设计的最大次数，则结束算法。

否则，选取下一个学习样本以及对应的输出期望，返回第三部，进入下一轮学习。

6 BP网络的设计

在进行BP网络的设计是，一般应从网络的层数、每层中的神经元个数和激活函数、初始值以及学习速率等几个方面来进行考虑，下面是一些选取的原则。

6.1 网络的层数

理论已经证明，具有偏差和至少一个S型隐层加上一个线性输出层的网络，能够逼近任何有理函数，增加层数可以进一步降低误差，提高精度，但同时也是网络复杂化。另外不能用仅具有非线性激活函数的单层网络来解决问题，因为能用单层网络解决的问题，用自适应线性网络也一定能解决，而且自适应线性网络的运算速度更快，而对于只能用非线性函数解决的问题，单层精度又不够高，也只有增加层数才能达到期望的结果。

6.2 隐层神经元的个数

网络训练精度的提高，可以通过采用一个隐含层，而增加其神经元个数的方法来获得，这在结构实现上要比增加网络层数简单得多。一般而言，我们用精度和训练网络的时间来恒量一个神经网络设计的好坏：

（1）神经元数太少时，网络不能很好的学习，训练迭代的次数也比较多，训练精度也不高。

（2）神经元数太多时，网络的功能越强大，精确度也更高，训练迭代的次数也大，可能会出现过拟合(over fitting)现象。

由此，我们得到神经网络隐层神经元个数的选取原则是：在能够解决问题的前提下，再加上一两个神经元，以加快误差下降速度即可。

6.3 初始权值的选取

一般初始权值是取值在(−1,1)之间的随机数。另外威得罗等人在分析了两层网络是如何对一个函数进行训练后，提出选择初始权值量级为s√r的策略，其中r为输入个数，s为第一层神经元个数。

6.4 学习速率

学习速率一般选取为0.01−0.8，大的学习速率可能导致系统的不稳定，但小的学习速率导致收敛太慢，需要较长的训练时间。对于较复杂的网络，在误差曲面的不同位置可能需要不同的学习速率，为了减少寻找学习速率的训练次数及时间，比较合适的方法是采用变化的自适应学习速率，使网络在不同的阶段设置不同大小的学习速率。

6.5 期望误差的选取

在设计网络的过程中，期望误差值也应当通过对比训练后确定一个合适的值，这个合适的值是相对于所需要的隐层节点数来确定的。一般情况下，可以同时对两个不同的期望误差值的网络进行训练，最后通过综合因素来确定其中一个网络。

7 BP网络的局限性

BP网络具有以下的几个问题：

(1)需要较长的训练时间：这主要是由于学习速率太小所造成的，可采用变化的或自适应的学习速率来加以改进。

(2)完全不能训练：这主要表现在网络的麻痹上，通常为了避免这种情况的产生，一是选取较小的初始权值，而是采用较小的学习速率。

(3)局部最小值：这里采用的梯度下降法可能收敛到局部最小值，采用多层网络或较多的神经元，有可能得到更好的结果。

8 BP网络的改进

P算法改进的主要目标是加快训练速度，避免陷入局部极小值等，常见的改进方法有带动量因子算法、自适应学习速率、变化的学习速率以及作用函数后缩法等。动量因子法的基本思想是在反向传播的基础上，在每一个权值的变化上加上一项正比于前次权值变化的值，并根据反向传播法来产生新的权值变化。而自适应学习速率的方法则是针对一些特定的问题的。改变学习速率的方法的原则是，若连续几次迭代中，若目标函数对某个权倒数的符号相同，则这个权的学习速率增加，反之若符号相反则减小它的学习速率。而作用函数后缩法则是将作用函数进行平移，即加上一个常数。

⛄四、部分源代码

% 清除环境变量

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 数组元素个数计数

ProblemParams.VarMin=repmat(ProblemParams.VarMin,1,ProblemParams.NPar); %复制矩阵，行数不变，仍然是roblemParams.VarMin，列数重复ProblemParams.NPar遍ProblemParams.VarMin

ProblemParams.VarMax=repmat(ProblemParams.VarMax,1,ProblemParams.NPar);

end

ProblemParams.SearchSpaceSize = ProblemParams.VarMax - ProblemParams.VarMin; %搜索区间

%% 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; %迭代次数

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”. 每次趋近的系数

AlgorithmParams.AssimilationAngleCoefficient = .5; % In the original paper assimilation angle coefficient is shown by “gama”. 夹角度数

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. 殖民者参数

if DisplayParams.PlotEmpires

DisplayParams.EmpiresFigureHandle = figure(‘Name’,‘Plot of Empires’,‘NumberTitle’,‘off’);

DisplayParams.EmpiresAxisHandle = axes;

end

DisplayParams.PlotCost = true; % “true” to plot. “false” 消耗参数

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);%建立国家子函数调用一个列向量 (AlgorithmParams.NumOfCountriesx1) 值为约束范围内的随机数

% 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. 按照相关耗费给国家排序

%stop 调试所用

Empires = CreateInitialEmpires(InitialCountries,InitialCost,AlgorithmParams, ProblemParams);%子函数调用得到帝国种群（殖民国加殖民地）

%% Main Loop

MinimumCost = repmat(nan,AlgorithmParams.NumOfDecades,1); %把nan复制了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

for Decade = 1:AlgorithmParams.NumOfDecades

AlgorithmParams.RevolutionRate = AlgorithmParams.DampRatio * AlgorithmParams.RevolutionRate;%进化率=0.99*原进化率

Remained = AlgorithmParams.NumOfDecades - Decadefor ii = 1:numel(Empires)%% Assimilation同化; Movement of Colonies Toward Imperialists (Assimilation Policy)Empires(ii) = AssimilateColonies(Empires(ii),AlgorithmParams,ProblemParams);%子函数调用%% Revolution; A Sudden Change in the Socio-Political Characteristics 有部分用重新生成的国家替代 Empires(ii) = RevolveColonies(Empires(ii),AlgorithmParams,ProblemParams);%子函数调用%% New Cost Evaluationif isempty(ProblemParams.CostFuncExtraParams)Empires(ii).ColoniesCost = feval(ProblemParams.CostFuncName,Empires(ii).ColoniesPosition);elseEmpires(ii).ColoniesCost = feval(ProblemParams.CostFuncName,Empires(ii).ColoniesPosition,ProblemParams.CostFuncExtraParams);end%% Empire Possession (****** Power Possession, Empire Possession)Empires(ii) = PossesEmpire(Empires(ii));%子函数调用%% Computation of Total Cost for EmpiresEmpires(ii).TotalCost = Empires(ii).ImperialistCost + AlgorithmParams.Zeta * mean(Empires(ii).ColoniesCost);endif numel(Empires) == 1 && AlgorithmParams.StopIfJustOneEmpire%如果只剩下一个殖民者breakend

%% 导出

if Decade == AlgorithmParams.NumOfDecades

positions(i,:)=[Empires(index).ImperialistPosition ];

results(i,:)=[MinimumCost(Decade)];

end

if DisplayParams.PlotCostrefreshdata(h_MinCostPlot);% refreshdata(h_MeanCostPlot);drawnow; grid on;hold on;xlabel('迭代代数/次');ylabel('目标函数值');pause(0.01); end

% if DisplayParams.PlotCost

% refreshdata(h_MinCostPlot);%刷新图片的数据

% refreshdata(h_MeanCostPlot);

% drawnow;grid on;hold on;

% xlabel(‘迭代代数/次’);ylabel(‘目标函数值’);

% pause(0.01);

% end

end % End of Algorithm

end

MinimumCost(end)

save Cost MinimumCost

save pos positions

function DisplayEmpires(Empires,AlgorithmParams,ProblemParams,DisplayParams)

if ~DisplayParams.PlotEmpiresreturn;endif (ProblemParams.NPar ~= 2) && (ProblemParams.NPar ~= 3)return;endif ~any(findall(0)==DisplayParams.EmpiresFigureHandle)return;endif ProblemParams.NPar == 2for ii = 1:numel(Empires)plot(DisplayParams.EmpiresAxisHandle,Empires(ii).ImperialistPosition(1),Empires(ii).ImperialistPosition(2),'p',...'MarkerEdgeColor','k',...'MarkerFaceColor',DisplayParams.ColorMatrix(ii,:),...'MarkerSize',70*numel(Empires(ii).ColoniesCost)/AlgorithmParams.NumOfAllColonies + 13);hold onplot(DisplayParams.EmpiresAxisHandle,Empires(ii).ColoniesPosition(:,1),Empires(ii).ColoniesPosition(:,2),'ok',...'MarkerEdgeColor','k',...'MarkerFaceColor',DisplayParams.ColorMatrix(ii,:),...'MarkerSize',8);endxlim([DisplayParams.AxisMargin.Min(1) DisplayParams.AxisMargin.Max(1)]);ylim([DisplayParams.AxisMargin.Min(2) DisplayParams.AxisMargin.Max(2)]);hold offendif ProblemParams.NPar == 3figure(1)for ii = 1:numel(Empires)plot3(DisplayParams.EmpiresAxisHandle,Empires(ii).ImperialistPosition(1),Empires(ii).ImperialistPosition(2),Empires(ii).ImperialistPosition(3),'p',...'MarkerEdgeColor','k',...'MarkerFaceColor',DisplayParams.ColorMatrix(ii,:),...'MarkerSize',70*numel(Empires(ii).ColoniesCost)/AlgorithmParams.NumOfAllColonies + 13);hold onplot3(DisplayParams.EmpiresAxisHandle,Empires(ii).ColoniesPosition(:,1),Empires(ii).ColoniesPosition(:,2),Empires(ii).ColoniesPosition(:,3),'ok',...'MarkerEdgeColor','k',...'MarkerFaceColor',DisplayParams.ColorMatrix(ii,:),...'MarkerSize',8);endxlim([DisplayParams.AxisMargin.Min(1) DisplayParams.AxisMargin.Max(1)]);ylim([DisplayParams.AxisMargin.Min(2) DisplayParams.AxisMargin.Max(2)]);zlim([DisplayParams.AxisMargin.Min(3) DisplayParams.AxisMargin.Max(3)]);hold offendpause(0.05);

end