Lesson_02.ipynb

Kernel: Python 3 (system-wide)

# EXECUTE FIRST

# computational imports
from pyomo.environ import *
import numpy as np
import pandas as pd

# plotting imports
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("darkgrid")

# NOTE - If the videos are too small, you can change the dimensions here and re-run this cell.
def play_video(vid_name, w=640, h=360):
    from IPython.display import display, IFrame
    from IPython.core.display import HTML
    vid_path = "https://media.uwex.edu/content/ds/ds775_r19/" + vid_name + "/index.html"
    print(vid_path)
    hlink = '<a href = ' + vid_path + ' target = """_blank""">Open video in new tab</a>'
    display(IFrame(vid_path, width=w, height=h))
    display(HTML(hlink))

Lesson 02: Sensitivity Analysis and Abstract Modeling

Abstract Modeling for LP

To tackle larger linear programs it isn't practical to type out all of the variables, the terms in the objective function, or the individual constraints. We need to move from a concrete model like this:

to something more abstract like this:

We'll still use the ConcreteModel object in Pyomo, but we'll start each problem by first laying out the data and then adding indexed variables to the model object. The big idea in writing models that are generalizable or "abstract" is to separate the data from the model so that swapping in different or larger data is simple.

Below we present several examples of generalizing models so that the model can be easily adapted to different problems. Study these examples and use them to guide you in the homework. Rather than give detailed videos for these examples we've tried to carefully write the necessary information. Brief videos are provided to highlight a point or two.

To make the code easier to read and debug we use lists of meaningful variable names instead of simply using integers to index our variables. We use those variable names as keys in dictionaries to look up values of model coefficients. At the very bottom of this notebook you'll find a bit of material about Python dictionaries. However, if you want additional help you should seek out a tutorial. Here is an example tutorial, but there are many others. Feel free to share any resources you find on Piazza.

Wyndor - Abstract Formulation

Problem Description

Again, here is the Wyndor model from the textbook. $Z$ is the profit in thousands of dollars. $d$ and $w$ are the batches of doors and windows, respectively. The constraints, in order, represent the production capacities of Plants 1, 2, and 3.

Maximize $Z = 3 d + 5 w$

Subject to:

$\begin{array}{ccccc} d & & & \leq & 4 \\ & & 2w & \leq & 12 \\ 3d & + & 2w & \leq & 18 \end{array}$

$d \geq 0$, $w \geq 0$

Mathematical Formulation

We want to move from the above, very concrete, implementation to a more abstract representation which can be generalized to larger problems.

The index sets:

Let $Pr$ represent the set of products. Symbolically we write $Pr = \{ \mbox{doors}, \mbox{windows} \}$.

Let $Pl$ be set of plants so $Pl = \{ \mbox{Plant1}, \mbox{Plant2}, \mbox{Plant3} \}.$

In Python we can represent these index sets as any iterable object. While it is possible to use Python sets here, we usually use Python lists because those are both iterable and subscriptable. In Python the index sets look like this:

products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']

The objective function:

For each product $pr$ let $x_{pr}$ be the number of batches of that product to produce and let $c_{pr}$ be the profit rate per batch of $pr$ that is produced. The objective function could be written as $$Z = c_{\mbox{doors}} x_{\mbox{doors}} + c_{\mbox{windows}} x_{\mbox{windows}}.$$

More generally we can write: $$Z = \sum_{pr \in Pr} c_{pr} x_{pr}.$$

In Python this will look like:

sum(profit_rate[pr] * model.weekly_prod[pr] for pr in products)

The constraints:

To start, we can make all the constraints look the same: $$\begin{array}{ccccc} 1 x_{\mbox{doors}} & + & 0 x_{\mbox{windows}} & \leq & 4 \\ 0 x_{\mbox{doors}} & + & 2 x_{\mbox{windows}} & \leq & 12 \\ 3 x_{\mbox{doors}} & + & 2 x_{\mbox{windows}} & \leq & 18 \end{array}$$

Now we have three constraints, one for each plant, that has the form: $$(\mbox{hours per batch of doors}) * x_{\mbox{doors}} + (\mbox{hours per batch of windows}) * x_{\mbox{windows}} \leq \mbox{hourly capacity}.$$

If we let $h_{pl,pr}$ represent the hours needed for a batch of product $pr$ at plant $pl$ and let $a_{pl}$ be the number of hours available at plant $pl$, then we can express the constraints as:

for each plant $pl \in Pl$ $$\sum_{pr \in Pr} h_{pl,pr} x_{pr} \leq a_{pl}.$$

In Python this will look like:

sum(hours_per_batch[pl][pr] * model.weekly_prod[pr] for pr in products) <= hours_available[pl]

For comparison, the concrete and abstract models are shown side-by-side:

Concrete Model	Abstract Model
Maximize $Z = 3 d + 5 w$	Maximize $Z = \displaystyle \sum_{pr \in Pr} c_{pr} x_{pr}$
Subject to:

$\begin{array}{ccccc} d & & & \leq & 4 \\ & & 2w & \leq & 12 \\ 3d & + & 2w & \leq & 18 \end{array}$

Subject to:

$\displaystyle \sum_{pr \in Pr} h_{pl,pr} x_{pr} \leq a_{pl}, \mbox{ for each } pl \in Pl$

$d \geq 0$, $w \geq 0$ $x_{pr} \geq 0, \mbox{ for each } pr \in Pr$

In the Pyomo solutions below we'll share a few different approaches for the data structures used to both store the problem data and to construct the model. In each case we share some pluses and minuses to each approach. We offer a video for the first solution and code for all the approaches. After the four solution approaches we discuss some advantages and disadvantages to each approach. If four approaches is too many, then just focus on the second as it will work for all the problems we encounter.

Explanation of Pyomo Solution 1 (video)

The video below explains some parts of the solution code. If you're content to study the code on your own, then you don't need to watch it. There is a bit in the middle that shows you how to view the abstractly constructed constraints that is generally useful.

# execute this cell for video
play_video("ds775_lesson2-wyndor-abstract")

https://media.uwex.edu/content/ds/ds775_r19/ds775_lesson2-wyndor-abstract/index.html

Pyomo Solution 1

This solution is the one referenced in the video. Here the data and the model are separated quite well, but the dictionaries for the model data are fully typed out which isn't practical except for tiny models. For larger models the data would be loaded into lists or arrays from a file or database then the data would be parsed into dictionaries or Pandas series and dataframes. The advantage to those data structures is that they allow us to reference the data using strings that label the variables and constraints which makes the models easier to read and debug.

Further below we present 3 other solution approaches that illustrate different approaches to storing the problem data and using it in the model construction. After each solution we'll both print out the model and also note some pluses and minuses for that approach.

from pyomo.environ import *
import numpy as np
import pandas as pd

# abstract Wyndor

### PROBLEM DATA ###
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
profit_rate = {'Doors': 3, 'Windows': 5} 
hours_available = {'Plant1': 4, 'Plant2': 12, 'Plant3': 18}
hours_per_batch = {
    'Plant1': {
        'Doors': 1,
        'Windows': 0
    },
    'Plant2': {
        'Doors': 0,
        'Windows': 2
    },
    'Plant3': {
        'Doors': 3,
        'Windows': 2
    }
}

