BIM production revisited

# install dependencies and select solver
%pip install -q amplpy

SOLVER = "highs"

from amplpy import AMPL, ampl_notebook

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

Problem description

We consider BIM raw material planning, but now with more sophisticated pricing and acquisition protocols. There are now three suppliers, each of which can deliver the following materials:

• A: silicon, germanium and plastic

• B: copper

• C: all of the above

For the suppliers, the following conditions apply. Copper should be acquired in multiples of 100 gram, since it is delivered in sheets of 100 gram. Unitary materials such as silicon, germanium and plastic may be acquired in any number, but the price is in batches of 100. Meaning that 30 units of silicon with 10 units of germanium and 50 units of plastic cost as much as 1 unit of silicon but half as much as 30 units of silicon with 30 units of germanium and 50 units of plastic. Furthermore, supplier C sells all materials and offers a discount if purchased together: 100 gram of copper and a batch of unitary material cost just 7. This set price is only applied to pairs, meaning that 100 gram of copper and 2 batches cost 13.

The summary of the prices in € is given in the following table:

Supplier	Copper per sheet of 100 gram	Batch of units	Together
A	-	5	-
B	3	-	-
C	4	6	7

Next, for stocked products inventory costs are incurred, whose summary is given in the following table:

Copper per 10 gram	Silicon per unit	Germanium per unit	Plastic per unit
0.1	0.02	0.02	0.02

The holding price of copper is per 10 gram and the copper stocked is rounded up to multiples of 10 grams, meaning that 12 grams pay for 20.

The capacity limitations of the warehouse allow for a maximum of \(10\) kilogram of copper in stock at any moment, but there are no practical limitations to the number of units of unitary products in stock.

Recall that BIM has the following stock at the beginning of the year:

Copper	Silicon	Germanium	Plastic
480	1000	1500	1750

The company would like to have at least the following stock at the end of the year:

Copper	Silicon	Germanium	Plastic
200	500	500	1000

The goal is to build an optimization model using the data above and solve it to minimize the acquisition and holding costs of the products while meeting the required quantities for production. The production is made-to-order, meaning that no inventory of chips is kept.

from io import StringIO
import pandas as pd

demand_data = """
chip, Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec
logic, 88, 125, 260, 217, 238, 286, 248, 238, 265, 293, 259, 244
memory, 47, 62, 81, 65, 95, 118, 86, 89, 82, 82, 84, 66
"""

demand_chips = pd.read_csv(StringIO(demand_data), index_col="chip")
display(demand_chips)

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
chip
logic	88	125	260	217	238	286	248	238	265	293	259	244
memory	47	62	81	65	95	118	86	89	82	82	84	66

use = dict()
use["logic"] = {"silicon": 1, "plastic": 1, "copper": 4}
use["memory"] = {"germanium": 1, "plastic": 1, "copper": 2}
use = pd.DataFrame.from_dict(use).fillna(0).astype(int)
display(use)

	logic	memory
silicon	1	0
plastic	1	1
copper	4	2
germanium	0	1

demand = use.dot(demand_chips)
display(demand)

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
silicon	88	125	260	217	238	286	248	238	265	293	259	244
plastic	135	187	341	282	333	404	334	327	347	375	343	310
copper	446	624	1202	998	1142	1380	1164	1130	1224	1336	1204	1108
germanium	47	62	81	65	95	118	86	89	82	82	84	66

%%writefile BIMproduction_v1.mod

set periods ordered;
set products;
set supplying_batches;
set supplying_copper;
set unitary_products;

param delta{products, periods}; # demand
param pi{supplying_batches};    # aquisition price
param kappa{supplying_copper};  # price of copper sheet
param beta;
param gamma{products};          # unitary holding costs
param Alpha{products};          # initial stock
param Omega{products};          # desired stock
param copper_bucket_size;
param stock_limit;
param batch_size;
param copper_sheet_mass;

var y{periods, supplying_copper} integer >= 0;
var x{unitary_products, periods, supplying_batches} >= 0;
var ss{products, periods} >= 0;
var uu{products, periods} >= 0;
var bb{periods, supplying_batches} integer >= 0;
var pp{periods} integer >= 0;
var rr{periods} integer >= 0;

s.t. units_in_batches {t in periods, s in supplying_batches}:
    sum{p in unitary_products} x[p,t,s] <= batch_size * bb[t,s];
s.t. copper_in_buckets {t in periods}:
    ss['copper', t] <= copper_bucket_size * rr[t];
