Friday, October 28, 2016

GAMS Modeling: Inventing Set Element Labels

I often see models with set elements 1,2,…. Here is an example:

image

Of course if we have real data often we can use real names. But in some cases we need to use numbers as IDs. This is also typically the case when we don’t have real data and have to invent things. BTW, inventing data for a model is actually a very useful exercise: it will force you to think about the model in a more structured way.

If I have to invent labels I tend not to use 1,2,3,… but rather prefix the number by a string indicating what it is, e.g.:

set
  i
/i1*i10/
  j
/j1*j10/
  w
'warehouse' /wh1*wh12/

;

The reason to use these longer labels is that the output of the model becomes easier to interpret. Below is an example:

Bare numbers Prefixed numbers
set
  p
'product' /1*10/
  f
'factory' /1*4/
  w
'warehouse' /1*3/
  t
'time period' /1*5/

;

positive variables x(p,f,w,t) 'product flow from factory to warehouse';

Declarations
 

set
  p
'product' /pr1*pr10/
  f
'factory' /fact1*fact4/
  w
'warehouse' /wh1*wh3/
  t
'time period' /t1*t5/

;

positive variables x(p,f,w,t) 'product flow from factory to warehouse';

---- VAR x  product flow from factory to warehouse

              LOWER        LEVEL        UPPER       MARGINAL

1 .1.1.1        .          20.0000      +INF           .         
1 .1.1.2        .          85.0000      +INF           .         
1 .1.1.3        .          60.0000      +INF           .         
1 .1.1.4        .          35.0000      +INF           .         
1 .1.1.5        .          30.0000      +INF           .         
1 .1.2.1        .          25.0000      +INF           .         
1 .1.2.2        .          35.0000      +INF           .         
1 .1.2.3        .          90.0000      +INF           .         
1 .1.2.4        .          10.0000      +INF           .         
1 .1.2.5        .          55.0000      +INF           .         
1 .1.3.1        .         100.0000      +INF           .         
...

First part of the solution listing. Without recalling the declaration of variable x this is difficult to understand.

 

Using the prefixes we do not have to remember the ordering of the indices for this variable

---- VAR x  product flow from factory to warehouse

                      LOWER       LEVEL       UPPER      MARGINAL

pr1 .fact1.wh1.t1       .         20.0000     +INF          .         
pr1 .fact1.wh1.t2       .         85.0000     +INF          .         
pr1 .fact1.wh1.t3       .         60.0000     +INF          .         
pr1 .fact1.wh1.t4       .         35.0000     +INF          .         
pr1 .fact1.wh1.t5       .         30.0000     +INF          .         
pr1 .fact1.wh2.t1       .         25.0000     +INF          .         
pr1 .fact1.wh2.t2       .         35.0000     +INF          .         
pr1 .fact1.wh2.t3       .         90.0000     +INF          .         
pr1 .fact1.wh2.t4       .         10.0000     +INF          .         
pr1 .fact1.wh2.t5       .         55.0000     +INF          .         
pr1 .fact1.wh3.t1       .        100.0000     +INF          .
         
...

----     26 VARIABLE x.L  product flow from factory to warehouse

INDEX 1 = 1

              1           2           3           4           5

1.1      20.000      85.000      60.000      35.000      30.000
1.2      25.000      35.000      90.000      10.000      55.000
1.3     100.000      60.000     100.000      80.000      15.000
2.1      65.000      20.000      30.000      70.000      45.000
2.2      40.000      40.000      15.000      20.000      60.000
2.3      85.000      25.000      70.000      80.000      35.000
3.1      15.000      55.000      20.000      90.000      30.000
3.2      30.000      60.000      75.000      65.000      50.000
3.3      45.000      15.000      35.000       5.000      35.000
4.1      20.000      65.000      60.000      80.000      30.000
4.2      70.000      80.000      65.000      30.000      10.000
4.3      15.000      65.000      55.000       5.000      80.000

INDEX 1 = 2

              1           2           3           4           5

1.1      10.000      20.000      55.000      80.000      20.000
1.2       5.000      60.000      65.000      40.000      40.000
1.3      25.000      25.000      15.000      95.000      40.000
2.1      80.000      35.000      15.000      75.000      10.000

...

Same with a default display output. Again only partially reproduced here.
 

----     24 VARIABLE x.L  product flow from factory to warehouse

INDEX 1 = pr1

                   t1          t2          t3          t4          t5

