Can we work more safely and healthily with robot partners? A human‐friendly robot–human‐coordinated order fulfillment scheme

Abstract

This work addresses the issue of intelligent robot–human‐coordinated parts‐to‐picker order fulfillment carried out in a human‐friendly manner. One unique feature of the proposed approach involves integrating a real‐time data‐driven stochastic‐dynamic model with a fatigue accumulation function. The optimal solutions help achieve coordination between human pickers and robots, such that robots can agilely adapt to the coordinated pickers’ efficiency and fatigue conditions. Specifically, the proposed method estimates human pickers’ instantaneous performance and robot queue lengths, which are then fed back in real time as indexes to adjust robots’ speeds of handling racks and moving them to human pickers. Using data that are provided by a giant electronic commerce (e‐commerce) company, our analyses demonstrate that the proposed robot–picker coordination system permits alleviating a picker's fatigue without much influence on picking efficiency. In particular, a picker's accumulated fatigue can be reduced by 53.74% at the expense of lowering picker efficiency by 14.79% if the proposed robot–picker coordination system is applied in the focal firm of the study case. Through our scenario design and sensitivity analysis, additional findings and managerial insights, including the rules of human‐friendly robot behaviors for coordination with human pickers in different operational scenarios are provided. They facilitate the development of “human‐friendly” intelligent robot‐human‐coordinated order‐fulfillment systems for intelligent logistics operations.

INTRODUCTION

Business operations have entered the Industry 4.0 era in which robotics are widely applied in the production, logistics, and e‐commerce companies (Choi et al., 2022; Loffler et al., 2022; Olsen & Tomlin 2020; Rai et al., 2021). This is evidenced by the claims made by Amazon.com:

It's better for everybody…Workers no longer would have to walk massive warehouse floors to find the right power drill–instead, robots would bring the drill directly to them… The hour or more it took to process a package had been shaved down to as little as 15 minutes…Robots, Amazon insists, are good for workers… “They make the job safer,”… (Evans, 2020)

However, picking three times faster also implies more wear and tear due to repetitive motion and working faster at lifting and handling products… So along with the drive to automate more warehouse tasks comes much higher expectations for workers…The robots have raised the average picker's productivity from around 100 items per hour to what Mr. Long and others have said is a target of around 300 or 400, though the numbers vary across teams and facilities… (Del Rey, 2019).

Indeed, there is no doubt that advancements in such disruptive technologies as robotics and related applications (Azadeh et al., 2019; Chung, 2021; Perera et al., 2020; Shi et al., 2021; Wang et al., 2022) have brought new benefits to electronic commerce (e‐commerce) companies (such as Amazon, Alibaba, and Jingdong) by increasing the efficiency of order fulfillment. However, they are simultaneously creating new safety challenges.

Order picking is a critical step in customer order fulfillment (Batt & Gallino, 2019; de Koster et al., 2007; de Vries et al., 2016; Frazelle, 2002). In a typical warehouse, it is not only costly and time‐consuming but also labor‐intensive. In general, the cost of order picking is estimated to account for more than 50% of the total warehousing cost (Tompkins et al., 2010). Furthermore, order picking typically relies on manual labor to retrieve multiple items from storage, sort them, and then package them to fulfill customer orders (C.‐M. Chen et al., 2010; de Koster et al., 2007; de Vries et al., 2016). Owing to the large variety of customer orders and stock‐keeping units (SKUs) in most e‐commerce operations, each order‐picking task for a customer order is difficult to replicate. The ordering picking's efficiency and accuracy depend mainly on the human pickers’ experience (Batt & Gallino, 2019) and performance (de Vries et al., 2016). Moreover, in e‐commerce, agile order picking can be more challenging when the ordered products are diverse, demand is fluctuating, customer returns are increasing and customers expect more flexible and expeditious logistics services (Frazelle, 2002). Briefly, order picking, which depends on intensive labor and pickers’ performance, remains crucial to the success of fulfilling customer orders with high service quality, regardless of the methods/technologies that are deployed in a warehouse.

Despite a variety of innovative automation technologies, such as autonomous mobile robots, have been increasingly introduced to carry out order fulfillment jointly with human pickers in warehouses (Banker, 2016; Tobe, 2018; Wang et al., 2022), “parts‐to‐picker” order fulfillment systems raise several new issues on humanity and safety. This creates a challenge on robot–human coordination. A certain number of robots that carry racks automatically from the storage area to the order picking area may be deployed to facilitate the picking tasks of human pickers. Poor robot–human coordination typically increases the number of robots that carry, queue and await handling by pickers in the order‐picking area. Note that a human picker's efficiency is not as controllable as that of a robot and may vary with time and psychophysical condition. Thus, an unexpected delay in picking items may inevitably occur if robots just continue to carry racks to her and are unable to respond to her performance. More seriously, recent reports have revealed a substantial rising injury rate in automated warehouses (Al Elew & Oh, 2020; Doll, 2020; Evans, 2020). For example, Amazon's warehouse injury rates have increased every year since 2016 after they began using robots in parts‐to‐picker order fulfillment. The serious injury rate in 2019 was 7.7 injuries per 100 employees, which was 33% higher than in 2016 and nearly double the corresponding recent industry standard (Evans, 2020). Order picking often requires pickers to perform repetitive tasks in awkward postures for a whole day, easily causing musculoskeletal disorders (Grosse et al., 2015; Lavender et al., 2012), which accounts for up to 33% of all injuries and illnesses in the United States in 2013 (Bureau of Labor Statistics, 2014).

Motivated by the aforementioned robot–human coordination issues in parts‐to‐picker order fulfillment practices, this paper aims to answer the following questions.

How should human pickers and coordinated robots react and adapt to each other to fulfill orders in a human‐friendly manner considering robot–human collaboration?

How does such a human‐friendly robot–human coordination scheme affect the performance of an order fulfillment system in a parts‐to‐picker order fulfillment center? Specifically, what are the critical factors and how do they influence the systems’ performance from humanity and safety perspectives?