### MODEL CONSTRUCTION ###
#Declaration
model = ConcreteModel()

#Decision Variables
model.weekly_prod = Var(products, domain=NonNegativeReals)

#Objective
model.profit = Objective(expr=sum(profit_rate[pr] * model.weekly_prod[pr]
                               for pr in products),
                      sense=maximize)

#Constraints
model.capacity = ConstraintList()
for pl in plants:
    model.capacity.add(
        sum(hours_per_batch[pl][pr] * model.weekly_prod[pr]
            for pr in products) <= hours_available[pl])

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)



### OUTPUT ###

# note that we're using f-strings for output here which is a little different and cleaner than in the video
print(f"Maximum Profit = ${1000*model.profit():,.2f}")
for j in products:
    print(f"Batches of {j} = {model.weekly_prod[j]():.1f}")

Maximum Profit = $36,000.00
Batches of Doors = 2.0
Batches of Windows = 6.0

Here is what the constraints and objective function look like with the labeled variables:

# the objective function:
model.profit.pprint()

# the constraints:
model.capacity.pprint()

profit : Size=1, Index=None, Active=True
    Key  : Active : Sense    : Expression
    None :   True : maximize : 3*weekly_prod[Doors] + 5*weekly_prod[Windows]
capacity : Size=3, Index=capacity_index, Active=True
    Key : Lower : Body                                          : Upper : Active
      1 :  -Inf :   weekly_prod[Doors] + 0*weekly_prod[Windows] :   4.0 :   True
      2 :  -Inf : 0*weekly_prod[Doors] + 2*weekly_prod[Windows] :  12.0 :   True
      3 :  -Inf : 3*weekly_prod[Doors] + 2*weekly_prod[Windows] :  18.0 :   True

Notice that using meaningful labels for the variables and the indices make it easier to read and debug the model which is a big plus. The downside to the data approach above is that we had to type out the dictionaries that store all the model parameters. That wouldn't be feasible to larger models. We explore a better approach in the next solution.

Pyomo Solution 2 (Solutions 2 and 3 are the preferred approaches for homework)

The only difference here is that we have the data loaded into lists that might be similar to how the data would appear after it is loaded from an external source. The data is parsed from the lists into dictionaries to facilitate the model construction. The model construction and solution are identical to those in Solution 1. At the end of this notebook you'll find several examples of constructing dictionaries from lists. We also use this approach throughout our examples and homework solutions. You could use a similar approach for constructing dictionaries from any array-like data type.

### PROBLEM DATA ###

# load data
products = ['doors','windows']
            
          
plants = ['Plant1', 'Plant2', 'Plant3']
profit_rate_list = [3,5]
hours_avail_list = [4,12,18]
hours_per_batch_list = [ [1,0], [0,2], [3,2] ]

# parse lists into dictionaries
profit_rate = dict( zip( products, profit_rate_list) )
hours_available = dict( zip( plants, hours_avail_list) )
hours_per_batch = { pl: { pr: hours_per_batch_list[i][j] for j,pr in enumerate(products)} for i,pl in enumerate(plants)}

### MODEL CONSTRUCTION ###
#Declaration
model = ConcreteModel()

#Decision Variables
model.weekly_prod = Var(products, domain=NonNegativeReals)

#Objective
model.profit = Objective(expr=sum(profit_rate[pr] * model.weekly_prod[pr]
                               for pr in products),
                      sense=maximize)

#Constraints
model.capacity = ConstraintList()
for pl in plants:
    model.capacity.add(
        sum(hours_per_batch[pl][pr] * model.weekly_prod[pr]
            for pr in products) <= hours_available[pl])
    
model.pprint()

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###

# note that we're using f-strings for output here which is a little different and cleaner than in the video
print(f"Maximum Profit = ${1000*model.profit():,.2f}")
for j in products:
    print(f"Batches of {j} = {model.weekly_prod[j]():.1f}")

2 Set Declarations
    capacity_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {1, 2, 3}
    weekly_prod_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {'doors', 'windows'}

1 Var Declarations
    weekly_prod : Size=2, Index=weekly_prod_index
        Key     : Lower : Value : Upper : Fixed : Stale : Domain
          doors :     0 :  None :  None : False :  True : NonNegativeReals
        windows :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    profit : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 3*weekly_prod[doors] + 5*weekly_prod[windows]

1 Constraint Declarations
    capacity : Size=3, Index=capacity_index, Active=True
        Key : Lower : Body                                          : Upper : Active
          1 :  -Inf :   weekly_prod[doors] + 0*weekly_prod[windows] :   4.0 :   True
          2 :  -Inf : 0*weekly_prod[doors] + 2*weekly_prod[windows] :  12.0 :   True
          3 :  -Inf : 3*weekly_prod[doors] + 2*weekly_prod[windows] :  18.0 :   True

5 Declarations: weekly_prod_index weekly_prod profit capacity_index capacity
Maximum Profit = $36,000.00
Batches of doors = 2.0
Batches of windows = 6.0

Pyomo Solution 3

Solution approaches 2 and 3 are what you should use in the homework. We almost always use approach 2 in the lessons and solutions.

This is similar to Solution 2 except the data is parsed from the lists into Pandas dataframes and series to facilitate the model construction. The advantage to Pandas objects is that we can use our lists of labels as indices. In the model construction we only had to change how the values in hours_per_batch are referenced. This approach will work on most problems, but we'll encounter a problem in the next lesson with three index variables for which you'll need to use dictionaries.

### PROBLEM DATA ###

# load data
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
profit_rate_list = [3,5]
hours_avail_list = [4,12,18]
hours_per_batch_list = [ [1,0], [0,2], [3,2] ]

# parse lists into Pandas
profit_rate = pd.Series( profit_rate_list, index = products)
hours_available = pd.Series( hours_avail_list, index = plants)
hours_per_batch = pd.DataFrame( hours_per_batch_list, index = plants, columns = products)

### MODEL CONSTRUCTION ###
#Declaration
model = ConcreteModel()

#Decision Variables
model.weekly_prod = Var(products, domain=NonNegativeReals)

#Objective
model.profit = Objective(expr=sum(profit_rate[pr] * model.weekly_prod[pr]
                               for pr in products),
                      sense=maximize)

#Constraints
model.capacity = ConstraintList()
for pl in plants:
    model.capacity.add(
        sum(hours_per_batch.loc[pl,pr] * model.weekly_prod[pr]
            for pr in products) <= hours_available[pl])

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###

# note that we're using f-strings for output here which is a little different and cleaner than in the video
print(f"Maximum Profit = ${1000*model.profit():,.2f}")
for j in products:
    print(f"Batches of {j} = {model.weekly_prod[j]():.1f}")

Maximum Profit = $36,000.00
Batches of Doors = 2.0
Batches of Windows = 6.0

Pyomo Solution 4

In this solution we loop over the lists directly in the model construction. This approach is fine for quick and dirty solutions, but when possible the model and code should be written using meaningful index variables as in Solutions 2 and 3.

