MATLAB 函数画图h=0.1,0.01,0.001

Exercise:

 结果呈现:

 方式:

g = colormap(lines);hold on;
for i = 1:3
    x = 0:power(10,-i):2*pi;
    y = (exp(-x)) .* (sin((x.^2)/2));
    m = diff(y)./diff(x);
    plot(x(1:end-1),m,'Color',g(i,:));
end
hold off;
set(gca,'YLim',[-0.25,0.25]);
set(gca,'XLim',[0, 2*pi]);
set(gca,'XTick',0:pi/2:2*pi);
set(gca,'YTick',-0.2:0.1:0.2);
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'});
h = legend('h=0.1','h=0.01','h=0.001');
box on;

编程两大法宝:

百度 + 再读一遍你写的代码

整理以下代码:%% 含电转气-碳捕集(P2G-CCS)耦合的综合能源系统低碳经济调度 clc clear %% 参数设定 P_C2P=sdpvar(1,24); %CCPP-P2G系统的总能耗 P_P2G=sdpvar(1,24); %P2G设备的能耗 P_CC=sdpvar(1,24); %碳捕集能耗 P_OP=sdpvar(1,24); %碳捕集运行能耗 P_GN=sdpvar(1,24); %碳捕集电厂净出力 P_G=sdpvar(1,24); %碳捕集电厂等效出力 P_GC=sdpvar(1,24); %碳捕集电厂提供的碳捕集能耗 P_Galpha=sdpvar(1,24); %碳捕集电厂的烟气处理能耗 Q_CC=sdpvar(1,24); %CCPP-P2G系统捕集的CO2的总量,单位为吨 Q_P2Gsum=sdpvar(1,24); %P2G消耗的CO2的总量,单位为吨 V_P2G=sdpvar(1,24); %P2G生成的天然气 P_alpha=sdpvar(1,24); %总烟气处理能耗 alpha_1=sdpvar(1,24); %t时刻进行烟气处理的烟气量中由 碳捕集电厂运行产生的烟气提供的部分 alpha_2=sdpvar(1,24); %流入烟气存储罐的气量 alpha_3=sdpvar(1,24); %t时刻进行烟气处理的烟气量中由储气装置提供的烟气量 P_WC=sdpvar(1,24); %风电机组提供的碳捕集能耗 P_VC=sdpvar(1,24); %光伏机组提供的碳捕集能耗 P_WIC=sdpvar(1,24); % 碳捕集电厂中火电部分提供的碳捕集能耗 P_Valpha=sdpvar(1,24); %光伏机组提供的烟气处理能耗 P_Walpha=sdpvar(1,24); %风力机组提供的烟气处理能耗 P_WIalpha=sdpvar(1,24); % 碳捕集电厂中火电部分提供的烟气处理能耗 P_WN=sdpvar(1,24); %风力机组的上网功率 P_VN=sdpvar(1,24); %光伏机组的上网功率 P_WIN=sdpvar(1,24); % 碳捕集电厂的上网功率 P_WI=sdpvar(1,24); % 碳捕集发电出力 Q_N=sdpvar(1,24); %碳捕集电厂的CO2的净排放量,单位为t/h P_PH=sdpvar(1,24); %CHP机组的输出总功率 P_CHP=sdpvar(1,24); %CHP机组的输出电功率 H_CHP=sdpvar(1,24); %CHP机组的输出热功率 V_CHP=sdpvar(1,24); %CHP机组消耗的天然气量 H_GB=sdpvar(1,24); %燃气锅炉的输出热功率 V_GB=sdpvar(1,24); %燃气锅炉消耗的天然气量 S_ES=sdpvar(1,24); %电储能在t时段末的蓄电量,单位为MW P_ESC=sdpvar(1,24); %电储能的充电功率 P_ESD=sdpvar(1,24); %电储能的放电功率 S_TS=sdpvar(1,24); %热储能在t时段末的蓄热量,单位为MW H_TSC=sdpvar(1,24); %热储能的充热功率 H_TSD=sdpvar(1,24); %热储能的放热功率 P_EM=sdpvar(1,24); %系统在电网的购电量 P_CUT=sdpvar(1,24); %各级中断负荷功率之和 P_Cmax=sdpvar(1,24); %碳捕集系统的运行能耗上限 miu_ESC=binvar(1,24); %充电的布尔变量 miu_ESD=binvar(1,24); %放电的布尔变量 miu_TSC=binvar(1,24); %充热的布尔变量 miu_TSD=binvar(1,24); %放热的布尔变量 V_WIalpha=sdpvar(1,24); %烟气存储罐储气量 Q_CS=sdpvar(1,24); %流入烟气存储罐的气量 Q_P2G=sdpvar(1,24); %流入烟气存储罐的气量 P_A=15*ones(1,24); %CCPP-P2G系统固定能耗(因占比较少设为定值),单位为MW %风电机组的预测出力 P_W=[232.75,247.44,219.09,188.78,239.58,232.84,188.52,159.84,111.45,51.23,119.88,137.29,141.39,115.78,135.24,143.44,151.64,195.69,159.70,180.94,203.38,193.64,155.32,247.43]; %光伏机组的预测出力 P_V=[0,0,0,0,0,22,63,97,110,118,128,132,133,136,131,133,120,85,37,0,0,0,0,0]; %电负荷 P_EL=[457,319,296,228,184,297,406,509,607,687,803,857,845,793,832,801,795,731,640,593,554,518,525,409]; %热功率 H_HL=[109,131,158,153,139,121,111,98,82,57,22,12,42,62,89,99,122,131,148,160,139,131,119,74]; %购电价 k_EM=[38.85,39.18,36.89,35.57,39.84,43.77,51.31,64.10,74.59,77.21,85.41,89.02,82.46,80.49,83.11,81.80,78.52,73.93,69.67,76.89,74.26,66.39,55.57,46.72]; S_ES_init=60; S_TS_init=30; P_WA=0.2*P_W; %t时段弃风功率 P_VA=0.2*P_V; %t时段弃光功率 %% 约束条件 C=[]; %CCPP-P2G系统能耗以及CCPP出力 for t=1:24 C=[C, P_C2P(t)==P_P2G(t)+P_CC(t), %CCPP-P2G系统总能耗约束 P_P2G(t)==P_WA(t)+P_VA(t), %P2G消纳的弃风、弃光量约束 P_CC(t)==P_A(t)+P_OP(t), %碳捕集能耗约束 P_GN(t)==P_G(t)-P_GC(t)-P_Galpha(t), %碳捕集电厂功率约束 ]; end %CCPP-P2G系统碳利用量以及天然气生成量 for t=1:24 C=[C, Q_CC(t)==P_OP(t)/0.269, %CCPP-P2G系统捕集的CO2的总量与耗能约束 P_OP是碳捕集运行能耗 Q_P2Gsum(t)==0.2*0.6*P_P2G(t), %P2G设备的消耗的CO2量和电功率约束 V_P2G(t)==3.6*0.6*P_P2G(t)/39, %P2G设备的生成天然气 ]; end % 碳捕集电厂烟气处理模型 for t=1:24 C=[C,P_alpha(t)==0.513*(alpha_1(t)+alpha_3(t)),]; %烟气处理系统能耗 end %碳捕集-风电-光伏联合运行策略 for t=1:24 C=[C, P_GC(t)+P_WC(t)+P_VC(t)+P_WIC(t)==P_CC(t), %碳捕集能耗等式约束 P_OP(t)==0.269*Q_CC(t), %CCPP-P2G系统捕集的CO2的总量与耗能约束 P_Valpha(t)+P_Walpha(t)+P_Galpha(t)+P_WIalpha(t)==P_alpha(t), %烟气处理能耗等式约束 P_WN(t)+P_WC(t)+P_Walpha(t)==P_W(t), %风力机组的出力约束 P_VN(t)+P_VC(t)+P_Valpha(t)==P_V(t), %光伏机组的出力约束 P_WIN(t)+P_WIC(t)+P_WIalpha(t)==P_WI(t), % 碳捕集电厂的出力约束 Q_N(t)==0.96*P_G(t)-Q_CC(t), %碳捕集电厂的碳排放约束 ]; end %CHP机组和燃气锅炉模型 for t=1:24 C=[C, P_PH(t)==P_CHP(t)+H_CHP(t), %CHP机组的输出功率约束 P_CHP(t)==V_CHP(t)*39*0.35, %CHP机组的输出电功率约束 H_CHP(t)==V_CHP(t)*39*0.40, %CHP机组的输出热功率约束 H_GB(t)==V_GB(t)*39*0.40, %CHP机组的输出热功率约束 ]; end %储能装置模型 for t=2:24 C=[C, S_ES(t)==S_ES(t-1)*(1-0.001)+0.95*P_ESC(t)-P_ESD(t)/0.95, %电储能运行约束 S_TS(t)==S_TS(t-1)*(1-0.01)+0.88*H_TSC(t)-H_TSD(t)/0.88, %热储能运行约束 ]; end %电功率和热功率平衡约束 for t=1:24 C=[C, P_GN(t)+P_WIN(t)+P_CHP(t)+P_WN(t)+P_VN(t)+P_ESD(t)+P_EM(t)==P_EL(t)+P_ESC(t), %电功率平衡约束 H_CHP(t)+H_GB(t)+H_TSD(t)==H_HL(t)+H_TSC(t), %热功率平衡约束 ]; end %碳捕集电厂约束 for t=1:24 C=[C, 100<=P_G(t)<=400, %碳捕集电厂出力上下限约束 0<=Q_CC(t)<=0.96*400, 15<=P_GC(t)+P_WC(t)+P_VC(t)+P_WIC(t)<=P_Cmax(t), %碳捕集系统的运行能耗上下限 P_Cmax(t)==0.269*0.96*P_G(t), %碳捕集系统的运行能耗上限赋值 ]; end for t=2:24 C=[C, -60<=P_G(t)-P_G(t-1)<=60, %碳捕集电厂出力爬坡速率约束 -65<=P_GC(t)+P_WC(t)+P_VC(t)+P_WIC(t)-P_GC(t-1)-P_WC(t-1)-P_VC(t-1)-P_WIC(t-1)<=65, %碳捕集电厂碳捕集能耗爬坡速率约束 ]; end %CHP机组电热出力以及爬坡约束 for t=1:24 C=[C, 0<=P_CHP(t)<=140, %CHP的电功率出力约束 0<=H_CHP(t)<=160, %CHP的热功率出力约束 ]; end for t=2:24 C=[C,-120<=P_PH(t)-P_PH(t-1)<=120,]; %CHP机组出力爬坡速率约束 end %燃气锅炉热出力以及爬坡约束 for t=1:24 C=[C,0<=H_GB(t)<=200,]; %燃气锅炉的热功率出力约束 end for t=2:24 C=[C,-100<=H_GB(t)-H_GB(t-1)<=100,]; %燃气锅炉的爬坡速率约束 end %P2G运行约束 C=[C,0<=P_P2G<=200,]; %P2G运行功率上下限 %电储能和热储能约束 for t=1:24 C=[C, 0<=P_ESC(t)<=miu_ESC(t)*80, %充电上下限约束 0<=P_ESD(t)<=miu_ESD(t)*80, %放电上下限约束 0<=miu_ESC(t)+miu_ESD(t)<=1,%不可能同时出现充放电 10<=S_ES(t)<=120, %储电量上下限约束 0<=H_TSC(t)<=miu_TSC(t)*40, %充热上下限约束 0<=H_TSD(t)<=miu_TSD(t)*40, %放热上下限约束 0<=miu_TSC(t)+miu_TSD(t)<=1,%不可能同时出现充放热 10<=S_TS(t)<=50, %储热量上下限约束 ]; end C=[C,S_ES_init==S_ES(24),S_ES(1)==S_ES_init*(1-0.001)+0.95*P_ESC(1)-P_ESD(1)/0.95,]; %电储能和热储能初始和结束的约束 C=[C,S_TS_init==S_TS(24),S_TS(1)==S_TS_init*(1-0.01)+0.88*H_TSC(1)-H_TSD(1)/0.88,]; % 碳捕集电厂约束 C=[C,sum(P_WI)<=1500,]; % 碳捕集电厂日总出力上限 for t=1:24 C=[C, 60<=P_WI(t)<=100, % 碳捕集电站的出力上下限 0.1*400<=V_WIalpha(t)<=0.9*400, %储气装置容量的上下限 0<=alpha_1(t)<=160, 0<=alpha_2(t)<=160, %流入储气装置的上下限 0<=alpha_3(t)<=160, %流出储气装置的上下限 ]; end for t=2:24 C=[C,-40<=P_WI(t)-P_WI(t-1)<=40,]; % 碳捕集电站的爬坡速率约束(文中未给参数,暂取40) end for t=2:24 C=[C,V_WIalpha(t)==V_WIalpha(t-1)+alpha_2(t)-alpha_3(t),]; %储气装置的容量变化约束 end C=[C,V_WIalpha(1)==0.4*400+alpha_2(1)-alpha_3(1),]; C=[C,Q_CC==Q_CS+Q_P2G,Q_CS>=0,Q_P2G>=0,0<=P_VC,0<=P_WIC,0<=P_GC,]; C=[C,P_EM>=0,Q_P2Gsum-Q_P2G>=0,P_WIN>=0,0<=P_Valpha,0<=P_Walpha,0<=P_Galpha,0<=P_WIalpha,0<=P_WC,]; %% 目标函数 C_F=200*24+17*sum(P_G)+0.04*P_G*(P_G');%碳捕集电厂燃料成本C_F C_WI=19.8*(0.96-0.76)*sum(P_WI); % 碳捕集电厂成本 V_BUY=V_CHP+V_GB-V_P2G; %天然气购买量 C_H=0.419*sum(V_BUY); %CHP机组和燃气锅炉成本 Q_BUY=Q_P2Gsum-Q_P2G; %购买的二氧化碳的量 C_PG=120*sum(Q_BUY)+20*sum(P_P2G); %P2G成本 C_CS=4.89*sum(Q_CS); %碳封存成本 C_W=21.4*sum(P_W)+14.2*sum(P_V); %系统运行维护成本 C_M=sum(k_EM.*P_EM); %电力市场购电成本 %碳排放权配额模型 E_e_buy=0.68*C_M; %购电配额 E_CHPGB=0.102*3.6*sum(V_BUY); %CHP和GB配额 Q_Q=0.76*P_GN; %碳捕集电厂的碳排放配额 z=sum(V_BUY); E_IES=E_e_buy+E_CHPGB+sum(Q_Q); %IES总碳排放配额 %实际碳排放模型 E_e_buy_a=0.711*C_M; E_CHPGB_a=0.249*3.6*sum(V_BUY); E_IES_a=E_e_buy_a+E_CHPGB_a+sum(Q_N); E=E_IES-E_IES_a %实际IES总碳排放 %阶梯碳交易成本(分段线性化) E_v=sdpvar(1,5)%每段区间内的长度,分为5段,每段长度是2000 lamda=19.8;%碳交易基价 C=[C, E==sum(E_v),%总长度等于E 0<=E_v(1:4)<=2000,%除了最后一段,每段区间长度小于等于1000 0<=E_v(5), ]; %碳交易成本 C_CO2=0; for v=1:5 C_CO2=C_CO2+(lamda+(v-1)*0.25*lamda)*E_v(v); end Obj=C_F-C_CO2+C_WI+C_H+C_PG+C_CS+C_W+C_M; %% 求解 ops=sdpsettings('solver','cplex','verbose',2,'usex0',0); ops.cplex.mip.tolerances.mipgap=1e-6; result=optimize(C,Obj,ops); if result.problem == 0 % problem =0 代表求解成功 else error('求解出错!'); end %% 数据分许与画图部分 %碳捕集能耗 P_CC=-P_CC; P_WC=double(P_WC); P_VC=double(P_VC); P_WIC=double(P_WIC); P_CC=double(P_CC); for t=1:24 Plot_p(1,t)=P_WC(t); Plot_p(2,t)=P_VC(t); Plot_p(3,t)=P_WIC(t); end Plot_p=Plot_p'; figure b=bar(Plot_p,'stacked'); b(1).FaceColor = [.9 .8 .1]; b(2).FaceColor = [.0 .9 .0]; b(3).FaceColor = [.9 .1 .3]; hold on bar(P_CC,'stacked'); xlabel('时段'); ylabel('功率/MW'); legend('风电提供能耗','光伏提供能耗',' 碳捕集电厂提供能耗','碳捕集总能耗'); legend('boxoff'); title('综合能源系统的碳捕集能耗'); box off %CO2排放量和处理量 Q_N=double(Q_N); Q_CC=double(Q_CC); figure bar([Q_N',Q_CC']) xlabel('时段'); ylabel('二氧化碳量/t'); legend('二氧化碳排放量','二氧化碳捕集量'); legend('boxoff'); title('CO2排放量和处理量'); box off % 碳捕集电厂出力和能量市场购电量 P_WI=double(P_WI); P_WIN=double(P_WIN); P_EM=double(P_EM); figure d=bar(P_EM,'stacked'); d(1).FaceColor = [.9 .1 .3]; hold on plot(P_WI,'b-.','LineWidth',1.5); hold on plot(P_WIN,'m-.o','LineWidth',1.5); xlabel('时段'); ylabel('功率/MW'); legend('综合能源系统购电量',' 碳捕集电厂等效出力',' 碳捕集电厂净出力'); legend('boxoff'); title('综合能源系统购电和电厂出力情况'); box off % 电负荷和各供电单元出力 % P_GN(t)+P_WIN(t)+P_CHP(t)+P_WN(t)+P_VN(t)+P_ESD(t)+P_EM(t)==P_EL(t)+P_ESC(t), %电功率平衡约束 P_GN=double(P_GN); P_WIN=double(P_WIN); P_CHP=double(P_CHP); P_WN=double(P_WN); P_VN=double(P_VN); P_ESD=double(P_ESD); P_EM=double(P_EM); P_EL=double(P_EL); P_ESC=double(P_ESC); for t=1:24 Plot_PLoad(1,t)=P_GN(t); Plot_PLoad(2,t)=P_WIN(t); Plot_PLoad(3,t)=P_CHP(t); Plot_PLoad(4,t)=P_WN(t); Plot_PLoad(5,t)=P_VN(t); Plot_PLoad(6,t)=P_ESD(t); Plot_PLoad(7,t)=P_EM(t); Plot_PLoad(8,t)=-1*P_ESC(t); Plot_PLoad(9,t)=-1*P_EL(t); end Plot_PLoad=Plot_PLoad'; figure f1=bar(Plot_PLoad,'stacked'); f1(1).FaceColor = [.9 .1 .3]; f1(2).FaceColor = [.2 .7 .2]; f1(3).FaceColor = [.2 .2 .7]; f1(4).FaceColor = [.9 .8 .1]; f1(5).FaceColor = [.6 .8 .1]; hold on plot(P_EL,'m-.s','LineWidth',2); xlabel('时间/h'); ylabel('功率/MW'); title('电负荷和各机组供电出力'); hold on legend('碳捕集电厂净出力','碳捕集电厂的出力','CHP机组电出力','风力机组出力',... '光伏机组出力','电储能的放电功率','电网的购电量','电储能的充电功率','电负荷功率'); legend('boxoff'); grid on box off % 热负荷和各供热单元出力 % H_CHP(t)+H_GB(t)+H_TSD(t)==H_HL(t)+H_TSC(t), %热功率平衡约束 H_CHP=double(H_CHP); H_GB=double(H_GB); H_TSD=double(H_TSD); H_TSC=double(H_TSC); P_CHP=double(P_CHP); for t=1:24 Plot_HeatLoad(1,t)=H_CHP(t); Plot_HeatLoad(2,t)=H_GB(t); Plot_HeatLoad(3,t)=H_TSD(t); Plot_HeatLoad(4,t)=-1*H_TSC(t); Plot_HeatLoad(5,t)=-1*H_HL(t); end Plot_HeatLoad=Plot_HeatLoad'; figure f=bar(Plot_HeatLoad,'stacked'); f(1).FaceColor = [.9 .1 .3]; f(2).FaceColor = [.2 .7 .2]; f(3).FaceColor = [.2 .2 .7]; f(4).FaceColor = [.9 .8 .1]; f(5).FaceColor = [.6 .8 .1]; hold on plot(H_HL,'m-.s','LineWidth',2); xlabel('时间/h'); ylabel('功率/MW'); title('热负荷和各机组供热出力'); hold on legend('CHP机组热出力','燃气锅炉热出力','热储能放热','热储能储热','热负荷功率'); legend('boxoff'); grid on box off % 电储能 P_ESD=-P_ESD; P_ESC1=double(P_ESC);P_ESD1=double(P_ESD);S_ES=double(S_ES); for t=1:24 P(1,t)=P_ESC1(t); P(2,t)=S_ES(t); end P=P'; figure a=bar(P,'stacked'); a(1).FaceColor = [.2 .7 .2]; a(2).FaceColor = [.9 .1 .3]; hold on bar(P_ESD1,'stacked'); hold on plot(S_ES,'c-.x','LineWidth',1.5); xlabel('时间/h'); ylabel('功率/MW'); title('IES电储能装置功率状态'); hold on legend('充电功率','电储能蓄电量','放电功率'); legend('boxoff'); grid on box off % 热储能 H_TSD=-H_TSD; S_TS=double(S_TS);H_TSC=double(H_TSC);H_TSD=double(H_TSD); for t=1:24 R(1,t)=H_TSC(t); R(2,t)=S_TS(t); end R=R'; figure bar(R,'stacked'); hold on a=bar(H_TSD,'stacked'); a(1).FaceColor = [.9 .1 .3]; hold on plot(S_TS,'m-.o','LineWidth',1.5); xlabel('时间/h'); ylabel('功率/MW'); title('IES热储能装置功率状态'); hold on legend('充热功率','放热功率','热储能储热量'); legend('boxoff'); grid on box off
09-02
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值