数学建模社区-数学中国

标题: [源码] 基于粒子群算法的电力系统机组组合问题的程序 [打印本页]

作者: wan121036604 时间: 2013-1-26 21:09
标题: [源码] 基于粒子群算法的电力系统机组组合问题的程序
机组组合是电力系统优化运行的一个重要方面，由于它能带来显著的经济效益，所以一直是现代电力系统每天运行计划主要的优化任务。从数学角度讲，机组组合问题是一个多约束 NP 难组合优化问题，很难得到理论上的最优解。到目前为止仍不存在既能全面考虑各种约束，又能获得理想的运算速度和精度的实用算法。如何提高求解机组组合问题的精度和速度对电力系统的经济运行仍具有重要的意义。
以下是一个算例，希望可以起到抛砖引玉的作用：
上传下载附件（相关数据）

以下是改编的程序：（程序来自《MATLAB智能算法30案例分析》（作者：史峰等））

%% 电力系统机组组合问题粒子群优化算法
%清空环境
clc
clear
tic

%% 初始参数
weight=500;%总重量限制
xs=[0.00156,0.00194,0.00482
7.92,7.94,7.97
561,7.85,78
300,310,150
0.0315,0.0420,0.0630
100,100,50];
%初始化程序
Dim=3;    %粒子维数
xSize=200; %种群个数
MaxIt=1000; %迭代次数
c1=0.8; %算法参数
c2=0.8; %算法参数
wmax=1.2;  %惯性因子
wmin=0.1;  %惯性因子
objnum=xSize; %类中物品个数？
x1=[600*rand(xSize,1),400*rand(xSize,1),200*rand(xSize,1)];
x2=[100*ones(xSize,1),100*ones(xSize,1),50*ones(xSize,1)];
x=x1;
x(x2>x1)=0;  %粒子初始化,随机产生小于四的组合
v=zeros(xSize,Dim);    %速度初始化

xbest=x;          %个体最佳值
gbest=x(1,:);    %粒子群最佳位置

% 粒子适应度值
px=zeros(1,xSize); %粒子价值目标
cx=zeros(1,xSize); %重量约束

% 最优值初始化
pxbest=zeros(1,xSize); %粒子最优价值目标
cxbest=zeros(1,xSize);  %记录重量，以求约束

% 上一次的值
pxPrior=zeros(1,xSize);%粒子价值目标
cxPrior=zeros(1,xSize);%记录重量，以求约束

%计算初始目标向量
for i=1:xSize%60  每一个粒子群的解集
for j=1:Dim %控制类别
      px(i) = px(i)+xs(1,j)*x(i,j)^2+xs(2,j)*x(i,j)+xs(3,j)+abs(xs(4,j)*sin(xs(5,j)*(x(i,j)-xs(6,j))));  %粒子价值
      cx(i) = cx(i)+x(i,j);  %粒子重量
end
end
% 粒子最优位置
pxbest=px;cxbest=cx;

%% 初始筛选非劣解
flj=[];
fljx=[];
fljNum=0;
%两个实数相等精度
for i=1:xSize
flag=0;  %支配标志
for j=1:xSize
      if j~=i
         if (px(i)<px(j))&&(cx(i)>=weight)
            flag=1;
            break;
         end
      end
end

%判断有无被支配
if flag==1
      fljNum=fljNum+1;
      % 记录非劣解
      flj(fljNum,1)=px(i);flj(fljNum,2)=cx(i);
      % 非劣解位置
      fljx(fljNum,:)=x(i,:);
end
end
%% 循环迭代
for iter=1:MaxIt

% 权值更新
w=wmax-(wmax-wmin)*iter/MaxIt;

%从非劣解中选择粒子作为全局最优解
s=size(fljx,1);
index=randi(s,1,1);
gbest=fljx(index,:);

