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.

KEYWORDS:

• High-throughput screening
• disturbance decoupling
• uncertain max-plus linear systems
• scheduling

1. 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 max-plus 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, 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.

2. Mathematical preliminaries of max-plus algebra

2.1. Idempotent semirings

Definition 1. A semiring is a set $\mathcal{S}$, equipped with two operations $\oplus$ and $\otimes$, such that $(\mathcal{S},\oplus )$ is a commutative monoid (the zero element will be denoted $\epsilon$), $(\mathcal{S},\otimes )$ is a monoid (the unit element will be denoted $e$), operation $\otimes$ is right and left distributive over $\oplus$, and $\epsilon$ is absorbing for the product (i.e. $\epsilon \otimes a=a\otimes \epsilon =\epsilon ,\mathrm{\forall }a$).

A semiring $\mathcal{S}$ is idempotent if $a\oplus a=a$ for all $a\in \mathcal{S}$. In an idempotent semiring $\mathcal{S}$, operation $\oplus$ induces a partial order relation:

(1) $\begin{array}{ccc}ab& \underset{\_}{\succ }\Leftarrow & a=a\oplus b,\hspace{0.17em}\mathrm{\forall }a,b\in \mathcal{S}.\end{array}$(1)

Then, $a\vee b=a\oplus b$. An idempotent semiring $\mathcal{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 $\mathrm{T}$ (for ‘top’). In this paper, we denote ${\bar{\mathbb{Z}}}_{max}=(\mathbb{Z}\cup \{-\mathrm{\infty },+\mathrm{\infty }\},max,-\mathrm{\infty },+,0)$ as the integer max-plus semiring, where $\epsilon =-\mathrm{\infty }$ is the neutral (zero) element for $max$ and $e=0$ is the unit element for $+$. A non-empty subset $\mathcal{B}$ of a semiring $\mathcal{S}$ is a sub-semiring of $\mathcal{S}$ if for all $a,b\in \mathcal{B}$ we have $a\oplus b\in \mathcal{B}$ and $a\otimes b\in \mathcal{B}$.

Definition 2. A mapping $f:\mathcal{S}\to \mathcal{S}$,where $\mathcal{S}$ is a complete idempotent semiring, is residuated if and only if $f(\epsilon )=\epsilon$ and $f$ is lower-semicontinuous, that is,

$f\left(\underset{i\in I}{\oplus }{a}_i\right)=\underset{i\in I}{\oplus }f\left(a_i\right),$

for any (finite or infinite) set $I$. The mapping $f$ is said to be residuated and $f^{\mathrm{\#}}$ is called its residual.

Theorem 1. (Baccelli et al. 1992)When $f$ is residuated, $f^{\mathrm{\#}}$ is the unique order preserving mapping such that

(2) $f\circ f^{\mathrm{\#}}\underset{\_}{\prec }\mathrm{I}\mathrm{d}and\hspace{0.17em}{f}^{\mathrm{\#}}\circ f\underset{\_}{\succ }\mathrm{I}\mathrm{d},$(2)

where $\mathrm{I}\mathrm{d}$ is the identity mapping from $\mathcal{S}$ to $\mathcal{S}$.

The operator $f\circ f^{\mathrm{\#}}$ denotes function composition, i.e. $f\circ f^{\mathrm{\#}}$ means that $f$ is a function of $f^{\mathrm{\#}}$.

It is straightforward that: $L_a:\mathcal{S}\to \mathcal{S},x\mapsto ax$ and $R_a:\mathcal{S}\to \mathcal{S},x\mapsto xa$ are lower semi-continuous. Therefore, these mappings are residuated i.e. $L_a(x)\underset{\_}{\prec }b$ (resp. $R_a(x)\underset{\_}{\prec }b$) admits a greatest solution, then the following notations are considered:

$\begin{array}{ccc}{L}_a^{\mathrm{♯}}(b)=a\text{invdiameter}b=\oplus \left\{x|axb\right\}\mathrm{a}\mathrm{n}\mathrm{d}& \hspace{0.17em}& \hspace{0.17em}\\ R_a^{\mathrm{♯}}(b)=b\text{diameter}a=\oplus \left\{x|xab\right\},\hspace{0.17em}\mathrm{\forall }a,b\in \mathcal{S},& \hspace{0.17em}& \hspace{0.17em}\end{array}$

where $L_a^{\mathrm{\#}}$ and $R_a^{\mathrm{\#}}(b)$ are the residual mappings.

Theorem 2. (Baccelli et al. 1992)Over a complete max-plus algebra, the implicit equation $x=ax\oplus b$ admits $x=a^{\ast }b$ as the least solution, where $a^{\ast }={\oplus }_{i\in \mathbb{N}}{a}^i$ ($(\cdot {)}^{\ast }$ is the Kleene Star operator) with $a^0=e$.

Definition 3 (Cohen et al. 1996, 1997, 2006). Let $\mathcal{S}$ be a complete idempotent semiring and let $C$ be a $n\times p$ matrix with entries in $\mathcal{S}$. We call null kernel of $C$ as the set of elements $x\in {\mathcal{S}}^p$ such that $Cx=\epsilon$, denoted as $\mathrm{k}\mathrm{e}\mathrm{r}C$. We call equivalence kernel of $L_C$ (denoted by ${ker}_{\mathrm{e}\mathrm{q}}C$), the subset of all pairs of elements of ${\mathcal{S}}^p$ whose components are both mapped by $L_C$ to the same element in ${\mathcal{S}}^n$,i.e. the following definition

(3) $\begin{array}{ccc}{ker}_{\mathrm{e}\mathrm{q}}C& :=& \left\{\left(s,s{ }^{\mathrm{\prime }}\right)\in {\left({\mathcal{S}}^p\right)}^2|Cs=Cs{ }^{\mathrm{\prime }}\right\}.\end{array}$(3)

Clearly ${ker}_{\mathrm{e}\mathrm{q}}C$, is an equivalence relation on ${\mathcal{S}}^p$, $i.e.$, $Cs=Cs{ }^{\mathrm{\prime }}\iff s{ }^{\mathrm{\prime }}\equiv s(\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{k}\mathrm{e}{\mathrm{r}}_{\mathrm{e}\mathrm{q}}C)$ and furthermore it is a congruence and then we can define the quotient ${\mathcal{S}}^p/{ker}_{eq}C$.

The subset of elements $s{ }^{\mathrm{\prime }}\in {\mathcal{S}}^p$ that are equivalent to $s$ modulo ${ker}_{eq}C$ is denoted as $[s{]}_C$, $i.e.$,

(4) $\begin{array}{ccc}{s}_C=\left\{s{ }^{\mathrm{\prime }}\in {\mathcal{S}}^p|s{ }^{\mathrm{\prime }}\equiv s(\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{k}\mathrm{e}{\mathrm{r}}_{\mathrm{e}\mathrm{q}}C)\right\}\subset {\mathcal{S}}^p.& \hspace{0.17em}& \hspace{0.17em}\end{array}$(4)

Definition 4 (Restricted map). Let $f:{\mathcal{S}}^p\to {\mathcal{S}}^n$ be a map and $\mathcal{A}\subseteq {\mathcal{S}}^p$. We will denote $f_{|\mathcal{A}}:\mathcal{A}\to {\mathcal{S}}^n$ the map defined by $f_{|\mathcal{A}}=f\circ \mathrm{I}{\mathrm{d}}_{|\mathcal{A}}$ where $\mathrm{I}{\mathrm{d}}_{|\mathcal{A}}:\mathcal{A}\to {\mathcal{S}}^p,x\mapsto x$ be the canonical injection. Identically, let $\mathcal{B}\subseteq {\mathcal{S}}^n$ with $\mathrm{I}\mathrm{m}f\subseteq \mathcal{B}$. Map ${ }_{\mathcal{B}|}f:{\mathcal{S}}^p\to \mathcal{B}$ is defined by $f=\mathrm{I}{\mathrm{d}}_{|\mathcal{B}}{\circ }_{\mathcal{B}|}f$,where $\mathrm{I}{\mathrm{d}}_{|\mathcal{B}}:\mathcal{B}\to {\mathcal{S}}^n,x\mapsto x$ be the canonical injection.

Definition 5 (Isotone map). A map $f:{\mathcal{S}}^p\to {\mathcal{S}}^p$ is said to be order preserving or isotone if the following property holds: $a\underset{\_}{\prec }b\Rightarrow f(a)\underset{\_}{\prec }f(b).$

Definition 6 (Closure map). An isotone map $f:{\mathcal{S}}^p\to {\mathcal{S}}^p$ is a closure map if $f\underset{\_}{\succ }\mathrm{I}{\mathrm{d}}_{{\mathcal{S}}^p}$ and $f\circ f=f$.

Proposition 1 (2003. A closure map $f:{\mathcal{S}}^p\to {\mathcal{S}}^p$ restricted to its image ${ }_{\mathrm{I}\mathrm{m}f|}f$ is a residuated map whose residual is the canonical injection $\mathrm{I}{\mathrm{d}}_{|\mathrm{I}\mathrm{m}f}:\mathrm{I}\mathrm{m}f\to {\mathcal{S}}^p$, $s\mapsto s$.

Corollary 1. Let $\mathcal{K}:{\mathcal{S}}^p\to {\mathcal{S}}^p,s\mapsto s^{\ast }$ be a map, where $s^{\ast }={\oplus }_{i\in \mathbb{N}}{s}^i$. The map ${ }_{\mathrm{I}\mathrm{m}\mathcal{K}|}\mathcal{K}$ is a residuated map whose residual is ${\left({ }_{\mathrm{I}\mathrm{m}\mathcal{K}|}\mathcal{K}\right)}^{\mathrm{\#}}=\mathrm{I}{\mathrm{d}}_{|\mathrm{I}\mathrm{m}\mathcal{K}}$. This means that $x=s^{\ast }$ is the greatest solution to inequality $x^{\ast }\underset{\_}{\prec }{s}^{\ast }$. 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,

(5) $a^{\ast }(b{a}^{\ast }{)}^{\ast }=(a\oplus b{)}^{\ast }=(a^{\ast }b{)}^{\ast }{a}^{\ast },$(5)

(6) $(ab{)}^{\ast }a=a(ba{)}^{\ast },$(6)

(7) $a(a\text{invdiameter}(ax))=ax,$(7)

(8) $(ab)\text{invdiameter}x=b\text{invdiameter}(a\circ \mathrm{\setminus }x).$(8)

2.2. Max-plus linear systems

A max-plus linear system is defined by the following equations:

$x(k)=Ax(k-1)\oplus Bu(k)$

(9) $y(k)=Cx(k),$(9)

where $x(k)\in {\bar{\mathbb{Z}}}_{max}^n$, $u(k)\in {\bar{\mathbb{Z}}}_{max}^p$, $y(k)\in {\bar{\mathbb{Z}}}_{max}^q$ and $k\in \mathbb{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)\in {\bar{\mathbb{Z}}}_{max}$ and predict the system evolution. For a state equation in (9), each increasing sequence $\left\{x(k)\right\}$, it is possible to define the transformation $X(\gamma )=\underset{k\in \mathbb{Z}}{\oplus }x(k){\gamma }^k$ where $\gamma$ is a backward shift operator in event domain (i.e. $Y(\gamma )=\gamma X(\gamma )\iff \left\{y(k)\right\}=\left\{x(k-1)\right\}$, (see, 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(\gamma )$ is a synthetic representation of the trajectory $x(k)$. The set of the formal power series in $\gamma$ is denoted by ${\bar{\mathbb{Z}}}_{max}[[\gamma ]]$ and constitutes an idempotent semiring. Therefore, the state equation in (9) becomes a polynomial equation or an event-domain representation,

$X(\gamma )=\bar{A}X(\gamma )\oplus BU(\gamma ),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\bar{A}=\gamma A,$

(10) $Y(\gamma )=CX(\gamma ),$(10)

where the state $X(\gamma )\in {\left({\bar{\mathbb{Z}}}_{max}[[\gamma ]]\right)}^n$, the output $Y(\gamma )\in {\left({\bar{\mathbb{Z}}}_{max}[[\gamma ]]\right)}^q$, and the input $U(\gamma )\in {\left({\bar{\mathbb{Z}}}_{max}[[\gamma ]]\right)}^p$, and matrices $\bar{A}\stackrel{\mathrm{\Delta }}{=}\gamma A\in {\left({\bar{\mathbb{Z}}}_{max}[[\gamma ]]\right)}^{n\times n}$, $B\in {\left({\bar{\mathbb{Z}}}_{max}[[\gamma ]]\right)}^{n\times p}$ and $C\in {\left({\bar{\mathbb{Z}}}_{max}[[\gamma ]]\right)}^{q\times n}$ represent the link between transitions. According to the state equation (10)(10) $Y(\gamma )=CX(\gamma ),$(10) , the evolution of the system is

$X(\gamma )={\bar{A}}^{\ast }BU(\gamma ),$

(11) $Y(\gamma )=C{\bar{A}}^{\ast }BU(\gamma ).$(11)

The trajectories $U(\gamma )$ and $Y(\gamma )$ can be related (Baccelli et al. (1992), p. 243) by the equation $Y(\gamma )=H(\gamma )U(\gamma )$, where $H(\gamma )=C{\bar{A}}^{\ast }B\in {\left({\bar{\mathbb{Z}}}_{max}[[\gamma ]]\right)}^{q\times p}$ is called the transfer matrix of the TEG. Entries of matrix $H$ are periodic series (Baccelli et al. (1992), p. 260) in the idempotent semiring, usually represented by $p(\gamma )\oplus q(\gamma )(\tau {\gamma }^{\nu }{)}^{\ast }$, where $p(\gamma )$ is a polynomial representing the transient behavior, $q(\gamma )$ is a polynomial corresponding to a pattern which is repeated periodically, the period being given by the monomial $(\tau {\gamma }^{\nu })$. For example, consider a transfer matrix $H(\gamma )=0{\gamma }^0\oplus 1\gamma \oplus 3{\gamma }^4\oplus (5{\gamma }^5\oplus 6{\gamma }^7)(3{\gamma }^4{)}^{\ast }$ represented in Figure 1. The transient behavior is given by $p(\gamma )=0{\gamma }^0\oplus 1\gamma \oplus 3{\gamma }^4$ while the periodic behavior is given by $q(\gamma )(\tau {\gamma }^{\nu }{)}^{\ast }=(5{\gamma }^5\oplus 6{\gamma }^7)(3{\gamma }^4{)}^{\ast }$.

Figure 1. Diagram showing a series with transient and periodic behavior (Hardouin et al., 2010).

3. 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.

Figure 2. Layout of an HTS system. Source: HTSC at Washington University in St. Louis.

Table 1. Resources and their associated tasks at the HTS plant at Washington University in St. Louis.

Resource	Device Type	Tasks
Cytomat2C	Microplate hotel	Feeds microplates into the system.
		Stores microplates after they are processed.
Teleshake	Incubator (4-shelf microplate shaker)	Allows the biochemical substances to incubate in a controlled environment.
Envision	Reader	Analyzes chemical changes occurring in biochemical substances.
SCARA	Robot arm (transportation device)	Transfers microplates between resources.
Multidrop-384	Liquid dispenser	Dispenses liquid substance into the microplates’ wells.

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 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$.

Figure 3. Diagram showing the start and release times of each activity of the assay process.

Figure 4. Gantt chart of a single batch of the assay process.

4. 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,

Figure 5. TEG of the assay process.

$x_2(k)=x_1(k)\otimes 32=32x_1(k),$

$x_{12}(k)=5\otimes x_{11}(k)\oplus 0\otimes x_9(k)\oplus 0\otimes x_{13}(k-4)$

$=5x_{11}(k)\oplus x_9(k)\oplus x_{13}(k-4).$

In ${\bar{\mathbb{Z}}}_{max}[[\gamma ]]$ domain, these equations become

$X_2(\gamma )=32X_1(\gamma ),$

$X_{12}(\gamma )=5X_{11}(\gamma )\oplus X_9(\gamma )\oplus {\gamma }^4{X}_{13}(\gamma ).$

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 $\epsilon$ are

$\begin{array}{cccccccc}{A}_{1,22}=\gamma & A_{2,1}=32& A_{3,2}=19& A_{4,3}=9& A_{5,4}=23& A_{6,5}=54& \hspace{0.17em}& \hspace{0.17em}\\ A_{6,8}=e& A_{7,10}=\gamma & A_{8,6}=18& A_{8,7}=1& A_{9,8}=18& A_{9,12}=e& \hspace{0.17em}& \hspace{0.17em}\\ A_{10,9}=1& A_{11,14}=\gamma & A_{12,11}=5& A_{12,9}=e& A_{12,13}={\gamma }^4& A_{13,12}=210& \hspace{0.17em}& \hspace{0.17em}\\ A_{13,16}=e& A_{14,13}=5& A_{15,18}=\gamma & A_{16,15}=1& A_{16,13}=e& A_{17,16}=18& \hspace{0.17em}& \hspace{0.17em}\\ A_{18,17}=1& A_{19,17}=e& A_{20,19}=54& A_{21,20}=20& A_{22,21}=34,& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}$

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 $\epsilon$ are

$B_{1,1}=B_{7,2}=B_{11,3}=B_{15,4}=C_{1,22}=e,$

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. (2013, 2016a).

4.1. 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

Figure 6. TEG of the assay process with disturbances.

$X(\gamma )=\bar{A}X(\gamma )\oplus BU(\gamma )\oplus SQ(\gamma ),$

$Y(\gamma )=CX(\gamma ),$

where matrices $A,B$ and $C$ are as earlier defined and the non- $\epsilon$ elements of the disturbance matrix are:

$S_{2,1}=S_{4,2}=S_{8,3}=S_{16,4}=S_{20,5}=e.$

Definition 7. (Shang, Hardouin, Lhommeau, Maia et al. 2016a) The system (12) is called disturbance decoupled by a state-feedback control $u(k)=Fx(k-1)\oplus v(k)$) (or an open-loop control $u(k)=Pv(k)$) if and only if any disturbance signal will not affect the system output $y(k)$ for all $k\in \mathbb{Z}$ and for any initial condition $x_0$.

Solving the DDP means that the control $U(\gamma )$ has to achieve

(12) $C{\bar{A}}^{\ast }BU(\gamma )\oplus C{\bar{A}}^{\ast }SQ(\gamma )=C{\bar{A}}^{\ast }BU(\gamma )$(12)

(13) $\iff C{\bar{A}}^{\ast }SQ(\gamma )\underset{\_}{\prec }C{\bar{A}}^{\ast }BU(\gamma ).$(13)

The control input is computed from an external input, $V(\gamma )$, and the disturbance, $Q(\gamma )$, which is measurable, as $U(\gamma )=PQ(\gamma )\oplus V(\gamma )$. For any external input $V(\gamma )$ and disturbance $Q(\gamma )$, the inequality above is equivalent to

(14) $C{\bar{A}}^{\ast }SQ(\gamma )C{\bar{A}}^{\ast }BPQ(\gamma )\oplus C{\bar{A}}^{\ast }BV(\gamma ).$(14)

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)\oplus 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(\gamma )$ 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(\gamma )$, such that the following equation holds,

(15) $C{\bar{A}}^{\ast }BU(\gamma )\oplus C{\bar{A}}^{\ast }SQ(\gamma )=C{\bar{A}}^{\ast }SQ(\gamma )$(15)

(16) $\iff C{\bar{A}}^{\ast }BU(\gamma )C{\bar{A}}^{\ast }SQ(\gamma ).$(16)

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(\gamma )=PQ(\gamma )$, 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

(17) $C{\bar{A}}^{\ast }BPQ(\gamma )\oplus C{\bar{A}}^{\ast }SQ(\gamma )=C{\bar{A}}^{\ast }SQ(\gamma ),\mathrm{\forall }Q(\gamma ).$(17)

The optimum prefilter matrix $P_{opt}$ that solves both the DDP and MDDP is given by

(18) $P_{opt}=\left(C{\bar{A}}^{\ast }B\right)\text{invdiameter}\left(C{\bar{A}}^{\ast }S\right).$(18)

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

$\begin{array}{ccc}{F}_{opt}& =& \left(C{\bar{A}}^{\ast }B\right)\text{invdiameter}\left(C{\bar{A}}^{\ast }\tilde{B}\right)\text{invdiameter}\left(C{\bar{A}}^{\ast }\tilde{B}\right)\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\tilde{B}=\left[B|S\right].\end{array}$

When simulated in MinMaxGD, a C++ toolbox for handling periodic series introduced in Cottenceau et al. (2000), we obtain

$P_{opt}=\left[\begin{array}{ccccc}459\gamma {(491\gamma )}^{\ast }& 431\gamma {(491\gamma )}^{\ast }& 354\gamma {(491\gamma )}^{\ast }& 126\gamma {(491\gamma )}^{\ast }& 54\gamma {(491\gamma )}^{\ast }\\ 104(491\gamma {)}^{\ast }& 76(491\gamma {)}^{\ast }& 490\gamma {(491\gamma )}^{\ast }& 262\gamma {(491\gamma )}^{\ast }& 190\gamma {(491\gamma )}^{\ast }\\ 118(491\gamma {)}^{\ast }& 90(491\gamma {)}^{\ast }& 13(491\gamma {)}^{\ast }& 276\gamma {(491\gamma )}^{\ast }& 204\gamma {(491\gamma )}^{\ast }\\ 332(491\gamma {)}^{\ast }& 304(491\gamma {)}^{\ast }& 227(491\gamma {)}^{\ast }& 490\gamma {(491\gamma )}^{\ast }& 418\gamma {(491\gamma )}^{\ast }\end{array}\right],$

and

$F_{opt}=\left[\begin{array}{c}\gamma {(491\gamma )}^{\ast }\\ 136\gamma {(491\gamma )}^{\ast }\\ 150\gamma {(491\gamma )}^{\ast }\\ 364\gamma {(491\gamma )}^{\ast }\end{array}\right].$

Each element in $P_{opt}$ and $F_{opt}$ is periodic with a period of $491$ time units. For instance,

$Popt(2,1)=104(491\gamma {)}^{\ast }$

$=104(e\oplus 491\gamma \oplus 982{\gamma }^2\oplus 1473{\gamma }^3\oplus 1964{\gamma }^4\oplus \cdots )$

$=104\oplus 595\gamma \oplus 1086{\gamma }^2\oplus 1577{\gamma }^3\oplus 2068{\gamma }^4\oplus \cdots .$

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

Figure 7. TEG model of the HTS system with integrated prefilter and the output controllers.

$V=P_{opt}Q,\mathrm{a}\mathrm{n}\mathrm{d}$

$U=\hat{F}CX\oplus V=\hat{F}Y\oplus V.$

For instance, the first external input, $V_1$, is gotten from the first row of $P_{opt}$ and the disturbance vector as follows:

$V_1=\left[\begin{array}{ccccc}459\gamma {(491\gamma )}^{\ast }& 431\gamma {(491\gamma )}^{\ast }& 354\gamma {(491\gamma )}^{\ast }& 126\gamma {(491\gamma )}^{\ast }& 54\gamma {(491\gamma )}^{\ast }\end{array}\right]\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ q_4\\ q_5\end{array}\right];$

and the first control input $U_1$ is

$U_1=\gamma (491\gamma {)}^{\ast }y\oplus V_1.$

In $\gamma (491\gamma {)}^{\ast }$, the cycle time is 491, and $\gamma$ represents a token in the cyclic mode.

5. 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 max-plus linear system.

An uncertain max-plus linear system is defined as

$x(k)=Ax(k-1)\oplus Bu(k)\oplus Sq(k),$

(19) $y(k)=Cx(k),$(19)

where system matrices lie in corresponding matrix intervals with known lower and upper bounds, specifically, $A\in \mathbf{A}=[A_l,A_u]\in I({\bar{\mathbb{Z}}}_{max}{)}^{n\times n}$, $B\in \mathbf{B}=[B_l,B_u]\in I({\bar{\mathbb{Z}}}_{max}{)}^{n\times p}$, $S\in \mathbf{S}=[S_l,S_u]\in I({\bar{\mathbb{Z}}}_{max}{)}^{n\times r}$, and $C\in \mathbf{C}=[C_l,C_u]\in I({\bar{\mathbb{Z}}}_{max}{)}^{q\times n}$. And the states are $x(k)\in X\cong {\bar{\mathbb{Z}}}_{max}^n$, the inputs are $u(k)\in U\cong {\bar{\mathbb{Z}}}_{max}^p$, the disturbances are $q(k)\in Q\cong {\bar{\mathbb{Z}}}_{max}^r$, and the outputs are $y(k)\in Y\cong {\bar{\mathbb{Z}}}_{max}^q$ and $k\in \mathbb{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.

5.1. 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 $\mathbf{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

(20) $\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{B}\mathbf{U}\bar{\oplus }\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{S}\mathbf{Q}=\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{B}\mathbf{U},$(20)

holds for intervals of matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\mathbf{S}$,and the disturbance interval $\mathbf{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 $\mathbf{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

(21) $\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{B}\mathbf{U}\bar{\oplus }\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{S}\mathbf{Q}=\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{S}\mathbf{Q},$(21)

is satisfied for intervals of matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\mathbf{S}$,and the disturbance interval $\mathbf{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 $\mathbf{U}$ which induces an output interval unchanged by the disturbance $Q(\gamma )$ in an interval $\mathbf{Q}$.

5.2. Open-loop controller interval solving WMDDP

Solving the WMDDP by an open-loop controller requires finding an optimal prefilter interval ${\mathbf{P}}_{\mathbf{o}\mathbf{p}\mathbf{t}}^{\mathbf{W}\mathbf{M}\mathbf{D}\mathbf{D}\mathbf{P}}=[P_{lopt},P_{uopt}]$ such that the control interval $\mathbf{U}=\mathbf{P}\mathbf{Q}$ solves the WMDDP. $\mathbf{U}=[U_l,U_u]$, $\mathbf{P}=[P_l,P_u]$ and $\mathbf{Q}$ is interval of the disturbances. ${\mathbf{P}}_{\mathbf{o}\mathbf{p}\mathbf{t}}^{\mathbf{W}\mathbf{M}\mathbf{D}\mathbf{D}\mathbf{P}}$ ensures the following equality holds

(22) $\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{B}{\mathbf{P}}_{\mathbf{o}\mathbf{p}\mathbf{t}}^{\mathbf{W}\mathbf{M}\mathbf{D}\mathbf{D}\mathbf{P}}\bar{\oplus }\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{S}=\mathbf{C}{\bar{\mathbf{A}}}^{\mathbf{\ast }}\mathbf{S},$(22)

$\iff \mathbf{C}{\bar{\mathbf{A}}}^{\ast }{\mathbf{B}\mathbf{P}}_{\mathbf{o}\mathbf{p}\mathbf{t}}^{\mathbf{W}\mathbf{M}\mathbf{D}\mathbf{D}\mathbf{P}}\underset{\_}{\prec }\mathbf{C}{\bar{\mathbf{A}}}^{\ast }\mathbf{S}.$

The optimal prefilter interval was found in Shang, Hardouin, Lhommeau, Maia et al. (2016a) as

(23) $P_{opt}^{WMDDP}\stackrel{\Delta }{=}C{\overline{A}}^{*}B\overline{\varphi }C{\overline{A}}^{*}S=\left[P_{lopt}^{WMDDP},P_{uopt}^{WMDDP}\right],$(23)

where

$P_{lopt}^{WMDDP}=C_l{\bar{A}}_l^{\ast }{B}_l\text{invdiameter}{C}_l{\bar{A}}_l^{\ast }{S}_l\wedge C_u{\bar{A}}_u^{\ast }{B}_u\text{invdiameter}{C}_u{\bar{A}}_u^{\ast }{S}_u,$

(24) $P_{uopt}^{WMDDP}=C_u{\bar{A}}_u^{\ast }{B}_u\text{invdiameter}{C}_u{\bar{A}}_u^{\ast }{S}_u.$(24)

5.3. 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 $\mathbf{U}=\mathbf{F}\mathbf{Y}\bar{\oplus }\mathbf{V}=\mathbf{F}\mathbf{C}\mathbf{X}\bar{\oplus }\mathbf{V}$. If we define $\tilde{H}=C{\bar{A}}^{\ast }\tilde{B}$, the controller $F$ has to satisfy

Figure 8. The output feedback structure.

(25) $H\left(e\epsilon \right)F*H\underset{\_}{\prec }H,$(25)

for all $\tilde{H}\in \tilde{\mathbf{H}}=[C_l{\bar{A}}_l^{\ast }{\tilde{B}}_l,C_u{\bar{A}}_u^{\ast }{\tilde{B}}_u]$, ${\tilde{B}}_l=\left[B_l|S_l\right]$ and ${\tilde{B}}_l=\left[B_u|S_u\right]$.

The interval of output feedback controls $\mathbf{F}=[F_l,F_u]$ is

$\stackrel{\mathbf{ˆ}}{\mathbf{F}}=\mathbf{H}\overline{\text{invdiameter}}\tilde{\mathbf{H}}\overline{\text{diameter}}\tilde{\mathbf{H}}=\left[{\hat{F}}_{lopt},{\hat{F}}_{uopt}\right]$

(26) $=\left[H_l\text{invdiameter}{\tilde{H}}_l\text{diameter}{\tilde{H}}_l\wedge H_u\text{invdiameter}{\tilde{H}}_u\text{diameter}{\tilde{H}}_u,H_u\text{invdiameter}{\tilde{H}}_u\text{diameter}{\tilde{H}}_u\right].$(26)

6. 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.

Figure 9. TEG of the assay process of the uncertain case with disturbances.

The controllers given in Eqs (24) and (26) are computed for the HTS system. The optimal lower and upper bounds of the prefilter are

$P_{lopt}^{WMDDP}=\left[\begin{array}{ccccc}459\gamma {(491\gamma )}^{\ast }& 431\gamma {(491\gamma )}^{\ast }& 354\gamma {(491\gamma )}^{\ast }& 126\gamma {(491\gamma )}^{\ast }& 54\gamma {(491\gamma )}^{\ast }\\ 104(491\gamma {)}^{\ast }& 76(491\gamma {)}^{\ast }& 490\gamma {(491\gamma )}^{\ast }& 262\gamma {(491\gamma )}^{\ast }& 190\gamma {(491\gamma )}^{\ast }\\ 118(491\gamma {)}^{\ast }& 90(491\gamma {)}^{\ast }& 13(491\gamma {)}^{\ast }& 276\gamma {(491\gamma )}^{\ast }& 204\gamma {(491\gamma )}^{\ast }\\ 332(491\gamma {)}^{\ast }& 304(491\gamma {)}^{\ast }& 227(491\gamma {)}^{\ast }& 490\gamma {(491\gamma )}^{\ast }& 418\gamma {(491\gamma )}^{\ast }\end{array}\right],$

and

$P_{uopt}^{WMDDP}=\left[\begin{array}{ccccc}484\gamma {(516\gamma )}^{\ast }& 456\gamma {(516\gamma )}^{\ast }& 379\gamma {(516\gamma )}^{\ast }& 126\gamma {(516\gamma )}^{\ast }& 54\gamma {(516\gamma )}^{\ast }\\ 104(516\gamma {)}^{\ast }& 76(516\gamma {)}^{\ast }& 515\gamma {(516\gamma )}^{\ast }& 262\gamma {(516\gamma )}^{\ast }& 190\gamma {(516\gamma )}^{\ast }\\ 118(516\gamma {)}^{\ast }& 90(516\gamma {)}^{\ast }& 13(516\gamma {)}^{\ast }& 276\gamma {(516\gamma )}^{\ast }& 204\gamma {(516\gamma )}^{\ast }\\ 357(516\gamma {)}^{\ast }& 329(516\gamma {)}^{\ast }& 252(516\gamma {)}^{\ast }& 515\gamma {(516\gamma )}^{\ast }& 443\gamma {(516\gamma )}^{\ast }\end{array}\right].$

The optimal lower and upper bounds of the output feedback controller are

${\hat{F}}_{lopt}=\left[\begin{array}{c}\gamma {(491\gamma )}^{\ast }\\ 136\gamma {(491\gamma )}^{\ast }\\ 150\gamma {(491\gamma )}^{\ast }\\ 364\gamma {(491\gamma )}^{\ast }\end{array}\right],$

and

${\hat{F}}_{uopt}=\left[\begin{array}{c}\gamma {(516\gamma )}^{\ast }\\ 136\gamma {(516\gamma )}^{\ast }\\ 150\gamma {(516\gamma )}^{\ast }\\ 389\gamma {(516\gamma )}^{\ast }\end{array}\right].$

In Figure 10, the upper-bound prefilter, $P_{uopt}^{WMDDP}$ and the upper-bound output feedback controller, ${\hat{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.

Figure 10. TEG model of the HTS system with integrated upper bound prefilter and the output controllers.

Figure 11. Gantt charts of the process with disturbance and with disturbance and prefilter.

Table 2. Activities and corresponding states.

Activity	States
1	$x_1-x_6$
2	$x_7-x_{10}$
3	$x_{11}-x_{14}$
4	$x_{15}-x_{18}$
5	$x_{19}-x_{22}$

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\mathrm{\%}$ improvement. Therefore, the prefilter ensures a just-in-time completion time.

Table 3. Comparison of activity duration without prefilter and with prefilter for $k=1$.

	Without prefilter			With prefilter
Activity	Start time (s)	Release time (s)	Duration (s)	Start time (s)	Release time (s)	Duration (s)
1	484	621	137	484	621	137
2	124	640	516	620	640	20
3	363	879	516	634	879	245
4	377	893	516	873	893	20
5	892	1000	108	892	1000	108

Table 4. Waiting time eliminated with prefilter.

Activity	Waiting time eliminated	$\mathrm{\%}$ Time saved
1	0	0
2	496	96
3	275	52
4	496	96
5	0	0

7. 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.

8. 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.

Acknowledgments

The authors would like to thank Dr Ma Xenia Ilagan, Associate Professor in the Department of Biochemistry and Molecular Biophysics, for providing access to the data from the High Throughput Screening Center at Washington University in St. Louis.

Disclosure statement

No potential conflict of interest was reported by the author(s).