After upgrading the gurobi version from 9.5 to 10.0.1, the model solution result is wrong?

[Electricman23]

Hi Johan and everyone.
In the process of running the above program, I encountered problems running on different computers with different results. Among them, in MATLAB 2018b, version: Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64), which can run the correct program; In Matlab2019b, version: Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (win64), some results cannot run the correct result, For example, the decision variable f1, whose value is the upper limit of my constraint, which is obviously incorrect. So I would like to ask you, is it an issue with my version mismatch, or something else.If for other reasons, please help answer, thank you.
The following is the log of calling gurobi version 10.0.1 in matlabR2019b to run the program：

Set parameter IntFeasTol to value 1e-06
Set parameter MIPGap to value 0.3
Set parameter ScaleFlag to value 1
Set parameter BarIterLimit to value 2000000000
Set parameter CrossoverBasis to value 0
Set parameter MinRelNodes to value 0
Set parameter NodeMethod to value 1
Set parameter PumpPasses to value 0
Set parameter ZeroObjNodes to value 0
Set parameter AggFill to value 10
Set parameter NonConvex to value 2
Set parameter PreSparsify to value 0
Set parameter TuneTrials to value 2
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (win64)

CPU model: Intel(R) Core(TM) i5-6300HQ CPU @ 2.30GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 4 physical cores, 4 logical processors, using up to 4 threads

Optimize a model with 29449 rows, 18322 columns and 128677 nonzeros
Model fingerprint: 0xfbaa4730
Variable types: 12905 continuous, 5417 integer (5417 binary)
Coefficient statistics:
  Matrix range     [1e-05, 4e+06]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e-05, 4e+06]
Presolve removed 24049 rows and 16197 columns
Presolve time: 0.42s
Presolved: 5400 rows, 2125 columns, 19895 nonzeros
Variable types: 1082 continuous, 1043 integer (1042 binary)

