前言
随着现代城市面临环境问题，原来燃煤的水和空间供暖设备已逐渐被电锅炉和热泵等电气设备所取代。此外，集中生产热能并通过管网分配热能的区域供暖系统，由于其更高的效率，在冬季漫长寒冷的国家和地区越来越受欢迎。供暖设备的电气化和区域供暖网络（DHN）的迅速发展，创造了综合能源分配系统，利用城市地区的多种能源。这种耦合基础设施的关键组成部分是能源枢纽（EH），它起着能源生产、转换和储存的作用。EH影响具有异构能源载体的不同物理系统的运行，并有可能同时参与电力、供暖和天然气市场。
市场出清模型
所谓出清，即确定市场中的零售商可以出售的电量和价格。由此扩展到电热互联综合能源系统中，该系统共有3个责任主体：电网PDN，热网DHN，系统运行商EH，电网（热网）将可出售的电量及价格（热量及价格）上报给EH，EH与其他网络进行招投标，以高价匹配高价的形式锁定交易方，从而降低本系统的成本。
通过利用电能存储单元（ESU）和储热单元（TSU）提供的灵活性以及电力、供暖和天然气市场之间的耦合，专注于EH的投标与定价。作出以下假设：
1）考虑未来一天的电力和供暖市场变化，对负荷进行预测。忽略模型的不确定性，在需要时使用场景方法进行补充。实时市场中平衡发电和需求之间的差异。市场由电力市场运营商（EMO）和供热市场运营商单独运营。
2）PDN和DHN通过EH连接。PDN中的潮流状态由线性支路流模型表示。DHN中的热状态通过具有恒定流量的液压和热组合公式来描述。联合市场将其电力（供暖）报价/投标价格和容量提交给EMO（HMO），EMO随后对电力（供热）市场出清，以最大限度地降低总生产成本，并根据报价进行付款，这两个市场的投标策略在空间和时间上是对应的。
3）电力和热力的价格受到约束，并不是价格越高越好，因为较较高的价格会导致较低的市场份额（EMO和HMO自发自用，而不是购买能源）。从电力市场的角度来看EH是一个生产消费者：它既可以消耗电力，也可以提供电力；从供暖市场的角度来看，它是一家生产商。
在PDN的决策问题中，不确定性通常源于市场价格和可再生能源发电。两个市场之间的能源价格取决于投标策略或双边协议，天然气价格由天然气零售商决定。鉴于目前的天然气市场，天然气价格保持不变。然而，PDN总线上的节点电价由上层电网决定。
由于储能的存在，在模型中不考虑不确定因素，可以减少可再生能源波动带来的影响。尽管如此，如果考虑不确定因素对经济性的影响，我们可以使用基于场景的随机规划方法，从而最大限度地减少EH的预期收益。
所提出的模型和方法可以在各种应用中提供有用的信息。能源枢纽所有者可以使用所提出的方法来确定电力和供暖市场的最优投标策略；投资者可以使用该模型来检验给定条件下能源枢纽的利润，帮助他选择更好的计划；政府机构可以使用该模型来调查能源枢纽的市场力量和综合能源系统中的均衡状态，从而确定关键的市场参数（如最高发行价格），维护市场公平。
系统结构如下：
程序介绍
电-热（EH）互联系统提高了系统的能源利用效率，其中能源母线在能源生产、转换和储存方面发挥着重要作用。在综合能源背景下，程序提出了考虑平衡约束数学（MPEC）模型，研究EH在电力和热力市场中如何保证运行的经济性。在上层模型中，EH向电力市场和热力市场提交价格和容量。在下层模型中，两个市场同时出清，EH和两个市场之间的电热交易被锁定。通过互补性和松弛性条件进行整数处理，并对双线性项进行二进制展开，得到混合整数线性规划的MPEC模型。程序中算例丰富，注释清晰，干货满满，创新性和可扩展性很高，足以撑起一篇高水平论文！下面对程序做简要介绍！
程序适用平台：Matlab+Yalmip+Cplex
参考文献：《Participation of an Energy Hub in Electricity and Heat Distribution Markets: An MPEC Approach》-IEEE Trans on Smart Grid；
程序结果
部分程序
%% 参数选择
BAU = 0; % 基准值 1, 0
Option_Price = 4; % 分时电价 1, 峰谷 2, 偶尔 3, LMP 4
Option_Gas = 1; % 便宜 1, 昂贵 2, 随机 3
Option_GenCostCap = 1; % 减少发电成本 & 增加 GTGB 容量
case 1 % 春季SeasonP = 1; SeasonD = 2;
case 2 % 夏季SeasonP = 2; SeasonD = 1;
case 3 % 秋季SeasonP = 3; SeasonD = 3;
case 4 % 冬季SeasonP = 4; SeasonD = 3;
%% 负荷参数
eval('SysPDN33'); % 负荷 0.35 ~0.5 p.u.eval('SysDHN32'); % 负荷 0.2 p.u.
%% 约束条件
% 上层约束为常数
Para_Price_P_in = 7758;Para_Price_P_out = 201314;
Para_Price_H = 52121;Para_Quan_P_in = 7474;
%% KKT条件
%%% 导出矩阵形式%%% 定义上层变量
clear Pr_E_in Pr_E_out p_hat_hub_in p_hat_hub_out
Pr_E_in = sdpvar(1,NT,'full');Pr_E_out = sdpvar(1,NT,'full');
p_hat_hub_in = sdpvar(1,NT,'full');p_hat_hub_out = sdpvar(1,NT,'full');
%%% 查找参数
Ind_Vara_P_out = getvariables(p_hub_out); % 无功约束
%%% 更新变量矩阵
x_kkt_PDN = sdpvar(length(recoverymodel_PDN.used_variables),1,'full');
big_M = 1e5;
%%% 变量 & 约束
Ind_Vara_Others = setdiff(1:size(x_kkt_PDN,1),[Ind_Vara_P_out, Ind_Vara_P_in]);
Ind_Cons_Others = setdiff(1:size(model_PDN.Aineq,1),[Ind_Cons_Quan_P_out, Ind_Cons_Quan_P_in]);
%% KKT约束
Obj_DOWN_DHN = (1+GenCost_Factor)*alphaGB(DataDHN.IndGB(1),:)*sum(Hg_GB(DataDHN.IndGB(1),:).^2) + betaGB(DataDHN.IndGB(1),:)*sum(Hg_GB(DataDHN.IndGB(1),:))+...
%% 上层变量
clear Pr_H h_hat_hubPr_H = sdpvar(1,NT,'full');h_hat_hub = sdpvar(1,NT,'full');
%%% 查找参数
Ind_Vara_H = getvariables(h_hub); % 无功变量
% !!!注意两次导出时变量的重叠!!!
[Index,ia_out,ib] = intersect(recoverymodel_DHN.used_variables,Ind_Vara_H);
%%% 更新矩阵
x_kkt_DHN = sdpvar(length(recoverymodel_DHN.used_variables),1,'full');
big_M = 1e5;
%% 市场博弈层F_PF_Hub_Bid = [];部分内容源自网络，侵权联系删除！
欢迎感兴趣的小伙伴关注并私信获取完整版代码，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！
