✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
????个人主页：Matlab科研工作室
????个人信条：格物致知。
更多Matlab仿真内容点击????
⛄ 内容介绍
随着能源市场由传统的垂直一体式结构向交互竞争型结构转变,社区综合能源系统的分布式特征愈发明显,传统的集中优化方法难以揭示多主体间的交互行为.该文提出一种基于主从博弈的社区综合能源系统分布式协同优化运行策略,将综合能源销售商作为领导者,新能源冷热电联供运营商和负荷聚合商作为跟随者,求解各方在追求目标最优时的交互策略.首先,介绍社区综合能源系统的交易模式及数学模型,并将其嵌入到主从博弈框架下,建立一主多从的分布式协同优化模型.其次,证明Stackelberg均衡的唯一性,并通过遗传算法和二次规划相结合的算法求解.最后,通过算例验证所提方法的有效性,供能侧的收益和用能侧的消费者剩余同时得到提升.
本文提出的 CIES 集成 IER、新能源 CCHP 系统以及可调节负荷于一体，以 IER 为纽带，新能源
CCHP 系统为基础，协同互联电力网络，实现经济、高效供能，科学、合理用能，具体架构如图 1 所示。本文中IER是基于电力市场中售电公司的概念提出的，在电能交易的基础上又考虑了热能交易，
满足用户的多样化需求。IER 作为源、荷之间的桥梁，基于供需关系，日前优化购入、售出的电价、
热价，从供能侧购买电、热等能源，并出售给用能侧，从中赚取收益。IER 这一模式的引入，能够提
供相比电网更加灵活的电价策略，对于引导分布式供能系统参与电力市场竞争、鼓励中小型社区用户
科学用能都具有积极作用。在能源交易过程中，IER同样需要承担因价格波动、供需不平衡而带来的风险。当 CCHP 输出电功率无法满足负荷需求时，IER必须高价从电网购电。考虑到热能转换、传输的损耗较大，所以热能主要由社区内的 CCHP 机组就地提供，为了避免出现供热中断，当 IER 不能满足热负荷需求时，需要支付一定的惩罚费用。
新能源CCHP系统将新能源发电与传统燃料发电优势互补，基于能量梯级利用的原则，同时满足
用户电、热、冷不同的能量需求，其结构示意图如图 2 所示。文中新能源包含风电、光伏等，并采用
最大化消纳原则。可控单元包括内燃发电机、燃气锅炉。内燃机发电的同时，缸套水和烟气中携带的
热量可以通过余热装置回收再利用，并与燃气锅炉产生的热量一起，在冬季经热交换器供热，或夏季
经吸收式制冷机转化为冷量为用户供冷。基于 IER的报价，运营商优化各设备的逐时出力，以获得更
高的收益。
⛄ 部分代码
figure
plot(trace(:,1),trace(:,2),'r-*','linewidth',1.5)
xlabel('迭代次数');
ylabel('目标函数值');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('DE上层目标函数曲线');
legend('能源运营商收益曲线')
% set(gca,'YLim',[-20000 20000]);
box off
figure
plot(trace(:,1),trace(:,4),'b-o','linewidth',1.5);
xlabel('迭代次数');
ylabel('目标函数值');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('DE下层目标函数曲线');
legend('用户收益曲线');%%用户收益曲线
% set(gca,'YLim',[0 30000]);
box off
x=1:24;
figure
% yyaxis left
plot(x,Ppv,'r-*','linewidth',1.5)
hold on
plot(x,Pwt,'g-*','linewidth',1.5)
grid
xlabel('时间/h');
ylabel('风光出力');
% yyaxis right
plot(x,de,'b-*',x,dh,'y-+','linewidth',1.5);
xlabel('时间/h');
ylabel('用户负荷曲线/kW');
% set(get(gca,'XLabel'),'Fontsize',14)
% set(get(gca,'YLabel'),'Fontsize',14)
% title('典型用户的负荷预测曲线');
legend('光伏预测曲线','风机预测曲线','电负荷曲线','热负荷曲线');
% set(gca,'XLim',[1 24]);%X轴的数据显示范围
% set(gca,'YLim',[0 2000]);
% legend('光伏预测曲线','风机预测曲线');
box off
xx=1:24;
figure
plot(xx,de,'-*',xx,de1+le,'-+','linewidth',1.5);
hold on
bar((de-de1),0.5)
grid
xlabel('时间/h');
ylabel('电负荷/kW');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('需求响应前后电负荷曲线');
legend('优化前电负荷','优化后电负荷','削减负荷量');
set(gca,'XLim',[0.5 24.]);%X轴的数据显示范围
set(gca,'YLim',[0 2000]);
box off
figure
plot(xx,dh,'-*',xx,dh+lh,'-+','linewidth',1.5);
hold on
grid
bar(lh,0.5)
xlabel('时间/h');
ylabel('热负荷/kW');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('需求响应前后热负荷曲线');
legend('优化前热负荷','优化后热负荷','削减负荷量');
set(gca,'XLim',[1 24]);%X轴的数据显示范围
set(gca,'YLim',[-500 1500]);
box off
figure
bar(Qgb,'stack')
hold on
grid
plot(xx,Pice,'r-*',xx,Pice*(1-n_ice)*n_ex/n_ice,'g-o','linewidth',1.5);
xlabel('时间/h');
ylabel('功率/kW');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('需求响应前后热负荷曲线');
legend('燃气轮机出力','燃气锅炉出力','余热锅炉出力');
set(gca,'XLim',[1 24]);%X轴的数据显示范围
box off
figure
bar(Qgb,'stack')
hold on
grid
plot(xx,Pice,'r-*',xx,Pice*(1-n_ice)*n_ex/n_ice,'g-o','linewidth',1.5);
xlabel('时间/h');
ylabel('功率/kW');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('需求响应前后热负荷曲线');
legend('燃气轮机出力','燃气锅炉出力','余热锅炉出力');
set(gca,'XLim',[1 24]);%X轴的数据显示范围
box off
figure
% plot(xx,Psell,'-*',xx,Pbuy,'-+');
bar(Psell','stack');
hold on
bar(Pbuy','stack')
xlabel('时间/h');
ylabel('功率/kW');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('市场交易结果');
legend('卖的电量','买的电量');
set(gca,'XLim',[0.5 24.5]);%X轴的数据显示范围
box off
figure
for xx=1:24
a1(xx)=best(xx);%ceb
a3(xx)=best(xx+48);%ces
end
xx=1:24;
stairs(xx,grid_fs,'r--','linewidth',1.5);
hold on
stairs(xx,grid_sw,'g--','linewidth',1.5);
hold on
stairs(xx,a1,'b-*','linewidth',1.5);
hold on
stairs(xx,a3,'c-+','linewidth',1.5);
xlabel('时间/h');
ylabel('电价/元');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('内部电价结果');
legend('电网分时电价','电网上网电价','IER购电电价','IER售电电价');
set(gca,'XLim',[1 24]);%X轴的数据显示范围
set(gca,'YLim',[0 2]);
box off
figure
for xx=1:24
a2(xx)=best(xx+24);%ceb
a4(xx)=best(xx+72);%ces
end
xx=1:24;
stairs(xx,ch_min,'r--','linewidth',1.2);
hold on
stairs(xx,ch_max,'g--','linewidth',1.2);
hold on
stairs(xx,a2,'b-*','linewidth',1.2);
hold on
stairs(xx,a4,'c-+','linewidth',1.2);
xlabel('时间/h');
ylabel('电价/元');
set(get(gca,'XLabel'),'Fontsize',14)
set(get(gca,'YLabel'),'Fontsize',14)
title('内部热价结果');
legend('热价下线','热价上限','IER购热价','IER售热价');
set(gca,'XLim',[1 24]);%X轴的数据显示范围
set(gca,'YLim',[0 1]);
box off
x=1:24;
PP=[Pbuy;Pdischarge;Pice;Pwt;Ppv];
PP1=[Psell;-Pcharge];
figure
bar(PP','stack');
% set(h,'Orientation','horizon')
hold on
bar(PP1','stack');
plot(x,value(de1+le),'g-*','linewidth',1.5);
legend('交换功率','蓄电池出力','燃气轮机出力','风电出力出力','光伏出力','售电电量','充电电量','负荷曲线');
xlabel('时段');ylabel('功率/kW');
hold off
x=1:24;
QQ=[Hto;Qgb*n_c;Pice*(1-n_ice)*n_ex/n_ice*n_c];
QQ1=[-Hti];
figure
bar(QQ',0.6,'stack');
hold on
bar(-value(dh+lh)',0.6,'stack');
bar(QQ1',0.6,'stack');
hold on
xlabel('时段');ylabel('功率/kW');
legend('蓄热槽出力','燃气锅炉出力','余热锅炉出力','总热负荷','充热量');
% set(h,'Orientation','horizon')
hold off
⛄ 运行结果
⛄ 参考文献
[1]王海洋, 李珂, 张承慧,等. 基于主从博弈的社区综合能源系统分布式协同优化运行策略[J]. 中国电机工程学报, 2020, 40(17):10.