# How do I solve a probability based knapsack problem?

I am looking to maximize my probability of reaching a given target value, by creating multiple groups of 6 items with different means and variances that stay within a weight limit. This is an example of the top of the data frame - the actual data frame is 300+ rows long.

index mu std Weight
item_a 73.450198 10.455766 111.
item_b 72.652175 9.477475 110.
item_c 73.033862 10.293721 108.
item_d 73.886648 10.426305 105.
item_e 68.409628 10.588617 103.

I want to maximise the probability that a group of 6 items exceeds a given value, whilst still under the weight limit (500).

It was pointed out to me that instead of calculating the exact probability I can minimize the z score for each group using the formula:

z_score = (target_score-mu_group) / sd_group


When I try to minimize this formula in CPLEX it minimizes the value for the sum of each z score instead of the z score of the group:

from docplex.mp.model import Model
target_score = 600
weight_limit = 500
w = df.weight
v = df.mu
y = df.std
N = len(w)

knapsack_model = Model('knapsack')

x = knapsack_model.binary_var_list(N, name= 'x')

knapsack_model.add_constraint(sum(w[i]*x[i] for i in range(N)) <= weight_limit)
knapsack_model.add_constraint(sum(x[i] for i in range(N)) == 6)

obj_fn = sum(((target_score-v[i])/((y[i]*y[i])) for i in range(N))
knapsack_model.set_objective('min',obj_fn)

knapsack_model.solve()


I am aware that this is because the binary variable used for each team is not within the z score formula, and is therefore only minimizing the sum of all z scores in the group instead of the z score of the group as a whole.

I imagine I will need a non-linear solver as the z score calculation is not linear, but I am not sure how to formulate this as I am reaching the edge of my understanding in this area.

The other non-trivial part of the problem I would like to solve is to maximize the probability of multiple groups reaching the target score - so create x amount of groups that maximize the probability of one of them reaching the target score.

I would also like to have control of the percentages of each item across these groups, for example being able to specify that item_a is present in no more than 10% of groups.

In my initial problem the items will be independent, but in the future I will be looking to incorporate a covariance matrix in to the solution.

I have been working primarily in python.

Notes: assume this problem is not solve-able using brute force. All items are normally distributed.

{'ID': {0: 23516794, 1: 23516795, 2: 23516796, 3: 23516797, 4: 23516798, 5: 23516799, 6: 23516800, 7: 23516801, 8: 23516802, 9: 23516803, 10: 23516804, 11: 23516805, 12: 23516806, 13: 23516807, 14: 23516808, 15: 23516809, 16: 23516810, 17: 23516811, 18: 23516812, 19: 23516813}, 'mu': {0: 77.22374694548756, 1: 84.05634969802539, 2: 80.96531404866622, 3: 87.57033049910319, 4: 74.78393207733146, 5: 79.82400732228508, 6: 78.16565147466483, 7: 78.14728328442523, 8: 75.98898560579204, 9: 78.8541978528528, 10: 75.07931255975794, 11: 85.11385652194893, 12: 76.72146255952251, 13: 80.29834170101859, 14: 77.86295417762425, 15: 74.70045311388368, 16: 81.24594041364621, 17: 76.21212603618336, 18: 73.42569984426565, 19: 76.1438949931049}, 'std': {0: 11.68598172116063, 1: 13.327661965136674, 2: 14.834799563943232, 3: 9.453573748770634, 4: 12.528214306044084, 5: 13.435066788790433, 6: 11.933781990021213, 7: 12.933356356518166, 8: 12.258175910542153, 9: 15.115553286074276, 10: 13.029399969876495, 11: 15.880786524364469, 12: 13.246754947691105, 13: 12.340407352537394, 14: 13.66745084735041, 15: 12.626895384909284, 16: 14.89605081596392, 17: 13.241317159371539, 18: 12.294456354577791, 19: 12.249112531502997}, 'Weight': {0: 96.0, 1: 109.0, 2: 105.0, 3: 86.0, 4: 87.0, 5: 107.0, 6: 81.0, 7: 82.0, 8: 106.0, 9: 84.0, 10: 83.0, 11: 85.0, 12: 98.0, 13: 94.0, 14: 92.0, 15: 103.0, 16: 80.0, 17: 88.0, 18: 89.0, 19: 90.0}}

{'ID': {0: 23516794, 1: 23516795, 2: 23516796, 3: 23516797, 4: 23516798, 5: 23516799, 6: 23516800, 7: 23516801, 8: 23516802, 9: 23516803, 10: 23516804, 11: 23516805, 12: 23516806, 13: 23516807, 14: 23516808, 15: 23516809, 16: 23516810, 17: 23516811, 18: 23516812, 19: 23516813}, 0: {0: 136.56216878730032, 1: -4.494339538364646, 2: -14.73409159687065, 3: -1.8369732959606888, 4: 1.5974817060365636, 5: -9.60875988350014, 6: -2.183995750123481, 7: 2.514233365391022, 8: -1.916442017276486, 9: -11.679165869066852, 10: -6.306976154792111, 11: -11.332834201301306, 12: -2.3812635674739857, 13: 1.4835959725328676, 14: -5.657222725636583, 15: -9.79722153520659, 16: -12.603609517512604, 17: -2.6398894441141043, 18: -2.291667714493622, 19: -7.3855786001326775}, 1: {0: -4.494339538364646, 1: 177.62657345695047, 2: -12.628111759388075, 3: -4.443689614302776, 4: -10.738219362736135, 5: -2.867861039491733, 6: -8.024258930122437, 7: -5.8343736141719615, 8: -5.0856685690709735, 9: -5.810127855540393, 10: -6.486789218975431, 11: -4.536185413209171, 12: -14.60542313812631, 13: -0.600165154002184, 14: -6.212283004005756, 15: -3.4578205514451863, 16: -6.836266779603481, 17: -9.031669667810993, 18: -9.651812989565943, 19: -0.3322923570894588}, 2: {0: -14.73409159687065, 1: -12.628111759388075, 2: 220.07127810237012, 3: -5.472041080378147, 4: -0.5530251090807978, 5: -4.843046719622863, 6: 2.652486193022895, 7: -5.864881828798867, 8: -13.817643791422551, 9: -10.526409041166408, 10: -19.840091014093485, 11: -10.619099101842991, 12: -7.225864283168083, 13: 0.5735926596177876, 14: 4.0005089456860965, 15: -8.240192896686342, 16: -2.7591820700695866, 17: -8.557254208071246, 18: 0.16726741873534237, 19: -1.0993500347207696}, 3: {0: -1.8369732959606888, 1: -4.443689614302776, 2: -5.472041080378147, 3: 89.3700566234453, 4: -8.848241985121726, 5: -6.625741125767621, 6: -0.21577281576776328, 7: -0.8512025218040024, 8: 0.6751643162844593, 9: 1.0524602516495232, 10: -1.7817610644992792, 11: -5.389103259494309, 12: -3.8824825750140897, 13: -7.0157826193021755, 14: -6.87163496978418, 15: -11.560853860022556, 16: 5.356758878641404, 17: -0.9098172019944819, 18: -1.9189977283684865, 19: -0.9362378417573064}, 4: {0: 1.5974817060365636, 1: -10.738219362736135, 2: -0.5530251090807978, 3: -8.848241985121726, 4: 156.95615369816767, 5: -8.768969357976829, 6: 0.6217679553017051, 7: -0.20504497477901446, 8: -1.0831551767658514, 9: -1.3066968240860621, 10: 5.490867046711869, 11: -2.8210361046394854, 12: 6.475286818155162, 13: -8.802027521591068, 14: -7.569425287856836, 15: -9.454334011368, 16: -9.679516491054699, 17: -9.165398431012328, 18: -4.73040679319504, 19: -7.686554609840943}, 5: {0: -9.60875988350014, 1: -2.867861039491733, 2: -4.843046719622863, 3: -6.625741125767621, 4: -8.768969357976829, 5: 180.50101961925955, 6: -7.054933927313439, 7: -11.358041935082174, 8: -8.629615651443382, 9: -7.1476470240365115, 10: -13.270380194728848, 11: 2.741649714411773, 12: -3.7257538015023197, 13: -4.8121739826093135, 14: 3.894065439133523, 15: 5.628402517781505, 16: -15.985098818113466, 17: -3.0115403619146024, 18: 1.12320328990123, 19: -1.025631802973544}, 6: {0: -2.183995750123481, 1: -8.024258930122437, 2: 2.652486193022895, 3: -0.21577281576776328, 4: 0.6217679553017051, 5: -7.054933927313439, 6: 142.4151525853546, 7: -0.1342070188433831, 8: -7.908278633887033, 9: -5.3017760571802794, 10: 4.9529559923762525, 11: -3.749047764361528, 12: -4.772361787449819, 13: -8.250291267512742, 14: -0.575220665125228, 15: -10.83153808101616, 16: -1.2168714342029199, 17: -3.8947200971928373, 18: -6.738190917666668, 19: -0.9565234278762215}, 7: {0: 2.514233365391022, 1: -5.8343736141719615, 2: -5.864881828798867, 3: -0.8512025218040024, 4: -0.20504497477901446, 5: -11.358041935082174, 6: -0.1342070188433831, 7: 167.27170664468886, 8: -1.7306718846505822, 9: -19.695139284680625, 10: -7.246281550693982, 11: -6.446113218776497, 12: -8.690130833341392, 13: -4.403524906950436, 14: -14.011147795374235, 15: -5.137164556199244, 16: -11.341072122630719, 17: 0.10147369388262845, 18: -1.7697234078007908, 19: 5.975718005677391}, 8: {0: -1.916442017276486, 1: -5.0856685690709735, 2: -13.817643791422551, 3: 0.6751643162844593, 4: -1.0831551767658514, 5: -8.629615651443382, 6: -7.908278633887033, 7: -1.7306718846505822, 8: 150.262876653796, 9: -15.072231128833636, 10: -0.1496075400717749, 11: -15.449852852260035, 12: -1.7637338713840054, 13: -2.0771821311264778, 14: -9.045513112177188, 15: -3.2033101826605255, 16: -7.9132725471029035, 17: -1.5278980319764106, 18: -6.262416358407891, 19: -12.651154101625123}, 9: {0: -11.679165869066852, 1: -5.810127855540393, 2: -10.526409041166408, 3: 1.0524602516495232, 4: -1.3066968240860621, 5: -7.1476470240365115, 6: -5.3017760571802794, 7: -19.695139284680625, 8: -15.072231128833636, 9: 228.47995114415113, 10: 1.5522799948287105, 11: -7.352331898572419, 12: -12.005452429573031, 13: -11.435036416340306, 14: 0.7097071851524513, 15: -6.134332961139446, 16: -14.090526602277816, 17: -2.2206619756610806, 18: -10.242162355051265, 19: -2.6407749930270024}, 10: {0: -6.306976154792111, 1: -6.486789218975431, 2: -19.840091014093485, 3: -1.7817610644992792, 4: 5.490867046711869, 5: -13.270380194728848, 6: 4.9529559923762525, 7: -7.246281550693982, 8: -0.1496075400717749, 9: 1.5522799948287105, 10: 169.76526357501763, 11: -6.376673487808323, 12: -10.590716933267549, 13: -3.4809085348779503, 14: -6.065589307674493, 15: -4.950793172098783, 16: 1.7465983331624262, 17: -12.462347651961771, 18: -6.79247415233182, 19: 6.0659190790992215}, 11: {0: -11.332834201301306, 1: -4.536185413209171, 2: -10.619099101842991, 3: -5.389103259494309, 4: -2.8210361046394854, 5: 2.741649714411773, 6: -3.749047764361528, 7: -6.446113218776497, 8: -15.449852852260035, 9: -7.352331898572419, 10: -6.376673487808323, 11: 252.19938063243615, 12: 5.204041430308693, 13: -23.284827959906316, 14: -6.531768966914284, 15: -3.494555758287847, 16: -14.535184801941293, 17: -4.55733965420517, 18: 4.111678333396179, 19: -2.316270358926943}, 12: {0: -2.3812635674739857, 1: -14.60542313812631, 2: -7.225864283168083, 3: -3.8824825750140897, 4: 6.475286818155162, 5: -3.7257538015023197, 6: -4.772361787449819, 7: -8.690130833341392, 8: -1.7637338713840054, 9: -12.005452429573031, 10: -10.590716933267549, 11: 5.204041430308693, 12: 175.47651664417913, 13: -12.99506383181884, 14: -10.512473183364055, 15: -4.949916172219996, 16: -6.5276526791083676, 17: -7.299556696505636, 18: 3.226564599579514, 19: 2.084662921550234}, 13: {0: 1.4835959725328676, 1: -0.600165154002184, 2: 0.5735926596177876, 3: -7.0157826193021755, 4: -8.802027521591068, 5: -4.8121739826093135, 6: -8.250291267512742, 7: -4.403524906950436, 8: -2.0771821311264778, 9: -11.435036416340306, 10: -3.4809085348779503, 11: -23.284827959906316, 12: -12.99506383181884, 13: 152.2856536265588, 14: -5.372661979592924, 15: -0.8370127370548941, 16: -9.78689911529412, 17: -3.1964142604668715, 18: 8.707018182079095, 19: -8.424931545799005}, 14: {0: -5.657222725636583, 1: -6.212283004005756, 2: 4.0005089456860965, 3: -6.87163496978418, 4: -7.569425287856836, 5: 3.894065439133523, 6: -0.575220665125228, 7: -14.011147795374235, 8: -9.045513112177188, 9: 0.7097071851524513, 10: -6.065589307674493, 11: -6.531768966914284, 12: -10.512473183364055, 13: -5.372661979592924, 14: 186.79921266473949, 15: -1.5788775100724504, 16: -16.952483639367813, 17: -14.518689150558647, 18: -11.119584217756573, 19: -10.17339135612405}, 15: {0: -9.79722153520659, 1: -3.4578205514451863, 2: -8.240192896686342, 3: -11.560853860022556, 4: -9.454334011368, 5: 5.628402517781505, 6: -10.83153808101616, 7: -5.137164556199244, 8: -3.2033101826605255, 9: -6.134332961139446, 10: -4.950793172098783, 11: -3.494555758287847, 12: -4.949916172219996, 13: -0.8370127370548941, 14: -1.5788775100724504, 15: 159.43848706144342, 16: -5.604885004509065, 17: -3.7491992092835695, 18: -2.4676427841801045, 19: -6.103848391110505}, 16: {0: -12.603609517512604, 1: -6.836266779603481, 2: -2.7591820700695866, 3: 5.356758878641404, 4: -9.679516491054699, 5: -15.985098818113466, 6: -1.2168714342029199, 7: -11.341072122630719, 8: -7.9132725471029035, 9: -14.090526602277816, 10: 1.7465983331624262, 11: -14.535184801941293, 12: -6.5276526791083676, 13: -9.78689911529412, 14: -16.952483639367813, 15: -5.604885004509065, 16: 221.89232991177963, 17: -5.850040172649914, 18: -15.64551677338866, 19: -11.237751321651324}, 17: {0: -2.6398894441141043, 1: -9.031669667810993, 2: -8.557254208071246, 3: -0.9098172019944819, 4: -9.165398431012328, 5: -3.0115403619146024, 6: -3.8947200971928373, 7: 0.10147369388262845, 8: -1.5278980319764106, 9: -2.2206619756610806, 10: -12.462347651961771, 11: -4.55733965420517, 12: -7.299556696505636, 13: -3.1964142604668715, 14: -14.518689150558647, 15: -3.7491992092835695, 16: -5.850040172649914, 17: 175.3324801150675, 18: -3.2911916920963638, 19: 1.033994997644956}, 18: {0: -2.291667714493622, 1: -9.651812989565943, 2: 0.16726741873534237, 3: -1.9189977283684865, 4: -4.73040679319504, 5: 1.12320328990123, 6: -6.738190917666668, 7: -1.7697234078007908, 8: -6.262416358407891, 9: -10.242162355051265, 10: -6.79247415233182, 11: 4.111678333396179, 12: 3.226564599579514, 13: 8.707018182079095, 14: -11.119584217756573, 15: -2.4676427841801045, 16: -15.64551677338866, 17: -3.2911916920963638, 18: 151.1536570546182, 19: 0.3524880702152624}, 19: {0: -7.3855786001326775, 1: -0.3322923570894588, 2: -1.0993500347207696, 3: -0.9362378417573064, 4: -7.686554609840943, 5: -1.025631802973544, 6: -0.9565234278762215, 7: 5.975718005677391, 8: -12.651154101625123, 9: -2.6407749930270024, 10: 6.0659190790992215, 11: -2.316270358926943, 12: 2.084662921550234, 13: -8.424931545799005, 14: -10.17339135612405, 15: -6.103848391110505, 16: -11.237751321651324, 17: 1.033994997644956, 18: 0.3524880702152624, 19: 150.04075780942372}}

• Are you assuming that item values are (a) normally distributed and (b) independent of each other?
– prubin
Jul 28, 2022 at 15:29
• yes that is correct. I am looking to include a covariance matrix at some point in the future - I will add that in to the problem to make it clear.
– will
Jul 28, 2022 at 15:44
• Are you able to share the full data? Jul 28, 2022 at 15:49
• @RobPratt what more information is needed from the data? Do you just want the data frame?
– will
Jul 28, 2022 at 15:54
• Yes, just the data frame. Jul 28, 2022 at 15:56

I'll tackle here the question of selecting a single team (of specified size with some upper limit on the weight) with maximum probability of success (which I'll operationalize by maximizing z-score (mu_team-target_score)/sd_team). Though at first blush this sounds like a very difficult optimization problem (a non-linear integer optimization problem), it turns out it can be effectively handled with branch and bound. I'll first explain the high-level idea and then share some working code.

