微网优化调度 | Matlab实现计及绿证及碳排放的含智能楼宇微网优化调度模型

算法如诗

已于 2023-04-24 09:11:22 修改

阅读量382

点赞数

分类专栏：微网优化调度文章标签：微网优化调度人工智能算法 matlab

于 2023-04-22 21:05:55 首次发布

本文链接：https://blog.csdn.net/qq_59771180/article/details/130310375

版权

微网优化调度专栏收录该内容

1 篇文章 0 订阅

订阅专栏

文章目录

效果一览

文章概述

微网优化调度 | Matlab实现计及绿证及碳排放的含智能楼宇微网优化调度模型
计及绿证交易及碳排放的含智能楼宇微网优化调度
代码主要做的是电热综合能源系统的优化调度问题，在传统的含风光储火的微网基础上，加入电动汽车以及智能楼宇单元，组成更加复杂的微网调度模型，并考虑微网参与碳市场以及绿色证书交易市场，并对交易结果进行了量化。结果非常清晰，出图效果也非常清楚，代码非常精品，注释几乎一行一注释；

程序设计

完整代码下载地址：Matlab实现计及绿证及碳排放的含智能楼宇微网优化调度模型

%% 加入绿色证书交易，消纳比例高！

%% 程序初始化
clc
clear
close all

%% 可转移负荷约束
for t=1:24
    C=[C,
    ushiftl(t)+ushiftq(t)<=1;
    shiftl(t)<= ushiftl(t)*pload(t)*0.1;
    shiftq(t)<=ushiftq(t)*pload(t)*0.1;
    0<=shiftl(t);
    0<=shiftq(t);
    ];
end

C=[C,sum(shiftl)==sum(shiftq)];%转移前后总负荷应当保持不变

%% 碳排放及碳交易约束
for t=1:24
    C=[C,
        eq(t)==nc*(pmt(t)+ppv(t)+wt(t)),%碳配额计算
        ];
end

for t=1:24
    C=[C,
        ep(t)==qg*(pmt(t)+pmgb(t)),%碳排放计算
        ];
end

%% 需求响应资源——可中断负荷约束
for m=1:3
    for t=1:24
        C=[C,
        0<=pil(m,t)<=cil(m)*pload(t),   %各级中断负荷约束
        ];
    end
end

for m=1:3
    for t=2:24
      C=[C,
    pil(m,t)+ pil(m,t-1)<=0.2*pload(t), %连续性中断负荷约束  
    ]  ;
    end
end

%% 燃气轮机出力约束
for t=1:24
    C=[C,
     xconv(t)*gtmin<=pmt(t)<=xconv(t)*gtmax ,%出力上下限约束
    ];
end
C=[C,pmt(1)<=ramp]; %初始爬坡率约束
C=[C,xconv(1)<=yconv(1)];%初始工作状态约束

for t=2:24
    C=[C,
    -ramp<=pmt(t)-pmt(t-1)<=ramp,%爬坡率约束
    xconv(t)-xconv(t-1)<=yconv(t), %工作状态约束 
    ];
end


%% 购售电量约束
for t=1:24
    C=[C,
       0<=umob(t)+umos(t)<=1,%状态互斥约束
       0<=pmgb(t)<=umob(t)*pmgmax,%购电量约束
       0<=pmgs(t)<=umos(t)*pmgmax,%售电量约束
      ];
end

%% 储能约束
for t=1:24
    C=[C,
        0<=gesc(t)<=gescmax,%充电功率约束
        0<=gesd(t)<=gesdmax,%放电功率约束
        0<=sess(t)<=sessmax,%蓄电量约束
        ];
end

C=[C,sess(1)==gesc(1)*uesc-gesd(1)/uesd];%初始蓄电量约束

for t=2:24
    C=[C,
       sess(t)==sess(t-1)+gesc(t)*uesc-gesd(t)/uesd,%蓄电量等式约束
       ];
end

