Formulation of Intlinprog problem with quadratic terms

Hi,

I want to minimize a function with quadratic terms. I started by linearizing the terms according to the suggestions mentioned in: http://orinanobworld.blogspot.de/2010/10/binary-variables-and-quadratic-terms.html. The idea behind this function is to apply a price from a list of prices (p1, p2) to each row so that the total sum is the minimum possible. I am thinking to the branch and bound algorithm. So each row will have only one price applied. The consumption should be split between the 2 rows so that the calculated price is the minimum possible.

Intlinprog is returning -2 and don't like the constraints. Is the formulation correct ? If yes How can debug/improve the constraints ? Is there a simpler way to formulate ? I could'nt find any example with variable products and constraints.

Thank you.

Here is The code with all the details corresponding to the formulation:

clear all;

%% Objectiv Function to minimize

% f = (x11 + x12) * (z11.p1 + z12.p2) + ...

% (x21 + x22) * (z21.p1 + z22.p2)

% f = x11.z11.p1 + x11.z12.p2 + x12.z11.p1 + x12.z12.p2 + ...

% x21.z21.p1 + x21.z22.p2 + x22.z21.p1 + x22.z22.p2

% The idea behind this function is that X are consumptions and p1 and p2 are exclusivs prices. % Either p1 or p2 apply for each of the rows X1n and X2n % Z variables are binary and therefore their sum is always 1 : z11 + z12 = 1

% z21 + z22 = 1 % % Y variables are introduced to linearize the function as suggested in http://orinanobworld.blogspot.de/2010/10/binary-variables-and-quadratic-terms.html % % f = y11.p1 + y12.p2 + y13.p1 + y14.p4 + ...

% y21.p1 + y22.p2 + y23.p1 + y24.p4

%% Inequality contsraints % Applied the suggestions from % http://orinanobworld.blogspot.de/2010/10/binary-variables-and-quadratic-terms.html % %Let z be binary and let L≤x≤U where L and U are bounds known in advance. %Introduce a new variable y. (Programming note: Regardless of whether x is real-valued, integer-valued or binary, y can be treated as real-valued.)

%Add the following four constraints:

% y ≤ Uz ==> y - Uz ≤ 0

% y ≥ Lz ==> y - Lz ≥ 0

% y ≤ x−L(1−z) ==> -x + y - Lz ≤ -L

% y ≥ x-U(1−z) ==> -x + y - Uz ≥ -U %

After some Modifications: % 0x + y - Uz ≤ 0

% 0x + -y + Lz ≤ 0

% -x + y - Lz ≤ -L

% x - y + Uz ≤ U

%% lower upper bounds

ub = [40; 40; 60; 60; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; 1; 1; 1; 1];

lb = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 0; 0; 0; 0 ];

% f = x11.z11.p1 + x11.z12.p2 + x12.z11.p1 + x12.z12.p2 + ...

% x21.z21.p1 + x21.z22.p2 + x22.z21.p1 + x22.z22.p2

% f = y11.p1 + y12.p2 + y13.p1 + y14.p4 + ...

% y21.p1 + y22.p2 + y23.p1 + y24.p4

% A = [ x11 x12 x21 x21 y11 y12 y13 y14 y21 y22 y23 y24 z11 z12 z21 z22]

