站内风机启动优化程序结果错误以及0"/>
matlab+yalmip+cplex/gurobi风电场站内风机启动优化程序结果错误以及0
目标函数：给定恢复时间内发电量最大
约束条件：稳态频率约束+惯量响应约束
风电场内有27台风机，共分为5组。
每个操作时步有两种可能：辅机启动、风机启动；
每组辅机或风机一起启动，且某组辅机启动后下一时步风机启动。
第一时步预估后，设置启动功率。
warning('off');
clear
clc
yalmip;%新能源场站内节点上的辅机及新能源机组信息
%第1列：节点号 第2列：辅机功率 第3列：额定容量 第4列：功频静特性系数（机组） 第5列：功频静特性系数（辅机） 第6列：惯性时间常数
Ebus=[
1,18000,2000000,25,1,10;
2,18000,2000000,25,1,10;
3,18000,2000000,25,1,10;
4,18000,2000000,25,1,10;
5,18000,2000000,25,1,10;
6,18000,2000000,25,1,10;
7,18000,2000000,25,1,10;
8,18000,2000000,25,1,10;
9,18000,2000000,25,1,10;
10,18000,2000000,25,1,10;
11,18000,2000000,25,1,10;
12,18000,2000000,25,1,10;
13,18000,2000000,25,1,10;
14,18000,2000000,25,1,10;
15,18000,2000000,25,1,10;
16,18000,2000000,25,1,10;
17,18000,2000000,25,1,10;
18,18000,2000000,25,1,10;
19,18000,2000000,25,1,10;
20,18000,2000000,25,1,10;
21,18000,2000000,25,1,10;
22,18000,2000000,25,1,10;
23,18000,2000000,25,1,10;
24,18000,2000000,25,1,10;
25,18000,2000000,25,1,10;
26,18000,2000000,25,1,10;
27,18000,2000000,25,1,10;
];
ND=length(Ebus);%有27个节点%新能源场站内线路信息
%第1列：线路号 第2列：首节点 第3列：末节点 第4列：rij 第5列：xij
Ebranch=[1,0,1;
2,1,2;
3,2,3;
4,3,4;
5,4,5;
6,0,6;
7,6,7;
8,7,8;
9,8,9;
10,9,10;
11,0,11;
12,11,12;
13,12,13;
14,13,14;
15,14,15;
16,0,16;
17,16,17;
18,17,18;
19,18,19;
20,19,20;
21,20,21;
22,0,22;
23,22,23;
24,23,24;
25,24,25;
26,25,26;
27,26,27;
];
NC=length(find(0==Ebranch(:,2)));%有5组%% 额定常数
NT=2*NC;%最大总时步为2倍分组数
fmax=0.2/50;%最大频率偏差；0.2/50P.U.
R=0.5/50;%频率变化量最大值 0.5HZ/s
SG=2500000;%外部带电域或储能的额定容量 2.5MW
KG=50;%外部带电域或储能的功频静特性系数 50p.u.
HG=5;%外部带电域或储能的惯性时间常数 5s
PL=Ebus(:,2);
SW=Ebus(:,3);
KW=Ebus(:,4);
KL=Ebus(:,5);
HW=Ebus(:,6);%% 决策变量
K=sdpvar(1,NT,'full');%每一个时步的功频静特性系数
H=sdpvar(1,NT,'full');%每一个时步的惯性时间常数
P1=sdpvar(1,NT,'full');%每一个时步辅机所需要的功率
P2=sdpvar(1,NT,'full');%每一个时步新能源场站的出力
T=sdpvar(NT,1,'full');%每一时步的时间
n=sdpvar(1,NT,'full');
k=sdpvar(1,NT,'full');
v=binvar(ND,NT,'full');%每一个时步每一台新能源机组辅机的运行状态
u=binvar(ND,NT,'full');%每一个时步每一台新能源机组的运行状态 0表示未稳定出力 1表示稳定出力
a=binvar(13,NT,'full');
b1=binvar(2,NT,'full');
b2=binvar(2,NT,'full');
b3=binvar(2,NT,'full');
b4=binvar(2,NT,'full');
b5=binvar(2,NT,'full');
b6=binvar(2,NT,'full');
b7=binvar(2,NT,'full');
b8=binvar(2,NT,'full');
b9=binvar(2,NT,'full');
b10=binvar(2,NT,'full');
S1=sdpvar(1,NT,'full');%% 目标函数
st=[];
Z=ones(1,NT);for m=1:NTf=-sum(P2(1,m)*T(m,1));end%% 约束条件
%根据约束对第一时步相关参数进行预设
for m=1for i=1:NDst=st+[S1(1,m)*(sum(v(i,m)*PL(i,:))+SG*Z(1,m))==1;K(1,m)==(sum(v(i,m)*KL(i,:)*PL(i,:))+KG*SG*Z(1,m))*S1(1,m);%./(sum(v(i,m).*PL(i,:))+SG*Z(1,m)H(1,m)==HG;P1(1,m)==90000;P2(1,m)==0;v(1:5,m)==1;v(6:27,m)==0;u(:,m)==0;T(m,1)==85
% P1(1,m)<=fmax*K(1,m)*SG;
% P1(1,m)<=2*R*H(1,m)*SG;];
% st=st+[implies(P1(1,m)>=90000,P1(1,m)==90000&P2(1,m)==0&v(1:5,m)==1&v(6:27,m)==0&u(:,m)==0&T(m,1)==85)];end
end%稳态频率约束
for m=2:NTfor i=1:NDst=st+[K(1,m)==(sum(u(i,m)*KW(i,:)*SW(i,:))+sum(v(i,m)*KL(i,:)*PL(i,:))+KG*SG*Z(1,m))/(sum(u(i,m)*SW(i,:))+sum(v(i,m)*PL(i,:))+SG*Z(1,m));P1(1,m)<=fmax*K(1,m)*(sum(u(i,m).*SW(i,:))+SG*Z(1,m));P2(1,m)<=fmax*K(1,m)*(sum(u(i,m).*SW(i,:)))];end
end%惯量响应约束
for m=2:NTfor i=1:NDst=st+[H(1,m)==(sum(u(i,m)*HW(i,:).*SW(i,:))+HG*SG*Z(1,m))/(sum(u(i,m)*SW(i,:))+SG*Z(1,m));P1(1,m)<=2*R*H(1,m)*(sum(u(i,m)*SW(i,:))+SG*Z(1,m));P2(1,m)<=2*R*H(1,m)*(sum(u(i,m)*SW(i,:)))];end
end%启动功率以及新能源机组出力等式关系
for m=2:NTfor i=1:NDst=st+[P1(1,m)==sum(v(i,m)*SW(i,:));0.1*sum(u(i,m)*SW(i,:))<=P2(1,m),P2(1,m)<=sum(u(i,m)*SW(i,:))];end
end%对每一时步的01决策变量进行计数
for m=1:NTfor i=1:NDst=st+[n(1,m)==sum(v(i,m));k(1,m)==sum(u(i,m))];end
end%时间约束 新能源辅机恢复时，新能源所在线路及变压器恢复考虑在这一时步（40+20+5*n）新能源机组恢复时（40+10）
for m=2:NTst=st+[sum(a(:,m))==1;implies(a(1,m),n(1,m)-n(1,m-1)==5,T(m,1)==85);implies(a(2,m),n(1,m)-n(1,m-1)==10,T(m,1)==85);implies(a(3,m),n(1,m)-n(1,m-1)==15,T(m,1)==85);implies(a(4,m),n(1,m)-n(1,m-1)==6,T(m,1)==90);implies(a(5,m),n(1,m)-n(1,m-1)==11,T(m,1)==90);implies(a(6,m),n(1,m)-n(1,m-1)==12,T(m,1)==90);implies(a(7,m),n(1,m)-n(1,m-1)==16,T(m,1)==90);implies(a(8,m),n(1,m)-n(1,m-1)==17,T(m,1)==90);implies(a(9,m),n(1,m)-n(1,m-1)==21,T(m,1)==90);implies(a(10,m),n(1,m)-n(1,m-1)==22,T(m,1)==90);implies(a(11,m),n(1,m)-n(1,m-1)==27,T(m,1)==90);%启动辅机时不同时步的时间implies(a(12,m),k(1,m)-k(1,m-1)>0,T(m,1)==50);%只要启动新能源机组，该时步时间为50simplies(a(13,m),n(1,m)-n(1,m-1)==0&k(1,m)-k(1,m-1)==0,T(m,1)==50);%启动完成后的时步皆为50s];
end%0/1决策变量约束
for m=2:NTst=st+[sum(b1(:,m))==1;implies(b1(2,m),v(1,m)==1,v(2:5,m)==1);implies(b1(2,m),v(1,m)==0,v(2:5,m)==0)]; st=st+[sum(b2(:,m))==1;implies(b2(2,m),v(6,m)==1,v(7:10,m)==1);implies(b2(2,m),v(6,m)==0,v(7:10,m)==0)]; st=st+[sum(b3(:,m))==1;implies(b3(2,m),v(11,m)==1,v(12:15,m)==1);implies(b3(2,m),v(11,m)==0,v(12:15,m)==0)]; st=st+[sum(b4(:,m))==1;implies(b4(2,m),v(16,m)==1,v(17:21,m)==1);implies(b4(2,m),v(16,m)==0,v(17:21,m)==0)]; st=st+[sum(b5(:,m))==1;implies(b5(2,m),v(22,m)==1,v(23:27,m)==1);implies(b5(2,m),v(22,m)==0,v(23:27,m)==0)]; %成组启动，当一组中第一台机组辅机启动，这一组中的机组辅机也启动% st=st+[sum(b6(:,m))==1;
% implies(b6(2,m),u(1,m)==1,u(2:5,m)==1);
% implies(b6(2,m),u(1,m)==0,u(2:5,m)==0)];
% st=st+[sum(b7(:,m))==1;
% implies(b7(2,m),u(6,m)==1,u(7:10,m)==1);
% implies(b7(2,m),u(6,m)==0,u(7:10,m)==0)];
% st=st+[sum(b8(:,m))==1;
% implies(b8(2,m),u(11,m)==1,u(12:15,m)==1);
% implies(b8(2,m),u(11,m)==0,u(12:15,m)==0)];
% st=st+[sum(b9(:,m))==1;
% implies(b9(2,m),u(16,m)==1,u(17:21,m)==1);
% implies(b9(2,m),u(16,m)==0,u(17:21,m)==0)];
% st=st+[sum(b10(:,m))==1;
% implies(b10(2,m),u(22,m)==1,u(23:27,m)==1);
% implies(b10(2,m),u(22,m)==0,u(23:27,m)==0)]; %成组启动，当一组中第一台机组启动，这一组中的机组也启动
end%辅机的状态量与新能源机组的状态量间的约束
for m=2:NTfor i=1:NDst=st+[v(i,m-1)<=v(i,m);%辅机启动后不会停止
% u(i,m-1)<=u(i,m);%新能源机组启动后也不会停止u(i,m)==v(i,m-1)];%辅机启动后下一时步新能源机组启动end
end%% 求解
ops=sdpsettings('solver', 'gurobi');
result=solvesdp(st,f);
double(f)%% 结论
A=double(P1);
B=double(P2);
C=double(u);
D=double(v);
F=double(K);
G=double(H);问题：
* Starting YALMIP integer branch & bound.
* Lower solver : fmincon-standard
* Upper solver : rounder
* Max iterations : 300
Warning : The continuous relaxation may be nonconvex. This means
that the branching process is not guaranteed to find a
globally optimal solution, since the lower bound can be
invalid. Hence, do not trust the bound or the gap...
Node Upper Gap(%) Lower Open
1 : Inf Inf NaN 0 Infeasible problem
+ 1 Finishing. Cost: Inf
ans =
0
更多推荐
matlab+yalmip+cplex/gurobi风电场站内风机启动优化程序结果错误以及0
发布评论