GAMS FAQ for keyword: MODELING
Missing trig functions (arccos, arcsin, tan)
> My model requires the tangent and inverse sine and cosine
> functions, but these functions are not built in to GAMS. What should
> I do?
All these functions can be constructed from the GAMS built-in
functions. For example, the definition of tan(x) is sin(x)/cos(x).
While computing arcsin and arccos from the GAMS built-in function
arctan is less trivial, the formulas, derivations, and some examples
are given in the ARCTRIG model.
Sensitivity analysis (parameter ranges) for LP
> How can I get sensitivity information (i.e. "parameter ranges") for
> the rhs and objective function coefficients in a linear program?
There are two types of sensitivity information that one can extract
from a solved GAMS model, information available from within GAMS
and information available only within the LP solver.
Upon solving an LP model, all GAMS solvers return 4 sets of values
to GAMS: level and marginal values for the primal variables, and
level and marginal values for the constraints. (N.B. The rhs term is
not included in the level values for the constraints.) From this
information, it is possible to compute, from within the GAMS model,
a range of parameter values for which the current solution (both
primal and dual variables) will not change. Such a range of parameter
values can be computed for all rhs coefficients having a zero marginal
value, where the range of parameter values depends on the current
level value for the constraint and the sense of the inequality in
question. Similarly, the set of obj. coefficients that can change
without changing the primal or dual level values coincides with the
set of primal variables currently at their bound. These coefficients
can move an infinite amount in one direction, and a finite amount
(determined by the marginal value or reduced cost) in the other.
A simple example helps to demonstrate how this is done.
Note that we cannot address the issue of changing level values
within the same basis in the GAMS model, since the GAMS model
itself does not store the current basis matrix.
To determine the range of values a particular parameter can take
without forcing a change in the current basis, but while allowing the
level values to change, one must have a representation of the current
basis. Since GAMS does not keep this information, these ranges
must be produced by the LP solver. Documentation on how this is done
is given elsewhere.
Documentation is also available for ANALYZE.
How large should BIG M be?
> How big should the constant M be in the equation, where Y(J) are
> binary variables:
>
> YBIN(J).. SUM(I,X(I,J)) =L= M*Y(J);
> I said in my model: SCALAR M /10000000000/ but that causes
> problems.
In general the answer should be as small as possible. Your model
with the M of 1.0E10 could not be solved by any of the MIP solvers
we tried out. It causes enormous scaling problems. However further
on in your model you have
CAP(J).. SUM(I,X(I,J)) =L= MAXCAP;
We can combine these two equations in one:
YBIN(J).. SUM(I,X(I,J)) =L= MAXCAP*Y(J);
Which both sets M to the smallest possible value and reduces the
number of constraints.
I have seen MIP models fail miserably when using a real big M.
Reusing optimal values
> How can I save the optimal values of a certain variable X that
> appears in one model, and reuse in another model that is otherwise
> unrelated, i.e. the save/restart facilities are not appropriate.
You may try the PUT statement to write out an include file for the
second model. These PUT statement could look like:
file pf /'inc.inc'/;
put pf;
put "PARAMETER P(I,J) /" /;
LOOP ((I,J),
PUT " ",I.TL,".",J.TL," ",X.L(I,J):20:10 /;
);
put "/;" /;
This would generate a file like:
PARAMETER P(I,J) /
SEATTLE .NEW-YORK 0.0000000000
SEATTLE .CHICAGO 300.0000000000
SEATTLE .TOPEKA 0.0000000000
SAN-DIEGO .NEW-YORK 325.0000000000
SAN-DIEGO .CHICAGO 0.0000000000
SAN-DIEGO .TOPEKA 275.0000000000
/;
The only thing to worry about is that the sets I and J are identical for
the first model and the second one. One way to make sure they are in
sync is to the put these set declarations in an include file.
How do I model piecewise linear functions?
> How do I model a piecewise linear function in GAMS. The
> following fragment is not accepted:
>
> Y =E= (a*X+b)$(X lt 0) + (c*X+d)$(X ge 0 and X lt 1) +
> (e*X+f)$(X ge 1);
That is a good question. With "if-then-else"'s within the
equations you make the model non-linear. Even non-linear in a nasty
sense: it can discontinuous. In your specific case it would make it
non-differentiable. The resulting model can not be solved using a
general LP or even NLP solver (NLP solvers like MINOS and
CONOPT like smooth continuous-differentiable functions).
In fact the piecewise linear function you describe can be modeled
with some MIP solvers using socalled SOS 2 variables, or with
general binary variables. Some textbooks like Nemhauser and
Wolsey, Integer and Combinatorial Optimization (page 10) and H.P.
Williams, Model Building in Mathematical Programming (section 7.3,
9.3).
Sometimes the following trick can be used: if you are minimizing Q
one can introduce the inequalities:
q >= aP+b q >= cP+d q >= eP+f
because of the form of the P-Q curve. [here a picture approximately
like:
|
|
\
\
-----
with two kinks].
(Note that the other way around does not work: q <= aP+b, y <=
cP+b, y <= eP+f would describe another feasible region!).
In most cases the best way would be to find a smooth approximation
for this function. One could fit for instance a polynomial or an
exponential. If needed you could construct a GAMS model for the
least squares fit!
Smooth approximations for MAX(X,0) and MIN(X,0)
> Do you know a smooth approximation for max(x,0), and min(x,0)?
This comes from Prof. Ignacio Grossmann (CMU):
Use the approximation
( sqrt( sqr(x) + sqr(epsilon) ) + x ) / 2
for max(x,0), where sqrt is the square root and sqr is the square.
The error err(x) in the above approximation is maximized at 0 (the
point of nondifferentiability), where err(0) = epsilon. As x goes to +/-
infinity, err(x) goes to 0. One can reduce the maximum error to
epsilon/2 by using the approximation given below. This provides a
better approximation near the point of nonsmoothness but is not so
accurate away from this point.
( sqrt( sqr(x) + sqr(epsilon) ) + x - epsilon ) / 2
Because min(x,0) = -max(-x,0), you can use the above
approximations for min(x,0) as well. Epsilon is a small positive
constant.
Modeling permutations
> How can I model in GAMS the following restriction? A vector of
> integer variables x(i) can assume values 1,2,...,card(i), and each x(i)
> should be different, i.e.
>
> x(i) <> x(j) for all i <> j
This looks like we need the help of a permutation matrix P=pij,
where there is exactly one 1 in each row and column, and the other
elements are zero. The identity matrix is a trivial example of a
permutation matrix. Other ones just have some rows and columns of
the identity matrix interchanged. To model this in a MIP use binary
variables P and add the constraints:
alias (i,j);
binary variables p(i,j);
equations rowsum(i) colsum(j) ;
rowsum(i).. sum(j, p(i,j)) =e= 1;
colsum(j).. sum(i, p(i,j)) =e= 1;
The different values for our vector x can now be calculated as:
x(i) =e= sum(j, ord(j)*p(i,j));
Because x is automatically integer if p is integer, we can declare x as
a normal continuous variable.
MIN function, don't use it
> I have the the following equations in my model:
>
> obj.. y=e=min(xa,xb,xc)+xd;
> model m /all/;
> solve m using dnlp maximizing y;
>
> MINOS was able to solve a very small instance of my model, but as
> the size increased MINOS was giving up.
The min function is a very dangerous one. To discourage users to
use this function we don't allow this in an NLP model, and you are
forced to declare the model as a DNLP model. The problem is that the
min function is not differentiable, just like the max and the abs. The
best advice I can give you is: stay away from them. In your case we
can reformulate the model quite easily to make it a normal NLP,
which solves without problems with MINOS:
obj1.. y=L=xa+xd ;
obj2.. y=L=xb+xd ;
obj3.. y=L=xc+xd ;
This is possible because we maximize y, which will assume
automatically in the optimal solution the minimum of the three
right-hand sides.