Imagine you're targeting a score of 550 (subject to some team size and weight requirements, which I'll assume are added as constraints in all optimization problems). Now imagine you solved two initial optimization problems:

• You optimize to maximize the mean of your team (ignoring its variance), and you get a team with a mean of 597.4 and a variance of 490.5.
• You optimize to minimize the variance of your team among any team whose mean is at least the target score of 550, and you learn that the minimum variance that can be achieved for an above-target team is 397.7.

At this point, we know two things:

1. You have already found a team with z-score (597.4-550)/sqrt(490.5) = 2.14, so the optimal solution must have a z-score at least that good.
2. Since every team with an above-target score has a variance of at least 397.7, no optimal solution can have a z-score of better than (597.4-550)/sqrt(397.7) = 2.38.

Now we have an upper and a lower bound -- the optimal solution's z-score falls between 2.14 and 2.38. The key question is how to proceed toward the optimal solution. The idea is an interesting one -- we pick a variance threshold of 444.1 (halfway between 397.7 and 490.5) and we optimize the team with the best mean subject to having a variance below 444.1. The optimization problem returns a team with mean 595.1 and a variance of 440.4. Now we can partition the space of solutions into two ranges:

• Region 1: For the set of teams with variance between 444.1 and 490.5, the best achievable mean is 597.4; the z-score of the best solution in the region must be at least (597.4-550)/sqrt(490.5) = 2.14 but can be no more than (597.4-550)/sqrt(444.1) = 2.25.
• Region 2: For the set of teams with variance between 397.7 and 444.1, the best achievable mean is 595.1; the z-score of the best solution in the region must be at least (595.1-550)/sqrt(440.4) = 2.15 but can be no more than (595.1-550)/sqrt(397.7) = 2.26.

