# Workload Balancing Among Heathcare Workers Under Uncertain Service Time
Using Distributionally Robust Optimization

Duy Anh Nguyen

###### Abstract.

Healthcare systems are facing serious challenges in balancing their human
resources to cope with volatile service demand, while at the same time
providing necessary job satisfaction to the healthcare workers. We propose in
this paper a distributionally robust optimization formulation to generate a
task assignment plan that promotes the fairness in allocation, attained by
reducing the difference in the total working time among workers, under
uncertain service time. The proposed joint chance constraint model is
conservatively approximated by a worst-case Conditional Value-at-Risk, and we
devise a sequential algorithm to solve the finite-dimensional reformulations
which are linear (mixed-binary) optimization problems. We also provide
explicit formula in the situation where the support set of the random vectors
is a hyperrectangle. The experiment with synthetic data suggests promising
results for our approach.

Keywords. Fairness in healthcare; task allocations; stochastic service time;
distributionally robust optimization; joint chance constraint.

The author is with the Hanoi Obstetrics and Gynecology Hospital, Vietnam.

##  1\. Introduction

The fast aging population and devastating pandemic outbreaks are exerting
extreme pressure on the worldwide healthcare systems. Apart from the lack of
medical devices, prescriptive drugs and the lack of developed infrastructure,
healthcare providers are in critical shortage of qualified workers, including
surgeons, medical doctors and nurses [1]. Many surveys indicate that a
majority of healthcare workers are not satisfied with their occupation due to
low salary, tremendous emotional stress and long shift working hours [12, 18,
32]. Because the healthcare workers hold direct responsibility to communicate
and provide service to patient, ensuring their job satisfaction is an
important factor to delivery the highest possible quality of treatment and
care to the community in need [9, 36, 24]. Hence, modern hospital management
needs to deal with numerous complex and conflicting criteria in the daily
operations [15].

Apart from the provision of necessary support and attractive compensation, a
reasonable work schedule is arguably one of the most critical factor to
promote job satisfaction in healthcare workers [10, 20]. Besides the number of
working hours per day/week and the characteristics of the work rotation, the
current research is focusing on promoting fairness in the scheduling or the
rostering of the workers. While there are many definitions for the notion of
fairness in the existing literature [11, 27, 33], in this paper, we focus on a
specific choice of fairness that ensure a balance workload among health
workers. Workload balancing has also been studied in the literature related to
scheduling, where the focus has been placed on the static setting [17, 5, 2].

In this paper, we use the duration of the task as a proxy to calculate the
workload of the tasks, and as such, other factors which may directly
contribute to the workload such as the complexity of the task are omitted.
Balancing the workload in this case is simplified to generating a shift
schedule in which all healthcare workers will spend a relatively similar
amount of time to finish the assigned jobs. In reality, the exact duration of
each task, also referred to as the service time, is rarely known in advanced,
and thus, it is difficult to correctly specify the true workload given to each
worker. If we use random vectors to model the service time, then in many
situations, we may not have access to the true distribution of these random
vectors. At best, we may have access to limited training or historical data,
from which we can try to infer the ambiguity set which contains the true
distribution with high confidence. Consequentially, workload balancing under
uncertainty remains a challenging and active area for both theoretical
research as well as operational implementation [19].

There are several attractive frameworks to deal with the inherent uncertain
nature of healthcare decision making at both the operational and the planning
level including stochastic programming [30] and robust optimization [3].
Distributionally robust optimization is an emerging framework for decision
support when the information about the distribution of the underlying
randomness is incorrect or incomplete [7, 14]. Distributionally robust
solutions have been shown to deliver superior performance in many healthcare
decision making tasks, including, but not limited to, appointment scheduling
[37, 13], nurse staffing [26], surgery block allocation [34], and location
planning of emergency service [16].

This paper aims to bridge the existing gap in the literature of fair
scheduling of healthcare workers by considering a flexible model to solve the
workload balancing under uncertain service time using the distributionally
robust optimization framework. This paper stems from one of the long-standing
efforts of the management team at our hospital111Hospital name is anonymized
to improve the scheduling and rostering of our healthcare workers, especially
in improving their job satisfaction. The contributions of this paper are
highlighted as follows.

  * •

We detail a joint chance constraint optimization problem that aims to produce
a fair allocation of tasks among a team of healthcare workers. This
optimization problem requires that the difference of the working time among
healthcare workers is bounded by a threshold with high probability. Based on
this formulation, we propose the distributionally robust optimization problem
with joint chance constraint to deal with situations when the true underlying
probability of the randomness is difficult to be identified, or when this
probability may entail the risk of misspecification.

  * •

We develop a conservative approximation of the joint chance constraint using
the worst-case Conditional Value-at-Risk when the ambiguity set contains the
expected values and support information of the true distribution. Furthermore,
we propose a sequential optimization algorithm to efficiently solve the
resulting reformulations, which are often linear (binary) optimization
problems. To facilitate practical implementation of the approach, we also
tailor the result to the special case where the random vector is supported on
a hyperrectangular set.

This paper unfolds as follows. Section 2 layouts our mathematical framework
for workload balancing under uncertain service time. Section 3 details the
approximation scheme, which will be tailored to the hyperrectangular support
set in Section 4. Section 5 reports the numerical results and Section 6
concludes the paper with future research directions.

Notations. At any given risk level $\varepsilon\in(0,1)$, the Conditional
Value-at-Risk of a measurable loss function $\ell(\xi)$ dependent on the
random vector $\xi$ whose distribution is governed by a probability measure
$\mathbb{Q}$ is defined as
$\mathbb{Q}\text{-CVaR}_{\varepsilon}(\ell(\xi))\triangleq\inf_{\tau}\left\\{\tau+\varepsilon^{-1}\mathds{E}_{\mathbb{Q}}[\max\\{0,\ell(\xi)-\tau\\}]\right\\}$,
where the operator $\mathds{E}_{\mathbb{Q}}$ denotes the integration under the
probability measure $\mathbb{Q}$. We use $|\mathcal{A}|$ to denote the
cardinality of the set $\mathcal{A}$. For any matrix (vector) $A$, we denote
by $A^{\top}$ its transpose. For two vectors $l$ and $u$ of the same
dimension, the inequalities $l\leq u$ and $l<u$ are understood in the element-
wise sense. For two integers $a\leq b$, we use $\llbracket a,b\rrbracket$ to
denote the set $\\{a,a+1,\ldots,b\\}$.

##  2\. Model Formulations

We consider a healthcare system with a collection $\mathcal{J}$ of healthcare
workers. Given a set of tasks $\mathcal{I}$, the planner aims to find the best
plan to assign each task $i\in\mathcal{I}$ to a worker $j\in\mathcal{J}$ to
maximize the system output, while at the same time ensures that the total
workload assigned to each worker $j\in\mathcal{J}$ is comparatively balanced.
We assume that task $i\in\mathcal{I}$ utilizes a random amount of resource
$\xi_{i}$, and assigning task $i\in\mathcal{I}$ to worker $j\in\mathcal{J}$
brings a deterministic reward $r_{ij}\in\mathbb{R}$. If we use the binary
decision variables $x_{ij}$ with

$x_{ij}=\begin{cases}1&\text{if task $i$ is assigned to worker $j$,}\\\
0&\text{otherwise.}\end{cases}$

to carry the information of the task assignment to each worker, then the total
time that takes worker $j$ to finish all given tasks can be written as the sum
$\sum_{i\in\mathcal{I}}\xi_{i}x_{ij}$. To promote the fairness in the
allocation of tasks to workers, it is plausible to impose that all workers
will require exactly the same total amount of time to finish the given tasks,
which can be equivalently translated to a mathematical constraint of the form

$\sum_{i\in\mathcal{I}}\xi_{i}x_{ij}=\sum_{i\in\mathcal{I}}\xi_{i}x_{ij^{\prime}}\quad\forall(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}.$

Unfortunately, because $\xi$ are random service time, the above constraint can
hardly be met in many settings. Indeed, if $\xi$ is governed by a probability
measure which is absolutely continuous with respect to a Lebesgue measure then
this constraint is met with probability 0 irrespective of any choice of
assignment $x$. A possible remedy is to constraint the time difference to be
lower than a positive threshold $\Delta\in\mathbb{R}_{+}$ as

$\Big{|}\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})\Big{|}\leq\Delta\quad\forall(j,j^{\prime})\in\mathcal{J}\times\mathcal{J},$

where the inequality should be understood as an almost-surely inequality,
which basically means it should hold with probability one. Unfortunately, this
almost-surely constraint is extremely pessimistic and oftentimes it leads to
the infeasibility of the optimization problem. Instead of imposing this
constraint, we seek to propose a relaxation using the joint chance constraint
reformulation. Thus, we propose the stochastic fairness-constrained allocation
problem which can be written as a joint chance constrained optimization
problem

(1)
$\begin{array}[]{cl}\max&f(x)\triangleq\displaystyle\sum_{(i,j)\in\mathcal{I}\times\mathcal{J}}r_{ij}x_{ij}\\\
\mathrm{s.t.}&x\in\mathbb{X}\\\
&\displaystyle\mathbb{P}\left(\Big{|}\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})\Big{|}\leq\Delta\;\forall(j,j^{\prime})\right)\geq
1-\varepsilon.\end{array}$

In problem (1), the set of feasible assignment $\mathbb{X}$ is defined as

$\mathbb{X}\triangleq\left\\{x\in\mathbb{\\{}0,1\\}^{|\mathcal{I}|\times|\mathcal{J}|}:\sum_{j}x_{ij}=1\;\forall
i\in\mathcal{I}\right\\},$

where the constraints defining $\mathbb{X}$ indicate that any task $i$ is
assigned to exactly one worker. The objective function $f(x)$ quantifies the
total reward of each allocation $x$ and it is assumed to be an affine function
of $x$. The joint chance constraint in (1) is prescribed using a balance
threshold $\Delta\in\mathbb{R}_{+}$ and a small tolerance
$\varepsilon\in(0,1)$. The sign $\forall(j,j^{\prime})$ in this constraint is
understood as for all pairs of workers
$(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}$. This joint chance constraint
dictates that the comparative workload among all workers is balanced up to a
threshold $\Delta$ with high probability $1-\varepsilon$.

Problem (1) is a stochastic program which requires the true probability
distribution $\mathbb{P}$ of the (joint) distribution of the random resource
utilization $\xi$ as input. Unfortunately, evaluating whether an allocation
$x$ is feasible for the joint chance constraint of (1) is in general #P-hard
[8]. Moreover, due to the complex nature of the tasks in the healthcare
system, full knowledge about $\mathbb{P}$ is rarely known in practice. At
best, the healthcare system may have access to some historical data which can
be used to infer $\mathbb{P}$. However, the available data may be very scarce
and may not represent well the distribution of certain task.

To alleviate the downsides of problem (1), we propose to approach problem (1)
using the lens of distributionally robust optimization. To this end, we define
the ambiguity set

