11institutetext: Australian Institute for Machine Learning, University of
Adelaide 22institutetext: Universidade Federal Rural de Pernambuco, Brazil
33institutetext: Universität Ulm, Germany 44institutetext: South Australian
Health and Medical Research Institute

# Self-supervised Mean Teacher for Semi-supervised Chest X-ray Classification

Fengbei Liu Yu Tian Filipe R. Cordeiro Vasileios Belagiannis
Ian Reid Gustavo Carneiro 1First two authors contributed equally to this
work.1First two authors contributed equally to this work.114††footnotemark:
4††footnotemark: 22331111

###### Abstract

The training of deep learning models generally requires a large amount of
annotated data for effective convergence and generalisation. However,
obtaining high-quality annotations is a laboursome and expensive process due
to the need of expert radiologists for the labelling task. The study of semi-
supervised learning in medical image analysis is then of crucial importance
given that it is much less expensive to obtain unlabelled images than to
acquire images labelled by expert radiologists. Essentially, semi-supervised
methods leverage large sets of unlabelled data to enable better training
convergence and generalisation than using only the small set of labelled
images. In this paper, we propose Self-supervised Mean Teacher for Semi-
supervised (S2MTS2) learning that combines self-supervised mean-teacher pre-
training with semi-supervised fine-tuning. The main innovation of S2MTS2 is
the self-supervised mean-teacher pre-training based on the joint contrastive
learning, which uses an infinite number of pairs of positive query and key
features to improve the mean-teacher representation. The model is then fine-
tuned using the exponential moving average teacher framework trained with
semi-supervised learning.We validate S2MTS2 on the thorax disease multi-label
classification problem from the dataset Chest X-ray14, [fengbei]: CheXpert and
ISIC2018, where we show that it outperforms the previous SOTA semi-supervised
learning methods by a large margin. Our code will be available at
https://github.com/fengbeiliu/semi-chest.

###### Keywords:

Semi-supervised learning Chest X-ray Self-supervised learning Multi-label
classification.

##  1 Introduction

Deep learning has shown outstanding results in medical image analysis problems
[lee2017deep, tian2020few, li2018thoracic, tian2019one, liu2020self,
jonmohamadi2020automatic]. However, this performance usually depends on the
availability of labelled datasets, which is expensive to obtain given that the
labelling process requires expert doctors. This limitation motivates the study
of semi-supervised learning (SSL) methods that train models with a small set
of labelled data and a large set of unlabelled data.

The current state-of-the-art (SOTA) SSL is based on pseudo-labelling methods
[rizve2021defense, lee2013pseudo], consistency-enforcing approaches
[berthelot2019remixmatch, tarvainen2017mean, laine2016temporal], self-
supervised and semi-supervised learning (S4L) [chen2020big, zhai2019s4l], and
graph-based label propagation [aviles2019graphx]. Pseudo-labelling is an
intuitive SSL technique, where confident predictions from the model are
transformed into pseudo-labels for the unlabelled data, which are then used to
re-train the model [lee2013pseudo]. Consistency-enforcing regularisation is
based on training for a consistent output given model [liu2020semi,
tarvainen2017mean] or input data [berthelot2019remixmatch, laine2016temporal]
perturbations. S4L methods are based on self-supervised pre-training [moco,
simclr], followed by supervised fine-tuning using few labelled samples
[chen2020big, zhai2019s4l]. Graph-based methods rely on label propagation on
graphs [aviles2019graphx]. Recently, Yang et al. [yang2020rethinking]
suggested that self-supervision pre-training provides better feature
representations than consistency-enforcing approaches in SSL. However,
previous S4L approaches use only the labelled data in the fine-tuning stage,
missing useful training information present in the unlabelled data.
Furthermore, self-supervised pre-training [simclr, moco] tends to use limited
amount of samples to represent each class, but recently, Cai et al.
[cai2020joint] showed that better representation can be obtained with an
infinite amount of samples. Also, recent research [rizve2021defense] suggests
that the student-teacher framework, such as the mean-teacher
[tarvainen2017mean], works better in multi-label semi-supervised tasks than
other SSL methods. We speculate that this is because other methods are usually
designed to work with softmax activation that only works in multi-class
problems, while mean-teacher [tarvainen2017mean] does not have this constraint
and can work in multi-label problems.

In this paper, we propose a self-supervised mean-teacher for semi-supervised
(S2MTS2) learning approach that combines S4L [chen2020big, zhai2019s4l] with
consistency-enforcing learning based on the mean-teacher algorithm
[tarvainen2017mean]. The main contribution of our method is the self-
supervised mean-teacher pre-training with the joint contrastive learning
[cai2020joint]. To the best of our knowledge, this is the first approach, in
our field, to train the mean teacher model with self-supervised learning. This
model is then fine-tuned with semi-supervised learning using the exponential
moving average teacher framework [tarvainen2017mean]. [fengbei]: We evaluate
proposed method on two thorax disease multi-label datasets (ChestX-ray 14
[wang2017chestx] and CheXpert [irvin2019chexpert]). We also test our method on
a multi-class skin condition dataset ISIC2018 [tschandl2018ham10000,
codella2019skin].we show that our method outperforms the SOTA on semi-
supervised learning [aviles2019graphx, liu2020semi, gyawali2020semi,
unnikrishnan2020semi]. Moreover, we investigate each component of our
framework for their contribution to the overall model in the ablation study.

Figure 1: Description of the proposed self-supervised mean-teacher for semi-
supervised (S2MTS2) learning. The main contribution of the paper resides in
the top part of the figure, with the self-supervised mean-teacher pre-training
based on joint contrastive learning, which uses an infinite number of pairs of
positive query and key features sampled from the unlabelled images to minimise
$\ell_{p}(.)$ in (1). This model is then fine-tuned with the exponential
moving average teacher in a semi-supervised learning framework that uses both
labelled and unlabelled sets to minimise $\ell_{cls}(.)$ and $\ell_{con}(.)$
in (2).

##  2 Related Works

SSL is a research topic that is gaining attention from the medical image
analysis community due to the expensive image annotation process
[cheplygina2019not] and the growing number of large-scale datasets available
in the field [wang2017chestx]. The current SOTA SSL methods are based on
consistency-enforcing approaches that leverage the unlabelled data to
regularise the model prediction consistency [tarvainen2017mean,
laine2016temporal]. Other related papers [cui2019semi] extend the mean teacher
[tarvainen2017mean] to encourage consistency between the prediction by the
student and teacher models for atrium and brain lesion segmentation. The SOTA
SSL method on Chest X-ray images [liu2020semi] exploits the consistency in the
relations between labelled and unlabelled data. None of these methods explores
a self-supervised consistency-enforcing method to pre-train an SSL model, as
we propose.

Self-supervised learning methods [simclr, moco] are also being widely
investigated in SSL because they can provide good representations
[chen2020big, zhai2019s4l]. However, these methods ignore the large amount of
unlabelled data to be used during SSL, which may lead to unsatisfactory
generalisation process. An important point in self-supervised learning is on
how to define the classes to be learned. In general, each class is composed of
a single pair of augmented images from the same image, and many pairs of
augmentations from different images [moco, simclr, cai2020joint, chen2020big].
The use of a single pair of images to form a class has been criticised by Cai
et al. [cai2020joint], who propose the joint contrastive learning (JCL), which
is an efficient way to form a class with an infinite number of augmented
images from the same image to leverage the statistical dependency between
different augmentations.

##  3 Method

In this section, we introduce our two-stage learning framework in detail (see
Fig. 1). We assume that we have a small labelled dataset, denoted by
$\mathcal{D}_{L}=\\{(\mathbf{x}_{i},\mathbf{y}_{i})\\}_{i=1}^{|\mathcal{D}_{L}|}$,
where the image is represented by
$\mathbf{x}\in\mathcal{X}\subset\mathbb{R}^{H\times W\times C}$, and class
$\mathbf{y}\in\\{0,1\\}^{|\mathcal{Y}|}$, where $\mathcal{Y}$ represents the
label set. We consider a multi-label problem and thus
$\sum_{c=1}^{|\mathcal{Y}|}\mathbf{y}_{i}(c)\in[0,|\mathcal{Y}|]$. The
unlabelled dataset is defined by
$\mathcal{D}_{U}=\\{\mathbf{x}_{i}\\}_{i=1}^{|\mathcal{D}_{U}|}$ with
$|\mathcal{D}_{L}|<<|\mathcal{D}_{U}|$.

Our model consists of a student and a teacher model [tarvainen2017mean],
denoted by parameters $\theta,\theta^{\prime}\in\mathbf{\Theta}$,
respectively, which parameterize the classifier
$f_{\theta}:\mathcal{X}\rightarrow[0,1]^{|\mathcal{Y}|}$. This classifier can
be decomposed as $f_{\theta}=h_{\theta_{1}}\circ g_{\theta_{2}}$, with
$g_{\theta_{2}}:\mathcal{X}\rightarrow\mathcal{Z}$ and
$h_{\theta_{1}}:\mathcal{Z}\rightarrow[0,1]^{|\mathcal{Y}|}$. The first stage
(top of Fig. 1) of the training consists of a self-supervised learning that
uses the images from $\mathcal{D}_{L}$ and $\mathcal{D}_{U}$, denoted by
$\mathcal{D}^{\mathcal{X}}=\\{\mathbf{x}_{i}|\mathbf{x}_{i}\in\mathcal{D}^{\mathcal{X}}_{L}\bigcup\mathcal{D}_{U}\\}_{i=1}^{|\mathcal{D}_{\mathcal{X}}|}$,
with $\mathcal{D}^{\mathcal{X}}_{L}$ representing the images from the set
$\mathcal{D}_{L}$, where our method minimises the joint contrastive learning
loss [cai2020joint], defined in (1). This means that during this first stage,
we only learn the parameters for $g_{\theta_{2}}$. The second stage (bottom of
Fig. 1) fine-tunes this pre-trained student-teacher model using the semi
supervised consistency loss defined in (2). Below we provide details on the
losses and training.

###  3.1 Joint Contrastive Learning to Self-supervise the Mean-teacher Pre-
training

The self-supervised pre-training of the mean-teacher using joint contrastive
learning (JCL) [cai2020joint], presented in this section, is the main
technical contribution of this paper. The teacher and student process an input
image to return the keys $\mathbf{k}\in\mathcal{Z}$ and the queries
$\mathbf{q}\in\mathcal{Z}$ with
$\mathbf{k}=g_{\theta^{\prime}_{2}}(\mathbf{x})$ and
$\mathbf{q}=g_{\theta_{2}}(\mathbf{x})$. We also assume that we have a set of
augmentation functions, i.e., random crop and resize, rotation and Gaussian
blur, denoted by
$\mathcal{A}=\\{a_{l}:\mathcal{X}\rightarrow\mathcal{X}\\}_{l=1}^{|\mathcal{A}|}$.
Then JCL minimises the following loss [cai2020joint]:

${\ell_{p}(\mathcal{D}^{\mathcal{X}},\theta_{2},\theta^{\prime}_{2})=-\frac{1}{|\mathcal{D}^{\mathcal{X}}|}\frac{1}{M}\sum_{i=1}^{|\mathcal{D}^{\mathcal{X}}|}\sum_{m=1}^{M}\left[\log\frac{\exp\left[\frac{1}{\tau}\mathbf{q}_{i}^{\top}\mathbf{k}_{i,m}^{+}\right]}{\exp\left[\frac{1}{\tau}\mathbf{q}_{i}^{\top}\mathbf{k}_{i,m}^{+}\right]+\sum_{j=1}^{K}\exp\left[\frac{1}{\tau}\mathbf{q}_{i}^{\top}\mathbf{k}_{i,j}^{-}\right]}\right],}$
(1)

where $\tau$ is the temperature hyper-parameter, the query
$\mathbf{q}_{i}=g_{\theta_{2}}(a(\mathbf{x}_{i}))$, with $a\in\mathcal{A}$.
the positive key $\mathbf{k}_{i,m}^{+}\sim p(\mathbf{k}_{i}^{+})$, with
$p(\mathbf{k}_{i}^{+})=\mathcal{N}(\mu_{\mathbf{k}_{i}},\Sigma_{\mathbf{k}_{i}})$
and $\mathbf{k}_{i}=g_{\theta^{\prime}_{2}}(a(\mathbf{x}_{i}))$ (i.e., a
sample from the data augmentation distribution for $\mathbf{x}$), and the
negative keys
$\mathbf{k}_{i,j}^{-}\in\\{\mu_{\mathbf{k}_{j}}\\}_{i,j\in\\{1,..,|\mathcal{D}^{\mathcal{X}}|\\},i\neq
j}$ represents a negative key for query $\mathbf{q}_{i}$. In (1), $M$ denotes
the number of positive keys, and Cai et al. [cai2020joint] describe a loss
that minimises a bound to (1) for $M\to\infty$ – below, the minimisation of
$\ell_{p}(.)$ in (1) is realised by the minimisation of this bound. As defined
above, the generative model $p(\mathbf{k}_{i}^{+})$ is denoted by the Gaussian
$\mathcal{N}(\mu_{\mathbf{k}_{i}},\Sigma_{\mathbf{k}_{i}})$, where the mean
$\mu_{\mathbf{k}_{i}}$ and covariance $\Sigma_{\mathbf{k}_{i}}$ are estimated
from a set of keys
$\\{\mathbf{k}^{+}_{i,l}=g_{\theta^{\prime}_{2}}(a_{l}(\mathbf{x}_{i}))\\}_{a_{l}\in\mathcal{A}}$
formed by different views of $\mathbf{x}_{i}$. The set of negative keys
$\\{\mu_{\mathbf{k}_{j}}\\}_{i,j\in\\{1,..,|\mathcal{D}^{\mathcal{X}}|\\},i\neq
j}$ is stored in a memory queue [moco] that is updated in a first-in-first-out
way, where the mean of the keys in
$\\{\mu_{\mathbf{k}_{i}}\\}_{i=1}^{|\mathcal{D}^{\mathcal{X}}|}$ are inserted
to the memory queue to replace the oldest key means from previous training
iterations. The memory queue has been designed to increase the number of
negative samples without sacrificing computation efficiency.

