[ Home | Downloads | Documentation | Solvers | APIs | Tools | Model Libraries | Resources | Sales | Support | Contact Us | Search ]

Special Language Features

# Introduction

This chapter introduces special features in GAMS that do not translate across solvers, or are specific to certain model types. These features can be extremely useful for relevant models, and are among the most widely used.

# Special MIP Features

Some special features have been added to GAMS to help in simplifying the modeling of MIP problems. Two special types of discrete variables are defined and discussed. Finally, creating priorities for the discrete variables is discussed. The solvers use this information when solving the problem.

## Types of Discrete Variables

The following types of discrete variables have been discussed so far in the book,

• binary variables These can take on values of 0 or 1 only.
• integer variables These can take on integer values between the defined bounds. The default lower and upper bounds are 0 and 100 respectively.

In addition to these two, two new types of discrete variables that are introduced in this section. Both these variables exploit special structures in MIP models during the solution phase. These are the following

• Special Ordered Sets (SOS) The precise definition of special ordered sets differ from one solver to another and the development of these features has been driven more by internal algorithmic consideration than by broader modeling concepts. GAMS offers sos1 and sos2 variables as two types of compromise features that model special ordered sets. Sections Special Order Sets of Type 1 (SOS1) and Special Order Sets of Type 2 (SOS2) discuss these two types of variables in greater detail.

The presence of any of the above types of discrete variables requires a mixed integer model and all the discreteness is handled by the branch and bound algorithm in the same way as binary and general integer variables are handled.

## Special Order Sets of Type 1 (SOS1)

At most one variable within a SOS1 set can have a non-zero value. This variable can take any positive value. Special ordered sets of type 1 are defined as follows,

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


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

Attention
• The default bounds for SOS1 variables are 0 to $$+\infty$$. As with any other variable, the user may set these bounds to whatever is required.
• The user can, in addition, explicitly provide whatever convexity row that the problem may need through an equation that requires the members of the SOS set to be less than a certain value. Any such convexity row would 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 non-zero values that one of the elements of the SOS1 variable s1 can take.

A special case of SOS1 variables is when exactly one of the elements of the set have to be non-zero. In this case, the defsoss1 equation will be

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


A common use of the use of this set is for the case where the non-zero value is 1. In such cases, the SOS1 variable behaves like 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 to model the case where at most one out of $$n$$ options can be selected. This is expressed as

sos1 variable x(i)
equation defx ;
defx.. sum(i,x(i)) =l= 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 assignment problem.

Attention
Not all MIP solvers allow SOS1 variables. Furthermore, among the solvers that allow their use, the precise definition can vary from solver to solver. Any model that contains these variables may not be transferable among solvers. Please verify how the solver you are interested in handles SOS1 variables by checking the relevant section of the Solver Manual.

## Special Order Sets of Type 2 (SOS2)

At most two variables within a SOS2 set can have non-zero values. The two non-zero values have to be adjacent. The most common use of SOS2 sets is to model piece-wise linear approximations to nonlinear functions.

Attention
The default bounds for SOS2 variables are 0 to $$+\infty$$. As with any other variable, the user may set these bounds to whatever is required.

Special ordered sets of type 2 are defined as follows,

sos2 Variable s2(i), t2(k,j), w2(i,j,k) ;


The members of the innermost (the right-most) index belongs to the same set. For example, in the sets defined above, s2 represents one special ordered set of type 2 with i elements, t2 defines k sets of j elements each, and w2 defines (i,j) sets with k elements each.

[PRODSCHX] shows SOS type formulations with binary, SOS1 and SOS2 sets. The default bounds for SOS variables are 0 to $$+\infty$$. As with any other variable, the user may set these bounds to whatever is required.

Attention
Not all MIP solvers allow SOS2 variables. Furthermore, among the solvers that allow their use, the precise definition can vary from solver to solver. Any model that contains these variables may not be transferable among solvers. Please verify how the solver you are interested in handles SOS2 variables by checking the relevant section of the Solver Manual.

## Semi-Continuous Variables

Semi-continuous variables are those whose values, if non-zero, must be above a given minimum level. This can be expressed algebraically as: Either $$x=0$$ or $$L\leq x\leq U$$.

By default, this lower bound ( $$L$$) is 1 and the upper bound ( $$U$$) is $$+\infty$$. The lower and upper bounds are set through .lo and .up. In GAMS , a semi-continuous variable is declared using the reserved phrase semicont variable. The following example illustrates its use.

