Python and GLPK to solve Linear programming problem (infeasible solution)

Question

I'm using Python and the GLPK Solver to find recommendations for an infeasible problem.

I'm optimizing for price, and this is the problem data:

Define problem data

foods = ['Food 0', 'Food 1', 'Food 2', 'Food 3', 'Food 4', 'Food 5', 'Food 6', 'Food 7'] nutrients = ['Nutrient 0', 'Nutrient 1', 'Nutrient 2', 'Nutrient 3', 'Nutrient 4', 'Nutrient 5', 'Nutrient 6', 'Nutrient 7'] food_costs = [2.10, 3.00, 1.50, 4.26, 0.50, 0.84, 2.40, 0.75]

nutrient_matrix = [ [0.10, 0.20, 0.15, 0.05, 0.05, 0.10, 0.15, 0.10], [0.05, 0.10, 0.20, 0.10, 0.15, 0.15, 0.15, 0.15], [0.15, 0.05, 0.10, 0.20, 0.25, 0.10, 0.10, 0.20], [0.05, 0.20, 0.05, 0.15, 0.25, 0.25, 0.20, 0.10], [0.10, 0.15, 0.25, 0.15, 0.05, 0.10, 0.05, 0.20], [0.30, 0.05, 0.15, 0.25, 0.10, 0.15, 0.05, 0.10], [0.20, 0.10, 0.10, 0.05, 0.20, 0.20, 0.15, 0.10], [0.15, 0.25, 0.15, 0.10, 0.15, 0.10, 0.10, 0.15], ] nutrient_bounds = [(0.1, 0.2), (0.05, 0.1), (0.1, 0.2), (0.1, 0.15), (0.1, 0.2), (0.1, 0.2), (0.1, 0.2), (0.1, 0.2)] food_bounds = [(0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3)] It's a nutrition problem where I have foods and nutrients, both having constraints.

Since the problem is infeasible, I have to make changes to the food constraints to make it feasible.

The feasible solution is found when Food 0 has its upper bound set to 0.4 and Food 2 has it's lower bound set to 0.

I'm trying to develop a code that gives recommendations to what Foods need their bounds changed and to what direction (lower / upper) when a problem is infeasible.

The result of the problem should be:

"Make changes to the upper bound of Food 0 ; Make changes to the lower bound of food 2."

I have tried everything from Farkas and Dual variables but I can't find the solution. Can anyone help?

Here is where I stopped on my code:

import cvxpy as cp import numpy as np

Define problem data

foods = ['Food 0', 'Food 1', 'Food 2', 'Food 3', 'Food 4', 'Food 5', 'Food 6', 'Food 7'] nutrients = ['Nutrient 0', 'Nutrient 1', 'Nutrient 2', 'Nutrient 3', 'Nutrient 4', 'Nutrient 5', 'Nutrient 6', 'Nutrient 7'] food_costs = [2.10, 3.00, 1.50, 4.26, 0.50, 0.84, 2.40, 0.75]

nutrient_matrix = [ [0.10, 0.20, 0.15, 0.05, 0.05, 0.10, 0.15, 0.10], [0.05, 0.10, 0.20, 0.10, 0.15, 0.15, 0.15, 0.15], [0.15, 0.05, 0.10, 0.20, 0.25, 0.10, 0.10, 0.20], [0.05, 0.20, 0.05, 0.15, 0.25, 0.25, 0.20, 0.10], [0.10, 0.15, 0.25, 0.15, 0.05, 0.10, 0.05, 0.20], [0.30, 0.05, 0.15, 0.25, 0.10, 0.15, 0.05, 0.10], [0.20, 0.10, 0.10, 0.05, 0.20, 0.20, 0.15, 0.10], [0.15, 0.25, 0.15, 0.10, 0.15, 0.10, 0.10, 0.15], ] nutrient_bounds = [(0.1, 0.2), (0.05, 0.1), (0.1, 0.2), (0.1, 0.15), (0.1, 0.2), (0.1, 0.2), (0.1, 0.2), (0.1, 0.2)] food_bounds = [(0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3), (0.1, 0.3)]

Define the function to solve the problem

def solve_problem(food_bounds, nutrient_bounds, nutrient_matrix, food_costs): # Create the variables food_vars = cp.Variable(len(foods), nonneg=True)

# Objective function: Minimize cost
objective = cp.Minimize(cp.sum(food_costs @ food_vars))

