# Knapsack Cover Inequalities

Consider the 0-1 knapsack Set 

$K := \Big\{\{0, 1\}^{n} : \sum_{j = 1}^n w_j x_j \le b \Big\}, w_j > 0, \forall j \in N = \{1, 2, ..., n\}, b > 0$

**Definition (Cover)** A set $C \subseteq N$ is a *cover* if $\sum_{j \in C} w_j > b$. A cover is minimal if $\sum_{j \in C \backslash \{j\}} w_j \le b, \forall j \in C$. 

**Prop** If $C \subseteq N$ is a cover for $K$, then $\sum_{j \in C} \le \vert C \vert - 1$ is valid for $K$.

The alternate reformulation of the Knapsack set is in terms of cover inequalities. In this tutorial, we will see how to generate a minimal cover, how to lift it to make it facet defining, and solve the resulting Knapsack problem.



```python
# First we import some packages
import numpy as np
from gurobipy import *
import time
```

### Generating problem instance


```python
def generateData(n): # n is the number of items
    np.random.seed(30) # Setting a seed value
    a = np.random.randint(1000, size=(n)) # Generating value of items 
    w = np.random.randint(1000, size=(n))# Generating weight of items 
    b = np.sum(w)*0.40 # Total weight of Knapsack
    return a, w, b
```

### Simply solving the Knapsack problem using Gurobi


```python
def solveKnapsackUsingGurobi(n, a, w, b, verbose = 1): # no. of items, value of items, weight of items, and Knapsack weight
    m = Model() # Preparing the model
    x = {i : m.addVar(lb=0, vtype=GRB.BINARY) for i in range(n)} # whether the item should be included in the Knapsack
    m.addConstr(sum([w[j] * x[j] for j in range(n)]) <= b) # Total weight should be less than the Knapsack capacity
    m.setObjective(sum([a[j] * x[j] for j in range(n)]), sense = GRB.MAXIMIZE) # Maximize the total value of the Knapsack
    m.Params.OutputFlag = 0; m.optimize()
    if m.status == 2:  
        if verbose == 1:            
            print('*********************** Gurobi ***************************')
            print('Items selected: ', {k for k in x if x[k].x > 0.4})
            print('Total value: ', m.objVal)
            print('**********************************************************')
        return {k for k in x if x[k].x > 0.4}, m.objVal
    else:
        print('Infeasible model')
        return {}, -float("inf")    
```

### Let's start preparing the initial model


```python
def preparePartialModel(n, a, w, b):
    m = Model()
    x = {i : m.addVar(lb=0.0, ub = 1.0, vtype=GRB.CONTINUOUS, name = str(i)) for i in range(n)}
    m.addConstr(sum([w[j] * x[j] for j in range(n)]) <= b)
    c, zeta = generateCover(n, a, w, b, {k:0 for k in range(n)})      
        
    m.addConstr(sum([x[j] for j in c]) <= len(c)-1)
    m.setObjective(sum([a[j] * x[j] for j in range(n)]), sense = GRB.MAXIMIZE)
    m.Params.OutputFlag = 0; m.Params.IntegralityFocus = 1
    return m

```

## Generating a cover/Seperation problem
To use cover inequalities in a cutting plane scheme, one is faced with the separation problem, that is, given a vector $\bar{x} \in \left[0, 1\right]^n$, find a cover inequality for $K$ that is violated by $\bar{x}$, or show that none exists. For this purpose, we solve the following problem:

$\zeta = \min\{\sum_{j \in C} (1-\bar{x}_j): C \text{ is a cover for } K\}$

- If $\zeta \ge 1$, then $\bar{x}$ satisfies all the cover inequalities for $K$. 
- If $\zeta < 1$, then an optimal over to the above problem yields the violated minimal cover.

One can also write the above prblem as the following binary program. Let $z_j = 1$, if $j$ is in the violated cover $C$.

