请你检查以下代码,指出不符合逻辑的地方: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_h = 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_hes = 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_hes(1,1) = H_SOC_min; % 热储能初始
for t=2:25 % 在一个周期内的充放热功率
Cons = [Cons,(H_SOC_hes(mod(t-1,24)+1)==(H_SOC_hes(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_hes(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_h(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
% 可平移热负荷量
Cons = [Cons,sum(Temp_H_shift(1,1:24)) == 3,sum(Temp_H_shift(1,5:21)) == 3];% 可平移热负荷 平移标志
for i=5:19 % 时段区间为5~21-3+1
Cons = [Cons,sum(Temp_H_shift(1,i:i+2)) >= 3*(Temp_H_shift(1,i)-Temp_H_shift(1,i-1))]; % 连续3个时段
end
for i=1:24
Cons = [Cons,HH_shift(1,i)== 45*Temp_H_shift(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
% 可削减热负荷
Cons=[Cons,sum(Temp_H_cut(1,1:24))==8,sum(Temp_H_cut(1,11:19))==8];
Cons=[Cons,2<=n2<=5];
for i=11:19-n2+1 % 时段区间为11~19-n2+1
Cons = [Cons,sum(Temp_H_cut(1,i:i+n1-1)) >= n1*(Temp_H_cut(1,i)-Temp_H_cut(1,i-1))];
end
for i=1:24
Cons = [Cons,Temp_H_cut(1,i)*0<=HH_cut(1,i)<=Temp_H_cut(1,i)*0.9*H_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_h(i)+0.83*P_mt_e(i)/0.45-H_HSS_cha(1,i)+H_HSS_dis(1,i)==H_fix(i)+H_cut(i)+HH_shift(i)-HH_cut(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_h(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.1*HH_shift(i)+0.3*PP_tran(i)+0.4*PP_cut(i)+0.2*HH_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_h(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 = f1+f2+f3+f4+f5+f6+f7;
%% CPLEX求解器
ops = sdpsettings('solver','cplex', 'verbose', 0); % 参数指定程序用cplex求解器
r = 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))]);
P_pv_e=value(P_pv_e);
P_wt_e=value(P_wt_e);
P_mt_e=value(P_mt_e);
P_GB_h=value(P_GB_h);
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=value([P_fix; P_cut; P_shift1; P_shift2; P_tran]);
bar(ee',1,'stack')
hold on
plot(P_fix+P_cut+P_shift1+P_shift2+P_tran,'g-*','linewidth',1)
hold on
grid on
xlabel('时间/h');
ylabel('电负荷功率/kW');
legend('基础电负荷','可消减电负荷','可平移电负荷1','可平移电负荷2','可转移电负荷','等效电负荷');
title('优化前用户侧柔性电负荷分布');
figure(2)
for i=1:24
PPP_cut(i) = P_cut(i)-PP_cut(i); % 所剩的可消减电负荷
end
ee1=value([P_fix; PPP_cut; PP_shift1; PP_shift2; PP_tran]);
bar(ee1','stack');
hold on
plot(P_fix+PPP_cut+PP_shift1+PP_shift2+PP_tran,'g-*','linewidth',1)
hold on
legend('基础电负荷','可消减电负荷','可平移电负荷1','可平移电负荷2','可转移电负荷','等效电负荷');
xlabel('时间/h');
ylabel('电负荷功率/kW');
title('优化后用户侧柔性电负荷分布');
figure(3)
hh=value([H_fix; H_cut; H_shift]);
bar(hh',1,'stack');
hold on
plot(H_fix+H_cut+H_shift,'c-*','linewidth',1)
hold on
grid on
legend('基础热负荷','可消减热负荷','可平移热负荷','等效热负荷');
xlabel('时间/h');
ylabel('热负荷功率/kW');
title('优化前用户侧柔性热负荷分布');
figure(4)
for i=1:24
HHH_cut(i) = H_cut(i)-HH_cut(i); % 所剩的可消减热负荷
end
hh1=value([H_fix;HHH_cut;HH_shift]);
bar(hh1','stack');
hold on
plot(H_fix+HHH_cut+HH_shift,'r-*','linewidth',1)
hold on
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
plot(P_load_e,'g-o');
hold on
grid on
plot(P_op_load_e,'g-*');
xlabel('时间/h');
ylabel('电负荷/kW');
title('需求响应前、后电负荷曲线');
legend('优化前电负荷','优化后电负荷');
figure(6)
for i=1:24
P_op_load_h(i)=H_fix(i)+H_cut(i)+HH_shift(i)-HH_cut(i);
end
plot(P_load_h,'r-o');
hold on
grid on
plot(P_op_load_h,'r-*');
xlabel('时间/h');
ylabel('热负荷/kW');
title('需求响应前、后热负荷曲线');
legend('优化前热负荷','优化后热负荷');
figure(7)
E_v=value(E_v);
bar(E_v,0.5)
xlabel('时间/4h');
ylabel('碳交易量/吨');
yyaxis right
ecc=(lamda+(v-1)*0.25*lamda)*E_v;
plot(ecc,'r-*')
ylabel('碳交易成本/元');
ylim([24 50]);
grid on
legend('阶梯式碳交易量','阶梯式碳交易成本');
title('碳排放轨迹追踪');
figure(8)
s=(((1303-798)*(P_buy_e+abs(P_sell_e))+(564.7-424)*(P_GB_h/9.7+P_mt_e/0.45/9.7)+...
(43-78)*P_wt_e+(154.5-78)*P_pv_e+91.3*(P_ESS_cha+P_ESS_dis)));
bar(s,0.5,'b')
ylabel('碳交易量/吨');
hold on
grid on
yyaxis right
plot(C_buy_price,'r-*');
xlabel('时间/h');
ylabel('市场电价/元');
title('碳交易总量');
legend('碳交易总量','市场电价');
figure(9)
ee_in = value([P_buy_e; P_ESS_dis; P_pv_e; P_wt_e; P_mt_e]);
ee_out = value([P_sell_e; -P_ESS_cha; -P_op_load_e]);
bar(ee_in','stack');
hold on
grid on
bar(ee_out','stack');
plot(P_op_load_e,'*-g');
legend('上网购电','蓄电池充电','光伏出力','风电出力','燃气轮机(电)','上网卖电','蓄电池放电','优化后电负荷');
title('电负荷平衡');
xlabel('时段');
ylabel('功率/kW');
figure(10)
hh_in = value([P_GB_h; H_HSS_dis; 0.83*P_mt_e/0.45]);
hh_out = value([-H_HSS_cha; -P_op_load_h]);
bar(hh_in','stack');
hold on
grid on
bar(hh_out','stack');
plot(P_op_load_h,'*-r');
legend('燃气锅炉','热储能充热','燃气轮机(热)','热储能放热','优化后热负荷');
title('热负荷平衡');
xlabel('时段');
ylabel('功率/kW');
最新发布