Patent Publication Number: US-11023816-B2

Title: Disruption forecasting in complex schedules

Description:
This application is a continuation of Ser. No. 14/657,378, filed Mar. 13, 2015, the entire disclosure of which is incorporated herein by reference. 
    
    
     BACKGROUND 
     The present invention relates to system and method for forecasting a probabilistic distribution of time delays added to a scheduled start time and a scheduled end time of a task. Specific applications include forecasting probabilistic distributions of arrival times and departure times for an item such as a transport or a work piece in a process. 
     Many types of processes involve the arrival and departure of a certain type of item such as a train at a railroad station or a work piece at a machining station. While an arrival time and departure time may be scheduled or planned, in certain cases the item may arrive later than scheduled and thus depart later than scheduled. These delays can have ramifications in the associated process such as causing other delays that can propagate throughout the process. The cost of some types of delays may be minimal while other types of delays may cause multiple delays throughout the process and thus be very costly. 
     Disclosed is a method for forecasting time delays added to a scheduled start time and a scheduled end time of a task. The method includes: generating a stochastic model of the task and resources affecting the task, the stochastic model comprising a reactionary delay component and a root cause delay component, the reactionary component being a function of previous task end times and the root cause delay component being an independent random process at a specific time; calculating a probability distribution of time delays added to the scheduled start time as a combination of the reactionary delay component and the root cause delay component using the stochastic model to provide a probability distribution of start times; and calculating a probability distribution of time delays added to the scheduled end time as a combination of the reactionary delay component and the root cause delay component using the stochastic model to provide a probability distribution of end times; wherein the generating, the calculating a probability distribution of time delays added to the scheduled start time, and the calculating a probability distribution of time delays added to the scheduled end time are implemented by a processor. 
     Also disclosed is a system method for forecasting time delays added to a scheduled start time and a scheduled end time of a task. The system includes a processor configured to: generate a stochastic model of the task and resources affecting the task, the stochastic model comprising a reactionary delay component and a root cause delay component, the reactionary component being a function of previous task end times and the root cause delay component being an independent random process at a specific time; calculate a probability distribution of time delays added to the scheduled start time as a combination of the reactionary delay component and the root cause delay component using the stochastic model to provide a probability distribution of start times; calculate a probability distribution of time delays added to the scheduled end time as a combination of the reactionary delay component and the root cause delay component using the stochastic model to provide a probability distribution of end times; transmit a signal comprising the probability distribution of start times and the probability distribution of end times to a signal receiving device. The system further includes a signal receiving device configured to receive the signal comprising the probability distribution of start times and the probability distribution of end times. 
     Additional features and advantages are realized through the techniques of the present invention. Other embodiments and aspects of the invention are described in detail herein and are considered a part of the claimed invention. For a better understanding of the invention with the advantages and the features, refer to the description and to the drawings. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The subject matter which is regarded as the invention is particularly pointed out and distinctly claimed in the claims at the conclusion of the specification. The forgoing and other features, and advantages of the invention are apparent from the following detailed description taken in conjunction with the accompanying drawings in which: 
         FIG. 1  is a schematic diagram depicting aspects of a train arriving and later departing a train station according to an embodiment of the invention; 
         FIG. 2  is a schematic diagram depicting aspects of a work piece arriving at a computerized numerical control (CNC) machine for machining and then later departing the CNC machine after being machined according to an embodiment of the invention; 
         FIG. 3  is a schematic diagram of a computer processing system for implementing a method for disruption forecasting in complex schedules; 
         FIG. 4  depicts aspects of a Markov model with three states; 
         FIG. 5  depicts aspects of a Markov model with three states and emissions from the three states; 
         FIG. 6  depicts aspects of data inputs and calculated outputs using a process model and forecasting; 
         FIG. 7  is an overview of predictive outputs of the process model and forecasting in one embodiment; 
         FIG. 8  is an overview of state-estimation outputs of the embodiment; 
         FIG. 9  is an overview of process components of the embodiment for a departure D i  and arrival A i  of a single train movement i; 
         FIG. 10  presents method steps of the prediction of the departure D i  and arrival A i  of a single train movement i, using only a single resource az; and 
         FIG. 11  depicts aspects of a Gantt chart for a machining application. 
     
    
    
     DETAILED DESCRIPTION 
     Embodiments of systems and methods described herein relate to calculating a probabilistic delay distribution for scheduled resources such as transports or work pieces in non-limiting embodiments. Rather than forecasting a single number for a delay time, the systems and methods give an explicit probability distribution over possible delays. Consider the following example. Suppose a railway operations controller is evaluating a particular train scheduled in the near future. It might be the case that a departure delay under 40 minutes for this train will incur relatively little cost to the railway, but a delay exceeding 40 minutes could be very costly due to missed passenger connections, missed departure slots at stations, crew overtime, etc. If the operations controller knows that there is a 35% chance of a delay over 40 minutes, he or she can weigh the potential cost and risk of not acting against the cost of some mitigating action, like swapping the rail equipment assigned to this train. A single point forecast that does not account for uncertainty, such as “the expected departure delay is 26 minutes” is far less useful for guiding this decision. 
     The systems and methods disclosed herein identify and characterize the behavior of “delay states” from historical data. These states capture the time-dependent nature of certain kinds of delays, which are useful for accurately forecasting future delays. As a railway example of a delay state, suppose that many trains are experiencing departure delays from Boston on a winter morning. A person would reason that “there is probably bad weather in Boston right now and that will probably continue to affect departures for at least the next several hours.” The systems and methods disclosed herein identify patterns indicating these types of time-dependent sources of delay from data, and characterize the magnitude of delays experienced in each state together with the usual duration of time spent in each state. 
     It can be appreciated that the systems and methods disclosed herein may be applied to various processes. One process for teaching purposes is a transport process such as a railway process for transporting passengers.  FIG. 1  is a schematic diagram depicting aspects of a train  90  arriving and later departing a train station  91  according to an embodiment of the invention. The train  90  has an arrival time, a task time such as unloading and loading passengers, and a departure time. It can be appreciated that the train  90  may be representative of a transport in other types of transportation systems. Another example of a process is a manufacturing process such as machining a work piece.  FIG. 2  is a schematic diagram depicting aspects of a work piece  92  arriving at a computerized numerical control (CNC) machine  93  for machining and then later departing the CNC machine after being machined according to an embodiment of the invention. 
       FIG. 3  depicts a block diagram of a computer system for implementing the teachings disclosed herein according to an embodiment. Referring now to  FIG. 3 , a block diagram of a computer system  10  suitable for providing communication over cross-coupled links between independently managed compute and storage networks according to exemplary embodiments is shown. Computer system  10  is only one example of a computer system and is not intended to suggest any limitation as to the scope of use or functionality of embodiments described herein. Regardless, computer system  10  is capable of being implemented and/or performing any of the functionality set forth hereinabove. 
     Computer system  10  is operational with numerous other general purpose or special purpose computing system environments or configurations. Examples of well-known computing systems, environments, and/or configurations that may be suitable for use with computer system  10  include, but are not limited to, personal computer systems, server computer systems, thin clients, thick clients, cellular telephones, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, set top boxes, programmable consumer electronics, network PCs, minicomputer systems, mainframe computer systems, and distributed cloud computing environments that include any of the above systems or devices, and the like. 
     Computer system  10  may be described in the general context of computer system-executable instructions, such as program modules, being executed by the computer system  10 . Generally, program modules may include routines, programs, objects, components, logic, data structures, and so on that perform particular tasks or implement particular abstract data types. Computer system  10  may be practiced in distributed cloud computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer system storage media including memory storage devices. 
     As shown in  FIG. 3 , computer system  10  is shown in the form of a general-purpose computing device, also referred to as a processing device. The components of computer system may include, but are not limited to, one or more processors or processing units  16 , a system memory  28 , and a bus  18  that couples various system components including system memory  28  to processor  16 . 
     Bus  18  represents one or more of any of several types of bus structures, including a memory bus or memory controller, a peripheral bus, an accelerated graphics port, and a processor or local bus using any of a variety of bus architectures. By way of example, and not limitation, such architectures include Industry Standard Architecture (ISA) bus, Micro Channel Architecture (MCA) bus, Enhanced ISA (EISA) bus, Video Electronics Standards Association (VESA) local bus, and Peripheral Component Interconnects (PCI) bus. 
     Computer system  10  may include a variety of computer system readable media. Such media may be any available media that is accessible by computer system/server  10 , and it includes both volatile and non-volatile media, removable and non-removable media. 
     System memory  28  can include computer system readable media in the form of volatile memory, such as random access memory (RAM)  30  and/or cache memory  32 . Computer system  10  may further include other removable/non-removable, volatile/non-volatile computer system storage media. By way of example only, storage system  34  can be provided for reading from and writing to a non-removable, non-volatile magnetic media (not shown and typically called a “hard drive”). Although not shown, a magnetic disk drive for reading from and writing to a removable, non-volatile magnetic disk (e.g., a “floppy disk”), and an optical disk drive for reading from or writing to a removable, non-volatile optical disk such as a CD-ROM, DVD-ROM or other optical media can be provided. In such instances, each can be connected to bus  18  by one or more data media interfaces. As will be further depicted and described below, memory  28  may include at least one program product having a set (e.g., at least one) of program modules that are configured to carry out the functions of embodiments of the disclosure. 
     Program/utility  40 , having a set (at least one) of program modules  42 , may be stored in memory  28  by way of example, and not limitation, as well as an operating system, one or more application programs, other program modules, and program data. Each of the operating system, one or more application programs, other program modules, and program data or some combination thereof, may include an implementation of a networking environment. Program modules  42  generally carry out the functions and/or methodologies of embodiments of the invention as described herein. 
     Computer system  10  may also communicate with one or more external devices  14  such as a keyboard, a pointing device, a display  24 , etc.; one or more devices that enable a user to interact with computer system/server  10 ; and/or any devices (e.g., network card, modem, etc.) that enable computer system/server  10  to communicate with one or more other computing devices. Such communication can occur via Input/Output (I/O) interfaces  22 . Still yet, computer system  10  can communicate with one or more networks such as a local area network (LAN), a general wide area network (WAN), and/or a public network (e.g., the Internet) via network adapter  20 . As depicted, network adapter  20  communicates with the other components of computer system  10  via bus  18 . It should be understood that although not shown, other hardware and/or software components could be used in conjunction with computer system  10 . Examples include, but are not limited to: microcode, device drivers, redundant processing units, external disk drive arrays, RAID systems, tape drives, and data archival storage systems, etc. 
     The systems and methods disclosed herein cover two use cases with two distinct input-output combinations, but having a common forecasting algorithm as discussed further below. In both cases, historical schedule and delay data, planned schedule and real-time delay data are input into a stochastic model. The forecasting algorithm then operates on the stochastic model. 
     In the next three paragraphs, a Markov model and a hidden Markov model for use herein are discussed in general. Markov models are one of the most basic stochastic models, which are fully determined by a set of states X={X i } transition probabilities p ij  from state X i  to state X j . The transition probabilities can be seen as a matrix p=(p ij ) or a probability distribution of the next state, conditional on the current state.  FIG. 4  illustrates one example of a Markov model. 
     Notice that the states are often not observable. Instead, “emissions” are observed within a set Y={y j }, where the emission probabilities r ij  of producing y j  in state x i  can be seen as a matrix r=(r ij ) or a probability distribution of the emission of y j , conditional on the current state X i .  FIG. 5  illustrates one example of a Markov model with three states and emissions. 
     The problem of estimating the most likely X, p, and r, given y, is known as the problem of training a Hidden Markov Model (HMM). Once a HMM is trained it becomes just a Markov model. A hidden Markov model (HMM) is a statistical model, which assumes the system is a Markov process with unobserved (hidden) states X, observations y, state transition probabilities p, and output probabilities r. 
     Training of an HMM model involves, given a history of observations y, estimating the most likely state transition probabilities p and output probabilities r. Application of the HMM model involves considering the state transition probabilities p and output probabilities r learned previously, using a history of observations y, to estimate the next state and observation. In a number of fields, including information theory, signal processing, and speech recognition, algorithms have been developed to compute the local maximum likelihood solution of the problem of training the HMM model in polynomial time. In particular, reference may be had to the algorithms of Viterbi, Baum-Welch, and Baldi-Chauvin. 
     The application of an HMM takes the input of: N—number of states; T—number of observations; y t =1 . . . T—observation at time t; θ i =1 . . . N—emission parameter associated with state I; φ i =1 . . . N, j=1 . . . N—probability of transition from state i to state j; x t =1 . . . T—state at time t; F(y|θ)—probability distribution of an observation, parametrized on θ and assuming x T +1˜Categorical(φ xT ), y T−1 ˜F(θ xT ) and produces: F(x T+1 |y, x)—probability distribution of the next state, parametrized on the above; and F(y T+1 |y, x)—probability distribution of the next observation, parametrized on the above, which can be done in time linear in the dimensions of the output, for a samples-based representation of the distributions. Notice that Categorical refers to a discrete random variable, i.e. a variable with a finite support set. 
     Iterated application of the HMM model involves considering the state transition probabilities p and output probabilities r learned previously, using a history of observations y, to estimate the next k states and k observations, in a recursive fashion. Notice that by iterated application of the HMM model, one can estimate the states and emissions arbitrarily far in the future. 
     The disclosure is discussed in a general form, before two specific examples in scheduling and railway operations are introduced. 
     Consider a problem where there is a set of tasks and a set of resources. A partial order O on the set of tasks, in which the tasks are to be performed, is known. Formally, a partial order is a relation, which is reflexive, antisymmetric, and transitive. The start time s i  when the task i should be performed is also known. One task can possibly include of multiple sub-tasks, which have to take place in a given order, one after another. Each task i can hence be associated with one start time s i  and multiple further times corresponding to the start-times of the sub-tasks. Multiple bipartite graphs may be used to capture the resource requirements required for tasks and sub-tasks. In each, one part of the vertices represents the tasks and one part of the vertices represents the resources, and there are no edges between two resources. Each task requires all resources it is adjacent to in one of the bipartite graphs, but the choice of the graph to use can vary task to task. 
     In the setting above, it is assumed that there is a state X r,t  of each resource r, which changes over time t, with the changes governed by a random process, but that the state is not directly observable. The random processes may or may not be independent across the resources. It is not assumed that the range is known, from which the values X r,t  are drawn, but it is assumed that its cardinality is known. 
     In the setting above, once a task i starts, one can obtain one observation of the actual start time, denoted a i , and multivariate, mixed-integer “additional information” c i , which are related to the state X r,ai  of all resources r required by task i at time a i . That is, the “additional information” may describe multiple categorical phenomena, which are encoded as integers, and multiple phenomena, which are described by real numbers. For example, the state of one resource may be, in part, described by the weather at the location of the resource, which can be described by the wind speed (real variable) and the information whether the sun is visible (binary, i.e. integer, variable). 
     In the setting above, at time t, one can describe the actual start times a i  and the additional information c i  for all the tasks i that have started by time t, as the “delay data”. 
     In the setting above, it is desired to forecast the actual start time of task i as a random variable A i  using a two-factor model, where A i  is a convolution of the two factors. One factor is the “reactionary delay component”, which is a function of the actual start times A j  for all (i, j) in the partial order O. The other factor is the “the root cause delay component,” which is based on the unobservable states of the resources required by the task in one of the bipartite graphs. Notice that at time t, one can use the observations of actual start times a j  for all (i, j) in the partial order O, where task j has started by time t, and an estimate of A j  for all (i, j) in the partial order O, where task j has not started by time t. Notice that at any time t, one may have to use an estimate of the state X r,t  of all resources r used by task i in one of the bipartite graphs, in the order to estimate the root cause delay component of the actual start time A i . 
     The system and method described herein output a stochastic model, which given the delay data up to time t, produces estimates of probability distributions of the actual start times A i  for each scheduled task i, which has not been started at time t, yet, as illustrated in  FIG. 6 , and given the same input, estimates of probability distributions of states X r,u  for all times u&gt;t and all resources r. Notice that the outputs of the disclosed system and method are functions. 
     The method operations above hence include two steps: obtaining the stochastic model (“training”) and applying the stochastic model (“forecasting”). 
     Subsequently to forecasting, the method operations above at time t may also include considering a risk measure of the random variables comprising times A i  for each scheduled task i, which has not been started at time t, yet, and states X r,u  for all times u&gt;t and all resources r (“risk estimation”). Formally, a risk measure is a mapping from a set of random variables to the real numbers. One may consider, e.g., the mode, or the most likely time. One may also consider the expectation of the random variable. These two would be “risk neutral” risk measures. One may also consider the value at risk with coefficient c of the random variable, which is the value of the cumulative distribution function at a certain probability c. One may also consider the conditional value at risk with coefficient c′ of the random variable, which is the value of the cumulative distribution function integrated over an interval [c′, 1] for continuous random variables. Those are “basic-risk averse” approaches. 
     Subsequently to risk estimation, the method operations above at time t may also include displaying the random variables comprising times A i  for each scheduled task i, which has not been started at time t, yet, and states X r,u  for all times u&gt;t and all resources r to a user approximately using a user interface (e.g., a display or printer). One may, e.g. update a graphical representation of the schedule in the form of a Gantt chart, with the blocks representing the tasks moved by the amount suggested by the risk measure, or vary the size and the transparency of the blocks, so as to suggest what is the probability of their starting at the given times. Considering there is a fixed resolution of the screen, one can capture the random variable only approximately, in effect producing a histogram approximation. One may also represent the possible evolutions of the state of each resource over time by plotting several sample paths in a plot, where the horizontal axis is the time and the vertical axis is the state or the value of some function thereof. 
     Subsequently to risk estimation, the method operations at time t above may also include storing the approximate representation of random variables comprising times A i  for each scheduled task i, which has not been started at time t, yet, and states X r,u  for all times u&gt;t and all resources r in memory or a storage medium. The representations include the histogram approximations and sample paths approximations, as above. 
     Subsequently to risk estimation, the method operations may also include updating a schedule to consider the times as random variables or any of their risk measures as discussed above. 
     Development of the stochastic model for scheduling a machining application job is now discussed. This gives the complete detail for a simple example. In the usual job scheduling problem, let M={M 1 , M 2 , . . . , M m } be a finite set of “machines” and J={J 1 , J 2 , . . . , J n } a finite set of “jobs”. For each 1≤j≤n, job J j  to be processed by a prescribed sequence of machines, S j ∈M×M× . . . M, 1≤1≤n. For example, if sequence S 1  for job J 1  is (M 3 , M 2 , M 1 ), the job has to be processed at machine M 3  first, then at M 2 , and M 1 . Each machine M z  can accept any job with sequence (M x , . . . , M y , M z , . . . ) that has been processed on machines M x , . . . , M y  and has not been processed on M z . The processing of a job j by a machine i is denoted by: (i,j). Further, a generalization of the job-scheduling concept towards resources is considered. (Similar generalizations underlie much of project scheduling, but the exact formalization varies, so one formalization is introduced here.) Let    i,j  ⊂ M×J, 1≤i≤m, 1≤j≤n denote the set of jobs, which need to be processed before job I j  can be processed at machine M i . Clearly for every job j with sequence S j =( . . . , M y , M z , . . . ). 
     With respect to processing and set-up times, two quantities are defined. The first quantity is time- and machine-dependent processing time p j,a,t =Y a,t +ƒ(j,t), which is the time it takes to process job J j  at machine α at time t. Notice that for a discretization of time, p j,a,z  could be seen as a rank-2 tensor, which can be factored into a random matrix Y a,z  and a deterministic function ƒ(j,t), representable by a matrix. The second quantity is sequence-dependent set-up time τ α,β +V α,β , 1≤m, n≤m, which is the time taken to move any job from machine α to machine β, or possibly, the time it takes to prepare machine α for processing any job, which has been previously processed at machine β. The sequence dependent set-up times is composed of a deterministic component, represented by a full m×m matrix T=(τ α,β ) and a noise term V α,β . Sequence-dependent set-up time τ α,β +V α,β , 1≤m, n≤m, which is the time taken to move any job from machine α to machine β, or possibly, the time it takes to prepare machine α for processing any job, which has been previously processed at machine β. The sequence dependent set-up times is composed of a deterministic component, represented by a full m×m matrix T=(τ α,β ) and a noise term V α,β . 
     Four representations of a “solution” of the scheduling problem are now introduced. “The Plan”: Let X denote all possible assignments of jobs to machines, such that each job is performed exactly once at each machine within its prescribed sequence. Elements x∈X may be written as n×m matrices, where column i represents jobs that machine M i  accepts, in the given order. For example, the sparse matrix 
             x   =     (         2       2       1           3       2       1           1       3                     )           
