Sustainable and Robust Home Healthcare Logistics: A Response to the COVID-19 Pandemic

Abstract

Today, research on healthcare logistics is an important challenge in developing and developed countries, especially when a pandemic such as COVID-19 occurs. The responses required during such a pandemic would benefit from an efficiently designed model for robust and sustainable healthcare logistics. In this study, we focus on home healthcare logistics and services for planning the routing and scheduling of caregivers to visit patients’ homes. Due to the need for social distancing during the COVID-19 pandemic, these services are highly applicable for reducing the growth of the epidemic. In addition to this challenge, home healthcare logistics and services must be redesigned to meet the standards of a triple bottom line approach based on sustainable development goals. A triple bottom line approach finds a balance between economic, environmental, and social criteria for making a sustainable decision. Although, recently, the concept of green home healthcare has been studied based on the total cost and green emissions of home healthcare logistics and services, as far as we know, no research has been conducted on the formulation of a triple bottom line approach for home healthcare logistics and services. To achieve social justice for caregivers, the goal of balancing working time is to find a balance between unemployment time and overtime. Another contribution of this research is to develop a scenario-based robust optimization approach to address the uncertainty of home healthcare logistics and services and to assist with making robust decisions for home healthcare planning. Since our multi-objective optimization model for sustainable and robust home healthcare logistics and services is more complex than other studies, the last novel contribution of this research is to establish an efficient heuristic algorithm based on the Lagrangian relaxation theory. An initial solution is found by defining three heuristic algorithms. Our heuristic algorithms use a symmetric pattern for allocating patients to pharmacies and planning the routing of caregivers. Then, a combination of the epsilon constraint method and the Lagrangian relaxation theory is proposed to generate high-quality Pareto-based solutions in a reasonable time period. Finally, an extensive analysis is done to show that our multi-objective optimization model and proposed heuristic algorithm are efficient and practical, as well as some sensitivities are studied to provide some managerial insights for achieving sustainable and robust home healthcare services in practice.

1. Introduction

The triple bottom line approach aims to redesign logistics and supply chains for implementation of sustainable development goals based on economic, environmental, and social criteria. Taking these factors into consideration creates a sustainable network design [1,2]. Based on international standards such as ISO (International Organization for Standardization (ISO).) 14000 and ISO 26000, respectively, for environmental sustainability and social responsibility [3,4], a sustainable integrated logistics network ensures that all decision-makers focus on economic, environmental, and social aspects simultaneously [5,6]. Simultaneous consideration of these factors is the formulation of sustainability for economic equity, environmental preservation, and social justice criteria. A long-term process to improve a logistics network relies on individuals’ abilities and their actions to perform recent advancements in global logistics sustainability [4,5]. To cope with this new trend in logistics, in this paper, we deploy sustainable home healthcare logistics and services that meet global logistics sustainability.

Nowadays, another main challenge is to create a suitable response and action to an epidemic virus, i.e., COVID-19, for all healthcare practitioners and international associations. Home healthcare services have a significant role in optimizing the cost of public health and reducing the rate of growth by social distancing [6]. Home healthcare services are performed at patients’ homes [7], starting from a pharmacy, then, visiting a group of patients, and then returning to a laboratory to check patient treatments to update patients’ health records. According to the European Commission reports [8], the population of people older than 60 years will increase to 32% in 2030 and reach around 54% by 2060. Therefore, the demand for home healthcare services is growing, as demonstrated in recent decades. It goes without saying that these services are more efficient and cheaper than the same services in hospitals and retirement homes [9]. The best way to reduce the pandemic growth is to perform the care services in homes instead of at public places. It is difficult to consider the social distancing in hospitals and public health services. Therefore, promotion of these services for all patients is an important goal to save patients and to address the COVID-19 pandemic.

To perform home healthcare services, caregivers face different challenges such as the unavailability of patients based on patients’ time windows. In this regard, travel and service times are clearly uncertain. A caregiver must meet a patient’s time window if they want to keep the patient’s appointment; otherwise, a new appointment for this patient is scheduled. There are many factors that impact the uncertainty of travel and service times as well as the time window for patients. For example, road conditions, traffic, driving skills, and caregiver skills can lead to uncertainty for caregivers’ availability and travel and service times [7,9,10]. To simulate this uncertainty, a scenario-based approach is one alternative to consider all pessimistic, optimistic, and realistic scenarios for these parameters. A scenario-based model can formulate all these uncertainties with probabilistic scenarios in order to make a robust decision based on the expected value and deviation of each scenario. To this end, in this study, we propose a multi-objective robust optimization approach for a sustainable home healthcare logistics network problem.

One significant contribution of this research is to shift from a green home healthcare to a sustainable home healthcare logistics network considering economic equity, environmental preservation, and social justice criteria. These factors are based on a triple bottom line approach for our home healthcare logistics network. Formulation of these factors requires a multi-objective robust optimization model for our proposed home healthcare logistics network. Figure 1 shows the goals of this study aimed at achieving sustainable home healthcare logistics and services. In addition to the total cost and environmental pollution, this study contributes to social justice for home healthcare services with working time balancing constraints to optimize both unemployment time and overtime. These objectives and constraints make our sustainable and robust home healthcare logistics network more complex than other studies and this fact highlights the need for an intelligent solution algorithm for this problem.

Home healthcare logistics and services are an extension to the problem of route optimization [11,12,13,14]. Classically, home healthcare logistics have been shown to be an NP-hard problem in the strong sense [15,16]. This study applies an efficient and strong heuristic algorithm based on the Lagrangian relaxation theory and Pareto-based analysis. Initial solutions are generated by heuristic algorithms. Then, the lower and upper bounds are updated by the Lagrangian relaxation theory. Based on the epsilon constraint method, Prato-based solutions are generated. To confirm the high performance of our Lagrangian relaxation-based heuristic algorithm, a comparison is done by the epsilon constraint method using the exact solver.

In conclusion, this study highlights the following contributions:

• A multi-objective robust optimization model for home healthcare logistics is developed.
• A triple bottom line approach is proposed to model sustainable home healthcare logistics based on (1) total cost as the economic factor, (2) environmental emissions of the home healthcare logistics for the environmental sustainability factor, and (3) caregivers’ unemployment times for balancing the working time of caregivers as the social factor.
• A scenario-based robust optimization approach is devoted to address uncertainty and a response to the COVID-19 pandemic.
• Three heuristic algorithms for decisions are proposed based on symmetric patterns for the allocation of patients and routing of caregivers.
• A Lagrangian relaxation-based heuristic algorithm based on a Pareto analysis is developed for the first time.

This paper is organized as follows: In Section 2, we study the background of this research topic and identify the research gaps to confirm our contributions; in Section 3, we propose our multi-objective robust optimization model for our sustainable home healthcare logistics; in Section 4, we create an innovative and intelligent solution algorithm for the proposed optimization model; in Section 5, we provide a set of simulated test studies for analyzing our optimization model and heuristic algorithm, and following comparisons and sensitivity analyses, we discuss and explain managerial insights from the results; finally, in Section 6, we provide a summary of our findings, as well as limitations and future research recommendations.

2. Literature Review

Since home healthcare services have recently become important as a solution for hospitalization problems, researchers have proposed optimization techniques to help practitioners in this industry [13,14,15,16]. Although there are some review papers on this research area [17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33], in this review, we provide the most important and recent studies which have made a significant contribution to home healthcare services and logistics.

Optimization models for home healthcare logistics have been studied from 1977 [17] to 2021 [18,19,20,21,22,23,24,25,26,27,28]. A special decision support system (SDSS) was proposed by Begur et al. [17]. Later, such techniques were considered by some other studies such as Eveborn et al. [19]. The impact of transportation costs on the economic performance of home healthcare companies is undeniable and important; the earliest studies such as Cheng and Rich [29] in 1998, for the first time, applied a routing optimization with time windows for the application of home healthcare logistics. There is no doubt that routing and scheduling optimization problems are difficult, and many studies have contributed various heuristic and metaheuristic algorithms [30].

One of these popular algorithms was the particle swarm optimization (PSO) algorithm. Akjiratikarl et al. [15] applied this swarm-based heuristic algorithm to find a model for efficient scheduling of caregivers in Ukraine. Trautsamwieser et al. [20] proposed a variable neighborhood search (VNS) algorithm to address home healthcare logistics faced with a flood in Austria. A fast heuristic algorithm based on the Lagrangian relaxation theory was first proposed by Kergosien et al. [31] for routing and scheduling in a home healthcare logistics network. Different heuristic approaches with the Tabu search (TS) algorithm and genetic algorithm (GA) were implemented in the model proposed in [29] by Liu et al. [11]. In another work, Liu et al. [12] also offered a hybrid metaheuristic algorithm as a variation of TS based on the strategies of local searches in another paper.

In 2012, a study by Rasmussen et al. [21] presented another single-depot, single-period home healthcare scheduling model to minimize the total cost of caregiver routes. They considered a branch-and-price algorithm to solve this complicated model. Nickel et al. [32], as one of the first studies, proposed a multi-period home healthcare model under uncertainty. To create a plan for caregivers’ routing and scheduling, they considered policies for master and operational scheduling based on a case study in Germany. In 2014, a model for a multi-service home healthcare problem was proposed by Mankowska et al. [22]. Their model, in a multi-period environment, provided the possibility of different caregivers visiting patients. They also assumed that each caregiver had two different skills simultaneously. In 2015, different types of vehicles such as cars and public transportation were proposed for a multi-service home healthcare system addressed by Hiermann et al. [26]. By solving a case study in Austria and a set of random tests, they offered a two-stage solution algorithm. The first step was a heuristic algorithm based on constraint programming and the second one was three capable metaheuristic algorithms, i.e., memetic algorithm (MA), simulated annealing (SA), and scatter search (SS). In another similar paper, Fikar and Hirsch [23] proposed a new home healthcare logistics model with different vehicles and time windows. They assumed that the caregivers were able to visit the patients by walking. Frifita et al. [27] introduced, for the first time, the term synchronized visits. A VNS algorithm was utilized to solve their home healthcare problem.

