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