$title Happy Milk example with concave cost functions

option limrow=0, limcol=0;

* abscissae and ordinates of the piecewise-linear cost function
set pieces /1*3/;

set      region /A,B/
         milkType /raw,hifat,lofat/
         costfuncDefs /grad1, grad2, break, maxmilk/
         productType /cream, milk/

table butterfat(region,milkType)
       raw        hifat    lofat
A      .05         .80      .01
B      .035        .80     .005;


* cost DECREASES after the breakpoint; cannot use an LP model

table costfunc(region,costfuncDefs)
      grad1      grad2   break    maxmilk
A       .56      .555      500     5000
B       .44       .43      700     5000;

* define the piecewise linear cost functions via abscissae and ordinates

parameter Xpurchase(region,pieces), Xcost(region,pieces);
loop(region,
 Xpurchase(region,'1') = 0; 
 Xcost(region,'1') = 0;
 Xpurchase(region,'2') = costfunc(region,'break');
 Xcost(region,'2') = costfunc(region,'break') * costfunc(region,'grad1');
 Xpurchase(region,'3') = costfunc(region,'maxmilk');
 Xcost(region,'3') = Xcost(region,'2') + 
    costfunc(region,'grad2')*(Xpurchase(region,'3') - Xpurchase(region,'2'));
)

display Xcost, Xpurchase;
 
parameter
         separationCost(region) / A  .05, B  .06 /
         maxDemand(productType) / cream     300, milk      2000 /
         minRichness(productType) / cream    .50, milk     .02 /
         prices(productType) / cream   4.50, milk      .80 /;

sos2 variables gamma(region,pieces);
variables
         purchase(region)
         amtsAfterSeparation(region,milkType)
	 amtsMixed(region,milkType,productType)
         amtProduct(productType),
         saleableProduct(ProductType),
         purchaseCost(region), 
	 totalCost, 
	 sales, 
	 profit;

positive variable purchase, 
         amtsAfterSeparation, amtsMixed, amtProduct, 
	 saleableProduct, purchaseCost, totalCost, sales;

equations
	 EQpiecewise(region)                   constraints on the SOS2
	 EQpurchase(region)                    figure purchase from the pieces
	 EQcost(region)                        figure purchase cost
         fatConservationSeparation(region)     fat conserved in separation
         volumeConservationSeparation(region)  volume conserved in separation
	 volumeMixed(region,milkType)          mix at most what's available
	 totalProduct(productType)             volume of each final product
         qualityProduct(productType)           min fat content of final product
         howMuchCanYouSell(productType)        amt sold may be <= amt produced
         productionCost                        total cost of production
         revenue                               sales revenue
         objective;

* constraints on the SOS2 variabes that model the piecewise linear function
EQpiecewise(region)..
	sum(pieces, gamma(region,pieces)) =e= 1;

EQpurchase(region)..
	purchase(region) =e= 
             sum(pieces, gamma(region,pieces)*Xpurchase(region,pieces));

EQcost(region)..
	purchaseCost(region) =e=
	     sum(pieces, gamma(region,pieces)*Xcost(region,pieces));

* ensure that the total amount of fat in the products of separation equals
* the total fat in the raw milk purchased

fatConservationSeparation(region)..
	purchase(region)*butterfat(region,"raw") =e= sum(milkType, 
   amtsAfterSeparation(region,milkType) * butterfat(region,milkType));

* ensure that volume is conserved in separation

volumeConservationSeparation(region)..
	purchase(region) =e= sum(milkType, amtsAfterSeparation(region,milkType));

* ensure that the amounts of separated product contributed to 
* final product satisfy a volume conservation constraint

volumeMixed(region,milkType)..
	sum(productType, amtsMixed(region,milkType,productType)) =l=
	   amtsAfterSeparation(region,milkType);

* calculate the amount of each final product

totalProduct(productType)..
         amtProduct(productType) =e= sum((region,milkType),
           amtsMixed(region,milkType,productType));

* ensure that each final product has the required content of fat

qualityProduct(productType)..
         sum((region,milkType),
         butterfat(region,milkType) * amtsMixed(region,milkType,productType))
         =g= minRichness(productType) * amtProduct(productType);

* from the amounts produced, extract the (possibly smaller) amount 
* that can actually be sold.

howMuchCanYouSell(productType)..
         saleableProduct(productType) =l= amtProduct(productType);

* get total cost by adding purchase cost to separation cost
productionCost..
	totalCost =e= 
  sum(region, purchaseCost(region)) +
  sum(region, (amtsAfterSeparation(region,"hifat")
             + amtsAfterSeparation(region,"lofat")) 
     * separationCost(region));

* calculate total sales revenue

revenue..
         sales =e= sum(productType,
           saleableProduct(productType) * prices(productType));

objective..
         profit =e= sales - totalCost;


model milky / all /;

* set upper bound on sales

saleableProduct.up(productType) = maxDemand(productType);

option mip=cplex;

* important to set the RELATIVE precision small (setting optca doesn't help!)
milky.optcr = 1.e-6;

solve milky using mip maximizing profit;

display purchase.l, amtsAfterSeparation.l, amtProduct.l, profit.l;