Troubles solving MIP using Gurobi

[Astghik]

Hello,

I'm trying to simulate scenario-based MPC problem using Mixed Integer Constraints. However, when I run the code, at some point Gurobi won't find a solution. And in case it finds a solution, the behavior of my object becomes really weird. The code is quite big, so I will just post the part of defining the objective and constraints.

Here I have 2 uncertain variables, w1 and w2 and N scenarios generated for each of them, K is my prediction horizon.

x = sdpvar(repmat(nx,1,K+1),ones(1,K+1));

y1 = sdpvar(repmat(ny,1,K),repmat(N,1,K));

y2 = sdpvar(repmat(ny,1,K),repmat(N,1,K));

r = sdpvar(repmat(ny,1,K+1),ones(1,K+1));

obs1=sdpvar(ny,2);

obs2=sdpvar(ny,2);

delta=0.04;

alpha=0.95;

cons=[];

objective = 0;

M=60;

za1=sdpvar(N,K); 

z1 = sdpvar(1,K);

za2=sdpvar(N,K); 

z2 = sdpvar(1,K);

dk1 = binvar(repmat(2*ny,1,K),repmat(N,1,K));

dk2 = binvar(repmat(2*ny,1,K),repmat(N,1,K));

e=0.01;

for k = 1:K

    for i=1:N

        cons=cons+[sum((C*x{k}-y1{k}(:,i)-sum(w1{i}(:,1:k),2)).^2)-z1(k)-za1(i,k)<=0];

        cons=cons+[sum((C*x{k}-y2{k}(:,i)-sum(w2{i}(:,1:k),2)).^2)-z2(k)-za2(i,k)<=0];

    end

    cons=cons+[y1{k}(1,:)-obs1(1,1)+e<=M*dk1{k}(1,:),...

               obs1(1,2)-y1{k}(1,:)+e<=M*dk1{k}(2,:),...

               y1{k}(2,:)-obs1(2,1)+e<=M*dk1{k}(3,:),...

               obs1(2,2)-y1{k}(2,:)+e<=M*dk1{k}(4,:),...

               y1{k}(3,:)-obs1(3,1)+e<=M*dk1{k}(5,:),...

               obs1(3,2)-y1{k}(3,:)+e<=M*dk1{k}(6,:),...

               sum(dk1{k})<=5];

    cons=cons+[y2{k}(1,:)-obs2(1,1)+e<=M*dk2{k}(1,:),...

               obs2(1,2)-y2{k}(1,:)+e<=M*dk2{k}(2,:),...

               y2{k}(2,:)-obs2(2,1)+e<=M*dk2{k}(3,:),...

               obs2(2,2)-y2{k}(2,:)+e<=M*dk2{k}(4,:),...

               y2{k}(3,:)-obs2(3,1)+e<=M*dk2{k}(5,:),...

               obs2(3,2)-y2{k}(3,:)+e<=M*dk2{k}(6,:),...

               sum(dk2{k})<=5];

    objective = objective + (C*x{k}-r{k})'*(C*x{k}-r{k})+u{k}'*R*u{k};

    cons = cons+[x{k+1} == A*x{k}+B*u{k},...

               1/N*(z1(k)+za1(:,k)/(1-alpha))<=delta, za1(:,k)>=0,  1/N*(z2(k)+za2(:,k)/(1-alpha))<=delta, za2(:,k)>=0];

end

 

objective = objective + (C*x{K+1}-r{K+1})'*(C*x{K+1}-r{K+1});

parameters_in = {x{1},[r{:}],obs1,obs2};

solutions_out = {[u{1}]};

opt=sdpsettings('solver','gurobi');

opt.showprogress=1;

opt.debug=1;

opt.verbose=2;

controller = optimizer(cons,objective,opt,parameters_in,solutions_out)

[Johan Löfberg]

impossible to respond with code that runs

code looks reasonable

[Astghik]

Thank you for the prompt response!

I've tried to make a reproducible code. You just need to run the main.m file.

[Johan Löfberg]

looks reasonable.

you would have to define "weird" etc. typical mistake is that you have a dynamical system and plan with a too short horizon/too aggressively, leading to the vehicle running into walls in the predictions etc leading to infeasibility. I note you don't have any Riccati-based terminal penalty, and a short horizon, which might easily lead to such short-sighted behavior. Big-m constant not properly selected is another typical problem

[Astghik]

