目录
主要内容
程序主要做的是一个虚拟电厂或者微网单元的日前优化调度模型，考虑了光伏出力和负荷功率的双重不确定性，采用随机规划法处理不确定性变量，构建了虚拟电厂随机优化调度模型。具体来看，首先是基于蒙特卡洛算法，对预测的光伏以及负荷曲线进行场景生成，然后基于快概率距离快速消除法进行削减，直至削减至5个场景，然后采用随机调度的方法，对多场景下的虚拟电厂调度策略进行优化调度。
程序中燃气轮机、储能部分模型以及随机优化算法也是和下述文档一致。
1.1 场景生成及缩减
clc
clear
close all
%% 光伏场景生成以及削减
%生成负荷场景并削减%
%负荷出力预测均值E
Ww=[0 0 0 0 0 0 3.8 3.9 4.5 5.2 6.5 7.3 7.4 7.2 7.4 6.5 5.5 4.8 0 0 0 0 0 0];
% Ww=[0,0,0,0,0,1,2.5,4,5,5.5,5.8,5.7,5.5,5.3,5.1,5,3.8,2.5,1.2,0,0,0,0,0];
W=0.3*Ww
%取标准差为负荷出力预测值E的5%-20%,这里x=E*10%
l=W*0.1;
Ws=[];
%生成一个负荷场景，E+x*randn(1,24)，其中randn(1,24)为生成随机数的标准正态分布
m=200; %生成m个场景
for i=1:m
s=W+l.*randn(1,24);
Ws=[Ws;s];
end
figure(1)
[ss,gg]=meshgrid(1:200,1:24 );
plot3(ss,gg,Ws,'-');
grid
xlabel('场景');
ylabel('时刻');
zlabel('负荷功率');
title('负荷场景生成图');
% legend('负荷曲线1','负荷曲线2 ','负荷曲线3 ','负荷曲线4 ')
Ws_d=Ws; %定义削减后的场景
%场景削减
pi=1/m*ones(m,1); %蒙特卡罗生成的场景为等概率场景，建立每个场景的概率向量
%计算负荷场景Ws中每对场景的几何距离x
x=zeros(m,m);
for i=1:m
for j=1:m
x(i,j)=sum(abs(Ws(i,:)-Ws(j,:)));
end
end
%计算每个场景与剩余场景的概率距离之和y
y=zeros(m,1);
for i=1:m
y(i)=1/m*sum(x(i,:));
end
k=length(y);
%不断削减场景，直到剩余5个场景
while(k>5)
d=find(y==min(y)); %选定与剩余场景的概率距离之和最小的场景
x_2=x+100*eye(k); %构造新的x，以便找出风电场景Ws中与场景d几何距离最小的场景r
r=find(x_2(d,:)==min(x_2(d,:)));
pi(r)=pi(r)+pi(d); %将d场景的概率加到r场景上
%在负荷场景中删除d场景
pi(d)=[];
Ws_d(d,:)=[];
x(d,:)=[];
x(:,d)=[];
y(d)=[];
k=length(y);
end
figure(2)
[ss,gg]=meshgrid(1:5,1:24 );
plot3(ss,gg,Ws_d,'-');
grid
xlabel('场景');
ylabel('时刻');
zlabel('风机出力值');
title('场景削减图');
得到场景生成和缩减结果：
1.2 随机优化调度
模型中最核心的是随机优化模型的构建和处理，采用非预期运行约束处理随机优化问题，其中传统机组随机优化模型如下：
非预期约束是如何处理随机优化问题的呢？国内一些文献也是采用该类方法，如下所示，非预期约束也就是为了限定在风电、光伏以及其他因素存在的情况下，有且只有一种状态是被用在日前调度模型里的。
%主要内容：考虑光伏、负荷的不确定性，实现虚拟电厂的随机优化
%模型中，考虑了燃气轮机、储能，光伏，负荷四种单元，并考虑可以并网
clc
clear
close all
ppv=[0,0,0,0,0,0,1.34060681352041,1.13391651200680,1.40252116109929,1.74124103883023,1.49174914115029,2.53240793479920,2.04065972732231,2.35828237099756,1.55809477900388,1.64805990357317,1.87064383698464,1.29603798576622,0,0,0,0,0,0;0,0,0,0,0,0,1.18358196815285,1.20530370391118,1.35854372362592,1.32447876876720,1.70829563301966,2.51396807020581,1.87219733041650,1.51045430771326,2.14995685962274,2.36901954872155,1.49939386631594,1.24439925735197,0,0,0,0,0,0;0,0,0,0,0,0,1.01175119002517,1.31035724158009,1.08516736148552,1.94136887715302,1.93148023426153,2.08839059797124,2.16179304729596,1.66329638481645,2.26989483226563,2.44265872929522,1.36505525184415,1.44002086391399,0,0,0,0,0,0;0,0,0,0,0,0,1.04539328775374,1.03998113428866,1.18246929195709,1.69424435782838,1.64461569874716,2.06924213909286,2.43963602014925,1.68729987832366,1.87805529528869,1.51992087170733,2.01158827571321,1.29724518700574,0,0,0,0,0,0;0,0,0,0,0,0,1.32311860610225,1.28695117820334,1.34520597073109,1.46724609252811,1.74247560818882,1.72276716702525,2.57775664860191,1.79578137388444,2.66126270645276,2.15099361545381,1.66252664174142,1.63885764533067,0,0,0,0,0,0];
pload=[[0.431807609584421,0.486148921643357,0.721272534818094,0.898346267079082,0.944425724763714,0.541378372225899,0.882373258746104,1.56831193839750,1.86281649971011,1.99584713547466,2.19733700575767,2.29773449997527,2.47099246563420,1.66872324677400,1.77591203172356,1.32477944622446,1.95414778150622,1.54974896023576,1.74359932022130,1.44846182755801,0.955252219886710,1.25160480095078,1.49222676285844,0.682744656361740;0.483635538027688,0.495506878187239,0.714763198577902,0.799435561018225,1.01952146263562,0.800604394686264,0.945961336073563,1.28885742430358,1.90848488028237,2.10268787259155,1.91839575715085,1.85220123265328,2.17764591245376,1.78073710542207,1.73239834633622,1.42351956423826,2.03383140201547,2.58777879153454,1.89810311996295,1.45544594829551,1.11733457030710,0.835193945440698,1.45888556170521,0.737906988149361;0.475566682884608,0.538938919866738,0.694384127185758,0.948064145000996,0.869283322295024,0.693421312226079,1.24531781162386,1.18155825364294,2.24893779918531,2.22929790054771,1.91247902153441,1.60109975699407,2.04025369591884,1.62927362526231,1.94798019608120,1.30811995960382,1.92258374438264,2.34731329057125,1.88682790373556,1.82091112812632,0.761274726486709,1.14465406229059,1.66887813069703,0.687295849991400;0.506197427493947,0.555237638826251,0.711932167351328,0.723600759144382,1.04613991322447,0.658159592550928,1.07636853878386,1.06165890483766,2.06902267439968,2.48488721178039,2.34092228143214,2.63450984929863,2.35352014698963,2.08958282274750,1.65408992214257,1.54410971082492,1.94697999777118,1.92406941886044,2.00023938323277,2.10124494602628,1.04036325929128,0.984006134878927,1.94126189741016,0.651136186886876;0.504831465058796,0.522880453611977,0.805416372257449,0.776866377863578,1.09776915278538,0.499517330314076,0.895214516450348,1.27340524098249,1.51733245769139,1.46096552221328,2.35562457523011,1.98343189442592,1.73365278725908,1.98502985411742,2.00458112949663,1.35148256752316,1.86202305372752,2.23007113968369,1.56324879332238,1.57649124469752,0.995615716738264,1.11287433329993,1.96306488103865,0.572126166745211]];
%% 定义变量
%%定义市场购电电价以及售电电价
xb=[630,630,630,630,630,630,1020,1020,1020,1520,1520,1520,1520,1520,1020,1020,1020,1520,1520,1520,1020,1020,630,630];
xs1=[100,100,100,100,100,100,380,380,380,800,800,800,800,800,380,380,380,800,800,800,380,380,100,100];
xs=1.05*xs1;
%% 定义燃气轮机参数
a=600;%固定开机费用
kcp=100;%分段线性化费用
sconv=100;%启停费用
gtmax=3.31;%出力上限
gtmin=1.3;%最小出力值
ramp=1.5;%爬坡率
%% 定义储能参数
gescmax=1;%充电功率上限
gesdmax=1;%放电功率上限
sessmax=4;%蓄电量上限
sessmin=0;%蓄电量最小值
uesc=0.95;%充电效率
uesd=0.95;%放电效率
kil=[500,700,800];%中断负荷补偿费用
%% 其他输入参数
pmgmax=20;%最大交易量
%负荷值
% pload=[1.5 1.8 2.3 2.8 3.2 2.1 3.3 4.2 5.9 6.8 7.5 7.2 7.1 6.5 5.9 4.8 5.6 6.8 6.8 6.2 3.3 3.6 5.4 2.4];
% %光伏出力
% ppv=[0 0 0 0 0 0 3.8 3.9 4.5 5.2 6.5 7.3 7.4 7.2 7.4 6.5 5.5 4.8 0 0 0 0 0 0];
%% 定义变量sdpvar/binvar
umob=binvar(1,24);%是否购电
umos=binvar(1,24);%是否售电
umospf=binvar(5,5,24)
umobpf=binvar(5,5,24)
pmgb=sdpvar(1,24);%市场购电量
pmgs=sdpvar(1,24);%市场售电量
pmgbpf=sdpvar(5,5,24);
pmgspf=sdpvar(5,5,24);
xconv=binvar(1,24);%燃气轮机工作状态变量
yconv=binvar(1,24);%燃气轮机启停状态变量
xconvpf=binvar(5,5,24);
yconvpf=binvar(5,5,24);
pmt=sdpvar(1,24);%燃气轮机出力
pmtpf=sdpvar(5,5,24);
gesc=sdpvar(1,24);%储能充电功率
gesd=sdpvar(1,24);%储能放电功率
sess=sdpvar(1,24);%蓄电池蓄电量
gescpf=sdpvar(5,5,24);
gesdpf=sdpvar(5,5,24);
sesspf=sdpvar(5,5,24);
%% 约束条件
C=[];%初始化约束
%% 燃气轮机出力约束
for p=1:5
for f=1:5
for t=1:24
C=[C,
xconvpf(p,f,t)*gtmin<=pmtpf(p,f,t)<=xconvpf(p,f,t)*gtmax ,%出力上下限约束
];
end
end
end
%
for p=1:5
for f=1:5
C=[C,
pmtpf(p,f,1)<=ramp,%初始爬坡约束
xconvpf(p,f,1)<=yconvpf(p,f,1),%初始启停约束
];
end
end
%
for p=1:5
for f=1:5
for t=2:24
C=[C,
-ramp<=pmtpf(p,f,t)-pmtpf(p,f,t-1)<=ramp,%爬坡率约束
xconvpf(p,f,t)-xconvpf(p,f,t-1)<=yconvpf(p,f,t), %工作状态约束
];
end
end
end