✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
逐步提升风电等可再生能源发电占比,并对火电机组进行低碳化改造,辅之以多类需求侧资源,是实现能源电力碳达峰,碳中和目标的重要手段.首先,挖掘源荷两侧低碳资源并分析其低碳特性,源侧在碳捕集电厂中装设烟气旁路系统与溶液存储器,形成碳捕集电厂综合灵活运行方式进而与风电协调配合;荷侧调用不同响应速度的价格型,激励型需求响应资源克服多时间尺度下碳捕集电厂综合灵活运行方式的局限,通过源荷资源协调优化,从而提高系统的低碳性能.其次,构建源荷协调的日前—日内—实时3阶段低碳经济调度模型,优化系统的负荷及旋转备用分配计划,并改善失负荷与弃风问题.最后,在改进的IEEE-39节点系统中进行算例分析,结果表明该调度方法能够利用源荷可调节资源的调度优势,实现电力系统低碳经济调度的目标.
⛄ 部分代码
%% 日前调度
clc
clear
close all
warning off % 关闭版本不兼容提示
%% 数据导入
% 调度周期
T=24;
% 风电预测功率
P_prew=1.4*[299.316845700000,267.373698500000,280.085373600000,284.772683200000,329.099586800000,289.710754300000,168.284358400000,159.064053800000,186.034614400000,198.132934700000,131.119251700000,85.7752555300000,58.5152130000000,91.2575791600000,116.950628000000,139.509689700000,149.377514500000,191.110973000000,267.212495800000,291.126903600000,300.994923900000,334.319139100000,353.357157400000,342.351496400000];
% 电负荷
EleLoad=0.74*[1209.12680000000,1391.49080000000,1417.72400000000,1431.18080000000,1398.85200000000,1322.37000000000,1410.94800000000,1492.84800000000,1554.33600000000,1390.15800000000,1770.80400000000,1786.93200000000,1512.25200000000,1305.10800000000,1210.39800000000,1228.51400000000,1440.85200000000,1770.42600000000,1969.12800000000,1685.12400000000,1488.74600000000,1445.37400000000,1286.13800000000,1377.11000000000];
% 热负荷
HeatLoad=0.7*[358.676650700000,358.415415000000,370.788488900000,369.577305000000,344.641167600000,321.889910800000,352.311998500000,326.544656500000,315.744330800000,265.474819400000,232.010437700000,227.021414800000,296.575979700000,292.170468200000,267.079475800000,273.994086500000,288.173414800000,270.458472200000,227.607740900000,345.550477800000,322.246141400000,342.399073100000,347.968033700000,370.337263500000];
% 分时电价
gama_elec=100*[0.29 0.29 0.29 0.29 0.29 0.69 0.69 0.69 0.84 0.84 0.84 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.64 0.64 0.55 0.55 0.29 0.29];
%% 相关参数
a=[0.00048,0.00031,0.0002]; % 成本系数a
b=[16.2,17.3,16.6]; % 成本系数b
c=[1000,970,700]; % 成本系数c
eg=[0.91,0.95,0.98]; % 碳排放系数
GT_P_max=[150,80];
E_bata=0.81; % 碳捕集效率tu
P_yita=[1.05,1.05,1.05]; % 再生塔和压缩机最大工作状态系数
P_lamda=[0.269,0.269,0.269]; % 单位捕碳量能耗
P_G_max=[400,455,200]; % 机组最大出力
P_G_min=[200,120,100]; % 机组最小出力
P_D1=0.5*[20,20,15]; % 碳捕集固定能耗
P_D=ones(3,24);
for i=1:3
P_D(i,:)=P_D1(i).*ones(1,24);
end
R_u=4*[50,50,25]; % 机组上爬坡
M_MEA=61.08; % M_MEA的摩尔质量
M_co2=44; % 二氧化碳的摩尔质量
sita=0.4; % 再生塔解析量
CR=30; % 醇胺溶液浓度(%)
rou_R=1.01; % 醇胺溶液密
V_CR=[2000,2000,2000]; % 溶液储液装置容量
K_R=1.17; % 乙醇胺溶剂成本系数
fai=1.5; % 溶剂运行损耗系数
% h
h_min=[200,120,100];
h_max=[401,456,201];
h_gn=5;
hl1=(h_max-h_min)./h_gn;
hl2=zeros(3,h_gn+1);
% k
k_min=[-400,-455,-200];
k_max=[-199,-119,-99];
k_gn=5;
kl1=(k_max-k_min)./k_gn;
kl2=zeros(3,k_gn+1);
for i=1:3
hl2(i,:)=h_min(i):hl1(i):h_max(i);
end
for i=1:3
kl2(i,:)=k_min(i):kl1(i):k_max(i);
end
%% 决策变量
% 电力源出力
GT_P = sdpvar(2,24,'full'); % 燃气轮机电出力
P_w = sdpvar(1,24,'full'); % 风电机组出力
P_G = sdpvar(3,24,'full'); % 火电机组出力
EB=sdpvar(2,24,'full'); % 电锅炉出力
% 热力源出力
GT_H = sdpvar(2,24,'full'); % 燃气轮机热出力
EB_H=sdpvar(2,24,'full'); % 电锅炉热出力
% 天然气
P_gas=sdpvar(2,24,'full'); % 天然气需求
% 碳捕集相关
E_G=sdpvar(3,24,'full'); % 碳捕集机组产生的总碳排放
E_total_co2=sdpvar(3,24,'full'); % 机组捕获的总碳排放
E_CG=sdpvar(3,24,'full'); % 储液装置提供的待捕集二氧化碳量
P_B=sdpvar(3,24,'full'); % 机组运行能耗
P_J=sdpvar(3,24,'full'); % 机组净出力
V_CA=sdpvar(3,24,'full'); % 机组净出力
V_FY=sdpvar(3,24,'full'); % 富液体积
V_PY=sdpvar(3,24,'full'); % 贫液体积
P_tran=sdpvar(1,24,'full'); % 系统可转移电负荷
P_cut=sdpvar(1,24,'full'); % 系统可削减电负荷
P_DE=sdpvar(1,24,'full'); % 系统经过过需求响应后的电负荷
H_tran=sdpvar(1,24,'full'); % 系统可转移热负荷
H_cut=sdpvar(1,24,'full'); % 系统可削减热负荷
H_DE=sdpvar(1,24,'full'); % 系统经过过需求响应后的热负荷
gn=5;
P_G_line= sdpvar(3,24,'full'); % 火电机组出力
w1=sdpvar(gn+1,T,'full');
w2=sdpvar(gn+1,T,'full');
w3=sdpvar(gn+1,T,'full');
w4=sdpvar(gn+1,T,'full');
w5=sdpvar(gn+1,T,'full');
z1=binvar(gn, T,'full');
z2=binvar(gn, T,'full');
z3=binvar(gn, T,'full');
z4=binvar(gn, T,'full');
z5=binvar(gn, T,'full');
y1= sdpvar(3,24,'full');
h_w1=sdpvar(h_gn+1,T,'full');
h_w2=sdpvar(h_gn+1,T,'full');
h_w3=sdpvar(h_gn+1,T,'full');
h_w4=sdpvar(h_gn+1,T,'full');
h_w5=sdpvar(h_gn+1,T,'full');
h_z1=binvar(h_gn,T,'full');
h_z2=binvar(h_gn,T,'full');
h_z3=binvar(h_gn,T,'full');
h_z4=binvar(h_gn,T,'full');
h_z5=binvar(h_gn,T,'full');
y2= sdpvar(3,24,'full');
k_w1=sdpvar(k_gn+1,T,'full');
k_w2=sdpvar(k_gn+1,T,'full');
k_w3=sdpvar(k_gn+1,T,'full');
k_w4=sdpvar(k_gn+1,T,'full');
k_w5=sdpvar(k_gn+1,T,'full');
k_z1=binvar(k_gn, T,'full');
k_z2=binvar(k_gn, T,'full');
k_z3=binvar(k_gn, T,'full');
k_z4=binvar(k_gn, T,'full');
k_z5=binvar(k_gn, T,'full');
%% 约束条件
C = []; %约束条件初始
for t=1:24
for i=1:3
C = [C,
0<=E_CG(i,t),
0<=P_B(i,t),
E_G(i,t)==eg(i)*P_G(i,t), % 碳捕集机组产生的总碳排放
E_total_co2(i,t)==E_CG(i,t)+0.25*E_bata*eg(i)*(y1(i,t)-y2(i,t)), % 机组捕获的二氧化碳总量
0<=E_total_co2(i,t)<=P_yita*E_bata*eg(i)*P_G_max(i),
P_B(i,t)==P_lamda(i)*E_total_co2(i,t), % 机组运行能耗
P_G(i,t)==P_J(i,t)+P_D(i,t)+P_B(i,t), % 机组输出总功率
P_G_min(i)-P_lamda(i)*P_yita*E_bata*eg(i)*P_G_max(i)-P_D(i)<=P_J(i,t)<=P_G_max(i)-P_D(i), % 碳捕集电厂净出力范围
0<= P_w(t)<= P_prew(t), % 风电出力区间约束
sum(EB(:,t))+P_w(t)<=P_prew(t);
P_G_min(i)<= P_G(i,t)<=P_G_max(i), % 火电机组出力约束
];
end
end
C=[C,min(sum(R_u),sum(P_G_max)-sum(P_G))>=0.01*max(P_DE),]; % 旋转备用约束
C = [C, y1(1,:)==hl2(1,:).^2*h_w1];
C = [C, y1(2,:)==hl2(2,:).^2*h_w2];
C = [C, y1(3,:)==hl2(3,:).^2*h_w3];
C = [C, h_w1(1,:)<=h_z1(1,:)];
for i=2:h_gn
C = [C, h_w1(i,:)<=h_z1(i-1,:)+h_z1(i,:)];
end
C = [C, h_w1(h_gn+1,:)<=h_z1(h_gn,:)];
C = [C, sum(h_z1)==ones(1,T)];
C = [C, h_w2(1,:)<=h_z2(1,:)];
for i=2:h_gn
C = [C, h_w2(i,:)<=h_z2(i-1,:)+h_z2(i,:)];
end
C = [C, h_w2(h_gn+1,:)<=h_z2(h_gn,:)];
C = [C, sum(h_z2)==ones(1,T)];
C = [C, h_w3(1,:)<=h_z3(1,:)];
for i=2:h_gn
C = [C, h_w3(i,:)<=h_z3(i-1,:)+h_z3(i,:)];
end
C = [C, h_w3(h_gn+1,:)<=h_z3(h_gn,:)];
C = [C, sum(h_z3)==ones(1,T)];
C = [C, h_w4(1,:)<=h_z4(1,:)];
for i=2:h_gn
C = [C, h_w4(i,:)<=h_z4(i-1,:)+h_z4(i,:)];
end
C = [C, h_w4(h_gn+1,:)<=h_z4(h_gn,:)];
C = [C, sum(h_z4)==ones(1,T)];
C = [C, h_w5(1,:)<=h_z5(1,:)];
for i=2:h_gn
C = [C, h_w5(i,:)<=h_z5(i-1,:)+h_z5(i,:)];
end
C = [C, h_w5(h_gn+1,:)<=h_z5(h_gn,:)];
C = [C, sum(h_z3)==ones(1,T)];
C = [C, y2(1,:)==kl2(1,:).^2*k_w1];
C = [C, y2(2,:)==kl2(2,:).^2*k_w2];
C = [C, y2(3,:)==kl2(3,:).^2*k_w3];
C = [C, k_w1(1,:)<=k_z1(1,:)];
for i=2:k_gn
C = [C, k_w1(i,:)<=k_z1(i-1,:)+k_z1(i,:)];
end
Plot_EleNet=Plot_EleNet';
figure
bar(Plot_EleNet,'stacked');
hold on
plot(EleLoad,'k-','LineWidth',1.5);
plot(P_DE,'r-o','LineWidth',1.5);
xlabel('时间/h');
ylabel('功率/kW');
title('电功率平衡');
hold on
legend('火电机组1','火电机组2','火电机组3','风电','燃气轮机1','燃气轮机2','电锅炉','电负荷','P-DE');
box off
% 热力系统
for t=1:24
Plot_HeatNet(1,t)=sum(GT_H(:,t));
Plot_HeatNet(2,t)=sum(EB_H(:,t));
end
Plot_HeatNet=Plot_HeatNet';
figure
bar(Plot_HeatNet,'stacked');
hold on
plot(HeatLoad,'k-','LineWidth',1.5);
plot(H_DE,'r-o','LineWidth',1.5);
xlabel('时间/h');
ylabel('功率/kW');
title('热功率平衡');
hold on
legend('燃气轮机','电锅炉','热负荷','H-DE');
box off
% 风电出力
figure
plot(P_w+sum(EB),'g-o','LineWidth',1.5)
hold on
plot(P_prew,'k-','LineWidth',1.5)
xlabel('时间/h');
ylabel('功率/kW');
title('风电功率');
% 碳捕集能耗
figure
plot(sum(P_B),'g-^','LineWidth',1.5)
hold on
xlabel('时间/h');
ylabel('功率/kW');
title('碳捕集能耗');
box off
% 火电机组净出力
figure
plot(sum(P_J),'r-o','LineWidth',1.5)
hold on
xlabel('时间/h');
ylabel('功率/kW');
title('火电机组净出力');
box off
% 二氧化碳捕集量
figure
plot(E_total_co2(1,:),'r-','LineWidth',1.5)
hold on
plot(E_total_co2(2,:),'g-^','LineWidth',1.5)
plot(E_total_co2(3,:),'b-*','LineWidth',1.5)
legend('火电机组1','火电机组2','火电机组3')
xlabel('时间/h');
ylabel('碳捕集量');
title('不同机组二氧化碳捕集量');
⛄ 运行结果
⛄ 参考文献
[1]崔杨, 邓贵波, 曾鹏,等. 计及碳捕集电厂低碳特性的含风电电力系统源-荷多时间尺度调度方法[J]. 中国电机工程学报, 2022, 42(16):18.
扫一扫
1840
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
