数学建模社区-数学中国

标题: 一种基于伸缩因子的基础PSO算法程序 [打印本页]

作者: 夜雨声烦 时间: 2016-4-26 21:41
标题: 一种基于伸缩因子的基础PSO算法程序
function PSOfirst()
%% 清空环境
clear;
clc;

%% 参数设置
w=0.9;%权值将影响PSO 的全局与局部搜优能力，值较大，全局搜优能力强，局部搜优能力弱;反之，则局部搜优能力增强，而全局搜优能力减弱。
c1=0.1;%加速度，影响收敛速度
c2=0.1;
dim=6;%6维，表示企业数量
swarmsize=100;%粒子群规模，表示有100个粒子
maxiter=200;%最大迭代次数，影响时间
minfit=0.001;%最小适应值
vmax=0.01;%最大速度
vmin=-0.01;%最小速度
ub=[0.2,0.2,0.2,0.2,0.2,0.2];%解向量的最大限制
lb=[0.01,0.01,0.01,0.01,0.01,0.01];%解向量的最小限制

%% 种群初始化
range=ones(swarmsize,1)*(ub-lb);%产生200个粒子的初始坐标，初始解位置
swarm=rand(swarmsize,dim).*range+ones(swarmsize,1)*lb;%粒子群位置矩阵,每行表示一组解
Y1=[33.08;
21.85;
6.19;
11.77;
9.96;
17.15;];
Y=Y1./100;%将百分数化为小数
[ym,yn]=size(Y);
for i=1:swarmsize %% YX的约束
s=swarm(i,

;
ss=s';
while sum(Y.*ss)<0.1*sum(Y)
ss=rand(dim,1).*((ub-lb)')+ones(dim,1).*((lb)');
end
swarm(i,

=ss';
end
vstep=rand(swarmsize,dim)*(vmax-vmin)+vmin;%粒子群速度矩阵
fswarm=zeros(swarmsize,1);%预设空矩阵，存放适应值
%% 计算初始种群适应度
for i=1:swarmsize
X=swarm(i,

;
[SUMG,G]=jn(X);
fswarm(i,

=SUMG;
%fswarm(i,

=feval(jn,swarm(i,

);%以粒子群位置的第i行为输入，求函数值，对应输出给适应值
end
fswarm

%% 个体极值和群体极值
[bestf,bestindex]=min(fswarm);%求得适应值中的最小适应值，和，其所在的序列
gbest=swarm;%暂时的个体最优解为自己
fgbest=fswarm;%暂时的个体最优适应值
zbest=swarm(bestindex,

;%所在序列的对应的解矩阵序列，全局最佳解
fzbest=bestf;%全局最优适应值

%% 迭代寻优
iter=0;
yfitness=zeros(1,maxiter);%1行100列矩阵，存放100个最优值的空间矩阵
x1=zeros(1,maxiter);%存放x的空间
x2=zeros(1,maxiter);
x3=zeros(1,maxiter);
x4=zeros(1,maxiter);
x5=zeros(1,maxiter);
x6=zeros(1,maxiter);
while((iter<maxiter)&&(fzbest>minfit))
for j=1:swarmsize
% 速度更新
vstep(j,

=w*vstep(j,

+c1*rand*(gbest(j,

-swarm(j,

)+c2*rand*(zbest-swarm(j,

);
if vstep(j,

>vmax
vstep(j,

=vmax;%速度限制
end
if vstep(j,

<vmin
vstep(j,

=vmin;
      end
      % 位置更新
      swarm(j,

=swarm(j,

+vstep(j,

;
      for k=1:dim
         if swarm(j,k)>ub(k)
            swarm(j,k)=ub(k);%位置限制
         end
         if swarm(j,k)<lb(k)
            swarm(j,k)=lb(k);
         end
      end

      % 适应值
      X=swarm(j,

;
[SUMG,G]=jn(X);
fswarm(j,

=SUMG;
      % 可在此处增加约束条件，若满足约束条件，则进行适应值计算

      %
      % 个体最优更新
      if fswarm(j)<fgbest(j) %如果当前的函数值比个体最优值小
         gbest(j,

=swarm(j,

;%个体最优解更新
         fgbest(j)=fswarm(j);%个体最优值更新
      end
      % 群体最优更新
      if fswarm(j)<fzbest%如果当前的函数值比群体最优值大
         zbest=swarm(j,

;%群体最优解更新
         fzbest=fswarm(j);%群体最优值更新
      end
end
iter=iter+1;
yfitness(1,iter)=fzbest;
x1(1,iter)=zbest(1);%将全局最优解的第1个元素，依次存储，共有MAXITER个
x2(1,iter)=zbest(2);
x3(1,iter)=zbest(3);
x4(1,iter)=zbest(4);
x5(1,iter)=zbest(5);
x6(1,iter)=zbest(6);
end
min(yfitness)
fzbest
zbest
X=zbest;
[SUMG,G]=jn(X);
GGbest=G;GGbest
%% 画图
figure(1)
plot(yfitness,'linewidth',2)
title('最优基尼系数优化曲线','fontsize',14);
xlabel('迭代次数','fontsize',14);
ylabel('基尼系数','fontsize',14);

figure(2)
plot(x1,'b')
hold on
plot(x2,'g')
hold on
plot(x3,'r')
hold on
plot(x4,'c')
hold on
plot(x5,'m')
hold on
plot(x6,'y')
title('x优化曲线','fontsize',14);
xlabel('迭代次数','fontsize',14);
ylabel('参数值','fontsize',14);
legend('x1','x2','x3','x4','x5','x6',88)

%% 适应度函数,即为目标函数，这里为基尼系数函数
function [SUMG,G]=jn(X)
%% 已知数据
% A矩阵，行表示企业编号，列表示员工、营业收入、税收总额，其中数据位百分数
A1=[ 30.8 59.2 39.92;
17.6 9.5  31.42;
13.6 7.1  6.62;
9.5  7 5.64;
23.8 5.8  4.79;
4.7  11.4 11.6;];
A=A1./100;%将百分数化为小数
[am,an]=size(A);%am=6;an=3
% Y矩阵，行表示企业编号，列表示二氧化碳百分比，其中为百分数
Y1=[33.08;
21.85;
6.19;
11.77;
9.96;
17.15;];
Y=Y1./100;%将百分数化为小数
[ym,yn]=size(Y);%ym=6;yn=1
%% 代入X解向量，X为1行6列向量
XX=X';%将矩阵转置
one=ones(ym,yn);
newx=one-XX;%1减去对应位置的解
%% 计算基尼系数G
G=zeros(an,1);%3行1列
for j=1:an
aj=A(:,j);
yx1=Y.*newx;
yx=yx1./sum(yx1);
ya=yx./aj;
compose=[ya,aj,yx;];
newm=sortrows(compose,1);%将ya矩阵从小到大升序排列；
ajnew=newm(:,2);
yxnew=newm(:,3);
yxnewsum=zeros(ym,yn);
for ii=1:ym
      yxnewsum(ii,yn)=sum(yxnew(1:ii));
end
yxnewsum2=zeros(ym,yn);
for iii=1:ym
      if iii==1
         yxnewsum2(iii,yn)=yxnewsum(iii,yn);
      else
      yxnewsum2(iii,yn)=yxnewsum(iii-1,yn)+yxnewsum(iii,yn);
      end
end
ay=ajnew.*yxnewsum2;
gj=1-sum(ay);
G(j)=gj;
end
GMAX=[0.3;0.3;0.2;];
if ((G(1)-GMAX(1)>0)||(G(2)-GMAX(2)>0)||(G(3)-GMAX(3)>0))
G=GMAX;
end
SUMG=0.61*G(1)+0.19*G(2)+0.2*G(3);
%输出G，基尼系数

作者: 夜雨声烦 时间: 2016-4-26 21:43
这好像有问题，帖子中笑脸应该是“

”的，大家使用时注意一下！

作者: 成哥cc 时间: 2016-4-30 20:18
00000000000000000

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)