$$
	\begin{align}
	& \zeta = \underset{z}{\text{minimize}}
	& & \sum_{j=1}^n (1-\bar{x}_j) z_j \\
	& \text{subject to}
	& & \sum_{j= 1}^n w_j z_j \ge b + 1\\
	& & & z \in \{0,1\}^n
	\end{align}
$$


```python
def generateCover(n, a, w, b, xbar): # no. of items, value of items, weight of items, Knapsack weight, and xbar
     m1 = Model()
     z = {i:m1.addVar(lb=0.0, vtype=GRB.BINARY) for i in range(n)}
     m1.addConstr(sum([w[j] * z[j] for j in range(n)]) >= b + 1)
     m1.setObjective(sum([(1-xbar[j]) * z[j] for j in range(n)]), sense = GRB.MINIMIZE)
     m1.Params.OutputFlag = 0; m1.optimize()
     return [j for j in range(n) if z[j].x > 0.4], m1.objVal # return the cover and the value of zeta
```

## Lifting 

This is esentially a procedure to convert a minimal cover inequality $\sum_{j \in C} x_j \le \vert C \vert - 1$ into a facet defining inequality for $K$.

### Sequential Lifting
- Choose an ordering $j_1, j_2, ..., j_l$ of the indices in $N \backslash C$. Set t = 1.
- The valid inequality $\sum_{i = 1}^{t-1} \alpha_{j_i} x_{j_i} + \sum_{j \in C} x_j \le \vert C \vert - 1$  has been obtained so far.
- To calculate the largest value of $\alpha_{j_t}$ for which the inequality
$\alpha_{j_t} x_t  + \sum_{i = 1}^{t-1} \alpha_{j_i} x_{j_i} + \sum_{j \in C} x_j \le \vert C \vert -1$ is valid, solve the following Knapsack problem:

$$
	\begin{align}
	& \xi_t = \underset{x}{\text{maximize}}
	& & \sum_{i = 1}^{t-1} \alpha_{j_i} x_{j_i} + \sum_{j \in C} x_j \\
	& \text{subject to}
	& & \sum_{i = 1}^{t-1} \alpha_{j_i} x_{j_i} + \sum_{j \in C} w_j x_j \le b - w_{j_t}\\
	& & & x \in \{0,1\}^{\vert C \vert + t - 1}
	\end{align}
$$
- Set  $\alpha_{j_t} = \vert C \vert - 1 - \xi_t$.
- Stop if $t = r$.

We remark that different orderings of $N \backslash C$ may produce different lifted inequalities. Furthermore, not all possible liftings can be derived from the above procedure.


```python
def sequentialLifting(c, n, a, w, b): # Cover set, no. of items, value of items, weight of items, and Knapsack weight.
    alpha = {k:0.0 for k in range(n) if k not in c} # coefficients of x_j, j not in cover
    m2 = Model()
    x_dash = {i:m2.addVar(lb=0.0, vtype=GRB.BINARY) for i in c}
    toBechecked = list(alpha.keys()); checked = [] # choosing the ordering of indices one-by-one and finding its coeffcient
    if len(toBechecked) != 0:        
        i = toBechecked.pop(0); 
        x_dash[i] = m2.addVar(lb=0.0, vtype=GRB.BINARY) # including a variable.
        m2.update()
        tmpc = sum([x_dash[j] * w[j] for j in c]) + sum([alpha[j]*x_dash[j] for j in checked]) # constraint of already in the Knapsack
        tmpo = sum([x_dash[j]  for j in c]) + sum([alpha[j]*x_dash[j] for j in checked])
        m2.setObjective(tmpo, sense = GRB.MAXIMIZE)
        m2.addConstr(tmpc <= b - w[i])
        m2.Params.OutputFlag = 0; m2.optimize()       
        if len(c) - 1 - m2.objVal >= 1: # Checking if the current item can be included in the cover.
            alpha[i] = len(c) - 1 - m2.objVal
            checked.append(i)
    return alpha # return coefficients of x_j in the lifted inequality.

```