To address the above research questions, we build a stochastic optimal control‐based model and develop a data‐driven solution approach. We propose a robot–picker‐coordinated order fulfillment mechanism under which the robot and picker efficiencies have the following characteristics: (i) The robot and picker efficiencies can be effectively controlled and coordinated if real‐time data concerning newly arriving racks are available and used. (ii) The robot efficiency heavily depends on real‐time data concerning newly arriving racks, whereas picker efficiency highly depends on the length of the queue of unprocessed racks. (iii) The picker‐to‐robot relative performance (RP) from the perspective of either efficiency or effectiveness during peak hours is less than that during either normal or off‐peak hours. The picker‐to‐robot RP is best in the “off‐peak” scenario. Moreover, we find that the robot efficiency is not the most important metric. Rather, human efficiency dominates the efficiency of a parts‐to‐picker order fulfillment system, particularly during peak hours. Accordingly, controlling the number of queued racks by slowing down their delivery by robots to the picking–packaging area is in fact a more efficient measure than slowing down the rack‐processing of a picker (in order to alleviate that picker's fatigue); this is especially prominent during a peak‐hour period. Last but not least, employing data that were provided by a well‐established e‐commerce company, our analyses demonstrate that a picker's accumulated fatigue can be reduced by over 50% at the expense of lowering picker efficiency by around 15% if the proposed robot–picker coordination system is applied.

The rest of this paper is organized as follows. Section 2 reviews the relevant literature to elucidate its contribution to the field of innovative logistics. Section 3 describes the problem and characterizes the proposed robot–picker‐coordinated order fulfillment system. Section 4 presents a proposed stochastic‐dynamic optimal control model. Section 5 develops a real‐time data‐driven approach to estimating and controlling the state variables of the proposed robot–picker‐coordinated order fulfillment system. Section 6 presents the computational results of a practice‐based real data analysis. Section 7 draws conclusions, provides managerial implications, and discusses future research. To enhance exposition, all technical proofs are put in an Online Appendix.

LITERATURE REVIEW

Recently, in the Industry 4.0 era, robotic technologies (such as Amazon robotics and Jingdong robots) that are used by e‐commerce companies (Luo & Choi, 2022) for parts‐to‐picker order fulfillment have emerged (Banker, 2016; Shi et al., 2021; Tam, 2014; Tobe, 2018). They have revolutionalized the way of order fulfillment operated in warehouses. Most of the extant literature on parts‐to‐picker order fulfillment systems focuses mainly on increasing the warehousing efficiency utilizing diverse methodologies. For example, Enright and Wurman (2011) identified several parts‐to‐picker order fulfillment‐related resource allocation problems, including order allocation and robot allocation problems. Recently, Boysen et al. (2017) proposed a mixed integer programming model to minimize the number of rack visits to a stationary picking station in a mobile robot‐based order‐picking problem. The authors solved the optimization problem by first decomposing the problem into two subproblems of rack sequencing (for a given order sequence) and order sequencing (for a given rack sequence). By contrast, Lamballais et al. (2017) proposed four queuing network models associated with various warehousing layouts and robot zoning strategies to analyze the performance of an automated storage and parts‐to‐picker system with the objectives of optimizing robot utilization, maximizing order throughput, and minimizing order cycle time. Different from the above optimization‐based analytical models, Bozer and Aldarondo (2018) utilized a simulation‐based method to evaluate and compare the performances of two types of parts‐to‐picker picking systems (called “miniload” and “Kiva” systems) in the aspects of expected throughputs and expected container retrieval times in order processing. Some of the subsequent literature, including Weidinger et al. (2018) and Yuan et al. (2019), aimed at exploring the problem of pod/rack storage, which is the antecedent of parts‐to‐picker order fulfillment, where the inventory of items is stored and spread over multiple mobile pods that are carried by robots moving between storage and picking (stowing) zones. More recently, Wang et al. (2022) utilized an “approximate dynamic programming” based branch‐and‐price approach to solve the optimal robot scheduling problem for parts‐to‐picker order fulfillment systems. Different from the above literature, which solely focuses on efficiency, Wang et al. (2022) incorporated fluctuations of the working states of human pickers into the model, and hence the human factor was explored. However, human safety was not yet examined by Wang et al. (2022). In addition, Boldrer et al. (2022) proposed a hierarchical framework to address the problem of multi‐robot navigation in human‐shared working environments. Even though the approach of Boldrer et al. (2022) does not focus on parts‐to‐picker order fulfillment systems, the authors’ idea of conceptualizing a safe and socially aware navigation in their proposed methodological framework is noteworthy and consistent with our research goal.

Despite remarkable advances made to increase efficiency in parts‐to‐picker order fulfillment, issues such as safety and harmony in the human–robot interactions remain under‐explored in relevant intralogistics operations and related areas (Azadeh et al., 2019; Chen et al., 2022; Enright & Wurman, 2011). As argued by Enright and Wurman (2011), systems optimization at the high level, considering dual objective functions for both workers and robots, remains challenging as these two objectives may not be compatible. Details of the operational features and challenges of parts‐to‐picker order fulfillment can also be found in Enright and Wurman (2011). Drawing from experimental results, Chen et al. (2022) further suggested the urgent necessity of developing novel robotic motion methods to ensure human safety for applications of human—robot coordination in warehouses. Note that Azadeh et al. (2019) comprehensively reviewed the literature on robotized and automated warehouse systems. They pointed out that issues related to human—machine interaction in automated warehousing environments have not been fully addressed.

Another stream of relevant literature concerns the association of human factors with order‐picking performance (Batt & Gallino, 2019; Grosse et al., 2015, 2017; Loske, 2022; de Vries et al., 2016). As argued by Grosse et al. (2017), human factors determine the performance of an order‐picking system, which is time‐consuming and labor‐intensive, but they are often ignored in management‐oriented research. Incorporating human factors into order‐picking models for intralogistics operations and warehouse management is, thus, indispensable—particularly with respect to production, operations, and logistics management (de Koster et al., 2007). Observe that the related literature can be divided into two categories: the first focuses on relevant aspects of human cognition (e.g., human learning) to improve productivity/efficiency, and the second cares about the welfare and health of workers. The corresponding literature reviews are detailed below.

Prior studies on cognition have investigated its impact on order‐picking performance. Some of them (Grosse & Glock, 2015) integrated the concept of human learning (characterized by learning curves) with mathematical models to study how pickers’ performance can be improved by learning from experience. de Vries et al. (2016) conducted a field experiment to examine three types of picker‐to‐parts order‐picking methods under two incentive schemes for pickers with different regulatory foci (“prevention focus” vs. “promotion focus”). The work of de Vries et al. (2016) is pioneering in “elaborately” aligning order‐picking methods, incentive mechanisms, and regulatory focus to bridge the gaps between organizational behavior theories and warehousing practices. Batt and Gallino (2019) empirically estimated the effects of several factors, including a picker's learning from the experience of walking and searching processes, on picking time for order fulfillment in online women's apparel retailing. One remarkable feature of their work is the empirical demonstration of the effects of pickers’ heterogeneous learning capabilities, characterized by different learning curves, on order‐picking performance. Rather recently, Loske (2022) combined parametric and nonparametric approaches to analyze how the interactions between humans, machines, and intelligent software impact human learning and perception of work characteristics in the transition to an automated order‐picking system. The author's empirical findings verified that the real‐time feedback provided by the order‐picking system can facilitate human learning by doing tasks in the perception–cognition–motor–action cycle, thus clarifying the need for a human‐centered work system design.

The literature on the welfare and health of order‐picking workers aims to develop integrated models that consider both ergonomic (such as human energy expenditure and fatigue) and economic (such as picking time) performances in order picking (Grosse et al., 2015, 2017). For example, Grosse et al. (2015) proposed a conceptual framework that considers human factors in four critical categories (viz., the “perceptual, mental, physical, and psychosocial”). The authors suggested that human factors can be incorporated into the planning models to improve the performance of order‐picking systems. More recently, Glock, Grosse, Abedinnia, and Emde (2019) proposed an integrated model to determine the optimal processing sequence for orders, pallet rotation, and picker routing to minimize total picking effort. The authors considered the spinal loads on pickers and consequent risks of injury. Similar efforts had been made previously (Glock, Grosse, Kim, et al., 2019; Larco et al., 2017).

Despite the fact that the idea of integrating human factors into order‐picking decision support systems is promising and necessary for human wellbeing, the above‐cited order‐picking studies on human factors are concerned mostly with manual order picking (picker‐to‐parts order picking) rather than parts‐to‐picker order‐picking systems, which are the focus of our study. To be specific, a parts‐to‐picker order‐picking system has the following operational features that differentiate it from a traditional picker‐to‐parts order‐picking system, which relies mainly on manpower.

First, a parts‐to‐picker ordering picking system relies substantially on the coordination between pickers and robots to fulfill customer orders. Mobile robots lift racks that store SKUs and transport them from the storage area to stationary pickers to facilitate subsequent picking and packaging tasks by those pickers. In such parts‐to‐picker picking environments, robots must move forward and backward between the storage area and stationary pickers, forming closed loops between the storage area and the locations of the coordinated pickers. Accordingly, the efficiency of order fulfillment in a parts‐to‐picker picking system is contingent jointly upon the coordination between robots and pickers in carrying out the two consecutive tasks of handling racks (by mobile robots) and packaging (by pickers). Nevertheless, human–robot coordination often generates coordinating complexity and uncertainty as a result of the difficulty of sharing information, communication, and mutual adaption by the coordinated dyad (Faraj & Xiao, 2006).

Second, from a managerial perspective, an environment that includes working robots/machines is likely to differ from the one that includes only human labor. Psychological effects (e.g., stress and strains) have attracted attention since computer‐integrated/computer‐aided/robot‐aided manufacturing technologies were introduced for factory automation (Karuppan & Schniederjans, 1995; Rosenthal, 1984). As argued by Rosenthal (1984), the most difficult problems in automating a factory are managerial rather than technical. Olsen and Tomlin (2020) further highlighted the issue of managing worker–machine interfaces when robotics and artificial intelligence were introduced to convert manual operations into lean, digitized, and highly automated operations. Therefore, the issues that are raised by working with robots/machines seem not to be limited to those around technology and productivity. Instead, in the post‐automation era, new psychological (Karuppan & Schniederjans, 1995), behavioral (Hinds, 2004), and organization‐behavioral (Barrett et al., 2012; Beane & Orlikowski, 2015) issues should be crucial for operations managers.

Considering the above research gaps in the literature on parts‐to‐picker order fulfillment and intralogistics operations, this work contributes to the field of intelligent logistics operations in the following three ways.

First, motivated by the sociotechnical theory (Cooper & Foster, 1971; Mumford, 2006), this work urges that the harmony of interrelations between human workers and robots in a parts‐to‐picker order fulfillment center should be considered in realizing intelligent logistics systems. In particular, this work promotes the development of adaptive human‐friendly robotics to assist pickers in fulfilling orders with the fewest human–robot coordination conflicts. On the way toward “perfect automation,” the coexistence of robots and humans in task environments is inevitable, requiring carefully planned mechanisms for agile robot—human coordination that are humanity‐oriented and machine‐assisted.

Second, the proposed methodology is novel in conceptualizing the aforementioned philosophy of machine‐assisted humanity orientation into modeling and analysis to address the issue of robot–human coordination for parts‐to‐picker order fulfillment systems. Methodologically, this work integrates a discrete‐time nonlinear dynamic stochastic model with a time‐varying fatigue accumulation function to characterize and estimate a human picker's fatigue‐dependent working state. Then, this work develops a stochastic optimal control‐based data‐driven solution approach, which permits not only estimating human pickers’ instantaneous performance and robot queue lengths but also feeding back the aforementioned human picker's fatigue‐dependent working states in real time to the coordinated robots to adjust their efficiencies in handling and distributing mobile racks to the picker. Such a real‐time data‐driven stochastic optimal control‐based approach has its unique features and relative merits in addressing the critical yet under‐explored issue of real‐time robot–human coordination for machine‐assisted order fulfillment systems.

Third, this work is innovative in the realization of ideas from occupational psychology (psychology), organizational behavior (sociology), and robotics (engineering) to address picker–robot coordination for intelligent logistics in contemporary operations (management). Psychological factors are considered in the development of a novel robot–human‐coordinated mechanism for use in an automated parts‐to‐picker order fulfillment system. Through normative analysis and real‐data‐based empirical studies, the above concept (i.e., machine‐assisted humanity orientation) and proposed robot–picker‐coordinated order fulfillment mechanism are confirmed to be theoretically reasonable, scientifically solid and demonstrated to be of practical value in resolving the issue of interest.

PROBLEM DESCRIPTION

Typically, with our discussions with managers from a well‐established e‐commerce company in mainland China and referring to the problem settings of Wang et al. (2022), the workflow of mobile robot‐based parts‐to‐picker order fulfillment operations investigated in this paper contains the following steps. First, SKUs are received in the receiving area of a warehouse/distribution center; meanwhile, mobile robots that carry empty racks are dispatched to load the received SKUs. This is followed by the step of storage assignment, where mobile robots move loaded racks to the storage area (Li et al., 2020). If the SKUs in the racks are required in an order processing list, then the mobile robots move the racks with those required SKUs from the storage area to the picking–packaging area (termed as the step of rack handling and distribution). Then, stationary pickers in the picking–packaging area pick out the SKUs that are identified in the picking list and package them into cardboard boxes to fulfill customer orders (termed as the step of picking and packaging). This step is followed by automatically moving the completed cardboard boxes by means of conveyers to the outbound area for vehicular loading and dispatching. Mobile robots that have completed the delivery of the racks with the required SKUs to the stationary pickers return from the picking area to the storage area for repeated execution of the step of rack handling and distribution, to facilitate pickers’ carrying out picking and packaging for order fulfillment. Similar mobile robot‐based parts‐to‐picker order fulfillment operations can also be readily found in numerous practical cases, including Amazon.com (Demaitre, 2019; Gharehgozli & Zaerpour, 2020; Garland, 2022).

This work focuses on coordinating the step of rack handling and distribution (executed by robots) as well as the step of picking and packaging (executed by pickers). Specifically, the performance of mobile robots in handling and moving racks (number of racks handled and moved per robot per time unit) significantly influences the associated picker's performance in the subsequent picking and packaging tasks. The effect of psychological and behavioral uncertainties, which may be evident in the interaction between pickers and robots, on an associated picker's performance is considered herein. Accordingly, the problem of interest and goals of this work are as follows.

Consider a few‐to‐one parts‐to‐picker picking system, in which a few mobile robots work collaboratively with a human picker. A set of C customer orders, involving S SKUs, and being stored in N racks must be handled and moved by M mobile robots from the storage area to the picking area. Then, they have to be picked and packaged by one associated picker in a given working period T (e.g., 3 hours) without the intervention of scheduled breaks. A shared storage (scattered storage) policy, which is typically applied to mobile robot‐based parts‐to‐picker picking systems (Weidinger et al., 2018), is adopted. Different SKUs may be stored in a rack, and each rack can be associated with one or more customer orders. Two goals must be achieved. The first is to complete the picking and packaging of the aforementioned S SKUs that are stored and scattered over R racks by smoothly coordinating the operations of M mobile robots and the associated picker during a working period T. This goal must be achieved under the condition that the picker's fatigue increment (

$\mathrm{\Delta }{\rho }_k$

) in any time interval k (

$k=\{1,2,3\text{…}K\}$

;

$T=K·\overline{t}$

) during

