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);
Best Answer