一、前言
在常规的电热微网优化中，可以得到蓄电池、外网交互、燃料电池、余热锅炉等设备的功率，但是常规模型已经比较基础了，在发文章或者做毕业论文的时候用太初级，工作量和深度难以满足要求，因此可以考虑不确定变量的鲁棒性、考虑机会约束等，本文用王锐的《含可再生能源的热电联供型微网经济运行优化》为例，分析机会约束规划理论建立模型和编程方法，同时给出常规粒子群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; %% Initialization empty_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 Position particle(i).Position=unifrnd(VarMin,VarMax,VarSize);%随机初始化变量 % Initialize Velocity particle(i).Velocity=zeros(VarSize);%速度 % Evaluation particle(i).Cost=CostFunction(particle(i).Position);%目标 % Update Personal Best particle(i).Best.Position=particle(i).Position;%更新个体最优 particle(i).Best.Cost=particle(i).Cost; % Update Global Best%更新全局最优 if particle(i).Best.Cost<GlobalBest.Cost GlobalBest=particle(i).Best; end end BestCost=zeros(MaxIt,1); %% PSO Main Loop for it=1:MaxIt for 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<VarMin | 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 Best if particle(i).Cost<particle(i).Best.Cost particle(i).Best.Position=particle(i).Position; particle(i).Best.Cost=particle(i).Cost; % Update Global Best更新全局最优 if particle(i).Best.Cost<GlobalBest.Cost GlobalBest=particle(i).Best; end end end BestCost(it)=GlobalBest.Cost; disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]); w=w*wdamp; end bestPosition = 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:24 fun = 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)<=0 h = [h, Pex(t)+Pfl(t)+Pwt(t)+Ppv(t)+Pbt(t)/eta_ch-Pel(t) ]; %=0 else h = [h, Pex(t)+Pfl(t)+Pwt(t)+Ppv(t)+Pbt(t)*eta_dis-Pel(t) ]; %=0 end 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:N G=G+max(0, g(n))^2; end H=0; for m=1:M H=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
3.置信水平检部分
该部分程序可详见评论区程序链接。
七、程序结果
1.不含随机变量的程序结果
2.含随机变量处理的程序结果
通过结果能够看出，该方法实现了所有约束，优化效果均较好。
八、程序链接
程序包括两部分，分别是不含随机变量和包含随机变量处理两个程序代码，可以详细分析两种方法的对比结果，需要的可私信，下面是程序讲解视频！
考虑源荷不确定性的微网经济运行优化matlab