(2)
$\mathcal{Q}\triangleq\left\\{\mathbb{Q}\in\mathcal{P}(\mathbb{R}^{|\mathcal{I}|}):\mathds{E}_{\mathbb{Q}}[\xi]=\mu,~{}\mathbb{Q}(\xi\in\Xi)=1\right\\},$

where $\mathcal{P}(\mathbb{R}^{|\mathcal{I}|})$ is the set of all probability
measures on $\mathbb{R}^{|\mathcal{I}|}$ and the support set $\Xi$ is a
compact polyhedron (polytope) that can be described effectively by $M$ linear
constraints of the form

(3) $\Xi=\left\\{\xi\in\mathbb{R}^{|\mathcal{I}|}:G\xi\leq h\right\\}$

for some matrix $G\in\mathbb{R}^{M\times|\mathcal{I}|}$ and vector
$h\in\mathbb{R}^{M}$. Furthermore, we assume that the mean vector $\mu$
belongs to the interior of the support $\Xi$, that means
$\mu\in\mathrm{int}(\Xi)$, or equivalently $G\mu<h$. Descriptively speaking,
the ambiguity set $\mathcal{Q}$ contains all probability measures under which
$\xi$ has a given mean vector $\mu$ and the support of $\Xi$ is contained in
the polytope $\Xi$. Using this ambiguity set $\mathcal{Q}$, we define the set
of distributionally robust fair allocations

$\displaystyle\mathcal{X}\triangleq\left\\{\begin{array}[]{l}x\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{J}|}:\\\
\min\limits_{\mathbb{Q}\in\mathcal{Q}}\displaystyle\mathbb{Q}\left(\Big{|}\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})\Big{|}\leq\Delta\;\forall(j,j^{\prime})\right)\geq
1-\varepsilon\end{array}\right\\}$

that contains allocations $x$ such that the workload among workers are
comparatively balanced for all probability measures $\mathbb{Q}$ in
$\mathcal{Q}$. Instead of solving the stochastic fairness-constrained
allocation problem (1), we propose to solve the following distributionally
robust fairness-constrained allocation problem

(4) $\max\left\\{f(x):x\in\mathbb{X}\cap\mathcal{X}\right\\},$

which can also be written in a more explicit form as

$\begin{array}[]{cl}\max&f(x)\\\ \mathrm{s.t.}&x\in\mathbb{X}\\\
&\min\limits_{\mathbb{Q}\in\mathcal{Q}}\displaystyle\mathbb{Q}\left(\Big{|}\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})\Big{|}\leq\Delta\;\forall(j,j^{\prime})\right)\geq
1-\varepsilon.\end{array}$

Before moving on to Section 3 where we provide a solution procedure to problem
(4), we would like to first discuss the components of (4) in further details.
While the set of allocations $\mathbb{X}$ has been defined in a minimalistic
fashion, we emphasize that there is a great flexibility to incorporate
additional operational constraints into $\mathbb{X}$ to promote a higher level
of fairness in the assignment. These constraints may include, but are not
restricted to, the following:

  * •

skillset restrictions: if a worker $j$ does not possess the skill required by
task $i$, then the linear constraint $x_{ij}=0$ can be added into
$\mathbb{X}$,

  * •

workers’ preferences: if a worker $j$ is not willing to perform task $i$, then
the linear constraint $x_{ij}=0$ can be added into $\mathbb{X}$,

  * •

pre-assignment: if task $i$ has to assigned to worker $j$, then the linear
constraint $x_{ij}=1$ can be added into $\mathbb{X}$,

  * •

knapsack constraints: if worker $j$ can perform at most $\theta_{j}$ tasks,
then the linear constraints $\sum_{i\in\mathcal{I}}x_{ij}\leq\theta_{j}$ can
be added into $\mathbb{X}$.

Similarly, there are various penalty terms that can be incorporated into the
loss function $f$ to reflect real-life incurred costs. However, these
modifications of the feasible set $\mathbb{X}$ and the objective function $f$
does not alter the reformulation or approximation of the distributionally
robust fair allocation set $\mathcal{X}$. Because our main goal in this paper
is to explore efficient methods to embody the set $\mathcal{X}$ into the
optimization phase, we thus keep $\mathbb{X}$ and $f$ as simple as possible to
avoid any unwanted effects on the quality of the fair allocation.

Returning to the joint chance constraint in (4), we would like to emphasize
that the choice of the joint chance constraint is imperatively preferable
compared to a set of individual chance constraints. Indeed, consider the
following set of individual chance constraints

(5)
$\left.\begin{array}[]{l}\displaystyle\mathbb{P}\left(\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})\leq\Delta\right)\geq
1-\epsilon_{jj^{\prime}}^{1}\\\
\displaystyle\mathbb{P}\left(\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij^{\prime}}-x_{ij})\leq\Delta\right)\geq
1-\epsilon_{jj^{\prime}}^{2}\end{array}\right\\}\quad\forall(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}$

for a collection of parameters
$(\epsilon_{jj^{\prime}}^{1},\epsilon_{jj^{\prime}}^{2})_{(jj^{\prime})\in\mathcal{J}\times\mathcal{J}}$.
Theoretically, an allocation $x$ which is feasible for the above set of
individual chance constraint will satisfy

$\mathbb{P}\left(\Big{|}\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})\Big{|}\leq\Delta\;\forall(j,j^{\prime})\right)\geq
1-\epsilon$

with $\epsilon=\sum(\epsilon_{jj^{\prime}}^{1}+\epsilon_{jj^{\prime}}^{2})$ by
the Bonferroni inequality [25, Section 6.1]. This implies that an allocation
that satisfies the set of individual chance constraint may not theoretically
deliver the necessary guarantee of fairness jointly over all pairs of workers
with high probability. Consequentially, to achieve a fixed tolerance
$\varepsilon$ in the joint chance constraint formulation (1), one may resort
to solve an approximation using individual chance constraints of the form (5)
with the parameters $\epsilon_{jj^{\prime}}^{1}$ and
$\epsilon_{jj^{\prime}}^{2}$ being set to a sufficiently small value,
typically by setting to a common value
$\epsilon_{jj^{\prime}}^{1}=\epsilon_{jj^{\prime}}^{2}=\varepsilon/(2|\mathcal{J}|^{2})$.
Unfortunately, for real-life applications with a large number of workers, this
choice of individual parameters tends to be over-conservative and may, in
specific cases, fail to even identify a feasible allocation to problem (1).
The incapability of the Bonferroni approach to deliver sufficient performance
guarantee for the joint chance constraint in the distributionally robust
optimization setting has been numerically studied in [23]. This argument
justifies the merit of embracing the joint chance constraint formulation in
(1) and in its distributionally robust counterpart (4). Finally, one can also
augment the threshold $\Delta$ to pair-specific threshold
$\Delta_{jj^{\prime}}$ for each $(j,j^{\prime})$ to capture more subtle
specification of the system. This parameters augmentation will not alter the
approach in this paper, thus, we keep a common threshold $\Delta$ for the sake
of simplicity.

##  3\. Tractable Approximations

Problem (4) is intuitively appealing, however, it is also notoriously
difficult to be solved. Indeed, evaluating whether an allocation $x$ belongs
to the set $\mathcal{X}$ involves solving an infinite dimensional optimization
problem over the space of probability measures. Leveraging on diverse results
from the field of distributionally robust optimization, we derive in this
section a safe and tractable approximation for problem (4). To this end, we
notice that the argument in the joint chance constraint

$\Big{|}\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})\Big{|}\leq\Delta\quad\forall(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}$

can be written explicitly as a collection of $2\times|\mathcal{J}|$ linear
constraints

$\left.\begin{array}[]{r}\displaystyle\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})-\Delta\leq
0\\\
\displaystyle-\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})-\Delta\leq
0\end{array}\right\\}\quad\forall(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}.$

Adopting a similar strategy as in [6], we associate the first and second set
of constraints with a set of strictly positive scaling factors
$\alpha_{jj^{\prime}}$ and $\beta_{jj^{\prime}}$ respectively, and the above
collection of linear constraints is equivalent to

$\left.\begin{array}[]{r}\displaystyle\alpha_{jj^{\prime}}\Big{(}\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})-\Delta\Big{)}\leq
0\\\
\displaystyle\beta_{jj^{\prime}}\Big{(}-\sum_{i\in\mathcal{I}}\xi_{i}(x_{ij}-x_{ij^{\prime}})-\Delta\Big{)}\leq
0\end{array}\right\\}\;\forall(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}.$

Thus, for any strictly positive scaling matrices
$\alpha,\beta\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$, the set
of probabilistically fair allocations $\mathcal{X}$ can be written as

$\displaystyle\mathcal{X}=\left\\{\begin{array}[]{l}x\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{J}|}:\\\
\inf\limits_{\mathbb{Q}\in\mathcal{Q}}\displaystyle\mathbb{Q}\big{(}\max_{k\in\llbracket
1,K\rrbracket}\\{\xi^{\top}a_{k}(x)+b_{k}\\}\leq 0\big{)}\geq
1-\varepsilon\end{array}\right\\},$

where $K=2|\mathcal{J}|^{2}$ is the total number of individual constraints in
the joint chance constraint. For any $k\in\llbracket 1,K\rrbracket$, the
$i$-th component of the vector $a_{k}(x)\in\mathbb{R}^{|\mathcal{I}|}$
dependent on $x$ is defined as

(6)
$\displaystyle[a_{k}(x)]_{i}=\begin{cases}\alpha_{jj^{\prime}}(x_{ij}-x_{ij^{\prime}})&\text{if
}k=(j-1)|\mathcal{J}|+j^{\prime},\\\
\beta_{jj^{\prime}}(x_{ij^{\prime}}-x_{ij})&\text{if
}k=|\mathcal{J}|^{2}+(j-1)|\mathcal{J}|+j^{\prime},\end{cases}$

for any $i\in\mathcal{I}$, and the scalar $b_{k}\in\mathbb{R}$ is defined as

(7) $b_{k}=\begin{cases}-\alpha_{jj^{\prime}}\Delta&\text{if
}k=(j-1)|\mathcal{J}|+j^{\prime},\\\ -\beta_{jj^{\prime}}\Delta&\text{if
}k=|\mathcal{J}|^{2}+(j-1)|\mathcal{J}|+j^{\prime}.\end{cases}$

We emphasize that the vectors $a_{k}$ and the scalars $b_{k}$ are dependent on
the specific choice of the scaling factors $\alpha$ and $\beta$, however, this
dependence is made implicit at this moment to avoid cluttered notations. Next,
we define the worst-case Conditional Value-at-Risk function

$\mathcal{R}_{\alpha,\beta}(x)\triangleq\sup\limits_{\mathbb{Q}\in\mathcal{Q}}\mathbb{Q}\text{-CVaR}_{\varepsilon}\left(\max_{k\in\llbracket
1,K\rrbracket}\\{\xi^{\top}a_{k}(x)+b_{k}\\}\right).$