fact1.wh1      20.000      85.000      60.000      35.000      30.000
fact1.wh2      25.000      35.000      90.000      10.000      55.000
fact1.wh3     100.000      60.000     100.000      80.000      15.000
fact2.wh1      65.000      20.000      30.000      70.000      45.000
fact2.wh2      40.000      40.000      15.000      20.000      60.000
fact2.wh3      85.000      25.000      70.000      80.000      35.000
fact3.wh1      15.000      55.000      20.000      90.000      30.000
fact3.wh2      30.000      60.000      75.000      65.000      50.000
fact3.wh3      45.000      15.000      35.000       5.000      35.000
fact4.wh1      20.000      65.000      60.000      80.000      30.000
fact4.wh2      70.000      80.000      65.000      30.000      10.000
fact4.wh3      15.000      65.000      55.000       5.000      80.000

INDEX 1 = pr2

                   t1          t2          t3          t4          t5

fact1.wh1      10.000      20.000      55.000      80.000      20.000
fact1.wh2       5.000      60.000      65.000      40.000      40.000
fact1.wh3      25.000      25.000      15.000      95.000      40.000
fact2.wh1      80.000      35.000      15.000      75.000      10.000
...

image With the GDX viewer things can become even more mysterious and difficult to decipher as we can flip dimensions around. You get a little bit of help with this:image
The 5th dimension is LO,L(evel),M,UP. We see we flipped dimensions 4 and 3 (t and w).
Here we see immediately what is going on. image

Notes

  1. It makes sense to use this approach even with real (not invented) numeric IDs. Once the data arrives in GAMS it is often not so easy to “repair” set element labels. If the data is sitting in a database it can be helpful to prefix numeric IDs in the SQL query:

    SELECT 'pr' || productid
    FROM
    table


    For a spreadsheet, often the easiest way is to insert a column or create a copy of the table with the new prefixed names.
  2. The construct i.val does not work with numbers that are prefixed. I usually work around this to create a parameter val(i) that contains the value of element i. If needed I create a superset of i that has an easy way to populate val(s). E.g.

    set
       s
    /i3*i10/
       i(s)
    /i3*i5,i8*i10/
    ;
    parameter val(s);
    val(s) =
    ord
    (s)+2;
    * now we can use val(i)

  3. It is not always possible to distinguish dimensions by their element names. E.g.

    set i 'nodes' /node1*node10/;
    alias
    (i,j);
    positive variable
    flow(i,j);


    Here, by nature, the two indices i and j have the same names because they have the same domain.

Thursday, October 27, 2016

Gurobi 7

Gurobi Optimizer 7.0 Now Available

This release offers significantly higher performance and powerful new modeling capabilities, among other enhancements. Specifically, v7.0 includes:
  • Significant performance improvements on real-world models
    • MIP: 19% faster overall and 30% faster on models that take >100s to solve
    • LP (concurrent): 10% faster overall and 25% faster on models that take >100s to solve
    • Significant QCP, MIQCP, and MIQP performance improvements
  • Python modeling enhancements — More easily translate mathematical models into code
  • Support for multiple objectives — Specify multiple objectives and their relative priorities
  • MIP solution pool support — Obtain more than just one optimal solution to a MIP model
  • New general constraints — Directly use MIN/MAX, ABS, AND/OR, and Indicator constraints
  • Expanded Python support for the Mac — Gurobi now supports Python 3.5 on the Mac
Additional enhancements include an even more flexible tuning tool, useful new parameters, and enhanced .NET property support. You can learn more on our What's new in Gurobi 7.0 page.

Significantly Enhanced Gurobi Instant Cloud

The Gurobi Instant Cloud has undergone a significant set of upgrades. These include enhanced API support for easier launching of cloud instances directly from Gurobi APIs, a more intuitive online interface to help users get up and running more quickly, and new machine pool support which simplifies the launching and management of cloud machines.
You can learn more about these enhancements, and our new cloud subscription plans, on our Cloud Overview page.

Sunday, October 23, 2016

Optimization of p*q: price and quantity

When doing revenue optimization we need to take into account the price elasticity of demand. Not sure if this publisher has a good handle on that.

image

 

References

  1. George Bittlingmayer, “The Elasticity of Demand for Books, Resale Price Maintenance and the Lerner Index“,  Journal of Institutional and Theoretical Economics (JITE) / Zeitschrift für die gesamte Staatswissenschaft, Vol. 148, No. 4 (December 1992), pp. 588-606
  2. http://yetanothermathprogrammingconsultant.blogspot.com/2016/09/getting-rid-of-fly-ash_12.html: another example where we use a price elasticity to optimize a firm’s profit.

Friday, October 21, 2016

ROML: R Optimization Modeling Language

From: https://rdrr.io/rforge/ROML/:

Utilizing the possibilities of R's lazy evaluation ROML provides an algebraic optimization modeling language. ROML is still under very active developments, so it's possible that there will be heavy changes and some parts are not so well tested as they should be.

Keep track of this.

References:

  1. http://yetanothermathprogrammingconsultant.blogspot.com/2016/09/lpmip-models-in-r.html

Thursday, October 20, 2016

Excel and very large data files

