请你检查以下代码:clc
clear
%% 综合能源系统
% 电系统:风机、光伏、CCHP机组(燃气轮机、余热锅炉、吸收式制冷机)
% 热系统:CCHP机组、燃气锅炉、蓄热罐
% 冷系统:CCHP机组、蓄冰槽
%% 导入电、 热、冷负荷
P_Load_e0 = [6150,5347,4737,4886,4960,4830,4998,6056,7004,9140,9938,10645,11331,11443,11499,11292,11212,11595,10822,9441,9092,7794,7174,6777];
P_Load_h0 = [5986,5952,6033,5722,5977,6060,5872,5815,5562,5424,5320,5124,4962,4986,4962,5019,5017,5179,5181,5305,5454,5643,5964,6033];
P_Load_c0 = [4626,4725,4947,4822,4977,5060,4872,4815,4662,4484,4290,4184,3992,3986,3923,4089,4017,4179,4261,4383,4458,4724,5012,5423];
%% 计算风机、光伏出力功率
% 光照辐射度
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];
P_pv_max = zeros(1,24);
for t=1:24
P_stc = 0.26; G_stc = 1; n = 0.9645; PV_k = 50;
P_pv_max = (P_stc * gz' / G_stc * 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];
P_wt_max = zeros(1,24);
for t=1:24
V_ci = 3; V_co = 25; V_ra = 11; Pwt_r = 10; WT_k = 50;
if fs(t) < V_ci || fs(t) > V_co
P_wt_max(t) = 0;
elseif fs(t) >= V_ci && fs(t) <= V_ra
P_wt_max(t) = Pwt_r * (fs(t) - V_ci) / (V_ra - V_ci) * WT_k;
else
P_wt_max(t) = 10 * WT_k;
end
end
%% 外部市场数据:常规能源市场、天然气市场
% 常规能源市场价格:
C_cg_ce = [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];
% 绿电能源市场价格:
C_ld_ge = [0.42,0.45,0.43,0.42,0.42,0.38,0.42,1.26,1.29,1.24,1.29,1.18,0.72,0.78,0.76,0.79,0.74,0.72,1.28,1.25,1.25,1.21,0.46,0.34];
% 天然气价格数据(元/m3),辽宁省是2.33元/m3
C_P_gas = [2.58*ones(1,8),2.18*ones(1,14),2.58*ones(1,1),2.18*ones(1,1)];
%% 决策变量初始化
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_h = sdpvar(1,24,'full'); % 燃气轮机的发出废热
% 余热锅炉(HB)
P_HB_h = sdpvar(1,24,'full'); % 余热锅炉的实际热出力值
P_fh_HB = sdpvar(1,24,'full'); % 余热锅炉吸收的废热
% 吸收式制冷机(AC)
P_AC_c = sdpvar(1,24,'full'); % 吸收式制冷机的实际冷出力值
P_fh_AC = sdpvar(1,24,'full'); % 吸收式制冷机吸收的废热
% 燃气锅炉(GB)
P_GB_h = sdpvar(1,24,'full');
%% 常规能源市场交易(绿电能源、常规能源)
P_e_net = sdpvar(1,24,'full'); % 电力市场净购电量=购电量-售电量
P_e_buy = sdpvar(1,24,'full'); % 电力市场购电量
Temp_net = binvar(1,24,'full'); % 电网购|售电标志
% 绿电能源日前市场
P_ld_buy = sdpvar(1,24,'full'); % 向外电网的购买的绿电功率
B_ld_buy = binvar(1,24,'full'); % 向外电网的购买的绿电标志,取1则买绿电,取0则卖绿电
% 常规能源日前市场
P_cg_buy = sdpvar(1,24,'full'); % 向外电网的购买的常规能源发电功率
B_cg_buy = binvar(1,24,'full'); % 向外电网的购买的常规能源发电标志,取1则买常规能源,取0则卖常规能源
%% 天然气交易市场(天然气)
Gas_all = sdpvar(1,24,'full'); % 系统的总耗气量
G_g_GT = sdpvar(1,24,'full'); % 燃气轮机的天然气耗气量
G_g_GB = sdpvar(1,24,'full'); % 燃气锅炉的天然气耗气量
G_g_buy = sdpvar(1,24,'full'); % 综合能源系统1购买天然气量
%% 设备数据参数
% CCHP机组:燃气轮机、余热锅炉、吸收式制冷机
k_eff_CCHP_e = 0.38;
k_eff_CCHP_h = 0.46;
k_eff_CCHP_c = 0.52;
% 燃气轮机(GT)
P_GT_max = 8000; % 燃气最大出力
rs_GT = 0.03; % 燃气轮机热损率
% 余热锅炉(HB)
P_HB_h_max = 8000;
% 吸收式制冷机(AC)
P_AC_c_max = 6000;
% 燃气锅炉(GB)
P_GB_max = 5000; % 燃气锅炉最大出力
k_eff_GB = 0.85; % 燃气锅炉效率
% 蓄电池
k_eff_ess_cha = 0.95; % 蓄电池充电效率
k_eff_ess_dis = 0.85; % 蓄电池放电效率
SOC_ESS(1) = 0.5;
SOC_Emax = 0.9;
P_Ess_cha_max = 4800;
P_Ess_dis_max = 4800;
Ebess = 5000;
% 蓄热罐
k_eff_hss_cha = 0.85; % 蓄热罐充热效率
k_eff_hss_dis = 0.85; % 蓄热罐放热效率
SOC_HSS(1) = 0.5;
SOC_Hmax = 0.9;
P_Hss_cha_max = 4500;
P_Hss_dis_max = 4500;
Hbess = 5000;
% 蓄冰槽
k_eff_css_cha = 0.85; % 蓄冰槽充冷效率
k_eff_css_dis = 0.85; % 蓄冰罐放冷效率
SOC_Cmax = 0.9;
SOC_CSS(1) = 0.5;
P_Css_cha_max = 4500;
P_Css_dis_max = 4500;
Cbess = 5000;
%% 约束条件
Cos = [];
% 电/ 热/ 冷负荷需求响应部分
for t=1:24
Cos=[Cos,
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
Cos=[Cos,sum(P_e_tra)==0,]; % 转移的电负荷总量为0约束
Cos=[Cos,sum(P_h_DR)==0,]; % 转移的热负荷总量为0约束
Cos=[Cos,sum(P_c_DR)==0,]; % 转移的冷负荷总量为0约束
% ========== CCHP机组(冷热电三联供) ==========
% (1)燃气轮机GT:
% 1)最大出力约束
for t=1:24
Cos=[Cos,0<=P_GT_e(t)<=P_GT_max,
P_GT_e(t)==G_g_GT(t)*k_eff_CCHP_e*9.88,
P_GT_h(t)==P_GT_e(t)*(1-k_eff_CCHP_e-rs_GT)/k_eff_CCHP_e];
end
% 2)爬坡约束
for t=1:23
Cos=[Cos,-0.2*2000<=P_GT_e(t+1)-P_GT_e(t)<=0.2*2000];
end
%(2)余热锅炉HB:
% 1)出力约束
for t=1:24
Cos=[Cos,P_HB_h(t)==P_fh_HB(t)*k_eff_CCHP_h,
0<=P_HB_h(t)<=P_HB_h_max,]
end
% 2)爬坡约束
Cos=[Cos,-0.2*1000 <=P_HB_h(2:24)-P_HB_h(1:23)<=0.2*1000];
%(3)吸收式制冷机AC
for t=1:24
Cos=[Cos, 0<=P_AC_c(t)<=P_AC_c_max,
P_AC_c(t)==P_fh_AC(t)*k_eff_CCHP_c];
end
% (4)CCHP机组废气约束
for t=1:24
Cos=[Cos,P_fh_HB(t)+0.56*P_fh_AC(t)<=P_GT_h(t)];
end
% ========== 燃气锅炉GB(热) ==========
Cos=[Cos,P_GB_h(t)==G_g_GB(t)*k_eff_GB*9.88,
0<=P_GB_h(t)<=P_GB_max,
-0.2*1200 <=P_GB_h(2:24)-P_GB_h(1:23)<=0.2*1200];
% ========== 蓄电池机组(电) ==========
% 1)蓄电池的约束条件:蓄电池功率约束在[min,max]内,且统一时刻只允许充电、放电、不充不放的三种状态中的一种
Cos=[Cos,P_ESS_cha>=0,P_ESS_cha<=P_Ess_cha_max*B_ESS_cha,
0<=P_ESS_dis,P_ESS_dis<=P_Ess_dis_max*B_ESS_dis,
B_ESS_cha+B_ESS_dis<=1]; %充放电状态唯一
% 2)蓄电池的约束条件:蓄电池功率与容量耦合约束,始末时刻负荷电量相等约束,与容量被限制在[min,max]
for i=1:23
Cos=[Cos,SOC_ESS(i+1)==SOC_ESS(i)+(P_ESS_cha(i)*k_eff_ess_cha-P_ESS_dis(i)/k_eff_ess_dis)/Ebess];
end
Cos=[Cos,SOC_ESS(24)==SOC_ESS(1),SOC_ESS>=0.2,SOC_ESS<=0.9];
% ========== 蓄热罐机组(热) ==========
% 1)蓄热罐的约束条件:蓄热罐功率约束在[min,max]内,且统一时刻只允许充热、放热、不充不放的三种状态中的一种
Cos=[Cos,P_HSS_cha>=0,P_HSS_cha<=P_Hss_cha_max*B_HSS_cha,
P_HSS_dis>=0,P_HSS_dis<=P_Hss_dis_max*B_HSS_dis,
B_HSS_cha+B_HSS_dis<=1];% 充放电状态唯一
% 2)蓄热罐的约束条件:蓄热罐功率与容量耦合约束,始末时刻负荷热能相等约束,与容量被限制在[min,max]
for i=1:23
Cos=[Cos,SOC_HSS(i+1)==SOC_HSS(i)+(P_HSS_cha(i)*k_eff_hss_cha-P_HSS_dis(i)/k_eff_hss_dis)/Hbess];
end
Cos=[Cos,SOC_HSS(24)==SOC_HSS(1),SOC_HSS>=0.2,SOC_HSS<=0.9];
% ========== 蓄冰槽机组(冷) ==========
% 1)蓄冰槽的约束条件:蓄冷罐功率约束在[min,max]内,且统一时刻只允许充冷、放冷、不充不放的三种状态中的一种
Cos=[Cos,P_CSS_cha>=0,P_CSS_cha<=P_Css_cha_max*B_CSS_cha,
P_CSS_dis>=0,P_CSS_dis<=P_Css_dis_max*B_CSS_dis,
B_CSS_cha+B_CSS_dis<=1];% 充放电状态唯一
% 2)蓄冰槽的约束条件:蓄冷罐功率与容量耦合约束,始末时刻负荷冷能相等约束,与容量被限制在[min,max]
for i=1:23
Cos=[Cos,SOC_CSS(i+1)==SOC_CSS(i)+(P_CSS_cha(i)*k_eff_css_cha-P_CSS_dis(i)/k_eff_css_dis)/Cbess];
end
Cos=[Cos,SOC_CSS(24)==SOC_CSS(1)];
Cos=[Cos,SOC_CSS>=0.2,SOC_CSS<=0.9];
% 天然气购买约束条件
for t=1:24
Cos=[Cos, % GB耗气量约束,CCHP发电效率=CCHPt时刻天然气消耗量*燃气轮机发电效率*天然气燃烧热值转换系数
Gas_all(t)==G_g_GT(t)+G_g_GB(t), % 总耗气量约束
0<=Gas_all(t)<=12000,
0<=G_g_GT(t)<=6000,
0<=G_g_GB(t)<=6000];
end
% 光伏、风机约束条件
for t=1:24
Cos=[Cos,
0<=P_wt_e(t)<=P_wt_max(t),
0<=P_pv_e(t)<=P_pv_max(t)];
end
%% 电力交易市场(常规能源市场)
for t=1:24
Cos=[Cos,P_e_buy(t)==0.78*P_ld_buy(t)+P_cg_buy(t),P_e_buy(t)>=0,P_e_buy(t)<=10000];
end
% 绿电日前市场预测购、售电功率
for t=1:24
Cos=[Cos,200<=P_ld_buy(t),P_ld_buy(t)<=8000*B_ld_buy(t)];
end
% 常规能源日前市场预测购、售电功率
for t=1:24
Cos=[Cos,0<=P_cg_buy(t),P_cg_buy(t)<=8000*B_cg_buy(t)];
end
% 电力净购买量,若P_E_net>0,则不需要买电;若P_E_net<0,需要买电
for t=1:24
Cos=[Cos,P_e_net(t)==P_wt_e(t)+P_pv_e(t)+P_GT_e(t)+P_ESS_dis(t)-P_ESS_cha(t)-P_Load_e(t)];
end
% 电力购售量
for t=1:24
Cos=[Cos,implies(Temp_net(1,t),[P_e_net(1,t)>=0,P_e_buy(1,t)==0]),
implies(Temp_net(1,t),[P_e_net(1,t)<=0,P_e_buy(1,t)==-P_e_net(1,t)])];
end
%% 阶梯碳交易成本
% 碳排放权配额模型
E_e_buy=0.728*sum(P_cg_buy); % 购电配额
E_CCHP=(sum(P_GT_e)/(0.375*9.88)+(sum(P_GB_h)/(0.82*9.88)))*0.385; % CCHP、GB配额
E_IES=E_e_buy+E_CCHP; % IES总碳排放配额
E1=0.728*(P_cg_buy)+0.102*3.6*((P_GT_e)/(0.375*9.88))+0.385*(P_GB_h)/(0.82*9.88);
% 实际碳排放模型
E_e_buy_a=1.08*sum(P_cg_buy);
E_CCHP_a=0.934*(sum(P_GT_e)/(0.375*9.88)+sum(P_GB_h)/(0.82*9.88));
E_IES_a=E_e_buy_a+E_CCHP_a;
E=E_IES_a-E_IES; % 实际IES总碳排放
% 阶梯碳交易成本(分段线性化)
E_v=sdpvar(1,7) % 每段区间内的长度,分为5段,每段长度是2000
lamda=0.250; % 碳交易基价
Cos=[Cos,
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
Cos=[Cos,P_wt_e(t)+P_pv_e(t)+P_GT_e(t)+P_ESS_dis(t)+P_e_buy(t)==P_Load_e(t)+P_ESS_cha(t), % 电功率平衡
P_HB_h(t)+P_GB_h(t)+P_HSS_dis(t)==P_Load_h(t)+P_HSS_cha(t), % 热功率平衡
P_AC_c(t)+P_CSS_dis(t)==P_Load_c(t)+P_CSS_cha(t), % 冷功率平衡
G_g_buy(t)== Gas_all(t)]
end
%% 目标函数(购能成本+运行成本+需求侧成本)
Cost_f1 = sum(P_ld_buy)*C_ld_ge*ones(1,24)'+sum(P_cg_buy)*C_cg_ce*ones(1,24)'+sum(Gas_all)*C_P_gas*ones(1,24)';
Cost_f2 = sum(P_GT_e)*0.04+sum(P_HB_h)*0.025+sum(P_AC_c)*0.025+0.024*sum(P_GB_h)+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+sum(P_ESS_cha)*0.018+sum(P_ESS_dis)*0.018;
Cost_f3 = 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_f4 = C_CO2;
Cost_all = Cost_f1+Cost_f2+Cost_f3+Cost_f4;
%% 求解器相关配置
ops=sdpsettings('solver','cplex','verbose',0,'showprogress',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
%% 进行求解计算
sol=optimize(Cos,Cost_all,ops);
if sol.problem==0
disp('>>求解成功!')
else
disp('>>求解失败,失败原因:')
disp(sol.info)
end
Cost_all = value(Cost_all);
display(['>>>园区综合能源系统1独立运行成本的全局最优规划值为 : ', num2str(Cost_all)]);
Cost_f1=value(Cost_f1);
display(['1、园区综合能源系统1的购能成本为 : ', num2str(Cost_f1)]);
Cost_f2=value(Cost_f2);
display(['2、园区综合能源系统1的运行成本为 : ', num2str(Cost_f2)]);
Cost_f3=value(Cost_f3);
display(['3、园区综合能源系统1的需求成本为 : ', num2str(Cost_f3)]);
Cost_f4=value(Cost_f4);
display(['4、园区综合能源系统1的环境成本为 : ', num2str(Cost_f4)]);
最新发布