phpstudy集成环境

 首页 > 这程序咋运行

这程序咋运行

function [ A b_L b_U ] = juzhen
%JUZHEN Summary of this function goes here
% Detailed explanation goes here

% clear
% clc
%光伏的概率
% ezplot(@(x) exp(1x)-1/0.864xexp(1x),[0,0.864]);
a=0.11-0.01
b=1.5+0.5
eta=a;
Q=b;
% Q=2.5;
Q=Q/2;
n=310;
a=0.2652;
Ac=700;
ita=0.087;
ku=0.864;
xishu=0.5;
w=0;
C=0.6251;
lamda=3.5507;
% lamda=4;
phi=0.409sin(2pi/365*n-1.39);
I0=1.1(sin(a)sin(phi)+cos(a)cos(phi)cos(w));
T=zeros(1,18);
Tdot=zeros(1,18);
TT=zeros(1,18);
TTdot=zeros(1,18);
%根据时角计算输出
T(9)=0.611;
Tdot(9)=0.016;
T(15)=T(9);
Tdot(15)=Tdot(9);

T(12)=(sin(a)sin(phi)+cos(a)cos(phi)cos(0))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*T(9);
Tdot(12)=(sin(a)sin(phi)+cos(a)cos(phi)cos(0))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*Tdot(9);

T(11)=(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/12))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*T(9);
Tdot(11)=(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/12))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*Tdot(9);
T(13)=T(11);
Tdot(13)=Tdot(11);

T(14)=(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/6))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*T(9);
Tdot(14)=(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/6))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*Tdot(9);
T(10)=T(14);
Tdot(10)=Tdot(14);

T(16)=(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/3))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*T(9);
Tdot(16)=(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/3))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))*Tdot(9);
T(8)=T(16);
Tdot(8)=Tdot(16);