Recently, in addition to the problem of total cost, new objectives have been considered in the literature [31,32,33,34,35,36]. As one of the earliest works, Braekers et al. [9] proposed a bi-objective multi-period home healthcare system with satisfaction levels for the patients. They considered a heuristic algorithm based on dynamic programming to solve this practical problem. Fikar et al. [36] developed a multi-objective multi-depot home healthcare logistics model. The total cost and time of the home healthcare services were their objectives. They showed that home healthcare decisions should be made in real time and based on this reason, a discreet event-driven metaheuristic algorithm was proposed. Except Fikar et al., [36], no paper [37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63] has applied a real-time optimization to monitor and control the home healthcare services and caregivers’ routing and scheduling decisions.

Recently, researchers have become interested in the uncertainties associated with home healthcare logistics [7,25,28,32,48,63]. In this regard, Shi et al. [7] proposed a fuzzy environment for home healthcare logistics. They applied a hybrid GA as their solution method. In another work, the same authors [28] considered the stochasticity of travel and service times. In addition to SA and the GA, they applied the bat algorithm (BA) and cuckoo search (CS) to do a comparison among different metaheuristic approaches. Next, [48] proposed a robust optimization model for single objective, single-depot and single-period home healthcare routing with time windows. SA and VNS were considered to solve their robust optimization model. [56] also added synchronized visits and green emissions in another study that used a hybrid of TS and SA to solve the model. In another recent work, Fathollahi-Fard et al. [59] proposed a robust optimization model and included working time balancing and continuity of care in their home healthcare model. They proposed different heuristic and metaheuristic algorithms including red deer algorithm, PSO, GA, and MA for solving their model. Shahnejat-Bushehri et al. [62] developed another robust optimization model for home healthcare routing and scheduling with temporal dependencies to visit the patients. They integrated a Mont Carlo simulation with SA and MA for solving their model.

Environmental sustainability is another home healthcare logistics research area, which was first considered by Fathollahi-Fard et al. [38]. They proposed a bi-objective, single-depot and single-period home healthcare logistics model. With four heuristic algorithms based on routing of the caregivers, they found a tradeoff between the total cost and environmental pollution. They also proposed a hybridization of SA and the salp swarm algorithm (SSA). In another work, the concept of green home healthcare logistics was studied using a location-routing model. The same authors [52] proposed clustering patients in their home healthcare logistics model. To find Pareto-based solutions, five modifications of SA were developed. The recommendations of both researchers encouraged studying a triple bottom line approach for home healthcare logistics in order to consider the economic, environmental, and social impacts in an uncertain environment, which had never been studied in the literature.

Recent studies continue to focus on the development of new solutions for the home healthcare logistics [37,38,39,40,41,42,43,44,45,46,47,48,49,50]. For example, Fathollahi-Fard et al. [40] proposed a solution, based on the Lagrangian relaxation theory, for home healthcare routing with a travel balancing supposition. Lin et al. [24] solved a multi-service home healthcare routing and scheduling model using a combination of a harmony search algorithm (HSA) with a GA. Decerle et al. [42] developed a hybrid MA and ant colony optimization (ACO) solution for home healthcare logistics that considered time windows, synchronization, and working time balancing. In another important recent study, three efficient heuristic algorithms based on a combination of VNS and SA algorithms with the Lagrangian relaxation theory were studied by Fathollahi-Fard et al. [45]. Their home healthcare logistics network was a single-depot and single-period routing optimization model with travel balancing using different types of vehicles. More recently, Frifita et al. [53] introduced a number of modified VNS algorithms for a multi-depot home healthcare routing model with working time balancing for the caregivers and time windows. In another interesting paper, an adaptive memory search was applied to a recently developed metaheuristic algorithm, the so-called social engineering optimizer, by Fathollahi-Fard et al. [54]. They considered a balance between the total cost and satisfaction of patients. Their model was formulated using a fuzzy algorithm. At last, but not least, Fathollahi-Fard et al. [58], in another work, applied bi-level programming for the first time in the area of home healthcare logistics. They considered the possibility of outsourcing home care services in a multi-depot home healthcare logistics network. Liu et al. [64] proposed a multi-period home healthcare routing and scheduling problem considering synchronized visits and satisfaction of patients. They proposed an efficient metaheuristic algorithm combining the dynamic programming method. It goes without saying that another bi-level optimization that considered patient satisfaction was developed by Chen et al. [65]. They proposed a three-stage metaheuristic approach that combined an iterated local search with a large-neighborhood search and heuristic algorithm for solving large-scale datasets.

We collected the most revenant studies by conducting a search in Scopus source to identify research gaps and the differences between our study as compared with the aforementioned studies, as shown in

Table 1. Literature review of home healthcare services and logistics.

Reference	Number of Objectives		Number of Depots		Number of Periods		Uncertainty			Suppositions
	Single Objective	Multi-Objective	Single Depot	Multi-Depot	Single Period	Multi-Period	Fuzzy	Stochastic	Robust	Synchronized Visits	Multi-Service	Satisfaction	Working Time Balancing	Patient Clustering	Green Emissions	Solution Algorithm
[17]	√	-	√	-	√	-	-	-	-	-	-	-	-	√	-	SDSS
[18]	√	-	√	-	√	-	-	-	-	-	-	-	-	√	-	Heuristic
[15]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	PSO
[19]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	Exact
[20]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	VNS
[21]	√	-	√	-	√	-	-	-	-	√	-	-	-	-	-	Exact
[11]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	GA and TS
[12]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	Heuristic
[22]	√	-	√	-	√	-	-	-	-	√	√	-	-	-	-	Heuristic
[23]	√	-	√	-	√	-	-	-	-	√	√	-	-	-	-	VNS
[9]	-	√	√	-	-	√	-	-	-	-	√	√	-	-	-	Heuristic
[7]	√	-	√	-	√	-	√	-	-	-	-	-	-	-	-	GA
[24]	√	-	√	-	√	-	-	-	-	√	√	-	-	-	-	GA and HSA
[25]	√	-	√	-	√	-	-	√	-	-	√	-	-	-	-	Heuristic
[26]	√	-	√	-	√	-	-	-	-	√	√	-	-	-	-	MA, VNS, SA
[27]	√	-	√	-	√	-	-	-	-	√	-	-	-	-	-	VNS
[28]	√	-	√	-	√	-	-	√	-	-	-	-	-	-	-	GA, SA, CS, BA
[29]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	Exact
[30]	√	-	√	√	-	√	-	-	-	-	-	-	-	√	-	DSS
[31]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	Exact
[32]	√	-	√	-	-	√	-	√	-	-	-	-	-	√	-	Heuristic
[34]	√	-	√	-	-	√	-	-	-	√	√	-	√	√	-	Exact
[35]	√	-	-	√	-	√	-	-	-	√	-	-	-	-	-	Exact
[36]	-	√	√	-	-	√	-	-	-	√	-	-	-	-	-	Event-driven-based metaheuristic
[37]	√	-	√	-	-	√	-	-	-	-	√	-	-	√	-	Decomposition method
[38]	-	√	√	-	√	-	-	-	-	-	-	-	-	-	√	Heuristics and metaheuristics
[39]	√	-	√	-	√	-	-	-	-	-	√	-	-	-	-	MA
[40]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	Lagrangian relaxation method
[41]	-	√	√	-	-	√	-	-	-	-	√	√	-	√	-	Epsilon constraint method
[42]	√	-	√	-	√	-	-	-	-	√	√	-	√	√	-	Hybrid of MA and ACO
[43]	-	√	√	-	√	-	-	-	-	-	√	√	-	-	-	NSGA-II
[44]	√	-	√	-	√	-	-	-	-	-	√	√	-	-	-	Multi-objective MA
[45]	√	-	√	-	√	-	-	-	-	-	-	-	-	-	-	Lagrangian method and heuristics
[46]	√	-	√	-	-	√	-	-	-	√	√	√	-	√	-	Lagrangian and VNS methods
[47]	√	-	-	√	-	√	-	√	-	-	√	-	-	√	-	Heuristics
[48]	√	-	√	-	√	-	-	-	√	-	-	-	-	-	-	Mont Carlo simulation and heuristics
[49]	√	-	√	-	-	√	-	-	-	√	√	-	-	√	-	Heuristic
[50]	√	-	-	√	√	-	-	-	-	-	-	-	-	√	-	Heuristics
[51]	√	-	-	√	-	√	-	-	-	√	√	-	-	√	-	Heuristics
[52]	-	√	-	√	-	√	-	-	-	-	-	-	-	√	√	SA
[53]	√	-	√	-	√	-	-	-	-	√	-	-	√	-	-	VNS
[54]	-	√	-	√	-	√	√	-	-	-	-	√	-	√	-	Metaheuristics
[55]	√	-	-	√	-	√	-	-	-	-	-	-	-	-	-	Benders’ decomposition
[56]	√	-	√	-	√	-	-	-	√	√	√	-	-	-	√	Hybrid of TS and SA
[57]	√	-	-	√	-	√	-	√	-	√	-	-	-	-	-	Branch and cut algorithm
[58]	√	-	-	√	-	√	-	-	-	-	√	-	-	√	-	Heuristics and metaheuristics
[59]	-	√	-	√	√	-	-	-	√	-	√	-	√	√	-	Red deer algorithm
[60]	√	-	-	√	-	√	-	√	-	√	-	√	-	√	-	Metaheuristic algorithm
[61]	-	√	-	√	-	√	-	-	-	-	√	-	√	-	-	Modified social engineering optimizer
[62]	√	-	-	√	√	-	-	-	√	-	-	√	-	-	-	Mont Carlo simulation with SA and MA
[63]	√	-	-	√	√	-	-	-	-	-	-	-	-	-	-	Hybrid GA
[64]	√	-	-	√	-	√	-	-	-	√	-	-	√	-	-	Metaheuristics
[65]	√	-	-	√	-	√	-	-	-	√	-	√	-	-	-	Hybrid metaheuristic
This paper	-	√	-	√	-	√	-	-	√	-	√	-	√	√	√	Heuristics with Lagrangian relaxation based on epsilon constraint

