请你检查以下代码,指出需要修改的地方:%% 基于能源集线器概念,结合需求侧柔性负荷的可平移、可转移、可削减特性,构建了含风光储、燃气轮机、柔性负荷等
% 综合考虑了系统运行成本和碳交易成本,建立了以总成本最低为优化目标的 IES 低碳经济调度模型,采用cplex求解器对算例进行求解。
% 场景2:只考虑包括可平移、可转移、可削减柔性电负荷,不考虑柔性热负荷参与系统优化调度的情况
clc
clear
%% 读取数据
% 电负荷、热负荷、光伏、风机、购电价、售电价
P_load_e = [160 150 140 140 130 135 150 180 215 250 275 320 335 290 260 275 270 280 320 360 345 310 220 160]; % 电负荷
P_load_h = [135 140 150 135 140 120 115 100 115 115 160 180 190 170 140 130 145 200 220 230 160 150 140 130]; % 热负荷
P_pv_max = [0 0 0 0 0 10 15 25 45 75 90 100 80 100 50 40 30 15 10 0 0 0 0 0 ]; % 光伏预测数据
P_wt_max = [60 65 70 75 80 85 90 100 125 150 130 110 100 120 125 130 140 160 180 200 175 160 155 150]; % 风机预测数据
C_buy_price = [0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.53 0.53 0.53 0.82 0.82 0.82 0.82 0.8 0.53 0.53 0.53 0.82 0.82 0.82 0.53 0.53 0.53]; % 购电价
C_sell_price = [0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.42 0.42 0.42 0.65 0.65 0.65 0.65 0.65 0.42 0.42 0.42 0.65 0.65 0.65 0.42 0.42 0.42]; % 售电价
%% 需求响应数据
P_cut = [10 10 10 10 10 10 15 15 25 50 50 50 50 50 50 50 50 50 50 50 40 40 15 10]; % 可削减电负荷
Temp_P_cut = binvar(1,24,'full'); % 电负荷削减标志
PP_cut = sdpvar(1,24,'full'); % 电负荷消减量
n1 = zeros(1,1); % 消减连续
H_cut = [25 25 25 25 25 25 25 25 30 40 40 40 40 40 40 40 40 40 50 50 30 30 20 15]; % 可削减热负荷
Temp_H_cut = binvar(1,24,'full'); % 热负荷削减标志
HH_cut = sdpvar(1,24,'full'); % 热负荷消减量
n2 = zeros(1,1); % 消减连续
P_tran = [0 0 0 0 0 0 0 0 0 0 0 0 25 25 25 25 0 0 0 0 0 0 0 0 ]; % 可转移电负荷
Temp_P_tran = binvar(1,24,'full'); % 可转移电负荷 转移标志
PP_tran = sdpvar(1,24,'full'); % 电负荷转移量
P_shift1 = [0 0 0 0 0 0 0 0 0 0 0 25 25 0 0 0 0 0 0 0 0 0 0 0 ]; % 可平移电负荷1
Temp_P_shift1 = binvar(1,24,'full'); % 可平移电负荷1 平移标志
PP_shift1 = sdpvar(1,24,'full'); % 可平移电负荷1量
P_shift2 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 25 25 25 0 0 ]; % 可平移电负荷2
Temp_P_shift2 = binvar(1,24,'full'); % 可平移电负荷2 平移标志
PP_shift2 = sdpvar(1,24,'full'); % 可平移电负荷2量
H_shift=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 45 45 45 0 0 0 0 ]; % 可平移热负荷
Temp_H_shift = binvar(1,24,'full'); % 可平移热负荷 平移标志
HH_shift=sdpvar(1,24,'full'); % 可平移热负荷量
for i=1:24
P_fix(i) = P_load_e(i)-P_shift1(i)-P_shift2(i)-P_tran(i)-P_cut(i); % 基础电负荷
end
for i=1:24
H_fix(i) = P_load_h(i)-H_shift(i)-H_cut(i); % 基础热负荷
end
%% 定义机组变量
P_pv_e = sdpvar(1,24,'full'); % 光伏电输出功率
P_wt_e = sdpvar(1,24,'full'); % 风机电输出功率
P_mt_e = sdpvar(1,24,'full'); % 燃气轮机电输出功率
P_GB_e = sdpvar(1,24,'full'); % 燃气锅炉输出热功率
P_buy_e = sdpvar(1,24,'full'); % 从电网购电电量
P_sell_e = sdpvar(1,24,'full'); % 向电网售电电量
P_net = sdpvar(1,24,'full'); % 与电网交换功率
Temp_net = binvar(1,24,'full'); % 购|售电标志
P_ESS_cha = sdpvar(1,24,'full'); % 充电功率
P_ESS_dis = sdpvar(1,24,'full'); % 放电功率
B_UP_cha = binvar(1,24,'full'); % 充电标志
B_UP_dis = binvar(1,24,'full'); % 放电标志
B_SOC_ess = sdpvar(1,24,'full');% 电储能余量
H_HSS_cha = sdpvar(1,24,'full'); % 储热充热
H_HSS_dis = sdpvar(1,24,'full'); % 储热放热
B_UH_cha = binvar(1,24,'full'); % 充热标志
B_UH_dis = binvar(1,24,'full'); % 放热标志
H_SOC_hss = sdpvar(1,24,'full'); % 热储能余量
%% 储能参数
% 电储能容量//自损/充电/放电
E_SOC_max = 0.95*100;
E_SOC_min = 0.4*100;
e_loss = 0.001;
e_cha = 0.9;
e_dis = 0.9;
% 热储能容量//自损/充热/放热
H_SOC_max = 0.95*100;
H_SOC_min = 0.4*100;
h_loss = 0.001;
h_cha = 0.9;
h_dis = 0.9;
%% 约束条件
Cons =[];
% 电储能容量约束、SOC约束、充电约束、放电约束、充放电状态约束、爬坡约束
B_SOC_ess(1,1) = E_SOC_min; % 电储能初始
for t=2:25 % 在一个周期内的充放电功率
Cons = [Cons,(B_SOC_ess(mod(t-1,24)+1)==(B_SOC_ess(mod(t-2,24)+1)*(1-e_loss)+(e_cha*P_ESS_cha(mod(t-2,24)+1)-(1/e_dis)*P_ESS_dis(mod(t-2,24)+1))))];
end
% 周期净交换功率为零, 初始功率相等即可
for i=1:24
Cons = [Cons,E_SOC_min<=B_SOC_ess(1,i)<=E_SOC_max]; % 容量约束限制
end
for i=1:24
Cons = [Cons,30*B_UP_cha(1,i)<=P_ESS_cha(1,i)<=40*B_UP_cha(1,i)]; % 电储能充电约束
Cons = [Cons,30*B_UP_dis(1,i)<=P_ESS_dis(1,i)<=40*B_UP_dis(1,i)]; % 电储能放电约束
end
% 蓄电池充放电约束
for i=1:24
Cons = [Cons,B_UP_cha(1,i)+B_UP_dis(1,i)<=1]; % 不同时充放电
end
Cons = [Cons,sum(B_UP_cha(1,1:24))+sum(B_UP_dis(1,1:24))==16]; % 使用寿命小于24
% 热储能容量约束、SOC约束、充热约束、放热约束、充放热状态约束
H_SOC_hss(1,1) = H_SOC_min; % 热储能初始
for t=2:25 % 在一个周期内的充放热功率
Cons = [Cons,(H_SOC_hss(mod(t-1,24)+1)==(H_SOC_hss(mod(t-2,24)+1)*(1-h_loss)+(h_cha*H_HSS_cha(mod(t-2,24)+1)-(1/h_dis)*H_HSS_dis(mod(t-2,24)+1))))];
end
% 全周期净交换功率为零, 初始功率相等即可
for i=1:24
Cons = [Cons,H_SOC_min<=H_SOC_hss(1,i)<=H_SOC_max]; % 容量约束限制
end
for i=1:24
Cons = [Cons,5*B_UH_cha(1,i)<=H_HSS_cha(1,i)<=30*B_UH_cha(1,i)]; % 热储能充电约束
Cons = [Cons,5*B_UH_dis(1,i)<=H_HSS_dis(1,i)<=30*B_UH_dis(1,i)]; % 热储能放电约束
end
% 蓄热池充放电约束
for i=1:24
Cons = [Cons,B_UH_cha(1,i)+B_UH_dis(1,i)<=1]; % 不同时充放热
end
Cons=[Cons,sum(B_UH_cha(1,1:24))+sum(B_UH_dis(1,1:24))==16]; % 使用寿命小于24
% 运行机组约束
for i=1:24
Cons = [Cons,0<=P_pv_e(i)<=P_pv_max(i)]; % 光伏上、下限约束
Cons = [Cons,0<=P_wt_e(i)<=P_wt_max(i)]; % 风机上、下限约束
Cons = [Cons,0<=P_mt_e(i)<=65]; % 燃气轮机上、下限约束
Cons = [Cons,0<=P_GB_e(i)<=160]; % 燃气锅炉上、下限约束
Cons = [Cons, -160<=P_net(i)<=160,0<=P_buy_e(i)<=160,-160<=P_sell_e(i)<=0]; % 主网功率交换约束
Cons = [Cons, implies(Temp_net(i),[P_net(i)>=0,P_buy_e(i)==P_net(i),P_sell_e(i)==0])]; % 购电情况约束
Cons = [Cons, implies(1-Temp_net(i),[P_net(i)<=0,P_sell_e(i)==P_net(i),P_buy_e(i)==0])]; % 售电情况约束
end
% 可平移电负荷1量
Cons = [Cons,sum(Temp_P_shift1(1,1:24)) == 2,sum(Temp_P_shift1(1,5:21)) == 2]; % 可平移电负荷1 平移标志
for i=5:20 % 时段区间为5~21-2+1
Cons = [Cons,sum(Temp_P_shift1(1,i:i+1)) >= 2*(Temp_P_shift1(1,i)-Temp_P_shift1(1,i-1))]; % 连续2个时段
end
for i=1:24
Cons = [Cons,PP_shift1(1,i) == 25*Temp_P_shift1(1,i)]; % 可平移电负荷1量
end
% 可平移电负荷2量
Cons = [Cons,sum(Temp_P_shift2(1,1:24)) == 3,sum(Temp_P_shift2(1,7:23)) == 3]; % 可平移电负荷2 平移标志
for i=7:21 % 时段区间为7~23-3+1
Cons = [Cons,sum(Temp_P_shift2(1,i:i+2)) >= 3*(Temp_P_shift2(1,i)-Temp_P_shift2(1,i-1)-Temp_P_shift2(1,i-2))]; % 连续3个时段
end
for i=1:24
Cons = [Cons,PP_shift2(1,i)== 25*Temp_P_shift2(1,i)]; %可平移电负荷2量
end
% 可转移电负荷(大于5自然会大于2)
for i=1:24
Cons = [Cons,Temp_P_tran(i)*8<=PP_tran(i)<=Temp_P_tran(i)*26.7 ]; % 可转移电负荷
end
Cons = [Cons,sum(Temp_P_tran(1,1:24)) == 5,sum(Temp_P_tran(1,4:22)) ==5]; % 可转移电负荷
Cons = [Cons,sum(Temp_P_tran(1,1:24)) ==5]; % 可转移电负荷
for i=4:18 % 时段区间为4~22-5+1
Cons = [Cons,sum(Temp_P_tran(1,i:i+4)) >= 5*(Temp_P_tran(1,i)-Temp_P_tran(1,i-1))];
end
% 可削减电负荷
Cons = [Cons,sum(Temp_P_cut)==8,sum(Temp_P_cut(1,5:22))==8];
Cons = [Cons,2<=n1<=5];
for i=5:22-n1+1 % 时段区间为5~22-n1+1
Cons = [Cons,sum(Temp_P_cut(1,i:i+n1-1)) >= n1*(Temp_P_cut(1,i)-Temp_P_cut(1,i-1))];
end
for i=1:24
Cons = [Cons,0<=PP_cut(1,i)<=Temp_P_cut(1,i)*0.9*P_cut(i)]; % 可消减电负荷
end
% 电、热平衡约束
for i=1:24
P_op_load_e(i) = P_fix(i)+P_cut(i)+PP_shift1(i)+PP_shift2(i)+PP_tran(i)-PP_cut(i);
end
for i=1:24
P_op_load_h(i) = H_fix(i)+H_cut(i)+HH_shift(i)-HH_cut(i);
end
for i=1:24
Cons = [Cons,P_mt_e(i)+P_pv_e(i)+P_wt_e(i)+P_net(i)-P_ESS_cha(1,i)+P_ESS_dis(1,i)==P_fix(i)+P_cut(i)+PP_shift1(i)+PP_shift2(i)+PP_tran(i)-PP_cut(i)]; % 电平衡约束
Cons = [Cons,P_GB_e(i)+0.83*P_mt_e(i)/0.45-H_HSS_cha(1,i)+H_HSS_dis(1,i)==P_load_h(i)]; % 热平衡约束
end
%% 目标函数
% 从大电网的购电成本
f1 = 0;
for i=1:24
f1 = f1+P_buy_e(i)*C_buy_price(i);
end
% 向大电网的售电成本
f2 = 0;
for i=1:24
f2 = f2+P_sell_e(i)*C_sell_price(i);
end
% 运维成本
f3 = 0;
for i=1:24
f3 = f3+0.72*P_pv_e(i)+0.52*P_wt_e(i); % 风机、光伏运维成本
end
% 燃料成本
f4 = 0;
for i=1:24
f4 = f4+2.5*P_GB_e(i)/9.7+2.5*P_mt_e(i)/0.45/9.7; % 耗气成本
end
% 储能成本
f5 = 0;
for i=1:24
f5 = f5+0.5*(P_ESS_cha(i)+P_ESS_dis(i)+H_HSS_cha(i)+H_HSS_dis(i)); % 储能成本
end
% 补偿成本
f6 = 0;
for i=1:24
f6 = f6+0.2*(PP_shift1(i)+PP_shift2(i))+0.3*PP_tran(i)+0.4*PP_cut(i);
end
% 碳交易成本
Q_carbon=0; % 碳排放量-碳配额量(克)
for i=1:24
Q_carbon=Q_carbon+(((1303-798)*(P_buy_e(i)+abs(P_sell_e(i)))+(564.7-424)*(P_GB_e(i)/9.7+P_mt_e(i)/0.45/9.7)+...
(43-78)*P_wt_e(i)+(154.5-78)*P_pv_e(i)+91.3*(P_ESS_cha(i)+P_ESS_dis(i))));
end
E_v=sdpvar(1,5); % 每段区间内的长度,分为5段,每段长度是2000
lamda=0.15*10^(-3); % 碳交易基价
Cons=[Cons,
Q_carbon==sum(E_v), % 总长度等于Q_carbon
0<=E_v(1:4)<=120000, % 除了最后一段,每段区间长度小于等于120000g
0<=E_v(5),
];
f7 = 0;
for v=1:5
f7 = f7+(lamda+(v-1)*0.25*lamda)*E_v(v);
end
% 目标函数
F = f3+f4+f1+f2+f5+f6+f7;
%% CPLEX求解器
ops = sdpsettings('solver','cplex', 'verbose',0); % 参数指定程序用cplex求解器
optimize(Cons,F,ops)
Cost_F = value(F);
display(['通过Yalmip求得的最优规划值为 : ', num2str(Cost_F)]);
display(['>>>购电成本为 : ', num2str(value(f1))]);
display(['>>>卖电收益为 : ', num2str(value(f2))]);
display(['>>>运维成本为 : ', num2str(value(f3))]);
display(['>>>燃料成本为 : ', num2str(value(f4))]);
display(['>>>储能成本为 : ', num2str(value(f5))]);
display(['>>>补偿成本为 : ', num2str(value(f6))]);
display(['>>>环境成本为 : ', num2str(value(f7))]);
P_pv_e=value(P_pv_e);
P_wt_e=value(P_wt_e);
P_mt_e=value(P_mt_e);
P_GB_e=value(P_GB_e);
P_ESS_cha=value(P_ESS_cha);
P_ESS_dis=value(P_ESS_dis);
H_HSS_cha=value(H_HSS_cha);
H_HSS_dis=value(H_HSS_dis);
P_buy_e=value(P_buy_e);
P_sell_e=value(P_sell_e);
PP_shift1=value(PP_shift1);
PP_shift2=value(PP_shift2);
PP_tran=value(PP_tran);
PP_cut=value(PP_cut);
HH_shift=value(HH_shift);
HH_cut=value(HH_cut);
%% 画图
figure(1)
ee_s=value([P_fix; P_cut; P_shift1; P_shift2; P_tran]);
bar(ee_s','stack');
grid on
hold on
plot(P_load_e,'g-*')
legend('基础电负荷','可消减电负荷','可平移电负荷1','可平移电负荷2','可转移电负荷','优化前电负荷');
xlabel('时间/h');
ylabel('电负荷功率/kW');
title('优化前:用户侧柔性电负荷分布');
figure(2)
for i=1:24
PPP_cut(i)=P_cut(i)-PP_cut(i); % 所剩的可消减电负荷
end
eee=value([P_fix; PPP_cut; PP_shift1; PP_shift2; PP_tran]);
bar(eee','stack');
grid on
hold on
plot(value(P_op_load_e),'g-*')
legend('基础电负荷','可消减电负荷','可平移电负荷1','可平移电负荷2','可转移电负荷','优化后电负荷');
xlabel('时间/h');
ylabel('电负荷功率/kW');
title('优化后:用户侧柔性电负荷分布');
figure(3)
hh=value([H_fix; H_cut; H_shift]);
bar(hh','stack');
grid on
hold on
plot(P_load_h,'r-*')
legend('基础热负荷','可消减热负荷','可平移热负荷','优化前热负荷');
xlabel('时间/h');
ylabel('热负荷功率/kW');
title('优化前:用户侧柔性热负荷分布');
figure(4)
for i=1:24
HHH_cut(i)=H_cut(i)-0; % 所剩的可消减热负荷
end
hhh_in=value([H_fix; HHH_cut; H_shift]);
bar(hhh_in','stack');
grid on
hold on
plot(value(P_load_h),'r-*')
legend('基础热负荷','可消减热负荷','可平移热负荷','优化后热负荷');
xlabel('时间/h');
ylabel('热负荷功率/kW');
title('优化后:用户侧柔性热负荷分布');
figure(5)
for i=1:24
P_op_load_e(i)=P_fix(i)+P_cut(i)+PP_shift1(i)+PP_shift2(i)+PP_tran(i)-PP_cut(i);
end
x=1:24;
plot(x,P_load_e,'-r*',x,P_op_load_e,'-bo');
grid on
xlabel('时间/h');
ylabel('电负荷/kW');
title('需求响应前、后电负荷曲线');
legend('优化前电负荷','优化后电负荷');
figure(6)
x=1:24;
plot(x,P_load_h,'-r*',x,P_load_h,'-bo');
grid on
xlabel('时间/h');
ylabel('热负荷/kW');
title('需求响应前后热负荷曲线');
legend('优化前热负荷','优化后热负荷');
figure(7)
stairs(x,C_buy_price,'-r')
hold on
stairs(x,C_sell_price,'-b')
hold on
grid on
xlabel('时间/h');
ylabel('价格/元');
title('市场分时电价曲线');
legend('购电价','售电价');
figure(8)
plot(x,P_load_e,'-o')
hold on
plot(x,P_load_h,'-s')
hold on
plot(x,P_pv_max,'-^')
hold on
plot(x,P_wt_max,'-p')
grid on
xlabel('时间/h');
ylabel('功率/kW');
title('源-荷动态预测曲线');
legend('电负荷','热负荷','光伏机组','风电机组');
figure(9)
eee_in=[P_buy_e; P_ESS_dis; P_pv_e; P_wt_e; P_mt_e];
eee_out=[P_sell_e; -P_ESS_cha];
bar(eee_in','stack');
hold on
bar(eee_out','stack');
grid on
plot(x,P_op_load_e,'-g*');
xlabel('时段');
ylabel('功率/kW');
legend('市场购电功率','蓄电池充电','光伏出力','风电出力','燃气轮机(电)','市场售电功率','蓄电池放电','优化后电负荷');
title('电负荷平衡分析');
figure(10)
hhh_in=[P_GB_e; H_HSS_dis; 0.83*P_mt_e/0.45];
hhh_out=[-H_HSS_cha];
bar(hhh_in','stack');
hold on
bar(hhh_out','stack');
grid on
plot(x,P_load_h,'-r*');
xlabel('时段');
ylabel('功率/kW');
legend('燃气锅炉','热储能放热','燃气轮机','热储能充热','优化后热负荷');
title('热负荷平衡分析');