How do I solve this non-linear optimisation problem based on simulations?

Question

I have an optimisation problem that is essentially a knapsack problem with a non-linear objective.

I have an input dataframe that contains a row for each item, each item has columns defining its mean value and weight. I also have a covariance matrix that contains the variances and covariances between each item.

I take this information and simulate the values for each item a number of times, and store that data.

Then i define my optimisation problem, where i want to select 20 groups of 6 items, so that i maximise the probability of any one of these groups reaching a target value.

Essentially, I am trying to maximise:

$$P(max(x) > T)$$

where x is the values of the groups of lineups, and T is the target score.

My current method:

I am using pyomo in python to solve this problem. First I define a binary decision variable and some constraints, such as making sure no groups are identical and that they are under the weight limit, then i define the objective rule.

I try and do this by summing the simulated scores for each group of items, doing this for each simulation, and keeping track of the number of times the value reaches the target value.

I then try to solve the problem using the solver.

However, when i run it, the model does not initialise values for the binary decision variable, i assume because it cannot find a solution. When i initialise the values myself with groups that make sense, it runs through and calculates the probability, but does not change any of the groups at all.

Here is my current attempt:

# Define input data
n_items = 120
m_groups = 20
target_value = 400
n_simulations = 10000

# Generate simulated scores for all players in the tournament
simulated_scores = multivariate_normal.rvs(mean=df['value'], cov=cov_mat, size=n_simulations)

# Define model
model = ConcreteModel()

# Define decision variables with initial values
model.x = Var(range(n_items), range(m_groups), within=Binary, initialize=0)

def objective_rule(model):
    prob = 0
    for j in range(m_groups):
        # Get indices of selected players in lineup j
        for s in range(n_simulations):
            sim_scores = sum(model.x[i,j] * simulated_scores[s, i] for i in range(n_items))
            if value(sim_scores) > target_value:
                prob += 1
    total_prob = prob / (n_simulations*20)   
    return total_prob

model.objective = Objective(rule=objective_rule, sense=maximize)

# Define constraints
model.players_per_lineup = ConstraintList()
for j in range(m_groups):
    model.players_per_lineup.add(sum(model.x[i,j] for i in range(n_items)) == 6)

model.group_per_item = ConstraintList()
for i in range(n_items):
    model.group_per_item.add(sum(model.x[i,j] for j in range(m_groups)) <= 3)

model.weight_limit = ConstraintList()
for j in range(m_groups):
    model.weight_limit.add(sum(model.x[i,j] * df.iloc[i]['weights'] for i in range(n_items)) <= 300)

# Solve model
solver = SolverFactory('glpk')

I assume this is not working as I am not fully understanding the extent of this optimisation problem, and taking a simplistic approach that cannot be solved using this method.

I'm looking for a simple way to produce a solution for this non-linear problem that produces the optimal groups of items. Ideally something that works with simulations that are not always normally distributed.