目录
1 主要内容
模型架构图
目标函数
非合作博弈流程
2 部分代码
3 程序结果
4 下载链接
1 主要内容
该程序复现《基于非合作博弈的风-光-氢微电网容量优化配置》，程序包含3种场景，场景1中包含风电、光伏和制氢-储氢-发电3种分布式电源，而场景2中仅包含风电和光伏，场景3中包含风电和制氢-储氢-发电。（通过设置相应的变量为0来实现对比场景）
以风电、光伏和制氢-储氢-发电系统 3 个投资方作为博弈主体，并以各投资方收益最大化为优化目标，建立基于非合作博弈的风-光-氢微电网容量优化配置模型。考虑各博弈参与者的投资成本、运维成本、购售电成本、弃风弃光惩罚费用和负荷中断惩罚费用等经济因素，利用粒子群算法对各博弈参与者的容量配置进行单独优化，确定各博弈参与者收益最大化的 Nash 均衡点。程序采用matlab编写，无需其他软件包，注释清楚，方便学习！
模型架构图
-
目标函数
1.投资成本
2.运维成本
3.购电成本
4.售电收益
5.弃风弃光成本
6.负荷中断补偿成本
7.补贴电价
%% 1.投资费用
f=(r*(1+r)^m)/(((1+r)^m)-1);
C_wt_inv=num_WT*P_wt*e_wt*f*(1/12); %风机的月投资费用
C_pv_inv=num_PV*P_pv*e_pv*f*(1/12); %光伏的月投资费用
C_h_inv=(num_FC*P_fc*e_fc*f+num_HST*V_hst*e_hst*f+num_EC*W_EC*e_ec*f)*(1/12); %制氢-储氢-发电系统的月投资费用
%% 2.运维成本
y_wt=0.01; %风电运维成本占总投资成本的比例
C_wt_run=y_wt*C_wt_inv; %风机的月运行维护成本
y_pv=0.001; %光伏运维成本占总投资成本的比例
C_pv_run=y_pv*C_pv_inv; %光伏的月运行维护成本
y_h=0.01;%制氢-储氢-发电系统运维成本占总投资成本的比例
C_h_run=y_h*C_h_inv;%制氢-储氢-发电系统的月运行维护成本
%% 3.购售电费用、弃风弃光费用、负荷中断补偿费用
p=Load-(num_WT*P_wind+num_PV*pv_end);
E_hst=zeros(1,721);
E_hst(1)=E_hst_init;
P_FC=zeros(1,720);
P_EC=zeros(1,720);
for i=1:720
if p(i)>0
if p(i)<(num_FC*P_FC_max) %判断p是否小于燃料电池的总功率
P_FC(i)=p(i);
elseif p(i)>(num_FC*P_FC_max)
P_FC(i)=(num_FC*P_FC_max);
end
E_hst(i+1)=E_hst(i)-(P_FC(i)*W_g2p)/(a*eff_g2p);
if E_hst(i+1)<=0 %判断储氢罐的储氢量是否被消耗完
P_FC(i)=0;
E_hst(i+1)=E_hst(i);%恢复到上一次的状态
end
elseif p(i)<0
if abs(p(i))<(num_EC*P_EC_max) %判断p是否小于电解槽的总功率
P_EC(i)=abs(p(i));
elseif abs(p(i))>(num_EC*P_EC_max)
P_EC(i)=(num_EC*P_EC_max);
end
E_hst(i+1)=E_hst(i)+P_EC(i)*a*eff_p2g*W_p2g;
if E_hst(i+1)>=E_hst_all %判断储氢罐的储氢量是否超过最大容量
P_EC(i)=0;
E_hst(i+1)=E_hst(i);%恢复到上一次的状态
end
end
end
P_H=P_FC-P_EC; %求取氢储系统的实际交换功率
P_EL=Load-(num_WT*P_wind+num_PV*pv_end+P_FC);%求取等效负荷
P_buy=zeros(1,720);
P_waste_wt=zeros(1,720);
P_waste_pv=zeros(1,720);
P_break=zeros(1,720);
P_exc=zeros(1,720);
P_wt_se=zeros(1,720);
P_pv_se=zeros(1,720);
P_max=zeros(1,720);
P_s=zeros(1,720);
for i=1:720
P_max(i)=P_EC(i)+Load(i)+P_grid;
P_exc(i)=num_WT*P_wind(i)+num_PV*pv_end(i)-P_max(i);
if P_EL(i)>0 %购电
if P_EL(i)< P_grid
P_buy(i)=P_EL(i);
elseif P_EL(i)> P_grid
P_buy(i)=P_grid;
P_break(i)=P_EL(i)-P_grid; %需要中断的负荷
end
elseif P_EL(i)<0 %售电
if P_exc(i)<=0
P_wt_se(i)=num_WT*P_wind(i);
P_pv_se(i)=num_PV*pv_end(i);
if (P_wt_se(i)+P_pv_se(i))>(P_EC(i)+Load(i))
P_s(i)=(P_wt_se(i)+P_pv_se(i))-(P_EC(i)+Load(i));
end
elseif P_exc(i)>0
P_wt_se(i)= (num_WT*P_wind(i)*P_max(i))/(num_WT*P_wind(i)+num_PV*pv_end(i));
P_pv_se(i)= (num_PV*pv_end(i)*P_max(i))/(num_WT*P_wind(i)+num_PV*pv_end(i));
P_waste_wt(i)= num_WT*P_wind(i)-P_wt_se(i);
P_waste_pv(i)= num_PV*pv_end(i)-P_pv_se(i);
P_s(i)=P_grid;
end
end
end
p_waste=zeros(1,720);
for i=1:720
if P_exc(i)>=0
p_waste(i)=P_exc(i);
end
end
% 1)购电成本
C_buy_wt=(c_buy*sum(P_buy)*num_WT*P_wt)/(num_WT*P_wt+num_PV*P_wt+num_FC*P_fc);%风电博弈参与者的购电成本
C_buy_pv=(c_buy*sum(P_buy)*num_PV*P_pv)/(num_WT*P_wt+num_PV*P_wt+num_FC*P_fc);%光伏博弈参与者的购电成本
C_buy_fc=(c_buy*sum(P_buy)*num_FC*P_fc)/(num_WT*P_wt+num_PV*P_wt+num_FC*P_fc);%燃料电池博弈参与者的购电成本
% 2)售电收益
C_sell_wt=c_buy*sum(P_wt_se);%风电博弈参与者的售电收益
C_sell_pv=c_buy*sum(P_pv_se);%光伏博弈参与者的售电收益
C_sell_fc=c_buy*sum(P_FC)+0.5*sum(P_EC)*eff_p2g*W_p2g*c_se_o2;%燃料电池博弈参与者的售电收益+卖氧气
% 3)弃风弃光成本
C_waste_pv=c_waste*sum(P_waste_pv);%光伏博弈参与者的弃风弃光费用
C_waste_wt=c_waste*sum(P_waste_wt);%风电博弈参与者的弃风弃光费用
% 4)负荷中断补偿成本
C_break_wt=(c_break*sum(P_break)*num_WT*P_wt)/(num_WT*P_wt+num_PV*P_wt+num_FC*P_fc);%风电博弈参与者的负荷中断补偿费用
C_break_pv=(c_break*sum(P_break)*num_PV*P_pv)/(num_WT*P_wt+num_PV*P_wt+num_FC*P_fc);%光伏博弈参与者的负荷中断补偿费用
C_break_fc=(c_break*sum(P_break)*num_FC*P_fc)/(num_WT*P_wt+num_PV*P_wt+num_FC*P_fc);%燃料电池参与者的负荷中断补偿费用
% 5)补贴电价
C_subsidy_wt=c_sub_wt*num_WT*sum(P_wind);
C_subsidy_pv=c_sub_pv*num_PV*sum(pv_end);
C_subsidy_h=c_sub_h*sum(P_EC);
-
非合作博弈流程
由此可以看出，非合作博弈模型并不复杂，寻找纳什均衡点的过程也很清晰，大家在做模型的过程中也可以通过设置策略提高深度，该文献就是一个很好的参考。
2 部分代码
%% %设定均衡点初值 PV0=0;%光伏组件的初始个数 WT0=7;%风机的初始台数 EC0=5;%电解槽的初始个数 HST0=6;%储氢罐的初始个数 FC0=9;%燃料电池的初始个数 %% %各博弈参与者或联盟依次进行独立优化决策 %此处选用粒子群算法作为优化算法 %此处以多目标优化问题为例，从博弈的角度来求解收益的最大值Uwt、Upv、Uh %参数初始化 c1=3; %个体学习因子 c2=3; %社会学习因子 w_max=0.9; %设置最大惯性权重为0.9 w_min=0.4; %设置最小惯性权重为0.4 max_die_dai=30; %迭代次数设置为100 size_zhong_qun=50; %粒子种群规模设置为100 PV_N=randi(2650,size_zhong_qun,1); %初始化粒子PV的位置 WT_N=randi(22,size_zhong_qun,1); %初始化粒子WT的位置 EC_N=randi(16,size_zhong_qun,1); %初始化粒子EC的位置 HST_N=randi(20,size_zhong_qun,1); %初始化粒子HST的位置 FC_N=randi(110,size_zhong_qun,1); %初始化粒子FC的位置 v_PV=3.*rands(size_zhong_qun,1); %初始化粒子PV的飞翔速度 v_WT=3.*rands(size_zhong_qun,1); %初始化粒子WT的飞翔速度 v_h=3.*rands(size_zhong_qun,3); %初始化粒子EC的飞翔速度 %定义适应度函数 for i=1:size_zhong_qun f_wt(i,:)=fitness(WT_N(i,:),PV0,EC0,HST0,FC0,1); %序号1代表fitness函数当前返回wt的收益 end for i=1:size_zhong_qun f_pv(i,:)=fitness(WT0,PV_N(i,:),EC0,HST0,FC0,2); %序号1代表fitness函数当前返回pv的收益 end for i=1:size_zhong_qun f_h(i,:)=fitness(WT0,PV0,EC_N(i,:),HST_N(i,:),FC_N(i,:),3); %序号1代表fitness函数当前返回h的收益 end %最大速度值过大，容易导致越过最优区域，值过小，容易陷入局部最优 v_PV_max=3; %速度最大值为3 v_PV_min=-3; %速度最小值为-3 v_WT_max=3; %速度最大值为3 v_WT_min=-3; %速度最小值为-3 v_H_max=3; %速度最大值为3 v_H_min=-3; %速度最小值为-3 PV_N_max=5000; %种群PV中个体位置的最大值 PV_N_min=10; %种群PV中个体位置的最小值 WT_N_max=50; %种群WT中个体位置的最大值 WT_N_min=1; %种群WT中个体位置的最小值 EC_N_max=30; %种群EC中个体位置的最大值 EC_N_min=2; %种群EC中个体位置的最小值 HST_N_max=40; %种群HST中个体位置的最大值 HST_N_min=1; %种群HST中个体位置的最小值 FC_N_max=200; %种群FC中个体位置的最大值 FC_N_min=4; %种群FC中个体位置的最小值 %群体中个体的最佳位置随着迭代的进行而不断发生变化 personal_best_local_x=WT_N; %种群wt中每个个体的最佳位置，随着迭代的进行，数值将会发生改变 personal_best_value_x=f_wt; %种群wt中每个个体的最佳适应度值 [global_best_value_x,i]=max(personal_best_value_x); %计算出当前的全局最佳值和最佳值所对应的位置 global_best_local_x=personal_best_local_x(i,:); %当前的全局最佳位置
3 程序结果
