HW 2¶

Please formulate the baa99 problem as an All-in-One LP model in Pyomo (see the second page in Lecture 3.2), and find the decisions. The first stage decisions are bounded within (0, 217).

The demand vector (d1, d2) has two scenarios, scenario1: (161.326406, 156.3195565); scenario 2: (52.29173734, 48.17833053).

By: Chengyi (Jeff) Chen

%load_ext autotime
%load_ext nb_black
%matplotlib inline

from IPython.display import Latex
import numpy as np
import cvxpy as cp
import pandas as pd
import scipy.stats
from scipy import optimize
import matplotlib.pyplot as plt
from matplotlib import rc

plt.rcParams["figure.dpi"] = 300
plt.rcParams["figure.figsize"] = (16, 12)

time: 970 ms (started: 2021-02-16 14:20:15 +08:00)

BAA99 Problem¶

Problem Setup¶

We are an electronics company that sells computer type 1 and type 2 for \$8 and \\$4 respectively. To build the computers, we manufacture two types of chips, quantity $x_1$ and $x_2$, which are \$4 and \\$2 respectively. Chip 1 can be used to make both Type 1 and Type 2 computer, but chip 2 can only be used to make Type 2 computer. Hence, the quantity of chip 1 we allocate to making type 1 computer is denoted as $w_{11}$, quantity of chip 1 to type 2 is $w_{12}$, and quantity of chip 2 to type 2 is $w_{22}$. Demand follows a Bernoulli distribution with $p=0.5$ such that high demand for type 1 computers in scenario 1 is $d^{s = 1}_{1} = 161.326406$ and demand for type 2 computers in scenario 1 is $d^{s = 1}_{2} = 156.3195565$ while low demand for type 1 computers in scenario 2 is $d^{s = 2}_{2} = 52.29173734$ and demand for type 2 computers in scenario 2 is $d^{s = 2}_{2} = 48.17833053$. Also, the backorder cost for insufficient supply of computers is the same for both computers at \$10 while holding cost for excess supply of computers produced is \\$0.20. We denote $u^s_i, v^s_i$ as the quantity of computers that are in insufficient or excess, $s = 1, 2$ and $i = 1, 2$. Notice that for each scenario $s$ if $u^s_i > 0 \implies v^s_i = 0$ and $v^s_i > 0 \implies u^s_i = 0$. Our optimization problem is then to formulate a way to maximize the net profit accounting for all scenarios:

\begin{aligned} \underset{x_1, x_2, w^s_{11}, w^{s}{12}, w^{s}{22}, u^s_i, v^s_i}{\text{Maximize }} &-4x_1 - 2x_2 + \frac{1}{n}\sum^{n=2}{s=1} \left(8 \left(w^s{11} - v^s_1\right) + 4\left(w^{s}{12} + w^{s}{22} - v^s_2\right) - \left( 0.2v^s_1 + 0.2v^s_2 \right) - \left( 10u^s_1 + 10u^s_2 \right) \right) \ \text{subject to } &x_1 - w^s_{11} - w^s_{12} = 0 \ &x_2 - w^s_{22} = 0 \ &w^s_{11} + u^s_1 - v^s_1 = d^s_1 \ &w^s_{12} + w^s_{22} + u^s_2 - v^s_2 = d^s_2 \ &s = 1, 2 \text{ and all variables }\geq 0, \in \mathcal{N}. \ \end{aligned}

# Create static optimization variables
x = cp.Variable((2,), integer=False)

# Create scenario dependent variables
num_scenarios = 2
scenarios = {}
for s in range(num_scenarios):
    scenarios[s] = {}
    scenarios[s]["w"] = cp.Variable((3,), integer=False)  # [w^s_11, w^s_12, w^s_22]
    scenarios[s]["u"] = cp.Variable((2,), integer=False)  # [u^s_1, u^s_2]
    scenarios[s]["v"] = cp.Variable((2,), integer=False)  # [v^s_1, v^s_2]
    scenarios[s]["d"] = (
        [161.326406, 156.3195565] if s == 0 else [52.29173734, 48.17833053]
    )  # [d^s_1, d^s_2]

