6. Inventory 2: Periodic Safety Stock¶

ISE 754, Fall 2024

Package Used: Functions from the following package are used in this notebook for the first time:

  • DataStructures: DataStructures.jl, Implements the priority queue data structure used to manage simulation events.

1. Base Stock with Lost Sales¶

Estimate Performance Measures¶

Estimate the out-of-stock $\pi_0$ and average-inventory $\overline{q}$ performance measures.

In [1]:
using Random, DataStructures, DataFrames, Statistics

qmax = 3                                             # Maximum inventory level
t_stop = 2.0                                         # Simulation stopping time
rₐ = 6.0                                             # Demand rate
tₑ = 1/8.0                                           # Order lead time

rng = Xoshiro(1234)                                  # Arrival time random number stream
sch = []                                             # Initialize scheduler array                    
push!(sch, t -> t + randexp(rng)/rₐ)                 # Time next demand event occurs
push!(sch, t -> t + tₑ)                              # Time next order arrives

E = PriorityQueue((:demand, 0) => 0.0)               # Event queue: intial event
S = Dict(:q => qmax)                                 # State: initial inventory
out = DataFrame(t=Float64[], evt=Symbol[], q=Int[])  # Log output

id = 0
while !isempty(E)
    (evt, idi), t = dequeue_pair!(E)                   # Get next event and remove from E
    
    t > t_stop && break                              # Stop simulation if t > t_stop

    if evt == :demand                                # Process next event
        id += 1
        enqueue!(E, (:demand, id), sch[1](t))        # Schedule next demand event
        if S[:q] > 0
            S[:q] -= 1
            id += 1
            enqueue!(E, (:order, id), sch[2](t))      # Schedule next order event
        end                                           # (only q > 0 => lost sales)
    elseif evt == :order
        S[:q] += 1
    end

    push!(out, (t, evt, S[:q]))                       # Log event and state

    println("\n$t, q = $(S[:q]), $evt-$idi")
    [println("    $t:\t$evt-$id") for ((evt, id), t) in E]
end
0.0, q = 2, demand-0
    0.09280528708070662:	demand-1
    0.125:	order-2

0.09280528708070662, q = 1, demand-1
    0.125:	order-2
    0.21780528708070662:	order-4
    0.2344699608729862:	demand-3

0.125, q = 2, order-2
    0.21780528708070662:	order-4
    0.2344699608729862:	demand-3

0.21780528708070662, q = 3, order-4
    0.2344699608729862:	demand-3

0.2344699608729862, q = 2, demand-3
    0.28535049374745236:	demand-5
    0.3594699608729862:	order-6

0.28535049374745236, q = 1, demand-5
    0.3594699608729862:	order-6
    0.41035049374745236:	order-8
    0.4205452126471838:	demand-7

0.3594699608729862, q = 2, order-6
    0.41035049374745236:	order-8
    0.4205452126471838:	demand-7

0.41035049374745236, q = 3, order-8
    0.4205452126471838:	demand-7

0.4205452126471838, q = 2, demand-7
    0.5395194331650394:	demand-9
    0.5455452126471838:	order-10

0.5395194331650394, q = 1, demand-9
    0.5455452126471838:	order-10
    0.5984192148045301:	demand-11
    0.6645194331650394:	order-12

0.5455452126471838, q = 2, order-10
    0.5984192148045301:	demand-11
    0.6645194331650394:	order-12

0.5984192148045301, q = 1, demand-11
    0.6645194331650394:	order-12
    0.7234192148045301:	order-14
    0.9195563434894847:	demand-13

0.6645194331650394, q = 2, order-12
    0.7234192148045301:	order-14
    0.9195563434894847:	demand-13

0.7234192148045301, q = 3, order-14
    0.9195563434894847:	demand-13

0.9195563434894847, q = 2, demand-13
    0.9580777815169099:	demand-15
    1.0445563434894847:	order-16

0.9580777815169099, q = 1, demand-15
    1.0061293889135055:	demand-17
    1.0445563434894847:	order-16
    1.08307778151691:	order-18