Notice that because $a_{k}(x)$ and $b_{k}$ are dependent on the scaling
matrices $\alpha$ and $\beta$, the worst-case CVaR also depends on these
scaling matrices and this dependence is made explicit. We now define the
following set

$\displaystyle\mathcal{X}^{\circ}(\alpha,\beta)=\left\\{x\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{J}|}:\mathcal{R}_{\alpha,\beta}(x)\leq
0\right\\}.$

We can show that $\mathcal{X}^{\circ}(\alpha,\beta)\subseteq\mathcal{X}$ for
any strictly positive scaling matrices
$\alpha,\beta\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$ [21]. As a
consequence, the optimization problem

(8) $\begin{array}[]{cl}\max&f(x)\\\
\mathrm{s.t.}&x\in\mathbb{X}\cap\mathcal{X}^{\circ}(\alpha,\beta)\end{array}$

constitutes a conservative approximation of the fairness-constrained
distributionally robust allocation problem (4).

######  Theorem 3.1 (Reformulation of worst-case CVaR).

For any fixed allocation $x\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{J}|}$,
scaling matrices $\alpha\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$
and $\beta\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$, the worst-
case CVaR $\mathcal{R}_{\alpha,\beta}(x)$ equals to the optimal value of the
linear optimization problem

$\begin{array}[]{cl}\min&\displaystyle\gamma+\mu^{\top}\lambda\\\\[8.61108pt]
\mathrm{s.t.}&\gamma\in\mathbb{R},\,\tau\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\,\eta_{k}\in\mathbb{R}_{+}^{M}\quad\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}b_{k}-(1-\varepsilon)\tau\leq\varepsilon(\gamma-h^{\top}\eta_{k})\\\
\varepsilon(G^{\top}\eta_{k}+\lambda)=a_{k}(x)\end{array}\right\\}\quad\forall
k\in\llbracket 1,K\rrbracket\\\ &\tau\leq\gamma-h^{\top}\eta_{0},\quad
G^{\top}\eta_{0}+\lambda=0.\end{array}$

###### Proof.

By denoting $(a)^{+}=\max\\{0,a\\}$, the worst-case CVaR can be expressed as

$\displaystyle\sup\limits_{\mathbb{Q}\in\mathcal{Q}}~{}\mathbb{Q}\text{-CVaR}_{\varepsilon}\left(\max_{k\in\llbracket
1,K\rrbracket}\\{\xi^{\top}a_{k}(x)+b_{k}\\}\right)$ (9a) $\displaystyle=$
$\displaystyle\sup\limits_{\mathbb{Q}\in\mathcal{Q}}\inf\limits_{\tau\in\mathbb{R}}\tau+\frac{1}{\varepsilon}\mathds{E}_{\mathbb{Q}}\left[\left(\max_{k\in\llbracket
1,K\rrbracket}\\{\xi^{\top}a_{k}(x)+b_{k}\\}-\tau\right)^{+}\right]$ (9b)
$\displaystyle=$
$\displaystyle\sup\limits_{\mathbb{Q}\in\mathcal{Q}}\inf\limits_{\tau\in\mathbb{R}}\mathds{E}_{\mathbb{Q}}\left[\max_{k\in\llbracket
0,K\rrbracket}\\{\xi^{\top}c_{k}(x)+d_{k}(\tau)\\}\right]$ (9c)
$\displaystyle=$
$\displaystyle\inf\limits_{\tau\in\mathbb{R}}\sup\limits_{\mathbb{Q}\in\mathcal{Q}}\mathds{E}_{\mathbb{Q}}\left[\max_{k\in\llbracket
0,K\rrbracket}\\{\xi^{\top}c_{k}(x)+d_{k}(\tau)\\}\right],$

where equality (9a) exploits the definition of CVaR. In (9b), we have defined

$\displaystyle c_{0}(x)=0,$ $\displaystyle d_{0}(\tau)=\tau,$ $\displaystyle
c_{k}(x)=\frac{a_{k}(x)}{\varepsilon},$ $\displaystyle
d_{k}(\tau)=\frac{b_{k}}{\varepsilon}+\left(1-\frac{1}{\varepsilon}\right)\tau\;\forall
k\in\llbracket 1,K\rrbracket.$

Finally, the equality in (9c) holds by Sion’s minimax theorem [31] because the
objective function is concave in $\mathbb{Q}$, convex in $\tau$, and the
ambiguity set $\mathcal{Q}$ is convex and weakly compact thanks to the
compactness of the support set $\Xi$. For any fixed value $x$ and for any
$\tau$, we can invoke Proposition 6.1 to rewrite the inner supremum problem in
(9c), we find

$\displaystyle\sup\limits_{\mathbb{Q}\in\mathcal{Q}}~{}\mathbb{Q}\text{-CVaR}_{\varepsilon}\left(\max_{k\in\llbracket
1,K\rrbracket}\\{\xi^{\top}a_{k}(x)+b_{k}\\}\right)$
$\displaystyle=\left\\{\begin{array}[]{cll}\inf&\displaystyle\gamma+\mu^{\top}\lambda\\\\[8.61108pt]
\mathrm{s.t.}&\gamma\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\,\eta_{k}\in\mathbb{R}_{+}^{M}&\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}h^{\top}\eta_{k}\leq\gamma-d_{k}(\tau)\\\
G^{\top}\eta_{k}=c_{k}(x)-\lambda\end{array}\right\\}&\forall k\in\llbracket
0,K\rrbracket.\end{array}\right.$

Substituting the expressions for $c_{k}(x)$ and $d_{k}(\tau)$ into the above
optimization problem and rearranging terms lead to the postulated linear
program. ∎

From what we have discussed so far, the optimization problem

(10) $\begin{array}[]{cl}\max&f(x)\\\ \mathrm{s.t.}&x\in\mathbb{X}\\\
&\gamma\in\mathbb{R},\,\tau\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\,\eta_{k}\in\mathbb{R}_{+}^{M}\quad\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}b_{k}-(1-\varepsilon)\tau\leq\varepsilon(\gamma-h^{\top}\eta_{k})\\\
\varepsilon(G^{\top}\eta_{k}+\lambda)=a_{k}(x)\end{array}\right\\}\quad\forall
k\in\llbracket 1,K\rrbracket\\\ &\tau\leq\gamma-h^{\top}\eta_{0},\quad
G^{\top}\eta_{0}+\lambda=0\\\ &\gamma+\mu^{\top}\lambda\leq 0\end{array}$

is a linear program and it constitutes a tractable conservative approximation
of (4) for any fixed value
$\alpha,\beta\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$. More
specifically, the optimal value of (10) is a lower bound on the optimal value
of (4) and the optimizer $x^{\star}$ in the variable $x$ of problem (10) is a
feasible solution to (4) for any fixed value of strictly positive scaling
parameters $\alpha,\beta$. One can systematically find the best lower bound by
jointly optimizing over $\alpha$ and $\beta$ in (10). Unfortunately,
optimizing jointly over the decision variables of (10) and $\alpha,\beta$ is a
non-convex optimization problem because it involves the bilinear terms between
$x$ and $\alpha,\beta$. Inspired by [38], we now devise an optimization
algorithm that sequentially optimizes over $(\alpha,\beta)$ and the allocation
$x$. To this end, we now consider for any fixed allocation $x\in\mathbb{X}$
the following optimization problem

(11)
$\min\left\\{\mathcal{R}_{\alpha,\beta}(x)~{}:~{}(\alpha,\beta)\in\mathcal{S}\right\\}$

that searches for $(\alpha,\beta)$ that minimizes the worst-case CVaR of a
given allocation $x$. To make problem (11) well-defined, we constrain
$(\alpha,\beta)$ to a compact set $\mathcal{S}$ defined as

(12)
$\mathcal{S}\triangleq\left\\{\begin{array}[]{l}(\alpha,\beta)\in\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}\times\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}:\\\
\displaystyle\sum_{(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}}~{}\alpha_{jj^{\prime}}=1,\sum_{(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}}~{}\beta_{jj^{\prime}}=1\\\
\alpha_{jj^{\prime}}\geq\delta,~{}\beta_{jj^{\prime}}\geq\delta\quad\forall(j,j^{\prime})\in\mathcal{J}\times\mathcal{J}\end{array}\right\\}$

for some strictly positive but small enough scalar $\delta$ so that
$\mathcal{S}$ is non-empty. We emphasize that the overall scale of $\alpha$
and $\beta$ is immaterial, and thus normalizing the sum of the elements in
$\alpha$ and $\beta$ to 1 in the definition of the set $\mathcal{S}$ does not
restrict any generality. The generic sequential algorithm is depicted in
Algorithm 1.

Algorithm 1 Generic sequential algorithm

0: Maximum iteration $T$, stopping tolerance $\theta$

Initialize $g_{0}\leftarrow+\infty$, $t\leftarrow 1$

Initialize $\alpha_{jj^{\prime}}\leftarrow 1/|\mathcal{J}|^{2}$,
$\beta_{jj^{\prime}}\leftarrow 1/|\mathcal{J}|^{2}\,\forall(j,j^{\prime})$

while $t\leq T$ do

Fix $(\alpha,\beta)$ and find a solution $(x_{t},v_{t})$ of  (13)
$g_{t}=\left\\{\begin{array}[]{cl}\max&f(x)-:\mathds{M}v\\\
\mathrm{s.t.}&x\in\mathbb{X},~{}v\in\mathbb{R}_{+},~{}\mathcal{R}_{\alpha,\beta}(x)\leq
v.\end{array}\right.$

if $|g_{t}-g_{t-1}|/g_{t}<\theta$ then break endif

Find $(\alpha,\beta)$ that solves (11) with $x$ being fixed to $x_{t}$

Set $t\leftarrow t+1$

end while

$x_{t}$

Problem (13) can be written explicitly as a linear program in the form of
problem (10) two minor modifications: (i) the right hand side of the last
constraint of (10) is replaced by the non-negative decision variable $v$, and
(ii) the auxiliary variables is penalized with a big-$\mathds{M}$ cost
parameter in the objective function. This additional decision variable $v$ is
necessary: indeed, it relaxes problem (10) and renders problem (10) feasible
even when the scaling parameters $(\alpha,\beta)$ are badly initialized.
Moreover, problem (11) can be expressed in the equivalent form as

(14) $\begin{array}[]{cl}\min&\gamma+\mu^{\top}\lambda\\\
\mathrm{s.t.}&(\alpha,\beta)\in\mathcal{S}\\\
&\gamma\in\mathbb{R},\,\tau\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\,\eta_{k}\in\mathbb{R}_{+}^{M}\quad\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}b_{k}(\alpha,\beta)-(1-\varepsilon)\tau\leq\varepsilon(\gamma-h^{\top}\eta_{k})\\\
\varepsilon(G^{\top}\eta_{k}+\lambda)=a_{k}(\alpha,\beta,x)\end{array}\right\\}\quad\forall
k\in\llbracket 1,K\rrbracket\\\ &\tau\leq\gamma-h^{\top}\eta_{0},\quad
G^{\top}\eta_{0}+\lambda=0,\end{array}$

