matlab-粒子群源码优化模糊隶属度函数值

Posted 2022-11-12 studyer_domi

1、内容简介

略
594-可以交流、咨询、答疑

粒子群源码优化模糊隶属度函数值

2、内容说明

定义
粒子群优化算法(Particle Swarm optimization,PSO)又翻译为粒子群算法、微粒群算法、或微粒群优化算法。是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。通常认为它是群集智能 (Swarm intelligence, SI) 的一种。它可以被纳入多主体优化系统(Multiagent Optimization System, MAOS)。粒子群优化算法是由Eberhart博士和kennedy博士发明。
模拟捕食
PSO模拟鸟群的捕食行为。一群鸟在随机搜索食物，在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻离食物最近的鸟的周围区域。
启示
PSO从这种模型中得到启示并用于解决优化问题。PSO中，每个优化问题的解都是搜索空间中的一只鸟。我们称之为“粒子”。所有的粒子都有一个由被优化的函数决定的适应值(fitnessvalue)，每个粒子还有一个速度决定他们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。
PSO初始化
PSO初始化为一群随机粒子(随机解)，然后通过迭代找到最优解，在每一次迭代中，粒子通过跟踪两个“极值”来更新自己。第一个就是粒子本身所找到的最优解，这个解叫做个体极值pBest，另一个极值是整个种群找到的最优解，这个极值是全局极值gBest。另外也可以不用整个种群而只是用其中一部分最优粒子的邻居，那么在所有邻居中的极值就是局部极值。

3、仿真分析

