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

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))