3.1 Data preprocessing and distribution analysis

Data preprocessing is a crucial step in ML, particularly when working with real-world medical datasets that often suffer from skewed distributions, noise, and sparsity. In our study, the target variable, DURATA, representing the total operative time in minutes, was derived from the difference between USCITASALA (exit from OR) and INGRESSOSALA (entry). Both fields were parsed as datetime objects, and duration was computed as the difference in minutes. Entries with negative or zero durations, likely due to data entry errors, were removed from the dataset to ensure data integrity and modeling accuracy. As the goal of our study is regression, addressing target distribution characteristics is particularly important. As shown in Figure 1, the original distribution of the intervention durations was highly skewed, with a large concentration of cases around 15 minutes and a long tail of less frequent, prolonged procedures. This distribution skewness can significantly affect model performance, leading to biased predictions and poor generalization, particularly for minority classes or rare cases. To address these issues, we employed the following preprocessing steps:

1. 1. Diagnoses that appeared only once within each department (REPARTO) were identified and grouped by K-Means clustering (with up to 3 clusters). This step reduced data sparsity and helped generalize rare cases by associating them with similar groups according to the department.

2. 2. We filtered out extreme values that lay outside the typical range of durations using the Interquartile Range (IQR) method, that is a statistical measure of dispersion, calculated as the difference between the third quartile (Q3) and the first quartile (Q1) of the data (i.e., IQR = Q3 – Q1) (Dekking et al. Reference Dekking, Kraaikamp, Lopuhaä and Meester2005). This improved the central tendency and distribution spread, enhancing the learning stability of regression models.

3. 3. To mitigate issues related to multicollinearity, we removed features that exhibited high pairwise correlation (above a 0.95 threshold). Redundant features can introduce noise, inflate model complexity, and impair generalization. By identifying and eliminating highly correlated variables, we reduced dimensionality and improved model robustness and interpretability. Indeed, by removing highly correlated features, we aimed to reduce redundancy in the dataset. When two or more features are strongly correlated, they tend to capture the same underlying information, making it difficult to isolate their individual contributions to the model. This redundancy complicates interpretation especially when assessing which variables are truly driving predictions, and reduces model stability. As discussed in Kuhn et al. (Reference Kuhn, Johnson, Kuhn and Johnson2013), eliminating one of several highly correlated features typically does not impair predictive performance and can lead to a simpler, more interpretable model.

In particular, Figure 2 illustrates the duration distribution after Step 2 of preprocessing. Compared to the raw dataset, the cleaned data shows a more compact distribution with fewer extreme outliers and a less pronounced right skew. Furthermore, the preprocessing performed in Step 3 resulted in a reduction of feature dimensionality from 32 to 23, reflecting the elimination of less informative features.

Fig. 1. Histogram of intervention durations (before preprocessing).

Fig. 2. Histogram of intervention durations (after preprocessing).

3.2 Predictive modeling

To predict surgical procedure durations, several ML regressors were implemented and compared. All models were trained and evaluated under identical preprocessing and evaluation pipelines to ensure fairness. The dataset was split into training and testing subsets (80 % and 20 %, respectively) using stratified sampling based on procedure duration. This strategy ensured that the distribution of durations, particularly for rare or long procedures, was preserved across both sets. In this way, we avoided scenarios where specific duration ranges (e.g., particularly long procedures) appeared only in the test set, which could otherwise result in biased or unreliable performance estimations. The following algorithms were explored:

• • Decision Tree Regressor (DT): a simple, interpretable model that recursively splits the data based on feature thresholds to minimize prediction error (Breiman et al. Reference Breiman, Friedman, Olshen and Stone1984).

• • Random Forest Regressor (RF): an ensemble method that constructs multiple decision trees using bootstrap sampling and random feature selection at each split (Breiman Reference Breiman2001). Final predictions are computed as the average of individual tree predictions.

• • Gradient Boosting Regressor (GB): a sequential ensemble technique where each new tree is trained to correct the residuals of the previous ensemble (Friedman Reference Friedman2001). It combines weak learners into a strong predictor and includes hyperparameters for controlling learning rate, tree depth, and regularization.

