Introduction to Linear and Integer Programming¶

Description: Basic introduction to linear programming and AMPL via a lemonade stand example

Tags: ampl-lecture, amplpy, ampl, introduction, linear programming, lemonade stand

Notebook author: Gyorgy Matyasfalvi <gyorgy@ampl.com>

References:

AMPL a Modeling Language for Mathematical Programming – Robert Fourer et al.
Linear Programming (Foundations and Extensions) – Robert J. Vanderbei
Introduction to Linear Optimization – Dimitris Bertsimas et al.
Convex Analysis and Optimization – Dimitri P. Bertsekas et al.

# Install dependencies
%pip install -q amplpy pandas numpy matplotlib

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["open", "gurobi"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

# Import all necessary libraries
import matplotlib.pyplot as plt
import numpy as np
from math import log, cos, sin, pi, sqrt

Linear problems¶

Convex (aka continuous) linear problems represent one of the more tractable classes of optimization problems. A plethora of algorithms and optimization techniques exist that adeptly tackle such challenges. Fortunately, many real-world situations can be effectively represented as continuous linear programs. Examples include specific production planning scenarios, multi-period electric power capacity planning, certain scheduling dilemmas, path selection in communication networks, pattern classification, and optimal control of linear systems. A linear programming problem (LP) in standard form is articulated as (all LPs can be expressed in the standard form):

$ \begin{equation} \begin{array}{lll} &\min & c^T x & \ &\textrm{s.t.} & Ax = b & \ & & x \geq 0, &
\end{array} \tag{4} \end{equation} $

where $c^T \in \mathbb{R}^n$, $b \in \mathbb{R}^m$, and $A \in \mathbb{R}^{m \times n}$.

Integer linear problems¶

While many real-world applications can be formulated as purely linear programs, practical scenarios often necessitate our decision variables to be limited to integer values. This is especially true when modeling decisions about quantities, such as the number of products to manufacture or the number of workers to assign to a particular task. Linear integer programs are commonly abbreviated as IPs. A pure integer program can be articulated as follows:

$ \begin{equation} \begin{array}{lll} &\min & c^T x & \ &\textrm{s.t.} & Ax = b & \ & & x \geq 0 & \ & & x \in \mathbb{Z}, &
\end{array} \tag{5} \end{equation} $

where $c^T \in \mathbb{R}^n$, $b \in \mathbb{R}^m$, and $A \in \mathbb{R}^{m \times n}$.

Solver Expectations for Linear Programming Problems (LPs)¶

In the realm of Linear Programming (LP), given feasible conditions and sufficient computation time, any solver is expected to find a globally optimal solution. This principle applies to non-linear solvers like Knitro as well, which are capable of efficiently solving LPs since linear problems are inherently a subset of non-linear problems. The preference for using dedicated linear solvers stems from their specially designed algorithms, which are highly optimized for speed when addressing LPs.

When we venture into the territory of non-linear problems, the proficiency of linear solvers is often constrained, particularly with non-linearities present in both the constraints and the objective function. However, if non-linearities are confined solely to the objective function, certain algorithms developed for linear optimization may be successfully adapted.

Among non-linear problems, only a select group, such as convex quadratic programs, are within the purview of advanced linear solvers. Notably, Gurobi, which was once confined to solving convex quadratically constrained problems, has broadened its horizons. Its latest release (still in beta) includes the functionality of a global solver, enabling it to address a wider spectrum of optimization challenges.

Although HiGHS was unable to resolve the ‘difficult’ optimization task, the new version of Gurobi will present a promising alternative.

EXERCISE 5¶

In the code cell provided below, adjust the settings of the solver to ensure that it can find the optimal solution. You will need to select a solver that can solve general nonlinear problems.
Modify the program such that if the optimal solution is successfully located, it should display the message ‘Success: Problem solved to optimality!’.

ampl = AMPL()  # Instantiate an AMPL object
ampl.eval(
    r"""
            var x;                                              # Declare our variable
            minimize difficult: (x-1)^2 - log(x+1);             # Define our objective function
          """
)

ampl.option["solver"] = "highs"  # Select the solver
ampl.solve()  # Attempt to solve
print(
    ampl.option["solve_result_table"]
)  # Print the solve result table, this will inform us of the various solution codes.
result_df = ampl.get_data("solve_result_num").to_pandas()  # Retrieve solution status
# Check if the problem was solved if not print warning
srn = result_df["solve_result_num"].iloc[0]
if srn != 0:
    print(f"Warning:\tProblem not solved to optimality!\n\t\tsolve_result_num = {srn}")

HiGHS 1.6.0: 
------------ WARNINGS ------------
WARNING:  "PLApprox"
  An expression of type 'LogConstraint' has been
piecewise-linearly approximated. Set cvt:plapprox:reltol
to control precision (currently 0.010000).
WARNING:  "PLApproxDomain"
  Argument domain of a 'LogConstraint'
has been reduced to [0.000001, 1000000.000000] for numerical reasons
(partially controlled by cvt:plapprox:domain.)
------------ WARNINGS ------------
WARNING:  "PLApprox"
  An expression of type 'LogConstraint' has been
piecewise-linearly approximated. Set cvt:plapprox:reltol
to control precision (currently 0.010000).
WARNING:  "PLApproxDomain"
  Argument domain of a 'LogConstraint'
has been reduced to [0.000001, 1000000.000000] for numerical reasons
(partially controlled by cvt:plapprox:domain.)
  Error -1 for call Highs_run(lp())

0	solved
100	solved?
200	infeasible
300	unbounded
400	limit
500	failure

Warning:	Problem not solved to optimality!
		solve_result_num = 500

SOLUTION¶

Pick a non-linear solver such as Knitro, Ipopt, or Lindoglobal. Because the problem is convex all non-linear solvers will find the optimal solution.
Modify the if-statement to print a “success” message if the problem was solved to optimality.

ampl = AMPL()  # Instantiate an AMPL object
ampl.eval(
    r"""
            var x;                                              # Declare our variable
            minimize difficult: (x-1)^2 - log(x+1);             # Define our objective function
          """
)

ampl.option["solver"] = "knitro"  # Select the solver
ampl.solve()  # Attempt to solve
print(
    ampl.option["solve_result_table"]
)  # Print the solve result table, this will inform us of the various solution codes.
result_df = ampl.get_data("solve_result_num").to_pandas()  # Retrieve solution status
# Check if the problem was solved if not print warning
srn = result_df["solve_result_num"].iloc[0]
if srn != 0:
    print(f"Warning:\tProblem not solved to optimality!\n\t\tsolve_result_num = {srn}")
else:
    print("Success:\tProblem solved to optimality!")

Cannot find "knitro"

0	solved
100	solved?
200	infeasible
300	unbounded
400	limit
500	failure

Warning:	Problem not solved to optimality!
		solve_result_num = 999

Solver Expectations for Integer Programming Problems¶

Solving integer programming problems is often exponentially more challenging than their linear programming counterparts. The most formidable among these are the Integer Non-Linear Programs (MINLPs), which can be exceedingly complex to model and solve—sometimes even involving the complex plane.

Even the linear versions of these problems, referred to as Mixed-Integer Linear Programs (MILPs), are significantly more difficult to solve compared to standard linear programs (LPs).

Most solvers utilize a method called branch-and-bound, which, in the worst-case scenario, may need to explore an exponentially large number of subproblems proportional to the number of integer variables. The branch-and-bound algorithm begins by solving the continuous relaxation of the integer problem. If this yields an integer solution, the problem is considered solved. However, should the solution contain non-integer values, the algorithm branches by introducing new constraints designed to eliminate the fractional values. For instance, if $ x = 2.5 $, the solver would create two new subproblems: one with the constraint $ x \leq 2 $ and another with $ x \geq 3 $. By doing so, the solver ensures that neither of these branches can return the fractional value $ x = 2.5 $. The process continues, solving each subproblem and branching further until an integer solution is found.

Considering that real-world problems can involve millions of variables, the enumeration tree (which contains all possible subproblem solutions) can become enormous and unwieldy.

Due to these complexities, optimization problems involving integers often require specialized, and typically commercial, solvers to tackle them efficiently.

Let’s now turn to formulating our first constrained optimization problem. We will start with a simple lemonade stand example and gradually evolve it into a more general production model.

Lemonade stand¶

Imagine you own a small lemonade stand and want to make the most profit possible. You offer two drinks: lemonade $(x_{\textrm{lemonade}})$, which requires two cups of water, one lemon and two cups of sugar per glass, and iced tea $(x_{\textrm{iced_tea}})$, which requires one tea bag, one cup of sugar, and two cups of water per glass. Thankfully, you have unlimited access to water for free.

However, you have limited “dry” ingredients to work with each day, including 10 lemons, 8 tea bags, and 20 cups of sugar. Fortunately, you receive a daily inheritance of ingredients from a relative, so you don’t have to worry about purchasing them.

Your goal is to decide how many glasses of each drink to make while staying within your daily ingredient limits to maximize your profit. You can sell each glass of lemonade for $1.50 and each glass of iced tea for $1.00. You also pay the city a daily permit fee of $2 to keep your stand open. Since you’re the only lemonade stand in town, you won’t be spending any money on iced tea or lemonade anywhere else, and will sell all that you make.

To formulate this as an integer programming problem, we can use the following model:

$ \begin{equation} \begin{array}{lcrcrl} &\textrm{Objective:} & \max & 1.5x_{\textrm{lemonade}} & + & x_{\textrm{iced_tea}} & - & 2 \ &\textrm{Subject to:} & & & & & & \ & & & x_{\textrm{lemonade}} & & & \leq & 10 \ & & & & & x_{\textrm{iced_tea}} & \leq & 8 \ & & & 2x_{\textrm{lemonade}} & + & x_{\textrm{iced_tea}} & \leq & 20 \ & & & & & & & \ & & & x_{\textrm{lemonade}} & & & \geq & 0 \ & & & & & x_{\textrm{iced_tea}} & \geq & 0 \ & & & x_{\textrm{lemonade}} & , & x_{\textrm{iced_tea}} & \in & \mathbb{Z} \end{array} \end{equation} $

EXERCISE 6¶

With the integer programming problem outlined above, your task is to construct an AMPL model. We’re familiar with declaring variables and objective functions, yet specifying additional attributes for variables, such as non-negativity or integrality, and defining constraints is new territory. Fear not, for integrating these aspects into AMPL is a straightforward process. Refer to the provided sections ‘The AMPL language’ and ‘Model formulation in AMPL’ below to guide you in crafting your first AMPL model.
Below, create a code cell and begin with the %%writefile lemonade.mod magic command. Follow with the formulation of your AMPL model.
Next, in a separate code cell, initialize your AMPL environment and load your model using the ampl.read("lemonade.mod") command.
Develop and execute the necessary code to calculate the optimal amounts of lemonade and iced tea to prepare in order to maximize profits.
Finally, output the quantities of lemonade and iced tea that should be produced, along with the projected profit.

The AMPL language¶

Section A.1 of the appendix of the AMPL book provides a detailed discussion of AMPL’s lexical rules. We will mention the most important below:

AMPL models involve variables, constraints, and objectives, expressed with the help of sets and parameters. These are called model entities, each of which has an alphanumeric name such as lemonade. Upper-case letters are distinct from lower-case so lemonade and Lemonade would correspond to two separate entities.

The most basic format for declaring any entity in AMPL consists of the entity type (variable, constraint, etc.), the name you give it, and a semicolon at the end of the statement. For example,
```
var lemonade;
```

AMPL is a “declarative” language, meaning that entities used in expressions must be declared before they can be used.

An expression is similar to a mathematical equation or formula and typically involves variables, constants, and mathematical operators such as addition, subtraction, multiplication, and division. Expressions can also include inequalities, logical operators, and other constructs that allow for complex modeling. In AMPL, expressions are used extensively in the definition of the objective function and constraints of a mathematical program. For example:
```
maximize profit: 1.5*lemonade + iced_tea - fee;
```

Comments in AMPL files begin with: #. For example,

# Define the objective function
maximize profit: 1.5*lemonade + iced_tea - fee;

AMPL ignores white space. For example,

var lemonade >= 0;
var iced_tea >= 0; 

is equivalent to

var lemonade >= 0; var iced_tea >= 0;

Each AMPL statement should end with a semicolon, with the exception of include commands, which we will discuss in more detail later.

Model formulation in AMPL¶

The following sections will cover the basic syntax in AMPL for declaring these entities. We will begin with an overview of how to declare variables, then move on to defining the objective function and constraints.

Variable declarations¶

To declare our variables in AMPL, we first use the var keyword followed by the variable name, optional attributes, and a semicolon.


    var name alias_opt indexing_opt attributes_opt ;

    attributes: 
        binary
        integer
        symbolic
        >= expr
        <= expr
        := expr
        default expr
        = expr
        coeff indexing_opt constraint expr
        cover indexing_opt constraint
        obj indexing_opt objective expr
        in sexpr
        suffix sufname expr

Optional clauses marked with opt subscript such as alias or indexing are discussed in detail in the AMPL book.

As a brief note, expr is shorthand for expression, and sexpr stands for set expression.

It is also worth mentioning that by default, AMPL initializes all variables to zero unless an initial value is provided via an attribute (e.g. var x := 1; in which case the initial value of x becomes 1).

Parameter declarations¶

To declare a parameter in AMPL, we use the keyword param, followed by the parameter name, an optional initial value assignment, and a semi-colon.


   param name alias_opt indexing_opt attributes_opt ;

    attributes: 
        binary
        integer
        symbolic
        relop expr
        default expr
        = expr
        in sexpr

    relop: 
        <  <=  =  ==  !=  <>  >  >=

Optional clauses marked with opt subscript such as alias or indexing are discussed in detail in the AMPL book.

As a brief note, expr is shorthand for expression, sexpr stands for set expression, and relop is short for relational operator.

Objective declarations¶

To specify an objective function in AMPL, we first have to specify the direction of optimization via the keywords maximize or minimize, which is followed by a name, a colon, an expression involving our variables, and a semicolon.


    maximize name alias_opt indexing_opt  
        : expression [ suffix-initialization ] ;


    minimize name alias_opt indexing_opt  
        : expression [ suffix-initialization ] ;

Optional clauses, such as suffix-initialization (marked with [ ]) and alias or indexing (marked with opt subscript), are discussed in detail in the AMPL book.

Reminder: Variables must be declared before they can be used in an expression.

Constraint declarations¶

To declare our constraints in AMPL, we first use the keywords subject to, which is followed by a name, a colon, an expression involving our variables, and a semicolon.

Just like the objective function, we can give constraints an informative name and then provide the appropriate expression to describe the constraint.


    [ subject to ] name alias_opt indexing_opt 
        [:=  initial_dual ] [ default initial_dual ] 
        : constraint expression [ suffix-initialization ] ;

Constraint expressions can take two forms one-sided and two-sided.

One-sided constraint expressions must be specified in one of the following forms:

$$ \begin{equation} \begin{array}{lcl} \textrm{expression} & <= & \textrm{expression} \ \textrm{expression} & = & \textrm{expression} \ \textrm{expression} & >= & \textrm{expression} \end{array} \end{equation} \tag{2} $$

Whereas two-sided constraint expressions must be declared in one of the following forms:

$$ \begin{equation} \begin{array}{lcccr} \textrm{constant expression} & <= & \textrm{expression} & <= & \textrm{constant expression} \ \textrm{constant expression} & >= & \textrm{expression} & >= & \textrm{constant expression} \ \end{array} \end{equation} \tag{3} $$

A constant expression is an expression that does not contain any variables.

Optional clauses, such as initial_dual or suffix-initialization (marked with [ ]) and alias or indexing (marked with opt subscript), are discussed in detail in the AMPL book.

The subject to keyword used to declare constraints can be shortened to s.t. or even omitted entirely (not recommended).

SOLUTION¶

The %%writefile lemonade.mod magic will create the file lemonade.mod in the directory this Jupyter Notebook is located in.

%%writefile lemonade.mod
# Define the decision variables
var lemonade integer, >= 0;
var iced_tea integer, >= 0;

# Define the parameter(s)
param fee = 2;

# Define the objective function
maximize profit: 1.5*lemonade + iced_tea - fee;

# Define the constraints
subject to lemon_constraint: lemonade <= 10;
subject to tea_bag_constraint: iced_tea <= 8;
subject to sugar_constraint: 2*lemonade + iced_tea <= 20;

Writing lemonade.mod

ampl = AMPL()  # Instantiate an AMPL object

ampl.read("lemonade.mod")  # Read the model file

ampl.option["solver"] = "gurobi"  # Select the solver

ampl.solve()  # Attempt to solve

df_sol = ampl.get_data("lemonade, iced_tea, profit").to_pandas()  # Get solution

print(df_sol.T)  # Print solution

Gurobi 11.0.0:Gurobi 11.0.0: optimal solution; objective 15
0 simplex iterations
 
           0
lemonade   6
iced_tea   8
profit    15

Mulled wine stand¶

Your lemonade stand works well during the summer months, but as winter approaches, you sense the need to switch your business model. Thankfully, your inheritance of ingredients can be adapted from lemons to wine and spices, and you decide to open a mulled wine and hot tea stand. As before, your goal is to maximize your daily profit by deciding how many cups of mulled wine and hot tea to sell, subject to the constraint that you have a limited amount of wine, spices, sugar, and tea bags available to you each day.

In this example you will offer two drinks: mulled wine $(x_{\textrm{mulled_wine}})$, which requires two cups of wine, two tablespoons of spice and four tablespoons of sugar per glass, and hot tea $(x_{\textrm{hot_tea}})$, which requires one tea bag, two tablespoons of sugar, and two cups of water per glass.

Thankfully, you still have access to unlimited water, however your daily inheritance of ingredients has changed to: 12 tablespoons of spices, 8 cups of tea, 30 tablespoons of sugar, and 15 cups of wine.

Your goal will again be to maximize profit by determining how many glasses of each drink to make given your limited daily ingredients. You can sell each glass of wine for $2.00 and each glass of tea for $1.50. Since you have cornered the market on drink stands in town, you can expect to sell all the drinks that you make and remember that your ingredient list is limited but costs nothing.

To formulate this as an integer programming problem, we can use the following model:

$$ \begin{equation} \begin{array}{lcrcrl} &\textrm{Objective:} & \max & 2x_{\textrm{mulled_wine}} & + & 1.5x_{\textrm{hot_tea}} & - & 2 \ &\textrm{Subject to:} & & & & & & \ & & & 2x_{\textrm{mulled_wine}} & & & \leq & 12 \ & & & & & x_{\textrm{hot_tea}} & \leq & 8 \ & & & 4x_{\textrm{mulled_wine}} & + & 2x_{\textrm{hot_tea}} & \leq & 30 \ & & & 2x_{\textrm{mulled_wine}} & & & \leq & 15 \ & & & & & & & \ & & & x_{\textrm{mulled_wine}} & & & \geq & 0 \ & & & & & x_{\textrm{hot_tea}} & \geq & 0 \ & & & x_{\textrm{mulled_wine}} & , & x_{\textrm{hot_tea}} & \in & \mathbb{Z} \end{array} \end{equation} \tag{4} $$

EXERCISE 7¶

With the integer programming problem outlined above, your task is to construct an AMPL model.
Below, create a code cell and begin with the %%writefile mulled_wine.mod magic command. Follow with the formulation of your AMPL model.
Next, in a separate code cell, initialize your AMPL environment and load your model using the ampl.read("mulled_wine.mod") command.
Develop and execute the necessary code to calculate the optimal amounts of mulled wine and hot tea to prepare in order to maximize profits.
Finally, output the quantities of mulled wine and hot tea that should be produced, along with the projected profit.

SOLUTION¶

%%writefile mulled_wine.mod
# Define the decision variables
var mulled_wine integer, >= 0;
var hot_tea integer, >= 0;

# Define the parameter(s)
param fee = 2;

# Define the objective function
maximize profit: 2*mulled_wine + 1.5*hot_tea - fee;

# Define the constraints
subject to spice_constraint: 2*mulled_wine <= 12;
subject to tea_bag_constraint: hot_tea <= 8;
subject to sugar_constraint: 4*mulled_wine + 2*hot_tea <= 30;
subject to wine_constraint: 2*mulled_wine <= 15;

Writing mulled_wine.mod

ampl = AMPL()  # Instantiate an AMPL object

ampl.read("mulled_wine.mod")  # Read the model file

ampl.option["solver"] = "highs"  # Select the solver

ampl.solve()  # Attempt to solve

df_sol = ampl.get_data("mulled_wine, hot_tea, profit").to_pandas()  # Get solution

print(df_sol.T)  # Print solution

HiGHS 1.6.0:HiGHS 1.6.0: optimal solution; objective 16.5
1 simplex iterations
1 branching nodes
 
                0
mulled_wine   4.0
hot_tea       7.0
profit       16.5

Hydrothermal Scheduling Problem with Conic Programming

Introduction to Mathematical Optimization