Tables with many rows (much more than the one million row limit in Excel) can be loaded in Excel’s Power Pivot Data Model. Here is an example from Statistics New Zealand with a CSV file of more than 23 million records.

image 

Actually Excel stores this data quite efficiently:

image

Note: wc reports one extra row because of the headers in line 1.

References
  1. http://yetanothermathprogrammingconsultant.blogspot.com/2016/07/excel-pivot-table-and-large-data-sets.html

Wednesday, October 19, 2016

Microsoft Open Sources Gradient Boosting Machine code

Microsoft abandoned their MS Solver Foundation framework, now focusing on machine learning (and R) ?

References

  1. Alexey Natekin, Alois Knoll, ”Gradient Boosting Machines, A Tutorial”, 1983, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3885826/pdf/fnbot-07-00021.pdf

Commercial Vehicle Routing

Source

An interesting approach to provide vehicle routing services is: https://routific.com/. You send a spreadsheet with the data (addresses) and get the ordering back,

  • The CEO has a degree from Erasmus University in Rotterdam
  • They got some funding to get started ($1.2M).
  • Article in Fast Company: https://www.fastcoexist.com/3062836/world-changing-ideas/this-bee-inspired-algorithm-helps-delivery-companies-plan-the-most-effi
  • One more article: http://www.digitaltrends.com/cool-tech/routific-bee-algorithm-delivery/
  • This is based on a “Bee Inspired” heuristic. There are quite a few of them:
    artificial bee colony (ABC) algorithm, honeybees mating optimization (HBMO) algorithm, artificial beehive algorithm (ABHA), bee colony optimization (BCO) algorithm, bee colony inspired algorithm (BCiA), bee swarm optimization (BSO) algorithm, bee system (BS) algorithm, BeeHive algorithm, bees algorithm (BeA), bees life algorithm (BLA), bumblebees algorithm, honeybee social foraging (HBSF) algorithm, OptBees algorithm, simulated bee colony (SBC) algorithm, virtual bees algorithm (VBA), and wasp swarm optimization (WSO) algorithm.  This list is from (1)
  • Looks like simple heuristics are preferred, compared to sophisticated mathematical programming based algorithms.

References

  1. Bo Xing, Wen-Jing Gao, “Bee Inspired Algorithms”, in Innovative Computational Intelligence: A Rough Guide to 134 Clever Algorithms, 2013.

Tuesday, October 18, 2016

MIP Modeling: if constraint

From http://cs.stackexchange.com/questions/51025/cast-to-boolean-for-integer-linear-programming/51089#51089

could you please explain the your logic with binary variable delta, for more general case, for instance, if y={a: x>0, b: x<0}.

If we assume \(y\) can be either \(a\) or \(b\) (constants) when \(x=0\), and \(x \in [L,U]\)  (with \(L<0\) and \(U>0\)) then we can write

\[\begin{matrix}
y = \begin{cases}a&x\ge 0\\b&x\le 0 \end{cases}
& \Longrightarrow
& \boxed{\begin{align}&L(1-\delta) \le x \le U \delta\\&y=a\delta + b (1-\delta)\\ &\delta \in \{0,1\} \end{align}}
\end{matrix} \]

I am not sure if the left part is formulated mathematically correctly. May be I should say:

\[y = \begin{cases} a&x> 0\\
                   b&x< 0 \\
                   a \text{ or } b & x=0
\end{cases}\]

If we want to exclude \(x=0\) then we need to add some small numbers:

\[\begin{matrix}
y = \begin{cases}a&x> 0\\b&x < 0 \end{cases}
& \Longrightarrow
& \boxed{\begin{align}&\varepsilon + (L-\varepsilon)(1-\delta) \le x \le –\varepsilon +(U+\varepsilon)\delta\\&y=a\delta + b (1-\delta)\\&\delta \in \{0,1\}\\
&\varepsilon=0.001\end{align}}
\end{matrix} \]
Of course in many cases we can approximate \(U\approx U+\varepsilon\) (similar for \(L\)), I was a bit precise here.

If \(a\) and \(b\) are variables (instead of constants) things become somewhat more complicated. I believe the following is correct:

\[\boxed{\begin{align}
&\varepsilon + (L-\varepsilon)(1-\delta) \le x \le –\varepsilon +(U+\varepsilon)\delta\\
&a+M_1(1-\delta)\le y \le a+ M_2(1-\delta)\\
&b+M_3\delta\le y \le b+M_4\delta\\
&\delta \in \{0,1\}\\
&\varepsilon=0.001\\
&M_1 = b^{LO}-a^{UP}\\
&M_2 = b^{UP}-a^{LO}\\
&M_3 = a^{LO}-b^{UP}\\
&M_4 = a^{UP}-b^{LO}
\end{align}}\]

I used here \(a \in [a^{LO},a^{UP}]\) and \(b \in [b^{LO},b^{UP}]\).