where $a_{k}(\alpha,\beta,x)$ and $b_{k}(\alpha,\beta)$ are defined as in (6)
and (7), respectively, with their dependence on $(\alpha,\beta)$ being made
explicit. For any allocation $x$, $a_{k}(\alpha,\beta,x)$ and
$b_{k}(\alpha,\beta)$ are affine functions of $\alpha$ and $\beta$, and hence
problem (14) is a linear program.

The sequence $\\{g_{t}\\}_{t\in\mathbb{N}}$ is non-decreasing and thus it is
guaranteed to converge. If Algorithm 1 terminates with $v_{t}=0$, then the
allocation $x_{t}$ is feasible for problem (4) and the terminal value of
$g_{t}$ is a lower bound for the objective value of (4).

##  4\. Tractable Approximations for Hyperrectangular Support Set

Based on the general results of the previous section, we consider now the
hyperrectangular support set, where $\Xi$ admits the following specific
representation

(15) $\Xi=\left\\{\xi\in\mathbb{R}^{|\mathcal{I}|}:l\leq\xi\leq u\right\\}$

for some lower bound vector $l\in\mathbb{R}^{|\mathcal{I}|}$ and upper bound
vector $u\in\mathbb{R}^{|\mathcal{I}|}$ satisfying $l<u$. This
hyperrectangular set (15) is useful and practical in many ways. First,
prescribing the set $\Xi$ in this case requires the estimation of the upper
bound $u_{i}$ and lower bound $l_{i}$ for each task $i\in\mathcal{I}$, which
can be done efficiently given any available dataset. Furthermore, this set
$\Xi$ is intuitive, and it is easy to be communicated and explained to the
stakeholders involved in the decision making process, including the healthcare
workers, the managers of the healthcare unit and the human resources
department. Finally, even though the set $\Xi$ may look over-simplified, this
choice of the support set has been used widely in the literature and has shown
promising performance in inventory management with uncertain demand [28], in
vehicle routing problem with uncertain demand or travelling time [22]. Our aim
in this section is to detail the algorithm and the involved optimization
problems tailored for the hyperrectangular support set to facilitate the
implementation of this model for practical purposes.

When $\Xi$ is a hyperrectangle, the value of the worst-case CVaR
$\mathcal{R}_{\alpha,\beta}(x)$ can be reformulated explicitly using the
following corollary.

######  Corollary 4.1 (Reformulation for hyperrectangular support set).

Suppose that the support set $\Xi$ is a hyperrectangle as in (15). For any
fixed allocation $x\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{J}|}$, scaling
parameters $\alpha\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$ and
$\beta\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$, the worst-case
CVaR value $\mathcal{R}_{\alpha,\beta}(x)$ equals to the optimal value of the
linear optimization problem

$\begin{array}[]{cl}\min&\displaystyle\gamma+\mu^{\top}\lambda\\\\[8.61108pt]
\mathrm{s.t.}&\gamma\in\mathbb{R},\,\tau\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\\\
&\eta_{1k}\in\mathbb{R}_{+}^{|\mathcal{I}|},\,\eta_{2k}\in\mathbb{R}_{+}^{|\mathcal{I}|}\quad\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}b_{k}-(1-\varepsilon)\tau\leq\varepsilon(\gamma-u^{\top}\eta_{1k}+l^{\top}\eta_{2k})\\\
\varepsilon(\eta_{1k}-\eta_{2k}+\lambda)=a_{k}(x)\end{array}\right\\}\forall
k\in\llbracket 1,K\rrbracket\\\
&\tau\leq\gamma-u^{\top}\eta_{10}+l^{\top}\eta_{20},\quad\eta_{10}-\eta_{20}+\lambda=0.\end{array}$

This corollary follows directly from Theorem 3.1 by noticing that the support
set $\Xi$ in the form (15) can be expressed in the general form (3) with
$M=2|\mathcal{I}|$ constraints using

$G=\begin{bmatrix}I\\\
-I\end{bmatrix}\in\mathbb{R}^{2|\mathcal{I}|\times|\mathcal{I}|},\quad
h=\begin{pmatrix}u\\\ -l\end{pmatrix}\in\mathbb{R}^{2|\mathcal{I}|}.$

Furthermore, when $\Xi$ is dictated by (15), then for any fixed scaling
parameters
$\alpha,\beta\in\mathbb{R}_{++}^{|\mathcal{J}|\times|\mathcal{J}|}$, the
optimal value of the program

(16) $\begin{array}[]{cl}\max&f(x)\\\
\mathrm{s.t.}&x\in\mathbb{X},\,\gamma\in\mathbb{R},\,\tau\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\\\
&\eta_{1k}\in\mathbb{R}_{+}^{|\mathcal{I}|},\,\eta_{2k}\in\mathbb{R}_{+}^{|\mathcal{I}|}\quad\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}b_{k}-(1-\varepsilon)\tau\leq\varepsilon(\gamma-u^{\top}\eta_{1k}+l^{\top}\eta_{2k})\\\
\varepsilon(\eta_{1k}-\eta_{2k}+\lambda)=a_{k}(x)\end{array}\right\\}\forall
k\in\llbracket 1,K\rrbracket\\\
&\tau\leq\gamma-u^{\top}\eta_{10}+l^{\top}\eta_{20},\quad\eta_{10}-\eta_{20}+\lambda=0\\\
&\gamma+\mu^{\top}\lambda\leq 0\end{array}$

provides a valid lower bound on the optimal value of problem (4). We are now
ready to adapt the generic procedure delineated in Algorithm 1 to sequentially
solve over the decision $x$ and the scaling parameters $\alpha$ and $\beta$ in
order to find a competitive solution for problem (4). Notice that
$K=2|\mathcal{J}|^{2}$, and the optimization problem (13) over the decision
variables $x$ with $(\alpha,\beta)$ being fixed can be written using the
summation notations as

(17) $\begin{array}[]{cl}\max&f(x)-\mathds{M}v\\\
\mathrm{s.t.}&x\in\mathbb{X},\,v\in\mathbb{R}_{+},\,\gamma\in\mathbb{R},\,\tau\in\mathbb{R}\\\
&\eta_{1jj^{\prime}i}^{-}\in\mathbb{R}_{+},\,\eta_{2jj^{\prime}i}^{-}\in\mathbb{R}_{+}\quad\forall(j,j^{\prime},i)\in\mathcal{J}\times\mathcal{J}\times\mathcal{I}\\\
&\eta_{1jj^{\prime}i}^{+}\in\mathbb{R}_{+},\,\eta_{2jj^{\prime}i}^{+}\in\mathbb{R}_{+}\quad\forall(j,j^{\prime},i)\in\mathcal{J}\times\mathcal{J}\times\mathcal{I}\\\
&\eta_{10i}\in\mathbb{R}_{+},\,\eta_{20i}\in\mathbb{R}_{+},\,\lambda_{i}\in\mathbb{R}\quad\forall
i\in\mathcal{I}\\\
&\left.\begin{array}[]{l}-(1-\varepsilon)\tau-\alpha_{jj^{\prime}}\Delta\leq\varepsilon\big{(}\gamma-\displaystyle\sum_{i}(u_{i}\eta_{1jj^{\prime}i}^{-}-l_{i}\eta_{2jj^{\prime}i}^{-})\big{)}\\\
-(1-\varepsilon)\tau-\beta_{jj^{\prime}}\Delta\leq\varepsilon\big{(}\gamma-\displaystyle\sum_{i}(u_{i}\eta_{1jj^{\prime}i}^{+}-l_{i}\eta_{2jj^{\prime}i}^{+})\big{)}\end{array}\right\\}\forall(j,j^{\prime})\\\
&\left.\begin{array}[]{l}\varepsilon(\eta_{1jj^{\prime}i}^{-}-\eta_{2jj^{\prime}i}^{-}+\lambda_{i})=\alpha_{jj^{\prime}}(x_{ij}-x_{ij^{\prime}})\\\
\varepsilon(\eta_{1jj^{\prime}i}^{+}-\eta_{2jj^{\prime}i}^{+}+\lambda_{i})=\beta_{jj^{\prime}}(x_{ij^{\prime}}-x_{ij})\end{array}\right\\}\forall(j,j^{\prime},i)\\\
&\tau\leq\gamma-\displaystyle\sum_{i}(u_{i}\eta_{10i}-l_{i}\eta_{20i})\\\
&\eta_{10i}-\eta_{20i}+\lambda_{i}=0\quad\forall i\in\mathcal{I}\\\
&\gamma+\displaystyle\sum_{i}\mu_{i}\lambda_{i}\leq v,\end{array}$

where the notation $\forall(j,j^{\prime},i)$ in the constraints means
$\forall(j,j^{\prime},i)\in\mathcal{J}\times\mathcal{J}\times\mathcal{I}$.
Problem (17) is a mixed binary optimization problem because of the binary
requirements of the assignment feasible set $\mathbb{X}$.

Similarly, the optimization problem (11) that solves over $\alpha$ and $\beta$
with the values of $x$ being fixed can be written using the summation
notations as

(18)
$\begin{array}[]{cl}\min&\gamma+\displaystyle\sum_{i}\mu_{i}\lambda_{i}\\\
\mathrm{s.t.}&(\alpha,\beta)\in\mathcal{S},\,\gamma\in\mathbb{R},\,\tau\in\mathbb{R}\\\
&\eta_{1jj^{\prime}i}^{-}\in\mathbb{R}_{+},\,\eta_{2jj^{\prime}i}^{-}\in\mathbb{R}_{+}\quad\forall(j,j^{\prime},i)\in\mathcal{J}\times\mathcal{J}\times\mathcal{I}\\\
&\eta_{1jj^{\prime}i}^{+}\in\mathbb{R}_{+},\,\eta_{2jj^{\prime}i}^{+}\in\mathbb{R}_{+}\quad\forall(j,j^{\prime},i)\in\mathcal{J}\times\mathcal{J}\times\mathcal{I}\\\
&\eta_{10i}\in\mathbb{R}_{+},\,\eta_{20i}\in\mathbb{R}_{+},\,\lambda_{i}\in\mathbb{R}\quad\forall
i\in\mathcal{I}\\\
&\left.\begin{array}[]{l}-(1-\varepsilon)\tau-\alpha_{jj^{\prime}}\Delta\leq\varepsilon\big{(}\gamma-\displaystyle\sum_{i}(u_{i}\eta_{1jj^{\prime}i}^{-}-l_{i}\eta_{2jj^{\prime}i}^{-})\big{)}\\\
-(1-\varepsilon)\tau-\beta_{jj^{\prime}}\Delta\leq\varepsilon\big{(}\gamma-\displaystyle\sum_{i}(u_{i}\eta_{1jj^{\prime}i}^{+}-l_{i}\eta_{2jj^{\prime}i}^{+})\big{)}\end{array}\right\\}\forall(j,j^{\prime})\\\
&\left.\begin{array}[]{l}\varepsilon(\eta_{1jj^{\prime}i}^{-}-\eta_{2jj^{\prime}i}^{-}+\lambda_{i})=\alpha_{jj^{\prime}}(x_{ij}-x_{ij^{\prime}})\\\
\varepsilon(\eta_{1jj^{\prime}i}^{+}-\eta_{2jj^{\prime}i}^{+}+\lambda_{i})=\beta_{jj^{\prime}}(x_{ij^{\prime}}-x_{ij})\end{array}\right\\}\forall(j,j^{\prime},i)\\\
&\tau\leq\gamma-\displaystyle\sum_{i}(u_{i}\eta_{10i}-l_{i}\eta_{20i})\\\
&\eta_{10i}-\eta_{20i}+\lambda_{i}=0\quad\forall i\in\mathcal{I},\end{array}$

