Error in 0-1 decision variable when solving model and the multiplication and division problem of the two decision variables.

[Electricman23]

When I compare the inequality of 0-1 decision variables, the result shows that it is 0.5. If there is a decision variable v(i,m) and the constraint v(i,m+1)>=v(i,m) is set, the ideal goal is not achieved.
The expected possibility is
step m step (m+1)
0 0
0 1
1 1
I assigned a value to V in the first step. E.g. v(1:5,1)=1,v(6:27,m)=0.

[johanlofberg]

You have to use code-blocks <> when you paste code into posts. As it is now, it cannot be copied to run as it is broken.

There are no official commands in yalmip called solvesdp or double https://yalmip.github.io/tutorial/basics/

It is unclear what your actual question is now.

[johanlofberg]

and code gets messy to read when you use & to separate constraints. in particular when you mix using & and ,

[Electricman23]

You have to use code-blocks <> when you paste code into posts. As it is now, it cannot be copied to run as it is broken.

There are no official commands in yalmip called solvesdp or double https://yalmip.github.io/tutorial/basics/

It is unclear what your actual question is now.

First of all, thank you for your reply!
I found that the running time was too long when I ran the program. Therefore, I reduced the time step to two time steps, knowing that the first time step was the value set for me, and after running, I found that the K(1,m) and H(1,m) of the second time step could not get the result.
I have tried two methods, the first is to multiply the denominator B in the right division to the left; the second is to add an auxiliary variable C so that BC = 1, and then substitute C into the original form.
My latest modified program is shown below and would appreciate your help.

[johanlofberg]

No, you haven't placed the code in a code block <> as it then would look like

x == 5*y
y == z*x

while yours is simply copy-pasted leading to the broken code

x == 5y
y == zx

[johanlofberg]

...and you are still using solvesdp and double

[Electricman23]

Sorry, because it is the first time to contact this community platform and the relationship between beginners, it may produce a lot of small mistakes. Double appears in my program because I modify it according to someone else's program to find the value in the decision variable. As for uploading code programming, I will systematically learn from this community platform to avoid duplication of effort, thank you.

[johanlofberg]

then use value/optimize as shown on the tutorial link

[Electricman23]

Dear johanlofberg! After studying, I make changes to the program according to your guidance.Hope to get your guidance.
As for the question, I would like to make the following statement:
I found that the running time was too long when I ran the program. Therefore, I reduced the time step to two time steps, knowing that the first time step was the value set for me, and after running, I found that the K(1,m) and H(1,m) of the second time step could not get the result.
I have tried two methods, the first is to multiply the denominator B in the right division to the left; the second is to add an auxiliary variable C so that BC = 1, and then substitute C into the original form.
The current procedure adopts the first approach, but fails to meet the desired requirements.

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.
R=0.5/50;%频率变化量最大值  0.5HZ/s
SG=2500000;%外部带电域或储能的额定容量  2.5MW
KG=50;%外部带电域或储能的功频静特性系数  50p.u.
HG=5;%外部带电域或储能的惯性时间常数  5s
PL=18000;
SW=2000000;
KW=25;
KL=1;
HW=10;

%% 决策变量
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');%每一时步的时间
n=sdpvar(NC,NT,'full');
k=sdpvar(NC,NT,'full');
v=binvar(NC,NT,'full');%每一个时步每一组新能源机组辅机的运行状态
u=binvar(NC,NT,'full');%每一个时步每一组新能源机组的运行状态   0表示未稳定出力  1表示稳定出力
a=binvar(13,NT,'full');
b1=binvar(3,NT,'full');
b2=binvar(3,NT,'full');
b3=binvar(3,NT,'full');
b4=binvar(3,NT,'full');
b5=binvar(3,NT,'full');
S1=sdpvar(1,NT,'full');
S2=sdpvar(1,NT,'full');
S3=sdpvar(1,NT,'full');


%% 目标函数
st=[];

for m=1:NT
                f=-sum(P2(1,m)*T(1,m));
end


%% 约束条件
%根据约束对第一时步相关参数进行预设
for m=1
        st=st+[v(1,m)==1;
               v(2:5,m)==0;
               u(1:5,m)==0;
               H(1,m)==HG;
               K(1,m)==KG;
               P1(1,m)==90000;
               P2(1,m)==0;
%                P1(1,m)<=fmax*K(1,m)*SG;
%                P1(1,m)<=2*R*H(1,m)*SG;
               ];
end

% for m=1
%     for i=1:ND
%         st=st+[S1(1,m)*(sum(PL(i,1)*v(i,m))+SG)==1;
%                K(1,m)==(sum(KL(i,1)*PL(i,1)*v(i,m))+KG*SG)*S1(1,m);
%                ];
%     end
% end

%0/1决策变量约束，辅机的状态量与新能源机组的状态量间的约束

for m=2
     for i=1:NC
        st=st+[v(i,m-1)==v(i,m)];
     end
end