### PROBLEM DATA ###

# load data
num_products = 2
num_plants = 3
profit_rate = [3,5]
hours_available = [4,12,18]
hours_per_batch = [[1,0],[0,2],[3,2]]

### MODEL CONSTRUCTION ###
#Declaration
model = ConcreteModel()

#Decision Variables
model.weekly_prod = Var(range(num_products), domain=NonNegativeReals)

#Objective
model.profit = Objective(expr=sum(profit_rate[pr] * model.weekly_prod[pr]
                               for pr in range(num_products)),
                      sense=maximize)

# Constraints
model.capacity = ConstraintList()
for pl in range(num_plants):
    model.capacity.add(
        sum(hours_per_batch[pl][pr] * model.weekly_prod[pr]
            for pr in range(num_products)) <= hours_available[pl])
    
model.pprint()

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Maximum Profit = ${1000*model.profit():,.2f}")
for pr in range(num_products):
    print(f"Batches of {pr} = {model.weekly_prod[pr]():.0f}")

2 Set Declarations
    capacity_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {1, 2, 3}
    weekly_prod_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {0, 1}

1 Var Declarations
    weekly_prod : Size=2, Index=weekly_prod_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :  None : False :  True : NonNegativeReals
          1 :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    profit : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 3*weekly_prod[0] + 5*weekly_prod[1]

1 Constraint Declarations
    capacity : Size=3, Index=capacity_index, Active=True
        Key : Lower : Body                                : Upper : Active
          1 :  -Inf :   weekly_prod[0] + 0*weekly_prod[1] :   4.0 :   True
          2 :  -Inf : 0*weekly_prod[0] + 2*weekly_prod[1] :  12.0 :   True
          3 :  -Inf : 3*weekly_prod[0] + 2*weekly_prod[1] :  18.0 :   True

5 Declarations: weekly_prod_index weekly_prod profit capacity_index capacity
Maximum Profit = $36,000.00
Batches of 0 = 2
Batches of 1 = 6

The model code really hasn't changed much other than our loops are indexed by integers instead of lists of labels. This makes the model output a little harder to read and debug because we don't know which index values correspond to Windows or Doors:

# the objective function:
model.profit.pprint()

# the constraints:
model.capacity.pprint()

profit : Size=1, Index=None, Active=True
    Key  : Active : Sense    : Expression
    None :   True : maximize : 3*weekly_prod[0] + 5*weekly_prod[1]
capacity : Size=3, Index=capacity_index, Active=True
    Key : Lower : Body                                : Upper : Active
      1 :  -Inf :   weekly_prod[0] + 0*weekly_prod[1] :   4.0 :   True
      2 :  -Inf : 0*weekly_prod[0] + 2*weekly_prod[1] :  12.0 :   True
      3 :  -Inf : 3*weekly_prod[0] + 2*weekly_prod[1] :  18.0 :   True

For a small model such as this one it's still pretty easy to read the output and sort out the decision variables, but using meaningful labels can be really helpful for larger problems so we ask you not to use this approach in your homework.

Advantages and Disadvantages for the Four Approaches

The models are essentially the same in all four approaches, but the data preparation and labeling is different.

• Solution 1 - Typing the dictionaries directly

  • Advantage: it's easy to type out the dictionaries for small problems

  • Disadvantage: typing out the dictionaries for large problems isn't feasible