T(17)=(sin(a)sin(phi)+cos(a)cos(phi)cos(5pi/12))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))T(9);
Tdot(17)=(sin(a)sin(phi)+cos(a)cos(phi)cos(5pi/12))/(sin(a)sin(phi)+cos(a)cos(phi)cos(pi/4))Tdot(9);
T(7)=T(17);
Tdot(7)=Tdot(17);
for i=8:18
PVmax(i)=Ac0.093(T(i-1)0.864-Tdot(i-1)0.746);
alfa(i)=T(i-1)/Tdot(i-1);
end
for i=8:18
TT(i)=T(i-1);
TTdot(i)=Tdot(i-1);
end
T=TT;
Tdot=TTdot;
% C=0.29;
% lamda=1.21;
% C=0.6251;
% lamda=3.5507;
% Ac=5;
% ezplot(@(pv) C(ku-0.5(alfa(9)-sqrt(alfa(9)^2-4pv/(0.093Tdot(9)Ac))))exp(0.5lamda(alfa(9)-sqrt(alfa(9)^2-4pv/(0.093Tdot(9)Ac))))/(kuAc0.093T(9)),[0,0.24]);
% prb=int('C(ku-xishu(alfa(9)-sqrt(alfa(9)^2-4pv/(itaTdot(9)Ac))))exp(xishulamda(alfa(9)-sqrt(alfa(9)^2-4pv/(itaTdot(9)Ac))))/(kuAcitaT(9))','pv','0','x');%符号积分
% ezplot(@(pv) C(ku-0.5(alfa(7)-sqrt(alfa(7)^2-4pv/(0.093Tdot(7)650))))exp(0.5lamda(alfa(7)-sqrt(alfa(7)^2-4pv/(0.093Tdot(7)650))))/(kuAc0.093T(7)),[0,8.4052]);%画图
%------------20151121--------------
%prb=int('C(ku-xishu(alfa(7)-sqrt(alfa(7)^2-4pv/(itaTdot(7)Ac))))exp(xishulamda(alfa(7)-sqrt(alfa(7)^2-4pv/(itaTdot(7)Ac))))/(kuAcitaT(7))','pv','0','x');%符号积分
pr=@(x) -1/2C(-exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)kulamda^2xishu(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)itaTdot(7)Ac-exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)kulamdaitaTdot(7)Ac-exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)xishu^2lamda^2alfa(7)^2itaTdot(7)Ac+4exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)xishu^2lamda^2x+exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)xishu^2lamda^2(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)alfa(7)itaTdot(7)Ac-2exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)itaTdot(7)Ac+exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)xishulamdaalfa(7)itaTdot(7)Ac-2exp(-xishulamda(-(-alfa(7)^2itaTdot(7)Ac+4x)/ita/Tdot(7)/Ac)^(1/2)+xishulamdaalfa(7))exp((alfa(7)^2)^(1/2)xishulamda)itaTdot(7)Ac+exp(xishulamdaalfa(7))kulamda^2xishu(alfa(7)^2)^(1/2)itaTdot(7)Ac-exp(xishulamdaalfa(7))xishu^2lamda^2(alfa(7)^2)^(1/2)alfa(7)itaTdot(7)Ac+exp(xishulamdaalfa(7))kulamdaitaTdot(7)Ac+exp(xishulamdaalfa(7))lamda^2xishu^2alfa(7)^2itaTdot(7)Ac-exp(xishulamdaalfa(7))xishualfa(7)lamdaitaTdot(7)Ac+2exp(xishulamdaalfa(7))itaTdot(7)Ac+2exp(xishulamdaalfa(7))(alfa(7)^2)^(1/2)xishulamdaitaTdot(7)Ac)/ita/Ac/exp((alfa(7)^2)^(1/2)xishulamda)/lamda^3/ku/T(7)/xishu^2;
%光伏离散化
PVXL=zeros(24,50);
lamdaa=zeros(1,24);
caa=zeros(1,24);
lamdaa=[3.5507 3.5507 3.5507 3.5507 3.5507 3.5507 3.5507 3.5707 3.5707 4.5 5.1 5.5 7 5.5 5.1 4.5 3.5507 3.5707 3.7507 3.5507 3.5507 3.5507 3.5507 3.5507];%改过，让概率逼近1，,35707，0.6351
caa=[0.6251 0.6251 0.6251 0.6251 0.6251 0.6251 0.6251 0.6351 0.6351 0.41 0.3 0.245 0.105 0.245 0.3 0.41 0.6351 0.6351 0.6251 0.6251 0.6251 0.6251 0.6251 0.6251];

                                                     %7点