. In total, 54 studies were reviewed. The studies were classified as: multi-objective or single objective, multi-depot or single-depot, and multi-period or single-period home healthcare models. If the studies considered uncertainty in their model, their model was classified to one of three groups based on fuzzy, stochastic, or robust optimization frameworks. The main new contributions of the aforementioned studies were classified according to the following six items: synchronized visits, multi-service, satisfaction levels, working time balancing, patient clusters, and green emissions of the home healthcare logistics. Finally, the solution algorithms, as one of the main novelties in the majority of the studies, are given in Table 1. On the basis of these criteria, the following findings are identified:

• Only nine papers ([9,36,38,41,43,52,54,59,61]) considered more than one objective function. However, they did not consider unemployment time and green emissions as their objectives, simultaneously.
• The simultaneous study of a multi-service, multi-period, and multi-depot home healthcare logistics model under uncertainty was rarely studied in the literature. A study by Fathollahi-Fard et al. [59] was the only study classified in this group.
• As summarized in Table 1, only nine studies included uncertainty in their models [7,9,25,28,32,47,48,54,56,57,59,60,62]. Among them, robust optimization was applied by four studies [48,56,59,62]. With the exception of the study by Shi et al. [56], no other study considered green emissions.
• The environmental sustainability for home healthcare logistics was considered by only two studies [38,52]; both of these studies did not consider uncertainty in their home healthcare logistics models.
• As summarized in Table 1, no study optimized green emissions, working time balancing, and multi-service suppositions, simultaneously.
• The majority of studies innovated new heuristic and metaheuristic algorithms; only a few of them considered the Lagrangian relaxation theory [40,45,46] or Benders’ decomposition [55]. However, no study applied an efficient Lagrangian relaxation-based heuristic approach based on the epsilon constraint method.

To fill these research gaps, in this study, we address an efficient robust optimization model for sustainable home healthcare logistics considering working time balancing and green emissions. To solve the proposed model, strong reformulations and heuristic algorithms based on Pareto-based Lagrangian relaxation are generated.

3. Proposed Model for the Problem

In this study, the problem was to find a robust and sustainable solution for the design of a home healthcare logistics and services network which included allocation, scheduling, and routing decisions. Here, the problem statement is introduced. Next, the main assumptions for the proposed problem are illustrated. Finally, a multi-objective robust optimization model to deal with uncertainty in home healthcare logistics and services based on a triple bottom line approach is developed.

3.1. Problem Statement

In line with other home healthcare studies [7,38,40,45,54], we consider a city with V patients, P pharmacies, and L laboratories. A general statement of the proposed problem is presented graphically in Figure 2. Regarding home healthcare logistics, the proposed model fits with a two-dimensional geographical plan. In this regard, the distances between the patients, pharmacies, and laboratories were defined. Before making routing and scheduling decisions, for all periods, first, the patients were assigned to pharmacies. This assignment creates patient clusters, which reduces transportation costs for the companies. This allocation is associated with the distance ($AC$), as noted earlier. Figure 2a shows an example of our allocation decision problem with 24 patients and three pharmacies. To address the needs of patients, each pharmacy employs $N_p$ caregivers composing four categories that include some doctors, nurses, physiotherapist, and nutritionist. Therefore, each caregiver is able to do one or more roles (${\beta }_{inp}^k$). These different caregivers create a multi-service home healthcare system. It goes without saying that this multi-service system confirms a full range of home care services. Accordingly, the patients need specific types of services per period (${\delta }_i^{kt}$). There is a variable cost for each home care service ($W{C}_{np}^{kt}$), while there is a fixed salary for each caregiver per period ($F{C}_{np}^{kt}$). On each day, each caregiver’s starting point is at a pharmacy and after some patient visits, the caregivers go back to a laboratory. As an example, these routing and scheduling activities are drawn in Figure 2b for period t and Figure 2c shows the period t + 1.

As our problem is an extension of routing optimization, each caregiver employs one vehicle. As such, each route requires a vehicle with a transportation cost per distance unit ($T{C}_n$) and, in this study, this vehicle has a specific fuel consumption rate (FER), i.e., 0.25 L and CO₂ emissions rate (CER), i.e., 2.61 kg CO₂/L. It should be noted that, since the vehicle load is the same in all nodes, these rates will not change. There is no additional passenger for each trip and the number of drugs for patients is very low and does not limit the capacity of the vehicles.

There is overtime for a caregiver, if the caregiver works over the maximum time in each period ($O{C}_{np}^{kt}$). The maximum working time of the caregivers in each period is fixed ($W_{\mathit{max}}$). As such, the caregivers may not be available for all periods (${\gamma }_{np}^t$). The main challenge of scheduling in home healthcare logistics is the uncertainty associated with travel and service times, as well as the time windows. With regards to the scenario s, service time is uncertain and depends on patients’ health records and types of diseases ($W_{ki}^{ks}$). The travel time per each scenario is varied ($T_{\mathit{ijs}}$). At last, but not least, a time window limits the availability of each patient based on the earliest ($E{P}_{is}^t$) and the latest time ($L{P}_{is}^t$).

Based on the aforementioned problem statement, a scenario-based robust optimization model is developed for a multi-objective home healthcare problem considering a triple bottom line approach.

3.2. Assumptions

The proposed problem follows these assumptions to build a multi-objective robust optimization model for sustainable home healthcare logistics.

• The model has three objectives, i.e., to optimize the total cost, environmental pollution and the total unemployment time, simultaneously.
• The proposed model is also classified as a multi-depot, multi-period and multi-service home healthcare problem.
• To deal with uncertainty, travel and service times in addition to time windows are uncertain. To address these uncertain parameters, a scenario-based robust optimization theory is employed.
• As a multi-service model, the caregivers include doctors, nurses, physiotherapists, and nutritionists. In a multi-service concept, each caregiver can support one or more roles.
• The model assumes that the locations of pharmacies and laboratories are fixed. Its number is also the same. Therefore, there is no facility location planning in our healthcare network.
• The demand of patients must be met in their specified time window. A shortage is not allowed for the home healthcare services.
• The caregivers are not available for all shifts (time periods). It means that each caregiver has a rest in some time periods.
• Each patient has a time window to show their availability in each period.
• The total cost includes a fixed fee for the employment of caregivers, a variable cost for each home care service, transportation costs, and overtime costs. These factors help our model to achieve economic sustainability.
• To achieve environmental sustainability, the emerged CO₂ for the transportation of caregivers is minimized by the second objective function.
• To achieve social sustainability, the model considers the working time balancing constraints and minimizes the unemployment time and overtime of the caregivers in each period.

3.3. Formulation

We have defined all notations in the Supplementary Materials to reduce the length of this paper. This study deploys a novel multi-objective scenario-based, robust optimization model. In this section, first, we introduce the concept of robust optimization. Then, we state our objectives and constraints for our sustainable home healthcare logistics network design problem.

• Robust optimization

In this study, we utilize the concept of robust optimization proposed firstly by Mulvey et al. [66]. To illustrate this method, assume a minimization model as follows:

$$min Z_s=fy+c_s{x}_s$$

(1)

where $Z_s$ is the total cost per scenario. For binary variables ($y$) as the fixed costs, we have defined the coefficient of $f$, while $c_s$ is the coefficients for continuous variables ($x_s$) which are dependent on the probabilistic scenario. Mulvey et al. [66] updated the model in Equation (1) to the define the robust optimization as follows:

$$Z=\mathit{min}(\lambda {\displaystyle \sum }_{s\in S}{\pi }_s{Z}_s+\left(1-\lambda \right){\displaystyle \sum }_{s\in S}{\pi }_s{(Z_s-{\displaystyle \sum }_{s^{\prime }\in S}{\pi }_{s^{\prime }}{Z}_{s^{\prime }})}^2)$$

(2)

The objective function is divided into two parts. The first part calculates the expected total cost based on the probability of each scenario (${\pi }_s$), while the second part calculates the variance of each scenario. To show the importance of each part, $\lambda$ is a weight factor. The second part of this model is nonlinear and makes it difficult. Mulvey et al. [66] called this model robust optimization and defined the following constraints which are the budget constraint set:

$$Ty+A{x}_s\ge b_s \forall s\in S$$

(3)

As given in Relation (3), for binary variables, we have defined the coefficient $T$, while $A$ is the coefficient of continuous variables. The budget as the right number is defined as $b_s$ which is dependent on each scenario. Similar to the concept of two-stage stochastic programming, binary variables are first-stage variables and continuous variables are second-stage variables in the robust optimization. This classification refers to the fact that the binary variables are not related to the scenarios. However, other variables with regards to probabilistic scenarios are estimated by the information about the uncertain parameters which are correlated to the occurrences of the scenarios.

One main disadvantage of the robust optimization proposed by Mulvey et al. [66], is the nonlinear term in the objective function. In this regard, Leung et al. [67] provided an update to the formulation of Mulvey and made it linear. They defined an objective function which was linear with one auxiliary variable. They found that their model was easier than the model by Mulvey et al. [66]. Their update to the robust optimization concept is:

$$Z=\mathit{min}({\displaystyle \sum }_{s\in S}{\pi }_s{Z}_s+\lambda {\displaystyle \sum }_{s\in S}{\pi }_s(Z_s-{\displaystyle \sum }_{s^{\prime }\in S}{\pi }_{s^{\prime }}{Z}_{s^{\prime }}+2{\theta }_s))$$

(4)

where ${\theta }_s$ represents an auxiliary variable.

Based on this new objective function, our constraints are:

$$T’y+A’x_s\ge b{’}_s \forall s\in S$$

(5)

$$Z_s-{\displaystyle \sum }_{s^{\prime }\in S}{\pi }_{s^{\prime }}{Z}_{s^{\prime }}+{\theta }_s\ge 0 \forall s\in S$$

(6)

Equation (5) is the budget constraint which is similar to the model by Mulvey et al. [66], while Equation (6) ensures that there is a non-negative standard deviation for each scenario. Using the above robust optimization framework, the proposed multi-objective model for sustainable home healthcare logistics is established as follows:

• Economic sustainability

Equation (7) is the first objective function for achieving economic sustainability. The goal is to minimize the total cost ($f_1$). This objective includes six terms; four terms are for the first-stage variables and two other terms are for the second-stage variables. The first and second terms are the assignment costs to generate the cluster for each pharmacy. The transportation cost of routing caregivers is given in the third term. The fourth term is the fixed cost of home healthcare services. The last term supports the variable cost for the working time of the caregivers and the overtime costs which are uncertain and estimated by the robust optimization concept.

