考虑源荷随机特征的热电联供微网优化（含matlab程序）

一、前言

二、含可再生能源的CHP型微网系统

三、CCP理论

四、具体模型

五、不含随机变量分析的matlab程序设计

1.粒子群寻优功能代码段

2.目标函数子程序

3.其他代码段

六、基于CCP的粒子群优化程序

1.含随机变量的约束条件处理

2.随机变量生成

3.置信水平检验部分

七、程序结果

1.不含随机变量的程序结果

2.含随机变量处理的程序结果

八、程序下载链接

一、前言

在常规的电热微网优化中，可以得到蓄电池、外网交互、燃料电池、余热锅炉等设备的功率，但是常规模型已经比较基础了，在发文章或者做毕业论文的时候用太初级，工作量和深度难以满足要求，因此可以考虑不确定变量的鲁棒性、考虑机会约束等，本文用王锐的《含可再生能源的热电联供型微网经济运行优化》为例，分析机会约束规划理论建立模型和编程方法，同时给出常规粒子群PSO和基于CCP理论的粒子群算法上面的区别。

二、含可再生能源的CHP型微网系统

三、CCP理论

四、具体模型

具体模型不在此赘述了，可以详见文献资料。

五、不含随机变量分析的matlab程序设计

1.粒子群寻优功能代码段

function [ bestPosition, fitValue ] = ...
PSOFUN( CostFun,nVar,VarMin,VarMax,MaxIt,nPop )
%% PSO Parameters 
CostFunction=@(x) CostFun(x);        % Cost Function
w=1;            % Inertia Weight
wdamp=0.99;     % Inertia Weight Damping Ratio
c1=1.5;         % Personal Learning Coefficient
c2=2.0;         % Global Learning Coefficient
VarSize=[1 nVar];   % Size of Decision Variables Matrix
% Velocity Limits
VelMax=0.1*(VarMax-VarMin);
VelMin=-VelMax;%% Initializationempty_particle.Position=[];
empty_particle.Cost=[];
empty_particle.Velocity=[];
empty_particle.Best.Position=[];
empty_particle.Best.Cost=[];particle=repmat(empty_particle,nPop,1);GlobalBest.Cost=inf;for i=1:nPop% Initialize Positionparticle(i).Position=unifrnd(VarMin,VarMax,VarSize);%随机初始化变量% Initialize Velocityparticle(i).Velocity=zeros(VarSize);%速度% Evaluationparticle(i).Cost=CostFunction(particle(i).Position);%目标% Update Personal Bestparticle(i).Best.Position=particle(i).Position;%更新个体最优particle(i).Best.Cost=particle(i).Cost;% Update Global Best%更新全局最优if particle(i).Best.CostendendBestCost=zeros(MaxIt,1);%% PSO Main Loopfor it=1:MaxItfor i=1:nPop% Update Velocity更新速度particle(i).Velocity = w*particle(i).Velocity ...+c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...+c2*rand(VarSize).*(GlobalBest.Position-particle(i).Position);% Apply Velocity Limits速度约束particle(i).Velocity = max(particle(i).Velocity,VelMin);particle(i).Velocity = min(particle(i).Velocity,VelMax);% Update Position更新位置particle(i).Position = particle(i).Position + particle(i).Velocity;% Velocity Mirror Effect变量越限镜像处理IsOutside=(particle(i).Position| particle(i).Position>VarMax);particle(i).Velocity(IsOutside)=-particle(i).Velocity(IsOutside);% Apply Position Limits%变量越限处理particle(i).Position = max(particle(i).Position,VarMin);particle(i).Position = min(particle(i).Position,VarMax);% Evaluation新目标particle(i).Cost = CostFunction(particle(i).Position);% Update Personal Bestif particle(i).Costif particle(i).Best.CostendendendBestCost(it)=GlobalBest.Cost;disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);w=w*wdamp;endbestPosition = GlobalBest.Position;%记录最优位置
fitValue = GlobalBest.Cost;%记录最优目标结果end

粒子群代码部分主要是实现初始粒子设置和迭代寻优功能，采用结构型变量实现简约、利于理解的代码功能，且注释比较清晰，方便学习！

2.目标函数子程序