Root relaxation: objective -2.394388e+08, 1769 iterations, 0.07 seconds (0.03 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -2.394e+08    0   24          - -2.394e+08      -     -    0s
     0     0 -2.394e+08    0   60          - -2.394e+08      -     -    0s
     0     0 -2.394e+08    0   52          - -2.394e+08      -     -    0s
     0     0 -2.394e+08    0    6          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0   10          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    3          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    6          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    3          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    3          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    3          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    4          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    4          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    4          - -2.394e+08      -     -    1s
     0     0 -2.394e+08    0    4          - -2.394e+08      -     -    1s
     0     2 -2.394e+08    0    4          - -2.394e+08      -     -    2s
*   49    80              16    -2.39439e+08 -2.394e+08  0.00%  41.4    2s

Cutting planes:
  Learned: 10
  Gomory: 5
  Cover: 35
  Implied bound: 13
  Clique: 22
  MIR: 35
  Flow cover: 21
  GUB cover: 1
  RLT: 1
  Relax-and-lift: 6

Explored 106 nodes (10142 simplex iterations) in 2.60 seconds (1.27 work units)
Thread count was 4 (of 4 available processors)

Solution count 1: -2.39439e+08 
No other solutions better than -2.39439e+08

Optimal solution found (tolerance 3.00e-01)
Warning: max constraint violation (9.9000e-06) exceeds tolerance
Best objective -2.394387793003e+08, best bound -2.394387793003e+08, gap 0.0000%

Attached is the program I have written for your reference.

yalmip('clear')
clear
clc;
yalmip;
t1=tic;

Ww=[5.8,6.7,5.8,5.1,6.3,5];%12小时预测值%,6.2,6,4.1,6,7,6.8,6.5,6.9,5,5.6,6,5.8,6.2,4.7,3.3,4.4,5.6,5%h：24小时预测值
Ww=5*Ww;%sum(Ww)=708
Wwd=Ww(1:length(Ww)-1);
Wwp=Ww(2:length(Ww));
Ww1=Wwp-Wwd;

Epart=[1,2,3,4,5,19,20,21,22,23,24,25];%
DSu=[0,0,0,1,0,0,1,1,0,0,0,0];%可中断负荷：0：不可中断；1：可中断
LJu=[0,0,0,0,0,0,0,0,1,1,0,0];%候选连接点：0：不可连接；1：可连接
PL=((1-DSu).*Epart)';
PL(find(PL==0))=[];%不可中断负荷
PR=(DSu.*Epart)';
PR(find(PR==0))=[];%可中断负荷
LJ=(LJu.*Epart)';
LJ(find(LJ==0))=[];
nN=length(PL);
nDS=length(PR);%sum(DSu);%可中断负荷数量
nLJ=length(LJ);
nMEPS=1;
Fre=1;%负荷可中断次数

Ebus11=[1,0,0,5;
2,100,60,2;
3,90,40,3;
4,120,80,2;
5,60,30,5;
19,90,40,5;
20,90,40,2;
21,90,40,2;
22,90,40,3; 
23,90,40,4; 
24,420,200,3;
25,420,200,3;
    ];
numbus=size(Ebus11,1);%包括与输电网连接的节点
cg=[1  0
    2  1
    3  3
    4  5
    5  10
    ];

Ebranch11=[1 ,1 ,2 ,0.0922,0.0407;
2 ,2 ,3 ,0.4930,0.2511;
3 ,3 ,4 ,0.3660,0.1864;
4 ,4 ,5 ,0.3811,0.1941;
18,2, 19,0.1640,0.1565;
19,19,20,1.5042,1.3554;
20,20,21,0.4095,0.4784;
21,21,22,0.7089,0.9373;
22,3, 23,0.4512,0.3083;
23,23,24,0.8980,0.7091;
24,24,25,0.8960,0.7011;
    ];
Ebranch11rij22=zeros();Ebranch11xij22=zeros();Ebr1=zeros();Ebr2=zeros();
numbr=size(Ebranch11,1);
for j=1:numbr
    Ebranch11rij22(j,1)=Ebranch11(j,4)^2;
    Ebranch11xij22(j,1)=Ebranch11(j,5)^2;
    Ebr1(j)=find(Ebus11(:,1)==Ebranch11(j,2));%线路首节点在节点矩阵的位置
    Ebr2(j)=find(Ebus11(:,1)==Ebranch11(j,3));%线路尾节点在节点矩阵的位置
end

P_S=zeros(numbus,1);
for i=1:numbus
      if Ebus11(i,1)==1
       P_S(i)=1;
    else
       P_S(i)=0;
      end
end


%假设间歇性能源连接在节点12
P_D=zeros(numbus,1);
for i=1:numbus
      if Ebus11(i,1)==Ebus11(12,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(Ebus11(i,1)==Ebranch11(j,2))  
       P_B(i,j)=1;
    elseif find(Ebus11(i,1)==Ebranch11(j,3))
       P_B(i,j)=-1;
    else
       P_B(i,j)=0;
    end
  end
end
lu=zeros(numbr);                   %对矩阵lu进行初始化
for i=1:numbr
    lu(i,1)=Ebranch11(i,1);                 %将路径矩阵lu的第一列储存为支路号
   j=Ebranch11(i,2);
        n1=i;
        k=2;
        while(j*P_S(find(Ebus11(:,1)==j))==0)            %一直搜寻到母节点0结束   %值运算：&&、||；向量运算：&、|（两边都会计算）
        for m=1:numbr
        if(Ebranch11(m,3)==j)
            %if(Ebranch11(m,3)==find(Ebus11(:,1)==j))          %在末节点定位j
            lu(n1,k)=Ebranch11(m,1);     
            j=Ebranch11(m,2);          %将首节点赋给j,然后再在末节点定位j
            k=k+1;
        end
        end
        end
   % end
end
EE1(1)=0;
EE1(2:numbus)=sum(lu(1:numbr,:)~=0,2);%length(lu(i,:)~=0);
toc(t1);
t2=tic;
T=10*60;%评估时长3h
UB=12.66;%基准电压
Eplcon=200;%设线路传输功率上限为600kw；限制了线路传输功率上线后，电网输入功率也被限制！
Esup=1000;%输电网总有功容量
Tprob=20;%ceil(Esup/(PEmax+min(Ww)));%向上取整
PEmax=Esup/Tprob;%输电母线单次输入最大功率
Tplus=5;%每段半小时，还有五段
P_MEPS=20;
t_MEPS=T*ones(1,Tprob);%[15  10];
t_MEPS(9)=10;t_MEPS(10)=8;
PMprice=0.8;
xt=0.1;%负荷恢复最小量
Eprice=0.8;%从电网购买电功率的价格为0.8元/kw/10000*10/60
thva=10^(-2);
risk=0.001;%线路恢复风险
Eneed=sum(Ebus11(:,2));
Ecould=sum(Ww)+Esup;
sdpvar tt;%tt=sdpvar(1,1,'full');%length(delta(1,:));%总时步数
delta1=sdpvar(numbus,Tprob,'full');

tDS1=sdpvar(nDS,Fre+1,'full');%各可中断负荷的中断开始时刻
tDS1(:,1)=PR;
PEbus1=sdpvar(numbus,Tprob,'full');%节点恢复负荷量，需求侧响应节点可以为负数
QEbus1=sdpvar(numbus,Tprob,'full');%节点恢复负荷量
PE1=PEmax*ones(1,Tprob);%电网输入功率
QE1=sdpvar(1,Tprob,'full');%电网输入功率
PW=sdpvar(1,Tprob,'full');%间歇性电源输入功率
PM1=sdpvar(numbus,Tprob,'full');%MEPS输入功率
V_121=sdpvar(numbus-1,Tprob,'full');%节点电压的平方
diV221=sdpvar(numbus,Tprob,'full');%节点电压的平方的倒数
V1=sdpvar(numbus,Tprob,'full');%节点电压的平方
V221=[UB*UB*ones(1,Tprob);V_121];
El1=binvar(numbr,Tprob,'full');%线路恢复状态
Er1=binvar(numbus,Tprob,'full');%线路恢复状态
Er_t1=binvar(numbus,Tprob,'full');%线路恢复状态
EW=binvar(1,Tprob,'full');%间歇性能源供电状态
EM1=binvar(numbus,Tprob,'full');%MEPS连接状态
EDS1=binvar(numbus,Tplus,'full');%需求响应标志：-1为响应，0为无动作，1为不响应
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');%恢复风险 %改之前为sdpvar(numbr,Tprob,'full');
kkp1=sdpvar(numbus,Tprob,'full');
Pload1=sdpvar(numbr,Tprob,'full');
Pload221=sdpvar(numbr,Tprob,'full');
Qload1=sdpvar(numbr,Tprob,'full');
Qload221=sdpvar(numbr,Tprob,'full');
Lijload221=sdpvar(numbr,Tprob,'full');
pritim1=sdpvar(1,Tprob,'full');%各操作时段的时长
moment1=sdpvar(numbus,Tprob,'full');%各节点的操作完成时刻

PEBUS1=sdpvar(numbus,Tplus,'full');%节点恢复负荷量，需求侧响应节点可以为负数
PLOAD1=sdpvar(numbr,Tplus,'full');
F1=sdpvar(numbus,Tplus,'full');%恢复收益


C=[];
%第一个时段结束后已经有负荷恢复
C = [C,(delta1(1,:)==0):'delta(1,:)==0',
        (V1(1,:)==UB):'V(1,:)',
        (Er1(1,:)==0):'Er(1,:)==0'
        (Er_t1(1,:)==0):'Er_t(1,:)==0'
        (PEbus1(1,:)==0):'Er(1,:)==0'
        (EDS1(1,:)==0):'Er(1,:)==0'
        (EM1(1:8,:)==0):'EM(1,:)==0'
        (EM1(11:12,:)==0):'EM(1,:)==0'
        (PM1(1:8,:)==0):'EM(1,:)==0'
        (PM1(11:12,:)==0):'EM(1,:)==0'
        ];
toc(t2);
t3=tic;

for i=2:numbus
      for p=1: EE1(i)
         %Ek为每个节点的路径所经过的线路在Ebranch11第一列中取的位置
          Ek(i,p)=find(Ebranch11(:,1)==lu(i-1,p)); %包括节点1的时候：Ek(i,p)=find(Ebranch11(:,1)==lu(i-1,p));
          guo(i,Ek(i,p))=1;
      end
     C = [C,(implies(Er1(i,1),delta1(i,1)==2*EE1(i)+1)):'delta(i,1)'];
     %负荷投入总量
     C = [C,(sum(PEbus1(i,:))<=Ebus11(i,2)):'sum(PEbus(i,:)',
         (Epro1(i,:)==PEbus1(i,:)/Ebus11(i,2)):'Epro',
         ];%+sum(DG(:,2),1)%     5%是什么的：机组装机容量？域内负荷？j
        C = [C,(implies(~El1(i-1,1),V221(Ebr2(i-1),1)==0)):'V22==0',
            (implies(V221(i,1)~=0&Er_t1(i,1),diV221(i,1)*V221(i,1)==1)):'diV22',
              (implies(V221(i,1)==0,diV221(i,1)==0)):'diV22',
              (V221(i,1)==V1(i,1)^2):'V22(i,1)==V(i,1)^2',%*V(i,t)
              ];
end
 C = [C,(Epro_t1==cumsum(Epro1,2)):'Epro'];

%0-1变量与连续变量的乘积线性化，不是MC包络：https://blog.csdn.net/zhangxd212489/article/details/122263798
%线路传输约束
 C = [C,[-Eplcon<=Pload1,Pload1<=Eplcon]:'Pload'];
 C = [C,[-Eplcon<=Qload1,Qload1<=Eplcon]:'Qload'];
 C = [C,[0<=Lijload221,Lijload221<=10^6]:'Lijload22'];
 C = [C,[-max(Ebus11(:,2))<=PEbus1,PEbus1<=max(Ebus11(:,2))]:'PEbus'];%负荷恢复量
 C = [C,[0<=QEbus1,QEbus1<=max(Ebus11(:,3))]:'QEbus'];
 C = [C,[0<=QE1,QE1<=2*Eplcon]:'QE'];
 C = [C,[0<=PW,PW<=max(Ww)]:'PW'];
 C = [C,[0<=PM1,PM1<=P_MEPS]:'PM1'];
 %电压约束
 C = [C,[0<=V_121,V_121<=1.1*1.1*UB*UB]:'V_12'];
 C = [C,[0<=diV221,diV221<=1.1*1.1*UB*UB]:'diV22'];%C = [C,(1/(UB*1.1*1.1)<=diV22<=1/(UB*0.9*0.9)):'V'];
 C = [C,[0<=V1,V1<=1.1*UB]:'V'];%C = [C,(UB*0.9<=V<=UB*1.1):'V'];%用有名值便于计算功率值V22
 %显式约束
 C = [C,[0<=kkp1,kkp1<=max(EE1)]:'kk1'];%要测试kk1是否在0和10之间，如果所有参数均为数值标量，则使用(0<=kk1)&(kk1<=10)，否则使用(0<=kk1)&&(kk1<=10)
 C = [C,[0<=delta1,delta1<=Tprob]:'delta'];
 C = [C,[0<=pritim1,pritim1<=max(EE1)+1]:'pritim'];
 C = [C,[0<=moment1,moment1<=Tprob*max(EE1)]:'moment'];
 C = [C,[-1<=Epro1,Epro1<=1]:'Epro'];
 C = [C,[-1<=Epro_t1,Epro_t1<=1]:'Epro_t'];
 C = [C,[0<=f1,f1<=10^6]:'f'];
 C = [C,[0<=c1,c1<=10^6]:'c'];
 C = [C,[0<=r1,r1<=10^4]:'r'];
 C = [C,[-Eplcon<=PLOAD1,PLOAD1<=Eplcon]:'Pload'];
 C = [C,[-PEmax<=PEBUS1,PEBUS1<=PEmax]:'PEbus'];%负荷恢复量
 C = [C,[-10^4<=F1,F1<=10^6]:'f'];

C = [C,(pritim1==max(delta1)):'tim'
%       (0.1<=tim):'时步不为0'
];%
 for t=2:Tprob
    for i=2:numbus%n_N%-n_DSu  %防止出现%Ek(1,1)=0的情况，已经将节点1的负荷删去：Er(1,:)=0
      %节点操作时序约束
      C = [C,(Er_t1(i,t-1)<=Er_t1(i,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)',
       (implies(~El1(i-1,t),V221(Ebr2(i-1),t)==0)):'V22(Ebr2(i-1),t)==0'];
      C = [C,(kkp1(i,t)==nnz(El1(Ek(i,1:EE1(i)),t-1))):'kk1(i,t-1)'%kk1：上一时段节点i的路径的恢复数量
];%kk1：上一时段节点i的路径的恢复数量
      C = [C,
          (EW(t-1)<=EW(t)):'Er(i,t-1)<=Er(i,t)'
          (EM1(i,t-1)<=EM1(i,t)):'Er(i,t-1)<=Er(i,t)'
          ];%不能用find！！！！diff(sign(diff(PEbus(i,a))))==2；a为sdpvar
    end
 end

for t=1:Tprob    
C = [C,(P_B*Pload1(:,t)+PEbus1(:,t)==P_S*PE1(:,t)+P_D*PW(:,t)+PM1(:,t)):'t*P_S*PE(:,t)'];%*Ebus11(:,2) +P_G*PDG
C = [C,(P_B*Qload1(:,t)+QEbus1(:,t)==P_S*QE1(:,t)):'t*P_S*QE(:,t)'];%*Ebus11(:,2) +P_G*QDG+P_D*QW(:,t)
    for i=2:numbus
          C = [C,(implies(V221(i,t)~=0&Er_t1(i,t),diV221(i,t)*V221(i,t)==1)):'diV22',
              (implies(V221(i,t)==0,diV221(i,t)==0)):'diV22',
              (V221(i,t)==V1(i,t)^2):'V22(i,t)==V(i,t)^2',];%*V(i,t)
           C = [C,(PEbus1(i,t)<=Er1(i,t)*Ebus11(i,2)):'PEbus',
              (implies(Er1(i,t),Er1(i,t)*xt*Ebus11(i,2)<=PEbus1(i,t))):'PEbus',%调用 "varargout" 时，未对一个或多个输出参数赋值。
              (implies(~Er1(i,t),PEbus1(i,t)==0)):'PEbus',
              (0<=sum(PEbus1(i,1:t))):'PEbus',%调用 "varargout" 时，未对一个或多个输出参数赋值。
           ]; 
      C = [C,(implies(sum(PE1(1:t))+sum(sum(PM1(:,1:t)))>=0.3*Eneed,EW(t)==1)):'Er(i,t-1)<=Er(i,t)'];
      C = [C,(implies(sum(PE1(1:t))+sum(sum(PM1(:,1:t)))<=(0.3-thva)*Eneed,EW(t)==0)):'Er(i,t-1)<=Er(i,t)'];
      C = [C,(implies(EW(t),PW(t)==Ww(1))):'Er(i,t-1)<=Er(i,t)'];
      C = [C,(implies(~EW(t),PW(t)==0)):'Er(i,t-1)<=Er(i,t)'];
           C = [C,(implies(EM1(i,t),t_MEPS(i)<=sum(pritim1(1:t)))):'EM'
          (sum(EM1(:,t))<=nMEPS):'EM'
          (implies(EM1(i,t),PM1(i,t)==P_MEPS)):'EM'
          (implies(~EM1(i,t),PM1(i,t)==0)):'EM'    %i=1没有定但没有报implies的警告
           ]; 
      C = [C,(implies(Epro_t1(i,t)==0,Er_t1(i,t)==0)):'Er_t(i,t)==0'
             (implies(Epro_t1(i,t)>=thva,Er_t1(i,t)==1)):'Er_t(i,t)==1'
             ];
      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))):'恢复时刻'
     
    end
    for j=1:numbr
   C = [C,(Pload221(j,t)==Pload1(j,t)^2):'Pload22',%*Pload(j,t)
       (Qload221(j,t)==Qload1(j,t)^2):'Qload22',%*Qload(j,t)
       (Lijload221(j,t)==(Pload221(j,t)+Qload221(j,t))*diV221(Ebr1(j),t)):'Lijload22',
       (implies(El1(j,t),V221(Ebr1(j),t)-V221(Ebr2(j),t)==2*(Pload1(j,t)*Ebranch11(j,4)+Qload1(j,t)*Ebranch11(j,5))...
       -(Ebranch11rij22(j,1)+Ebranch11xij22(j,1))*Lijload221(j,t))):'V22(Ebr1(j),t)-V22(Ebr2(j),t)',
       ];
       %线路操作时序约束
   C = [C, %(implies(El(j,t),El(j,t:tt)==1)):'El(j,t+1:end)',%El(find(El(:,t)~=1),t)==0,
       (Pload1(j,t)<=El1(j,t)*Eplcon):'Pload(j,t)<=',
           (-El1(j,t)*Eplcon<=Pload1(j,t)):'<=Pload(j,t)',
       (Qload1(j,t)<=El1(j,t)*Eplcon):'Pload(j,t)<=',
           (-El1(j,t)*Eplcon<=Qload1(j,t)):'<=Pload(j,t)',
           ];
   C = [C,(implies(Pload1(j,t)==0,El1(j,t)==0)):'El(j,t)==0'%不能改成implies(~El(j,t),Pload(j,t)==0)，否则El和Pload不相符
       ];%implies(Pload(j,t)==0,El(j,t)==0)%implies(~El(j,t),Pload(j,t)==0)
    end