Please see attached the video of complete simulation, which is what I call weird. My agent avoids the obstacles but doesn't follow the reference path after it.

Also, I can't use bigger prediction horizons as Gurobi doesn't find a solution for it as shown in the attached screenshot.

[Johan Löfberg]

you talked about infeasibility, but then how can you even finish the simulation ?

you should use a relevant terminal penalty to balance the too short horizon

does it work with no obstacles. 1 obstacle. longer horizon.  long horizon and slacks. other solver

[Astghik]

The part I can't understand is that sometimes the solver returns infeasible, sometimes it works till the end, and sometimes it just gets stuck at some point.

In case of no obstacles, it works fine, and even in the case of 1 obstacle, it works. Before using Gurobi I used Cplex, but it was even worse. As I understood, you think that the problem is associated with the tuning parameters?

I just added the Riccati-based terminal penalty as per your advice. But the simulations take quite a long time to complete. Will update you as soon as I get a result.

Thanks!!

[Johan Löfberg]

since you have random stuff, isn't it rather reasonable that it might become infeasible for some setups

[Astghik]

I don't see any reason why it should become infeasible for some setups, because obstacle avoidance constraints are given in such a way, that they can be violated to some extent.

After adding the terminal penalty and increasing the prediction horizon, the controller works a lot better, but it became much slower and at some point becomes infeasible.

[Johan Löfberg]

recursive feasibility. with a too short horizon you put too much momentum into the vehicle, and then boom you don't have the possibility to break/steer to avoid an obstacle

[Astghik]

Professor Lofberg,

Thank you for your advice! I increased the horizon, as well as added a terminal cost as you recommended. However, the issue I got was because of a small mistake in the coding.

But still, it takes about 2-3 hours to complete one simulation. Is there a way to speed up the process?

Thank you!

[Johan Löfberg]

There is no simple trick. A clever reformulation might render it trivial to solve. Solving MISOCPs is typically much worse than solving MILP or MIQP. Tweaking options can make huge differences. Isn't sum(dk1{k})==5 a sufficient model (forcing to 1 doesn't change anything as it becomes redundant, but ==5 instead of <=5 is combinatorially much smaller as it has 6 possibilities vs 63)

Note that your code is a bit dangerous as you will get symmetric obs1 etc if ny=2. You should use 'full' flag for safety. Same with dk which will be symmetric if 2*ny = N

[Astghik]

So in my case, I have an obstacle with obs1=[x_min x_max; y_min y_max; z_min z_max] coordinates. I want my y1 to be outside of the obs1, that's why I put sum(dk1{k})<=5.

[Johan Löfberg]

Such a disjunction should say that you must violate 1 out of 6 constraints (since a box has 6 sides), so it should be possible to represent with a model which says sum(violation) == 1

implies(violation1, ontop) + implies(violation2, onbottom) + implies(violation3, totheleft) + ...+sum(violation == 1)

[Astghik]

I want "at least" one constraint to be violated. Then, in that case, it will be sum(violation)>=1, which is the same as sum(dk1{k})<=5.

[Johan Löfberg]

As I said, you have to think about the combinatorial complexity of how you model things

f1(x) >=0 or f2(x)>= 0 is most efficiently represented as 

f1(x) >= M*d1

f2(x) >= M*d2

d1 + d2 == 1

With that model, you will definitely have that at-least 1 of them is violated. It could be that both are violated (i.e. one is violated for free so to speak). The number of combinatorial cases in this model is 2. Had you used the model d1 + d2 >= 1, you would have had 3 combinatorial cases, and thus a worse model. For n disjunctions, you get n combinatorial cases, whereas a sum(d)>=1 would lead to 2^n-1 possibilities

[frankshi]

Hi Professor, I am also studying a MISOCP problem with large-scale constraints and variables. The modeling process (adding the constraints) takes a very long time with about several hours. Is it possible to use other modeling languages such as GAMS to get a fast modeling process? Thank you very much.

[Johan Löfberg]

Your question is weird. You are asking on the YALMIP forum if it is possible to use GAMS. It's like asking on the Toyota forum if it is possible to drive a Ford. Of course, but you will not get any help on that forum.

Slow modelling can typically be fixed by vectorizing poorly vectorized code. In SOCP case often even more improvement by using the low-level cone operator

[frankshi]

Hi Professor,

Thanks. My question is actually how to speed up the modeling process. I am currently using the norm command to model the SOCP constraints. I will study the cone operator to check the improvement. Thanks very much.