How do I include covariance in PulP?

Viewed 46

I have a data frame of ID's, weights and variances (input as dict in code segment below)

I currently have some working code that maximises the sum of 6 variances within the weight limit, but I also want to include covariance in this calculation.

I have a dataframe which I am using as a covariance matrix, which has the variance and covariance relationship between each ID (also included as dict below). I want to include this matrix in the calculation for the variance, by summing each covariance of each pair in the group together, and then summing this with the total variance, to end up maximising the total variance+covariance of the group of 6.

Is this possible to do in Pulp? I do not understand how exactly pulp iterates through the potentially combinations so I struggle to see how I can make this work.

This is my current code:

data = data.set_index('ID')

limit = 50000

w = data.Weight
v = data.Variance

items = list(sorted(data.index))
# Create model
m = LpProblem("Knapsack", LpMaximize)

# Variables
x = LpVariable.dicts('p', items, lowBound=0, upBound=1, cat=LpInteger)

# Objective
m += sum(v[i]*x[i] for i in items)

# Constraint
m += sum(w[i]*x[i] for i in items) <= limit

m += sum(x[i] for i in items) == 6

# Optimize
m.solve()    


{'Weight': {23486388: 6400, 23486470: 4600, 23486060: 5300, 23486018: 6700, 23485910: 7000, 23486424: 5400, 23487054: 3100, 23486516: 3900, 23486026: 6200, 23486412: 5600}, 'Variance': {23486388: 93.44613056000001, 23486470: 61.323284000000015, 23486060: 69.48265599999999, 23486018: 61.02131600000001, 23485910: 68.400996, 23486424: 56.798176000000005, 23487054: 38.09551020408164, 23486516: 48.381836159999985, 23486026: 27.320515999999998, 23486412: 54.67685376}}

{(23485910,): {23485910: 68.400996, 23486018: 9.114639999999998, 23486026: 6.949136326530611, 23486060: 8.937408163265308, 23486412: 11.643567346938775, 23486424: -0.9125706122448977, 23486388: 0.0, 23486470: 0.0, 23487054: 0.0, 23486516: 0.0}, (23486018,): {23485910: 9.114639999999998, 23486018: 61.02131600000001, 23486026: 0.4844571428571431, 23486060: 5.883240816326529, 23486412: -3.7593159183673452, 23486424: 1.5452489795918374, 23486388: 0.0, 23486470: 0.0, 23487054: 0.0, 23486516: 0.0}, (23486026,): {23485910: 6.949136326530611, 23486018: 0.4844571428571431, 23486026: 27.320515999999998, 23486060: 2.548579591836733, 23486412: 9.487351836734694, 23486424: -2.6440285714285716, 23486388: 0.0, 23486470: 0.0, 23487054: 0.0, 23486516: 0.0}, (23486060,): {23485910: 8.937408163265308, 23486018: 5.883240816326529, 23486026: 2.548579591836733, 23486060: 69.48265599999999, 23486412: -9.567245714285711, 23486424: -0.5396938775510204, 23486388: 0.0, 23486470: 0.0, 23487054: 0.0, 23486516: 0.0}, (23486412,): {23485910: 11.643567346938775, 23486018: -3.7593159183673452, 23486026: 9.487351836734694, 23486060: -9.567245714285711, 23486412: 54.67685376, 23486424: -9.41125714285714, 23486388: 0.0, 23486470: 0.0, 23487054: 0.0, 23486516: 0.0}, (23486424,): {23485910: -0.9125706122448977, 23486018: 1.5452489795918374, 23486026: -2.6440285714285716, 23486060: -0.5396938775510204, 23486412: -9.41125714285714, 23486424: 56.798176000000005, 23486388: 0.0, 23486470: 0.0, 23487054: 0.0, 23486516: 0.0}, (23486388,): {23485910: 0.0, 23486018: 0.0, 23486026: 0.0, 23486060: 0.0, 23486412: 0.0, 23486424: 0.0, 23486388: 93.44613056000001, 23486470: 0.0, 23487054: 0.0, 23486516: 0.0}, (23486470,): {23485910: 0.0, 23486018: 0.0, 23486026: 0.0, 23486060: 0.0, 23486412: 0.0, 23486424: 0.0, 23486388: 0.0, 23486470: 61.323284000000015, 23487054: 0.0, 23486516: 0.0}, (23487054,): {23485910: 0.0, 23486018: 0.0, 23486026: 0.0, 23486060: 0.0, 23486412: 0.0, 23486424: 0.0, 23486388: 0.0, 23486470: 0.0, 23487054: 38.09551020408164, 23486516: 0.0}, (23486516,): {23485910: 0.0, 23486018: 0.0, 23486026: 0.0, 23486060: 0.0, 23486412: 0.0, 23486424: 0.0, 23486388: 0.0, 23486470: 0.0, 23487054: 0.0, 23486516: 48.381836159999985}}

I want to end up with an objective function like this:

m += sum((v[i]+covar[i])*x[i] for i in items)

or at least something that achieves the same goal.

1 Answers

Here is an example that I think will help. Realize the comment is correct that you need to consider all of the possible pairs amongst the selected items to include the covariance, so if you are selecting 6 items out of N then you will have choose(6,2) covariance terms, or 15 covariance terms in the result.

You will also need a "helper variable" for each pair, so the problem will get kinda big quick. If you have 100 items to choose from, then you will need (100^2)/2 - 100 of these "pair" variables, which is the number of pairs in the upper triangular (minus the diagonal) of the covariance matrix.

Here is a toy example with 4 items. Note if you change the constraint to allow <= 4 selections, it will only take 3 because of the negative cov in the matrix. If you force 4, it will give an (obviously) inferior solution to 3, and you will see all 6 pairs (choose(4,2)) pairs used.

Code:

# knapsack maximizing variance

import pulp
from itertools import product

v = [1.2, 0.4, 1.3, 0.8]

cv = [[ 1.0, 0.5, 0.4, 0.9 ],
      [ 0.5, 1.0, -0.5, -0.7],
      [ 0.4, -0.5, 1.0, 0.7],
      [ 0.9, -0.7, 0.7, 1.0]]

I = list(range(len(v)))
II = [(i, j) for (i,j) in list(product(I,I)) if j>i]  # upper triangular to prevent dupes
print(II)

m = pulp.LpProblem("var_knapsack", pulp.LpMaximize)

x = pulp.LpVariable.dicts('select', I, cat=pulp.LpBinary)

x_pair = pulp.LpVariable.dicts('pair_select', II, cat=pulp.LpBinary)

# Objective
m += 2*sum(x_pair[i, j] * cv[i][j] for i, j in II) + sum(x[i] * v[i] for i in I)

# Constraints
# link the pair selection to the singles
for i, j in II:
    m += x_pair[i, j] <= (x[i] + x[j]) / 2  
    m += x_pair[i, j] >= x[i] + x[j] - 1

# limit to 3 selections
m += sum(x[i] for i in I) == 3

m.solve()

for i in I:
    if x[i].varValue: print(f'select {i}')

for ii in II:
    if x_pair[ii].varValue: print(f'paired: {ii}')

Yields:

Result - Optimal solution found

Objective value:                7.30000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

select 0
select 2
select 3
paired: (0, 2)
paired: (0, 3)
paired: (2, 3)
Related