The training of the student-teacher model [moco, zhou2020comparing,
tarvainen2017mean] is achieved by updating the student parameter using the
loss in (1), as in
$\theta_{2}(t)=\theta_{2}(t-1)-\nabla_{\theta_{2}}\ell_{p}(\mathcal{D}^{\mathcal{X}},\theta_{2},\theta^{\prime}_{2})$,
where $t$ is the training iteration. The teacher model parameter is updated
with exponential moving average (EMA) with
$\theta^{\prime}_{2}(t)=\alpha\theta^{\prime}_{2}(t-1)+(1-\alpha)\theta_{2}(t)$,
with $\alpha\in[0,1]$. For this pre-training stage, we notice that training
for more epochs always improve the model regularisation given that it is
difficult to overfit the training set with the loss in (1). Hence, we select
the last epoch student model $g_{\theta_{2}}(.)$ to initialise the fine-tuning
stage, defined below in Sec. 3.2.

###  3.2 Fine-tuning the Mean Teacher

To fine tune the mean teacher, we follow the approach in [moco,
tarvainen2017mean] using the following loss to train the student model:

$\ell_{t}(\mathcal{D}_{L},\mathcal{D}_{U},\theta,\theta^{\prime})=\frac{1}{|\mathcal{D}_{L}|}\sum_{(\mathbf{x}_{i},\mathbf{y}_{i})\in\mathcal{D}_{L}}\ell_{cls}(\mathbf{y}_{i},f_{\theta}(\mathbf{x}_{i}))+\frac{1}{|\mathcal{D}|}\sum_{\mathbf{x}_{i}\in\mathcal{D}}\ell_{con}(f_{\theta}(\mathbf{x}_{i}),f_{\theta^{\prime}}(\mathbf{x}_{i})),$
(2)

where
$\ell_{cls}(\mathbf{y}_{i},f_{\theta}(\mathbf{x}_{i}))=-\mathbf{y}_{i}^{\top}\log(f_{\theta}(\mathbf{x}_{i}))$,
$\ell_{con}(f_{\theta}(\mathbf{x}_{i}),f_{\theta^{\prime}}(\mathbf{x}_{i}))=\|f_{\theta}(\mathbf{x}_{i})-f_{\theta^{\prime}}(\mathbf{x}_{i})\|^{2}$,
and $\mathcal{D}=\mathcal{D}_{U}\bigcup\mathcal{D}^{\mathcal{X}}_{L}$. The
training of the student-teacher model [moco, tarvainen2017mean,
zhou2020comparing] is achieved by updating the student parameter using the
loss in (2), as in
$\theta(t)=\theta(t-1)-\nabla_{\theta}\ell_{t}(\mathcal{D}_{L},\mathcal{D}_{U},\theta,\theta^{\prime})$,
where $t$ is the training iteration. The teacher model parameter is updated
with exponential moving average (EMA) with
$\theta^{\prime}(t)=\alpha\theta^{\prime}(t-1)+(1-\alpha)\theta(t)$, with
$\alpha\in[0,1]$. After finishing the fine-tuning stage, we select the teacher
model $f_{\theta^{\prime}}(.)$ to estimate the multi-label classification for
test images.

##  4 Experiment

###  4.1 Dataset Setup

[fengbei]: We use Chest X-ray14 [wang2017chestx], CheXpert [irvin2019chexpert]
and ISIC2018 [tschandl2018ham10000, codella2019skin] to evaluate our method.
Chest X-ray14 The dataset contains 112,120 chest x-ray images from 30,805
different patients. There are 14 different labels (each label represents a
disease) in the dataset, where each patient can have multiple diseases at the
same time, forming a multi-label classification problem. To compare with
previous papers [aviles2019graphx, liu2020semi], we adopt the official
train/test data split. For the self-supervised pre-training of the mean
teacher, we used all the unlabelled images (86k samples) from the training
set. For the semi-supervised fine-tuning of the mean teacher, we follow the
papers [aviles2019graphx, liu2020semi] and experiment with training sets
containing different proportions of labelled data (2%,5%,10%,15%,20%). We
report the classification result on the official test set (26k samples) using
area under the ROC curve (AUC).
CheXpert The dataset contains 220K+ images with 14 different diseases. For
pre-processing, we remove all lateral view images and treat uncertain label as
negative labels. We follow [gyawali2020semi] semi-supervised setup and
experiment with 100/200/300/400/500K number of labeled subset. We report
result on official test set using AUC.
ISIC2018 ISIC 2018 skin condition dataset consists 10,015 images with seven
different labels. Each image is associated with one of the seven labels,
forming a multi-class classification problem. We follow [liu2020semi]
train/test split for fair comparison. We report AUC, F1 score and Sensitivity
compare with baselines.

###  4.2 Implementation Details

[fengbei]: For all dataset, we use DenseNet121 [huang2017densely] as our
backbone model. For self-supervised pre-training, we follow [chen2020big] and
replace the projectionn head from a two-layer multi-layer perceptron (MLP) to
a three-layer MLP. For dataset pre-processing, we resized Chest X-ray14 images
from 1024 $\times$ 1024 to 512 $\times$ 512 for faster processing and CheXpert
and ISIC2018 to 224 $\times$ 224 for fair comparison. We use the data
augmentation proposed in [simclr], consisting of random resize and crop,
random rotation, random horizonntal flipping, except for random grayscale
because X-ray images are originally in grayscale. The batch size is 128 for
Chest X-ray14 and 256 for CheXpert and ISIC2018, learning rate is 0.05. For
finee-tuning stage, we use batch size 32 with 16 labeled and 16 unlabeled. We
trained w with learning rate deecayed by 0.1 at 15 and 25 epochs for all
dataset.

For all dataset, We use the DenseNet-121 [huang2017densely] as our backbone
model. We use Adam optimiser and learning rate of 0.03. Furthermore, we follow
Chen et al.’s discovery [chen2020big] and replace the projection head from a
two-layer multi-layer perceptron (MLP) [simclr] to a three-layer MLP.
similarly to [hermoza2020region, wang2017chestx]. We resized the original
images from 1024 $\times$ 1024 to 512 $\times$ 512 for faster processing. For
the self-supervised pre-training of the mean teacher, we use the data
augmentation proposed in [simclr], consisting of random resize and crop,
random rotation, and random horizontal flipping, except for the random
grayscale because the X-ray images are originally in grayscale. We use Adam
optimiser and learning rate of 0.03. Furthermore, we follow Chen et al.’s
discovery [chen2020big] and replace the projection head from a two-layer
multi-layer perceptron (MLP) [simclr] to a three-layer MLP. The pre-training
stage takes 100 epochs,

the batch size is set to 16, and learning rate is 0.05. The fine-tuning takes
30 epochs with learning rate decayed by 0.1 at 15, 25 epochs. The whole code
is written in Pytorch. We use 4 Nvidia Volta-100 for the self-supervised stage
and 1 Nvidia RTX 2080ti for fine-tuning.

###  4.3 Experimental Results

In this section, we show the semi-supervised experimental results on ChestX-
ray14. We evaluate our approach based on different percentage of labelled data
(i.e., 2$\%$, 5$\%$, 10$\%$, 15$\%$, 20$\%$) on the official test set, as
shown in Table 6. The set of labelled data used for each percentage above
follows the same strategy of previous works [aviles2019graphx, liu2020semi].
Our S4L achieves the SOTA AUC results on all different percentages of labels.
Our model surpasses the previous SOTA SRC-MT [liu2020semi] by a large margin
of 8.7% and 6.8% AUC gap for the 2% and 5% labelled set cases, respectively,
where we use a backbone architecture of lower complexity (Densenet121 instead
of the DenseNet169 of [liu2020semi]). Using the same Densenet121 backbone,
GraphXnet [aviles2019graphx] fails to classify precisely for the 2% and 5%
labelled set cases. Our method surpasses GraphXnet by more than 20% AUC gap in
both cases. Furthermore, we achieve the SOTA results of the field for the 10%,
15% and 20% labelled set cases, outperforming all previous semi-supervised
methods [liu2020semi, aviles2019graphx]. It is worth noting that our model
trained with 5% of the labelled set achieves better results than SRC-MT with
15% of labelled.We also compare with a recently proposed self-supervised
methods, MoCo V2 [chen2020improved], adapted to our semi-supervised task,
using MoCo V2 to pretrain the network followed by the fine-tuning stage using
different percentages of labelled data. Our method outperforms MoCo V2 by
almost 10% AUC when using 2% of labelled set, and almost 3% AUC for 10% of
labelled set. Our result for 20% labelled set achieves comparable 81.06% AUC
performance as the supervised learning approaches – 81.20% from MoCo V2
(Densenet 121) and 81.75% from SRC-MT (Densenet 169) using 100% of the
labelled samples. Such result indicates the effectiveness of our proposed
S2MTS2 in SSL benchmark problems.

We also show the class-level performance using 20% of the labelled data and
compare with other SOTA methods in Tab. 2. We compare with the previous
baselines, namely original mean teacher (MT) with Densenet169, SRC-MT with
Densenet169, MoCo V2, and GraphXNet with Densenet121. We also train a baseline
Densenet121 model with 20% labelled data using Imagenet pre-trained model. Our
method achieves the best results on nine classes, surpassing the original MT
[tarvainen2017mean] and its extension SRC-MT [liu2020semi] by a large margin,
demonstrating the effectiveness of our self-supervised learning.

Furthermore, we compare our approach on the fully-supervised Chest X-ray14
benchmark in Tab. 3. To the best of our knowledge, Hermoza et al.
[hermoza2020region] has the SOTA supervised classification method containing a
complex structure (relying on the weakly-supervised localisation of lesions)
with a mean AUC of 82.1% (over the 14 classes), while ours reports a mean AUC
of 82.5%. Hence, our model, using the whole labelled set, achieves the SOTA
performance on 8 classes and an average that surpasses the previous supervised
methods by a minimum of 0.4% and a maximum of 8% AUC. This indicates that our
method not only performs well on semi-supervised tasks but also achieves SOTA
results on supervised learning. In Fig. 2, we show some examples of
predictions by our model (trained with 20% labelled data) on the test set of
Chest X-ray14 [wang2017chestx].

Label Percentage | 2% | 5% | 10% | 15% | 20% | 100%
---|---|---|---|---|---|---
Graph XNet* [aviles2019graphx] | 53.00 | 58.00 | 63.00 | 68.00 | 78.00 | N/A
SRC-MT* [liu2020semi] | 66.95 | 72.29 | 75.28 | 77.76 | 79.23 | 81.75
NoTeacher [unnikrishnan2020semi] | 72.60 | 77.04 | 77.61 | N/A | 79.49 | N/A
MOCO V2 [chen2020improved] | 65.97 | 73.84 | 77.07 | 79.37 | 80.17 | 81.20
Ours | 75.69 | 78.96 | 79.90 | 80.31 | 81.06 | 82.51

Table 1: Mean AUC result (over the 14 disease classes) for different label set
training percentages. * indicates the methods that use Densenet169 as backbone
architecture. [fengbei]: I have updated NoTeacher result Figure 2: Examples of
classification by our model (trained with 20% labelled data) on Chest X-ray14
[wang2017chestx] images. Red indicates the ground truth label.[fengbei]: Can
we remove these figures, we have to add multiple new tables and there are no
space left

Method | Densenet-121 | GraphXNet [aviles2019graphx] | MOCO V2 [chen2020improved] | MT [liu2020semi] * | SRC-MT [liu2020semi] * | Ours
---|---|---|---|---|---|---
Atelectasis | 75.75 | 71.89 | 77.21 | 75.12 | 75.38 | 78.57
Cardiomegaly | 80.71 | 87.99 | 85.84 | 87.37 | 87.7 | 88.08
Effusion | 79.87 | 79.2 | 81.62 | 80.81 | 81.58 | 82.87
Infiltration | 69.16 | 72.05 | 70.91 | 70.67 | 70.4 | 70.68
Mass | 78.40 | 80.9 | 81.71 | 77.72 | 78.03 | 82.57
Nodule | 74.49 | 71.13 | 76.72 | 73.27 | 73.64 | 76.60
Pneumonia | 69.55 | 76.64 | 71.08 | 69.17 | 69.27 | 72.25
Pneumothorax | 84.70 | 83.7 | 85.92 | 85.63 | 86.12 | 86.55
Consolidation | 71.85 | 73.36 | 74.47 | 72.51 | 73.11 | 75.47
Edema | 81.61 | 80.2 | 83.57 | 82.72 | 82.94 | 84.83
Emphysema | 89.75 | 84.07 | 91.10 | 88.16 | 88.98 | 91.88
Fibrosis | 79.30 | 80.34 | 80.96 | 78.24 | 79.22 | 81.73
Pleural Thicken | 73.46 | 75.7 | 75.65 | 74.43 | 75.63 | 76.86
Hernia | 86.05 | 87.22 | 85.62 | 87.74 | 87.27 | 85.98
Mean | 78.19 | 78.88 | 80.17 | 78.83 | 79.23 | 81.06