A =   [ 0  0  0  0         1   0   0  0      0   0  0  0     -ub(1) 0   0   0;
        0  0  0  0        -1   0   0  0      0   0  0  0      1  0   0   0;
       -1  0  0  0         1   0   0  0      0   0  0  0     -1  0   0   0;
        1  0  0  0        -1   0   0  0      0   0  0  0     ub(1)  0   0   0;

          0  0  0  0        0   1   0 0        0   0  0  0     -ub(1) 0   0   0;
          0  0  0  0        0  -1   0 0        0   0  0 0       1  0   0   0;
         -1  0  0  0        0   1   0 0        0   0  0 0      -1  0   0   0;
          1  0  0  0        0  -1   0 0        0   0  0 0     ub(1)  0   0   0;          

          0  0  0  0        0   0   1 0        0   0  0 0       0 -ub(2)  0 0;
          0  0  0  0        0   0  -1 0        0   0  0 0       0  1   0 0;
          0 -1  0  0        0   0   1 0        0   0  0 0       0 -1   0 0;
          0  1  0  0        0   0  -1 0        0   0  0 0       0 ub(2)   0 0;

          0  0  0  0        0   0   0  1       0   0  0 0       0 -ub(2) 0 0;
          0  0  0  0        0   0   0 -1       0   0  0 0       0  1  0 0;
          0 -1  0  0        0   0   0  1       0   0  0 0       0 -1  0 0;
          0  1  0  0        0   0   0 -1       0   0  0 0       0 ub(2)  0 0;

          0  0  0  0        0   0   0 0        1  0  0  0        0 0 -ub(3) 0;
          0  0  0  0        0   0   0 0       -1  0  0  0        0 0  1  0;
          0  0 -1  0        0   0   0 0        1  0  0  0        0 0 -1  0;
          0  0  1  0        0   0   0 0       -1  0  0  0        0 0 ub(3)  0;

          0  0  0  0        0   0   0 0        0  1  0  0        0 0 -ub(3) 0;
          0  0  0  0        0   0   0 0        0 -1  0  0        0 0  1  0;
          0  0 -1  0        0   0   0 0        0  1  0  0        0 0 -1  0;
          0  0  1  0        0   0   0 0        0 -1  0  0        0 0 ub(3)  0;

          0  0  0  0        0   0   0 0        0  0  1  0        0 0  0 -ub(4);
          0  0  0  0        0   0   0 0        0  0 -1  0        0 0  0  1 ;
          0  0  0 -1        0   0   0 0        0  0  1  0        0 0  0 -1 ;
          0  0  0  1        0   0   0 0        0  0 -1  0        0 0  0 ub(4) ;          

          0  0  0  0        0   0   0 0        0  0  0  1        0  0 0 -ub(4);
          0  0  0  0        0   0   0 0        0  0  0 -1        0  0 0  1 ;
          0  0  0 -1        0   0   0 0        0  0  0  1        0  0 0 -1 ;
          0  0  0  1        0   0   0 0        0  0  0 -1        0  0 0 ub(4) ;          

          0  0  0  0        0   0   0 0        0  0  0  0        1  1 1 1 ];

 B = [0; 0; -1; ub(1);
      0; 0; -1; ub(1);
      0; 0; -1; ub(2);
      0; 0; -1; ub(2);
      0; 0; -1; ub(3);
      0; 0; -1; ub(3);
      0; 0; -1; ub(4);
      0; 0; -1; ub(4);
      1];

%% equality Constraints

% x11 + x21 = 10

% x12 + x22 = 16

% z11 + z12 = 1

% z21 + z22 = 1

Beq = [10; 30; 1; 1];

Aeq = [1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1];

%% Prices

p1 = 2 ; p2 = 4 ; Prices = [p1; p2]; intcon = [13, 14, 15, 16];

%% Objectiv Function % f = x11.z11.p1 + x11.z12.p2 + x12.z11.p1 + x12.z12.p2 + ...

% x21.z21.p1 + x21.z22.p2 + x22.z21.p1 + x22.z22.p2

% f = y11.p1 + y12.p2 + y13.p1 + y14.p2 + ...

% y21.p1 + y22.p2 + y23.p1 + y24.p2

f = [0; 0; 0; 0; p1; p2; p1; p2; p1; p2; p1; p2; 0; 0; 0; 0];

%% Executer optimisation

options = optimoptions('intlinprog','Display','iter'); [x,fval,exitflag,output] = intlinprog(f,intcon,A,B,Aeq,Beq,lb,ub,options);

Answer 1

[Alan Weiss]

Without reading your problem in detail, I think that this example is relevant in that it shows a correct formulation of a MIQP solution via MILP.

Alan Weiss

MATLAB mathematical toolbox documentation

[Mohammed]

Thank you. It contains alot of valuable informations and explains many thnigs that I need to do. I still need to do some modif. on the problem formulation.