## Thursday, September 29, 2016

### Matlab vs GAMS: Linear Programming

Regularly I am asked “What is the difference between Matlab and GAMS when modeling optimization problems?”. I’ll try to answer this question with some examples.

#### Matlab

In Matlab LP and MIP models are modeled using a matrix data structure. Basically a linear optimization model is viewed as:

 \boxed{\begin{align}\min\>&c^Tx\\&A_1 x\le b_1\\&A_2 x=b_2\\&\ell\le x \le u\end{align}}

We need to populate the vectors $$c$$, $$b_1$$, $$b_2$$, $$\ell$$ and $$u$$ and the matrices $$A_1$$, $$A_2$$ and then call linprog or intlinprog.

Some notes:

• The matrices $$A_1$$, $$A_2$$ can become very, very large and sparse for larger models. In that case it may help to use a sparse matrix storage scheme, preventing all the zeros being stored. We’ll see an example of this later on in a subsequent post.
• There are alternative modeling frameworks available for use with Matlab, including Tomlab and CVX. These systems take a step toward a somewhat more equation based modeling paradigm. In some sense one could think about them as positioned between pure Matlab and a modeling system like AMPL and GAMS.
• Instead of the Matlab optimization toolbox solvers linprog or intlinprog, you can also use high-performance solvers such as Cplex and Gurobi. They have a similar matrix oriented modeling interface as the Matlab solvers.
• Matlab’s quadratic solver quadprog also has a similar matrix oriented interface.
• General NLP (nonlinear programming) problems have a totally different interface in Matlab. They require the objective and constraints to be modeled as functions. Introducing nonlinearities to an LP model will force you to rewrite the the model code.
• In practice Matlab’s matrix approach works best for highly structured models, and is not as convenient for more complex and messy models.
• Many LP/MIP models in R are using an matrix-based modeling paradigm. Much of what I am saying here also applies to R.
• Indexing in Matlab is done by integers, and array indexing (vectors, matrices) is one-based.
• Matlab default variable bounds are $$-\infty \le x \le \infty$$. Most LP and MIP solvers use a different default: if a variable has no bounds it is considered positive.

#### GAMS

GAMS is an equation based modeling system. This means an optimization model is collection of algebraic equations. For instance a capacity constraint can look like:

 Mathematics:$\sum_j x_{i,j} \le cap_{i} \>\> \forall i$GAMS:Capacity(i)..    sum(j, x(i,j)) =l= cap(i);

Of course we can write a whole LP as just one or two equations, but that is almost never a good modeling approach:

 objective..   z =e= sum(j, c(j)*x(j));e(i)..   sum(j, a(i,j)*x(j)) =e= b(i);

Notes:

• GAMS is similar to other modeling systems such as AMPL. These notes largely apply to those modeling systems as well.
• Equation based modeling allows an easy transition to non-linear programming. We don’t have to change the modeling paradigm when the model becomes non-linear.
• A properly designed model will is both compact and allows us to think and reason about the model. Often we don’t need to maintain a separate document that states the mathematical model.
• A modeling system helps to be efficient in writing and maintaining complex models.
• GAMS is somewhat unusual in not having an objective function but rather an objective variable. The objective function is formulated as a normal constraint.
• By default GAMS variables are free variables. As said before most other modeling software assume by default that variables are non-negative.
• Indexing in GAMS is done using strings.

#### A simple LP transportation problem

Let’s start with an extremely simple transportation problem from the GAMS model library (1) and originally from (2) (though with slightly different data, making the model degenerate).

The mathematical model looks like:

 \boxed{\begin{align}\min\>&\sum_{i,j} c_{i,j} x_{i,j}\\&\sum_j x_{i,j} \le a_i\>\>\forall i\\&\sum_i x_{i,j} \ge b_j\>\>\forall j\\&x_{i,j}\ge 0\end{align}}

Below is a side-by-side comparison of the Matlab code and the GAMS model. As you will see the Matlab matrix approach is very different from an equation based language.

Matlab code GAMS code

NumPlants = 2;

NumMarkets = 3;

