This example was developed by Fabricio Oliveira at 13/03/2016. Main reference is: Bertsimas, D. and Sim, M., 2004. The price of robustness. Operations research, 52(1), pp.35-53.

# Deterministic Model¶

Let us consider a knapsack problem. In this problem, we have to make a selection of a given number of items that are available. Each item $j$ has a value $c_j$ associated to it and a weigth $b_j$. We want to make an optimal selection that has maximised value, given that we have a limited capacity $b$.

Let $x_j$ be the binary variable representing the decision of selecting $(x_j = 1)$ or not a given item. Thus, our optimisation problem can be posed as follows.

\begin{align} \max_x &\sum_j c_jx_j \\ \text{s.t.: } &\sum_j a_jx_j \leq b\\ &x_j \in \lbrace 0, 1\rbrace \end{align}

A bit of coding to solve the deterministic version of this problem would look like this...

In [1]:
import pulp as pp
import numpy as n
import random as r
import matplotlib.pyplot as plot

#Some initial declarations
items = range(0,50)
c = n.zeros(50)
a = n.zeros(50)
b = 500

#Creating instances randomly
for j in items:
c[j] = r.randrange(10,200)
a[j] = r.randrange(10,100)

#Creating the problem instance
detProb = pp.LpProblem("Knapsack Deterministic", pp.LpMaximize)

#Declaring our x variables
x = pp.LpVariable.dicts("x", items, 0 , 1, 'Integer')

#Including OF
detProb += sum(c[j]*x[j] for j in items)

#Including constraint
detProb += sum(a[j]*x[j] for j in items) <= b

#Solve the model
detProb.solve()

#Printing out the solution
print("Optimal selection: ")
selection = ""
for j in items:
if x[j].value() > 0:
if selection == "":
selection += str(j)
else:
selection += ", " + str(j)
print(selection)


Optimal selection:
3, 7, 9, 11, 19, 22, 23, 25, 32, 38, 45, 47, 48, 49


Now, what happens if we simulate that selection? In other words, what are the chances that we end up having a knapsack that is actually heavier than my weight limit $b$?

Suppose that the weights $a_j$ follow a Gaussian distribution with average $\overline{a}_j$ and deviation $0.5\overline{a}$. A simple Monte Carlo simulation show us the following.

In [2]:
simulations = range(10000)
simulResult = n.zeros(len(simulations))
aTilde = n.zeros(len(items))

for s in simulations:
for j in items:
aTilde[j] = r.gauss(a[j], 0.5*a[j]) #sampling an observed \tilde{a}

if sum(aTilde[j]*x[j].value() for j in items) <= b: #registering if the experiment is feasible
simulResult[s] =  1

print("Chance of feasibility.: " + str(sum(simulResult)*100/len(simulations)) + "%")

Chance of feasibility.: 50.51%


# Robust Model¶

Using the RO philosophy, let's assume that the weights $a_j$ are uncertain and take value in an interval which is defined as $\left[\overline{a}_j - \hat{a}_j, \overline{a}_j + \hat{a}_j \right]$. Let us set $\hat{a} = 0.5\overline{a}$. By doing so, it lead us to following Robust Counterpart

\begin{align} \max_x &\sum_j c_jx_j \\ \text{s.t.: } &\sum_j a_jx_j + \Gamma\pi + \sum_j p_j \leq b\\ &\pi + p_j \geq \hat{a}x_j, \ \forall j\\ &x_j \in \lbrace 0, 1\rbrace \\ &p_j \geq 0 ,\forall j, \pi \geq 0 \end{align}