Npv=[1 1 1 1 1 1 1 3 8 7 9 10 10 10 9 7 8 3 1 1 1 1 1 1];
for i=8:18
Npv(i)=fix(PVmax(i)/Q)+1; %离散步长设置处
end
for i=8:18%从0：00开始的，7点才有功率
C=caa(i);
lamda=lamdaa(i);
N=Npv(i)-1;
%------------20151121--------------
% prb=int('C(ku-xishu(alfa(i)-sqrt(alfa(i)^2-4pv/(itaTdot(i)Ac))))exp(xishulamda(alfa(i)-sqrt(alfa(i)^2-4pv/(itaTdot(i)Ac))))/(kuAcitaT(i))','pv','0','x');%符号积分
pr=@(x) -1/2C(-exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)kulamda^2xishu(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)itaTdot(i)Ac-exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)kulamdaitaTdot(i)Ac-exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)xishu^2lamda^2alfa(i)^2itaTdot(i)Ac+4exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)xishu^2lamda^2x+exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)xishu^2lamda^2(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)alfa(i)itaTdot(i)Ac-2exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)itaTdot(i)Ac+exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)xishulamdaalfa(i)itaTdot(i)Ac-2exp(-xishulamda(-(-alfa(i)^2itaTdot(i)Ac+4x)/ita/Tdot(i)/Ac)^(1/2)+xishulamdaalfa(i))exp((alfa(i)^2)^(1/2)xishulamda)itaTdot(i)Ac+exp(xishulamdaalfa(i))kulamda^2xishu(alfa(i)^2)^(1/2)itaTdot(i)Ac-exp(xishulamdaalfa(i))lamda^2xishu^2(alfa(i)^2)^(1/2)alfa(i)itaTdot(i)Ac+exp(xishulamdaalfa(i))kulamdaitaTdot(i)Ac+exp(xishulamdaalfa(i))lamda^2xishu^2alfa(i)^2itaTdot(i)Ac-exp(xishulamdaalfa(i))xishualfa(i)lamdaitaTdot(i)Ac+2exp(xishulamdaalfa(i))itaTdot(i)Ac+2exp(xishulamdaalfa(i))(alfa(i)^2)^(1/2)xishulamdaitaTdot(i)Ac)/ita/Ac/exp((alfa(i)^2)^(1/2)xishulamda)/lamda^3/ku/T(i)/xishu^2;
PVXL(i,1)=pr(PVmax(i)0.51/N);
for j=2:N
PVXL(i,j)=pr(PVmax(i)(0.51/N+(j-1)/N))-pr(PVmax(i)(0.51/N+(j-2)/N));
end
PVXL(i,N+1)=pr(PVmax(i))-pr(PVmax(i)(0.51/N+(N-1)/N));
end
%风电概率
% v_in = 5; %测试点
% v_r = 17;
v_in = 4;
v_r = 15;
v_out = 26;
WXL=[];
ca= zeros(0); %平均风速 1x24
k0a= zeros(0); %形状系数 1x24
Na=zeros(0); %离散份数 1x24
Vavr=[13.5 13.35 12.98 10.65 10.38 9.43 10.13 10.12 10.27 10.73 12.36 11.73 10.23 9.43 8.1 7.45 6.38 7.35 8.15 10.15 14.03 14.98 14.58 14.54];
% k0a=[2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ];
k0a=[2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.01 1.98 1.97 2.15 1.89 1.89 1.89 1.92 1.94 1.98 1.99 2.0 2.0 2.4 2.0 2.0 2.0 ];
for i=1:24
ca(i)=Vavr(i)/gamma(1+1/k0a(i));
end

for i=1:24
Nw(i)=fix(30/Q)+1;
end
for i=1:24
c=ca(i);
k0=k0a(i);
N=Nw(i)-1;%项数换成长度
Pro0_p =@(pwind) 1-exp(-((v_in+(v_r-v_in)*pwind)/c).^k0)+exp(-(v_out/c).^k0);%风电概率分布函数 pwind标幺值
WXL(i,1)=Pro0_p(0.5*1/N);
for j=2:N
WXL(i,j)=Pro0_p(0.51/N+(j-1)/N)-Pro0_p(0.51/N+(j-2)/N);
end
WXL(i,N+1)=Pro0_p(0.99999)-Pro0_p(0.5*1/N+(N-1)/N)+exp(-(v_r/c)^k0)-exp(-(v_out/c)^k0);
end
%离散卷积
probility=sparse(zeros(24,300));
for i=8:18%第8个开始

for p=1:Nw(i)
    for q=1:Npv(i)
        probility(i,p+q)=probility(i,p+q)+WXL(i,p)*PVXL(i,q);
    end
end

end
%##########################
for i=8:18 %如何处理概率的问题初始元素位置麻烦
for j=1:Nw(i)+Npv(i)-1

   Probability(i,j)= probility(i,j+1);

end
end
%##########################
% delta=[Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q];
% Pg=zeros(1,24);
% for i=8:18
% for j=2:Nw(i)+Npv(i)
% Pinterval(i,j)=delta(i)(j-2)probility(i,j);
% Pg(i)=Pg(i)+Pinterval(i,j);
% end
% end
% for i=1:7
% for j=2:Nw(i)
% Pinterval(i,j)=delta(i)(j-1)WXL(i,j);
% Pg(i)=Pg(i)+Pinterval(i,j);
% end
% end
% for i=19:24
% for j=2:Nw(i)
% Pinterval(i,j)=delta(i)(j-1)WXL(i,j);
% Pg(i)=Pg(i)+Pinterval(i,j);
% end
% end
%##########################
for i=1:7
for j=1:Nw(i)

   Probability(i,j)= WXL(i,j);

end
end
for i=19:24
for j=1:Nw(i)

   Probability(i,j)= WXL(i,j);