% for m=3
%      for i=1:ND
%         st=st+[implies(v(i,m-1)==1,v(i,m)==1)];%加入这一步后决策变量变为0.5(与矩阵的乘除法有关系)
%      end
% end
% 
% for m=4
%      for i=1:ND
%         st=st+[v(i,m-1)==v(i,m)];%加入这一步后决策变量变为0.5(与矩阵的乘除法有关系)
%      end
% end

for m=2:NT
     for i=1:NC
        st=st+[u(i,m)==v(i,m-1)];%辅机启动后下一时步新能源机组启动
     end
end

%对每一时步的01决策变量进行计数
for m=1:NT
      for i=1:NC
            st=st+[n(i,m)==v(i,m)*NE(i,1);
                   k(i,m)==u(i,m)*NE(i,1)];
      end
end

% 时间约束    新能源辅机恢复时，新能源所在线路及变压器恢复考虑在这一时步（40+20+5*n）新能源机组恢复时（40+10）
for m=1
    for i=1:NC
            st=st+[T(1,m)==max(60+5*n(i,m))];
    end
end


% 稳态频率约束
for m=2:NT
    for i=1:NC
        st=st+[
               (SW*sum(k(i,m))+PL*sum(n(i,m))+SG)*K(1,m)==KW*SW*sum(k(i,m))+KL*PL*sum(n(i,m))+KG*SG;
               P1(1,m)<=fmax*K(1,m)*(SW*sum(k(i,m))+SG);
               P2(1,m)<=fmax*K(1,m)*SW*sum(k(i,m))];
    end
end

% 惯量响应约束
for m=2:NT
    for i=1:NC
        st=st+[(SW*sum(k(i,m))+SG)*H(1,m)==HW*SW*sum(k(i,m))+HG*SG;
               P1(1,m)<=2*R*H(1,m)*(SW*sum(k(i,m))+SG);
               P2(1,m)<=2*R*H(1,m)*SW*sum(k(i,m))];
    end
end

% 启动功率以及新能源机组出力等式关系
for m=2:NT
    for i=1:NC
        st=st+[P1(1,m)==PL*sum(n(i,m));
               0.1*SW*sum(k(i,m))<=P2(1,m)<=SW*sum(k(i,m))];
    end
end

%% 求解
gurobi_stop_gap = 10^(-4); % 设置求解混合整数规划的 Gap 为10^(-4)
longest_solve_time = 300; % 设置最长求解时间为 300 秒
NumericFocus = 2;%
options = sdpsettings('verbose',2,'solver','GUROBI', 'debug', 1, 'gurobi.NonConvex', 2,'gurobi.timelimit',longest_solve_time,...
    'gurobi.MIPGap', gurobi_stop_gap);%cplex%,'gurobi.NumericFocus',NumericFocus
result=optimize(st,f);
value(f)

%% 结论
AP1=value(P1);
AP2=value(P2);
Bv=value(v);
Bu=value(u);
BK=value(K);
BH=value(H);
CT=value(T);
Cn=value(n);
Ck=value(k);

[Electricman23]

Thankfully, I might have been programming the wrong way without your reminder, and I'd like to see if gurobi can solve the objective function as a value multiplied by two decision variables in a constraint. At the same time, I would also like to know how to view the diagnostic reports like the one you showed, because I can't find it on the Internet and want you to teach.
Second, I hope you'll help see if these things that appear after the modifications indicate that although my program is not necessarily correct, the solver used is gurobi and not something else.

Set parameter TimeLimit to value 300
Set parameter IntFeasTol to value 1e-06
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.2 build v9.5.2rc0 (win64)
Thread count: 4 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 119 rows, 81 columns and 360 nonzeros
Model fingerprint: 0x74d9283f
Model has 1 quadratic objective term
Variable types: 39 continuous, 42 integer (42 binary)
Coefficient statistics:
  Matrix range     [1e+00, 2e+06]
  Objective range  [0e+00, 0e+00]
  QObjective range [2e+00, 2e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 9e+04]
Presolve removed 119 rows and 79 columns
Presolve time: 0.00s
Presolved: 2 rows, 4 columns, 4 nonzeros
Presolved model has 1 bilinear constraint(s)
Variable types: 4 continuous, 0 integer (0 binary)

Root relaxation: unbounded, 1 iterations, 0.00 seconds (0.00 work units)

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

     0     0  postponed    0               -          -      -     -    0s
     0     0  postponed    0               -          -      -     -    0s
H    0     0                    -9.37129e+14          -      -     -    0s
     0     2  postponed    0      -9.371e+14          -      -     -    0s

Explored 3 nodes (5 simplex iterations) in 0.13 seconds (0.00 work units)
Thread count was 4 (of 4 available processors)

Solution count 1: -9.37129e+14 
No other solutions better than -9.37129e+14

Model is unbounded
Best objective -9.371290370337e+14, best bound -9.371290370337e+14, gap 0.0000%

ans =

     0

[johanlofberg]

It's the variable you call result. You don't display it as you have a ; after the call to optimize (but gurobi is writing it clearly in the log too)

[Electricman23]

I already know. But when I have a diagnostic report like you showed earlier, what should I do to make my model viable. According to the page, should I hide the objective function to see if it is feasible. But my goal is to multiply the two decision variables in the constraint, whether this is one of the reasons why my goal is unbounded, I hope you can help explain.