end

for t=1:Tprob    
  C = [C,(c1(t)==(PE1(t)*Eprice+sum(PM1(:,t))*PMprice)*(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(Ebus11(i,4),2)):'f'];%PEbus：分别为每个时段恢复量 
   C = [C,(r1(i,t)==(1-(1-risk)^EE1(i))*f1(i,t)):'r'];%
    end
end
    C = [C,(P_B*PLOAD1+PEBUS1==P_D*Ww1):'P_B*PLOAD1(:,t)+PEBUS1(:,t)'];%*Ebus11(:,2) +P_G*PDG+PM1(:,t)

for i=1:numbus
              
    C = [C,(EDS1(i,:)<=DSu(i)):'EDS(i,:)<=DSu(i)',];%*Ebus11(:,2) +P_G*PDG+PM1(:,t)
      for t=1:Tplus
    C = [C,(sum(PEbus1(i,:))+sum(PEBUS1(i,1:t))<=Ebus11(i,2)):'P_B*PLOAD1(:,t)+PEBUS1(:,t)'];%*Ebus11(:,2) +P_G*PDG+PM1(:,t)
    C = [C,(0<=sum(PEbus1(i,:))+sum(PEBUS1(i,1:t))):'P_B*PLOAD1(:,t)+PEBUS1(:,t)'];%*Ebus11(:,2) +P_G*PDG+PM1(:,t)
    C = [C,
        (implies(EDS1(i,t),PEBUS1(i,t)<=-thva)):'PEbus',
        (implies(~EDS1(i,t),0<=PEBUS1(i,t))):'PEbus',%调用 "varargout" 时，未对一个或多个输出参数赋值。
        ];%*Ebus11(:,2) +P_G*PDG+PM1(:,t)
   C = [C,(F1(i,t)==(PEBUS1(i,t))*(T/6)*cg(Ebus11(i,4),2)):'F'];%PEbus：分别为每个时段恢复量 
      end
end

toc(t3);
t4=tic;
 %目标函数

Objective=sum(sum(c1))+sum(sum(r1))-sum(sum(f1))-sum(sum(F1));%净收益的负值-sum(sum(faft))+sum(sum(R1))
%%%%%调用求解器CPLEX%%%%%%%%
gurobi_stop_gap = 0.15; % 设置求解混合整数规划的 Gap 为10^(-4)
longest_solve_time = 600; % 设置最长求解时间为 300 秒
NumericFocus = 2;
options = sdpsettings('verbose',2,'solver','GUROBI+', 'debug', 1, 'gurobi.NonConvex', 2,...
    'gurobi.MIPGap', gurobi_stop_gap);%cplex%,'gurobi.NumericFocus',NumericFocus,'gurobi.timelimit',longest_solve_time
%%%%%%%%%%%优化%%%%%%%%%%%
sol = optimize(C,Objective,options);%Constraints为约束；Objective为目标函数；options为属性。
%用来检查约束
%sol = optimize(C(1:9178));if sol.problem==0;display('Feasible');else;display('Infeasible');end
%第一个是1017(120:6997)

toc(t4);
t5=tic;

tt=double(tt);
Ebus11=double(Ebus11);
Ebranch11=double(Ebranch11);
kkp1=double(kkp1);
delta1=double(delta1);
PE1=double(PE1);
QE1=double(QE1);
PW=double(PW);
PM1=double(PM1);
diV221=double(diV221);
V221=double(V221);
V1=double(V1);
Lijload221=double(Lijload221);
Epro1=double(Epro1);
Epro_t1=double(Epro_t1);
Pload1=double(Pload1);
Pload221=double(Pload221);
Qload1=double(Qload1);
Qload221=double(Qload221);
PEbus1=double(PEbus1);
QEbus1=double(QEbus1);
El1=double(El1);
Er1=double(Er1);
Er_t1=double(Er_t1);
EW=double(EW);
EM1=double(EM1);
EDS1=double(EDS1);
tDS1=double(tDS1);
pritim1=double(pritim1);
moment1=double(moment1);
f1=double(f1);
c1=double(c1);
r1=double(r1);
PEBUS1=double(PEBUS1);
PLOAD1=double(PLOAD1);
F1=double(F1);
Objective=double(Objective);
toc(t5);

[Electricman23]

Hello, according to the process of looking for the problem, I changed the solver version, replaced the original 10.0.1 with the 9.5.1 version, and found that the result is still not wrong.
Also, I use MATLAB2018b and solver gurobi9.5.1 on another computer, and the calendar after running the program is listed below for your reference:

Set parameter MIPGap to value 0.15
Set parameter CrossoverBasis to value 0
Set parameter NodefileDir to value ""
Set parameter NonConvex to value 2
Set parameter PreSOS2BigM to value 0
Set parameter TuneTrials to value 3
Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 28636 rows, 10599 columns and 119138 nonzeros
Model fingerprint: 0xd3d5c8f6
Model has 1832 quadratic constraints
Variable types: 4515 continuous, 6084 integer (6084 binary)
Coefficient statistics:
  Matrix range     [1e-05, 4e+06]
  QMatrix range    [8e-01, 1e+01]
  QLMatrix range   [1e+00, 4e+04]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e-05, 4e+06]
  QRHS range       [2e+04, 4e+04]
