NETW 707 Modeling and Simulation Amr El Mougy Maggie Mashaly

1. NETW 707 Modeling and Simulation Amr El Mougy Maggie Mashaly

2. Lecture(3) Simulation Examples – Cont’d

3. Simulation of Queuing Systems

4. Single Server Queue • Dynamic, event-based model • A queuing system is described by

5. Simple Single Channel Queuing System Calling Population Server Waiting Line

6. Simple Single Channel Queuing System Assumptions Calling Population Server Waiting Line • Once units join the queue they will eventually be served • Service times are of random length according to a probability distribution that does not change with time • The system capacity has no limit, i.e. the queue can be infinite • Unlimited potential calling population • Constant arrival rate, i.e. no matter how many units arrived or are still in the system, inter-arrival times follow the same statistical distribution • Arrivals for service occur one at a time in a random manner • Units are served in the order of their arrival, i.e. First In First Out (FIFO) or First Come First Served (FCFS) • Single server or multiple parallel servers

7. Important Note • For any single or multi-channel queue, the overall effective arrival rate must be lower than the service rate. Otherwise, the waiting line will grow without bound • When queues grow without bound they are termed “explosive” or “unstable”

8. System State, Events and Simulation Clock for a Queuing Model • System state: • Number of units in the system • Server status: busy or idle • Events: • Arrival event • Service beginning event • Departure event, i.e. service completion • The simulation clock is used to track simulation time

9. Arrival Event Flow Diagram Arrival event Schedule the next arrival event No Server Busy? Yes Add 1 to the number in queue Set delay = 0 for this customer and gather statistics Queue full? Yes Add 1 to number of customers delayed Write error msg and stop sim No Make server busy Stop time of arrival of this customer Schedule a departure event for this customer Arrival event

10. Departure Event Flow Diagram Departure event No Yes Subtract 1 from the number in queue Is queue empty? Make server idle Compute delay of customer entering service and gather statistics Eliminate departure event from consideration Add 1 to the number of customers delayed Schedule a departure event for this customer Move each customer in queue (if any) up one place Return

11. Potential Unit Actions Upon Arrival

12. Server Outcomes after the Completion of a Service

13. Example of Single Server Queue Simulation • A small grocery store has one check out counter. Customers arrive at the check counter at random times that range from 1 to 8 minutes apart. Assume that inter-arrival times are integer valued, with each of the 8 values having equal probability. The service times vary from 1 to 6 minutes – also integer valued – with probabilities shown in the table below • Analyze the system by simulating the arrival and departure of 100 customers • Compute the measures of performance of the queuing model

14. Model Responses and Simulation Table for the Grocery Store Simulation

15. Typical Performance Measures of a Queuing System = = 0.255

16. Typical Performance Measures of a Queuing System =

17. Typical Performance Measures of a Queuing System Average time a customer spends in the system = average time a customer waits in queue + average time a customer spends in service = 1.32 + 3.27 = 4.59 min

18. Typical Performance Measures of a Queuing System

19. Typical Performance Measures of a Queuing System The experiment was run 50 times, each trial represents one day

20. Start Implementation using C Initialization routine Main program Timing routine 1 Set simulation clock = 0 Initialize system state and statistical counters Initialize event list 0 0. Invoke the initialization routine Invoke the timing routine Invoke event routing i Determine the next event type i Advance the simulation clock i Repeatedly 2 Event routine i Library routines Update system state Update statistical counters Generate future events and add to event list Generate random variables Simulation over? No Yes Report generator Compute estimates of interest Write report Stop

21. Implementation using C /*External definitions for single server queuing system */ #include <stdio.h> #include <math.h> #include <lcgrand.h> /*header file for RNG .*/ #define Q_LIMIT 100 /*limit on queue length*/ #define BUSY 1 /*Mnemonics for server being busy*/ #define IDLE 0 /*and idle */ intnext_event_type, num_custs_delayed, num_delays_required, num_events, num_in_q, server status; float area_num_in_q, area_server_status, mean_interarrival, mean_service, sim_time, time_arrival[Q_LIMIT+1], time_last_event, time_next_event[3], total_of_delays; FILE *infile, *outfile; void initialize(void); void timing(void); void arrive(void); void depart(void); void report(void); void update_time_avg_stats(void); float expon(float mean);

