Mathematical Programming (Python)

# Import the Pyomo library, which can be installed via pip, conda, or from GitHub: https://github.com/Pyomo/pyomo
import pyomo.environ as pyo
# Import Matplotlib to generate plots
import matplotlib.pyplot as plt
# Import numpy and scipy for certain numerical calculations below
import numpy as np
from scipy.special import gamma
import math
import os

# Generate the feasible region plot of this problem

# Define meshgrid for feasible region
d = np.linspace(-0.5,3.5,300)
x1,x2 = np.meshgrid(d,d)

# Define the lines for the constraints
x = np.linspace(x1.min(), x1.max(), 2000)
# x2 <= x1 + 2
x21 = x + 2
# 8*x1 + 2*x2 <= 17
x22 = (17-8*x)/2.0
# Objective: max 5.5*x1 + 2.1*x2
# generate heatmap from objective function
objective_heatmap = np.fromfunction(lambda i, j:5.5 * i + 2.1 * j,(300,300), dtype=float)

# Plot feasible region
fig, ax = plt.subplots()
feas_reg = ax.imshow( (
    (x1>=0) &           # Bound 1 
    (x2>=0) &           # Bound 2
    (x2 <= x1 + 2) &    # Constraint 1
    (8*x1 + 2*x2 <= 17) # Constraint 2
    ).astype(int) 
    * objective_heatmap, # objective function
                extent=(x1.min(),x1.max(),x2.min(),x2.max()),origin="lower", cmap="OrRd", alpha = 1)

# Make plots of constraints
ax.plot(x, x21, label=r'$x_2 \leq x_1 + 2$')
ax.plot(x, x22, label=r'$8x_1 + 2x_2 \leq 17$')

# Nonnegativitivy constraints
plt.plot(x, np.zeros_like(x), label=r'$x_2 \geq 0$')
plt.plot(np.zeros_like(x), x, label=r'$x_1 \geq 0$')

# plot formatting
plt.title("LP Feasible Region \n max $z = 5.5x_{1} + 2.1x_{2}$",fontsize=10,y=1)
plt.xlim(x1.min(),x1.max())
plt.ylim(x2.min(),x2.max())
plt.legend(loc='upper right', prop={'size': 6})
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.grid(alpha=0.2)

# Define the model
model = pyo.ConcreteModel(name='Simple example LP')
# Define the variables
model.x = pyo.Var([1,2], domain=pyo.NonNegativeReals)
# Define the objective function
def _obj(m):
    return 5.5*m.x[1] + 2.1*m.x[2]

model.obj = pyo.Objective(rule = _obj, sense=pyo.maximize)
# Define the constraints
def _constraint1(m):
    return m.x[2] <= m.x[1] + 2

def _constraint2(m):
    return 8*m.x[1] + 2*m.x[2] <= 17
    
model.Constraint1 = pyo.Constraint(rule = _constraint1)

model.Constraint2 = pyo.Constraint(rule = _constraint2)
# Print the model
model.pprint()

# Define the solvers GLPK and CBC
if IN_COLAB:
    opt_glpk = pyo.SolverFactory('glpk', executable='/usr/bin/glpsol')
    opt_cbc = pyo.SolverFactory('cbc', executable='/usr/bin/cbc')
else:
    opt_glpk = pyo.SolverFactory('glpk')
    opt_cbc = pyo.SolverFactory('cbc')

for solver_name, solver in [('GLPK', opt_glpk), ('CBC', opt_cbc)]:
    if not solver.available(exception_flag=False):
        raise RuntimeError(f"{solver_name} not found. Run the GLPK/CBC install cell above.")
    print(f"{solver_name}: OK")
# Here we could use another solver, e.g. gurobi or cplex
# opt_gurobi = pyo.SolverFactory('gurobi')

# Each solve checks termination and feasibility, and compares the objective
# with the tutorial's known optimum when supplied (14.08 for this LP).
# A failed check stops execution. BONMIN's nonconvex example omits the
# objective check because it cannot certify a global optimum.
def check_solution(model, results, expected_objective=None):
    """Check termination, feasibility, and the tutorial's known objective.

    Absolute tolerance 1e-6 in each model expression's units accommodates
    solver feasibility slack without hiding a wrong solution.
    Omit expected_objective for BONMIN on the nonconvex model, where it
    cannot certify a global optimum.
    """
    pyo.assert_optimal_termination(results)
    tolerance = 1e-6
    for variable in model.component_data_objects(pyo.Var, active=True):
        value = pyo.value(variable)
        assert math.isfinite(value), "Non-finite decision variable"
        if variable.lb is not None:
            assert value >= variable.lb - tolerance, "Variable below its lower bound"
        if variable.ub is not None:
            assert value <= variable.ub + tolerance, "Variable above its upper bound"
        if variable.is_integer():
            assert abs(value - round(value)) <= tolerance, "Non-integer decision"
    for constraint in model.component_data_objects(pyo.Constraint, active=True):
        value = pyo.value(constraint.body)
        assert math.isfinite(value), "Non-finite constraint residual"
        if constraint.has_lb():
            assert value >= pyo.value(constraint.lower) - tolerance, "Constraint lower bound violated"
        if constraint.has_ub():
            assert value <= pyo.value(constraint.upper) + tolerance, "Constraint upper bound violated"
    if expected_objective is not None:
        assert math.isclose(pyo.value(model.obj), expected_objective,
                            rel_tol=0, abs_tol=tolerance), "Unexpected objective"