Presolve removed 20431 rows and 5836 columns
Presolve time: 0.53s
Presolved: 22847 rows, 7827 columns, 63853 nonzeros
Presolved model has 209 quadratic constraint(s)
Presolved model has 3528 bilinear constraint(s)
Variable types: 5865 continuous, 1962 integer (1936 binary)

Root relaxation: objective -2.666593e+06, 5311 iterations, 0.47 seconds (0.32 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -2666593.1    0 2064          - -2666593.1      -     -    2s
     0     0 -2639548.4    0 2126          - -2639548.4      -     -    5s
     0     0 -2639125.8    0 2096          - -2639125.8      -     -    5s
     0     0 -2629372.1    0 2302          - -2629372.1      -     -    7s
     0     0 -2628719.6    0 2316          - -2628719.6      -     -    7s
     0     0 -2620742.0    0 2633          - -2620742.0      -     -    9s
     0     0 -2620562.2    0 2708          - -2620562.2      -     -    9s
     0     0 -2613245.4    0 2388          - -2613245.4      -     -   10s
     0     0 -2613119.9    0 2501          - -2613119.9      -     -   11s

I would like to ask if there is a mismatch between yalmip and the gurobi version.

[johanlofberg]

You haven't showed finished logs from computers, so we don't know what differs

Also, you haven't confirmed that you are running the same yalmip version on both (installed from scratch so you haven't edited something in some of them)

