Thursday, April 29, 2021

What is this BRATIO option in GAMS?

This is a note about BRATIO, a very exotic GAMS option related to providing an advanced basis to a Simplex-based LP (or NLP) solver. 


Marginals


For every variable and equation used in a model, GAMS stores the level (.L) and marginal (.M). The level is just its value. I.e. to print results, we can write

display x.l(i,j);

The marginals correspond to duals for equations and reduced costs for variables. Printing is similar:

display capacity.m(wh); 

The levels and marginals are typically populated (or updated) by the solver. If not set, their default values are zero. This makes sense as GAMS stores everything in sparse data structures. That means a variable with default values will not use any memory. This is a property we use extensively in large models: only a small fraction of variables 

production(product,variant,facility,time)

maybe actually used e.g. because a facility can only produce certain products/variants.


Nonlinear models: initial values  


For nonlinear models, it is often important to provide a reasonable starting point. We can do this by assigning values:

x.l(i,j) = 100;

A special note: if the initial value is outside the bounds, GAMS will project the value to the closest bound. I.e. if the lower bound is 110, the actual initial value passed on to the solver will not be 100 but rather 110. 

Besides setting initial values, they can also come from a model solved earlier. This also means: solving two identical models in sequence will make things very easy for the second solve: 

option nlp=ipopt;
solve chain using nlp minimizing energy;
solve chain using nlp minimizing energy;

If we run this we see:

               S O L V E      S U M M A R Y 

     MODEL   chain               OBJECTIVE  energy
     TYPE    NLP                 DIRECTION  MINIMIZE
     SOLVER  IPOPT               FROM LINE  74

**** SOLVER STATUS     1 Normal Completion
**** MODEL STATUS      2 Locally Optimal
**** OBJECTIVE VALUE                5.0686

 RESOURCE USAGE, LIMIT         29.531 10000000000.000
 ITERATION COUNT, LIMIT       362    2147483647
 EVALUATION ERRORS              0             0

               S O L V E      S U M M A R Y

     MODEL   chain               OBJECTIVE  energy
     TYPE    NLP                 DIRECTION  MINIMIZE
     SOLVER  IPOPT               FROM LINE  75

**** SOLVER STATUS     1 Normal Completion
**** MODEL STATUS      2 Locally Optimal
**** OBJECTIVE VALUE                5.0686

 RESOURCE USAGE, LIMIT          0.078 10000000000.000
 ITERATION COUNT, LIMIT         0    2147483647
 EVALUATION ERRORS              0             0

We see the first solve takes 30 seconds and 362 iterations. The second solve took 0.1 seconds and 0 iterations. IPOPT immediately sees we are giving it an optimal solution.

This scheme also works when we mix solvers (e..g the first solve with CONOPT), or when we make some changes in the model.


Linear programming: advanced basis


LP solvers using a simplex method (either a primal or dual simplex algorithm), don't really use the level values, but rather basis information. This data tells the solver which variables and equations are basic (to be determined by the solver by solving a linear system) or non-basic (temporarily fixed at one of the bounds). GAMS will pass on, besides levels, marginal indicators for the rows and columns. They look like:  
  • x.m (or e.m) is zero: the variable or equation is basic
  • x.m (or e.m) is non-zero or EPS: the variable or equation is non-basic.
For an optimal or just feasible solution, all basic variables and equations are typically between their bounds (they can be at bound, which is called a degenerate solution). The non-basic variables and equations are typically pegged at their bound (or one of the bounds). A special case is a free variable. They are usually basic, but sometimes they can be non-basic (in that case they are usually zero). For more information on how to find out how the basis looks like, see [1].

 Detail: the levels may be used to determine if a non-basic variable/equation is "at-lower" or "at-upper". The levels are not used for basic variables.

A basic solution provided by a solver always has:
  • The number of basic variables and equations is equal to the number of equations.
  • The number of non-basic variables and equations is equal to the number of variables in the model.
Let's try again two solves in a row:  

option lp=cbc;
solve one minimizing cost using lp;
solve one minimizing cost using lp;               

Unfortunately, the listing file does not display the correct iteration count for the first solve:

               S O L V E      S U M M A R Y 

     MODEL   one                 OBJECTIVE  cost
     TYPE    LP                  DIRECTION  MINIMIZE
     SOLVER  CBC                 FROM LINE  1090 

**** SOLVER STATUS     1 Normal Completion
**** MODEL STATUS      1 Optimal
**** OBJECTIVE VALUE            27569.5989

 RESOURCE USAGE, LIMIT          0.033 10000000000.000
 ITERATION COUNT, LIMIT         0    2147483647