Great news -- after this second round we now know that the overall optimal solution must have a z-score of at least max(2.14, 2.15) = 2.15, but can't have a z-score any bigger than max(2.25, 2.26) = 2.26. Both our upper and lower bounds have improved!

Branch and bound now proceeds by iteratively identifying the region with the largest upper bound (region 2 at the moment) and splits it as we did above, always improving the upper bound at least a little bit and sometimes improving the lower bound as well. You keep the process going until the upper and lower bounds are sufficiently close.

This is quite straightforward to do with any optimization solver. I'll use the python pulp package here:

import math
import pulp

# Optimize a model
def opt(optVar, numSelect, weight_limit, varMin=None, varMax=None, muMin=None):
allId = [x for x in df["ID"]]
x = pulp.LpVariable.dicts("x", allId, 0, 1, pulp.LpBinary)

# Set objective to either maximize the mean or minimize the variance
if optVar == "mean":
prob = pulp.LpProblem(optVar, pulp.LpMaximize)
prob += pulp.lpSum([df["mu"][idx] * x[idx] for idx in allId])
else:
if optVar == "minVar":
prob = pulp.LpProblem(optVar, pulp.LpMinimize)
else:
prob = pulp.LpProblem(optVar, pulp.LpMaximize)
prob += pulp.lpSum([df["std"][idx]*df["std"][idx]*x[idx] for idx in allId])