1.0061293889135055, q = 0, demand-17
    1.0445563434894847:	order-16
    1.077386895781228:	demand-19
    1.08307778151691:	order-18
    1.1311293889135055:	order-20

1.0445563434894847, q = 1, order-16
    1.077386895781228:	demand-19
    1.08307778151691:	order-18
    1.1311293889135055:	order-20

1.077386895781228, q = 0, demand-19
    1.08307778151691:	order-18
    1.1311293889135055:	order-20
    1.202386895781228:	order-22
    1.4765925327935197:	demand-21

1.08307778151691, q = 1, order-18
    1.1311293889135055:	order-20
    1.202386895781228:	order-22
    1.4765925327935197:	demand-21

1.1311293889135055, q = 2, order-20
    1.202386895781228:	order-22
    1.4765925327935197:	demand-21

1.202386895781228, q = 3, order-22
    1.4765925327935197:	demand-21

1.4765925327935197, q = 2, demand-21
    1.6015925327935197:	order-24
    1.7365716632050692:	demand-23

1.6015925327935197, q = 3, order-24
    1.7365716632050692:	demand-23

1.7365716632050692, q = 2, demand-23
    1.7379633350028847:	demand-25
    1.8615716632050692:	order-26

1.7379633350028847, q = 1, demand-25
    1.7526948698151428:	demand-27
    1.8615716632050692:	order-26
    1.8629633350028847:	order-28

1.7526948698151428, q = 0, demand-27
    1.8615716632050692:	order-26
    1.8629633350028847:	order-28
    1.8776948698151428:	order-30
    1.9154229178912858:	demand-29

1.8615716632050692, q = 1, order-26
    1.8629633350028847:	order-28
    1.8776948698151428:	order-30
    1.9154229178912858:	demand-29

1.8629633350028847, q = 2, order-28
    1.8776948698151428:	order-30
    1.9154229178912858:	demand-29

1.8776948698151428, q = 3, order-30
    1.9154229178912858:	demand-29

1.9154229178912858, q = 2, demand-29
    2.040422917891286:	order-32
    2.0768947413354795:	demand-31
In [2]:
first(out, 5)
Out[2]:
5×3 DataFrame
Rowtevtq
Float64SymbolInt64
10.0demand2
20.0928053demand1
30.125order2
40.217805order3
50.23447demand2
In [3]:
π₀ = sum(out.q .== 0)/nrow(out)
q̄ = mean(out.q)
π₀, q̄
Out[3]:
(0.0967741935483871, 1.7096774193548387)
In [4]:
using CairoMakie

lines(out.t, out.q, axis=(xlabel="Time", ylabel="Inventory"))
Out[4]:
No description has been provided for this image
In [5]:
using Random, DataStructures, DataFrames, Statistics

function base_stock(qmax, sch, t_stop)

    E = PriorityQueue((:demand, 0) => 0.0)               # Event queue: intial event
    S = Dict(:q => qmax)                                 # State: initial inventory
    out = DataFrame(t=Float64[], evt=Symbol[], q=Int[])  # Output log

    id = 0
    while !isempty(E)
        (evt, _), t = dequeue_pair!(E)                   # Get next event, remove from E
        
        t > t_stop && break                              # Stop simulation if t > t_stop
    
        if evt == :demand                                # Process next event
            id += 1
            enqueue!(E, (:demand, id), sch[1](t))
            if S[:q] > 0
                S[:q] -= 1
                id += 1
                enqueue!(E, (:order, id), sch[2](t))
            end
        elseif evt == :order
            S[:q] += 1
        end
    
        push!(out, (t, evt, S[:q]))                       # Log event and state
    end
    
    return sum(out.q .== 0)/nrow(out), mean(out.q), out
end
Out[5]:
base_stock (generic function with 1 method)
In [6]:
qmax, t_stop = 3, 2.0
rₐ, tₑ = 6.0, 1/8.0

sch, rng = [], Xoshiro(1234)                                          
push!(sch, t -> t + randexp(rng)/rₐ)
push!(sch, t -> t + tₑ)

