Abstract:
A method and system for scheduling a plurality of concurrent job sequences with constraints of precedence and mutual exclusion is described. In particular, a scheduler for a system performing biological analyses on clinical samples including non-preemptible jobs which must use a set of resources (machines) and having constraints on release times and execution times is described. The scheduling method uses on-purpose algorithmic solution that extracts the core of DBM data structures and Floyd-Warshall algorithm from the engine of a full-fledged Model Checking tool and tailors it to the specific requirements of the case.

Description:
FIELD OF TECHNOLOGY 
       [0001]    The present invention relates to a method and system for scheduling a plurality of concurrent sequences of timed jobs that are subject to constraints of mutual exclusion and mutual delay. In particular, the present invention relates to a scheduler for a system that performs multiple concurrent biological analyses that apply different analysis protocols on different samples while sharing a set of electromechanical devices. 
       BACKGROUND 
       [0002]    While the problem of scheduling a plurality of concurrent and timed job sequences with constraints of mutual delay and mutual exclusion with allowed bounded time slacks can be applied to a large variety of contexts, in the present application, without loss of generality, we address the problem with reference to a concrete example of a run-time scheduler for an electromechanical system for biological analyses, such as for example the VIDAS system engineered by bioMérieux 
         [0003]    The system of the present example is able to run a group of multiple concurrent analyses, each alternating steps of biological reaction of the sample and transfer operations carried out by shared electro-mechanical devices; the duration of each reaction step is determined by the biological protocol of the analysis, but bounded waiting times can be added after the completion of each reaction step. This gives raise to the scheduling problem of determining waiting times after each step so as to comply with maximum allowed delays and to avoid conflicts in the usage of shared resources, while maximizing the number of concurrent analyses and minimizing their overall completion time. 
         [0004]    With more detail, each type of biological analysis is comprised by a pretreatment phase and an analytic protocol. At the beginning of the analysis, the sample is contained in a test-tube uniquely identified by a barcode, while other tubes contain dilution and incubation fluids and other tubes are empty. During the pretreatment phase, an automatic pipettor repeatedly pours off the sample among the various tubes; each pipettor action lasts for a deterministic amount of time. Waiting times are allowed between successive actions, but they are constrained to range between minimum and maximum values determined by incubation/reaction periods and sample deterioration properties. After completion of the pretreatment phase, the analytic protocol follows a fixed sequence of steps during which the sample is combined with reagents and multiple measurements are taken by a read-head whose operations take deterministic durations (see  FIG. 1 ). 
         [0005]    For efficient exploitation of electromechanical components, multiple analyses, possibly of different types, are carried out concurrently. To this end, the system composes multiple sections, one for each analysis. Each section is in turn divided in slots, carrying different samples that can be subject to different pretreatments and can operate on samples of different subjects. However, since the read-head is designed to take measures on a whole section at once, all slots in the same section are constrained to run the same analytic protocol (see  FIG. 2 ). 
         [0006]    The pipettor and the read-head are shared among different slots and sections and cannot be used simultaneously by two different analyses. This raises a scheduling problem that must determine feasible values of the initial delay of each pretreatment and waiting times between events of the same analysis. Such values must avoid the need of contemporary actions of the pipettor and the read-head, avoid waiting times beyond allowed bounds and keep the overall time span of the group of analyses as short as possible. 
       Problem Formulation 
       [0007]    The problem can be formally cast in the general notation of job shop scheduling as a set of n independent and non-preemptible jobs J 1  . . . J n , with release times r 1  . . . r n  and deterministic Execution Times e 1  . . . e n , each job being statically assigned to one out of m independent machines M 1  . . . M m , subject to constraints of precedence and mutual exclusion:
       for any two jobs J i  and J k , a precedence constraint may require that the delay between the completion of J k  and the start of J i  range within a minimum d −   ik  and a maximum d +   ik :       
 
         [0000]        d   −   ik   ≦r   i −( r   k   +e   k ) d+   ik  
       for any two jobs J i  and J k , a mutual exclusion constraint may require that the execution periods of J k  and J i  do not overlap; this may be achieved with either of the following two determinations, which guarantee that J i  starts after completion of J k  or viceversa:       
 
         [0000]        r   i   ≧r   k   +e   k  OR  r   i   +e   i   ≦r   k    
         [0010]    The scheduling problem amounts to determine the release times r 1  . . . r n , subject to all prescribed precedence and mutual exclusion constraints, so as to minimize the completion time of the overall plurality of job sequences.  FIG. 3   a  illustrates an instance of the problem: Jobs J 1 , J 2 , J 3 , and J 4  are allocated on Machines M 1  and M 2 ; Job J 2  has a mutual exclusion with respect to J 4  and a precedence constraint bounding its start time with respect to the completion time of J 1  and J 3 . A precedence constraints also bounds the delay of the start time of J 4  with respect to the completion of J 3 .  FIG. 3   b  shows a solution for the problem. 
         [0011]    The scheduling problem defined above is characterized by the intertwining of non-deterministic release times with mutual exclusion and precedence constraints, and it can thus be conveniently formulated on the ground of non-deterministic formalisms such as Time Petri Nets (TPN) or Timed Automata (TA) through a direct translation of scheduling concepts like jobs, constraints and resources, into places/locations, transitions, pre- and post-conditions. This opens the way to the application of modeling and verification tools that manage the combination of concurrency constraints and non-deterministic timers through symbolic analysis of the state space based on Regions (see R. Alur, D. L. Dill: Automata for modeling real-time systems. In Proc. Of 17th ICALP, 1990) or Difference Bounds Matrix (DBM) zones (see D. Dill. Timing Assumptions and Verification of Finite-State Concurrent Systems. In Proc. Workshop on Computer Aided Verification Methods for Finite State Systems, 1989). In this perspective, the scheduling problem can be formulated as the identification of a witness in the symbolic state space that satisfies a Real Time Model Checking formula (see e.g. A. Fehnker. Scheduling a steel plant with timed automata. In Real-Time Computing Systems and Applications, 1999. RTCSA &#39;99. Sixth International Conference on, pages 280-286, 1999) that selects “all behaviors that reach a completion state within a maximum time T, after traversing only states without conflicts in the usage of resources”. The optimal value of T can be obtained through polynomial iteration, which is supported in various tools and notably in the Uppaal environment (see e.g. K. Larsen. Resource-efficient scheduling for real time systems. In Rajeev Alur and Insup Lee, editors, Embedded Software, volume 2855 of Lecture Notes in Computer Science, pages 16-19. Springer Berlin-Heidelberg, 2003.; or I. AlAttili, F. Houben, G. Igna, S. Michels, F. Zhu, and F. W. Vaandrager. Adaptive scheduling of data paths using uppaal tiga. In QFM, pages 1-11, 2009). In the Oris tool (G. Bucci and L. Carnevali and L. Ridi and E. Vicario. Oris: a Tool for Modeling, Verification and Evaluation of Real-Time Systems. International Journal of Software Tools for Technology Transfer, 12(5):391-403, 2010), the optimum value T is derived in direct manner by evaluating the minimum and maximum time spent along each sequence of transitions that satisfy the expected sequencing constraints (E. Vicario. Static Analysis and Dynamic Steering of Time Dependent Systems Using Time Petri Nets. IEEE Trans. on SW Eng., 27(1):728-748, August 2001). 
         [0012]    Unfortunately, straight application of this type of symbolic state space analysis techniques faces the problems of state space explosion (as the number of states grows exponentially with the number of jobs) and curse of dimensionality (as the size of data structures encoding each state grows quadratically with the number of sequences). This results in memory and time complexity that are by far not affordable for an on-line scheduler, even if computation is performed on a conventional PC, which is often not the case in industrial applications. 
         [0013]    Just to figure out the order of magnitude of practical complexity, state space analysis based on DBM zones with the Oris tool was experimented on the TPN representation of a simplistic case made of two Sections, three Slots per Section, and deterministic delays: the TPN enumeration procedure would produce a huge state space of 255128 DBM state classes, which requires about 20 minutes on a dual-core 2 GHz desktop computer, with 4216 possible paths from the initial state to the end terminal node. 
         [0014]    Even applying state of the art techniques for reducing complexity through model simplification, end even partitioning the state space through DBM Zones (which are much more compact than Regions), processing times and memory demand of this class of state space analysis techniques are not compatible with the intended use within an industrial application as the one described above: this requires that the scheduling problem be solved on-line in a near-real-time approach, and that its execution be affordable within an embedded component that runs on limited processing and memory resources. Hardware resources normally available in this kind of biological testing machines (but more generally in industrial electromechanical machines) is very limited and not compatible with the huge demand of such a complex mathematical algorithm implementation. 
         [0015]    A solution that could trade generality of analysis for efficiency in the specific problem so as to run faster and with lower demand of memory and processing power would be therefore highly desirable. 
       OBJECTS OF THE DISCLOSURE 
       [0016]    It is an object of the present disclosure to overcome at least some of the problems associated with the prior art. 
         [0017]    A further object of the present invention is to provide a method that requires limited memory and processing time for computing an acceptable and efficient schedule for a plurality of concurrent sequences of timed steps subject to constraints of mutual exclusion and mutual delay and allowed to be delayed by bounded waiting times. 
         [0018]    This is represented as a scheduling problem arising in a real industrial electromechanical system, which includes finding feasible timings for the actions performed by a set of shared electromechanical components so as to complete a set of concurrent biological analyses in the minimum possible timespan while satisfying mutual exclusion and mutual delay constraints. In the specific example here described, the scheduling algorithm is required to run on board of the hardware shipped with an electromechanical system presently equipped with a 200 MHz CPU and with less than 10 MB of memory. Moreover, in order to guarantee an adequate responsiveness for the operators interacting with the system, the algorithm is allowed to find sub-optimal yet efficient solutions but it is required to end conclusively within a short amount of time in the order of few seconds. 
         [0019]    Another object of the invention is to provide an electromechanical system that executes a plurality of concurrent biological analyses, each alternating periods of reaction and transfer operations performed by shared devices under timing and mutual delay constraints determined by the biological analysis protocol, and under mutual exclusion constraints determined by the exclusive usage of shared devices. 
       SUMMARY 
       [0020]    The present disclosure provides a scheduling method as set out in the accompanying claims. 
         [0021]    According to a first embodiment, the present disclosure provides a scheduling method, in a system adapted to execute a plurality of non-preemptible jobs J 1  . . . J n , with deterministic execution times e 1  . . . e n , each job being statically assigned to one of m independent machines M 1  . . . M m , for determining a release time r 1  . . . r n  for each job J 1  . . . J n  so as to minimize the completion time of overall set of jobs, while respecting the following precedence constraint:
       for a predetermined set of pairs J x  and J y  of the plurality of jobs, the delay between the completion of J y  and the start of J x  range within a minimum delay d −   xy  and a maximum delay d+ xy :       
 
         [0000]        d   −   xy   ≦r   x −( r   y   +e   y )≦ d+   xy  
 
         [0000]    and the following exclusion constraints:
       for a predetermined set of pairs J i  and J k  of the plurality of jobs, one of the following two determinations of exclusion constraint is fulfilled:       
 
         [0000]        r   i   ≧r   k   +e   k  OR  r   i   +e   i   ≦r   k          so that the execution periods of J k  and J i  do not overlap,
 
the method comprising the steps of: A) building a frontier set including a Difference Bound Matrix (DBM) zone that collects release times satisfying the precedence constraints; B) selecting one of the DBM zones and removing it from the frontier set: C) responsive to the selected DBM zone not satisfying at least one exclusion constraint of one pair J i  and J k  of the predetermined set of pairs, building one restricted zone of the selected DBM zone per each of the two determinations resolving the conflict (r i ≧r k +e k  or r i +e i ≦r k ); D) checking if the restricted DBM zones are non-empty and adding the non-empty restricted DBM zones to the frontier set; E) repeating the steps B to D until a selected DBM zone satisfies all the exclusion constraints.
         
         [0025]    The method of the present invention allows scheduling of concurrent sequences of timed jobs with constraints of mutual exclusion and delay, by capturing sets of equivalent possible solutions through Difference Bound Matrix (DBM) zones and by organizing the identification of a feasible solution as a search in a graph-theoretical structure. This enables a series of algorithmic optimizations and heuristic approaches to speed the search of solutions while keeping an exact (i.e. not quantized nor approximated) representation of the set of possible timing configurations. 
         [0026]    According to another embodiment the step of repeating is performed until all DBM zones in the frontier set satisfies all the exclusion constraints. Alternatively the step of repeating is stopped when a predetermined maximum search time has been reached. In another embodiment an optimal solution (schedule) is selected within all DBM zones satisfying all exclusion constraints. 
         [0027]    According to another embodiment, the selection of step B is based on an heuristic evaluation, including the so called longest-interval heuristic (safe) or the so called shortest-interval heuristic (greedy). 
         [0028]    According to another embodiment, a computer program for performing the method is provided. 
         [0029]    According to a further embodiment of the invention, a system implementing the method above is provided. Such system can be used for example in an electromechanical apparatus for performing biological analyses. 
         [0030]    The present invention offers a number of benefits. One of the advantages of a preferred embodiment of the present invention is that a complex scheduling operation can be performed with limited time and hardware resources. In principle the problem can be formulated as a case of Real Time Model Checking and solved through well-established state space analysis approaches. However, such enumerative approaches are not viable in the specific setting, due to hardware restrictions imposed by the expected embedded platform. With the method of the present invention a direct on-purpose adaptation of DBM theory largely reduces the complexity of the algorithm, also by allowing the adoption of heuristics that drive the search toward feasible solutions. In a preferred embodiment of the present invention, subsequent algorithmic refinements permit to further reduce computation times, attaining the expected performance requirements. 
         [0031]    An additional advantage of the method and system according to an embodiment of the present invention is that of allowing the case where one or more analyses are added to an ongoing schedule by reclaiming unused resources without jeopardizing the execution of already scheduled analyses. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0032]    Reference will now be made, by way of example, to the accompanying drawings, in which: 
           [0033]      FIG. 1  is a schematic representation of the physical architecture of system components running a single analysis in an electromechanical apparatus for performing biological tests; 
           [0034]      FIG. 2  is a schematic representation, of the physical architecture of the overall system composing multiple analyses in slots and sections in an electromechanical apparatus for performing biological tests according to an embodiment of the present invention; 
           [0035]      FIG. 3   a  is a schematic representation of an instance of the scheduling problem of determining the release times time r 1  . . . r n , subject to all prescribed precedence and mutual exclusion constraints, so as to minimize the completion time of overall set of jobs.  FIG. 3   b  is a possible solution that determines the value of release times; 
           [0036]      FIG. 4   a  shows a 2-dimensional DBM encoding the solution set for a pair of variables t i  and t j .  FIG. 4   b  shows a 2-dimensional DBM highlighting the presence of a non-effective constraint.  FIG. 4   c  shows the graph-theoretical representation of a DBM involving the three variables t i , t j  and t k . 
           [0037]      FIG. 5  shows an example of the binary tree that is constituted by subsequent choices for the exclusion constraints of the scheduling method according to a preferred embodiment of the present invention; 
           [0038]      FIG. 6  shows the results of performance testing experiments with two different heuristics, according to a preferred embodiment of the present invention. 
       
    
    
     DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 
       [0039]    To the purposes of the subsequent treatment, we introduce here the main concepts on which the scheduling solution relies: Difference Bound Matrices (DBMs) and Floyd-Warshall algorithm. 
         [0040]    A DBM encodes a system of inequalities identifying a set of solutions for a vector of N variables r=         t 1  . . . t N              
         [0000]    
       
         
           
             D 
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           τ 
                           i 
                         
                         - 
                         
                           τ 
                           j 
                         
                       
                       ≤ 
                       
                         b 
                         ij 
                       
                     
                   
                   
                     
                       i 
                       , 
                       
                         j 
                         ∈ 
                         
                           
                             [ 
                             
                               1 
                                
                               
                                   
                               
                                
                               … 
                                
                               
                                   
                               
                                
                               N 
                             
                             ] 
                           
                           ⋃ 
                           * 
                         
                       
                       , 
                       
                         
                           b 
                           ij 
                         
                         ∈ 
                         ℚ 
                       
                     
                   
                 
                 
                   
                     
                       
                         τ 
                         * 
                       
                       = 
                       0 
                     
                   
                   
                     □ 
                   
                 
               
             
           
         
       
     
         [0041]    Any 2-dimensional projection of a DBM results in a polyhedron of the type represented in  FIG. 4   a . A DBM may be not completely formed when one or more constraints are not effective (see  FIG. 4   b ); in this case an infinite number of values for these constraints yield the same set of solutions. The limiting value, such that any further restriction would change the solution set, is proven to exist and be unique if and only if the following triangular inequality holds: 
         [0000]    
       
      
       b 
       ij 
       ≦b 
       ik 
       +b 
       kj 
       ∀i,j,k  
      
     
         [0042]    This condition identifies the concept of normal form, that is defined as the representation that satisfies the triangular inequality. 
         [0043]    The set of constraints of a DBM can be conveniently cast into a graph-theoretical formulation (see  FIG. 4   c ) in which nodes of the graph represent variables, and edges are labeled with values of DBM coefficients (i.e. the edge that connects nodes corresponding to variables i and j is labeled with coefficient b ij ): in this perspective, the normal form of a DBM becomes an all-shortest-path problem where b ij  coefficients are computed as the shortest path from node i to j. 
         [0044]    It is easily proven that the assumption that DBM is not empty excludes the existence of negative cycles (i.e. cycles in the graph whose coefficients sum up to a negative value). Thus, the problem can be solved by N repeated applications of the well-known Dijkstra algorithm, with total complexity O(N·N 2 ). Floyd-Warshall algorithm provides a direct solution that still runs in O(N 3 ), but it performs much simpler operations than Dijkstra within its inner cycle and is preferable when dealing with dense adjacency matrices, like those corresponding to DBM domains (see E. Vicario. Static Analysis and Dynamic Steering of Time Dependent Systems Using Time Petri Nets. IEEE Trans. on SW Eng., 27(1):728-748, 2001). 
         [0045]    We report in Alg. 1 the Floyd Warshall pseudo-code. Its core invariant is variable b ij   k , which encodes “the shortest path from j to i under the restriction that only vertices with index higher than k can be visited as intermediary nodes”. Under this statement, it is clear that b ij   N  can be given the initial values of coefficients b ij , which gives meaning to the first for-loop in the code. It is also clear that b ij   0  comprises the final solution of the problem, which in the code is attained by decreasing k from N to 1 through the second for-loop. The real core of the algorithm is thus restrained to the way how the nested double loop on i, j is able to derive b ij   k-1  from b ij   k . 
         [0000]    
       
         
               
             
               
               
             
           
               
                   
               
               
                 ALGORITHM 1: FLOYD-WARSHALL ALGORITHM 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 1  begin 
               
               
                   
                 2   for [i, j] ∈ [0, N − 1] × [0, N − 1] do 
               
               
                   
                 3    b ij   N−1  = b ij ; 
               
               
                   
                 4   end 
               
               
                   
                 5   for k = N − 1 → 0 do 
               
               
                   
                 6    for [i, j] ∈ [0, N − 1] × [0, N − 1] do 
               
               
                   
                 7     b ij   k−1  = min{b ij   k ,b ik   k  + b kj   k }; 
               
               
                   
                 8    end 
               
               
                   
                 9   end 
               
               
                   
                 10 end 
               
               
                   
                   
               
             
          
         
       
     
         [0046]    Difference Bound Matrix data structures and Floyd-Warshall algorithm constitute the basic ground for general-purpose Verification and Model Checking tools, which in principle could provide the functionality to solve the scheduling problem. In a preferred embodiment of the present invention, in order to restrain the complexity of state space analysis of a general model, we adapted the algorithmic core of Difference Bound Matrix (DBM) Zones to the identification of timings that result from different orderings between actions subject to mutual exclusion constraints by leveraging on the specific regular structure of scheduled sequences. This resulted in various optimized adaptations of the Floyd-Warshall algorithm which is used to maintain DBM Zones in normal form. 
         [0047]    We then organized DBMs corresponding to different determinations of different sets of mutual exclusion constraints in a graph-theoretical structure allowing incremental traversal. In broad terms, the algorithm iteratively constructs candidate solutions by incrementally adding exclusion constraints and verifying the conformance of temporal parameters with the constraints imposed by the problem; in this setting, DBM structures are used to encode the set of constraints that determine precedences and exclusions among jobs, while Floyd-Warshall algorithm is repeatedly invoked to keep DBMs in their normal form, in order to allow (efficient) manipulation and comparison of obtained results. 
         [0048]    The algorithm initially builds a DBM zone D prec  representing all precedence constraints of the problem. Then, it repeatedly restricts the zone D prec  by adding each exclusion constraint in its two different determinations, and discarding restricted zones that turn out to be empty. Any non-empty zone that has been restricted with some determination for all exclusion constraints comprises a set of solutions for the scheduling problem. This solution univocally corresponds to a single determination of exclusion constraints, and its DBM zone contains all and only the timings that are consistent with these determinations. The optimal timing can be derived in straightforward and direct manner from the normal form of the DBM zone. 
         [0049]    In this perspective, our scheduling algorithm can be regarded as the traversal of a graph-theoretical structure, where each node of the structure is associated with a DBM Zone. In particular: the root is the zone D prec , that satisfies all precedences but no exclusions; intermediate nodes restrict D prec  with a subset of exclusion constraints in some determination; leaf nodes restrict D prec  with all exclusion constraints in some determination; each edge of the graph-theoretical structure represents the choice for some determination of one or more exclusion constraints. Without loss of generality, Alg. 2 reports the algorithm pseudo-code of an implementation according to a preferred embodiment of the present invention: in the specific case each edge adds a single mutual exclusion constraint, so that the graph-theoretical structure becomes a binary tree. 
         [0000]    
       
         
               
             
               
               
             
           
               
                   
               
               
                 ALGORITHM 2: SCHEDULING ALGORITHM 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 1  begin 
               
               
                   
                 2   Dprec = DBM encoding all precedence constraints; 
               
               
                   
                 3   Normalize(Dprec); 
               
               
                   
                 4   Enqueue(&lt;Dprec, 0&gt;); 
               
               
                   
                 5   Dsol = dummy solution; 
               
               
                   
                 6   while QueueIsNotEmpty( ) and time limit is not exceeded do 
               
               
                   
                 7    &lt;D, c&gt; = Dequeue( ); 
               
               
                   
                 8    Dr = AddRightConstraint(D,c + 1); 
               
               
                   
                 9   Normalize(Dr); 
               
               
                   
                 10  if (IsNotEmpty(Dr)) then 
               
               
                   
                 11   if (c + 1 == C) then 
               
               
                   
                 12    Dsol = GetBestSolutiort(D r , Dsol); 
               
               
                   
                 13   else 
               
               
                   
                 14    Enqueue(&lt;D r , c + 1&gt;); 
               
               
                   
                 15  Dl = AddLeftConstraint(D, c + 1); 
               
               
                   
                 16  Normalize(Dl); 
               
               
                   
                 17  if (IsNotEmpty(D l )) then 
               
               
                   
                 18   if (c + 1 == C) then 
               
               
                   
                 19    Dsol = GetBestSolution(Dl , Dsol); 
               
               
                   
                 20   else 
               
               
                   
                 21    Enqueue(&lt;Dl, c + 1&gt;); 
               
               
                   
                 22  end 
               
               
                   
                 23 end 
               
               
                   
                 24 return Dsol ; 
               
               
                   
                 25 end 
               
               
                   
                   
               
             
          
         
       
     
         [0050]    The implementation relies on some data structures and invariants:
       The index c represents the number of exclusion constraints added to the restriction, which also takes the meaning of depth in the binary tree; besides, C is the number of all exclusion constraints, so that a node at depth C is a leaf.   D sol  represents the best identified solution at any point of the computation. It is initialized with some dummy valuation, which might for instance be a trivial and inefficient schedule or some symbolic value that can be regarded as worse than any other one.   The algorithm relies on a queue that contains nodes in the tree, each identified by a DBM zone D and a depth index c.   Exclusion constraints are assumed to be sorted so that an index c uniquely identifies an exclusion constraint.       
 
         [0055]    The algorithm also relies on various functions:
       Normalize(D) applies Floyd-Warshall algorithm to normalize the DBM D taken as parameter;   Enqueue(&lt;D, c&gt;) and Dequeue( ) perform respectively addition and extraction operations on the queue; QueueIsNotEmpty( ) checks if the queue contains at least one element;   AddRightConstraint(D, c) restricts D by adding the constraint of order c in its right determination (i.e. if the constraint is between jobs Ji and Jj, then the added constraint is ri−rj≧ej); in the same manner, AddLeftConstraint(D, c) restricts D by adding the constraint of order c in its left determination (i.e. if the constraint is between jobs Ji and Jj, then the added constraint is r j −r i ≧e i ).   GetBestSolution(D′,D″) receives two DBM zones and selects the one that includes the best timing, e.g. the timing that minimizes the overall completion time for all job sequences.       
 
         [0060]    The above described example is based on the choice between a pair of determinations for each mutual exclusion constraint for each pair of jobs, so that the graph-theoretical structure becomes a binary tree. Those skilled in the art will appreciate that other implementations are possible, e.g. if the algorithm could choose among a set of determinations allowing for the presence of minimum slack delays between conflicting jobs (i.e. constraints of the type r i −r j ≧e j +d ij , where d ij  is a positive quantity), then the resulting structure would become an k-ary tree, with k equal to the number of possible slack delays associated with each mutual exclusion constraint. 
         [0061]    The organization of the algorithm gives space to various possible optimizations that significantly affect performance, outlined in the following. 
       Heuristic Choice Between Exclusion Constraints Determinations 
       [0062]    Depending on the actual implementation of Enqueue(&lt;D, c&gt;) and Dequeue( ) operations, various policies can be enforced for the visit of the search tree (i.e. breadth-first, depth-first, etc.). We adopted a depth-first visit, which is the simplest natural choice as the objective of the algorithm is to find, as soon as possible, a leaf associated to a solution for the scheduling problem. Enumeration queue Q is thus implemented as a LIFO queue. 
         [0063]    Algorithm behaviour is also determined by the order in which D r  and D, are added to the queue (i.e. the order of lines 14 and 21 of Alg. 2): for instance, if D, is always added first, a preorder visit is performed; vice versa, a postorder visit is obtained by always adding D r  first. In our case, a good heuristics should always make choices that most likely lead to feasible solutions without adding too much overhead to the search algorithm. We report here two simple heuristics that influence algorithm behaviour according to the conflicting values of tightness and feasibility:
       feasibility represents the property of exploring first loose configurations that most likely can be extended up to completion without encountering unfeasible constraints. Intuitively, this can be achieved by selecting the exclusion constraint that leaves more degree of freedom to constrained jobs J i  and J j  in the resulting configuration, i.e. the one for which the distance between release times |r i −r j | is free to range in the largest interval (we called it “safe” or “longest-interval” heuristics);   tightness property aims at the opposite goal of obtaining compact schedules, even at the risk of encountering more unfeasible solutions. As a consequence, the resulting heuristics (which we called “greedy” or “shortest-interval” heuristics) selects the constraint that narrows the most the distance between r i  and r j .       
 
       Problem Simplification Through Time Intervals Analysis 
       [0066]    The order given to exclusion constraints impacts on the performance of the algorithm conditioning the early discard of determinations that do not lead to complete schedules. Moreover, some refinements aim to the reduction of the number of conflicting job pairs C; although this kind of optimizations do not impact on the theoretical complexity of the algorithm, in most practical cases they are able to drastically reduce computation time. Some of the choices on which exclusion constraint should be added at each step of the algorithm are already resolved by constraints of D prec ; for instance, if two events requiring the same resource are explicitly constrained to occur in sequence, the exclusion constraint is redundant and it can be ignored. In the specific application this determines a strong speed-up, as jobs represent system actions which are naturally structured along linear sequences. A more subtle but highly effective reduction can be achieved considering those event pairs that, even if not being explicitly constrained, are assured to be non-overlapping by delays configuration. The simplest example of such a condition is constituted by three jobs J i  J j  and J k  such that J i  and J k  are both constrained to start after the completion of J j , and d ij   + +e i &lt;d kj   − . The satisfaction of these three conditions, in fact, guarantees that J i  and J k  are not-overlapping even if no explicit constraint exists between them. More complex interactions are possible, requiring an higher effort to be identified; moreover, as constraints are iteratively added during the search procedure, similar conditions can dynamically occur in the construction of zones D r  and D l  (lines 8 and 15 of Alg. 2), even if they are not verified in the initial D prec  zone; however, a straightforward detection of such conditions can be obtained through direct observation of the normal form of zone D extracted from the queue (line 7 of the Alg. 2). In turn, this can be translated into an optimization of the search algorithm that enables early conflicts detection, preventing the unnecessary construction and normalization of zones D r  and D l . 
       Reducing Complexity Through Warm-Restart Floyd-Warshall Algorithm 
       [0067]    A second optimization manages to reduce the asymptotical time complexity of the procedure, as it operates on the structure of its core engine, the Floyd-Warshall algorithm. 
         [0068]    The standard formulation of Floyd-Warshall algorithm used to find the normal form of DBMs (Alg. 1) runs in time O(N 3 ) with respect to the dimensionality of the zone (in our case, with respect to the number of jobs). However, the normalization invoked at lines 9 and 16 of Alg. 2 presents a peculiarity enabling a specific algorithmic refinement; in fact, zones D r  and D, are obtained adding a single constraint to zone D extracted from the queue at step 7, which is in normal form. This implies that only some of the coefficients of D r  and D l  are affected by the normalization; according to this, Floyd-Warshall algorithm can be rewritten to run in time O(N 2 ) (see Alg. 3). 
         [0000]    
       
         
               
             
               
               
             
           
               
                   
               
               
                 ALGORITHM 3: OPTIMIZED FLOYD-WARSHALL 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 1  begin 
               
               
                   
                 2   for k = 0, . . . , N − 1 do 
               
               
                   
                 3    if b kj  &gt; b ki  + b ij  then 
               
               
                   
                 4     b kj  = b ki  + b ij ; 
               
               
                   
                 5     for l = 0, . . . , N − 1 do 
               
               
                   
                 6      b kl  = min{b kl , b ki  + b ij  + b ji }; 
               
               
                   
                 7     end 
               
               
                   
                 8    end 
               
               
                   
                 9   end 
               
               
                   
                 10  for h = 0, . . . , N − 1 do 
               
               
                   
                 11   b ih  = min{b ih , b ij  + b jh }; 
               
               
                   
                 12  end 
               
               
                   
                 13 end 
               
               
                   
                   
               
             
          
         
       
     
       Problem Simplification Through Events Aggregation 
       [0069]    A third optimization has been introduced to reduce the size of DBM matrix N when the problem comprises jobs with immediate delays, i.e. jobs that are constrained to run immediately one after the other. This can happen if and only if a precedence constraint is set between the two jobs J i  and J j  such that d ij =d ij   + =d ij   − =0. 
         [0070]    In this case jobs J i  and J j  can be merged into a single job J ij  with Release Time r ij =r f  and Execution Time e ij =e i +e j . This results in a DBM matrix reduced by one row and one column. 
       Low-Level Code Optimizations 
       [0071]    The core of the search algorithm in Alg. 2 amounts to just some dozens of lines of code; in order to successfully deploy a production version on the target platform, a lightweight ANSI-C library has been implemented, providing accessory functionalities like input-output management, execution time control and on-line switching among different heuristics, for a total of about 2600 lines of code. Due to the severe hardware limitations, several code optimizations has been adopted in addition to the algorithmic refinements, so as to further reduce space and time cost. The most important among them involves the actual encoding of intermediate solutions. In the pseudo-code of Alg. 2 a queue is used to spool intermediate solutions D as they are generated; although this is the standard approach adopted in search algorithms, in the present case it would easily saturate the available memory. As an alternative solution, the implemented algorithm entirely memorizes only the initial zone D prec  and the actual intermediate zone D, using a stack to encode a compact representation of exclusion constraints chosen until then. Whenever a rollback is due, the algorithm rebuilds a new solution from D prec  according to the informations found into the stack; in this case, since multiple coefficients are modified, a full normalization (i.e. the one of Alg. 1) is required, but this small overhead is largely justified by the strong reduction in space occupation. 
         [0072]      FIG. 5  is a schematic representation of the solution according to a preferred embodiment of the present invention. The search is formulated as a tree traversal which performs a series of successive binary decision on exclusion constraints. The search algorithm aims at finding zones that produce feasible solutions and selecting among such solutions the one corresponding to the optimal solution. The key features of the present method is that of enabling heuristic selection of “reasonable” paths according to predetermined criteria so that a solution can be found without exploring all possible cases. In the present example we consider two possible heuristics: “greedy” (or shortest-interval heuristic) and “safe” (or longest-interval heuristic). The crosses indicates a node in the tree from where no possible solutions can be found, therefore all nodes and leaves departing from such node can be discarded without further calculations. 
         [0073]    The problem simplifications mentioned above reached the results of speeding the computation by reducing complexity. In particular, the time interval analysis minimizes the number of independent constraints; the warm restart reduces complexity of the core step from O(N 3 ) to O(N 2 ); while event aggregation reduces the number of variables. 
         [0074]    Table A and B of  FIG. 6  show the results of a performance test with a comparison of the present method using two different implementations, respectively with a “greedy” heuristic and with a “safe” one. 
         [0075]    Experimental results were measured on a desktop computer (with dual-core 2 GHz processor). This is believed to be around 10 times faster than a real environment on a electromechanical machine as the one where the present algorithm should be run (i.e. an execution on a target platform would take about 10 times more than the same execution on the desktop computer of the present test). 
         [0076]    We ran our search algorithm on a test suite composed by three Analytical Protocols (AP1-AP3), which are realistic examples taken from real-world biological analyses practice: 
         [0000]    AP1: 2 sections, 1 to 3 pipettor actions per slot, 4 read-head actions per section, non-deterministic delays;
 