presents machine M 1  performs three jobs J 2 , J 3 , J 1  in the order J 2 , J 3 , J 1 , while machine M 2  will accept the jobs in the order J 1 , J 2 , J 3  and machine M 3  will accept the jobs in the order J 1 , J 3 . “The Schedule”: for each 1≤i≤m and 1≤j≤n, s i,j  denotes the scheduled time of processing job J ƒ  on machine M i , with s i, j  undefined, if S j  does not include M i . Unless stated otherwise, we assume s i,j  is a deterministic non-negative integer. The schedule can be represented as a sparse m×n matrix  =s i,j . “The Actual Times”: for each 1≤i≤m and 1≤j≤n, D i,j  and A i,j  are random variables representing the actual time of starting and finishing job J j  on machine M i , with A i,j , D i,j  undefined, if S j  does not include M i . “The Observations”: for each 1≤i≤m and 1≤j≤n, y i,j  represents one sample of the random variable Y i,t  drawn by processing job J j  on machine M i  at time A i,j . Observations of the actual times could be introduced as well, but those will not be necessary for the present discussion.
 
     In scheduling, it is customary to represent problem instances by graphs and combinations of instances and “the plan” by planar embeddings of such graphs, sometimes related to the so called “Gantt charts” as illustrated in  FIG. 11 . In such graphs and embeddings, oriented edges (i.e., horizontal) represent the sequences, in which each job should be accepted by machines, while non-oriented (i.e., diagonal) edges represent the exclusion constraints: the fact two jobs may not run at the same machine at one time. (Non-oriented edges represent further requirements within    i,j .) In the planar embeddings, each job is represented by one row, with time plotted on the horizontal axis. In the example in  FIG. 11 , the second job in the second row cannot be performed at the same time as the second job in the third row, the first job in the second row cannot be performed at the same time as the first job in the third row or the third job in the first row, and the first job in the first row cannot be performed at the same time as the third job in the second row. In practice, “machines” may represent actual CNC machines and “jobs” may represent the physical artifacts or work pieces to be machined. “Machines” may also represent workers in a factory and “jobs” may represent their work. In the context of railways discussed above, “machines” may represent stations in a railway network and “jobs” may represent train runs. 
     Throughout, it is assumed that the following are given: “The Time”: t∈N; “The Instance” given by sets M, J and resource requirements    i,j  ⊂ M×J, 1≤i≤m, 1≤j≤n, and the input to scheduling algorithms; “The Plan” X, which is generally the output of scheduling algorithms; “The Schedule” S i,j  for each 1≤i≤m and 1≤j≤n, which is generally the output of scheduling algorithms; and “The Observations” y i,j  for each 1≤i≤m and 1≤j≤n such that J ƒ  on machine M i  was started at time A i,j ≤t. In the training phase, an internal representation of a certain Markov chain is produced. In the forecasting phase, this Markov chain is used to derive distributional forecasts of D i,j ,A i,j , where D i,j ,A i,j &gt;t. 
     As suggested above, the following model is assumed for processing of job J j  at machine α in sequence S j =( . . . , α, β, . . . ). 
                     D     α   ,   j       =       s     α   ,   j       +     f   ⁡     (       s     α   ,   j       -       max       (     i   ,   k     )     ∈       α   ,   j           ⁢     {     A     i   ,   k       }         )       +     Y     α   ,     s     i   ,   j                     (   1   )                 A     α   ,   j       =       D     α   ,   j       +     τ     α   ,   β       +     V     α   ,   β                 (   2   )               