π₀, q̄, out = base_stock(qmax, sch, t_stop)
π₀, q̄
Out[6]:
(0.0967741935483871, 1.7096774193548387)
In [7]:
t_stop = 20.0
π₀, q̄, out = base_stock(qmax, sch, t_stop)
π₀, q̄
Out[7]:
(0.05286343612334802, 1.907488986784141)
In [8]:
lines(out.t, out.q, axis=(xlabel="Time", ylabel="Inventory"))
Out[8]:
No description has been provided for this image
In [9]:
π₀, q̄, out = base_stock(2qmax, sch, t_stop)
π₀, q̄
Out[9]:
(0.0, 4.804)
In [10]:
lines(out.t, out.q, axis=(xlabel="Time", ylabel="Inventory"))
Out[10]:
No description has been provided for this image

Optimal Base-Stock Policy¶

In [11]:
p, c, h = 6, 3, .3             # Unit price and cost, and carrying rate
rₐ, tₑ = 6.0, 1/8.0            # Demand rate and order time
t_stop = 2000.0                # Simulation stopping time

TPh(π₀, q̄) = (p - c)*(1 - π₀)*rₐ - c*h*q̄    # Total profit function

UBqmax = 2round(Int, rₐ)                   # Upper bound for qmax

TPᵒ, qmaxᵒ = -Inf, -1
for i in 0:UBqmax                          # Run base_stock simulation
    sch, rng = [], Xoshiro(1234)
    push!(sch, t -> t + randexp(rng)/rₐ)
    push!(sch, t -> t + tₑ)

    π₀, q̄, _ = base_stock(i, sch, t_stop)
    TP = TPh(π₀, q̄)
    if TP > TPᵒ
        qmaxᵒ, TPᵒ = i, TP
    end
    println(i, "\t", TP, "\t", π₀, "\t", q̄)
end

qmaxᵒ, TPᵒ  # Optimal qmax and total profit
0	0.0	1.0	0.0
1	6.21531863389991	0.6365310740409409	0.3634689259590591
2	12.18713873473918	0.2734295227524972	0.9901442841287459
3	14.84114828513786	0.08439811701412239	1.8218728984532615
4	15.157059921210863	0.01961434791623471	2.766535351441012
5	14.55543424317618	0.0037220843672456576	3.7528535980148883
6	13.713550493985366	0.0006200653135463603	4.750320367078666
7	12.82358630952381	8.267195767195767e-5	5.749917328042328
8	11.925074404761904	0.0	6.749917328042328
9	11.025074404761906	0.0	7.749917328042328
10	10.125074404761907	0.0	8.749917328042327
11	9.225074404761907	0.0	9.749917328042327
12	8.325074404761907	0.0	10.749917328042327
Out[11]:
(4, 15.157059921210863)

Random Order Lead Times¶

In [12]:
TPᵒ, qmaxᵒ = -Inf, -1
for i in 0:UBqmax                          # Run base_stock simulation
    sch, rng1, rng2 = [], Xoshiro(1234), Xoshiro(5678)
    push!(sch, t -> t + randexp(rng1)/rₐ)
    push!(sch, t -> t + randexp(rng2)*tₑ)

    π₀, q̄, _ = base_stock(i, sch, t_stop)
    TP = TPh(π₀, q̄)
    if TP > TPᵒ
        qmaxᵒ, TPᵒ = i, TP
    end
    println(i, "\t", TP, "\t", π₀, "\t", q̄)
end

qmaxᵒ, TPᵒ  # Optimal qmax and total profit
0	0.0	1.0	0.0
1	6.188330871491876	0.638109305760709	0.361890694239291
2	12.140508015453618	0.2763000133220836	0.9845463830543097
3	14.809035840484604	0.08640417297661114	1.817432273262662
4	15.141027130175063	0.020700240603999005	2.7626317099477307
5	14.56536839493922	0.0035557760688001323	3.7451418175804183
6	13.723809523809521	0.00045469576719576717	4.742228835978836
7	12.831994047619048	0.0	5.742228835978836
8	11.931994047619048	0.0	6.742228835978836
9	11.031994047619047	0.0	7.742228835978836
10	10.131994047619049	0.0	8.742228835978835
11	9.231994047619049	0.0	9.742228835978835
12	8.33199404761905	0.0	10.742228835978835
Out[12]:
(4, 15.141027130175063)

