This file shows how to solve Example 4.4.1 (Drug Running) from Van Roy and Maison in Pyomo.

In [1]:
from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")

Given are the probabilities of interception for the various ports.

In [2]:
IP = {('SD'):1/3,
      ('LA'):1/2,
      ('SF'):3/5}

We extract the list of ports as usual.

In [3]:
Cities = list(IP.keys())

Decision variables are the $p_i$, the probability of sending a boat to port $i$, and the overall probability of getting through without interception $s$.

In [4]:
model = ConcreteModel()
model.p = Var(Cities, within=NonNegativeReals)
model.s = Var(within=Reals)

Now the model. According to the minimax principle, the optimal solution is given by maximizing $s$ subject to $$ \sum_{i} p_i = 1 $$ and $$ s \leq \sum_{i} p_i \, s_{ij} $$ where $s_ij$ is the probability of success at port $i$ assuming the coast guard vessel is located at port $j$. Thus, $s_{ij}=1$ if $i \neq j$ and $s_{ii} = 1 - \text{IP}_i$.

In [5]:
def intercept_rule(model, j):
    return model.s <= sum(model.p[i] for i in Cities) - IP[j]*model.p[j]
    
model.intercept = Constraint(Cities, rule=intercept_rule)

model.probabilities = Constraint(expr = sum(model.p[i] for i in Cities) == 1)

model.success = Objective(expr = model.s, sense=maximize)

We solve and output the optimal landing probabilities.

In [6]:
model.dual = Suffix(direction=Suffix.IMPORT)
results = opt.solve(model)
model.p.get_values()
Out[6]:
{'LA': 0.3, 'SD': 0.45, 'SF': 0.25}

The dual variables are the corresponding probabilities that the coast guard should choose to guard the different ports. Since the pay-off matrix for this problem is symmetric, these dual probabilities come out identical to the primal one. In more general two-player zero-sum games, they may be different.