input.init.saved_veh_file='BATTERY_CAPACITOR_defaults_in';
[a,b]=adv_no_gui('initialize',input);
%%
c1=1.5;%学习因子
c2=1.5;%学习因子
maxgen=100; %进化次数
sizepop=40; %种群规模
%粒子更新速度
Vmax=0.8;%最大前进速度
Vmin=-0.8; %最小前进速度
elpha=1.5;
%随机确定种群
popmax=[ 0.25 0.55 0.15 0.5 0.15 0.6 0.15 0.8 0.75 1 0.45 0.75 0.35 0.75 0.75 0.9 0.45 0.8 0.4 0.85 0.8 1 0.1 0.25 0.1 0.4 0.1 0.6 0.1 0.8 0.9 1 0.25 0.7 0.25 0.65 0.7 1 0.4 0.8 0.4 0.75 0.8 0.9 0.35 0.55 0.3 0.55 0.3 0.8 0.85 1 0.2 0.3 0.25 0.55 0.2 0.85 0.9 1];
popmin=[0 0.35 0 0.2 0 0.4 0 0.6 0.5 0.7 0.22 0.45 0.15 0.35 0.4 0.65 0.2 0.45 0.15 0.45 0.45 0.75 0 0.15 0 0.35 0 0.45 0 0.65 0.7 0.78 0 0.28 0 0.35 0.4 0.7 0.25 0.4 0 0.49 0.5 0.85 0.2 0.25 0 0.36 0 0.55 0.62 0.83 0 0.15 0 0.30 0 0.5 0.55 0.85];
par_num=60;
%产生初始粒子和速度
Range=repmat((popmax-popmin),[sizepop 1]);%remat 是复制和平铺矩阵
Lower=repmat(popmin, [sizepop 1]); pop=rand(sizepop,par_num).*Range+Lower;
V=rands(sizepop,par_num); %初始化速度
for i=1:sizepop
    %随机产生一个种群，并计算适应度
    x1=pop(i,1);x2=pop(i,2);x3=pop(i,3);x4=pop(i,4);x5=pop(i,5);x6=pop(i,6);x7=pop(i,7);x8=pop(i,8);x9=pop(i,9); x10=pop(i,10);
    y1=pop(i,11);y2=pop(i,12);y3=pop(i,13);y4=pop(i,14);y5=pop(i,15);y6=pop(i,16);
    z1=pop(i,17);z2=pop(i,18);z3=pop(i,19);z4=pop(i,20);z5=pop(i,21);z6=pop(i,22);
    a1=pop(i,23);a2=pop(i,24);a3=pop(i,25);a4=pop(i,26);a5=pop(i,27);a6=pop(i,28);a7=pop(i,29);a8=pop(i,30);a9=pop(i,31);a10=pop(i,32);
    b1=pop(i,33);b2=pop(i,34);b3=pop(i,35);b4=pop(i,36);b5=pop(i,37);b6=pop(i,38);
    c1=pop(i,39);c2=pop(i,40);c3=pop(i,41);c4=pop(i,42);c5=pop(i,43);c6=pop(i,44);
    d1=pop(i,45);d2=pop(i,46);d3=pop(i,47);d4=pop(i,48);d5=pop(i,49);d6=pop(i,50);d7=pop(i,51);d8=pop(i,52);
    %输入驱动功率 Preq
    a=addvar(a,'input','Preq',[0,1]); %给FIS添加新的变量 Preq
    a=addmf(a,'input',1,'TS','trap',[x1,x2]); %将隶属度函数TS添加到 FIS
    a=addmf(a,'input',1,'S','gaussmf',[x3,x4]);
    a=addmf(a,'input',1,'M','gaussmf',[x5,x6]);
    a=addmf(a,'input',1,'B','gaussmf',[x7,x8]);
    a=addmf(a,'input',1,'TB','trap',[x9,x10]);
    %输入电池电量 SOCbat
    a=addvar(a,'input','SOCbat',[0.2,0.9]);
    a=addmf(a,'input',2,'L','zmf',[y1,y2]);
    a=addmf(a,'input',2,'M','gaussmf',[y3,y4]);
    a=addmf(a,'input',2,'H','smf',[y5,y6]); 
    %输入电容电量 SOCuc 
    a=addvar(a,'input','SOCuc',[0.2,1]);
    a=addmf(a,'input',3,'L','zmf',[z1,z2]);
    a=addmf(a,'input',3,'M','gaussmf',[z3,z4]);
    a=addmf(a,'input',3,'H','smf',[z5,z6]);
    %输出比例因子 Kuc
    a=addvar(a,'output','Kuc',[0,1]); 
    a=addmf(a,'output',1,'TS','trap',[a1,a2]);
    a=addmf(a,'output',1,'S','gaussmf',[a3,p4]);
    a=addmf(a,'output',1,'M','gaussmf',[a5,a6]);
    a=addmf(a,'output',1,'B','gaussmf',[a7,a8]); 
    a=addmf(a,'output',1,'TB','trap',[a9,a10]);
    %输入 rulelist 模糊规则
    rulelist=[1 1 2 4 1 1;3 1 1 4 2 1;2 1 1 4 1 1;5 1 1 2 1 1;3 1 1 2 1 1;1 2 1 4 1 1;4 2 1 4 1 1;2 2 1 4 1 1;
              5 2 1 4 1 1;3 2 1 2 1 1;1 3 1 4 1 1;4 3 1 4 1 1;2 3 1 4 1 1;5 3 1 4 1 1;3 3 1 4 1 1;1 1 2 4 1 1;
              4 1 2 2 1 1;2 1 2 1 1 1;5 1 2 3 1 1;3 1 2 3 1 1;1 2 2 4 1 1;4 2 2 4 1 1;2 2 2 2 1 1;5 2 2 1 1 1;
              3 2 2 3 1 1;1 3 2 4 1 1;4 3 2 4 1 1;2 3 2 4 1 1;5 3 2 2 1 1;3 3 2 1 1 1;1 1 3 2 1 1;4 1 3 1 1 1;
              2 1 3 3 1 1;5 1 3 5 1 1;3 1 3 5 1 1;1 2 3 4 1 1;4 2 3 2 1 1;2 2 3 1 1 1;5 2 3 3 1 1;3 2 3 5 1 1;
              1 3 3 4 1 1;4 3 3 4 1 1;2 3 3 2 1 1;5 3 3 1 1 1;3 3 3 3 1 1];
    a=addrule(a,rulelist);%FIS 中添加模糊规则
    a1=setfis(a,'DefuzzMethod','mom'); %采用极大平均法进行反模糊化
    writefis(a1,'ll');
    A1=readfis('ll');
    a=newfis('JJ');
    %输入电池电量 SOCbat
    a=addvar(a,'input',SOCbat',[0.2,0.9]);
    a=addmf(a,'input',1,'L','zmf',[b1,b2]);
    a=addmf(a,'input',1,'M','gaussmf',[b3,b4]);
    a=addmf(a,'input',1,'H','smf',[b5,b6]);
    %输入电容电量 SOCuc
    a=addvar(a,'input','SOCuc',[0.2,1]);
    a=addmf(a,'input',2,'L','gaussmf',[c1,c2]);
    a=addmf(a,'input',2,'M','gaussmf',[c3,c4]);
    a=addmf(a,'input',2,'H','smf',[c5,c6]);
    a=addmf(a,'input',2,'TL','zmf',[c7,c8]);
    %输出比例因子 Kuc
    a=addvar(a,'output','Kuc',[0,1]);
    a=addmf(a,'output',1,'TS','gaussmf',[d1,d2]);
    a=addmf(a,'output',1,'S','gaussmf',[d3,d4]);
    a=addmf(a,'output',1,'M','gaussmf',[d5,d6]);
    a=addmf(a,'output',1,'B','smf',[d7,d8]);
    %输入 rulelist
    rulelist=[2 2 2 4 1 1;1 2 2 1 1 1;3 2 2 2 1 1;2 3 2 4 1 1;1 3 2 4 1 1;3 3 2 1 1 1;2 1 3 3 1 1;
        1 1 3 3 1 1;3 1 3 3 1 1;2 2 3 3 1 1;1 2 3 2 1 1;3 2 3 2 1 1];
    a=addrule(a,rulelist);
    a1=setfis(a,'DefuzzMethod','mom');
    writefis(a1,'JJ');
    A2=readfis('JJ');
    input.cycle.param='cycle.name','cycle.soc','cyc.number';
    input.cycle.value='CYC_UDDS','off',1%输入工况 NEDC
    [a,b]=adv_no_gui('drive_cycle',input)%工况驱动
    input.resp.param='distance','ess_used_all';%行驶里程和总体能量消耗
    [a,b]=adv_no_gui('other_info',input)
    D=max(b.other.value1);
    L=max(b.other.value2); fitness(i)=L/D; %当前最小适应度
end
%寻找最佳的适应度值
[bestfitness bestindex]=min(fitness);
zbest=pop(bestindex,:); %全局最佳
gbest=pop; %个体最佳
fitnessgbest=fitness; %个体最佳适应度值
fitnesszbest=bestfitness; %全局最佳适应度值
%迭代寻优
for ii=1:maxgen
for j=1:sizepop poop=pop(j,:);
    k=fix(rand*sizepop)+1;
    while(j==k)
        k=fix(rand*sizepop)+1;
    end
    jj=fix(rand*par_num)+1;
    %进行邻域搜索
    poop(1,jj)=pop(j,jj)+(rand-0.5)*elpha*(pop(j,jj)-pop(k,jj));
    V(j,:) = V(j,:)+c1*rand*(zbest(1,:)-poop(1,:))+c2*rand*(gbest(j,:)-poop(1,:));
    V(j,find(V(j,:)>Vmax))=Vmax;
    V(j,find(V(j,:)<Vmin))=Vmin;
    %种群进行更新
    poop(1,:)=poop(1,:)+0.5*V(j,:);
    if poop(1,jj)>popmax(jj)
        poop(1,jj)=popmax(jj);
    end
    if poop(1,jj)<popmin(jj)
        poop(1,jj)=popmin(jj);
    end
    if rand>0.82
        k=ceil(par_num*rand);
        pop(j,jj)=rand;
    end
    %适应度值
    x1=pop(i,1);x2=pop(i,2);x3=pop(i,3);x4=pop(i,4);x5=pop(i,5);x6=pop(i,6);x7=pop(i,7);x8=pop(i,8);x9=pop(i,9); x10=pop(i,10);
    y1=pop(i,11);y2=pop(i,12);y3=pop(i,13);y4=pop(i,14);y5=pop(i,15);y6=pop(i,16);
    z1=pop(i,17);z2=pop(i,18);z3=pop(i,19);z4=pop(i,20);z5=pop(i,21);z6=pop(i,22);
    a1=pop(i,23);a2=pop(i,24);a3=pop(i,25);a4=pop(i,26);a5=pop(i,27);a6=pop(i,28);
    a7=pop(i,29);a8=pop(i,30);a9=pop(i,31);a10=pop(i,32);
    b1=pop(i,33);b2=pop(i,34);b3=pop(i,35);b4=pop(i,36);b5=pop(i,37);b6=pop(i,38);
    c1=pop(i,39);c2=pop(i,40);c3=pop(i,41);c4=pop(i,42);c5=pop(i,43);c6=pop(i,44);
    d1=pop(i,45);d2=pop(i,46);d3=pop(i,47);d4=pop(i,48);d5=pop(i,49);d6=pop(i,50);
    d7=pop(i,51);d8=pop(i,52);
    a=newfis('ll');
    %输入驱动功率 Preq
    a=addvar(a,'input','Preq',[0,1]);
    a=addmf(a,'input',1,'TS','trap',[x1,x2]);
    a=addmf(a,'input',1,'S','gaussmf',[x3,x4]);
    a=addmf(a,'input',1,'M','gaussmf',[x5,x6]);
    a=addmf(a,'input',1,'B','gaussmf',[x7,x8]);
    a=addmf(a,'input',1,'TB','trap',[x9,x10]);
    %输入电池电量 SOCbat
    a=addvar(a,'input','SOCbat',[0.2,0.9]); 
    a=addmf(a,'input',2,'L','zmf',[y1,y2]);
    a=addmf(a,'input',2,'M','gaussmf',[y3,y4]);
    a=addmf(a,'input',2,'H','smf',[y5,y6]);
    %输入电容电量 SOCuc
    a=addvar(a,'input','SOCuc',[0.2,1]);
    a=addmf(a,'input',3,'L','zmf',[z1,z2]);
    a=addmf(a,'input',3,'M','gaussmf',[z3,z4]); 
    a=addmf(a,'input',3,'H','smf',[z5,z6]);
    %输出比例因子 Kuc
    a=addvar(a,'output','Kuc',[0,1]);
    a=addmf(a,'output',1,'TS','trap',[p1,p2]);
    a=addmf(a,'output',1,'S','gaussmf',[p3,p4]);
    a=addmf(a,'output',1,'M','gaussmf',[p5,p6]);
    a=addmf(a,'output',1,'B','gaussmf',[p7,p8]);
    a=addmf(a,'output',1,'TB','trap',[p9,p10]);
    %输入 rulelist
    rulelist=[1 1 2 4 1 1;3 1 1 4 2 1;2 1 1 4 1 1;5 1 1 2 1 1;3 1 1 2 1 1;1 2 1 4 1 1;4 2 1 4 1 1;2 2 1 4 1 1;
              5 2 1 4 1 1;3 2 1 2 1 1;1 3 1 4 1 1;4 3 1 4 1 1;2 3 1 4 1 1;5 3 1 4 1 1;3 3 1 4 1 1;1 1 2 4 1 1;
              4 1 2 2 1 1;2 1 2 1 1 1;5 1 2 3 1 1;3 1 2 3 1 1;1 2 2 4 1 1;4 2 2 4 1 1;2 2 2 2 1 1;5 2 2 1 1 1;
              3 2 2 3 1 1;1 3 2 4 1 1;4 3 2 4 1 1;2 3 2 4 1 1;5 3 2 2 1 1;3 3 2 1 1 1;1 1 3 2 1 1;4 1 3 1 1 1;
              2 1 3 3 1 1;5 1 3 5 1 1;3 1 3 5 1 1;1 2 3 4 1 1;4 2 3 2 1 1;2 2 3 1 1 1;5 2 3 3 1 1;3 2 3 5 1 1;
              1 3 3 4 1 1;4 3 3 4 1 1;2 3 3 2 1 1;5 3 3 1 1 1;3 3 3 3 1 1];
    a=addrule(a,rulelist);
    a1=setfis(a,'DefuzzMethod','mom');
    writefis(a1,'ll');
    A1=readfis('ll');
    a=newfis('JJ');
    %输入电池电量 SOCbat
    a=addvar(a,'input','SOCbat',[0.2,0.9]);
    a=addmf(a,'input',1,'L','zmf',[b1,b2]);
    a=addmf(a,'input',1,'M','gaussmf',[b3,b4]);
    a=addmf(a,'input',1,'H','smf',[b5,b6]);
    %输入电容电量 SOCuc
    a=addvar(a,'input','SOCuc',[0.2,1]); 
    a=addmf(a,'input',2,'L','gaussmf',[c1,c2]);
    a=addmf(a,'input',2,'M','gaussmf',[c3,c4]);
    a=addmf(a,'input',2,'H','smf',[c5,c6]);
    a=addmf(a,'input',2,'TL','zmf',[c7,c8]);
    %输出比例因子 Kuc
    a=addvar(a,'output','Kuc',[0,1]);
    a=addmf(a,'output',1,'TS','gaussmf',[d1,d2]);
    a=addmf(a,'output',1,'S','gaussmf',[d3,d4]);
    a=addmf(a,'output',1,'M','gaussmf',[d5,d6]);
    a=addmf(a,'output',1,'B','smf',[d7,d8]);
    %输入 rulelist
    rulelist=[2 2 2 4 1 1;1 2 2 1 1 1;3 2 2 2 1 1;2 3 2 4 1 1;1 3 2 4 1 1;3 3 2 1 1 1;2 1 3 3 1
    1;1 1 3 3 1 1;3 1 3 3 1 1;2 2 3 3 1 1;1 2 3 2 1 1;3 2 3 2 1 1];
    a=addrule(a,rulelist);
    a1=setfis(a,'DefuzzMethod','mom');
    writefis(a1,'JJ');
    A2=readfis('JJ');
    input.cycle.param='cycle.name','cycle.soc','cyc.number';
    input.cycle.value='CYC_UDDS','off',1
    [a,b]=adv_no_gui('drive_cycle',input)
    input.resp.param='distance','ess_used_all';
    [a,b]=adv_no_gui('other_info',input) 
    D=max(b.other.value1);
    L=max(b.other.value2);
if (length(ess_used_all) ==1380)
fitness(i)=max(ess_used_all);
else fitness(i)=80000000000000000000;
end
fitness=L/D; %当前最小适应度
fitness(i)=((1*cos((1+1)*pop(i,1)+1))+(2*cos((2+1)*+2))+(3*cos((3+1)*pop(i,1)+3))+(4*cos((4+1)*pop(i,1)+4))+(5*cos((5+1)*pop(i,1)+5))).*((1*cos((1+1)*pop(i,2)+1))+(2*cos((2+1)*pop(i,2)+2))+(3*cos((3+1)*pop(i,2)+3))+(4*cos((4+1)*pop(i,2)+4))+( 5*cos((5+1)*pop(i,2)+5)))
end
end
%个体最优粒子更新
if fitness<fitnessgbest(j)
gbest(j,:)=pop(j,:);
fitnessgbest(j)=fitness;
end
%群体最优粒子更新
if fitness<fitnesszbest
zbest=pop(j,:);
fitnesszbest=fitness;
end

yy(ii)=fitnesszbest;

plot(yy,'linewidth',2);
title(['适应度度曲线' '终止代数= ' num2str(maxgen)]); xlabel('进化代数');ylabel('适应度');
grid on

4、参考论文

略