本程序参考SCI期刊论文《Optimal dispatch of zero-carbon-emission micro Energy Internet integrated with non-supplementary fired compressed air energy storage system》，程序中有详细的热网模型，温度控制模块，压缩机模块，综合能源系统规模较大，非常具有可扩展性，算例丰富，注释清晰，干货满满，下面对文章和程序做简要介绍！
为了以清洁和集成的方式利用热和电，通过引入非补充燃烧压缩空气储能（NSF-CAES）技术，提出了零碳排放综合能源系统（ZCE-MEI）架构。本文考虑了一种典型的ZCE-MEI，它将配电网（PDN）和区域供热网（DHN）与NSF-CAES相结合。NSF-CAES架构的制定考虑了热动态和压力行为，以增强调度灵活性。利用改进的DistFlow模型允许几个离散和连续的无功功率补偿器来维持PDN的电压质量。首先将ZCE-MEI的最优操作建模为混合整数非线性规划（MINLP）。通过多次变换和简化，将该问题转化为混合整数线性规划（MILP），利用CPLEX高效求解。采用由NSF-CAES母线、33节点PDN和8节点DHN组成的典型测试系统来验证所提出的ZCE-MEI在降低运营成本和弃风方面的有效性。
文章算例架构：
文章结果：
程序结果：程序图片共25张，包括文章中所展示的，限于篇幅，仅列出部分图片。该程序稍加修改就可以出新的算例文章！！！
部分程序：
clc
close all
clear all
format short
NT = 24; % total dispatch period
Nc = 2; % total stages of compressor
Ne = 2; % total stages of turbine
Sb = 10; % Base power MW
Vb= 12.66; % Base Voltage kV
Zb = Vb^2/Sb; % Base impedence
Ib = Sb/(sqrt(3)*Vb); % kA
%% Power Bus
% No.(1)|Type(2)|Pd(3)|Qd(4)
powerbus = [
1 3 0 0; % MW MVar
2 1 0.100 0.060;
3 1 0.090 0.040;
4 1 0.120 0.080;
5 1 0.060 0.030;
6 1 0.060 0.020;
7 1 0.200 0.100;
8 1 0.200 0.100;
9 1 0.060 0.020;
10 1 0.060 0.020;
11 1 0.045 0.030;
12 1 0.060 0.035;
13 1 0.060 0.035;
14 1 0.120 0.080;
15 1 0.060 0.010;
16 1 0.060 0.020;
17 1 0.060 0.020;
18 1 0.090 0.040;
19 1 0.090 0.040;
20 1 0.090 0.040;
21 1 0.090 0.040;
22 1 0.090 0.040;
23 1 0.090 0.050;
24 1 0.420 0.200;
25 1 0.420 0.200;
26 1 0.060 0.025;
27 1 0.060 0.025;
28 1 0.060 0.020;
29 1 0.120 0.070;
30 1 0.200 0.100;
31 1 0.150 0.070;
32 1 0.210 0.100;
33 1 0.060 0.040;
];
N_bus1 = size(powerbus,1);
Pd_ratio = powerbus(:,3)/sum(powerbus(:,3)); % Active load ratio
Qd_ratio = powerbus(:,4)/sum(powerbus(:,4)); % Reactive load ratio
Pd0 = [63 62 60 58 59 65 72 85 95 99 100 99 93 92 90 88 90 92 96 98 96 90 80 70]/15-1.5;% MW
Qd0 = [18 16 15 14 15.5 15 16 17 18 19 20 20.5 21 20.5 21 19.5 20 20 19.5 19.5 18.5 18.5 18 18]/10; % system load MVar
Pd = Pd_ratio * Pd0;
Qd = Qd_ratio * Qd0;
Pd = Pd/Sb;
Qd = Qd/Sb; % p.u
U2_min = 0.95^2;
U2_max = 1.05^2;
%% Compensator
% Location(1)|Max(2)|Min(3)|Step(4)
ComCap = [
5 0.2 0 0.05;
10 0.2 0 0.05;
13 0.2 0 0.05;
17 0.2 0 0.05;
20 0.2 0 0.05;
23 0.2 0 0.05;
30 0.2 0 0.05;
];
v = 2; % Step number for linearization
N_ComCap = size(ComCap,1); % Number of compensator
Ind_ComCap = ComCap(:,1);
S = ComCap(:,4);
%% SVG
SVG = [ % Mvar % Location(1)|Max(2)|Min(3)
4 0.1 0;
9 0.1 0;
14 0.1 0;
];
Ind_SVG = SVG(:,1);
SVG(:,2:3) = SVG(:,2:3)/Sb; % SVG p.u
%% Heat Node
% No(1)|Hd(2)|Pr_SR_min(3)|tao_S_max(4)|tao_S_min(5)|tao_R_max(6)|tao_R_min(7)|mass flow(8)
heatnode = [ %kW %par %℃
1 0 50000 120 90 80 60 0;
2 0 50000 120 90 80 60 0;
3 0 50000 120 90 80 60 0;
4 0 50000 120 90 80 60 0;
5 250 50000 120 90 80 60 2;
6 250 50000 120 90 80 60 2;
7 250 50000 120 90 80 60 2;
8 500 50000 120 90 80 60 4;
];
N_bus2 = size(heatnode,1);
H_ratio = heatnode(:,2)/sum(heatnode(:,2));
H_hd0 = [1250*ones(1,4), 1150*ones(1,4), 1000*ones(1,4), 800*ones(1,4), 1150*ones(1,4), 1250*ones(1,4)]; % kW
H_Hd = H_ratio * H_hd0;
Nd_Hd = find(heatnode(:,2)>0);
m_Hd = heatnode(Nd_Hd,8); % Mass flow ratio of heat load
tao_NS_max = heatnode(:,4);
tao_NS_min = heatnode(:,5);
tao_NR_max = heatnode(:,6);
tao_NR_min = heatnode(:,7);
%% Power Gen
Ind_gen = [2 7 19 26];
% Wg1 = [11.7 11.3 11.3 12.3 13.5 14.9 16.4 17.2 17.7 18 17.9 17.4 ...
% 16.3 16.1 16.2 16.6 16.8 16.9 16.8 16.6 16.4 16.5 16.6 16.8]/10;% Wind Gen #2(MW) 1.8 MW
% Wg1 = [58.27 82.12 89.22 84.73 77.25 65.13 75.91 71.55 73.40 49.11 30.71 18.09 ...
% 10.80 12.50 15.00 21.62 15.00 10.88 14.50 12.54 16.00 28.41 30.34 37.10]/50; % Wind Gen #1(MW)
% Wg1 = [58.27 82.12 89.22 84.73 77.25 75.13 75.91 71.55 73.40 49.11 30.71 28.09 ...
% 20.80 22.50 25.00 21.62 25.00 20.88 24.50 22.54 26.00 28.41 40.34 47.10]/30; % Wind Gen #1(MW)
Wg1 = [6.88 7.08 7.20 7.16 6.96 6.52 6.44 5.98 5.72 5.54 5.36 5.12 ...
4.64 4.56 4.60 4.64 4.52 4.52 4.92 5.40 5.96 6.56 6.68 6.72]-4.2;% Wind Gen #1(MW) MW 3MW
Wg2 = [16.6 16.4 16.5 16.6 16.8 11.7 11.3 11.3 12.3 13.5 14.9 16.4 17.2 17.7 18 17.9 17.4 ...
16.3 16.1 16.2 16.6 16.8 16.9 16.8]/40;% Wind Gen #2(MW) 0.6 MW
Wg3 = Wg2;
Wg4 = Wg2;
Wg1 = Wg1/Sb; % p.u
Wg2 = Wg2/Sb;
Wg3 = Wg3/Sb;
Wg4 = Wg4/Sb;
Wg = [Wg1;Wg2;Wg3;Wg4];
%% Heat Gen
% Location(1)|Hg_max(2)|Hg_min(3)|C_A(4)|C_B(5)|C_(6)|Mass flow(7) %
heatgen = [% kW % kg/s
2 2500 0 0.05 20 0 10
];
N_gen = size(heatgen,1);
Nd_HS = heatgen(:,1); % Node of heat station
m_HS = heatgen(:,7);
Hg_min = heatgen(:,3);
Hg_max = heatgen(:,2);
%% CAES Hub
yita_comp = [0.80, 0.75];
yita_turb = [0.86, 0.86];
beta_comp = [11.6,8.15];
pi_turb = [8.9,8.9];
Pcomp_min = zeros(Nc,1);
Pcomp_max = 500*ones(Nc,1); %kW
Pturb_min = zeros(Ne,1);
Pturb_max = 1000*ones(Ne,1); %kW
Vst = 2000; % m^3
% Pcomp_min = zeros(Nc,1);
% Pcomp_max = 1500*ones(Nc,1); %kW
% Pturb_min = zeros(Ne,1);
% Pturb_max = 3000*ones(Ne,1); %kW
%
% Vst = 10000; % m^3
k = 1.4; % adiabatic exponent
Rg = 0.297; % KJ/(kg.K)
cp_a = 1.007; % KJ/(kg.K) 25℃
cp_w = 4.2; % KJ/(kg.K) 25℃
cp_s = 2.5; % KJ/(kg.K) 25℃
tao_am = 15; % ambient temperature
tao_K = 273.15;
tao_am = tao_am + tao_K;
tao_str = 40;
tao_str = tao_str + tao_K;
tao_salt_min = 60 + tao_K;
tao_salt_max = 320 + tao_K;
pr_am = 0.101*1e3; % Kpa ambient pressure
pr_st_min = 8.4*1e3; % Kpa
pr_st_max = 9.0*1e3; % Kpa
qm_comp_min = 0;
qm_comp_max = 2.306/3.6; % kg/s 1MW CAES
qm_turb_min = 0;
qm_turb_max = 8.869/3.6; % kg/s 1MW CAES
H_str_min = 0.2*1e3; %kW
H_str_max = 3.0*1e3; %kW
% qm_comp_min = 0;
% qm_comp_max = 9.3/3.6; % kg/s
% qm_turb_min = 0;
% qm_turb_max = 43.7/3.6; % kg/s
% H_str_min = 1e3; %kW.h
% H_str_max = 20*1e3;%kW.h
pr_comp_in1 = pr_am*ones(1,NT); % Fix pressure
pr_comp_out1 = beta_comp(1)*pr_comp_in1;
pr_comp_in2 = pr_comp_out1;
pr_comp_out2 = beta_comp(2)*pr_comp_in2;
y_comp1 = (beta_comp(1))^((k-1)/k);
y_comp2 = (beta_comp(2))^((k-1)/k);
pr_turb_in1 = pr_st_min*ones(1,NT);
pr_turb_out1 = pr_turb_in1/pi_turb(1);
pr_turb_in2 = pr_turb_out1;
pr_turb_out2 = pr_turb_in2/pi_turb(2);
y_turb1 = (pi_turb(1))^(-(k-1)/k);
y_turb2 = (pi_turb(2))^(-(k-1)/k);
tao_comp_in1 = tao_am*ones(1,NT); % Fix pressure
tao_comp_in2 = (40 + tao_K)*ones(1,NT);
tao_comp_out1 = tao_comp_in1/yita_comp(1).*(y_comp1-1+yita_comp(1));
tao_comp_out2 = tao_comp_in2/yita_comp(2).*(y_comp2-1+yita_comp(2));
tao_cold_s_in1 = tao_comp_out1;
tao_cold_s_out1 = 90 + tao_K ; % Fix salt heat exchanger output temperature 90℃
tao_cold_w_in1 = tao_cold_s_out1;
tao_cold_w_out1 = 40 + tao_K ; % water heat exchanger output temperature 40℃
tao_cold_s_in2 = tao_comp_out2;
tao_cold_s_out2 = 90 + tao_K;
tao_cold_w_in2 = tao_cold_s_out2;
tao_cold_w_out2 = 40 + tao_K ;
tao_turb_in1 = (280 + tao_K)*ones(1,NT);
tao_turb_in2 = (280 + tao_K)*ones(1,NT);
tao_turb_out1 = tao_turb_in1*yita_turb(1).*(y_turb1-1+1/yita_turb(1));
tao_turb_out2 = tao_turb_in2*yita_turb(2).*(y_turb2-1+1/yita_turb(2));
tao_heat_in1 = tao_str;
tao_heat_in2 = tao_turb_out1;
tao_heat_out1 = tao_turb_in1;
tao_heat_out2 = tao_turb_in2;
CAES_ind = 2;
%% Power Line
% No.(1)|From bus(2)|To bus(3)|r(4)|x(5)|P_line_max(6)|P_line_min(7) %
branch = [
1 1 2 0.0922 0.0470 9.9 0;
2 2 3 0.4930 0.2512 9.9 0;
3 3 4 0.3661 0.1864 9.9 0;
4 4 5 0.3811 0.1941 9.9 0;
5 5 6 0.8190 0.7070 9.9 0;
6 6 7 0.1872 0.6188 9.9 0;
7 7 8 0.7115 0.2351 9.9 0;
8 8 9 1.0299 0.7400 9.9 0;
9 9 10 1.0440 0.7400 9.9 0;
10 10 11 0.1967 0.0651 9.9 0;
11 11 12 0.3744 0.1298 9.9 0;
12 12 13 1.4680 1.1549 9.9 0;
13 13 14 0.5416 0.7129 9.9 0;
14 14 15 0.5909 0.5260 9.9 0;
15 15 16 0.7462 0.5449 9.9 0;
16 16 17 1.2889 1.7210 9.9 0;
17 17 18 0.7320 0.5739 9.9 0;
18 2 19 0.1640 0.1565 9.9 0;
19 19 20 1.5042 1.3555 9.9 0;
20 20 21 0.4095 0.4784 9.9 0;
21 21 22 0.7089 0.9373 9.9 0;
22 3 23 0.4512 0.3084 9.9 0;
23 23 24 0.8980 0.7091 9.9 0;
24 24 25 0.8959 0.7071 9.9 0;
25 6 26 0.2031 0.1034 9.9 0;
26 26 27 0.2842 0.1447 9.9 0;
27 27 28 1.0589 0.9338 9.9 0;
28 28 29 0.8043 0.7006 9.9 0;
29 29 30 0.5074 0.2585 9.9 0;
30 30 31 0.9745 0.9629 9.9 0;
31 31 32 0.3105 0.3619 9.9 0;
32 32 33 0.3411 0.5302 9.9 0;
];
N_line = size(branch,1);
line_i = branch(:,2);
line_j = branch(:,3);
r = branch(:,4)/Zb;
x = branch(:,5)/Zb;
Pmax = branch(:,6)/Sb;
Pmin = branch(:,7)/Sb;
%% OLTC
% Line No.(1)|K_max(2)|K_min(3)|K_Step(4)|%
OLTC = [
1 1.05 0.95 0.01;
18 1.05 0.95 0.01;
22 1.05 0.95 0.01;
25 1.05 0.95 0.01;
];
t_OLTC = 0.95:0.01:1.05; % available tap value
T_OLTC = repmat(t_OLTC',1,NT);
n_OLTC= length(t_OLTC); % num of total tap value
N_OLTC = size(OLTC,1);% num of OLTC
Ind_OLTC = OLTC(:,1);
Ind_subline = zeros(N_line,2); % Index of Children line of power netwrok
for i = 1: N_line
temp = find(line_i == line_j(i));
if ~isempty(temp)
Ind_subline(i,1:length(temp)) = temp;
end
end
Pg_min = zeros(N_bus1,NT);
Pg_max = zeros(N_bus1,NT);
以上就是本次介绍的主要内容，小编非常推荐该程序代码，该模型对综合能源系统方向和压缩空气储能方向的小伙伴有很大帮助，代码的可扩展性强，欢关注下方公众号获取完整版代码，小编会继续推送更有质量的学习资料、文章和程序代码！