Table 2: Class-level AUC comparison between our S2MTS2 and other semi-
supervised SOTA approaches trained with 20% of labelled data. * denotes the
methods that use Densenet-169 as backbone.

Method | Wang et al. [wang2017chestx] | Li et al. [li2018thoracic] | CheXNet [rajpurkar2017chexnet] | CRAL [guan2020multi] | Ma et al. [ma2019multi] | Hermoza et al. [hermoza2020region] | Ours
---|---|---|---|---|---|---|---
Atelectasis | 70 | 72.9 | 75.5 | 78.1 | 77.7 | 77.5 | 78.7
Cardiomegaly | 81 | 84.6 | 86.7 | 88.3 | 89.4 | 88.1 | 87.4
Effusion | 75.9 | 78.1 | 81.5 | 83.1 | 82.9 | 83.1 | 83.8
Infiltration | 66.1 | 67.3 | 69.4 | 69.7 | 69.6 | 69.5 | 70.9
Mass | 69.3 | 74.3 | 80.2 | 83 | 83.8 | 82.6 | 83.3
Nodule | 66.9 | 75.8 | 73.5 | 76.4 | 77.1 | 78.9 | 79.9
Pneumonia | 65.8 | 63.3 | 69.8 | 72.5 | 72.2 | 74.1 | 73.9
Pneumothorax | 79.9 | 79.3 | 82.8 | 86.6 | 86.2 | 87.9 | 87.1
Consolidation | 70.3 | 72 | 72.2 | 75.8 | 75 | 74.7 | 75.9
Edema | 80.5 | 71 | 83.5 | 85.3 | 84.6 | 84.6 | 84.5
Emphysema | 83.3 | 75.1 | 85.6 | 91.1 | 90.8 | 93.6 | 93.7
Fibrosis | 78.6 | 76.1 | 80.3 | 82.6 | 82.7 | 83.3 | 83.4
Pleural Thicken | 68.4 | 73 | 74.9 | 78 | 77.9 | 79.3 | 79.3
Hernia | 87.2 | 66.8 | 89.4 | 91.8 | 93.4 | 91.7 | 93.3
Mean | 74.5 | 73.9 | 78.9 | 81.6 | 81.7 | 82.1 | 82.5

Table 3: Class-level AUC comparison between our S2MTS2 and other supervised
SOTA approaches trained with 100% of labelled data. [fengbei]: Remove fully
supervised detail table? Not enough space

###  4.4 Ablation Study

We study the impact of different components of our proposed S2MTS2 in Tab. 4.
Using the proposed self-supervised learning with just the student model, our
model achieves at least 72.95% mean AUC on various percentages of labelled
training data. Adding the JCL component improves the baseline by around 1%
mean AUC on each training percentage. Adding the mean teacher boosts the
result by 1.5% to 2% mean AUC on each training percentage. [fengbei]: We also
conducted experiments on ISIC2018 and CheXpert datasets to show that our model
works in other datasets. We used the same hyperparameters listed in our paper
except for general parameters (batch size, learning etc…) according to image
size. On ISIC2018, we follow the SOTA method SRC-MT [liu2020semi] setup with
20% labeled and 80% unlabeled split. On SRC-MT’s test set, our method
outperform SRC-MT in all metrics. For CheXpert, we compared our method with
LatentMixing[gyawali2020semi] setup with 100/200/300/400/500k number of
labeled data. Our result is higher than LatentMixing in all labeled subset.
The combination of all our proposed three components achieves SOTA performance
on semi-supervised task.

Self-supervised | JCL | MT | AUC (2%) | AUC (5%) | AUC (10%) | AUC (15%) | AUC (20%)
---|---|---|---|---|---|---|---
✓ |  |  | 72.95 | 76.82 | 78.54 | 79.28 | 80.14
✓ | ✓ |  | 73.60 | 77.46 | 79.18 | 79.83 | 80.62
✓ |  | ✓ | 74.80 | 77.66 | 79.08 | 79.70 | 80.57
✓ | ✓ | ✓ | 75.69 | 78.96 | 79.90 | 80.31 | 81.06

Table 4: Ablation studies of our method with different components. ”Self-
supervised” indicates the traditional self-supervised learning with
contrastive loss [moco]. ”JCL” replaces contrastive loss with (1), ”MT” stands
for fine-tuned with student-teacher learning instead only fine-tuned on only
labelled samples.

Label Percentage (k) | 100 | 200 | 300 | 400 | 500
---|---|---|---|---|---
LatentMixing [gyawali2020semi] | 65.12 | 66.41 | 67.39 | 67.96 | 68.47
Ours | 66.15 | 67.85 | 70.83 | 71.37 | 71.58

Table 5: Mean AUC result (over the 14 disease classes) for different label set
training percentages.

Method | AUC | Sensitivity | F1
---|---|---|---
Supervised | 90.15 | 65.50 | 52.03
MT | 92.96 | 69.75 | 59.10
SRC-MT | 93.58 | 71.47 | 60.68
Ours | 94.71 | 72.14 | 62.67

Table 6: AUC, Sensitivity and F1 result for ISIC 2018 with 20% label and 80%
unlabel.

##  5 Conclusion

In this paper, we presented a novel semi-supervised framework, the Self-
supervised Mean Teacher for Semi-supervised (S2MTS2) learning. The main
contribution of S2MTS2 is the self-supervised mean teacher pre-trained based
on joint contrastive learning [cai2020joint], using an infinite number of
pairs of positive query and key features. This model is then fine-tuned with
the exponential moving average teacher framework. S2MTS2 is validated on the
thorax disease multi-label classification problem from the dataset Chest
X-ray14 [wang2017chestx]. The experiments show that our method outperforms the
previous SOTA semi-supervised learning methods by a large margin in all
benchmarks containing a varying percentage of labelled data. We also show that
the method holds the SOTA results on Chest X-ray14 [wang2017chestx] even for
the fully-supervised problem. The ablation study shows the importance of three
main components of the method, namely self-supervised learning, JCL, and the
mean-teacher model. We will investigate the performance of our method on other
semi-supervised medical imaging benchmarks in the future.


2022

[1]Roel Brouwer

[1]Department of Information and Computing Sciences, Utrecht University,
Princetonplein 5, 3584 CC Utrecht, The Netherlands

# A hybrid local search algorithm for the Continuous Energy-Constrained
Scheduling Problem

<EMAIL_ADDRESS>Marjan van den Akker<EMAIL_ADDRESS>Han Hoogeveen
<EMAIL_ADDRESS>*

###### Abstract

We consider the Continuous Energy-Constrained Scheduling Problem (CECSP). A
set of jobs has to be processed on a continuous, shared resource. A schedule
for a job consists of a start time, completion time, and a resource
consumption profile. We want to find a schedule such that: each job does not
start before its release time, is completed before its deadline, satisfies its
full resource requirement, and respects its lower and upper bounds on resource
consumption during processing. Our objective is to minimize the total weighted
completion time. We present a hybrid local search approach, using simulated
annealing and linear programming, and compare it to a mixed-integer linear
programming (MILP) formulation. We show that the hybrid local search approach
matches the MILP formulation in solution quality for small instances, and is
able to find a feasible solution for larger instances in reasonable time.

###### keywords:

continuous scheduling, resource-constrained scheduling, mixed-integer linear
programming, simulated annealing

##  1 Introduction

When considering scheduling problems with resource constraints, jobs are often
assumed to have a fixed duration or are not allowed to change the amount of
resource they consume over time. In certain applications, however, these
assumptions are too limiting. Typically, such cases can be found in areas
where a (cumulative) amount of work (or resource) is required to complete a
job, but the rate of consumption is not necessarily fixed, such as energy
related applications. Think of demand-response in electricity consumption and
electric vehicle charging, for example. To find schedules that allow for this,
we consider an extension of the traditional scheduling problem, previously
introduced as the Continuous Energy-Constrained Scheduling Problem (CECSP) by
Nattaf . (2014).

This problem is described as follows. A set $\\{J_{1},\ldots,J_{n}\\}$ of jobs
has to be processed on a continuous resource $P$. This means that, at any
time, multiple jobs can be processed simultaneously and at different rates, as
long as their total consumption does not exceed $P$. Each job $J_{j}$ requires
a total amount of resource equal to $E_{j}$. We want to find a schedule such
that: each job $J_{j}$ does not start before its release time $r_{j}$, is
completed before its deadline $\bar{d}_{j}$, and respects its lower and upper
bounds ($P^{-}_{j},P^{+}_{j}$) on the resource consumption rate during
processing. Preemption is not allowed: from its start until its completion,
each job must consume resources at a rate of $P^{-}_{j}$ units. Our objective
is to minimize the total weighted completion time. We look at the case where
both the resource and time are continuous. We assume that there is no
efficiency function influencing resource consumption, and no explicit
precedence relations exist between jobs. It should be noted that the presented
approach can easily be extended to deal with (piece-wise) linear efficiency
functions and precedence relations. The present work, however, does not
include these extensions.

Our contribution. We present a hybrid local search approach (Section 6), using
simulated annealing and linear programming, to solve instances of the CECSP.
In addition, we provide a mixed-integer linear programming (MILP) formulation
(Section 5) and compare the performance of our hybrid approach to it in terms
of quality and runtime. Finally, we found that the feasibility problem is
solvable in polynomial time if we drop the lower bounds on resource
consumption from the problem. We formulate a flow-based solution algorithm for
this case (Section 7.2), that we use to select problem instances for the
computational tests.

As far as we are aware, we are the first to propose a decomposition and try a
combination of local search and mathematical programming techniques for CECSP.
We show that it can compete with exact approaches on small instances, and has
the potential to find good solutions for larger ones. Furthermore, the
approach can be adapted to be applied to related (energy-constrained)
scheduling problems.

The rest of this work is structured as follows. In Section 2, we discuss
related work, followed by a detailed problem description in Section 3. We
explain our event-based modeling approach in Section 4. From there, we build a
MILP formulation in Section 5, and introduce our hybrid local search approach
in Section 6. We discuss how our test instances are generated in Section 7,
which includes a description of our flow-based algorithm for the feasibility
problem without lower bounds. In Section 8 we discuss the results, followed by
the conclusions and suggestions for future work in Section 9.

##  2 Related Work

The CECSP can be seen as a variant of the Resource Constrained Project
Scheduling Problem (RCPSP). This problem concerns the scheduling of activities
subject to precedence and resource constraints. The surveys by Hartmann
Briskorn (2010, 2022) provide a good overview. Baptiste . (1999) defined the
Cumulative Scheduling Problem (CuSP) as a subproblem of the RCPSP. In the case
of CuSP, we have a single resource with a given capacity, and we need to
schedule a number of activities without exceeding this capacity. Each activity
has a release time, deadline, (fixed) processing time and (constant) resource
capacity requirement. More recent work further defined the CECSP as a
generalization of CuSP (Nattaf ., 2014) and provided a hybrid branch-and-bound
algorithm (Nattaf ., 2017) and mixed-integer linear program (Nattaf ., 2019)
to find exact solutions for small instances. In the CECSP, the resource
capacity requirement is a range with a lower and upper bound, rather than a
fixed value, and the consumption rate can vary during the execution of the
job. As a result, the processing time is no longer fixed, but depends on the
consumption rate during its execution. Through its relation to CuSP, Nattaf .
(2016) proved that finding a feasible solution for CECSP is NP-complete.

The CECSP is also closely related to the scheduling of malleable jobs (as
introduced by Turek . (1992)) on parallel machines, which involves the
scheduling of jobs on $P$ machines, while the number of machines assigned to a
job can change during its execution.

Table 1: Example instance of the CECSP with three jobs $j$ | $E_{j}$ | $r_{j}$ | $\bar{d}_{j}$ | $P^{-}_{j}$ | $P^{+}_{j}$ | $w_{j}$ | $B_{j}$
---|---|---|---|---|---|---|---
1 | 70.0 | 0.0 | 3.0 | 10.0 | 30.0 | 1.0 | 0.0
2 | 20.0 | 1.5 | 3.0 | 10.0 | 40.0 | 3.5 | 0.0
3 | 45.0 | 2.5 | 4.0 | 10.0 | 50.0 | 5.0 | 0.0
|  |  |  |  |  |  |
Figure 1: Optimal schedule for the example instance with three jobs

##  3 Problem description

An instance of the CECSP problem is defined by the following properties:

  * •

Available resource: constant $P$;

  * •

For each job $J_{j},$ $j\in\\{1,...,n\\}$:

    * –

Resource requirement $E_{j}$;

    * –

Release time $r_{j}$;

    * –

Deadline $\bar{d}_{j}$;

    * –

Lower bound $P^{-}_{j}$;

    * –

Upper bound $P^{+}_{j}$;

    * –

Linear cost function $w_{j}\cdot C_{j}+B_{j}$.

For each job, we need to determine a start time $S_{j}$, a completion time
$C_{j}$ and a resource consumption profile $p_{j}(t)$ such that
$\sum_{j}w_{j}C_{j}+B_{j}$ is minimal, while the following constraints are
respected:

  * C1

The amount of resource a job $j$ consumes is exactly $E_{j}$;

  * C2

A job $j$ does not start before its release time $r_{j}$, i.e. $S_{j}\geq
r_{j}$;

  * C3

A job $j$ completes no later than its deadline $\bar{d}_{j}$, i.e.
$C_{j}\leq\bar{d}_{j}$;

  * C4