function fun = fun_objective(x)
%% 准备工作 
parameter; %输入所有的数据 
% 各个决策变量的含义 
Pfl = x(1:24);      % 燃料电池出力 
Pbt = x(25:48);  % 蓄电池出力 
Pex = x(49:72); % 交互功率
Pgb = x(73:96); % 锅炉出力
fun =0; 
%% 书写目标函数
for t=1:24fun = fun + 1/2*(Cph+Cse)*Pex(t) + 1/2*(Cph-Cse)*abs( Pex(t) ) ...+ Cgas*(Pfl(t)/eta_fl+ Pgb(t)/eta_gb ) + Pfl(t)*Cfl_om + ...Pfl(t)*r_fl*eta_hrbl*Cbl_om + abs(Pbt(t))*Cbt_om+ ...Pgb(t)*Cgb_om + Pwt(t)*Cwt_om + Ppv(t)*Cpv_om;
end
%% 书写约束
% ******************* 等式约束****************************
h=[];
for t=1:24 % (1) 电能平衡约束if Pbt(t)<=0h = [h, Pex(t)+Pfl(t)+Pwt(t)+Ppv(t)+Pbt(t)/eta_ch-Pel(t) ]; %=0elseh = [h, Pex(t)+Pfl(t)+Pwt(t)+Ppv(t)+Pbt(t)*eta_dis-Pel(t) ]; %=0end
end
for t=1:24 % (2) 热能平衡约束h = [h, Pgb(t)+Pfl(t)*r_fl*eta_hrbl-Pth(t) ]; %=0
end% (3) 电池储能初始和最终状态相等约束h = [h, sum(Pbt) ]; %=0
% ******************* 不等式约束 ***************************
g=[];
for t=2:24 % (1) 燃料电池爬坡约束g=[g, Pfl(t)-Pfl(t-1)-deltaP_up] ; % <=0    g=[g, -( Pfl(t)-Pfl(t-1)-deltaP_down ) ] ; % <=0    
end
for t=1:24 % (2) 余热锅炉约束g=[g, Pfl(t)*r_fl*eta_hrbl-Pbl_max ] ; % <=0    g=[g, - ( Pfl(t)*r_fl*eta_hrbl-Pbl_min ) ] ; % <=0    
end
for t=1:24 % (3) 蓄电池约束g=[g, Wbt_init-sum(Pbt(1:t))-Wbt_max  ] ; % <=0    g=[g, -( Wbt_init-sum(Pbt(1:t))-Wbt_min )  ] ; % <=0    
end
%**********************罚函数处理*************************
Big=100000;
small=0.01;
N=length(g);
M=length(h);
G=0;
for n=1:NG=G+max(0, g(n))^2;
end
H=0;
for m=1:MH=H+max(  0, abs(h(m))-small  )^2;
end
%*******************加入罚函数后的目标函数******************
fun=fun+Big*(H+G);
end

目标函数代码完美复刻了文献中的目标函数和约束条件，约束部分采用清晰简明的等式和不等式部分，方便理解，采用罚函数的形式形成最终目标函数值。

3.其他代码段

其他主要就是主函数调用和参数定义代码段，由于功能性单一，在此不再列出，可以通过评论区链接进行购买下载。

六、基于CCP的粒子群优化程序

1.含随机变量的约束条件处理

%等式约束随机变量处理
Pex = Pel - Pfl - max(Pbt,0)*eta_dis - min(Pbt,0)/eta_ch - Pwt - Ppv;
Pgb = Pth - Pfl*r_fl*eta_hrbl

2.随机变量生成

考虑源荷的随机特征，具体参数如下：

for i=1:24%风电、光伏服从预测值为均值，0.1为方差的均匀分布Pwtmax(i)=Pwt(i)+1.2^0.5/2;Pwtmin(i)=Pwt(i)-1.2^0.5/2;Ppvmax(i)=Ppv(i)+1.2^0.5/2;Ppvmin(i)=Ppv(i)-1.2^0.5/2;%热电负荷服从预测值偏差+-10%的均匀分布Pelmax(i)=1.1*Pel(i);Pelmin(i)=0.9*Pel(i);Pthmax(i)=1.1*Pth(i);Pthmin(i)=0.9*Pth(i);%形成新的随机变量Pwt1(i)=Pwtmin(i)+(Pwtmax(i)-Pwtmin(i))*rand;Ppv1(i)=Ppvmin(i)+(Ppvmax(i)-Ppvmin(i))*rand;Pel1(i)=Pelmin(i)+(Pelmax(i)-Pelmin(i))*rand;Pth1(i)=Pthmin(i)+(Pthmax(i)-Pthmin(i))*rand;
end