# Weight, variance, and mean limits as well as cardinality:
prob += pulp.lpSum([df["Weight"][idx]*x[idx] for idx in allId]) <= weight_limit
prob += pulp.lpSum(x) == numSelect
if varMin is not None:
prob += pulp.lpSum([df["std"][idx]*df["std"][idx]*x[idx] for idx in allId]) >= varMin
if varMax is not None:
prob += pulp.lpSum([df["std"][idx]*df["std"][idx]*x[idx] for idx in allId]) <= varMax
if muMin is not None:
prob += pulp.lpSum([df["mu"][idx]*x[idx] for idx in allId]) >= muMin

# Solve and return solution
prob.solve(pulp.PULP_CBC_CMD(msg=0))
selected = [idx for idx in allId if x[idx].value() > 0.999]
mu = sum([df["mu"][idx]*x[idx].value() for idx in allId])
var = sum([df["std"][idx]*df["std"][idx]*x[idx].value() for idx in allId])
return mu, var, selected

# Optimize a team with the specified parameters
def optTeam(numSelect, target_score, weight_limit, tolerance, verbose=False):
M, Mvar, best = opt("mean", numSelect, weight_limit, None, None, None)
if M == target_score:
return 0, 0, best
elif M > target_score:
positiveZ = True
_, Vbest, _ = opt("minVar", numSelect, weight_limit, None, None, target_score)
else:
positiveZ = False
_, Vbest, _ = opt("maxVar", numSelect, weight_limit, None, None, None)