A job $j$ only consumes resources between its start and completion time, i.e.
$p_{j}(t)=0$ for $t<S_{j}$ or $t\geq C_{j}$;

  * C5

While active, the resource consumption of a job $j$ never drops below
$P_{j}^{-}$ or rises above $P_{j}^{+}$, i.e. $P_{j}^{-}\leq p_{j}(t)\leq
P_{j}^{+}$ for $S_{j}\leq t\leq C_{j}$;

  * C6

The total amount of resource consumed by all jobs together at any given time
can never exceed $P$, i.e. $\sum_{j}p_{j}(t)\leq P$ for all possible values of
$t$.

An example instance with three jobs and $P=50.00$ is provided in Table 1.
Figure 1 gives a visual representation of the optimal schedule for this
instance.

Comparing the problem description above to the version originally introduced
by Nattaf . (2014), we notice two differences:

  1. 1.

Our version does not consider efficiency functions (i.e., the amount of work
done is equal to the amount of resource consumed), while in the work of Nattaf
et al., (piece-wise) linear efficiency functions are considered (i.e., the
amount of work done is a linear function of the amount of resource consumed).

  2. 2.

Our objective function minimizes the total weighted completion time, while the
objective in the work of Nattaf et al. is to minimize the amount of resource
consumed.

##  4 Modeling: event-based approach

In this work, we will use an event-based model of the problem. This means that
we define a schedule by determining the timing of events, and the activity of
jobs during the intervals that are bounded by these events.

We consider two types of events, associated with a job $J_{j}$: its start
$S_{j}$ and completion $C_{j}$. Now let us consider a list $\mathcal{E}$, that
lists all events in the order that they occur over time.

For a given order of events, we can divide the schedule into $2n-1$ intervals,
where each interval is bounded by two consecutive events. In the following,
when we refer to an interval, we always mean the span of time in between two
consecutive events. So, no event ever occurs during an interval. Also note
that an interval can be empty, if it is bounded by two events that happen at
the same time.

During these intervals, the set of jobs that are being processed does not
change. We can show that there always exists an optimal solution in which
$p_{j}(t)$ remains constant during each interval. We can therefore limit
ourselves to such solutions.

######  Theorem 1.

For any feasible schedule $\mathcal{S}$ (with start times $S_{j}$, completion
times $C_{j}$ and resource consumption profiles $p_{j}(t)$ for all jobs
$J_{j},j\in\\{1,...,n\\}$) that follows a given event order $\mathcal{E}$, a
feasible schedule $\mathcal{S}^{\prime}$ exists with the same objective value
where the resource consumption $p_{j}^{\prime}(t)$ of all jobs remains
constant during each interval.

Proof: Consider any interval between two events, $i$ and $i^{\prime}$. The
interval starts at $t_{i}$ and ends at $t_{i^{\prime}}$.

A job $J_{j}$ is inactive during this interval if $t_{i}\geq C_{j}$ or
$t_{i^{\prime}}\leq S_{j}$ in the schedule, and active otherwise.

The resource consumption profile of any active job $J_{j}$ in $\mathcal{S}$
within this interval does not violate the lower and upper bounds for this job:
$P^{-}_{j}\leq p_{j}(t)\leq P^{+}_{j}\ \forall t,t_{i}\leq t\leq
t_{i^{\prime}}$.

Therefore, the average consumption of job $J_{j}$ within this interval,
$\bar{p}_{j}=\frac{\int_{t_{i}}^{t_{i^{\prime}}}p_{j}(t)\
dt}{t_{i^{\prime}}-t_{i}}$, does not exceed the bounds on consumption for job
$J_{j}$: $P^{-}_{j}\leq\bar{p}_{j}\leq P^{+}_{j}$.

In this interval, all active jobs in $\mathcal{S}$ together consume at most
($t_{i^{\prime}}-t_{i})P$ of the resource, i.e.:

$\sum_{j}\bar{p}_{j}(t_{i^{\prime}}-t_{i})\leq(t_{i^{\prime}}-t_{i})P$

From this, it follows that $\sum_{j}\bar{p}_{j}\leq P$, i.e., the sum of the
average consumption of all jobs does not exceed the available amount of
resource.

We obtain $\mathcal{S}^{\prime}$, from $\mathcal{S}$ by replacing all resource
consumption profiles $p_{j}(t)$ with a step-wise constant function
$p_{j}^{\prime}(t)$ that remains constant during any interval (in between two
consecutive events), and takes the value $\bar{p}_{j}$, i.e. the average
resource consumption of job $J_{j}$ in $\mathcal{S}$, during that interval.

We have shown that the schedule $\mathcal{S}^{\prime}$ is indeed feasible. For
each job, the resulting consumption rate is within the bounds allowed for that
job and the total amount of available resource is not exceeded either. No
other constraints are affected.

It remains to show that $\mathcal{S}^{\prime}$ has the same objective value.
The objective depends only on the value of $C_{j}$. As the timing of all
events in $\mathcal{S}^{\prime}$ is identical to that in $\mathcal{S}$, the
objective is not affected by the modifications to the schedule.

######  Corollary 2.

For any optimal schedule $\mathcal{S}^{*}$, there exists an optimal schedule
${\mathcal{S}^{*}}^{\prime}$ where the resource consumption
${p_{j}^{*}}^{\prime}(t)$ of all jobs remains constant during each interval.

This means that we can ignore solutions with a non-constant resource
consumption profile for any job during any interval. Therefore, we can fully
represent a resource consumption profile of job $J_{j}$ with at most $2n-1$
values $p_{j,i}$, one for each interval in which $J_{j}$ is being processed,
indicating the total amount of resource consumed during that interval.

In the following, we will maintain a link between events and jobs by their
index. Every job has two associated events. For a job with index $j$, we will
say that its start event has index $i=2j-1$, and its completion event has
index $i+1=2j$. This allows us to uniquely identify the start and completion
events of any job easily in notation. _The index of an event $i$ does not
indicate its position in the event order._ We will maintain this indexing
scheme throughout this work.

Furthermore, intervals are defined by the event at the start of that interval.
For example, let $i$ and $i^{\prime}$ be the indices of the first two events.
The first interval starts when the first event happens: at $t_{i}$; and ends
when the next event starts: at $t_{i^{\prime}}$. Then, $p_{j,i}$ represents
the amount of resource that job $J_{j}$ consumes in the interval associated
with event $i$, spanning the time between $t_{i}$ and $t_{i^{\prime}}$.

##  5 Mixed-integer linear programming formulation

Using the event-based approach, we developed a MILP formulation that will find
an optimal schedule. In this section, we gradually build this formulation.

Although we arrived at it independently, our formulation is similar to the
event based formulations presented by Nattaf . (2016). The models share their
focus on intervals, delimited by events, but differ in how the order of these
events is expressed. Nattaf . (2016) use assignment variables for this: binary
variables that indicate whether events occur at a given position in the order,
or whether a job is active during a given interval. Our model uses binary
ordering variables, to indicate the relative order between pairs of events.

First, we define the decision variables of our model. For each of the $2n$
events $i$, $t_{i}\geq 0$ will denote the time at which it takes place. For
each job $J_{j}$ and event $i$, $p_{j,i}\geq 0$ is the amount of resource that
$J_{j}$ consumes during the interval defined by event $i$. We know the
consumption of a job $J_{j}$ during the interval it completes to be 0 and
therefore set $p_{j,2j}=0$. Note that we do not know beforehand what event
will be the last in the event order. The interval that is defined by that
event is, theoretically, open-ended. However, in the formulations below you
will note that any $p_{j,i}$ associated with this interval is forced to 0, as
all jobs complete in or before the interval.

Additionally, we define binary variables $a_{i,i^{\prime}}$ and
$b_{i,i^{\prime}}$ to determine the order of events. These can be viewed as
helper variables. For two events $i$ and $i^{\prime}$, $a_{i,i^{\prime}}$
represents the order they occur in. If $a_{i,i^{\prime}}=1$, then $i$ occurs
before $i^{\prime}$ in the event order. On top of that, $b_{i,i^{\prime}}$
indicates whether an event is an immediate successor of another:
$b_{i,i^{\prime}}=1$ if $i^{\prime}$ occurs immediately after $i$. In our
formulation, we will use $a_{i,i^{\prime}}$ where possible and
$b_{i,i^{\prime}}$ if we have to: exclusively for the enforcement of lower
bounds.

The objective is to minimize the combined cost functions of the jobs:

$\min\sum\limits_{j=1}^{n}\left(w_{j}t_{2j}+B_{j}\right)$ (1)

Now, let us go through the constraints we listed in Section 3 above.

First, we ensure that the total amount of resource consumed by each job
$J_{j}$ equals $E_{j}$ (C1), $\forall j\in\\{1,...,n\\}$:

$\sum\limits_{i=1}^{2n}p_{j,i}=E_{j}$ (2)

Next, we enforce the release time (C2) and deadline (C3) for each job $J_{j}$.
Since we can identify the start and completion event of a job by its index
$i$, we obtain that $\forall j\in\\{1,...,n\\}$:

$t_{2j-1}\geq r_{j}$ (3) $t_{2j}\leq\bar{d}_{j}$ (4)

We ensure that $p_{j,i}=0$ for any interval after the completion event of
$J_{j}$ (with index $2j$) has happened or before its start event (with index
$2j-1$) has taken place (C4). We use a big $M$ constraint with the order
variables $a_{i,i^{\prime}}$ to enforce this exactly for those intervals. So
$\forall j\in\\{1,...,n\\},i\in\\{1,...,2n\\}$, and sufficiently large $M$
($M=E_{j}$ will suffice) we have:

$p_{j,i}\leq a_{i,2j}M$ (5) $p_{j,i}\leq a_{2j-1,i}M$ (6)

In between the start and completion events of a job $J_{j}$, we need to choose
values for $p_{j,i}$ that respect the lower ($P^{-}_{j}$) and upper
($P^{+}_{j}$) bounds on resource consumption (C5). The upper bound can be
enforced in between any two events. We use a big $M$ constraint with order
variable $a_{i,i^{\prime}}$ to activate the constraint only for events that
happen in the right chronological order, hence $\forall
j\in\\{1,...,n\\},i,i^{\prime}\in\\{1,...,2n\\},i\neq i^{\prime}$, and
sufficiently large $M$:

$p_{j,i}\leq P_{j}^{+}(t_{i^{\prime}}-t_{i})+Ma_{i^{\prime},i}$ (7)

The lower bound only applies for consecutive events, which both fall within
the execution window of job $J_{j}$. We use a big $M$ term with the
$b_{i,i^{\prime}}$ helper variable to activate the constraint only for events
that are immediate successors. Besides that, we include the variables
$a_{i^{\prime},2j-1}$ and $a_{2j,i}$ to make sure the constraint is only
active if $i$ and $i^{\prime}$ are between the start event of job $J_{j}$
(index $2j-1$) and its completion (index $2j$). Hence $\forall
j\in\\{1,...,n\\},i,i^{\prime}\in\\{1,...,2n\\},i\neq i^{\prime}$, and
sufficiently large $M$:

$\displaystyle p_{j,i}\geq$ $\displaystyle P_{j}^{-}(t_{i^{\prime}}-t_{i})$
(8) $\displaystyle-(1-b_{i,i^{\prime}}+a_{i^{\prime},2j-1}+a_{2j,i})M$

During an interval the total resource consumption of all busy jobs cannot
exceed $P$ at any time. This has to be enforced for all intervals, $\forall
i,i^{\prime}\in\\{1,...,2n\\},i\neq i^{\prime}$, and sufficiently large $M$:

$\sum\limits_{j=1}^{n}p_{j,i}\leq P(t_{i^{\prime}}-t_{i})+Ma_{i^{\prime},i}$
(9)

We need to make sure that the event times are consistent with the variables
$a_{i,i^{\prime}}$ that control the order of events. Therefore, $\forall
i,i^{\prime}\in\\{1,...,2n\\},i\neq i^{\prime}$, and sufficiently large $M$,
we get:

$t_{i}\leq t_{i^{\prime}}+Ma_{i^{\prime},i}$ (10)

We ensure that either $a_{i,i^{\prime}}=1$ or $a_{i^{\prime},i}=1$, but not
both, for any pair of distinct events. We could, of course, get rid of half of
our binary variables by expressing $a_{i^{\prime},i}$ in terms of
$a_{i,i^{\prime}}$ for all $i^{\prime}>i$. We choose not to do so for ease of
notation, therefore we need to ensure that $\forall
i,i^{\prime}\in\\{1,...,2n\\},i<i^{\prime}$:

$a_{i,i^{\prime}}+a_{i^{\prime},i}=1$ (11)

We can note here that we also already know: $\forall j\in\\{1,...,n\\}$,
$a_{2j-1,2j}=1$.

Finally, we need some constraints to make sure the $b_{i,i^{\prime}}$
variables indicate events that are immediate successors. First making sure
that $b_{i,i^{\prime}}$ can only be 1 if $i^{\prime}$ occurs immediately after
$i$, $\forall i,i^{\prime}\in\\{1,...,2n\\},i\neq i^{\prime}$, and
sufficiently large $M$:

$\sum\limits_{i^{\prime\prime}=1}^{2n}a_{i,i^{\prime\prime}}-\sum\limits_{i^{\prime\prime}=1}^{2n}a_{i^{\prime},i^{\prime\prime}}\leq
1+(1-b_{i,i^{\prime}})M$ (12)
$\sum\limits_{i^{\prime\prime}=1}^{2n}a_{i,i^{\prime\prime}}-\sum\limits_{i^{\prime\prime}=1}^{2n}a_{i^{\prime},i^{\prime\prime}}\geq
1-(1-b_{i,i^{\prime}})M$ (13)

