clear all
%% 输入原始数据 光照强度：w/m^2 风速：m/s 电热负荷：kw等
Solar=[0 0 0 0 0 9 115 282 454 620 744 801 793 709 593 457 297 127 15 0 0 0 0 0];%光照强度
Wind=[7.14 7.20 7.22 6.89 7.18 7.07 6.43 5.97 6.33 6.56 6.63 6.70 6.65 6.92 6.86 6.85 6.89 7.01 7.02 6.97 7.05 7.06 7.03 7.10];%风速
Load_E =[30420.49 32008.64 32673.34 33327.69 29716.69 29540.05 27917.86 37536.115 44387.93 42669.255 44470.04 32343.06 34827.06 41421.16 41072.135 43739.79 38737.75 33883.05 36011.1 34328.19 34977.48 32621.59 32920.36 32306.72];%电负荷
Load_H =[27394.52 35116.88 36361.76 37704.68 38214.26 40975.34 42859.76 50462.42 55585.58 58076.48 59895.92 56950.16 57761.84 58141.46 57643.28 59874.26 55563.92 55509.2 48967.88 44603.96 43942.76 36773.13 34553.72 30860.12];%热负荷
Price_E=[0.35 0.35 0.35 0.35 0.35 0.35 0.35 0.61 0.97 0.97 0.97 0.61 0.35 0.61 0.97 0.97 0.97 0.97 0.97 0.97 0.97 0.61 0.61 0.35];%分时电价，元/kwh
Price_Gas=2.7;%天然气价格，元/m^3
a_Gas=35600;%天然气热值，35800kJ/m^3
alpha=0.08;%碳捕集设备选择捕集产生的二氧化碳比例，燃烧产生的二氧化碳不一定要全部捕集，这里0.08=8%，至少捕集8%产生的二氧化碳
%% 输入设备模型
%光伏
S_PV=40000;%光伏面积
Q_PV=0.15*S_PV;%一平米折合约0.15kw
Pe_PV=0.001*0.157*S_PV*Solar;%光伏功率，效率0.157，0.001转换成kw/m2,光照强度已载入
PV_omfix=0.02*Q_PV;%固定运维成本,按天结算，如设备清灰这样的人工费
PV_omwork=0.039*sum(Pe_PV);%可变运维成本，与发电量有关
PV_om=PV_omfix+PV_omwork;%总运维成本
%风电WP
N_WP=4;%总台数
Q_sWP=2000;%单台额定容量2MW
for i=1:24 %风机，单台额定2MW，切入风速3，额定12，切出20
if 3<=Wind(i)&&Wind(i)<=12
Pe_swind(i)=222.2*(Wind(i)-3);
elseif 12<Wind(i)&&Wind(i)<=20
Pe_swind(i)=Q_sWP;
else
Pe_swind(i)=0;
end
end
Pe_WP=N_WP*Pe_swind;%总出力
WP_om=0.07*sum(Pe_WP);%风电运维
%热电联产CHP
Q_CHP=25000;%CHP容量，kw
Qg_CHP= sdpvar(1, 24);%CHP耗气速率，m3/s
Prg_CHP=Qg_CHP*a_Gas;%CHP燃烧天然气产生功率
Pe_CHP=0.31*Prg_CHP;%CHP输出电功率，产电效率0.31
Ph_CHP=0.5*(1-0.31)*Prg_CHP;%CHP产热功率;余热回收效率0.5
CHP_omfix=0.02*Q_CHP;
CHP_omwork=0.03*sum(Pe_CHP);
CHP_om=CHP_omfix+CHP_omwork;
%燃气锅炉GB
Q_GB=18000;%燃气锅炉容量，kw
Qg_GB= sdpvar(1, 24);%燃气锅炉耗气速率，m3/s
Ph_GB=0.9*Qg_GB*a_Gas;%燃气锅炉产热,效率0.9
GB_omfix=0.03*Q_GB;
GB_omwork=0.02*sum(Ph_GB);
GB_om=GB_omfix+GB_omwork;
%电锅炉
Q_EB=16000;%电锅炉容量，kw
Pe_EB=sdpvar(1,24);%电锅炉电功率
Ph_EB=0.95*Pe_EB;%电热锅炉产热，效率0.95
EB_omfix=0.03*Q_EB;
EB_omwork=0.04*sum(Ph_EB);
EB_om=EB_omfix+EB_omwork;
%电储能，EES
Q_EES=5000;%蓄电池容量kwh
Qe_EES= sdpvar(1, 25);%电池容量
Pe_EESc= sdpvar(1, 24);%充电功率
Pe_EESdisc= sdpvar(1, 24);%放电功率
EES_omfix=0.03*Q_EES;
EES_omwork=0.02*sum(Pe_EESc+Pe_EESdisc);
EES_om=EES_omfix+EES_omwork;
%碳捕集设备CCS
Q_CCS=6000;% 最大电功耗，除以900为最大捕集能力 kg/s（CCS工作电耗：0.25kwh/kg）
CCS_max=Q_CCS/900;%每秒最大捕集能力，kg
CCS_maxyear=CCS_max*24*3600*365*0.001;%年最大捕集量，t/a
Yita_CCS=0.9;%碳捕集效率
Q_CO2in=sdpvar(1, 24);%输入碳捕集设备的CO2量可控,kg/s
Q_CO2product=1.964*(Qg_CHP+Qg_GB);%CO2直接产生量，1m^3天然气燃烧产生1.964kgCO2，kg/s
Q_CO2capture=Yita_CCS*Q_CO2in;%CO2捕集量，kg/s
Pe_CCS=900*Q_CO2in;%CCS工作电耗，0.25kwh/kg,0.25*3600=900
Ph_CCS=3500*Q_CO2capture;%CCS工作热耗，3.5GJ/t，3500kJ/kg
CCS_om=0.14*CCS_maxyear;%每天的运行维护成本与年最大捕集量有关
Pe_buy=sdpvar(1,24);%电网购电
%% 约束条件
Cons=[
Pe_buy+Pe_PV+Pe_WP+Pe_CHP+Pe_EESdisc-Pe_EESc-Pe_EB-Pe_CCS==Load_E;%电功率平衡
Ph_CHP+Ph_GB+Ph_EB-Ph_CCS==Load_H; %热功率平衡
0<=Pe_buy;%保证电力不外送，政策要求小机组就地消纳不上网
0.2*Q_CHP<=Pe_CHP<=Q_CHP; %CHP容量约束
0.1*Q_EB<=Ph_EB<=Q_EB;%电锅炉容量约束
0.2*Q_GB<=Ph_GB<=Q_GB;%燃气炉容量约束
0<=Qe_EES<=Q_EES;
Qe_EES(1)==0.2*Q_EES;%给定一个初始容量
Qe_EES(25)==0.2*Q_EES;
0<=Pe_EESc<=0.2*Q_EES;%%储能爬坡约束，因为调度时长1h比较长，其他设备不考虑爬坡约束
0<=Pe_EESdisc<=0.2*Q_EES;
0<=Q_CO2capture<=CCS_max;%碳捕集约束，不超过最大捕集能力
alpha*Q_CO2product<=Q_CO2in<=Q_CO2product;%碳捕集待捕集比例
];
for i=1:24
Cons = Cons +[
Qe_EES(i+1)==0.98*Qe_EES(i)+(Pe_EESc(i)*0.9- Pe_EESdisc(i)/0.9);%蓄电池储能
];
end
%% 目标函数 系统总运行成本最小 包括购电成本、购气成本、碳交易成本、运维成本
Cost_Ebuy=sum(Price_E.*Pe_buy); %购电成本，元
Cost_Gasbuy=Price_Gas*sum(3600*(Qg_CHP+Qg_GB)); %购气成本，元 注意：Qg_CHP+Qg_GB 是每秒的耗气量，要乘以3600是每个小时的 再sum是一天24小时的
%计算碳交易成本，日结 参考《上海市2021年碳排放配额分配方案》 1kwh=0.0036GJ
CHP_quota=0.06885*0.0036*(sum(Ph_CHP)+sum(Pe_CHP));%CHP配额，t 热电比大于1按供热分配:单位配额0.06885t/GJ 年度综合供热量=年度实际供热量+年度发电折算供热量
GB_quota=0.06233*0.0036*sum(Ph_GB);%GB产热配额，t 单位配额0.06233t/GJ
CO2_quota=GB_quota+CHP_quota;%总配额
CO2_emission=0.001*3600*sum(Q_CO2product-Q_CO2capture);%总排放量
Cost_CO2trade=80*(CO2_emission-CO2_quota);%碳交易成本，碳交易价格80元/t，如果实际排放量小于配额，表现为收益
Cost_om=PV_om+WP_om+CHP_om+GB_om+EB_om+EES_om+CCS_om;%所有设备的运维成本
Cost_total=Cost_Ebuy+Cost_Gasbuy+Cost_CO2trade+Cost_om;%总运行成本
%% 开始求解
options=sdpsettings( 'verbose', 2,'fmincon.maxiter',100000,'fmincon.MaxFunEvals',100000);%sdpsettings的设置可以去CSDN上查阅，这样就够用，一般默认选择Cplex求解
sol = optimize(Cons,Cost_total,options);%在这里可以更换“目标函数” 格式：optimize(约束条件，目标函数，求解要求)
if sol.problem == 0
%外购电力功率
Pe_buy=value(Pe_buy);
%CHP
Qg_CHP=value(Qg_CHP);
Ph_CHP=value(Ph_CHP);
Pe_CHP=value(Pe_CHP);
Prg_CHP=value(Prg_CHP);
CHP_omwork=value(CHP_omwork);
CHP_om=value(CHP_om);
%电锅炉
Pe_EB=value(Pe_EB);
Ph_EB=value(Ph_EB);
EB_omwork=value(EB_omwork);
EB_om=value(EB_om);
%燃气锅炉
Qg_GB=value(Qg_GB);
Ph_GB=value(Ph_GB);
GB_omwork=value(GB_omwork);
GB_om=value(GB_om);
%蓄电池
Qe_EES=value(Qe_EES);
Pe_EESc=value(Pe_EESc);
Pe_EESdisc=value(Pe_EESdisc);
Real_EESdis=Pe_EESdisc- Pe_EESc;
EES_om=value(EES_om);
EES_omwork=value(EES_omwork);
%碳捕集设备
Pe_CCS=value(Pe_CCS);
Ph_CCS=value(Ph_CCS);
Q_CO2capture=value(Q_CO2capture);
Q_CO2in=value(Q_CO2in);
Q_CO2product=value(Q_CO2product);
CO2_emission=value(CO2_emission);
%目标函数部分转化成值
Cost_Ebuy=value(Cost_Ebuy)
Cost_Gasbuy=value(Cost_Gasbuy)
Cost_CO2trade=value(Cost_CO2trade)
Cost_om=value(Cost_om)
Cost_total=value(Cost_total)
X=['计算完成'];
disp(X)
else
Y=['出错了'];
disp(Y)
end
%% 开始画图
%先把耗电耗热的取反
Load_E=-Load_E;
Load_H=-Load_H;
Pe_EB=-Pe_EB;
Pe_CCS=-Pe_CCS;
Ph_CCS=-Ph_CCS;
Pe_EESc=-Pe_EESc;
figure(1);
x=1:24;
plot(x,Pe_PV','-r*')
hold on
plot(x,Pe_WP','-bo')
xlabel('时间/h'); ylabel('功率/kw');title('新能源出力曲线');
legend('光伏','风电','Location','East')
figure(2);%电功率平衡
Figure_E_in=[ Pe_buy' Pe_PV' Pe_WP' Pe_CHP' Pe_EESdisc'];
bar(Figure_E_in,'stacked');
hold on
Figure_E_out=[Pe_EESc' Pe_EB' Pe_CCS' Load_E'];
bar(Figure_E_out,'stacked');
xlabel('时间/h'); ylabel('功率/kw');title('电功率平衡');
legend('外购电','光伏','风电','热电联产','储能放电','储能充电','电锅炉耗电','碳捕集耗电','电负荷','Location','NorthWest');
grid on;axis([0 25,-inf,inf])
figure(3);%热功率平衡
Figure_H_in=[Ph_CHP' Ph_GB' Ph_EB' ];
bar(Figure_H_in,'stacked');
hold on
Figure_H_out=[Ph_CCS' Load_H'];
bar(Figure_H_out,'stacked');
xlabel('时间/h'); ylabel('功率/kw');title('热功率平衡');
legend('热电联产','燃气锅炉','电锅炉','碳捕集耗热','热负荷','Location','NorthWest');
grid on;axis([0 25,-inf,inf])
figure(4);
subplot(1,2,1);
pie([Cost_Ebuy,Cost_Gasbuy,Cost_om,Cost_CO2trade]);
title('各部分成本比例');
legend('购电成本','购气成本','运维成本','碳交易成本','Location','South');
subplot(1,2,2);
Cos