22. Implementation using C main() /*main function*/ { /*open input and output files*/ infile = fopen(“mm1.in”, “r”); outfile = fopen(“mm1.out”, “w”); /*specify the number of events for the timing function*/ num_events = 2; /*read input parameters*/ fscanf(infile, “%f %f %d”, &mean_interarrival, &mean_service, &num_delays_required); /*write report heading and input parameters*/ fprintf(outfile, “Single server queuing system\n\n”); fprintf(outfile, “Mean interarrival time%11.3f minutes\n\n”, mean_interarrival); fprintf(outfile, “Mean service time%16.3f minutes\n\n”, mean_service); fprintf(outfile, “Number of customers%14d\n\n”, num_delays_required); /*Initialize the simulation*/ Initialize(); /*Run the simulation while more delays are still needed*/ While(num_custs_delays < num_delays_required) { /*Determine the next event*/ timing(); /*Update time-average statistical accumulators*/ update_time_avg_stats(); /*Invoke the appropriate event function*/ switch(next_event_type) { case 1: arrive(); break; case 2: depart(); break; } } /*Invoke the report generator and end the simulation*/ report() fclose(infile); fclose(outfile); return0; }

23. Implementation using C void initialize(void) /*Initialize function*/ { /*Initialize the simulation clock*/ sim_time = 0.0; /*Initialize the state variables*/ server_status = IDLE; num_in_q = 0; time_last_event = 0.0; /*Initialize the statistical counters*/ num_custs_delayed = 0; total_of_delays = 0.0; area_num_in_q = 0.0; area_server_status = 0.0; /* Initialize event list. Since no customers are present, the departure (service completion) event is eliminated from consideration*/ time_next_event[1] = sim_time + expon(mean_interarrival); time_next_event[2] = 1.0e+30; }

24. Implementation using C void timing(void) /*Timing function*/ { inti; float min_time_next_event = 1.0e+29; next_event_type = 0; /*Determine the event type of the next event to occur*/ for (i = 1; i <= num_events; ++i) if (time_next_event[i] < min_time_next_event) { min_time_next_event = time_next_event[i]; next_event_type = i; } /*Check to see whether the event list is empty*/ if (next_event_type == 0) /*The event list is empty, so stop the simulation*/ fprintf(outfile, “\nEvent list empty at time %f”, sim_time); exit(1); } /*The event list is not empty, so advance the simulation clock*/ sim_time = min_time_next_event; }