# Here we solve the optimization problem, the option tee=True prints the solver output
result_obj = opt_glpk.solve(model, tee=True)
check_solution(model, result_obj, expected_objective=14.08)

# Display solution of the problem
model.display()

# Optimal solution LP
ax.scatter(1.3,3.3,color='r', label='optimal solution LP')
ax.get_legend().remove()
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
fig.canvas.draw()
fig

# We obtain the same solution with CBC
result_obj = opt_cbc.solve(model, tee=False)
check_solution(model, result_obj, expected_objective=14.08)
model.display()

# Here we solve the optimization problem, the option tee=True prints the solver output
result_obj_ipopt = opt_ipopt.solve(model, tee=True)
check_solution(model, result_obj_ipopt, expected_objective=14.08)
model.display()

# Define grid for integer points
x1_int, x2_int = np.meshgrid(range(math.ceil(x1.max())), range(math.ceil(x2.max())))
idx = ((x1_int>=0) & (x2_int <= x1_int + 2) & (8*x1_int + 2*x2_int <= 17) & (x2_int>=0))
x1_int, x2_int = x1_int[idx], x2_int[idx]
ax.scatter(x1_int,x2_int,color='k', label='integer points')

# Plotting optimal solution IP
# plt.title("LP Feasible Region \n max $z = 5.5x_{1} + 2.1x_{2}$",fontsize=10,y=1)
ax.set_title("ILP Feasible Region \n max $z = 5.5x_{1} + 2.1x_{2}$")
ax.get_legend().remove()
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
fig.canvas.draw()
fig

# Define the integer model
model_ilp = pyo.ConcreteModel(name='Simple example IP, 47-779 QuIP')
#Define the variables
model_ilp.x = pyo.Var([1,2], domain=pyo.Integers)
# Define the objective function
model_ilp.obj = pyo.Objective(rule = _obj, sense=pyo.maximize)
# Define the constraints
model_ilp.Constraint1 = pyo.Constraint(rule = _constraint1)

model_ilp.Constraint2 = pyo.Constraint(rule = _constraint2)
# Print the model
model_ilp.pprint()

# Here we solve the optimization problem, the option tee=True prints the solver output
result_obj_ilp = opt_cbc.solve(model_ilp, tee=False)
check_solution(model_ilp, result_obj_ilp, expected_objective=11.8)

model_ilp.display()

# Optimal solution ILP
ax.scatter(1,3,color='c', label='optimal solution ILP')
ax.get_legend().remove()
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
fig.canvas.draw()
fig

fig2, ax2 = plt.subplots()
n = np.arange(1,100,1)
ax2.plot(n,np.exp(n*np.log(2)), label=r'$2^n$')
ax2.plot(n,gamma(n), label=r'$n!$')

ax2.plot(n,3.154E16*np.ones_like(n), 'g--', label=r'ns in a year')
ax2.plot(n,4.3E26*np.ones_like(n), 'k--', label=r'age of the universe in ns')
ax2.plot(n,6E79*np.ones_like(n), 'r--', label=r'atoms in the universe')
plt.yscale('log')

plt.legend()
plt.xlabel(r'$n$')
plt.ylabel('Possible solutions')

# Define grid for integer points
# Remove the previous heatmap if this cell is re-run out of order.
try:
    feas_reg.remove()
except (ValueError, AttributeError, NameError):
    pass  # already removed or not yet created
feas_reg = ax.imshow( (
    (x1>=0) &           # Bound 1 
    (x2>=0) &           # Bound 2
    (x2 <= x1 + 2) &    # Constraint 1
    (8*x1 + 2*x2 <= 17) & # Constraint 2
    ((x2-1)**2 <= 2-x1) # Nonlinear constraint
    ).astype(int) 
    * objective_heatmap, # objective function, 
                extent=(x1.min(),x1.max(),x2.min(),x2.max()),origin="lower", cmap="OrRd", alpha = .8)


x1nl = 2- (x - 1)**2
# Nonlinear constraint
nl_const = ax.plot(x1nl, x, label=r'$(x_2-1)^2 \leq 2-x_1$')

# Plotting optimal solution INLP
ax.set_title("INLP Feasible Region \n max $z = 5.5x_{1} + 2.1x_{2}$")
ax.get_legend().remove()
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
fig.canvas.draw()
fig

# Define the integer model
model_cinlp = pyo.ConcreteModel(name='Simple example convex INLP, 47-779 QuIP')
#Define the variables
model_cinlp.x = pyo.Var([1,2], domain=pyo.Integers)
# Define the objective function
model_cinlp.obj = pyo.Objective(rule = _obj, sense=pyo.maximize)
# Define the constraints
model_cinlp.Constraint1 = pyo.Constraint(rule = _constraint1)

model_cinlp.Constraint2 = pyo.Constraint(rule = _constraint2)

