【EI复现】多区域综合能源系统热网建模及系统运行优化（Matlab代码实现）-平芜编程栈

💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。

💥1 概述

文献来源：

目前，关于 IES 的研究多集中于 CCHP 系统建模、规划与经济调度方面[16-21]。文献[1]对 CCHP 微网建模、规划与能量管理的研究现状进行了综述，总结了当前研究中存在的不足；文献[2]提出了

CCHP 微网母线式结构，给出了 CCHP 微网的一般结构，利用该结构建立了 CCHP 微网的日前经济调度的混合整数线性规划 (mixed integer linearprogramming，MILP)模型；文献[22]采用矩阵形式，构建了 CCHP 系统的输入输出模型，并建立了系统评价模型，确定了系统热电联供机组的最优容量；文献[9]以德黑兰东部某商住综合区为例，建立了含有 DHCN 的区域 CCHP 系统规划模型，研究了系统设备的最优配置问题，但文中的 DHCN 模型过于简单，缺乏相应的理论依据；文献[23]根据热力学定律，推导了能量在线路中的普遍化传递方程，并基于广义基尔霍夫定律建立了能量网络方程，为能量网络的建模与分析奠定了理论基础，但该模型过于抽象，在工程应用中存在一定困难。针对上述问题，本文建立了区域热网能量传输模型和含有热网的多区域 IES 运行优化模型。首先，基于传热学基本原理与管网基本理论建立热网能量传输的通用模型，该模型以热网热媒流量、温度作为优化变量，能够较为精确的描述热网状态，但由于含有非线性项，计算量大，求解较为困难。其次，根据热网通用模型推导出热网热损方程，将其线性化得到热网能量流模型，模型仅含有热网传输热功率变量，且为混合整数线性模型，给系统优化求解带来极大方便。再次，推导热网流量–温度方程，该方程可以通过热网功率流求解网络状态参数(热媒流量、温度)。最后，将 CCHP 运行优化模型与热网模型相结合，建立含有热网的多区域 IES混合整数线性规划模型。仿真结果表明，多区域 IES协同运行能够带来可观的经济效益，同时验证本文建立的热网模型的有效性。

本文研究的多区域 IES 结构如图 4 所示，该系统包含若干个不同区域的 CCHP 系统、热网、电网、水网和燃气网。当允许向电网售电时，CCHP 系统可将多余的电能出售给电网。各个 CCHP 系统通过热网进行热功率交换，热能供应过剩的区域向热网注入热能，热能供应不足的区域从热网得到热能。由于 CCHP 系统与热网热能交互具有双向性，从而其耦合环节应具有方向选择性。如图 5 所示。

算例：

这项研究工作复现了《多区域综合能源系统热网建模及系统运行优化》模型，并对该模型进行了进一步扩展和优化。通过基于传热学基本原理的建模，成功建立了区域热网能量传输通用模型，实现了对热网热损方程的线性化处理，从而实现了热网能量流的建模。在已有冷热电联供系统运行优化模型的基础上，新建了一个包含热网的多区域综合能源系统（IES）优化混合整数线性规划模型。该模型以最小化电网购电费用、最大化向电网售电所得费用、最小化燃气费、弃光成本和热网运行费用之和为优化目标，同时考虑了CCHP系统、热网和各种设备之间的耦合关系约束条件。

为了更好地模拟多样化的可再生能源出力情景，研究采用了蒙特卡洛和kmeans聚类技术生成了可再生能源出力的随机优化场景。此外，研究还对原先的确定性模型进行了调整，转变为日前确定性和实时鲁棒优化模型，以更好地适应实际运行中的波动情况。

该程序采用了MATLAB编程结合Gurobi求解器进行实现。最近对程序进行了修正和优化，为对综合能源系统方向感兴趣的同学提供了学习和参考的机会。欢迎有兴趣的同学下载学习，以深入了解多区域热网系统建模与优化的相关技朮。

这次针对程序的修正和优化工作涵盖了三个主要方面，具体如下所述：

1. 针对程序的求解速度进行了优化，通过对程序代码的精心调整和优化，成功缩短了求解时间，目前程序的求解时间仅需大约5分钟左右，大大提升了程序的效率和实用性。

2. 在这次修正中，成功解决了热功率失衡的问题，这对于热网系统的稳定运行至关重要。通过调整算法和参数，有效地处理了热功率失衡情况，提升了系统的整体性能和可靠性，确保了热功率分配的均衡性。

3. 另外，在程序优化的过程中，引入了 addpath 命令，解决了程序对文件夹依赖的问题。该改进使得用户下载程序后能够直接运行，无需手动调整文件夹路径或名称，提升了程序的易用性和操作便捷性，为用户提供了更加流畅的体验和高效的工作环境。这一改进为用户节省了大量时间和精力，使得程序的部署和应用更加便利和无障碍。

📚2 运行结果

由于运行结果图比较多，这里仅展示主场景的运行结果图：

其余运行结果为了美观，去掉Matlab图框。

部分代码：

% 鲁棒优化

% CCHP约束-----日前
% [MPconstrains,func] = MP2(StateParams,StateTemData,Params,NumOfDistract); % MP问题
%%
HasBattery = [0,1,0,1];
constrains = [];
% Func= [];

for i = 1:NumOfDistract
cons = CCHPConstrains2(StateParams(i),Params.(['CCHP',num2str(i)]),HasBattery(i));
cons = [ cons;StateParams(i).Hhe * Params.(['CCHP',num2str(i)]).HE.n + StateParams(i).Hex == Params.(['CCHP',num2str(i)]).HeatLoad'];% Hex 与热网交换（**********）
constrains = [constrains;cons];
Func(i) = Obj2(StateParams(i),Params.(['CCHP',num2str(i)]),Params,i);

end

TemConstrains = [];
for t = 1:24
Hex= [StateParams(1).Hex(t);StateParams(2).Hex(t);StateParams(3).Hex(t);StateParams(4).Hex(t)];
cons = HeatingNetworkConstraints11(StateTemData{t},Params,Hex);

TemConstrains = [TemConstrains;cons];

% 热网耦合约束（*）
end
constrains = [constrains;
TemConstrains;];

func = 0;
for i = 1:NumOfDistract
func = func + Obj2(StateParams(i),Params.(['CCHP',num2str(i)]),Params,i);
end
MPconstrains=constrains;
%%
% 鲁棒约束----实时
theta = sdpvar(1);
NowRobusCost = sdpvar(1,NumOfScence);
NowRobustConstrains = cell(1,NumOfScence);
Cconstrains = MPconstrains;
for i = 1:NumOfScence
[NowConstrainss,NowCost] = SPSingleRobustTest(DifferentScence{i},DifferentTemData{i},Params,StateParams,StateTemData,NumOfDistract,i);
NowRobustConstrains{i} = NowConstrainss;
NowRobusCost(i) = NowCost;
Cconstrains = [Cconstrains;
NowConstrainss;
theta >= NowCost;];
end

opt = sdpsettings('verbose',1,'solver','gurobi');
opt.gurobi.MIPGap=0.1;
result = optimize(Cconstrains,func+theta,opt); % 主问题求解
a1 = value(func+theta)

%% 打印结果
% 日前调度
mkdir('./优化结果')
mkdir('./优化结果/主场景')
ResultPrint(StateParams,StateTemData,NumOfDistract,Params,'主场景');
% 实时结果

for i = 1:NumOfScence
mkdir(['./优化结果/场景',num2str(i)])
ResultPrint(DifferentScence{i},DifferentTemData{i},NumOfDistract,Params,['场景',num2str(i)]);
end