s.t. inventory_capacity {t in periods}:
    ss['copper', t] <= stock_limit;
s.t. pairs_in_batches {t in periods}:
    pp[t] <= bb[t,'C'];
s.t. pairs_in_sheets {t in periods}:
    pp[t] <= y[t,'C'];
s.t. bought {t in periods, p in products}:
    (if p == 'copper' then
        copper_sheet_mass * sum{s in supplying_copper} y[t,s]
    else
        sum{s in supplying_batches} x[p,t,s])
    == uu[p,t];

var acquisition_cost =
    sum{t in periods}(
        sum{s in supplying_batches} pi[s] * bb[t,s] +
        sum{s in supplying_copper} kappa[s] * y[t,s] -
        beta * pp[t]
    );
var inventory_cost =
    sum{t in periods}(
        gamma['copper'] * rr[t] + 
        sum{p in unitary_products} gamma[p] * ss[p,t]
    );

minimize total_cost: acquisition_cost + inventory_cost;
   
s.t. balance {p in products, t in periods}:
    (if t == first(periods) then
        Alpha[p]
    else
        ss[p,prev(t)]) +
    uu[p,t] == delta[p,t] + ss[p,t];
s.t. finish {p in products}:
    ss[p, last(periods)] >= Omega[p];

Overwriting BIMproduction_v1.mod

def BIMproduction_v1(
    demand,
    existing,
    desired,
    stock_limit,
    supplying_copper,
    supplying_batches,
    price_copper_sheet,
    price_batch,
    discounted_price,
    batch_size,
    copper_sheet_mass,
    copper_bucket_size,
    unitary_products,
    unitary_holding_costs,
):
    m = AMPL()
    m.read("BIMproduction_v1.mod")

    m.set["periods"] = demand.columns
    m.set["products"] = demand.index
    m.param["delta"] = demand

    m.set["supplying_batches"] = supplying_batches
    m.param["pi"] = price_batch

    m.set["supplying_copper"] = supplying_copper
    m.param["kappa"] = price_copper_sheet

    m.param["beta"] = price_batch["C"] + price_copper_sheet["C"] - discounted_price

    m.set["unitary_products"] = unitary_products
    m.param["gamma"] = unitary_holding_costs

    m.param["Alpha"] = existing
    m.param["Omega"] = desired

    m.param["batch_size"] = batch_size
    m.param["copper_bucket_size"] = copper_bucket_size
    m.param["stock_limit"] = stock_limit
    m.param["copper_sheet_mass"] = copper_sheet_mass

    return m

m1 = BIMproduction_v1(
    demand=demand,
    existing={"silicon": 1000, "germanium": 1500, "plastic": 1750, "copper": 4800},
    desired={"silicon": 500, "germanium": 500, "plastic": 1000, "copper": 2000},
    stock_limit=10000,
    supplying_copper=["B", "C"],
    supplying_batches=["A", "C"],
    price_copper_sheet={"B": 300, "C": 400},
    price_batch={"A": 500, "C": 600},
    discounted_price=700,
    batch_size=100,
    copper_sheet_mass=100,
    copper_bucket_size=10,
    unitary_products=["silicon", "germanium", "plastic"],
    unitary_holding_costs={"copper": 10, "silicon": 2, "germanium": 2, "plastic": 2},
)

m1.option["solver"] = SOLVER
m1.option["highs_options"] = "outlev=1"
m1.solve()

HiGHS 1.5.1: tech:outlev=1
Running HiGHS 1.5.1 [date: 2023-06-22, git hash: 93f1876]
Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
152 rows, 236 cols, 444 nonzeros
113 rows, 197 cols, 366 nonzeros
96 rows, 156 cols, 285 nonzeros

Solving MIP model with:
   96 rows
   156 cols (0 binary, 72 integer, 23 implied int., 61 continuous)
   285 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   -inf            inf                  inf        0      0      0         0     0.0s
 R       0       0         0   0.00%   109104          114244             4.50%        0      0      0        67     0.1s
 L       0       0         0   0.00%   109957.085557   110216             0.23%      926     61      0       331     0.9s

2.8% inactive integer columns, restarting
Model after restart has 92 rows, 147 cols (11 bin., 58 int., 21 impl., 57 cont.), and 270 nonzeros

         0       0         0   0.00%   109957.143248   110216             0.23%       32      0      0       466     0.9s
         0       0         0   0.00%   110065.335725   110216             0.14%       32     15      2       500     0.9s
 L       0       0         0   0.00%   110065.431752   110216             0.14%       48     16      2       503     1.0s