### Balas' lifiting of minimal cover inequalities
Let $C$ be a minimal cover for $K$, and let 
$\sum_{j \in C} x_j + \sum_{j \in N \backslash C} w_j x_j \le \vert C \vert - 1$ be a lifting of the cover inequality associated with $C$. Up to permuting the indices, assume that C = \{1, ..., t\} and $a_1 \ge a_2 \ge ... \ge a_t$. 

Let $\mu_0 =  0$, and $\mu_h = \sum_{l = 0}^h w_j,$ for $h = 1, ..., t$	. Let $\lambda = \mu_t - b$ (note that $\lambda > 0$).

If above inequality defines a facet of $conv(K)$, then the following hold for every $j \in N \backslash C$.
- If, for some $h$, $\mu_h ≤ w_j ≤ \mu_{h+1} − \lambda$, then $\alpha_j = h$.
- If, for some $h$, $\mu_{h+1} − \lambda < w_j < \mu_{h+1}$, then $h ≤ \alpha_j ≤ h + 1$ (i.e., select $\alpha_j = h+1$).


```python
def BalasLifting(c, n, a, w, b): # Cover set, no. of items, value of items, weight of items, and Knapsack weight.
    orderCoverWeightsList = sorted(c, key=lambda item: w[item], reverse =True)
    alpha = {k:0.0 for k in range(n) if k not in c} # coefficients of x_j, j not in cover
    mu = {0:0.0}   
    ind = 0; cum = 0
    for k in orderCoverWeightsList:
        if ind == 0:
            mu[ind+1] = w[k]
            ind += 1; cum += w[k]
        else:
            mu[ind+1] = w[k] + cum
            ind += 1; cum += w[k]
    lamb = mu[len(mu)-1] - b
    
    
    for j in alpha:
        h = [h for h in mu if h < len(mu)-1 and mu[h] <= w[j] <= mu[h+1] - lamb]
        if len(h) != 0:
            alpha[j] = h[0]
        else:
            h = [h for h in mu if h < len(mu)-1 and mu[h] -lamb <= w[j] <= mu[h]]
            alpha[j] = h[-1]
    return alpha # return coefficients of x_j in the lifted inequality.
```

### Solving Knapsack problem using cover inequalities
Putting above peices together.


```python
def solveKnapsackUsingCoverInequalities(n, a, w, b, lifting = 'sequential', verbose = 1):
    m = preparePartialModel(n, a, w, b)
    m.optimize()
    if m.status == 2:
        x = {i:m.getVarByName(str(i)).x for i in range(n)}   
        c, zeta = generateCover(n, a, w, b, x)        
    else:    
        print('Infeasible model')       
    alpha = []
   
    while zeta < 1:
        m.addConstr(sum([m.getVarByName(str(i)) for i in c]) + sum([m.getVarByName(str(i)) * alpha[i] for i in alpha]) <= len(c) - 1)
        m.optimize()
        x = {i:m.getVarByName(str(i)).x for i in range(n)}   
        c, zeta = generateCover(n, a, w, b, x) 
        #print('Cover', sum([w[k] for k in c])-min([w[k] for k in c]) <= b)
        if lifting == 'sequential':
            alpha = sequentialLifting(c, n, a, w, b)
        elif lifting == 'Balas':
            alpha = BalasLifting(c, n, a, w, b)
        elif lifting == 'None':
            alpha = []
        else:
            print('lifting is unknown.. terminating')
            return 'NA'
        zeta = round(zeta, 4)
        #print(c, zeta, alpha)
        
    if verbose == 1:
        print('******************* Cover Inequality *********************')
        print('Items selected: ', {i:m.getVarByName(str(i)).x for i in range(n) if m.getVarByName(str(i)).x > 0.4} )
        print('Total value: ', sum([a[k] for k in range(n) if round(m.getVarByName(str(k)).x) > 0.4]))
        print('**********************************************************')
    return {i:m.getVarByName(str(i)).x for i in range(n) if m.getVarByName(str(i)).x > 0.4}, sum([a[k] for k in range(n) if round(m.getVarByName(str(k)).x) > 0.4]), len(m.getConstrs()) - 1
             
```