$T$

should not exceed a predetermined fatigue degree (

$\overline{\rho }$

), such that the accumulated fatigue degree (

${\rho }_K$

) in the last time interval K does not exceed its upper bound (ρ_max), where

$\overline{t}$

is the length of a time interval. Meanwhile, in the few‐to‐one parts‐to‐picker picking system, the second goal is that the processed racks in the picking area are efficiently returned to the storage area by mobile robots for the next round of order processing.

Based on the above problem description and defined goals, this work proposes a real‐time robot–picker‐coordinated order fulfillment system that comprises the following four subsystems: (1) robot‐based rack handling (Subsystem 1), (2) mobile robot‐carried rack movement (Subsystem 2), (3) human‐based picking and packaging (Subsystem 3), and (4) robot‐based rack return (Subsystem 4), which are depicted in Figure 1. The number of newly arriving (

$A_k$

) racks that are loaded with SKUs from the receiving area to the storage area is treated as an exogenous variable, which is given in each time interval k (

$\forall k$

). Subsystems 2 and 4 are virtual as they are specified to determine the time‐varying distribution rates (

${\mu }_k$

and

${\nu }_k$

) of forward and reverse logistics for mobile racks that move between the storage (Subsystem 1) and picking–packaging (Subsystem 3) areas in any time interval k (

$\forall k$

). Subsystems 1 and 3 correspond to task operations that are executed by robots and associated pickers in the storage and picking–packaging areas, respectively. The outputs of Subsystems 1 and 3 are the number of racks (

$X_k$

) that are handled by mobile robots and the number of racks (

$Y_k$

) that are processed by a picker in each time interval k (

$\forall k$

). Both

$X_k$

and

$Y_k$

can be treated as metrics that are evaluated by contemporary detection technologies, such as image processing‐based detectors.

FIGURE 1

The proposed real‐time robot–picker‐coordinated order fulfillment system

The proposed robot–picker‐coordinated order fulfillment system has three distinctive features. First, the rack‐based order processing cycle, which integrates the forward and reverse logistical flows of mobile racks between the storage and picking–packaging areas, is considered. A closed‐loop rack reuse process facilitates the corresponding logistical resource management in the system. Second, the interactions between robots and associated pickers are considered, and they are characterized by the instantaneous inputs (

$Q_k$

and

$Q_{k|k}$

), which are defined as the perceived queue lengths of unprocessed racks at the beginning and end of time interval k, respectively, and output (

$Y_k$

, which is the number of racks that are processed by a picker in time interval k).

$Q_k$

and

$Q_{k|k}$

are given by:

1

where

$Q_{k-1|k-1}$

(

$Q_{k-1|k-1}=Q_{k-1}-Y_{k-1}$

) represents the number of unprocessed racks in the queue in the picking–packaging area (Subsystem 3) at the end of time interval

$k-1$

; and

${\stackrel{\sim }{X}}_k$

represents the number of unprocessed racks that arrive at the picking–packaging area (Subsystem 3) in time interval k. Third, the picker's psychological and physical conditions, characterized by the accumulated fatigue degree (

${\rho }_k$

) and its association with

$Y_k$

and

$Q_k$

, are considered and fed back to adjust the handling, transportation, and return decisions of the robots that are executed in Subsystems 1, 2, and 4, respectively. Hence, the picking and packaging tasks can be completed in a manner that includes critical human factors.

MODELING

This section proposes a discrete‐time nonlinear dynamic stochastic model to formulate and solve the above problem. Specifically, Section 4.1 defines the variables of system states (including state variables and control variables), and Section 4.2 explicates the proposed dynamic stochastic model.

Specification of system states

The state variables of Subsystem 3 are specified to define the expected performance of a picker in coordination with mobile robots. Then, the state variables for Subsystem 1 are specified to characterize the performance of the coordinated robots. Thereafter, the state variables of Subsystems 2 and 4 are specified to characterize the forward and reverse logistical flows of mobile racks that are required for coordination between Subsystems 1 and 3 (robots’ and human pickers’ actions).

States of human picker

Consider the states that are ideally associated with a picker's performance (Subsystem 3). Ideally, the picking and packaging tasks of a picker are carried out efficiently by seamlessly receiving and processing the loaded racks that are delivered by mobile robots. However, such an ideal state of Subsystem 3 is obtained at the cost of the picker's accumulated fatigue (

${\rho }_k$

). Worker fatigue is multidimensional and often increases with time (Glock, Grosse, Kim, et al., 2019; Jaber et al., 2013), stress, and workload (Do et al., 2018; Kc & Terwiesch, 2009; MacDonald, 2003), particularly when a picker must perform repetitive picking and packaging tasks. Based on the learning–forgetting–fatigue–recovery models established in the ergonomics‐related literature (Glock, Grosse, Kim, et al., 2019; Jaber et al., 2013), the commonly used deterministic exponential form of accumulated fatigue is extended herein into a dynamic stochastic form. Specifically, the fatigue growth rate (

${\lambda }_k$

) that determines the speed of fatigue accumulation is associated with the perceived queue length of unprocessed racks (

$Q_k$

) and work rate (

$Y_k$

). Therein,

$Q_k$

and

$Y_k$

can be treated, respectively, as the instantaneous mental and physical workloads that have been empirically verified to have an adverse impact on fatigue (Tan & Netessine, 2014). This model utilizes

$Q_k$

and

$Y_k$

to capture how both the instantaneous mental and physical workloads affect a picker's accumulated fatigue in any given time interval k (

$\forall k$

). Thus, the picker's accumulated fatigue (

${\rho }_k$

) in time interval k (

$\forall k$

) varies with time and can be expressed in a recursive form given by:

$${\rho }_k={\rho }_{k-1}+\mathrm{\Delta }{\rho }_k,\forall k,$$

2

where

${\rho }_{k-1}$

is the picker's accumulated fatigue in time interval

$k-1$

, and

$\mathrm{\Delta }{\rho }_k$

represents the fatigue increment in time interval k and equals

$1-e^{-{\lambda }_k}$

(

$\mathrm{\Delta }{\rho }_k=1-e^{-{\lambda }_k}$

). According to the empirical results of MacDonald (2003) from multiple regression analyses, both workloads (

$Q_k$

) and work rates (

$Y_k$

) have amplified impact on fatigue, inferring that

${\lambda }_k$

can be characterized by

$Q_k$

and

$Y_k$

in a simple linear manner as:

$${\lambda }_k={\beta }_0+{\beta }_1{Y}_k+{\beta }_2{Q}_k,\forall k,$$

3

where β₀ (

${\beta }_0\ge 0$

) is the intrinsic fatigue growth rate; β₁ (

${\beta }_1\ge 0$

) and β₂ (

${\beta }_2\ge 0$

) are the coefficients that capture the associations of

$Y_k$

and

$Q_k$

with

${\lambda }_k$

, respectively. The following condition on the picker's fatigue that accumulates in any time interval (Equation 4) must hold to ensure the wellbeing of the human worker):

$${\rho }_k\le {\rho }_{\max },\forall k.$$

4

The accumulated fatigue (

${\rho }_k$

) in time interval k may have a negative ergonomic effect on a picker's work rate (

$Y_{k+1}$

) in the next time interval

$k+1$

(Jaber et al., 2013; MacDonald, 2003): a higher

${\rho }_k$

often corresponds to a lower

$Y_{k+1}$

. Similarly, as shown by Kc and Terwiesch (2009), the overwork that is associated with accumulated fatigue has been demonstrated to affect a worker's performance (e.g., service rate). Do et al. (2018) assumed that the service rate has a multiplicative form under the effect of overwork. With reference to the fatigue‐recovery model of Jaber et al. (2013) and the above literature, in this work, a simple negative exponential function (

$f({\rho }_k)$

,

$f({\rho }_k)\equiv e^{-{\rho }_k}$

) is used to capture the moderating effect of

${\rho }_k$

, which can be expressed in a simple multiplicative form as shown in Equation (5):

