johanlofberg
commented
1 year ago You have an old version of YALMIP installed on the machine where you are running gurobi 10.
The option nonconvex is not something you have to bother with, YALMIP automatically sets this flag since several years back.
Hello, sir! When I call the gurobi10.0.1 solver using matlab2019a, the report says Warning: Solver not applicable (gurobi). However, using MATLAB 2018b on another computer to call gurobi 9.5.1 to run the program, there is no error, is it my gurobi academic application problem or program error, or is there a version incompatibility problem between matlab and gurobi? Below is my program for your reference. ‘’‘ yalmip( 'clear' ) clear clc; yalmip; t1 = tic; Ww = [ 6.8, 6.5, 6.9, 5, 5.6, 6, 5.8, 6.2, 4.7, 4.4];%风电出力数据 Ww = 5*Ww; Wwd = Ww( 1:length( Ww ) - 1 ); Wwp = Ww( 2:length( Ww ) ); Ww1 = Wwp - Wwd;
Epart2 = [ 6, 7, 8, 26, 27, 28, 29, 30, 31, 32, 33 ];%分区内节点号
%节点号、有功功率、无功功率、重要度 Ebus22=[ 6,60,20,3; 7,200,100,4; 8,200,100,2; 26,60,25,4; 27,60,25,3; 28,60,20,3; 29,120,70,2; 30,200,600,3; 31,150,70,3; 32,210,100,4; 33,60,40,2; ]; numbus = size( Ebus22, 1 );%节点数 cg = [ 1, 0; 2, 1; 3, 3; 4, 5; 5, 10 ];%不同重要度对应的单位恢复收益 cg(:,2)=20cg(:,2); %线路号、线路首节点、尾节点、rij、xij Ebranch22=[ 6 ,6 ,7 ,0.1872,0.6188; 7 ,7 ,8 ,0.7144,0.2351; 25,6, 26,0.2030,0.1034; 26,26,27,0.2842,0.1447; 27,27,28,1.0590,0.9337; 28,28,29,0.8042,0.7006; 29,29,30,0.5075,0.2585; 30,30,31,0.9744,0.9630; 31,31,32,0.3105,0.3619; 32,32,33,0.3410,0.5302; ]; Ebranch22rij22 = zeros();%rij^2 Ebranch22xij22 = zeros();%xij^2 Ebr1 = zeros();%线路首节点的位置 Ebr2 = zeros();%线路尾节点的位置 numbr = size( Ebranch22, 1 ); for j = 1:numbr Ebranch22rij22( j, 1 ) = Ebranch22( j, 4 )^2; Ebranch22xij22( j, 1 ) = Ebranch22( j, 5 )^2; Ebr1( j ) = find( Ebus22( :, 1 )==Ebranch22( j, 2 ) ); Ebr2( j ) = find( Ebus22( :, 1 )==Ebranch22( j, 3 ) ); end P_S = zeros( numbus, 1 );%常规电源连接节点 for i = 1:numbus if Ebus22( i, 1 )==6 P_S( i ) = 1; else P_S( i ) = 0; end end P_D = zeros( numbus, 1 );%风电机组连接节点 for i = 1:numbus if Ebus22( i, 1 )==Ebus22( 10, 1 ) P_D( i ) = 1; else P_D( i ) = 0; end end P_B = zeros( numbus, numbr );%线路-节点关联矩阵 for i = 1:numbus for j = 1:numbr if find( Ebus22( i, 1 )==Ebranch22( j, 2 ) ) P_B( i, j ) = 1; elseif find( Ebus22( i, 1 )==Ebranch22( j, 3 ) ) P_B( i, j ) = -1; else P_B( i, j ) = 0; end end end lu = zeros( numbr );%各节点到电源的路径 for i = 1:numbr lu( i, 1 ) = Ebranch22( i, 1 ); j = Ebranch22( i, 2 ); n1 = i; k = 2; while (jP_S( find( Ebus22( :, 1 )==j ) )==0) for m = 1:numbr if (Ebranch22( m, 3 )==j) lu( n1, k ) = Ebranch22( m, 1 ); j = Ebranch22( m, 2 ); k = k + 1; end end end end EE1( 1 ) = 0;%各节点到电源经过的线路数 EE1( 2:numbus ) = sum( lu( 1:numbr, : )~=0, 2 ); toc( t1 ); t2 = tic; T = 10;%观察时间，单位：h UB = 12.66;%基准电压 Eplcon = 200;%线路传输上限 Esup = 1000;%常规机组总出力 Tprob = 20;%恢复步数 PEmax = Esup/Tprob;%单步恢复最大有功 Tplus = length(Ww)-1;%风机波动次数 xt = 0.1;%最小负荷恢复比例 thva = 10^(-3);%足够小的正数 risk = 0.001;%负荷恢复风险 % PMprice = 0.8;%移动应急资源供电成本 Eprice = 0.8;%常规机组供电成本 Eneed = sum( Ebus22( :, 2 ) );%分区内总待恢复负荷量 Ecould = sum( Ww ) + Esup;%分区内总固定电源有功 PEbus1 = sdpvar( numbus, Tprob, 'full' );%单步节点负荷有功恢复量 QEbus1 = sdpvar( numbus, Tprob, 'full' );%单步节点负荷无功恢复量 delta1= sdpvar( numbus, Tprob, 'full' );%计算恢复时间 pritim1= sdpvar( 1, Tprob, 'full' );%每个恢复时段时长 moment1= sdpvar( numbus, Tprob, 'full' );%每个负荷的恢复时刻 PE1 = PEmax*ones( 1, Tprob );%单步常规机组有功出力 QE1 = sdpvar( 1, Tprob, 'full' );%单步常规机组无功出力 PW = sdpvar( 1, Tprob, 'full' );%单步风电机组有功出力 El1=binvar(numbr,Tprob,'full');%线路恢复状态 kkp1=sdpvar(numbus,Tprob,'full'); Er1 = binvar( numbus, Tprob, 'full' );%单步节点负荷恢复状态：1为恢复；0为不恢复 Er_t1 = binvar( numbus, Tprob, 'full' );%节点负荷恢复状态：1为已恢复；0为从未恢复过 EW = binvar( 1, Tprob, 'full' );%风电出力状态； Epro1 = sdpvar( numbus, Tprob, 'full' );%单步节点负荷有功恢复比例 Epro_t1 = sdpvar( numbus, Tprob, 'full' );%节点负荷有功恢复比例
%下面的是需求响应用来表示分段的，可以删掉
f1 = sdpvar( numbus, Tprob, 'full' );%单步负荷恢复收益 c1 = sdpvar( 1, Tprob, 'full' );%单步负荷恢复成本 r1 = sdpvar( numbus, Tprob, 'full' );%单步负荷恢复风险 Pload1 = sdpvar( numbr, Tprob, 'full' );%单步线路传输有功 Qload1 = sdpvar( numbr, Tprob, 'full' );%单步线路传输无功 PEBUS1 = sdpvar( numbus, Tplus, 'full' );%单步风电波动导致的负荷变化量，为正表示负荷恢复量，为负表示负荷需求响应量 PLOAD1 = sdpvar( numbr, Tplus, 'full' );%单步风电波动导致的线路传输有功波动 F1 = sdpvar( numbus, Tplus, 'full' );%单步风电波动导致的负荷恢复收益变化量，为正表示负荷恢复收益，为负表示负荷需求响应补偿
C = [ ]; toc( t2 ); t3 = tic; for i = 2:numbus for p = 1:EE1( i ) Ek( i, p ) = find( Ebranch22( :, 1 )==lu( i - 1, p ) );%各节点到常规电源经过的节点 guo( i, Ek( i, p ) ) = 1; end C = [ C, (sum( PEbus1( i, : ) )<=Ebus22( i, 2 )):'sum(PEbus(i,:)';%有功总恢复量≤额定有功 (Epro1( i, : )==PEbus1( i, : )/Ebus22( i, 2 )):'Epro' %计算各节点不同时段恢复负荷的比例 ]; end C = [ C, (Epro_t1==cumsum( Epro1, 2 )):'Epro' ];%计算各节点不同时段已恢复负荷比例 %显式约束：上下限约束 C = [ C, [ -Eplcon<=Pload1, Pload1<=Eplcon ]:'Pload' ]; C = [ C, [ -Eplcon<=Qload1, Qload1<=Eplcon ]:'Qload' ]; C = [ C, [0<=PEbus1, PEbus1<=max( Ebus22( :, 2 ) ) ]:'PEbus' ]; C = [ C, [ 0<=QEbus1, QEbus1<=max( Ebus22( :, 3 ) ) ]:'QEbus' ]; C = [ C, [ 0<=QE1, QE1<=2Eplcon ]:'QE' ]; C = [ C, [ 0<=PW, PW<=max( Ww ) ]:'PW' ]; C = [ C, [ -1<=Epro1, Epro1<=1 ]:'Epro' ]; C = [ C, [ -1<=Epro_t1, Epro_t1<=1 ]:'Epro_t' ]; C = [ C, [ 0<=f1, f1<=10^6 ]:'f1' ]; C = [ C, [ 0<=c1, c1<=10^6 ]:'c1' ]; C = [ C, [ 0<=r1, r1<=10^4 ]:'r1' ]; C = [ C, [ -Eplcon<=PLOAD1, PLOAD1<=Eplcon ]:'Pload' ]; C = [ C, [ -max( Ebus22( :, 2 ) ) <=PEBUS1, PEBUS1<=max( Ebus22( :, 2 ) ) ]:'PEbus' ];%需求响应 C = [C,(pritim1==max(delta1)):'tim']; for t = 2:Tprob for i = 2:numbus C = [ C, (Er_t1( i, t - 1 )<=Er_t1( i, t )):'Er(i,t-1)<=Er(i,t)' ];%节点一旦恢复，后面时间都为已恢复 C = [ C, (EW( t - 1 )<=EW( t )):'Er(i,t-1)<=Er(i,t)'; ]; C = [C,(El1(i-1,t-1)<=El1(i-1,t)):'El(i,t-1)<=El(i,t)']; C = [C,(kkp1(i,t)==nnz(El1(Ek(i,1:EE1(i)),t-1))):'kk1(i,t-1)'];%kk1：上一时段节点i的路径的恢复数量 C = [C,(implies(~Er_t1(i,t-1)&Er1(i,t),delta1(i,t)==2(EE1(i)-kkp1(i,t))+1)):'~Er(i,t-1)&Er(i,t)'];%length(lu(i,:)-intersect(El(:,t).Ebranch11(:,1),lu(i,:))) % C = [C,(implies(delta(i,t)==2(EE(i)-kk1(i,t))+1,[Er_t(i,t-1)==0,Er(i,t)])):'~Er(i,t-1)&Er(i,t)'];%length(lu(i,:)-intersect(El(:,t).Ebranch11(:,1),lu(i,:))) C = [C,(implies(Er_t1(i,t-1)&Er1(i,t),delta1(i,t)==0.1)):'Er(i,t-1)&Er(i,t)'];%length(lu(i,:)-intersect(El(:,t).Ebranch11(:,1),lu(i,:))) end end C = [ C, (PW( 1 )==Ww( 1 )):'Er(i,t-1)<=Er(i,t)' ];%风电机组第一时步输出有功 C = [ C, (PW( 2:Tprob )==0):'Er(i,t-1)<=Er(i,t)' ];%风电机组一小时波动一次，因此前20分钟无波动 for t = 1:Tprob C = [ C, (P_BPload1( :, t ) + PEbus1( :, t )==P_SPE1( :, t ) + P_DPW( :, t ) ):'tP_SPE(:,t)' ];%节点有功平衡 C = [ C, (P_BQload1( :, t ) + QEbus1( :, t )==P_SQE1( :, t )):'tP_SQE(:,t)' ];%节点无功平衡 for i = 1:numbus C = [C,(implies(~Er1(i,t),delta1(i,t)==0)):'delta(i,t)==0', (moment1(i,t)==Er1(i,t)(sum(pritim1(1,1:t))-pritim1(1,t)+delta1(i,t))):'moment' ];%(1<=sum(Er(i,1))):'恢复时刻' C = [ C, (PEbus1( i, t )<=Er1( i, t )Ebus22( i, 2 )):'PEbus'; %恢复量与恢复状态关联 (Er1( i, t )xtEbus22( i, 2 )<=PEbus1( i, t ) ):'PEbus';%恢复量与恢复状态关联 ]; C = [ C, (implies( Epro_t1( i, t )==0, Er_t1( i, t )==0 )):'Er_t(i,t)==0';%有功恢复量为0，从未恢复过此节点 (implies( Epro_t1( i, t )>=thva, Er_t1( i, t )==1 )):'Er_t(i,t)==1' %有功恢复量不为0，此节点曾得到恢复 %不用>号的原因：yalmip不支持 ]; end end for t = 1:Tprob%假设间隔为1分钟 C = [ C, (c1( t )==(PE1( t )Eprice )(T-sum(pritim1(1,1:t))+pritim1(1,t))):'c' ];%成本=（常规机组有功出力增量电价+移动应急资源有功出力增量电价）出力时间 for i = 1:numbus C = [ C, (f1( i, t )==PEbus1( i, t )(T-sum(pritim1(1,1:t))+pritim1(1,t)-delta1(i,t))cg( Ebus22( i, 4 ), 2 )):'f' ];%收益=各节点有功出力增量时间单位恢复收益 C = [ C, (r1( i, t )==(1 - (1 - risk)^EE1( i ))f1( i, t )):'r' ];%风险=收益线路恢复风险^恢复节点所需经过的线路数 end end C = [ C, (P_BPLOAD1 + PEBUS1==P_DWw1):'P_BPLOAD1(:,t)+PEBUS1(:,t)' ];%风电出力波动时，节点有功平衡 for i = 1:numbus for t = 1:Tplus C = [ C, (sum( PEbus1( i, : ) ) + sum( PEBUS1( i, 1:t ) )<=Ebus22( i, 2 )):'sum<=Ebus11( i, 2 )' ];%节点的有功负荷总是≤额定有功 C = [ C, (0<=sum( PEbus1( i, : ) ) + sum( PEBUS1( i, 1:t ) )):'0<=sum' ];%0≤节点的有功负荷总是 C = [ C, (F1( i, t )==(PEBUS1( i, t ))(T - t)cg( Ebus22( i, 4 ), 2 )):'F' ];%若节点不参与需求响应，收益=有功负荷增量时间*负荷单位恢复收益 end end C = [ C, (implies(-thva<=PEBUS1<=thva,PEBUS1==0)):',PEBUS1==0' ];%若有功负荷的恢复量或需求响应量很小，就将其视为0
toc( t3 ); t4 = tic; Objective = sum( sum( c1 ) ) + sum( sum( r1 ) ) - sum( sum( f1 ) ) - sum( sum( F1 ) ) ;%+ c_M*sum( GM( :, time ) )/G_MEGS gurobi_stop_gap = 0.15; longest_solve_time = 600; NumericFocus = 2; options = sdpsettings( 'verbose', 2, 'solver', 'GUROBI', 'debug', 1, 'gurobi.NonConvex', 2, 'gurobi.MIPGap', gurobi_stop_gap ); sol = optimize( C, Objective, options ); %sol = optimize( C(1:3593) );if sol.problem==0;display('Feasible');else;display('Infeasible');end %用来查约束，若运行时出现警告或报错导致不能进行，用上面式子表示：删去目标函数及其他设定，仅检查约束条件，若结果为Feasible，说明约束条件没有问题，需要检查目标函数或者数据；否则约束条件有问题，检查约束条件 %检查约束条件时可以检查部分，比如：C(1:3593)为检查约束条件中1:3593条，其余不检查，约束条件中紫色的文字可以帮助确定有问题的约束条件的位置 toc( t4 ); t5 = tic; Ebus22 = double( Ebus22 ); Ebranch22 = double( Ebranch22 ); PE1 = double( PE1 ); QE1 = double( QE1 ); PW = double( PW ); % PM1 = double( PM1 ); Epro1 = double( Epro1 ); Epro_t1 = double( Epro_t1 ); Pload1 = double( Pload1 ); Qload1 = double( Qload1 ); PEbus1 = double( PEbus1 ); QEbus1 = double( QEbus1 ); Er1 = double( Er1 ); Er_t1 = double( Er_t1 ); EW = double( EW ); f1 = double( f1 ); c1 = double( c1 ); r1 = double( r1 ); PEBUS1 = double( PEBUS1 ); PLOAD1 = double( PLOAD1 ); F1 = double( F1 ); Objective = double( Objective ); toc( t5 ); ’‘’ I hope you will not hesitate to teach and reply to the students' confusion.