## Solve the models and print results


```python
import pandas as pd
n = 1000 # No. of items
results = []; results.append(['Total items', n])
a, w, b = generateData(n)
start = time.time()
xg, objg = solveKnapsackUsingGurobi(n, a, w, b, verbose = 0)
#print('Gurobi model took', time.time()-start); start = time.time()
timeF = time.time()-start
results.append(['No. of items selected', len(xg)])
results.append(['Value of items selected', objg])
results.append(['No. of covers generated', float("NAN")])
results.append(['Computational time', timeF])

for lift in ['None', 'sequential', 'Balas']:    
    start = time.time()
    xc, objc, nc = solveKnapsackUsingCoverInequalities(n, a, w, b, lifting = lift, verbose = 0)
    timeF = time.time()-start
    #print('CG model took', time.time()- start)
    results[0].append(n)
    results[1].append(len(xc))
    results[2].append(objc)
    results[3].append(nc)
    results[4].append(timeF)

# Create the pandas DataFrame
pd.options.display.float_format = "{:,.2f}".format
df = pd.DataFrame(results, columns = ['Particulars', 'IP with Gurobi', 'Knapsack cover with no lifting', 'Knapsack cover with sequential lifting', 'Knapsack cover with Balas lifting'])
display(df)
```

    Academic license - for non-commercial use only - expires 2021-09-17
    Using license file C:\Users\Pramesh Kumar\gurobi.lic



<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Particulars</th>
      <th>IP with Gurobi</th>
      <th>Knapsack cover with no lifting</th>
      <th>Knapsack cover with sequential lifting</th>
      <th>Knapsack cover with Balas lifting</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>Total items</td>
      <td>1,000.00</td>
      <td>1,000.00</td>
      <td>1,000.00</td>
      <td>1,000.00</td>
    </tr>
    <tr>
      <th>1</th>
      <td>No. of items selected</td>
      <td>561.00</td>
      <td>560.00</td>
      <td>560.00</td>
      <td>560.00</td>
    </tr>
    <tr>
      <th>2</th>
      <td>Value of items selected</td>
      <td>374,128.00</td>
      <td>374,076.00</td>
      <td>374,076.00</td>
      <td>374,076.00</td>
    </tr>
    <tr>
      <th>3</th>
      <td>No. of covers generated</td>
      <td>nan</td>
      <td>6.00</td>
      <td>6.00</td>
      <td>11.00</td>
    </tr>
    <tr>
      <th>4</th>
      <td>Computational time</td>
      <td>0.12</td>
      <td>0.42</td>
      <td>0.66</td>
      <td>5.07</td>
    </tr>
  </tbody>
</table>
</div>


8b22ad0ac33ee7fae305741ee5d7afea45c49e74
ipynb
Jupyter Notebook
Knapsack Cover Inequalities/Knapsack Cover Inequalities notebook.ipynb
prameshk/Discrete-Optimization
7cfb1e0977ab071ba5faf2b515e7838ff09fe579
MIT
2021-06-16T07:04:08.000Z
2021-11-08T21:35:01.000Z
Knapsack Cover Inequalities/Knapsack Cover Inequalities notebook.ipynb
prameshk/Discrete-Optimization
7cfb1e0977ab071ba5faf2b515e7838ff09fe579
MIT
Knapsack Cover Inequalities/Knapsack Cover Inequalities notebook.ipynb
prameshk/Discrete-Optimization
7cfb1e0977ab071ba5faf2b515e7838ff09fe579
MIT
2021-11-08T21:51:30.000Z
2021-11-08T21:51:30.000Z
Qwen/Qwen-72B
1. YES
2. YES