where: s α,j  is the scheduled time of starting processing job J j  on machine α; D α,j  is the actual start of processing job J j  on machine α;    α,j  is the resource requirements, which need to finish before job J j  can be processed on machine α; τ α, β  is the sequence dependent set-up time, e.g. scheduled duration of movement from machine α to β; A α,j  is the actual availability of job J j  to be accepted on the machine β following machine α in sequence S j =( . . . , α, β, . . . ); V α,β  is a noise term for moving any job J j  from machine β following machine α in sequence S j =( . . . , α, β, . . . ); and Y α,t  is a noise term for processing job J j  at machine α at time t.
 
     Notice that in this model, the crucial random variable is Y. In particular delays (D α,j −s α,j ) + , where (⋅) +  denotes the positive part, are explained by two factors in this model. The first factor is the propagation of delays, i.e. the availability of resources required to process the job. This “reactionary component” is, however, rather mechanistic and is represented as: ƒ(x)=c−x, if x≤c, or =0, otherwise, where c is a constant in the chosen units of time, to capture the impact that late arrival of resources has on start of processing. The second factor is the random variables Y α,t , whose observations are y α,t . To simplify notation, y α,t  is replaced for processing job J j  if at machine α at time t, such as s i,j , by simply y j . In particular, the constant c in the expression may be the minimum amount of time, required by the machining of the part at the given machine, under the best possible conditions. 
     For any j and k&gt;j, it is desired to know the conditional density ƒ(y k |(y 0 ,c 0 ), . . . , (y j , c j )). 
     A HMM is used to capture this conditional density: 
               f   ⁢     {         y   k     |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )       =       ∑     x     t   k         ⁢       f   ⁡     (       y   k     ,       x     t   k       |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )       .             