%% 群体更新
for i=1:xSize
      %速度更新
      v(i,:)=w*v(i,:)+c1*rand(1,1)*(xbest(i,:)-x(i,:))+c2*rand(1,1)*(gbest-x(i,:));

      %位置更新
      x(i,:)=x(i,:)+v(i,:);%?????????????????????????????????????????????
      x(i,:) = rem(x(i,:),objnum)/double(objnum);
      index1=find(x(i,:)<=0);
      if ~isempty(index1)
         x(i,index1)=rand(size(index1));
      end
      x(i,:)=ceil(3*x(i,:));
end

%% 计算个体适应度
pxPrior(:)=0;
cxPrior(:)=0;
for i=1:xSize%60  每一个粒子群的解集
for j=1:Dim %控制类别
      pxPrior(i) = pxPrior(i)+xs(1,j)*x(i,j)^2+xs(2,j)*x(i,j)+xs(3,j)+abs(xs(4,j)*sin(xs(5,j)*(x(i,j)-xs(6,j))));  %粒子价值
      cxPrior(i) = cxPrior(i)+x(i,j);  %粒子重量
end
end
%% 更新粒子历史最佳
for i=1:xSize
      %现在的支配原有的，替代原有的
      if (px(i)>pxPrior(i))&&(cxPrior(i)>=weight)
            xbest(i,:)=x(i,:);%没有记录目标值
            pxbest(i)=pxPrior(i);cxbest(i)=cxPrior(i);
      end

      %彼此不受支配，随机决定
      if ~(px(i)>pxPrior(i))&&(cxPrior(i)>=weight)
         if rand(1,1)<0.5
            xbest(i,:)=x(i,:);
               pxbest(i)=pxPrior(i);cxbest(i)=cxPrior(i);
         end
      end
end

%% 更新非劣解集合
px=pxPrior;
cx=cxPrior;
%更新升级非劣解集合
s=size(flj,1);%目前非劣解集合中元素个数

%先将非劣解集合和xbest合并
pppx=zeros(1,s+xSize);
cccx=zeros(1,s+xSize);
pppx(1:xSize)=pxbest;pppx(xSize+1:end)=flj(:,1)';
cccx(1:xSize)=cxbest;cccx(xSize+1:end)=flj(:,2)';
xxbest=zeros(s+xSize,Dim);
xxbest(1:xSize,:)=xbest;
xxbest(xSize+1:end,:)=fljx;

%筛选非劣解
flj=[];
fljx=[];
k=0;

for i=1:xSize+s
      flag=0;%没有被支配
      %判断该点是否非劣
      for j=1:xSize+s
         if j~=i
            if ((pppx(i)<pppx(j)) &&  (cccx(i)>=weight)) %有一次被支配
                  flag=1;
                  break;
            end
         end
      end

      %判断有无被支配
      if flag==1
         k=k+1;
         flj(k,1)=pppx(i);flj(k,2)=cccx(i);%记录非劣解
         fljx(k,:)=xxbest(i,:);%非劣解位置
      end
end

%    去掉重复粒子
repflag=0; %重复标志
k=1;       %不同非劣解粒子数
flj2=[];    %存储不同非劣解
fljx2=[]; %存储不同非劣解粒子位置
%    flj21=[];
%    fljx21=[];
flj2(k,:)=flj(1,:);
fljx2(k,:)=fljx(1,:);
for j=2:size(flj,1)
      repflag=0;  %重复标志
      for i=1:size(flj2,1)
         result=(fljx(j,:)==fljx2(i,:));
         if length(find(result==1))==Dim
            repflag=1;%有重复
         end
      end
%       粒子不同，存储
      if repflag==0
         k=k+1;
         flj2(k,:)=flj(j,:);
         fljx2(k,:)=fljx(j,:);
      end

end

%    非劣解更新
flj=flj2;
fljx=fljx2;

end

%绘制非劣解分布
plot(flj(:,1),flj(:,2),'o')
xlabel('P')
ylabel('R')
title('最终非劣解在目标空间分布')
disp('非劣解flj中两列依次为P，C')
flj
fljx
toc

360截图20130123161902665.jpg (186.72 KB, 下载次数: 579)