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.

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)
```

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)) + "%")
```

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()
```