MIP 4: Production-Inventory Planning

Packages Used: No new packages are used in this notebook.

1. Inventory Costs in Production-Inventory Planning

In the production-inventory planning model presented in Ex 7 of Net Mod 3, the inventory cost for each stage i was given, $c_1^i = 5$ and $c_2^i = 25$ \$/ton:

It is more convenient to use an annual inventory carrying rate, $h$, multiplied by the cummulative production cost:

$ \begin{eqnarray*} \quad c_i^i &=& h\sum_{k=1}^{i} c_k^p \\ h &=& \mbox{per-period inventory carrying rate (\$/\$-period) .} \\ \end{eqnarray*} $

For example, in Ex 7, assuming an annual inventory carrying rate of 0.3 \$/\\$-yr, the monthly rate is $ h = 0.3/12$. The inventory cost per ton at stage i is then the monthly rate times the cumulative cost to produce each ton from stage 1 through i:

$ \begin{eqnarray*} \quad c_1^i &=& \dfrac{0.3}{12} 200 = 5 \\ c_2^i &=& \dfrac{0.3}{12} \bigl( 200 + 800 \bigr) = 25 \end{eqnarray*} $.

2. Single-Product Production-Inventory Planning with Fixed Cost

In the production-inventory planning model presented in Net Mod 3, the production capacity available at each stage for each period was given; as a result, the model could be solved as an LP. Now, the model will include making a decision at each stage for each period regarding how much, if any, capacity is available. The capacity is assumed to be available in discrete units that are either available at 100% of the unit's capacity or not available at all. Binary decision variables are used to make these decisions. This internal capacity decision requires the model to be solved as a MILP.

Fixed Production Cost and Economies of Scale

Cost data from existing facilities can be used to fit linear estimate. Economies of scale in production $ \implies k > 0 \mbox{ and } \beta < 1 $

MILP Model

$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{i=1}^{m} \sum_{j=1}^{t} c_i^p x_{ij} + \sum_{i=1}^{m} \sum_{j=2}^{t+1} c_i^i y_{ij} &+& \sum_{i=1}^{m} \sum_{j=1}^{t} c_i^k z_{ij} \\ \quad \mbox{subject to} \quad -x_{ij} + x_{(i+1)j} - y_{ij} + y_{i(j+1)} &=& 0 \mbox{,} &\quad i = 1, \dots, m-1; j = 1, \dots, t &\quad (a) \\ -x_{mj} - y_{mj} + y_{m(j+1)} &=& d_j \mbox{,} &\quad j = 1, \dots, t &\quad (b) \\ x_{ij} &\le& u_{i} z_{ij} \mbox{,} &\quad i = 1, \dots, m; j = 1, \dots, t &\quad (c) \\ y_{i,1} &=& y_i^0 \mbox{,} &\quad i = 1, \dots, m \\ y_{i(t+1)} &=& y_i^{t+1} \mbox{,} &\quad i = 1, \dots, m \\ x, y \ge 0; z &\in& \bigl\{ 0, 1 \bigr\} \mbox{,} \\ \end{eqnarray*} $

where,

$ \begin{eqnarray*} \quad m &=& \mbox{number of production stages} \\ t &=& \mbox{number of periods of production} \\ c_i^p &=& \mbox{variable production cost (dollar/ton) in stage } i \\ x_{ij} &=& \mbox{production (ton) at stage } i \mbox{ in period } j \\ c_i^i &=& \mbox{inventory cost (dollar/ton) in stage } i \\ y_{ij} &=& \mbox{stage-} i \mbox{ inventory (ton) from period } j-1 \mbox{ to } j \\ c_i^k &=& \mbox{fixed production cost (dollar/period) in stage } i \\ z_{ij} &=& \mbox{production indicator at stage } i \mbox{ in period } j \\ d_j &=& \mbox{demand (ton) in period } j \\ u_{i} &=& \mbox{production capacity (ton) of stage } i \mbox{ in period } j \\ y_i^0 &=& \mbox{initial inventory (ton) of stage } i \\ y_i^{t+1} &=& \mbox{final inventory (ton) of stage } i \mbox{.} \\ \end{eqnarray*} $

In the objective function, inventory costs are summed from 2 to t + 1 because the period 1 initial inventory cost was accounted for in period 0 prior to period 1 (rolling horizon). Constraints (a) and (b) in the model are the flow-balance equations for each node. Constraint (c) represents available production capacity and uses the binary variables z to indicate whether production occurs.

