Weird results with 3-dim binary matrix

MIDCP
1

Hello there,
I am trying to run the code below. Yet it produces unexpected results. Under all constraints and the data given we would expect exactly one item of u(i,:,: ) to be set to 1 . Yet all u’s are always set to 0.
For information - one item u(i,h1,: ) can be set to 1 if and only if y(i,h1) has been set to 1, and only one index of u can be set to 1 in the chosen row.

Could you think of a possible reason? I am trying to find a bug in implementation yet I totally fail. All seems good but the results are bad and instead of getting some u’s=1 I am always getting zeros. What is more, I have tried to solve it with both Gurobi and Mosek, always receiving the same result no matter the solver.

------- My code: ----

clear all
n_hnon = 3; 
n_href = 2; 
n_h = n_hnon + n_href; 
n_p = 15;
s_href = 0.8 + (0.93-0.8).* rand(n_href,1); 
s_hnon = 0.1 + (0.3-0.1).*rand(n_hnon,1); 
s_hs=[s_hnon;s_href];
H_max_non=100*ones(n_hnon,1);
H_max_ref=100*ones(n_href,1);
eta=0.2;
%% Optimization
cvx_begin
variable s_p2(n_p)
variable s_term1(n_p)
variable s_term2(n_p)
variable s_term3(n_p)
variable s_term4(n_p)
variable w(n_p,n_hnon)
variable b(n_p,n_href) binary
variable y(n_p,n_hnon) binary
variable u(n_hnon,n_href,n_p) binary

maximize sum(s_p2)
subject to
for p1=1:n_p
    s_term1(p1) == eta*sum(s_hnon'.*y(p1,:))
end
for p2=1:n_p
   s_term2(p2) == (1-eta)*sum(s_hnon'.*w(p2,:))
end
for p3=1:n_p
    s_term3(p3) == sum(s_href'.*b(p3,:))
end
for p4=1:n_p
    for h1 = 1:n_hnon
        s_term4(p4) == sum((100.*squeeze(u(h1,:,p4))))
    end    
end
s_p2==s_term1+s_term2+s_term3+s_term4
%Product of binaries linearisation
w<=y
for p5=1:n_p
    for h11=1:n_hnon
        w(p5,h11)<=1-sum(u(h11,:,p5))
        w(p5,h11)>=y(p5,h11)+(1-sum(u(h11,:,p5)))-1
    end
end
for h12 = 1:n_hnon
     sum(y(:,h12))<=H_max_non(h12)
end
for h2 = 1:n_href
    sum(b(:,h2))+sum(sum(u(:,h2,:)))<=H_max_ref(h2)
end
%choose one
for p7=1:n_p
    for h13=1:n_hnon
          sum(u(h13,:,p7)) - y(p7,h13) <= 0           
    end
end

for p9=1:n_p
    sum(y(p9,:))+sum(b(p9,:))==1
end

b==0

cvx_end
2

The problem is infeasible if the constraint sum(u(:)) >= 1 is added, so of course the optimal solution has all elements of u equal 0. However, if that constraint is added and the 'choose one" double for loop constraints are removed, the problem is feasible, the solution has all elements of u equal 1, and the optimal objective value is far greater (better) than the original problem.

I will leave further diagnosis as to what is 'wrong" with your model to you, since I have no idea what it is supposed to model or do. You may find https://yalmip.github.io/debugginginfeasible , all of which except section 1 also applies to CVX, to be helpful.

3

Thanks Mark. I also found that the problem is linked to the double “choose one” loop. Yet, it was supposed to be a fairly simple constraint in the form: sum(u(i,:,p)) <= y(p,i) \forall p an i
Which should not cause any issues. There should exist a solution where exactly one element of u in dim1,dim2,dim3 is set to 1. Basically I am trying to model the situation where exactly one column element of u at index i,…,p can be set to 1 if and only if y(p,i) has been set to 1 . If y(p,i) is 0 all elements of u(i,:,p) should be 0.

Frankly I fail to understand why this constr causes infeasibiloty. There is no other contradicting contr (or at least I am not seeing it,) Do you have any tips why this is not OK and where to dig further?

Maybe I am making a mistake in indexing of the 3-dim array?