$$\begin{array}{c} min f_1={\displaystyle {\displaystyle \sum }_{l\in l}}{\displaystyle {\displaystyle \sum }_{p\in P}^{ }}AC.D_{pl}^L.Y_{pl}^L\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{i\in V}^{ }}{\displaystyle {\displaystyle \sum }_{p\in P}^{ }}AC.D_{ip}^P.Y_{ip}^P\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{p\in P}^{ }}{\displaystyle {\displaystyle \sum }_{l\in L}^{ }}{\displaystyle {\displaystyle \sum }_{i\in V}^{ }}{\displaystyle {\displaystyle \sum }_{j\in V}^{ }}{\displaystyle {\displaystyle \sum }_{n\in N_p}^{ }}{\displaystyle {\displaystyle \sum }_{t\in T}^{ }}{\displaystyle {\displaystyle \sum }_{k\in K}^{ }}T{C}_{np}.D_{ij}^V.X_{ijnpl}^{tk}+{\displaystyle {\displaystyle \sum }_{p\in P}^{ }}{\displaystyle {\displaystyle \sum }_{i\in V}^{ }}{\displaystyle {\displaystyle \sum }_{n\in N_p}^{ }}{\displaystyle {\displaystyle \sum }_{t\in T}^{ }}{\displaystyle {\displaystyle \sum }_{k\in K}^{ }}F{C}_{np}^{kt}.Y_{inp}^{kt}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\displaystyle \sum }_{s\in S}^{ }}{\pi }_s{f}_{1s}+\lambda {\displaystyle {\displaystyle \sum }_{s\in S}^{ }}{\pi }_s\left(f_{1,s}-{\displaystyle {\displaystyle \sum }_{s^{\prime }\in S}^{ }}{\pi }_{s^{\prime }}{f}_{1,s’}+2{\theta }_{1,s}\right)\hfill \end{array}$$

(7)

• Environmental sustainability

Sustainable transportation refers to vehicles which have a low impact on the environment [38,52]. Since transportation is the heart of a home healthcare logistics network, environmental pollution should be minimized by our optimization model. In this regard, we focus on the minimization of CO₂ emissions. $G{E}_{ij}$ is the amount of CO₂ emissions for traveling from node i to node j. The calculation of CO₂ emissions $\left(G{E}_{ij}=CER\times FCR\times D_{\mathit{ij}}^V\right)$ is related to three factors, i.e., the type of fuel, the rate of fuel consumption in a vehicle, and the travel distance between two nodes. Therefore, the second objective function ($f_2$), which is not related to uncertainty, is to minimize the environmental pollution for the transportation in our home healthcare logistics network, which is given in Equation (8):

$$min f_2={\displaystyle \sum }_{p\in P}{\displaystyle \sum }_{l\in L}{\displaystyle \sum }_{i\in V}{\displaystyle \sum }_{j\in V}{\displaystyle \sum }_{n\in N_p}{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}G{E}_{ij}.X_{ijnpl}^{tk}$$

(8)

• Social sustainability

Social sustainability for home healthcare services refers to the satisfaction of caregivers [59]. One way to improve the satisfaction of caregivers is to achieve working time balancing for the caregivers. Therefore, the third objective ($f_3$), as given in Equation (9), minimizes the unemployment time of each caregiver which is uncertain and, accordingly, the robust optimization is adopted as follows:

$$min f_3={\displaystyle \sum }_{s\in S}{\pi }_s{f}_{3s}+\lambda {\displaystyle \sum }_{s\in S}{\pi }_s\left(f_{3,s}-{\displaystyle \sum }_{s^{\prime }\in S}{\pi }_{s^{\prime }}{f}_{3,s’}+2{\theta }_{3,s}\right)$$

(9)

• Constraints

Equations (10)–(30) provide the constraints for our multi-objective robust optimization model as follows:

$$f_{1,s}-{\displaystyle \sum }_{s^{\prime }\in S}{\pi }_{s^{\prime }}{f}_{1,s^{\prime }}+{\theta }_{1,s}\ge 0 \forall s\in S$$

(10)

$$f_{1,s}={\displaystyle \sum }_{i\in P}{\displaystyle \sum }_{n\in N_p}{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}W{C}_{np}^t.W_{is}^{kt}.Y_{inp}^{kt}+{\displaystyle \sum }_{p\in P}{\displaystyle \sum }_{n\in N_p}{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}O{C}_{np}^{kt}{O}_{nps}^{kt} \forall s\in S$$

(11)

$$f_{3,s}-{\displaystyle \sum }_{s^{\prime }\in S}{\pi }_{s^{\prime }}{f}_{3,s^{\prime }}+{\theta }_{2,s}\ge 0 \forall s\in S$$

(12)

