1 / 76

Chapter 10 Sample path sensitivity analysis

Chapter 10 Sample path sensitivity analysis. Learning objectives : Understand the mathematics behind discrete event simulation Learn how to derive gradient information during a simulation Able to verify the correctness of sample gradients Textbook :

yair
Download Presentation

Chapter 10 Sample path sensitivity analysis

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Chapter 10 Sample path sensitivity analysis • Learning objectives : • Understand the mathematics behind discrete event simulation • Learn how to derive gradient information during a simulation • Able to verify the correctness of sample gradients • Textbook : • C. Cassandras and S. Lafortune, Introduction to Discrete Event Systems, Springer, 2007

  2. Plan • Introduction • Perturbation analysis by example • GSMP - Generalised Semi-Markovian Process • Some mathematical basis of discrete event simulation • Perturbation analysis of GSMP • Correctness of gradient estimators • Applications to manufacturing systems • Sensitivity estimation revisited • Extensions of perturbation analysis

  3. Introduction

  4. Sensitivity measures • Two types of sensitivitymeasures: • Derivatives : d J(q)/dq • Finitedifference : J(q+D) - J(q) • where J(q) is the performance funtion of a system atparameterq • Sensitivitymeasures are for • determination of the direction of change to take • optimization of a system design • performing "what-if" analysis

  5. Sensitivity measures • Difficulties of sensitivityanalysis: • Unknown performance function J(q) • Estimation by simulation in most cases • Main ideas of perturbation analysis • Prediction of the behavior of perturbed system q+Dalong the simulation of the system q • Main steps of perturbation analysis • Perturbation generation: changes in quantitiesdirectlylinked to q • Perturbation propagation: changes in otherquantitiesgenerated by perturbation generation

  6. Stochastic approximation • Stochastic approximation (gradient basedminimization) : • x(n+1) = x(n) - a(n)g(n) • where • x(n) : value of x atiteration n • a(n) > 0 : stepsize • g(n) : gradient estimation of a performance function f(x) • Round-Robin Theorem: Assume that f(x) iscontinuouslydifferentiable and thereexists a single x* suchthatdf(x*)/dx = 0. If • E[g(n)] = df(x)/dx • and • limn a(n)= 0, Sn a(n) = ∞ • thenx(n) converges withprobability 1 to x*.

  7. Perturbation analysis by example

  8. An inspection machine • Products arrive randomly to an inspection machine • Each product requires a random number of tests • The inspection of a production requires a setup of time g • Each test takes q time units • Products are inspected in FIFO order • The buffers are of unlimited capacity • The system is empty initially • Performance measure = T(q,g), average time at the inspection station.

  9. A G/G/1 queueing model Service time X = g + Lq Inter-arrival time A

  10. Sample path Definition : A SamplePathis the sequence (event, state, time) of a simulation or experimentation. Example : (Start, 0, 0), (arrival, 1, A1), (arrival, 2, A1+A2), (arrival, 3, A1+A2+A3), (departure, 2, A1+X1), … Property : The samplepathcontains all information of the simulation.

  11. An inspection machine • Each simulation run or a samplepathischaracterized by : • Li : number of tests of product i • Xi = g + Liq : service time of product i • Ai : inter-arrival time betweenproduct i and product i-1 • ti : time at the inspection station of product i (alsocalled system time) • Attention ti ≠ Xi

  12. An inspection machine The best estimator of the average time at the inspection station: Under some regularity conditions Goal of perturbation analysis Evaluate ∂T(q, g)/∂ g and ∂T(q, g)/∂ q

  13. Event-based system dynamics {e1, e2, ..., en, ...} : sequence of events with en{a, d} tn : occurrence time of en with tn = 0 qn : number of parts in the system at time tn+ with qn = 0 Dn : time in state qn En : set of active events at tn+ ren : remaining time till occurrence of event e at time tn+

  14. Event-based system dynamics en : sequence of events tn : occurrence time of en qn : number of parts at time tn+ Dn : time in state qn En : set of active events at tn+ ren : remaining time of event e Time to nextevent :Dn= Min (ren , eEn) Nextevent : en+1 = Argmin (ren , eEn) New state qn+1 = qn +1 if en+1 = a qn+1 = qn -1 if en+1 = d Update eventclocks ren+1 = ren - Dn En+1 = En -{en+1} Create new event Case e = a and lessthan N arrivals En+1 = En+1 +{a} ran+1 = An+1 Case e = d and qn+1 > 0 AND Case e = a and qn = 0 En+1 = En+1 +{d} ran+1 = g+Ln+1q

  15. Event-based system dynamics en : sequence of events tn : occurrence time of en qn : number of parts at time tn+ Dn : time in state qn En : set of active events at tn+ ren : remaining time of event e Performance evaluation : Sn = total time of all parts in the system up to event n Sn+1 = Sn + qnDn T(q,g) = SV /N where eVis the departureevent of N-th part

  16. Event-based system dynamics Simulation algorithm 1/ Initialization : q = 0, t = 0, E = {a}, generater.v. A, ra = A Na = Nd = S = 0 2/ Nextevent: e= Argmin (re, eE), D= re Statistics : Na = Na+ 1(e=a), Nd = Nd+ 1(e=d), S = S + qD 3/ State update : q = q + 1(e = a)-1(e= d) 4/ Update eventclocks: re = re - D, E = E -{e} 5/ Generate new events If e=a  Na < N, E = E +{a}, generater.v. A, ra = A If (e=a  q= 1) or (e=d  Nd < N), E = E +{d}, generater.v. L, rd = g+Lq

  17. Perturbation analysis • Typical question (What-If) : • Whatwouldbe the average system time • if the experimentationwereperformedwithparameterq+Dq? • Definitions : • Nominal system : the system withparameterq of the simulation • Perturbed system : the system withparameterq+Dq of the What-If question? • Nominal samplepath: samplepath of the nominal system • Perturbedsamplepath: samplepath of the perturbed system (preferablyunder the samerandom condition as the nominal samplepath

  18. Perturbation analysis • Typical perturbation analysisanswer: • Perturbation generation • The perturbation of parameterqgeneratesdirectly a perturbation ∆Xi = Li∆q to eachproduct service time (recallthat Xi = g + Liq) • Perturbation propagation • The perturbation of the service time of a productwillbepropagated to otherproducts • Fundamental PA condition : the sequence of eventsdoes not change for smallenough perturbation Dq

  19. Perturbation analysis

  20. Event-based perturbation analysis Nextevent : en+1 =e*= Argmin (ren , eEn) Time to nextevent : Dn= re*n , Dn/q= re*n /q New state qn+1 Update eventclocks ren+1 = ren - Dn , ren+1/q= ren+1 /q - Dn/q En+1 = En -{en+1} Create new event New arrival : En+1 = En+1 +{a}, ran+1 = An+1, ran+1/q= 0 New depart : En+1 = En+1 +{d}, rdn+1 = g + Ln+1q, rdn+1/q= Ln+1 Statistics : Sn+1 = Sn + qnDn,  Sn+1 /q=  Sn/q+ qnDn/q

  21. Event-based system dynamics Simulation algorithm 1/ Initialization : q = 0, t = 0, E = {a}, generater.v. A, ra = A Na = Nd = S = 0, grad_ra = 0 2/ Nextevent: e= Argmin (re, eE), D= re , grad_D = grad_re Statistics : Na = Na+ 1(e=a), Nd = Nd+ 1(e=d), S = S + qD, grad_S = grad_S + q*grad_D 3/ State update : q = q + 1(e = a)-1(e= d) 4/ Update eventclocks: E = E -{e}, re = re - D, grad_re= grad_re- grad_D 5/ Generate new events If e=a  Na < N, E = E +{a}, generate A, ra = A, grad_ra = 0 If (e=a  q= 1) or (e=d  Nd < N), E = E +{d}, generate L, rd = g+Lq, grad_rd = L 6/ If Nd < N, go to 2; otherwise, T = S/N, grad_T = grad_S /N, END

  22. Busy period • Definition : • A busyperiod (BP) is a period of time duringwhich the inspection station isneveridle. • For a product i of the 1stbusyperiod (BP1): • arrival date : A1+A2 + ...+Ai • departure time: A1+ X1+X2 + ...+Xi • system time : ti = Xi + (X1+X2 + ...+Xi-1) - (A2 + ...+Ai) • Total system time of products of BP1

  23. Busy period • EachbusyperiodBPmisinitiated by the arrival of a productthatfinds an empty system • Let km +1be the productthatinitiatesbusyperiodBPm • System time of the i-th product of BPm • Total system time of products of BPm • where nmis the number of productsinspected in BPm

  24. Busy period • Total system time of products of BPm • where M is the total number of busyperiodsobservedduring the simulation • M is a random variable thatdepends on the simulation realization.

  25. Busy period-based Perturbation analysis • Busyperiodsremain the samednder the fundamental condition. • System time in the perturbed system • As a result,

  26. Busy period-based Perturbation analysis • To summarize, for Xi = g + Liq,

  27. Busy period-based Perturbation analysis Perturbation analysisalgorithmduring the simulation 0) Initialization: gradXq = 0, gradtq = 0, gradTq = 0 gradXg = 0, gradtg = 0, gradTg = 0 1) During the simulation : If the end of service of productwith L tests, 1.1) Perturbation generation: gradXq = L and gradXg = 1 1.2) Perturbation propagation (system time) gradtq = gradtq+ gradXq (Total system time) gradTq = gradTq + gradtq gradtq = gradtg+ gradXg gradTg = gradTg + gradtg 1.3) If the system isempty, then gradtq = 0 and gradtg = 0 2) At the end of the simulation, ∂T/∂q = gradTq/N and ∂T/∂q = gradTq/N

  28. Validation of the results Unbiasedness Consistency

  29. GSMP - Generalised Semi-Markovian Process A general framework for representation of Discrete event systems

  30. Basic concepts • States : a state is a possible system configuration • Ex : number of products to inspect • Events : the system state changes only upon the occurrence of events • Ex : arrival and end of inspection of a product • Transition probabilities : which rule the change of the state when an event occurs • Ex : routing probabilities among multiple alternatives • Clock : which indicate the remaining time to occurrence of events. • A clock is set according to some random variable when it is activated • It then decrements as long as it is active. • An event occurs when its clock reaches zero. • Attention : Remaining processing times of products are associated with clocks of related events and not in the definition of the state

  31. GSMP • A GSMP is defined by : • S : set of states • E : set of events • E(s) : set of events that are possible at state s • p(s'; s, e) : probability of jumping to state s' when event e occurs at state s • Fe : probability distribution of a new clock of event e • Example (Station d'inspection) : • S = {0, 1, …} : number of products waiting for inspection • E = {a, d} with a = product arrival, d = product departure • E(s) = {a, d} if s > 0 and E(0) = {a} • p(s+1; s, a) = 1, p(s-1; s, d) d 1 • Fa = distribution of inter-arrival time, Fd = distribution of service time

  32. GSMP Non-Interruption (NI) condition : " s, s' Œ S and e Œ E(s), p(s'; s, e) > 0  E(s') E(s) - {e} The NI condition implies that, once activated, an event remain active till its occurrence.

  33. Simulation algorithm of an GSMP 1. Initialization : • Number of events : n = 0 • Set initial state : S0 • Set clocks of activated events : C0(e) = X(e, 1), "e  Œ E(S0) • Set event counter : N(e, 0) = 0, " e Œ E 2. Determine the next event : 3. Sojourn time at Sn : tn+1 = Cn(en+1) 4. Update the event counters :

  34. Simulation algorithm of an GSMP • 5. Update the system state : • Sn+1 = F(Sn , en+1 , U(en+1, N(en+1, n+1))) • where • Fis a functionsuchthat P{F(s, e, U) = s'} = p(s'; s, e) • U is a random variable withuniform distribution on [0,1]. • Update existingclocks : • Cn+1(e) = Cn+1(e) - tn+1, "e  E(sn) -{en+1} • Generate new clocks : • Cn+1(e) = X(e, N(e, n) +1), "e  E(sn+1) - (E(sn) -{en+1}) • 8. If the end of simulation, stop. Otherwise, n:= n+1 and go to step 2

  35. Simulation algorithm of an GSMP Routing mechanisms • Each active eventeisassociated a uniformrandom variable U(e), calledroutingindicator • Step 5 of the GSMP simulation algorithmisdecomposeintotwosteps: • 5.1. Generation of the routingindicator : U(en+1,N(en+1, n+1)) • 5.2. Determination of the new state : • Sn+1 = F(Sn , en+1 , U(en+1,N(en+1, n+1)) ) • whereFissuchthat P{F(s, e, U) = s'} = p(s'; s, e). • Example of F : • Let S = { s1, s2, … , sm }. Sn+1 := sisuchthat :

  36. Simulation algorithm of an GSMP Performance measures • Let • Zt : be the state of the system at time t and • f(Zt) : the costassociatedwith state Zt. • The following performance measures are considered • where T is a given constant • where N is a giveninteger and tNis occurrence time of N-th event, i.e. tN = t1 + t2 + ... + tN ; • where T(a,k) is the time of k-th occurrence of event a .

  37. Some mathematical basis of discrete event simulation

  38. Random Number Generation The linear congruential technique • Goal: • Generate samples of U[0, 1] distribution • The linear congruential technique : • Xk+1 = (aXk + c) mod m, k = 0, 1, ... • where • a = multiplier, c = increment, m = modulus • X0 = seed • Quality strongly depends on the selection of parameters a, c, m. • The "pseudo" random samples have cycle of length at most m. • m = 2n , as large as possible with 231 for 32-bit computer, is often used to maximize the cycle and to ease to the modula arithmetic

  39. Random variate generation The inverse transform technique • Goal: • Generate samples of probability distribution F(x) • The inverse tranform technique : • Generate a sample U of distribution U[0, 1] • Transform U into X with X = F-1(U) • Online proof • Examples : • EXP(l) with F = 1 - exp(-lx) • Weibull with F = 1 - exp (-(l(x-g))b) • Drawback: • The computation of the inverse transform can be time consuming.

  40. Random variate generation Acceptance-rejection technique • Majorizingfunction: • g(x)suchthat g(x) ≥ f(x) • Densityfunctionrelated to g(x) • Acceptance-rejection technique : • Select a majorizingfunction g(x) and itsrelateddensityfunction h(x) • Generate a randomvariate Y withpdf h(x) • Generate a randomnumber U[0,1] independentlyfrom Y • If U ≤ f(Y)/g(Y), • then set X = Y. Otherwise, repeatstep 2. • Criteria for selection of majorizingfunction: • Ease of generation of distribution h(x) • Small rejection points, i.e. small surface between g(x) and f(x)

  41. Terminating simulation output analysis Point estimation • Observations: • X1, X2, ..., Xn : simulation observations of true but unknown performance q • Sample mean • Unbiasedness: • For iid observations X1, X2, ..., Xn, • Sample variance: strong law of large number

  42. Terminating simulation output analysis Interval estimation Central limit theorem tends to a standard normal distribution as n goes to infinity. (1-a)-percent confidence interval (why): • where • za/2is a defined on a standard normal distribution F(z). • s2isevaluated by sample variance • Student's t distribution isalsoused to avoidchecking how large n is for normal distribution approximation

  43. Nonterminating simulation output analysis Problem setting • Observations: • X1, X2, ... : simulation observations of a nonterminating simulation • Unknown steady state performance : • Time average: • Under ergodicity condition

  44. Nonterminating simulation output analysis Replications with deletions • I. • Simulation for a time horizon of m+r steps • Ignore data of the first r steps • Estimation with data of the remaining m steps • II. • Perform K replications of I • Use point estimation methods to evaluate the performance • Key questions : • How long should the warmup interval r be?

  45. Nonterminating simulation output analysis Batch means • Based on a single but long simulation run • Observed data X1, X2, ... are grouped into Batches • Assume batches of equal length m with a warmup period of length r • Batch mean Bj is estimated • Batch means are then used in point estimation

  46. Nonterminating simulation output analysis Regenerative simulation • Assumption: • There exists some regenerative state (regeneration point) after which the system is independent of past history • Regeneration cycles • R1, R2, ... : regeneration times • [Rj, Rj+1) : regenerative cycle or period • Cycle total (Lj) and cycle length (Nj) of cycle j • Steady state mean

  47. Perturbation analysis of GSMP

  48. Functional GSMP • System paramets : • Assume that the probability distribution of new clocksdepend on a continuousparameterq. • Let : • Fe(x; q) : distribution of new clocks of evente • Xq(e, n) : n-th clock of evente • Hypothesis (H1) : • "e and q, Fe(x; q) iscontinuous in x and Fe(x; q) = 0. • Hypothesis (H2) : • "e and q, Xq(e, n)iscontinuousdifferentiable in q. • Remark: • Under H1, the probabilitythattwoclocksreachzeroat the same time isnull

  49. Perturbation generation Remark : U = Fe(Xq(e, n); q) est une v.a. uniforme dans [0, 1] The inverse transform technique for generation of eventclocks X : Generate a random variable U uniformlydistributed on [0, 1] Xq(e, n) = Fe-1(U, q) Property : When the inverse tranform technique isused for generation of eventclocks, the followinginfinitesimal perturbations are generated :

  50. Perturbation generation Special cases : • scaleparameters : Fe(x; q) = Ge(x/q) dXq(e, n) / dq = Xq(e, n) / q Ex : a exponentiallydistributedrandom variable withmeanq • location parameters : Fe(x; q) = Ge(x - q) dXq(e, n) / dq = 1 Ex : X = g + L qwheregis a location parameter

More Related