%% 功率平衡约束
for t=1:24
    C=[C,
      gesc(t)+pload(t)-shiftl(t)+shiftq(t)-sum(pil(:,t))+pmgs(t)+gcvb(t)*200+gcvr(t)*200==wt(t)+gesd(t)+ppv(t)+pmgb(t)+pmt(t)+gdvr(t)*200+gdvb(t)*200  
    ];
end

for t=1:24
    C=[C,
     cshift(t)==50*shiftl(t)+30*shiftq(t) %计算转移负荷的补偿费用
    ];
end

%% 绿色证书交易
for t=1:24
    C=[C,
    etotal(t)==pload(t)-shiftl(t)+shiftq(t)-gesd(t)+gesc(t)+gcvb(t)*200+gcvr(t)*200-gdvr(t)*200-gdvb(t)*200;  
    Ere(t)==wt(t)+ppv(t);
    Ere(t)+Qgc(t)>=0.125*etotal(t);
    Cgc(t)==Qgc(t)*116.6
    ];
end

%% 电动汽车约束-比亚迪电动汽车以及日产电动汽车
for t=1:24
    C=[C,
    57/1000*0.15<=svb(t)<=57/1000*0.95,%蓄电池容量约束
    0<=gcvb(t)<=57/1000*0.2*ucvb(t),%充电功率约束
    0<=gdvb(t)<=57/1000*0.2*udvb(t),%放电功率约束
    ucvb(t)+udvb(t)<=1,%状态变量约束
    
    24/1000*0.15<=svr(t)<=24/1000*0.95,%同上
    0<=gcvr(t)<=24/1000*0.2*ucvr(t),%同上
    0<=gdvr(t)<=24/1000*0.2*udvr(t),%同上
    ucvr(t)+udvr(t)<=1,%同上
    ];
end
%电动汽车蓄电池蓄电量等式约束
for t=2:23
    C=[C,
    svr(t)==svr(t-1)+0.9*gcvr(t)-gdvr(t)/0.9-dvr(t)*0.228/1000  ;  
    svb(t)==svb(t-1)+0.9*gcvb(t)-gdvb(t)/0.9-dvb(t)*0.229/1000;
    ];
end
%初始蓄电量约束
C=[C,svb(1)==36.42/1000+0.9*gcvb(1)-gdvb(1)/0.9-dvr(t)*0.228/1000];
C=[C,svr(1)==4.68/1000+0.9*gcvr(1)-gdvr(1)/0.9-dvr(1)*0.228/1000];

%末期蓄电量约束
C=[C,svb(24)==10.75/1000];
C=[C,svr(24)==5.68/1000];

%% 费用计算    
F2=sum(a*xconv+kcp*pmt+yconv*sconv)%燃气轮机费用
F1=sum(pmgb.*xb-pmgs.*xs);%购售电费用
%电动汽车电池损耗费用
F4=1000*sum(2*(gdvr/0.9+0.228/1000*dvr))+1000*sum(2*(gdvb/0.9+0.228/1000*dvb))
%需求响应负荷费用
F3=sum(pil(1,:)*kil(1))+sum(pil(2,:)*kil(2))+sum(pil(3,:)*kil(3));
F5=cc*sum(ep-eq);%碳市场交易成本
F6=sum(cshift)%可转移负荷补偿费用
F7=sum(Cgc);%绿证交易收入



