Diagnose infeasibility#
Description: This notebook demonstrates how to deal with infeasible models.
Tags: tutorials, infeasibility, mp
Notebook author: Marcos Dominguez Velad <marcos@ampl.com>
Model author: N/A
References:
AMPL a Modeling Language for Mathematical Programming – Robert Fourer et al.
# Install dependencies
%pip install -q amplpy pandas
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["gurobi"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
For a linear program that has no feasible solution, you can ask solvers like Gurobi or CPLEX to find an irreducible infeasible subset (or IIS): a collection of constraints and variable bounds that is infeasible but that becomes feasible when any one constraint or bound is removed. If a small IIS exists and can be found, it can provide valuable clues as to the source of the infeasibility. You turn on the IIS finder by changing the iisfind directive from its default value of 0 to either 1 (for a relatively fast version) or 2 (for a slower version that tends to find a smaller IIS).
The following example shows how IIS finding might be applied to the infeasible diet problem.
%%writefile diet.mod
set NUTR;
set FOOD;
param cost {FOOD} > 0;
param f_min {FOOD} >= 0;
param f_max {j in FOOD} >= f_min[j];
param n_min {NUTR} >= 0;
param n_max {i in NUTR} >= n_min[i];
param amt {NUTR,FOOD} >= 0;
var Buy {j in FOOD} >= f_min[j], <= f_max[j];
minimize Total_Cost: sum {j in FOOD} cost[j] * Buy[j];
subject to Diet {i in NUTR}:
n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];
Writing diet.mod
import pandas as pd
import numpy as np
# Simply load data for the problem
def load_data(model):
# First generate the data in Python
foods = ["BEEF", "CHK", "FISH", "HAM", "MCH", "MTL", "SPG", "TUR"]
nutrs = ["A", "B1", "B2", "C", "NA", "CAL"]
food_df = pd.DataFrame(
[
("BEEF", 3.19, 2, 10),
("CHK", 2.59, 2, 10),
("FISH", 2.29, 2, 10),
("HAM", 2.89, 2, 10),
("MCH", 1.89, 2, 10),
("MTL", 1.99, 2, 10),
("SPG", 1.99, 2, 10),
("TUR", 2.49, 2, 10),
],
columns=["FOOD", "cost", "f_min", "f_max"],
).set_index("FOOD")
nutr_df = pd.DataFrame(
[
("A", 700, 20000),
("C", 700, 20000),
("B1", 700, 20000),
("B2", 700, 20000),
("NA", 0, 40000),
("CAL", 16000, 24000),
],
columns=["NUTR", "n_min", "n_max"],
).set_index("NUTR")
amt_df = pd.DataFrame(
np.array(
[
[60, 8, 8, 40, 15, 70, 25, 60],
[20, 0, 10, 40, 35, 30, 50, 20],
[10, 20, 15, 35, 15, 15, 25, 15],
[15, 20, 10, 10, 15, 15, 15, 10],
[938, 2180, 945, 278, 1182, 896, 1329, 1397],
[295, 770, 440, 430, 315, 400, 370, 450],
]
),
columns=food_df.index.to_list(),
index=nutr_df.index.to_list(),
)
# Load the data into Ampl:
ampl.set["FOOD"] = foods
ampl.set["NUTR"] = nutrs
# Load data related to foods:
ampl.set_data(food_df, "FOOD")
# Load data related to nutrients:
ampl.set_data(nutr_df, "NUTR")
ampl.param["amt"] = amt_df
After solve detects that there is no feasible solution, it is repeated with the directive iisfind 1:
ampl.read("diet.mod")
load_data(ampl)
ampl.solve(solver="gurobi")
Gurobi 12.0.3: Gurobi 12.0.3: infeasible problem
2 simplex iterations
suffix dunbdd OUT;
ampl.solve(solver="gurobi", mp_options="iisfind=1")
Gurobi 12.0.3: alg:iisfind = 1
Gurobi 12.0.3: infeasible problem
0 simplex iterations
suffix iis symbolic OUT;
option iis_table '\
0 non not in the iis\
1 low lower bound in the iis\
2 fix both bounds in the iis\
3 upp upper bound in the iis\
4 mem member\
5 pmem possible member\
6 plow possibly lower bound\
7 pupp possibly upper bound\
8 bug\
';
Again, AMPL shows any suffix statement that has been executed automatically. Our interest is in the new suffix named .iis. An associated option iis_table, also set up by the solver driver and displayed automatically by solve, shows the strings that may be associated with .iis and gives brief descriptions of what they mean. You can use display to look at the .iis values that have been returned:
You can display the .iis values that have been returned:
vars_iis_df = ampl.get_data("_varname, _var.iis").to_pandas()
cons_iis_df = ampl.get_data("_conname, _con.iis").to_pandas()
print(vars_iis_df)
print(cons_iis_df)
_varname _var.iis
1 Buy['BEEF'] low
2 Buy['CHK'] low
3 Buy['FISH'] low
4 Buy['HAM'] upp
5 Buy['MCH'] low
6 Buy['MTL'] non
7 Buy['SPG'] low
8 Buy['TUR'] low
_conname _con.iis
1 Diet['A'] non
2 Diet['C'] non
3 Diet['B1'] non
4 Diet['B2'] mem
5 Diet['NA'] mem
6 Diet['CAL'] non
This information indicates that the IIS consists of five lower and two upper bounds on the variables, plus the constraints providing the lower bound on B2 and the upper bound on NA in the diet. Together these restrictions have no feasible solution, but dropping any one of them will permit a solution to be found to the remaining ones. If dropping the bounds is not of interest, then you may want to list only the constraints in the IIS. A print statement produces a concise listing:
cons_in_iis = cons_iis_df[cons_iis_df["_con.iis"] != "non"]
print(cons_in_iis)
_conname _con.iis
4 Diet['B2'] mem
5 Diet['NA'] mem
You could conclude in this case that, to avoid violating the bounds on amounts purchased, you might need to accept either less vitamin B2 or more sodium, or both, in the diet. Further experimentation would be necessary to determine how much less or more, however, and what other changes you might need to accept in order to gain feasibility. (A linear program can have several irreducible infeasible subsets, but Gurobi’s IIS finding algorithm detects only one IIS at a time.)