25. Implementation using C void arrive(void) /*Arrival event function*/ { float delay; /*Schedule next arrival*/ time_next_event[1] = sim_time + expon(mean_interarrival); /*Check to see whether the server is busy*/ if (server_status == BUSY) { /*Server is busy, so increment number of customers in queue*/ { ++num_in_q; /*Check for overflow*/ if (num_in_q > Q_LIMIT) { /*The queue has overflowed, so stop the simulation*/ fprintf(outfile, “\nOverflow of the array time arrival at”); fprintf(outfile, “ time %f”, sim_time); exit(2); } /*There is still room in the queue, so store the time of arrival of the arriving customer at the new end of time_arrival*/ time_arrival[num_in_q] = sim_time; } else { /*Server is idle, so arriving customer has a delay of zero. (The following two statements are for program clarity and do not affect the results of the simulation*/ delay = 0; total_of_delays += delay; /*Increment the number of customers delayed, and make server busy*/ ++num_custs_delayed; server_status = BUSY; /*Schedule a departure (service completion)*/ time_next_event[2] = sim_time + expon(mean_service); } }

26. Implementation using C void depart(void) /*Departure event function*/ { inti; float delay; /*Check to see whether the queue is empty*/ if (num_in_q == 0) { /*The queue is empty so make the server idle and eliminate departure (service completion) event from consideration*/ server_status = IDLE; time_next_event[2] = 1.0e+30; } else { /*The queue is nonempty, so decrement the number of customers in the queue*/ --num_in_q; /*Compute the delay of the customer who is beginning service and update the total delay accumulator*/ delay = sim_time – time_arrival[1]; total_of_delays += delay; /*Increment the number of customers delayed, and schedule departure*/ ++num_custs_delayed; time_next_event[2] = sim_time + expon(mean_service); /*Move each customer in queue (if any) up one placed*/ for (i = 1, i <= num_in_q, ++i) time_arrival[i] = time_arrival[i + 1]; } }

27. Implementation using C void report(void) /*Report generator function*/ { /*Compute and write estimates of the desired measures of performance*/ fprintf(outfile, “\nAverage delay in queue%11.3f minutes\n\n”, total_of_delays / num_custs_delayed); fprintf(outfile, “Average number in queue%10.3f\n\n”, area_num_in_q / sim_time); fprintf(outfile, “Server utilization%15.3f\n\n”, area_server_status / sim_time); fprintf(outfile, “Time simulation ended%12.3f minutes”, sim_time); } void update_time_avg_stats(void) /*Update area accumulators for time-average statistics*/ { float time_since_last_event; /*Compute time since last event, and update last event time marker*/ time_since_last_event = sim_time – time_last_event; time_last_event = sim_time; /*Update area under numer_in_qeue function*/ area_num_in_q += num_in_q *time_since_last_event; /*Update area under server-busy indicator function*/ area_server_status += server_status *time_since_last_event; } float expon(float mean) /*Exponential variate generation function*/ { /*Return an exponential random variate with mean “mean”*/ return –mean * log(lcgrand(1)); }

28. Simulation Output Single-server queuing system Mean interarrival time 1.000 minutes Mean service time 0.500 minutes Number of customers 1000 Average delay in queue 0.430 minutes Average number in queue 0.418 Server utilization 0.460 Time simulation ended 1027.915 minutes

29. Simulation of Inventory Systems

30. The Simple Inventory System • The system is called an (M, N) system • The maximum inventory level is (M units) • The review time is N (weeks, days, hrs, etc.) I M Q1 Q3 Amount in Inventory Q2 T Time N N N

31. Basic Concepts of Inventory Systems • The lead time is the length of time between placement and receipt of an order • Demands are usually not known with certainty, they are modeled with a probability distribution. In the figure, demand as shown as being uniform over the review period. In practice, demands fluctuate over time • In the second cycle, demand drops below zero. These units are back-ordered, i.e. these units will be satisfied first when a new order arrives • In some inventory systems, a portion of sales may be lost rather than back-ordered when inventory runs out

32. Basic Concepts of Inventory Systems • The total cost (or total profit) of an inventory system is the measure of performance • Inventory system parameters: • Note: Lead time may be random if it cannot be controlled

33. System costs • Ordering cost = K + iZ K = setup cost, Z = S – I if I < s 0 if I ≥ s • Holding cost h (of every item held in inventory per month) • Backlog cost Pi (cost of every extra record of every item in the backlog)

34. Example: Refrigerator Inventory Problem • A company selling refrigerators maintains inventory by conducting a review every 1 month and making a decision about the quantity to order. The factory is trying to evaluate the following ordering policies • The demand, i.e. the number of refrigerators purchased by customers, is randomly distributed as shown in Table (1). The times between demands are exponentially distributed with mean 0.1 month After the company places an order to replenish their supply, the lead time uniformly distributed between 0.5 and 1 month • Assuming the ending inventory of the last review cycle is 60 refrigerators, K = 32.0, i = 3.0, h = 1.0, pi= 5.0: construct a simulation table and calculate the measures of performance Table (1)

35. Implementation using C Demand Evaluate Order Arrival End Simulation

36. Implementation using C Order Arrival Event Demand event Generate the size of this demand Increment the inventory level by the amount previously ordered Eliminate order-arrival event from consideration Decrement the inventory level by the demand size Schedule the next demand event Return Order arrival routine Return Demand routine

37. Implementation using C Inventory evaluation event Yes No Is I(t) < s? Determine amount to be ordered [S – I(t)] Incur ordering costs and gather statistics Schedule order arrival event for this order Schedule the next inventory-evaluation event Inventory evaluation routine Return

38. Implementation using C Update time avg stats Was I(t) during the previous interval negative, zero, or positive? Negative Positive Update backlog costs Update holding costs Zero Return Updating statistical accumulators

39. Implementation using C /*External definitions for inventory system*/ #include <stdio.h> #include <math.h> #include “lcgrand.h” /*header file for RNG int amount, bigs, initial_inv_level, inv_level, next_event_type, num_events, num_months, num_values_demand, smalls; float area_holding, area_storage, holding_cost, incremental_cost, maxlag, mean_interdemand, minlag, prob_distrib_demand[26], setup_cost, shortage_cost, sim_time, time_last_event, time_next_event[5], total_ordering_cost; FILE *infile, *outfile void initialize(void); void timing(void); void order_arrival(void); void demand(void); void evaluate(void); void report(void); void update_time_avg_stats(void); float expon(float mean); Intrandom_integer(float prob_distrib [ ]); float uniform(float a, float b);

40. Implementation using C main( ) /*Main function*/ { inti, num_policies; /*Open input and output files*/ infile = fopen(“inv.in”, “r”); outfile = fopen(“inv.out”, “w”); /*Specify the number of events for the timing function*/ num_events = 4; /*Read input parameters*/ fscanf(infile, “%d %d %d %d %f %f %f %f %f %f %f”, &initial_inv_level, &num_months, &num_policies, &num_values_demand, &mean_interdemand, &setup_cost, &incremental cost, &holding_cost, &shortage_cost, &minlag, &maxlag); for (i = 1; i <= num_values_demand; ++i) fscanf(infile, “%f”, &prob_distrib_demand[i]); /*Write report heading and input parameters*/ fprintf(outfile, “single-product inventory system\n\n”); fprintf(outfile, “Initial inventory level%24d items\n\n”, initial_inv_level); fprintf(outfile, “Number of demand sizes%25d\n\n”, num_values_demand); fprintf(outfile, “Distribution function of demand sizes ”); for (i = 1; i <= num_values_demand; ++i) fprintf(outfile, “%8.3f”, prob_distrib_de,amd[i]); fprintf(outfile, “\n\nMeaninterdemand time%26.2f\n\n”, mean_interdemand); fprintf(outfile, “Delivery lag range%29.2f to%10.2f months\n\n”, minlag, maxlag); fprintf(outfile, “Length of the simulation%23d months\n\n”, num_months); fprintf(outfile, “K = %6.1f i = %6.1f h = %6.1f pi = %6.1f\n\n”, setup_cost, incremental_cost, holding_cost, shortage_cost); fprintf(outfile, “Number of policies%29d\n\n”, num_policies); fprintf(outfile, “ Average Average”); fprintf(outfile, “ Average Average\n”); fprintf(outfile, “ Policy total_costordering_cost”); fprintf(outfile, “ holding cost shortage cost”);

41. Implementation using C – main() cont’d /*Run the simulation varying the inventory policy*/ for (i = 1, i =< num_policies, ++i) { /*Read the inventory policy, and initialize the simulation*/ fscanf(infile, “%d %d”, &smalls, &bigs); initialize(); /*Run the simulation until it terminates after an end-simulation event (type 3) occurs*/ do { /*Determine the next event*/ timing(); /*Update time-average statistical accumulators*/ update_time_avg_stats(); /*Invoke the appropriate event function*/ switch(next_event_type) { case 1: order_arrival(); break; case 2: demand(); break; case 4: evaluate(); break; case 3: report(); break; } /*If the event just executes was not the end-simulation event (type 3), continue simulating. Otherwise, end the simulation for the current (s, S) pair and go on to the next pair (if any)*/ } while (next_event_type != 3); } /*End the simulation*/ Fclose(infile); Fclose(outfile); Return 0; }

42. Implementation using C void initialize(void) /*Initialize function*/ { /*Initialize the simulation clock*/ sim_time = 0.0; /*Initialize the state variables*/ inv_level = initial_inv_level; time_last_event = 0.0; /*Initialize the statiastical counters*/ total_ordering_cost = 0.0; area_holding = 0.0; area_shortage = 0.0; /*Initialize the event list. Since no order is outstanding, the order-arrival event is eliminated from consideration*/ time_next_event[1] = 1.0e+30; time_next_event[2] = sim_time + expon(mean_interdemand); time_next_event[3] = num_months; time_next_event[4] = 0.0; }

43. Implementation using C void order_arrival(void) /*Order arrival event function*/ { /*Increment the inventory level by the amount ordered*/ inv_level += amount; /*Since no order is now outstanding, eliminate the order-arrival event from consideration*/ time_next_event[1] = 1.0e+30; } void demand(void) /*Demand event function*/ { /*Decrement the inventory level by a generated demand size inv_level -= random_integer(prob_distrib_demand); /*Schedule the time of the next demand*/ time_next_event[2] = sim_time + expon(mean_interdemand);

44. Implementation using C void evaluate(void) /*Inventory-evaluation event function*/ { /*check whether the inventory level is less than smalls*/ if (inv_level < smalls) { /*The inventory level is less than smalls, so place an order for the appropriate amount*/ amount = bigs – inv_level; total_ordering_cost = setup_cost + incremental_cost * amount; /*Schedule the arrival of the order*/ time_next_event[1] = sim_time + uniform(minlag, maxlag); } /*Regardless of the place-order decision, schedule the next inventory evaluation*/ time_next_event[4] = sim_time + 1.0; } Void report(void) /*Report generator function*/ { /*Compute and write estimates pf desired measures of performance*/ float avg_holding_cost, avg_ordering_cost, avg_shortage_cost; avg_ordering_cost = total_ordering_cost / num_months; avg_holding_cost = holding_cost * area_holding / num_months; avg_shortage_cost = shortage_cost * area_shortage / num_months; fprintf(outfile, “\n\n(%3d, %3d)%15.2f%15.2f%15.2f%15.2f”, smalls, bigs, avg_ordering_cost + aveg_holding_cost + avg_shortage_cost, avg_ordering_cost, avg_holding_cost, avg_shortage_cost); }

45. Implementation using C void update_time_avg_stats(void) /*Update area accumulators for time-average statistics*/ { float time_since_last_event; /*Compute time since last event, and update last-event-time marker*/ time_since_last_event = sim_time – time_last_event; time_last_event = sim_time; /*Determine the status of the inventory level during the previous interval. If the inventory level during the previous interval was negative, update area_shortage. If it was positive, update area_holding. If it was zero, no update is needed */ if (inv_level < 0) area_shortage -= inv_level * time_since_last_event; else if (inv_level >0) area_holding += inv_level *time_since_last_event; } intrandom_integer(float prob_distrib[ ]) /* Random integer generator function*/ { int I; float u; /*Generate a U(0, 1) random variate*/ u = lcgrand(1); /*Return a random integer in accordance with the (cumulative) distribution function prob_distrib*/ for (I = 1; u >= prob_distrib[i]; ++i) ; return i; } float uniform(float a, float b) /*Uniform random variate function*/ { /*Return a U(a, b) random variate*/ return a + lcgrand(1) * (b – a); }

46. Simulation Output Single-product inventory system Initial inventory level 60 times Number of demand size 4 Distribution function of demand sizes 0.167 0.500 0.833 1.00 Mean interdemand time 0.10 months Delivery lag range 0.50 to 1.00 months Length of the simulation 120 months K = 32.0 i = 3.0 h = 1.0 pi = 5.0 Number of policies 9