[johanlofberg]

If the problem is infeasible, you have to find what makes it infeasible. Once you fixed that, you add the objective and if it turns out to be unbounded, you have to figure out why it is unbounded. There are strategies described for both debugging adventures.

[Electricman23]

Okay, I'll follow the tutorial in https://yalmip.github.io/debuggingunbounded to check the errors, and I'll ask you more questions later.

[Electricman23]

Dear Johanlofberg, I read the website tutorial carefully and wanted to see the feasible problem by hiding the objective function. But there are very few tutorials on the web, so I failed to achieve the goal of hiding the objective function. I hereby ask you how I should do it, whether I need to modify the relevant content of the solver, if so, please point to it in the solution module below me.

gurobi_stop_gap = 10^(-4); 
longest_solve_time = 300; 
NumericFocus = 2;%
options=sdpsettings('verbose',2,'solver','gurobi+', 'debug', 1, 'gurobi.NonConvex', 2,'gurobi.timelimit',longest_solve_time,...
    'gurobi.MIPGap', gurobi_stop_gap);%cplex%,'gurobi.NumericFocus',NumericFocus
result=optimize(st,f,options)
value(f);

[johanlofberg]

Well, now you used bnb instead again as you didn't give the options structure, but I'm not sure what you are saying or asking at all. It is feasible (as I already told you above), and thus unboundedness is the issue to fix

[Electricman23]

Okay, does this mean that my constraints are fine, the problem is that my objective function has no boundaries. But doesn't the gurubi solver solve the nonlinear problem of multiplying two continuous decision variables? I'll follow the site tutorial below.

[johanlofberg]

The latest model you posted is obviously feasible. Gurobi supports quadratic models.

[Electricman23]

I want to explain why I didn't give the option structure for constrained validation, because when I add options, it will report an error, the error is as follows:

错误使用 solvesdp (line 89)
Second argument (the objective function h) should be an sdpvar or logdet object (or empty). It appears as if you sent an options structure in
the second argument.

出错 optimize (line 31)
[varargout{1:nargout}] = solvesdp(varargin{:});

出错 CZYH0810LY (line 294)
result=optimize(st,options)

[johanlofberg]

options goes in the third argument, you have called optimize with it as the second argument

[Electricman23]

Dear Johanlofberg, Hello! After yesterday's exchange, I have a preliminary grasp of the method of examining constraints to determine whether the model is feasible. So I added a new module to the previous model, but unfortunately, the model started a new non-feasible problem. I would like to ask you to help me see if there is a logical problem with my new module.
This is the constraint description, when the step m is odd, there are the following constraints for v of the current time step:
1.v(i,m-1)=1→ v(i,m)=1
that is, when the 0-1 decision variable for the m-1 step is 1, and the mth step is also 1.
2.v(i,m-1)=0→ v(i,m)=0/1
i.e. when step v is 0 at m-1, v of step m is judged by the following constraints as to whether it is 0 or 1.

 b1=binvar(3,1);
 for m=3
        for i=1:NC
        st=st+[
                    sum(b1)==1,
                    implies(b1(1),[v(i,m-1)==0,v(i,m)==0]);
                    implies(b1(2),[v(i,m-1)==0,v(i,m)==1]);
                    implies(b1(3),[v(i,m-1)==1,v(i,m)==1]);
                  ];
      end
end

[Electricman23]

I'm sorry, I thought too much, and added some meaningless constraints to it, thanks for the teaching. At the same time, I would like to ask you what you just said about the option in the third argument, then how I should modify it in optimize so that the program avoids errors after hiding the objective function f.

gurobi_stop_gap = 10^(-4); % 设置求解混合整数规划的 Gap 为10^(-4)
longest_solve_time = 300; % 设置最长求解时间为 300 秒
NumericFocus = 2;%
options=sdpsettings('verbose',2,'solver','gurobi+', 'debug', 1, 'gurobi.NonConvex', 2,'gurobi.timelimit',longest_solve_time,...
    'gurobi.MIPGap', gurobi_stop_gap);%cplex%,'gurobi.NumericFocus',NumericFocus
result=optimize(st,f,options)

[johanlofberg]

use [] or 0 as objective

[Electricman23]

Okay, thanks, I have another question that has been confused for a long time and would like to ask you. When I hide the objective function to determine the feasibility of the model, the result obtained after running it can not be regarded as a correct result, but only proves the feasibility of my model, and only when I add the objective function will I optimize to get the final result.

[johanlofberg]

Of course. To get the optimal solution you must solve the actual problem you are interested in. However, your model is flawed as it is unbounded. Removing the objective is only something you do when you have feasibility problems and want to figure out what causes infeasibility (your initial model you posted was infeasible, the model you later posted was unbounded, two completely different issues you are dealing with)

[Electricman23]

I can understand that you mean that my original model was not feasible, but after I modified it, the current model is already feasible, but when I added the objective function, my model created a new problem, that is, the goal is infinite. Therefore, the next thing I should do is to modify the target according to the tutorial above the website.