# Create constraints.
constraints = [
    x >= 0,
    x <= 217,
]
for s in range(num_scenarios):
    constraints += [
        x[0] - scenarios[s]["w"][0] - scenarios[s]["w"][1] >= 0,
        x[0] - scenarios[s]["w"][0] - scenarios[s]["w"][1] <= 0,
        x[1] - scenarios[s]["w"][2] >= 0,
        x[1] - scenarios[s]["w"][2] <= 0,
        scenarios[s]["w"][0] + scenarios[s]["u"][0] - scenarios[s]["v"][0]
        <= scenarios[s]["d"][0],
        scenarios[s]["w"][0] + scenarios[s]["u"][0] - scenarios[s]["v"][0]
        >= scenarios[s]["d"][0],
        scenarios[s]["w"][1]
        + scenarios[s]["w"][2]
        + scenarios[s]["u"][1]
        - scenarios[s]["v"][1]
        <= scenarios[s]["d"][1],
        scenarios[s]["w"][1]
        + scenarios[s]["w"][2]
        + scenarios[s]["u"][1]
        - scenarios[s]["v"][1]
        >= scenarios[s]["d"][1],
        scenarios[s]["w"] >= 0,
        scenarios[s]["u"] >= 0,
        scenarios[s]["v"] >= 0,
    ]

# Form objective.
obj = cp.Maximize(
    np.array([-4, -2]) @ x
    + np.mean(
        [
            8 * (scenarios[s]["w"][0] - scenarios[s]["v"][0])
            + 4 * (scenarios[s]["w"][1] + scenarios[s]["w"][2] - scenarios[s]["v"][1])
            - np.array([0.2, 0.2]) @ scenarios[s]["v"]
            - np.array([10, 10]) @ scenarios[s]["u"]
            for s in range(num_scenarios)
        ]
    )
)

# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()

print("Linear Programming Solution")
print("=" * 40)
print(f"Status: {prob.status}")
print(f"The optimal value is: {prob.value}")
print(f"The optimal solution is: x = {x.value}")

Linear Programming Solution
========================================
Status: optimal
The optimal value is: 283.8060209104777
The optimal solution is: x = [161.326406   156.31955651]
time: 53.5 ms (started: 2021-02-16 14:20:16 +08:00)

Scenario 1: High Demand¶

pd.Series(
    {
        r"$w^1_{11}$": str(round(scenarios[0]["w"][0].value, 2)),
        r"$w^1_{12}$": str(round(scenarios[0]["w"][1].value, 2)),
        r"$w^1_{22}$": str(round(scenarios[0]["w"][2].value, 2)),
        r"$u^1_{1}$": str(round(scenarios[0]["u"][0].value, 2)),
        r"$u^1_{2}$": str(round(scenarios[0]["u"][0].value, 2)),
        r"$v^1_{1}$": str(round(scenarios[0]["v"][0].value, 2)),
        r"$v^1_{2}$": str(round(scenarios[0]["v"][0].value, 2)),
    },
    name="value",
).to_frame()

	value
$w^1_{11}$	161.33
$w^1_{12}$	0.0
$w^1_{22}$	156.32
$u^1_{1}$	0.0
$u^1_{2}$	0.0
$v^1_{1}$	-0.0
$v^1_{2}$	-0.0

time: 9.09 ms (started: 2021-02-16 14:20:16 +08:00)

Scenario 2: Low Demand¶

pd.Series(
    {
        r"$w^2_{11}$": str(round(scenarios[1]["w"][0].value, 2)),
        r"$w^2_{12}$": str(round(scenarios[1]["w"][1].value, 2)),
        r"$w^2_{22}$": str(round(scenarios[1]["w"][2].value, 2)),
        r"$u^2_{1}$": str(round(scenarios[1]["u"][0].value, 2)),
        r"$u^2_{2}$": str(round(scenarios[1]["u"][0].value, 2)),
        r"$v^2_{1}$": str(round(scenarios[1]["v"][0].value, 2)),
        r"$v^2_{2}$": str(round(scenarios[1]["v"][0].value, 2)),
    },
    name="value",
).to_frame()

	value
$w^2_{11}$	114.09
$w^2_{12}$	47.23
$w^2_{22}$	156.32
$u^2_{1}$	-0.0
$u^2_{2}$	-0.0
$v^2_{1}$	61.8
$v^2_{2}$	61.8

time: 8.26 ms (started: 2021-02-16 14:20:16 +08:00)

ΨΦ

HW 2¶

BAA99 Problem¶

Problem Setup¶

Scenario 1: High Demand¶

Scenario 2: Low Demand¶