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')
ax = fig.add_subplot(111)

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