clc
clear
%% DIES:谈判破裂点模型
% 电系统:风机、光伏、燃气轮机
% 热系统:余热锅炉、蓄热罐
% 冷系统:吸收式制冷机、蓄冰槽
% 气系统:P2G机组、CCS机组
%% 导入电/热/冷负荷
P_Load_e0 = [8800,8600,8500,8500,8600,9250,12000,12100,12000,11000,9500,9000,8500,8500,8750,9400,10600,11900,12100,12000,8900,8800,8750,8500];
P_Load_h0 = [5681,5800,5854,6312,6512,6706,8258,7312,7823,8234,8432,7680,7600,7480,7520,7480,7320,6800,6012,6523,6125,4902,4612,4426];
P_Load_c0 = [1987,2053,2234,2342,2553,2712,3915,4124,4351,4612,4612,4203,4200,4350,4500,4305,4245,4050,3750,3000,2250,2000,1750,1500];
% 电力市场分时价格
C_buy_e = [0.30,0.32,0.33,0.35,0.35,0.34,0.36,1.16,1.19,1.14,1.16,1.15,0.62,0.73,0.70,0.64,0.71,0.65,1.17,1.12,1.11,1.12,0.34,0.28];
% 天然气价格数据(元/m3)
C_buy_g = [2.58*ones(1,8),2.18*ones(1,14),2.58*ones(1,1),2.18*ones(1,1)];
%% 计算风机、光伏出力功率
% 光照辐射度
gz = [0;0;0;0;0.0708300000000000;0.170690000000000;0.325740000000000;0.380270000000000;0.340240000000000;0.453630000000000;0.855530000000000;0.742260000000000;0.666500000000000;0.699100000000000;0.535630000000000;0.227290000000000;0.179330000000000;0.0335400000000000;0;0;0;0;0;0];
for t=1:24
Pstc=0.26; Gstc=1; n=0.9645; PV_k=5000;
P_pv_max(t)=(Pstc*gz(t)/Gstc*n)*PV_k;
end
% 风速大小
fs=[10.7391000000000;11.6160000000000;9.85170000000000;9.98270000000000;10.5166000000000;5.91910000000000;4.21130000000000;4.73580000000000;2.89860000000000;5.54320000000000;4.93960000000000;2.99680000000000;3.64480000000000;5.17060000000000;2.45990000000000;5.43100000000000;6.83090000000000;10.5836000000000;10.8847000000000;12.3629000000000;14.4617000000000;12.1296000000000;10.5531000000000;12.8199000000000];
for t=1:24
Vci=3; Vco=25; Vr=11; Pwt_r=10; WT_k=50;
if all(fs(t)<Vci)||all(fs(t)>Vco)
P_wt_max(t)=0*WT_k;
elseif all(fs(t)>=Vci)&&all(fs(t)<=Vr)
P_wt_max(t)=[Pwt_r*(fs(t)-Vci)/(Vr-Vci)]*WT_k;
elseif all(fs(t)>Vr)&&all(fs(t)<=Vco)
P_wt_max(t)=10*WT_k;
end
end
%% 决策变量初始化
% 需求响应初始化
P_Load_e = sdpvar(1,24,'full'); % 经过需求响应后实际的电负荷
P_Load_h = sdpvar(1,24,'full'); % 经过需求响应后实际的热负荷
P_Load_c = sdpvar(1,24,'full'); % 经过需求响应后实际的热负荷
P_e_cut = sdpvar(1,24,'full'); % 可削减电负荷
P_e_tra = sdpvar(1,24,'full'); % 可转移电负荷
P_h_DR = sdpvar(1,24,'full'); % 可削减热负荷
P_c_DR = sdpvar(1,24,'full'); % 可削减冷负荷
%% 设备变量初始化
% 蓄电池
SOC_ESS = sdpvar(1,24,'full'); % 储电设备的储电余量
P_ESS_cha = sdpvar(1,24,'full'); % 储电设备的充电功率
P_ESS_dis = sdpvar(1,24,'full'); % 储电设备的放电功率
B_ESS_cha = binvar(1,24,'full'); % 储电设备的放电状态位,取1时为放电,0为未放电
B_ESS_dis = binvar(1,24,'full'); % 储电设备的充电状态位,取1时为充电,0为未充电
% 蓄热罐
SOC_HSS = sdpvar(1,24,'full'); % 储热设备的储电余量
P_HSS_cha = sdpvar(1,24,'full'); % 储热设备的充热功率
P_HSS_dis = sdpvar(1,24,'full'); % 储热设备的放热功率
B_HSS_cha = binvar(1,24,'full'); % 储热设备的放热状态位,取1时为放热,0为未放热
B_HSS_dis = binvar(1,24,'full'); % 储热设备的充热状态位,取1时为充热,0为未充热
% 蓄冰槽
SOC_CSS = sdpvar(1,24,'full'); % 储冷设备的储电余量
P_CSS_cha = sdpvar(1,24,'full'); % 储冷设备的充冷功率
P_CSS_dis = sdpvar(1,24,'full'); % 储冷设备的放冷功率
B_CSS_cha = binvar(1,24,'full'); % 储冷设备的放冷状态位,取1时为放冷,0为未放冷
B_CSS_dis = binvar(1,24,'full'); % 储冷设备的充冷状态位,取1时为充冷,0为未充冷
% 风机、光伏
P_wt_e = sdpvar(1,24,'full'); % 风力的实际出力值
P_pv_e = sdpvar(1,24,'full'); % 光伏的实际出力值
% CCHP系统:
% 燃气轮机(GT)
P_GT_e = sdpvar(1,24,'full'); % 燃气轮机的实际电出力值
P_GT_fh = sdpvar(1,24,'full'); % 燃气轮机的废热
% 余热锅炉(WB)
P_WB_h = sdpvar(1,24,'full'); % 余热锅炉的实际热出力值
P_fh_WB = sdpvar(1,24,'full'); % 余热锅炉吸收的废热
% 吸收式制冷机(AC)
P_AC_c = sdpvar(1,24,'full'); % 吸收式制冷机的实际冷出力值
P_fh_AC=sdpvar(1,24,'full'); % 吸收式制冷机吸收的废热
% P2G电转气设备(电解槽EL→甲烷反应器MR)
% EL:电解槽(电解槽EL将电能转化为氢能,氢能经由甲烷反应器MR,进一步转化为天然气)
P_e_EL = sdpvar(1,24,'full'); % 燃气轮机供给电解槽EL的电出力
P_EL_H = sdpvar(1,24,'full'); % 电解槽生产H2
% MR:甲烷反应器
P_C_MR = sdpvar(1,24,'full'); % 输入MR设备的CO2
P_H_MR = sdpvar(1,24,'full'); % 输入MR设备的H2
P_MR_g = sdpvar(1,24,'full'); % MR设备生产天然气
P_C_p2g = sdpvar(1,24); % P2G所用的二氧化碳量
% CCS碳捕集机组
P_e_CCS = sdpvar(1,24,'full'); % 碳捕集机组产生的总能耗
E_total_co2 = sdpvar(1,24,'full'); % 机组捕获的总碳排放
P_CCS_B = sdpvar(1,24,'full'); % 机组运行能耗
P_CCS_G = 2000; % 碳捕集固定能耗为2000kwh
Soc_CO2 = sdpvar(1,24,'full'); % CCS的碳捕集量/P2G所用的二氧化碳量
% 电力市场相关
P_buy_e = sdpvar(1,24,'full'); % 电力市场购电量
P_buy_g = sdpvar(1,24,'full'); % 购买的天然气量
%% 设备数据
% CCHP:燃气轮机、余热锅炉、吸收式制冷机
k_eff_CCHP_e = 0.38; k_eff_CCHP_h = 0.46;k_eff_CCHP_c = 0.52;
% 燃气轮机(GT)
%(1)燃气最大出力
P_GT_max=8000;
%(2)燃气轮机热损率
rs_GT=0.03;
% 余热锅炉(WGB)
P_h_WB_max=7000;
% 吸收式制冷机(AC)
P_c_AC_max=7000;
% 储能电池
k_eff_ebess_cha=0.95; % 蓄电池充电效率
k_eff_ebess_dis=0.95; % 蓄电池放电效率
SOC_ESS(1)=0.5;
SOC_Emax=0.9;
P_Ebess_cha_max=5800;
P_Ebess_dis_max=5800;
Ebess=6000;
% 蓄热罐
k_eff_hbess_cha=0.85; % 蓄电池充电效率
k_eff_hbess_dis=0.85; % 蓄电池放电效率
SOC_HSS(1)=0.5;
P_Hbess_cha_max=3500;
P_Hbess_dis_max=3500;
Hbess=4000;
% 蓄冷槽
k_eff_cbess_cha=0.85; % 蓄冰槽充冷效率
k_eff_cbess_dis=0.85; % 蓄冰罐放冷效率
SOC_CSS(1)=0.5;
P_Cbess_cha_max=3500;
P_Cbess_dis_max=3500;
Cbess=4000;
% P2G电转气装置
ngas_G=0.35; %气转电效率
ngas_h=0.9; %气转热效率
ngas_c=0.9; %气转冷效率
%8.CCS碳捕集装备
P_CCS_G=2000; % 固定能耗2000kwh
%% 导入约束条件
Cons=[];
% 需求响应约束条件
for t=1:24
Cons=[Cons,
P_Load_e(t)==P_Load_e0(t)-P_e_cut(t)-P_e_tra(t), % 电负荷功率平衡约束
P_Load_h(t)==P_Load_h0(t)-P_h_DR(t), % 热负荷功率平衡约束
P_Load_c(t)==P_Load_c0(t)-P_c_DR(t), % 冷负荷功率平衡约束
0 <=P_e_cut(t)<= 0.05*P_Load_e0(t), % 可削减电功率上下 限约束
-0.1*P_Load_e0(t)<=P_e_tra(t)<=0.1*P_Load_e0(t), % 可转移电功率上下限约束
-0.1*P_Load_h0(t)<=P_h_DR(t) <=0.1*P_Load_h0(t), % 可削减热功率上下限约束
-0.1*P_Load_c0(t)<=P_c_DR(t) <=0.1*P_Load_c0(t)]; % 可削减热功率上下限约束
end
Cons=[Cons,sum(P_e_tra)==0,]; % 转移的电负荷总量为0约束
Cons=[Cons,sum(P_h_DR)==0,]; % 转移的热负荷总量为0约束
Cons=[Cons,sum(P_c_DR)==0,]; % 转移的冷负荷总量为0约束
% 设备约束条件
% CCHP(冷热电三联供)
% (1)燃气轮机GT:
% 最大出力约束
for t=1:24
Cons=[Cons,0<=P_GT_e(t)<=P_GT_max,
P_GT_e(t)==Gas_g_GT(t)*k_eff_CCHP_e*9.88,
P_GT_fh(t)==P_GT_e(t)*(1-k_eff_CCHP_e-rs_GT)/k_eff_CCHP_e,
P_GT_e(t)==P_e_EL(t)+P_e_CCS(t)+P_GT_e(t),
0<=P_e_EL(t),0<=P_e_CCS(t),0<=P_GT_e(t)];
end
% 爬坡约束
for t=1:23
Cons=[Cons,-0.2*2000<=P_GT_e(t+1)-P_GT_e(t)<=0.2*2000];
end
%(2)余热锅炉WB:
% 出力约束
for t=1:24
Cons=[Cons,P_WB_h(t)==P_fh_WB(t)*k_eff_CCHP_h,0<=P_WB_h(t)<=P_h_WB_max,]
end
% 爬坡约束
Cons=[Cons,-0.2*1000 <=P_WB_h(2:24)-P_WB_h(1:23)<=0.2*1000];
%(3)吸收式制冷机AC:
for t=1:24
Cons=[Cons, 0<=P_c_AC(t)<=P_c_AC_max,
P_c_AC(t)==P_fh_AC(t)*k_eff_CCHP_c];
end
Cons=[Cons,-0.2*2000 <=P_c_AC(2:24)-P_c_AC(1:23)<=0.2*2000];
%(4)CCHP废气约束
for t=1:24
Cons=[Cons,P_fh_WB(t)+P_fh_AC(t)<=P_GT_fh(t)];
end
%% 蓄电池约束条件
% 1)蓄电池的约束条件:蓄电池功率约束在[min,max]内,且统一时刻只允许充电、放电、不充不放的三种状态中的一种
% 储能电站的充放电功率约束,Big-M法进行线性化处理
M=1000; % 这里的M是个很大的数
for t=1:24
Cons=[Cons,
0<=P_ESS_cha(t)<=P_Ebess_cha_max,
0<=P_ESS_cha(t)<=B_ESS_cha(t)*M,
0<=P_ESS_dis(t)<=P_Ebess_dis_max,
0<=P_ESS_dis(t)<=B_ESS_dis(t)*M,
B_ESS_dis(t)+B_ESS_cha(t)<=1];
end
% 2)蓄电池的约束条件:蓄电池功率与容量耦合约束,始末时刻负荷电量相等约束,与容量被限制在[min,max]
for i=1:23
Cons=[Cons,SOC_ESS(i+1)==SOC_ESS(i)+(P_ESS_cha(i)*k_eff_ebess_cha-P_ESS_dis(i)/k_eff_ebess_dis)/Ebess];
end
Cons=[Cons,SOC_ESS(24)==SOC_ESS(1),SOC_ESS>=0.2,SOC_ESS<=0.9];
%% 蓄热罐约束条件
% 1)蓄热罐的约束条件:蓄热罐功率约束在[min,max]内,且统一时刻只允许充热、放热、不充不放的三种状态中的一种
M=1000; % 这里的M是个很大的数
for t=1:24
Cons=[Cons,
0<=P_HSS_cha(t)<=P_Hbess_cha_max,
0<=P_HSS_cha(t)<=B_HSS_cha(t)*M,
0<=P_HSS_dis(t)<=P_Hbess_dis_max,
0<=P_HSS_dis(t)<=B_HSS_dis(t)*M,
B_HSS_dis(t)+B_HSS_cha(t)<=1];
end
% 2)蓄热罐的约束条件:蓄热罐功率与容量耦合约束,始末时刻负荷电量相等约束,与容量被限制在[min,max]
for t=1:23
Cons=[Cons,SOC_HSS(t+1)==SOC_HSS(t)+(P_HSS_cha(t)*k_eff_hbess_cha-P_HSS_dis(t)/k_eff_hbess_dis)/Hbess];
end
Cons=[Cons,SOC_HSS(24)==SOC_HSS(1)];
Cons=[SOC_HSS>=0.2,SOC_HSS<=0.9];
%% 蓄冰槽约束条件
% 1)蓄冰槽的约束条件:蓄电池功率约束在[min,max]内,且统一时刻只允许充电、放电、不充不放的三种状态中的一种
M=1000; % 这里的M是个很大的数
for t=1:24
Cons=[Cons,
0<=P_CSS_cha(t)<=P_Cbess_cha_max,
0<=P_CSS_cha(t)<=B_CSS_cha(t)*M,
0<=P_CSS_dis(t)<=P_Cbess_dis_max,
0<=P_CSS_dis(t)<=B_CSS_dis(t)*M,
B_CSS_dis(t)+B_CSS_cha(t)<=1];
end
% 2)蓄冰槽的约束条件:蓄电池功率与容量耦合约束,始末时刻负荷电量相等约束,与容量被限制在[min,max]
for i=1:23
Cons=[Cons,SOC_CSS(i+1)==SOC_CSS(i)+(P_CSS_cha(i)*k_eff_cbess_cha-P_CSS_dis(i)/k_eff_cbess_dis)/Cbess];
end
Cons=[Cons,SOC_CSS(24)==SOC_CSS(1)];
Cons=[Cons,SOC_CSS>=0.2,SOC_CSS<=0.9];
%% GT—P2G-CCS约束
for t=1:24
Cons=[Cons,
P_EL_H(t)==P_e_EL(t)*0.85, % 电解槽EL在t时刻的制氢气功率=燃气轮机向电解槽供给的电出力*能量转化效率0.85
0-P_e_EL(t)-P_e_CCS(t)<=P_GT_e(t)<=6000-P_e_EL(t)-P_e_CCS(t),
0<=P_e_CCS(t)<=3000,
0<=P_e_EL(t)<=3000,
0<=P_GT_e<=6000,
max((0-0.15*P_WB_h(t)-P_e_EL(t)-P_e_CCS(t)),(0.85*(P_WB_h(t))-P_e_CCS(t)-P_e_EL(t)))<=P_GT_e(t)<=6000-0.20*P_WB_h(t)-P_e_EL(t)-P_e_CCS(t), % CCHP的热电耦合约束
max((0-0.15*P_WB_h(t)),(0.85*(P_WB_h(t)-50)-600-1200))<=P_GT_e(t)<=6000-0.20*P_WB_h(t)-0-0, % 考虑P2G和CCS后的CCHP的热电耦合约束
(0.55/(1+0.5*1.02))*max((0-0.15*P_WB_h(t)-P_GT_e(t)),(0.85*(P_WB_h(t)-50)-P_GT_e(t)))<=P_MR_g(t)<=(0.55/(1+0.5*1.02))*(6000-0.20*P_WB_h(t)-P_GT_e(t)), % 产气功率上下限约束
-1000<=(P_GT_e(2:24)+P_e_EL(2:24)+P_e_CCS(2:24))-(P_GT_e(1:23)+P_e_EL(1:23)+P_e_CCS(1:23))<=1000, % CCHP的爬坡约束
P_MR_g(t)==0.55*P_e_EL(t), % P2G产气功率与耗电量约束
E_total_co2(t)==P_e_CCS(t)/0.55, % CCS的耗电量与碳捕集量约束
P_C_p2g==1.02*P_e_EL(t), % P2G运行所需要的二氧化碳量与电功率约束
Soc_CO2(1) == 2000+0.97*E_total_co2(1)-P_C_p2g(1)/0.97 ,
Soc_CO2(24) == 2000,
1000<=Soc_CO2(t)<= 3000];
end
for t=2:24
Cons=[Cons,Soc_CO2(t) == Soc_CO2(t-1) + 0.97*E_total_co2(t)-P_C_p2g(t)/0.97,]; % 储电设备容量变化约束
end
% 碳捕集相关
Cons=[Cons,0<=E_total_co2<=0.55*(P_e_EL+P_e_CCS+P_GT_e+0.15*P_WB_h)];
% 天然气交易
for t=1:24
Cons=[Cons, % GB耗气量约束,CCHP发电效率=CCHPt时刻天然气消耗量*燃气轮机发电效率*天然气燃烧热值转换系数
0<=Gas_all(t)<=8000,
0<=Gas_g_GT(t)<=6000]; % 总耗气量约束
end
% 电力交易市场
for t=1:24
Cons=[Cons,P_buy_e(t)>=0,P_buy_e(t)<=8000];
end
%% 阶梯碳交易成本
% 阶梯碳交易成本(分段线性化)
E_v=sdpvar(1,7); % 每段区间内的长度,分为5段,每段长度是2000
lamda=0.250; % 碳交易基价
Cons=[Cons,
E==sum(E_v), % 总长度等于E
0<=E_v(1:6)<=2000, % 除了最后一段,每段区间长度小于等于2000
0<=E_v(7)];
% 碳交易成本
C_CO2=0;
for v=1:7
C_CO2=C_CO2+(lamda+(v-1)*0.25*lamda)*E_v(v);
end
%% 电平衡约束、热平衡约束、冷平衡约束
for t=1:24
Cons=[Cons,P_wt_e(t)+P_pv_e(t)+P_GT_e(t)+P_ESS_dis(t)+P_buy_e(t)==P_Load_e(t)+P_ESS_cha(t), % 电平衡
P_WB_h(t)+P_HSS_dis(t)==P_Load_h(t)+P_HSS_cha(t) , % 热平衡
P_c_AC(t)+P_CSS_dis(t)==P_Load_c(t)+P_CSS_cha(t), % 冷平衡
P_buy_g(t)+P_MR_g(t)== Gas_all(t),
P_GT_e(t)==P_e_EL(t)+P_e_CCS(t)+P_GT_e(t)]
end
%% 设备运行成本约束
Cel=sdpvar(1,1,'full');
Cmr=sdpvar(1,1,'full');
Cccs=sdpvar(1,1,'full');
for t=1:24
Cons=[Cons,Cel==(sum(P_e_EL)+sum(P_e_CCS))*0.234,
Cccs==sum(P_e_CCS)*0.138;
P_GT_e(t)==P_e_EL(t)+P_e_CCS(t)+P_GT_e(t),
0<=P_GT_e(t)];
end
%% 目标函数(购能成本+机组运维成本+需求侧响应成本+环境成本)
Cost_f1 = sum(P_buy_e)*C_buy_e*ones(1,24)';
Cost_f2 = sum(P_buy_g)*C_buy_g*ones(1,24)'
Cost_f3 = sum(P_GT_e)*0.04+sum(P_WB_h)*0.025+sum(P_c_AC)*0.025+sum(P_wt_e)*0.016+sum(P_pv_e)*0.018+sum(P_HSS_cha)*0.016+sum(P_HSS_dis)*0.016+sum(P_CSS_cha)*0.017+sum(P_CSS_dis)*0.017+Cel+Cmr+Cccs+sum(P_ESS_cha)*0.018+sum(P_ESS_dis)*0.018
Cost_f4 = 0.3*sum(abs(P_e_tra))+0.3*sum(abs(P_e_cut))+ +0.016*sum(abs(P_h_DR))+0.016*sum(abs(P_c_DR));
Cost_f5 = C_CO2;
F_all = Cost_f1+Cost_f2+Cost_f3+Cost_f4+Cost_f5;
%% 求解器相关配置
ops=sdpsettings('solver','cplex','verbose',0,'showprogress',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
% 进行求解计算:
sol=optimize(Cons,Cost_all,ops);
if sol.problem==0
disp('>>求解成功!')
else
disp('>>求解失败,失败原因:')
disp(sol.info)
end
% 输出求解结果:
F_all = value(F_all);
display(['>>>园区综合能源系统的全局最优规划值为 : ', num2str(F_all)]);
Cost_f1=value(Cost_f1);
display(['1、园区综合能源系统的购电成本为 : ', num2str(Cost_f1)]);
Cost_f2=value(Cost_f2);
display(['2、园区综合能源系统的购气成本为 : ', num2str(Cost_f2)]);
Cost_f3=value(Cost_f3);
display(['3、园区综合能源系统的运维成本为 : ', num2str(Cost_f3)]);
Cost_f4=value(Cost_f4);
display(['4、园区综合能源系统的需求成本为 : ', num2str(Cost_f4)]);
Cost_f5=value(Cost_f5);
display(['5、园区综合能源系统的环境成本为 : ', num2str(Cost_f5)]);
最新发布