This seems to be a bug. We can look at the log instead:

Calling CBC main solution routine...
0  Obj 1894.9246 Primal inf 5281.3095 (85) Dual inf 1497.1161 (26)
68  Obj 4884.1415 Primal inf 9801.2782 (117)
120  Obj 11984.03 Primal inf 7311.583 (114)
193  Obj 18810.234 Primal inf 1723.8549 (97)
250  Obj 24693.927 Primal inf 2654.5824 (76)
323  Obj 27569.599
Optimal - objective value 27569.599
Optimal objective 27569.5989 - 323 iterations time 0.032, Presolve 0.01

Calling CBC main solution routine...
Optimal - objective value 27569.599
Optimal objective 27569.5989 - 0 iterations time 0.012, Presolve 0.01


CBC has the most crowded and ugly log, but at least it shows the last iteration. The first solve takes 323 iterations while the restart takes 0 iterations.

When we use Cplex (I used 1 thread which means dual simplex is used), things are a bit more complicated. Again we use:

option lp=cplex
solve one minimizing cost using lp;
solve one minimizing cost using lp;               

The first solve uses 250 iterations, and the second solve 2 iterations. The reason is that Cplex option advind=2 is used for the second solve. This uses the presolve which messes things up just a little bit so that the solver needs to do 2 (degenerate) iterations. If we use advind=1 (no presolve), we see better behavior: 0 iterations.

Use cases


Here we used two identical models in a row, which in practice is not so useful. A more practical use case is to solve a base case followed by a number of scenarios that change the model a little bit. The changes can be a combination of

  • changing some coefficients in the model
  • adding some constraints
  • deleting some constraints
and also
  • solving different models that share a portion of the equations and variables
  • setting the marginals directly [2,3]
As GAMS passes on the primal and dual information when available, we can even change the model type. The first model can be an LP, while we have some non-linearities to the second model.

Accept or reject a basis 


One question that needs to be addressed is: when do we decide a previous model is too different from the current one so that using an advanced basis does not make sense anymore? GAMS uses a somewhat unique, heuristic approach. Here comes the BRATIO option in play. Let \(m\) be the number of equations in the model, and let \(nbr\) the number of nonzero marginals for these equations (that is the number of non-basic rows), then we accept the basis if: \[nbr \gt \mathit{BRATIO} \times m\] The default value of BRATIO is 0.25. This makes some sense: if the number of basic rows is too large (with zero row marginals), we likely have changed the model too much (remember that the default value for a marginal is zero). 

I usually leave this default in place, with a few exceptions. We have the following special cases:
  • option BRATIO=1; means: ignore the basis and let the solver start from scratch. Hopefully, the solver can "crash" a good basis [4]. 
  • option BRATIO=0; indicates: always accept the basis, and pass it on to the solver [2,3]. 
The GAMS documentation states:

The value specified for bRatio will cause a basis to be discarded if the number of basic variables is smaller than bRatio times the number of equations.

This description is somewhat incorrect: the BRATIO option only looks at row marginals. For a proper basis (number of basics equal to the number of rows), we have that the number of non-basic rows is equal to the number of basic columns. However, for an incomplete or overspecified basis, this is not the case.

Notes

  • When solving multiple large LP scenarios, it may make sense to use an interior point (barrier) algorithm for the first solve, followed by a number of simplex solves.
  • Similarly, if the scenarios are NLPs, it may help to use an interior point solver for the first scenario or base case (e.g. using IPOPT or KNITRO), followed by an active set method for the other scenarios (e.g. using CONOPT or SNOPT).
  • I usually don't worry about BRATIO, but in some cases, I set it to zero (accept basis) or one (ignore basis).
  • For MIP models, BRATIO is usually not applicable or useful. Most MIP models spend a small fraction in the first LP and a MIP presolve will likely destroy a basis.
  • BRATIO is very unique to GAMS. To my knowledge no other tool has something like this.

References


  1. What is this EPS in a GAMS solution? https://yetanothermathprogrammingconsultant.blogspot.com/2017/09/what-is-this-eps-in-gams-solution.html
  2. Solving sparse systems of equations with an optimizer: a basis trick, https://yetanothermathprogrammingconsultant.blogspot.com/2018/02/solving-sparse-systems-of-equations.html
  3. Inverting a matrix by LP, https://yetanothermathprogrammingconsultant.blogspot.com/2021/04/inverting-matrix-by-lp.html
  4. B.A. Murtagh, Advanced Linear Programming (McGraw-Hill, New York, 1981).

No comments:

Post a Comment