%并网模式下，碳排放最优调度模型
clc;clear;close all;% 程序初始化
%% 读取数据
shuju=xlsread('shuju数据.xlsx'); %把一天划分为24小时
load_e=shuju(2,:); %电负荷
load_g=shuju(3,:); %气负荷
P_PV=shuju(4,:); %光电出力预期值
P_WT=shuju(5,:); %风电出力预期值
T_out=shuju(6,:); %室外温度
%% 各变量及常量定义
P_G3=sdpvar(1,24,'full');%微燃气轮机电功率出力
e_G3=0.26;%燃气轮机电效率
h_G3=0.68;%燃气轮机热效率
P_EH=sdpvar(1,24,'full');%余热锅炉输出热功率
EH=0.6;%余热回收效率
P_GH=sdpvar(1,24,'full');%燃气锅炉输出热功率
GH=0.85;%燃气锅炉效率
P_AC=sdpvar(1,24,'full');%吸收式制冷机输出冷功率
COP_AC=0.8;%吸收式制冷机制冷系数
P_EC=sdpvar(1,24,'full');%电制冷机输出冷功率
COP_EC=3;%电制冷机冷系数
P_EG=sdpvar(1,24,'full');%P2G设备输出气功率
EG=0.6;%P2G设备综合转换效率
Pbuy=sdpvar(1,24,'full');%从电网购电电量
Psell=sdpvar(1,24,'full');%向电网售电电量
Pnet=sdpvar(1,24,'full');%交换功率
Temp_net=binvar(1,24,'full'); % 购|售电标志
Gbuy=sdpvar(1,24,'full');%从气网购气量
R=0.93;cc=0.54; %房间热阻和热容，这里分析改变等效热阻Ｒ时对冷负荷的影响
load_h=sdpvar(1,24);
T_hui=sdpvar(1,24);
%% 约束条件
C =[];
%% 热冷负荷
%对于热负荷，利用PMV值得到供热时室内温度
T_out_hot=sdpvar(1,26); %供热时室外温度
for i=3:26
C=[C,T_out_hot(1)==0;T_out_hot(2)==0];%-1、0时刻的供热室外温度
C=[C,T_out_hot(i)==T_out(i-2)];
end
T_in_hot=sdpvar(1,26); %供热时室内温度
T_gong=sdpvar(1,26); %供水温度
PMV=1;%PMV值，公式4,这里分析改变PMV指标时对热负荷的影响
for i=3:26
C=[C,T_in_hot(1)==25;T_in_hot(2)==25];%-1、0时刻的供热室内温度
C=[C,-PMV<=(0.303*exp(-0.036*70)+0.028)*(70-0-3.05*(10^(-3))*(5733-6.99*70-2000)-0.42*(70-58.15)-1.7*(10^(-5))*70*(5867-2000)-1.4*(10^(-3))*70*(34-T_in_hot(i)-3.96*(10^(-8))*1.15*((32+273)^4-(29.7+273)^4))-1.15*4.7*(32-T_in_hot(i)))<=PMV];
end
%公式1，2
for i=3:26
C=[C,T_gong(1)==75;T_gong(2)==75];%热网供水初始温度
C=[C,T_in_hot(i)==0.6991*T_in_hot(i-1)+0.1011*T_gong(i-1)+0.1998*T_out_hot(i-1)]; %公式1
C=[C,T_hui(i-2)<=T_gong(i)<=120];%热网供水温度
C=[C,T_hui(i-2)==0.5721*T_in_hot(i-1)+0.0607*T_in_hot(i-2)+0.2112*T_gong(i)-0.0243*T_gong(i-1)-0.0104*T_gong(i-2)+0.3317*T_out_hot(i)-0.3169*T_out_hot(i-1)+0.1741*T_out_hot(i-2)];%公式2
C=[C,5<=load_h(i-2)];%热负荷下限
C=[C,load_h(i-2)==0.63*(T_gong(i)-T_hui(i-2))]; %公式21
end
% %对于冷负荷，已知条件为
%根据公式3得到冷负荷，
load_c=sdpvar(1,24);%冷负荷
T_in_cold=sdpvar(1,25);%供冷时室内温度，初始为-15
for i=2:25
C=[C,T_in_cold(1)==-15];
C=[C,-20<=T_in_cold(i)<=-15];
C=[C,load_c(i-1)*(1-exp(-1/(R*cc)))==-T_in_cold(i)*R+T_out(i-1)*(1-exp(-1/(R*cc)))+T_in_cold(i-1)*exp(-1/(R*cc))];
end
%% 机组约束
for i=1:24
C=[C,9<=P_G3(i)<=18];%燃气轮机上下限约束
C=[C,0<= P_EH(i)<=43];%余热锅炉上下限约束
C=[C,0<=P_GH(i)<=36];%锅炉上下限约束
C=[C,0<=P_AC(i)<=20];%吸收式制冷机出力上下限约束
C=[C,0<=P_EC(i)<=30];%电制冷机出力上下限约束
C=[C,0<=P_EG(i)<=30];%P2G出力上下限约束
C = [C, -15<=Pnet(i)<=15,0<=Pbuy(i)<=15, -15<=Psell(i)<=0]; %主网功率交换约束
C = [C, implies(Temp_net(i),[Pnet(i)>=0,Pbuy(i)==Pnet(i),Psell(i)==0])]; %购电情况约束
C = [C, implies(1-Temp_net(i),[Pnet(i)<=0,Psell(i)==Pnet(i),Pbuy(i)==0])]; %售电情况约束
C= [C,P_EH(i)==P_WT(i)/e_G3*h_G3*GH]; %余热回收约束
end
for i=1:24
C = [C,P_WT(i)+P_PV(i)+Pnet(i)+P_G3(i)-P_EC(i)/COP_EC-P_EG(i)/EG==load_e(i)]; %电平衡
C = [C,P_EH(i)+P_GH(i)-P_AC(i)/COP_AC==load_h(i)];%热平衡约束
C=[C,P_EC(i)+P_AC(i)==load_c(i)];%冷平衡约束
C=[C,Gbuy(i)+P_EG(i)==P_WT(i)/e_G3+P_GH(i)/GH+load_g(i)];%气平衡约束
end
%% 目标函数
%------------------碳排最小--------------------%
%购天然气成本
C_Ng=0;
C_gas=0.25; %气价参数，设置不同的气价参数进行灵敏度分析
for i=1:24
C_Ng=C_Ng+C_gas*Gbuy(i);
end
%卖电收益
C_sell=0;
C_sell_ele=0.2;%大电网卖电价参数，设置不同的卖电价参数进行灵敏度分析
for i=1:24
C_sell=C_sell+Psell(i)*C_sell_ele;
end
%买电成本
C_buy=0;
C_buy_ele=0.7;%大电网买电价参数，设置不同的买电价参数进行灵敏度分析
for i=1:24
C_buy=C_buy+C_buy_ele*Pbuy(i);
end
%碳排放量
C_carbon_emission=0;
for i=1:24
C_carbon_emission=C_carbon_emission+0.4483*Gbuy(i)+0.805*Pbuy(i); %天然气碳排，热力碳排
end
%目标函数
F=C_carbon_emission; %此时是以碳排放最优的调度模型
ops = sdpsettings('solver','cplex', 'verbose', 2);%通过将verbose设置为0，解决者将以最小的显示运行。通过增加值，显示级别被控制（通常1提供了适度的显示级别，而2提供了大量的信息）。
optimize(C,F,ops)
value(F)%费用
C_Ng=value(C_Ng)
C_sell=value(C_sell)
C_buy=value(C_buy)
C_carbon_emission=value(C_carbon_emission)
load_h=value(load_h)
load_c=value(load_c)
%% 画图
x=1:24;
figure
plot(x,P_PV,'-r')
hold on
plot(x,P_WT,'--b')
legend('光伏预测曲线','风机预测曲线');
T_h=value(T_in_hot(1,3:26));
figure
plot(x,T_h,'-r',x,T_out,'-b')
hold on
legend('室内温度','室外温度');
figure
plot(x,load_e,'g-*',x,value(load_h),'r-o',x,value(load_c),'b-s',x,load_g,'m-*');
xlabel('时间/h');
ylabel('负荷/MW');
title('负荷曲线');
legend('电负荷','热负荷','冷负荷','气负荷');
tt=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
%电平衡
x=1:24;
PP=value([Pbuy;P_G3;P_WT;P_PV;tt;tt]);
PP1=value([Psell;tt;tt;tt;-P_EC/COP_EC;-P_EG/EG]);
figure
bar(PP','stack');
legend('电网','燃气轮机出力','风电出力','光伏出力','电制冷机耗电');
hold on
bar(PP1','stack');
plot(x,load_e,'g-*');
xlabel('时段');ylabel('功率/MW');
title('电网络平衡');
%热平衡
x=1:24;
QQ=value([P_GH;P_EH;tt]);
gg1=value([tt;tt;-P_AC/COP_AC]);
figure
bar(QQ','stack');
legend('燃气锅炉出力','余热锅炉出力','吸收式制冷机出力');
hold on
bar(gg1','stack');
plot(x,value(load_h),'r-o');
xlabel('时段');ylabel('功率/MW');
title('热网络平衡');
%冷平衡
x=1:24;
LL=value([P_EC;P_AC]);
figure
bar(LL','stack');
legend('电制冷机','吸收式制冷机');
hold on
plot(x,value(load_c),'b-s');
xlabel('时段');ylabel('功率/MW');
title('冷网络平衡');
%气平衡
x=1:24;
gg=value([Gbuy;P_EG;tt;tt]);
gg1=value([tt;tt;-P_WT/e_G3;-P_GH/GH]);
figure
bar(gg','stack');
legend('气网','P2G','燃气轮机耗气','燃气锅炉耗气');
hold on
bar(gg1','stack');
plot(x,load_g,'m-*');
xlabel('时段');ylabel('功率/MW');
title('气网络平衡');
