目录

1主要内容

1.1 场景生成及缩减

1.2 随机优化调度

2 程序链接

1主要内容

程序主要做的是一个虚拟电厂或者微网单元的日前优化调度模型，考虑了光伏出力和负荷功率的双重不确定性，采用随机规划法处理不确定性变量，构建了虚拟电厂随机优化调度模型。具体来看，首先是基于蒙特卡洛算法，对预测的光伏以及负荷曲线进行场景生成，然后基于快概率距离快速消除法进行削减，直至削减至5个场景，然后采用随机调度的方法，对多场景下的虚拟电厂调度策略进行优化调度。

程序中燃气轮机、储能部分模型以及随机优化算法也是和下述文档一致。

1.1 场景生成及缩减

clcclearclose all​%% 光伏场景生成以及削减%生成负荷场景并削减%%负荷出力预测均值EWw=[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:ms=W+l.*randn(1,24);Ws=[Ws;s];end​figure(1)[ss,gg]=meshgrid(1:200,1:24 );plot3(ss,gg,Ws,'-');gridxlabel('场景');ylabel('时刻');zlabel('负荷功率');title('负荷场景生成图');% legend('负荷曲线1','负荷曲线2 ','负荷曲线3 ','负荷曲线4 ')​​Ws_d=Ws; %定义削减后的场景%场景削减pi=1/m*ones(m,1); %蒙特卡罗生成的场景为等概率场景，建立每个场景的概率向量%计算负荷场景Ws中每对场景的几何距离xx=zeros(m,m); for i=1:m    for j=1:m            x(i,j)=sum(abs(Ws(i,:)-Ws(j,:)));    endend​%计算每个场景与剩余场景的概率距离之和yy=zeros(m,1);for i=1:my(i)=1/m*sum(x(i,:));endk=length(y);​%不断削减场景，直到剩余5个场景while(k>5)d=find(y==min(y)); %选定与剩余场景的概率距离之和最小的场景x_2=x+100*eye(k); %构造新的x，以便找出风电场景Ws中与场景d几何距离最小的场景rr=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,'-');gridxlabel('场景');ylabel('时刻');zlabel('风机出力值');title('场景削减图');

得到场景生成和缩减结果：

1.2 随机优化调度

模型中最核心的是随机优化模型的构建和处理，采用非预期运行约束处理随机优化问题，其中传统机组随机优化模型如下：

非预期约束是如何处理随机优化问题的呢？国内一些文献也是采用该类方法，如下所示，非预期约束也就是为了限定在风电、光伏以及其他因素存在的情况下，有且只有一种状态是被用在日前调度模型里的。

摘自《交直流混合配电网多阶段随机优化调度模型_裴蕾》

%主要内容：考虑光伏、负荷的不确定性，实现虚拟电厂的随机优化%模型中，考虑了燃气轮机、储能，光伏，负荷四种单元，并考虑可以并网clcclearclose allppv=[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/binvarumob=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    endend% for p=1:5    for f=1:5       C=[C,           pmtpf(p,f,1)<=ramp,%初始爬坡约束           xconvpf(p,f,1)<=yconvpf(p,f,1),%初始启停约束       ];     endend​% 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    endend

程序结果：