2. Order Point with Lost Sales¶

Estimate Performance Measures¶

Estimate the out-of-stock $\pi_0$, average-inventory $\overline{q}$, and average number of orders $\overline{n}_o$ performance measures.

In [13]:
using Random, DataStructures, DataFrames, Statistics

function order_point(qmax, qmin, sch, t_stop)

    E = PriorityQueue((:demand, 0, 1) => 0.0)            # Event queue: intial demand
    S = Dict(:q => qmax, :p => qmax)                     # State: on-hand and position
    out = DataFrame(t=Float64[], evt=Symbol[], q=Int[], p=Int[], qₒ=Int[])
    
    id, nₒ = 0, 0
    while !isempty(E)
        (evt, _, qₒ), t = dequeue_pair!(E)               # Get next event, remove from E
        
        t > t_stop && break                              # Stop simulation if t > t_stop
    
        if evt == :demand                                # Process next event
            id += 1
            enqueue!(E, (:demand, id, 1), sch[1](t))
            if S[:q] > 0
                S[:q] -= 1
                S[:p] -= 1
                id += 1
                if S[:p] < qmin
                    enqueue!(E, (:order, id, qmax - S[:p]), sch[2](t))
                    S[:p] = qmax
                end
            end
        elseif evt == :order
            S[:q] += qₒ                                  # Order size                             
            nₒ += 1
        end
    
        push!(out, (t, evt, S[:q], S[:p], qₒ))           # Log event and state
    end
    
    return sum(out.q .== 0)/nrow(out), mean(out.q), nₒ/t_stop, out
end
Out[13]:
order_point (generic function with 1 method)
In [14]:
qmax, qmin, t_stop = 5, 2, 2.0              # Max and min inv level and sim stop time
rₐ, tₑ = 6.0, 1/8.0                         # Demand rate and order time

sch, rng = [], Xoshiro(1234)
push!(sch, t -> t + randexp(rng)/rₐ)
push!(sch, t -> t + tₑ)

π₀, q̄, n̄ₒ, out = order_point(qmax, qmin, sch, t_stop)
π₀, q̄, n̄ₒ
Out[14]:
(0.21052631578947367, 2.210526315789474, 1.5)
In [15]:
out
Out[15]:
19×5 DataFrame
Rowtevtqpqₒ
Float64SymbolInt64Int64Int64
10.0demand441
20.0928053demand331
30.23447demand221
40.28535demand151
50.41035order554
60.420545demand441
70.539519demand331
80.598419demand221
90.919556demand151
100.958078demand041
111.00613demand041
121.04456order444
131.07739demand331
141.47659demand221
151.73657demand151
161.73796demand041
171.75269demand041
181.86157order444
191.91542demand331
In [16]:
using CairoMakie

dcf() = display(current_figure())

lines(out.t, out.q, axis=(xlabel="Time", ylabel="Inventory"), label="Inv On Hand")
lines!(out.t, out.p, label="Inv Position")

axislegend(position=:lb)
dcf();
No description has been provided for this image

Optimal Order-Point Policy¶

In [17]:
p, c, h, cₒ = 6, 5, .1, 2       # Unit price and cost, carrying rate, and order cost
rₐ, tₑ = 6.0, 1/8.0             # Demand rate and order time
t_stop = 2000.0                 # Sim stop time

TPh(π₀, q̄, n̄ₒ) = (p - c)*(1 - π₀)*rₐ - c*h*q̄ - cₒ*n̄ₒ   # Total profit function

UBqmax = 2round(Int, rₐ)                   # Upper bound for qmax

TPᵒ, qmaxᵒ, qminᵒ = -Inf, -1, -1
for i in 0:UBqmax
    for j in 0:i
        sch, rng = [], Xoshiro(1234)
        push!(sch, t -> t + randexp(rng)/rₐ)
        push!(sch, t -> t + tₑ)
        
        π₀, q̄, n̄ₒ, _ = order_point(i, j, sch, t_stop)
        TP = TPh(π₀, q̄, n̄ₒ)
        if TP > TPᵒ
            qmaxᵒ, qminᵒ, TPᵒ = i, j, TP
        end
        println(i, "\t", j, "\t", TP, "\t", π₀, "\t", q̄, "\t", n̄ₒ)
    end
