Introduction

In this chapter we will introduce special GAMS features that are useful for specific model types. The features include model scaling, conic programming and features that facilitate mixed integer as well as indicator constraints, a feature that does not translate across solvers.

Special Mixed Integer Programming (MIP) Features

Some special features have been added to GAMS to help simplifying the modeling of Mixed Integer Programming (MIP) problems. In GAMS MIP is the model type for mixed integer linear programs, this section used MIP more generally, we consider model with discrete variables, including non-linear expressions and pure discrete problem. We will first present details on discrete variables in GAMS, then we will discuss how to customize priorities for the branching process. Next, we will cover the model attributes that are important for MIPs and we will conclude with some hints that will make mixed integer programming with GAMS easier.

Types of Discrete Variables

Variables and variable types are introduced in chapter Variables. GAMS provides six discrete variable types: binary, integer, sos1, sos2, semicont and semiint. In the following subsections we will present details and examples for each of these discrete variable types. Note that if any discrete variables feature in a model, it has to be a mixed integer model or one of the related model types, like MINLP or MIQCP. See section Classification of Models for a full listing of all GAMS model types.

Binary Variables

Binary variables can take values of 0 (zero) and 1 (one) only. They are declared as follows:

Binary Variable var_name [(index_list)] [text];

The keyword binary indicates that this is a binary variable and then the usual conventions for variable declarations are followed. Alternatively, the variable may be declared first and specified as binary later. Consider the following code snippets from the orthogonal Latin Square model [LATIN]:

Sets k  "rows"      / row1*row4 /
     l  "columns"   / col1*col4 /
     v  "values"    / val1*val4 /;
alias (i,j,v);

Variables x(i,j,k,l) "pairs (i,j) allocated to cell(k,l)"
          z          "some objective";
Binary Variable x;

Equations c1(i,j)    "for each cell pick only one item pair";
c1(i,j).. sum((k,l), x(i,j,k,l)) =e= 1;

Note that the binary variable x is used in equation c1 to model the restriction that in each cell only one item pair is allowed. Binary variables are often used to model logical conditions such as imposing mutual exclusivity or complementarity.

Note that the default lower bound is 0 (zero) and the default upper bound is 1 (one). If the relaxed versions of the discrete models is solved, binary variables are treated like positive variables with the upper bound of 1. In addition, an infinite priority may be used to override binary specifications, see section Setting Priorities for Branching below for more information.

Even though the only possible values are 0 and 1, a solver might return a value for binary variable that is only close to 0 or 1. Every solver works with tolerances and also uses a tolerance to determine if a value is close enough to an integer values. So it is unwise to use code as a(i)$(b.l(i)=1) = yes; because one will potentially miss some elements. A safe way to write such code is: a(i)$(b.l(i)>0.5) = yes;. Rounding the level of a binary variable after the solve is also possible, but it is not done by the solver or the solver link because even small rounding can lead to infeasibilities.

A binary variable can also have a truely fractional value after a solver if the model status does not indicate a feasible integer solution (model status 1 or 8).

Integer Variables

Integer variables are discrete variables that can take only values between their bounds. The user may change both bounds from the default value. The default lower bound is 0 (zero) and the default upper bound inside GAMS is +inf, and the same upper bound is passed on to the solver.

Note that in relaxed model types the integrality requirement is relaxed. In addition, an infinite priority may be used to override integer specifications, see section Setting Priorities for Branching below for more information. Integer variables are declared as follows:

Integer Variable var_name [(index_list)] [text];

The keyword integer indicates that this is an integer variable and then the usual conventions for variable declarations are followed. Alternatively, the variable may be declared first and specified as integer later. Consider the following code snippets from the power scheduling model [MAGIC]:

Sets t  "demand blocks" / 12pm-6am, 6am-9am, 9am-3pm, 3pm-6pm, 6pm-12pm /
     g  "generators"    / type-1, type-2, type-3 /;

Variables  x(g,t)  "generator output (1000mw)"
           n(g,t)  "number of generators in use"
           cost    "total operating cost (l)";
Integer Variable n;

The integer variable n models the number of generators of various types that are in use at any of the time blocks.

Special Order Sets of Type 1 (SOS1)

SOS1 variables are a set of variables, such that at most one variable within the group may have a nonzero value. This variable may take any positive value. Special ordered sets of type 1 are defined as follows:

SOS1 Variable var_name (index_list) [text];

The keyword SOS1 indicates that this is a SOS1 variable and then the usual conventions for variable declarations are followed. Alternatively, the variable may be declared first and specified as SOS1 later. Consider the following example:

SOS1 Variable s1(i), t1(k,j), w1(i,j,k) ;

Note that the members of the innermost (the right-most) index belong to the same SOS set. For example in the sets defined above, s1 represents one special ordered set of type 1 with i elements, t1 defines k sets with j elements each and w1 defines (i,j) sets with k elements each.

The default bounds for SOS1 variables are zero and +inf. As with any other variable, the user may change these bounds. Further, the user may explicitly provide whatever convexity row that the problem may need through an equation that requires the members of the SOS1 set to be less than a certain value. Any such convexity row will implicitly define bounds on each of the variables.

Consider the following example:

SOS1 Variable s1(i); Equation defsoss1; 
defsoss1.. sum(i,s1(i)) =l= 3.5 ;

The equation defsoss1 implicitly defines the nonzero value that one of the elements of the SOS1 variable s1 may take as equal to or smaller than 3.5. Note that it is also possible that all variables s1 equal zero.

A special case arises when one of the elements of the set has to be nonzero and equal to a number, say 3.5. In this case equation defsoss1 will be:

defsoss1.. sum(i,s1(i)) =e= 3.5 ;

Frequently the nonzero value equals 1. As a result, the SOS1 variable is effectively a binary variable. It is only treated differently by the solver at the level of the branch and bound algorithm. For example, consider the following example where we want to model that one out of n options has to be selected. This is expressed as:

SOS1 Variable x(i); Equation defx;
defx.. sum(i, x(i)) =e= 1 ;

The variable x can be made binary without any change in meaning and the solution provided by the solver will be indistinguishable from the SOS1 case.

The use of special ordered sets may not always improve the performance of the branch and bound algorithm. If there is no natural order the use of binary variables may be a better choice. A good example of this is the classical assignment problem (see [H.P. Williams (2013) Model Building in Mathematical Programming], Wiley, Section 9.3.

Note that any model with SOS1 variables requires a MIP solver, because the solution process needs to impose the restrictions of at most one nonzero level values may be present.

For an example where SOS1 variables are used, see the production scheduling model [PRODSCHX].

Special Order Sets of Type 2 (SOS2)

SOS2 variables are a set of variables, such that at most two variables within the set may have nonzero values and these variables have to be adjacent. This requirement implies that the set is ordered, see chapter Sets as Sequences: Ordered Sets for details on ordered sets in GAMS. Note that the nonzero variables may take any positive value. Special ordered sets of type 2 are defined as follows:

SOS2 Variable var_name [(index_list)] [text];

The keyword SOS2 indicates that this is a SOS2 variable and then the usual conventions for variable declarations are followed. Alternatively, the variable may be declared first and specified as SOS2 later. Consider the following example:

Set i  / i1*i5 /;
SOS2 Variable s2(i), t2(k,j), w2(i,j,k);

The members of the innermost (the right-most) index belong to the same set. For example, in the sets defined above, s2 represents one special ordered set of type 2 with elements for each member of the set i. At most two variables s2 may be nonzero and they must reference adjacent elements of the set i. Note that the variables s2('i1') and s2('i2') are adjacent, but the variables s2('i1') and s2('i3') are not. Further, t2 defines k sets of SOS2 variables with j elements each and the adjacency requirement refers to the set j which must be ordered. Similarly, w2 defines (i,j) sets with k elements each and the adjacency requirement refers to the set k which must be ordered.

The default bounds for SOS2 variables are zero and +inf. As with any other variable, the user may change these bounds. SOS2 variables are most often used to model piece-wise linear approximations to nonlinear functions. The production scheduling model [PRODSCHX] shows SOS type formulations with binary, SOS1 and SOS2 sets.

Note that any model with SOS2 variables requires a MIP solver, because the solution process needs to impose the restrictions of adjacency and that no more than two nonzero level values may be present.

Semi-Continuous Variables

Semi-continuous variables are either zero or above a given minimum level. This can be expressed algebraically as: either $x=0$ or $L\leq x\leq U$. By default, the lower bound $L$ is 1 and the upper bound $U$ is +inf. As usual, these bounds may be changed with the variable attributes .lo and .up. Semi-continuous variables are defined as follows:

SemiCont Variable var_name [(index_list)] [text];

The keyword semicont indicates that this is a semi-continuous variable and then the usual conventions for variable declarations are followed. Alternatively, the variable may be declared first and specified as semicont later. Consider the following example:

SemiCont Variable x;
x.lo = 1.5;  x.up = 23.1;

The slice of code above declares the variable x to be a semi-continuous variable that may either be zero or behave as a continuous variable between 1.5 and 23.1.

Note that any model with semi-continuous variables requires a MIP solver, because the solution process needs to impose the discontinuous jump between zero and the threshold value.

Note

Not all MIP solvers allow semi-continuous variables. We recommend users to verify how the solver they are interested in handles semi-continuous variables by checking the relevant section of the respective solver manual.
The lower bound has to be less than the upper bound, and both bounds have to be greater than zero, otherwise GAMS will report an error.
The variable solution listing might show the level outside the lower and upper bound which for other variables indicates an infeasible variable, but not so for semi-continuous variables.
Semi-continuous variables are especially helpful if the upper bound is +inf and no implicit bound can be easily derived. If a finite upper bound is available it can be computational more efficient to replace the semi-continuous variable sc with lower bound scLow by a continuous variable x and binary variable b and the following equations:
```
Equation xForceLowerBnd   "Force x to be greater than scLow if b is 1"
         xForceZero       "Force x to be zero if b is zero";
xForceLowerBnd.. x =g= scLow*b;
xForceZero..     x =l= x.up*b;
```

Semi-Integer Variables

Semi-integer variables are either zero or integer and above a given minimum value. This can be expressed algebraically as: either $x = 0$ or $x \in \{L,\ldots,U\}$. By default, the lower bound $L$ is 1 and the upper bound $U$ inside GAMS is +inf and the same values are passed on to the solver. As usual, these default bounds may be changed with the variable attributes .lo and .up. Note that in relaxed model types the integrality requirement is relaxed. In addition, an infinite priority may be used to override integer specifications, see section Setting Priorities for Branching below for more information.

Semi-integer variables are defined as follows:

SemiInt Variable var_name [(index_list)] [text];

The keyword semiint indicates that this is a semi-integer variable and then the usual conventions for variable declarations are followed. Alternatively, the variable may be declared first and specified as semiint later. Consider the following example:

SemiInt Variable x;
x.lo = 2; x.up = 25;

The slice of code above declares the variable x to be a semi-integer variable that may either be zero or take any integer value between 2 and 25. Note that the bounds for semiint variables have to take integer values, otherwise GAMS will flag an error during model generation. Note further, that any model with semi-integer variables requires a MIP solver.

Note

Not all MIP solvers allow semi-integer variables. We recommend users to verify how the solver they are interested in handles semi-integer variables by checking the relevant section of the respective solver manual.
The lower bound has to be less than the upper bound, and both bounds have to be greater than zero, otherwise GAMS will report an error.
The variable solution listing might show the level outside the lower and upper bound which for other variables indicates an infeasible variable, but not so for semi-integer variables.
Semi-integer variables are especially helpful if the upper bound is +inf and no implicit bound can be easily derived (together with the appropriate IntVarUp setting). If a finite upper bound is available, it can be computationally more efficient to replace the semi-integer variable si, with lower bound siLow, by an integer variable i and a binary variable b and the following equations:
```
Equation iForceLowerBnd   "Force i to be greater than siLow if b is 1"
         iForceZero       "Force i to be zero if b is zero";
iForceLowerBnd.. i =g= siLow*b;
iForceZero..     i =l= i.up*b;
```

Setting Priorities for Branching

By setting priorities users may specify an order for choosing variables to branch on during a branch and bound search for MIP models. Without priorities the MIP algorithm will internally determine which variable is the most suitable to branch on. Priorities for individual variables may be used only if the model attribute .prioropt is set to 1; the respective GAMS statement is:

mymodel.prioropt = 1;

Here mymodel is the name of the model specified in the model statement. The default value is NA.

If the model attribute .prioropt is set to 1, the variable attribute .prior may be used to set the priorities of individual discrete variables. Note that there is one .prior value for each individual component of a multidimensional variable. Priorities may be set to any real value; the default value is 1. As a general rule, the most important variables should be given the highest priority. The highest priority is denoted by the lowest nonzero value in the .prior attribute. Functionally, the attribute .prior establishes in what order variables are to be branched on in the branch-and-bound algorithm while searching for a solution. Variables with a specific .prior value will branched on earlier until all fractional variables with higher .prior values have been branched on.

Note: The variable attribute .prior of a discrete variable may be used to relax the discrete restriction on that variable: setting the .prior value to +inf will relax a variable permanently (or until .prior gets a finite value assigned). This relaxation is done independently of the model attribute .prioropt.

Consider the following example:

z.prior(i,'small')   = 3;
z.prior(i,'medium')  = 2;
z.prior(i,'large')   = 1;

In this example the variables z(i,'large') are branched on before the variables z(i,'medium'), which in turn are branched on before the variables z(i,'small').

Note that knowledge about the problem may help to determine which variables should be considered first. For example, consider a problem with a binary variable u representing a yes/no decision whether to build a factory and other binary variables representing equipment selections within that factory. We would naturally want to explore whether or not the factory should be built before considering what specific equipment to be purchased within the factory. Therefore we would set the priority values lower for u. By assigning a higher priority - a lower value of the attribute .prior - to the build/nobuild decision variable u, we can force this logic into the tree search and thus speed up computation time since uninteresting portions of the tree are left unexplored.

Note

The lower the value given to the .prior suffix, the higher the priority for branching.
All members of any SOS1 or SOS2 set should be given the same priority value since it is the set itself which is branched upon rather than the individual members of the set.
While any value is accepted for .prior many solvers scale all giving priorities in the integer range of ${0,...,1000}$.
Branching priorities were a very important feature in the early days of mixed integer programming. Nowadays, it is not easy to find branching priorities that improve on the solvers default selection.
Global non-linear optimization solvers branch on continuous variables too (see, for example BARON). In GAMS one cannot set the branching priority of a continuous variable. Such branching priorities need to be communicated via a solver option file.

Miscellaneous Hints

We will conclude the discussion of mixed integer models in GAMS with this section where we offer a variety of hints that are meant to make special facilities of mixed integer programming solvers more accessible.

Model Attributes for Mixed Integer Programming in GAMS

GAMS offers several model attributes that may be used to influence MIP solver performance or report on results of MIPs. These model attributes include Cheat, CutOff, NodLim, ObjEst, OptCA, OptCR, PriorOpt and TryInt.

The Branch and Cut and Heuristic Facility

Hard MIP problems can be solved faster with the help of user supplied routines that generate cutting planes and good integer feasible solutions. The GAMS Branch-and-Cut-and-Heuristic (BCH) automates all major steps necessary to define, execute and control the use of user defined routines within the framework of general purpose branch-and-cut codes. It is documented in Branch-and-Cut-and-Heuristic Facility (BCH).

Branch and Bound Output

While the log output for each solver differs, some key figures are usually displayed for branch-and-bound based solvers. For example, solving a linear mixed integer model with CPLEX will yield output like the following:

       Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                            0.0000   5.05646e+08              --- 
Found incumbent of value 0.000000 after 0.01 sec. (0.73 ticks)
      0     0   2.23031e+07    12        0.0000   2.23031e+07       17     --- 
      0     0   2.23031e+07     7        0.0000       Cuts: 8       23     --- 
*     0+    0                       2.23031e+07   2.23031e+07             0.00%
Found incumbent of value 2.2303094e+07 after 0.02 sec. (1.08 ticks)
 ...
 
Fixing integer variables, and solving final LP...
...

Solution satisfies tolerances.

MIP Solution:     22303093.628684    (23 iterations, 0 nodes)
Final Solve:      22303093.628684    (0 iterations)

Best possible:    22303113.950091
Absolute gap:           20.321407
Relative gap:            0.000001

A brief explanation of the columns follows:

Node is the number of branch and bound nodes so far.
Nodes Left is the number of problems created during the branching process that are yet to be examined.
Objective gives the current objective function value of the relaxes node problem.
IInf gives the number of discrete variables with fractional solution levels.
Best Integer gives the incumbent solution. Note the last solution in that column is not necessarily the global best solution.
Cuts/Best Bound gives the current lower bound on the solution.
ItCnt gives the accumulated LP iteration count
Gap gives the maximum percentage difference from the theoretical optimum.

Note that it is common that solves of mixed integer models end with a gap between the solution found and the best possible solution. This may be controlled by limits (e.g. time), solver options and model attributes like .OptCR and .OptCA.

Nonlinear MIPs

Modelers may wish to impose integer restrictions on nonlinear formulations combining two hard model types: MIP and NLP. Such MINLP models can be solved with a selection of solvers. Many solvers, e.g. DICOPT and SBB, provide a local optimum where others, e.g. ANTIGONE and BARON provide a global optimum. In most cases both types of MINLP solver make use of MIP and NLP solvers to calculate a solution. Such subsolvers need to be licensed for the solver to succeed.

Model Termination Conditions and Recommended Actions

Recall that the termination condition of the model after the solution process has been completed is stored in the model attribute .modelStat. A list of all possible model statuses is given in section Model Status. We can easily check for the existence of a feasible solution (status 1 and 8). Linear models and MINLP problems solved with global solvers can achieve the "OPTIMAL" status, given sufficient resources (time) and a setting of OptCR and OptCA to 0. All other cases do not yield a feasible integer solution. If a problem is reported as infeasible (status 4,5,10, and 19), it might be a good idea to see if the relaxed version of the model is already infeasible. Debugging models that are relaxed feasible but integer infeasible is very difficult.

Frequent Problems

There are some problems users frequently encounter either due to GAMS settings or problem characteristics:

Default bounds

One needs to be aware that while the GAMS upper bound for integer and semi-integer variables is +inf, the bound that is passed to the solver can be different, namely 100 (see the discussion about integer variables). This can lead to unexpected results (e.g. infeasibilities or suboptimal solutions declared as optimal).

The nonending quest

Integer programming is a quite desirable formulation technique. However, integer problems are theoretically hard and the solution process (in the worst case) of exponential complexity. There are many ways that focus on improving the solution time of the solver. As with all models scaling is important (especially when using bigM formulations). For particular problems many different formulations exist and the literature about a particular problem together with the ability of GAMS to rapidly prototype and experiment is the best constellation to get the best results for the problem at hand. For fine tuning, some MIP solvers provide automated tuning tools (see e.g. Cplex tuning) that tweak the solver options to get the best performance.

Model Scaling - The Scale Option

The rules for good scaling are exclusively based on algorithmic needs. GAMS has been developed to increase the efficiency of modelers, and one of the best ways seems to be to encourage modelers to write their models using a notation that is as natural as possible. The units of measurement are one part of this natural notation. However, there is a potential conflict between what the modeler thinks is a good unit and what constitutes a well-scaled model.

The Scale Option

To facilitate the translation between a natural model and a well scaled model, GAMS has introduced the concept of a scale factor, both for variables and equations. The notations and definitions are quite simple. Scaling is turned off by default. Setting the model attribute .scaleopt to 1 turns on the scaling feature. For example,

model mymodel /all/ ;
mymodel.scaleopt = 1 ;
solve mymodel using nlp maximizing dollars ;

The statement should be inserted somewhere after the model statement and before the solve statement. To turn scaling off again, mymodel.scaleopt has to be set to zero before the next solve.

In most respects GAMS scaling is hidden from the user. The solution values reported back from a solution algorithm are always reported in the notation of the user. The algorithm's internal representation of the equations and variables are only reflected in the derivatives in the equation and column listings in the GAMS output if the values of the options limrow and limcol are positive. In addition, the internal representations will appear in the debugging output from the solution algorithm if the option sysout is set to on.

Scaling Variables

The scale factor of a variable is defined using the variable attribute .scale in the following way:

myvar.scale(i,j) = c;

The scale factor c is a number or a numerical expression that evaluates to a number. Note that the default scale factor is 1.

Note that there is one scale value for each individual component of a multidimensional variable.

Assume that $c$ is the scale factor of a variable $V_u$. Assume further, that the variable seen by the algorithm is $V_a$. Then we have: $ \mathbf{V_a =V_u/c}$. This means that each variable as seen by the user is divided by the scale factor.

For example, consider the following code snippet:

Positive Variables x1, x2;
Equation eq;
eq.. 200*x1 + 0.5*x2 =l= 5;
x1.up = 0.01;
x2.up = 10;
x1.scale = 0.01;
x2.scale = 10;

By setting x1.scale to 0.01 and x2.scale to 10, the model seen by the solver is:

Positive Variables xPrime1, xPrime2;
Equation eq;
eq.. 2*xPrime1 + 5*xPrime2 =l= 5;
xPrime1.up = 1;
xPrime2.up = 1;

Note that the solver does not see the variables x1 or x2, but rather the scaled (and better-behaved) variables xPrime1 and xPrime2. Note further, that upper and lower bounds on variables are automatically scaled in the same way as the variable itself.

Attention

Discrete variables cannot be scaled.
Expert Note. Internally, GAMS stores with each variable and equation one additional attribute or field (besides fields for level, marginal and lower and upper bound). Depending on the type of variable and sometimes even model type or solver, this field has different names in the GAMS language. For continuous variables, the field is called scale, while for discrete variable it is called prior. For stochastic 2-stage linear program models solved with DECIS, this field is called stage. The field .stage can lead to confusing results. Consider the following example:
```
Variable x;
x.scale = 0.1;
display x.stage;
```
The output is:
```
----      3 VARIABLE x.scale = 0.100
```
The field .scale has to be in a certain range ( >1e-20 and no special value), but this is only checked at model generation time. The field .prior can be any number and even +inf (but no other special values). For further information on .prior, see section Setting Priorities for Branching. For an introduction to variable and equation fields, see sections Variable Attributes and Equation Attributes respectively.

Scaling Equations

The scale factor of an equation is defined using the equation attribute .scale in the following way:

mzeqn.scale(i,j) = d;

The scale factor d is a number or a numerical expression that evaluates to a number. Note that the default scale factor is 1.

Assume that $d$ is the scale factor of an equation $G_u$. Assume further, that the equation seen by the algorithm is $G_a$. Then we have: $ \mathbf{G_a =G_u/d}$. This means that each equation as seen by the user is divided by the scale factor.

For example, consider the following equations:

Positive Variables y1, y2;
Equations eq1, eq2;
eq1..   200*y1 + 100*y2 =l= 500;
eq2..       3*y1 - 4*y2 =g= 6;

By setting eq1.scale to 100, the model seen by the solver is:

Positive Variables y1, y2;
Equations eqPrime1, eq2;
eqprime1..   2*y1 + 1*y2 =l= 5;
eq2..        3*y1 - 4*y2 =g= 6;

Note: The user may have to perform a combination of equation and variable scaling to obtain a well-scaled model.

Consider the following example:

Positive variables x1, x2;
Equations  eq1, eq2;
eq1..   100*x1 + 5*x2 =g= 20;
eq2..   50*x1 - 10*x2 =l= 5;
x1.up = 0.2;
x2.up = 1.5;

Setting the following scale values:

x1.scale  = 0.1;
eq1.scale = 5;
eq2.scale = 5;

will result in the solver seeing the following well-scaled model:

Positive Variables xPrime1, x2;
Equations eqPrime1, eqPrime2;
eqPrime1..   2*xPrime1 + x2 =g= 4;
eqPrime2..   xPrime1 - 2*x2 =l= 1;
xPrime1.up = 2;
x2.up      = 1.5;

Scaling Derivatives

In nonlinear models the derivatives also need to be well-scaled. Assume that the derivatives in the model of the user are denoted by $d(G_u)/d(V_u)$. Assume further, that the derivatives in the scaled model seen by the algorithm are denoted by $d(G_a)/d(V_a)$. Then we have: $\mathbf{d(G_a)/d(V_a) = d(G_u)/d(V_u) \cdot c/d}$, where $c$ is the scale factor for the variable and $d$ is the scale factor for the equation.

The user may affect the scaling of derivatives by scaling both the equation and variable involved.

Scaling Data

Scaling input data is independent of the model attribute .scaleopt and may contribute considerably towards achieving a well-scaled model. We recommend users to try to define the units of the input data such that the largest values expected for decision variables and their marginals is under a million, if possible.

For example, in US agriculture about 325 million acres are cropped and the corn crop is 9-10 billion bushels per year. When defining production data, we could enter land in 1000's of acres and all other resources in 1000's of units. We could also define the corn crop in millions of bushels. The data will be simultaneously scaled, hence if resource endowments are quoted in 1000's, corn yields are divided by 1000. This scaling results in a corn production variable in the units of millions. Consumption statistics would need to be scaled accordingly. Money units could also be in millions or billions of dollars. Such data scaling generally greatly reduces the disparity of coefficients in the model.

Conic Programming in GAMS

Conic programming models minimize a linear function over the intersection of an affine set and the product of nonlinear cones. The problem class involving second order (quadratic) cones is known as Second Order Cone Programs (SOCP). These are nonlinear convex problems that include linear and (convex) quadratic programs as special cases.

Conic programs allow the formulation of a wide variety of application models, including problems in engineering and financial management. Examples are portfolio optimization, Truss topology design in structural engineering, Finite Impulse Response (FIR) filter design and signal processing, antenna array weight design, grasping force optimization, quadratic programming, robust linear programming and norm minimization problems.

For more information, see References and Links at the end of this section.

Introduction to Conic Programming

Conic programs can be thought of as generalized linear programs with the additional nonlinear constraint $x \in C$, where $C$ is required to be a convex cone. The resulting class of problems is known as conic optimization and has the following form:

\[ \begin{array}{rl} \text{minimize} & c^Tx \\ \text{subject to} & Ax \le r^c, \\ & x \in [l^x, u^x] \\ & x \in C \\ \end{array} \]

where $A\in \mathbb{R}^{m \times n}$ is the constraint matrix, $x \in \mathbb{R}^n$ the decision variable and $c \in \mathbb{R}^n$ the objective function cost coefficients. The vector $r^c \in \mathbb{R}^m$ represents the right-hand side and the vectors $l^x, u^x \in \mathbb{R}^n$ are lower and upper bounds on the decision variable $x$.

Now partition the set of decision variables $x$ into sets $S^t, t=1,...,k$, such that each decision variables $x$ is a member of at most one set $S^t$. For example, we could have

\[ S^1 = (x_1, x_4, x_7) \quad \text{and} \quad S^2 = (x_6, x_5, x_3, x_2). \]

Let $x_{S^t}$ denote the variables $x$ belonging to set $S^t$. Then define

\begin{equation*} C := \left \{ x \in \mathbb{R}^n : x_{S^t} \in C_t, t=1,...,k \right \}, \end{equation*}

where $C_t$ must have one of the following forms:

Quadratic cone (also referred to as Lorentz or ice cream cone):
\[ C_t = \left \{ x \in \mathbb{R}^{n^t} : x_1 \ge \sqrt{\sum_{j=2}^{n^t}x_j^2} \right \}. \]
Rotated quadratic cone (also referred to as hyperbolic constraints):
\[ C_t = \left \{ x \in \mathbb{R}^{n^t} : 2x_1x_2 \ge \sum_{j=3}^{n^t}x_j^2, ~x_1,x_2 \ge 0 \right \}. \]

These two types of cones allow the formulation of quadratic, quadratically constrained and many other classes of nonlinear convex optimization problems.

Implementation of Conic Constraints in GAMS

The recommended way to write conic constraints is by using a quadratic formulation. Many solvers have the capability to identify the conic constraints in a QCP model even if it is not in perfect form but can be easily reformulated to fit in the described form. However, some solvers (namely MOSEK) expect the conic constraints to be precisely in the form given above. In earlier versions this form was enforced by a special type of equation, the =c= equation type. Moreover, such solvers have other requirements (e.g. disjunctive cones) that can be easily fulfilled by simple reformulation steps. Much progress is expected on the solver side in the coming years, so we don't go into much detail here.

Observe that we could formulate conic problems as regular NLPs using the following constraints:

Quadratic cone:

x('1') =g= sqrt[ sum(i$[not sameas(i,'1')], sqr[x(i)]) ];

Rotated quadratic cone:

2*x('1')*x('2') =g= sum(i$[not sameas(i,'1') and not sameas(i,'2')], sqr[x(i)]);

Here x('1') and x('2') are positive variables.

The following example illustrates the different formulations for conic programming problems. Note that a conic optimizer usually outperforms a general NLP method for the reformulated (NLP) cone problems.

Example

Consider the following example, which illustrates the use of rotated conic constraints. We will give reformulations of the original problem in regular NLP form and in conic form (with conic constraints).

The original problem is:

\begin{align} \text{minimize} \; & \sum_{i=1}^n \frac{d_i}{x_i} \\ \text{subject to}\; & a\,x \le b \\ & x_i \in [l_i,u_i], & i=1,\ldots,n \end{align}

where $x \in \mathbb{R}^n$ is the decision variable, $d, a, l, u \in \mathbb{R}^n$ are parameters with $l_i>0$ and $d_i \ge 0$ and $b \in \mathbb{R}$ is a scalar parameter. The original model may be written in GAMS using the equations:

defobj..  sum(n, d(n)/x(n)) =e= obj;
e1..      sum(n, a(n)*x(n)) =l= b;
Model orig /defobj, e1/;
x.lo(n) = l(n);
x.up(n) = u(n);

We can write an equivalent QCP formulation by using the substitution $t_i=1/x_i$ in the objective function and adding a constraint. As we are dealing with a minimization problem, $d_i \ge 0$ and $x_i \ge l_i > 0$, we can relax the equality $t_ix_i=1$ into an inequality $t_ix_i \ge 1$, which results in an equivalent problem with a convex feasible set:

\begin{align} \text{minimize} \; & \sum_{i=1}^n d_i t_i \\ \text{subject to}\; & a\,x \le b \\ & t_i x_i \ge 1, & i=1,\ldots,n \\ & x \in [l,u], \\ & t \ge 0, \\ \end{align}

where $t \in \mathbb{R}^n$ is a new decision variable. The GAMS formulation of this QCP is:

defobjc..     sum(n, d(n)*t(n)) =e= obj;
e1..          sum(n, a(n)*x(n)) =l= b;
coneqcp(n)..  t(n)*x(n)         =g= 1;

Model cqcp /defobjc, e1, coneqcp/;
t.lo(n) = 0;
x.lo(n) = l(n);
x.up(n) = u(n);

Note that the constraints $t_i x_i \ge 1$ are almost in rotated conic form. If we introduce a variable $z \in \mathbb{R}^n$ with $z_i = \sqrt{2}$, then we can reformulate the problem using conic constraints as:

\begin{align} \text{minimize} \; & \sum_{i=1}^n d_i t_i \\ \text{subject to}\; & a\,x \le b \\ & z_i = \sqrt{2}, & i=1,\ldots,n \\ & 2 t_i x_i \ge z_i^2, & i=1,\ldots,n \\ & x \in [l,u],\\ & t \ge 0, \\ \end{align}

The GAMS formulation using conic equations is:

defobjc..        sum(n, d(n)*t(n)) =e= obj;
e1..             sum(n, a(n)*x(n)) =l= b;
e2(n)..          z(n) =e= sqrt(2);
coneperfect(n).. 2*x(n)*t(n) =g= sqr(z(n));

Model cperfect  /defobjc, e1, e2, coneperfect/;
t.lo(n) = 0;
x.lo(n) = l(n);
x.up(n) = u(n);

The complete model is listed below:

Set       n / n1*n10 /;
Parameter d(n), a(n), l(n), u(n);
Scalar    b;

d(n) = uniform(1,2);
a(n) = uniform (10,50);
l(n) = uniform(0.1,10);
u(n) = l(n) + uniform(0,12-l(n));

Variables x(n);
x.l(n) = uniform(l(n), u(n));
b      = sum(n, x.l(n)*a(n));

Variables t(n), z(n), obj;
Equations defobjc, defobj, e1, e2(n), coneqcp(n), coneperfect(n), conenlp(n);

defobjc..          sum(n, d(n)*t(n)) =e= obj;
defobj..           sum(n, d(n)/x(n)) =e= obj;
e1..               sum(n, a(n)*x(n)) =l= b;
coneqcp(n)..       t(n)*x(n) =g= 1;
e2(n)..            z(n) =e= sqrt(2);
coneperfect(n)..   2*x(n)*t(n) =g= sqr(z(n));

Model cqcp     /defobjc, e1, coneqcp/;
Model cperfect /defobjc, e1, e2, coneperfect/;
Model orig     /defobj, e1/;

t.lo(n) = 0;
x.lo(n) = l(n);
x.up(n) = u(n);

Option qcp=cplex;
Solve cqcp min obj using qcp;
Option qcp=mosek;
Solve cperfect min obj using qcp;
Solve orig min obj using nlp;

Sample Conic Models in GAMS

Conic models in the GAMS model library include:

[EMFL]: A multiple facility location problem,
[FDESIGN]: Linear Phase Lowpass Filter Design,
[IMMUN]: Financial Optimization: Risk Management,
[PMEANVAR]: Mean-Variance Models with variable upper and lower Bounds,
[QP7]: A portfolio investment model using rotated quadratic cones (quadratic program using a Markowitz model),
[ROBUSTLP]: Robust linear programming as an SOCP,
[SPRINGCHAIN]: Equilibrium of System with Piecewise Linear Spring,
[TRUSSM]: Truss Toplogy Design with Multiple Loads

References and Links

A. Ben-Tal and A. Nemirovski, Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications, MPS/SIAM Series on Optimization, SIAM Press, 2001.
M. Lobo, L. Vandenberghe, S. Boyd and H. Lebret, Applications of Second-Order Cone Programming, Linear Algebra and its Applications, 284:193-228, November 1998, Special Issue on Linear Algebra in Control, Signals and Image Processing.
MOSEK ApS, MOSEK Modeling Cookbook, 2015.
G. Pataki G and S. Schmieta, The DIMACS Library of Semidefinite-Quadratic-Linear Programs. Tech. rep., Computational Optimization Research Center, Columbia University, 2002.
Seventh Dimacs Implementation Challenge on Semidefinite and Related Optimization Problems.

Indicator Constraints

An indicator constraint is a way of expressing relationships between variables by specifying a binary variable to control whether or not a constraint takes effect. For example, indicator constraints are useful in problems where there are fixed charges to express only if a given variable comes into play.

So-called Big M formulations require to estimate an upper bound of an expression in a model. In most cases the model data can be used to determine relatively small number for such coefficients. In some cases it is not possible to find small Big M values and a resulting solution may exhibit trickle flow and have other unwanted side-effects. The main purpose of indicator constraints is to overcome these limitations of Big M formulations. Generally, the use of indicator constraints is not warranted when the unwanted side-effects of Big M formulations are not present.

Consider the following example:

constr01.. x1 + x2 + x3 =l= 1e+9*y; // may cause problems

Here we use a Big M formulation, that relies on the $x$ values to sum to less than one billion. Note that this formulation may cause numeric instability or undesirable solutions in some situations. Alternatively, we could use an indicator constraint:

constr01$(y=0).. x1 + x2 + x3 =l= 0; // alternative

Note that $y$ is a binary variable; the logical condition makes sure that the constraint is only active if $y=0$. Unfortunately, the $() expressions do not allow endogenous variables and we need a different way to specify such implications. Indicators are supported by CPLEX, GUROBI, SCIP, and XPRESS.

We should mention a formulation of indicators without the explicit indicator constraints. We can use a SOS1 constraint to express that a constraint holds or not based on a binary variable by having a SOS1 set of a slack variable and the binary variable SOS1(slack,1-y):

Positive variable slack;
constr01.. x1 + x2 + x3 =l= slack; // alternative

Recall that only one of the variables in a SOS1 constraint can be non-zero, so if 1-y is non-zero (i.e. y=0) then slack must be zero, i.e. the constraint holds. If 1-y=0 is non-zero, i.e. y=1 then slack can be used to make the constraint feasible for any setting of x1, x2, and x3. Here two unrelated variable, slack and y are part of an SOS1 constraint and some tricks are required to formulate this properly in GAMS:

Set oo / slack, yExpr /;
SOS1 variable indic(oo);
Equation iOn, iOff;
Positive variable slack;
constr01.. x1 + x2 + x3 =l= slack;
iOn..      indic('slack') =e= slack;
iOff..     indic('yExpr') =e= 1-y;

Indicator Constraints with GAMS

In the remainder of this section we will describe how to specify indicator constraints when using GAMS. Please consult the corresponding solver manual for information on whether indicator constraints are supported and possible differences in their specification.

The example from above will be implemented with a constraint in combination with additional information in a solver option file, e.g. CPLEX: cplex.opt:

constr01.. x1 + x2 + x3 =l= 0;

and the following entry in the option file:

indic constr01$y 0

This has the following effect: equation constr01 will become an indicator constraint and becomes active in a solution where the binary variable takes the value 0. If the value of $y$ in a solution is 1, the constraint is not active.

Note that this way of entering an indicator constraint is dangerous since the option files changes the model (usually an option file has some effect on the performance of the algorithm). Therefore, the solver will abort if there is a problem processing the indicator options in the solver option file.

Attention: If the model is given to a solver without the option file containing the indicator mapping (or to a solver that does not understand the indic keyword, a very different model will be solved. The current implementation of indicator constraints requires a significant amount of caution from the user.

There are two ways of entering the equation/binary variable mapping in a solver option file: with an indexed format and using labels.

The indexed format is a convenient shorthand notation which borrows its syntax from the GAMS syntax. It requires that the indices for the binary variable are already present in the index set of the equation.

Consider the following example of an invalid GAMS syntax with an endogenous variable in the dollar condition:

equ1(i,j,k)$(ord(i) < ord(j) and bin1(i,k)=1).. lhs =l= rhs;

This may be specified with the following equation in the GAMS file:

equ1(i,j,k)$(ord(i) < ord(j)).. lhs =l= rhs;

plus a solver option file with the following entry:

indic equ1(i,j,k)$bin1(i,k) 1

The label format is used in cases where the binary variable indices are not present in the equation indices or the binary variable is adjusted with lags or leads. In these cases the mapping of all individual equations and variables of the indicator constraints need to be specified. An example follows.

Set             i /i1*i3/, j /j1*j2/;
Binary variable bin1(j);
Equation        equ1(i,j);

equ1(i,j)$(bin1(j++1)=0)..   lhs =e= 0;

Note that the example above is not valid GAMS code. Instead, we will combine a valid GAMS equation and a solver option file using the label format as follows:

equ1(i,j).. lhs =e= 0;

and the solver option file

indic equ1('i1','j1')$bin1('j2') 0
indic equ1('i1','j2')$bin1('j1') 0
indic equ1('i2','j1')$bin1('j2') 0
indic equ1('i2','j2')$bin1('j1') 0
indic equ1('i3','j1')$bin1('j2') 0
indic equ1('i3','j2')$bin1('j1') 0

Note that the lines in such option files need not be entered manually. They may be easily generated using the GAMS The Put Writing Facility. For example, the lines above may be generated as follows:

file fcpx / cplex.opt /;
fcpx.pc=8;
loop((i,j), put fcpx 'indic' equ1(i,j) '$' bin1(j++1) '0'/);
putclose fcpx;

Attention: Mixing and matching between indexed format and label format is not possible. While expressions like eq(i,j,k)$x(i,k) and eq('i1','j2','k3')$x('i1','k6') are valid, combinations such as eq(i,'j2',k)$x(i,k) are not supported.

There are situations where the indicator binary variable exist in the indicator constraint only and hence will not be generated by GAMS to be passed on to the solver. In such cases the solver will issue the following error message:

Error: Column is not of type Variable

There is an easy way to fix this problem: adding the binary indicator variable artificially to the model. For example, it may be added with the coefficient eps to the objective:

defobj.. z =e= ... + eps*sum(j, bin1(j));

An Example for Indicator Constraints with GAMS

In this subsection we will comment on parts of the model [TRNSINDIC]. This model uses big M formulations and indicator constraints to solve the same problem. In addition, a formulation that makes it easy to switch between these two is presented. It is a fixed-charge network example based on the well-known model [TRNSPORT].

Recall that i is the set of canning plants and j is the set of markets where cases of some product are to be shipped. First, the basic model is reformulated to a MIP by introducing the binary variable use(i,j) and two new equations:

Binary Variable use(i,j)     is 1 if arc is used in solution;

Equations       minship(i,j) ensure minimum shipping
                maxship(i,j) ensure zero shipping if use variable is 0;

minship(i,j)..  x(i,j) =g= minshipping*use(i,j);
maxship(i,j)..  x(i,j) =l= bigM*use(i,j);

Note that minshipping is a scalar denoting the minimum amount of cases that may be shipped and bigM is a sufficiently large number, as usual.

Next, the same problem is solved with indicator constraints: the two new equations are reformulated and a CPLEX option file with information on indicator constraints is added:

Equations  iminship(i,j) ensure minimum shipping using indicator constraints
           imaxship(i,j) ensure zero shipping if use variable is 0 using indicator constraints;

iminship(i,j)..   x(i,j) =g= minshipping;
imaxship(i,j)..   x(i,j) =e= 0;

Model indicatorModel /cost, supply, demand, iminship, imaxship/ ;

file fcpx Cplex Option file / cplex.opt /;
putclose fcpx 'indic iminship(i,j)$use(i,j) 1' / 'indic imaxship(i,j)$use(i,j) 0';
indicatorModel.optfile = 1;

Solve indicatorModel using mip minimizing z ;

Note that the option file contains an entry for each of the two equations. Note further, that the binary variable use moved from the equations to the option file. However, it also features in the objective equation of the model, therefore this is not problematic. Observe that the indexed format for the equation/binary variable mapping was used in the option file. Alternatively, the label format may be used:

loop((i,j),
     put fcpx 'indic ' iminship.tn(i,j) '$' use.tn(i,j) yes
            / 'indic ' imaxship.tn(i,j) '$' use.tn(i,j) no / );
  putclose fcpx;

In a final step the model is reformulated again such that the same problem may be solved with and without indicator constraints. This can be especially useful for debugging a model with indicator constraints. For details, see [TRNSINDIC].

Another example of a model with indicator constraints is the model [BILINEAR] where various formulations to represent bilinear product terms are demonstrated.

Table of Contents

Introduction

Special Mixed Integer Programming (MIP) Features

Types of Discrete Variables

Binary Variables

Integer Variables

Special Order Sets of Type 1 (SOS1)

Special Order Sets of Type 2 (SOS2)

Semi-Continuous Variables

Semi-Integer Variables

Setting Priorities for Branching

Miscellaneous Hints

Model Attributes for Mixed Integer Programming in GAMS

The Branch and Cut and Heuristic Facility

Branch and Bound Output

Nonlinear MIPs

Model Termination Conditions and Recommended Actions

Frequent Problems

Model Scaling - The Scale Option

The Scale Option

Scaling Variables

Scaling Equations

Scaling Derivatives

Scaling Data

Conic Programming in GAMS

Introduction to Conic Programming

Implementation of Conic Constraints in GAMS

Example

Sample Conic Models in GAMS

References and Links

Indicator Constraints

Indicator Constraints with GAMS

An Example for Indicator Constraints with GAMS