F=F1+F2+F3+F4+F5+F6+F7 ; 
ops=sdpsettings('solver','cplex','verbose',2,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
result=optimize(C,F,ops);

if result.problem == 0 % problem =0 代表求解成功
else
    error('求解出错');
end

%% 结果数值读取
F1=value(F1);
F5=value(F5);
F7=value(F7);

pmgb=value(pmgb);
pmgs=value(pmgs);
gesc=value(gesc);
gesd=value(gesd);

pmt=value(pmt);
pil=value(pil);
gcvb=value(gcvb);
gdvb=value(gdvb);
gcvr=value(gcvr);
gdvr=value(gdvr);
svb=value(svb);
svr=value(svr);
ep=value(ep);
eq=value(eq);
eqp=eq-ep;
shiftl=value(shiftl);
shiftq=value(shiftq);

etotal=value(etotal);
Ere=value(Ere);
Qgc=value(Qgc);
Cgc=value(Cgc);

%% 结果展示
%1、各机组出力结果展示
figure
bar(gesc,'stacked');
hold on
bar(-gesd,'stacked')
hold on
plot(pload,'r-o','LineWidth',1.5);
hold on
plot(-ppv,'g-*','LineWidth',1.5);
hold on
plot(-pmt,'b-d','LineWidth',1.5);
hold on
plot(-wt,'y-o','LineWidth',1.5);
hold on
xlabel('时间/h');
ylabel('功率/MW');
title('聚合单元的基本调度结果');
legend('储能充电','储能放电','负荷','光伏','燃气轮机','风机出力');
legend('boxoff');
box off

% 可转移负荷转移特性
figure
plot(pload-shiftl+shiftq,'r-*','LineWidth',1.5)
hold on
plot(pload,'g--*','LineWidth',1.5)
xlabel('时间/h');
ylabel('负荷功率/MW');
title('可转移负荷响应结果');
legend('转移后的负荷曲线','转移前的负荷曲线');
legend('boxoff');
box off

% 转移负荷结果
figure
bar(shiftl-shiftq,'stacked')
xlabel('时间/h');
ylabel('转移的负荷功率/MW');
yyaxis right
plot(xs,'r-o','LineWidth',1.5);
hold on
ylabel('购电电价/（元/MWh）');
title('可转移负荷响应结果');
legend('负荷转移值','分时电价值');
legend('boxoff');
box off

%每一辆电动汽车调控结果-比亚迪
figure
bar(gcvb,'stacked');
hold on
bar(-gdvb,'stacked')
hold on
ylabel('功率/MW');
yyaxis right
plot(svb,'r-o','LineWidth',1.5);
xlabel('时间/h');
ylabel('蓄电量/MWh');
title('比亚迪电动汽车调度结果');
legend('电动汽车充电','电动汽车放电','蓄电量曲线');
legend('boxoff');
box off

%日产电动汽车调度结果
figure
bar(gcvr,'stacked');
hold on
bar(-gdvr,'stacked')
hold on
ylabel('功率/MW');
yyaxis right
plot(svr,'r-o','LineWidth',1.5);
xlabel('时间/h');
ylabel('蓄电量/MWh');
title('日产电动汽车调度结果');
legend('电动汽车充电','电动汽车放电','蓄电量曲线');
legend('boxoff');
box off

%中断负荷调度结果
figure
bar(pil',0.6,'stacked');
hold on
ylabel('功率/MW');
yyaxis right
plot(xb,'r-o','LineWidth',1.5);
hold on
ylabel('购电电价/（元/MWh）');
xlabel('时间/h');
title('可中断负荷调度结果');
legend('一级负荷','二级负荷','三级负荷','市场电价');
legend('boxoff');
box off

%储能分时电价调度结果
figure
bar(gesc,0.6,'stacked');
hold on
bar(-gesd,0.6,'stacked');
ylabel('功率/MW');
yyaxis right
plot(xb,'r-o','LineWidth',1.5);
hold on
ylabel('购电电价/（元/MWh）');
xlabel('时间/h');
title('分时电价下储能优化结果');
legend('储能充电','充能放电','电价');
legend('boxoff');
box off

figure
bar(ep/qg,0.6,'stacked');
hold on
plot(eq/nc,'r-o','LineWidth',1.5);
hold on
ylabel('额度');
xlabel('时间/h');
title('碳配额与碳排放曲线');
legend('碳排放曲线','碳配额曲线');
legend('boxoff');
box off

figure
bar(etotal,'c')
hold on
plot(Ere,'r-o','LineWidth',1.5)
hold on
plot(Qgc,'b--','LineWidth',1.5)
xlabel('时间')
ylabel('功率')
legend('总功率','风光实际消纳量','绿证交易量')
legend('boxoff');
box off

figure
Ff=-[F1;F5;F7];
bar(Ff,0.4,'g')
xlabel('项目')
ylabel('费用')
box off