where the feasible set $\mathcal{S}$ is defined as in (12). Furthermore,
problem (18) is a linear continuous optimization problem, and techniques to
warm-start problem (18) can be applied to reduce the computational solution
time in the iterations of Algorithm 2.

Algorithm 2 Specific sequential algorithm for hypercube support set

0: Maximum iteration $T$, stopping tolerance $\theta$

Initialize $g_{0}\leftarrow+\infty$, $t\leftarrow 1$

Initialize $\alpha_{jj^{\prime}}\leftarrow 1/|\mathcal{J}|^{2}$,
$\beta_{jj^{\prime}}\leftarrow 1/|\mathcal{J}|^{2}\,\forall(j,j^{\prime})$

while $t\leq T$ do

Fix $(\alpha,\beta)$ and find a solution $(x_{t},v_{t})$ of problem (17)

Set $g_{t}$ to the optimal value of problem (17)

if $|g_{t}-g_{t-1}|/g_{t}<\theta$ then break endif

Find $(\alpha,\beta)$ that solves (18) with $x$ being fixed to $x_{t}$

Set $t\leftarrow t+1$

end while

$x_{t}$

##  5\. Numerical Experiment

We now showcase the prowess of the distributionally robust optimization
framework to promote fair assignments using numerical experiments. To reduce
the complexity of the experiments, we concentrate on the hyperrectangular
support set as presented in Section 4, and find the fair assignment using
Algorithm 2. All codes are implemented in MATLAB and are made available upon
request, and the linear optimization problems are solved using Gurobi 9. Due
to the nondisclosure agreement, we are unable to report the performance with
real-life data, thus we report in this paper only results using synthetic
data.

As a benchmark, we propose to compare the solution obtained from Algorithm 2
against the optimal solution of the optimization problem using the expected
values as input

(19) $\begin{array}[]{cl}\max&f(x)-\mathds{M}v\\\
\mathrm{s.t.}&x\in\mathbb{X},\,v\in\mathbb{R}_{+}\\\
&\displaystyle-\Delta-v\leq\sum_{i}\mu_{i}(x_{ij}-x_{ij^{\prime}})\leq\Delta+v\quad\forall(j,j^{\prime}),\end{array}$

where, once again, the variable $v$ is added to ensure the feasibility of the
optimization problem with a penalization parameter big-$\mathds{M}$ in the
objective function. Note that the random quantities $\xi$ are replaced by the
mean values $\mu$ in the constraints of problem (19).

The numerical experiment is conducted as follows. We set $|\mathcal{I}|=20$
tasks with $|\mathcal{J}|=5$ workers. We generate the vector
$\mu\in\mathbb{R}^{20}$ where each element is independent and identically
distributed (i.i.d.) with a uniform distribution in $[0,100]$, and the lower
bound $l=\mu-r$ and the upper bound $u=\mu+r$, where
$r=\min\\{\mu,r^{\prime}\\}$ and $r^{\prime}$ is a random vector whose
elements are i.i.d. with a uniform distribution on $[0,3]$. The minimum
operator in the definition of $r$ is required to make the values of the lower
bound $l$ non-negative. The reward $r_{ij}$ is generated i.i.d. from a uniform
distribution on $[0,100]$ for any $(i,j)\in\mathcal{I}\times\mathcal{J}$. We
set the threshold $\Delta=5$ and the tolerance $\varepsilon=5\%$, and then
invoke Algorithm 2 to obtain the distributionally robust assignment with
$T=40$ and $\theta=10^{-4}$, and solve problem (19) to obtain the assignment
using the expected values.

In the test process, we generate 10,000 samples $\\{\widehat{\xi}_{k}\\}$ from
a uniform distribution on $[l,u]$, where each element is drawn independently
of each other. From the $k$-th generated sample $\widehat{\xi}_{k}$, we record
the realized total time that worker $j$ needs to finish all the assigned tasks

$\widehat{T}_{j}=\sum_{i\in\mathcal{I}}\widehat{\xi}_{ki}x_{ij},$

where $x$ is either the distributionally robust assignment or the assignment
using mean values. The spread in the total time, which is a measurement of the
degree of fairness in the assignment, is computed as

(20) $s_{k}=\max\limits_{j}\widehat{T}_{j}-\min\limits_{j}\widehat{T}_{j}.$

If this spread exceeds the threshold, that is, $s_{k}>\Delta$, then the
realization $\widehat{\xi}_{k}$ generates a failure to the task assignment
plan. The probability that an assignment fails to satisfy the fairness
constraint is computed empirical as the total number of failure scenarios over
10,000 overall samples.

The aforementioned procedure is replicated 500 times to obtain 500 values of
the violation probability for the two approaches. The histogram comparing the
probability of violation is presented in Figure 1. One can observe that the
distributionally robust optimization approach successfully promotes the
fairness in the assignment by significantly reducing the probability that the
total time spread is bigger than the threshold $\Delta$. Averaging over 500
replications, the probability of violation is 1.37% for the distributionally
robust optimization and is 61.73% for the optimization using the expected
values approach. Notice that the average probability of 1.37% is significantly
lower than the prescribed tolerance of $\varepsilon=5\%$. This effect stems
from the combination of the distributionally robust approach and the
conservativeness of the joint chance constraint approximation.

Figure 2 compares the spreads of one particular replication using 10,000
samples. It is obvious that the spreads is smaller for the distributionally
robust approach compared to the optimization using the expected values
approach. Furthermore, from Figure 2a, the empirical spread is kept below the
threshold $\Delta=5$ with high probability, which demonstrates the
effectiveness of the distributionally robust optimization approach in
satisfying the joint chance constraint.

We emphasize that guaranteeing the fairness in the task assignment does not
come for free. Indeed, a fair assignment which is obtained by solving the
distributionally robust optimization problem often results in a lower reward
compared to the optimization with the expected values approach. Figure 3
depicts the rewards from each approach taken over 500 replications. The mean
reward for the distributionally robust optimization approach is 527.62, and
this quantity for the optimization using the expected values is 626.90.

(a) Distributionally robust optimization approach.

(b) Optimization with expected values approach.

Figure 1. Probability of violation taken over 500 replications for both
approaches.

(a) Distributionally robust optimization approach.

(b) Optimization with expected values approach.

Figure 2. Spreads, calculated in unit time as in (20) of a particular
replication with 10,000 samples.

(a) Distributionally robust optimization approach.

(b) Optimization with expected values approach.

Figure 3. Rewards taken over 500 replications for both approaches.

##  6\. Concluding Remarks

In this paper, we explore how the methodology of distributionally robust
optimization can be applied to generate fair task assignment. The notion of
fairness is measured by the difference in the time the workers are required to
finish the assigned tasks. We induce a fair assignment by imposing a joint
chance constraint that bounds the time difference by a certain threshold with
high probability, and we formulate the distributionally robust optimization
problem with joint chance constraint to alleviate the depedence on the
identification of the true underlying distribution, which is sometimes even
impossible. We develop a conservative approximation of the joint chance
constraint using Conditional Value-at-Risk, and propose a sequential
optimization algorithm to efficiently solve the resulting reformulations.
Finally, we tailor the result to a hyperrectangular support set and report
promising numerical results using synthetic datasets. The results in this
paper are a substantial part of the efforts from the hospital management team
of the Hanoi Obstetrics and Gynecology Hospital to improve the job
satisfaction of our healthcare workers, which in consequence aims to improve
the quality of our healthcare service.

The results in this paper reveal several directions for future research. While
this paper focuses on the mean-support ambiguity set, many other types of
ambiguity sets are available and can be used to promote fairness in the
allocation of tasks to healthcare workers. These new ambiguity sets include,
but is not restricted to, the ambiguity sets prescribed using the
$\phi$-divergence [4] and the Wasserstein distance [14]. Moreover, as we have
noted in Section 5, the assignment using the distributionally robust
optimization approach usually have lower rewards compared to other approaches
which is less stringent on the fairness criterion. However, this reduction in
the reward can be relieved in a certain extent by injecting more flexibility
in the recourse actions to adapt to the realization of the uncertain quantity
as the plan rolls over. This approach can be implemented using a re-
optimization on the second-stage decisions, which has been first demonstrated
in the power systems scheduling setting [23] and can be potentially applied to
this problem. We leave these ideas for future research.

## Proofs

######  Proposition 6.1 (Strong duality).

Suppose that the ambiguity set $\mathcal{Q}$ is defined as in (2) and that the
loss function $\ell(\xi)$ is a convex, piecewise affine function of the form
$\ell(\xi)=\max_{k\in\llbracket 0,K\rrbracket}\\{c_{k}^{\top}\xi+d_{k}\\}$ for
some $c_{k}\in\mathbb{R}^{|\mathcal{I}|}$ and $d_{k}\in\mathbb{R}~{}\forall
k\in\llbracket 0,K\rrbracket$. The worst-case expected loss can be
reformulated as

$\displaystyle\max\limits_{\mathbb{Q}\in\mathcal{Q}}~{}\mathds{E}_{\mathbb{Q}}[\ell(\xi)]=\left\\{\begin{array}[]{cll}\min&\displaystyle\gamma+\mu^{\top}\lambda\\\\[8.61108pt]
\mathrm{s.t.}&\gamma\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\,\eta_{k}\in\mathbb{R}_{+}^{M}&\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}h^{\top}\eta_{k}\leq\gamma-d_{k}\\\
G^{\top}\eta_{k}=c_{k}-\lambda\end{array}\right\\}&\forall k\in\llbracket
0,K\rrbracket.\end{array}\right.$

The result presented in Proposition 6.1 can be recovered as a special case of
[35, Theorem 1]. An elementary proof is presented here for completeness.

###### Proof.

By exploiting the definition of the ambiguity set $\mathcal{Q}$, the worst-
case expected loss can be written as an infinite dimensional optimization
problem