$$f_{3,s}={\displaystyle \sum }_{p\in P}{\displaystyle \sum }_{n\in N_p}{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{Z}_{nps}^{kt} \forall s\in S$$

(13)

$$S{T}_{inps}^t+W_{is}^{kt}+T_{ijs}\le S{T}_{jnps}^t+M\left(1-X_{ijnpl}^{tk}\right) \forall i,j\in V,i\ne j,p\in P,l\in L,n\in N_p,t\in T,k\in Ks\in S$$

(14)

$$E{P}_{is}^t\le S{T}_{inps}^t\le L{P}_{is }^t\forall i\in V,p\in P,n\in N_p,k\in K,t\in T ,s\in S$$

(15)

$${{\displaystyle \sum }}_{k\in K}{O}_{nps}^{kt}\ge \mathit{max}\left\{{{\displaystyle \sum }}_{i\in V}{{\displaystyle \sum }}_{k\in K}{Y}_{inp}^{kt}\cdot W_{is}^{kt}\cdot {\beta }_{inp}^k\cdot {\delta }_i^{kt}-W_{\mathit{max}},0\right\}\hspace{1em}\hspace{1em}\forall p\in P,n\in N_p,t\in T,s\in S$$

(16)

$${{\displaystyle \sum }}_{k\in K}{Z}_{nps}^{kt}\le \mathit{max}\left\{W_{\mathit{max}}-{{\displaystyle \sum }}_{i\in P}{{\displaystyle \sum }}_{k\in K}{Y}_{inp}^{kt}\cdot W_{is}^{kt}\cdot {\beta }_{inp}^k\cdot {\delta }_i^{kt},0\right\} \forall p\in P,n\in N_p,t\in T,s\in S$$

(17)

$${\displaystyle \sum _{i\in V}{\displaystyle \sum _{j\in V}{\displaystyle \sum _{n\in N_p}{\displaystyle \sum _{k\in K}{\displaystyle \sum _{t\in T}{X}_{ijnpl}^{tk}}}}}}\le M\cdot Y_{pl}^L\hspace{1em}\hspace{1em}\forall p\in P,l\in L$$

(18)

$${\displaystyle \sum _{l\in L}{Y}_{pl}^L}=1\hspace{1em}\hspace{1em}\forall p\in P$$

(19)

$${\displaystyle \sum _{p\in P}{Y}_{pl}^L}=1\hspace{1em}\hspace{1em}\forall l\in L$$

(20)

$${\displaystyle \sum _{\begin{array}{l}j\in V\\ i\ne j\end{array}}{\displaystyle \sum _{n\in N_p}{\displaystyle \sum _{k\in K}{\displaystyle \sum _{l\in L}{\displaystyle \sum _{t\in T}{X}_{ijnpl}^{tk}}}}}}\le M\cdot Y_{ip}^P\hspace{1em}\hspace{1em}\forall p\in P,i\in V$$

(21)

$${\displaystyle \sum _{p\in P}{Y}_{ip}^P}=1\hspace{1em}\hspace{1em}\forall i\in V$$

(22)

$${\displaystyle \sum _{\begin{array}{l}j\in V\\ i\ne j\end{array}}{\displaystyle \sum _{l\in L}{X}_{ijnpl}^{tk}}}=Y_{inp}^{tk}\hspace{1em}\hspace{1em}\forall i\in V,n\in N_p,p\in P,t\in T,k\in K$$

(23)

$${\displaystyle \sum _{p\in P}{\displaystyle \sum _{n\in N_p}{Y}_{inp}^{kt}\cdot {\beta }_{inp}^k}}={\delta }_i^{kt}\hspace{1em}\hspace{1em}\forall i\in V,k\in K,t\in T$$

(24)

$${\displaystyle \sum _{\begin{array}{l}j\in V\\ i\ne j\end{array}}{X}_{ijnpl}^{tk}}-{\displaystyle \sum _{\begin{array}{l}j\in V\\ i\ne j\end{array}}{X}_{jinpl}^{tk}}=0\hspace{1em}\hspace{1em}\forall i\in V,p\in P,l\in L,n\in N_p,k\in K,t\in T$$

(25)

$${\displaystyle \sum _{i\in V}{X}_{i0npl}^{tk}}=1\hspace{1em}\hspace{1em}\forall j\in V,p\in P,l\in L,n\in N_p,k\in K,t\in T$$

(26)

$${\displaystyle \sum _{j\in V}{X}_{0jnpl}^{tk}}=1\hspace{1em}\hspace{1em}\forall i\in V,p\in P,l\in L,n\in N_p,k\in K,t\in T$$

(27)

$${\displaystyle \sum _{\begin{array}{l}i\in V\\ i\ne j\end{array}}{\displaystyle \sum _{j\in V}{X}_{ijnpl}^{tk}}}\le M\cdot {\gamma }_{np}^t\hspace{1em}\hspace{1em}\forall l\in L,p\in P,n\in N_p,k\in K,t\in T$$

(28)

$$X_{ijnpl}^{tk},Y_{inp}^{kt},Y_{ip}^P,Y_{pl}^L\in \left\{0,1\right\}\hspace{1em}\hspace{1em}\forall i,j\in V,p\in P,l\in L,n\in N_p,k\in K,t\in T$$

(29)

$$S{T}_{inps}^t,O_{nps}^{kt},Z_{nps}^{kt},{\theta }_{1,s},{\theta }_{2,s}\ge 0 \hspace{1em}\hspace{1em}\forall i,j\in V,p\in P,l\in L,n\in N_p,k\in K,t\in T,s\in S$$

(30)

To summarize the role of each constraint, Equations (10) and (11) are the robust optimality limitations for the first objective function. As such, Equations (12) and (13) are the robust optimality limitations for the third objective. Equation (14) shows the relation of routing and scheduling decisions for the caregivers regarding their starting time to visit the patients and the availability of the patients. Equation (15) confirms that each patient has a time window limitation. Equations (16) and (17) are the working time balancing constraints to calculate the overtime and unemployment time of the caregivers, respectively. Equation (18) provides a link between the allocation decisions with routing and scheduling decisions. Equations (19) and (20) state that each pharmacy works with only one laboratory. Equation (21) is similar to Equation (18) and shows that a caregiver for pharmacy p can visit the patient i based on the assignment of patients and pharmacies. Similar to Equations (19) and (20), Equation (22) means that each patient should be assigned to only one pharmacy. Equation (23) shows that the assignment of patients and the routing and scheduling decisions of the caregivers are correlated. Equation (24) ensures that there is no shortage in the model and the demand of each caregiver for the type of the homecare services must be met. Equation (25) ensures that each caregiver visits a patient, and then leaves this patient. Equations (26) and (27) state that, for each caregiver, the starting point is a pharmacy and the end point is a laboratory, respectively. Equation (28) states that the availability of caregivers has a direct impact on the routing and scheduling decisions of the caregivers. Finally, Equation (29) supports the binary variables and Equation (30) shows the non-negative continuous decision variables.

In conclusion, for the first time, we study a multi-objective robust optimization model for a sustainable home healthcare logistics network. This model has not been studied before and, in the following, we provide a plan to solve it mathematically.

4. Proposed Solution

The proposed solution for the model presented in Section 3 is based on a combination of the epsilon constraint method to generate Pareto-based solutions and the Lagrangian relaxation theory with three heuristic algorithms to find an efficient solution in reasonable time. One benefit of our solution algorithm is the opportunity to find an optimal solution which is very close to the global solution. It should be noted that in the case of multi-objective optimization, there is no global solution and there is a set of non-dominated solutions which can dominate other solutions.

To define the concept of multi-objective optimization, assume that we only consider two objectives for our model including f₁ and f₂. Next, note that there are two solutions for our model: Solution A and Solution B. Let $Z_1^A$ and $Z_2^A$ be the values for the first (f₁) and second (f₂) objective functions in Solution A and, as such, Solution B has $Z_1^B$ and $Z_2^B$. In the case of minimization, for both objectives, Solution A dominates Solution B when $\left(\begin{array}{c}{Z}_1^A\\ Z_2^A\end{array}\right)\le \left(\begin{array}{c}{Z}_1^B\\ Z_2^B\end{array}\right)$ and $\left(\begin{array}{c}{Z}_1^A\\ Z_2^A\end{array}\right)\ne \left(\begin{array}{c}{Z}_1^B\\ Z_2^B\end{array}\right)$ [59,61].

Here, first, we introduce the epsilon constraint method to show how our proposed solution algorithm performs the Pareto-based analysis. Then, the Lagrangian relaxation method [68] is applied to our optimization model given in Section 3.3. Next, we illustrate our procedures for the generation of initial solutions by heuristic algorithms. Finally, the main loop of our algorithm is explained to apply all the epsilon constraint, Lagrangian relaxation, and heuristic methods in an integrated way.

4.1. Epsilon Constraint Method

The epsilon constraint method was first proposed by Haimes et al. [69]. On the basis of this method, we consider one objective as the main objective and other objectives are limited by favorable bounds. The model given in Equations (7)–(30) can be transformed into a single objective model as follows (Equations (10)–(30)):

$$\begin{array}{c}{\mathit{min} f}_1\hfill \\ s.t.\hfill \\ f_2\le e_1\hfill \\ f_3\le e_2\hfill \\ f_{2\mathit{min}}\le e_1\le f_{2\mathit{max}}\hfill \\ f_{3min}\le e_2\le f_{3max}\hfill \end{array}$$

(31)

where $e_1$ and $e_2$ are the limit bounds for the second and the third objective functions. These bounds are limited by the positive and negative ideal solutions for each objective function. To achieve the positive ideal solution for the second objective function ($f_{2min}$), we only optimize the second objective function without the limitations of other objective functions. Regarding the negative ideal solution ($f_{2max}$), the model given in Equation (31) without the bound limits should be run. As such, the limitations for the bound of the third objective function, ($f_{3min}, f_{3max}$), are calculated in the same way as we calculated the positive and negative ideal solutions for the second objective function.

In each run, the bounds are updated, and the solutions are noted. For our model, we have considered one cut for the bounds as the average of the positive and negative ideal solutions. Therefore, the epsilon constraint method, in this case, can generate nine solutions maximally (if all cases are feasible). These cases are given in

Table 2. All possible cases for the Pareto-based analysis by the epsilon constraint method for our optimization model.

No. of Solutions	${\mathit{e}}_1$	${\mathit{e}}_2$
S1	$f_{2\mathit{min}}$	$f_{3\mathit{max}}$
S2	$f_{2\mathit{max}}$	$f_{3\mathit{min}}$
S3	$f_{2\mathit{min}}$	$f_{3\mathit{min}}$
S4	$f_{2\mathit{max}}$	$f_{3\mathit{max}}$
S5	$f_{2\mathit{min}}$	$(f_{3\mathit{max}}+f_{3\mathit{min}})/2$
S6	$f_{2\mathit{max}}$	$(f_{3\mathit{max}}+f_{3\mathit{min}})/2$
S7	$(f_{2\mathit{max}}+f_{2\mathit{min}})/2$	$f_{3\mathit{min}}$
S8	$(f_{2\mathit{max}}+f_{2\mathit{min}})/2$	$f_{3\mathit{max}}$
S9	$(f_{2\mathit{max}}+f_{2\mathit{min}})/2$	$(f_{3\mathit{max}}+f_{3\mathit{min}})/2$

.

After an assessment of these solutions based on the definition of the multi-objective optimization concept, the non-dominated solutions would be selected to provide a full set of alternatives for decision-makers.

4.2. Lagrangian Relaxation Theory

This study uses the concept of the Lagrangian relaxation theory for solving the proposed routing and scheduling optimization model. According to the Lagrangian relaxation theory, first, we generate an efficient reformulation for the main model, and then, we apply an iterative algorithm for improving lower and upper bounds for the main optimization model. The applied iterative algorithm can be classified as a sub-gradient algorithm for optimization problems [68]. In the area of home healthcare logistics optimization, there are several studies which have considered the Lagrangian relaxation theory [40,45,46]. The proposed hybrid algorithm differs from previous works [40,45,46,61,62] as it is combined with the epsilon constraint method and three heuristic algorithms.

The Lagrangian relaxation theory uses a lower bound and an upper bound. These two bounds have been updated in each iteration. First of all, in this method, the main optimization model should be relaxed by removing a number of very difficult constraints from it. The Lagrangian relaxation reformulation is evaluated by two main criteria, i.e., the computational time for this reformulation and the quality of objective function based on its deviation from the original model.

These relaxed constraints are multiplied by the Lagrange multiplier in the objective function. For the case of minimization, after solving the relaxed model, a lower bound would be generated. This solution may be infeasible but optimal. To check the quality of this lower bound, we need to define an upper bound for our Lagrangian relaxation reformulation. In the majority of Lagrangian relaxation-based algorithms, this upper bound is selected randomly which is inefficient and may increase the computational cost [40,45,62]. However, the proposed algorithm uses three fast heuristic algorithms to find an upper bound to define an initial upper bound which is feasible but not optimal.

To select the best reformulation for our optimization model, we conducted some experiments which are not reported here. From our tests to find the most difficult constraint set, Equation (14) as a big-M type constraint [70] makes the model more complex and increases the computational time significantly. Therefore, the proposed reformulation as our relaxed model is (Equations (10)–(13); Equations (15)–(30)):

$$\begin{array}{c}{\mathit{min} f}_1^{’}=f_1+{\displaystyle {\displaystyle \sum }_{i\in V}}{\displaystyle {\displaystyle \sum }_{j\in V}}{\displaystyle {\displaystyle \sum }_{s\in S}}{\mathsf{\pi }}_{\mathit{ijs}}\left({\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{l\in L}}{\displaystyle {\displaystyle \sum }_{t\in T}}{\displaystyle {\displaystyle \sum }_{k\in K}}{\displaystyle {\displaystyle \sum }_{n\in N_p}}{\mathit{ST}}_{\mathit{inps}}^t+W_{\mathit{is}}^{\mathit{kt}}+T_{\mathit{ijs}}-{\mathit{ST}}_{\mathit{jnps}}^t+M\left(1-X_{\mathit{ijnpl}}^{\mathrm{tk}}\right)\right)\hfill \\ s.t. \hfill \\ f_2\le e_1\hfill \\ f_3\le e_2\hfill \\ f_{2\mathit{min}}\le e_1\le f_{2\mathit{max}}\hfill \\ f_{3\mathit{min}}\le e_2\le f_{3\mathit{max}}\hfill \end{array}$$

(32)

where ${\pi }_{ijs}$ are Lagrange multipliers. After solving this relaxed model, a lower bound is found. Now, we should check the quality of this lower bound by an upper bound.

4.3. Initialization

Solving Equation (32) provides a lower bound for the original model given in Equation (31). Although this lower bound is most probably infeasible, it is an optimal solution. Now, we want to define the search space for finding a feasible solution and consider a greedy algorithm based on different heuristic algorithms for decisions.

As reviewed in Table 1, many studies have contributed heuristic algorithms (greedy search) for solving home healthcare optimization models [9,22,32,38,45,49,51]. In this study, we extend the heuristic algorithms proposed by Fathollahi-Fard et al. [45] for our multi-objective robust and sustainable home healthcare logistics network model. The search space is based on the main decision variables for our optimization model including $X_{ijnpl}^{tk}$, $Y_{pl}^L$, and $Y_{ip}^P$. Among these variables, the assignment of pharmacies to laboratories ($Y_{pl}^L$) and patients to pharmacies ($Y_{ip}^P$) are the same for all periods and probabilistic scenarios. If we assume predefined values for $Y_{pl}^L={\overline{Y}}_{pl}^L$ and $Y_{ip}^P={\overline{Y}}_{ip}^P$, the main part of the search space is defined by feasible sets for our routing decision variable ($X_{ijnpl}^{tk}$). In this regard, we have defined three heuristic algorithms to define the relations in this variable. These relations are our search space to define a feasible solution which may be optimal.

To define $Y_{pl}^L={\overline{Y}}_{pl}^L$, each pharmacy is assigned to one laboratory regarding the minimum array for the following matrix:

$$DCS{L}_{pl}=D_{pl}^L\times AC;\forall p\in P;l\in L$$

(33)

where $DCS{L}_{pl}$ is the assignment cost for both the pharmacies and the laboratories. To allocate the patients to pharmacies ($Y_{ip}^P={\overline{Y}}_{ip}^P$), the proposed strategy considers the distance of patients to pharmacies ($D_{ip}^P$) and the distance of patients to laboratories ($D_{il}^{VL}$). To generate the cluster for each pharmacy, the average distance is considered. Therefore, the minimum array of each pharmacy is selected by:

$$DISP{L}_{ip}={\displaystyle \sum _{l\in L}(D_{ip}^P\times D_{il}^{VL}\times Y_{pl}^L)/2};\forall i\in V,p\in P$$

(34)

where $DISP{L}_{ip}$ is the matrix for the allocation of patients to form a pharmacy cluster. Based on Equations (33) and (34), the allocation decisions have been made. The routing and scheduling decisions ($X_{ijnpl}^{tk}$) are the difference of each heuristic. In this regard, the following matrix is generated regarding the caregivers and patients:

$$DT{C}_{ijnp}=D_{ij}\times T{C}_{np};\forall i,j\in V;p\in P,n\in N_p$$

(35)

where $DT{C}_{ijnp}$ is the total transportation cost.

We call our heuristic algorithms LGEC1, LGEC2, and LGEC3. More details about these algorithms are given in the Supplementary Materials file. All these heuristic algorithms follow a symmetric pattern for decisions on the allocation of patients and caregivers’ routing.

The difference of these algorithms is graphically shown in Figure 3 based on the symmetric pattern of patients. In this example, we have five patients, one caregiver and one home care service for all these patients, is requested. With regards to the first algorithm, i.e., LGEC1, the route for the caregiver is $\{3\to 2\to 4\to 1\to 5\}$, in which Patient 3 is the first node for the caregiver which is the closest patient to the pharmacy (Figure 3a). LGEC2 finds the optimal route of $\{2\to 4\to 1\to 3\to 5\}$ in which Patient 2 is selected based on the mean transportation cost to the closest patient (Figure 3b). The last route generated by the last version of our hybrid algorithm, i.e., LGEC3 is $\{1\to 4\to 2\to 3\to 5\}$. In this solution, Patient 1 who is the furthest patient to the laboratory, is selected (Figure 3c).

To show how we define a feasible solution in our search space for the routing decision variables ($X_{ijnpl}^{tk}$), consider an example with one pharmacy, one laboratory and 10 patients. This example with our computations, is shown in Figure 4. We have assumed that in our multi-service, only nurses and physiotherapists are available in this example. In this regard, we have one tour for the nurse and another tour for the physiotherapist. The transportation cost per unit is set to two. The first heuristic algorithm is LGEC. After running this algorithm, the nurse’s route starts with Patient 3. Next, the nurse does the home healthcare service for Patient 7, which is the patient closest to Patient 3. After the completion of this tour for visiting all patients needing a nurse, the tour is $0\to 3\to 7\to 1\to 8\to 4\to 0$. In a similar procedure, the physiotherapist’ route starts with Patient 5. The patient closest to Patient 5 is Patient 9. After visiting all patients, the complete route is $0\to 5\to 9\to 10\to 2\to 6\to 0$. The second heuristic algorithm for the routing decision, so-called LGEC2, computes the mean transportation cost for each patient. Based on this criterion, first, the nurse visits Patient 8. The route of the physiotherapist starts with Patient 10. Other visits are based on the lowest traveling distance. After the completion of this cycle for the nurse, the route is $0\to 8\to 4\to 3\to 7\to 1\to 0$. As such, the route for the physiotherapist is $0\to 10\to 9\to 2\to 6\to 5\to 0$. Finally, we run the third heuristic algorithm, so-called LGEC3. The first patient to be visited by the nurse is based on the criterion of furthest patient to the laboratory. In this regard, the route for the nurse is exactly the same as LGEC1 in this example. However, the physiotherapist starts with Patient 2. Hence, the complete cycle for the physiotherapist is $0\to 2\to 9\to 10\to 6\to 5\to 0$.

4.4. Proposed Algorithm

The main notations of our algorithm are as follows:

t	Number of iterations
${\mathit{LB}}_t$	Lower bound at iteration t
${\mathit{UB}}_t$	Upper bound at iteration t
${\mathsf{\pi }}_{\mathit{ijs}}^t$	Lagrange multiplier at iteration t
$f_1^t$	Optimal solution found in Equation (31) at iteration t
$f_t$	Updater for the Lagrange multiplier per iteration
α	Reduction rate of the updater
Maxt	Maximum number of iterations

Generally, the flowchart of the proposed hybrid heuristic algorithm is shown in Figure 5. The main loop is based on the following steps:

Step 0: Initialize the bounds of other objectives (e₁ and e₂) and Lagrange multiplier ${\pi }_{ijs}^0$ and set t = 0;

Step 1: Let ${\pi }_{ijs}={\pi }_{ijs}^t$ and solve the relaxed problem given in Equation (32). Replace this solution for the main model given in Equation (31). Then, update the lower bound as follows:

$$L{B}_t=\mathit{max}\left(L{B}_t,f_1^t\right)$$

(36)

Step 2: Select one of heuristic algorithms, i.e., LGEC1, LGEC2 and LGEC3. Then, find the upper bound.

Step 3: Update the Lagrange multiplier as follows:

$${\pi }_{ijs}^{t+1}=\mathit{max}\left({\pi }_{ijs}^t+{\mu }^t\times \left({\displaystyle \sum }_{p\in P}{\displaystyle \sum }_{l\in L}{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{k\in K}{\displaystyle \sum }_{n\in N_p}S{T}_{inps}^t+W_{is}^{kt}+T_{ijs}-S{T}_{jnps}^t+M\left(1-X_{ijnpl}^{tk}\right)\right),0\right)$$

(37)

where ${\mu }^t=f_t\times \left|\frac{L{B}_t-L{B}_{t+1}}{{\left(U{B}_{t+1}-L{B}_{t+1}\right)}^2}\right|$ and being $f_t$ a number distributed by $U(0,2)$ in the first iteration. Then, it would be reduced per iteration by $f_{t+1}=f_t\times \alpha$ with no improvement. Notably, $\alpha$ would be tuned by a range from 0.5 to 1;

Step 4: t = t + 1;

Step 5: If a feasible lower bound reaches or t satisfies the maximum number of iteration (Maxt) then, note this solution and go to Step 6. Otherwise, go to Step 1;

Step 6: If there is no possible case from Table 2, generate the best non-dominated solutions. Otherwise, go to Step 0.

In conclusion, the proposed hybrid algorithm has four input parameters including $f_t, \alpha , {\pi }_{ijs}^0$, and Maxt. These parameters would be tuned regarding the characteristics of the simulated test studies. The proposed hybrid heuristic algorithm is in three versions, so-called LGEC1, LGEC2 and LGEC3. The main innovation of the proposed hybrid heuristic algorithm is to use a near-optimal solution as the upper bound and the use of epsilon constraint method to generate the Pareto-based solutions.

5. Computational Results

Here, the proposed hybrid heuristic algorithm, i.e., LGEC1, LGEC2 and LGEC3, are applied. The calibrations of the model’s parameters with the three classifications of complexity levels and the details of these algorithms are given in the first subsection. Next, the validation and comparison of the results along with the sensitivity analyses are provided and, finally, a comprehensive discussion about the findings and managerial solutions is provided. It should be noted that all experiments for the coding of the algorithms and model were done in GAMS 24.7.3 software and MATLAB R2013a software on a laptop with Core^TM i5, 2.40 GHz, and RAM 4 GB and the use of the Windows 8 operating system.

5.1. Simulated Test Studies

As there is no available data base to match with our novel multi-objective robust optimization model for sustainable home healthcare logistics, we have generated our benchmarks based on the recent studies in the literature review [38,40,45,50,52]. In three classifications of the complexity levels, 12 simulated test studies are generated as given in

Table 3. Size of the simulated test studies.

Classification	Test Study	${\mathit{N}}_{\mathit{p}}$	$\mathit{P} \mathrm{and} \mathit{L}$	$\mathit{V}$	$\mathit{T}$	Maxt
Small	SP1	2	2	10	2	10
	SP2	3	2	25	4	10
	SP3	4	3	40	6	20
	SP4	4	3	65	8	20
Medium	MP5	6	3	80	14	30
	MP6	6	4	85	18	30
	MP7	6	5	95	24	30
	MP8	6	5	100	28	30
Large	LP9	8	6	120	32	40
	LP10	8	6	150	36	40
	LP11	8	7	160	40	50
	LP12	8	8	200	42	50

. For the calibration of the heuristic algorithms, the maximum iteration is tuned in Table 3 and its amount is 10 iterations for very small sizes and 50 iterations for the very large sizes. Other parameters of the algorithm, including $f_t, \alpha , \mathrm{and} {\pi }_{ijs}^0$, are calibrated as 2, 0.99, and 0.3, respectively. For the calibration of the robust optimization, the coefficient of the robust optimality ($\lambda$) is set as 0.5.

The range of deterministic and the scenario-based parameters is given in

Table 4. The range of deterministic parameters.

Parameters	Range
$(x_i,y_i)$	$1000\times (U(0,1),U(0,1))$
$(x_j,y_j)$	$1000\times (U(0,1),U(0,1))$
$(x_p,y_p)$	$1000\times (U(0,1),U(0,1))$
$(x_l,y_l)$	$1000\times (U(0,1),U(0,1))$
$D_{ij}^V$	$\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$ (meter)
$D_{ip}^P$	$\sqrt{{(x_i-x_p)}^2+{(y_i-y_p)}^2}$ (meter)
$D_{pl}^L$	$\sqrt{{(x_p-x_l)}^2+{(y_p-y_l)}^2}$ (meter)
$D_{il}^{VL}$	$\sqrt{{(x_i-x_l)}^2+{(y_i-y_l)}^2}$ (meter)
$T{C}_n$	$rand\left\{2,3,4,5\right\}$ ($)
$F{C}_{np}^{kt}$	$rand\left\{8,10,12,14,16\right\}$ ($)
$W{C}_{np}^t$	$rand\left\{0.4,0.5,0.6,0.7,0.8\right\}$ ($)
$O{C}_{np}^{kt}$	$rand\left\{1,2,3\right\}$ ($)
${\gamma }_{np}^t$	$rand\left\{0,1\right\}$
${\beta }_{inp}^k$	$rand\left\{0,1\right\}$
${\delta }_i^{kt}$	$rand\left\{0,1\right\}$
$W_{\mathit{max}}$	5 × 60 × 60 (min)
$AC$	2 ($)
FER	0.25 (L)
CER	2.61 (kgCO2/L)

and

Table 5. The range of uncertain parameters.

Parameter	General Scenarios
	Optimistic	Realistic	Pessimistic
${\mathit{EP}}_{\mathit{is}}^t$	rand{100, 351, …, 450} (min)	rand{450, 451, …, 600} (min)	rand{600, 601, …, 800} (min)
${\mathit{LP}}_{\mathit{is}}^t$	rand{500, 501, …, 850} (min)	rand{851, 851, …, 1100} (min)	rand{1100, 1101, …, 1320} (min)
$T_{\mathit{ijs}}$	$\frac{D_{\mathit{ij}}^V}{{{\displaystyle \sum }}_{i\in V}{{\displaystyle \sum }}_{j\in V}{D}_{\mathit{ij}}^V}\times 70$ (min)	$\frac{D_{\mathit{ij}}^V}{{{\displaystyle \sum }}_{i\in V}{{\displaystyle \sum }}_{j\in V}{D}_{\mathit{ij}}^V}\times 90$ (min)	$\frac{D_{\mathit{ij}}^V}{{{\displaystyle \sum }}_{i\in V}{{\displaystyle \sum }}_{j\in V}{D}_{\mathit{ij}}^V}\times 110$ (min)
$W_{\mathit{is}}^{\mathit{kt}}$	$U\left(10,20\right)$ (min)	$U\left(15,25\right)$ (min)	$U\left(20,30\right)$ (min)

, respectively. It should be noted that our planning horizon starts from zero minute and ends from one day to 21 days. From per period, eight hours are limited. Each caregiver may work each day in two periods, i.e., morning and evening.

5.2. Extensive Analyses and Sensitivities

In our experiments, the proposed heuristic algorithms are validated by the epsilon constraint method in small sizes. The results of SP1 as an example are given in

Table 6. Results of SP1.

Epsilon Constraint			LGEC1			LGEC2			LGEC3
${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$	${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$	${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$	${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$
7136.2	225.4	6.5	7235.81	226.84	6.75	7462.1	206.28	6.45	7392.2	226.84	6.5
7185.7	314.2	2.15	7256.42	314.72	2.15	7456.2	308.2	2.09	7416.8	312.6	2.1
7919.5	226.2	2.1	8546.17	221.4	2.1	8219.3	218.8	2.12	7682.5	224.1	2.14
7056.8	314.2	6.7	7182.81	314.72	6.75	7367.8	309.15	6.66	7543.2	312.8	6.12
7208.2	223.8	4	7365.6	226.4	4.1	7513.6	216.13	4.35	7481.5	226.84	3.85
7286.6	318.72	4.2	7482.54	314.72	4.35	7455.8	314.72	4.25	7518.3	309.72	3.96
7438.2	265.8	2.06	7518.25	260.7	2.15	7677.4	268.44	2.08	7602.4	256.8	2.11
7808.4	259.92	6.34	7648.32	265.8	6.75	7523.6	254.24	6.56	7416.2	249.71	6.34

Note: The best values are shown in bold.

. The non-dominated solutions are depicted in Figure 6. Among nine possible solutions, as given in Table 2, eight solutions are considered to be the non-dominated solutions for each algorithm in Table 6. It should be noted that the units of the objective functions are $, kg CO₂, and minutes, respectively.

As given in Figure 6, the solutions of LGEC2 outperform the other algorithms’ solutions, as all objectives are in minimization form. The solutions of LGEC1 and LGEC3 are almost the same.

For each solution of the heuristic algorithms, the gap between the lower bound and the upper bound is computed as follows:

$$\mathrm{Gap}=\frac{upper bound-lower bound }{lower bound }\times 100$$

(38)

Clearly, this gap is for the first objective function, and we provided an average for all non-dominated solutions from our heuristic algorithms, i.e., LGEC1, LGEC2, and LGEC3. We have also provided the average of the objectives in

Table 7. Results of algorithms (CPU = computational time (seconds).

Algorithms		SP1	SP2	SP3	SP4	MP5	MP6	MP7	MP8	LP9	LP10	LP11	LP12
LGEC1	$f_1$	7529.49	20,987.2	33,822	71,300	124,528.6	167,349.6	217,454.4	294,784	327,136	414,576	510,480	600,516
	$f_2$	270.8	754.8099	1216.417	2564.322	4478.702	6018.77	7820.802	10,601.98	11,765.53	14,910.33	18,359.54	21,597.71
	$f_3$	4.45	13.5	32.43	55.43	83.95	79.61	106.02	127.25	119.45	155.85	169.06	167.040
	Gap	0.078	0.0892	0.1205	0.1403	0.1662	0.0992	0.1205	0.1853	0.1102	0.2081	0.2382	0.2891
	CPU	98.1036	115.0362	236.4876	712.467	1472.815	1692.641	4346.857	6942.648	12,344.1	15,450.8	15,919.45	18,827.49
LGEC2	$f_1$	7584.475	20,632.8	32,875.8	68,304.8	114,875.6	147,252.6	219,043.2	279,784.4	342,624	421,524	499,760	559,776
	$f_2$	261.995	712.731	1135.6482	2359.4931	3968.2157	5086.6335	7566.5386	9664.7578	11,835.463	14,560.953	17,263.505	19,336.673
	$f_3$	4.32	15.5	30.86	58.67	79.44	78.39	118.45	121.21	130.27	149.05	163.67	176.74
	Gap	0.054	0.0652	0.0924	0.1105	0.1402	0.1482	0.1682	0.1402	0.1654	0.1572	0.2103	0.2571
	CPU	98.1684	114.1407	234.3996	711.5886	1472.855	1691.582	4353.564	6938.969	12,339.14	15,445.88	15,923.21	18,831.67
LGEC3	$f_1$	7506.638	20810	33822	69,802.4	124,528.6	15,7301.1	217454.4	287,284.2	334,880	418,050	505,120	580,146
	$f_2$	264.926	734.43132	1193.6538	2463.4824	4394.8921	5551.5067	7674.4508	10,138.9	11,818.662	14,753.917	17,826.812	20,474.646
	$f_3$	4.14	14.2	32.43	57.81	83.95	78.04	106.02	125.63	124.82	151.17	165.19	172.81
	Gap	0.0852	0.1206	0.1043	0.1402	0.1892	0.0865	0.1205	0.1782	0.1273	0.1672	0.2981	0.3452
	CPU	98.9874	115.7796	236.565	713.2833	1475.178	1695.827	4354.003	6939.268	12,339.1	15,454.67	15,922.76	18,830.81

Note: The best values are shown in bold.

as well as the average gap between the lower bound and the upper bound and the computational time. The best values in this table are shown in bold.

As given in Table 7, LGEC2, in most of the test problems, outperforms the other algorithms. From the average of the solutions, most of the LGEC2′s solutions are the best. The gap between the lower bound and the upper bound for each solution of the heuristic algorithms is also important to confirm the performance of the heuristic algorithms; most of the results of LGEC2 are shown in bold as the best gap. Based on the criterion of CPU time, in half of test studies, LGEC2 is the best and in the others, LGEC1 is the best point.

The behavior of the algorithms’ gaps is depicted in Figure 7. The gaps of the heuristic algorithms are the same and there are significant competition. However, in most of the test studies in small sizes and large ones, LGEC2 is the best and outperforms other algorithms. After this algorithm, LGEC3 is the best, except in one test in which LGEC1 is better than LGEC3 (i.e., LP9).

The computational time behavior of the algorithms is depicted in Figure 8. The first issue is that the CPU time of the hybrid heuristic algorithm is lower than the epsilon constraint method. In addition, all the heuristic algorithms have similar behavior based on this criterion. However, as previously noted in Table 7, both LGEC1 and LGEC2 confirm the best behavior with regards to the solution time.

Finally, some sensitivity analyses are performed to evaluate the key parameters of he model to gain practical insights. To this end, the best heuristic algorithm in our study, i.e., LGEC2, is selected. For all sensitivity analyses, SP2 is selected as a small test study. The behavior of the objectives, on average, of all non-dominated solutions is reported in the results.

First, we analyze the impact of robust optimality $(\lambda$) on the objectives. This parameter is set as 0, 0.5, 0.75, 1, and 1.5 in our sensitivity analyses. The results are given in

Table 8. Sensitivity analysis of robust optimality.

$\mathit{\lambda }$	0	0.5	0.75	1	1.5
$f_1$	15,640.4	20,632.8	23,129	25,625.5	30,617.6
$f_2$	712.731	712.731	712.731	712.731	712.731
$f_3$	11.5	15.5	17.5	19.5	21.5

and the behavior of the objective functions based on their normalized values is shown in Figure 9. These analyses show the impact of uncertainty on our robust optimization model. With an increase in robust optimality, the amount of the objective functions increased, except for the second objective function which is not linked with uncertainty.

More sensitivity analyses are performed using four parameters which include transportation cost ($T{C}_n$), fixed cost of caregiver employment ($F{C}_{np}^{kt}$), service cost of the caregivers ($W{C}_{np}^{kt}$), overtime cost of the caregivers ($O{C}_{np}^{kt}$), and the maximum time of working for each caregiver ($W_{\mathit{max}}$). All of these parameters have a high impact on the behavior of the objectives, and they are very important for managers of home healthcare logistics networks.

Without a doubt, transportation is one of the main sectors for a sustainable home healthcare logistics network. Based on the results given in

Table 9. Sensitivity analysis of transportation cost.

$\mathit{T}{\mathit{C}}_{\mathit{n}}$	3	4	5	6	7
$f_1$	15,999.14	20,632.8	25,266.45	29,900.11	34,533.76
$f_2$	712.731	712.731	712.731	712.731	712.731
$f_3$	15.5	15.5	15.5	15.5	15.5

and the behavior of the normalized values of the objectives depicted in Figure 10, there are no changes in the amount of environmental emissions and unemployment time of the caregivers, while the transportation cost increases. This parameter has a high impact on the financial costs of the home healthcare logistics model and reduces the routes of the caregivers to visit the patients.

A sensitivity analysis of the fixed cost of caregivers’ employment is done by an increase in this parameter and the behavior of the objective functions is analyzed. The results are given in

Table 10. Sensitivity analysis of the fixed employment cost of the caregivers.

$\mathit{F}{\mathit{C}}_{\mathit{n}\mathit{p}}^{\mathit{k}\mathit{t}}$	8	9	10	11	12
$f_1$	21,054.93	20,632.8	20,680.62	21,674.92	20,789.45
$f_2$	712.731	712.731	731.824	744.82	751.835
$f_3$	15.5	15.5	14	13	12.5

and Figure 11. This parameter results in some changes in the behavior of all objectives in the majority of the sensitivity analyses. An increase in the cost of fixed employment of the caregivers generally increases the total cost. However, it also effects the routes of the caregivers and the number of employed caregivers is reduced. This reduction increases the carbon emissions of the home healthcare logistics network as the route of each caregiver has been generally increased. This factor leads to a reduction in the unemployment time of the caregivers as they must work more to satisfy the demands of the patients.

The service cost of home healthcare is analyzed using some sensitivities. The results are given in

Table 11. Sensitivity analysis of the service cost of the caregivers.

$\mathit{W}{\mathit{C}}_{\mathit{n}\mathit{p}}^{\mathit{k}\mathit{t}}$	0.3	0.4	0.5	0.6	0.7
$f_1$	22,420.38	20,632.8	21,705.58	21,919.39	22,214.07
$f_2$	712.731	712.731	712.731	712.731	712.731
$f_3$	15.5	15.5	15.5	15.5	15.5

and the behavior of the objectives in terms of normalized values is depicted in Figure 12. Service cost is a variable cost for caregivers associated with each home care service. Similar to transportation cost, service cost is only linked with financial issues and there is no impact on environmental emissions and unemployment time, as confirmed by the results. Regarding financial costs, if the service cost fits the model’s parameters, it can reduce the total cost. However, the general behavior of the first objective function is to be increased while the service cost increases.

The sensitivity analysis of the overtime cost is another important factor which is linked with uncertainty. The results are given in

Table 12. Sensitivity analysis of the overtime cost.

$\mathit{O}{\mathit{C}}_{\mathit{n}\mathit{p}}^{\mathit{k}\mathit{t}}$	2	3	4	5	6
$f_1$	20,492.09	20,632.8	20,890.04	20,791.55	20,531.98
$f_2$	712.731	712.731	712.731	712.731	712.731
$f_3$	15.5	15.5	15.5	14	14.5

and the behavior of the objectives is depicted in Figure 13, based on the normalized values of the results. Overtime cost generally increases the total cost. However, an increase in this factor finally leads to a well-tuned level of this parameter to reduce the overtime of the caregivers. After this calibration, unemployment time is also reduced. It shows that the model achieves working time balancing, aimed at social justice for the caregivers. Although this parameter has a significant impact on the economic and social objectives, there is no change in the carbon emissions of the home healthcare logistics model based on the results.

Finally, the maximum working time in each period is another important factor to study the working time balancing of the caregivers. A sensitivity analysis of this parameter was conducted and the results are reported in

Table 13. Sensitivity analysis of the maximum time of working.

${\mathit{W}}_{\mathbf{max}}$	25	50	75	100	125
$f_1$	21,671.19	20,632.8	22,595.39	22,080.3	22,080.3
$f_2$	712.731	712.731	712.731	682.916	654.185
$f_3$	0	15.5	18.7	42.711	67.562

. The behavior of sustainability dimensions is also shown in Figure 14. The results confirm that this factor has a high impact on a sustainable home healthcare logistics model with economic, environmental, and social goals. When this parameter increases, the third objective regarding unemployment time increases. It shows that a proper amount of this factor leads to a well-tuned level of unemployment time as its ideal value. However, this parameter slightly reduces green emissions. When unemployment time increases, overtime reduces. It may also reduce a caregiver’s route. The behavior of the financial cost is varied and, in general, it is increased. In conclusion, a well-tuned level of this parameter is useful to achieve a sustainable home healthcare logistics model.

5.3. Managerial Implications

As the ageing problem worsens, home healthcare is becoming increasingly important in the functioning of sustainable development goals, especially a triple bottom line approach for sustainability. The services provided by a home healthcare logistics network include housekeeping, nursing, drug deliveries, and cleaning, as well as different types of the caregivers have been studied economically and environmentally [38,52], but are yet to be improved on the social dimension (Table 1).

In this study, we proposed a sustainable home healthcare logistics model based on a triple bottom line approach to optimize financial, environmental, and social goals, simultaneously (Figure 1). Since some parameters linked with time and availability are uncertain, a scenario-based robust optimization approach was developed. This multi-objective robust optimization model characterized with real-life constraints such as clustering of patients, working time balancing, multi-services, and time windows, was solved using a hybrid heuristic algorithm that combined the epsilon constraint method with the Lagrangian relaxation theory and used three heuristic algorithms (Figure 5).

Without a doubt, this study was more complex than the majority of studies on home healthcare logistics. In addition to the total cost and carbon emissions of transportation, unemployment time was minimized as the third objective function. With the use of a multi-objective robust optimization model, working time balancing was considered to find a tradeoff between unemployment time and overtime and multi-service home healthcare was applied to provide a full range of care by using doctors, nurses, physiotherapists and nutritionist. The validation of the developed heuristic algorithms, i.e., LGEC1, LGEC2 and LGEC3, as shown in Table 6, showed that the second version of the heuristic algorithm, i.e., LGEC2 outperformed the other algorithms based on the non-dominated solutions (Figure 6). The performance and applicability of this algorithm, based on the criteria of quality (Figure 7) and time (Figure 8), are shown in Table 7, which lends support for the development of our Pareto-based Lagrangian relaxation-based heuristic algorithm in other optimization problems.

One significant managerial implication is that this study conceptually shifts from green home healthcare logistics management to sustainable home healthcare logistics management based on working time balancing and environmental pollution of transportation. Other practical solutions are taken from the results of the sensitivity analyses. The first point is that the proposed robust optimization model is very sensitive to the robust optimality factor (Figure 9). This factor supports both the favorable and the unfavorable scenarios to find a more feasible solution. This reason satisfies the optimality of the deterministic solution as compared with a robust model.

The transportation cost (Figure 10) and the variable cost of home care services (Figure 12) are important to achieve economic sustainability. In addition, overtime cost effects economic as well as social sustainability factors (Figure 13). However, the fixed employment cost of the caregivers to change the number of caregivers in each period (Figure 11) and the maximum time of working (Figure 14) play key roles to achieve a triple bottom line approach of sustainability. Therefore, it is very challenging for managers of home healthcare logistics to find a well-tuned amount of these parameters with regards to other real-life limitations to achieve sustainability for a home healthcare logistics network.

6. Conclusions and Future Remarks

In this paper, a multi-objective robust optimization model was applied to a practical routing and scheduling optimization problem for the application of a home healthcare logistics network. The needs and benefits derived from an efficient design for a robust and sustainable healthcare logistics network are highly significant to create a response or action to the COVID-19 pandemic. This study focused on home healthcare logistics and services for planning the routing and scheduling of caregivers to visit patients’ homes. Due to the need for social distancing during the COVID-19 pandemic, these services are highly applicable to reduce the growth of the epidemic. In this regard, in this study, we addressed the challenge of the COVID-19 pandemic for healthcare systems. In addition, the home healthcare logistics network was redesigned to meet the standards of a triple bottom line approach based on sustainable development goals. Hence, this model, for the first time, formulated economic, environmental, and social criteria simultaneously to address a triple bottom line approach in home healthcare logistics and services.

The total cost and the environmental pollution of the home healthcare logistics and services were optimized in addition to the unemployment time of the caregivers. Three heuristic algorithms based on the epsilon constraint method and the Lagrangian relaxation theory were developed to address the proposed model. Extensive validation, comparison, and sensitivity analyses were performed and highlighted that the second version of our heuristic algorithms, so-called LGEC2, was very efficient and that the parameters of this model, especially the fixed employment cost and maximum working time of the caregivers, were significant to achieve a sustainable home healthcare logistics network based on economic, environmental, and social goals.

In this study, for the first time, although. we conceptualized a sustainable home healthcare logistics network, the proposed multi-objective robust optimization model was still very general and many other suppositions could be added to enhance its applicability. Considering Industry 4.0 and industrial informatics for monitoring patients using advanced technologies have rarely been modeled in this research area [71]. Patients’ satisfaction and care continuity are two important factors to increase the quality of home care services. An introduction of fuzzy logic with the possibility of programming instead of a robust optimization makes another contribution for future work [72]. The proposed heuristic algorithms were another significant contribution of this study and, therefore, the application of some novel and state of the art metaheuristic algorithms such as the social engineering optimizer [73], is very interesting to solve a very large real dataset for the proposed model as compared with our heuristic algorithms. Without a doubt, the proposed heuristic algorithms with the use of the Lagrangian relaxation theory and the epsilon constraint method, provide an introduction for solving other multi-objective routing and scheduling optimization problems.