109号资源-源程序：论文可在知网下载《综合能源系统低碳优化调度》本人博客有解读

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 E=0; 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将光照强度：w/m^2转换成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); %2000=2220*（12-3）因为上限容量为2MW，所以风速大于12时满额输出 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 sdpvar-创建实数型决策变量，等待求解 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);%电锅炉电功率 例如:x = sdpvar(1,1); %x为1×1的实数变量 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 1h=3600s 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小时的 %计算碳交易成本，日结 参考《上海市2022年碳排放配额分配方案》 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求解 %verbose：展示求解细节的设置。0表示完全不显示，1表示适度显示，2则是完全显示。 % MaxFunEvals:允许进行函数评价的最大次数,取值为正整数. Maxiter: 允许进行迭代的最大次数,取值为正整数 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','-ro','MarkerFaceColor','r') %光伏出力 -表示线条 r红色 *点的形状 hold on plot(x,Pe_WP','-bo') %风电出力 xlabel('时间/h'); ylabel('功率/kw');title('24h风电和光伏出力曲线'); legend('光伏','风电','Location','East') %图例 grid on figure(2); %电功率平衡 Figure_E_in=[ Pe_buy' Pe_PV' Pe_WP' Pe_CHP' Pe_EESdisc']; %产电量 bar(Figure_E_in,'stacked'); %stack - 将来自输入表或时间表的数据堆叠到输出表或时间表的一个变量中 hold on Figure_E_out=[Pe_EESc' Pe_EB' Pe_CCS' Load_E' ]; %耗电量 a=bar(Figure_E_out,'stacked'); a(3).FaceColor = [.1 .7 .3]; a(1).FaceColor = [.0 .0 .0]; a(4).