一、核心思路
针对 V2G 技术这一应用及背景，设计了三种改进后的粒子群算法：第一种是基于模拟退火粒子群算法的电动汽车充放电调度策略，第二种是基于莱维飞行退火粒子群的电动汽车充放电调度策略，第三种方法是基于混沌变异小生镜退火粒子群的电动汽车充放电调度策略，该三种方法均可避免粒子的早熟现象，最后通过实验验证了基于混沌变异小生镜退火粒子群的车辆调度算法拥有更佳的电动汽车调度性能，可降低电动汽车用户成本，对电网负荷削峰填谷的作用显著。
针对传统电动汽车充放电调度成本较高问题，本文在标准粒子群算法基础上设计了一种模拟退火粒子群算法的电动汽车充放电调度策略。引入了分时电价约束，并在目标函数中加入充电成本，解决了用户成本较高问题。仿真实验表明，标准粒子群算法与模拟退火粒子群算法调度策略可降低用户成本更小，且模拟退火粒子群算法更优。
针对传统模拟退火粒子群算法在充放电调度中陷入局部最优，导致电网波动与成本较大问题，本文提出两种调度策略：基于改进粒子群算法、混沌变异小生镜粒子群的电动汽车充放电策略。首先介绍了粒子群算法的衍生算法，分析其特点，并加入单双号限行约束，同时优化了电动汽车多目标模型。将其莱维飞行、模拟退火的思想引入到粒子群算法中。将其混沌变异小生境、模拟退火的思想引入其粒子群算法中。分别设计出改进粒子群算法、混沌变异小生镜粒子群算法，并对其进行充放电调度仿真，仿真实验表明，改进粒子群算法、混沌变异小生镜粒子群调度策略可显著提高电网稳定性和用户参与调度积极性，对比分析，基于混沌变异小生镜退火粒子群算法更为有效，大大降低了电动汽车用户成本与了电网负荷压力。
二、代码与结果
function fv = fitness2(x,~,Varin) res1=0; res2=0; for i=1:6 res1 = res1+Varin.a(i)*x(i)*x(i)+Varin.b(i)*x(i)+Varin.c(i); res2 = res2+Varin.Ea(i)*x(i)*x(i)+Varin.Eb(i)*x(i)+Varin.Ec(i)+Varin.Ed(i)*exp(Varin.Ee(i)*x(i)); end fv(1) = res1; fv(2) = res2; end
if x(i,j)<Varin.pMin(j) if(randi([0,0],1)==0) x(i,j)=Varin.pMin(j); v(i,j)=-v(i,j)*unifit; if(unifit>0.2) unifit = unifit -0.1; end else x(i,j)=Varin.pMin(j)+(Varin.pMax(j)-Varin.pMin(j))*rand; %随机初始化位置 v(i,j)=(Varin.pMax(j)-Varin.pMin(j))*rand*0.5; end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 改进的多目标粒子群算法，包括多个测试函数
% 对程序中的部分参数进行修改将更好地求解某些函数
%% 主函数
function []=main()
Varin = load('mydata.mat');%导入了约束函数的参数
ZDTFV=cell(1,50); %// 创建元胞数组
ZDT=zeros(1,50); %//0矩阵
funcname = 'ZDT1';
times = 10;%相当于独立运行十次程序
M = 100;%MOPSO中的迭代次数
for i=1:times %//循环10次，做以下的迭代
tic; %//计时开始
%[np,nprule,dnp,fv,goals,pbest]=ParticleSwarmOpt(funcname,100,200,2.0,1.0,0.5,M,30,Varin);%--ZDT3 zeros(1,9)-5-》zeros(1,29)
[np,nprule,dnp,fv,goals,pbest]=ParticleSwarmOpt(funcname,100,200,2.0,1.0,0.4,M,10,[0,zeros(1,9)],[1,zeros(1,9)+5],Varin);%--ZDT4
elapsedTime=toc; %//计时结束
ZDTFV(i)={fv};
ZDT(i)=elapsedTime;
display(strcat('迭代次数为',num2str(i)));
end
zdt1fv=cell2mat(ZDTFV');
for i =1:times
display(strcat(i,':'));
disp(ZDT(i));%也就是用时啦
end
disp(zdt1fv);
disp('接下来，更新后的适应度值为：');
zdt1fv=GetLeastFunctionValue(zdt1fv);
disp(zdt1fv);
figure(9)
plot(zdt1fv(:,1),zdt1fv(:,2),'k*');
%以下设置坐标轴的字体形式和大小
xlabel('$f_1$','interpreter','latex','FontSize',25);
ylabel('$f_2$','interpreter','latex','FontSize',25);
set(gca,'FontName','Times New Roman','FontSize',25)%设置坐标轴字体和刻度的大小,get current axes返回当前坐标轴对象的句柄值
%if(strcmp(funcname,'ZDT3'))
axis([0 1 0 1]);
%else
% axis([20555 20790 190.4 190.7]);
%end
% 以下求收敛度
% 先均匀取真实Pareto最优解上的点，再求其两目标函数值
p_true=0:0.002:1;
pf_true1=p_true;
pf_true2=1-p_true.^2;
r=size(zdt1fv,1);
for i=1:r%对每一个非劣解，
for j=1:501
d(i,j)=sqrt((zdt1fv(i,1)-pf_true1(j))^2+(zdt1fv(i,2)-pf_true2(j))^2);
end
end
%下面对d中每行求最小值，即为第i个点与各点的最小距离
for i=1:r
dmin(i)=min(d(i,:));
end
Cmean=mean(dmin);
%Cfangcha=var(dmin)% 各数与均值的差的平方，除以个数-1，这是样本方差
disp('样本均值:');
disp(Cmean);
Cvariance=var(dmin,1);% 各数与均值的差的平方，除以个数,这是数学上方差的定义
disp('样本方差:');
disp(Cvariance)
% 下面求多样性delta
% 先对zdt1fv按第一列升序排序，即按横坐标（第一个目标值）从左到右
zdt1fv=sortrows(zdt1fv);%按第一列升序排序结束，下面计算df(与左边的极值解前沿的距离)和dl（与最右边的极值解前沿的距离）
df=sqrt((zdt1fv(1,1)-0)^2+(zdt1fv(1,2)-1)^2);
r=length(zdt1fv);
dl=sqrt((zdt1fv(r,1)-1)^2+(zdt1fv(r,2)-0)^2);
for i=1:r-1
%第i个和第i+1个解前沿之间的距离为d(i)
c(i)=sqrt((zdt1fv(i,1)-zdt1fv(i+1,1))^2+(zdt1fv(i,2)-zdt1fv(i+1,2))^2);
end
%下面求d的均值
meanNum=mean(c);
%代入公式计算delta的值
%先求和号的部分
sum=0;
for i=1:r-1
sum=sum+abs(c(i)-meanNum);
end
delta=(df+dl+sum)/(df+dl+(r-1)*meanNum);
disp('多样性为：');
disp(delta);
end
%% MOPSO函数定义
%function [np,nprule,dnp,fv,goals,pbest] = ParticleSwarmOpt(funcname,N,Nnp,cmax,cmin,w,M,D,Varin)
function [np,nprule,dnp,fv,goals,pbest] = ParticleSwarmOpt(funcname,N,Nnp,cmax,cmin,w,M,D,lb,ub,Varin)
%待优化的目标函数:fitness
%电厂约束函数：fitness2
%内部种群(粒子数目)：N
%外部种群(非劣解集):Nnp
%适应度参数
%学习因子1:cmax
%学习因子2:cmin
%惯性权重:w
%最大迭代次数：M
%问题的维数：D
%目标函数取最小值时的自变量值:xm
%目标函数的最小值:fv
%迭代次数:cv
%非劣检查:flag
%自适应度参数:unifit:1->0.1
format long;
unifit = 1;
flag = 0;
NP=[];%非劣解集
Dnp=[];%非劣解集距离
params = struct('isfmopso',true,'istargetdis',false,'stopatborder',true);%ZTD2->isfmopso(false->true)改了一下 ZTD3问题时应为true
%x0=lb+(ub-lb).*rand([1,D]);
%T=size(fitness(x0,funcname),2);
T = 2;
goals=zeros(M,N,T);%记下N个粒子M次迭代T维目标变化
% %----初始化种群的个体--------///第1步///
% %x(1,:)=x0;
% %v(1,:)=(ub-lb).*rand([1,D])*0.5;
x = zeros(N,D);
v = zeros(N,D);
% for i=1:N
% for j=1:D
% x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %随机初始化位置
% v(i,j)=(ub(j)-lb(j))*rand*0.5; %随机初始化速度
% end
% end
% %----计算目标向量----------
% %---速度控制
% %vmax=(ub-lb)*0.5;
%vmin = -vmax;
%----初始化种群的个体--------///第1步///
for i=1:N
for j=1:D
x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %随机初始化位置
v(i,j)=(ub(j)-lb(j))*rand*0.5; %随机初始化速度
end
end
%----计算目标向量----------
%---速度控制
vmax=(ub-lb)*0.5;
vmin= -vmax;
%-----求出初始NP-----------第2步///
NP(1,:)=x(1,:);%第一个默认加入（非劣解集），粒子为一行不确定列的，列数表示决策变量数，即问题的维数
NPRule=[0,0,0];%非劣解集参数
Dnp(1,1)=0;
for i=2:N
[NP,NPRule,Dnp,flag] = compare(flag,x(i,:),NP,NPRule,Dnp,Nnp,funcname,params,Varin);
end
%-----初始自身最好位置------///第3步
pbest = x;%自身最优解
%-----在确定每个粒子所对就的目标方格-------//第4步///
%------进入主要循环，按照公式依次迭代------------
for t=1:M
c = cmax - (cmax - cmin)*t/M;
w1=w-(w-0.3)*t/M;
for i=1:N
%-----获得全局最优-------/第5步/
[gbest,NPRule] = GetGlobalBest(NP,NPRule,Dnp);
v(i,:)=w1*v(i,:)+c*rand*(pbest(i,:)-x(i,:))+c*rand*(gbest-x(i,:));
for j=1:D
if v(i,j)>vmax(j)
v(i,j)=vmax(j);
elseif v(i,j)<vmin(j)
v(i,j)=vmin(j);
end
end
x(i,:)=x(i,:)+v(i,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------采取措施，避免粒子飞出空间----------第7步/
%速度位置钳制
if(params.stopatborder)%粒子随机停留在边界
if x(i,1)>ub(1)
x(i,1)=ub(1);
v(i,1)=-v(i,1);
end
if x(i,1)<lb(1)
x(i,1)=lb(1)+(ub(1)-lb(1))*rand; %随机初始化位置
v(i,1)=(ub(1)-lb(1))*rand*0.5;
end
for j=2:D
if x(i,j)>ub(j)
if(randi([0,2],1)==0)%改了0->1
x(i,j)=ub(j);
v(i,j)=-v(i,j);
else
x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %随机初始化位置
v(i,j)=(ub(j)-lb(j))*rand*1.5;
end
end
if x(i,j)<lb(j)
if(randi([0,0],1)==0)
x(i,j)=lb(j);
v(i,j)=-v(i,j)*unifit;
if(unifit>0.2)
unifit = unifit -0.1;
end
else
x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %随机初始化位置
v(i,j)=(ub(j)-lb(j))*rand*0.5;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------每个粒子的目标向量-----------------//第8步///
goals(t,i,:)=fitness(x(i,:),funcname,Varin);
%----------------调整自身---------------------------//第9步/
domiRel = DominateRel(pbest(i,:),x(i,:),funcname,params,Varin);%x,y的支配关系
if domiRel==1%pbest支配新解
continue;
else
if domiRel==-1%新解支配pbest
pbest(i,:) = x(i,:);
elseif(rand*2<1)%新解与pbest互相不支配
pbest(i,:) = x(i,:);
end
%-----------------对NP进行更新和维护-----------------//第10步
[NP,NPRule,Dnp,flag] = compare(flag,x(i,:),NP,NPRule,Dnp,Nnp,funcname,params,Varin);
if flag==1%为克服本算法易陷入局部最优的问题，引入非劣排查监测机制
[NP,flag,x,v] = fresh(NP,flag,x,v);
end
end
end
end
np = NP;%非劣解
nprule=NPRule;
dnp = Dnp;%非劣解之间的距离
r=size(np,1);
fv=zeros(r,T);
for i=1:r
fv(i,:)=fitness(np(i,:),funcname,Varin);
end
end
%%%%%%%%%%%%%%%--------------主函数结束--------------%%%%%%%%%%%%%%%%%
%% 将粒子维护到外部种群
function [np_out,nprule_out,dnp_out,flag] = compare(flag,x,np,nprule,dnp,nnp,funcname,params,Varin)
%np:现有非劣解
%x:需要比较的量
Nnp = nnp;%非劣解集空间
r=size(np,1);%非劣解的个数
np_out=np;%非劣解复本
nprule_out = nprule;
dnp_out = dnp;%非劣解集点之间距离
if r==0
return;
end
for i=r:-1:1
domiRel=DominateRel(x,np(i,:),funcname,params,Varin);
if domiRel==1 %NP(i)被x支配
np_out(i,:)=[];%非劣解剔除该解
nprule_out(i,:)=[];
dnp_out(i,:)=[];
if ~isempty(dnp_out)
dnp_out(:,i)=[];
end
elseif domiRel==-1 %x被NP(i)支配,返回不再比较
return;
end
end
r1=size(np_out,1);%现有非劣解的行列
np_out(r1+1,:)=x;%与所有非支配集粒子比较均占优或不可比较，则NP中加入x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
迭代次数
结果
针对传统模拟退火粒子群算法在充放电调度中陷入局部最优，导致电网波动与成本较大问题，本文在其基础上提出两种调度策略：基于改进粒子群算法与混沌变异小生镜粒子群的电动汽车充放电策略，引入节点功率平衡与单双号限行约束，保证电网稳定性，并将用户放电收入和电网负荷峰谷差作为寻优对象，其中改进粒子群算法为避免传统模拟退火粒子群算法易陷入局部最优问题，简化了粒子群算法速度项，并采用动态惯性权重与莱维飞行策略，而混沌变异与小生镜操作可更好促进种群进化。并将几种算法通过设计调度实验，得出几种混合算法能在相同场景下，对比了几种算法的特点。仿真实验表明，改进粒子群算法、混沌变异小生镜粒子群调度策略均可显著提高电网稳定性和用户参与调度积极性，与传统的模拟退火粒子群算法方法相比，所提方法可以提升原算法的搜索效率，优化了系统负荷与了使用成本。且混沌变异小生镜粒子群更好满足用户侧与配电网侧要求。
691
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
