using Random, DataStructures, DataFrames, Statistics

function base_stock_two_echelon(qmax_r, qmax_dc, sch_r, sch_dc, t_stop)
    n = length(qmax_r)                  # Number of retailers

    E = PriorityQueue{Tuple{Symbol, Int, Int}, Float64}()  # Event queue
    for i in 1:n                        # Initial demand events for each retailer
        enqueue!(E, (:demand, i, 0), 0.0)
    end

    # State variables
    S = Dict(
        :q_r => copy(qmax_r),            # Inventory levels at retailers (vector)
        :q_dc => qmax_dc,                # Inventory level at DCr
    )

    # Output log
    out = DataFrame(t=Float64[], evt=Symbol[], loc=Int[], q_r=Vector{Int}[], q_dc=Int[])

    id = 0
    while !isempty(E)
        (evt, i, _), t = dequeue_pair!(E)  # Get next event and time

        t > t_stop && break              # Stop simulation if t > t_stop

        if evt == :demand                # Schedule next demand event for retailer i
            id += 1
            enqueue!(E, (:demand, i, id), sch_r[i][1](t))

            if S[:q_r][i] > 0            # Inventory available, fulfill demand
                S[:q_r][i] -= 1
                id += 1
                enqueue!(E, (:ret2dc_order, i, id), t)  # Immediate order placement
            end

        elseif evt == :ret2dc_order      # Retailer i places order to DC
            if S[:q_dc] > 0              # DC has inventory, fulfill retailer's order
                S[:q_dc] -= 1

                id += 1                  # Schedule arrival of inventory at retailer i
                enqueue!(E, (:dc2ret_ship, i, id), sch_r[i][2](t))

                id += 1                  # DC places immediate replenishment order
                enqueue!(E, (:sup2dc_ship, 0, id), sch_dc(t))
            else                         # DC is out of stock
                id += 1                  # Schedule sup2dc2ret crossdock for backorder
                enqueue!(E, (:sup2dc2ret_ship, i, id), sch_r[i][2](sch_dc(t)))
            end

        elseif evt == :dc2ret_ship        # Inventory arrives at retailer i from DC
            S[:q_r][i] += 1

        elseif evt == :sup2dc_ship        # DC's replenishment arrives from supplier
            S[:q_dc] += 1
        
        elseif evt == :sup2dc2ret_ship    # Backorder arrives at retailer i
            S[:q_r][i] += 1
        end

        push!(out, (t, evt, i, copy(S[:q_r]), S[:q_dc]))   # Log state
    end

    # Calculate metrics
    π₀_r = zeros(n)
    for i in 1:n
        demands = [row for row in eachrow(out) if row.evt == :demand && row.loc == i]
        stockouts = [row for row in demands if row.q_r[i] == 0]
        π₀_r[i] = length(stockouts) / length(demands)
    end

    # Average inventory levels
    q̄_r = [mean([row.q_r[i] for row in eachrow(out)]) for i in 1:n]
    q̄_dc = mean([row.q_dc for row in eachrow(out)])

    return π₀_r, q̄_r, q̄_dc, out
end

base_stock_two_echelon (generic function with 1 method)

qmax_r, qmax_dc = [8, 8], 30 
rₐ = [3.0, 3.0]
t_dc2ret = [1.0, 1.0]
t_sup2dc = 7.0

n = length(qmax_r)
sch_r = [[] for _ in 1:n]
for i = 1:n
    sch_r[i], rng = [], Xoshiro(1234+i)                                          
    push!(sch_r[i], t -> t + randexp(rng)/rₐ[i])
    push!(sch_r[i], t -> t + t_dc2ret[i])
end                                         
sch_dc = t -> t + t_sup2dc

t_stop = 1000.0
π₀_r, q̄_r, q̄_dc, out = base_stock_two_echelon(qmax_r, qmax_dc, sch_r, sch_dc, t_stop)
@show π₀_r
@show q̄_r
@show q̄_dc;