$$Y_{k+1}=f\left({\rho }_k\right)Y_k=e^{-{\rho }_k}{Y}_k,\forall k.$$

5

Accordingly, Equation (5) can be used to elucidate how the accumulated fatigue in any given time interval k influences a picker's instantaneous picking and packaging performance in the next time interval

$k+1$

. Moreover, the system states and control variables of robots in response to a picker's fatigue‐dependent performance can also be determined.

In addition to the above concern about accumulated fatigue, the efficiency of picking and packaging tasks at each time interval is an important issue. Ideally, a picker processes all loaded racks as they arrive to the picking–packaging area from the storage area so that no racks remain in a queue in any time interval (

$Q_k=Y_k$

,

$\forall k$

). Let

${\gamma }_k$

be the picker‐efficiency state, which is defined as

${\gamma }_k\equiv \frac{Y_k}{Q_k}=\frac{Y_k}{Q_{k-1}+{\stackrel{\sim }{X}}_k}$

(Equation 1), and

$\overline{\gamma }$

(

$\overline{\gamma }=1,\forall k$

) be the ideal state of

${\gamma }_k$

. Then,

${\gamma }_k$

is expected to remain stable in each time interval (

${\gamma }_{k+1}={\gamma }_k$

,

$\forall k$

) under the ideal condition. If it does so, then

${\gamma }_{k+1}={\gamma }_k=\overline{\gamma }=1$

in all time intervals (

$\forall k$

), including the initial time interval (

$k=0$

). In practice,

${\gamma }_k$

may change over time and exhibit stochastic features under the influence of other system states (e.g., accumulated fatigue

${\rho }_k$

) and exogenous variables (e.g.,

$A_k$

,

$\forall k$

). According to Equation (5),

$Y_{k+1}$

, which is one component of

${\gamma }_{k+1}$

, can be influenced by

$f({\rho }_k)$

owing to the moderating effect of accumulated fatigue. Given

$Q_k$

and

${\stackrel{\sim }{X}}_{k+1}$

, which are exogenous variables to this subsystem (Subsystem 3),

${\gamma }_{k+1}$

can be associated with

${\gamma }_k$

in a simple deterministic form as

${\gamma }_{k+1}=e^{-{\rho }_k}{\gamma }_k$

to indicate the effect of accumulated fatigue on worker performance (Jaber et al., 2013). If the stochastic features of the state variables are further considered, the above deterministic form can be extended into a recursive form with a Gaussian white noise term (

$w_{{\gamma }_k}$

), as shown in Equation (6), to capture its dynamic and stochastic features:

6

where

$w_{{\gamma }_k}$

is a Gaussian white noise term of

${\gamma }_k$

. Then, the ideal state of

${\gamma }_k$

(

$\overline{\gamma }$

) is

$\overline{\gamma }=1$

(

$\forall k$

), which means that no racks are queuing in the picking–packaging area in any time interval.

States of mobile robots

In the ideal state, robots handle SKUs and racks in the storage area (Subsystem 1) in a way that maintains equal volumes of inbound and outbound rack flows of Subsystem 1 (

$X_k=A_k+{\stackrel{\sim }{Y}}_k$

,

$\forall k$

). Let

${\delta }_k$

be the robot‐efficiency rate, which is defined as

${\delta }_k\equiv \frac{X_k}{A_k+{\stackrel{\sim }{Y}}_k}$

, and

$\overline{\delta }$

(

$\overline{\delta }=1$

,

$\forall k$

) be the ideal state of

${\delta }_k$

. Then,

${\delta }_k$

should remain the same in all time intervals under the ideal condition, such that

${\delta }_{k+1}={\delta }_k=\overline{\delta }=1$

(

$\forall k$

). Despite the fact that fatigue accumulation does not apply to Subsystem 1, involving robots’ operations, a robot‐efficiency state (

${\delta }_k$

) can be influenced by the exogenous variables (e.g.,

$A_k$

,

$\forall k$

) and

${\stackrel{\sim }{Y}}_k,$

which is contingent upon the picker's performance. Like

${\gamma }_k$

,

${\delta }_k$

has dynamic and stochastic features. Thus,

${\delta }_k$

can be expressed in the recursive form with a Gaussian white noise term (

$w_{{\delta }_k}$

) as follows: ¹

7

Control variables for forward and reverse logistics flows of racks

Subsystems 2 and 4 are designed to regulate the forward and reverse logistical distribution flows of racks between the storage area and picking–packaging area. Let

${\mu }_k$

and

${\nu }_k$

be the rack distribution rates that are associated with Subsystems 2 and 4 (Figure 1) and defined as

${\mu }_k\equiv \frac{{\stackrel{\sim }{X}}_k}{X_k}$

and

$v_k\equiv \frac{{\stackrel{\sim }{Y}}_k}{Y_k}$

, respectively. Then,

${\mu }_k$

and

${\nu }_k$

can be used as control variables to regulate the numbers of unprocessed (

${\stackrel{\sim }{X}}_k$

) and processed (

${\stackrel{\sim }{Y}}_k$

) racks that arrive at Subsystems 3 and 1 in each given time interval k by setting

${\stackrel{\sim }{X}}_k=u_k{X}_k$

and

${\stackrel{\sim }{Y}}_k={\nu }_k{Y}_k$

, respectively, such that the system states (

${\gamma }_k$

and

${\delta }_k$

) that are associated with the human picker and robots become:

$${\gamma }_k=\frac{Y_k}{Q_{k-1}+{\mu }_k{X}_k},\forall k,$$

8

$${\delta }_k=\frac{X_k}{A_k+{\nu }_k{Y}_k},\forall k.$$

9

Based on the above defined efficiency‐related states (

${\gamma }_k$

and

${\delta }_k$

) and control variables (

${\mu }_k$

and

${\nu }_k$

), and given the number of newly arriving racks (

$A_k$

) and the measured previous queue length of racks (

$Q_{k-1}$

), the following analytical results (Theorem 4.1) regarding the outputs

${\widehat{X}}_k$

and

${\widehat{Y}}_k$

(estimates of instantaneous efficiency rates associated with coordinated robots and picker) from Subsystems 1 and 3, respectively, can be derived. Note that all technical proofs are provided in Online Appendix A. Theorem 4.1

Given

${\gamma }_k$

,

${\delta }_k$

,

${\mu }_k$

,

${\nu }_k$

,

$A_k$

,

$Q_{k-1}$

and

${\mathbf{\Omega }}_k\equiv {\delta }_k{\mu }_k{\nu }_k{\gamma }_k$

, the estimates of instantaneous projections of

$X_k$

and

$Y_k$

(denoted by

${\widehat{X}}_k$

and

${\widehat{Y}}_k$

) are obtained as:

$${\widehat{X}}_k=\frac{1}{1-{\mathbf{\Omega }}_k}\left({\delta }_k{A}_k+\frac{{\mathbf{\Omega }}_k{Q}_{k-1}}{{\mu }_k}\right),\forall k,$$

10

$${\widehat{Y}}_k=\frac{1}{1-{\mathbf{\Omega }}_k}\left(\frac{{\mathbf{\Omega }}_k{A}_k}{{\nu }_k}+{\gamma }_k{Q}_{k-1}\right),\forall k.$$

11

In Equations (10) and (11) of Theorem 4.1,

${\mathbf{\Omega }}_k$

is the product of the subsystem states (

${\gamma }_k$

,

${\delta }_k$

) and control variables (

${\mu }_k$

,

${\nu }_k$

);

$A_k$

can be treated as an exogenous variable of the system; and

$Q_{k-1}$

refers to queuing racks in the previous time interval

$k-1$

. Both

$A_k$

and

$Q_{k-1}$

are predetermined; however,

${\gamma }_k$

,

${\delta }_k$

,

${\mu }_k$

, and

${\nu }_k$

(

$\forall k$

) are unknown and must be estimated (as detailed in Section 4.2).

Based on Theorem 4.1, Corollary 4.1 is also used to characterize and estimate the time‐varying fatigue growth rate (

${\lambda }_k$

) and fatigue increment (

$\mathrm{\Delta }{\rho }_k$

) in any time interval k (

$\forall k$

). Corollary 4.1

Given Theorem 4.1 , estimates of

${\lambda }_k$

and

$\mathrm{\Delta }{\rho }_k$

(denoted by

${\widehat{\lambda }}_k$

and

$\mathrm{\Delta }{\widehat{\rho }}_k$

) are obtained as:

$${\widehat{\lambda }}_k={\beta }_0+({\beta }_1{\gamma }_k+{\beta }_2)\left(\frac{{\delta }_k{\mu }_k{A}_k+Q_{k-1}}{1-{\mathbf{\Omega }}_k}\right),\forall k,$$

12

$$\mathrm{\Delta }{\widehat{\rho }}_k=1-e^{-\left({\beta }_0+\left({\beta }_1{\gamma }_k+{\beta }_2\right)\left(\frac{{\delta }_k{\mu }_k{A}_k+Q_{k-1}}{1-{\mathbf{\Omega }}_k}\right)\right)},\forall k.$$

13

Based on Corollary 4.1, substituting Equation (13) into Equation (2) yields

$${\rho }_k={\rho }_{k-1}+1-e^{-\left({\beta }_0+({\beta }_1{\gamma }_k+{\beta }_2)\left(\frac{{\delta }_k{\mu }_k{A}_k+Q_{k-1}}{1-{\mathbf{\Omega }}_k}\right)\right)},\forall k.$$

14

Using the correlations between the system states and control variables specified above, the human picker and robots can be readily coordinated for superior performance in an order fulfillment center once the optimal values of the state and control variables can be determined. It is noted that the proposed model is data‐driven in which the input data required (including

$Y_k$

and

$A_k$

) are collected in each time interval k (

$\forall k$

) for estimating efficiency‐related states (

${\gamma }_k$

and

${\delta }_k$

) and control variables (

${\mu }_k$

and

${\nu }_k$

) such that the state variables (

${\gamma }_k$

and

${\delta }_k$

) can move toward the associated ideal values. Therein, the data used are collected in each time interval, which should be expressed in a discrete‐time manner rather than a continuous time manner for dynamically estimating and controlling state variables. Thus, the problem is formulated in a discrete‐time rather than a continuous‐time manner as presented in the next Section 4.2.