[Electricman23]

Hello, do you mean can I understand that the yalmip version will have an effect on the result of running the program, or that a low version of yalmip will cause incorrect results in the model solved with the gurobi solver.
As for the yalmip version you said, I am looking for a method inquiry and will reply to you later.
Below I will give the full log of both computers after running.
First, the complete log of the first computer after running the program using matalbR2019b called gurobi9.5.1:

Set parameter IntFeasTol to value 1e-06
Set parameter MIPGap to value 0.3
Set parameter ScaleFlag to value 1
Set parameter BarIterLimit to value 2000000000
Set parameter CrossoverBasis to value 0
Set parameter MinRelNodes to value 0
Set parameter NodeMethod to value 1
Set parameter PumpPasses to value 0
Set parameter ZeroObjNodes to value 0
Set parameter AggFill to value 10
Set parameter NonConvex to value 2
Set parameter PreSparsify to value 0
Set parameter TuneTrials to value 2
Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)
Thread count: 4 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 29449 rows, 18322 columns and 128677 nonzeros
Model fingerprint: 0xfbaa4730
Variable types: 12905 continuous, 5417 integer (5417 binary)
Coefficient statistics:
  Matrix range     [1e-05, 4e+06]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e-05, 4e+06]
Presolve removed 24049 rows and 16197 columns
Presolve time: 0.43s
Presolved: 5400 rows, 2125 columns, 19895 nonzeros
Variable types: 1082 continuous, 1043 integer (1042 binary)