13.0% inactive integer columns, restarting
Model after restart has 71 rows, 118 cols (12 bin., 46 int., 7 impl., 53 cont.), and 218 nonzeros

         0       0         0   0.00%   110065.431752   110216             0.14%       16      0      0       695     1.0s
         0       0         0   0.00%   110065.431752   110216             0.14%       16     15      0       718     1.0s

13.8% inactive integer columns, restarting
Model after restart has 57 rows, 96 cols (3 bin., 42 int., 6 impl., 45 cont.), and 173 nonzeros

         0       0         0   0.00%   110091.208095   110216             0.11%       11      0      0       783     1.1s
         0       0         0   0.00%   110091.208095   110216             0.11%       11     11      0       796     1.1s

Solving report
  Status            Optimal
  Primal bound      110216
  Dual bound        110206.854428
  Gap               0.0083% (tolerance: 0.01%)
  Solution status   feasible
                    110216 (objective)
                    0 (bound viol.)
                    2.30926389122e-14 (int. viol.)
                    0 (row viol.)
  Timing            1.31 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             11
  LP iterations     1510 (total)
                    453 (strong br.)
                    347 (separation)
                    466 (heuristics)
HiGHS 1.5.1: optimal solution; objective 110216
1510 simplex iterations
11 branching nodes
absmipgap=9.14557, relmipgap=8.29786e-05

stock = m1.var["ss"].to_pandas().unstack()
stock.columns = stock.columns.get_level_values(1)
stock = stock.reindex(index=demand.index, columns=demand.columns)
display(stock)

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
silicon	912.0	787.0	527.0	310.0	72.0	15.0	0.0	0.0	0.0	56.0	54.0	500.0
plastic	1615.0	1428.0	1087.0	805.0	472.0	68.0	1.0	36.0	24.0	0.0	0.0	1000.0
copper	4354.0	3730.0	2528.0	1530.0	388.0	8.0	44.0	14.0	90.0	54.0	50.0	2042.0
germanium	1453.0	1391.0	1310.0	1245.0	1150.0	1032.0	946.0	857.0	775.0	693.0	609.0	543.0

import matplotlib.pyplot as plt, numpy as np

stock.T.plot(drawstyle="steps-mid", grid=True, figsize=(13, 4))
plt.xticks(np.arange(len(stock.columns)), stock.columns)
plt.show()

m1.var["pp"].to_pandas().T

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
pp.val	0.0	0.0	0.0	0.0	0.0	3.0	5.0	6.0	6.0	7.0	6.0	20.0

df = m1.var["uu"].to_pandas().unstack()
df.columns = df.columns.get_level_values(1)
df = df.reindex(index=demand.index, columns=demand.columns)
display(df)

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
silicon	0.0	0.0	0.0	0.0	0.0	229.0	233.0	238.0	265.0	349.0	257.0	690.0
plastic	0.0	0.0	0.0	0.0	0.0	0.0	267.0	362.0	335.0	351.0	343.0	1310.0
copper	0.0	0.0	0.0	0.0	0.0	1000.0	1200.0	1100.0	1300.0	1300.0	1200.0	3100.0
germanium	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

b = m1.var["bb"].to_pandas().unstack().T
b.index = b.index.get_level_values(1)
b = b.reindex(columns=demand.columns)
b.index.names = [None]
display(b)

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
A	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
C	0.0	0.0	0.0	0.0	0.0	3.0	5.0	6.0	6.0	7.0	6.0	20.0

y = m1.var["y"].to_pandas().unstack().T
y.index = y.index.get_level_values(1)
y = y.reindex(columns=demand.columns)
y.index.names = [None]
display(y)

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
B	0.0	0.0	0.0	0.0	0.0	7.0	7.0	5.0	7.0	6.0	6.0	11.0
C	0.0	0.0	0.0	0.0	0.0	3.0	5.0	6.0	6.0	7.0	6.0	20.0

x = m1.var["x"].to_pandas().unstack(1)
x.columns = x.columns.get_level_values(1)
x = x.reindex(columns=demand.columns)
x.index.names = ["materials", "supplier"]
display(x)

		Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
materials	supplier
germanium	A	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	C	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
plastic	A	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	C	0.0	0.0	0.0	0.0	0.0	0.0	267.0	362.0	335.0	351.0	343.0	1310.0
silicon	A	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	C	0.0	0.0	0.0	0.0	0.0	229.0	233.0	238.0	265.0	349.0	257.0	690.0