(21)
$\max\limits_{\mathbb{Q}\in\mathcal{Q}}~{}\mathds{E}_{\mathbb{Q}}[\ell(\xi)]=\left\\{\begin{array}[]{cl}\max&\displaystyle\int_{\mathbb{R}^{|\mathcal{I}|}}\ell(\xi)~{}\mathbb{Q}(\mathrm{d}\xi)\\\
\mathrm{s.t.}&\mathbb{Q}\in\mathcal{M}(\mathbb{R}^{|\mathcal{I}|})\\\
&\displaystyle\int_{\mathbb{R}^{|\mathcal{I}|}}\mathbbm{1}_{\Xi}(\xi)~{}\mathbb{Q}(\mathrm{d}\xi)=1\\\
&\displaystyle\int_{\mathbb{R}^{|\mathcal{I}|}}\xi~{}\mathbb{Q}(\mathrm{d}\xi)=\mu,\end{array}\right.$

where $\mathcal{M}(\mathbb{R}^{|\mathcal{I}|})$ denotes the set of all
positive measures on $\mathbb{R}^{|\mathcal{I}|}$ and $\mathbbm{1}_{\Xi}$ is
the indicator of the set $\Xi$, that is,

$\mathbbm{1}_{\Xi}(\xi)=\begin{cases}1&\text{if }\xi\in\Xi,\\\
0&\text{otherwise}.\end{cases}$

Because $\mu$ lies in the interior of the support set $\Xi$ by assumption, the
Slater type condition holds for problem (21). As a consequence, strong duality
holds by [29, Proposition 3.4], and the optimal value of problem (21) equals
to the optimal value of the following semi-infinite optimization problem

(22)
$\begin{array}[]{cl}\min&\displaystyle\gamma+\mu^{\top}\lambda\\\\[8.61108pt]
\mathrm{s.t.}&\lambda\in\mathbb{R}^{|\mathcal{I}|},\,\gamma\in\mathbb{R}\\\
&\displaystyle\gamma+\xi^{\top}\lambda\geq\ell(\xi)\quad\forall\xi\in\Xi.\end{array}$

By exploiting the piecewise affine form of the loss function $\ell$, the
constraint of problem (22) can be re-expressed as a system of $K$ semi-
infinite linear constraints

$\displaystyle\gamma+\xi^{\top}\lambda\geq\xi^{\top}c_{k}+d_{k}\quad\forall\xi\in\Xi\quad\forall
k\in\llbracket 0,K\rrbracket,$

which is further equivalent to a system of $K$ robust linear constraints

$\displaystyle\max_{\xi\in\Xi}~{}\xi^{\top}(c_{k}-\lambda)\leq\gamma-
d_{k}\quad\forall k\in\llbracket 0,K\rrbracket$

Formulating the dual linear program of the supremum problem on the left hand
side of the above constraint, we find that problem (22) is equivalent to

$\begin{array}[]{cll}\min&\displaystyle\gamma+\mu^{\top}\lambda\\\\[8.61108pt]
\mathrm{s.t.}&\gamma\in\mathbb{R},\,\lambda\in\mathbb{R}^{|\mathcal{I}|},\,\eta_{k}\in\mathbb{R}_{+}^{M}&\forall
k\in\llbracket 0,K\rrbracket\\\
&\left.\begin{array}[]{l}h^{\top}\eta_{k}\leq\gamma-d_{k}\\\
G^{\top}\eta_{k}=c_{k}-\lambda\end{array}\right\\}&\forall k\in\llbracket
0,K\rrbracket.\end{array}$

This observation completes the proof. ∎

## References

  * [1] C. Aluttis, T. Bishaw, and M. W. Frank, The workforce for health in a globalized context – global shortages and international migration, (2014).
  * [2] R. Bai, E. K. Burke, G. Kendall, J. Li, and B. McCollum, A hybrid evolutionary approach to the nurse rostering problem, IEEE Transactions on Evolutionary Computation, 14 (2010), pp. 580–590.
  * [3] A. Ben-Tal, L. El Ghaoui, and A. Nemirovski, Robust Optimization, Princeton Series in Applied Mathematics, Princeton University Press, October 2009\.
  * [4] A. Ben-Tal, D. D. Hertog, and J. Vial, Deriving robust counterparts of nonlinear uncertain inequalities, Mathematical Programming, 149 (2015), pp. 265–299.
  * [5] E. K. Burke, J. Li, and R. Qu, A hybrid model of integer programming and variable neighbourhood search for highly-constrained nurse rostering problems, European Journal of Operational Research, 203 (2010), pp. 484 – 493\.
  * [6] W. Chen, M. Sim, J. Sun, and C.-P. Teo, From CVaR to uncertainty set: Implications in joint chance-constrained optimization, Operations Research, 58 (2010), pp. 470–485.
  * [7] E. Delage and Y. Ye, Distributionally robust optimization under moment uncertainty with application to data-driven problems, Operations Research, 58 (2010), pp. 595–612.
  * [8] M. Dyer and A. Frieze, On the complexity of computing the volume of a polyhedron, SIAM Journal on Computing, 17 (1988), pp. 967–974.
  * [9] A. Gabrani, A. Hoxha, J. G. (CYCO), E. P. (ZAIMI), E. Zaimi, and E. Avdullari, Perceived organizational commitment and job satisfaction among nurses in albanian public hospitals: A cross-sectional study, International Journal of Healthcare Management, 9 (2016), pp. 110–118.
  * [10] C. A. Glass and R. A. Knight, The nurse rostering problem: A critical appraisal of the problem structure, European Journal of Operational Research, 202 (2010), pp. 379 – 389.
  * [11] M. H. A. Ibrahim, R. Ahmad, N. K. Ibrahim, S. Chuprat, and H. Haron, Nurse scheduling with fairness criteria for public hospital, in 2011 IEEE International Conference on Computer Applications and Industrial Electronics (ICCAIE), Dec 2011, pp. 91–95.
  * [12] A. Khani, M. Jaafarpour, and A. Dyrekvandmogadam, Quality of nursing work life, Journal of Clinical and Diagnostic Research, 2 (2008), pp. 1169–1174.
  * [13] Q. Kong, S. Li, N. Liu, C.-P. Teo, and Z. Yan, Appointment scheduling under time-dependent patient no-show behavior, Forthcoming in Management Science, (2020).
  * [14] D. Kuhn, P. M. Esfahani, V. A. Nguyen, and S. Shafieezadeh-Abadeh, Wasserstein distributionally robust optimization: Theory and applications in machine learning, INFORMS TutORials in Operations Research, (2019), pp. 130–166.
  * [15] P. Landa, M. Sonnessa, E. Tanfani, and A. Testi, Multiobjective bed management considering emergency and elective patient flows, International Transactions in Operational Research, 25 (2018), pp. 91–110.
  * [16] K. Liu, Q. Li, and Z.-H. Zhang, Distributionally robust optimization of an emergency medical service station location and sizing problem with joint chance constraints, Transportation Research Part B: Methodological, 119 (2019), pp. 79 – 101.
  * [17] B. Maenhout and M. Vanhoucke, Branching strategies in a branch-and-price approach for a multiple objective nurse scheduling problem, Journal of Scheduling, 13 (2010), pp. 77–93.
  * [18] T. Moradi, F. Maghaminejad, and I. Azizi-Fini, Quality of working life of nurses and its related factors, Nursing and Midwifery Studies, 3 (2014).
  * [19] M. Mousazadeh, S. A. Torabi, M. S. Pishvaee, and F. Abolhassani, Health service network design: a robust possibilistic approach, International Transactions in Operational Research, 25 (2018), pp. 337–373.
  * [20] J. W. Nelson and A. M. Cavanagh, Development of an international tool to measure nurse job satisfaction by testing the healthcare environment survey beyond jamaica and the united states to scotland: A cross sectional study utilizing exploratory factor analysis, International Journal of Healthcare Management, 11 (2018), pp. 314–318.
  * [21] A. Nemirovski and A. Shapiro, Convex approximations of chance constrained programs, SIAM Journal on Optimization, 17 (2007), pp. 969–996.
  * [22] V. Nguyen, J. Jiang, K. Ng, and K. Teo, Satisficing measure approach for vehicle routing problem with time windows under uncertainty, European Journal of Operational Research, 248 (2016), pp. 404 – 414.
  * [23] C. Ordoudis, V. Nguyen, D. Kuhn, and P. Pinson, Energy and reserve dispatch with distributionally robust joint chance constraints, Available on Optimization Online, (2018).
  * [24] P. Y. Permarupan, A. A. Mamun, N. Hayat, R. A. Saufi, and N. K. Samy, Nursing management challenges: Effect of quality of work life on depersonalization, International Journal of Healthcare Management, 0 (2020), pp. 1–10.
  * [25] A. Prékopa, Stochastic Programming, Springer, 1995.
  * [26] M. Ryu and R. Jiang, Nurse staffing under absenteeism: A distributionally robust optimization approach, arXiv preprint arXiv:1909.09875, (2019).
  * [27] J. Sangai and A. Bellabdaoui, Workload balancing in nurse scheduling problem models and discussion, in 2017 International Colloquium on Logistics and Supply Chain Management (LOGISTIQUA), 2017, pp. 82–87.
  * [28] C.-T. See and M. Sim, Robust approximation to multiperiod inventory management, Operations Research, 58 (2010), pp. 583–594.
  * [29] A. Shapiro, On duality theory of conic linear problems, in Semi-Infinite Programming, Kluwer Academic Publishers, 2001, pp. 135–165.
  * [30] A. Shapiro, D. Dentcheva, and A. Ruszczynski, Lectures on Stochastic Programming: Modeling and Theory, SIAM, 2014.
  * [31] M. Sion, On general minimax theorems, Pacific Journal of Mathematics, 8 (1958), pp. 171–176.
  * [32] K. Suleiman, Z. Hijazi, M. A. Kalaldeh, and L. A. Sharour, Quality of nursing work life and related factors among emergency nurses in Jordan, Journal of Occupational Health, 61 (2019), pp. 398–406.
  * [33] A. Uhde, N. Schlicker, D. P. Wallach, and M. Hassenzahl, Fairness and decision-making in collaborative shift scheduling systems, arXiv preprint arXiv:2001.09755v2, (2020).
  * [34] Y. Wang, Y. Zhang, and J. Tang, A distributionally robust optimization approach for surgery block allocation, European Journal of Operational Research, 273 (2019), pp. 740 – 753.
  * [35] W. Wiesemann, D. Kuhn, and M. Sim, Distributionally robust convex optimization, Operations Research, 62 (2014), pp. 1358–1376.
  * [36] S. Yanık, J. Kalcsics, S. Nickel, and B. Bozkaya, A multi-period multi-criteria districting problem applied to primary care scheme with gradual assignment, International Transactions in Operational Research, 26 (2019), pp. 1676–1697.
  * [37] Y. Zhang, S. Shen, and S. A. Erdogan, Distributionally robust appointment scheduling with moment-based ambiguity set, Operations Research Letters, 45 (2017), p. 139–144.
  * [38] S. Zymler, D. Kuhn, and B. Rustem, Distributionally robust joint chance constraints with second-order moment information, Mathematical Programming, 137 (2013), pp. 167–198.