• • Extreme Gradient Boosting (XGBoost): a highly optimized and regularized gradient boosting framework (Chen and Guestrin Reference Chen and Guestrin2016), known for its scalability and efficiency. XGBoost introduces system-level optimizations (e.g., parallelization) and algorithmic enhancements like shrinkage and sparsity-aware split finding.

• • K-Nearest Neighbors Regressor (KNN): a non-parametric method that predicts the target as the average of the $k$ nearest neighbors in the feature space (Altman Reference Altman1992).

• • Support Vector Regressor (SVR): a margin-based regression technique that aims to find a function within an $\epsilon$ -tube from the true outputs, penalizing predictions only when errors exceed $\epsilon$ (Drucker et al. Reference Drucker, Burges, Kaufman, Smola and Vapnik1996).

3.2.1 Deep learning (DL) models

DL models are generally not well-suited for tabular data, as highlighted by Shwartz-Ziv and Armon (Reference Shwartz-Ziv and Armon2022) and by Borisov et al. (Reference Borisov, Leemann, Seßler, Haug, Pawelczyk and Kasneci2024). Nevertheless, some of the recent DL approaches obtained good performance on tabular data, as the one proposed by Arik and Pfister (Reference Arik and Pfister2021), called TabNet. Indeed, TabNet is a deep neural architecture tailored for tabular data that uses sequential attention to select features at each decision step, enabling both high performance and interpretability. Motivated by this promising finding, we tested the performance of TabNet on our dataset and we adapted and evaluated the classical DL models Multi-Layer Perceptron (MLP) (Almeida Reference Almeida2020) and 1D Convolutional Neural Network (1D-CNN) (Ige and Sibiya Reference Ige and Sibiya2024). Each model was optimized via grid search and the best configuration is reported in Table 2. In more detail, MLP was implemented as a feedforward neural network with the Rectified Linear Unit (ReLU), a nonlinear activation function that introduces nonlinearity into the model and helps mitigate the vanishing gradient problem. It outputs the input directly if it is positive and zero otherwise, effectively retaining only the positive part of its argument. We also used mean squared error loss and adaptive moment estimation (Adam) optimizer to make training faster and more stable, especially on noisy or sparse data (Kingma and Ba Reference Kingma and Ba2015). Our implementation followed a standard configuration with two hidden layers and 64 units per layer. 1D-CNN was adapted to tabular input by reshaping the feature vectors as sequences, applying a single convolutional layer with 32 filters and kernel size 2, followed by a fully connected hidden layer and regression output. TabNet was tested using a flexible wrapper that allowed us to tune key architectural parameters such as the size of internal layers, the number of steps in the attention process, and the strength of feature sparsity. Unlike standard neural networks, TabNet learns which features to use at each step, making it more efficient for tabular tasks.

Table 2. Best hyperparameter configuration for each algorithm

Table 3. Model performance using best parameters. Best results are in bold

3.3 Evaluation metrics and experiments

