Handling of concave in modelling.

[Rupaj Kumar Nayak]

Whether Yalmip handles the following?

x(i), x(j), x(k) are binary variables.

X(i,j) ,X(i,k), X(j,k) are variables.

constraint is X(i,j)-X(i,k)<= min(x(i),x(j))- min(x(i),x(k)).

It is to be noted that R.H.S. is concave and L.H.S is affine and CVX is unable to handle it.

[Johan Löfberg]

Yes, you can write it as exactly that. YALMIP will recognize the nonconvexity of the constraint, and switch to a MILP representation of the second min

Note that this model requires all involved variables to explicitly bounded in your model, i.e., you must have explicit upper and lower bounds on x(i), x(j) and x(k) (which they are by definition since they are binary, so most likely you do not have to add it. Try it and you will see. If bounds are missing, YALMIP will scream about lousy bounds in big-M model)

[Johan Löfberg]

Since the involved variables are binary in the complicating inner min, the exact model is simple to derive. Indeed, YALMIP will recognize and exploit the structure.

YALMIP will replace min(x(i),x(k)) with a new variable z(i,k) which satisfies

z(i,k) <= x(i)
z(i,k) <= x(k)
z(i,k) <= x(i)+x(k)-1

[Rupaj Kumar Nayak]

Many Thanks Johan.

I am new to Yalmip and like to be clarified on the following.

I am solving the problem as follows:

min x^TA0x+b0^Tx
%         s.t.  x^TAix+bi^Tx+ci <= 0, i=1,...,nc, convex, i=1:nnc nonconvex
%                a^Tx+b <= 0,
%               x\in { 0, 1}^n integer

I was using CVX which requires convex formulation of constraint, whereas my constraint is non-convex in nature.

I am solving the above QCQP by SDP relaxation.

Would you like to guide me how I can use Yalmip to solve the above QCQP and solve by SDP relaxation by calling Sedumi.

[Johan Löfberg]

Are you solving the SDP relaxation by addding quadratic constraints of the type x*(x-1)=0? In the end, I will tell you how to solve the problem, an SDP relaxation is most likely a poor choice.

Define the model (depends on how you store data of course)