end
end
%##########################
avr=zeros(24,1);
for i=1:24

 for j=1:Nw(i)+Npv(i)-1
  avr(i)=avr(i)+ Probability(i,j);
 end

end
for i=1:24

 for j=1:Nw(i)+Npv(i)-1
     Probability(i,j)=Probability(i,j)/avr(i);
 end

end
delta=[Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q];
Pg=zeros(1,24);
for i=8:18
for j=2:Nw(i)+Npv(i)-1

Pinterval(i,j)=delta(i)*(j-1)*Probability(i,j);
Pg(i)=Pg(i)+Pinterval(i,j);

end
end
for i=1:7

for j=2:Nw(i)
  Pinterval(i,j)=delta(i)*(j-1)*Probability(i,j);  
  Pg(i)=Pg(i)+Pinterval(i,j);

end
end
for i=19:24

for j=2:Nw(i)
  Pinterval(i,j)=delta(i)*(j-1)*Probability(i,j);  
  Pg(i)=Pg(i)+Pinterval(i,j);

end
end

Pg=Pg*2;%平均出力太小所以乘2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%UC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Name='UNIT COMMITMNET';
% eta=0.08;
% eta=1;
L=400;
IntVars=[];
c=[];
x_L=[];
x_U=[];
b_U=[];
b_L=[];
%x=[r1 r2 r3 p1 p2 p3 i1 i2 i3 pin pout su1 su2 su3 sd1 sd2 sd3 soc Cres uload];201在2024个变量后，继续每个时段至少加11个变量
c1=[0.04 0.04 0.04 0.35 0.35 0.26 1.2 1.2 1 0 0.0001 1.6 1.6 3.5 0 0 0 0 0.0001 0]';%1*20
x_L1=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 162 0 0]';%201
x_U1=[25 25 55 30 30 65 1 1 1 220/0.9 220/0.9 1 1 1 1 1 1 280 220/0.9 200]';%20*1
PD=[30 31 32 35 32 50 75 100 120 130 150 150 145 135 135 120 140 110 105 80 60 50 30 30];
%生成A阵
A1=[0 0 0 1 1 1 0 0 0 -1 0.9 0 0 0 0 0 0 0 0 -1;
0 0 0 0 0 0 0 0 0 0.9 -1 0 0 0 0 0 0 -1 0 0; %17*20 %行数加一

0 0 0 1 0 0 -30 0 0  0 0 0 0 0 0 0 0 0 0 0;

0 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 -30 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 -65 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 -1 0 0 10 0 0 0 0 0 0 0 0 0 0 0;
1 0 0 1 0 0 -30 0 0 0 0 0 0 0 0 0 0 0 0 0;
-1 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 1 0 0 -30 0 0 0 0 0 0 0 0 0 0 0 0;
0 -1 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 1 0 0 -65 0 0 0 0 0 0 0 0 0 0 0;
0 0 -1 0 0 -1 0 0 10 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 0 0;
];