Model

This subsection proposes a discrete‐time nonlinear dynamic stochastic model that is characterized by recursive and measurement equations subject to boundary constraints (Anderson, Jr. et al., 2006; Santina et al., 1994; Sheu, 2002). Specifically, the time‐varying relationships among state variables, control variables, and measurement data defined in the previous section (Section 3), can be characterized in three forms: (1) recursive equations, (2) measurement equations, and (3) boundary constraints, as follows.

The recursive equations specify time‐varying relationships between the states of a dynamic and stochastic system in the next and current time intervals (

$k+1$

and k,

$\forall k$

) . Among such system states,

${\gamma }_k$

, and

${\delta }_k$

are treated as independent states, whereas

${\rho }_k$

is a state that depends on

${\gamma }_k$

and

${\delta }_k$

, according to Corollary 4.1 and Equation (14). In the proposed model, only the time‐varying relationships of independent states (

${\gamma }_k$

, and

${\delta }_k$

) must be formulated recursively; then, the dependent state

${\rho }_k$

can be readily determined. Based on Equations (6) and (7), the recursive equations that are associated with

${\gamma }_k$

and

${\delta }_k$

can be expressed in a generalized vector form as in Equation (15):

$${\mathbf{\Phi }}_{k+1}={\mathbf{F}}_k+{\mathbf{W}}_k,\forall k,$$

15

where

${\mathbf{\Phi }}_{k+1}$

,

${\mathbf{F}}_k$

and

${\mathbf{W}}_k$

are 2 × 1 time‐varying state vectors and can be further expressed as:

16

17

18

where

${\rho }_k={\rho }_{k-1}+1-e^{-({\beta }_0+({\beta }_1{\gamma }_k+{\beta }_2)(\frac{{\delta }_k{\mu }_k{A}_k+Q_{k-1}}{1-{\mathbf{\Omega }}_k}))}$

as in Equation (14) of Corollary 4.1.

Accordingly, the recursive equations (Equations 15 –18) reveal that the system states

${\gamma }_k$

and

${\delta }_k$

change over time, following the standard Gaussian–Markov processes and can be used to project

${\gamma }_{k+1}$

, and

${\delta }_{k+1}$

one‐time interval ahead, based on

${\gamma }_k$

and

${\delta }_k$

in any time interval k(

$\forall k$

).

The measurement equations characterize time‐varying relationships between measurements and system states. In this work,

$Y_k$

(the number of mobile racks that are processed by a picker in each given time interval k) is used as the measurement variable because it can be detected/measured readily using various detection technologies (such as image processing and automated counters). Based on Equation (11) in Theorem 4.1, we can further postulate the deterministic form of time‐varying relationships between the measurement (

$Y_k$

) and system states as

$Y_k=\frac{1}{1-{\mathbf{\Omega }}_k}(\frac{{\mathbf{\Omega }}_k{A}_k}{{\nu }_k}+{\gamma }_k{Q}_{k-1})$

. Owing to the measurement error in

$Y_k$

in any given time interval k (

$\forall k$

), a Gaussian white noise term

${\epsilon }_k$

is used to characterize this property. Since only one measurement variable (

$Y_k$

) is involved, one measurement equation is associated with

$Y_k$

in the proposed model. The generalized form of the measurement equation is given by:

$${\mathbf{Y}}_k={\mathbf{H}}_k+{\mathit{\epsilon }}_k,\forall k,$$

19

where

${\mathbf{Y}}_k$

,

${\mathbf{H}}_k$

and

${\mathit{\epsilon }}_k$

are 1 × 1 time‐varying state vectors, which can be further expressed as:

$${\mathbf{Y}}_k=Y_k,\forall k,$$

20

$${\mathbf{H}}_k=\left[\frac{1}{1-{\mathbf{\Omega }}_k}\left(\frac{{\mathbf{\Omega }}_k{A}_k}{{\nu }_k}+{\gamma }_k{Q}_{k-1}\right)\right],\forall k,$$

21

$${\mathit{\epsilon }}_k={\epsilon }_k,\forall k.$$

22

Additionally, the picker's accumulated fatigue (

${\rho }_k$

,

$\forall k$

) in any time interval should not exceed a preset upper bound (

${\rho }_{\max }$

), as indicated in Equation (4). This constraint is incorporated into the proposed model.

REAL‐TIME DATA‐DRIVEN STATE ESTIMATION

This section presents a stochastic optimal control‐based method for the real‐time estimation of system states (

${\gamma }_k$

and

${\delta }_k$

) and control variables (

${\mu }_k$

and

${\nu }_k$

), given information about the number of newly arriving racks (

$A_k$

) and the number of racks processed by the picker,

$Y_k$

, in each time interval (

$\forall k$

). Based on the principles of optimality that are applied in stochastic optimal control theory (Santina et al., 1994), the proposed algorithm searches for the optimal solutions for

${\gamma }_k$

,

${\delta }_k$

,

${\mu }_k$

and

${\nu }_k,$

which are updated using measurements of

$Y_k$

and the measurement equation (Equation 19) in each time interval k (

$\forall k$

) and then fed back as inputs into the recursive equation (Equation 15) such that the objective function ξ (Equation 23) is minimized:

23

where

${\widehat{\mathbf{\Phi }}}_k$

and

${\widehat{\mathit{U}}}_k$

are 2 × 1 vectors that contain the estimated system states (

${\gamma }_k$

and

${\delta }_k$

) and estimated control variables (

${\mu }_k$

and

${\nu }_k$

) in time interval k, respectively;

${\mathbf{\Phi }}_k^{\ast }$

and

${\mathit{U}}_k^{\ast }$

are 2 × 1 vectors that contain the ideal values of system states and control variables that are associated with

${\widehat{\mathbf{\Phi }}}_k$

and

${\widehat{\mathbf{U}}}_k$

, respectively; and

${\mathbf{\Psi }}_k^{\mathrm{\Phi }}$

and

${\mathbf{\Psi }}_k^U$

represent 2 × 2 time‐varying diagonal, positive‐definite weighting matrixes that are associated with

${\widehat{\mathbf{\Phi }}}_k$

and

${\widehat{\mathbf{U}}}_k$

, respectively. The objective function ξ in Equation (23) is the cost function, which is a scalar quadratic performance measure that indicates the deviation of the estimated system states and control variables from their ideal values in the order processing period.

To estimate the aforementioned system states and control variables in real time, a stochastic optimal control‐based method is developed using an extended Kalman filter. Kalman filtering is a well‐known statistical method for the linear‐quadratic estimation of system states of dynamic and stochastic systems. It has been successfully and extensively utilized in solving diverse estimation and control problems in various fields, such as the tracking and navigation of different sorts of vehicles and traffic signal control (Sheu, 2002). Whereas the basic Kalman filter applies to linear problems, the extended Kalman filter has been developed particularly for nonlinear problems that are characterized by either recursive or measurement equations in dynamic and stochastic models (Santina et al., 1994). The primary computational steps in the proposed estimation method are summarized below.

Initialization

This step initializes system states and all of the inputs that are required to trigger the subsequent computational steps. Let

$k=0$

and

$\mathrm{\Delta }{\mathbf{\Phi }}_k\equiv {\mathbf{\Phi }}_k-{\widehat{\mathbf{\Phi }}}_k$

(

$\forall k$

). Then, preset the initial estimates of the 2 × 1 state vector (

${\widehat{\mathbf{\Phi }}}_0$

), the 2 × 1 control variable vector (

${\widehat{\mathit{U}}}_0$

), the number of queued racks (Q ₀), and the 2 × 2 covariance matrix of the state estimation error (

${\widehat{\mathit{P}}}_{0|0}$

), where

$Q_0=0$

;

${\widehat{\mathbf{\Phi }}}_0$

,

${\widehat{\mathit{U}}}_0$

and

${\widehat{\mathit{P}}}_{0|0}$

are given by:

24

$${\widehat{\mathit{P}}}_{0|0}\equiv E\left[\mathrm{\Delta }{\mathbf{\Phi }}_0\mathrm{\Delta }{\mathbf{\Phi }}_0^{\mathrm{T}}\right]={\mathbf{P}}_0,(\mathrm{where}k=0).$$

25

Prior prediction of system states

In this step, the state vector

${\mathbf{\Phi }}_{k+1}$

and covariance matrix of the state estimation error

${\mathbf{P}}_{k+1}$

in time interval

$k+1$

are predicted using the estimates of system states and measurements made in time interval k (

$\forall k$

). The prior prediction of

${\mathbf{\Phi }}_{k+1}$

(

${\widehat{\mathbf{\Phi }}}_{k+1|k}$

) will be updated for the linear minimum mean square (LMMS) estimation of

${\mathbf{\Phi }}_{k+1}$

using the extended Kalman filter in the next stage; and the prior prediction of

${\mathbf{P}}_{k+1}$

(

${\widehat{\mathit{P}}}_{k+1|k}$

) is a prerequisite for the aforementioned state estimation. Using Equation (15), Corollaries 5.1 and 5.2 are provided to facilitate the prior predictions of

${\widehat{\mathbf{\Phi }}}_{k+1|k}$

and

${\widehat{\mathit{P}}}_{k+1|k}$

, respectively. Corollary 5.1

Given Equation ( 15 ) (recursive equation) and its characteristics, let

, where

${\widehat{\gamma }}_k$

and

${\widehat{\delta }}_k$

are the LMMS estimates of

${\gamma }_k$

and

${\delta }_k$

that are based on the data regarding