# Classification of chemical compounds based on the correlation between in
vitro gene expression profiles††thanks: This study is supported by JSPS
KAKENSHI Grant Numbers JP20K12203.

Jun-ichi Takeshita Research Institute of Science for Safety and
Sustainability, National Institute of Advanced Industrial Science and
Technology (AIST), Tsukuba, Japan<EMAIL_ADDRESS>Akinobu Toyoda
School of Science and Engineering, University of Tsukuba, Tsukuba, Japan.
Hidenori Tani Environmental Management Research Institute, National Institute
of Advanced Industrial Science and Technology (AIST), Tsukuba, Japan.
Yasunori Endo Faculty of Systems and Information Engineering, University of
Tsukuba, Tsukuba, Japan. Sadaaki Miyamoto Faculty of Systems and Information
Engineering, University of Tsukuba, Tsukuba, Japan.

###### Abstract

Toxicity evaluation of chemical compounds has traditionally relied on animal
experiments; however, the demand for non-animal-based prediction methods for
toxicology of compounds is increasing worldwide. Our aim was to provide a
classification method for compounds based on in vitro gene expression
profiles. The in vitro gene expression data analyzed in the present study was
obtained from our previous study. The data concerned nine compounds typically
employed in chemical management. We used agglomerative hierarchical clustering
to classify the compounds; however, there was a statistical difficulty to be
overcome. We needed to properly extract RNAs for clustering from more than
30,000 RNAs. In order to overcome this difficulty, we introduced a
combinatorial optimization problem with respect to both gene expression levels
and the correlation between gene expression profiles. Then, the simulated
annealing algorithm was used to obtain a good solution for the problem. As a
result, the nine compounds were divided into two groups using 1,000 extracted
RNAs. Our proposed methodology enables read-across, one of the frameworks for
predicting toxicology, based on in vitro gene expression profiles.

keywords: Statistical classification, Multiobjective combinatorial
optimization, Chemical toxicity, Alternatives to animal experiments, In vitro
gene expression, Mathematical formulation

##  1 Introduction

Traditionally, toxicity evaluation of chemical compounds has relied on animal
experiments [[2]]. However, in terms of time, cost efficiency, and animal
welfare concerns, there is an increasing demand for the development of non-
animal-based methodologies for predicting chemical toxicity. Recently, several
elements have been proposed for use in toxicity prediction methods, such as
(quantitative) structure-activity relationships ((Q)SAR), quantitative
activity-activity relationships (QAAR), and read-across [[9]].

Read-across is a method whereby the toxicity of a given compound is predicted
without the use of animal test data, but instead based on the animal toxicity
data for similar compounds. One study, for example, conducted read-across
based on similarity in chemical structure and toxicology expert judgment
[[7]]. For read-across, it is important to properly group compounds using data
that does not rely on animal experiments. In order to facilitate the grouping
of compounds, the Organisation for Economic Co-operation and Development
(OECD) published the OECD QSAR Toolbox [[6, Section 4.3]], [[10]]. In
addition, the Japanese government and participating academic institutes
developed the Hazard Evaluation Support System Integrated Platform (HESS)
[[11]], [[6, Section 4.4]], [[8]]. These platforms group compounds based on in
silico parameters, that is, chemical structures and essential physicochemical
parameters; and by employing them, we can successfully group compounds based
on these parameters.

In order to predict toxicity, it is useful to use not only in silico
parameters but also in vitro parameters, because the latter can reflect
certain biological characteristics of compounds. In fact, in the case of
predicting hepatotoxicity, one of the most prevalent forms of toxicity, two
studies have reported that using in vitro parameters increased the accuracy of
discriminative models for predicting the presence or absence of
hepatotoxicity, compared to using in silico parameters [[4]], [[5]]. These
studies suggest the hypothesis that grouping compounds based on in vitro
parameters would increase the accuracy of read-across approaches for
predicting chemical toxicity, compared to grouping based on in silico
parameters.

The present study proposes a methodology for grouping compounds using in vitro
gene expression data. In order to group compounds, agglomerative hierarchical
clustering was applied; however, a statistical difficulty appeared. The gene
expression data, from which it was necessary to properly extract a limited
number of RNAs for clustering compounds, included more than 30,000 RNAs. In
order to overcome this difficulty, we introduced a multiobjective
combinatorial optimization problem with respect to both gene expression levels
and the correlation between gene expression profiles. Then, we applied the
simulated annealing algorithm, a metaheuristic algorithm, to obtain a good
solution for the multiobjective combinatorial optimization problem.

##  2 Methods and materials

###  2.1 Gene expression data

The present study used the data reported in [12]. In that study, cell-based
assays were conducted, in duplicate, using mouse embryonic stem cells, for
nine compounds typically employed in chemical management: bis-phthalate,
$p$-dicholorobenzene, phenol, trichloroethylene, benzene, chloroform,
$p$-cresol, and tri-$n$-butyl-phosphate. In other words, each compound had two
samples, and each group of nine samples had the same control condition. Then,
the gene expression levels were quantified using the fragments per kilobase of
exon per million mapped fragments (FPKMs), and reported for a total of 32,586
RNAs. We used these RNAs and their FPKMs in the present study’s analysis.

###  2.2 Gene expression ratio

For any compound-treated group, the following gene expression ratio was used
for each gene:

$f(x,y)=\log_{2}\left(\frac{y}{x}\right),$

where $x$ and $y$ denote the gene expression levels of the control and
compound-treated groups, respectively. Note that if $x$ or $y$ was zero, then
the next smallest value in the respective group was used.

###  2.3 Clustering

In order to group compounds, agglomerative hierarchical clustering (the
average linkage between the merged groups) was used, because this method can
be used for any dissimilarity measures. The following dissimilarity measure
was used for the hierarchical clustering: for any two compounds, $x$ and $y$,
the dissimilarity measure, or distance, between $x$ and $y$, say, $d(x,y)$ was
defined by

$d(x,y)=\frac{1-\text{corr}(x,y)}{2},$

where $\mbox{corr}(x,y)$ is the correlation coefficient between the respective
FPKM vectors of Compound $x$ and Compound $y$. The dissimilarity measure takes
a value between $0$ and $1$.

###  2.4 Selection of RNAs

In order to extract RNAs that clearly revealed the difference between
compounds, the present study introduced the following combinatorial
optimization problem to extract $n$ RNAs for a given natural number $n$:

objective function $\displaystyle\quad U=(1-\alpha)U_{1}+\alpha U_{2}\
(0\leq\alpha\leq 1),$ $\displaystyle\quad
U_{1}=\frac{1}{\text{Count}(w)}\times$ $\displaystyle\qquad\sum_{x,y\in
E,x\neq
y}w_{x,y}\frac{\left|(1/n)\sum_{i=1}^{n}(x_{i}-\overline{x})(y_{i}-\overline{y})\right|}{\left\\{(1/n)\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}\right\\}^{1/2}\left\\{(1/n)\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}\right\\}^{1/2}},$
$\displaystyle\quad
U_{2}=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{\|r_{i}\|}{\max_{j}\|r_{j}\|}\right),$
$\displaystyle\mbox{subject to\quad}\\{x_{i}\\},\\{y_{i}\\}\subset\Gamma\mbox{
and }\\#\\{x_{i}\\}=\\#\\{y_{i}\\}=n,$

where $E$ and $\Gamma$ denote the set of compound-treated groups and the set
of all the RNAs (32,586 RNAs), respectively; $x_{i}$ and $y_{i}$ denote the
expression ratios of RNAi for $x$ and $y\in E$, respectively; $\overline{x}$
and $\overline{y}$ denote the means of $x_{i}$ and $y_{i}$, respectively;
$w_{x,y}$ denotes weights and takes a value in $\\{-1,0,1\\}$; and
$\text{Count}(w)$ denotes the number of $w_{x,y}$ taking a number of $1$. Note
that, for some $x$ and $y$, if $w_{x,y}=1$, we can extract RNAs that increase
the correlation between Compound $x$ and Compound $y$; and if $w_{x,y}=-1$, we
can extract RNAs that decrease the correlation between Compound $x$ and
Compound $y$; otherwise, we are not interested in the correlation between the
two compounds. In addition, $r_{i}$ designates the FPKM vector for RNAi, and
$\|\cdot\|$ represents the Euclidean norm.

The function $U_{1}$ takes a value in $[0,1]$, since the sigma component is
divided by $\text{Count}(w)$. The function $U_{2}$ also takes a value in
$[0,1]$. Thus, the objective function $U$ takes a value in $[0,1]$, since
$\alpha$ is a parameter taking a value in $[0,1]$. Roughly speaking, the
function $U_{1}$ describes the strength of the correlation between two
compound-treated groups, $x$ and $y$, based on extracted RNAs. The function
$U_{2}$ describes the gene expression level of extracted RNAs. Then, the
function $U$ is a linear combination of $U_{1}$ and $U_{2}$.

In the present study, the simulated annealing algorithm, which was originally
introduced by [3] and [1], was used to obtain a good solution for the
combinatorial optimization problem. Let $n$ be the number of RNAs we want to
extract, $T$ be the initial temperature, $T_{t}\ (0<T_{t}<T)$ be the final
temperature, and $\gamma\ (0<\gamma<1)$ be the cooling rate. Then, the
following is the algorithm to extract RNAs.

Step 1

Choose $n$ RNAs randomly as the initial state. Let $R$ be the set of the
chosen $n$ RNAs, and then calculate the objective function $U$ using the set
$R$. In addition, set the initial temperature $T$.

Step 2

As a neighbor solution, generate a set $R^{\prime}$, in which random elements
of $R$ and $\overline{R}$ (the complement set of $R$) are exchanged. Then,
calculate the objective function $U^{\prime}$ using the set $R^{\prime}$.

Step 3

If $U<U^{\prime}$, then $R:=R^{\prime}$; otherwise $R:=R^{\prime}$ with the
following probability $P$:

$P=\exp\left(-\frac{|U^{\prime}-U|}{T}\right).$

Step 4

$T:=\gamma T.$ If $T\geq T_{t}$, output the set $R$ as the final state;
otherwise, go to Step 2.

##  3 Results and discussion

###  3.1 Grouping compounds using all the RNAs

Figure 1 shows a dendrogram obtained by applying agglomerative hierarchical
clustering (the average linkage between the merged groups) to the data set of
gene expression ratios, and Figure 2 a dendogram obtained by applying similar
clustering to the set of gene expression levels, for the nine compounds in
duplicate and all the RNAs (32,586 RNAs). In both figures, the $y$-axes show
the dissimilarity measures where two clusters were merged, while the $x$-axes
show the distribution of compounds. Note that subscripts 1 and 2 refer to the
sample numbers; that is, if any two compounds have the same subscript number,
the control conditions for the two compounds are identical.