AP2: 2 sections, 1 to 3 pipettor actions per slot, 4 read-head actions per section, non-deterministic delays;
 
AP3: 2 sections, 1 pipettor action per slot, 4 read-head actions per section, non-deterministic delays;
 
         [0077]    Experiments has been performed by incrementally adding proposed refinements to the search algorithm, which is run adopting both shortest-interval and longest-interval heuristics. Results are shown in Tables A and B of  FIG. 6  respectively; it is particularly worth noting how they compare to the 20-minutes computation required by the Model Checking approach. As expected, successive refinements of the algorithm produce better results: in particular, an higher depth of the search tree gives major relevance to the pruning performed through time interval analysis and to the speed-up of Floyd-Warshall algorithm. Another interesting result comes from the comparison between longest-interval and shortest-interval heuristics; on realistic test cases, the latter has shown better performances than the former, both in performing a full search and in converging to the best solution (on the selected cases, shortest-interval heuristics always finds the best solution as the first one). 
         [0078]    For each tested biological protocol (AP1-AP3) and each optimization level (Basic, TIA, TIA+WR, TIA+WR+EA), performance is characterized as a triple first/best/total where first is the time at which the first solution has been found, best is the time at which the best solution has been found and total is the total execution time. TIA is an acronym for Time Interval Analysis, WR is an acronym for Warm Restart and EA is an acronym for Event Aggregation. On the last two columns, we also report the time span of the first and the best solutions found. 
         [0079]    It will be appreciated that alterations and modifications may be made to the above without departing from the scope of the disclosure. Naturally, in order to satisfy local and specific requirements, a person skilled in the art may apply to the solution described above many modifications and alterations. Particularly, although the present disclosure has been described with a certain degree of particularity with reference to preferred embodiment(s) thereof, it should be understood that various omissions, substitutions and changes in the form and details as well as other embodiments are possible; moreover, it is expressly intended that specific elements described in connection with any disclosed embodiment of the disclosure may be incorporated in any other embodiment as a general matter of design choice. 
         [0080]    For example, similar considerations apply if the components have different structure or include equivalent units. Also, it should be easily understood by those skilled in the art that the present invention may be extended to.