Scheduling and control of high throughput screening systems with uncertainties and disturbances

ABSTRACT Based on the existing research work on the modeling and control of cyclically operated high-throughput screening systems, this paper presents the scheduling and control of high-throughput screening systems with uncertainties and disturbances. Different definitions of disturbance decoupling problems are considered in order to achieve the optimal scheduling of high-throughput screening systems. Online scheduling was achieved by computing optimal prefilters and output feedback controllers, which can be computed using residuation theory. The control strategies achieved by the controllers are optimal in the sense of the just-in-time criterion commonly used in industrial schedule practice.


Introduction
A successful drug discovery is an extremely time-consuming procedure, including initial target identification and validation, pre-clinical trials on animals, regulatory approval to start trials in humans, clinical trials, submission of marketing and manufacturing authorization, licensing review, product sale, and post-marketing surveillance (Major, 1998;Mayer et al., 2008;Mayr & Fuerst, 2008;Noah, 2010;Pereira & Williams, 2007). With the development of robotics and high-speed computing technology, it is feasible to develop automatic systems that can screen a large number of biochemical compounds in a short period. Such an automatic compound screening and analyzing process is called high-throughput screening (HTS) in drug discovery of pharmaceutical industries. HTS is a standard technology routinely employed in the pharmaceutical industry for drug discovery processes. It is used for initial screening in the process of drug discovery to reduce what is an almost infinite number of possible combinations of compounds to a reasonably few enough possibilities on which further testing can be carried out. Many HTS systems consist of several activities on several different resources. An HTS operation can incorporate multiple batches, each with hundreds of events, where a batch is a combination of all the operations to be performed on a set of substances for complete analysis. HTS provides a practical and efficient method to test a large number of synthetic compounds in miniaturized in vitro assays to identify hit targets of interest. Then, the chemical compounds that have therapeutic and useful pharmacological or biological activities, called leads, are evaluated and undergo lead optimization to identify promising lead compounds. Followed by the initial synthesis and animal testing in preclinical trials and three phases of clinical trials on humans, a drug can be put on the market after the Food and Drug Administration (FDA) approval.
It is necessary to develop a schedule for an HTS system that ensures that analysis is completed in the shortest possible time, i.e. increase the throughput or reduce the makespan (makespan is the completion time of the last activity of a process). The optimal schedule must also ensure consistency in time spent on each activity for every batch. This means that the operation of the HTS system must be cyclic. There have been various methods developed to obtain an optimal schedule for an HTS system and other cyclic systems. The problem of obtaining a time-optimal schedule for an HTS system, as an example of a strictly cyclic system, was formulated and solved as mixed integer linear programming (MILP) problem by Mayer and Raisch (2004). This MILP approach is useful for obtaining a static schedule, which works well so long as there are no changes to the cyclic design. A time-optimal schedule was also obtained by Mayer and Raisch (2003) for cyclic systems with multicapacity resources. HTS processes with hierarchical nested cycles were considered by Mayer et al. (2008), where mixed integer optimization was adopted to minimize the overall cycle time. In these existing works, the schedules obtained are static as they do not tolerate any deviation from the cyclic processes.
Starting from a static schedule, a max-plus algebra model of an HTS system was obtained in Brunsch and Raisch (2009). In a max-plus linear system, the addition and multiplication operations in the classical control systems are replaced by max and addition, respectively. The real number is also replaced by a special algebra called maxplus algebra with the operations changed accordingly (Baccelli et al., 1992;Golan, 1999). The scheduling problem can be transformed into the problem of computing a feedback controller that generates possible control actions in the presence of disturbances for such a max-plus linear system model. The controller generates a control input that gives an optimal schedule in spite of system malfunctions and delays. Brunsch et al. (2012) presented an optimal feedback controller. There have been investigations in this vein into scheduling and control of dioid models of discrete event systems in the presence of disturbances (Shang et al. (2013(Shang et al. ( , p. 2016a), all of which can be applied to HTS systems. In Shang, et al. (2016b), the concept of disturbance decoupling problem (DDP) for uncertain max-plus linear systems whose system matrices vary between intervals was introduced.
In this paper, an HTS process is first modeled in max-plus algebra and a static schedule is developed for the process based on the model. Online scheduling was then considered by computing controllers that ensure an optimal schedule when there are disturbances. The result is extended to the uncertain case where system matrices lie in intervals. The rest of this paper is organized as follows. Section 2 gives the mathematical fundamentals of the max-plus algebra that is to be used in describing the HTS system. In Section 3, the specifications of the HTS system considered are described. Sections 4 and 5 detail online scheduling of certain and uncertain HTS systems, respectively. Simulation results are presented in Section 6 and a discussion of the result and contribution of the study is presented in Section 7. Section 8 concludes this paper with future research directions.

Idempotent semirings
Definition 1. A semiring is a set S, equipped with two operations � and � , such that ðS; �Þ is a commutative monoid (the zero element will be denoted ε), ðS; �Þ is a monoid (the unit element will be denoted e), operation � is right and left distributive over � , and ε is absorbing for the product (i.e. ε � a ¼ a � ε ¼ ε; "a).
A semiring S is idempotent if a � a ¼ a for all a 2 S. In an idempotent semiring S, operation � induces a partial order relation: (1) An idempotent semiring S is complete if sums of infinite numbers of terms are always defined, and if multiplication distributes over infinite sums too. In particular, the sum of all the elements of the idempotent semiring is denoted T (for 'top').
In this paper, we denote Z max ¼ ðZ [ fÀ 1; þ1g; max; À 1; þ; 0Þ as the integer maxplus semiring, where ε ¼ À 1 is the neutral (zero) element for max and e ¼ 0 is the unit element for þ . A non-empty subset B of a semiring S is a sub-semiring of S if for all a; b 2 B we have a � b 2 B and a � b 2 B.
for any (finite or infinite) set I. The mapping f is said to be residuated and f # is called its residual.
Theorem 1. (Baccelli et al. 1992)When f is residuated, f # is the unique order preserving mapping such that (2) where Id is the identity mapping from S to S.
The operator f � f # denotes function composition, i.e. f � f # means that f is a function of f # .
It is straightforward that: L a : S ! S; x7 !ax and R a : S ! S; x7 !xa are lower semicontinuous. Therefore, these mappings are residuated i.e. L a ðxÞ � b (resp. R a ðxÞ � b) admits a greatest solution, then the following notations are considered: where L # a and R # a ðbÞ are the residual mappings.
Theorem 2. (Baccelli et al. 1992)Over a complete max-plus algebra, the implicit equation Definition 3 (Cohen et al. 1996(Cohen et al. , 1997(Cohen et al. , 2006. Let S be a complete idempotent semiring and let C be a n � p matrix with entries in S. We call null kernel of C as the set of elements x 2 S p such that Cx ¼ ε, denoted as ker C. We call equivalence kernel of L C (denoted by ker eq C), the subset of all pairs of elements of S p whose components are both mapped by L C to the same element in S n ,i.e. the following definition Clearly ker eq C, is an equivalence relation on S p , i:e:, Cs ¼ Cs 0 , s 0 ;sðmodker eq CÞ and furthermore it is a congruence and then we can define the quotient S p = ker eq C.
The subset of elements s 0 2 S p that are equivalent to s modulo ker eq C is denoted as ½s� C , i:e:, Definition 4 (Restricted map). Let f : S p ! S n be a map and A � S p . We will denote f jA : A ! S n the map defined by Definition 5 (Isotone map). A map f : S p ! S p is said to be order preserving or isotone if the following property holds: a � b ) f ðaÞ � f ðbÞ: Actually, the greatest solution achieves equality.
There are basic properties for star and residuation operations in the residuation theory (Baccelli et al. 1992), for example,

Max-plus linear systems
A max-plus linear system is defined by the following equations: where xðkÞ 2 Z n max , uðkÞ 2 Z p max , yðkÞ 2 Z q max and k 2 Z. This kind of system makes it possible to describe the behaviors of TEGs, by associating to each transition a firing date sequence x i ðkÞ 2 Z max and predict the system evolution. For a state equation in (9), each increasing sequence xðkÞ f g, it is possible to define the transformation Baccelli et al. (1992), p. 228). This transformation is analogous to the z-transform used in discrete-time classical control theory and the formal series XðγÞ is a synthetic representation of the trajectory xðkÞ. The set of the formal power series in γ is denoted by Z max ½½γ�� and constitutes an idempotent semiring. Therefore, the state equation in (9) becomes a polynomial equation or an event-domain representation, where the state XðγÞ 2 Z max ½½γ�� À � n , the output YðγÞ 2 Z max ½½γ�� À � q , and the input � q�n represent the link between transitions. According to the state equation (10), the evolution of the system is The trajectories UðγÞ and YðγÞ can be related (Baccelli et al. (1992) (Baccelli et al. (1992), p. 260) in the idempotent semiring, usually represented by pðγÞ � qðγÞðτγ ν Þ � , where pðγÞ is a polynomial representing the transient behavior, qðγÞ is a polynomial corresponding to a pattern which is repeated periodically, the period being given by the monomial ðτγ ν Þ.

Modeling of HTS systems
We consider an HTS system from the High-Throughput Screening Center (HTSC) at Washington University in St. Louis whose layout is shown in Figure 2. The system consists of different resources as identified in Figure 2 on which different activities are  carried out. Table 1 lists the resources used in the assay and their associated tasks.
Cytomat2C is the microplate hotel where all the microplates are kept and where they are returned after all the processes are completed. The microplates have either 384 wells or 96 wells each. Teleshake is the station where the incubation occurs while reading is done in Envision. SCARA is the robot doing all of the transfer activities between resources while necessary biochemical substances are added to the microplates in the Multidrop.
An enzymatic assay process that consists of 11 activities was carried out on the HTS system of Figure 2. To develop a schedule for the process, we are interested in the start and release times of each activity of the process. Figure 3 shows the start and release times of each of the activities. o i is the start time of activity i and r i its release time. It should be noted that for this particular process, the start time of activity i corresponds to the release time of activity ði À 1Þ. It can be seen that the 11 activities of the process are carried out on 5 resources. Activities 1 and 11 are carried out on the Cytomat2C; activities 2; 4; 6; 8; and 10 are transfer activities carried out by the SCARA robot; activities 5 and 9 are done on the Envision. The Multidrop and the Teleshake are single activity resources responsible for activity 3 and acticvity 7, respectively. All the resources are single capacity resources, except for the Teleshake, which has a capacity of 4. he assay process utilized full 384 -well microplates. If the microplates were not full or 96 -well microplates were used, the duration of activities 5 and 9, which are done on the Envision, might be reduced. For this process, the duration of other activities that include transfer, liquid dispensing and incubaTtion would remain unchanged because they are independent of the number of wells. The Gantt chart for one batch of the assay process is shown in Figure 4. We see that the makespan of a single process is 491 s. It can also be seen from the Gantt chart that there is some overlap between some Transfers microplates between resources.

Multidrop-384
Liquid dispenser Dispenses liquid substance into the microplates' wells. activities during a batch in which case a batch occupies more than one resource at the same time. Such is the case during transfer activities 6 and 8 and during the incubation, activity 7.

Scheduling of HTS systems by max-plus algebra
A timed event graph (TEG), shown in Figure 5 was generated to capture the information provided in Figure 3. The TEG has 12 transitions representing the start and release events of the different activities. The 4 tokens in the TEG indicate that the Teleshake has a capacity of 4. Each transition can be described by max-plus algebra. For instance, x 2 ðkÞ ¼ x 1 ðkÞ � 32 ¼ 32x 1 ðkÞ; x 12 ðkÞ ¼ 5 � x 11 ðkÞ � 0 � x 9 ðkÞ � 0 � x 13 ðk À 4Þ ¼ 5x 11 ðkÞ � x 9 ðkÞ � x 13 ðk À 4Þ: In Z max ½½γ�� domain, these equations become X 2 ðγÞ ¼ 32X 1 ðγÞ; X 12 ðγÞ ¼ 5X 11 ðγÞ � X 9 ðγÞ � γ 4 X 13 ðγÞ:  The HTS system can be modeled as a max-plus linear system as in Eq. (10). The elements of matrix A are the states of the system and are given by the minimal time of each activity as shown in Figure 5. Matrices B and C are the input and output matrices respectively. The output signal is the release time of the last activity. The elements of matrix A that are not ε are where A i;j gives the element at row i and column j in matrix A.
The elements of matrices B and C that are not ε are where B i;j and C i;j give the element at row i and column j in matrices B and C. It is trivial to obtain an off-line schedule for the system by computing its output trajectory as given in Eq. (11). To obtain an on-line schedule in the presence of disturbances, we need to find an optimal controller that solves the disturbance decoupling problem as explained in Shang et al. (2013Shang et al. ( , 2016a.

Online scheduling for HTS system
Disturbances are unforeseen disruptions to the normal evolution of an HTS system where a disturbance can be a breakdown of the robot arm or an obstruction in its way. Disturbances are usually modelled as uncontrollable inputs in a system. For our running example, we assume that there is a disruption in the movement of the robot arm. These disturbances are shown as inputs q 1 to q 5 in Figure 6. The system equations become XðγÞ ¼ AXðγÞ � BUðγÞ � SQðγÞ; YðγÞ ¼ CXðγÞ; where matrices A; B and C are as earlier defined and the non-ε elements of the disturbance matrix are: Solving the DDP means that the control UðγÞ has to achieve , CA � SQðγÞ � CA � BUðγÞ: The control input is computed from an external input, VðγÞ, and the disturbance, QðγÞ, which is measurable, as UðγÞ ¼ PQðγÞ � VðγÞ. For any external input VðγÞ and disturbance QðγÞ, the inequality above is equivalent to Definition 8. (Shang, Hardouin, Lhommeau, Maia et al. 2016a) The max-plus linear system described in (12) is called modified disturbance decoupled by a state feedback control uðkÞ ¼ Fxðk À 1Þ � vðkÞ (or an open-loop control uðkÞ ¼ PvðkÞ) if and only if the system output signals will not be disturbed more than the output signals influenced by the disturbances.
The objective of the modified disturbance decoupling problem (MDDP) is to find the greatest open-loop or output feedback control UðγÞ such that the output trajectories will not be delayed more than the disturbance signals have acted on the system. This means finding the greatest control, UðγÞ, such that the following equation holds, , CA � BUðγÞCA � SQðγÞ: The scheduling problem can be transformed into a DDP or MDDP in max-plus linear systems. This is done by finding a prefilter or feedback controller that generates a control input to solve the DDP or MDDP. Such a control input generates an optimal schedule in spite of system malfunctions and delays, based on the just-in-time criterion according to the industry standard. If the control UðγÞ ¼ PQðγÞ, where P is a prefilter which generates the control by taking the measured disturbances into account, then, solving the MDDP is equivalent to finding a prefilter P satisfying CA � BPQðγÞ � CA � SQðγÞ ¼ CA � SQðγÞ; "QðγÞ: The optimum prefilter matrix P opt that solves both the DDP and MDDP is given by According to Shang et al.(2016a), the optimal feedback controller F opt that generates a schedule according to the just-in-time criterion is given by When simulated in MinMaxGD, a C++ toolbox for handling periodic series introduced in Cottenceau et al. (2000), we obtain Each element in P opt and F opt is periodic with a period of 491 time units. For instance, A realization of the TEG together with the controllers is shown in Figure 7. It can be seen that the disturbances pass through the prefilter before going into the system. We see that For instance, the first external input, V 1 , is gotten from the first row of P opt and the disturbance vector as follows: ; and the first control input U 1 is In γð491γÞ � , the cycle time is 491, and γ represents a token in the cyclic mode.

Online scheduling of uncertain HTS systems
For the HTS process we have been considering, the time delays indicated in Figure 5 and Figure 6 give the minimum amount of time to be spent on each activity. In some cases, certain activities have a lower limit and an upper limit of the time they take rather than a minimal time. In such a situation where the times for some activities are not fixed, we have what we call an uncertain HTS system, which can be described by uncertain maxplus linear system. An uncertain max-plus linear system is defined as xðkÞ ¼ Axðk À 1Þ � BuðkÞ � SqðkÞ; yðkÞ ¼ CxðkÞ; where system matrices lie in corresponding matrix intervals with known lower and upper bounds, specifically, A 2 A ¼ ½A l ; A u � 2 IðZ max Þ n�n , B 2 B ¼ ½B l ; B u � 2 IðZ max Þ n�p , S 2 S ¼ ½S l ; S u � 2 IðZ max Þ n�r , and C 2 C ¼ ½C l ; C u � 2 IðZ max Þ q�n . And the states are xðkÞ 2 X ffi Z n max , the inputs are uðkÞ 2 U ffi Z p max , the disturbances are qðkÞ 2 Q ffi Z r max , and the outputs are yðkÞ 2 Y ffi Z q max and k 2 Z. An uncertain HTS system with disturbances can be disturbance decoupled by different control inputs in much the same way as for systems without uncertainties. However, the controls have to be selected very restrictively in order to solve the DDP and the MDDP for all the variations of the system matrices, as well as for arbitrary disturbances. Sometimes, such non-trivial controls do not even exist for all variations of system matrices. Therefore, we relax the definitions of the DDP and MDDP in order to achieve more general results.

Weakly DDP and weakly MDDP
Definition 9. (WDDP) The uncertain max-plus linear system in (19) is called weakly disturbance decoupled if and only if there exists a control U belonging to an interval U such that the output trajectory interval, generated by the disturbances and the uncertain matrices, belongs to the output trajectory interval, generated by the control U. In other words, solving WDDP means that the following equality holds for intervals of matrices A, B, C, S,and the disturbance interval Q.

Definition 10. (WMDDP) The max-plus linear system described in (19) is called weakly modified disturbance decoupled if and only if it exists a control U belonging to an interval
U such that the output trajectories belong to the interval generated by any disturbances and uncertain matrices. In other words, solving WMDDP means that the following equality is satisfied for intervals of matrices A, B, C, S,and the disturbance interval Q.
The controllers solving the DDP and MDDP have to solve them for any arbitrary variations of system matrices, as well as the disturbances. On the other hand, the WDDP and WMDDP only requires that there exists an interval of control U which induces an output interval unchanged by the disturbance QðγÞ in an interval Q.

Open-loop controller interval solving WMDDP
Solving the WMDDP by an open-loop controller requires finding an optimal prefilter interval P WMDDP opt ¼ ½P lopt ; P uopt � such that the control interval U ¼ PQ solves the WMDDP. U ¼ ½U l ; U u �, P ¼ ½P l ; P u � and Q is interval of the disturbances. P WMDDP opt ensures the following equality holds , CA � BP WMDDP opt � CA � S: The optimal prefilter interval was found in Shang, Hardouin, Lhommeau, Maia et al. (2016a) as where

Feedback controller interval solving WMDDP
As introduced in Shang, Hardouin, Lhommeau, Maia et al. (2016a), the output feedback controller structure solving the WMDDP is shown in Figure 8. The output feedback control interval that solves the WMDDP is for all e H 2 e The interval of output feedback controls

Simulation results
The uncertain case of the HTS assay process is considered to simulate controllers that preserve open-loop behavior in the presence of disturbance. The incubation time in the uncertain HTS process has an upper and a lower bound; it takes up to 210 s and no more than 235 s as shown in Figure 9. The controllers given in Eqs (24) and (26) are computed for the HTS system. The optimal lower and upper bounds of the prefilter are  In Figure 10, the upper-bound prefilter, P WMDDP uopt and the upper-bound output feedback controller, F uopt are integrated into the TEG of the HTS process. The Gantt chart of the process is shown in Figure 11. The states corresponding to the activities marked in the Gantt chart are given in Table 2. It is seen from the Gantt chart that the prefilter eliminates unnecessary wait time caused by the disturbance. With the prefilter, the completion time of each activity satisfies the just-in-time criterion. This is better seen in Tables 3 and 4. Table 3 shows the start times, release times and durations of each activity with and without disturbance and when the prefilter is included. We see, for example, that activity 2 has a duration of 516 seconds without prefilter. This is reduced to just 20 s, which is the minimum required time, by the prefilter. Table 4 shows that this is a 96% improvement. Therefore, the prefilter ensures a just-in-time completion time.

Discussion
This study has successfully obtained a max-plus algebraic model for an HTS system that is operated cyclically. We take further the method introduced in Brunsch and Raisch (2009) for obtaining an online schedule for cyclically operated HTS systems with disturbances by obtaining a schedule for an HTS system with disturbances and uncertainties. This study also presented optimal controllers that solve the DDP, and thereby solving the scheduling problem for an HTS system with disturbances and uncertainties. Whereas in Shang et al. (2016b), controllers that solve DDP for uncertain max-plus linear systems whose system matrices vary between intervals was introduced, we show that the controllers solving DDP also solve the scheduling problem in cyclically operated discrete event systems. The supervisory control proposed in this study will be useful for generating static schedules for HTS systems and other discrete event systems. More importantly, it can be used to generate online schedule, which is important for when there are disturbances or deviations from a predetermined schedule during operation.
In summation, an HTS process is first modeled in max-plus algebra and a static schedule is developed for the process based on the model. Online scheduling was then considered by computing controllers that ensure an optimal schedule when there are disturbances. The result is extended to the uncertain case where system matrices lie in intervals. This paper  presents the solution of the scheduling problem for HTS systems, with and without uncertainties, by solving different types of disturbance decoupling problems of max-plus linear systems. Open-loop controls and the output feedback controls solving the DDP, which ultimately solves an online scheduling problem, were constructed using the residuation the ory and proved to be optimal schedules for the just-in-time control criterion.

Future research
In future research, the method employed here will be extended to other high-throughput processes such as High-Throughput in vitro Absorption, Distribution, Metabolism and Excretion (HT-ADME) and High-Throughput Mass Spectrometry (HT-MS) by modeling these processes as a max-plus linear systems. Of particular interest in future work will be HTS systems with nested cycles. The authors are working on other optimization models such as mixed integer programming and queueing theory to obtain additional results and to compare with the results obtained here. x 1 À x 6 2 x 7 À x 10 3 x 11 À x 14 4 x 15 À x 18 5 x 19 À x 22