# Nutrient constraints
nutrient_constraints = [
    cp.sum(np.array(nutrient_matrix)[i] @
           food_vars) >= nutrient_bounds[i][0]
    for i in range(len(nutrients))
] + [
    cp.sum(np.array(nutrient_matrix)[i] @
           food_vars) <= nutrient_bounds[i][1]
    for i in range(len(nutrients))
]

# Food constraints
food_constraints = [
    food_vars[i] >= food_bounds[i][0] for i in range(len(foods))
] + [
    food_vars[i] <= food_bounds[i][1] for i in range(len(foods))
]

# Sum of food percentages constraint
total_percentage_constraint = [cp.sum(food_vars) == 1]

# Combine all constraints
constraints = nutrient_constraints + \
    food_constraints + total_percentage_constraint

# Create the problem
prob = cp.Problem(objective, constraints)

# Solve the problem
prob.solve(solver=cp.GLPK)

return prob, constraints, food_vars

def compute_dual_values(original_prob, constraints): dual_vars = [cp.Variable(nonneg=True) for _ in constraints] slack_vars = [cp.Variable() for _ in constraints]

dual_constraints = []
for dual_var, constraint, slack_var in zip(dual_vars, constraints, slack_vars):
    A = constraint.args[0].T  # Transpose of the constraint matrix
    x = constraint.args[1]
    dual_constraints.append(slack_var == A @ x - dual_var)

dual_prob = cp.Problem(cp.Minimize(cp.sum(slack_vars)), dual_constraints)
dual_prob.solve(solver=cp.GLPK)

return [dual_var.value for dual_var in dual_vars]

def analyze_food_constraints(constraints, dual_vars): food_constraint_indices = list( range(len(constraints) - len(foods) * 2, len(constraints)))

largest_duals = sorted(food_constraint_indices,
                       key=lambda i: abs(dual_vars[i]), reverse=True)

recommendations = []
for i in largest_duals:
    if abs(dual_vars[i]) > 1e-6:
        food_index = i - len(constraints) + len(foods)
        if food_index < len(foods):
            bound_type = "lower"
        else:
            food_index -= len(foods)
            bound_type = "upper"

        recommendations.append(
            f"Recommend changing the {bound_type} bound of Food {food_index}.")

return recommendations

Call the function to solve the original problem

original_prob, constraints, food_vars = solve_problem( food_bounds, nutrient_bounds, nutrient_matrix, food_costs)

Check if the original problem is infeasible

if original_prob.status == cp.INFEASIBLE: print("The original problem is infeasible. Analyzing constraints...")

# Obtain the dual variables (shadow prices)
dual_vars = compute_dual_values(original_prob, constraints)

# Analyze food constraints
recommendations = analyze_food_constraints(constraints, dual_vars)

# Print recommendations
for recommendation in recommendations:
    print(recommendation)

else: print("Status:", original_prob.status)

# Print the optimal food quantities and cost if the problem is feasible
print("Optimal food quantities (percentage of meal):")
for i, food in enumerate(foods):
    print(f"{food}: {food_vars.value[i] * 100:.2f}%")

print("Optimal meal cost per kilogram: ${:.2f}".format(
    original_prob.value))

Answer 1

The first step is deciding what constraints you are willing to relax. It appears that you have chosen to relax the bounds on food quantities (as opposed to nutrient requirements).

The next step is to decide how to "price" the relaxations. Is relaxing the upper bound on food #1 by 1 kg. better than, the same as or worse than relaxing the lower bound on food #2 by 1 kg.? This is a question that probably needs to be answered by whoever instituted those bounds in the first place.

A possible third step is to decide how much you are willing to relax each bound. That might be an individual decision or a "group" decision. For instance, if there are three foods that are considered "pasta", you might say that you are willing to increase total pasta by 10 kg. and leave it to the solver to decide how to split the increase among the three foods.

Once you know all this, you can add variables to the model representing the increase or decrease in each relevant bound, change the objective function to minimizing the total cost/penalty for those changes, and solve the modified problem. The solution will tell you how to adjust the bounds to make the problem feasible, after which you can solve the original model with the modified bounds.

Besides the somewhat subjective aspects of the questions above, there is one other caveat. If you make the minimum changes necessary to achieve feasibility, you might wind up with a model with a very small feasible region (perhaps a single point), which might produce a very unsatisfying solution to the original problem. An alternative approach is to insert the new variables (sometimes described as "artificial" variables) and optimize a weighted combination of the original objective and a cost/penalty for the bounds changes. The tricky part of this approach is assigning a cost for the changes (e.g., increasing the bound on food 1 by 1 kg. "costs" the same as buying 3 kg. of food 2).