Ex 1: Single-Product Production-Inventory Planning with Fixed Cost

Question: A single product is produced in a two-stage production process. Stage one has a capacity to produce 40 tons of the product per month, and there is a fixed cost of \$6000 per month and a variable cost of \\$200 for each ton of product produced. For stage two, capacity is 80 tons per month, and there is a fixed cost of \$18,000 and a variable cost of \\$800. The fixed costs are only incurred when any amount of production occurs at that stage for that month. Demand for the next six months is estimated to be 25, 15, 10, 50, 25, and 15 tons per month. Currently, five tons of material are in storage at stage one and 35 tons at stage two; at the end of the planning period, ten tons of material should be in storage at stage one and 15 at stage two. Up to 30 tons of material per month can be stored at each stage. Assuming that the annual inventory carrying rate is $0.4 \mbox{ \$/\$-yr }$, determine the production plan that minimizes total costs over the planning horizon and what those total costs will be.

Solution: Can extend the MCNF model from Net Mod 3 by using the production indicator binary variable to turn capacity on or off at each stage for each period. One modeling issue is how to determine inventory cost on a per-ton basis if the production cost includes a fixed component. An approximation is to allocate the fixed cost based on capacity:

$ \quad c_i^i = h\sum_{j=1}^{i} \left( c_j^p + \dfrac{c_j^k}{u_j} \right) $

The maximum inventory storage constraint can be modeled as an upper bound on y.

u = [40, 80]                       # capacity of each stage (ton)
d = [25, 15, 10, 50, 25, 15]       # demand in each period (ton)
cp = [200, 800]                    # variable production cost for each stage ($/ton)
ck = [6000, 18000]                 # fixed production cost for each stage ($/month)
h = 0.4/12                         # monthly inventory carrying rate (1/month)
ci = cumsum(cp + ck./u)*h          # inventory cost for each stage ($/ton)
y0 = [5, 25]                       # initial inventory of each stage (ton)
yt = [10, 15]                      # final inventory of each stage (ton)
m = length(u)                      # number of production stages
t = length(d)                      # number of periods of production
Ylb = [y0 zeros(m,t-1) yt]         # same LB and UB to fix inital and final inventory levels
Yub = [y0 repeat([30],m,t-1) yt]   # max inventory storage level

2×7 Matrix{Int64}:
  5  30  30  30  30  30  10
 25  30  30  30  30  30  15

using JuMP, Cbc
model = Model(Cbc.Optimizer)
@variable(model, 0 <= X[i=1:m, j=1:t] )
@variable(model, Ylb[i,j] <= Y[i=1:m, j=1:t+1] <= Yub[i,j] )
@variable(model, Z[i=1:m, j=1:t], Bin )
@objective(model, Min, sum(cp[i]*X[i,j] for i=1:m, j=1:t ) + 
    sum(ci[i]*Y[i,j] for i=1:m, j=2:t+1 ) + sum(ck[i]*Z[i,j] for i=1:m, j=1:t ) )
for j = 1:t
    for i = 1:m-1                                      # (a) flow-balance all but last stage
        @constraint(model, -X[i,j] + X[i+1,j] - Y[i,j] + Y[i,j+1] == 0 )   
    end                                                # (b) flow-balance for last stage
    @constraint(model, -X[m,j] - Y[m,j] + Y[m,j+1] == -d[j] )              
    for i = 1:m                                        # (c) production capacity
        @constraint(model, X[i,j] <= u[i]*Z[i,j] )                         
    end
end
set_optimizer_attribute(model, "logLevel", 0)
set_optimizer_attribute(model, "seconds", 30.0)
optimize!(model)
println(solution_summary(model).termination_status)
TCᵒ, Xᵒ, Yᵒ, Zᵒ = objective_value(model), value.(X), value.(Y), round.(value.(Z))

OPTIMAL

(212275.0, [0.0 24.999999999999986 … 0.0 40.0; 0.0 29.999999999999986 … 0.0 29.999999999999996], [5.0 5.0 … 0.0 10.0; 25.0 0.0 … 0.0 15.0], [0.0 1.0 … 0.0 1.0; 0.0 1.0 … 0.0 1.0])

Optimal production flows (and required demands):