end

qmaxᵒ, qminᵒ, TPᵒ  # Optimal qmax. qmin, and total profit
0	0	0.0	1.0	0.0	0.0
1	0	0.0	1.0	0.0	0.0
1	1	-4.907920907225175	0.6365310740409409	0.3634689259590591	3.4535
2	0	0.00045469576719587646	0.999917328042328	8.267195767195767e-5	0.0
2	1	-1.5947231926249397	0.4671882581271228	0.7991872877244056	2.196
2	2	-6.564649278579356	0.2734295227524972	0.9901442841287459	5.2145
3	0	0.0008680555555557741	0.999834656084656	0.000248015873015873	0.0
3	1	-0.06678875905737947	0.36875775181147596	1.262484496377048	1.6115
3	2	-1.5329270447420917	0.23328583642063266	1.3684240524365916	2.7245
3	3	-7.113325151311365	0.08439811701412239	1.8218728984532615	5.848
4	0	0.0012400793650796929	0.9997519841269841	0.000496031746031746	0.0
4	1	0.757801885503484	0.3053012706653914	1.7367809810083346	1.271
4	2	0.24040081993062135	0.17798801639861242	1.8653421633554084	1.8795
4	3	-1.4740277469478356	0.06576026637069922	2.310932297447281	2.962
4	4	-7.519953763217915	0.01961434791623471	2.766535351441012	6.0095
5	0	0.0015707671957669667	0.9996693121693122	0.0008267195767195767	0.0
5	1	1.2269314914011837	0.2624753312658585	2.21243304200733	1.046
5	2	1.0970909151681263	0.1419881008088776	2.3759609599572165	1.4315
5	3	0.3200133831310308	0.05365701836290072	2.778089013383131	1.9845
5	4	-1.7411028924707441	0.015842349304482226	3.2520975932877016	3.01
5	5	-7.982759305210918	0.0037220843672456576	3.7528535980148883	6.042
6	0	0.0018601190476189278	0.9995866402116402	0.001240079365079365	0.0
6	1	1.501567623983301	0.22385373929316923	2.7166198805153674	0.8985
6	2	1.5241850926504266	0.12228468318412103	2.858213616489694	1.1565
6	3	1.1138992772362575	0.0438962933492474	3.2754459253365162	1.4925
6	4	0.026382477041449803	0.011106974435343758	3.773951352692976	2.01
6	5	-2.189369108648917	0.0025908163827793396	4.257648420704482	3.0225
6	6	-8.473880575420612	0.0006200653135463603	4.750320367078666	6.0475
7	0	0.00210813492063491	0.9995039682539683	0.001736111111111111	0.0
7	1	1.6394100461301901	0.1998974884674526	3.2004100461301896	0.7805
7	2	1.7549336371404163	0.10139561378524636	3.377385360296212	0.974
7	3	1.497882738629305	0.03768376009386431	3.766029401615018	1.1965
7	4	0.7987167813658025	0.010322922181048173	4.246691371095818	1.508
7	5	-0.41568150812352744	0.0020463847203274215	4.746806399603125	2.015
7	6	-2.676361186132392	0.0004960590861489279	5.252769663230998	3.0235
7	7	-8.971454695767196	8.267195767195767e-5	5.749917328042328	6.048
8	0	0.0023148148148149136	0.9994212962962963	0.0023148148148148147	0.0
8	1	1.6950885888113965	0.17940347232527082	3.6929811544739577	0.691
8	2	1.8461905626134298	0.09067150635208711	3.861560798548094	0.8395
8	3	1.6828532917139607	0.03213677639046538	4.248652099886493	1.0
8	4	1.1668344704017635	0.007304803252704844	4.748673420164014	1.2075
8	5	0.3367790859183808	0.0016535485151134334	5.260599245981877	1.5115
8	6	-0.9004425497612703	0.00018602343895330811	5.7366528182551	2.0155
8	7	-3.1714292328042335	5.511463844797178e-5	6.24619708994709	3.024
8	8	-9.470958664021165	0.0	6.749917328042328	6.048
9	0	0.002480158730158938	0.9993386243386243	0.002976190476190476	0.0
9	1	1.6890641871625667	0.16309238152369526	4.184763047390522	0.62
9	2	1.8534219141719512	0.08457454652706091	4.346261613331367	0.733
9	3	1.7325966690803765	0.02896451846488052	4.75923244026068	0.857
9	4	1.3134182024383323	0.0070172951516869865	5.264956053303091	1.006
9	5	0.6995072343943782	0.00137797988149373	5.74844977263332	1.209
9	6	-0.16024452675441525	0.00026456776241814934	6.271314240359812	1.5115
9	7	-1.403031746031746	6.200396825396825e-5	6.741319444444445	2.016
9	8	-3.671098544973545	0.0	7.24619708994709	3.024
9	9	-9.970958664021165	0.0	7.749917328042328	6.048
10	0	0.002604166666666984	0.9992559523809523	0.003720238095238095	0.0
10	1	1.6393509833585473	0.14969742813918305	4.676928895612708	0.562
10	2	1.816915572792363	0.07331443914081145	4.862395584725537	0.656
10	3	1.7102966095462238	0.025005515922629992	5.277340589835993	0.7505
10	4	1.3574900882650853	0.0057878744031254524	5.763565330632325	0.863
10	5	0.8627910997732422	0.0010629251700680273	6.229662698412699	1.008
10	6	0.21362142611091972	0.00013778849466069583	6.733103685842232	1.2095
10	7	-0.6500582010582017	6.613756613756614e-5	7.251322751322752	1.512
10	8	-1.9026597222222223	0.0	7.741319444444445	2.016
10	9	-4.171098544973545	0.0	8.24619708994709	3.024
10	10	-10.470958664021165	0.0	8.749917328042327	6.048
11	0	0.002686838624338385	0.9991732804232805	0.004546957671957672	0.0
11	1	1.5600502818832849	0.13652293158616485	5.181624257199451	0.515
11	2	1.7350512048192768	0.06573795180722891	5.373042168674699	0.592
11	3	1.643594296776115	0.025165661529297893	5.740823468096195	0.6675
11	4	1.3212876672056446	0.0058797589298838744	6.266867558430104	0.755
11	5	0.8841553208420747	0.0012298343340808797	6.76293134630688	0.8635
11	6	0.3691686507936507	0.0	7.229662698412699	1.008
11	7	-0.2862752325180855	6.889424733034792e-5	7.733723734068206	1.2095
11	8	-1.1496613756613758	0.0	8.251322751322752	1.512
11	9	-2.4026597222222223	0.0	8.741319444444445	2.016
11	10	-4.671098544973545	0.0	9.24619708994709	3.024
11	11	-10.970958664021165	0.0	9.749917328042327	6.048
12	0	0.002728174603174473	0.9990906084656085	0.005456349206349206	0.0
12	1	1.4490361879935598	0.12865138388407574	5.664111017403972	0.4735
12	2	1.6223188097768328	0.059738879611355704	5.882495825110065	0.539
12	3	1.5385563909774438	0.01819548872180451	6.296541353383459	0.602
12	4	1.2616240511980945	0.0043161184700104185	6.740958475963685	0.671
12	5	0.8599120305724999	0.0011758653634158889	7.244065554494011	0.7555
12	6	0.3881747685185186	0.00021701388888888888	7.765046296296297	0.864
12	7	-0.13083134920634887	0.0	8.229662698412698	1.008
12	8	-0.7858618670341024	0.0	8.733723734068205	1.2095
12	9	-1.6496613756613758	0.0	9.251322751322752	1.512
12	10	-2.9026597222222223	0.0	9.741319444444445	2.016
12	11	-5.171098544973545	0.0	10.24619708994709	3.024
12	12	-11.470958664021165	0.0	10.749917328042327	6.048
Out[17]:
(9, 2, 1.8534219141719512)