nodes = [[Vbest, Mvar, M, best]]

while True:
allLB = [(x[2]-target_score)/math.sqrt(x[1]) for x in nodes]
currLB, currLBPos = max([(val, idx) for idx, val in enumerate(allLB)])
allUB = [(x[2]-target_score)/math.sqrt(x[0]) for x in nodes]
currUB, currUBPos = max([(val, idx) for idx, val in enumerate(allUB)])
if verbose:
print(len(nodes), currLB, currUB)

if (positiveZ and currUB / currLB < tolerance) or (not positiveZ and currLB / currUB < tolerance):
break

# Split the node with the largest upper bound
mid = (nodes[currUBPos][0] + nodes[currUBPos][1]) / 2.0
if positiveZ:
newMu, newVar, newBest = opt("mean", numSelect, weight_limit, None, mid, None)
else:
newMu, newVar, newBest = opt("mean", numSelect, weight_limit, mid, None, None)
nodes.append([nodes[currUBPos][0], newVar, newMu, newBest])
nodes[currUBPos][0] = mid

return currLB, currUB, nodes[currLBPos][3]

LB, UB, sol = optTeam(6, 600, 500, 1.0001, True)
print(sol)


The approach outputs the lower and upper bound from each step of the branch and bound procedure. Here I optimize with target score 600 and weight limit 500, which is done in under a second of computation time with the default CBC solver:

1 -7.1048510789624855 -6.060011148175446
2 -6.611541985047641 -6.172741906265993
3 -6.355377019877131 -6.250015803032588
4 -6.355377019877131 -6.3020359480846775
5 -6.355377019877131 -6.328537892630254
6 -6.355377019877131 -6.341914862728255
7 -6.355377019877131 -6.3486352364871745
8 -6.355377019877131 -6.352003444878156
9 -6.355377019877131 -6.353689560661244
10 -6.355377019877131 -6.354533122228607
11 -6.355377019877131 -6.354955029028778
[10, 12, 15, 23, 38, 55]


Clearly the optimal team will rarely break score 600, with probability of just 1.04e-10. Lower target scores are much more achievable -- the optimal team with target score 500 has probability 0.55%, the optimal team with target score 450 has probability 30.4%, and the optimal team with target score 400 has probability 95.7%.

• This is great, thanks! I really like the approach. Would a similar idea be possible while including a covariance matrix?
– will
Aug 2, 2022 at 20:26
• @will If you're referring selecting a single team where there are covariances between the scores of those team members, that should not be too tough. You would update the opt function to compute the variance of a team as the sum of the variances of each selected player plus the sum of the covariances of each pair of selected players. You'll need decision variable products x_i * x_j for each pair (i, j) with non-zero covariance to do this. You can search this site for how to model the product of binary variables in a formulation. Aug 3, 2022 at 17:04
• @will I'd be glad to update the code myself to capture covariances, but I'd ask you to first update your data example to include some example covariances. Aug 4, 2022 at 16:41
• I have updated the data frames to include a covariance matrix. I had to reduce the size of the data as the whole covariance matrix for 120 items wouldnt fit in this post. Thanks for your help! @josliber
– will
Aug 10, 2022 at 16:41
• Hey @josliber, do you need a different format for the covariance dataframe to make this possible? Happy to edit the input data as needed.
– will
Aug 31, 2022 at 8:40