semicont variable x ;
x.lo = 1.5 ; x.up = 23.1 ;


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

Attention
• Not all MIP solvers allow semi-continuous variables. Please verify that the solver you are interested in can handle semi-continuous variables by checking the relevant section of the Solver Manual.
• The lower bound has to be less than the upper bound, and both bounds have to be greater than 0. GAMS will flag an error if it finds that this is not the case.

## Semi-Integer Variables

Semi-integer variables are those whose values, if non-zero, must be integral above a given minimum value. This can be expressed algebraically as: Either $$x = 0$$ or $$x \in \{L,\ldots,U\}$$

By default, this lower bound ( $$L$$) is 1 and the upper bound ( $$U$$) is 100. The lower and upper bounds are set through .lo and .up. In GAMS , a semi-integer variable is declared using the reserved phrase semiint variable. The following example illustrates its use.

semiint variable x ;
x.lo = 2 ; x.up = 25 ;


The above slice of code declares the variable x to be semi-continuous variable that can either be 0, or can take any integer value between 2 and 25.

Attention
• Not all MIP solvers allow semi-integer variables. Please verify that the solver you are interested in can handle semi-integer variables by checking the relevant section of the Solver Manual.
• The lower bound ( $$L$$) has to be less than the upper bound ( $$U$$), and both bounds have to be greater than 0. GAMS will flag an error during model generation if it finds that this is not the case.
• The bounds for semiint variables have to take integer values. GAMS will flag an error during model generation if it finds that this is not the case.

## Setting Priorities for Branching

The user can specify an order for picking variables to branch on during a branch and bound search for MIP models through the use of priorities. Without priorities, the MIP algorithm will determine which variable is the most suitable to branch on. The GAMS statement to use priorities for branching during the branch and bound search is:

mymodel.prioropt = 1 ;


where mymodel is the name of the model specified in the model statement. The default value is 0 in which case priorities will not be used.

Using the .prior suffix sets the priorities of the individual variables. Note that there is one prior value for each individual component of a multidimensional variable. Priorities can be set to any real value. The default value is 1. As a general rule of thumb, the most important variables should be given the highest priority.

The following example illustrates its use,

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


In the above example, z(i,'large') variables are branched on before z(i, 'small') variables.

Attention
• 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.

# 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, and there is unfortunately 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 suffix .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. In order to turn scaling off again, set the model.scaleopt parameter to 0 before the next solve.

The scale factor of a variable or an equation is referenced with the suffix .scale, i.e. the scale factor of variable x(i)  is referenced as x.scale(i). Note that there is one scale value for each individual component of a multidimensional variable or equation. Scale factors can be defined using assignment statements. The default scale factor is always 1.

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

## Variable Scaling

The scale factor on a variable, $$V_s$$, is used to relate the variable as seen by the user, $$V_u$$, to the variable as seen by the algorithm, $$V_a$$, as follows: $$V_a =V_u/V_s$$

For example, consider the following equation,

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.

Attention
• Upper and lower bounds on variables are automatically scaled in the same way as the variable itself.
• Integer and binary variables cannot be scaled.

## Equation Scaling

Similarly, the scale factor on an equation, $$G_s$$, is used to relate the equation as seen by the user, $$G_u$$, to the equation as seen by the algorithm, $$G_a$$, as follows: $$G_a =G_u/G_s$$

For example, consider the following equations,

positive variables y1,y2 ;
equation 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 ;
equation eqprime1, eq2 ;
eqprime1.. 2*y1 + 1*y2 =l= 5 ;
eq2.. 3*y1 - 4*y2 =g= 6 ;

Attention
The user may have to perform a combination of equation and variable scaling until a well-scaled model is obtained.

Consider the following example,

positive variables x1,x2 ;
equation 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 ;
equation eqprime1, eqprime2 ;
eqprime1.. 2*xprime1 + x2 =g= 4 ;
eqprime2.. xprime1 - 2*xprime2 =l= 1 ;
xprime1.up = 2 ; x2.up = 1.5 ;


## Scaling of Derivatives

For nonlinear models, the derivatives also need to be well scaled. The derivatives in the scaled model seen by the algorithm, i.e. $$d(G_a)/d(V_a)$$ are related to the derivatives in the user's model, $$d(G_u)/d(V_u)$$ through the formula: $$d(G_a)/d(V_a) = d(G_u)/d(V_u) \cdot V_s/G_s$$.

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