which, once again, can be modelled as follows (note that we have embedded the model into a function so we can use in the following experiment.

In [3]:
def solveForGamma(Gamma):

#defining it for the robust formulation
aHat = 0.5*a

#Creating the problem instance
robProb = pp.LpProblem("Knapsack Robust", pp.LpMaximize)

#Declaring the selection variables
x = pp.LpVariable.dicts("x", items, 0 , 1, 'Integer')

#... and the additional dual variables
pi = pp.LpVariable("pi", 0 , 5000,'Continuous')
p = pp.LpVariable.dicts("p", items, 0 , 5000,'Continuous')

#Including OF
robProb += sum(c[j]*x[j] for j in items)

#Including constraints
robProb += sum(a[j]*x[j] for j in items) + Gamma*pi + sum(p[j] for j in items) <= b

for j in items:
robProb += pi + p[j] >= aHat[j]*x[j]

#Solve the model
robProb.solve()

#...and simulate on-the-fly (note to self: stop being lazy and use a simulation procedure as well)
simulations = range(10000)

#simulation elements
simulResult = n.zeros(len(simulations))
aTilde = n.zeros(len(items))

#Performing Monte Carlo simulation...
for s in simulations:
for j in items:
aTilde[j] = r.gauss(a[j], 0.5*a[j])

if sum(aTilde[j]*x[j].value() for j in items) <= b:
simulResult[s] =  1

obsProb = sum(simulResult)*100/len(simulations)
optimalValue = robProb.objective.value()

#Printing the output
print("Optimal selection: ")
selection = ""
for j in items:
if x[j].value() > 0:
if selection == "":
selection += str(j)
else:
selection += ", " + str(j)
print(selection)

return [obsProb, optimalValue]


Let's do a similar experiment, but now different values of $\Gamma$ to see what happens as we increase it.

In [4]:
totalGammas = 15
obsProb = n.zeros(totalGammas)
optimalValue = n.zeros(totalGammas)

#Let's repeat the experiment for several Gamma values and store observed probabilities and costs
for gamma in range(totalGammas):
print("Done for Gamma = " + str(gamma))
[obsProb[gamma], optimalValue[gamma]] = solveForGamma(gamma)

#Boring bit to get a pretty graph...
%matplotlib inline
fig = plot.figure()
fig.set_size_inches(18.5, 10.5, forward=True)
plot.title('The price of robustness')

def theoryBound(gamma): #this is the theoretical bound from Li et al. 2012, thm. 3.1.
return 100*(1 - n.exp(-(gamma**2)/(2*len(items))))

gammaRange = n.arange(0, totalGammas, 1)

#Plot on regular axis
ax.plot(obsProb, '-', label = 'Chance of feas.')
ax.plot(theoryBound(gammaRange), '-g', label = 'Theoretical bound.')

#Plot on the reverse axis
ax2 = ax.twinx()
ax2.plot(optimalValue, '-r', label = 'Obj. Function Values')

#get labels for legend
lines, labels = ax.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2, loc='center right')

ax.grid()
ax.set_xlabel("Gamma")
ax.set_ylabel("Probabilities")
ax2.set_ylabel("Obj. Function Values")

plot.show()


Done for Gamma = 0
Optimal selection:
3, 7, 9, 11, 19, 22, 23, 25, 32, 38, 45, 47, 48, 49
Done for Gamma = 1
Optimal selection:
3, 7, 9, 11, 14, 22, 23, 25, 32, 38, 40, 45, 47, 48, 49
Done for Gamma = 2
Optimal selection:
3, 7, 9, 11, 13, 14, 23, 25, 31, 32, 45, 47, 48, 49
Done for Gamma = 3
Optimal selection:
3, 7, 9, 11, 13, 22, 23, 25, 32, 40, 45, 47, 48, 49
Done for Gamma = 4
Optimal selection:
3, 7, 9, 11, 13, 14, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 5
Optimal selection:
3, 7, 9, 11, 14, 22, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 6
Optimal selection:
3, 7, 9, 11, 22, 23, 25, 32, 40, 45, 47, 48, 49
Done for Gamma = 7
Optimal selection:
3, 7, 9, 11, 22, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 8
Optimal selection:
7, 9, 11, 14, 22, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 9
Optimal selection:
3, 7, 9, 11, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 10
Optimal selection:
3, 7, 9, 11, 14, 22, 23, 25, 32, 45, 48, 49
Done for Gamma = 11
Optimal selection:
7, 9, 11, 22, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 12
Optimal selection:
7, 9, 11, 22, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 13
Optimal selection:
7, 9, 11, 22, 23, 25, 32, 45, 47, 48, 49
Done for Gamma = 14
Optimal selection:
7, 9, 11, 22, 23, 25, 32, 45, 47, 48, 49