## 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

----     66 VARIABLE z.L        =      153.675  total transportation costs in \$1000

----     66 VARIABLE x.L  shipment quantities in cases

new-york     chicago      topeka

seattle        50.000     300.000
san-diego     275.000                 275.000

This was a very simple model model, about the simplest one can imagine. But even with this simple model we can illustrate the differences between a matrix oriented model system like Matlab and an equation based system like GAMS. Already in the Matlab code, I had to apply some careful matrix indexing tricks. The main reason is we need to map from variables $$x_{i,j}$$ to matrix columns (each variable is a column in the coefficient matrices $$A_1$$ and $$A_2$$), i.e. we have to map from a matrix to a vector. If we have multiple variables and/or variables that have an irregular shape, things become even more dicey in the Matlab approach.

A comparison of a slightly more complicated MIP model is shown in (4).

#### References

1. http://www.gams.com/modlib/libhtml/trnsport.htm
2. G.B.Dantzig, Linear Programming and Extensions, Princeton University Press, Princeton, New Jersey, 1963.
3. Robert Fourer, “Modeling languages versus matrix generators for linear programming”, ACM Transactions on Mathematical Software, Volume 9, Issue 2, June 1983, Pages 143-183.
4. http://yetanothermathprogrammingconsultant.blogspot.com/2016/10/matlab-vs-gams-integer-programming.html shows a Matlab example of a slightly more complex MIP where Matlabs matrix approach becomes somewhat nightmarish.