Towers of Hanoi: inventory and network formulation

A standard problem with 4 disks requires 15 moves

The towers of Hanoi problem [1] is a famous puzzle demonstrating recursion. The task is to move a stack of disks from one peg to another. We can only move one disk at a time, and we need to obey the rule that larger disks can never be on top of a smaller disk. The moves for the 3 disk problem are:

Click on the picture to enlarge.

Wolf, goat and cabbage problem: MIP and network model

The Wolf, goat, and cabbage problem can be stated as [1]:

A farmer with a wolf, a goat, and a cabbage must cross a river by boat. The boat can carry only the farmer and a single item. If left unattended together, the wolf would eat the goat, or the goat would eat the cabbage. How can they cross the river without anything being eaten?

A long transcontinental flight was a good opportunity to try to attack this problem. Here are some approaches to model this. Not at all a very useful or practical model, but still interesting, I think (although I may be in a small minority here). 

Inventory model

In this model, we keep track of the inventory of items (wolf, goat, cabbage). We assume the starting inventory is: all items are on the left bank of the river. The final inventory should be: all items are on the right bank. In this model, I assume the following numbering scheme:

----     31 SET trip  trips

trip1 ,    trip2 ,    trip3 ,    trip4 ,    trip5 ,    trip6 ,    trip7 ,    trip8 ,    trip9 ,    trip10


----     31 SET dir  direction of trip

L->R,    R->L


----     31 SET tripDir  trip direction combos

              L->R        R->L

trip1          YES
trip2                      YES
trip3          YES
trip4                      YES
trip5          YES
trip6                      YES
trip7          YES
trip8                      YES
trip9          YES
trip10                     YES

Small MIP, proving optimality is difficult

This is a simple, smallish MIP model that is difficult to solve to proven optimality. The problem is stated as [1]:

I have grid of dimensions H and W of boolean variables. The only constraint is that if a variable is false then at least one of the adjacent variables (top, right, left, bottom, diagonals don't count) must be true. The goal is to minimize the number of true values in the grid.

A high-level mathematical representation of this can be:

High-level model
\[\begin{align}\min&\sum_{i,j}\color{darkred}x_{i,j}\\ & \color{darkred}x_{i,j}=0 \implies \color{darkred}x_{i-1,j} + \color{darkred}x_{i+1,j} + \color{darkred}x_{i,j-1} + \color{darkred}x_{i,j+1} \ge 1 \\ & \color{darkred}x_{i,j}\in \{0,1\} \end{align}\]

GAMS 48 tests

Some minor quibbles. 

gdx2sqlite

The latest version of GAMS contains a replacement of gdx2sqlite. This dumps a GDX file into a SQLite database. It is a tool I use a lot. Here is a comparison using the indus89 model in the GAMS model library:

Prevent Loops in GAMS

 

This book [1] on DEA models has an accompanying website with all the GAMS models [2]. 

Of course, I'll be doing some nitpicking on the GAMS code. 

Solving DEA Models with GAMS

Data Envelopment Analysis (DEA) models are somewhat special. They typically consist of small LPs, of which a whole bunch have to be solved. The CCR formulation (after [1]), for the \(i\)-th DMU (Decision Making Unit), can be stated as [2]:

CCR LP Model

\[\begin{align} \max \>& \color{darkred}{\mathit{efficiency}}_i=\sum_{\mathit{outp}} \color{darkred}u_{{\mathit{outp}}} \cdot \color{darkblue}y_{i,{\mathit{outp}}} \\ & \sum_{\mathit{inp}} \color{darkred}v_{{\mathit{inp}}} \cdot \color{darkblue}x_{i,{\mathit{inp}}} = 1 \\ & \sum_{\mathit{outp}} \color{darkred}u_{{\mathit{outp}}} \cdot \color{darkblue}y_{j,{\mathit{outp}}} \le \color{darkred}v_{{\mathit{inp}}} \cdot \color{darkblue}x_{j,{\mathit{inp}}} && \forall j \\ & \color{darkred}u_{{\mathit{outp}}} \ge 0, \color{darkred}v_{{\mathit{inp}}} \ge 0 \end{align}\]

Multiple Solutions in Minimum Spanning Tree example

In [1], I discussed some LP and MIP formulations for the Minimum Spanning Tree (MST) problem. 

Minimum Spanning Tree visualized through Google Maps

Here, I focus on two formulations: a multicommodity network approach (this can be solved as a large LP) and a MIP formulation based on techniques we know from the Traveling Salesman Problem (TSP). The main issue I want to discuss is the presence of multiple optimal solutions.

Circle Packing and HTML reporting

Little example. Here, we try to pack \(n\) circles with a given radius \(r_i\) into a larger disc with an unknown radius \(R\). The goal is to minimize \(R\). The underlying model is simple: 

Packing of Circles
\[\begin{align} \min\> & \color{darkred}R \\ & \sum_c \left(\color{darkred}p_{i,c}-\color{darkred}p_{j,c}\right)^2 \ge \left(\color{darkblue}r_i+\color{darkblue}r_j\right)^2 & \forall i\lt j \\ & \sum_c \color{darkred}p_{i,c}^2 \le \left(\color{darkred}R-\color{darkblue}r_i\right)^2 & \forall i \\ & \color{darkred}R \ge 0\\ & c \in \{x,y\} \\ \end{align}\]

Revised Simplex LP Solver written in GAMS

I am teaching some GAMS classes, and a question arose: "How does the Simplex method work?" It's not easy to answer in a few sentences, but I want to touch upon the concept of a basis anyway. Once you have a good intuition of what a basis is, a simple Simplex method is not so far-fetched. I find the tableau presentation somewhat confusing and far removed from what actual Simplex solvers do. I strongly prefer the Revised Simplex Method in matrix notation. 

Minor rant: I just don't understand the appeal of the tableau method. It looks to me like an invention for torturing undergrad students. Most of all, it is not very structure-revealing; it does not help you understand the underlying concepts. But about 100% of the LP textbooks insist we should learn that first.

As a gimmick, I implemented a simplified version in the GAMS language. This reminds me that someone spent the effort writing a Basic interpreter in TeX [1]. This is probably just as useful.

Small non-convex MINLP: Pyomo vs GAMS

 In [1], the following Pyomo model (Python fragment) is presented:

model.x = Var(name="Number of batches", domain=NonNegativeIntegers, initialize=10)                    
model.a = Var(name="Batch Size", domain=NonNegativeIntegers, bounds=(5,20))

# Objective function
def total_production(model):
    return model.x * model.a
model.total_production = Objective(rule=total_production, sense=minimize)

# Constraints
# Minimum production of the two output products
def first_material_constraint_rule(model):
    return sum(0.2 * model.a * i for i in range(1, value(model.x)+1)) >= 70
model.first_material_constraint = Constraint(rule=first_material_constraint_rule)

def second_material_constraint_rule(model):
    return sum(0.8 * model.a * i for i in range(1, value(model.x)+1)) >= 90
model.second_material_constraint = Constraint(rule=second_material_constraint_rule)

# At least one production run
def min_production_rule(model):
    return model.x >= 1
model.min_production = Constraint(rule=min_production_rule)

GAMS listing file: missing Unicode support

Newer versions of GAMS allow UTF-8 encoded strings as labels. That is very welcome, as these labels may come from data sources that just use Unicode characters. However, when printing to the listing file, we miss proper Unicode support. At first, I thought, "OK, just a few misaligned tables. No big deal." Here is a constructed example showing this may be a bit more problematic.

GAMS: SMAX and sparsity

This is a discussion about the SMAX function in GAMS and how it behaves for sparse data.

The data structure we were facing was something like:

set

i 'cases' /case1*case100000/

j 'attribute' /j1*j25/

k 'attribute' /k1*k25/

t 'type' /typ1*typ2/

;

parameter p(i,j,k,t) 'positive numbers';

* note: for each i we have only one (j,k)

Plotting NUTS-2 maps from GAMS

 NUTS-2 regions are statistical subnational regions (often provinces), mainly for the EU and UK [1]. 

NUTS hierarchy (from [1])

In [2] we can find mapping information in the form of Shapefiles[3] and related formats. I used the GeoJSON[4] format, and created a Python notebook script to extract a GAMS set from that file. The file is reproduced in the appendix below. The NUTS-2 codes form the set elements, and the name is stored as explanatory text. There is an option to generate Latin names instead of using the native alphabet. The Latin names are also inside the geojson file. E.g. we have: 

EL65  'Πελοπόννησος'

which is in the Greek alphabet. Using the Latin name, this would look like:

EL65  'Peloponnisos'

Some GAMS embedded Python notes

Here are two issues you may want to be aware of. The discussion below is relevant for Windows and not so much for Unix variants.

Using Python Raw strings for directories

Here we use some Python inside an otherwise empty GAMS model:

$onEmbeddedCode Python: dir = r"%system.fp%" print(dir) $offEmbeddedCode

Some approaches for moving data between MS Access and GAMS

Moving data between different environments is always more difficult than we hope. Here I list some approaches and actually try them out on a small dataset. We hit some bugs along the way and also a few conceptual stumbling blocks (even for this stylized example). We had some issues with Access as well as GAMS Connect. 

This question came up in an actual project. My idea was: "Let me show you how this can be done". I am afraid, I got carried away a bit. But it demonstrates that we should not underestimate these, at first sight, menial tasks. When the data set becomes larger, the problems compound. We can't eyeball the data, and statistically, it is more likely we encounter some problems.

Export GAMS GDX file to different Python formats (CSV,Feather,Pickle)

A GAMS GDX file is a binary file with a number of GAMS symbols with their data. When collaborating with colleagues that work with Python, it may be useful to convert a GAMS GDX with relevant data to data files that Python can consume. 

The idea is here that a GAMS user can run the GAMS script, without knowing Python or having Python installed. (GAMS comes with its own somewhat barebones Python, which is used inside the GAMS script). On the other hand, the Python user will not need to have GAMS installed to read the generated data files. This is opposed to Python packages and tools that can interact with GDX files directly. You can see a few in the PyPI directory [1]. Those will require both a GAMS system and a Python environment.

A network model: Pyomo vs GAMS

Network models are an important class of optimization models, in itself but also as part of larger models. In [1] a nice presentation is given on how to set up a network model in the Python-based modeling tool Pyomo.

 

The simple shortest path problem is reproduced here:

GAMS Singleton Set Issue

There is a nasty problem with the singleton set in GAMS. A singleton set is a set that can contain only zero or one element. Adding more elements implies replacement. I don't use it very often. Here is a reason why. 

Consider the small example:

set i /i1*i20/; parameter p(i); p(i) = uniform(0,10); * find last element >= 5 scalar threshold /5/; singleton set k(i); loop(i,   k(i)$(p(i)>=threshold) = yes; ); display p,k;

I would expect to see as output for k the last element in p(i) that is larger than 5. However, we see:

----     12 PARAMETER p  

i1  1.717,    i2  8.433,    i3  5.504,    i4  3.011,    i5  2.922,    i6  2.241,    i7  3.498,    i8  8.563
i9  0.671,    i10 5.002,    i11 9.981,    i12 5.787,    i13 9.911,    i14 7.623,    i15 1.307,    i16 6.397
i17 1.595,    i18 2.501,    i19 6.689,    i20 4.354


----     12 SET k  

                                                      ( EMPTY )

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. 

Eight soldiers lining up: a very large permutation problem

In [1] an intriguing problem is posted:

There are 8 soldiers, gathering and lining up every morning for their military service. The commander at the head of these soldiers demands that the morning lineup of these soldiers be arranged differently for every next day according to the following rule:

        Any three soldiers cannot be lined up next to each other in the same order for others days.

For example; If ABCDEFGH is the first arrangement for day 1, on the other days, ABC,BCD, CDE, DEF, EFG and FGH cannot be lined up next to each other in the same order any more, but ACB arrangement is okay for other days until used once since they are not in the same order. 

What is the maximum number of days can this happen?

 

Let's first make the problem a bit smaller. Say we have just 4 soldiers. This gives us  \(4!=24\) permutations or line-ups. Let's have a look at the following sets: 

• \(p\) is the set of permutations
• \(t\) is the set of substrings of length 3 that we can form from \(p\)
• \(\color{darkblue}pt(p,t)\) is the mapping between \(p\) and \(t\)

----     53 SET p  permutations or line ups

ABCD,    ABDC,    ACBD,    ACDB,    ADBC,    ADCB,    BACD,    BADC,    BCAD,    BCDA,    BDAC,    BDCA,    CABD
CADB,    CBAD,    CBDA,    CDAB,    CDBA,    DABC,    DACB,    DBAC,    DBCA,    DCAB,    DCBA


----     53 SET t  substrings of length 3

ABC,    BCD,    ABD,    BDC,    ACB,    CBD,    ACD,    CDB,    ADB,    DBC,    ADC,    DCB,    BAC,    BAD,    BCA
CAD,    CDA,    BDA,    DAC,    DCA,    CAB,    CBA,    DAB,    DBA


----     55 SET pt  mapping (computed in Python)

ABCD.ABC,    ABCD.BCD,    ABDC.ABD,    ABDC.BDC,    ACBD.ACB,    ACBD.CBD,    ACDB.ACD
ACDB.CDB,    ADBC.ADB,    ADBC.DBC,    ADCB.ADC,    ADCB.DCB,    BACD.ACD,    BACD.BAC
BADC.ADC,    BADC.BAD,    BCAD.BCA,    BCAD.CAD,    BCDA.BCD,    BCDA.CDA,    BDAC.BDA
BDAC.DAC,    BDCA.BDC,    BDCA.DCA,    CABD.ABD,    CABD.CAB,    CADB.ADB,    CADB.CAD
CBAD.BAD,    CBAD.CBA,    CBDA.CBD,    CBDA.BDA,    CDAB.CDA,    CDAB.DAB,    CDBA.CDB
CDBA.DBA,    DABC.ABC,    DABC.DAB,    DACB.ACB,    DACB.DAC,    DBAC.BAC,    DBAC.DBA
DBCA.DBC,    DBCA.BCA,    DCAB.DCA,    DCAB.CAB,    DCBA.DCB,    DCBA.CBA


----     66 PARAMETER np                   =           24  card(p)
            PARAMETER nt                   =           24  card(t)
            PARAMETER npt                  =           48  card(pt)

In [1], Rob Pratt proposes a MIP model for this:

Model A
\[\begin{align}\max\> & \color{darkred}z = \sum_p \color{darkred}x_p \\ &\sum_{p|\color{darkblue}{pt}(p,t)} \color{darkred}x_p \le 1&& \forall t \\ & \color{darkred}x_p \in \{0,1\}\end{align} \]