$A_k$

and

$Y_k$

collected in time interval k . Then, the prior prediction of the LMMS estimate of

${\mathbf{\Phi }}_{k+1}$

(

${\widehat{\mathbf{\Phi }}}_{k+1|k}$

) that is based on the data regarding

$A_k$

and

$Y_k$

that were collected in time interval k is given by:

$${\widehat{\mathbf{\Phi }}}_{k+1|k}={\widehat{\mathbf{F}}}_k,\forall k.$$

26

The prior prediction (

${\widehat{\mathit{P}}}_{k+1|k}$

) of the covariance matrix of the state estimation error

${\mathbf{P}}_{k+1}$

as well as

${\widehat{\mathbf{\Phi }}}_{k+1|k}$

must be determined in this stage. Using Equations (15) and (26) (Corollary 5.1), Corollary 5.2 is obtained to determine

${\widehat{\mathit{P}}}_{k+1|k}$

. Corollary 5.2

Given Equation ( 15 ) (the recursive equations) and Corollary 5.1 , let

${{\widehat{\mathit{F}}}^{\prime }}_k={\frac{\partial {\mathbf{F}}_k}{\partial {\mathbf{\Phi }}_k}|}_{{\mathbf{\Phi }}_k={\widehat{\mathbf{\Phi }}}_k}$

based on the data regarding

$A_k$

and

$Y_k$

that were collected in time interval k, and let

${\mathbf{R}}_k\equiv E[{\mathbf{W}}_k{\mathbf{W}}_k^{\mathrm{T}}]$

. Then, the prior prediction of

${\mathbf{P}}_{k+1}$

(

${\widehat{\mathit{P}}}_{k+1|k}$

) that is based on the data regarding

$A_k$

and

$Y_k$

that were collected in time interval k is given by:

$${\widehat{\mathit{P}}}_{k+1|k}={{\widehat{\mathit{F}}}^{\prime }}_k{\widehat{\mathit{P}}}_{k|k}{{\widehat{\mathit{F}}}^{\prime \mathrm{T}}}_k+{\mathbf{R}}_k,\forall k.$$

27

Correction of prior predictions

In this stage, the prior predictions

${\widehat{\mathbf{\Phi }}}_{k+1|k}$

and

${\widehat{\mathit{P}}}_{k+1|k}$

are corrected using the Kalman gain (

${\mathbf{G}}_{k+1}$

) and the measurements made through time interval

$k+1$

(

$\forall k$

). To achieve this purpose, Corollary 5.3, which associates the corrected state vector (

${\widehat{\mathbf{\Phi }}}_{k+1}$

) with its prior prediction (

${\widehat{\mathbf{\Phi }}}_{k+1|k}$

), the Kalman gain (

${\mathbf{G}}_{k+1}$

) and the measurements made through time interval

$k+1$

(

$\forall k$

), is provided. Corollary 5.3

Given Equations ( 5 ) (fatigue prediction function) and ( 19 ) (measurement equation), let

$\mathrm{\Delta }{\mathbf{y}}_{k+1|k}$

be the measurement residual in time interval

$k+1$

, defined as

$\mathrm{\Delta }{\mathbf{y}}_{k+1|k}\equiv {\mathbf{Y}}_{k+1}-{\widehat{\mathit{Y}}}_{k+1|k}={\mathbf{Y}}_{k+1}-e^{-{\widehat{\rho }}_k}{\widehat{\mathit{H}}}_k$

(

$\forall k$

), where

${\widehat{\mathit{H}}}_k=[\frac{1}{1-{\widehat{\mathbf{\Omega }}}_k}(\frac{{\widehat{\mathbf{\Omega }}}_k{A}_k}{{\widehat{\nu }}_k}+{\widehat{\gamma }}_k{Q}_{k-1})]$

(by Equation 21 ), and let

$\mathrm{\Delta }{\mathbf{Y}}_k$

be a

$k\times 1$

residual vector of the measurement residuals through time interval k that is defined as

. Then, the LMMS estimate of

${\mathbf{\Phi }}_{k+1}$

(

${\widehat{\mathbf{\Phi }}}_{k+1}$

), which is based on the data that were collected through time interval

$k+1$

is given by:

$${\widehat{\mathbf{\Phi }}}_{k+1}={\widehat{\mathbf{\Phi }}}_{k+1|k}+{\mathbf{G}}_{k+1}·\mathrm{\Delta }{\mathbf{y}}_{k+1|k}+{\mathbf{m}}_{\mathbf{\Phi }},\forall k,$$

28

where m _Φ is the mean of

${\mathbf{\Phi }}_{k+1}$

, and

${\mathit{G}}_{k+1}$

is the Kalman gain, which is defined as:

$${\mathbf{G}}_{k+1}\equiv E\left[{\mathbf{\Phi }}_{k+1}\mathrm{\Delta }{\mathbf{y}}_{k+1|k}^{\mathrm{T}}\right]·{\left\{E[\mathrm{\Delta }{\mathbf{y}}_{k+1|k}\mathrm{\Delta }{\mathbf{y}}_{k+1|k}^{\mathrm{T}}]\right\}}^{-1},\forall k.$$

29

As seen in Equations (28) and (29) of Corollary 5.3, the LMMS estimate of

${\mathbf{\Phi }}_{k+1}$

(

${\widehat{\mathbf{\Phi }}}_{k+1}$

) given the measurements made through time interval

$k+1$

can be derived following the determination of the Kalman gain

${\mathbf{G}}_{k+1}$

as in Equation (29). Therefore, Corollary 5.4 is provided to determine the Kalman gain

${\mathbf{G}}_{k+1}$

. Corollary 5.4

Given Corollaries 5.2 and 5.3, and Equation ( 19 ) (measurement equation), let

${{\widehat{\mathit{H}}}^{\prime }}_{k+1}={\frac{\partial {\mathbf{H}}_{k+1}}{\partial {\mathbf{\Phi }}_{k+1}}|}_{{\mathbf{\Phi }}_{k+1}={\widehat{\mathbf{\Phi }}}_{k+1|k}}$

and

${\mathbf{\Lambda }}_{k+1}\equiv E[{\mathit{\epsilon }}_k{\mathit{\epsilon }}_k^{\mathrm{T}}]$

. Then, the Kalman gain

${\mathbf{G}}_{k+1}$

is given by:

$$\begin{array}{c}\hfill {\mathbf{G}}_{k+1}={\widehat{\mathit{P}}}_{k+1|k}{{\widehat{\mathit{H}}}^{\prime \mathrm{T}}}_{k+1}{\left({{\widehat{\mathit{H}}}^{\prime }}_{k+1}{\widehat{\mathit{P}}}_{k+1|k}{{\widehat{\mathit{H}}}^{\prime \mathrm{T}}}_{k+1}+{\mathbf{\Lambda }}_{k+1}\right)}^{-1},\forall k.\end{array}$$

30

Based on Corollaries 5.3 and 5.4, the corrected state vector

${\widehat{\mathbf{\Phi }}}_{k+1}$

can be obtained using Equations (28) and (30). Alternatively,

${\widehat{\mathbf{\Phi }}}_{k+1}$

can also be expressed as:

31

Finally, the recursive equation for correcting the prior prediction of the covariance matrix of state estimation error

${\widehat{\mathit{P}}}_{k+1|k}$

into

${\widehat{\mathit{P}}}_{k+1|k+1}$

must be derived for the recursive estimation of system states in the next time interval

$k+1$

. To achieve this purpose, Corollary 5.5 is provided. Corollary 5.5

Given Corollaries 5.2 to 5.4, let the corrected covariance matrix of state estimation error

${\widehat{\mathit{P}}}_{k+1|k+1}$

be defined as

${\widehat{\mathit{P}}}_{k+1|k+1}\equiv E[\mathrm{\Delta }{\mathbf{\Phi }}_{k+1|k+1}\mathrm{\Delta }{\mathbf{\Phi }}_{k+1|k+1}^{\mathrm{T}}]$

, where

$\mathrm{\Delta }{\mathbf{\Phi }}_{k+1|k+1}\equiv {\mathbf{\Phi }}_{k+1}-{\widehat{\mathbf{\Phi }}}_{k+1}$

. Then,

${\widehat{\mathit{P}}}_{k+1|k+1}$

is given by:

$${\widehat{\mathit{P}}}_{k+1|k+1}=(\mathbf{I}-{\mathbf{G}}_{k+1}{{\widehat{\mathit{H}}}^{\prime }}_{k+1}){\widehat{\mathit{P}}}_{k+1|k}+{\mathbf{m}}_{\mathbf{\Phi }}{\mathbf{m}}_{\mathbf{\Phi }}^{\mathrm{T}},\forall k.$$

32

The result of Corollary 5.5 in (32) is important for the estimation of the control variable vector as shown in the next subsection.

Estimation of control variable vector

Using the corrected state vector

${\widehat{\mathbf{\Phi }}}_{k+1}$

(Equation 31) and the covariance matrix of the state estimation error

${\widehat{\mathit{P}}}_{k+1|k+1}$

(Equation 32), the control variable vector

${\widehat{\mathit{U}}}_{k+1}$

(

,

$\forall k$

) can be estimated. The principles of stochastic optimal control theory and existing algorithms (Santina et al., 1994; Sheu, 2002) are applied here to estimate

${\widehat{\mathit{U}}}_{k+1}$

. Consistent with the principles of stochastic optimal control, the corrected state vector

${\widehat{\mathbf{\Phi }}}_{k+1}$

is fed back through the optimal control gain matrix

${\mathbf{\Theta }}_{k+1}$

to determine the control variable vector

${\widehat{\mathit{U}}}_{k+1}$

:

$${\widehat{\mathit{U}}}_{k+1}=-{\mathbf{\Theta }}_{k+1}{\widehat{\mathbf{\Phi }}}_{k+1}+{\mathit{\eta }}_{k+1},\forall k.$$