And second, by requiring that exactly (or equivalently: at least) $2n-1$
variables $b_{i,i^{\prime}}$ actually take the value of 1:

$\sum_{i=1}^{2n}\sum_{i^{\prime}\neq i}b_{i,i^{\prime}}=2n-1\\\ $ (14)

This completes the initial MILP formulation. However, we have added two types
of valid inequalities to speed up the process of solving the model. We know
that the range of feasible processing times for a job is limited by the lower
and upper bound on its resource consumption, $\forall j\in\\{1,...,n\\}$:

$t_{2j}-t_{2j-1}\leq E_{j}/P^{-}_{j}$ (15) $t_{2j}-t_{2j-1}\geq E_{j}/P+_{j}$
(16)

Indeed, the addition of these inequalities has a significant impact on the
time needed to find an optimal solution. The MILP formulation that our hybrid
local search approach will be tested against consists of all equations (1) -
(16).

##  6 Hybrid local search

In the MILP formulation described above, we can identify two separable parts
to the problem. The binary variables $a_{i,i^{\prime}}$ (and
$b_{i,i^{\prime}}$) fully determine the order of events $\mathcal{E}$. Then,
the continuous variables $t_{i}$ and $p_{j,i}$ determine the exact schedule
for this order.

We can decompose the problem into these two parts. Once we have an order of
events $\mathcal{E}$, determining an optimal schedule boils down to solving an
LP, as all binary variables disappear from the formulation once the event
order is fixed.

A global overview of our approach is given in Algorithm 1. The elements in
this overview will be discussed in more detail in the following sections.

Algorithm 1 Hybrid local search approach

1:$prec\leftarrow$ ImplicitPrecedences($J$)

2:$\mathcal{E},best\\_score\leftarrow$ InitialSolution($J$, $P$, $prec$)

3:$T\leftarrow T_{init}$

4:for $iter\in\\{1,...,max\\_iter\\}$ do

5: $\mathcal{E},score\leftarrow$ GetNeighbor($\mathcal{E}$, $T$, $J$, $P$,
$prec$)

6: if $score<best\\_score$ then

7: $\mathcal{E}_{best},best\\_score\leftarrow\mathcal{E},score$

8: end if

9: if $iter\%\alpha_{period}=0$ then

10: $T\leftarrow T\cdot\alpha$

11: end if

12:end for

13:return $\mathcal{E}_{best},best\\_score$

###  6.1 Finding an optimal schedule for a given order: linear program

We will adapt the MILP formulation described above to solve only the
(sub)problem of finding an optimal schedule for a given event order.

The LP will consist of Equation (1) - (4) and (7) - (10), eliminating any part
of the inequalities that contain the variables $a_{i,i^{\prime}}$ and
$b_{i,i^{\prime}}$.

Instead of enforcing $p_{j,i}$ to be 0 outside of $J_{j}$’s processing window
(as in Equation (5) and (6)), we only take into account $p_{j,i}$ for $i$
within the processing window, which can easily be done as we already know the
order of events. Similarly, upper and lower bounds are only enforced on
intervals within the processing window.

For ease of notation, let us define two functions:

  * •

$I(e)$ gives the event index $i$ of the event in position $e$ of the event
order;

  * •

$E(i)$ gives the position in the event order $e$ of the event with index $i$.

Note that these functions are known already before we attempt to solve the LP,
as they only depend on the event order.

So, for example the inequality enforcing the order of events (Equation (10) in
the MILP formulation) now becomes, $\forall e\in\\{1,...,2n-1\\}$:

$t_{I(e)}\leq t_{I(e+1)}$ (17)

And the inequality controlling the total work of a job (Equation (2) in the
MILP formulation) now becomes, $\forall j\in\\{1,...,n\\}$:

$\sum\limits_{e\in\\{E(2j-1),...,E(2j)-1\\}}p_{j,I(e)}=E_{j}$ (18)

Note that $\\{E(2j-1),...,E(2j)-1\\}$ represents the set of positions in the
event order between the start and completion event of a job $J_{j}$, including
the former and excluding the latter.

Finally, to evaluate the relative quality of infeasible event orders, we allow
for the violation of lower bounds, upper bounds, and/or resource availability
constraints, with appropriate cost. For this purpose, we introduce three types
of slack variables:

  * •

$s^{-}_{j,i}$ allowing the resource assignment to job $J_{j}$ in interval $i$
to go below its lower bound $P^{-}_{j}$;

  * •

$s^{+}_{j,i}$ allowing the resource assignment to job $J_{j}$ in interval $i$
to go over its upper bound $P^{+}_{j}$;

  * •

$s^{t}_{i}$ allowing for increasing the amount of available resource during an
interval $i$.

This affects the inequalities enforcing these bounds in the LP (Equation (7) -
(9) in the MILP formulation), $\forall
j\in\\{1,...,n\\},e\in\\{E(2j-1),...,E(2j)-1\\}$:

$p_{j,I(e)}\geq P_{j}^{-}(t_{I(e+1)}-t_{I(e)})-s^{-}_{j,I(e)}$ (19)
$p_{j,I(e)}\leq P_{j}^{+}(t_{I(e+1)}-t_{I(e)})+s^{+}_{j,I(e)}$ (20)

And, $\forall e\in\\{1,...,2n-1\\}$:

$\sum\limits_{j=1}^{n}p_{j,I(e)}\leq P(t_{I(e+1)}-t_{I(e)})+s^{t}_{I(e)}$ (21)

These slack variables are penalized in the objective function, using $L^{B}$
for violation of upper and lower bounds and $L^{R}$ for excessive resource
usage. The slack variables are non-negative, but do not have an enforced upper
bound. This adds the following term to the objective:

$\sum\limits_{i=1}^{2n-1}\left(L^{R}s^{t}_{i}+\sum\limits_{j=1}^{n}L^{B}(s^{-}_{j,i}+s^{+}_{j,i})\right)$
(22)

This allows us to get an idea of the quality of infeasible event orders that
we can use to guide the local search towards solutions that are closer to
feasibility. For many event orders, no feasible schedule exists. Scoring the
‘degree of infeasibility’, helps us to get from an infeasible order to a
feasible one.

###  6.2 Finding a good event order: local search

It remains to find good event orders to feed to the LP. For this, we use a
local search algorithm: simulated annealing. Below, we will describe the core
of the local search, and the neighborhood operators that are used.

We start our local search with the event order that results from the greedy
algorithm (see below). Every iteration, we modify our current solution using
one of the neighborhood operators. Improvements are always accepted, while
solutions that are worse than the current one are accepted with a probability
determined by the current value of the temperature $T$. The precise settings
for the simulated annealing will be discussed in Section 8.

If a candidate solution is rejected, we keep trying with the same neighborhood
operator. We do this by generating a random permutation of the event list each
time we start with a neighborhood operator, and apply the operator to the
events in that order, until a candidate solution is accepted. If this is
unsuccessful we move on to the next neighborhood operator. If no solution is
accepted after looping through all three neighborhood operators, we terminate
the search and report the best seen solution up to that point.

We have three neighborhood operators that we use to modify the event order and
obtain new candidate solutions:

  * •

Swap two adjacent events. We take the next event from the random permutation,
at position $e\in\\{1,...,2n-1\\}$ and its successor at position $e+1$ and
reverse their order. This means that the event originally at position $e$ now
is at position $e+1$ and vice versa.

  * •

Move a single event. We take the next event from the random permutation, at
position $e\in\\{1,...,2n\\}$. We then determine its range of movement. On
both sides, it is limited by the first event we encounter with which the event
has a precedence relation. How we determine these precedence relations is
discussed in Section 6.3. From this range, we then select a new position,
where the probability of a position being selected is inversely proportional
to its distance to the original position of the event. That is, smaller
displacements are more likely to be selected than larger ones. The event is
then moved to its new position, with the other events in between the current
and former position of the event shifted accordingly.

  * •

Move both events associated with a single job. We take the next job from the
random permutation, job $J_{j}$ $j\in\\{1,...,n\\}$. We then look at both of
its associated events, and determine their combined range of movement. As
before, we determine the number of positions both can move in either direction
without crossing over an event with which they have a precedence relation (see
Section 6.3). We then take the minimum of these values to get a combined
range. From this range, we select a new position with uniform probability. The
events are then moved to their new position, both offset from their original
position by the same amount, with the other events in between the current and
former position of the events shifted accordingly.

Each of these neighborhood operators results in a small modification of the
event order, necessitating an update of the LP for evaluation. We modify only
those parts of the LP that are affected, and resolve the model to obtain the
score of the new candidate solution, starting from the current solution.

###  6.3 Preprocessing and initial solution: greedy algorithm

We obtain an initial event order using a greedy algorithm. The algorithm
assigns resource to jobs that have been released and have not been fully
served, and extracts the order of events from the resulting schedule. In the
construction of the initial order, we ignore lower bounds on resource
consumption and do not consider the cost function at all. Keep in mind that we
are only interested in the order of events at this point, not in the schedule
itself.

We take the release times $r_{j}$ and deadlines $\bar{d}_{j}$ of all jobs, and
sort them in ascending order. In this way, we obtain a sequence of time
periods during which the set of jobs available for processing remains the
same.

We loop over these time periods, keeping track of the list of available jobs
during the period. Then, we decide the assignment of resources during each
time period in the following way:

  1. 1.

We assign all available jobs their minimal amount of resource required. If
this exceeds the total amount of resource available, we stop after this step
and continue to the next time period;

  2. 2.

In order of increasing deadline, we assign as much resource to each job as it
can consume. This may be bounded by the upper bound on resource consumption of
the job, or by the amount of resource that is left to be assigned at this
point. We continue until the resource in the time period is fully used or no
other jobs are available.

The minimal amount of resource a job has to consume is given by
$\max\\{0,\tilde{E}_{j}-(\bar{d}_{j}-t_{end})P^{+}_{j}\\}$, while its maximum
is given by $\min\\{\tilde{E}_{j},P^{+}_{j}(t_{end}-t_{start})\\}$. Here
$t_{start}$, $t_{end}$ are the start and end of the current time period and
$\tilde{E}_{j}$ is the _remaining_ resource requirement of the job.

At any time during this procedure:

  * •

If a job is assigned resources for the first time, we append the start event
for that job to our current list;

  * •

If a job’s full resource requirement is satisfied, we append the completion
event for that job to our current list.

At the end of the procedure we have a complete order of events, which will be
the starting solution for our local search. We determine the quality of the
optimal schedule for this initial event order by solving the LP above. This
schedule is likely to assign resources differently from the greedy algorithm
we used to find the event order in the first place. Note that the schedule may
have positive penalty terms, i.e. it does not have to meet the available
resource bound $P$ or the job lower or upper bounds $P_{j}^{-}$ and
$P_{j}^{+}$.

To complete the description of our approach, we will discuss the preprocessing
steps that we use in our algorithm. We deduce a number of implicit precedence
constraints that we use in the neighborhood operators:

  1. 1.

A start event of a job $S_{j}$ always has to occur before its corresponding
completion event $C_{j}$;

  2. 2.

Let $u_{j}=E_{j}/P^{+}_{j}$ be the minimal processing time required to
complete a job. Then the possible time ranges of a start and completion event
are $[r_{j},\bar{d}_{j}-u_{j}]$ and $[r_{j}+u_{j},\bar{d}_{j}]$, respectively.
For any two events, if the last possible time within the range of event $i$
lies before the first possible time within the range of event $i^{\prime}$, we
know that $i$ has to come before $i^{\prime}$ in any event order.

##  7 Test instances

We generated a number of instances to evaluate our approach. We consider three
important characteristics to group these instances: the number of jobs ($n$),
the amount of resource ($P$) and whether an instance is adversarial ($a$).

An instance is said to be adversarial, when the ordering of weights is such
that the weight of jobs monotonically increases with increasing deadline. This
is meant as a counter to the greedy algorithm that determines a starting
solution, as this assigns resources to jobs in order of increasing deadline.

By varying $P$, we also affect the ratios of $P/P^{-}_{j}$, $P/P^{+}_{j}$ and
$P/E_{j}$, as the sampling of values for $E_{j}$, $P^{+}_{j}$ and $P^{-}_{j}$
is (largely) unaffected by the chosen values of $n$, $P$ and $a$. Each of
these ratios are more interesting properties of an instance than $P$ in
itself.

Table 2 lists the values of $n$, $P$ and $a$ for which instances were
generated. Each combination of $n$, $P$ and $a$ was used to generate four
unique instances.

Below we describe the procedure for generating instances.

Table 2: Overview of values of $n$, $P$ and $a$ used for generating instances $n$ | 5, 10, 15, 20, 30, 50
---|---
$P$ | 25.00, 50.00, 100.00, 200.00
$a$ | True, False
|

###  7.1 Instance generation method

The instance generation is based on the three parameters discussed above ($n$,
$P$ and $a$) as well as four scaling parameters: $a_{\mathrm{maxlow}}$,
$a_{\mathrm{minupp}}$, $a_{\mathrm{rshift}}$ and $a_{\mathrm{pws}}$. For each
job we sample seven values, that together fully determine the six properties
of a job as described in Section 3.

  * •

A resource requirement $E_{j}$ from $U(10,100)$.

  * •