A2=[0 0 0 1 1 1 0 0 0 -1 0.9 0 0 0 0 0 0 0 0 -1;%20*20

0 0 0 0 0 0 0 0 0 0.9 -1 0 0 0 0 0 0 -1 0 0;
0 0 0 1 0 0 -30 0 0  0 0 0 0 0 0 0 0 0 0 0;

0 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 -30 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 -65 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 -1 0 0 10 0 0 0 0 0 0 0 0 0 0 0;
1 0 0 1 0 0 -30 0 0 0 0 0 0 0 0 0 0 0 0 0;
-1 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 1 0 0 -30 0 0 0 0 0 0 0 0 0 0 0 0;%x=[r1 r2 r3 p1 p2 p3 i1 i2 i3 pin pout su1 su2 su3 sd1 sd2 sd3 soc];1*18
0 -1 0 0 -1 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 1 0 0 -65 0 0 0 0 0 0 0 0 0 0 0;
0 0 -1 0 0 -1 0 0 10 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0;
0 0 0 0 0 0 -1 0 0 0 0 1 0 0 -1 0 0 0 0 0;
0 0 0 0 0 0 0 -1 0 0 0 0 1 0 0 -1 0 0 0 0; %停开机0-1变量
0 0 0 0 0 0 0 0 -1 0 0 0 0 1 0 0 -1 0 0 0; %i(t-1)没考虑
];
B1=[A1 zeros(17,20*23)];
B2=[zeros(20,20) A2 zeros(20,20*22)];
B3=[zeros(20,202) A2 zeros(20,2021)];
B4=[zeros(20,203) A2 zeros(20,2020)];
B5=[zeros(20,204) A2 zeros(20,2019)];
B6=[zeros(20,205) A2 zeros(20,2018)];
B7=[zeros(20,206) A2 zeros(20,2017)];
B8=[zeros(20,207) A2 zeros(20,2016)];
B9=[zeros(20,208) A2 zeros(20,2015)];
B10=[zeros(20,209) A2 zeros(20,2014)];
B11=[zeros(20,2010) A2 zeros(20,2013)];
B12=[zeros(20,2011) A2 zeros(20,2012)];
B13=[zeros(20,2012) A2 zeros(20,2011)];
B14=[zeros(20,2013) A2 zeros(20,2010)];
B15=[zeros(20,2014) A2 zeros(20,209)];
B16=[zeros(20,2015) A2 zeros(20,208)];
B17=[zeros(20,2016) A2 zeros(20,207)];
B18=[zeros(20,2017) A2 zeros(20,206)];
B19=[zeros(20,2018) A2 zeros(20,205)];
B20=[zeros(20,2019) A2 zeros(20,204)];
B21=[zeros(20,2020) A2 zeros(20,203)];
B22=[zeros(20,2021) A2 zeros(20,202)];
B23=[zeros(20,2022) A2 zeros(20,201)];
B24=[zeros(20,20*23) A2];
A=[B1;B2;B3;B4;B5;B6;B7;B8;B9;B10;B11;B12;B13;B14;B15;B16;B17;B18;B19;B20;B21;B22;B23;B24];
[m,n] = size(A);

for i=1:23 %Soc(t-1)+Pin(t)-Pout(t)=Soc(t);

A(17+20*(i-1)+2,20*i-2)=1;

end
for k=1:23

i=17+(k-1)*20;%su-sd-i(t)+i(t-1)=0;
j=(k-1)*20;
A(i+18,j+7)=1;
A(i+19,j+8)=1;
 A(i+20,j+9)=1;

end
i=0;
j=0;
for i=1:24

j=j+Npv(i)+Nw(i)-1;

end
pronum=j;
% S=zeros(823,730);
S=[A,zeros(477,j)];
p=0;
q=0;%新加的元素坐标
k=[];
for p=1:24

k(p)=0;
for q=1:p
    k(p)=k(p)+Npv(q)+Nw(q)-1;%减过1
end

end
p=1;
% V=[];

for q=1:pronum
        if(q>k(p))
          p=p+1;  
        end
    S(477+q,1+(p-1)*20)=-1/L;%备用设置较为复杂,储能还未考虑
    S(477+q,2+(p-1)*20)=-1/L;
    S(477+q,3+(p-1)*20)=-1/L;
    S(477+q,20*p-1)=-0.9*1/L;
    S(477+q,24*20+q)=1;      
end

S=sparse(S);

p=0;
 for i=1:24
     for j=1:Npv(i)+Nw(i)-1
         if i~=1
     S(477+pronum+i,24*20+j+k(i-1))=Probability(i,j);%注意这里probility未改动，后期需要将风电概率加进来
         else
         S(477+pronum+i,24*20+j)=Probability(i,j);
         end
     end
 end

for i=1:24 %储能备用容量限制 0.9Cres+0.9Pout<=20 Cres变量定义在最后
S(477+pronum+24+i,20*i-1)=0.9;
S(477+pronum+24+i,11+(i-1)*20)=0.9;

S(477+pronum+24+24+i,20i-1)=0.9;% 0.9Cres<=0.9*（Soc(t)-Socmin)
S(477+pronum+24+24+i,20*i-2)=-0.9;
end
S=sparse(S);

x_U2=[25 25 55 30 30 65 1 1 1 220/0.9 220/0.9 1 1 1 1 1 1 216 220/0.9 300]';
for i=1:24 %传统24小时组合