π₀_r = [0.3802485723883104, 0.4130360205831904]
q̄_r = [2.578978399528327, 2.5653642064640616]
q̄_dc = 4.166532668703436

out

filter(r -> r.evt == :sup2dc2ret_ship, out)

Row	t	evt	loc	q_r	q_dc
	Float64	Symbol	Int64	Array…	Int64
1	0.0	demand	1	[7, 8]	30
2	0.0	demand	2	[7, 7]	30
3	0.0	ret2dc_order	1	[7, 7]	29
4	0.0	ret2dc_order	2	[7, 7]	28
5	0.456446	demand	1	[6, 7]	28
6	0.456446	ret2dc_order	1	[6, 7]	27
7	0.572171	demand	1	[5, 7]	27
8	0.572171	ret2dc_order	1	[5, 7]	26
9	0.955705	demand	1	[4, 7]	26
10	0.955705	ret2dc_order	1	[4, 7]	25
11	1.0	dc2ret_ship	2	[4, 8]	25
12	1.0	dc2ret_ship	1	[5, 8]	25
13	1.11036	demand	1	[4, 8]	25
⋮	⋮	⋮	⋮	⋮	⋮
18646	999.23	sup2dc_ship	0	[0, 2]	6
18647	999.355	dc2ret_ship	1	[1, 2]	6
18648	999.427	sup2dc_ship	0	[1, 2]	7
18649	999.494	sup2dc_ship	0	[1, 2]	8
18650	999.551	demand	1	[0, 2]	8
18651	999.551	ret2dc_order	1	[0, 2]	7
18652	999.632	demand	2	[0, 1]	7
18653	999.632	ret2dc_order	2	[0, 1]	6
18654	999.692	dc2ret_ship	2	[0, 2]	6
18655	999.718	demand	1	[0, 2]	6
18656	999.735	sup2dc_ship	0	[0, 2]	7
18657	999.858	dc2ret_ship	1	[1, 2]	7

Row	t	evt	loc	q_r	q_dc
	Float64	Symbol	Int64	Array…	Int64
1	13.431	sup2dc2ret_ship	1	[2, 1]	13
2	13.47	sup2dc2ret_ship	1	[3, 1]	13
3	13.4785	sup2dc2ret_ship	1	[4, 1]	13
4	13.6191	sup2dc2ret_ship	1	[4, 0]	11
5	13.6471	sup2dc2ret_ship	2	[3, 1]	10
6	14.0319	sup2dc2ret_ship	2	[3, 3]	10
7	14.0383	sup2dc2ret_ship	2	[3, 4]	10
8	14.4642	sup2dc2ret_ship	1	[3, 3]	9
9	14.465	sup2dc2ret_ship	2	[3, 4]	9
10	14.8347	sup2dc2ret_ship	1	[5, 4]	9
11	15.515	sup2dc2ret_ship	1	[3, 4]	3
12	15.6855	sup2dc2ret_ship	2	[3, 6]	4
13	25.1261	sup2dc2ret_ship	2	[1, 1]	16
⋮	⋮	⋮	⋮	⋮	⋮
864	990.305	sup2dc2ret_ship	1	[1, 2]	14
865	990.387	sup2dc2ret_ship	1	[2, 2]	14
866	990.495	sup2dc2ret_ship	2	[2, 3]	15
867	990.56	sup2dc2ret_ship	2	[2, 4]	15
868	990.682	sup2dc2ret_ship	1	[2, 4]	14
869	990.683	sup2dc2ret_ship	1	[3, 4]	14
870	990.895	sup2dc2ret_ship	1	[3, 5]	13
871	991.048	sup2dc2ret_ship	1	[2, 5]	10
872	991.057	sup2dc2ret_ship	1	[3, 5]	10
873	991.263	sup2dc2ret_ship	1	[4, 5]	11
874	991.269	sup2dc2ret_ship	2	[4, 6]	11
875	991.689	sup2dc2ret_ship	2	[1, 7]	7

6. Inventory 3: Multi-Echelon Inventory Systems¶

Ex: Two-Echelon Supply Chain¶

Estimate Performance Measures¶