matlab+gurobi风电场站内按照发电量最大的风机恢复顺序问题

 风电场站内风机成组恢复时，直到最后一步才会将所有风机恢复，而中间过程不进行其他操作。

风电场10个时步，奇数时步恢复辅机，偶数时步恢复风机。

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组
NE=[5;5;5;6;6];

%% 额定常数
NT=2*NC;%最大总时步为2倍分组数
fmax=0.2/50;%最大稳态频率偏差；0.2/50P.U.
ft=0.5;%最大暂态频率偏差；
R=0.5/50;%频率变化量最大值  0.5HZ/s
SG=2500000;%外部带电域或储能的额定容量  2.5MW
KG=100;%外部带电域或储能的功频静特性系数  50p.u.
HG=25;%外部带电域或储能的惯性时间常数  5s
dfG=2;
PL=18000;
SW=2000000;
KW=25;
KL=1;
HW=10;
dfW=4;

%% 决策变量
K=sdpvar(1,NT,'full');%每一个时步的功频静特性系数
H=sdpvar(1,NT,'full');%每一个时步的惯性时间常数
P1=sdpvar(1,NT,'full');%每一个时步辅机所需要的功率
P2=sdpvar(1,NT,'full');%每一个时步新能源场站的出力
T=sdpvar(1,NT,'full');%每一时步的时间
t=sdpvar(NC,NT,'full');%每一组辅机或者风机恢复需要的时间
n=sdpvar(NC,NT,'full');%每一组已恢复的辅机数
k=sdpvar(NC,NT,'full');%每一组已恢复的风机数
Nn=sdpvar(1,NT,'full');
Nk=sdpvar(1,NT,'full');
v=binvar(NC,NT,'full');%每一个时步每一组新能源机组辅机的运行状态
u=binvar(NC,NT,'full');%每一个时步每一组新能源机组的运行状态   0表示未稳定出力  1表示稳定出力
a=binvar(11,NT,'full');


%% 目标函数
st=[];
% f=[];
for m=1:NT
                f=-sum(P2(1,m)*T(1,m));
end

%为解决目标无界性，对所有决策变量进行约束
% for m=1:NT
%     for i=1:NC
        st=st+[0<=t<=100;
               0<=n<=10;
               0<=k<=10;];
%     end
% end

% for m=1:NT
        st=st+[0<=K<=100;
               0<=H<=100;
               0<=P1<=27*PL;
               0<=P2<=27*SW;
               40<=T<=100;
               0<=Nn<=27;
               0<=Nk<=27;
               ];
% end

%% 约束条件
%根据约束对第一时步相关参数进行预