model_cinlp.Constraint3 = pyo.Constraint(expr = (model_cinlp.x[2]-1)**2 <= 2 - model_cinlp.x[1])


# Print the model
model_cinlp.pprint()

# Define the solver BONMIN
if IN_COLAB:
    opt_bonmin = pyo.SolverFactory('bonmin', executable='/content/solvers/bonmin')
else:
    opt_bonmin = pyo.SolverFactory('bonmin')
if not opt_bonmin.available(exception_flag=False):
    raise RuntimeError("BONMIN not found. Run the IDAES solver setup cell near the top of this notebook.")

# Here we solve the optimization problem, the option tee=True prints the solver output
result_obj_cinlp = opt_bonmin.solve(model_cinlp, tee=False)
check_solution(model_cinlp, result_obj_cinlp, expected_objective=9.7)
model_cinlp.display()

ax.scatter(1,2,color='m', label='optimal solution convex INLP')
ax.get_legend().remove()
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
fig.canvas.draw()
fig

# Define grid for integer points
# Remove the previous heatmap if this cell is re-run out of order.
try:
    feas_reg.remove()
except (ValueError, AttributeError, NameError):
    pass  # already removed or not yet created
feas_reg = ax.imshow( (
    (x1>=0) &           # Bound 1 
    (x2>=0) &           # Bound 2
    (x2 <= x1 + 2) &    # Constraint 1
    (8*x1 + 2*x2 <= 17) & # Constraint 2
    ((x2-1)**2 <= 2-x1) & # Nonlinear constraint 1
    ((x2-1)**2 >= x1+0.5) # Nonlinear constraint 2
    ).astype(int) 
    * objective_heatmap,  # objective function,
                extent=(x1.min(),x1.max(),x2.min(),x2.max()),origin="lower", cmap="OrRd", alpha = .8)


x1nl = -1/2 + (x - 1)**2
# Nonlinear constraint
nl_const = ax.plot(x1nl, x, label=r'$(x_2-1)^2 \geq x_1 + 1/2$')

ax.get_legend().remove()
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
fig.canvas.draw()
fig

# Define the integer model
model_ncinlp = pyo.ConcreteModel(name='Simple example non-convex INP, 47-779 QuIP')
#Define the variables
model_ncinlp.x = pyo.Var([1,2], domain=pyo.Integers)
# Define the objective function
model_ncinlp.obj = pyo.Objective(rule = _obj, sense=pyo.maximize)
# Define the constraints
model_ncinlp.Constraint1 = pyo.Constraint(rule = _constraint1)

model_ncinlp.Constraint2 = pyo.Constraint(rule = _constraint2)

model_ncinlp.Constraint3 = pyo.Constraint(expr = (model_ncinlp.x[2]-1)**2 <= 2 - model_ncinlp.x[1])

model_ncinlp.Constraint4 = pyo.Constraint(expr = (model_ncinlp.x[2]-1)**2 >= 1/2 + model_ncinlp.x[1])


# Print the model
model_ncinlp.pprint()

# Trying to solve the problem with BONMIN we might obtain the optimal solution, but we have no guarantees
result_obj_ncinlp = opt_bonmin.solve(model_ncinlp, tee=False)
check_solution(model_ncinlp, result_obj_ncinlp)
model_ncinlp.display()

# Define the solver COUENNE
if IN_COLAB:
    opt_couenne = pyo.SolverFactory('couenne', executable='/content/solvers/couenne')
else:
    opt_couenne = pyo.SolverFactory('couenne')
if not opt_couenne.available(exception_flag=False):
    raise RuntimeError("COUENNE not found. Run the IDAES solver setup cell near the top of this notebook.")

# Trying to solve the problem with global MINLP solver COUENNE
result_obj_ncinlp = opt_couenne.solve(model_ncinlp, tee=False)
check_solution(model_ncinlp, result_obj_ncinlp, expected_objective=4.2)
model_ncinlp.display()

ax.scatter(0,2,color='g', label='optimal solution nonconvex INLP')
ax.get_legend().remove()
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
fig.canvas.draw()
fig

# =============================================================================
# EXERCISE 1: Change the objective sense
# =============================================================================
# Create a scratch version of the linear model that minimizes the same linear expression instead of maximizing it. Compare the new solution with the original optimum and identify which active constraints changed.
#
# Hint: Keep the feasible region unchanged so only the objective direction changes.
#
# Your code here:
# =============================================================================

# =============================================================================
# EXERCISE 2: Add an integrality constraint
# =============================================================================
# Choose one continuous decision variable and restrict it to integer values. Re-solve the model and describe whether the objective value changes because the feasible set became smaller.
#
# Hint: Use a copied model or a fresh variable declaration so the original demonstration remains reproducible.
#
# Your code here:
# =============================================================================

# =============================================================================
# EXERCISE 3: Explain a new constraint geometrically
# =============================================================================
# Add one additional linear constraint that cuts off part of the original feasible region. Sketch or describe which side of the line remains feasible and predict the effect before solving.
#
# Hint: Pick a simple bound such as a weighted sum of the two plotted variables.
#
# Your code here:
# =============================================================================