请你检查以下代码，说明需要优化或修改的地方：%% 考虑零碳排放的综合能源系统优化调度方法（电网+热网）
clc
close
%% 一、控制参数设置
N_Time = 24; % 调度周期（小时）
S_base = 10; % 基准功率 MW
V_base = 12.66; % 基准电压 kV
Z_base = V_base^2/S_base; % 基准阻抗
%% ==================== 风电机组参数配置 ====================
Ind_gen = [2, 7, 19, 26]; % 风机并网的节点
P_WT_g1 = [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; % 风电机组 #1(MW) MW 3MW
P_WT_g2 = [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; % 风电机组 #2(MW) 0.6 MW
P_WT_g3 = P_WT_g2;
P_WT_g4 = P_WT_g2;
P_WT_g1 = P_WT_g1/S_base; % 标幺化处理
P_WT_g2 = P_WT_g2/S_base;
P_WT_g3 = P_WT_g3/S_base;
P_WT_g4 = P_WT_g4/S_base;
P_WT_g = [P_WT_g1;P_WT_g2;P_WT_g3;P_WT_g4]; % 4台风机
%% ==================== 电网节点参数配置 ====================
% 电源母线参数 (MW、MVar）
% No.(1) |Type(2)| Pd(3)| Qd(4)
power_bus = [
1 3 0 0;
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_bus_1 = size(power_bus,1); % 33节点网络
Pd_ratio = power_bus(:,3)/sum(power_bus(:,3)); % 有功负载比
Qd_ratio = power_bus(:,4)/sum(power_bus(:,4)); % 无功负载比
Pd_0 = [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
Qd_0 = [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; % 基准无功功率 MVar
Pd = Pd_ratio * Pd_0; % 节点电负荷计算
Qd = Qd_ratio * Qd_0;
Pd = Pd/S_base; % 标幺化有功负载
Qd = Qd/S_base; % 标幺化无功负载
U2_min = 0.95^2; % 电压边界设置
U2_max = 1.05^2;
%% ==================== 电网支路参数配置 ====================
% 电力线支路参数
% 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);
B_line_i = branch(:,2);
B_line_j = branch(:,3);
r_pu = branch(:,4)/Z_base; % 电阻标幺值
x_pu = branch(:,5)/Z_base; % 电抗标幺值
P_max_pu = branch(:,6)/S_base; % 功率限值标幺化
P_min_pu = branch(:,7)/S_base;
%% ==================== 无功补偿设备配置 ====================
% 电容无功补偿器参数
% Location(1) |Max(2) |Min(3) |Step(4)
Com_Cap = [
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;
N_Com_Cap = size(Com_Cap,1); % 补偿器数量
Ind_Com_Cap = Com_Cap(:,1); % 补偿器位置
Com_Cap_Step = Com_Cap(:,4); % 电容步长调整
% 静止无功补偿器参数
SVG = [ % Location(1) |Max(2)| Min(3)
4 0.1 0;
9 0.1 0;
14 0.1 0;
];
Ind_SVG = SVG(:,1); % SVG数量
SVG(:, 2:3) = SVG(:, 2:3) / S_base; % 标幺化
%% ==================== 可调变压器配置 ====================
% 可调变压器参数
% 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]; % 可用抽头值
T_OLTC = repmat(t_OLTC',1,N_Time);
n_OLTC= length(t_OLTC); % 总损失价值
N_OLTC = size(OLTC,1); % 变压器数量
Ind_OLTC = OLTC(:,1); % 变压器所在线路索引
Ind_sub_line = zeros(N_line,2); % 电网子线路索引
for i = 1: N_line
temp = find(B_line_i == B_line_j(i));
if ~isempty(temp)
Ind_sub_line(i,1:length(temp)) = temp;
end
end
Pg_min = zeros(N_bus_1,N_Time);
Pg_max = zeros(N_bus_1,N_Time);
%% ==================== 热网节点参数配置 ====================
% 热网节点参数（kW、par、℃）
% 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)
heat_node = [
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_bus_2 = size(heat_node,1);
H_ratio = heat_node(:,2)/sum(heat_node(:,2));
H_hd_0 = [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_hd_0; % 节点热负荷计算
Nd_Hd = find(heat_node(:,2)>0);
m_Hd = heat_node(Nd_Hd,8); % 热负荷质量流量比
tao_NS_max = heat_node(:,4); % 温度边界设置
tao_NS_min = heat_node(:,5);
tao_NR_max = heat_node(:,6); % 温度边界设置
tao_NR_min = heat_node(:,7);
%% ==================== 热力发电参数配置 ====================
% 热力发电参数（kW、kg/s）
% Location(1) |Hg_max(2) |Hg_min(3) |C_A(4) |C_B(5) |C_(6) |Mass flow(7)
heat_gen = [
2 2500 0 0.05 20 0 10
];
N_gen = size(heat_gen,1); % 热力机组数量
Nd_HS = heat_gen(:,1); % 热力站的节点位置
mf_HS = heat_gen(:,7); % 热力机组质量流量
Hg_min = heat_gen(:, 3); % 最小出力
Hg_max = heat_gen(:, 2); % 最大出力
%% ==================== 热网管道参数配置 ====================
% 热网传输通道
% No.(1) |From node(2) |To node(3) |L(4) |u_p(5) |u_T(6) |ms_max(7) |ms_min(8) |mr_max(9) |mr_min(10)
pipe = [
1 1 2 1000 0.04 0.5*1e-3 1e6 0 10 10;
2 2 3 1000 0.04 0.5*1e-3 1e6 0 8 8;
3 3 4 1000 0.04 0.5*1e-3 1e6 0 6 6;
4 2 5 1000 0.04 0.5*1e-3 1e6 0 2 2;
5 3 6 1000 0.04 0.5*1e-3 1e6 0 2 2;
6 4 7 1000 0.04 0.5*1e-3 1e6 0 2 2;
7 4 8 1000 0.04 0.5*1e-3 1e6 0 4 4;
];
N_pipe = size(pipe,1); % 获取管道数量
pipe_i = pipe(:,2); % 提取节点i首端
pipe_j = pipe(:,3); % 提取节点j末端
L_pipe = pipe(:,4); % 管道长度
lamada_pipe = pipe(:,6); % 热传导系数
m_pipe_max = pipe(:, 7); % 最大流量
m_pipe_min = pipe(:, 8); % 最小流量
ms_pipe = pipe(:, 9); % 发送端流量
mr_pipe = pipe(:, 10); % 接收端流量
S_pipe_F = zeros(N_bus_2,2); % 发送端拓扑矩阵
S_pipe_T = zeros(N_bus_2,2); % 接收端拓扑矩阵
%% ==================== 压缩空气储能参数配置 ====================
% 集线器参数 CAES (kW)
N_ct = 2; % 压缩机数量
N_et = 2; % 涡轮机数量
% 运行机组参数
yita_comp = [0.80, 0.75]; % 压缩机等熵效率
yita_turb = [0.86, 0.86]; % 涡轮机等熵效率
beta_comp = [11.6, 8.15]; % 压缩机压力比
pi_turb = [8.9, 8.9];
P_comp_min = zeros(N_ct, 1);
P_comp_max = 500 * ones(N_ct, 1);
P_turb_min = zeros(N_et, 1);
P_turb_max = 1000 * ones(N_et, 1);
% 物理常量参数
Vst = 2000; % m^3
k = 1.4; % 绝热指数
Rg = 0.297; % KJ/(kg.K)
cp_a = 1.007; % KJ/(kg.K)
cp_w = 4.2; % KJ/(kg.K)
cp_s = 2.5; % KJ/(kg.K)
% 环境温度参数
tao_K = 273.15; % 摄氏转开尔文偏移量
tao_am = 15 + tao_K; % 直接计算开尔文温度
tao_str = 40 + tao_K; % 介质初始温度
tao_salt_min = 60 + tao_K;
tao_salt_max = 320 + tao_K;
% 环境压力参数
pr_am = 0.101 * 1e3; % kPa
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;
qm_turb_min = 0;
qm_turb_max = 8.869 / 3.6;
% 热功率上、下限参数
H_str_min = 0.2*1e3; % kW
H_str_max = 3.0*1e3; % kW
% ---------- 压力计算（向量化）----------
% 预压缩计算
k_ratio = (k - 1) / k; % 复用指数计算
% 压缩机压力计算
pr_comp_in_1 = pr_am * ones(1, N_Time);
pr_comp_out_1 = beta_comp(1) * pr_comp_in_1;
pr_comp_in_2 = pr_comp_out_1;
pr_comp_out_2 = beta_comp(2) * pr_comp_in_2;
% 涡轮机压力计算
pr_turb_in_1 = pr_st_min * ones(1, N_Time);
pr_turb_out_1 = pr_turb_in_1 / pi_turb(1);
pr_turb_in_2 = pr_turb_out_1;
pr_turb_out_2 = pr_turb_in_2 / pi_turb(2);
% ---------- 温度计算（复用计算因子）----------
% 压缩过程温度因子计算
y_comp = beta_comp .^ k_ratio; % 向量化计算
tao_comp_in_1 = tao_am * ones(1, N_Time);
tao_comp_in_2 = (40 + tao_K) * ones(1, N_Time);
tao_comp_out_1 = tao_comp_in_1 ./ yita_comp(1) .* (y_comp(1) - 1 + yita_comp(1));
tao_comp_out_2 = tao_comp_in_2 ./ yita_comp(2) .* (y_comp(2) - 1 + yita_comp(2));
% 热交换温度计算
tao_cold_s_in_1 = tao_comp_out_1;
tao_cold_s_out_1 = 90 + tao_K; % 直接计算开尔文温度
tao_cold_w_in_1 = tao_cold_s_out_1;
tao_cold_w_out_1 = 40 + tao_K;
tao_cold_s_in_2 = tao_comp_out_2;
tao_cold_s_out_2 = 90 + tao_K;
tao_cold_w_in_2 = tao_cold_s_out_2;
tao_cold_w_out_2 = 40 + tao_K;
% 涡轮温度计算（复用计算因子）
y_turb = pi_turb .^ (-k_ratio); % 向量化计算
tao_turb_in_1 = (280 + tao_K) * ones(1, N_Time);
tao_turb_in_2 = (280 + tao_K) * ones(1, N_Time);
tao_turb_out_1 = tao_turb_in_1 .* yita_turb(1) .* (y_turb(1) - 1 + 1./yita_turb(1));
tao_turb_out_2 = tao_turb_in_2 .* yita_turb(2) .* (y_turb(2) - 1 + 1./yita_turb(2));
% 热能传递计算
tao_heat_in_1 = tao_str;
tao_heat_in_2 = tao_turb_out_1;
tao_heat_out_1 = tao_turb_in_1;
tao_heat_out_2 = tao_turb_in_2;
% 系统标识
CAES_ind = 2;
%% 二、运行变量设置
%% ==================== 电网运行变量 ====================
P_line_load = sdpvar(N_line,N_Time,'full'); % 每条线路的有功功率
Q_line_load = sdpvar(N_line,N_Time,'full'); % 每条线路的无功功率
P_sub = sdpvar(N_line,N_Time,'full');
Q_sub = sdpvar(N_line,N_Time,'full');
U2_node = sdpvar(N_bus_1,N_Time,'full'); % 母线电压幅值的平方
Pg_node = sdpvar(N_bus_1,N_Time,'full'); % 向每条母线注入发电机有功功率
Qc_node = sdpvar(N_bus_1,N_Time,'full'); % 无功补偿器SVG
P_grid = sdpvar(1,N_Time); % 上网购电有功
Q_grid = sdpvar(1,N_Time); % 上网购电无功
%% ==================== 热网运行变量 ====================
tao_PS_src = sdpvar(N_pipe,N_Time,'full'); % 供水管“源”侧温度
tao_PS_load = sdpvar(N_pipe,N_Time,'full'); % 供水管“荷”侧温度
tao_PR_src = sdpvar(N_pipe,N_Time,'full'); % 回流管“源”侧温度
tao_PR_load = sdpvar(N_pipe,N_Time,'full'); % 回流管“荷”侧温度
tao_NS_node = sdpvar(N_bus_2,N_Time,'full'); % 供电网络节点温度
tao_NR_node = sdpvar(N_bus_2,N_Time,'full'); % 回流网络节点温度
Hg_HP = sdpvar(N_bus_2,N_Time,'full'); % 配备CAES的热泵的热功率
%% ==================== CAES系统运行变量 ====================
on_comp = binvar(1,N_Time); % 压缩机 comp 的开/关
on_turb = binvar(1,N_Time); % 涡轮机 turb 的开/关
P_comp_1 = sdpvar(1,N_Time); % 压缩机功率 comp1
P_comp_2 = sdpvar(1,N_Time); % 压缩机功率 comp2
P_turb_1 = sdpvar(1,N_Time); % 涡轮机功率 turb1
P_turb_2 = sdpvar(1,N_Time); % 涡轮机功率 turb2
P_caes_d = sdpvar(1,N_Time);
P_caes_g = sdpvar(1,N_Time);
%% ==================== 线性化辅助变量 ====================
pr_st = sdpvar(1,N_Time); % 储气罐压强
pr_st_0 = sdpvar(1,1);
qm_comp = sdpvar(1,N_Time);
qm_turb = sdpvar(1,N_Time);
y1 = sdpvar(1,N_Time); % 线性化储气室压强约束
y2 = sdpvar(1,N_Time);
h1 = sdpvar(1,N_Time); % 线性化储热系统SOC
h2 = sdpvar(1,N_Time);
HM = 1e7; % 大M法常数
H_coll_s1 = sdpvar(1,N_Time);
H_coll_s2 = sdpvar(1,N_Time);
H_cons1 = sdpvar(1,N_Time);
H_cons2 = sdpvar(1,N_Time);
H_coll_sum = sdpvar(1,N_Time);
H_cons_sum = sdpvar(1,N_Time);
H_str = sdpvar(1,N_Time); % 储热罐中存贮的热量
H_str_0 = sdpvar(1,1);
Hg_CAES = sdpvar(1,N_Time); % 蓄热环节可供热负荷
% 无功补偿器分段线性化
for i = 1:N_Com_Cap % 电容器
xd_cap{i} = binvar(v+1,N_Time);
delta_cap{i} = sdpvar(v+1,N_Time);
end
for i = 1:N_OLTC % 变压器
rd_tap{i} = binvar(n_OLTC,N_Time);
h_tap{i} = sdpvar(n_OLTC,N_Time);
end
%% 三、约束条件设置
%% ==================== CAES系统约束 ====================
F_turb = []; % 每级功率定义
F_comp = []; % 每级功率定义
F_oper = []; % 运行约束
F_power = []; % 功率平衡约束
F_air_str = []; % 储气罐压强动态约束
F_cold = [];
F_heat = [];
F_heat_str = [];
% 压缩机约束
F_comp = [F_comp, P_comp_1 == 1/yita_comp(1)*k/(k-1)*Rg*qm_comp.*tao_comp_in_1.*(y_comp1-1)];
F_comp = [F_comp, P_comp_2 == 1/yita_comp(2)*k/(k-1)*Rg*qm_comp.*tao_comp_in_2.*(y_comp2-1)];
F_comp = [F_comp, P_comp_min(1)*on_comp <= P_comp_1 <= P_comp_max(1)*on_comp ]; % 每级消耗的功率约束
F_comp = [F_comp, P_comp_min(2)*on_comp <= P_comp_2 <= P_comp_max(2)*on_comp ];
F_comp = [F_comp, P_caes_d == P_comp_1 + P_comp_2]; % 总功率定义
% 涡轮机约束
F_turb = [F_turb, P_turb_1 == yita_turb(1)*k/(k-1)*Rg*qm_turb.*tao_turb_in_1.*(1-y_turb1)];
F_turb = [F_turb, P_turb_2 == yita_turb(2)*k/(k-1)*Rg*qm_turb.*tao_turb_in_2.*(1-y_turb2)];
F_turb = [F_turb, P_turb_min(1)*on_turb <= P_turb_1 <= P_turb_max(1)*on_turb]; % 每级发出功率约束
F_turb = [F_turb, P_turb_min(2)*on_turb <= P_turb_2 <= P_turb_max(2)*on_turb];
F_turb = [F_turb, P_caes_g == P_turb_1 + P_turb_2];
% 机组运行约束
F_oper = [F_oper, 0 <= on_comp + on_turb <= 1]; % 充放电不能同时进行
F_oper = [F_oper, qm_comp_min*on_comp <= qm_comp <= qm_comp_max*on_comp]; % 质量流量非否约束
F_oper = [F_oper, qm_turb_min*on_turb <= qm_turb <= qm_turb_max*on_turb];
% 边界约束和质量流量约束
F_air_str = [F_air_str, pr_st(1) == pr_st_0 + 1/Vst * Rg * tao_str * 3600*(y1(1)- y2(1))];
F_air_str = [F_air_str, pr_st(1,2:N_Time) == pr_st(1,1:N_Time-1) + 1/Vst * Rg * tao_str * 3600*(y1(2:N_Time)- y2(2:N_Time))];
F_air_str = [F_air_str, pr_st(1,N_Time) == pr_st_0];
F_air_str = [F_air_str, qm_comp_min*on_comp <= y1 <= qm_comp_max*on_comp];
F_air_str = [F_air_str, qm_turb_min*on_turb <= y2 <= qm_turb_max*on_turb];
F_air_str = [F_air_str, qm_comp_min*(1-on_comp) <= qm_comp - y1 <= qm_comp_max*(1-on_comp)];
F_air_str = [F_air_str, qm_turb_min*(1-on_turb) <= qm_turb - y2 <= qm_turb_max*(1-on_turb)];
F_air_str = [F_air_str, pr_st_min <= pr_st <= pr_st_max];
F_air_str = [F_air_str, pr_st_min <= pr_st_0 <= pr_st_max];
% 热回收系统约束 - 压缩热回收
F_cold = [F_cold, H_coll_s1 == cp_a*qm_comp.*(tao_cold_s_in_1 - tao_cold_s_out_1)]; % 回热系统热量约束，冷却环节
F_cold = [F_cold, H_coll_s2 == cp_a*qm_comp.*(tao_cold_s_in_2 - tao_cold_s_out_2)]; % 第1,2级导热油回收的热量
F_cold = [F_cold, H_coll_sum == H_coll_s1 + H_coll_s2]; % 只计及导热油收集的热量，回收的总热量
% 热回收系统约束 - 膨胀热消耗
F_heat = [F_heat, H_cons1 == cp_a*qm_turb.*(tao_heat_out_1 - tao_heat_in_1)]; % 第1-2级消耗的热量
F_heat = [F_heat, H_cons2 == cp_a*qm_turb.*(tao_heat_out_2 - tao_heat_in_2)];
F_heat = [F_heat, H_cons_sum == H_cons1 + H_cons2]; % 加热环节消耗的总热量
% 线性化约束和边界约束
F_heat_str = [F_heat_str, H_str(1,1) == H_str_0 + h1(1,1) - h2(1,1) - Hg_CAES(1,1)]; % 回热系统SOC% 高温储热罐中存贮的热量
F_heat_str = [F_heat_str, H_str(1,2:N_Time) == H_str(1,1:N_Time-1) + h1(1,2:N_Time) - h2(1,2:N_Time) - Hg_CAES(1,2:N_Time)];
F_heat_str = [F_heat_str, -HM*on_comp <= h1 <= HM*on_comp];
F_heat_str = [F_heat_str, -HM*on_turb <= h2 <= HM*on_turb];
F_heat_str = [F_heat_str, -HM*(1-on_comp) <= H_coll_sum - h1 <= HM*(1-on_comp)];
F_heat_str = [F_heat_str, -HM*(1-on_turb) <= H_cons_sum - h2 <= HM*(1-on_turb)];
F_heat_str = [F_heat_str, H_str_min <= H_str <= H_str_max]; % Box 约束
F_heat_str = [F_heat_str, H_str_min <= H_str_0 <= H_str_max];
F_heat_str = [F_heat_str, H_str(1,N_Time) == H_str_0];
%% ==================== 电力网络约束 ====================
F_P = [];
F_Q = [];
F_U = [];
% 节点功率约束
for t = 1:N_Time
for i = 1:N_bus_1 % 设置未放置发电机的母线节点的注入(P，Q)功率及a,b,c为 0， 发电机与SVG节点约束
if isempty(find(Ind_gen == i))
F_P = [F_P, Pg_node(i,t) == 0];
else
Pg_max(i,:) = P_WT_g(find(Ind_gen == i),:);
end
% 判断是否装有SVG
if isempty(find(Ind_SVG ==i)) % 未安装SVG
F_Q = [F_Q, Qc_node(i,t) == 0];
else % 安装SVG
temp = find(Ind_SVG == i); % 索引
F_Q = [F_Q, SVG(temp,3) <= Qc_node(i,t) <= SVG(temp,2)];
end
end
end
% 线路子线路潮流
for t = 1:N_Time
for i = 1:N_line
num_temp = size(find(Ind_sub_line(i,:) == 0),2);
if num_temp == 1
F_P = [F_P, P_sub(i,t) == P_line_load(Ind_sub_line(i,1),t)];
F_Q = [F_Q, Q_sub(i,t) == Q_line_load(Ind_sub_line(i,1),t)];
elseif num_temp == 2
F_P = [F_P, P_sub(i,t) == 0];
F_Q = [F_Q, Q_sub(i,t) == 0];
else
F_P = [F_P, P_sub(i,t) == P_line_load(Ind_sub_line(i,1),t) + P_line_load(Ind_sub_line(i,2),t)];
F_Q = [F_Q, Q_sub(i,t) == Q_line_load(Ind_sub_line(i,1),t) + Q_line_load(Ind_sub_line(i,2),t)];
end
end
end
% 初始化计数器
OTLC_count = zeros(1,N_Time);
ComCap_count = zeros(1,N_Time);
M = 1000;
% 基准节点约束
F_P = [F_P, P_line_load(1,:) == P_grid(1,:)];
F_Q = [F_Q, Q_line_load(1,:) == 0];
F_U = [F_U, U2_node(1,:) == 1.05^2];
Vs1 = 1.05^2 * ones(1,N_Time);
yite = 0.8;
% 母线功率平衡与设备约束
for t = 1:N_Time
for i = 1:N_line
if B_line_j(i) == CAES_ind
CAES_ind;
F_P = [F_P, P_line_load(i,t) + Pg_node(B_line_j(i),t) + P_caes_g(t)/1e3/S_base == P_sub(i,t) + Pd(B_line_j(i),t) + P_caes_d(t)/1e3/S_base];
else
F_P = [F_P, P_line_load(i,t) + Pg_node(B_line_j(i),t) == P_sub(i,t) + Pd(B_line_j(i),t)];
end
% 判断是否装有并联补偿电容
if ~isempty(find(Ind_Com_Cap == B_line_j(i))) % 安装补偿电容
ComCap_count(t) = ComCap_count(t) + 1;
F_Q = [F_Q, Q_line_load(i,t) + 0.5*(U2_node(B_line_j(i),t)*Com_Cap(ComCap_count(t),3) + Com_Cap_Step(ComCap_count(t))*(2^0*delta_cap{ComCap_count(t)}(1,t) + ...
2^1*delta_cap{ComCap_count(t)}(2,t) + 2^2*delta_cap{ComCap_count(t)}(2,t)))+ Qc_node(B_line_j(i),t) == Q_sub(i,t) + Qd(B_line_j(i),t)];
for m = 1:v+1
F_U = [F_U, U2_node(B_line_j(i),t) - M*(1-xd_cap{ComCap_count(t)}(m,t)) <= delta_cap{ComCap_count(t)}(m,t) <= U2_node(B_line_j(i),t) + M*(1-xd_cap{ComCap_count(t)}(m,t))];
F_U = [F_U, -M*xd_cap{ComCap_count(t)}(m,t) <= delta_cap{ComCap_count(t)}(m,t) <= M*xd_cap{ComCap_count(t)}(m,t)];
end
F_U = [F_U, 0 <= 2^0*xd_cap{ComCap_count(t)}(1,t)+ 2^1*xd_cap{ComCap_count(t)}(2,t) + 2^2*xd_cap{ComCap_count(t)}(3,t) <= (Com_Cap(ComCap_count(t),2) - ...
Com_Cap(ComCap_count(t),3))/Com_Cap(ComCap_count(t),4)];
else % 未安装补偿电容
F_Q = [F_Q, Q_line_load(i,t) + Qc_node(B_line_j(i),t) == Q_sub(i,t) + Qd(B_line_j(i),t)];
end
if ~isempty(find(Ind_OLTC == i)) % 含OTLC的支路
OTLC_count(t) = OTLC_count(t)+1;
F_U = [F_U, sum(h_tap{OTLC_count(t)}(:,t)./T_OLTC(:,t).^2,1) == U2_node(B_line_i(i),t)-(r_pu(i)*P_line_load(i,t)+x_pu(i)*Q_line_load(i,t))/Vs1(1,t)];
for k = 1:n_OLTC
F_U = [F_U, -M*(1-rd_tap{OTLC_count(t)}(k,t)) + U2_node(B_line_j(i),t) <= h_tap{OTLC_count(t)}(k,t) <= U2_node(B_line_j(i),t) + M*(1-rd_tap{OTLC_count(t)}(k,t))];
F_U = [F_U, -M*rd_tap{OTLC_count(t)}(k,t) <= h_tap{OTLC_count(t)}(k,t)<= M*rd_tap{OTLC_count(t)}(k,t)];
end
F_U = [F_U,sum(rd_tap{OTLC_count(t)}(:,t),1) == 1];
else % 不含OTLC的支路
F_U = [F_U, U2_node(B_line_j(i),t)== U2_node(B_line_i(i),t)-(r_pu(i)*P_line_load(i,t)+x_pu(i)*Q_line_load(i,t))/Vs1(1,t)];
end
% 线路潮流约束
F_P = [F_P, P_min_pu(i) <= P_line_load(i,t) <= P_max_pu(i)];
end
for i = 1:N_bus_1
F_P = [F_P,Pg_min(i,t) <= Pg_node(i,t) <= Pg_max(i,t)];
F_U = [F_U, U2_min <= U2_node(i,t) <= U2_max];
end
end
F_P = [F_P, P_grid >= 0];
%% ==================== 热力站约束 ====================
F_H = [];
Count_HS = 0;
Count_Hd = 0;
for j = 1:N_bus_2
if ~isempty(find(Nd_HS == j)) % 若该节点设置供热机组
Count_HS = Count_HS + 1;
F_H = [F_H, Hg_HP(j,:) + Hg_CAES(1,:) == cp_w*mf_HS(Count_HS)*(tao_NS_node(j,:) - tao_NR_node(j,:))]; % 增加CAES产生的热能
F_H = [F_H, Hg_min <= Hg_HP(j,:) <= Hg_max];
elseif ~isempty(find(Nd_Hd == j)) % 若该节点设置供热负荷
Count_Hd = Count_Hd + 1;
F_H = [F_H, Hg_HP(j,:) == zeros(1,N_Time)];
F_H = [F_H, cp_w*m_Hd(Count_Hd,:)*(tao_NS_node(j,:) - tao_NR_node(j,:)) == H_Hd(j,:)];
else % 其他节点
F_H = [F_H, Hg_HP(j,:) == zeros(1,N_Time)];
end
F_H = [F_H, tao_NS_min(j) <= tao_NS_node(j,:) <= tao_NS_max(j)];
F_H = [F_H, tao_NR_min(j) <= tao_NR_node(j,:)<= tao_NR_max(j)];
end
F_H = [F_H, 0 <= Hg_CAES <= 0.2*H_str];
%% ==================== 供热网络约束 ====================
F_PH = [];
for i = 1:N_bus_2
temp1 = find(pipe_i == i); % 管道首节点为i
temp2 = find(pipe_j == i); % 管道末节点为i
S_pipe_F(i,1:length(temp1)) = temp1; % 首节点为i的管道集合
S_pipe_T(i,1:length(temp2)) = temp2; % 末节点为i的管道集合
end
for i = 1:N_bus_2
% 供水管道温度节点混合
num_temp1 = size(find(S_pipe_T(i,:) == 0),2);
if num_temp1 == 1 % 以节点i末端的管道数为1
b = S_pipe_T(i,1); % 管道编号
F_PH = [F_PH, ms_pipe(b)* tao_PS_load(b,:) == ms_pipe(b) * tao_NS_node(i,:)];
F_PH = [F_PH, tao_PR_src(b,:) == tao_NR_node(i,:)];
elseif num_temp1 == 0 % 以节点i末端的管道数为2
b = S_pipe_T(i,:); % 管道编号
F_PH = [F_PH, ms_pipe(b(1))* tao_PS_load(b(1),:) + ms_pipe(b(2))* tao_PS_load(b(2),:) == (ms_pipe(b(1))+ms_pipe(b(2))) * tao_NS_node(i,:)];
F_PH = [F_PH, tao_PR_src(b(1),:) == tao_NR_node(i,:)];
F_PH = [F_PH, tao_PR_src(b(2),:) == tao_NR_node(i,:)];
else
% 以节点i末端的管道数为0 ,首节点
end
% 回水管道温度节点混合
num_temp2 = size(find(S_pipe_F(i,:) == 0),2);
if num_temp2 == 1 % 以节点i首端的管道数为1
b = S_pipe_F(i,1); % 管道编号
F_PH = [F_PH, mr_pipe(b)* tao_PR_load(b,:) == mr_pipe(b) * tao_NR_node(i,:)];
F_PH = [F_PH, tao_PS_src(b,:) == tao_NS_node(i,:)];
elseif num_temp2 == 0 % 以节点i首端的管道数为2
b = S_pipe_F(i,:); % 管道编号
F_PH = [F_PH, mr_pipe(b(1))* tao_PR_load(b(1),:) + mr_pipe(b(2))* tao_PR_load(b(2),:) == (mr_pipe(b(1))+mr_pipe(b(2))) * tao_NR_node(i,:)];
F_PH = [F_PH, tao_PS_src(b(1),:) == tao_NS_node(i,:)];
F_PH = [F_PH, tao_PS_src(b(2),:) == tao_NS_node(i,:)];
else
% 以节点i末端的管道数为0 末节点
end
end
for i = 1:N_pipe
% 温度变化方程
F_PH = [F_PH, tao_PS_load(i,:) == (tao_PS_src(i,:) - (tao_am-tao_K)).*exp(-lamada_pipe(i)*L_pipe(i)/(cp_w*ms_pipe(i))) + (tao_am-tao_K)];
F_PH = [F_PH, tao_PR_load(i,:) == (tao_PR_src(i,:) - (tao_am-tao_K)).*exp(-lamada_pipe(i)*L_pipe(i)/(cp_w*mr_pipe(i))) + (tao_am-tao_K)];
end
最新发布
扫一扫
816
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