x = sdpvar(n,1);
objective = x'*A0*x + b0^T*x;
Constraints = a'*x + b <= 0;
for i = 1:nc+nnc
Constraints = [Constraints, x'*A{i}*x + bi{i}^T*x + c(i) <= 0];
end


This model is used for any solver.

YALMIP has a built-in moment relaxation framework (i.e., semidefinite relaxations). The standard first-order relaxation, which probably is what you refer to, is solved by default
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.solvemoment

solvemoment([Constraints, x.*(x-1) == 0],objective)

Stronger relaxations can be derived and computed too

solvemoment([Constraints, x.(x-1) == 0],objective,[],2)

Alternatively, use standard solver call, but specify moment as solver

solvesdp([Constraints, x.*(x-1) == 0],objective,sdpsettings('moment.order',2))

In this format, you can switch to the semidefinite relaxation framework sparsepop if you want (solves sparse polynomial relaxations)

solvesdp([Constraints, x.*(x-1) == 0],objective,sdpsettings('solver','sparsepop'))

The first-order relaxation can be encoded semi-manually, by telling YALMIP to relax all monomials to linear variables. With that, a semidefinite relaxation is equivalent to

solvesdp([Constraints, x.*(x-1) == 0,[1;x]*[1 x'] >= 0],objective,sdpsettings('relax',1))

However, I would not use semidefinite relaxations for this problem. The problem can be solved using a MILP solver or a convex MISOCP solver. A MILP solver can be used by recognizing that a bilinear term x(i)*x(j) can be handled by introducing a new variable w(i,j) which models x(i) AND x(j), which can be done using three linear constraints. This works in low dimensions, or if the quadratics are sparse so you don't have to introduce that many new variables.

Alternatively, exploit the fact that a nonconvex quadratic function in binary variables always can be rewritten as a convex function, by exploiting the fact that x(i)^2 = x(i). Hence, x'*Q*x = x'*(Q+r*I)*x-r*sum(x) for any scalar r. Pick r large enough so that matrix is positive semidefinite. Convex MISOCPs are solved reasonably effciently by solvers gurobi, cplex and mosek.

Here, I ensure all quadratic functions have matrices with smallest eigenvalue larger than 1.

x = binvar(n,1);
m = min(eig(A0));
r = 1-m;
objective = x'*(A0+r*eye(n))*x-r*sum(x)+b0^T*x;
Constraints = a'*x + b <= 0;
for i = 1:nc+nnc
m = min(eig(A{i}));
r = 1-m;
Constraints = [Constraints, x'*(A{i}+m*eye(n))*x-m*sum(x)+bi{i}^T*x + c(i)<= 0];
end
solvesdp(Constraints,objective)


Hopefully easily solved if you have any the mentioned solvers installed (if you don't, it will use the internal bnb solver, which is vaaaastly inferior)

[Johan Löfberg]

With the code I have shown above, it should be trivial. The only difference is the declaration of stuff X and Y

X =sdpvar(n);
Y =[1 x';x X];
Constraints =[Y>=0]

but as I said, a mixed-integer SDP based relaxation is a terrible choice with 99% certainty. Sure, as a fun test, but you must compare it with a convex MISOCP model.

Go from there, and if you run into problems, let me know where

[Rupaj Kumar Nayak]

Please check the following code and correct me. The optimal solution is -3.327.
yalmip('clear');
Q0=[21 17 0;17 -24 0;0 0 0];
c0=[2;-14;0];
Q1=[2 2 0;2 2 0;0 0 0];
Q2=[-5 -4 0;-4 -5 0;0 0 0];
Q3=[0 0 0; 0 0 0 ;0 0 0];
Q4=[1 0 0;0 0 0;0 0 0];
Q5=[0 0 0;0 1 0; 0 0 0];
Q6=[0 0 0; 0 0 0; 0 0 1];
Q=cat(7,Q1,Q2,Q3,Q4,Q5,Q6);
a=[8 -4 1 -1 0 0 ; 6 4 2 0 -1  0;  0 0 0 0 0 0];
b=[-9 4 -2 0 0 0 ];
m=6;n=3;

tic


   x = binvar(n,1);
    X =sdpvar(n,n);
   
objective = minimize(trace(Q0*X) + c0'*x )

 for i=1:m
 constraints=trace(Q(:,:,i)*X)+a(:,i)'*x+b(i)<=0;
end
  
     diag(X)==x;
    %Adding RLT to SDP
    for i=1:n
      for j=i:n
          if (i<=j)
          x(i)+x(j)-X(i,j)<=1;
          X(i,j)-x(i)<=0;
          X(i,j)-x(j)<=0;
           X(i,j)>=0;
          end
          
           
      end
    end
   
 % Adding traingle inequality
if(n>=3)
    for i=1:n
        for j=i:n
            for k=j:n
              
                if(i<=j && j<=k)
                   X(i,j)-X(i,k)-X(j,k)<=min(x(i),x(j))-min(x(j),x(k));
                     
               end
            end
        end
    end
end
 
 
    % finally, the constraints

    Y =[1 x';x X];

 Constraints =[Y>=0];
ops = sdpsettings('solver','sedumi');
solvesdp(constraints,objective)

toc

[Johan Löfberg]

Your problem is infeasible, both the exact convex MISOCP model, and the (unnecessary) semidefinite relaxation

Also, you seem to have missed the introductory examples on the Wiki, describing how to define constraints.

Here is the convexified MISOCP model. Infeasible

x = binvar(n,1);
r = min(eig(Q0))-1;
objective = x'*(Q0-r*eye(n))*x+r*sum(x)+c0'*x;
constraints=[];
for i=1:m
    r = min(eig(Q(:,:,i)))-1;
    constraints=[constraints,
x'*(Q(:,:,i)-r*eye(n))*x+a(:,i)'*x+r*sum(x)+b(i)<=0];
end
solvesdp(constraints,objective);

[Johan Löfberg]

Your two first quadratic inequalities are

2*x(1)^2+4*x(1)*x(2)+2*x(2)^2+8*x(1)+6*x(2) <= 9
-4*x(1)+4*x(2)-5*x(1)^2-8*x(1)*x(2)-5*x(2)^2 <= -4

The only feasible choice for the first is x1=0 and x2=1. However, this is infeasible in the second. Hence, problem is infeasible

Your data is probably wrong, or you are missing a constant somewhere