33

In Equation (33),

${\mathbf{\Theta }}_{k+1}$

and

${\mathit{\eta }}_{k+1}$

are given by:

$$\begin{array}{c}\hfill {\mathbf{\Theta }}_{k+1}={\left[{\mathbf{B}}_{k+1}^{\mathrm{T}}{\mathbf{\Gamma }}_{k+2}{\mathbf{B}}_{k+1}+{\mathbf{\Psi }}_{k+1}^{\mathrm{U}}\right]}^{-1}{\mathbf{B}}_{k+1}^{\mathrm{T}}{\mathbf{\Gamma }}_{k+2}{\widehat{{\mathit{F}}^{\prime }}}_{k+1},\forall k,\end{array}$$

34

35

where

${\mathbf{B}}_{k+1}={\frac{\partial {\mathbf{F}}_{k+1}}{\partial {\mathbf{U}}_{k+1}}|}_{{\mathbf{U}}_{k+1}={\widehat{\mathit{U}}}_{k|k}}$

(

$\forall k$

), and matrix

${\mathbf{\Gamma }}_{k+2}$

should satisfy the Riccati equation as follows:

$$\begin{array}{c}\hfill {\mathbf{\Gamma }}_{k+1}={\mathbf{\Psi }}_{k+1}^{\mathrm{\Phi }}+{{\mathbf{F}}^{\prime }}_{k+1}^{\mathrm{T}}{\mathbf{\Gamma }}_{k+2}{{\mathbf{F}}^{\prime }}_{k+1}-{{\mathbf{F}}^{\prime }}_{k+1}^{\mathrm{T}}{\mathbf{\Gamma }}_{k+2}{\mathbf{B}}_{k+1}{\mathbf{\Theta }}_{k+1},\forall k.\end{array}$$

36

All of the system states that are characterized by the corrected state vector (

${\widehat{\mathbf{\Phi }}}_{k+1}$

) and the control variable vector (

${\widehat{\mathit{U}}}_{k+1}$

) in time interval

$k+1$

can be estimated using the aforementioned approach and the data that are collected through time interval

$k+1$

.

Utilizing the estimates of

${\widehat{\mathbf{\Phi }}}_{k+1}$

and

${\widehat{\mathbf{U}}}_{k+1}$

, the state‐dependent variables, including

${\widehat{X}}_{k+1}$

and

${\widehat{Y}}_{k+1}$

(by Theorem 4.1), together with

$Q_{k+1}$

(by Equation 1), can be updated and used as indices of the instantaneous performance of the system in terms of the efficiencies of the coordinated robots and the picker. Moreover, a picker's fatigue degree (

${\widehat{\rho }}_{k+1}$

) accumulated to time interval

$k+1$

can be estimated (by Corollary 4.1 and Equation 14) in real time to determine whether it satisfies the boundary constraint (Equation 4) that protects the safety of the workers.

The proposed data‐driven real‐time state estimation‐control approach can be coded in any appropriate computer programming language. This can not only facilitate the recursive calculation and optimal control of system states but also monitor system performance in real time. Figure 2 summarizes the primary algorithmic computational procedures based on the proposed approach.

FIGURE 2

Algorithmic procedures for data‐driven real‐time state estimation‐control

REAL‐PRACTICE‐BASED ANALYSIS

To demonstrate the validity and applicability of the proposed real‐time data driven state estimation and control approach for robot–picker‐coordinated order fulfillment, data were collected from a dominant e‐commerce company in mainland China that uses robots for parts‐to‐picker order fulfillment. The data were obtained from a warehouse of that company, which is a 2500m² robot‐carried mobile‐rack warehouse that stores small cosmetic items, in which 11 human pickers work with mobile robots to pick out customer orders. The raw data provide details about the pickers’ picking and packaging tasks at the study site over 6 days in March 2019. Data concerning the number of racks (

$A_k$

) that carry newly arriving SKUs to the storage area and the number of racks (

$Y_k$

) that are processed by each picker in the picking–packaging area at each time interval are used.

A total of 21 datasets that contain a total of 756 data points are generated from the collected raw data to facilitate real‐time state estimation and control in this empirical study. Each dataset contains 18 pairs of data points

$A_k$

and

$Y_k$

, measured every 10 min in a 3‐hour (3‐h) task period. The data points

$A_k$

are inputs to and

$Y_k$

are the measurements for the proposed approach. The generated 21 datasets are classified into three groups, associated with peak‐hour, normal, and off‐peak operational scenarios to demonstrate the applicability of the proposed approach to diverse operational scenarios in order fulfillment centers. Specifically, seven, eight, and six datasets are obtained for Scenarios I (peak hour), II (normal), and III (off‐peak), respectively. Using the 21 datasets, the robot–picker‐coordinated order fulfillment performance and accumulated fatigue of the pickers in peak‐hour, normal, and off‐peak operational scenarios (Scenarios I, II, and III) are comparatively analyzed.

Given the aforementioned 21 datasets (756 data points), the real‐practice‐based empirical analysis is carried out through two phases to demonstrate the relative performance (RP) of the proposed robot–picker coordination mechanism for parts‐to‐picker order fulfillment.

The first phase of the empirical analysis aims to compare the results generated using the proposed approach with the measurements gained from the 21 datasets in the aspects of picker efficiency (the average of

$Y_k$

values) and accumulated fatigue (the average of ρ _K values). The comparative results are presented in Figure 3, which help assess the RP of the proposed approach against the current parts‐to‐picker order fulfillment operations of the focal firm in the aspects of picker efficiency and fatigue alleviation, respectively.

FIGURE 3

Empirical results of phase I (robot–picker coordination relative to real operations). Scenarios: I, II, and III represent “peak,” “normal,” and “off‐peak” scenarios; Indexes: a and b represent the scenarios “with” and “without” considering robot–picker coordination; ARP, average of relative performance; RP, relative performance

Overall, the empirical results yielded in the first phase provide the following findings as summarized in Observation 6.1. Observation 6.1

The proposed robot–picker coordination system permits alleviating a picker's fatigue in order fulfillment without much influence on picking efficiency, compared with the current performance of the focal firm in the study case—As can be seen in Figure 3, a picker's accumulated fatigue can be reduced by 53.74% if the proposed robot–picker coordination system is applied in the focal firm of the study case. Such a fatigue alleviation, however, is carried out at the expense of lowering picker efficiency by 14.79%. Nevertheless, the above observation is encouraging for the applicability of the proposed method in practical cases of parts‐to‐picker order fulfillment systems as it implies that a more harmonious and enjoyable robot–human co‐working environment can be created to facilitate the development of intelligent and human‐friendly robot–human‐coordinated order fulfillment systems.

In the second phase of the empirical analysis, this work aims to demonstrate the RP of robot–picker coordination in contrast with the case without considering the robot–picker coordination for gaining more managerial insights. Therein, empirical results that are obtained with and without consideration of the robot–picker coordination mechanism are compared. Without robot–picker coordination, both control variables

${\mu }_k$

and

${\nu }_k$

are set to 1 (

${\mu }_k={\nu }_k=1$

) to mimic the current parts‐to‐picker order fulfillment operations of the focal firm of the case studied.

Figure 4 plots the empirical results concerning practice‐based empirical analysis, including the average values of scenario‐based state variables, control variables, system performance, and associated RP, as determined by comparing the outputs with and without robot–picker coordination. The findings and managerial implications of Figure 4 are discussed below.

FIGURE 4

Empirical results of phase II (robot–picker coordination performance). Scenarios: I, II, and III represent “peak,” “normal,” and “off‐peak” scenarios; indexes: a and b represent the scenarios “with” and “without” considering robot–picker coordination; RP, relative performance

Observation 6.2

The proposed robot–picker coordination system outperforms the system without coordination not only in efficiency but also the fatigue of human pickers—As seen in Figure 4, the RP values that are associated with

${\gamma }_k$

in all three scenarios (peak‐hour, normal, and off‐peak scenarios) are positive, indicating that the pickers’ working performance with robot–picker coordination is more efficient than without coordination. This higher efficiency is attributed to shorter queue racks (

$Q_k$

), owing to adjustments of the number of mobile racks that arrive at the picking and packaging area in each time interval based on estimates of the control variable μ _k in real time using the proposed approach. Additionally, the values of accumulated fatigue ( ρ _k ) that are found in all three scenarios with robot–picker coordination are lower than the scenarios without coordination. This observation is highly encouraging as it implies that customer orders can be efficiently fulfilled without the need for pickers to expend more energy and make more efforts while operating under the proposed robot–picker coordination, relative to the situation without it.

Observation 6.3

Robot–picker coordination improves system performance more in the “normal” scenario (Scenario II) than in the “peak‐hour” and “off‐peak” scenarios (Scenarios I and III)—This inference is drawn by comparing the RP values, particularly those associated with

$Q_k$

and

${\rho }_k$

in the three scenarios. The “minus” signs of the RP values that are associated with

$Q_k$

and

${\rho }_k$

imply that the proposed robot–picker coordination system alleviates pickers’ accumulated fatigue (

${\rho }_k\downarrow$

) by reducing the number of racks that queue (

$Q_k\downarrow$

) in the picking–packaging area. Such an effect of robot–picker coordination on either

$Q_k$

or

${\rho }_k$

is highly significant under normal conditions (Scenario II) and is reasonable as the length of queues of racks (

$Q_k$

) a key factor that is related to the stress and workloads of pickers as captured by Equation (3) in the proposed model.

Observation 6.4

(Human‐friendly robot behavior for coordination with human pickers in different operational scenarios). Under the proposed robot–picker coordination scheme, robots can exhibit different coordination behaviors in response to the instantaneous performance of coordinated pickers in different scenarios, as follows.

Under peak‐hour conditions (Scenario I)—Robots tend to handle racks with high efficiency (by increasing

$X_k$

value) in the storage area, and this efficiency almost matches the coordinated picker's efficiency in processing racks (