To evaluate the performance of the ML and DL models in predicting surgical procedure durations, we employed three standard regression metrics. The first one is the mean absolute error (MAE) that reflects the average magnitude of prediction errors, and is computed as the average of the absolute differences between the predicted values $\hat {y}_i$ and the true values $y_i$ , that is $\text{MAE} = \frac {1}{n} \sum _{i=1}^{n} |y_i - \hat {y}_i|$ . The second metric is the root mean squared Error (RMSE), which penalizes larger errors more than MAE and reflects the model’s sensitivity to outliers. RMSE is calculated as the square root of the mean squared error, that is $\text{RMSE} = \sqrt { \frac {1}{n} \sum _{i=1}^{n} (y_i - \hat {y}_i)^2 }$ . The last one is the coefficient of determination (R²), which measures the proportion of variance in the target variable explained by the model, that is $\textrm {R}^{2}$ $ = 1 - \frac { \sum _{i=1}^{n} (y_i - \hat {y}_i)^2 }{ \sum _{i=1}^{n} (y_i - \bar {y})^2 }$ , where $\bar {y}$ is the mean of the observed values. In addition to model evaluation, we conducted hyperparameter tuning using grid search with cross-validation. The final configuration for each model was selected based on the combination of hyperparameters that achieved the lowest MAE and RMSE averaged across the validation folds during cross-validation, as shown in Table 2 (for details about the hyperparameters see https://github.com/DeMaCS-UNICAL/ML4ORS). Moreover, Table 3 highlights the performance of each ML model in predicting surgical procedure durations. Among all candidates, the XGBoost Regressor achieved the best results, with the lowest MAE (12.36 minutes), lowest RMSE (18.20 minutes), and highest $\textrm {R}^{2}$ score (0.79). RF and GB followed closely behind, also demonstrating strong predictive performance with slightly higher error metrics. In contrast, simpler models such as DT, KNN, and SVR showed substantially lower performance. The DT and KNN regressors yielded higher MAE values (16.99 and 15.69 minutes, respectively) and lower $\textrm {R}^{2}$ scores (0.55 and 0.64), indicating limited generalization. The SVR model performed comparably to the DT, with the highest MAE (17.22 minutes) and RMSE (26.87 minutes), and an $\textrm {R}^{2}$ of 0.55, suggesting both struggled to capture the underlying structure of the data.

These discrepancies in performance can be partially explained by the distributional characteristics of the target variable. As shown in Figure 2, the duration data remains skewed even after preprocessing. Most surgeries are short-clustered around 15–20 minutes, but a long tail of high-duration outliers exists. This distribution skewness makes the prediction task inherently challenging, particularly for models like SVR and KNN, which are more sensitive to variance and do not adapt well to skewed or noisy distributions. In contrast, ensemble methods like XGBoost and GB are well-suited to handle skewed data, thanks to their iterative learning approach and ability to capture complex, nonlinear relationships. It is important to observe that the overall error is still non-negligible (as it is around 12 minutes on average). Nevertheless, this level of accuracy, combined with the confidence estimation framework (where over 65 % of predictions were classified as High or Moderate confidence) makes the models suitable for potential deployment in time-sensitive hospital operations and to be employed to improve the robustness of the schedules, as we will show in Section 4.3.

Concerning DL approaches, it is possible to observe that they do not outperform XGBoost in this task. However, TabNet demonstrates improved performance compared to other DL models, with an MAE of 14.45 minutes, RMSE of 22.38 minutes, and $\textrm {R}^{2}$ of 0.69, which is better than both the 1D-CNN (MAE: 15.24, RMSE: 22.57, $\textrm {R}^{2}$ : 0.68) and MLP (MAE: 15.73, RMSE: 23.28, $\textrm {R}^{2}$ : 0.66). Nevertheless, despite its competitive performance, TabNet still remains slightly behind XGBoost (MAE: 12.36, RMSE: 18.20, $\textrm {R}^{2}$ : 0.79).

Fig. 3. SHAP summary plot for the best-performing regression model (XGBoost).

Moreover, to interpret the contribution of each feature to the model’s predictions, we employed SHAP (SHapley Additive exPlanations) that assigns each feature an importance value for a particular prediction (Lundberg and Lee Reference Lundberg and Lee2017). Figure 3 shows the SHAP summary plot for the best-performing model (XGBoost), where features are ranked by their overall contribution to the predictions. The horizontal spread reflects the variability in feature impact across all samples, while the color gradient indicates the feature value (from low in black to high in red). As an example, REGRICOVERO, representing the type of hospital admission, shows a wide range of SHAP values, highlighting its strong and variable influence on predicted operative duration. Higher values of this feature (in red) tend to increase the predicted duration, whereas lower values (in blue) are associated with shorter procedures. This suggests the model has learned to associate certain admission types (e.g., complex or urgent cases) with longer surgery times.

3.4 Confidence estimation

To assess the precision of individual predictions, we computed the absolute percentage error (APE) for each instance. This metric quantifies the relative difference between the predicted value $\hat {y}$ and the true target value $y$ :

\begin{equation*} \text{APE} = \frac {|\hat {y} - y|}{y} \cdot 100. \end{equation*}

In particular, the absolute value ensures that the error is always positive, and the percentage form allows for intuitive interpretation across varying scales.

Based on the APE, we categorized each prediction into one of four confidence levels, namely high confidence when APE is less than 10 %, moderate confidence when APE is between 10 % and 25 %, low confidence when APE is between 25 % and 50 %, and very low confidence when APE is greater than or equal to 50 %.

4.1 ASP-based ORS solution

In this section, we briefly describe an existing ASP-based encoding for the ORS problem. We begin by describing the problem’s input and expected output and then show the core rules that model the scheduling constraints and objectives within the ASP framework. This encoding serves as the baseline for the extensions introduced in the following sections.

4.1.1 Data model

The input data is specified by means of the following constants and atoms. Instances of registration(ID,P,SP,DUR) represent the registration of the patient identified by an ID (ID) with priority level (P), the requested specialty (SP), and the expected duration of the surgery (DUR). Instances of mss(OR,SP,SHIFT,DAY) represent which specialty (SP) is assigned to an OR (OR) in a shift (SHIFT) on a day (DAY). Instances of shift(SHIFT,DURATION) indicate that the total duration of all surgeries scheduled in shift SHIFT must not exceed DURATION (expressed in time slots). The output is represented by atoms of the form x(ID,P,OR,DAY,SHIFT), whose intuitive meaning is that the patient identified by an ID (ID) having a priority (P) is assigned to the OR (OR) in the shift (SHIFT) on the day (DAY).

Listing 1. ASP encoding for the ORS problem.

4.1.2 Encoding

The related encoding is shown in Listing 1. The choice rule assigns an OR, a day, a shift to each registration. The second rule ensures that each registration is assigned at most once. The third rule ensures that the total duration of the assigned registration does not exceed the length of the shift. The fourth rule ensures that every registration with priority $p_1$ is assigned to some OR on a day, shift. The weak constraint then minimizes the number of unassigned registrations with priority $p_2$ , $p_3$ , and $p_4$ . Finally, the last constraint is added to ensure that the specific OR “OR A” is assigned to at most one patient, as it was reserved for emergencies and used in this limited way in the original data.

Listing 2. Additional rules to take into account confidence.

4.2 Integration of the ML predictions in ASP

As discussed before, the ORS solutions described by Dodaro et al. (Reference Dodaro, Galatà, Gebser, Maratea, Marte, Mochi and Scanu2024) can only verify whether the encoding aligns with the actual data and, at most, suggest alternative schedules that could have been computed. Thus, it is not possible to generate provisional schedules. Furthermore, the resulting schedules are not always robust. In this section, we describe how this encoding has been extended to incorporate predictive information obtained through ML techniques. The core idea is to exploit the confidence scores produced by the predictive models (specifically, the XGBoost algorithm, which achieved the best results in our evaluation, as shown in Section 3.3) in order to guide the scheduling decisions made by the ASP program.

To this end, we first update the parameter DUR of the registrations. Indeed, now the duration is predicted by the ML algorithm. Then, we also introduce atoms of the form confidence(ID,L) that collect the confidence associated to each registration, where L is 1 if the confidence is High, 2 if the confidence is Moderate, 3 if the confidence is Low, and 4 if the confidence is Very Low.

Finally, we handle these two changes by adding a new set of rules, reported in Listing 2. The key idea is to prefer balanced and high-confidence scheduling decisions. Indeed, the first rule derives the sum of all the confidences assigned to a particular day. The subsequent two rules derive the values corresponding to the maximum and the minimum sum of confidences in each day, respectively. Then, the first weak constraint penalizes solutions where the maximum sum of confidence scores assigned to a day and OR is high, whereas the second one promotes balance across ORs and days by penalizing the difference between the maximum and minimum confidence sums. In this way, we aim at distributing equally among the ORs and days the patients with high confidence value.

4.3 Experiments

This section presents the experimental evaluation conducted to assess the effectiveness of the proposed neuro-symbolic approach.

First of all, it is important to observe that, in the context of ORS, it is not possible to determine the exact duration of a surgical procedure in advance. Typically, the parameter DUR of the atoms registration(ID,P,SP,DUR) is assigned by hospital staff based on prior experience. In most cases, it corresponds either to the average duration of surgeries of the same type or to the average duration of surgeries performed within the same department. Therefore, for the purposes of the following experimental analysis, we compare these two traditional approaches (average by procedure type and average by department) against ML-based strategies described in the previous sections. Specifically, we replace the DUR value in the ASP encoding with the respective predicted duration, compute the surgical schedule using ASP, and then evaluate the quality of the schedules based on the actual durations of the procedures. The real durations, available only after the scheduling has been executed, are crucial for accurately estimating how effectively the ORs were used on each day.

We evaluate all the approaches on real-world data provided by ASL1 Liguria and already used by Dodaro et al. (Reference Dodaro, Galatà, Gebser, Maratea, Marte, Mochi and Scanu2024). The evaluation focuses on key metrics such as the OR percentage occupancy (mean, standard deviation, minimum, and maximum), and the number of times in which an OR has been underbooked (an OR is used less than 80 % of the time) or overbooked (an OR is used for more than 100 % of the time). It is important to note that, with respect to these parameters, ASP solutions (even within the same solver) may exhibit unintuitive behavior. For instance, a solution that is preferable in terms of optimality (e.g., lower cost) may still yield worse overbooking and underbooking values, since these aspects are evaluated only in a post-processing phase. To mitigate this issue, in our encoding, the weak constraints related to confidence are assigned a lower priority level than the original ones. This ensures that performance differences can be attributed primarily to the introduction of the new constraints.

The comparison has been carried out on an Apple M1 CPU machine with 8 GB of physical RAM and a time limit of 60 s per run. As ASP system, we used clingo v. 5.6.2, configured with parameters --restart-on-model and --parallel-mode = 6. These parameters have been found to be effective in a preliminary analysis we performed with several options.

The results are presented in Table 4. The column VBA refers to the virtually best approach, in which the durations are set to their actual values, known only after the surgeries have been performed. This serves as a reference for the optimal performance achievable by predictive methods. The column Conf. refers to the method in which the parameter DUR is set using the durations predicted by the ML model XGBoost, with the ASP encoding also taking into account the associated confidence information. The column Pred. corresponds to the method where DUR is set using the XGBoost predictions without considering the confidence scores. The columns Dep. and Surg. represent the methods where DUR is assigned based on the average duration per department and per surgical procedure, respectively.

Table 4. Comparison of the different methods

For the evaluation of the scheduling results, it is important to note that the ideal situation is to maximize the usage of ORS without exceeding their available capacity. Thus, an OR occupancy rate close to, but not exceeding, 100 % is preferred. Schedules that lead to overbooking (i.e., planned durations exceeding the available OR time) are undesirable, as they can cause operational disruptions and delays. Similarly, underbooking (i.e., leaving significant unused OR time, i.e. below 80 %) should be minimized, although it is generally less critical than overbooking. Therefore, among the generated schedules, the most desirable are those that maintain high, balanced occupancy rates without causing overbooking and with minimal underbooking.

By analyzing the results, we observe that the method Conf. (XGBoost predictions with confidence information) achieves the best overall performance, followed by Pred. (XGBoost predictions without confidence). In Bordighera, the mean OR occupancy with Conf. is 96 %, which is very close to the ideal 100 % threshold, outperforming Dep. (101 %) and Surg. (88 %). Additionally, the standard deviation is reasonably low, indicating consistent scheduling performance across different days. The number of overbooked cases is slightly lower or comparable to the other methods, and underbooking is almost negligible. Also Pred. obtains a good performance being quite close to Conf.. In Imperia, all methods reach a mean occupancy around 100 %, but Conf. maintains better control over the standard deviation compared to Dep. and Surg., and slightly fewer overbooked cases are observed. Similarly, in Sanremo, the Conf. and Pred. methods result in a mean occupancy closer to 100 % (101 %), whereas methods like Dep., and Surg. tend to overbook, with mean occupancies exceeding 103 %. Furthermore, Conf. achieves the lowest number of overbooking (5) compared to all the other methods, and avoids underbooking altogether.

Overall, integrating confidence information into the scheduling process improves the average OR usage and also helps limiting extreme cases of overbooking and underbooking, leading to more balanced and operationally feasible schedules.