Root relaxation: objective -2.394388e+08, 1801 iterations, 0.07 seconds (0.03 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -2.394e+08    0   34          - -2.394e+08      -     -    0s
     0     0 -2.394e+08    0  163          - -2.394e+08      -     -    0s
     0     0 -2.394e+08    0  132          - -2.394e+08      -     -    1s
H    0     0                    -2.39439e+08 -2.394e+08  0.00%     -    1s
     0     0 -2.394e+08    0   20 -2.394e+08 -2.394e+08  0.00%     -    1s

Cutting planes:
  Learned: 51
  Gomory: 29
  Cover: 40
  Implied bound: 11
  Clique: 28
  MIR: 66
  StrongCG: 4
  Flow cover: 68
  RLT: 22
  Relax-and-lift: 8
  BQP: 2

Explored 1 nodes (7754 simplex iterations) in 1.33 seconds (0.58 work units)
Thread count was 4 (of 4 available processors)

Solution count 1: -2.39439e+08 
No other solutions better than -2.39439e+08

Optimal solution found (tolerance 3.00e-01)
Warning: max constraint violation (1.6364e-06) exceeds tolerance
Best objective -2.394387793003e+08, best bound -2.394387793003e+08, gap 0.0000%

As for the logs where I behaved in the recovery that another computer used matlabr2018b to call gurobi 5.3.1 to get the correct result, I need to ask the owner of the other computer if you need it and send it to you later some time.

[johanlofberg]

You have now shown two logs where it converges to -2.394e8,, and one where you have no solution in the log but a lower bound in the iteration which is consistent, so not sure what you have shown

You cannot use YALMIP version A with gurobi version B on MATLAB version C on computer D and compare with YALMIP version E with gurobi version F on MATLAB version G on computer H, and try to understand what is causing a difference. THere coulöd be a bug in any of the codes, or just an artifact from the fact that you are solving a nonconvex problem, so an epsilon difference anywhere can cause a completely different solution path and cause failure or trigger bugs etc

You can export the YALMIP model on the two setups and compare these (and thus conclude if YALMIP compilation differs, and/or use that exported model to call gurobi on the other machine to see if the same model causes different solutions)

[Electricman23]

Hello, there is also a problem that although the first program can produce results after running, it is obvious that the conclusion is wrong. For example, the decision variable f1, I set its upper limit to 10e+6, and all f1 after the program runs is 10e+6.
So, is this a model problem or something else. What should I do about YALMIP model export?

[johanlofberg]

https://yalmip.github.io/command/export/

[Electricman23]

Hello, I don't know if it has anything to do with the solver, when I use the computer that can run the results, when I am in the workspace yalmiptest, I find that there is GUROBI NONCONVEX this solver.

|        GUROBI|           GUROBI|       found|
|        GUROBI|        NONCONVEX|       found|
|        GUROBI|              MEX|   not found|

[johanlofberg]

not sure what you say, but this is what it is supposed to look like

[Electricman23]

I want to make it clear that the computer that can't run the results doesn't have this solver.
And then I want to ask about the role of this solver, GUROBI NONCONVEX。
Thank you！

[johanlofberg]

it means it is running a very old version of yalmip

[Electricman23]

If I run the results in the old version, but I can't run the results in the latest version of yalmip, does that mean that my model is still problematic?

[johanlofberg]

if the only thing that differs is the YALMIP version, and you get the wrong solution in the latest version but correct in an older, then obviously something is wrong in the latest version of YALMIP, and you would have to reveal which old version you are using

[Electricman23]

Very sorry, I may not have made it clear.
First, the working version is Yalmip'20210331'+Gurobi9.5.1.
Secondly, I replaced the version of gurobi to 10.0.1 and it failed.
At the same time, my initial version of yalmip did not know what version it was, but running the program also failed.

[johanlofberg]

The my guess is that gurobi simply is buggy. Most common bugs in gurobi is presolve errors, hence turn off presolve in gutobi

[Electricman23]

What I know about preslove is the role of algorithmic simplification models, so what you mean is that I simplified my model when using other versions of gurobi, causing my model to not run results.
One more thing, please ask how to turn off Preslove.

[johanlofberg]

What I am saying is that many solvers have bugs in their presolve routine, so it is always the first time I tweak when trying to understand why a solver computes the wrong thing. Set gurobi.presolve to 0

[Electricman23]

Okay, thank you very much for taking the time to help solve the problem, I will first modify it according to your suggestions, and if there are still problems in the future, I hope you will advise.

[johanlofberg]

After testing various combinations

Gurobi 9.5, YALMIP Develop branch, Your options: Returns solution objective -2.56 * 10^6
Gurobi 9.5, YALMIP Develop branch, No options: Returns solution objective -2.56 * 10^6
Gurobi 10.01, YALMIP Develop branch, Your options: claims infeasible at point where lower bound is -2.6 * 10^6
Gurobi 10.01, YALMIP Develop branch, Presolve turned off: no solution found after 1000s, lower bound around -2.6 * 10^6
Gurobi 10.01, YALMIP Develop branch, Cuts turned off: solution with objective -2.55 * 10^6 found

Hence, my conclussion is that 9.5 returns a correct solution with objective -2.56 * 10^6, while 10.01 has a bug in some of the cut generators causing to cut off solution space incorrectly

The logs you show above where objective is -2.4 * 10^8 I cannot reproduce and I have no idea how you got them.

[wnn666666]

When I met the same problem, how did you finally solve it? I also had a situation where the decision variable was the upper limit，My yalmip is 2021.03.31 and Gurobi is 10.0.3. Would you like to correct one of these versions.