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);
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
%% 个体极值和群体极值
[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
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);
% Y矩阵，行表示企业编号，列表示二氧化碳百分比，其中为百分数
Y1=[33.08;
21.85;
6.19;
11.77;
9.96;
17.15;];
Y=Y1./100;%将百分数化为小数
[ym,yn]=size(Y);
%% 代入X解向量，X
XX=X';%将矩阵转置
one=ones(ym,yn);
newx=one-XX;
%% 计算基尼系数G
G=zeros(an,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，基尼系数
评论30
最新资源