• Solution 2 - Building dictionaries from lists (Preferred)

  • Advantage: if the data is loaded in to lists (arrays or numpy arrays would also work well) then dictionaries of any size can be built; more levels of nesting can be used to accommodate variables with more than two indices (we'll meet an example in Lesson 3)

  • Disadvantage: The code is harder to write and requires some mastery of dictionaries.

• Solution 3 - Building Pandas data frames and series from lists or other array types (Preferred)

  • Advantage - accommodates large problems easily; data frames and series are easier to construct than dictionaries

  • Disadvantage - limited to problems in which the variables have one or two indices.

• Solution 4 - Using integers to loop directly over lists or arrays

  • Advantage - requires very little effort to prepare the data other than loading into lists or arrays

  • Disadvantages - code may be less readable and model output is more difficult to interpret

Goals and Guidelines for Abstract Modeling

• Separate the data from the model. You can type out lists or arrays with the raw data, but you should parse the coefficients into dictionaries or Pandas objects to be used in the model. In general you should pretend that the data is being loaded into lists, dataframes, or arrays and then processed before model construction.

• Suppose you want to write all of your constraints as "$\leq$" but you have some that go the wrong direction. For example $$3x - 2y \geq 5.$$ If you multiply both sides by $(-1)$ it reverses the inequality to yield $$-3x + 2y \leq -5.$$

• Getting all the inequality constraints in the same direction helps because the coefficients can be stored and used to construct the inequalities via a loop as we did with the Wyndor model above.

• Going forward, you'll lose points on your homework if you type out coefficients instead of using the abstract approach introduced here. It's more difficult, but it's essential for working with large problems. It's also great practice toward writing reusable code.

• If you want additional practice, many of the self-assessment problems below include abstract formulations.

Self Assessment: Investment Allocation

Textbook Problem 3.2-3. This is your lucky day. You have just won a $20,000 prize. You are setting aside $8,000 for taxes and partying expenses, but you have decided to invest the other $12,000. Upon hearing this news, two different friends have offered you an opportunity to become a partner in two different entrepreneurial ventures, one planned by each friend. In both cases, this investment would involve expending some of your time next summer as well as putting up cash. Becoming a full partner in the first friend’s venture would require an investment of $10,000 and 400 hours, and your estimated profit (ignoring the value of your time) would be $9,000. The corresponding figures for the second friend’s venture are $8,000 and 500 hours, with an estimated profit to you of $9,000. However, both friends are flexible and would allow you to come in at any fraction of a full partnership you would like. If you choose a fraction of a full partnership, all the above figures given for a full partnership (money investment, time investment, and your profit) would be multiplied by this same fraction. Because you were looking for an interesting summer job anyway (maximum of 600 hours), you have decided to participate in one or both friends’ ventures in whichever combination would maximize your total estimated profit. You now need to solve the problem of finding the best combination.

This self-assessment problem was also in Lesson 1 where you formulated and solved a linear model using the graphical method. Now copy and paste the code from the Wyndor problem above into a new cell and adapt it to solve this problem.

Web Mercantile - Abstract Formulation

We first saw this problem in Lesson 1 where it was solved by introducing individual variables for each possible month and lease duration combination. Such an approach was possible over a five month span, but imagine typing all of the variables and constraints for a 24 month, or longer, time period. In the solution below, using indexed variables, the code can be easily adapted to any time period.

The primary difference between this and the Wyndor example above is that the decision variables depend on both the month and the duration so they form an array. Study this example to see

Problem Description

This is problem 3.4-9, page 85, from the textbook.

Web Mercantile sells many household products through an online catalog. The company needs substantial warehouse space for storing its goods. Plans now are being made for leasing warehouse storage space over the next 5 months. Just how much space will be required in each of these months is known. However, since these space requirements are quite different, it may be most economical to lease only the amount needed each month on a month-by-month basis. On the other hand, the additional cost for leasing space for additional months is much less than for the first month, so it may be less expensive to lease the maximum amount required for the entire 5 months. Another option is the intermediate approach of changing the total amount of space leased (by adding a new lease and/or having an old lease expire) at least once but not every month.

The space requirement and the leasing costs for the various leasing periods are as follows:

Mathematical Formulation

Let $x_{m,d}$ represent the number of square feet to lease for $d$ months at the beginning of month $m$.

Let $r_d$ be the cost per square foot of leasing for a duration of $d$ months.

Let $s_m$ be the number of square feet required in month $m$.

Let $M$ be the set of months, e.g. $M = \{1,2,3,4,5\}.$

Let $D$ be the set of possible durations, e.g. $\{1,2,3,4,5\}.$

The total cost of leasing, to be minimized is:

$$Z = \sum_{m \in M} \sum_{d \in D} r_d x_{m,d}.$$

Note, that $x_{m,d}$ will be 0 for some combinations of $m$ and $d$. For instance, we'll have a constraint that $x_{5,2} = 0$ since in month 5 we cannot lease office space for 2 months. In fact, if $m+d > 6$ we must have $x_{m,d} = 0$. If we let $n_{mon}$ be the total number of months (e.g. $n_{mon} = 5$), then the constraint is $$x_{m,d} = 0 \mbox{ if } m + d > n_{mon}$$

For the final constraint we need to make sure we have adequate space leased in each month. We'll have to take into consideration that, for example, a lease made for 2 months at the beginning of month 1 will be available in the second month, but not in the third month. For each month we'll have add up all of the square feet that are leased in both the current month and previous months that are still available. Mathematically, we can add conditions to our sum like this:

for month $m$

$$\sum_{i \in M \mbox{ if } i \leq m} \sum_{\mbox{ }d \in D \mbox{ if } i + d > m} x_{i,d} \geq s_m.$$

That sum looks pretty complicated, but the first sum says to add up the leased square feet for all months up to and including the current month. The second sum says to add only the leased square feet for unexpired leases. If we treat the indices as integers we could also write it this way, for each month $m$:

$$\sum_{i = 1}^m \sum_{d = m - i + 1}^{m} x_{m,d} \geq s_m.$$

Still confused? Look at the data frame below where we have months as rows and durations as columns:

# details not important, dataframe for example
df = pd.DataFrame([[10, 11, 12, 13, 14],[15,16,17,18,0],[19,20,21,0,0],[22,23,0,0,0],[24,0,0,0,0]], 
                  index = [1,2,3,4,5], columns = [1,2,3,4,5] )
df

In month 1 the number of square feet leased is $10 + 11 + 12 + 13 + 14 = 60$.

In month 2 the number of square feet leased is $11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 = 116.$ Notice that the sum includes all the rentals from months 1 and 2 except for the 10 square feet rented at the beginning of month 1 for a duration of 1 month. That term is not included in the sum because $i + d = 1 + 1 = 2$ is not greater than 2 indicating that the 10 square feet is not available in month 2.

How many square feet are available in month 3? In month 4? In month 5? (answers are 150, 148, and 100 respectively)

Explanation of Pyomo Solution (video)

The video below explains some parts of the solution code. If you're content to study the code on you own, then you don't need to watch it.

# execute this cell for video
play_video("ds775_lesson2-webmercantile-abstract")

https://media.uwex.edu/content/ds/ds775_r19/ds775_lesson2-webmercantile-abstract/index.html

Pyomo Solution

# abstract Web Mercantile

### Problem Data ###
num_months = 5
max_duration = 5

months = range(1, num_months + 1) # generates 1, 2, ..., num_months
durations = range(1, max_duration + 1) # generates 1, 2, ..., max_duration

rent = dict(zip(durations, [65, 100, 135, 160, 190]))
space = dict(zip(months, [30000, 20000, 40000, 10000, 50000]))

### Pyomo Model ###

# Concrete Model
model = ConcreteModel(name="WebMerc2")

# Decision Variables
model.x_sqft = Var(months, durations, domain=NonNegativeReals)

# Objective (minimize is the default so it's not stated explicitly here)
model.obj = Objective(expr=sum(rent[d] * model.x_sqft[m, d] for m in months
                               for d in durations))

# Constraints
model.space_ct = ConstraintList()
for month in months:
    model.space_ct.add(
        sum(model.x_sqft[m, d] for m in months for d in durations
            if m <= month and m + d > month) >= space[month])

model.space_ct.pprint()

model.time_rule_ct = ConstraintList()
for m in months:
    for d in durations:
        if m + d > num_months + 1:
            model.time_rule_ct.add(model.x_sqft[m, d] == 0)

### Solution ###

solver = SolverFactory('glpk')
solver.solve(model)

### Display ###
# note that output is formatted with f-strings unlike the video

print(f"Total Cost = ${model.obj():,.2f}")

print("\nHere are the amounts to lease by month and duration:")
for m in months:
    for d in durations:
        if model.x_sqft[m, d]() > 0:
            print(f"Lease {model.x_sqft[m, d]():.0f} sq ft in month {m:d} for {d:d} months")

print("\nHere are the amounts needed and the total amount needed in each month:")
for m in months:
    amount_leased = sum(model.x_sqft[i, d]() for d in durations for i in months if i <= m and i + d > m)
    print(f"In month {m:d}, {space[m]:.0f} square feet are needed and {amount_leased:.0f} square feet are leased")

space_ct : Size=5, Index=space_ct_index, Active=True
    Key : Lower   : Body                                                                                                                                                                                                            : Upper : Active
      1 : 30000.0 :                                                                                                                                             x_sqft[1,1] + x_sqft[1,2] + x_sqft[1,3] + x_sqft[1,4] + x_sqft[1,5] :  +Inf :   True
      2 : 20000.0 :                                                                                     x_sqft[1,2] + x_sqft[1,3] + x_sqft[1,4] + x_sqft[1,5] + x_sqft[2,1] + x_sqft[2,2] + x_sqft[2,3] + x_sqft[2,4] + x_sqft[2,5] :  +Inf :   True
      3 : 40000.0 :                                           x_sqft[1,3] + x_sqft[1,4] + x_sqft[1,5] + x_sqft[2,2] + x_sqft[2,3] + x_sqft[2,4] + x_sqft[2,5] + x_sqft[3,1] + x_sqft[3,2] + x_sqft[3,3] + x_sqft[3,4] + x_sqft[3,5] :  +Inf :   True
      4 : 10000.0 :               x_sqft[1,4] + x_sqft[1,5] + x_sqft[2,3] + x_sqft[2,4] + x_sqft[2,5] + x_sqft[3,2] + x_sqft[3,3] + x_sqft[3,4] + x_sqft[3,5] + x_sqft[4,1] + x_sqft[4,2] + x_sqft[4,3] + x_sqft[4,4] + x_sqft[4,5] :  +Inf :   True
      5 : 50000.0 : x_sqft[1,5] + x_sqft[2,4] + x_sqft[2,5] + x_sqft[3,3] + x_sqft[3,4] + x_sqft[3,5] + x_sqft[4,2] + x_sqft[4,3] + x_sqft[4,4] + x_sqft[4,5] + x_sqft[5,1] + x_sqft[5,2] + x_sqft[5,3] + x_sqft[5,4] + x_sqft[5,5] :  +Inf :   True
Total Cost = $7,650,000.00

Here are the amounts to lease by month and duration:
Lease 30000 sq ft in month 1 for 5 months
Lease 10000 sq ft in month 3 for 1 months
Lease 20000 sq ft in month 5 for 1 months

Here are the amounts needed and the total amount needed in each month:
In month 1, 30000 square feet are needed and 30000 square feet are leased
In month 2, 20000 square feet are needed and 30000 square feet are leased
In month 3, 40000 square feet are needed and 40000 square feet are leased
In month 4, 10000 square feet are needed and 30000 square feet are leased
In month 5, 50000 square feet are needed and 50000 square feet are leased

Self Assessment: A Holiday Factory

A company is planning to design and manufacture children's toys over an 11-month span prior to a major holiday. In the beginning, the company won't need much space to design and plan the toy, but the square footage required will grow as they stockpile parts for the toys and start manufacturing. Toward the end of the 11-month span the square feet required will start to diminish as toys are shipped to stores and distribution centers. The square feet needed in each month is as follows:

Month	Square Feet	Month	Square Feet
1	2000	7	10000
2	2000	8	10000
3	3000	9	9000
4	4000	10	7000
5	6000	11	5000
6	10000

The rent per square foot starts at $20 per month, but decreases for leases of longer duration as follows.

Duration (months)	Rent ($ per square foot)
1	20
2	20 + 19 = 39
3	20 + 19 + 18 = 57
$\vdots$	$\vdots$
11	165

Learn how easy it is to reuse models by figuring out the pattern for durations and then adapting the example code above to find the minimum cost for leasing the required square footage. Try to find a way to compute the rent for each duration using a sum or a formula instead of just computing each by "hand".

Sausages Blending - Abstract Formulation

This example, first presented in Lesson 1, is adapted from here: http://benalexkeen.com/linear-programming-with-python-and-pulp-part-4/.

In Lesson 1 this problem was formulated and solved by introducing individual variables for each ingredient and sausage type. Here, we'll use an array of variables and configure the objective function and constraints in a way that can be more easily generalized. The second will be to write the model in a way that can be easily extended.

Problem Description

We're going to make sausages by blending pork, wheat, and starch. Our objective is to minimize the cost of making the sausages. The table below shows the ingredients available, the cost, and the amount of each ingredient available from our supplier:

Ingredient	Cost ($/kg)	Amount (kg)
Pork	4.32	7
Wheat	2.46	20.0
Starch	1.86	17

Additionally, we have 23 kg of pork on hand that we must use in the sausages.

We want to make 2 types of sausage:

• Economy ( > 40% pork)

• Premium ( > 60% pork)

Each sausage is 50 grams (0.05 kg).

According to government regulations, the most starch we can use in our sausages is 25% by weight.

We have a demand for 350 economy sausages and 500 premium sausages.

Mathematical formulation

We'll write things generically so that the model can be used for any numbers of ingredients and types.

Let $I$ represent the set of ingredients and let $T$ represent the set of types of sausage.

For each $i \in I$ and $t \in T$

Quantity	Description
$x_{i,t}$	kg of ingredient $i$ to use in sausage of type $t$ (decision variable)
$c_i$	cost per kg of ingredient $i$ in dollars
$pmin_{i,t}$	minimum proportion of ingredient $i$ in sausage of type $t$, use 0 for no minimum
$pmax_{i,t}$	maximum proportion of ingredient $i$ in sausage of type $t$, use 1 for no maximum
$imax_{i}$	maximum kilograms of ingredient $i$ that are available
$imin_{i}$	minimum kilograms of ingredient $i$ that must be used
$d_{t}$	total kg of sausage of type $t$ that are demanded

To write the cost function we need to add up the quantity of each ingredient used across all the sausage types, so for each ingredient $i$ this is $$\sum_{t \in T} x_{i,t}.$$ Now the cost of ingredient $i$ will be $$c_i \left( \sum_{t \in T} x_{i,t} \right).$$ Finally, the objective function is the total cost for all ingredients and types: $$\mbox{Cost} = \sum_{i \in I} c_i \left( \sum_{t \in T} x_{i,t} \right).$$

Now let's turn to the constraints. First we need the total kg of each type of sausage made to equal the demand for that type.

$$\sum_{i \in I} x_{i,t} = d_t \mbox{ for each } t \in T.$$

The sum on the left represents the total kg of ingredients used in sausage $j$.

Now we need the minimum proportion constraints: $$x_{i,t} \geq pmin_{i,t} \sum_{k \in I} x_{k,t}.$$

That sum on the right now uses the index $k$ which serves as a dummy variable to add up the total kg of each type of sausage.

The maximum proportion constraints are similar: $$x_{i,t} \leq pmax_{i,t} \sum_{k \in I} x_{k,t}.$$

The availability constraints enforce the maximum amount of ingredient $i$ that is available: $$\sum_{t \in T} x_{i,t} \leq imax_i.$$ The sum on the left here is the total kg of ingredient $i$ used over all types of sausages.

The "must use" constraints are similar: $$\sum_{t \in T} x_{i,t} \geq imin_i.$$

Here are the concrete and abstract model comparisons side-by-side for comparison:

Description	Concrete Model	Abstract Model
Minimize Cost =	$4.32 ( p_e + p_p) + 2.46( w_e + w_p) + 1.86 (s_e + s_p)$	$\displaystyle \sum_{i \in I} c_i \left(\sum_{t \in T} x_{i,t} \right)$
Demand Constraints	$p_e + w_e + s_e = 350 \times 0.05$ $p_p + w_p + s_p = 500 \times 0.05$	$\displaystyle \sum_{i \in I} x_{i,t} = d_t \mbox{ for each } t \in T$
Minimimum Proportion Constraints	$p_e \geq 0.4 (p_e + w_e + s_e)$ $p_p \geq 0.6 (p_p + w_p + s_p)$	$\displaystyle x_{i,t} \geq pmin_{i,t} \sum_{k \in I} x_{k,t}$
Maximum Proportion Constraints	$s_e \leq 0.25 (p_e + w_e + s_e)$ $s_p \leq 0.25 (p_p + w_p + s_p)$	$\displaystyle x_{i,t} \leq pmax_{i,t} \sum_{k \in I} x_{k,t}$
Maximum Ingredient Constraints	$p_e + p_p \leq 30$ $w_e + w_p \leq 20$ $s_e + s_p \leq 17$	$\displaystyle \sum_{t \in T} x_{i,t} \leq imax_i$
Minimum Ingredient Constraints	$p_e + p_p \geq 23$	$\displaystyle \sum_{t \in T} x_{i,t} \geq imin_i$

Explanation of Pyomo Solution (video)

The video below explains some parts of the solution code. If you're content to study the code on you own, then you don't need to watch it.

# execute this cell for video
play_video("ds775_lesson2-blending-abstract")

Pyomo Solution

# abstract Sausage Factory

### Problem Data ###

types = ['economy','premium']
ingredients = ['pork', 'wheat', 'starch']

cost = dict( zip( ingredients, [4.32, 2.46, 1.86] ) )

kg_per_sausage = 0.05
number_each_type = dict( zip( types, [350, 500] ) )

mnpi = [[.4,.6],[0,0],[0,0]] # min proportions
min_prop_ing = { ingredients[i]:{ types[j]:mnpi[i][j] for j in range(len(types)) } for i in range(len(ingredients)) }
mxpi = [[1,1],[1,1],[.25,.25]] # max proportions
max_prop_ing = { ingredients[i]:{ types[j]:mxpi[i][j] for j in range(len(types)) } for i in range(len(ingredients)) }

max_ingredient = dict( zip( ingredients, [30, 20, 17] ) )
min_ingredient = dict( zip( ingredients, [23,  0,  0] ) )

### Pyomo Model ###

# Concrete Model
M = ConcreteModel(name = "Sausages")

# Decision Variables
M.amount = Var(ingredients, types, domain = NonNegativeReals)

# Objective
M.cost = Objective( expr = sum( cost[i] * sum(M.amount[i,t] for t in types) 
                               for i in ingredients), sense = minimize )

M.tot_sausages_ct = ConstraintList()
for t in types:
    M.tot_sausages_ct.add( sum( M.amount[i,t] for i in ingredients ) 
                         == kg_per_sausage * number_each_type[t] )

M.min_prop_ct = ConstraintList()
for i in ingredients:
    for t in types:
        M.min_prop_ct.add( M.amount[i,t] >= min_prop_ing[i][t] *
                         sum( M.amount[k,t] for k in ingredients ) )

M.max_prop_ct = ConstraintList()
for i in ingredients:
    for t in types:
        M.max_prop_ct.add( M.amount[i,t] <= max_prop_ing[i][t] * 
                         sum( M.amount[k, t] for k in ingredients ) )
        
M.max_ingredient_ct = ConstraintList()
for i in ingredients:
    M.max_ingredient_ct.add( sum( M.amount[ i, t] for t in types ) <= 
                           max_ingredient[i] )
    
M.min_ingredient_ct = ConstraintList()
for i in ingredients:
    M.min_ingredient_ct.add( sum( M.amount[ i, t] for t in types ) >=
                           min_ingredient[i] )

### Solution ###
solver = SolverFactory('glpk')
solver.solve(M)

### Output ###
print(f"Total Cost = ${M.cost():,.2f}")

# put amounts in dataframe for nicer display
import pandas as pd
dvars = pd.DataFrame( [ [M.amount[i,t]() for t in types] for i in ingredients ],
                    index = ['Pork','Wheat','Starch'],
                    columns = ['Economy','Premium'])
print("Kilograms of each ingredient in each type of sausage:")
dvars

Total Cost = $140.96
Kilograms of each ingredient in each type of sausage:

Self Assessment: Supply and Demand Problem

Use what you've learned in the examples above to formulate and solve a linear programming model for textbook problem 3.4-11.

The Medequip Company produces precision medical diagnostic equipment at two factories. Three medical centers have placed orders for this month’s production output. The table below shows what the cost would be for shipping each unit from each factory to each of these customers. Also shown are the number of units that will be produced at each factory and the number of units ordered by each customer.

	Customer 1	Customer 2	Customer 3	Output
Factory 1	\$600	\$800	\$700	400 units
Factory 2	\$400	\$900	\$600	500 units
Order size	300 units	200 units	400 units

A decision now needs to be made about the shipping plan for how many units to ship from each factory to each customer.

(a) Formulate a linear programming model for this problem.

(b) Solve this model using Pyomo.

Sensitivity Analysis

What is Sensitivity Analysis?

We want to know what happens when the model coefficients change. What effect do changes have on the optimal solution? We say a bit more about the need for Sensitivity Analysis in the next video:

# execute cell for video
play_video("ds775_lesson2_what-is-sensitivity-analysis")

Common Types of Sensitivity Analysis for Linear Programs

Below are four videos. Each video will motivate one common sensitivity analysis question, using Desmos, and will also show you where to find that information in a sensitivity report. We'll show you how to produce a sensitivity report using Pyomo and the GLPK solver further below. Here is the Desmos we used in the videos.

Video 1 - Active Constraints

• An active or binding constraint is one for which equality holds at the optimal values of the decision variables. Changes to the coefficients in an active constraint can change the optimal solution.

# execute cell for video
play_video("ds775_lesson2_active-vs-inactive-constraints")

Video 2 - Shadow Prices

• A shadow price is the rate at which $Z$ changes for each unit increase in the amount of resource $b_i$ (when $b_i$ remains in its allowable range and none of the other parameters are allowed to change). We used Desmos to produce the graphs in the video. Here is the Desmos we used in the videos.

# execute cell for video
play_video("ds775_lesson2_shadow-prices")

Video 3 - Allowable Range of Constraint Right Side

• The allowable range of constraint right side, $b_i$ is the range values over which the corresponding constraint remains active when no other parameters are allowed to change. Outside this range the constraint is no longer active and is not involved in the optimal solution. Note, as $b_i$ changes within its allowable range the optimal objective function value and decision variables values will change, but will be determined by the same set of active constraints. Here is the Desmos we used in the videos.

# new video for allowable range of constraint right side
play_video("ds775_lesson2_allowable-range-constraint-bounds")

Video 4 - Allowable Range of Objective Function Coefficients

• Over what range of values can we change a single objective coefficient without changing the location of the optimal corner point feasible solution? Note that the optimal objective function value will change but not the optimal decision variable values. Here is the Desmos we used in the videos.

# new video for obj function coef allowable range
play_video("ds775_lesson2_allowable-range-obj-function-coefficients")

Self Assessment: Positive Shadow Price

True or False: A positive shadow price indicates that the right-hand side of that constraint is a sensitive parameter.

Self Assessment: Allowable Range (Objective Coef)

True or False: For the objective coefficients, the allowable range is the range of values over which the current optimal solution remains optimal, assuming no change in the other coefficients.

Self Assessment: Changing Parameters

True or False: Any change in any parameter will necessarily change the optimal solution.

Self Assessment: Graphical Exploration of Sensitivity

Use DESMOS to answer the following questions about this model formulation:

Maximize $Z = c_1 x + c_2 y = 2 x + 5 y$

subject to

$x + 2y \leq 10$ (resource 1, $b_1$)

$x + 3y \leq 12$ (resource 2, $b_2$)

$x_1, x_2 \geq 0$

(a) Set up this model in DESMOS and show the optimal value of $Z$ and the coordinates of $x$ and $y$ where the optimal solution occurs. Provide a screenshot of the graph you used to find this answer with the optimal $Z$ and coordinates of $x$ and $y$ displayed.

(b) With a slider for $Z$ (and $b_2$ if you like), find the shadow price for resource 2 ($b_2$). Provide a screenshot of the graph with the new optimal $Z$ and coordinates of $x$ and $y$ at this point displayed.

(c) With a slider for $b_2$, find the allowable range for resource 2 ($b_2$). Provide screenshots of the graphs you used to find this answer, one for the lower bound and one for the upper bound of $b_2$.

(d) With sliders for $Z$ and $c_2$, find the allowable range for the unit profit of activity 2 ($c_2$). Provide screenshots of the graphs you used to find this answer, one for the lower bound and one for the upper bound of $c_2$.

Generating a Sensitivity Report in Python

Pyomo doesn't automatically generate a sensitivity analysis, but it's possible to write the model to a standardized file (called an LP file) and then run a solver to produce a sensitivity report. For directions about interpreting the sensitivity report refer to the last few videos or look at the annotated image provided further below.

First, run the following cell to load up the abstract Wyndor model from above.

### PROBLEM DATA ###

# load data
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
profit_rate_list = [3,5]
hours_avail_list = [4,12,18]
hours_per_batch_list = [ [1,0], [0,2], [3,2] ]

# parse lists into dictionaries
profit_rate = dict( zip( products, profit_rate_list) )
hours_available = dict( zip( plants, hours_avail_list) )
hours_per_batch = { pl: { pr: hours_per_batch_list[i][j] for j,pr in enumerate(products)} for i,pl in enumerate(plants)}

### MODEL CONSTRUCTION ###
#Declaration
model = ConcreteModel()

#Decision Variables
model.weekly_prod = Var(products, domain=NonNegativeReals)

#Objective
model.profit = Objective(expr=sum(profit_rate[pr] * model.weekly_prod[pr]
                               for pr in products),
                      sense=maximize)

#Constraints
model.capacity = ConstraintList()
for pl in plants:
    model.capacity.add(
        sum(hours_per_batch[pl][pr] * model.weekly_prod[pr]
            for pr in products) <= hours_available[pl])

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###

# note that we're using f-strings for output here which is a little different and cleaner than in the video
print(f"Maximum Profit = ${1000*model.profit():,.2f}")
for j in products:
    print(f"Batches of {j} = {model.weekly_prod[j]():.1f}")
    
model.capacity.display()

Maximum Profit = $36,000.00
Batches of Doors = 2.0
Batches of Windows = 6.0
capacity : Size=3
    Key : Lower : Body : Upper
      1 :  None :  2.0 :   4.0
      2 :  None : 12.0 :  12.0
      3 :  None : 18.0 :  18.0

Now run the cell below to write the model to an LP file to create a sensitivity report. Note: the solver produces some output that we do not need. You can just ignore the output here.

model.write('model.lp', io_options={'symbolic_solver_labels': True})

('model.lp', 140368987244864)

!glpsol -m model.lp --lp --ranges sensit.sen

GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 -m model.lp --lp --ranges sensit.sen
Reading problem data from 'model.lp'...
3 rows, 2 columns, 4 non-zeros
26 lines were read
GLPK Simplex Optimizer 5.0
3 rows, 2 columns, 4 non-zeros
Preprocessing...
1 row, 2 columns, 2 non-zeros
Scaling...
 A: min|aij| =  2.000e+00  max|aij| =  3.000e+00  ratio =  1.500e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 1
*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (2)
*     2: obj =   3.600000000e+01 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.0 Mb (40400 bytes)
Write sensitivity analysis report to 'sensit.sen'...

Run the following cell to print the report. Alternately, you can open the report separately using anything that will display plain text. We've found that using "Print" in CoCalc to print to a pdf file in landscape produces a readable document.

# print sensitivity report
np.set_printoptions(linewidth=110)
f = open('sensit.sen', 'r')
file_contents = f.read()
print(file_contents)
f.close()

GLPK 5.0  - SENSITIVITY ANALYSIS REPORT                                                                         Page   1

Problem:    
Objective:  profit = 36 (MAXimum)

   No. Row name     St      Activity         Slack   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 c_u_capacity(1)_
                    BS       2.00000       2.00000          -Inf         .           -3.00000      30.00000 c_u_capacity(3)_
                                            .            4.00000        6.00000       4.50000      45.00000 c_u_capacity(2)_

     2 c_u_capacity(2)_
                    NU      12.00000        .               -Inf        6.00000      -1.50000      27.00000 c_u_capacity(1)_
                                           1.50000      12.00000       18.00000          +Inf      45.00000 weekly_prod(Doors)

     3 c_u_capacity(3)_
                    NU      18.00000        .               -Inf       12.00000      -1.00000      30.00000 weekly_prod(Doors)
                                           1.00000      18.00000       24.00000          +Inf      42.00000 c_u_capacity(1)_

GLPK 5.0  - SENSITIVITY ANALYSIS REPORT                                                                         Page   2

Problem:    
Objective:  profit = 36 (MAXimum)

   No. Column name  St      Activity      Obj coef   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 weekly_prod(Doors)
                    BS       2.00000       3.00000        .                -Inf        .           30.00000 c_u_capacity(3)_
                                            .               +Inf        4.00000       7.50000      45.00000 c_u_capacity(2)_

     2 weekly_prod(Windows)
                    BS       6.00000       5.00000        .             3.00000       2.00000      18.00000 c_u_capacity(2)_
                                            .               +Inf        6.00000          +Inf          +Inf

End of report

Self-Assessment: Solve and Perform Sensitivity

Textbook Problem 4.7-6 (c, b)

Consider the following problem:

Maximize $Z = 5 x_1 + 4x_2 - x_3 + 3 x_4$

subject to

$3x_1 + 2x_2 - 3x_3 + x_4 \leq 24 \mbox{ (resource 1) }$

$3x_1 + 3x_2 + x_3 + 3x_4 \leq 36 \mbox{ (resource 2) }$

and $x_1 \geq 0, x_2 \geq 0, x_3 \geq 0, x_4 \geq 0.$

Use Pyomo to solve the problem and then generate sensitivity information. Use this information to identify the shadow price for each resource, the allowable range for each objective function coefficient, and the allowable range for each right-hand side. You should practice writing abstract code with the model and data separated.

Self-Assessment: Formulate, Solve, and Perform Sensitivity #1

Textbook Problem 7.3-4 (a & f with Pyomo)

One of the products of the G.A. Tanner Company is a special kind of toy that provides an estimated unit profit of $3. Because of a large demand for this toy, management would like to increase its production rate from the current level of 1,000 per day. However, a limited supply of two subassemblies (A and B) from vendors make this difficult. Each toy requires two subassemblies of type A, but the vendor providing these subassemblies would only be able to increase its supply rate from the current 2,000 per day to a maximum of 3,000 per day. Each toy requires only one subassembly of type B, but the vendor providing these subassemblies would be unable to increase its supply rate above the current level of 1,000 per day. Because no other vendors currently are available to provide these subassemblies, management is considering initiating a new production process internally that would simultaneously produce an equal number of subassemblies of the two types to supplement the supply from the two vendors. It is estimated that the company’s cost for producing one subassembly of each type would be $2.50 more than the cost of purchasing these subassemblies from the two vendors. Management wants to determine both the production rate of the toy and the production rate of each pair of subassemblies (one A and one B) that would maximize the total profit. The following table summarizes the data for the problem.

(a) Formulate the mathematical model for this problem and solve it using Pyomo in Python.

(f) Generate a sensitivity report to find the allowable range for the unit profit of each activity (toys and subassemblies).

Self-Assessment: Formulate, Solve, and Perform Sensitivity #2

Textbook Problem 7.3-5 (a, b, f)

Reconsider Problem 7.3-4. After further negotiations with each vendor, management of the G.A. Tanner Co. has learned that either of them would be willing to consider increasing their supply of their respective subassemblies over the previously stated maxima (3,000 subassemblies of type A per day and 1,000 of type B per day) if the company would pay a small premium over the regular price for the extra subassemblies. The size of the premium for each type of subassembly remains to be negotiated. The demand for the toy being produced is sufficiently high so that 2,500 per day could be sold if the supply of subassemblies could be increased enough to support this production rate. Assume that the original estimates of unit profits given in 7.3-4 are accurate.

(a) Formulate the mathematical model for this problem and with the original maximum supply levels and the additional constraint that no more than 2,500 toys should be produced per day. Solve it using Pyomo in Python.

(b) Without considering the premium, use Pyomo to determine the shadow price for the subassembly A constraint by solving the model again after increasing the maximum supply by 1. Use this shadow price to determine the maximum premium that the company should be willing to pay for each subassembly of this type.

(f) Use the sensitivity report to determine the shadow price for each of the subassembly constraints and the allowable range for the right-hand side of each of these constraints.

Constructing Dictionaries from Lists

Here are some examples to help with learning how to construct dictionaries from lists. We don't provide explanations for the code, but study and use any of the constructions as needed.

A single dictionary

First, let's look at creating a dictionary from two separate lists of keys and values. We generally prefer "method 7" for its compactness and simplicity. Study only "method 7" if there are too many options to absorb.

plants = ['Plant1', 'Plant2', 'Plant3']
hours_avail = [4, 12, 18]

# method 0 - type it out (This method should be avoided where possible.)
hours_avail_dict = { 'Plant1':4, 'Plant2':12, 'Plant3':18}
print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

# method 1 - use a for loop
hours_avail_dict = {}
for i in range( len(plants) ):
    hours_avail_dict[ plants[i] ] = hours_avail[i]

print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

# method 2 - use a for loop, v2 with zip
hours_avail_dict = {}
for p,h in zip(plants,hours_avail):
    hours_avail_dict[ p ] = h

print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

We like to use enumerate because it makes it easy to index both by position in the list and by the list elements. If you'd like to review the enumerate command this is a pretty good tutorial.

# method 3 - use a for loop, v3 with enumerate
hours_avail_dict = {}
for i, pl in enumerate(plants):
    hours_avail_dict[ pl ] = hours_avail[i]

print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

# method 4 - use a comprehension (compare to method 1)
hours_avail_dict = { plants[i]:hours_avail[i] for i in range( len(plants) )}
print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

# method 5 - use a comprehension with zip (compare to method 2)
hours_avail_dict = { p:h for p,h in zip(plants,hours_avail)}
print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

# method 6 - use a comprehension with enumerate (compare to method 3)
hours_avail_dict = { pl:hours_avail[i] for i,pl in enumerate(plants)}
print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

# method 7 - use the dictionary constructor along with zip (most 'Pythonic' approach)
hours_avail_dict = dict( zip( plants, hours_avail) )
print(hours_avail_dict)

{'Plant1': 4, 'Plant2': 12, 'Plant3': 18}

Nested Lists to Nested Dictionaries

For converting nested lists to nested dictionaries we generally prefer Methods 13 and 14 below. We use Method 13 in most of our solutions. Method 15 is the most 'Pythonic' but is not as easy to follow as some other methods.

# method 10 - nested dictionary from nested for loop
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
hours_per_batch = [ [1, 0], [0, 2], [3, 2] ]

hours_per_batch_dict = {}
for i in range( len(plants) ):
    hours_per_batch_dict[plants[i]] = {}
    for j in range( len(products) ):
        hours_per_batch_dict[plants[i]][products[j]] = hours_per_batch[i][j]
        
print(hours_per_batch_dict)

{'Plant1': {'Doors': 1, 'Windows': 0}, 'Plant2': {'Doors': 0, 'Windows': 2}, 'Plant3': {'Doors': 3, 'Windows': 2}}

# method 11 - nested dictionary from nested for loop with enumerate
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
hours_per_batch = [ [1, 0], [0, 2], [3, 2] ]

hours_per_batch_dict = {}
for i,pl in enumerate(plants):
    hours_per_batch_dict[ pl ]= {}
    for j,pr in enumerate(products):
        hours_per_batch_dict[ pl ][ pr ] = hours_per_batch[ i ][ j ]
        
print(hours_per_batch_dict)

{'Plant1': {'Doors': 1, 'Windows': 0}, 'Plant2': {'Doors': 0, 'Windows': 2}, 'Plant3': {'Doors': 3, 'Windows': 2}}

# method 12 - comprehension without enumerate
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
hours_per_batch = [ [1, 0], [0, 2], [3, 2] ]

hours_per_batch_dict = { plants[i]: { products[j]: hours_per_batch[i][j] for j in range( len(products) )} for i in range( len(plants) )}
        
print(hours_per_batch_dict)

{'Plant1': {'Doors': 1, 'Windows': 0}, 'Plant2': {'Doors': 0, 'Windows': 2}, 'Plant3': {'Doors': 3, 'Windows': 2}}

# method 13 - comprehension with enumerate
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
hours_per_batch = [ [1, 0], [0, 2], [3, 2] ]

hours_per_batch_dict = { pl: { pr: hours_per_batch[i][j] for j,pr in enumerate(products)} for i,pl in enumerate(plants)}
        
print(hours_per_batch_dict)

{'Plant1': {'Doors': 1, 'Windows': 0}, 'Plant2': {'Doors': 0, 'Windows': 2}, 'Plant3': {'Doors': 3, 'Windows': 2}}

# method 14 - comprehension with enumerate, replace inner with dict+zip construction
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
hours_per_batch = [ [1, 0], [0, 2], [3, 2] ]

hours_per_batch_dict = { pl: dict( zip( products, hours_per_batch[i]) ) for i,pl in enumerate(plants)}
        
print(hours_per_batch_dict)

{'Plant1': {'Doors': 1, 'Windows': 0}, 'Plant2': {'Doors': 0, 'Windows': 2}, 'Plant3': {'Doors': 3, 'Windows': 2}}

# method 15 - using dict+zip without enumerate
products = ['Doors', 'Windows']
plants = ['Plant1', 'Plant2', 'Plant3']
hours_per_batch = [ [1, 0], [0, 2], [3, 2] ]

hours_per_batch_dict = dict( zip( plants, [dict( zip( products, row)) for row in hours_per_batch] ))
        
print(hours_per_batch_dict)

{'Plant1': {'Doors': 1, 'Windows': 0}, 'Plant2': {'Doors': 0, 'Windows': 2}, 'Plant3': {'Doors': 3, 'Windows': 2}}