$X_k\approx Y_k$

). Moreover, robots tend to stabilize the volume of racks that are delivered to the coordinated pickers (by stabilizing

${\mu }_k$

value); meanwhile, the rack volume that is returned to the storage area is adjusted (by reducing

${\nu }_k$

value) based on the coordinated pickers’ instantaneous efficiency (by estimating

$Y_k$

value). Consequently, the accumulated number of queuing racks in the picking—packaging area can be well controlled (by reducing the

$Q_k$

value) to mitigate the accumulated stress and fatigue of coordinated human pickers during peak hours.

Under normal conditions (Scenario II)—Robots adjust their efficiency of handling racks in the storage area to keep the number of racks that are moved to the picking–packaging area almost equal to the number of racks that are processed by the coordinated pickers (

$X_k\approx Y_k$

) per unit time such that the accumulated number of racks that queue in the picking–packaging area remains stable. (

$Q_k$

remains approximately double the value of

$Y_k$

in this case). Meanwhile, the volumes of mobile racks (carried by robots) that move between the storage and picking–packaging areas are controlled to be almost equal to each other (

${\mu }_k\approx {\nu }_k$

).

Under off‐peak conditions (Scenario III)—Robots handle racks in the storage area more slowly (by reducing

$X_k$

value) to keep their efficiency a little lower than the coordinated pickers’ efficiency in processing racks (

$X_k\le Y_k$

). Meanwhile, the volumes of racks that are moved between the storage area and the picking–packaging area are adjusted in response to the coordinated picker's processing of racks. Specifically, robots tend to move racks forward fast (by increasing

${\mu }_k$

value) and to return racks from the picking–packaging area more slowly (by reducing

${\nu }_k$

value) during off‐peak hours.

Three empirical examples, based on three selected datasets, are further used to present graphically the aforementioned characteristics of robot–picker coordination behavior that is typically exhibited during “peak,” “normal,” and “off‐peak” hours as shown in Figures 5 to 7. Briefly, the coordinated robots adjust their efficiency (

$X_k$

) in handling racks in the storage area and controlling the volumes of racks (by

${\mu }_k$

and

${\nu }_k$

) that move between the storage area and the picking–packaging area in response to the coordinated pickers’ efficiency in processing racks (by

$Y_k$

) and accumulated fatigue (

${\rho }_k$

). Under such a proposed robot–picker coordination scheme, the lengths of the queues of racks in the picking–packaging area can also be controlled.

FIGURE 5

Empirical example (1)—system output during peak hours (Scenario I)

FIGURE 6

Empirical example (2)—system output during normal hours (Scenario II)

FIGURE 7

Empirical example (3)—system output during off‐peak hours (Scenario III)

The results of empirical example 1 (Figure 5) reveal a large difference between accumulated fatigue (

${\rho }_k$

) in a 3‐h peak period in Scenario I and those shown in the other two scenarios. Thus, Scenario I is used to conduct sensitivity analysis on two key parameters (β₁, and β₂), which determine the effects of picker efficiency (

$Y_k$

) and rack queue length (

$Q_k$

), respectively, on the fatigue growth rate (

${\lambda }_k$

) as captured by Equation (3). The parameter settings (

${\beta }_1=0.00005$

and

${\beta }_2=0.00005$

) and dataset that are used in empirical example 1 provide the baseline in the sensitivity analysis. Figure 8 plots the empirical results of the sensitivity analysis of the relationships among

${\rho }_k$

, β₁, and β₂, where the value of ρ is the average of the

${\rho }_k$

values that are estimated during the 3‐h peak period. Figure 8 has managerial implications concerning the psychological and physical effects of workload on human pickers, which are determined by the perceived rack queue lengths and the picking efficiency of pickers in coordination with robots during peak hours.

FIGURE 8

Sensitivity analysis of ρ (with respect to β₁ and β₂) in peak‐hour scenario

Observation 6.5

A picker's accumulated fatigue (ρ) increases significantly with either the picker's efficiency (

$Y_k$

) or perceived rack queue length (

$Q_k$

) during peak hours; the increase in ρ with

$Q_k$

exceeds that with

$Y_k$

.

During peak hours, a human picker is likely to be more sensitive to the perceived length of the queue of unprocessed mobile racks than to the number of racks processed. Therefore, perceiving any anomalous increase in the length of the queue of mobile racks would easily have a negative effect on the picker's psychological and behavioral responses (such as stress, anxiety, and panic).

CONCLUSIONS, MANAGERIAL IMPLICATIONS, AND RECOMMENDATIONS

Conclusion

Motivated by the introduction of robots and advanced technologies into intelligent logistics operations and the impact thereof on worker safety and wellbeing, this work has presented a real‐time data‐driven stochastic optimal control‐based approach to address the challenges of robot–picker‐coordinated order fulfillment. Specifically, the proposed approach combines a discrete‐time nonlinear stochastic‐dynamic model with a fatigue accumulation function (Section 4) into a real‐time system state estimation and control method (Section 5), to find the optimal solutions for the coordination between human pickers and robots. Theoretically and in practice, the coordinated robots can react and adapt to picker efficiency in real time, such that the coordinated robots and pickers can jointly fulfill orders in collaboration and in a human‐friendly manner. Sections 4 and 5 have presented the methodological characteristics and novelty of the proposed approach, which help address the first research question raised in Section 1. To answer the second question (Section 1), an empirical analysis has been performed using actual data that were collected from a dominant e‐commerce company in China that uses robots for “parts‐to‐picker” order fulfillment. Our empirical study has been conducted to generate various insights (called “observations”), which help answer the second question, concerning how the proposed human‐friendly, robot–human coordination system improves the performance of a parts‐to‐picker order fulfillment center in terms of both safety and human wellbeing.

Managerial implications

The analytical results that have been obtained using the proposed model (Sections 4 and 5) and the empirical results of the empirical analysis of actual data (Section 6) provide the following managerial implications.

Implication 1. Under the proposed robot–picker‐coordinated order fulfillment mechanism, the robot and picker efficiencies (

$X_k$

and

$Y_k$

) have the following characteristics.

• The robot and picker efficiencies (

  $X_k$
  and
  $Y_k$
  ) can be effectively controlled and coordinated if real‐time data concerning newly arriving racks (
  $A_k$
  ) that are loaded with SKUs are obtained.

• Robot efficiency (

  $X_k$
  ) depends strongly on
  $A_k$
  , whereas picker efficiency (
  $Y_k$
  ) depends strongly on the length of the queue of unprocessed racks (
  $Q_k$
  ).

The above managerial insights are derived from Theorem 4.1 and Observations 6.2 and 6.5. The effect of the perceived length of the queue (

$Q_k$

) of unprocessed racks on systems performance is emphasized and has implications for operations managers (e.g., in warehousing), enabling them better to manage the number of racks in queues by appropriately controlling the speeds with which robots deliver racks to the picking–packaging area.

Implication 2. Under the proposed robot–picker‐coordinated order fulfillment mechanism, the picker‐to‐robot RP in the aspect of either efficiency (

$Y_k/\phantom{Y_k{X}_k}{X}_k$

) or effectiveness (

${\gamma }_k/\phantom{{\gamma }_k{\delta }_k}{\delta }_k$

) in peak hours is less than that in either normal or off‐peak hours. The picker‐to‐robot RP is greatest in the “off‐peak” scenario. This finding is drawn from a comparison of empirical results (Figure 4) for

$Y_k/\phantom{Y_k{X}_k}{X}_k$

and

${\gamma }_k/\phantom{{\gamma }_k{\delta }_k}{\delta }_k$

among the three scenarios. It may help operations managers to determine appropriate criteria and standards for evaluating the efficiency and effectiveness of pickers in different scenarios.

Implication 3. Robot efficiency is not the most important metric. Rather, human efficiency dominates the efficiency of a parts‐to‐picker order fulfillment system, particularly during peak hours. Accordingly, controlling the number of queued racks by slowing down their delivery by robots to the picking–packaging area is more efficient than slowing down the rack‐processing of a picker to alleviate that picker's fatigue, particularly during a peak‐hour period. This implication is drawn from the empirical results (Figure 8) and is consistent with Observation 6.5 and Implication 1(b). It suggests the importance of accounting for human wellbeing and safety in the design and implementation of a human‐friendly intelligent robot–human‐coordinated order fulfillment system. In such a system, humans are the key element, and robots should be sufficiently intelligent and accommodating to assist humans in the joint intelligent logistics task of order fulfillment.

Recommendations for future research

Despite the novelty of the proposed methodology and its advantages in realizing a human‐friendly intelligent robot–human‐coordinated order fulfillment system, this research has several limitations. First, the proposed model has not incorporated factors related to worker psychology (e.g., perceived job worth, job attitude, and emotional involvement) and behavior (e.g., body posture) toward coordination and cooperation with robots, and these factors should be further explored and empirically verified in the future (MacDonald, 2003; Tan & Netessine, 2014; Tan et al., 2022). For example, body posture is regarded as a kind of physical demand, which may contribute to psychosocial hazards (i.e., negative effects on employees such as stress and health problems), as claimed in MacDonald (2003). Thus, body posture can be another issue to improve employees’ stress and health problems in human–robot collaborative work settings. Other related tasks that are undertaken in inbound logistics (such as storage assignment) and outbound logistics (such as resource/fleet management, vehicular loading, dispatching, and routing) have not yet been considered. They may deserve further investigation. Extending the proposed model for the cases of multiple working periods intervened by scheduled breaks is analytically challenging yet meaningful for future research. Last but not least, considering how robotics can be used to help cope with capacity bottlenecks and supply chain disruptions due to manpower shortage under/after the Coronavirus disease (COVID‐19) pandemic (Gupta et al., 2022) will also be interesting to study in the future.