c=[c;c1];
x_L=[x_L;x_L1];
 if i==24
x_U=[x_U;x_U2];
else
x_U=[x_U;x_U1];
end

end
for i=1:pronum %概率0-1变量

 c=[c;0];
 x_L=[x_L;0];
 x_U=[x_U;1];

end
% for i=1:24 %为了备用限制
% c=[c;0];
% x_L=[x_L;0];
% x_U=[x_U;20];
% end

b_U1=[PD(1)-Pg(1);-32;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;0; 0; 0 ];%1是功率平衡，2是储能平衡，3~8为出力限制，9~14为备用限制；15~17为开机限制
b_L1=[PD(1)-Pg(1);-32;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ];%SOC初始值为32.
b_U2=[PD(2)-Pg(2);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];%1是功率平衡，2是储能平衡，3~8为出力限制，9~14为备用限制；15~20为开关机限制
b_L2=[PD(2)-Pg(2);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U3=[PD(3)-Pg(3);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L3=[PD(3)-Pg(3);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U4=[PD(4)-Pg(4);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L4=[PD(4)-Pg(4);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U5=[PD(5)-Pg(5);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L5=[PD(5)-Pg(5);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U6=[PD(6)-Pg(6);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L6=[PD(6)-Pg(6);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U7=[PD(7)-Pg(7);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L7=[PD(7)-Pg(7);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U8=[PD(8)-Pg(8);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L8=[PD(8)-Pg(8);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U9=[PD(9)-Pg(9);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L9=[PD(9)-Pg(9);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U10=[PD(10)-Pg(10);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L10=[PD(10)-Pg(10);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U11=[PD(11)-Pg(11);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L11=[PD(11)-Pg(11);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U12=[PD(12)-Pg(12);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L12=[PD(12)-Pg(12);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U13=[PD(13)-Pg(13);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L13=[PD(13)-Pg(13);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U14=[PD(14)-Pg(14);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L14=[PD(14)-Pg(14);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U15=[PD(15)-Pg(15);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L15=[PD(15)-Pg(15);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U16=[PD(16)-Pg(16);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L16=[PD(16)-Pg(16);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U17=[PD(17)-Pg(17);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L17=[PD(17)-Pg(17);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U18=[PD(18)-Pg(18);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L18=[PD(18)-Pg(18);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U19=[PD(19)-Pg(19);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L19=[PD(19)-Pg(19);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U20=[PD(20)-Pg(20);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L20=[PD(20)-Pg(20);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U21=[PD(21)-Pg(21);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L21=[PD(21)-Pg(21);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U22=[PD(22)-Pg(22);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L22=[PD(22)-Pg(22);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U23=[PD(23)-Pg(23);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L23=[PD(23)-Pg(23);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];
b_U24=[PD(24)-Pg(24);0;0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0 ;1; 1; 1 ;0;0;0];
b_L24=[PD(24)-Pg(24);0;-200; -200; -200 ;-200; -200; -200; -200 ;-200 ;-200; -200; -200;-200;0; 0; 0 ;0;0;0];

b_U=[b_U1;b_U2;b_U3;b_U4;b_U5;b_U6;b_U7;b_U8;b_U9;b_U10;b_U11;b_U12;b_U13;b_U14;b_U15;b_U16;b_U17;b_U18;b_U19;b_U20;b_U21;b_U22;b_U23;b_U24];
b_L=[b_L1;b_L2;b_L3;b_L4;b_L5;b_L6;b_L7;b_L8;b_L9;b_L10;b_L11;b_L12;b_L13;b_L14;b_L15;b_L16;b_L17;b_L18;b_L19;b_L20;b_L21;b_L22;b_L23;b_L24];

p=1;
for i=1:24

for j=1:Nw(i)+Npv(i)-1

b_L=[b_L;(-Pg(i)+(j-1)2delta(i))/L];
b_U=[b_U;1+(-Pg(i)+(j-1)2delta(i))/L]; %需要考虑储能

end

end
for i=1:24 %保证备用概率的关系
b_L=[b_L;1-eta];
b_U=[b_U;1.02];%光伏那有点问题
end

for i=1:24 %储能备用限制1
b_L=[b_L;0];
b_U=[b_U;20*2];
end
for i=1:24 %储能备用限制2
b_L=[b_L;-200];
b_U=[b_U;-0.9162];
end

f_opt = 141278;
[M,N]=size(S);
x_0 = zeros(N,1);
for i=1:n

if mod(i,20)==7 |mod(i,20)==8 | mod(i,20)==9|mod(i,20)==12 |mod(i,20)==13 | mod(i,20)==14|mod(i,20)==15 |mod(i,20)==16 | mod(i,20)==17
   IntVars=[IntVars i];
end

end
for i=2024+1:2024+pronum

IntVars=[IntVars i];
end

for i=1:20*24
varible(i)=Result.x_k(i);
end
reshape(varible,20,24)
for i=20*24+1:size(Result.x_k)
zm(i-20*24)=Result.x_k(i);
end
zm
RS=0;
for i=1:n

if mod(i,20)==1 |mod(i,20)==2 | mod(i,20)==3%备用费用
  RS=RS+0.02*(Result.x_k(i));
end

end
SUFee=0;
for i=1:n
if mod(i,20)==12 |mod(i,20)==13 | mod(i,20)==14%备用费用

  SUFee=SUFee+1.5*(Result.x_k(i));
end

end
%#########
VV=S(478:918,24*20+1:825);
%#####################画图
%x=[r1 r2 r3 p1 p2 p3 i1 i2 i3 pin pout su1 su2 su3 sd1 sd2 sd3 soc Cres uload];201在2024个变量后，继续每个时段至少加11个变量
battryin=[];
battryout=[];
MT1=zeros(24,1);
MT2=zeros(24,1);
MT3=zeros(24,1);
reserve=zeros(24,1);
% for i=1:480
% j=fix(i/20)+1;
% if mod(i,20)==1|mod(i,20)==2|mod(i,20)==3
% reserve(j)=reserve(j)+Result.x_k(i);
% end
% end
% for i=1:480
% j=fix(i/20)+1;
% if mod(i,20)==19
% reserve(j)=reserve(j)+Result.x_k(i)*0.9;
% end
% end
% if eta==0.1
% plot(1:24,reserve(1:24),'bx-.','LineWidth',2);
% end
% if eta==0.05
% plot(1:24,reserve(1:24),'gd-.','LineWidth',2);
% end
% if eta==0.01
% plot(1:24,reserve(1:24),'r*:','LineWidth',2);
% end
% if eta==0
% plot(1:24,reserve(1:24),'kd--','LineWidth',2);
% end
% hold on
% xlabel('旋转备用 /KW');
% ylabel('时间/h');

for i=1:480

 j=fix(i/20)+1;
if mod(i,20)==4
    MT1(j)=MT1(j)+Result.x_k(i);
end

end
plot(1:24,MT1(1:24),'cs-','LineWidth',2);
hold on
for i=1:480

 j=fix(i/20)+1;
if mod(i,20)==5
    MT2(j)=MT2(j)+Result.x_k(i);
end

end
plot(1:24,MT2(1:24),'rd:','LineWidth',2);
hold on
for i=1:480

 j=fix(i/20)+1;
if mod(i,20)==6
    MT3(j)=MT3(j)+Result.x_k(i);
end

end
plot(1:24,MT3(1:24),'yv-.','LineWidth',2);
hold on

plot(1:24,Pg(1:24),'k*:','LineWidth',2);
hold on
% i=1:1:24;
% stairs(i,Pg(i));
% hold on
i=1:1:24;
stairs(i,PD(i),'LineWidth',2);
hold on
j=1;
k=1;
for i=1:480

if mod(i,20)==10

battryin(j)=Result.x_k(i);
j=j+1;

end

if mod(i,20)==11
battryout(k)=Result.x_k(i);
k=k+1;
end
end
bar(0.9*(battryout-battryin),0.6,'stacked');
xlabel('功率 /KW');
ylabel('时间/h');
f=Result.f_k;
return
end

【热门文章】
【热门文章】