A lower bound $P^{-}_{j}$ from
$U(0,\min\\{a_{\mathrm{maxlow}}P,a_{\mathrm{minupp}}E_{j}\\})$. The parameters
$a_{\mathrm{maxlow}}$ and $a_{\mathrm{minupp}}$ can be used to control the
largest possible lower bound. Here $a_{\mathrm{maxlow}}$ makes sure that the
lower bound on the resource consumption per time unit never exceeds a fraction
($a_{\mathrm{maxlow}}$) of what is available during a unit of time, leading to
instances where the lower bounds do not restrict the schedule to only
processing a single job at a time. Then, $a_{\mathrm{minupp}}$ has a similar
purpose, but as a fraction of the total resource requirement of a job,
ensuring that the lower bound does not restrict the possible duration of a job
too much.

  * •

An upper bound $P^{+}_{j}$ from $U(a_{\mathrm{minupp}}E_{j},E_{j})$. The
$a_{\mathrm{minupp}}$ parameter can be used to make the upper bound more
likely to be more restrictive, while also ensuring it does not drop below the
lower bound. This distribution also ensures that no job can be completed in
less than a single time unit.

  * •

A release time $r_{j}$ from
$U\left(-a_{\mathrm{rshift}}\frac{\sum_{j}E_{j}}{P},(1-a_{\mathrm{rshift}})\frac{\sum_{j}E_{j}}{P}\right)$.
Any negative values are treated as zero. The parameter $a_{\mathrm{rshift}}$
shifts the window for generating release times, such that a fraction of the
generated release times equal to $a_{\mathrm{rshift}}$ is expected to be a
negative number. In this way, one can control the expected number of jobs that
are available for processing immediately at time 0. The largest possible
release time is limited by the expected total processing time (at full
utilization).

  * •

a deadline $\bar{d}_{j}$:
$r_{j}+U\left(\frac{E_{j}}{\min\\{P,P^{+}_{j}\\}},a_{\mathrm{pws}}\frac{\sum_{j}E_{j}}{P}\right)$.
The parameter $a_{\mathrm{pws}}$ sets the latest possible deadline in terms of
the expected total processing time (at full utilization). The earliest
possible deadline is the release time plus the minimum required processing
time (its requirement divided by the minimum of its upper bound or the amount
of resource available per unit of time), avoiding trivially infeasible
instances.

  * •

A weight $w_{j}$ from $U(0,5)$ and an optional objective constant $B_{j}$ from
$U(0,10)$. Together they determine the cost function of job $J_{j}$:
$w_{j}\cdot C_{j}+B_{j}$.

All numbers are rounded to two decimal places.

This instance generation method does not entirely prevent that there might be
a shortage of available resource during a certain period of time, depending on
the sampled release times and deadlines. The odds of this happening can be
reduced by increasing the range of release times and/or processing windows,
using $a_{\mathrm{pws}}$ for example. We describe a method that gives a strong
indication of the feasibility of in instance in the following section. In
Section 7.3 we describe our choice of parameters that lead to mostly feasible
instances.

###  7.2 Checking feasibility

To get an indication of the (in)feasibility of generated instances, we use a
flow-inspired LP that solves the feasibility problem for the instance,
ignoring lower bounds ($P^{-}_{j}$).

The lower bounds might still render an instance infeasible, even if the LP
indicates that a feasible flow exists. However, it provides a good and quick
indication.

Figure 2: Example with three jobs and four intervals

We will formulate the problem as a max flow problem. For this, we construct a
bipartite graph, with one node for each job on one side, and one node for each
time interval on the other. To obtain relevant intervals, we sort the set of
all release times and deadlines and get an interval for every two consecutive
distinct values. This results in at most $2n-1$ intervals.

There is an arc between a job node ($j$) and a time interval node ($t$), if
the job can be processed during the interval $t=[t_{s},t_{e})$ (i.e.
$r_{j}\leq t_{s}$ and $t_{e}\leq\bar{d}_{j}$). The capacity of that arc is
equal to $(t_{e}-t_{s})P_{j}^{+}$, the maximum amount of power that job $j$
can receive during that interval.

We add a source that has an arc going to every job node, with $E_{j}$ (the
total amount of work for that job) as its capacity, and a sink with an
incoming arc from each time interval node with a capacity of $(t_{e}-t_{s})P$
(the available power in that interval).

A flow through this graph now will represent a power assignment to the jobs.
We use an LP analogue to solve this flow problem. If a max flow of
$\sum_{j}E_{j}$ is found, we know the problem without lower bounds is
feasible. If the max flow has a smaller value, we know that this instance has
no feasible solution.

An example with three jobs and four time intervals is illustrated in Figure 2.

###  7.3 Choice of parameters for instance generation

We performed a grid search on a number of possible combinations of parameter
settings, generating 100 random instances for each setting and using the flow-
inspired LP described above to check the feasibility.

This analysis was performed for different values of $n$, while $P$ was fixed
at $50.00$. A parameter setting was selected for each value of $n$ such that
for 99% of generated instances a feasible flow exists in the relaxed problem
(without lower bounds). The aim was to end up with mostly feasible instances,
with the occasional infeasible instance mixed in.

Note, however, that this analysis was performed for $P=50.00$ only. While the
same parameter settings were used for other values of $P$, the instances
generated are not necessarily equally likely to be feasible. We observe in
practice, however, that the likelihood of feasibility is very similar.

The parameter settings used for each value of $n$ are listed in Table 3.

Table 3: Parameter settings used for generating instances of different sizes | $a_{\mathrm{maxlow}}$ | $a_{\mathrm{minupp}}$ | $a_{\mathrm{rshift}}$ | $a_{\mathrm{pws}}$
---|---|---|---|---
$n=5,10$ | 0.25 | 0.25 | 0.125 | 2
$n=15,20,30,50$ | 0.25 | 0.25 | 0.125 | 1.5
|  |  |  |

##  8 Computational results

We ran computational experiments, measuring both solution time and the
objective value of the best found solution for the MILP formulation and the
hybrid local search algorithm. The instances we generated and used to run our
tests can be found here: https://github.com/RoelBrouwer/2022-cecsp-
data/releases/tag/jos and the code used to run the experiments is available
here: https://github.com/RoelBrouwer/continuousresource/releases/tag/jos.

###  8.1 Parameter settings

To run our simulated annealing algorithm, we need to determine the value for
the following parameters:

  * •

Initial temperature $T_{init}$;

  * •

Number of iterations between temperature updates $\alpha_{period}$;

  * •

Multiplication factor for updating the temperature $\alpha$;

  * •

Penalty for using slack variables $L^{B},L^{R}$;

  * •

Probability of selecting each neighborhood operator:

    * –

Swap $p_{s}$;

    * –

Single move $p_{m}$;

    * –

Paired move $p_{p}$.

For smaller values of $n$, we ran a (limited) grid search for a number of
combinations of parameter values, selecting values that performed well in
these instances, and extrapolated this to larger instance sizes. The selected
values are presented in Table 4.

In addition, since we consider the hybrid local search to be most relevant for
the instances with $n=50$, we allow restarts under certain conditions for
these larger instances. If a run completes before 1800 seconds have passed, we
perform 100 random swaps on the current event order, and restart the simulated
annealing with the obtained order as the initial solution.

Table 4: Parameter settings for simulated annealing used for computational experiments $T_{init}$ | $\alpha_{period}$ | $\alpha$ | $L^{B}$ | $L^{R}$ | $p_{s}$ | $p_{m}$ | $p_{p}$
---|---|---|---|---|---|---|---
$n$ | $(2n-1)\cdot 4$ | 0.95 | 5 | 5 | 0.75 | 0.15 | 0.1
|  |  |  |  |  |  |

For solving both the MILP formulation and the LP subproblem of the simulated
annealing algorithm we used CPLEX Studio 22.1, with a time limit of 3600
seconds, and limited to the use of a single thread. Everything else has been
implemented in the Python programming language. The processor of the system
used to run the tests on was an Intel(R) Xeon(R) Gold 6130 CPU @ 2.10GHz.

Table 5: Comparison of MILP formulation and simulated annealing algorithm on small instances | MILP | SA
---|---|---
$n$ | Avg. time | Avg. diff | Avg. time | Avg. diff
5 | 2.23 | 0.00 | 85.63 | 0.00
10 | 2300.45 | 0.00 | 280.18 | 0.00
15 | 3521.22 | 0.03 | 532.74 | 0.00
|  |  |  |

###  8.2 Results

The full results of these runs can be found in Tables 7 \- 12.

In these result tables, the first four columns uniquely describe the instance:
the number of jobs $n$, resource availability $P$, whether the weights are
sorted in ascending order (adversarial). As four instances of each type were
generated, the fourth column identifies each instance with a number 0-3. We
then list whether our feasibility test was passed by the instance or not (as
described in Section 7.2).

This is followed by a column indicating the best objective value that we have
encountered in any run of any of our algorithms. Note that this include some
extremely long runs done in the trail phase. Next we list the results of the
MILP and hybrid local search algorithm (SA). For each we present the run time
and best found objective value, while we also add the objective value of the
initial solution (found by the greedy algorithm described in Section 6.3) for
our hybrid local search approach.

Any cursive values indicate objective values corresponding to infeasible
solutions. In the case of a solution found by our hybrid local search, this
means that some slack variables have a non-zero value. An objective value of a
solution found by the MILP or hybrid local search is put in boldface if it is
equal to the best known objective value. Finally, we put LIMIT in the time
column under MILP if the solver did not terminate within 3600 seconds, meaning
it was unable to prove optimality of the current best found solution.

Note that the results for the hybrid local search in these tables are of a
single run of each of the algorithms, they are not averages or best-of
results.

In Table 5 we compare the performance of the MILP formulation and the
simulated annealing on small instances ($n=5,10,15$) in terms of run time and
solution quality. This analysis is based on the results in Tables 7 \- 9.
Infeasible instances were left out of the analysis. The average runtime and
difference are computed only over those instances where both approaches come
up with a feasible solution. Note that the average runtimes of the MILP for
$n=10,15$ are somewhat misleading, as many runs have been cut off at 3600
seconds, resulting in a runtime equal to the cut-off time for these runs. The
difference is defined to be the relative difference between the score of the
best known solution and the reported solution.

It is relevant to point out that for both $n=10$ and $n=15$ there are two
instances where the MILP is able to find a feasible solution, and the
simulated annealing is not, while for $n=15$ the simulated annealing finds a
feasible solution for five instances where the MILP fails.

We can see that the runtime of the hybrid local search approach scales better
than that of the MILP. Our approach is competitive with the MILP on smaller
instances: it finds equally good solutions, only taking more time for the
smallest instances. On larger instances, the hybrid local search finds better
solutions than the MILP in less time, if the MILP is able to find any feasible
solution at all. The difference is likely to become even clearer if we would
perform a number of reruns of the local search and only keep the best result,
as we already did for the larger instances ($n=50$). For $n=15$, for example,
about 6-7 reruns can be performed without exceeding the time limit of 3600
that the MILP also has to respect.

We do see, however, that the larger instances are much more challenging. In
Table 6, we show the number of instances for which we summarize the number of
instances for which the (flow) feasibility test predicted a feasible solution
exists, and present the number of times a feasible solution and the best known
solution were found by our hybrid approach and the MILP formulation. This
analysis is based on the results in Tables 7 \- 12. By allowing restarts, we
are able to find feasible solutions in the majority of cases, even for $n=50$,
but tackling larger instance sizes remains a challenge.

Note that the best found solution for $n=30$ is always the same as the
reported solution for the hybrid local search algorithm. We added the
instances of this size at a later stage to gain better insight, as the sizes
of $n=20$ and $n=50$ differ strongly in their results. We did not do trail
runs for these instances, and the reported solutions are therefore by default
always the best known. For the other large instances ($n=20,50$), we observe
that the number of times that our approach finds the best known solution
strongly decreases. To put this in perspective, we note that the best known
solution for these larger instances was often found during a trail run that
took many hours. These do not offer a realistic alternative to the presented
approach.

The initial solution that the greedy algorithm from Section 6.3 finds,
provides a good starting point for our local search. For $n=15$ it
occasionally even comes up with a better solution than the MILP finds in an
hour. In addition to that, the local search strongly improves the initial
solution for all instance sizes.

##  9 Conclusion and future work

Our hybrid local search approach matches the MILP formulation in solution
quality for small instances, and is able to find a feasible solution for
larger instances in reasonable time. We apply a decomposition of the problem
in two parts, where a local search algorithm is used to find event orders and
an LP is used to find an optimal schedule for a given event order. This
approach opens the door for finding solutions to larger instances. With
further tuning and refinement of the current approach it is promising for
finding solutions to instances of larger sizes.

One of the avenues for future work would be exactly this, further tuning the
approach to deal with larger instances. The addition of more preprocessing
steps (e.g. identifying more implicit precedences) and the experimentation
with different search strategies (e.g. adding some form of random restart) are
promising directions for further investigation. The MILP formulation can be
further strengthened by adding valid inequalities.

Finally, one could look at further extensions of the problem. Among those of
interest are the addition of precedence relations and changing $P$ from a
constant into a variable resource availability function $P(t)$.