# 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, for example:

• Portfolio Optimization
• Finite Impulse Response (FIR) Filter Design and Signal Processing
• Antenna Array Weight Design
• Grasping Force Optimization
• Robust linear programming
• Norm Minimization Problems

## 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 hyberbolic 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. However, some solvers accepts the conic constraints to be expressed with the =C= equation type. In this case, the conic equations are written as:

x('1') =C= sum(i$[not sameas(i,'1')], x(i));  • Rotated quadratic cone: x('1') + x('2') =C= sum(i$[not sameas(i,'1') and not sameas(i,'2')], x(i));


Note that the resulting nonlinear conic constraints result in "linear" constraints in GAMS. Thus the original nonlinear formulation is in fact a linear model in GAMS. We remark that we could formulate conic problems as regular NLP using constraints:

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)]);

where 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.

## Examples

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

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$$ parameters with $$l_i>0$$ and $$d_i \ge 0$$, and $$b \in \mathbb{R}$$ a scalar parameter. The original model can 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 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 (model cqcp) 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 =C= 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);
cone(n).. x(n) + t(n) =C= z(n);
Model clp /defobjc, e1, e2, cone/;
t.lo(n) = 0;
x.lo(n) = l(n);
x.up(n) = u(n);

Note that this formulation is a linear program in GAMS, although the constraints cone(n) represent the nonlinear rotated quadratic cone constraints.

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), cone(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);
cone(n).. x(n) + t(n) =C= z(n);
Model cqcp /defobjc, e1, coneqcp/;
Model clp /defobjc, e1, e2, cone/;
Model orig /defobj, e1/;
t.lo(n) = 0;
x.lo(n) = l(n);
x.up(n) = u(n);
Option qcp=cplexd;
Option lp=mosek;
Solve cqcp min obj using qcp;
Solve clp min obj using lp;
Solve orig min obj using nlp;

## Sample Conic Models in GAMS:

• emfl_socp.gms: Multiple facility location problem.
• fir_socp.gms: Linear phase lowpass filter design model
• qp7.gms: Portfolio investment model using rotated quadratic cones (quadratic program using a Markowitz model)
• springs.gms: Equilibrium of system with piecewise linear springs
• trussm.gms: Truss topology design problem with multiple loads

# Indicator Constraints

An indicator constraint is a way of expressing relationships among 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 often exhibit trickle flow, and at times, they behave unstably. The main purpose of indicator constraints is to avoid the unwanted side-effects 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.

For example, instead of the following Big M formulation, which relies on the $$x$$ values summing to less than one billion (a formulation that can cause numeric instability or undesirable solutions in some situations):

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


a natural way of entering an indicator constraint, where $$y$$ is a binary variable, would look like this:

constr01$(y=0).. x1 + x2 + x3 =l= 0; // alternative  In the following, we will describe how to specify indicator constraints when using GAMS/CPLEX. For other solver interfaces that support this construct, the syntax is similar. Please consult the corresponding solver manual for information on whether indicator constraints are supported and possible differences in their specification. The example from above would be implemented in the following way: constr01.. x1 + x2 + x3 =l= 0;  plus additional information in the Cplex option file cplex.opt: indic constr01$y 0


This has the following effect: Equation constr01 will become an indicator constraint and becomes effective 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, GAMS/CPLEX will abort if there is a problem processing the indicator options in the Cplex option file. Moreover, if the model is given to a different solver or to CPLEX without the option file containing the indicator mapping, a very different model is 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 Cplex option file:

1. using an indexed format,
2. using labels.

The indexed format is a convenient short hand notation which borrows its syntax from the GAMS syntax. It requires that the indexes for the binary variable are already present in the index set of the equation. For example, the following invalid GAMS syntax with an endogenous variable in the $-condition equ1(i,j,k)$(ord(i)<ord(j) and bin1(i,k)=1).. lhs =l= rhs;


would be specified like this in the GAMS file

equ1(i,j,k)$(ord(i)<ord(j)).. lhs =l= rhs;  plus the Cplex option file indic equ1(i,j,k)$bin1(i,k) 1


