电动汽车有序充放电

clc; clear; %% 抽样 load('Parameters.mat'); interval=zeros(1,96); EVNoCount=0; %可提供充电的总数目 EVstarttime=zeros(1,SampleNo); EVssoc=zeros(1,SampleNo); EVendtime=zeros(1,SampleNo); startt=12.25; for t=1:96 %96个时段 从中午12点开始计算，也就是从第50个时段开始计算, 12.25为第一个时间段 starttime=startt+0.25*(t-1); %%对电动汽车的到达时间及车辆数进行进行抽样 if starttime>24 starttime=starttime-24; end o=3.41;u=17.47; if starttime<=(u-12) liangshu=@(x)(1./(o.*((2.*pi).^0.5))).*exp(-((x+24-u).^2)./(2.*(o.^2)));%每个刷新时段内的入网量数的概率密度函数 else liangshu=@(x)(1./(o.*((2.*pi).^0.5))).*exp(-((x-u).^2)./(2.*(o.^2)));%每个刷新时段内的入网量数的概率密度函数 end s1=quadl(liangshu,starttime-0.25,starttime);%0.25代表刷新时段为0.25小时 intervalNo=round(s1*SampleNo); %该段时间的接入的车辆数目 interval(t)=intervalNo; %%对初始SOC进行抽样 reasonnumber=0; for i=1:intervalNo MU=2.98; SIGMA=1.14; R = lognrnd(MU,SIGMA); socs= 0.9- 15*R/(100*Ci); while socs<0.1 || socs>=0.3 R = lognrnd(MU,SIGMA); socs= 0.9- 15*R/(100*Ci); end socs=fix(socs*100)/100; endtimeone=normrnd(8.92,3.24); while endtimeone>=12 || endtimeone<=2 endtimeone=normrnd(8.92,3.24); end endtimeone=ceil((endtimeone*60)/15)+48; if (starttime>12)||(starttime<endtimeone) maxJ=endtimeone-t; if (0.25*maxJ*Pi>((0.9-socs)*Ci)) reasonnumber=reasonnumber+1; EVstarttime(1,EVNoCount+reasonnumber)=t; EVssoc(1,EVNoCount+reasonnumber)=socs; EVendtime(1,EVNoCount+reasonnumber)=endtimeone; end end end EVNoCount=EVNoCount+reasonnumber; end RealEVstarttime=zeros(1,EVNoCount); RealEVssoc=zeros(1,EVNoCount); RealEVendtime=zeros(1,EVNoCount); for i=1:EVNoCount RealEVstarttime(1,i)= EVstarttime(1,i); RealEVssoc(1,i)=EVssoc(1,i); RealEVendtime(1,i)=EVendtime(1,i); end clear EVstarttime EVssoc EVendtime; %%全局优化 dt=0.25;%时间间隔 Snj=binvar(EVNoCount,96); %控制变量，0-1型 Transformer=ones(1,96); Aj=ones(96,1); re=0.9*ones(EVNoCount,1); Cend=ones(EVNoCount,1); Price=[ 0.687 0.687 0.687 0.687 0.687 0.869 0.869 0.869 0.869 0.687 0.687 0.687 0.365 0.365 0.365 0.365 0.365 0.365 0.365 0.365 0.869 0.869 0.869 0.869]; pricekuo=zeros(1,96); for i=1:95 pricekuo(i)=Price(floor(i/4)+1); if(floor((i-1)/4)<floor(i/4)) pricekuo(i)=Price(floor(i/4)); end end pricekuo(96)=Price(24); %f=Pni*dt*sum(Snj*pricekuo'); f=max(fh+Pni*sum(Snj))-min(fh+Pni*sum(Snj)); F=[]; F=set((Pni*sum(Snj))<= (Brl*Transformer-fh));%构建约束条件 F=F+set((Pi*dt*Snj*Aj+Ci*(RealEVssoc'))>=(Ci*re))+set((Pi*dt*Snj*Aj+Ci*(RealEVssoc'))<=(Ci*Cend));%多个约束条件 for i = 1:EVNoCount for k=1:96 if (k<RealEVstarttime(i))|| (k> RealEVendtime(i)) F = [F,Snj(i,k)<1]; end end end options = sdpsettings('solver','cplex','verbose',1);%'solver' 参数指定程序用cplex求解器（如果已经安装，否则会报错），如果不指定 ‘solver’ 参数，他会根据决策变量类型自动挑选已安装的、最适合的求解器；'verbose' 指定显示冗余度（冗余度越大，你就可以看到越详细的求解过程信息）。 %求解 tic; sol = solvesdp(F,f,options);%函数solvesdp用来求解优化问题，这里解得是F条件下目标函数f的最小值，要求最大值f前面加个负号；s=solvesdp(F);%求解可行解问题；s=solvesdp(F, f);%求解一般优化问题，其中f为目标函数；s=solvesdp(F, f, options); %设定选项，比如选择算法等 if sol.problem == 0 solution = double(Snj);%提取解矩阵 else solution = double(Snj); display('Hmm, something went wrong!'); sol.info yalmiperror(sol.problem) end toc zongfh300=fh+Pni*sum(solution);%日常负荷+充电负荷 fh300=fh;%日常负荷 display(fh300); %保存参数数据 save('wx300','fh300','zongfh300' );