where t k  represents the scheduled time of departure k and X t     k    is a latent state variable at that time. It is assumed that the observations at time t are conditionally independent of all other observations given x t . So,
 
                     f   ⁡     (       y   k     ,       x     t   k       |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )       =       ⁢     f   ⁡     (         y   k     |     x     t   k         ,     (       y   0     ,     c   0       )     ,   …   ,     (       y   j     ,     c   j       )       )                       ⁢     p   ⁡     (         x     t   k       |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )                   =       ⁢       f   ⁡     (       y   k     |     x     t   k         )       ⁢       p   ⁡     (         x     t   k       |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )       .                   
It is assumed that the states follow a Markov chain, so
 
     
       
         
           
             
               p 
               ⁡ 
               
                 ( 
                 
                   
                     x 
                     
                       t 
                       k 
                     
                   
                   , 
                   
                     ( 
                     
                       
                         y 
                         0 
                       
                       , 
                       
                         c 
                         0 
                       
                     
                     ) 
                   
                   , 
                   … 
                   , 
                   
                     ( 
                     
                       
                         y 
                         j 
                       
                       , 
                       
                         c 
                         j 
                       
                     
                     ) 
                   
                 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 
                   
                     x 
                     0 
                   
                   , 
                   … 
                   , 
                   
                     x 
                     j 
                   
                 
               
               ⁢ 
               
                 
                   
                     ρ 
                     0 
                   
                   ⁡ 
                   
                     ( 
                     
                       x 
                       0 
                     
                     ) 
                   
                 
                 ⁢ 
                 
                   p 
                   ⁡ 
                   
                     ( 
                     
                       
                         x 
                         
                           t 
                           0 
                         
                       
                       | 
                       
                         x 
                         0 
                       
                     
                     ) 
                   
                 
                 ⁢ 
                 
                   ⋯p 
                   ⁡ 
                   
                     ( 
                     
                       
                         x 
                         
                           t 
                           k 
                         
                       
                       | 
                       
                         x 
                         
                           t 
                           j 
                         
                       
                     
                     ) 
                   
                 
                 ⁢ 
                 
                   q 
                   ⁡ 
                   
                     ( 
                     
                       
                         y 
                         0 
                       
                       , 
                       
                         
                           c 
                           0 
                         
                         | 
                         
                           x 
                           
                             t 
                             o 
                           
                         
                       
                     
                     ) 
                   
                 
                 ⁢ 
                 
                   
                     ⋯q 
                     ⁡ 
                     
                       ( 
                       
                         
                           y 
                           j 
                         
                         , 
                         
                           
                             c 
                             j 
                           
                           | 
                           
                             x 
                             
                               t 
                               j 
                             
                           
                         
                       
                       ) 
                     
                   
                   . 
                 
               
             
           
         
       
     
     Once the stochastic model is trained, the following factors are known: ρ 0 (x) is the probability mass function (PMF, for specifying a distribution of a discrete random variable) of initial states of the Markov chain; p(x t |x t-1 ) is the transition matrix of the Markov chain; q(c|x) is the conditional PMF of delay code given the state; and r (y|c, x) is the conditional distribution of delay times given the delay code and state. Once these factors are known, probabilities of states x of trains stations in the railway context (or machines in the machine context) can be iteratively forecast using p and probabilities of actual delays y can be iteratively forecast using r. See  FIG. 10  for example iterations. 
     HMMs in a railway context are now discussed. A simplified model of railway operations includes: a schedule as a collection of train movements with planned departure and arrival times, where a train movement can only commence when certain other train movements (“predecessor train movements” or “predecessors”) are complete; delays that are unanticipated changes in the start time or duration of a train movement; reactionary delays result from delays in predecessors (e.g., a train is delayed because it is behind a delayed train and must maintain a certain minimum headway; root delays that occur even when all predecessors are complete (e.g., construction on a section of track reduces speed below schedule, many boarding passengers delay departure from a station, mechanical problems affecting a signal). 
     The hidden states are time and location specific, and could be inferred from observations, e.g., reduced speed on a section of track, high volume of passengers at a station, faulty equipment at a station. In this embodiment, a method for modeling probabilities of states the stations are in and probabilities of delays as a function of predecessors&#39; delays and the states of the stations are presented. For example: several consecutive trains are delayed in departing a station; it might be inferred that this station is in a state of high passenger volume or a state related to faulty equipment or people on the railway tracks; based on historical data, how long this state will persist can be estimated; and based on this state estimate, the probabilities of future schedule disruptions and their reactionary impacts can be estimated. More complex embodiments are also contemplated. 
     In more detail: departure times are modeled as scheduled start time+reactionary delay+station-state-specific term; arrival times are modeled as start time+scheduled duration+train-movement-specific term; and models for the station-state-specific and train-movement-specific terms are learned from past data. In formal notation:
 
 D   i   =s   i +ƒ( s   i −max{ A   j } for  j∈R   i )+ Y   αi,S     i    and  A   i   =D   i +τ i   +V   i  
 
where: si is the scheduled departure time of train movement i; D i  is the actual departure time of train movement i; R i  is the set of “predecessor” movements for movement i; τ i  is the scheduled duration of train movement i; A i  is the actual arrival time of train movement i; α i  is the departure station of train movement i; Y αi,Si  is the noise term for train movement i departure time; and V i  is the noise term for train movement i.
 
     In this model, departure delays are explained by two factors. The first factor relates to the predecessors, i.e., the availability of resources required by the train movement prior to departure. In this model, the following equation is used to capture the impact that late arrival of resources has on departure time:
 
ƒ( x )=2 −x , if  x≤ 2, or =0, otherwise.
 
     The second factor relates to the random variables Y α,S . Departure delays are expected to be correlated at a given station throughout the day. One can interpret the ƒ as follows: Should the train arrive in the station and have the crew on-board more than three minutes ahead of the scheduled departure, it can leave on time at s_i, independent of how early it has arrived in the station. If the train was late and there were no passengers boarding the train, the train still needs to arrive at the station, board the crew, and leave the doors open for 3 minutes, such that the passengers who may be interested in boarding the train, may board the train. 
     As illustrated in  FIG. 7 , in this use case, system and method outputs include probability distributions over current and future states of the resources as a function of delay data. Outside of a function that provides probability distributions of the departure D i  and arrival A i  of the current train movement, probability distributions are also provided of the state X α, t  of each resource α at each time t. Notice that each train movement may depend on multiple resources, as is illustrated in  FIG. 8 . Multiple resources related to the departure of train movement i, e.g., crew α i  and outbound track β i  within a station, are modelled with multiple hidden Markov models (HMMs). The schedule data and noise-terms Y related to each resource are convolved (represented by asterisk within circle) with a function of the maximum of arrival times of the predecessors to obtain the departure time. From the departure time, the arrival time is computed, as is illustrated next in  FIG. 9 , which is used in computation of departure time, for train movements, for which the current train movement is a predecessor. 
       FIG. 10  presents a process flow diagram of method operations  100  for prediction of the departure D i  and arrival A i  of a single train movement i, using only a single resource α i . In operation 1, for each resource, the initial state distribution ρ 0 (x) is sampled to estimate the previous state x t-1 . This operation involves machine learning using stochastic model and previously obtained historical data. In operation 2, for each resource, the current state x t  conditioned on x t-1  is obtained using conditional probability distribution p(x t |x t-1 ). This operation involves machine learning using stochastic model and previously obtained historical data. In operation 3, for each resource, the delay code c i  conditioned on x t-1  is obtained using conditional probability distribution q(c t |x t ). This operation involves machine learning using stochastic model and previously obtained historical data. In operation 4, for each resource, the departure-related noise-term Y α     i     ,i  conditioned on x t  and c t  is obtained using conditional probability distribution r(y t |c t , x t ). This operation involves machine learning using stochastic model and previously obtained historical data. In operation 5, for each train movement, the departure time is estimated as a convolution (asterisk in circle) of Y α     i     ,i  the scheduled departure time s i , and a function of the maximum of arrival times {A j } of predecessors. This operation involves calculations using a processor. In operation 6, for each train movement, the arrival time is estimated as a convolution (asterisk in circle) of the departure time D i , the travel time τ i , and the arrival-related noise-term V i . In a simple embodiment, the travel time can be a deterministic value of the scheduled travel time, while in more elaborate embodiments, the travel time can be obtained by an independent prediction module. The noise term is a random variable that can delay the arrival of a train and involves machine learning using a stochastic model and historical data. In more elaborate embodiments, this may be Y_{\beta, i+t} for \beta being the next station and t being a travel time. In operation 7, for each train movement, the arrival time is propagated to train movements, for which the current train movement is a predecessor. This operation involves processing using a processor. In operation 8, for each resource, the state is propagated to the next iteration, increment the iteration counter t and repeat from operation 2 through to operation 8. The iterations terminate at the end of the point where forecasts are no longer necessary, e.g., at the end of the day, or at the end of the following day. In operations 1-4 and 7, the stochastic model can be a Markov model, which had been previously hidden, has been estimated in the training phase, and is now being applied. The data s i  in operation 5 and τ i  in operation 6 are schedule data. The method operations  100  may also include adding the probability distribution of time delays for the start time to the scheduled start time to provide a probability distribution of forecast start times. Similarly, the method operations  100  may also include adding the probability distribution of time delays for the end time to provide a probability distribution of forecast end times. 
     One example of training HMMs in the railway context is now discussed. For a given station, on a given day, n departures are observed. Departure i occurs at time D i  and reports a delay code of C i . The realization of Y αi,Si  for this train movement is:
 
 y   αi,si   =d   i   −s   i −ƒ( s   i −max{ a   j } for  j∈R   i ).
 
     At a given time of day, previous observations of the station-specific delay variables convey some information about future station-specific delays. To simplify notation, replace y αi,si  is simply replace by y i . For any j and k&gt;j, the conditional density
 
ƒ( y   k |( y   0   ,c   0 ), . . . ,( y   j   ,c   j )) is desired to be known.
 
     A HMM is used in this example to capture this conditional density:
 
ƒ( y   k |( y   0   ,c   0 ), . . . ,( y   j   ,c   j ))=Σƒ( y   k   x   tk |( y   0   ,c   0 ), . . . ,( y   j   ,c   j )) for summation over  x   tk  
 
where t k  represents the scheduled time of departure k and X tk  is a latent state variable at that time. It is assumed that the observations at time t are conditionally independent of all other observations given x t . Thus,
 
                     f   ⁡     (       y   k     ,       x   tk     |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )       =       ⁢     f   ⁡     (         y   k     |     x   tk       ,     (       y   0     ,     c   0       )     ,   …   ,     (       y   j     ,     c   j       )       )                       ⁢     p   ⁡     (         x   tk     |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )                   =       ⁢       f   ⁡     (       y   k     |     x   tk       )       ⁢     p   ⁡     (         x   tk     |     (       y   0     ,     c   0       )       ,   …   ,     (       y   j     ,     c   j       )       )                     
It is assumed that the states follow a Markov chain, thus
 
     p(x tk , (y 0 , c 0 ), . . . , (y j , c j ))=Σρ 0 (x 0 )p(x t0 |x 0 ) . . . p(x tk |x tj )q(y 0 , c 0 |x t0 ) . . . q(y j , c j |x tj ) for summation over x 0 , . . . , x tj . Above, it is assumed that the Markov chain governing state transitions and the conditional distributions governing observations are time-invariant, but these transitions and distributions may also be implemented as varying in time. Subsequently, algorithms may be used to find a local maximum likelihood solution in polynomial time, e.g., of Viterbi, Baum-Welch, or Baldi-Chauvin, as discussed above. 
     The two examples above are intentionally simple for teaching purposes and do not fully demonstrate the potential of the invention. Notice, in particular, that only two sub-tasks per task are considered in the scheduling example, where the first sub-task is the true start of the task and the other sub-task starts with the availability of the task at the following machine, and two sub-tasks per task in the railways example, where the first sub-task is the departure and the second is the arrival. Clearly, there could be multiple sub-tasks introduced into each of the examples, e.g., the arrival of a train at a station could consist of seven of more sub-tasks, including entering the area of the approach track, start of breaking, stop moving in a station, open a door, close a door, start the engine, leave the area of the station. In the railway example, the arrival time is considered A i =D i +τ i +V i , rather than a convolution of D i  with a more complex model of the travel time. One could employ an independent distributional forecasting model for the travel time. Further, a single delay integral code is considered in the railway example, rather multi-variate additional information, which may capture, e.g. weather at the railway station, numbers of ticket purchases over the past 1, 2, 3, . . . , minutes, people passing through the entrance barriers over the past 1, 2, 3, . . . , minutes, and similar. Even further, a single most recent observed delay and a single most recent observed delay code are considered in the example, rather than k most recent delays and delay codes, which may be useful in modelling a limited amount of history, in a so called windowing, at the price of increasing the number of states exponentially in the size of the window. Even further, only the fixed dependence resource is considered, whereas one could consider multiple alternatives for the dependence. That is, for at most k alternatives per task, one could have k bipartite graphs, and there could be a choice of all resources from any single of the k bipartite graphs for any single task. 
     The present invention provides several advantages. One advantage is that it takes into account that the causes of delays may be complex and are not precisely known. Hence, a probability distribution of delays of start and end times can provide a user such as an operations controller or manager the information needed to make decisions regarding the delays. Another advantage is that it allows published schedules to be updated in real time based on the latest received information. While the invention has been discussed above concerning advantages with respect to the transportation and manufacturing industries, the advantages are also applicable to other industries involved with tasks having scheduled start and end times for the tasks and, in particularly, where delays such as reactionary delays may propagate through a system and affect multiple tasks. 
     The present invention may be a system, a method, and/or a computer program product. The computer program product may include a non-transitory computer readable storage medium (or media) having computer readable program instructions thereon for causing a processor to carry out aspects of the present invention. 
     The computer readable storage medium can be a tangible device that can retain and store instructions for use by an instruction execution device. The computer readable storage medium may be, for example, but is not limited to, an electronic storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. A non-exhaustive list of more specific examples of the computer readable storage medium includes the following: a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), a static random access memory (SRAM), a portable compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanically encoded device such as punch-cards or raised structures in a groove having instructions recorded thereon, and any suitable combination of the foregoing. A computer readable storage medium, as used herein, is not to be construed as being transitory signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission media (e.g., light pulses passing through a fiber-optic cable), or electrical signals transmitted through a wire. 
     Computer readable program instructions described herein can be downloaded to respective computing/processing devices from a computer readable storage medium or to an external computer or external storage device via a network, for example, the Internet, a local area network, a wide area network and/or a wireless network. The network may comprise copper transmission cables, optical transmission fibers, wireless transmission, routers, firewalls, switches, gateway computers and/or edge servers. A network adapter card or network interface in each computing/processing device receives computer readable program instructions from the network and forwards the computer readable program instructions for storage in a computer readable storage medium within the respective computing/processing device. 
     Computer readable program instructions for carrying out operations of the present invention may be assembler instructions, instruction-set-architecture (ISA) instructions, machine instructions, machine dependent instructions, microcode, firmware instructions, state-setting data, or either source code or object code written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C++ or the like, and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The computer readable program instructions may execute entirely on the user&#39;s computer, partly on the user&#39;s computer, as a stand-alone software package, partly on the user&#39;s computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user&#39;s computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider). In some embodiments, electronic circuitry including, for example, programmable logic circuitry, field-programmable gate arrays (FPGA), or programmable logic arrays (PLA) may execute the computer readable program instructions by utilizing state information of the computer readable program instructions to personalize the electronic circuitry, in order to perform aspects of the present invention. 
     Aspects of the present invention are described herein with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer readable program instructions. 
     These computer readable program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. These computer readable program instructions may also be stored in a computer readable storage medium that can direct a computer, a programmable data processing apparatus, and/or other devices to function in a particular manner, such that the computer readable storage medium having instructions stored therein comprises an article of manufacture including instructions which implement aspects of the function/act specified in the flowchart and/or block diagram block or blocks. 
     The computer readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable apparatus or other device to produce a computer implemented process, such that the instructions which execute on the computer, other programmable apparatus, or other device implement the functions/acts specified in the flowchart and/or block diagram block or blocks. 
     The flowchart and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of instructions, which comprises one or more executable instructions for implementing the specified logical function(s). In some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts or carry out combinations of special purpose hardware and computer instructions. 
     The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one more other features, integers, steps, operations, element components, and/or groups thereof. 
     The description of the present invention has been presented for purposes of illustration and description, but is not intended to be exhaustive or limited to the invention in the form disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the invention. The embodiment was chosen and described in order to best explain the principles of the invention and the practical application, and to enable others of ordinary skill in the art to understand the invention for various embodiments with various modifications as are suited to the particular use contemplated 
     The flow diagrams depicted herein are just one example. There may be many variations to this diagram or the steps (or operations) described therein without departing from the spirit of the invention. For instance, the steps may be performed in a differing order or steps may be added, deleted or modified. All of these variations are considered a part of the claimed invention. 
     While the preferred embodiment to the invention had been described, it will be understood that those skilled in the art, both now and in the future, may make various improvements and enhancements which fall within the scope of the claims which follow. These claims should be construed to maintain the proper protection for the invention first described.