In Figure 1, we can see that there are two robust clusters, and each cluster
consists of the nine compounds with the same subscript number. This implies
that the two compound groups strongly depend on the gene expression levels in
the control conditions. In Figure 2, the two control conditions are closely
proximate, but there is no significant difference between the 18 samples (the
nine compounds in duplicate). These results indicate that there is no
significant difference between the gene expression patterns of the two control
conditions. Therefore, we may infer that the significant difference between
the two control conditions in Figure 1 was the result of the incremental
accumulation of difference among the RNAs, as more than 30,000 RNAs were in
the data set. These results indicate that it is unreliable to use the
information of all the RNAs for classifying compounds, but instead we must
extract a limited number of RNAs in order to properly classify compounds using
in vitro gene expression data.

###  3.2 Grouping compounds using extracted RNAs

We used the parameters $n=3,000,1,000,100$, and $\alpha=0.0,0.1,0.2,0.3$, for
the combinatorial optimization problem to extract RNAs. Figures 3, 4, and 5
show dendrograms obtained by agglomerative hierarchical clustering (the
average linkage between the merged groups) applied to the data set with the
nine compounds in duplicate, and $3,000$, $1,000$, and $100$ extracted RNAs,
respectively. Each figure has four panels. The upper-left (a), upper-right
(b), lower-left (c), and lower-right (d) panes correspond to the cases of
$\alpha=0.0,0.1,0.2$, and $0.3$, respectively. In each panel in these figures,
the $y$-axis shows the dissimilarity measures where two clusters were merged,
while the $x$-axis shows the distribution of compounds. The meaning and
function of the subscripts are as in Figures 1 and 2.

In Figures 3, 4, and 5, it is clear that the two samples are initially merged
for all the compounds, and then different compounds are merged. However, when
$3,000$ RNAs were extracted (Figure 3), there are roughly $0.05$ dissimilarity
measures between the two samples of each compound, although the maximum ranges
are roughly $0.20$. This result suggests that $3,000$ RNAs are not suitable
for classifying compounds because the dissimilarity measures between the two
samples of each compound are not small enough, compared to the dissimilarity
measures between the compounds. On the other hand, when $1,000$ or $100$ RNAs
were extracted (Figures 4 and 5, respectively), and $\alpha=0.1$ or $0.2$, the
two samples are merged with sufficiently small dissimilarity measures, and we
can clearly see the difference between compounds. These results demonstrate
the necessity of limiting the number of RNAs when classifying compounds to
assess the effects of compounds on gene expression patterns.

Next, we compare the cases when $1,000$ and $100$ RNAs were extracted (Figures
4 and 5, respectively). When $\alpha=0.1$, neither case shows clear cluster
structures; however, when $\alpha=0.2$, with $1,000$ extracted RNAs, the nine
compounds are divided into two groups, one consisting of three compounds (bis-
phthalate, trichloroethylene, and tri-n-butyl-phosphate) and the other of the
remaining six compounds. These results indicate that not only limiting the
number of RNAs but also reflecting the sizes of RNAs is effective for
revealing the difference between the compounds. Figures 6 ($\alpha=0$) and 7
($\alpha=0.2$) show the scatter plots of all $32,586$ RNAs between Samples $1$
and $2$, in the case of bis-phthalate. The red plots indicate the $1,000$
extracted RNAs, and the blue plots the remaining RNAs. When $\alpha=0.0$, it
is clear that, for the most part, only RNAs with near-zero FPKMs are
extracted; whereas, when $\alpha=0.2$, there is an increased number of RNAs of
significant size. Thus, using $n=1,000$ and $\alpha=0.2$ in the combinatorial
optimization problem would be best for extracting RNAs in order to classify
these compounds.

##  4 Conclusion

The present study has been applied agglomerative hierarchical clustering
methods and a multiobjective combinatorial optimization problem to classify
chemical compounds, based on in vitro gene expression data. This approach
enables read-across based on in vitro parameters, and the prediction of
chemical toxicity. However, even if we can properly classify compounds using
in vitro parameters, the results may not show the same similarity between
compounds as in in vivo experiments. There is a gap in the respective
similarities based on in vitro parameters and in vivo experiments. Therefore,
we must develop a means of extracting RNAs that reflects the similarity
determined by in vivo experiments, and this requires collaboration with
toxicology experts. This is a consideration for the future, but our approach
would aid in such analysis.

## acknowledgement

The authors would like to express appreciation to Mr. Ryosuke Abe, Dr. Hiroshi
Aoki, Dr. Hiroaki Sato, Dr. Masaki Torimura, and Dr. Masashi Gamo for thier
fruitful discussions for the first version of manuscript.

## References

  * [1] Cerny, V., Thermodynamical approach to the traveling salesman problem : An efficient simulation algorithm, J. Optim. Theory Appl., 45 (1985), 41-51.
  * [2] Eaton, D.L. and Gilbert, S.G., Principles of toxicology, In: Klaassen, C.D. and Watkins, III J.B. (eds.) Casarett Doull’s Essentials Toxicology, Third Edition, 5-20, McGraw-Hill Education, NY, USA. 2015
  * [3] Kirkpatrick, S., Gelatt, C.D., and Vecchi, M.P., Optimization by simulated annealing, Science, 220 (1983), 671-680.
  * [4] Liu, J., Mansouri, K., Judson, R.S., Martin, M.T., Hong, H., Chen, M., Xu, X., Thomas, R.S., and Shah, I., Predicting hepatotoxicity using ToxCast in vitro bioactivity and chemical structure, Chem. Res. Toxicol., 28 (2015), 738-751.
  * [5] Low, Y., Uehara, T., Minowa, Y., Yamada, H., Ohno, Y., Urushidani, T., Sedykh, A., Muratov, E., Kuzmin, V., Fourches, D., Zhu, H., Rusyn, I., and Tropsha, A., Predicting drug-induced hepatotoxicity using QSAR and toxicogenomics approaches, Chem. Res. Toxicol., 24 (2011), 1251-1262.
  * [6] Madden, J.C., Tools for grouping chemicals and forming categories, In: Cronin, M. et al. (eds) Chemical Toxicity Prediction: Category Formation and Read-Across, 72-97, Royal Society of Chemical Publishing, London. (2013)
  * [7] Mellor, C.L., Schultz, T.W., Przybylak, K.R., Richarz, A.N., Bradbury, S.P., and Cronin, M.T.D. Read-across for rat oral gavage repeated-dose toxicity for short-chain mono-alkylphenols: A case study, Comput. Toxicol., 2 (2017), 1-11.
  * [8] NITE: Hazard Evaluation Support System Integrated Platform (HESS), http://www.nite.go.jp/en/chem/qsar/hess-e.html, Accessed 5 January 2021.
  * [9] OECD, Guidance on grouping of chemicals, Second edition, Number 194, ENV/JM/MONO(2014)4, Paris, France. 2014
  * [10] OECD, The OECD QSAR Toolbox. http://www.oecd.org/env/ehs/risk-assessment/theoecdqsartoolbox.htm, Accessed 5 January 2021
  * [11] Sakuratani, Y., Zhang, H.Q., Nishikawa, S., Yamazaki, K., Yamada, T., Yamada, J., Gerova, K., Chankov, G., Mekenyan, O., and Hayashi, M., Hazard Evaluation Support System (HESS) for predicting repeated dose toxicity using toxicological categories, SAR QSAR Environ. Res., 24 (2013), 351-63.
  * [12] Tani, H., Takeshita, J., Aoki, H., Nakmura, K., Abe, R., Toyoda, A., Endo, Y., Miyamoto, S., Gamo, M., Sato, H., and Torimura, M., Identification of RNA biomarkers for chemical safety screening in mouse embryonic stem cells using RNA-seq, Biochem. Biophys. Res. Commun., 512 (2019), 641-646.

## Figures

Figure 1: A dendrogram obtained by applying aggregative hierarchical
clustering (the average linkage between the merged groups) to the data of the
gene expression ratios for the nine compounds and 32,586 RNAs. The $y$-axis
marks the dissimilarity measures at which the clusters merge, and the $x$-axis
the distribution of the nine compounds in duplicate. Figure 2: A dendrogram
obtained by applying aggregative hierarchical clustering (the average linkage
between the merged groups) to the data of the gene expression levels for the
nine compounds and 32,586 RNAs. The $y$-axis marks the dissimilarity measures
at which the clusters merge, and the $x$-axis the distribution of the nine
compounds in duplicate. Figure 3: Four dendrograms obtained by applying
aggregative hierarchical clustering (the average linkage between the merged
groups) to the data of the gene expression ratios for the nine compounds and
$3,000$ extracted RNAs. The upper-left (a), upper-right (b), lower-left (c)
and lower-right (d) panels are the cases of $\alpha=0.0,0.1,0.2$, and $0.3$,
respectively. In each panel, the $y$-axis marks the dissimilarity measures at
which the clusters merge, and the $x$-axis the distribution of the nine
compounds in duplicate. Figure 4: Four dendrograms obtained by applying
aggregative hierarchical clustering (the average linkage between the merged
groups) to the data of the gene expression ratios for the nine compounds and
$1,000$ extracted RNAs. The upper-left (a), upper-right (b), lower-left (c)
and lower-right (d) panels are the cases of $\alpha=0.0,0.1,0.2$, and $0.3$,
respectively. In each panel, the $y$-axis marks the dissimilarity measures at
which the clusters merge, and the $x$-axis the distribution of the nine
compounds in duplicate. Figure 5: Four dendrograms obtained by applying
aggregative hierarchical clustering (the average linkage between the merged
groups) to the data of the gene expression ratios for the nine compounds and
$100$ extracted RNAs. The upper-left (a), upper-right (b), lower-left (c) and
lower-right (d) panels are the cases of $\alpha=0.0,0.1,0.2$, and $0.3$,
respectively. In each panel, the $y$-axis marks the dissimilarity measures at
which the clusters merge, and the $x$-axis the distribution of the nine
compounds in duplicate. Figure 6: Scatter plot of all the RNAs ($32,586$ RNAs)
between the sample 1 and 2 in case of bis-phthalate. The red plots indicate
the $1,000$ extracted RNAs in case of $\alpha=0.0$, and the blue plots are the
rest RNAs. The sizes of the extracted RNAs are almost zeros. Figure 7: Scatter
plot of all the RNAs ($32,586$ RNAs) between the sample 1 and 2 in case of
bis-phthalate. The red plots indicate the $1,000$ extracted RNAs in case of
$\alpha=0.2$, and the blue plots are the rest RNAs. The number of RNAs whose
sizes are not zeros increase, compared to the case of $\alpha=0.0$.