Table 6: Number of feasible solutions found for each instance size | Flow | SA | MIP
---|---|---|---
$n$ | Feas. | Feas. | Best | Feas. | Best
5 | 28 | 28 | 28 | 28 | 28
10 | 32 | 30 | 26 | 32 | 23
15 | 32 | 30 | 15 | 27 | 4
20 | 31 | 30 | 8 | - | -
30 | 32 | 28 | 32 | - | -
50 | 32 | 21 | 9 | - | -
|  |  |  |  |
Table 7: Computational results for $n=5$ | Flow | Best | MILP | SA
---|---|---|---|---
$n$ | $P$ | adv.? | # | Feas.? | known | time | obj | time | obj | init
5 | 25 | 0 | 0 | yes | 163.58 | 2.04 | 163.58 | 24.18 | 163.58 | 177.77
5 | 25 | 0 | 1 | yes | 165.72 | 1.58 | 165.72 | 189.19 | 165.72 | 212.33
5 | 25 | 0 | 2 | yes | 99.42 | 18.87 | 99.42 | 59.68 | 99.42 | 104.79
5 | 25 | 0 | 3 | yes | 78.70 | 1.92 | 78.70 | 80.97 | 78.70 | 97.96
5 | 25 | 1 | 0 | yes | 113.21 | 1.96 | 113.21 | 104.72 | 113.21 | 132.90
5 | 25 | 1 | 1 | yes | 61.94 | 6.30 | 61.94 | 91.03 | 61.94 | 68.71
5 | 25 | 1 | 2 | yes | 73.06 | 1.83 | 73.06 | 111.09 | 73.06 | 88.02
5 | 25 | 1 | 3 | yes | 96.81 | 6.36 | 96.81 | 95.73 | 96.81 | 98.34
5 | 50 | 0 | 0 | yes | 72.39 | 1.56 | 72.39 | 75.05 | 72.39 | 76.19
5 | 50 | 0 | 1 | yes | 88.61 | 0.87 | 88.61 | 103.08 | 88.61 | 94.87
5 | 50 | 0 | 2 | yes | 95.38 | 1.13 | 95.38 | 100.95 | 95.38 | 101.10
5 | 50 | 0 | 3 | yes | 80.81 | 1.38 | 80.81 | 130.12 | 80.81 | 92.27
5 | 50 | 1 | 0 | yes | 98.25 | 4.23 | 98.25 | 121.34 | 98.25 | 112.77
5 | 50 | 1 | 1 | yes | 74.93 | 2.24 | 74.93 | 101.55 | 74.93 | 82.42
5 | 50 | 1 | 2 | yes | 83.88 | 0.94 | 83.88 | 92.08 | 83.88 | 99.66
5 | 50 | 1 | 3 | yes | 102.10 | 2.08 | 102.10 | 91.60 | 102.10 | 102.95
5 | 100 | 0 | 0 | yes | 75.25 | 0.66 | 75.25 | 91.27 | 75.25 | 91.20
5 | 100 | 0 | 1 | yes | 77.71 | 0.51 | 77.71 | 36.22 | 77.71 | 175.64
5 | 100 | 0 | 2 | yes | 49.32 | 0.95 | 49.32 | 62.31 | 49.32 | 56.07
5 | 100 | 0 | 3 | yes | 51.97 | 0.54 | 51.97 | 64.66 | 51.97 | 111.72
5 | 100 | 1 | 0 | yes | 53.80 | 0.45 | 53.80 | 123.77 | 53.80 | 54.48
5 | 100 | 1 | 1 | yes | 69.92 | 0.71 | 69.92 | 43.58 | 69.92 | 76.33
5 | 100 | 1 | 2 | yes | 93.13 | 0.92 | 93.13 | 42.37 | 93.13 | 100.24
5 | 100 | 1 | 3 | yes | 53.79 | 0.64 | 53.79 | 39.44 | 53.79 | 55.46
5 | 200 | 0 | 0 | yes | 67.13 | 0.38 | 67.13 | 89.15 | 67.13 | 137.92
5 | 200 | 0 | 1 | no | 120.94 | 0.43 | -1.00 | 0.24 | 169.91 | 169.91
5 | 200 | 0 | 2 | no | 74.41 | 0.45 | -1.00 | 1.78 | 74.44 | 146.13
5 | 200 | 0 | 3 | yes | 57.02 | 0.48 | 57.02 | 86.05 | 57.02 | 59.86
5 | 200 | 1 | 0 | yes | 56.35 | 0.39 | 56.35 | 67.70 | 56.35 | 57.73
5 | 200 | 1 | 1 | no | 81.35 | 0.45 | -1.00 | 1.62 | 81.61 | 148.08
5 | 200 | 1 | 2 | yes | 67.19 | 0.61 | 67.19 | 78.81 | 67.19 | 122.79
5 | 200 | 1 | 3 | no | 229.58 | 0.41 | -1.00 | 1.16 | 231.03 | 330.74
|  |  |  |  |  |  |  |  |  |
Table 8: Computational results for $n=10$ | Flow | Best | MILP | SA
---|---|---|---|---
$n$ | $P$ | adv.? | # | Feas.? | known | time | obj | time | obj | init
10 | 25 | 0 | 0 | yes | 359.47 | LIMIT | 360.49 | 193.25 | 359.47 | 424.51
10 | 25 | 0 | 1 | yes | 334.63 | LIMIT | 337.72 | 206.95 | 334.63 | 424.35
10 | 25 | 0 | 2 | yes | 225.99 | LIMIT | 226.76 | 204.47 | 232.05 | 250.28
10 | 25 | 0 | 3 | yes | 300.08 | LIMIT | 300.08 | 243.97 | 300.08 | 392.84
10 | 25 | 1 | 0 | yes | 345.71 | 525.62 | 345.71 | 277.95 | 345.71 | 401.68
10 | 25 | 1 | 1 | yes | 394.69 | 3228.46 | 394.69 | 218.89 | 394.69 | 446.65
10 | 25 | 1 | 2 | yes | 399.06 | LIMIT | 399.06 | 245.50 | 399.06 | 493.58
10 | 25 | 1 | 3 | yes | 378.15 | LIMIT | 378.15 | 280.72 | 378.15 | 400.00
10 | 50 | 0 | 0 | yes | 193.57 | LIMIT | 193.57 | 293.98 | 193.57 | 207.17
10 | 50 | 0 | 1 | yes | 220.46 | LIMIT | 221.18 | 283.44 | 220.87 | 239.74
10 | 50 | 0 | 2 | yes | 254.38 | 1131.58 | 254.38 | 332.78 | 254.38 | 289.27
10 | 50 | 0 | 3 | yes | 191.05 | LIMIT | 191.05 | 248.19 | 191.05 | 213.67
10 | 50 | 1 | 0 | yes | 162.90 | LIMIT | 162.90 | 281.60 | 162.90 | 179.64
10 | 50 | 1 | 1 | yes | 171.90 | LIMIT | 172.04 | 298.19 | 171.90 | 174.68
10 | 50 | 1 | 2 | yes | 194.43 | LIMIT | 194.56 | 1.17 | 552.53 | 595.79
10 | 50 | 1 | 3 | yes | 289.80 | LIMIT | 289.80 | 256.29 | 303.90 | 388.97
10 | 100 | 0 | 0 | yes | 155.15 | 29.48 | 155.15 | 312.78 | 155.15 | 162.38
10 | 100 | 0 | 1 | yes | 139.62 | LIMIT | 139.62 | 308.31 | 139.62 | 148.31
10 | 100 | 0 | 2 | yes | 203.37 | LIMIT | 203.37 | 284.11 | 203.37 | 211.90
10 | 100 | 0 | 3 | yes | 168.99 | LIMIT | 171.55 | 306.39 | 168.99 | 194.37
10 | 100 | 1 | 0 | yes | 205.44 | LIMIT | 205.44 | 348.42 | 205.44 | 219.86
10 | 100 | 1 | 1 | yes | 171.26 | LIMIT | 171.42 | 393.40 | 171.26 | 225.93
10 | 100 | 1 | 2 | yes | 143.83 | LIMIT | 143.83 | 253.98 | 143.83 | 163.53
10 | 100 | 1 | 3 | yes | 145.66 | LIMIT | 148.32 | 282.35 | 148.32 | 162.14
10 | 200 | 0 | 0 | yes | 125.38 | 181.10 | 125.38 | 221.56 | 125.38 | 130.17
10 | 200 | 0 | 1 | yes | 138.32 | 4.58 | 138.32 | 309.86 | 138.32 | 185.09
10 | 200 | 0 | 2 | yes | 107.74 | 33.39 | 107.74 | 276.53 | 107.74 | 113.15
10 | 200 | 0 | 3 | yes | 100.16 | 4.42 | 100.16 | 314.35 | 100.16 | 104.39
10 | 200 | 1 | 0 | yes | 158.68 | 18.55 | 158.68 | 307.19 | 158.68 | 182.99
10 | 200 | 1 | 1 | yes | 135.62 | 13.38 | 135.62 | 443.79 | 135.62 | 136.47
10 | 200 | 1 | 2 | yes | 108.72 | 7.00 | 108.72 | 193.36 | 108.72 | 109.42
10 | 200 | 1 | 3 | yes | 101.19 | LIMIT | 101.19 | 233.21 | 101.19 | 106.54
|  |  |  |  |  |  |  |  |  |
Table 9: Computational results for $n=15$ | Flow | Best | MILP | SA
---|---|---|---|---
$n$ | $P$ | adv.? | # | Feas.? | known | time | obj | time | obj | init
15 | 25 | 0 | 0 | yes | 831.86 | LIMIT | 897.55 | 538.08 | 831.86 | 871.10
15 | 25 | 0 | 1 | yes | 587.90 | LIMIT | -1.00 | 390.93 | 587.96 | 752.17
15 | 25 | 0 | 2 | yes | 779.20 | LIMIT | 790.29 | 399.92 | 779.20 | 985.93
15 | 25 | 0 | 3 | yes | 800.74 | LIMIT | 825.81 | 417.82 | 800.74 | 997.23
15 | 25 | 1 | 0 | yes | 733.29 | LIMIT | 770.22 | 326.79 | 733.29 | 882.03
15 | 25 | 1 | 1 | yes | 805.43 | LIMIT | 898.95 | 536.96 | 811.87 | 953.44
15 | 25 | 1 | 2 | yes | 649.83 | LIMIT | -1.00 | 563.08 | 649.83 | 847.58
15 | 25 | 1 | 3 | yes | 561.05 | LIMIT | 572.22 | 402.05 | 561.98 | 662.60
15 | 50 | 0 | 0 | yes | 542.53 | LIMIT | -1.00 | 448.71 | 542.57 | 634.80
15 | 50 | 0 | 1 | yes | 309.19 | LIMIT | 311.46 | 544.40 | 309.19 | 363.38
15 | 50 | 0 | 2 | yes | 353.07 | LIMIT | 358.26 | 411.76 | 353.08 | 382.82
15 | 50 | 0 | 3 | yes | 391.34 | LIMIT | -1.00 | 462.41 | 391.48 | 611.09
15 | 50 | 1 | 0 | yes | 322.26 | LIMIT | 325.22 | 520.40 | 322.26 | 349.78
15 | 50 | 1 | 1 | yes | 475.77 | LIMIT | 492.32 | 2.18 | 793.14 | 793.14
15 | 50 | 1 | 2 | yes | 589.86 | LIMIT | -1.00 | 436.84 | 589.86 | 626.66
15 | 50 | 1 | 3 | yes | 347.26 | LIMIT | 349.15 | 392.01 | 347.26 | 399.45
15 | 100 | 0 | 0 | yes | 235.75 | LIMIT | 246.03 | 427.83 | 239.26 | 250.06
15 | 100 | 0 | 1 | yes | 330.64 | LIMIT | 337.54 | 507.91 | 331.21 | 338.12
15 | 100 | 0 | 2 | yes | 314.33 | LIMIT | 321.52 | 643.93 | 314.33 | 348.71
15 | 100 | 0 | 3 | yes | 288.97 | LIMIT | 289.10 | 479.21 | 288.97 | 342.44
15 | 100 | 1 | 0 | yes | 279.48 | LIMIT | 291.05 | 542.24 | 279.58 | 299.08
15 | 100 | 1 | 1 | yes | 341.46 | LIMIT | 345.56 | 671.86 | 341.46 | 408.68
15 | 100 | 1 | 2 | yes | 215.79 | LIMIT | 230.71 | 2.30 | 520.72 | 520.72
15 | 100 | 1 | 3 | yes | 259.86 | LIMIT | 265.08 | 546.11 | 259.86 | 296.06
15 | 200 | 0 | 0 | yes | 168.23 | LIMIT | 190.76 | 579.06 | 168.24 | 181.16
15 | 200 | 0 | 1 | yes | 253.88 | LIMIT | 261.80 | 881.63 | 253.88 | 259.09
15 | 200 | 0 | 2 | yes | 143.38 | LIMIT | 144.71 | 702.94 | 143.51 | 151.02
15 | 200 | 0 | 3 | yes | 237.71 | 1496.19 | 237.71 | 549.86 | 237.71 | 243.81
15 | 200 | 1 | 0 | yes | 179.79 | LIMIT | 179.86 | 646.97 | 179.85 | 199.15
15 | 200 | 1 | 1 | yes | 198.08 | LIMIT | 198.08 | 526.01 | 198.10 | 207.30
15 | 200 | 1 | 2 | yes | 160.21 | LIMIT | 160.21 | 565.79 | 160.34 | 165.41
15 | 200 | 1 | 3 | yes | 244.38 | LIMIT | 244.38 | 556.86 | 244.41 | 254.76
|  |  |  |  |  |  |  |  |  |
Table 10: Computational results for $n=20$ | Flow | Best | SA
---|---|---|---
$n$ | $P$ | adv.? | # | Feas.? | known | time | obj | init
20 | 25 | 0 | 0 | yes | 1122.08 | 862.19 | 1125.55 | 1643.66
20 | 25 | 0 | 1 | yes | 1343.89 | 685.22 | 1343.89 | 1771.98
20 | 25 | 0 | 2 | yes | 1012.18 | 877.71 | 1012.18 | 1197.75
20 | 25 | 0 | 3 | no | 1445.11 | 732.78 | 1518.36 | 1784.42
20 | 25 | 1 | 0 | yes | 991.01 | 776.51 | 995.64 | 1195.65
20 | 25 | 1 | 1 | yes | 1316.27 | 797.72 | 1316.27 | 1768.37
20 | 25 | 1 | 2 | yes | 1094.52 | 549.41 | 1094.52 | 1399.41
20 | 25 | 1 | 3 | yes | 991.20 | 702.07 | 991.20 | 1123.76
20 | 50 | 0 | 0 | yes | 657.76 | 752.77 | 657.76 | 827.27
20 | 50 | 0 | 1 | yes | 765.34 | 712.16 | 766.59 | 951.32
20 | 50 | 0 | 2 | yes | 590.19 | 809.64 | 590.19 | 620.57
20 | 50 | 0 | 3 | yes | 729.51 | 710.34 | 729.51 | 1049.81
20 | 50 | 1 | 0 | yes | 571.33 | 678.59 | 571.48 | 635.55
20 | 50 | 1 | 1 | yes | 625.52 | 674.42 | 628.49 | 755.16
20 | 50 | 1 | 2 | yes | 548.63 | 565.51 | 548.74 | 737.98
20 | 50 | 1 | 3 | yes | 551.01 | 603.31 | 554.73 | 610.62
20 | 100 | 0 | 0 | yes | 301.95 | 989.00 | 302.12 | 314.61
20 | 100 | 0 | 1 | yes | 476.50 | 810.92 | 476.71 | 527.38
20 | 100 | 0 | 2 | yes | 363.01 | 902.16 | 363.12 | 443.13
20 | 100 | 0 | 3 | yes | 488.63 | 1461.01 | 488.68 | 506.45
20 | 100 | 1 | 0 | yes | 365.24 | 980.08 | 365.31 | 501.94
20 | 100 | 1 | 1 | yes | 493.23 | 1072.18 | 493.28 | 584.07
20 | 100 | 1 | 2 | yes | 415.11 | 1185.18 | 415.24 | 454.11
20 | 100 | 1 | 3 | yes | 359.19 | 664.77 | 359.25 | 482.25
20 | 200 | 0 | 0 | yes | 303.41 | 818.06 | 303.78 | 309.97
20 | 200 | 0 | 1 | yes | 310.99 | 883.86 | 311.31 | 364.42
20 | 200 | 0 | 2 | yes | 345.24 | 1064.08 | 345.28 | 367.84
20 | 200 | 0 | 3 | yes | 285.33 | 796.83 | 285.70 | 290.13
20 | 200 | 1 | 0 | yes | 305.00 | 1095.02 | 305.06 | 334.49
20 | 200 | 1 | 1 | yes | 337.89 | 753.72 | 337.97 | 346.15
20 | 200 | 1 | 2 | yes | 313.49 | 828.79 | 313.72 | 332.67
20 | 200 | 1 | 3 | yes | 318.78 | 1098.68 | 318.96 | 365.87
|  |  |  |  |  |  |  |
Table 11: Computational results for $n=30$ | Flow | Best | SA
---|---|---|---
$n$ | $P$ | adv.? | # | Feas.? | known | time | obj | init
30 | 25 | 0 | 0 | yes | 2088.94 | 1502.01 | 2088.94 | 2599.09
30 | 25 | 0 | 1 | yes | 3134.71 | 5.58 | 3134.71 | 3211.39
30 | 25 | 0 | 2 | yes | 2291.27 | 1406.69 | 2291.27 | 2744.97
30 | 25 | 0 | 3 | yes | 2082.82 | 1424.99 | 2082.82 | 2860.31
30 | 25 | 1 | 0 | yes | 2920.08 | 1142.54 | 2920.08 | 3380.82
30 | 25 | 1 | 1 | yes | 2963.77 | 1483.27 | 2963.77 | 3870.67
30 | 25 | 1 | 2 | yes | 3436.78 | 29.96 | 3436.78 | 3622.45
30 | 25 | 1 | 3 | yes | 3237.29 | 1760.39 | 3237.29 | 4677.28
30 | 50 | 0 | 0 | yes | 2426.46 | 5.62 | 2426.46 | 2476.53
30 | 50 | 0 | 1 | yes | 2538.60 | 5.75 | 2538.60 | 2589.17
30 | 50 | 0 | 2 | yes | 1052.26 | 1728.24 | 1052.26 | 1254.98
30 | 50 | 0 | 3 | yes | 2019.35 | 5.69 | 2019.35 | 2020.03
30 | 50 | 1 | 0 | yes | 1660.95 | 1601.76 | 1660.95 | 2041.24
30 | 50 | 1 | 1 | yes | 1732.74 | 1614.14 | 1732.74 | 2391.34
30 | 50 | 1 | 2 | yes | 1778.19 | 1381.26 | 1778.19 | 2190.11
30 | 50 | 1 | 3 | yes | 1703.09 | 308.90 | 1703.09 | 2080.97
30 | 100 | 0 | 0 | yes | 747.01 | 1591.51 | 747.01 | 1482.49
30 | 100 | 0 | 1 | yes | 904.57 | 1628.79 | 904.57 | 1262.40
30 | 100 | 0 | 2 | yes | 779.62 | 2140.44 | 779.62 | 929.26
30 | 100 | 0 | 3 | yes | 865.12 | 1481.99 | 865.12 | 1160.84
30 | 100 | 1 | 0 | yes | 836.81 | 2110.24 | 836.81 | 1118.03
30 | 100 | 1 | 1 | yes | 976.40 | 1890.41 | 976.40 | 1419.84
30 | 100 | 1 | 2 | yes | 754.56 | 1861.07 | 754.56 | 1085.14
30 | 100 | 1 | 3 | yes | 695.52 | 1690.72 | 695.52 | 902.68
30 | 200 | 0 | 0 | yes | 497.55 | 2059.61 | 497.55 | 549.87
30 | 200 | 0 | 1 | yes | 486.30 | 1826.58 | 486.30 | 513.94
30 | 200 | 0 | 2 | yes | 509.32 | 1685.78 | 509.32 | 552.84
30 | 200 | 0 | 3 | yes | 538.68 | 2378.16 | 538.68 | 604.61
30 | 200 | 1 | 0 | yes | 468.25 | 2631.23 | 468.25 | 623.84
30 | 200 | 1 | 1 | yes | 586.33 | 1823.30 | 586.33 | 715.84
30 | 200 | 1 | 2 | yes | 619.83 | 1826.30 | 619.83 | 645.95
30 | 200 | 1 | 3 | yes | 568.19 | 2356.85 | 568.19 | 687.19
|  |  |  |  |  |  |  |
Table 12: Computational results for $n=50$ | Flow | Best | SA
---|---|---|---
$n$ | $P$ | adv.? | # | Feas.? | known | time | obj | init
50 | 25 | 0 | 0 | yes | 6002.37 | 4293.20 | 6777.69 | 7958.75
50 | 25 | 0 | 1 | yes | 6573.68 | 4665.56 | 7245.46 | 7509.22
50 | 25 | 0 | 2 | yes | 6423.54 | 5003.16 | 6423.54 | 7689.15
50 | 25 | 0 | 3 | yes | 7027.25 | 5916.94 | 7027.25 | 8145.16
50 | 25 | 1 | 0 | yes | 6562.65 | 3602.98 | 7196.96 | 7217.22
50 | 25 | 1 | 1 | yes | 6257.97 | 5594.27 | 7278.21 | 7520.76
50 | 25 | 1 | 2 | yes | 7239.58 | 6254.32 | 7239.58 | 8007.63
50 | 25 | 1 | 3 | yes | 7911.54 | 5082.28 | 7925.13 | 8873.06
50 | 50 | 0 | 0 | yes | 4176.25 | 5303.19 | 4469.83 | 5420.81
50 | 50 | 0 | 1 | yes | 5837.52 | 3731.27 | 5857.73 | 5871.31
50 | 50 | 0 | 2 | yes | 3167.27 | 2318.80 | 3756.45 | 4459.79
50 | 50 | 0 | 3 | yes | 3429.76 | 3484.34 | 4080.28 | 4731.48
50 | 50 | 1 | 0 | yes | 3503.36 | 6897.78 | 3516.58 | 4041.55
50 | 50 | 1 | 1 | yes | 3081.82 | 6070.92 | 3081.82 | 4950.52
50 | 50 | 1 | 2 | yes | 4326.41 | 2764.69 | 4499.39 | 4509.67
50 | 50 | 1 | 3 | yes | 3548.99 | 3541.65 | 4359.87 | 4750.75
50 | 100 | 0 | 0 | yes | 1609.52 | 8184.16 | 1609.52 | 1928.63
50 | 100 | 0 | 1 | yes | 2205.73 | 8174.82 | 2205.73 | 2588.95
50 | 100 | 0 | 2 | yes | 1658.84 | 3838.16 | 2074.60 | 2127.92
50 | 100 | 0 | 3 | yes | 1947.59 | 2484.76 | 2282.36 | 2622.42
50 | 100 | 1 | 0 | yes | 1628.10 | 4110.64 | 1933.85 | 2240.80
50 | 100 | 1 | 1 | yes | 1716.27 | 9738.11 | 1716.27 | 2096.80
50 | 100 | 1 | 2 | yes | 1637.18 | 6344.71 | 1637.18 | 2402.85
50 | 100 | 1 | 3 | yes | 1880.46 | 6946.31 | 1880.46 | 2358.75
50 | 200 | 0 | 0 | yes | 1139.72 | 2123.86 | 1384.30 | 1415.43
50 | 200 | 0 | 1 | yes | 1237.01 | 2861.69 | 1352.59 | 1352.59
50 | 200 | 0 | 2 | yes | 1298.67 | 12194.00 | 1306.81 | 1439.92
50 | 200 | 0 | 3 | yes | 1211.63 | 9607.54 | 1219.67 | 1442.92
50 | 200 | 1 | 0 | yes | 1052.00 | 13363.99 | 1060.47 | 1169.24
50 | 200 | 1 | 1 | yes | 1265.67 | 2003.09 | 1458.85 | 1459.74
50 | 200 | 1 | 2 | yes | 1060.45 | 9351.82 | 1070.15 | 1223.95
50 | 200 | 1 | 3 | yes | 1095.45 | 8013.34 | 1103.34 | 1216.94
|  |  |  |  |  |  |  |

