$title Purchasing with price breaks
$ontext

$offtext

$set NB 480
$set NS 140

set BREAK0 'price break points' /0*%NB%/;
set BREAK1(BREAK0) 'price break curves' /1*%NB%/;
set SUPPL 'Suppliers' /1*%NS%/;

alias(s,SUPPL,i);
alias(b,b2,BREAK0);


parameters
  BASECOST(s)  'base Unit cost'
  BASEBREAK(s) 'base break level'
  COST(s,b)    'Unit costs'
  BR(s,b)    'Breakpoints (quantities at which unit cost changes)'
;

scalar BSIZE;
BSIZE = round(uniform(1,4));
option seed = 666 ;
BASECOST(s) = uniform(10,12);
BASEBREAK(s) = round(uniform(50,150));
COST(s,b) = BASECOST(s) * power(uniform(0.85,0.95), ord(b)-1) ;
BR(s,b)$(ord(b)>=2) = BASEBREAK(s)*BSIZE*(ord(b)-1)*uniform(0.9,1.1);

display COST;
display BR;


parameter CBR(SUPPL,BREAK0) 'Total cost at break points';
parameter MAXPERC(SUPPL) 'Maximum percentages for each supplier';

MAXPERC(s) = 10*round(uniform(100,250)/card(s)) ;
display MAXPERC;



scalar REQ 'Total quantity required';

REQ = sum((s,b)$(ord(b)=round(card(b)/2)), BR(s,b))/card(s) ;
display REQ;



positive variables buy(SUPPL) 'Quantity to purchase from supplier s';
sos2 variables w(SUPPL,BREAK0) 'Weights at breakpoint k for supplier s';
variables MinCost;
equations defobj, defbuy(s), OneW(s), Demand;

* Calculate total cost at breakpoints
CBR(s,'0') = 0;
loop(b$BREAK1(b),
  CBR(s,b)= CBR(s,b-1) + COST(s,b)*(BR(s,b)-BR(s,b-1));
);

* Objective: sum of costs*weights
defobj..
  MinCost =e= sum((s,b), CBR(s,b) * w(s,b));

* Define buy and also order the weight variables by breakpoint quantities
defbuy(s)..
  sum(b, BR(s,b) * w(s,b)) =e= buy(s);

* The convexity row (w sum to 1)
OneW(s)..
  sum(b, w(s,b)) =e= 1;

* The minimum quantity that must be bought
Demand..
  sum(s, buy(s)) =g= REQ;

* No more than the maximum percentage from each supplier
buy.up(s) = MAXPERC(s)*REQ/100.0;

model purchase /all/;
purchase.optcr = 1e-6;
solve purchase using mip minimizing MinCost;

parameter sos2_time;
sos2_time = purchase.resusd;


* Try to model it with binary variables.  Define all new variables
positive variables x(i), lambda(i,b);
binary variables z(i,b);

equations
  defobj_2, defbuy_2(s), demand2, convexity(i), sosone(i), lamz(i,b)
;

defobj_2..
  MinCost =E= sum((s,b), CBR(s,b) * lambda(s,b));

defbuy_2(s)..
  sum(b, BR(s,b) * lambda(s,b)) =e= x(s);

convexity(i)..
  sum(b, lambda(i,b)) =E= 1;

demand2..
  sum(i, x(i)) =G= REQ ;

x.up(i) = MAXPERC(i)*REQ/100.0;

sosone(i)..
  sum(b,z(i,b)) =E= 1;

lamz(i,b)..
  lambda(i,b) =L= z(i,b-1) + z(i,b);

model purchase_sos1 / defobj_2, defbuy_2, demand2, convexity, sosone, lamz/ ;
purchase_sos1.optcr = 1e-6;
solve purchase_sos1 using mip minimizing MinCost;

parameter sos1_time;
sos1_time = purchase_sos1.resusd;

display sos1_time, sos2_time;