MatLab uses indexing by integers.
GAMS is set based with set elements
represented by strings
sets
i
canning plants /seattle, san-diego /
j
markets        /new-york, chicago, topeka /
;

a = [350 600]';     % capacity

b = [325 300 275]'; % demand

d = [2.5 1.7 1.8;   % distance

2.5 1.8 1.4];

f = 90;             % unit freight cost

c = f * d / 1000;   % freight

Ordering is important. We need to remember index=1 means Seattle or New-York depending on the situation.
Ordering is not important. parameters
a(i)
capacity of plant i in cases

/  seattle     350

san-diego   600  /

b(j)
demand at market j in cases

/  new-york    325

chicago     300

topeka      275  / ;

table d(i,j)  distance in thousands of miles

new-york   chicago   topeka

seattle        2.5       1.7       1.8

san-diego      2.5       1.8       1.4  ;

scalar f  freight in dollars per case per 1000 miles /90/
;

parameter c(i,j)  transport cost in $1000 per case ; c(i,j) = f * d(i,j) / 1000 ; % allocate matrices Asupply = zeros(NumPlants,NumPlants*NumMarkets); Ademand = zeros(NumMarkets,NumPlants*NumMarkets); % populate submatrices for i=1:NumPlants for j=1:NumMarkets k = (i-1)*NumMarkets + j; Asupply(i, k) = 1; Ademand(j, k) = 1; end end % populate matrix A and vectors B,C A = [Asupply;-Ademand]; B = [a;-b]; C = reshape(c',NumPlants*NumMarkets,1); % positive variables lb = zeros(NumPlants*NumMarkets,1); Lots of index calculations. Note that we flip the ≥ demand constraint to a ≤ constraint. The matrix A looks like: A = 1 1 1 0 0 0 0 0 0 1 1 1 -1 0 0 -1 0 0 0 -1 0 0 -1 0 0 0 -1 0 0 -1 The GAMS code for the equations is much closer to the mathematical model. variables x(i,j) shipment quantities in cases z total transportation costs in$1000 ;

positive variable
x ;

equations

cost
define objective function
supply(i)
observe supply limit at plant i
demand(j)
satisfy demand at market j ;

cost..        z =e=
sum
((i,j), c(i,j)*x(i,j)) ;

supply(i)..
sum
(j, x(i,j)) =l= a(i) ;

demand(j)..
sum
(i, x(i,j)) =g= b(j) ;

options = optimoptions('linprog','Algorithm','dual-simplex');

[x,z] = linprog(C,A,B,[],[],lb,[],[],options);

x=reshape(x,NumMarkets,NumPlants)';

z,x

model transport /all/ ;
solve
transport using lp minimizing z ;
display
z.l,x.l;

The output of the two models is similar. As the LP model has multiple optimal solutions, we actually find two different solutions.

Matlab output GAMS output

z =

153.6750

x =

0   300     0
325     0   275

##### Alternative

If we can make the function $$f(x,y)$$ separable, we may be able to use a different scheme. E.g. assume $$f(x,y)=x\cdot y$$ with $$x,y>0$$. We can make this separable by applying a log transformation: $$\log(z)=\log(x)+\log(y)$$. We can implement each of the functions $$z’=\log(z),x’=\log(x),y’=\log(y)$$ using standard 1D piecewise linear formulations. Then we just add the constraint: $$z’=x’+y’$$.

Of course another alternative is using a nonlinear solver.

##### References
1. H.P.Williams, “Model Building in Mathematical Programming”, 5th edition, Wiley, 2013.
2. Misener, R. & Floudas, C.A., “Piecewise-Linear Approximations of Multidimensional Functions”, J. Optim. Theory Appl. (2010), pp.145-120.
In this paper SOS1 sets are used for the $$\delta$$ and $$\gamma$$ variables and SOS2 sets for the $$\beta$$ variables. Sadly, Chris Floudas passed away last month while vacationing in Greece.
3. GAMS: Piecewise linear functions with SOS2 variables. Basic GAMS syntax and rules for SOS2 sets.
4. 2D Interpolation With SOS2 variables. Example GAMS model.