## References

  *   * Baptiste . (1999) Baptiste1999Baptiste, P., Pape, C.L. Nuijten, W.  1999\.  Satisfiability tests and time-bound adjustments for cumulative scheduling problems Satisfiability tests and time-bound adjustments for cumulative scheduling problems.  Annals of Operations Research92305–333.  https://doi.org/10.1023/A:1018995000688
  * Hartmann Briskorn (2010) Hartmann2010Hartmann, S. Briskorn, D.  2010\.  A survey of variants and extensions of the resource-constrained project scheduling problem A survey of variants and extensions of the resource-constrained project scheduling problem.  European Journal of Operational Research2071–14.  https://doi.org/10.1016/j.ejor.2009.11.005
  * Hartmann Briskorn (2022) Hartmann2022Hartmann, S. Briskorn, D.  2022\.  An updated survey of variants and extensions of the resource-constrained project scheduling problem An updated survey of variants and extensions of the resource-constrained project scheduling problem.  European Journal of Operational Research2971–14.  https://doi.org/10.1016/j.ejor.2021.05.004
  * Nattaf . (2014) Nattaf2014Nattaf, M., Artigues, C. Lopez, P.  2014\.  A polynomial satisfiability test using energetic reasoning for energy-constraint scheduling A polynomial satisfiability test using energetic reasoning for energy-constraint scheduling.  International Conference on Project Management and Scheduling (PMS 2014) International conference on project management and scheduling (pms 2014) ( 169–172).  https://hal.archives-ouvertes.fr/hal-00978706
  * Nattaf . (2017) Nattaf2017Nattaf, M., Artigues, C. Lopez, P.  201710\.  Cumulative scheduling with variable task profiles and concave piecewise linear processing rate functions Cumulative scheduling with variable task profiles and concave piecewise linear processing rate functions.  Constraints22530–547.  https://doi.org/10.1007/S10601-017-9271-4
  * Nattaf . (2016) Nattaf2016Nattaf, M., Artigues, C., Lopez, P. Rivreau, D.  20163\.  Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions.  OR Spectrum38459–492.  https://doi.org/10.1007/S00291-015-0423-X
  * Nattaf . (2019) Nattaf2019Nattaf, M., Horváth, M., Kis, T., Artigues, C. Lopez, P.  2019\.  Polyhedral results and valid inequalities for the continuous energy-constrained scheduling problem Polyhedral results and valid inequalities for the continuous energy-constrained scheduling problem.  Discrete Applied Mathematics258188–203.  https://doi.org/10.1016/j.dam.2018.11.008
  * Turek . (1992) Turek1992Turek, J., Wolf, J.L. Yu, P.S.  1992\.  Approximate algorithms for scheduling parallelizable tasks Approximate algorithms for scheduling parallelizable tasks.  Proceedings of the 4th Annual ACM Symposium on Parallel Algorithms and Architectures Proceedings of the 4th annual acm symposium on parallel algorithms and architectures ( 323–332).  https://doi.org/10.1145/140901.141909