Matrix([Xᵒ; d'])

3×6 Matrix{Float64}:
  0.0  25.0  30.0  40.0   0.0  40.0
  0.0  30.0   0.0  70.0   0.0  30.0
 25.0  15.0  10.0  50.0  25.0  15.0

Optimal inventory flows:

Matrix(Yᵒ)

2×7 Matrix{Float64}:
  5.0  5.0   0.0  30.0   0.0  0.0  10.0
 25.0  0.0  15.0   5.0  25.0  0.0  15.0

Optimal production indicators:

Matrix(Zᵒ)

2×6 Matrix{Float64}:
 0.0  1.0  1.0  1.0  0.0  1.0
 0.0  1.0  0.0  1.0  0.0  1.0

Calculate total variable production, inventory, and fixed production cost components:

TCp, TCi, TCk = sum(sum(cp.*Xᵒ)), sum(sum(ci.*Yᵒ)), sum(sum(ck.*Zᵒ))

(130999.99999999999, 4479.166666666666, 78000.0)

TCp, TCi, TCk = sum((cp.*Xᵒ)[:]), sum((ci.*Yᵒ)[:]), sum((ck.*Zᵒ)[:])

(130999.99999999999, 4479.166666666666, 78000.0)

The sum of component costs exceeds total cost from solver because it includes first-period (time 0) inventory costs:

TCp+TCi+TCk, TCᵒ

(213479.16666666666, 212275.0)

Exclude first period from inventory costs:

TCi = sum(sum(ci.*Yᵒ[:,2:end]))

3274.9999999999986

Total component costs now match total cost from solver:

TCp+TCi+TCk, TCᵒ

(212274.99999999997, 212275.0)

Ex 2: Truckload Transport

A plant in Austin, TX, uses raw material from a supplier in Escanaba, MI. The raw material is transported via truckload, and the cost for each truckload is \$3 per mile, which is independent of how much is transported. Each truckload can carry up to 25 tons and 2750 cubic feet of material. Each unit of raw material costs \$2500, occupies 60 cubic feet, and weighs 800 pounds. A planning horizon of 13 weeks is used and is based on the following forecasted demand (in units):

21, 23, 9, 19, 41, 15, 62, 27, 11, 29, 23, 22, and 17

A maximum of 100 units can be stored at the plant. Currently, 13 units are in storage, and the same number should be in storage at the end of the planning horizon. Assuming that the annual inventory carrying rate is 0.4, determine the best transportation plan over the planning horizon.

d = [21, 23, 9, 19, 41, 15, 62, 27, 11, 29, 23, 22, 17]
dist = 1451                         # shortest-distance route via Google Maps
rTL = 3                             # TL cost per mile
cTL = rTL*dist                      # cost per TL
nCube = floor(2750/60)              # max cube units
nWt = floor(25*2000/800)            # max weight units
nTL = min(nCube,nWt)                # max unit is min of max cube and max weight

45.0

u = [nTL]                           # capacity of each stage (unit)
cp = [2500]                         # variable production cost for each stage ($/unit)
ck = [cTL]                          # fixed production cost for each stage ($)
h = 0.4/(365.25/7)                  # weekly inventory carrying rate (1/week)
ci = cumsum(cp + ck./u)*h           # inventory cost for each stage ($/unit)
y0 = [13]                           # initial inventory of each stage (unit)
yt = y0                             # final inventory of each stage (unit)
m = length(u)                       # number of production stages
t = length(d)                       # number of periods of production
Ylb = [y0 zeros(m,t-1) yt]          # same LB and UB to fix inital and final inventory levels
Yub = [y0 repeat([100],m,t-1) yt]   # max inventory storage level

1×14 Matrix{Int64}:
 13  100  100  100  100  100  100  100  100  100  100  100  100  13

using JuMP, Cbc
model = Model(Cbc.Optimizer)
@variable(model, 0 <= X[i=1:m, j=1:t] )
@variable(model, Ylb[i,j] <= Y[i=1:m, j=1:t+1] <= Yub[i,j] )
@variable(model, Z[i=1:m, j=1:t], Int )   # Int instead of Bin since can have multiple TL/wk
# Note: Unit prod cost not included in obj since does not impact soln, used only for inv cost
@objective(model, Min, 
    sum(ci[i]*Y[i,j] for i=1:m, j=2:t+1 ) + sum(ck[i]*Z[i,j] for i=1:m, j=1:t ) )
for j = 1:t
    for i = 1:m-1                                      # (a) flow-balance all but last stage
        @constraint(model, -X[i,j] + X[i+1,j] - Y[i,j] + Y[i,j+1] == 0 )   
    end                                                # (b) flow-balance for last stage
    @constraint(model, -X[m,j] - Y[m,j] + Y[m,j+1] == -d[j] )              
    for i = 1:m                                        # (c) production capacity
        @constraint(model, X[i,j] <= u[i]*Z[i,j] )                         
    end
end
set_optimizer_attribute(model, "logLevel", 0)
set_optimizer_attribute(model, "seconds", 30.0)
optimize!(model)
println(solution_summary(model).termination_status)
TCᵒ, Xᵒ, Yᵒ, Zᵒ = objective_value(model), value.(X), value.(Y), round.(value.(Z))

OPTIMAL

(37969.228820442615, [40.00000000000001 0.0 … 0.0 30.000000000000007], [13.0 32.00000000000001 … 0.0 13.0], [1.0 0.0 … 0.0 1.0])

# Report results
println(stdout,"\nTotal transport and inventory cost: ", TCᵒ)
println(stdout,"\nUnits transported/Demand: "), flush(stdout)
display(Matrix([Xᵒ; d']))
println(stdout,"\nInventory level: "), flush(stdout)
display(Matrix(Yᵒ))
println(stdout,"\nTL per week: "), flush(stdout)
display(Int.(Matrix(Zᵒ)))

Total transport and inventory cost: 37969.228820442615

Units transported/Demand: 

2×13 Matrix{Float64}:
 40.0   0.0  0.0  30.0  45.0   0.0  89.0   0.0  40.0   0.0  45.0   0.0  30.0
 21.0  23.0  9.0  19.0  41.0  15.0  62.0  27.0  11.0  29.0  23.0  22.0  17.0

Inventory level: 

1×14 Matrix{Float64}:
 13.0  32.0  9.0  0.0  11.0  15.0  0.0  27.0  0.0  29.0  0.0  22.0  0.0  13.0

TL per week: 

1×13 Matrix{Int64}:
 1  0  0  1  1  0  2  0  1  0  1  0  1

3. Multiple-product Production-Inventory Planning with Setup Cost

If multiple products are produced using shared resources, there is often a cost associated with reconfiguring or setting up the resource to produce a product if it was previously producing a different product; examples of such setup costs would be the material and labor cost associated with re-tooling/re-calibration of machines, cleaning/sterilization of equipment, and the lost production capacity while a resource is not available during set up. In the following model, set up cost is separate from production cost in the objective function:

$ \begin{eqnarray*} \quad \mbox{Minimize} \quad \sum_{i=1}^{m} \sum_{j=1}^{t} \sum_{k=1}^{g} c_{ik}^p x_{ijk} + \sum_{i=1}^{m} \sum_{j=2}^{t+1} \sum_{k=1}^{g} c_{ik}^i y_{ijk} &+& \sum_{i=1}^{m} \sum_{j=1}^{t} \sum_{k=1}^{g} c_{ik}^s v_{ijk} \\ \quad \mbox{subject to} \quad -x_{ijk} + x_{(i+1)jk} - y_{ijk} + y_{i(j+1)k} &=& 0 \mbox{,} & i = 1, \dots, m-1; j = 1, \dots, t; k = 1, \dots, g \;\;& (a) \\ -x_{mjk} - y_{mjk} + y_{m(j+1)k} &=& d_{jk} \mbox{,} & j = 1, \dots, t; k = 1, \dots, g & (b) \\ x_{ijk} &\le& u_{ik} z_{ijk} \mbox{,} & i = 1, \dots, m; j = 1, \dots, t; k = 1, \dots, g & (c) \\ -v_{i1k} + z_{i1k} &\le& z_{ik}^0 \mbox{,} & i = 1, \dots, m; k = 1, \dots, g & (d) \\ -v_{ijk} + z_{ijk} - z_{i(j-1)k} &\le& 0 \mbox{,} & i = 1, \dots, m; j = 2, \dots, t; k = 1, \dots, g & (e) \\ \sum_{k=1}^g z_{ijk} &=& 1 \mbox{,} & i = 1, \dots, m; j = 1, \dots, t & (f) \\ y_{i1k} &=& y_i^0 \mbox{,} & i = 1, \dots, m; k = 1, \dots, g \\ y_{i(t+1)k} &=& y_{ik}^{t+1} \mbox{,} & i = 1, \dots, m; k = 1, \dots, g \\ x, y \ge 0; z, v &\in& \bigl\{ 0, 1 \bigr\} \mbox{,} \\ \end{eqnarray*} $

where,

$ \begin{eqnarray*} \quad m &=& \mbox{number of production stages} \\ t &=& \mbox{number of periods of production} \\ g &=& \mbox{number of products} \\ c_{ik}^p &=& \mbox{production cost (dollar/ton) in stage } i \mbox{ for product } k \\ x_{ijk} &=& \mbox{production (ton) at stage } i \mbox{ in period } j \mbox{ for product } k \\ c_{ik}^i &=& \mbox{inventory cost (dollar/ton) in stage } i \mbox{ for product } k \\ y_{ijk} &=& \mbox{stage-} i \mbox{ inventory (ton) from period } j-1 \mbox{ to } j \mbox{ for product } k \\ c_{ik}^s &=& \mbox{setup cost (dollar) in stage } i \mbox{ for product } k \\ v_{ijk} &=& \mbox{setup indicator at stage } i \mbox{ in period } j \mbox{ for product } k \\ z_{ijk} &=& \mbox{production indicator at stage } i \mbox{ in period } j \mbox{ for product } k \\ d_{jk} &=& \mbox{demand (ton) in period } j \mbox{ for product } k \\ u_{ik} &=& \mbox{production capacity (ton) of stage } i \mbox{ in period } j \mbox{ for product } k \\ z_{ik}^0 &=& \mbox{initial setup at stage } i \mbox{ for product } k \\ y_{ik}^0 &=& \mbox{initial inventory (ton) of stage } i \mbox{ for product } k \\ y_{ik}^{t+1} &=& \mbox{final inventory (ton) of stage } i \mbox{ for product } k \mbox{.} \\ \end{eqnarray*} $

Constraints (a) and (b) are for flow balance, (c) for production capacity, (d) and (e) to force setups to occur, and (f) are linking constraints that ensure that only one product is produced at each stage during each period. Constraints (a)-(e) all refer to each product separately, and only (f) refer to multiple products together for each stage and period.

In order to ensure that a setup occurs only the first period that production switched to a different product, the constraint $ -v_t + z_t - z_{t-1} \le 0 $ can used. There are eight possible combinations of values that the three binary variables can take. In the table below, right, the combination $ [ 0 \; 1 \; 0 = 1 \nleq 0 ] $ is not feasible (i.e., allowing a product switch without incurring a setup), and this constraint ensures that this cannot happen. The combination $ [ 1 \; 1 \; 0 = 0 \leq 0 ] $ is feasible can will occur; the other three combinations where a setup is indicated ($ v_t = 1 $) are feasible but would never occur because they would only add an additional unnecessary setup cost and the objective is to minimize cost.

T = [0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1]
TT(f) = [T [f(row...) for row in eachrow(T)] [f(row...).<=0 for row in eachrow(T)]]
TT((vt,zt,zt_1) -> -vt + zt - zt_1)

8×5 Matrix{Int64}:
 0  0  0   0  1
 0  0  1  -1  1
 0  1  0   1  0
 0  1  1   0  1
 1  0  0  -1  1
 1  0  1  -2  1
 1  1  0   0  1
 1  1  1  -1  1

Ex 3: Multiple-Product Production-Inventory Planning with Setup Cost

Question: A plant can use a three-stage process to produce two products. A 13-week rolling horizon is used for planning production. The products’ forecasted demand, in tons, is the following, where D[j,k] is the demand in period j for product k:

D = [
  6  17
 40  14
 33  17
 24   4
 29   7
 24  10
 16  12
 13  18
 14  13
 17  23
 18   8
 43  17
 25   3
]

The plant can produce up to 60, 55, and 50 and 50, 45, and 35 tons per week for each stage of each product, respectively; has production costs of \$12, \\$75, and \$35 and \\$20, \$130, and \\$60 per ton for each stage of each product, respectively; and has costs of \$400, \\$90, and \$50 and \\$600, \$110, and \\$60 per set up for each stage of each product, respectively. The annual inventory carrying rate is 0.3. The plant is currently set up to produce each stage of the first product, and the first two weeks' demand for the second product is being held as finished product inventory. The final inventory for all stages of both products is zero. Determine the amount of each product that should be produced in order to minimize total costs over the planning horizon.

D = [6 40 33 24 29 24 16 13 14 17 18 43 25
    17 14 17  4  7 10 12 18 13 23  8 17  3]'  # demand in each period, each product (ton)
t = size(D,1)                         # number of periods
U = [60 50; 55 45; 50 35]             # capacity of each stage for each product (ton)
m = size(U,1)                         # number of production stages
g = size(U,2)                         # number of products
Cp = [12 20; 75 130; 35 60]           # production cost for each stage, each product ($/ton)
h = 0.3/(365.25/7)                    # weekly inventory carrying rate (1/wk)
Ci = cumsum(Cp,dims=1)*h              # inventory cost for each stage, each product ($/ton)
Cs = [400 600; 90 110; 50 70]         # setup cost for each stage for each product ($)
Y0 = [0 0; 0 0; 0 sum(D[[1 2],2])]    # initial inventory of each stage, each product (ton)

3×2 Matrix{Int64}:
 0   0
 0   0
 0  31

Yt = zeros(m,g)                       # final inventory of each stage for each product (ton)
Z0 = [1 0; 1 0; 1 0];                 # initial setup at each stage for each product
                                      # convert Y0 and Yt from 2D (m x g) to 3D (m x t x g)
Ylb = [reshape(Y0,m,1,g) zeros(m,t-1,g) reshape(Yt,m,1,g)]  
Yub = [reshape(Y0,m,1,g) repeat([Inf],m,t-1,g) reshape(Yt,m,1,g)]

3×14×2 Array{Float64, 3}:
[:, :, 1] =
 0.0  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  0.0
 0.0  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  0.0
 0.0  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  0.0

[:, :, 2] =
  0.0  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  0.0
  0.0  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  0.0
 31.0  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  Inf  0.0

using JuMP, Cbc
model = Model(Cbc.Optimizer)
@variable(model, 0 <= X[i=1:m, j=1:t, k=1:g] )
@variable(model, Ylb[i,j,k] <= Y[i=1:m, j=1:t+1, k=1:g] <= Yub[i,j,k] )
@variable(model, Z[i=1:m, j=1:t, k=1:g], Bin )
@variable(model, V[i=1:m, j=1:t, k=1:g], Bin )
@objective(model, Min, 
    sum(Cp[i,k]*X[i,j,k] for i=1:m, j=1:t, k=1:g ) + 
    sum(Ci[i,k]*Y[i,j,k] for i=1:m, j=2:t+1, k=1:g ) + 
    sum(Cs[i,k]*V[i,j,k] for i=1:m, j=1:t, k=1:g ) )
for k = 1:g
    for j = 1:t
        for i = 1:m-1                                  # (a) flow-balance all but last stage
            @constraint(model, -X[i,j,k] + X[i+1,j,k] - Y[i,j,k] + Y[i,j+1,k] == 0 )    
        end                                            # (b) flow-balance for last stage
        @constraint(model, -X[m,j,k] - Y[m,j,k] + Y[m,j+1,k] == -D[j,k] )               
        for i = 1:m                                    # (c) production capacity
            @constraint(model, X[i,j,k] <= U[i,k]*Z[i,j,k] )
        end
    end
    for i = 1:m                                         # (d) setup for first stage
        @constraint(model, -V[i,1,k] + Z[i,1,k] <= Z0[i,k] )
        for j = 2:t                                     # (e) setup all but first stage
            @constraint(model, -V[i,j,k] + Z[i,j,k] - Z[i,j-1,k] <= 0 )
        end
    end
end
for i = 1:m
    for j = 1:t                                         # (f) linking
        @constraint(model, sum(Z[i,j,k] for k=1:g) == 1 )              
    end
end
set_optimizer_attribute(model, "logLevel", 0)
set_optimizer_attribute(model, "seconds", 30.0)
optimize!(model)
println(solution_summary(model).termination_status)
TCᵒ,Xᵒ,Yᵒ,Zᵒ,Vᵒ = 
objective_value(model), value.(X), value.(Y), round.(value.(Z)), round.(value.(V))

OPTIMAL

(67314.92484599589, [35.000000000000036 50.0 … 0.0 0.0; 35.000000000000036 50.0 … 43.000000000000014 25.000000000000007; 35.000000000000036 49.99999999999999 … 43.00000000000001 25.000000000000007]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 25.000000000000007 0.0; 0.0 0.0 … 0.0 0.0; 0.0 29.000000000000036 … 0.0 0.0]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; 31.0 13.999999999999998 … 3.0 0.0], [1.0 1.0 … 0.0 0.0; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]

[0.0 0.0 … 1.0 1.0; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 1.0 0.0]

[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Matrix([Xᵒ[:,:,1]; D[:,1]'])   # Product 1 production level/Demand

4×13 Matrix{Float64}:
 35.0  50.0   0.0  50.0  50.0  57.0  60.0   0.0   0.0   0.0   0.0   0.0   0.0
 35.0  50.0   0.0  50.0  50.0  49.0   0.0   0.0   0.0   0.0   0.0  43.0  25.0
 35.0  50.0   0.0  50.0  50.0   0.0   0.0   0.0  49.0   0.0   0.0  43.0  25.0
  6.0  40.0  33.0  24.0  29.0  24.0  16.0  13.0  14.0  17.0  18.0  43.0  25.0

Matrix([Xᵒ[:,:,2]; D[:,2]'])   # Product 2 production level/Demand

4×13 Matrix{Float64}:
  0.0   0.0  50.0  0.0  0.0   0.0   0.0  31.0   0.0  23.0  28.0   0.0  0.0
  0.0   0.0  38.0  0.0  0.0   0.0  12.0  31.0   0.0  23.0  28.0   0.0  0.0
  0.0   0.0  28.0  0.0  0.0  10.0  12.0  31.0   0.0  23.0  28.0   0.0  0.0
 17.0  14.0  17.0  4.0  7.0  10.0  12.0  18.0  13.0  23.0   8.0  17.0  3.0

Matrix(Yᵒ[:,:,1])              # Product 1 inventory level

3×14 Matrix{Float64}:
 0.0   0.0   0.0  0.0   0.0   0.0   8.0  …  68.0  68.0  68.0  68.0  25.0  0.0
 0.0   0.0   0.0  0.0   0.0   0.0  49.0     49.0   0.0   0.0   0.0   0.0  0.0
 0.0  29.0  39.0  6.0  32.0  53.0  29.0      0.0  35.0  18.0   0.0   0.0  0.0

Matrix(Yᵒ[:,:,2])              # Product 2 inventory level

3×14 Matrix{Float64}:
  0.0   0.0  0.0  12.0  12.0  12.0  12.0  0.0   0.0  0.0  0.0   0.0  0.0  0.0
  0.0   0.0  0.0  10.0  10.0  10.0   0.0  0.0   0.0  0.0  0.0   0.0  0.0  0.0
 31.0  14.0  0.0  11.0   7.0   0.0   0.0  0.0  13.0  0.0  0.0  20.0  3.0  0.0

Matrix(Zᵒ[:,:,1])              # Product 1 production indicator

3×13 Matrix{Float64}:
 1.0  1.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  1.0  0.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0
 1.0  1.0  0.0  1.0  1.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0  1.0

Matrix(Vᵒ[:,:,1])              # Product 1 setup indicator

3×13 Matrix{Float64}:
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0  0.0

Matrix(Zᵒ[:,:,2])              # Product 2 production indicator

3×13 Matrix{Float64}:
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  1.0  1.0  1.0  0.0  1.0  1.0  0.0  0.0

Matrix(Vᵒ[:,:,2])              # Product 2 setup indicator

3×13 Matrix{Float64}:
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  1.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0

size(Cp),size(Xᵒ)   # Cp not compatible size with Xᵒ, need to reshape to (3, 1, 2)

((3, 2), (3, 13, 2))

TCp = sum((reshape(Cp,m,1,g).*Xᵒ)[:])   # Reshaped Cp broadcast for multiplication with Xᵒ
TCi = sum((reshape(Ci,m,1,g).*Yᵒ[:,2:end,:])[:])   # First period excluded
TCs = sum((reshape(Cs,m,1,g).*Vᵒ)[:])
TCp, TCi, TCs

(64564.00000000001, 390.92484599589324, 2360.0)

TCp+TCi+TCs, TCᵒ   # Check that production, inventory, and setup costs match total costs

(67314.9248459959, 67314.92484599589)