I coded the initial version of the problem (finding one group of six items) in R using a permutation-based genetic algorithm. Since this is a metaheuristic, (a) you are not guaranteed an optimal solution and (b) you can futz with all sorts of parameters (including but not limited to run length/time limits), as well as trying different random number seeds (or just rerunning), and get different results. As to how good the solution is, that is an open question until you come up with a benchmark.

As a caveat, I had to reduce your target score. As best I can tell, no combination of six items has even a remote chance of producing a score of 600 or more (meaning the highest probability the GA could get was around $$10^{-12}$$). The problem with that is that improvements are so tiny the GA code thinks it has stalled and quits early. With a target score of 450, the GA was able to find a combination with better than 15% chance of success.

The general idea is that the GA plays with candidate solutions that are permutations of the indices of your items (1 to 108 for your sample data). Given a permutation of items, a somewhat cumbersome function goes down the list and selects the six items with highest rank (earliest in the list) that collectively stay within the budget, and then calculates the likelihood of meeting or exceed the target assuming their payoffs are i.i.d. normal. There is a bunch of sorting involved in the move from permutation to item list, but the GA ran in seconds nonetheless.

Accommodating a covariance matrix should not be too hard. All it changes is the probability calculation given a choice of items. Accommodating disjoint sets of items (e.g., select four groups of 10 with no overlap) would be a somewhat straightforward extension. Selecting multiple groups with overlap might be trickier. Calculating the probability of a win would be fairly straightforward for disjoint groups, but with overlapping groups I think it would get a bit convoluted.