In cases where the binary variable indexes are not present in the equation indexes (or is adjusted using lags or leads), one needs to specify the mapping of all individual equations and variables of the indicator constraints. For example,

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;  This would be specified using the label format of indicator constraint as follows. The GAMS file equ1(i,j).. lhs =e= 0;  plus the Cplex 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  Such option files can be easily generated using The Put Writing Facility : file fcpx / cplex.opt /; fcpx.pc=8; loop((i,j), put fcpx 'indic' equ1(i,j) '$' bin1(j++1) '0'/); putclose fcpx;


There are situation where the binary variable shows up the in indicator constraints only and hence will not be generated by GAMS to be passed on to Cplex. In such cases Cplex will issue an error message:

Error: Column is not of type Variable


In such cases the binary variable has to be added artificially to the model, e.g. but adding them with cost EPS to the objective:

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


The GAMS model libraries contain a few examples that use indicator constraints. For example, check out the model [BILINEAR] in the GAMS Model Library.

Finally, we want to present a fixed-charge network example based on the [TRNSPORT] model from the GAMS Model library that used indicator constraints, big M forumulations and a formulation that makes it easy to switch between these two.

$Title Fixed Charge Transportation Problem with Indicator Constraints Sets i canning plants / seattle, san-diego / j markets / new-york, chicago, topeka / ; 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 thousand miles /90/ ; Parameter c(i,j) transport cost in thousands of dollars per case ; c(i,j) = f * d(i,j) / 1000 ; Parameter fixcost(i,j) fixed cost in thousands of dollars ; fixcost(i,j) = 10*d(i,j) / 1000 ; Scalar minshipping minimum shipping of cases /100/; Scalar bigM sufficiently large number; bigM = smax(i, a(i)); Variables x(i,j) shipment quantities in cases use(i,j) is 1 if arc is used in solution z total transportation costs in thousands of dollars ; Positive Variable x ; Binary Variable use; Equations cost define objective function supply(i) observe supply limit at plant i demand(j) satisfy demand at market j minship(i,j) ensure minimum shipping maxship(i,j) ensure zero shipping if use variable is 0; cost .. z =e= sum((i,j), c(i,j)*x(i,j) + fixcost(i,j)*use(i,j)) ; supply(i) .. sum(j, x(i,j)) =l= a(i) ; demand(j) .. sum(i, x(i,j)) =g= b(j) ; minship(i,j).. x(i,j) =g= minshipping*use(i,j); maxship(i,j).. x(i,j) =l= bigM*use(i,j); Option limrow=0, limcol=0, optcr=0, mip=cplex; Model bigMModel /all/ ; Solve bigMModel using mip minimizing z ; * Now let's build a model for the same problem using indicator constraints 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 ; * Let's do the same using indicator option with labels 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; Solve indicatorModel using mip minimizing z ; * Now let's build a model for the same problem that can be used with * and without indicator constraints. This can become handy when * debugging a model with indicator constraints Positive Variable minslack(i,j), maxslack(i,j); Equations xminship(i,j) ensure minimum shipping using indicator constraints and bigM xmaxship(i,j) ensure zero shipping ig use variable is 0 using indicator constraints and bigM bndminslack(i,j) ensure minslack is zero if use variable is 1 bndmaxslack(i,j) ensure maxslack is zero if use variable is 0; xminship(i,j).. x(i,j) =g= minshipping - minslack(i,j); xmaxship(i,j).. x(i,j) =e= 0 + maxslack(i,j); bndminslack(i,j).. minslack(i,j) =l= bigM*(1-use(i,j)); bndmaxslack(i,j).. maxslack(i,j) =l= bigM*use(i,j); Model indicatorbigMModel /cost, supply, demand, xminship, xmaxship, bndminslack, bndmaxslack/ ; * Let's first solve this without use of indicators indicatorbigMModel.optfile = 0; Solve indicatorbigMModel using mip minimizing z ; * Now we will use indicators and therefore we don't need the slacks putclose fcpx 'indic xminship(i,j)$use(i,j) 1' / 'indic xmaxship(i,j)$use(i,j) 0'; minslack.fx(i,j) = 0; maxslack.fx(i,j) = 0; indicatorbigMModel.optfile = 1; Solve indicatorbigMModel using mip minimizing z ; * We can also mix and match bigM with indicator constraints putclose fcpx 'indic xminship(i,j)$use(i,j) 1';

minslack.fx(i,j) = 0;   maxslack.up(i,j) = inf;
indicatorbigMModel.optfile = 1; Solve indicatorbigMModel using mip minimizing z ;