检查以下代码,指出需要优化的地方:clc
clear
%% 设置条件
N=6;
T=24;
P_DJC_max=80;
P_DJC_b=10;
n_DJC=0.87;
N_DJC=0;
%% 决策变量初始化
P_CCHP_e=sdpvar(1,T); %CCHP的输出电功率
P_CCHP_h=sdpvar(1,T); %CCHP的输出热功率
P_CCHP_c=sdpvar(1,T); %CCHP的输出冷功率
P_g_CCHP=sdpvar(1,T); %CCHP消耗天然气功率
P_e_DJC=sdpvar(N,T); %电解槽的用电功率
P_DJC_H2=sdpvar(N,T); %电解槽的产氢功率
I=binvar(N,T); %停机
S=binvar(N,T); %冷待机
L=binvar(N,T); %变载运行
R=binvar(N,T); %过载运行
V=binvar(N,T); %低载运行
Y=binvar(N,T); %启动间隔
Z=binvar(N,T); %关机间隔
W=binvar(N,T); %制氢惩罚
X=binvar(N,T); %从停机到完全启动的标志
P_H_MR=sdpvar(1,T); %输入MR设备的氢能功率
P_MR_g=sdpvar(1,T); %MR设备输出的天然气功率
P_H_HFC=sdpvar(1,T); %输入HFC设备的氢能功率
P_HFC_e=sdpvar(1,T); %HFC设备输出的电功率
P_HFC_h=sdpvar(1,T); %HFC设备输出的热功率
P_PV=sdpvar(1,T); %光伏消纳功率
P_e_GB=sdpvar(1,T); %输入GB设备的天然气功率
P_GB_h=sdpvar(1,T); %GB设备输出的热功率
%储能部分(电ES1、热ES2、气ES3、氢ES4)
P_ES1_cha=sdpvar(1,T);P_ES2_cha=sdpvar(1,T);P_ES3_cha=sdpvar(1,T);P_ES4_cha=sdpvar(1,T); %充放功率
P_ES1_dis=sdpvar(1,T);P_ES2_dis=sdpvar(1,T);P_ES3_dis=sdpvar(1,T);P_ES4_dis=sdpvar(1,T);
S_1=sdpvar(1,T);S_2=sdpvar(1,T);S_3=sdpvar(1,T);S_4=sdpvar(1,T); %各储能的实时容量状态
%引入充放标志二进制变量
B_ES1_cha=binvar(1,T);B_ES2_cha=binvar(1,T);B_ES3_cha=binvar(1,T);B_ES4_cha=binvar(1,T); %充标志
B_ES1_dis=binvar(1,T);B_ES2_dis=binvar(1,T);B_ES3_dis=binvar(1,T);B_ES4_dis=binvar(1,T); %放标志
P_e_buy=sdpvar(1,T); %购电功率
P_g_buy=sdpvar(1,T); %购气功率
%% 导入常数参数
%购电价格
jf=0.7; %峰段电价
jg=0.25; %谷段电价
jp=0.6487; %平段电价
c_e_buy=[jg.*ones(1,7) jf.*ones(1,5) jp.*ones(1,5) jf.*ones(1,4) jp.*ones(1,3)]; %电价情况
%电负荷
P_e_load=[1632.49275595938,1572.19275595938,1548.99275595938,1609.29275595938,1660.29275595938,1957.29275595938,2244.89275595938,1623.55187999109,1878.75187999109,2022.55187999109,2124.65187999109,2138.55187999109,2418.68481725823,2418.68481725823,2460.38481725823,2441.88481725823,2529.98481725823,2249.95187999109,2180.35187999109,2004.05187999109,1748.85187999109,2019.68481725823,1374.68481725823,1203.08481725823];
%热负荷
P_h_load=[1015.4 917.9 919.6 885.9 1192.2 1661.6 1603.7 1497.3 1570.9 1681.9 1378.7 1412.1 1512.3 1435.7 1476.2 1560.3 1666.7 1746 1780.5 2031.2 2152.3 1905.2 1705.6 1153.8];
%冷负荷
P_c_load=[229.916897506925;224.376731301939;216.066481994460;221.606648199446;224.376731301939;252.077562326870;268.698060941828;288.088642659280;299.168975069252;288.088642659280;293.628808864266;282.548476454294;279.778393351801;271.468144044321;271.468144044321;268.698060941828;277.008310249307;293.628808864266;307.479224376731;304.709141274238;293.628808864266;285.318559556787;277.008310249307;265.927977839335]';
%光伏出力最大值
P_PV_max=[1870+100 1760+100 1362.5+250 1500.1 1220.2+150 1360.7 1394.7 1067.08+200 1370.8 1192.7 1321.4 1340 1280 1428.2 1560.7 1600 1572 1739 1660.1 1677 1800 1756.8 1790.1 1856.6 ];
%购气价格
c_g_buy=0.35*ones(1,T);
%% 导入约束条件
C=[];
C=C+[P_DJC_H2==0.9*(P_e_DJC-P_DJC_b*S)]; %电氢转换关系
C=C+[P_DJC_b*S+0.3*P_DJC_max*L+P_DJC_max*R+0.2*P_DJC_max*V<=P_e_DJC<=0.3*P_DJC_max*V+1.5*P_DJC_max*R+P_DJC_b*S+P_DJC_max*L]; %PEM功率限制
C=C+[-I(:,1:T-2)+I(:,2:T-1)-I(:,3:T)<=0]; %停机时间约束
C=C+[L+S+I+R+V==1]; %状态互斥约束
C=C+[V(:,2:T)+R(:,2:T)+L(:,2:T)+S(:,2:T)+I(:,1 :T-1)-1<=Y(:,2:T)]; %启停间隔约束
C=C+[R(:,1:T-1)+L(:,1:T-1)+S(:,1:T-1)+I(:,1:T-1)-1<=Z(:,2:T)]; %启停间隔约束
C=C+[W(:,2:T)<=S(:,1:T-1)];
C=C+[W<=L+R+V];
C=C+[W(:,2:T)>=S(:,1:T-1)+L(:,2:T)+R(:,2:T)+V(:,2:T)-1];
%过载、低载时间限制
C=[C,sum(R)<=3];
for i=1:T-2
C=C+[R(:,i)+R(:,i+1)+R(:,i+2)<=2];
end
for i=1:T-2
C=C+[V(:,i)+V(:,i+1)+V(:,i+2)<=2];
end
C=[C,
%CHP的电-气能量转换约束
0<=P_g_CCHP<=1200, %CHP消耗的气功率上下限约束
P_CCHP_e==0.35*P_g_CCHP,
P_CCHP_h==0.45*P_g_CCHP,
-500<=P_g_CCHP(2:T)-P_g_CCHP(1:23)<=500, %CHP的爬坡约束(1-24时段)
];
C=[C,
P_MR_g==0.6*P_H_MR, %MR(甲烷反应器)的气-氢能量转换约束
0<=P_H_MR<=300, %MR消耗的氢功率的上下限约束
-0.2*300<=P_H_MR(2:T)-P_H_MR(1:23)<=0.5*300, %MR的爬坡约束(1-24时段)
];
C=[C,
P_HFC_e+P_HFC_h==0.95*P_H_HFC, %HFC(氢燃料电池)的电-氢能量转换约束
0<=P_H_HFC<=500, %HFC消耗的氢功率上下限约束
0.5*P_HFC_e<=P_HFC_h, %HFC的热电比上下限
P_HFC_h<=2.1*P_HFC_e,%HFC的热电比上下限
-0.2*500<=P_H_HFC(2:T)-P_H_HFC(1:23)<=0.2*500, %HFC的爬坡约束(1-24时段)
];
C=[C,
0<=P_PV<=P_PV_max, %风电出力约束
P_GB_h==0.90*P_g_GB, %GB的热-气能量转换约束
0<=P_g_GB<=1200, %GB的出力上下限约束
-500<=P_g_GB(2:T)-P_g_GB(1:23)<=500, %GB的爬坡约束(1-24时段)
];
C=[C,
0<=P_ES1_cha<=B_ES1_cha*0.5*1000, %储电设备的最大充电功率约束
0<=P_ES2_cha<=B_ES2_cha*0.5*1000, %储热设备的最大充热功率约束
0<=P_ES3_cha<=B_ES3_cha*0.5*1000, %储气设备的最大充气功率约束
0<=P_ES4_cha<=B_ES4_cha*0.5*400, %储氢设备的最大充氢功率约束
0<=P_ES1_dis<=B_ES1_dis*0.5*1000, %储电设备的最大放电功率约束
0<=P_ES2_dis<=B_ES2_dis*0.5*1000, %储热设备的最大放热功率约束
0<=P_ES3_dis<=B_ES3_dis*0.5*1000, %储气设备的最大放气功率约束
0<=P_ES4_dis<=B_ES4_dis*0.5*400, %储氢设备的最大放氢功率约束
S_1(1)==0.1*1000, %储电设备的初始容量
S_2(1)==0.1*1000, %储热设备的初始容量
S_3(1)==0.1*1000, %储气设备的初始容量
S_4(1)==80, %储氢设备的初始容量
%储能爬坡约束
-0.2*1000<=P_ES1_cha(2:T)-P_ES1_cha(1:23)<=0.2*1000
-0.2*1000<=P_ES1_dis(2:T)-P_ES1_dis(1:23)<=0.2*1000
-0.2*1000<=P_ES2_cha(2:T)-P_ES2_cha(1:23)<=0.2*1000
-0.2*1000<=P_ES2_dis(2:T)-P_ES2_dis(1:23)<=0.2*1000
-0.2*1000<=P_ES3_cha(2:T)-P_ES3_cha(1:23)<=0.2*1000
-0.2*1000<=P_ES3_dis(2:T)-P_ES3_dis(1:23)<=0.2*1000
-80<=P_ES4_cha(2:T)-P_ES4_cha(1:23)<=80
-80<=P_ES4_dis(2:T)-P_ES4_dis(1:23)<=80
%始末状态守恒约束
S_1(T)==S_1(1),
S_2(T)==S_2(1),
S_3(T)==S_3(1),
S_4(T)==S_4(1),
%充放状态唯一
B_ES1_cha+B_ES1_dis<=1,
B_ES2_cha+B_ES2_dis<=1,
B_ES3_cha+B_ES3_dis<=1,
B_ES4_cha+B_ES4_dis<=1,
%储能容量上下限约束
0.1*1000<=S_1<=0.9*1000,
0.1*1000<=S_2<=0.9*1000,
0.1*1000<=S_3<=0.9*1000,
80<=S_4<=400,
%储能容量变化约束
S_1(2:T)==S_1(1:23)+0.95*P_ES1_cha(2:T)-P_ES1_dis(2:T)/0.95,
S_2(2:T)==S_2(1:23)+0.95*P_ES2_cha(2:T)-P_ES2_dis(2:T)/0.95,
S_3(2:T)==S_3(1:23)+0.95*P_ES3_cha(2:T)-P_ES3_dis(2:T)/0.95,
S_4(2:T)==S_4(1:23)+0.95*P_ES4_cha(2:T)-P_ES4_dis(2:T)/0.95,
];
C=[C, sum(0.95*P_ES1_cha)-sum(P_ES1_dis/0.95)==0];
C=[C,sum(0.95*P_ES2_cha)-sum(P_ES2_dis/0.95)==0];
P_PEM_H2_N=sum(P_DJC_H2);
P_PEM_E_N=sum(P_e_DJC);
%% 功率平衡条件
C=[C,
P_e_buy==P_e_load+P_PEM_E_N+P_ES1_cha-P_ES1_dis-P_PV-P_CCHP_e-P_HFC_e, %电功率平衡约束
P_HFC_h+P_CCHP_h+P_GB_h==P_h_load+P_ES2_cha-P_ES2_dis, %热功率平衡约束
P_g_buy==P_c_load+P_ES3_cha-P_ES3_dis+P_g_CCHP+P_g_GB-P_MR_g, %气功率平衡约束
P_PEM_H2_N==P_H_MR+P_H_HFC+P_ES4_cha-P_ES4_dis, %氢功率平衡约束
0<=P_e_buy<=5000, %购电功率约束
0<=P_g_buy<=5000, %购气功率约束
];
%% 阶梯碳交易
%碳排放权配额模型
E_e_buy=0.728*sum(P_e_buy); %购电配额
E_CHP=0.57*sum(P_CCHP_h+6/3.6*P_CCHP_e); %CHP配额
E_GB=0.57*sum(P_GB_h); %GB配额
E_IES=E_e_buy+E_CHP+E_GB; %IES总碳排放配额
%实际碳排放模型
E_e_buy_a=1.08*sum(P_e_buy);
E_CHP_a=0.6101*sum(P_CCHP_h+6/3.6*P_CCHP_e);
E_GB_a=0.6101*sum(P_GB_h);
E_MR_a=0.5*sum(P_MR_g); %实际MR减少的碳排放
E_IES_a=E_e_buy_a+E_CHP_a+E_GB_a;
E=E_IES_a; %实际IES总碳排放
%阶梯碳交易成本(分段线性化)
E_v=sdpvar(1,5) %每段区间内的长度,分为5段,每段长度是2000
lamda=0.267; %碳交易基价
E_C= E-E_IES;
C=[C,
E_C==sum(E_v), %总长度等于E
0<=E_v(1:4)<=1000, %除了最后一段,每段区间长度小于等于2000
0<=E_v(5),
];
%碳交易成本
C_CO2=0;
for v=1:5
C_CO2=C_CO2+(lamda+(v-1)*0.25*lamda)*E_v(v);
end
%% 总目标函数
CY=30*sum(sum(Y));
obj=c_e_buy*P_e_buy'+c_g_buy*P_g_buy'+0.2*sum(P_PV_max-P_PV)+C_CO2+CY; %购能成本+弃风成本+碳交易成本+启停成本
C_E_BUY=c_e_buy*P_e_buy';
C_G_BUY=c_g_buy*P_g_buy';
C_wind=0.2*sum(P_PV_max-P_PV);
%% 求解器相关配置
ops=sdpsettings('solver','cplex','verbose',2,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
%% 进行求解计算
result=optimize(C,obj,ops);
if result.problem==0
% problem =0 代表求解成功
% do nothing!
else
error('求解出错');
end