- Email: [email protected]
Order release planning by iterative simulation and linear programming: Theoretical foundation and analysis of its shortcomings
Accepted Manuscript
Order release planning by iterative simulation and linear programming: Theoretical foundation and analysis of its shortcomings Hubert Missbauer PII: DOI: Reference:
S0377-2217(19)30594-6 https://doi.org/10.1016/j.ejor.2019.07.030 EOR 15931
To appear in:
European Journal of Operational Research
Received date: Accepted date:
16 July 2018 10 July 2019
Please cite this article as: Hubert Missbauer , Order release planning by iterative simulation and linear programming: Theoretical foundation and analysis of its shortcomings, European Journal of Operational Research (2019), doi: https://doi.org/10.1016/j.ejor.2019.07.030
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights Iterating between production planning and scheduling level is a promising technique.
This has been applied especially for order release planning.
Iterative order release mechanisms often fail to converge.
The paper identifies the reason for this failure.
This is done by analyzing the theoretical underpinning of this iterative procedure.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Order release planning by iterative simulation and linear programming: Theoretical foundation and analysis of its shortcomings Hubert Missbauer
CR IP T
University of Innsbruck Department of Information Systems, Production and Logistics Management Universitaetsstrasse 15 6020 Innsbruck, Austria Email: [email protected]
Abstract
AN US
Planning future order releases for complex manufacturing systems with substantial lead times must consider limited capacities and the resulting production smoothing problem as well as the highly nonlinear relationship between capacity utilization, work-in-process and flow times. An important solution approach that has been proposed for this problem is to iterate between an order release model with fixed lead times and a lead time estimation (usually simulation) model that estimates the flow
M
times for given order releases, providing the lead times for the next iteration. However, the convergence of this iterative procedure is highly unpredictable, limiting its practical use.
ED
The iterative mechanism is analyzed analytically for simplified formulations of the order release and lead time estimation model. We show that an order release procedure of this type that iterates on the lead times is a dual (price) coordination mechanism whose design does not meet the theoretical
PT
requirements, and there is no straightforward way to overcome this. The analysis also provides insights into the results of several numerical studies from the literature and suggests a possible
CE
research direction to improve the method. Order release mechanisms of this type are a special case of a broader class of production planning methods that iterate between the production planning and the production scheduling level in
AC
order to provide realistic values for lead times and planned capacities. Providing theoretical underpinning for this type of production planning methods is an important research objective, and the paper pursues this direction for the special case of order release.
Keywords: Manufacturing; Production planning & control; Order release; Simulation
2
ACCEPTED MANUSCRIPT
1 Problem description Manufacturing planning and control (MPC) is a complex task encompassing decisions with different time horizons and impacts on the manufacturing or supply chain processes, which must include mechanisms that compensate for changes in planning data as well as unplanned events. Due to this inherent complexity the overall MPC problem is usually decomposed hierarchically into several
CR IP T
interdependent subproblems, such as master planning/master production scheduling, MRP, lot sizing, and detailed scheduling, including capacity planning at the respective levels.
We consider one specific decomposition, namely the structuring of the MPC system into the planning and the scheduling level. While the planning level usually determines quantities for production, transportation, sales, etc., for discrete planning periods, that is, it determines flows, the
AN US
scheduling level is event-oriented and determines assignments of operations to resources, sequencing of operations and their exact timing; hence planning and scheduling are very different in nature (Dauzère–Pérès and Lasserre, 2002).
Given this hierarchical decomposition, the planning level requires information on operational characteristics of the manufacturing system, most important the effective work center capacities (net
M
of idle times) and flow times or their probability distributions that are the basis for determining planned lead times. Since these operational characteristics are strongly influenced by scheduling decisions, the planning level must anticipate the behavior of the manufacturing system controlled by
ED
scheduling using an anticipation function as discussed by Schneeweiss (2003). Reliable estimates of effective capacities and flow times are necessary to determine appropriate values of the parameters
PT
describing capacities and planned lead times in the planning model. Determining planned lead times has been widely explored especially in the MRP literature (see
CE
Milne et al., 2015). Lead time setting must consider the high variance of flow time distributions in practice as well as the highly nonlinear relationship between mean flow times and resource utilization well-known in queueing theory. This complex behaviour results from the discrete and stochastic
AC
nature of the material flow and the autonomy of the scheduling level with respect to detailed scheduling decisions, and limits our ability to accurately estimate the lead times or, more generally, to identify feasible combinations of output quantity and lead times. It also means that lead times should be adjusted to capacity utilization, which requires time-dependent lead times if resource utilization varies over time due to time-varying demand. Iterative coordination of planning and scheduling is thus an obvious approach: The planning level uses a set of parameter values representing capacities and lead times to generate a production plan that is executed by a scheduling model that serves as an anticipation function. The resulting work center
3
ACCEPTED MANUSCRIPT
utilizations and/or flow times are used to update the capacity and/or lead time parameters in the planning model, which produces a revised production plan, until – hopefully – the algorithm converges, that is, the planned work center capacities and/or lead times used in the planning model and the actual utilization levels and/or lead times resulting from the scheduling model agree within a sufficiently small tolerance (we shall specify this more precisely later). An overview of production planning mechanisms of this kind is given in Section 2.
CR IP T
Most research on iterative planning and scheduling algorithms has been performed on order release planning. Order release is a decision problem that results from the interpretation of the planning – scheduling hierarchy described above as an organisational structure, which in many practical settings is a reasonable structure for the MPC system: The planning level coordinates the material flow over the entire logistic chain by setting coordinated production targets for production units that perform well-defined portions of the manufacturing process (e.g., a component
AN US
manufacturing shop or a manufacturing cell) and are autonomous with respect to detailed scheduling. The production targets set by the planning level are the work orders and their required due dates and consider adequate capacity constraints and planned lead times for the production units (for a thorough description of this decision structure, see de Kok and Fransoo, 2003; Bertrand, Wortmann, and Wijngaard, 1990). Detailed scheduling is performed locally within the production units. In this MPC
M
structure, order release, defined as the transfer of the control over the work orders under consideration from the planning to the scheduling level, is an important decision problem. Controlled order release
ED
keeps the work in process (WIP) level and – via Little’s law – the average flow time under control, while timely release of the orders should guarantee high due date performance. The flow time of a work order is the time elapsed from release until its completion, that is, the release time is the arrival
PT
time of the order at the first work center on its routing. Within this decision structure order release controls the arrival times of the work orders at the
CE
production unit which, given the behavior of the production unit that can be modelled as a queueing network, largely determines average flow time and output. Anticipating these effects of order release
AC
decisions is difficult, especially if the target values for WIP and flow times vary over time due to time-varying demand. Coordinating order release planning and detailed scheduling, the latter represented by a simulation or transient queueing model, by iterating between these levels and adjusting the planning parameters based on the simulation results is an important research direction. If simulation is used this can be considered as a hybrid analytical/simulation method (for reviews, see Shanthikumar and Sargent, 1983; Figueira and Almada-Lobo, 2014). The potential advantages of this approach are obvious: A very complex order release problem including stochastic elements can be solved by combining two well-known techniques (linear programming and simulation) using
4
ACCEPTED MANUSCRIPT
commercial software, and the techniques, especially in the case of simulation, can be explained to practitioners relatively easily. However, these iterative mechanisms lack any guarantee of optimality. Worse, in some cases they fail to converge and thus do not provide a feasible solution. This largely precludes practical application. Recently the performance of these iterative mechanisms under various conditions has been explored carefully by means of numerical studies, but there are virtually no theoretical insights
CR IP T
into how these mechanisms work and the nature of the underlying coordination mechanism. Therefore, no satisfactory explanation of the convergence problems has been provided. This reveals a major research gap: Solving the overall planning and scheduling problem by iteratively solving the planning and scheduling subproblems is a decomposition and coordination approach, but there has hardly been any attempt to integrate these iterative mechanisms into the theory of mathematical decomposition. Once the formal decomposition and coordination method implied by a specific
AN US
iterative planning/scheduling algorithm is identified, the relevant mathematical theory can, at least in principle, be applied to identify whether or not this iterative mechanism meets the requirements derived from that theory.
In the following this line of research is pursued for order release algorithms that iterate between an order release model with fixed lead times and a simulation or transient queueing model that
M
estimates the flow times.
The paper is organized as follows: In Section 2 the relevant literature is reviewed and the iterative
ED
order release mechanisms under study are described using Hung and Leachman (1996) as the main reference. Section 3 presents a simplified order release model encompassing only one work center as an instrument for the subsequent analysis. The possibilities of iteratively solving this model are
PT
studied in Section 4. Section 5 illustrates the insights by numerical experiments, while Section 6
2
CE
presents the conclusions and possible future research directions for this work.
Related literature
AC
Production planning methods that iterate between the planning and the scheduling level in order to determine correct values for parameters of the production planning model have been developed for various production planning problems. We first give a more general overview and then focus on iterative linear programming (LP) – simulation algorithms for order release planning. Dauzère-Pérès and Lasserre (2002) present an integrated model for lot sizing and scheduling and an algorithm that iterates between a lot sizing module that assumes a fixed production sequence and a scheduling module that sequences the given lots. Wolosewicz, Dauzère-Pérès, and Aggoune (2015) use Lagrangian relaxation to coordinate lot sizing and scheduling for a fixed sequence of operations
5
ACCEPTED MANUSCRIPT
on the machines. In both papers the feedback from the scheduling level represents the effective capacity available. Negenman (2000) presents an algorithm that iterates between an LP model that calculates the production plan for a production network and a flexible flow shop scheduling model that is solved by a heuristic. The feedback information provided by the scheduling level is the completion times of the orders, which are used to update the available capacities of the respective resources. Albey and Bilge (2011) present a hierarchical production planning and control system
CR IP T
framework for a Flexible Manufacturing System that consists of three levels: aggregate planning, loading and detailed planning. The behavior of the shop floor for a given production plan is anticipated by means of simulation. The simulation result is used to update the coefficients that represent capacities in the upper-level modules. Convergence is not analyzed in detail, but the authors indicate that capacity updating is complex due to the special problem and decision structure. Bilgen and Çelebi (2013) present an integrated production scheduling and distribution planning model where
AN US
the operation times are updated by iterating with a simulation model. Convergence of the operation times seems to be satisfactory. Since both preventive and corrective maintenance determines capacity losses, iteratively adjusting production planning and maintenance has been proposed in Zhao, Wang, and Zheng (2014). Note that all these algorithms essentially iterate on the capacities. Most iterative planning – scheduling algorithms have been developed for order release planning
M
where flow times are load-dependent and can be controlled by order release decisions. Multi-period models for order release planning usually represent the production unit as a network of work centers
ED
m=1,..., M. The planning horizon is divided into discrete planning periods t=1,..., T. The material flow is represented by inventory balance equations for work in process (WIP) at each work center and for final products, usually distinguishing different products or families of products with similar routing.
PT
The resulting network flow model is a fluid model in discrete time. Since the effect of order releases on lead times must be anticipated in order to determine optimal release quantities over time, the
CE
highly nonlinear relationships among WIP, flow time and throughput must be represented in the model. If fixed, exogenous lead times are assumed the model can be formulated as a linear program
AC
(LP) with the release quantities per period as the main decision variables. This model structure reflects the material flow structure in a variety of manufacturing systems for discrete manufacturing, especially in component manufacturing shops. Much of the research in this field was motivated by semiconductor manufacturing, and the simulation models used represent a wafer fab. The iterative order release algorithms analyzed in the following consist of an order release model (usually linear programming) that calculates optimal order releases per period assuming timedependent lead time parameters, and a lead time estimation (usually simulation) model. The simulation receives the release quantities as input and yields the predicted flow times and/or actual
6
ACCEPTED MANUSCRIPT
production quantities over time which are used to update the capacities and/or lead times of the release model for the next iteration. Zäpfel (1984) proposed an order release mechanism that iterates between an LP release model with fixed lead times and a simulation model for lead time estimation. The release model is formulated without capacity constraints since the effects of overloading (congestion) are included in the flow times that are the feedback from simulation. No numerical tests are provided.
CR IP T
The first to provide numerical tests of this type of order release planning were Hung and Leachman (1996), who developed an iterative release mechanism and provided tests using a model of a semiconductor wafer fab. Since we consider this as a reference work, the procedure is described in more detail.
Hung and Leachman (1996) use a release model similar to the step-separated formulation in Leachman and Carmon (1992). The lead times are fixed parameters and are defined as the time
AN US
required for a lot of product i to reach operation j after being released into the plant. The lead times are epoch-based, that is, they are defined for the period boundaries and thus can vary within the period. Hung and Leachman use backward lead times at the period boundaries, defined as the lead time of an order that is finished (or arrives at a certain operation) at the start/end of a specified period.
AC
CE
PT
ED
M
This is depicted in Figure 1.
Figure 1: Epoch-based backward lead times to relate output to wafer releases. Source: Hung and Leachman (1996, p. 260)
In the Hung & Leachman model the releases and backward lead times determine both the arrival times at the work centers and the periods in which capacity at the work centers is required for
7
ACCEPTED MANUSCRIPT
production; the lead times are not modeled as time windows within which the production periods can be chosen, as in de Kok and Fransoo (2003) and Pürgstaller and Missbauer (2012). The simulation model collects statistics on the release times of the manufacturing lots (equivalent to work orders in the above description) and the time an operation on a lot is initiated. The time elapsing between these two events gives the flow time of this operation, as defined in the stepseparated formulation. The flow times at the period boundaries are interpolated from the flow times of
CR IP T
lots that initiate processing immediately before and immediately after the boundary epoch (Hung and Leachman, 1996, p. 262).
Hung and Leachman (1996) examine the rate of convergence of the flow time estimates to the actual flow times in the simulation and find that convergence to the correct expected flow time values can be quite rapid, but that the procedure can fail to converge in some cases which are not fully understood. Irdem, Kacar, and Uzsoy (2008) provide detailed numerical tests and report that for high
AN US
utilization of the bottleneck facilities the results of the Hung/Leachman procedure “are
…..
discouraging in terms of convergence”. For the low utilization case the results are classified as “rather discouraging”, but “when the bottleneck utilization is low, the value of the iterative procedure is very limited in terms of its ability to obtain improved system performance.” Reducing the period length from seven days to one day did not lead to a fundamental change of this pattern (Irdem, Kacar, and
M
Uzsoy, 2010, p. 9); the authors state that “…LP model and simulation model do not represent each other sufficiently.”
ED
Riaño (2003) uses a similar iterative procedure where the LP release model expresses the lead times as weights wkt that denote the fraction of input released in period k that will emerge as output in
PT
period t. The wkt are estimated using a transient queuing model. The algorithm iterates on the weights rather than on the lead times as in Hung and Leachmann (1996). Riaño (2003) reports that the algorithm does not always converge and that convergence may be affected by the choice of the initial
CE
condition. “However, when it does converge, it converges to the same solution regardless of the initial solution used. When it did not converge, it would cycle through two or three solutions” (p. 83).
AC
Hung and Hou (2001) use the same basic procedure as Hung and Leachman (1996) but replace
the simulation model with a function relating expected flow times to the loading rate at individual work centers, similar to that used by Voss and Woodruff (2003). The loading rate is defined as the machine hours of all orders arriving at the work center in the respective period divided by the work center capacity, its functional relationship to the expected flow time is estimated by simulation. Their attempt to use an M/M/s queueing model proved unsuccessful. They report short computation times and good convergence for short sub-periods, but do not explicitly comment on whether there were convergence problems in some cases.
8
ACCEPTED MANUSCRIPT
Byrne and Bakir (1999) use a simulation model to determine realistic estimates for the available capacities, Byrne and Hossain (2005) present an extended production planning model within this framework. Bang and Kim (2010) and Kim and Lee (2016) present refined methods for semiconductor production planning that iterate with a priority-based scheduling or simulation model, respectively. Kim and Kim (2001) simultaneously update the lead times and available capacities assumed in the
CR IP T
LP model using this iteration scheme. Simulation is used to anticipate the flow times and the utilization of the work centers. As in Riaño (2003), the lead times are expressed as loading ratios, defined as the proportion of the start quantity released in a period p that contributes to the workload at machine k in period t (Kim and Kim, 2001, p. 168).
Kim and Kim (2001) do not report convergence problems in their numerical tests. Irdem, Kacar, and Uzsoy (2010) report good convergence of this procedure for both high and low utilization and
AN US
conclude “that the convergence behavior of the KK [Kim & Kim, 2001] procedure is qualitatively different from that of the HL [Hung & Leachman, 1996] procedure” (p. 452f.). They conjecture that this difference might be due to the differences in the estimation of the loading ratios which are estimated directly from the simulation in the Kim & Kim procedure whereas the Hung & Leachman procedure estimates the flow times of individual lots at specific epochs from the simulation, and uses
M
these to calculate an estimate of the loading ratios. (p. 453). However, Albey, Bilge, and Uzsoy (2012) perform numerical studies and find that this good convergence performance of the Kim & Kim
ED
procedure holds only for aggregate convergence criteria, such as objective function value and total production, but “achieving convergence at the level of production levels of each product in each period proves to be difficult” (p.11). They conclude that the ability of these approaches to allocate
PT
resources efficiently to products across multiple time periods is questionable. The body of research we have reviewed yields the following insights: (1) Convergence of iterative
CE
LP and lead time estimation procedures is a difficult issue that substantially hinders their practical implementation. However, including the capacities into the iterative mechanism as in Kim and Kim
AC
(2001) seems to mitigate the convergence problems substantially; some papers on iterative planning scheduling algorithms that essentially iterate on the capacities report satisfactory convergence. (2) The iterative methods for order release planning are not derived from a well-defined decomposition and coordination method. This is in contrast to several methods that iterate between two production planning levels, such as the Lagrangian methods presented in Graves (1982) or the fix-and-optimize approach of Bahl and Ritzman (1984). Since the underlying theory is unclear the convergence behaviour has been treated as a purely empirical question. (3) The convergence problems of iterative methods for order release planning that iterate solely on the lead times seem to be independent of specific details of the methods and their implementation. Riaño (2003) encountered convergence
9
ACCEPTED MANUSCRIPT
problems even for single-stage systems, so we expect that the complex material flow structure in the lead time estimation model (e.g., wafer fab with re-entrant material flows) is not the decisive factor. Convergence problems occur in order release models that include or discard WIP holding costs or capacity constraints, respectively. Likewise, iterating on the lead times (Hung and Leachman, 1996) and the loading ratios that are derived from the observed flow times (Riaño 2003) both produce this phenomenon.
CR IP T
Therefore, we hypothesize that the convergence problem is caused by the basic mechanism of this iterative approach, and inferring the theoretical basis from the proposed iterative technique should allow us to identify this basic mechanism. We pursue this idea in the following analysis.
3
The simplified order release model
AN US
In this section we formulate an order release model for a single production unit that retains the essential elements of the iterative mechanisms while remaining sufficiently tractable to obtain the necessary insights. We prove analytically that iterative order release mechanisms that iterate on the lead times imply a well-defined decomposition method, namely price coordination based on Lagrangian relaxation, but implement it in an incomplete manner that does not meet the theoretical
M
requirements especially for convergence. Therefore, they cannot be expected to converge to the optimum since the optimum is not a fixed point of the iterations. We argue that this behavior is inherent to the structure of the procedure rather than the complexity of real-life situations.
ED
The order release model for a production unit within the iterative approach is characterized by the following assumptions: The production unit consists of a single (bottleneck) work center that
PT
produces a single product. Time is divided into planning periods t=1,…., T. The orders, expressed as hours of work, are released and arrive at the WIP inventory, are processed and after processing arrive
CE
at the finished goods inventory. Demand for each period is known and deterministic. For each period t there is a fixed, integer backward lead time vˆtb ( vˆtb is an integer multiple of the period length) that denotes the lead time of products that are finished in period t. The maximum capacity per period,
AC
denoted by C, is constant over the planning horizon. Later we extend the model in order to capture the emergence of the lead times. The decision variables (all measured in hours of work) are the work released in period t, denoted
Rt , and the output (actual production) X t in period t, the WIP Wt at the end of period t and the finished goods inventory I t at the end of period t. The parameter Dt denotes the demand in period t (measured in hours of work), while h and l denote the holding costs per unit and period for WIP and final products, respectively. The initial conditions are represented by zero initial WIP ( W0 0 ) and a
10
ACCEPTED MANUSCRIPT
finished goods inventory level I 0 . We assume that these initial conditions, together with the capacity per period C, allow a feasible solution. The order release model OR_ITER can be formulated as follows:
h Wt l I t t
subject to
(1)
t t
t
k 1
k 1
Wt Rk X k t
t
k 1
k 1
I t I 0 X k Dk t vˆtb
t
R X Xt C
k 1
k
X t , Rt , I t ,Wt 0
(2)
for all t
(3)
for t vˆtb 1,..., T
(4)
for all t
(5)
AN US
k 1
k
for all t
CR IP T
min
for all t
(6)
The objective function (1) minimizes the costs for WIP and finished goods inventory (FGI). Equations (2) and (3) are the respective inventory balance equations. For exposition we assume zero initial WIP implying Rt 0 for t 0 (which can be relaxed by defining fixed releases for past
M
periods). Equation (4) states that for each period t the cumulative output at the end of period t must be
ED
equal to the cumulative releases at the end of period t vˆtb , implementing the specified backward lead times vˆtb . For given vˆtb as assumed here this is a set of linear equations that relate the Rt and the Xt.
PT
The vˆtb must satisfy a no-passing property, that is, work finished in period t+1 cannot be released earlier than work finished in period t (see Carey, 1992, for the no-passing conditions on the lead
CE
times). If vˆtb1 vˆtb 1 , the output in period t and in period t+1 are released in the same period t vˆtb . Therefore, the constraints are formulated as inequalities. The WIP holding costs in the objective function guarantee that the constraints hold at equality except in this special case, and that the
AC
specified backward lead times are met. (Alternatively, Equation (4) can be formulated using loading ratios as described above.) Equation (5) shows the capacity constraints. All variables are non-negative (Equation (6)). The model exhibits the essential structure of the order release model within the iterative mechanisms: In order to be finished at a certain work center in period t, the work is released in period
t vˆtb . The lead times can be achieved as long as the capacities are not overloaded. Production capacity is consumed at the end of the lead time. The lead times vˆtb applied in iteration k are derived
11
ACCEPTED MANUSCRIPT
from the flow times obtained from the simulation run in iteration k-1 and are parameters in the release model. For the first iteration appropriate initial values for the lead times must be determined. The flow times are deterministic since OR_ITER is a fluid model with first-come-first-served sequencing. Thus, the problem of deriving deterministic lead time parameters from observed flow time distributions does not occur here; this is actually an additional degree of freedom. Holding WIP inventory incurs costs. In the Hung & Leachman model the costs of releases and the revenues from
CR IP T
completions are discounted which leads to an increase in the objective function as WIP is held. The subsequent analysis is performed on this model. The iterative release mechanisms developed in the literature often are based on less restrictive assumptions and allow, e.g., multiple work centers, fractional lead times (Hung and Leachman 1996) or lead time distributions (Riaño 2003). Our reasoning is as follows: If it can be proven that even in the simplified case represented by model OR_ITER the iterative mechanisms do not lead to the optimum, this is sufficient to disprove
sufficient for falsification). We offer the following essential insight.
AN US
the general claim that the iterative mechanisms lead to the optimal solution (one counterexample is
Theorem 1: The order release model OR_ITER with given time-varying backward lead times vˆtb for orders (or work) behaves as an order release model with zero lead time and time-varying
M
production costs (that is, it yields the same optimal solution except time-shifted Rt ) as long as the
ED
constraint that the release periods t vˆtb must be positive is not binding. The WIP cost term in the objective function (1) is
h Wt h vˆtb X t .
(7)
t
PT
t
Hence the WIP holding costs for backward lead times vˆtb are equivalent to production costs of
CE
h vˆtb for each unit of work that is produced and occupies capacity in period t. Proof: See Appendix 1.
AC
The intuition of the proof is the following: Finishing an order in period t implies that its lead time is vˆtb , which is also the time it spends in WIP. Thus a unit of output completed in period t will incur a WIP holding cost of h times vˆtb . Since the WIP holding costs are assumed to be proportional to the production quantity they can be viewed as production costs; the impact is the same. Note that the planned release volumes Rt have no further influence on the objective function and, as long as there are no further constraints on the releases, are not relevant to the result. The optimal
Rt can be obtained by offsetting the optimal X t by the backward lead times vˆtb . This does not hold if 12
ACCEPTED MANUSCRIPT
release periods obtained in this way would be before period 1 since the release quantities for the past periods are fixed. Therefore, it is assumed that the initial conditions that are determined by the past releases are such that this case does not occur. Corollary 1: If value is added during production, which means h l , and production maintains the no-passing property, which implies vˆtb1 vˆtb 1 , then the optimal production quantities X t do not depend on the backward lead times vˆtb . The X t that are optimal for instantaneous production
CR IP T
( vˆtb 0 for all t) remain optimal for all vˆtb 0 that satisfy the no-passing property.
Proof: Model OR_ITER with vˆtb 0 for all t is an LP production planning model with instantaneous production that produces output as late as possible. Demand in some period t is produced in an earlier period s
AN US
production from period t to period s
vˆtb vˆsb
l (t s ) . h
(8)
M
Since vˆtb vˆsb t s due to the no-passing property and l h 1 due to the assumption that value is added, this condition can never be satisfied.
□
ED
Theorem 1 provides an important insight: The iterative mechanism coordinates the order release model and the lead time estimation model using a price coordination mechanism. Thus we can apply the theory of price coordination (dual coordination; see Dirickx and Jennergren, 1979) to evaluate the
PT
theoretical consistency of this iterative mechanism. We must check, firstly, whether the order release and the scheduling problem can be coordinated by this price coordination mechanism and, secondly,
CE
whether the iterative mechanism is a valid price adjustment procedure, that is, whether the lead times
vˆtb derived from the lead time estimation model lead to correct prices (steering costs) for optimizing
AC
the release and output volumes Rt and X t , respectively. This requires formulating the monolithic model that jointly determines the optimal release volumes and lead times and which is solved by means of the iterative mechanism. In order to formulate the monolithic model implied by the iterative mechanism the assumptions of model OR_ITER are complemented/replaced by the following assumptions:
The backward lead times vˆtb are decision variables that result from the releases Rt and the outputs X t of the relevant periods.
13
ACCEPTED MANUSCRIPT
The capacity constraints (5) are not sufficient to describe the possible output. The dynamics of the material flow (variability effects) can lead to idleness in a period t even though the planned available work (termed load) Wt 1 Rt is higher than the capacity. Hence there is an upper bound on the output that results from the dynamics of the material flow and can be calculated by the lead time estimation model. In the following we assume that the upper bound on the output is determined by the system dynamics prior
CR IP T
to the period under consideration (history of the process). (Including expectations about future periods, e.g., dispatching rules that anticipate future events, is not considered here.) We do not specify the logic by which the (deterministic) upper bound to the output is derived from the stochasticity of the material flow. It can be the conditional expectation or some quantile of the conditional output distribution given a certain history of the
AN US
process. This is analogous to estimating clearing functions from simulation data which can be done by ordinary least squares or quantile regression (see Asmundsson, Rardin, Turkseven, and Uzsoy (2009)).
Including the simulation or transient queueing model for output and lead time estimation into the order release model renders the resulting “monolithic” model unsolvable by mathematical programming methods. The reasoning for this is given below. Therefore,
The model is modified as follows:
M
iterative methods are applied for solving the model.
ED
Upper bound on the output
The following constraints are added:
for all t
(9)
PT
X t ( R [ 1,..., t ], X [ 1,..., t 1])
The output in a period t is a function of the releases and output in the previous periods and of the
CE
releases in period t. The function represents all effects of the dynamics, including the stochasticity, of the material flow within the production unit and is expressed as a simulation or transient queueing model that predicts the output in a period t from the releases until period t and the history of the
AC
output quantities. At this point we do not specify the mathematical properties of this function. The constraint is formulated as an equality constraint since simulation models usually do not keep resources idle if work is available. Deduction of the backward lead times from WIP and output For given releases Rt and outputs X t the lead times vˆtb can be calculated from Equation (4) provided they are integer. Non-integer lead times can be calculated as follows (for a comprehensive treatment, see Hackman, 2008, p. 320 ff.): We define X tcum , Rtcum as the cumulative output and release
14
ACCEPTED MANUSCRIPT
quantities, respectively, at the end of period t. X cum (t ), R cum (t ) are defined analogously in between the period boundaries, assuming uniform releases and output within the period: t
Rtcum Rk
for t 1,2,..., T
(10)
for 0 t T and t non-integer
(11)
k 1
Rcum (t ) Rcum Rt 1 (t t ) t
CR IP T
with t the integer part of the real-valued number t (largest integer less than or equal to t). The definitions of X tcum and X cum (t ) are analogous. The backward lead time vˆtb for orders that are finished at the end of period t is defined by
for t 1,2,..., T
W0 Rcum (t vˆtb ) X tcum .
(12)
In the following the analysis is performed assuming integer lead times.
AN US
Reformulation of the objective function
The objective function (1) can (but need not) be reformulated by substituting (7) for the WIP cost term, which yields:
h vˆtb X t l I t .
min
t
t
M
The monolithic model OR_MON that is solved by the iterative mechanisms is then: Objective function:
h Wt l I t
or, equivalently:
t
ED
min
(13)
t
h vˆtb X t l I t
PT
min
t
(14)
t
AC
CE
Inventory balance equations
t
t
k 1
k 1
Wt Rk X k t
t
k 1
k 1
I t I 0 X k Dk
for all t
(15)
for all t
(16)
for all t
(17)
for all t
(18)
Relationship between release and production volumes t vˆtb
t
R X k 1
k
k 1
k
Output anticipation
X t ( R [ 1,..., t ], X [ 1,..., t 1])
Capacity constraints (redundant if Equation (18) is correctly specified)
15
ACCEPTED MANUSCRIPT
Xt C
for all t
(19)
X t , Rt , I t ,Wt , vˆtb 0
for all t
(20)
If the objective function (14) is used, the mechanics of the model OR_MON are as follows: The release and output variables (Rt and Xt, resp.) lead to a certain WIP level (Equation (15)) and – by Equation (17) – to certain lead times that, in turn, determine the first term in the objective function. Note that in this lead time formulation of the model, Equation (17) is an integral part of the model. In
CR IP T
the general case this formulation cannot be solved by mathematical programming techniques for two reasons: Firstly, the function is implicit in a simulation or transient queueing model and exhibits a complex structure even if approximations are used (e.g., integer and nonlinear as shown in Missbauer, 2009). Secondly, the lead time estimation constraints (17) contain variables that are indices of other variables ( vˆtb in the index of Rk), which renders the model unsolvable with mathematical
AN US
programming methods. This structure occurs because the lead times are modeled correctly as a result of WIP and output (see Hackman, 2008, p. 338 ff.). If lead times are modeled without explicit reference to WIP (as in Voss and Woodruff, 2006) this difficulty can be avoided at the price of imprecision in modeling the dynamics of the lead times.
The second difficulty described above (variables as indices of other variables) can be avoided
M
provided that there are no cut-off effects at the start of the planning horizon; that is, extending the planning horizon infinitely into the past does not change the optimal solution. In this case (which is
ED
assumed in the following) there are no relevant constraints on the lead times. Using the original formulation of the objective function (Equation (13)) changes the mechanics of the model, yielding the WIP-oriented formulation: The release and output variables (Rt and Xt, resp.) lead to a certain WIP
PT
level (see Equation (15)) that determines the first term in the objective function (13). The lead time variables are not required for the optimization, but are deduced from the optimal solution by Equation
CE
(17), which poses no problem. Of course the difficulty with the structure of Equation (18) remains. The lead time oriented and WIP-oriented formulation of the model are equivalent although they
AC
exhibit different structures for the optimization. Hence in the following we use the WIP-oriented formulation which is easier to handle and deduct the lead times if necessary.
16
ACCEPTED MANUSCRIPT
4
Analysis of the iterative solution of the model
4.1 Method of the analysis We can now ask under what conditions the model OR_ITER, together with a lead time updating
CR IP T
procedure, can achieve the optimal solution to the monolithic model OR_MON. A necessary condition for this is that there exists a fixed point of the iterative mechanism and the optimal solution is a fixed point. A fixed point of the iterative mechanism is a solution where the time-dependent lead time vector that leads to a certain order release plan is equal to the lead time vector that results from these order releases and from the subsequent lead time estimation (e.g., simulation). Thus if the iterative mechanism reaches a fixed point the solution will remain there. Defining the desired property
AN US
of the optimal solution as being a fixed point avoids ambiguity in the definition of convergence: If an iterative order release mechanism converges in the objective function value, it need not converge in terms of the production quantities of the products (see Section 2).
Per the Banach Fixed Point Theorem (O'Regan, Meehan, and Agarwal, 2001) convergence to
that yields the lead time vector (from the lead
M
a fixed point is guaranteed when the function g vˆ b
time estimation, e.g., simulation model) as a function of the lead time vector vˆ b used in the order
ED
release model is a contraction mapping. Due to the complexity of the function g (.) that consists of an LP order release and a lead time estimation model we do not pursue this direction. Instead we analyze whether the optimal solution to the monolithic model OR_MON is a fixed point of the iterative
PT
mechanism that consists of model OR_ITER and an appropriate lead time estimation model. This is
the case if g vˆ b * vˆ b * with vˆ b * the optimal lead time vector, that is, when the release model
CE
with optimal lead times as parameters leads to a release plan that, executed by the lead time estimation model, again yields the lead time vector vˆ b * .
AC
This implies the two questions outlined in Section 3: The first of these is whether there exist lead
time parameters in the release model OR_ITER that lead to the optimal solution of the monolithic model OR_MON; this refers to the coordinability (Mesarovic, Macko, and Takahara 1970) of order release and lead time estimation problem by this iterative mechanism. The second is under what conditions – if at all – the lead time estimation model can yield these lead time values, which refers to the validity of the price adjustment procedure. We start with the first issue. For the following analyses we assume that the constraint that the release periods t vˆtb must be positive is not binding.
17
ACCEPTED MANUSCRIPT
In order to make model OR_MON solvable we consider only a special case of the output constraint (18): The output in period t is assumed to depend only on the load (available work) in period t ( Wt 1 Rt ) (see Equation (21)). Again we argue that if it can be proven that even in this special case the iterative mechanism does not lead to the optimum, this disproves the general claim that the iterative mechanism leads to the optimal solution. Furthermore, it is formulated as an inequality constraint such that the load determines the maximum output.
CR IP T
Thus we obtain a conventional (one-dimensional) clearing function (see Equation (22)):
( R [ 1,..., t ], X [ 1,..., t 1]) f (Wt 1 Rt )
(21)
X t f (Wt 1 Rt ) 0
(22)
for all t
Clearing function models with a concave, saturating clearing function were introduced by Karmarkar (1989) and have been described extensively in the literature (e.g., Missbauer and Uzsoy,
AN US
2011). The feasible set is convex (see Carey, 1987, which analyzes a model with this structure for traffic flow control), and the model can be solved either by nonlinear programming or by linearizing the clearing function and solving the resulting linear program. Carey (1987) proves that if certain assumptions are satisfied, the output is always at its upper bound, which in our context means that the clearing function constraints (22) hold at equality. One of these assumptions, termed the cost function
M
assumption, requires that “the marginal cost per unit volume, per unit time should not increase on moving from one arc to the next” (p. 63f.). This assumption is not satisfied here because the FGI holding costs l usually are higher than the WIP holding costs h. However, in our single-stage model
ED
the only purpose of the WIP level is to enable the desired output, and X t f (Wt 1 Rt ) 0 only occurs if demand falls very sharply within a short period of time (Kefeli, Uzsoy, Fathi, & Kay, 2011,
PT
p. 387). This issue will be considered below when necessary. In the course of the analysis we will also benefit from some structural insights into this type of models that are available in the literature. We
CE
use the WIP-oriented formulation of the objective function; the lead time anticipation equation (17) is not required for the optimization, it just calculates the lead times from WIP and output.
AC
4.2 Analysis for arbitrary utilization The iterative release mechanism anticipates the dynamics of the production unit external to the
optimization model. Hence the output constraint (22) is not a part of the optimization model OR_ITER; capacity is only constrained by Equation (19) or (5), respectively. As stated in Theorem 1, the lead times define the prices for producing items (or for occupying capacity, respectively) in certain periods. Thus in model OR_ITER the relevant capacity constraints of the monolithic model (Equation (22); note that the capacity constraints (19) are redundant) are relaxed and replaced by
18
ACCEPTED MANUSCRIPT
defining prices for the capacity. This is a Lagrangian relaxation of the constraints (22), and the optimality conditions for the prices (and hence for the lead times) can be derived from standard theory. The Lagrangian function obtained by relaxing the clearing function constraints (22) of model OR_MON is
£( X t , ut ) h Wt l I t ut X t f (Wt 1 Rt ) t
(23)
t
or
CR IP T
t
£( X t , ut ) h Wt l I t ut X t ut f (Wt 1 Rt ) , t
t
t
t
(24)
where ut denotes the dual variables associated with the clearing function constraints (22) in period t.
AN US
The Lagrangian function (24) is the objective function of model OR_ITER (Equation (1)) extended by a price term for occupying capacity and by the term
u f (W t
t 1
Rt ) . Model
t
OR_MON with the simplified output constraints (22) satisfies Slater’s condition since in general there exists at least one solution with positive decision variables where the clearing function constraints
M
hold with strict inequalities (Boyd and Vandenberghe, 2009, p. 226), and due to its convexity the duality gap is zero. Maximizing Equation (24) over ut and minimizing it over the primal variables
ED
(subject to the remaining constraints of model OR_MON) leads to the optimal ut (partial duality; see Luenberger, 1989, p. 401). The implications of this structure for the validity of the iterative order release mechanism are given in the following theorems.
PT
Theorem 2: Minimizing the Lagrangian function (24) over the primal variables for given optimal dual variables (dual prices) ut* subject to the remaining constraints (15), (16), (19), (20) of model
CE
OR_MON need not lead to a feasible primal solution that satisfies the Karush-Kuhn-Tucker (KKT) optimality conditions due to multiple optima. Thus, in the general case no choice of the ut will lead to
AC
the optimal solution of model OR_MON. Proof: The structure of the multiple optima is shown in Appendix 2. (For the KKT conditions for
an extended clearing function model, see Kim and Uzsoy, 2008.)
□
Theorem 3: Solving model OR_ITER to optimality in general will not yield the optimal solution of model OR_MON for any choice of the lead times vˆtb . Proof: As shown above, model OR_ITER implies that the clearing function constraints are relaxed. In this case a solution to model OR_MON with optimal objective function value can only be
19
ACCEPTED MANUSCRIPT
obtained by minimizing the Lagrangian function (24) for optimal dual prices ut* . (24) is a nonlinear function due to the clearing function related term
u f (W t
t 1
Rt ) . The objective function (1) of
t
model OR_ITER does not include this term and is linear independent of the lead times vˆtb , that is, □
model OR_ITER does not minimize the Lagrangian function (24).
The intuition of these theorems is the following: Per Theorem 1 the lead times act as prices for
CR IP T
production, and it has been shown that this implies a price coordination mechanism that coordinates the order release and the lead time estimation model. Our mathematical analysis of the integrated model OR_MON and its Lagrangian decomposition shows that a necessary condition for this price coordination mechanism to work is that the Lagrangian function (24) is minimized (given the dual prices ut* ). However, model OR_ITER cannot do this because it minimizes the wrong function as
AN US
shown in Theorem 3. In addition, even if the Lagrangian function (24) is minimized for optimal dual prices the solution can be infeasible per Theorem 2.
There is no obvious way to overcome this fundamental problem since in general f (Wt 1 Rt ) represents a simulation model in practical applications, and integrating
u f (W t
t 1
Rt ) into
t
M
Equation (1) would require integrating a function coded in a simulation model into the mathematical programming model OR_ITER. Note that this structural problem of model OR_ITER – discarding the
ED
respective price term in objective function (1) – is not specifically related to clearing functions; it also holds for the general formulation of the output constraint (Equation (18)) although in this case the structure of the solution space is unknown unless function is specified.
PT
These insights lead to the following
Conclusion: It cannot be expected that any choice of the backward lead times vˆtb in model
CE
OR_ITER can lead to the optimal solution to model OR_MON. Thus, the optimal solution to model OR_ITER cannot be expected to be a fixed point.
AC
This does not exclude the possibility that there exists a fixed point under specific assumptions. Especially, if the clearing function and/or the load situation allow independence between order releases and lead times there is a fixed point. This is formulated in *
Theorem 4: If the lead times are independent of the order release quantities (vˆtb ) , then there
exists a fixed point and the algorithm will converge to this point. Proof: If the system is initialized with an arbitrary lead time vector vˆtb for all t, the resulting Rt(1) *
vector will lead to the lead time vector vˆtb by assumption. This leads to a second vector of order
20
ACCEPTED MANUSCRIPT
*
*
release quantities ( Rt(2) ) . Rt(2) and vˆtb , t=1, …. T, are a fixed point because (vˆtb ) can never change (by assumption) and leads to vector Rt(2) (by definition of Rt(2) ).
□
Theorem 4 refers to a situation where congestion effects due to limited capacities are not relevant. This cannot be assumed in general for capacitated production resources, but it can be a good approximation if capacity utilization is low.
CR IP T
However, it can be argued that, despite these structural deficiencies, model OR_ITER might lead to a reasonable heuristic solution to model OR_MON in certain cases provided that the prices for production implied by the objective function (Equation (1)) h vˆtb are equal (or sufficiently close) to the optimal prices that result from the Lagrangian relaxation of the clearing function constraint (22) in model OR_MON. This might be the case for high utilization, which is analyzed next.
AN US
4.3 The special case of high utilization
Although normally the clearing function depends on the values of the decision variables, the optimal solution might be in a region where the maximum output is close to the capacity, that is,
f (Wt 1 Rt ) C especially for bottleneck work centers that are of primary interest. In this case the
M
clearing function related term in Equation (24) is approximately an additive constant. Furthermore, the relaxed constraints (22) can be violated only to a minor extent due to the capacity constraints (19).
ED
We now analyze whether in this case the iterative release mechanism can lead to backward lead times
vˆtb that have the same impact as optimal dual variables (dual prices) ut* of the constraints (22). Since there is no duality gap and the possible infeasibility may be very small, this might lead to a reasonable
PT
solution to model OR_MON.
CE
For f (Wt 1 Rt ) C the Lagrangian function (24) is approximately:
AC
Substituting h
£( X t , ut ) h Wt l I t ut X t utC . t
W hvˆ X b t
t
t
t
t
t
(25)
t
(from Equation (7)) and collecting the Xt coefficients we get
t
£( X t , ut ) hvˆtb ut X t l I t utC , t
t
(26)
t
where the last term is an additive constant. Obviously, if correct lead times are assumed the coefficients of the Xt in the objective function of model OR_ITER ((1) or (14), respectively) are too small; the dual variable ut is added to the correct WIP holding cost rate hvˆtb in the Lagrangian function (26). However, this need not preclude a reasonable solution: The WIP holding cost
21
ACCEPTED MANUSCRIPT
coefficient h is difficult to estimate; production prior to demand triggers various processes in the company with very different cost structures. Proportionality of holding costs to quantity and time is a strong assumption, and holding cost coefficients can be argued to be theoretical constructs rather than observable values. Therefore, it is worth exploring under which assumptions an iterative procedure that implicitly assumes ut* 0 (as in objective function (1) of model OR_ITER) can lead to an
hvˆtb ut* . h vˆtb If ut* vˆtb this simplifies to
h h .
CR IP T
optimal solution for some modified holding cost coefficient h . This would require hvˆtb hvˆtb ut* or (27)
(28)
AN US
Hence if the backward lead times and the dual prices of the clearing function constraints are proportional, then ignoring the respective price terms in the optimization, as done in model OR_ITER, can be compensated by increasing the holding cost coefficient for WIP by a specified value. In other words, using a modified holding cost coefficient h h in the objective function (1) is consistent with the Lagrangian function (26) if h exceeds the actual value h by the value (from Equation (28)).
M
With this modification the Lagrangian function (26) can be minimized by minimizing objective function (1).
ED
In the following we analyze under what assumptions the condition ut* vˆtb that provides the possibility to compensate the structural deficiencies of objective function (1) by increasing the WIP
PT
holding costs is satisfied. Since we still analyse the special case f (Wt 1 Rt ) C , the workload in period t t Wt 1 Rt
CE
calculations.
is used as a proxy for the lead time vˆtb in the subsequent analytical
Kefeli, Uzsoy, Fathi, and Kay (2011) analyze order release models with a saturating clearing function approximated by a set of linear functions (tangents), distinguishing congested periods
AC
(Wt>0), non-congested periods (Wt=0 and Xt>0) and idle periods (Wt=0 and Xt=0). They find that in a congested interval, as long as work is released and FGI is carried in consecutive periods, the sum of the nonzero dual prices of capacity (more precisely: of the linear segments of the clearing function that are tight in a certain period) increases by l (the holding cost coefficient of FGI) for each consecutive time period the system holds FGI (p. 388). Since shifting the segments of the linearized clearing function that are tight upward (that is, increasing the intercept of both segments by 1) increases capacity by 1, this theorem holds for the dual price of the clearing function constraint. The proof extends trivially to nonlinear, saturating clearing functions by letting the number of tangents
22
ACCEPTED MANUSCRIPT
approach infinity and the difference of their slopes approach zero. The intuition behind this is the same as in production planning models with zero lead time: Increasing capacity by one unit in a period t of an interval where the capacity constraints are tight and FGI is held reduces FGI by the respective amount until period t. Hence the later the additional capacity is added, the larger the savings in FGI holding costs. This relationship was recognized very early from marginal cost analysis (Modigliani & Hohn, 1955).
CR IP T
Since in clearing function models additional capacity can be generated by higher WIP, there is a relationship with the WIP level: The independent variable of the clearing function, in our case the load, is only increased if the marginal benefit of the additional capacity that is generated (determined by the first derivative of the clearing function) is higher than the marginal cost of the higher WIP. Kim and Uzsoy (2008) analyze a nonlinear, concave (saturating) clearing function with Wt as
AN US
explanatory variable and obtain the KKT condition (Equation (11) on p. 394; we omit the possibility of engineering lots here)
h t 1 t
df ut dWt
(29)
with t the dual variable of the WIP balance constraint in period t. They show that for Wt>0 and
M
Rt>0 Equation (29) holds at equality and the term t 1 t vanishes. Assuming this and using the
ED
load t Wt 1 Rt as explanatory variable as in our model, this leads to
df ( t ) h . d ( t ) ut*
(30)
PT
Defining f '( t ) df ( t ) d ( t ) and a to be the intercept of the linear increase of the ut* as a function of t within the congested interval under consideration, that is, ut* a l t , we can write
CE
Equation (30) as:
f '(t )
h . a l t
(31)
AC
Since the clearing function is concave and strictly monotone in t , it has a well-defined inverse
function f '( 1) . Therefore, from Equation (31) we obtain
h t f '( 1) a l t
(32)
This can be interpreted as follows: For a given intercept a (which is determined by the initial conditions of the congested interval under consideration) the argument of f '( 1) in Equation (32) is determined for all periods of the congested interval and is decreasing over time (since l>0). The
23
ACCEPTED MANUSCRIPT
function f reflects the operational characteristics of the work center and is also given. Hence the load
t for all periods of the congested interval under consideration is given by Equation (32). Although we can assume that the function f is concave and saturating, the first derivative f '(t ) strongly depends on the characteristics of the arrival and departure process at the work center and on the ratio of period length and average service time. Hence we cannot draw any general conclusion on the
CR IP T
evolution of the relationship between ut* and t over time within a congested interval beyond the monotone increase of t (and the linear increase of ut* ) with t. In particular, we cannot assume that
ut* and t increase proportionally over time. Numerical examples lead to the following
Observation 1: Even within a congested interval the values of ut* and backward lead times vˆtb over time in model OR_MON can be far from proportional depending on the shape of the clearing
AN US
function. For general situations where congested and non-congested periods occur alternately there is no correspondence between ut* and vˆtb over time; in this case not even monotonicity can be assumed. For this situation to the best of our knowledge no analytical results for the dynamic behaviour of the dual prices are available.
The numerical experiments are demonstrated in Section 5 below.
M
However, Equations (30) and (31), respectively, indicate that there should be a functional form of the clearing function that leads to a linear increase of the load t over time within a congested
ED
interval and thus to a proportional relationship of load, which we use as a proxy for the lead time vˆtb , and the dual prices within a congested interval. This is the case if the clearing function exhibits the
PT
logarithmic shape f ( t ) K h Log t with K an integration constant. This is shown in Appendix 3. There is no justification for this functional form in queueing theory, and it is inconsistent
CE
with common assumptions about service processes since it predicts negative output if the load is below a certain positive value. Furthermore, it is not consistent with the condition f (Wt 1 Rt ) C
AC
for high load which provided the justification for this analysis. The above analysis reveals a second limitation formulated below. Theorem 5: We consider a congested interval with ut* vˆtb for all t . For h 0 the
condition h h (Equation (28)) implies that h l . Proof: ut* vˆtb , together with the linear increase of ut* over time within a congested interval with slope l , that is, ut*1 ut* l for all but the last t , means vˆtb1 vˆtb l . Solving for leads to
24
ACCEPTED MANUSCRIPT
l . vˆ vˆtb b t 1
The denominator vˆtb1 vˆtb 1 due to the no-passing condition, hence l and h h l for h 0 .
□
Theorem 5 means that the attempt to compensate for the difference between the Lagrangian function (26) and the objective function (1) of OR_ITER by modifying the WIP holding cost
CR IP T
coefficient, even under the best condition ut* vˆtb , requires a ratio between WIP and FGI holding cost coefficient that cannot result from standard cost accounting logic and fundamentally changes the nature of the release model OR_ITER.
From this analysis we can conclude that in general the backward lead times over time cannot be assumed to be proportional to the dual prices of the clearing function constraints. Even if they are for
AN US
a certain interval, model OR_ITER with conventional parameter settings respecting cost accounting practice does not minimize the objective function that is required by theory. The iterative mechanism in general will not lead to an optimal solution even if the value of the clearing function is (nearly) independent of the load in the relevant region. Therefore, iterating over the lead times cannot be expected to coordinate the release planning and the lead time estimation model even in the special
M
case of high utilization where the maximum output is approximately fixed.
ED
5 Numerical experiments
The insights derived above are illustrated by numerical examples. In Example 1 the nonlinear model
models, given by
PT
OR_MON is evaluated using the clearing function derived in Missbauer (2002) for discrete-time
(33)
CE
1 X t [C k t - C 2 2Ck k 2 2C t 2k t t2 ] 2
with t Wt 1 Rt , C=950, k=200. In an M/M/1 system this parameter setting corresponds to a
AC
service rate of 4.75 orders per period (see Missbauer 2002). The holding cost coefficients are h=0.5, l=1. The initial values for WIP and FGI are zero, the planning horizon T=20 periods. Demand is given, the optimal releases Rt , production quantities X t , dual prices ut* and lead times vˆtb are given in Table 1. The model was solved using the NLP feature of LINGO 13.0, and the vˆtb are obtained from Equation (12). Ideally the iterative mechanism should converge to this solution.
25
ACCEPTED MANUSCRIPT
Rt
Xt
vˆtb
0 0 0 0 0 900 900 900 900 900 900 900 900 900 600 400 200 100 0 0
0.00 0.00 0.00 186.84 1085.57 982.67 966.55 960.42 957.33 955.51 954.34 953.53 952.95 289.75 0.00 0.00 107.88 70.20 0.00 0.00
0.00 0.00 0.00 149.49 662.36 739.71 776.29 798.67 814.16 825.70 834.72 842.02 848.10 808.78 600.00 400.00 200.00 100.00 0.00 0.00
0.000 0.000 0.000 0.200 0.424 0.716 0.925 1.098 1.251 1.390 1.517 1.635 1.746 1.898 2.268 2.500 0.494 0.335 0.000 0.000
ut* 0.000 0.000 0.605 0.648 1.648 2.648 3.648 4.648 5.648 6.648 7.648 8.648 9.648 11.213 0.000 0.708 0.669 1.894 0.000 0.000
AN US
CR IP T
Dt
M
t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Table 1: Optimal solution of model OR_MON
ED
The vˆtb are fractional lead times defined at the period boundaries. As shown in Appendix 1, the
PT
WIP holding costs are (from Equation (37)):
X T 1 X X t 1 h vˆtb t vˆTb T . 2 2 t 1
CE
Substituting this for the term h
W
t
(34)
in the objective function (1), omitting the WIP balance
t
AC
equations (2) (and hence the variables Rt ) and solving model OR_ITER with these modifications yields the production quantities X t shown in Table 2. Since demand never exceeds capacity, output equals demand in all periods which follows from Corollary 1. This solution is different from the optimum of OR_MON. Due to the fractional lead times the cumulative releases R cum (t ) (in continuous time) are calculated from X t and vˆtb by Equation (12) as depicted in Figure 2, the releases are Rt Rtcum Rtcum 1 . The actual (simulated) output that results from these releases Rt is obtained from the lead time estimation model that - consistent with model OR_MON - consists of the
26
ACCEPTED MANUSCRIPT
WIP balance equations (15) and the clearing function (33) (as an equality constraint). This also yields the lead times as above (Equation (12)).
Rt
Xt
0 0 0 0 0 900 900 900 900 900 900 900 900 900 600 400 200 100 0 0
0.00 0.00 0.00 0.00 538.98 1175.68 1089.47 1058.24 1039.69 1025.49 1015.31 1048.89 847.84 293.68 66.53 66.53 76.29 57.38 0.00 0.00
0.00 0.00 0.00 0.00 0.00 900.00 900.00 900.00 900.00 900.00 900.00 900.00 900.00 900.00 600.00 400.00 200.00 100.00 0.00 0.00
vˆ
0.000 0.000 0.000 0.200 0.424 0.716 0.925 1.098 1.251 1.390 1.517 1.635 1.746 1.898 2.268 2.500 0.494 0.335 0.000 0.000
Estimation model Xt vˆb t
0.00 0.00 0.00 0.00 396.01 713.90 779.55 810.54 829.20 841.83 851.08 859.94 859.44 828.46 737.16 518.70 246.67 103.16 20.05 3.56
AN US
Dt
M
t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
b t
0.000 0.000 0.000 0.000 0.265 0.514 0.840 1.096 1.315 1.511 1.688 1.848 2.001 2.211 2.392 2.366 0.920 0.425 0.000 0.000
CR IP T
OR_ITER
ED
Table 2: Optimal X t of model OR_ITER, the resulting Rt and the simulated X t and vˆtb These results demonstrate that the solution diverges from the optimum of OR_MON and the
PT
optimum is not a fixed point.
Table 1 also demonstrates that the dual prices ut* and the estimated backward lead times vˆtb over
CE
time are not proportional even in the congested interval (periods 4 to 13). This is to be expected since the first derivative of the clearing function does not satisfy the condition given in Equation (32).
AC
In Example 2 we demonstrate the effect of high utilization. The parameter k of the clearing function (33) is set to 30 which makes the clearing function rather flat (close to the capacity) for moderate to high WIP levels. Demand exceeds capacity over several consecutive periods which forces the model to operate in this high-utilization regime. All other data remain unchanged. The left-hand side of Table 3 shows the results of OR_MON. The right-had side presents the relevant results of model OR_ITER using ut* and vˆtb resulting from OR_MON. The objective function is the Lagrangian function (25) discarding the additive constant
uC t
t
and substituting (34) for the WIP holding cost
term as above.
27
ACCEPTED MANUSCRIPT
Rt
Xt
vˆtb
ut*
0 0.00 0 0.00 0 0.00 0 10.06 0 1041.93 1000 964.26 1000 956.90 1000 954.26 1000 952.96 1000 952.21 1000 951.73 1000 951.40 1000 951.17 900 753.12 600 111.43 400 370.39 200 186.18 100 95.53 0 0.00 0 0.00
0.00 0.00 0.00 9.75 831.58 865.93 881.26 890.43 896.70 901.33 904.93 907.84 910.24 900.00 600.00 400.00 200.00 100.00 0.00 0.00
0.000 0.000 0.000 0.031 0.202 0.320 0.402 0.470 0.530 0.583 0.633 0.679 0.722 0.717 0.462 0.059 0.043 0.037 0.000 0.000
0.000 0.000 0.000 0.516 1.516 2.516 3.516 4.516 5.516 6.516 7.516 8.516 9.516 6.200 0.616 0.547 0.525 1.559 0.000 0.000
Dt
0.00 0.00 0.00 950.00 950.00 950.00 950.00 950.00 950.00 950.00 950.00 400.00 0.00 900.00 600.00 400.00 300.00 0.00 0.00 0.00
AN US
t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
OR_ITER with Lagrangian objective function Capacity Xt Rt dual price 0.00 0.00 0.00 985.53 1132.05 1063.58 1028.61 1015.24 1006.68 1001.32 632.41 134.59 646.57 596.22 388.93 285.70 282.58 0.00 0.00 0.00
0.000 0.000 0.000 0.320 0.270 0.198 0.148 0.110 0.078 0.050 0.024 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
CR IP T
OR_MON
M
Table 3: Optimal solution for high utilization As expected the result for OR_ITER is remarkably different from the optimum. The dual prices of
ED
the capacity constraints (5) are close to zero even in highly utilized periods, X13 is zero at reduced cost of 0.0222 which indicates that the difference to alternative solutions is very small. Further numerical studies reveal that if demand oscillates between two values that are below the
PT
capacity, both ut* and vˆtb exhibit an alternating pattern that seems to depend on the parameter setting
CE
of the clearing function. Proportionality between ut* and vˆtb cannot be expected for these irregular demand patterns.
Conclusions and implications for future research
AC
6
This paper analyzed the theoretical basis of a subset of production planning algorithms that iterate between the planning and the scheduling level, namely iterative LP – simulation algorithms for order release planning. The formal analysis is limited to algorithms that iterate solely on the lead times (not on the capacities) and seeks to identify the cause of the convergence problems that are observed frequently. From this analysis we can draw the following conclusions: (1) Adjusting the planned lead times for order release planning by iterating between an order release and a lead time estimation
28
ACCEPTED MANUSCRIPT
model is a dual (price) coordination mechanism although no prices are defined explicitly. Except for special cases this procedure will not lead to the optimal solution since the optimal solution is not a fixed point. (2) There seems to be no straightforward way to overcome this problem within this iterative framework: Calculating the optimal solution to the monolithic model OR_MON by an iterative price coordination mechanism requires minimization of the Lagrangian function (24) which in turn requires estimating the maximum output in the periods for a given trajectory of system states
CR IP T
prior to the period under consideration, and in our context the output estimation is done by simulation. It is unclear whether a simulation model can be integrated into the Lagrangian function in a way that leads to a tractable optimization model. (3) Including capacity constraints in the release model can partly compensate for these deficiencies. Thus, formulating the release model without capacity constraints cannot be recommended although this has been plausibly argued in the literature.
The question is how, if at all, the iterative order release mechanisms can be improved based on
AN US
these insights. This question is beyond the frontier of today’s research, so we can only present some preliminary insights.
An alternative to improving the price (dual) coordination mechanism is to replace/complement it by primal coordination, that is, to limit/assign the amount of resources available. The insight from Corollary 1 that, given realistic values for WIP holding costs and lead times, the lead time updating
M
procedure does not change the planned production quantities strongly supports this way. The iterative order release mechanism in Kim and Kim (2001) that iterates both on the lead times and on the
ED
maximum output seem to exhibit a much better convergence behaviour compared to mechanisms that only iterate on the lead times (see Section 2), and virtually all iterative production planning – scheduling algorithms outside the order release domain iterate on some representation of capacity. It
PT
seems reasonable to assume that iteratively updating the available capacities in the release model acts as a primal coordination mechanism and thus can partly compensate for the weaknesses of the dual
CE
coordination mechanism that is implied by lead time updating. Working out the underlying theory is a topic for future research.
AC
Virtually all iterative order release mechanisms use specific outcomes of the scheduling level (flow times and/or capacity utilization) as feedback information to the release model. This is a remarkable difference from both fix-and-optimize heuristics that use the decision variables directly and simulation-based optimization of the releases where the feedback is derived from the objective function value, such as a gradient in gradient-based algorithms or a fitness value in genetic algorithms. Future research must clarify the optimal design of this feedback mechanism and, since the scheduling level is often represented by a simulation model, the optimal positioning of simulation in hybrid analytical/simulation methods that solve order release and other production planning problems.
29
ACCEPTED MANUSCRIPT
Appendix 1: Proof of Theorem 1 Period-based integer lead times If lead times are integer and defined for each production period as in model OR_ITER, WIP holding costs of h are incurred each time a unit of the product is carried over from one period to the next. Thus the total WIP holding costs are h multiplied by the number of times this happens within the planning
CR IP T
horizon (analogous to the “parts times periods” in lot sizing). Since a WIP inventory at the end of period t implies that Wt product units are carried over from period t to period t+1, the total number of “parts times periods” in the WIP inventory are
W
and the total WIP holding costs are h
t
t
W . t
t
In the optimal solution to model OR_ITER the final WIP WT is zero since there is no incentive to
AN US
release work that is not completed and WIP inventory is penalized in the objective function. Hence all “parts times periods” can be assigned to a certain production period of the respective product units. The number of ”parts times periods” that can be assigned to production in period t is vˆtb X t by the fixed backward lead time assumption, so the number of “parts times periods” is
M
immediately leads to
hWt h vˆtb X t
b t
X t , which
t
(35)
t
ED
t
vˆ
Epoch-based, fractional lead times
Hung and Leachman (1996) assume that the backward lead times are defined at the period boundaries
PT
and can take fractional values at any time. This assumption is different from model OR_ITER. In the following we analyze whether Theorem 1 which states that time-varying lead times have the same
CE
effect as time-dependent production costs extends to the case of fractional, epoch-based lead times assuming uniform output rate within the periods. We denote by Pt the WIP inventory area (parts times
AC
periods) of the work finished in period t. Figure 2 exhibits the cumulative release and output and the WIP inventory areas Pt for assumed values of the backward lead times. For exposition we assume W0 = 0 which implies v0b 0 .
30
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 2: Relationship between cumulative releases, cumulative output, Pt and vˆtb
M
for five periods with vˆtb = 0.75, 1.5, 2, 1.5, 1 for t=1 to 5.
ED
The WIP inventory areas Pt are trapezoids and can be calculated as:
Pt
1 b vˆt 1 vˆtb X t 2
for t=1,.., T
(36)
CE
PT
For all t 1,..., T 1 each vˆtb is multiplied by X t 2 and X t 1 2 . Therefore: T
T 1
t 1
t 1
Pt vˆtb
X t X t 1 X vˆTb T 2 2
(37)
Equation (37) means that for all periods except period T the WIP holding costs per unit of work
AC
( h vˆtb ) are charged to the average production quantity of the periods t and t+1. Hence the insight that the backward lead times act as prices for production is still valid for epoch-based, fractional lead times, but the assignment of the prices to production periods is slightly different. A positive initial WIP W0 0 would require a modification of the calculation of Pt for the first period since for this period the Pt area is no longer a trapezoid (see Figure 2). This does not provide any additional relevant insight here.
31
ACCEPTED MANUSCRIPT
Appendix 2: Multiple optima in model OR_MON with relaxed and dualized clearing function We denote *t and X t* the load t Wt 1 Rt and the output, respectively, in period t at the optimum of OR_MON. The multiple optima maintain the total output
X t
* t
and the values of *t ,
Wt X t *t ,
CR IP T
that is, (38)
which makes the clearing function term in the Lagrangian function (24) an additive constant for given
ut . Substituting Wt *t X t into (24) and eliminating the clearing function term leads to the objective function
h *t h X t l I t ut X t Min ! t
t
t
AN US
t
(39)
The first term is an additive constant for given *t . This also holds for the second term (total output) if production is only shifted between periods. Given these conditions the relevant objective function is just
l I t ut X t Min !
(40)
t
M
t
If, starting from feasible optimal output quantities, production is shifted to the previous period, that is,
ED
X t X t* 1 and X t 1 X t*1 1 for some t, the first term of (40) increases by l and the second term decreases by ut ut 1 . If this change is made in a congested interval the objective function value
PT
remains unchanged since the dual prices satisfy ut* ut*1 l (see Section 4.3). This change is locally feasible since the relevant capacity constraint (the clearing function) is relaxed and the resulting
CE
decrease of Wt 1 due to the unchanged load (Equation (38)) is feasible since the optimal Wt is
AC
positive in a congested interval (see Section 4.3).
32
ACCEPTED MANUSCRIPT
Appendix 3: Proportionality of ut* and t in model OR_MON The shape of the clearing function for which ut* and t are proportional over time within a congested interval can be derived as follows: From Equation (31) and the proportionality condition
t ut* a l t from which t can be obtained as
t a l
we obtain by substituting t from (41) into (31)
f '( t )
a
h . t a
(41)
CR IP T
t
(42)
AN US
The solution of this differential equation leads to the clearing function f ( t ) K h Log t
M
with K the integration constant.
(43)
Acknowledgments
ED
The author wishes to thank Reha Uzsoy for valuable comments on earlier versions of the paper.
Declarations of interest
AC
CE
PT
None
33
ACCEPTED MANUSCRIPT
References Albey, E., & Bilge, Ü. (2011). A hierarchical approach to FMS planning and control with simulationbased capacity anticipation. International Journal of Production Research, 49, 3319-3342. Albey, E., Bilge, U., & Uzsoy, R. (2012). Performance of an iterative production planning approach for production systems with manufacturing flexibility. Pre-Prints of the 17th International Working Seminar on Production Economics 4, Innsbruck, Austria, 1-12.
CR IP T
Asmundsson, J., Rardin, R. L., Turkseven, C. H., & Uzsoy, R. (2009). Production planning with resources subject to congestion. Naval Research Logistics, 56, 142-157.
Bahl, H. C., & Ritzman, L. P. (1984). An integrated model for master scheduling, lot sizing and capacity requirements planning. The Journal of the Operational Research Society, 35, 389-399. Bang, J.-Y., & Kim, Y.-D. (2010). Hierarchical production planning for semiconductor wafer
AN US
fabrication based on linear programming and discrete-event simulation. IEEE Transactions on Automation Science & Engineering, 7, 326-336.
Bertrand, J. W. M., Wortmann, J. C., & Wijngaard, J. (1990). Production Control: A Structural and Design Oriented Approach. Amsterdam: Elsevier.
Bilgen, B., & Çelebi, Y. (2013). Integrated production scheduling and distribution planning in dairy supply chain by hybrid modelling. Annals of Operations Research, 211, 55-82.
M
Boyd, S., & Vandenberghe, L. (2004). Convex Optimization. Cambridge, UK: Cambridge University Press.
ED
Byrne, M. D., & Bakir, M. A. (1999). Production planning using a hybrid simulation-analytical approach. International Journal of Production Economics, 59, 305-311. Byrne, M. D., & Hossain, M. M. (2005). Production planning: an improved hybrid approach.
PT
International Journal of Production Economics, 93-94, 225-229.
69.
CE
Carey, M. (1987). Optimal time-varying flows on congested networks. Operations Research,. 35, 58-
Carey, M. (1992). Nonconvexity of the dynamic traffic assignment problem. Transportation Research
AC
Part B: Methodological, 26, 127-133. Dauzère–Pérès, S., & Lasserre, J. (2002). On the importance of sequencing decisions in production planning and scheduling. International Transactions in Operational Research, 9, 779-793.
de Kok, A. G., & Fransoo, J. C. (2003). Planning supply chain operations: definition and comparison of planning concepts. In A. G. de Kok & S. C. Graves (Eds.), OR Handbook on Supply Chain Management. Amsterdam: Elsevier, 597-675. Dirickx, Y. M. I., & Jennergren, L. P. (1979). Systems Analysis by Multilevel Methods: With Applications to Economics and Management. Chichester et al.: John Wiley & Sons.
34
ACCEPTED MANUSCRIPT
Figueira, G., & Almada-Lobo, B. (2014). Hybrid simulation–optimization methods: A taxonomy and discussion. Simulation Modelling Practice & Theory, 46, 118-134. Graves, S. C. (1982). Using Langrangean techniques to solve hierarchical production planning problems. Management Science, 28, 260-275. Hackman, S. T. (2008). Production Economics. Berlin: Springer. Hung, Y. F., & Hou, M. C. (2001). A Production Planning Approach Based on Iterations of Linear
Industrial Engineers, 18, 55-67.
CR IP T
Programming Optimization and Flow Time Prediction. Journal of the Chinese Institute of
Hung, Y. F., & Leachman, R. C. (1996). A production planning methodology for semiconductor manufacturing based on iterative simulation and linear programming calculations. IEEE Transactions on Semiconductor Manufacturing, 9, 257-269.
Irdem, D. F., Kacar, N. B., & Uzsoy, R. (2008). An experimental study of an iterative simulation-
AN US
optimization algorithm for production planning. 2008 Winter Simulation Conference. IEEE, Austin. 2176-2184.
Irdem, D. F., Kacar, N. B., & Uzsoy, R. (2010). An exploratory analysis of two iterative linear programming - simulation approaches for production planning. IEEE Transactions on Semiconductor Manufacturing, 23, 442-455.
M
Karmarkar, U. S. (1989). Capacity loading and release planning with work-in-progress (WIP) and lead-times. Journal of Manufacturing and Operations Management, 2, 105-123.
ED
Kefeli, A., Uzsoy, R., Fathi, Y., & Kay, M. (2011). Using a mathematical programming model to examine the marginal price of capacitated resources. International Journal of Production Economics, 131, 383-391.
PT
Kim, B., & Kim, S. (2001). Extended model for a hybrid production planning approach. International Journal of Production Economics, 73, 165-173.
CE
Kim, S. H., & Lee, Y. H. (2016). Synchronized production planning and scheduling in semiconductor fabrication. Computers and Industrial Engineering, 96, 72-85.
AC
Kim, S., & Uzsoy, R. (2008). Integrated planning of production and engineering process improvement. IEEE Transactions on Semiconductor Manufacturing, 21, 390-398.
Leachman, R. C., & Carmon, T. F. (1992). On capacity modeling for production planning with alternative machine types. IIE Transactions, 24, 62-72.
Luenberger, D. G. (1989). Linear and nonlinear programming. Massachusetts: Addison-Wesley Reading. Mesarovic, M. D., Macko, D., & Takahara, Y. (1970). Theory of Hierarchical, Multilevel, Systems: New York: Academic Press.
35
ACCEPTED MANUSCRIPT
Milne, R. J., Mahapatra, S., & Wang, C.-T. (2015). Optimizing Planned Lead Times for Enhancing Performance of MRP Systems. International Journal of Production Economics, 167, 220-231. Missbauer, H. (2002). Aggregate order release planning for time-varying demand. International Journal of Production Research, 40, 688-718. Missbauer, H. (2009). Models of the transient behaviour of production units to optimize the aggregate material flow. International Journal of Production Economics, 118, 387-397.
CR IP T
Missbauer, H., & Uzsoy, R. (2011). Optimization models of production planning problems. In K. G. Kempf, P. Keskinocak, & R. Uzsoy (Eds.), Planning Production and Inventories in the Extended Enterprise: A State of the Art Handbook. Berlin: Springer, 437-508.
Modigliani, F., & Hohn, F. E. (1955). Production planning over time and the nature of the expectation and planning horizon. Econometrica, 23, 46-66.
of Technology, Eindhoven, Netherlands.
AN US
Negenman, E. G. (2000). Material Coordination under Capacity Constraints. Eindhoven University
O'Regan, D., Meehan, M., & Agarwal, R. P. (2001). Contractions. In Fixed Point Theory and Applications (pp. 1-11). Cambridge: Cambridge University Press.
Puergstaller, P., & Missbauer H. (2012). Rule-based vs. optimization-based order release in workload
Economics, 140, 670-680.
M
control: A simulation study of a MTO manufacturer. International Journal of Production
Riaño, G. (2003). Transient Behavior of Stochastic Networks: Application to Production Planning
ED
with Load-Dependent Lead Times. Ph.D. dissertation, Georgia Institute of Technology, School of Industrial and Systems Engineering, Atlanta, GA. Schneeweiss, C. (2003). Distributed Decision Making. Berlin: Springer.
PT
Shanthikumar, J. G., & Sargent, R. G. (1983). A unifying view of hybrid simulation/analytic models and modeling. Operations Research, 31, 1030-1052.
CE
Voss, S., & Woodruff, D. L. (2006). Introduction to Computational Optimization Models for Production Planning in a Supply Chain. (2nd ed.). Berlin: Springer.
AC
Wolosewicz, C., Dauzère-Pérès, S., & Aggoune, R. (2015). A Lagrangian heuristic for an integrated lot-sizing and fixed scheduling problem. European Journal of Operational Research, 244, 3-12.
Zäpfel, G. (1984). Systemanalytische Konzeption der Produktionsplanung und –steuerung für Betriebe der Fertigungsindustrie. In C. Zink (Ed.), Sozio-technologische Systemgestaltung als Zukunftsaufgabe. München (Munich): Carl Hanser Verlag (in German). Zhao, S., Wang, L., & Zheng, Y. (2014). Integrating production planning and maintenance: an iterative method. Industrial Management & Data Systems, 114, 162-182.
36
Order release planning by iterative simulation and linear programming: Theoretical foundation and analysis of its shortcomings Hubert Missbauer PII: DOI: Reference:
S0377-2217(19)30594-6 https://doi.org/10.1016/j.ejor.2019.07.030 EOR 15931
To appear in:
European Journal of Operational Research
Received date: Accepted date:
16 July 2018 10 July 2019
Please cite this article as: Hubert Missbauer , Order release planning by iterative simulation and linear programming: Theoretical foundation and analysis of its shortcomings, European Journal of Operational Research (2019), doi: https://doi.org/10.1016/j.ejor.2019.07.030
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights Iterating between production planning and scheduling level is a promising technique.
This has been applied especially for order release planning.
Iterative order release mechanisms often fail to converge.
The paper identifies the reason for this failure.
This is done by analyzing the theoretical underpinning of this iterative procedure.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Order release planning by iterative simulation and linear programming: Theoretical foundation and analysis of its shortcomings Hubert Missbauer
CR IP T
University of Innsbruck Department of Information Systems, Production and Logistics Management Universitaetsstrasse 15 6020 Innsbruck, Austria Email: [email protected]
Abstract
AN US
Planning future order releases for complex manufacturing systems with substantial lead times must consider limited capacities and the resulting production smoothing problem as well as the highly nonlinear relationship between capacity utilization, work-in-process and flow times. An important solution approach that has been proposed for this problem is to iterate between an order release model with fixed lead times and a lead time estimation (usually simulation) model that estimates the flow
M
times for given order releases, providing the lead times for the next iteration. However, the convergence of this iterative procedure is highly unpredictable, limiting its practical use.
ED
The iterative mechanism is analyzed analytically for simplified formulations of the order release and lead time estimation model. We show that an order release procedure of this type that iterates on the lead times is a dual (price) coordination mechanism whose design does not meet the theoretical
PT
requirements, and there is no straightforward way to overcome this. The analysis also provides insights into the results of several numerical studies from the literature and suggests a possible
CE
research direction to improve the method. Order release mechanisms of this type are a special case of a broader class of production planning methods that iterate between the production planning and the production scheduling level in
AC
order to provide realistic values for lead times and planned capacities. Providing theoretical underpinning for this type of production planning methods is an important research objective, and the paper pursues this direction for the special case of order release.
Keywords: Manufacturing; Production planning & control; Order release; Simulation
2
ACCEPTED MANUSCRIPT
1 Problem description Manufacturing planning and control (MPC) is a complex task encompassing decisions with different time horizons and impacts on the manufacturing or supply chain processes, which must include mechanisms that compensate for changes in planning data as well as unplanned events. Due to this inherent complexity the overall MPC problem is usually decomposed hierarchically into several
CR IP T
interdependent subproblems, such as master planning/master production scheduling, MRP, lot sizing, and detailed scheduling, including capacity planning at the respective levels.
We consider one specific decomposition, namely the structuring of the MPC system into the planning and the scheduling level. While the planning level usually determines quantities for production, transportation, sales, etc., for discrete planning periods, that is, it determines flows, the
AN US
scheduling level is event-oriented and determines assignments of operations to resources, sequencing of operations and their exact timing; hence planning and scheduling are very different in nature (Dauzère–Pérès and Lasserre, 2002).
Given this hierarchical decomposition, the planning level requires information on operational characteristics of the manufacturing system, most important the effective work center capacities (net
M
of idle times) and flow times or their probability distributions that are the basis for determining planned lead times. Since these operational characteristics are strongly influenced by scheduling decisions, the planning level must anticipate the behavior of the manufacturing system controlled by
ED
scheduling using an anticipation function as discussed by Schneeweiss (2003). Reliable estimates of effective capacities and flow times are necessary to determine appropriate values of the parameters
PT
describing capacities and planned lead times in the planning model. Determining planned lead times has been widely explored especially in the MRP literature (see
CE
Milne et al., 2015). Lead time setting must consider the high variance of flow time distributions in practice as well as the highly nonlinear relationship between mean flow times and resource utilization well-known in queueing theory. This complex behaviour results from the discrete and stochastic
AC
nature of the material flow and the autonomy of the scheduling level with respect to detailed scheduling decisions, and limits our ability to accurately estimate the lead times or, more generally, to identify feasible combinations of output quantity and lead times. It also means that lead times should be adjusted to capacity utilization, which requires time-dependent lead times if resource utilization varies over time due to time-varying demand. Iterative coordination of planning and scheduling is thus an obvious approach: The planning level uses a set of parameter values representing capacities and lead times to generate a production plan that is executed by a scheduling model that serves as an anticipation function. The resulting work center
3
ACCEPTED MANUSCRIPT
utilizations and/or flow times are used to update the capacity and/or lead time parameters in the planning model, which produces a revised production plan, until – hopefully – the algorithm converges, that is, the planned work center capacities and/or lead times used in the planning model and the actual utilization levels and/or lead times resulting from the scheduling model agree within a sufficiently small tolerance (we shall specify this more precisely later). An overview of production planning mechanisms of this kind is given in Section 2.
CR IP T
Most research on iterative planning and scheduling algorithms has been performed on order release planning. Order release is a decision problem that results from the interpretation of the planning – scheduling hierarchy described above as an organisational structure, which in many practical settings is a reasonable structure for the MPC system: The planning level coordinates the material flow over the entire logistic chain by setting coordinated production targets for production units that perform well-defined portions of the manufacturing process (e.g., a component
AN US
manufacturing shop or a manufacturing cell) and are autonomous with respect to detailed scheduling. The production targets set by the planning level are the work orders and their required due dates and consider adequate capacity constraints and planned lead times for the production units (for a thorough description of this decision structure, see de Kok and Fransoo, 2003; Bertrand, Wortmann, and Wijngaard, 1990). Detailed scheduling is performed locally within the production units. In this MPC
M
structure, order release, defined as the transfer of the control over the work orders under consideration from the planning to the scheduling level, is an important decision problem. Controlled order release
ED
keeps the work in process (WIP) level and – via Little’s law – the average flow time under control, while timely release of the orders should guarantee high due date performance. The flow time of a work order is the time elapsed from release until its completion, that is, the release time is the arrival
PT
time of the order at the first work center on its routing. Within this decision structure order release controls the arrival times of the work orders at the
CE
production unit which, given the behavior of the production unit that can be modelled as a queueing network, largely determines average flow time and output. Anticipating these effects of order release
AC
decisions is difficult, especially if the target values for WIP and flow times vary over time due to time-varying demand. Coordinating order release planning and detailed scheduling, the latter represented by a simulation or transient queueing model, by iterating between these levels and adjusting the planning parameters based on the simulation results is an important research direction. If simulation is used this can be considered as a hybrid analytical/simulation method (for reviews, see Shanthikumar and Sargent, 1983; Figueira and Almada-Lobo, 2014). The potential advantages of this approach are obvious: A very complex order release problem including stochastic elements can be solved by combining two well-known techniques (linear programming and simulation) using
4
ACCEPTED MANUSCRIPT
commercial software, and the techniques, especially in the case of simulation, can be explained to practitioners relatively easily. However, these iterative mechanisms lack any guarantee of optimality. Worse, in some cases they fail to converge and thus do not provide a feasible solution. This largely precludes practical application. Recently the performance of these iterative mechanisms under various conditions has been explored carefully by means of numerical studies, but there are virtually no theoretical insights
CR IP T
into how these mechanisms work and the nature of the underlying coordination mechanism. Therefore, no satisfactory explanation of the convergence problems has been provided. This reveals a major research gap: Solving the overall planning and scheduling problem by iteratively solving the planning and scheduling subproblems is a decomposition and coordination approach, but there has hardly been any attempt to integrate these iterative mechanisms into the theory of mathematical decomposition. Once the formal decomposition and coordination method implied by a specific
AN US
iterative planning/scheduling algorithm is identified, the relevant mathematical theory can, at least in principle, be applied to identify whether or not this iterative mechanism meets the requirements derived from that theory.
In the following this line of research is pursued for order release algorithms that iterate between an order release model with fixed lead times and a simulation or transient queueing model that
M
estimates the flow times.
The paper is organized as follows: In Section 2 the relevant literature is reviewed and the iterative
ED
order release mechanisms under study are described using Hung and Leachman (1996) as the main reference. Section 3 presents a simplified order release model encompassing only one work center as an instrument for the subsequent analysis. The possibilities of iteratively solving this model are
PT
studied in Section 4. Section 5 illustrates the insights by numerical experiments, while Section 6
2
CE
presents the conclusions and possible future research directions for this work.
Related literature
AC
Production planning methods that iterate between the planning and the scheduling level in order to determine correct values for parameters of the production planning model have been developed for various production planning problems. We first give a more general overview and then focus on iterative linear programming (LP) – simulation algorithms for order release planning. Dauzère-Pérès and Lasserre (2002) present an integrated model for lot sizing and scheduling and an algorithm that iterates between a lot sizing module that assumes a fixed production sequence and a scheduling module that sequences the given lots. Wolosewicz, Dauzère-Pérès, and Aggoune (2015) use Lagrangian relaxation to coordinate lot sizing and scheduling for a fixed sequence of operations
5
ACCEPTED MANUSCRIPT
on the machines. In both papers the feedback from the scheduling level represents the effective capacity available. Negenman (2000) presents an algorithm that iterates between an LP model that calculates the production plan for a production network and a flexible flow shop scheduling model that is solved by a heuristic. The feedback information provided by the scheduling level is the completion times of the orders, which are used to update the available capacities of the respective resources. Albey and Bilge (2011) present a hierarchical production planning and control system
CR IP T
framework for a Flexible Manufacturing System that consists of three levels: aggregate planning, loading and detailed planning. The behavior of the shop floor for a given production plan is anticipated by means of simulation. The simulation result is used to update the coefficients that represent capacities in the upper-level modules. Convergence is not analyzed in detail, but the authors indicate that capacity updating is complex due to the special problem and decision structure. Bilgen and Çelebi (2013) present an integrated production scheduling and distribution planning model where
AN US
the operation times are updated by iterating with a simulation model. Convergence of the operation times seems to be satisfactory. Since both preventive and corrective maintenance determines capacity losses, iteratively adjusting production planning and maintenance has been proposed in Zhao, Wang, and Zheng (2014). Note that all these algorithms essentially iterate on the capacities. Most iterative planning – scheduling algorithms have been developed for order release planning
M
where flow times are load-dependent and can be controlled by order release decisions. Multi-period models for order release planning usually represent the production unit as a network of work centers
ED
m=1,..., M. The planning horizon is divided into discrete planning periods t=1,..., T. The material flow is represented by inventory balance equations for work in process (WIP) at each work center and for final products, usually distinguishing different products or families of products with similar routing.
PT
The resulting network flow model is a fluid model in discrete time. Since the effect of order releases on lead times must be anticipated in order to determine optimal release quantities over time, the
CE
highly nonlinear relationships among WIP, flow time and throughput must be represented in the model. If fixed, exogenous lead times are assumed the model can be formulated as a linear program
AC
(LP) with the release quantities per period as the main decision variables. This model structure reflects the material flow structure in a variety of manufacturing systems for discrete manufacturing, especially in component manufacturing shops. Much of the research in this field was motivated by semiconductor manufacturing, and the simulation models used represent a wafer fab. The iterative order release algorithms analyzed in the following consist of an order release model (usually linear programming) that calculates optimal order releases per period assuming timedependent lead time parameters, and a lead time estimation (usually simulation) model. The simulation receives the release quantities as input and yields the predicted flow times and/or actual
6
ACCEPTED MANUSCRIPT
production quantities over time which are used to update the capacities and/or lead times of the release model for the next iteration. Zäpfel (1984) proposed an order release mechanism that iterates between an LP release model with fixed lead times and a simulation model for lead time estimation. The release model is formulated without capacity constraints since the effects of overloading (congestion) are included in the flow times that are the feedback from simulation. No numerical tests are provided.
CR IP T
The first to provide numerical tests of this type of order release planning were Hung and Leachman (1996), who developed an iterative release mechanism and provided tests using a model of a semiconductor wafer fab. Since we consider this as a reference work, the procedure is described in more detail.
Hung and Leachman (1996) use a release model similar to the step-separated formulation in Leachman and Carmon (1992). The lead times are fixed parameters and are defined as the time
AN US
required for a lot of product i to reach operation j after being released into the plant. The lead times are epoch-based, that is, they are defined for the period boundaries and thus can vary within the period. Hung and Leachman use backward lead times at the period boundaries, defined as the lead time of an order that is finished (or arrives at a certain operation) at the start/end of a specified period.
AC
CE
PT
ED
M
This is depicted in Figure 1.
Figure 1: Epoch-based backward lead times to relate output to wafer releases. Source: Hung and Leachman (1996, p. 260)
In the Hung & Leachman model the releases and backward lead times determine both the arrival times at the work centers and the periods in which capacity at the work centers is required for
7
ACCEPTED MANUSCRIPT
production; the lead times are not modeled as time windows within which the production periods can be chosen, as in de Kok and Fransoo (2003) and Pürgstaller and Missbauer (2012). The simulation model collects statistics on the release times of the manufacturing lots (equivalent to work orders in the above description) and the time an operation on a lot is initiated. The time elapsing between these two events gives the flow time of this operation, as defined in the stepseparated formulation. The flow times at the period boundaries are interpolated from the flow times of
CR IP T
lots that initiate processing immediately before and immediately after the boundary epoch (Hung and Leachman, 1996, p. 262).
Hung and Leachman (1996) examine the rate of convergence of the flow time estimates to the actual flow times in the simulation and find that convergence to the correct expected flow time values can be quite rapid, but that the procedure can fail to converge in some cases which are not fully understood. Irdem, Kacar, and Uzsoy (2008) provide detailed numerical tests and report that for high
AN US
utilization of the bottleneck facilities the results of the Hung/Leachman procedure “are
…..
discouraging in terms of convergence”. For the low utilization case the results are classified as “rather discouraging”, but “when the bottleneck utilization is low, the value of the iterative procedure is very limited in terms of its ability to obtain improved system performance.” Reducing the period length from seven days to one day did not lead to a fundamental change of this pattern (Irdem, Kacar, and
M
Uzsoy, 2010, p. 9); the authors state that “…LP model and simulation model do not represent each other sufficiently.”
ED
Riaño (2003) uses a similar iterative procedure where the LP release model expresses the lead times as weights wkt that denote the fraction of input released in period k that will emerge as output in
PT
period t. The wkt are estimated using a transient queuing model. The algorithm iterates on the weights rather than on the lead times as in Hung and Leachmann (1996). Riaño (2003) reports that the algorithm does not always converge and that convergence may be affected by the choice of the initial
CE
condition. “However, when it does converge, it converges to the same solution regardless of the initial solution used. When it did not converge, it would cycle through two or three solutions” (p. 83).
AC
Hung and Hou (2001) use the same basic procedure as Hung and Leachman (1996) but replace
the simulation model with a function relating expected flow times to the loading rate at individual work centers, similar to that used by Voss and Woodruff (2003). The loading rate is defined as the machine hours of all orders arriving at the work center in the respective period divided by the work center capacity, its functional relationship to the expected flow time is estimated by simulation. Their attempt to use an M/M/s queueing model proved unsuccessful. They report short computation times and good convergence for short sub-periods, but do not explicitly comment on whether there were convergence problems in some cases.
8
ACCEPTED MANUSCRIPT
Byrne and Bakir (1999) use a simulation model to determine realistic estimates for the available capacities, Byrne and Hossain (2005) present an extended production planning model within this framework. Bang and Kim (2010) and Kim and Lee (2016) present refined methods for semiconductor production planning that iterate with a priority-based scheduling or simulation model, respectively. Kim and Kim (2001) simultaneously update the lead times and available capacities assumed in the
CR IP T
LP model using this iteration scheme. Simulation is used to anticipate the flow times and the utilization of the work centers. As in Riaño (2003), the lead times are expressed as loading ratios, defined as the proportion of the start quantity released in a period p that contributes to the workload at machine k in period t (Kim and Kim, 2001, p. 168).
Kim and Kim (2001) do not report convergence problems in their numerical tests. Irdem, Kacar, and Uzsoy (2010) report good convergence of this procedure for both high and low utilization and
AN US
conclude “that the convergence behavior of the KK [Kim & Kim, 2001] procedure is qualitatively different from that of the HL [Hung & Leachman, 1996] procedure” (p. 452f.). They conjecture that this difference might be due to the differences in the estimation of the loading ratios which are estimated directly from the simulation in the Kim & Kim procedure whereas the Hung & Leachman procedure estimates the flow times of individual lots at specific epochs from the simulation, and uses
M
these to calculate an estimate of the loading ratios. (p. 453). However, Albey, Bilge, and Uzsoy (2012) perform numerical studies and find that this good convergence performance of the Kim & Kim
ED
procedure holds only for aggregate convergence criteria, such as objective function value and total production, but “achieving convergence at the level of production levels of each product in each period proves to be difficult” (p.11). They conclude that the ability of these approaches to allocate
PT
resources efficiently to products across multiple time periods is questionable. The body of research we have reviewed yields the following insights: (1) Convergence of iterative
CE
LP and lead time estimation procedures is a difficult issue that substantially hinders their practical implementation. However, including the capacities into the iterative mechanism as in Kim and Kim
AC
(2001) seems to mitigate the convergence problems substantially; some papers on iterative planning scheduling algorithms that essentially iterate on the capacities report satisfactory convergence. (2) The iterative methods for order release planning are not derived from a well-defined decomposition and coordination method. This is in contrast to several methods that iterate between two production planning levels, such as the Lagrangian methods presented in Graves (1982) or the fix-and-optimize approach of Bahl and Ritzman (1984). Since the underlying theory is unclear the convergence behaviour has been treated as a purely empirical question. (3) The convergence problems of iterative methods for order release planning that iterate solely on the lead times seem to be independent of specific details of the methods and their implementation. Riaño (2003) encountered convergence
9
ACCEPTED MANUSCRIPT
problems even for single-stage systems, so we expect that the complex material flow structure in the lead time estimation model (e.g., wafer fab with re-entrant material flows) is not the decisive factor. Convergence problems occur in order release models that include or discard WIP holding costs or capacity constraints, respectively. Likewise, iterating on the lead times (Hung and Leachman, 1996) and the loading ratios that are derived from the observed flow times (Riaño 2003) both produce this phenomenon.
CR IP T
Therefore, we hypothesize that the convergence problem is caused by the basic mechanism of this iterative approach, and inferring the theoretical basis from the proposed iterative technique should allow us to identify this basic mechanism. We pursue this idea in the following analysis.
3
The simplified order release model
AN US
In this section we formulate an order release model for a single production unit that retains the essential elements of the iterative mechanisms while remaining sufficiently tractable to obtain the necessary insights. We prove analytically that iterative order release mechanisms that iterate on the lead times imply a well-defined decomposition method, namely price coordination based on Lagrangian relaxation, but implement it in an incomplete manner that does not meet the theoretical
M
requirements especially for convergence. Therefore, they cannot be expected to converge to the optimum since the optimum is not a fixed point of the iterations. We argue that this behavior is inherent to the structure of the procedure rather than the complexity of real-life situations.
ED
The order release model for a production unit within the iterative approach is characterized by the following assumptions: The production unit consists of a single (bottleneck) work center that
PT
produces a single product. Time is divided into planning periods t=1,…., T. The orders, expressed as hours of work, are released and arrive at the WIP inventory, are processed and after processing arrive
CE
at the finished goods inventory. Demand for each period is known and deterministic. For each period t there is a fixed, integer backward lead time vˆtb ( vˆtb is an integer multiple of the period length) that denotes the lead time of products that are finished in period t. The maximum capacity per period,
AC
denoted by C, is constant over the planning horizon. Later we extend the model in order to capture the emergence of the lead times. The decision variables (all measured in hours of work) are the work released in period t, denoted
Rt , and the output (actual production) X t in period t, the WIP Wt at the end of period t and the finished goods inventory I t at the end of period t. The parameter Dt denotes the demand in period t (measured in hours of work), while h and l denote the holding costs per unit and period for WIP and final products, respectively. The initial conditions are represented by zero initial WIP ( W0 0 ) and a
10
ACCEPTED MANUSCRIPT
finished goods inventory level I 0 . We assume that these initial conditions, together with the capacity per period C, allow a feasible solution. The order release model OR_ITER can be formulated as follows:
h Wt l I t t
subject to
(1)
t t
t
k 1
k 1
Wt Rk X k t
t
k 1
k 1
I t I 0 X k Dk t vˆtb
t
R X Xt C
k 1
k
X t , Rt , I t ,Wt 0
(2)
for all t
(3)
for t vˆtb 1,..., T
(4)
for all t
(5)
AN US
k 1
k
for all t
CR IP T
min
for all t
(6)
The objective function (1) minimizes the costs for WIP and finished goods inventory (FGI). Equations (2) and (3) are the respective inventory balance equations. For exposition we assume zero initial WIP implying Rt 0 for t 0 (which can be relaxed by defining fixed releases for past
M
periods). Equation (4) states that for each period t the cumulative output at the end of period t must be
ED
equal to the cumulative releases at the end of period t vˆtb , implementing the specified backward lead times vˆtb . For given vˆtb as assumed here this is a set of linear equations that relate the Rt and the Xt.
PT
The vˆtb must satisfy a no-passing property, that is, work finished in period t+1 cannot be released earlier than work finished in period t (see Carey, 1992, for the no-passing conditions on the lead
CE
times). If vˆtb1 vˆtb 1 , the output in period t and in period t+1 are released in the same period t vˆtb . Therefore, the constraints are formulated as inequalities. The WIP holding costs in the objective function guarantee that the constraints hold at equality except in this special case, and that the
AC
specified backward lead times are met. (Alternatively, Equation (4) can be formulated using loading ratios as described above.) Equation (5) shows the capacity constraints. All variables are non-negative (Equation (6)). The model exhibits the essential structure of the order release model within the iterative mechanisms: In order to be finished at a certain work center in period t, the work is released in period
t vˆtb . The lead times can be achieved as long as the capacities are not overloaded. Production capacity is consumed at the end of the lead time. The lead times vˆtb applied in iteration k are derived
11
ACCEPTED MANUSCRIPT
from the flow times obtained from the simulation run in iteration k-1 and are parameters in the release model. For the first iteration appropriate initial values for the lead times must be determined. The flow times are deterministic since OR_ITER is a fluid model with first-come-first-served sequencing. Thus, the problem of deriving deterministic lead time parameters from observed flow time distributions does not occur here; this is actually an additional degree of freedom. Holding WIP inventory incurs costs. In the Hung & Leachman model the costs of releases and the revenues from
CR IP T
completions are discounted which leads to an increase in the objective function as WIP is held. The subsequent analysis is performed on this model. The iterative release mechanisms developed in the literature often are based on less restrictive assumptions and allow, e.g., multiple work centers, fractional lead times (Hung and Leachman 1996) or lead time distributions (Riaño 2003). Our reasoning is as follows: If it can be proven that even in the simplified case represented by model OR_ITER the iterative mechanisms do not lead to the optimum, this is sufficient to disprove
sufficient for falsification). We offer the following essential insight.
AN US
the general claim that the iterative mechanisms lead to the optimal solution (one counterexample is
Theorem 1: The order release model OR_ITER with given time-varying backward lead times vˆtb for orders (or work) behaves as an order release model with zero lead time and time-varying
M
production costs (that is, it yields the same optimal solution except time-shifted Rt ) as long as the
ED
constraint that the release periods t vˆtb must be positive is not binding. The WIP cost term in the objective function (1) is
h Wt h vˆtb X t .
(7)
t
PT
t
Hence the WIP holding costs for backward lead times vˆtb are equivalent to production costs of
CE
h vˆtb for each unit of work that is produced and occupies capacity in period t. Proof: See Appendix 1.
AC
The intuition of the proof is the following: Finishing an order in period t implies that its lead time is vˆtb , which is also the time it spends in WIP. Thus a unit of output completed in period t will incur a WIP holding cost of h times vˆtb . Since the WIP holding costs are assumed to be proportional to the production quantity they can be viewed as production costs; the impact is the same. Note that the planned release volumes Rt have no further influence on the objective function and, as long as there are no further constraints on the releases, are not relevant to the result. The optimal
Rt can be obtained by offsetting the optimal X t by the backward lead times vˆtb . This does not hold if 12
ACCEPTED MANUSCRIPT
release periods obtained in this way would be before period 1 since the release quantities for the past periods are fixed. Therefore, it is assumed that the initial conditions that are determined by the past releases are such that this case does not occur. Corollary 1: If value is added during production, which means h l , and production maintains the no-passing property, which implies vˆtb1 vˆtb 1 , then the optimal production quantities X t do not depend on the backward lead times vˆtb . The X t that are optimal for instantaneous production
CR IP T
( vˆtb 0 for all t) remain optimal for all vˆtb 0 that satisfy the no-passing property.
Proof: Model OR_ITER with vˆtb 0 for all t is an LP production planning model with instantaneous production that produces output as late as possible. Demand in some period t is produced in an earlier period s
AN US
production from period t to period s
vˆtb vˆsb
l (t s ) . h
(8)
M
Since vˆtb vˆsb t s due to the no-passing property and l h 1 due to the assumption that value is added, this condition can never be satisfied.
□
ED
Theorem 1 provides an important insight: The iterative mechanism coordinates the order release model and the lead time estimation model using a price coordination mechanism. Thus we can apply the theory of price coordination (dual coordination; see Dirickx and Jennergren, 1979) to evaluate the
PT
theoretical consistency of this iterative mechanism. We must check, firstly, whether the order release and the scheduling problem can be coordinated by this price coordination mechanism and, secondly,
CE
whether the iterative mechanism is a valid price adjustment procedure, that is, whether the lead times
vˆtb derived from the lead time estimation model lead to correct prices (steering costs) for optimizing
AC
the release and output volumes Rt and X t , respectively. This requires formulating the monolithic model that jointly determines the optimal release volumes and lead times and which is solved by means of the iterative mechanism. In order to formulate the monolithic model implied by the iterative mechanism the assumptions of model OR_ITER are complemented/replaced by the following assumptions:
The backward lead times vˆtb are decision variables that result from the releases Rt and the outputs X t of the relevant periods.
13
ACCEPTED MANUSCRIPT
The capacity constraints (5) are not sufficient to describe the possible output. The dynamics of the material flow (variability effects) can lead to idleness in a period t even though the planned available work (termed load) Wt 1 Rt is higher than the capacity. Hence there is an upper bound on the output that results from the dynamics of the material flow and can be calculated by the lead time estimation model. In the following we assume that the upper bound on the output is determined by the system dynamics prior
CR IP T
to the period under consideration (history of the process). (Including expectations about future periods, e.g., dispatching rules that anticipate future events, is not considered here.) We do not specify the logic by which the (deterministic) upper bound to the output is derived from the stochasticity of the material flow. It can be the conditional expectation or some quantile of the conditional output distribution given a certain history of the
AN US
process. This is analogous to estimating clearing functions from simulation data which can be done by ordinary least squares or quantile regression (see Asmundsson, Rardin, Turkseven, and Uzsoy (2009)).
Including the simulation or transient queueing model for output and lead time estimation into the order release model renders the resulting “monolithic” model unsolvable by mathematical programming methods. The reasoning for this is given below. Therefore,
The model is modified as follows:
M
iterative methods are applied for solving the model.
ED
Upper bound on the output
The following constraints are added:
for all t
(9)
PT
X t ( R [ 1,..., t ], X [ 1,..., t 1])
The output in a period t is a function of the releases and output in the previous periods and of the
CE
releases in period t. The function represents all effects of the dynamics, including the stochasticity, of the material flow within the production unit and is expressed as a simulation or transient queueing model that predicts the output in a period t from the releases until period t and the history of the
AC
output quantities. At this point we do not specify the mathematical properties of this function. The constraint is formulated as an equality constraint since simulation models usually do not keep resources idle if work is available. Deduction of the backward lead times from WIP and output For given releases Rt and outputs X t the lead times vˆtb can be calculated from Equation (4) provided they are integer. Non-integer lead times can be calculated as follows (for a comprehensive treatment, see Hackman, 2008, p. 320 ff.): We define X tcum , Rtcum as the cumulative output and release
14
ACCEPTED MANUSCRIPT
quantities, respectively, at the end of period t. X cum (t ), R cum (t ) are defined analogously in between the period boundaries, assuming uniform releases and output within the period: t
Rtcum Rk
for t 1,2,..., T
(10)
for 0 t T and t non-integer
(11)
k 1
Rcum (t ) Rcum Rt 1 (t t ) t
CR IP T
with t the integer part of the real-valued number t (largest integer less than or equal to t). The definitions of X tcum and X cum (t ) are analogous. The backward lead time vˆtb for orders that are finished at the end of period t is defined by
for t 1,2,..., T
W0 Rcum (t vˆtb ) X tcum .
(12)
In the following the analysis is performed assuming integer lead times.
AN US
Reformulation of the objective function
The objective function (1) can (but need not) be reformulated by substituting (7) for the WIP cost term, which yields:
h vˆtb X t l I t .
min
t
t
M
The monolithic model OR_MON that is solved by the iterative mechanisms is then: Objective function:
h Wt l I t
or, equivalently:
t
ED
min
(13)
t
h vˆtb X t l I t
PT
min
t
(14)
t
AC
CE
Inventory balance equations
t
t
k 1
k 1
Wt Rk X k t
t
k 1
k 1
I t I 0 X k Dk
for all t
(15)
for all t
(16)
for all t
(17)
for all t
(18)
Relationship between release and production volumes t vˆtb
t
R X k 1
k
k 1
k
Output anticipation
X t ( R [ 1,..., t ], X [ 1,..., t 1])
Capacity constraints (redundant if Equation (18) is correctly specified)
15
ACCEPTED MANUSCRIPT
Xt C
for all t
(19)
X t , Rt , I t ,Wt , vˆtb 0
for all t
(20)
If the objective function (14) is used, the mechanics of the model OR_MON are as follows: The release and output variables (Rt and Xt, resp.) lead to a certain WIP level (Equation (15)) and – by Equation (17) – to certain lead times that, in turn, determine the first term in the objective function. Note that in this lead time formulation of the model, Equation (17) is an integral part of the model. In
CR IP T
the general case this formulation cannot be solved by mathematical programming techniques for two reasons: Firstly, the function is implicit in a simulation or transient queueing model and exhibits a complex structure even if approximations are used (e.g., integer and nonlinear as shown in Missbauer, 2009). Secondly, the lead time estimation constraints (17) contain variables that are indices of other variables ( vˆtb in the index of Rk), which renders the model unsolvable with mathematical
AN US
programming methods. This structure occurs because the lead times are modeled correctly as a result of WIP and output (see Hackman, 2008, p. 338 ff.). If lead times are modeled without explicit reference to WIP (as in Voss and Woodruff, 2006) this difficulty can be avoided at the price of imprecision in modeling the dynamics of the lead times.
The second difficulty described above (variables as indices of other variables) can be avoided
M
provided that there are no cut-off effects at the start of the planning horizon; that is, extending the planning horizon infinitely into the past does not change the optimal solution. In this case (which is
ED
assumed in the following) there are no relevant constraints on the lead times. Using the original formulation of the objective function (Equation (13)) changes the mechanics of the model, yielding the WIP-oriented formulation: The release and output variables (Rt and Xt, resp.) lead to a certain WIP
PT
level (see Equation (15)) that determines the first term in the objective function (13). The lead time variables are not required for the optimization, but are deduced from the optimal solution by Equation
CE
(17), which poses no problem. Of course the difficulty with the structure of Equation (18) remains. The lead time oriented and WIP-oriented formulation of the model are equivalent although they
AC
exhibit different structures for the optimization. Hence in the following we use the WIP-oriented formulation which is easier to handle and deduct the lead times if necessary.
16
ACCEPTED MANUSCRIPT
4
Analysis of the iterative solution of the model
4.1 Method of the analysis We can now ask under what conditions the model OR_ITER, together with a lead time updating
CR IP T
procedure, can achieve the optimal solution to the monolithic model OR_MON. A necessary condition for this is that there exists a fixed point of the iterative mechanism and the optimal solution is a fixed point. A fixed point of the iterative mechanism is a solution where the time-dependent lead time vector that leads to a certain order release plan is equal to the lead time vector that results from these order releases and from the subsequent lead time estimation (e.g., simulation). Thus if the iterative mechanism reaches a fixed point the solution will remain there. Defining the desired property
AN US
of the optimal solution as being a fixed point avoids ambiguity in the definition of convergence: If an iterative order release mechanism converges in the objective function value, it need not converge in terms of the production quantities of the products (see Section 2).
Per the Banach Fixed Point Theorem (O'Regan, Meehan, and Agarwal, 2001) convergence to
that yields the lead time vector (from the lead
M
a fixed point is guaranteed when the function g vˆ b
time estimation, e.g., simulation model) as a function of the lead time vector vˆ b used in the order
ED
release model is a contraction mapping. Due to the complexity of the function g (.) that consists of an LP order release and a lead time estimation model we do not pursue this direction. Instead we analyze whether the optimal solution to the monolithic model OR_MON is a fixed point of the iterative
PT
mechanism that consists of model OR_ITER and an appropriate lead time estimation model. This is
the case if g vˆ b * vˆ b * with vˆ b * the optimal lead time vector, that is, when the release model
CE
with optimal lead times as parameters leads to a release plan that, executed by the lead time estimation model, again yields the lead time vector vˆ b * .
AC
This implies the two questions outlined in Section 3: The first of these is whether there exist lead
time parameters in the release model OR_ITER that lead to the optimal solution of the monolithic model OR_MON; this refers to the coordinability (Mesarovic, Macko, and Takahara 1970) of order release and lead time estimation problem by this iterative mechanism. The second is under what conditions – if at all – the lead time estimation model can yield these lead time values, which refers to the validity of the price adjustment procedure. We start with the first issue. For the following analyses we assume that the constraint that the release periods t vˆtb must be positive is not binding.
17
ACCEPTED MANUSCRIPT
In order to make model OR_MON solvable we consider only a special case of the output constraint (18): The output in period t is assumed to depend only on the load (available work) in period t ( Wt 1 Rt ) (see Equation (21)). Again we argue that if it can be proven that even in this special case the iterative mechanism does not lead to the optimum, this disproves the general claim that the iterative mechanism leads to the optimal solution. Furthermore, it is formulated as an inequality constraint such that the load determines the maximum output.
CR IP T
Thus we obtain a conventional (one-dimensional) clearing function (see Equation (22)):
( R [ 1,..., t ], X [ 1,..., t 1]) f (Wt 1 Rt )
(21)
X t f (Wt 1 Rt ) 0
(22)
for all t
Clearing function models with a concave, saturating clearing function were introduced by Karmarkar (1989) and have been described extensively in the literature (e.g., Missbauer and Uzsoy,
AN US
2011). The feasible set is convex (see Carey, 1987, which analyzes a model with this structure for traffic flow control), and the model can be solved either by nonlinear programming or by linearizing the clearing function and solving the resulting linear program. Carey (1987) proves that if certain assumptions are satisfied, the output is always at its upper bound, which in our context means that the clearing function constraints (22) hold at equality. One of these assumptions, termed the cost function
M
assumption, requires that “the marginal cost per unit volume, per unit time should not increase on moving from one arc to the next” (p. 63f.). This assumption is not satisfied here because the FGI holding costs l usually are higher than the WIP holding costs h. However, in our single-stage model
ED
the only purpose of the WIP level is to enable the desired output, and X t f (Wt 1 Rt ) 0 only occurs if demand falls very sharply within a short period of time (Kefeli, Uzsoy, Fathi, & Kay, 2011,
PT
p. 387). This issue will be considered below when necessary. In the course of the analysis we will also benefit from some structural insights into this type of models that are available in the literature. We
CE
use the WIP-oriented formulation of the objective function; the lead time anticipation equation (17) is not required for the optimization, it just calculates the lead times from WIP and output.
AC
4.2 Analysis for arbitrary utilization The iterative release mechanism anticipates the dynamics of the production unit external to the
optimization model. Hence the output constraint (22) is not a part of the optimization model OR_ITER; capacity is only constrained by Equation (19) or (5), respectively. As stated in Theorem 1, the lead times define the prices for producing items (or for occupying capacity, respectively) in certain periods. Thus in model OR_ITER the relevant capacity constraints of the monolithic model (Equation (22); note that the capacity constraints (19) are redundant) are relaxed and replaced by
18
ACCEPTED MANUSCRIPT
defining prices for the capacity. This is a Lagrangian relaxation of the constraints (22), and the optimality conditions for the prices (and hence for the lead times) can be derived from standard theory. The Lagrangian function obtained by relaxing the clearing function constraints (22) of model OR_MON is
£( X t , ut ) h Wt l I t ut X t f (Wt 1 Rt ) t
(23)
t
or
CR IP T
t
£( X t , ut ) h Wt l I t ut X t ut f (Wt 1 Rt ) , t
t
t
t
(24)
where ut denotes the dual variables associated with the clearing function constraints (22) in period t.
AN US
The Lagrangian function (24) is the objective function of model OR_ITER (Equation (1)) extended by a price term for occupying capacity and by the term
u f (W t
t 1
Rt ) . Model
t
OR_MON with the simplified output constraints (22) satisfies Slater’s condition since in general there exists at least one solution with positive decision variables where the clearing function constraints
M
hold with strict inequalities (Boyd and Vandenberghe, 2009, p. 226), and due to its convexity the duality gap is zero. Maximizing Equation (24) over ut and minimizing it over the primal variables
ED
(subject to the remaining constraints of model OR_MON) leads to the optimal ut (partial duality; see Luenberger, 1989, p. 401). The implications of this structure for the validity of the iterative order release mechanism are given in the following theorems.
PT
Theorem 2: Minimizing the Lagrangian function (24) over the primal variables for given optimal dual variables (dual prices) ut* subject to the remaining constraints (15), (16), (19), (20) of model
CE
OR_MON need not lead to a feasible primal solution that satisfies the Karush-Kuhn-Tucker (KKT) optimality conditions due to multiple optima. Thus, in the general case no choice of the ut will lead to
AC
the optimal solution of model OR_MON. Proof: The structure of the multiple optima is shown in Appendix 2. (For the KKT conditions for
an extended clearing function model, see Kim and Uzsoy, 2008.)
□
Theorem 3: Solving model OR_ITER to optimality in general will not yield the optimal solution of model OR_MON for any choice of the lead times vˆtb . Proof: As shown above, model OR_ITER implies that the clearing function constraints are relaxed. In this case a solution to model OR_MON with optimal objective function value can only be
19
ACCEPTED MANUSCRIPT
obtained by minimizing the Lagrangian function (24) for optimal dual prices ut* . (24) is a nonlinear function due to the clearing function related term
u f (W t
t 1
Rt ) . The objective function (1) of
t
model OR_ITER does not include this term and is linear independent of the lead times vˆtb , that is, □
model OR_ITER does not minimize the Lagrangian function (24).
The intuition of these theorems is the following: Per Theorem 1 the lead times act as prices for
CR IP T
production, and it has been shown that this implies a price coordination mechanism that coordinates the order release and the lead time estimation model. Our mathematical analysis of the integrated model OR_MON and its Lagrangian decomposition shows that a necessary condition for this price coordination mechanism to work is that the Lagrangian function (24) is minimized (given the dual prices ut* ). However, model OR_ITER cannot do this because it minimizes the wrong function as
AN US
shown in Theorem 3. In addition, even if the Lagrangian function (24) is minimized for optimal dual prices the solution can be infeasible per Theorem 2.
There is no obvious way to overcome this fundamental problem since in general f (Wt 1 Rt ) represents a simulation model in practical applications, and integrating
u f (W t
t 1
Rt ) into
t
M
Equation (1) would require integrating a function coded in a simulation model into the mathematical programming model OR_ITER. Note that this structural problem of model OR_ITER – discarding the
ED
respective price term in objective function (1) – is not specifically related to clearing functions; it also holds for the general formulation of the output constraint (Equation (18)) although in this case the structure of the solution space is unknown unless function is specified.
PT
These insights lead to the following
Conclusion: It cannot be expected that any choice of the backward lead times vˆtb in model
CE
OR_ITER can lead to the optimal solution to model OR_MON. Thus, the optimal solution to model OR_ITER cannot be expected to be a fixed point.
AC
This does not exclude the possibility that there exists a fixed point under specific assumptions. Especially, if the clearing function and/or the load situation allow independence between order releases and lead times there is a fixed point. This is formulated in *
Theorem 4: If the lead times are independent of the order release quantities (vˆtb ) , then there
exists a fixed point and the algorithm will converge to this point. Proof: If the system is initialized with an arbitrary lead time vector vˆtb for all t, the resulting Rt(1) *
vector will lead to the lead time vector vˆtb by assumption. This leads to a second vector of order
20
ACCEPTED MANUSCRIPT
*
*
release quantities ( Rt(2) ) . Rt(2) and vˆtb , t=1, …. T, are a fixed point because (vˆtb ) can never change (by assumption) and leads to vector Rt(2) (by definition of Rt(2) ).
□
Theorem 4 refers to a situation where congestion effects due to limited capacities are not relevant. This cannot be assumed in general for capacitated production resources, but it can be a good approximation if capacity utilization is low.
CR IP T
However, it can be argued that, despite these structural deficiencies, model OR_ITER might lead to a reasonable heuristic solution to model OR_MON in certain cases provided that the prices for production implied by the objective function (Equation (1)) h vˆtb are equal (or sufficiently close) to the optimal prices that result from the Lagrangian relaxation of the clearing function constraint (22) in model OR_MON. This might be the case for high utilization, which is analyzed next.
AN US
4.3 The special case of high utilization
Although normally the clearing function depends on the values of the decision variables, the optimal solution might be in a region where the maximum output is close to the capacity, that is,
f (Wt 1 Rt ) C especially for bottleneck work centers that are of primary interest. In this case the
M
clearing function related term in Equation (24) is approximately an additive constant. Furthermore, the relaxed constraints (22) can be violated only to a minor extent due to the capacity constraints (19).
ED
We now analyze whether in this case the iterative release mechanism can lead to backward lead times
vˆtb that have the same impact as optimal dual variables (dual prices) ut* of the constraints (22). Since there is no duality gap and the possible infeasibility may be very small, this might lead to a reasonable
PT
solution to model OR_MON.
CE
For f (Wt 1 Rt ) C the Lagrangian function (24) is approximately:
AC
Substituting h
£( X t , ut ) h Wt l I t ut X t utC . t
W hvˆ X b t
t
t
t
t
t
(25)
t
(from Equation (7)) and collecting the Xt coefficients we get
t
£( X t , ut ) hvˆtb ut X t l I t utC , t
t
(26)
t
where the last term is an additive constant. Obviously, if correct lead times are assumed the coefficients of the Xt in the objective function of model OR_ITER ((1) or (14), respectively) are too small; the dual variable ut is added to the correct WIP holding cost rate hvˆtb in the Lagrangian function (26). However, this need not preclude a reasonable solution: The WIP holding cost
21
ACCEPTED MANUSCRIPT
coefficient h is difficult to estimate; production prior to demand triggers various processes in the company with very different cost structures. Proportionality of holding costs to quantity and time is a strong assumption, and holding cost coefficients can be argued to be theoretical constructs rather than observable values. Therefore, it is worth exploring under which assumptions an iterative procedure that implicitly assumes ut* 0 (as in objective function (1) of model OR_ITER) can lead to an
hvˆtb ut* . h vˆtb If ut* vˆtb this simplifies to
h h .
CR IP T
optimal solution for some modified holding cost coefficient h . This would require hvˆtb hvˆtb ut* or (27)
(28)
AN US
Hence if the backward lead times and the dual prices of the clearing function constraints are proportional, then ignoring the respective price terms in the optimization, as done in model OR_ITER, can be compensated by increasing the holding cost coefficient for WIP by a specified value. In other words, using a modified holding cost coefficient h h in the objective function (1) is consistent with the Lagrangian function (26) if h exceeds the actual value h by the value (from Equation (28)).
M
With this modification the Lagrangian function (26) can be minimized by minimizing objective function (1).
ED
In the following we analyze under what assumptions the condition ut* vˆtb that provides the possibility to compensate the structural deficiencies of objective function (1) by increasing the WIP
PT
holding costs is satisfied. Since we still analyse the special case f (Wt 1 Rt ) C , the workload in period t t Wt 1 Rt
CE
calculations.
is used as a proxy for the lead time vˆtb in the subsequent analytical
Kefeli, Uzsoy, Fathi, and Kay (2011) analyze order release models with a saturating clearing function approximated by a set of linear functions (tangents), distinguishing congested periods
AC
(Wt>0), non-congested periods (Wt=0 and Xt>0) and idle periods (Wt=0 and Xt=0). They find that in a congested interval, as long as work is released and FGI is carried in consecutive periods, the sum of the nonzero dual prices of capacity (more precisely: of the linear segments of the clearing function that are tight in a certain period) increases by l (the holding cost coefficient of FGI) for each consecutive time period the system holds FGI (p. 388). Since shifting the segments of the linearized clearing function that are tight upward (that is, increasing the intercept of both segments by 1) increases capacity by 1, this theorem holds for the dual price of the clearing function constraint. The proof extends trivially to nonlinear, saturating clearing functions by letting the number of tangents
22
ACCEPTED MANUSCRIPT
approach infinity and the difference of their slopes approach zero. The intuition behind this is the same as in production planning models with zero lead time: Increasing capacity by one unit in a period t of an interval where the capacity constraints are tight and FGI is held reduces FGI by the respective amount until period t. Hence the later the additional capacity is added, the larger the savings in FGI holding costs. This relationship was recognized very early from marginal cost analysis (Modigliani & Hohn, 1955).
CR IP T
Since in clearing function models additional capacity can be generated by higher WIP, there is a relationship with the WIP level: The independent variable of the clearing function, in our case the load, is only increased if the marginal benefit of the additional capacity that is generated (determined by the first derivative of the clearing function) is higher than the marginal cost of the higher WIP. Kim and Uzsoy (2008) analyze a nonlinear, concave (saturating) clearing function with Wt as
AN US
explanatory variable and obtain the KKT condition (Equation (11) on p. 394; we omit the possibility of engineering lots here)
h t 1 t
df ut dWt
(29)
with t the dual variable of the WIP balance constraint in period t. They show that for Wt>0 and
M
Rt>0 Equation (29) holds at equality and the term t 1 t vanishes. Assuming this and using the
ED
load t Wt 1 Rt as explanatory variable as in our model, this leads to
df ( t ) h . d ( t ) ut*
(30)
PT
Defining f '( t ) df ( t ) d ( t ) and a to be the intercept of the linear increase of the ut* as a function of t within the congested interval under consideration, that is, ut* a l t , we can write
CE
Equation (30) as:
f '(t )
h . a l t
(31)
AC
Since the clearing function is concave and strictly monotone in t , it has a well-defined inverse
function f '( 1) . Therefore, from Equation (31) we obtain
h t f '( 1) a l t
(32)
This can be interpreted as follows: For a given intercept a (which is determined by the initial conditions of the congested interval under consideration) the argument of f '( 1) in Equation (32) is determined for all periods of the congested interval and is decreasing over time (since l>0). The
23
ACCEPTED MANUSCRIPT
function f reflects the operational characteristics of the work center and is also given. Hence the load
t for all periods of the congested interval under consideration is given by Equation (32). Although we can assume that the function f is concave and saturating, the first derivative f '(t ) strongly depends on the characteristics of the arrival and departure process at the work center and on the ratio of period length and average service time. Hence we cannot draw any general conclusion on the
CR IP T
evolution of the relationship between ut* and t over time within a congested interval beyond the monotone increase of t (and the linear increase of ut* ) with t. In particular, we cannot assume that
ut* and t increase proportionally over time. Numerical examples lead to the following
Observation 1: Even within a congested interval the values of ut* and backward lead times vˆtb over time in model OR_MON can be far from proportional depending on the shape of the clearing
AN US
function. For general situations where congested and non-congested periods occur alternately there is no correspondence between ut* and vˆtb over time; in this case not even monotonicity can be assumed. For this situation to the best of our knowledge no analytical results for the dynamic behaviour of the dual prices are available.
The numerical experiments are demonstrated in Section 5 below.
M
However, Equations (30) and (31), respectively, indicate that there should be a functional form of the clearing function that leads to a linear increase of the load t over time within a congested
ED
interval and thus to a proportional relationship of load, which we use as a proxy for the lead time vˆtb , and the dual prices within a congested interval. This is the case if the clearing function exhibits the
PT
logarithmic shape f ( t ) K h Log t with K an integration constant. This is shown in Appendix 3. There is no justification for this functional form in queueing theory, and it is inconsistent
CE
with common assumptions about service processes since it predicts negative output if the load is below a certain positive value. Furthermore, it is not consistent with the condition f (Wt 1 Rt ) C
AC
for high load which provided the justification for this analysis. The above analysis reveals a second limitation formulated below. Theorem 5: We consider a congested interval with ut* vˆtb for all t . For h 0 the
condition h h (Equation (28)) implies that h l . Proof: ut* vˆtb , together with the linear increase of ut* over time within a congested interval with slope l , that is, ut*1 ut* l for all but the last t , means vˆtb1 vˆtb l . Solving for leads to
24
ACCEPTED MANUSCRIPT
l . vˆ vˆtb b t 1
The denominator vˆtb1 vˆtb 1 due to the no-passing condition, hence l and h h l for h 0 .
□
Theorem 5 means that the attempt to compensate for the difference between the Lagrangian function (26) and the objective function (1) of OR_ITER by modifying the WIP holding cost
CR IP T
coefficient, even under the best condition ut* vˆtb , requires a ratio between WIP and FGI holding cost coefficient that cannot result from standard cost accounting logic and fundamentally changes the nature of the release model OR_ITER.
From this analysis we can conclude that in general the backward lead times over time cannot be assumed to be proportional to the dual prices of the clearing function constraints. Even if they are for
AN US
a certain interval, model OR_ITER with conventional parameter settings respecting cost accounting practice does not minimize the objective function that is required by theory. The iterative mechanism in general will not lead to an optimal solution even if the value of the clearing function is (nearly) independent of the load in the relevant region. Therefore, iterating over the lead times cannot be expected to coordinate the release planning and the lead time estimation model even in the special
M
case of high utilization where the maximum output is approximately fixed.
ED
5 Numerical experiments
The insights derived above are illustrated by numerical examples. In Example 1 the nonlinear model
models, given by
PT
OR_MON is evaluated using the clearing function derived in Missbauer (2002) for discrete-time
(33)
CE
1 X t [C k t - C 2 2Ck k 2 2C t 2k t t2 ] 2
with t Wt 1 Rt , C=950, k=200. In an M/M/1 system this parameter setting corresponds to a
AC
service rate of 4.75 orders per period (see Missbauer 2002). The holding cost coefficients are h=0.5, l=1. The initial values for WIP and FGI are zero, the planning horizon T=20 periods. Demand is given, the optimal releases Rt , production quantities X t , dual prices ut* and lead times vˆtb are given in Table 1. The model was solved using the NLP feature of LINGO 13.0, and the vˆtb are obtained from Equation (12). Ideally the iterative mechanism should converge to this solution.
25
ACCEPTED MANUSCRIPT
Rt
Xt
vˆtb
0 0 0 0 0 900 900 900 900 900 900 900 900 900 600 400 200 100 0 0
0.00 0.00 0.00 186.84 1085.57 982.67 966.55 960.42 957.33 955.51 954.34 953.53 952.95 289.75 0.00 0.00 107.88 70.20 0.00 0.00
0.00 0.00 0.00 149.49 662.36 739.71 776.29 798.67 814.16 825.70 834.72 842.02 848.10 808.78 600.00 400.00 200.00 100.00 0.00 0.00
0.000 0.000 0.000 0.200 0.424 0.716 0.925 1.098 1.251 1.390 1.517 1.635 1.746 1.898 2.268 2.500 0.494 0.335 0.000 0.000
ut* 0.000 0.000 0.605 0.648 1.648 2.648 3.648 4.648 5.648 6.648 7.648 8.648 9.648 11.213 0.000 0.708 0.669 1.894 0.000 0.000
AN US
CR IP T
Dt
M
t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Table 1: Optimal solution of model OR_MON
ED
The vˆtb are fractional lead times defined at the period boundaries. As shown in Appendix 1, the
PT
WIP holding costs are (from Equation (37)):
X T 1 X X t 1 h vˆtb t vˆTb T . 2 2 t 1
CE
Substituting this for the term h
W
t
(34)
in the objective function (1), omitting the WIP balance
t
AC
equations (2) (and hence the variables Rt ) and solving model OR_ITER with these modifications yields the production quantities X t shown in Table 2. Since demand never exceeds capacity, output equals demand in all periods which follows from Corollary 1. This solution is different from the optimum of OR_MON. Due to the fractional lead times the cumulative releases R cum (t ) (in continuous time) are calculated from X t and vˆtb by Equation (12) as depicted in Figure 2, the releases are Rt Rtcum Rtcum 1 . The actual (simulated) output that results from these releases Rt is obtained from the lead time estimation model that - consistent with model OR_MON - consists of the
26
ACCEPTED MANUSCRIPT
WIP balance equations (15) and the clearing function (33) (as an equality constraint). This also yields the lead times as above (Equation (12)).
Rt
Xt
0 0 0 0 0 900 900 900 900 900 900 900 900 900 600 400 200 100 0 0
0.00 0.00 0.00 0.00 538.98 1175.68 1089.47 1058.24 1039.69 1025.49 1015.31 1048.89 847.84 293.68 66.53 66.53 76.29 57.38 0.00 0.00
0.00 0.00 0.00 0.00 0.00 900.00 900.00 900.00 900.00 900.00 900.00 900.00 900.00 900.00 600.00 400.00 200.00 100.00 0.00 0.00
vˆ
0.000 0.000 0.000 0.200 0.424 0.716 0.925 1.098 1.251 1.390 1.517 1.635 1.746 1.898 2.268 2.500 0.494 0.335 0.000 0.000
Estimation model Xt vˆb t
0.00 0.00 0.00 0.00 396.01 713.90 779.55 810.54 829.20 841.83 851.08 859.94 859.44 828.46 737.16 518.70 246.67 103.16 20.05 3.56
AN US
Dt
M
t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
b t
0.000 0.000 0.000 0.000 0.265 0.514 0.840 1.096 1.315 1.511 1.688 1.848 2.001 2.211 2.392 2.366 0.920 0.425 0.000 0.000
CR IP T
OR_ITER
ED
Table 2: Optimal X t of model OR_ITER, the resulting Rt and the simulated X t and vˆtb These results demonstrate that the solution diverges from the optimum of OR_MON and the
PT
optimum is not a fixed point.
Table 1 also demonstrates that the dual prices ut* and the estimated backward lead times vˆtb over
CE
time are not proportional even in the congested interval (periods 4 to 13). This is to be expected since the first derivative of the clearing function does not satisfy the condition given in Equation (32).
AC
In Example 2 we demonstrate the effect of high utilization. The parameter k of the clearing function (33) is set to 30 which makes the clearing function rather flat (close to the capacity) for moderate to high WIP levels. Demand exceeds capacity over several consecutive periods which forces the model to operate in this high-utilization regime. All other data remain unchanged. The left-hand side of Table 3 shows the results of OR_MON. The right-had side presents the relevant results of model OR_ITER using ut* and vˆtb resulting from OR_MON. The objective function is the Lagrangian function (25) discarding the additive constant
uC t
t
and substituting (34) for the WIP holding cost
term as above.
27
ACCEPTED MANUSCRIPT
Rt
Xt
vˆtb
ut*
0 0.00 0 0.00 0 0.00 0 10.06 0 1041.93 1000 964.26 1000 956.90 1000 954.26 1000 952.96 1000 952.21 1000 951.73 1000 951.40 1000 951.17 900 753.12 600 111.43 400 370.39 200 186.18 100 95.53 0 0.00 0 0.00
0.00 0.00 0.00 9.75 831.58 865.93 881.26 890.43 896.70 901.33 904.93 907.84 910.24 900.00 600.00 400.00 200.00 100.00 0.00 0.00
0.000 0.000 0.000 0.031 0.202 0.320 0.402 0.470 0.530 0.583 0.633 0.679 0.722 0.717 0.462 0.059 0.043 0.037 0.000 0.000
0.000 0.000 0.000 0.516 1.516 2.516 3.516 4.516 5.516 6.516 7.516 8.516 9.516 6.200 0.616 0.547 0.525 1.559 0.000 0.000
Dt
0.00 0.00 0.00 950.00 950.00 950.00 950.00 950.00 950.00 950.00 950.00 400.00 0.00 900.00 600.00 400.00 300.00 0.00 0.00 0.00
AN US
t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
OR_ITER with Lagrangian objective function Capacity Xt Rt dual price 0.00 0.00 0.00 985.53 1132.05 1063.58 1028.61 1015.24 1006.68 1001.32 632.41 134.59 646.57 596.22 388.93 285.70 282.58 0.00 0.00 0.00
0.000 0.000 0.000 0.320 0.270 0.198 0.148 0.110 0.078 0.050 0.024 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
CR IP T
OR_MON
M
Table 3: Optimal solution for high utilization As expected the result for OR_ITER is remarkably different from the optimum. The dual prices of
ED
the capacity constraints (5) are close to zero even in highly utilized periods, X13 is zero at reduced cost of 0.0222 which indicates that the difference to alternative solutions is very small. Further numerical studies reveal that if demand oscillates between two values that are below the
PT
capacity, both ut* and vˆtb exhibit an alternating pattern that seems to depend on the parameter setting
CE
of the clearing function. Proportionality between ut* and vˆtb cannot be expected for these irregular demand patterns.
Conclusions and implications for future research
AC
6
This paper analyzed the theoretical basis of a subset of production planning algorithms that iterate between the planning and the scheduling level, namely iterative LP – simulation algorithms for order release planning. The formal analysis is limited to algorithms that iterate solely on the lead times (not on the capacities) and seeks to identify the cause of the convergence problems that are observed frequently. From this analysis we can draw the following conclusions: (1) Adjusting the planned lead times for order release planning by iterating between an order release and a lead time estimation
28
ACCEPTED MANUSCRIPT
model is a dual (price) coordination mechanism although no prices are defined explicitly. Except for special cases this procedure will not lead to the optimal solution since the optimal solution is not a fixed point. (2) There seems to be no straightforward way to overcome this problem within this iterative framework: Calculating the optimal solution to the monolithic model OR_MON by an iterative price coordination mechanism requires minimization of the Lagrangian function (24) which in turn requires estimating the maximum output in the periods for a given trajectory of system states
CR IP T
prior to the period under consideration, and in our context the output estimation is done by simulation. It is unclear whether a simulation model can be integrated into the Lagrangian function in a way that leads to a tractable optimization model. (3) Including capacity constraints in the release model can partly compensate for these deficiencies. Thus, formulating the release model without capacity constraints cannot be recommended although this has been plausibly argued in the literature.
The question is how, if at all, the iterative order release mechanisms can be improved based on
AN US
these insights. This question is beyond the frontier of today’s research, so we can only present some preliminary insights.
An alternative to improving the price (dual) coordination mechanism is to replace/complement it by primal coordination, that is, to limit/assign the amount of resources available. The insight from Corollary 1 that, given realistic values for WIP holding costs and lead times, the lead time updating
M
procedure does not change the planned production quantities strongly supports this way. The iterative order release mechanism in Kim and Kim (2001) that iterates both on the lead times and on the
ED
maximum output seem to exhibit a much better convergence behaviour compared to mechanisms that only iterate on the lead times (see Section 2), and virtually all iterative production planning – scheduling algorithms outside the order release domain iterate on some representation of capacity. It
PT
seems reasonable to assume that iteratively updating the available capacities in the release model acts as a primal coordination mechanism and thus can partly compensate for the weaknesses of the dual
CE
coordination mechanism that is implied by lead time updating. Working out the underlying theory is a topic for future research.
AC
Virtually all iterative order release mechanisms use specific outcomes of the scheduling level (flow times and/or capacity utilization) as feedback information to the release model. This is a remarkable difference from both fix-and-optimize heuristics that use the decision variables directly and simulation-based optimization of the releases where the feedback is derived from the objective function value, such as a gradient in gradient-based algorithms or a fitness value in genetic algorithms. Future research must clarify the optimal design of this feedback mechanism and, since the scheduling level is often represented by a simulation model, the optimal positioning of simulation in hybrid analytical/simulation methods that solve order release and other production planning problems.
29
ACCEPTED MANUSCRIPT
Appendix 1: Proof of Theorem 1 Period-based integer lead times If lead times are integer and defined for each production period as in model OR_ITER, WIP holding costs of h are incurred each time a unit of the product is carried over from one period to the next. Thus the total WIP holding costs are h multiplied by the number of times this happens within the planning
CR IP T
horizon (analogous to the “parts times periods” in lot sizing). Since a WIP inventory at the end of period t implies that Wt product units are carried over from period t to period t+1, the total number of “parts times periods” in the WIP inventory are
W
and the total WIP holding costs are h
t
t
W . t
t
In the optimal solution to model OR_ITER the final WIP WT is zero since there is no incentive to
AN US
release work that is not completed and WIP inventory is penalized in the objective function. Hence all “parts times periods” can be assigned to a certain production period of the respective product units. The number of ”parts times periods” that can be assigned to production in period t is vˆtb X t by the fixed backward lead time assumption, so the number of “parts times periods” is
M
immediately leads to
hWt h vˆtb X t
b t
X t , which
t
(35)
t
ED
t
vˆ
Epoch-based, fractional lead times
Hung and Leachman (1996) assume that the backward lead times are defined at the period boundaries
PT
and can take fractional values at any time. This assumption is different from model OR_ITER. In the following we analyze whether Theorem 1 which states that time-varying lead times have the same
CE
effect as time-dependent production costs extends to the case of fractional, epoch-based lead times assuming uniform output rate within the periods. We denote by Pt the WIP inventory area (parts times
AC
periods) of the work finished in period t. Figure 2 exhibits the cumulative release and output and the WIP inventory areas Pt for assumed values of the backward lead times. For exposition we assume W0 = 0 which implies v0b 0 .
30
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 2: Relationship between cumulative releases, cumulative output, Pt and vˆtb
M
for five periods with vˆtb = 0.75, 1.5, 2, 1.5, 1 for t=1 to 5.
ED
The WIP inventory areas Pt are trapezoids and can be calculated as:
Pt
1 b vˆt 1 vˆtb X t 2
for t=1,.., T
(36)
CE
PT
For all t 1,..., T 1 each vˆtb is multiplied by X t 2 and X t 1 2 . Therefore: T
T 1
t 1
t 1
Pt vˆtb
X t X t 1 X vˆTb T 2 2
(37)
Equation (37) means that for all periods except period T the WIP holding costs per unit of work
AC
( h vˆtb ) are charged to the average production quantity of the periods t and t+1. Hence the insight that the backward lead times act as prices for production is still valid for epoch-based, fractional lead times, but the assignment of the prices to production periods is slightly different. A positive initial WIP W0 0 would require a modification of the calculation of Pt for the first period since for this period the Pt area is no longer a trapezoid (see Figure 2). This does not provide any additional relevant insight here.
31
ACCEPTED MANUSCRIPT
Appendix 2: Multiple optima in model OR_MON with relaxed and dualized clearing function We denote *t and X t* the load t Wt 1 Rt and the output, respectively, in period t at the optimum of OR_MON. The multiple optima maintain the total output
X t
* t
and the values of *t ,
Wt X t *t ,
CR IP T
that is, (38)
which makes the clearing function term in the Lagrangian function (24) an additive constant for given
ut . Substituting Wt *t X t into (24) and eliminating the clearing function term leads to the objective function
h *t h X t l I t ut X t Min ! t
t
t
AN US
t
(39)
The first term is an additive constant for given *t . This also holds for the second term (total output) if production is only shifted between periods. Given these conditions the relevant objective function is just
l I t ut X t Min !
(40)
t
M
t
If, starting from feasible optimal output quantities, production is shifted to the previous period, that is,
ED
X t X t* 1 and X t 1 X t*1 1 for some t, the first term of (40) increases by l and the second term decreases by ut ut 1 . If this change is made in a congested interval the objective function value
PT
remains unchanged since the dual prices satisfy ut* ut*1 l (see Section 4.3). This change is locally feasible since the relevant capacity constraint (the clearing function) is relaxed and the resulting
CE
decrease of Wt 1 due to the unchanged load (Equation (38)) is feasible since the optimal Wt is
AC
positive in a congested interval (see Section 4.3).
32
ACCEPTED MANUSCRIPT
Appendix 3: Proportionality of ut* and t in model OR_MON The shape of the clearing function for which ut* and t are proportional over time within a congested interval can be derived as follows: From Equation (31) and the proportionality condition
t ut* a l t from which t can be obtained as
t a l
we obtain by substituting t from (41) into (31)
f '( t )
a
h . t a
(41)
CR IP T
t
(42)
AN US
The solution of this differential equation leads to the clearing function f ( t ) K h Log t
M
with K the integration constant.
(43)
Acknowledgments
ED
The author wishes to thank Reha Uzsoy for valuable comments on earlier versions of the paper.
Declarations of interest
AC
CE
PT
None
33
ACCEPTED MANUSCRIPT
References Albey, E., & Bilge, Ü. (2011). A hierarchical approach to FMS planning and control with simulationbased capacity anticipation. International Journal of Production Research, 49, 3319-3342. Albey, E., Bilge, U., & Uzsoy, R. (2012). Performance of an iterative production planning approach for production systems with manufacturing flexibility. Pre-Prints of the 17th International Working Seminar on Production Economics 4, Innsbruck, Austria, 1-12.
CR IP T
Asmundsson, J., Rardin, R. L., Turkseven, C. H., & Uzsoy, R. (2009). Production planning with resources subject to congestion. Naval Research Logistics, 56, 142-157.
Bahl, H. C., & Ritzman, L. P. (1984). An integrated model for master scheduling, lot sizing and capacity requirements planning. The Journal of the Operational Research Society, 35, 389-399. Bang, J.-Y., & Kim, Y.-D. (2010). Hierarchical production planning for semiconductor wafer
AN US
fabrication based on linear programming and discrete-event simulation. IEEE Transactions on Automation Science & Engineering, 7, 326-336.
Bertrand, J. W. M., Wortmann, J. C., & Wijngaard, J. (1990). Production Control: A Structural and Design Oriented Approach. Amsterdam: Elsevier.
Bilgen, B., & Çelebi, Y. (2013). Integrated production scheduling and distribution planning in dairy supply chain by hybrid modelling. Annals of Operations Research, 211, 55-82.
M
Boyd, S., & Vandenberghe, L. (2004). Convex Optimization. Cambridge, UK: Cambridge University Press.
ED
Byrne, M. D., & Bakir, M. A. (1999). Production planning using a hybrid simulation-analytical approach. International Journal of Production Economics, 59, 305-311. Byrne, M. D., & Hossain, M. M. (2005). Production planning: an improved hybrid approach.
PT
International Journal of Production Economics, 93-94, 225-229.
69.
CE
Carey, M. (1987). Optimal time-varying flows on congested networks. Operations Research,. 35, 58-
Carey, M. (1992). Nonconvexity of the dynamic traffic assignment problem. Transportation Research
AC
Part B: Methodological, 26, 127-133. Dauzère–Pérès, S., & Lasserre, J. (2002). On the importance of sequencing decisions in production planning and scheduling. International Transactions in Operational Research, 9, 779-793.
de Kok, A. G., & Fransoo, J. C. (2003). Planning supply chain operations: definition and comparison of planning concepts. In A. G. de Kok & S. C. Graves (Eds.), OR Handbook on Supply Chain Management. Amsterdam: Elsevier, 597-675. Dirickx, Y. M. I., & Jennergren, L. P. (1979). Systems Analysis by Multilevel Methods: With Applications to Economics and Management. Chichester et al.: John Wiley & Sons.
34
ACCEPTED MANUSCRIPT
Figueira, G., & Almada-Lobo, B. (2014). Hybrid simulation–optimization methods: A taxonomy and discussion. Simulation Modelling Practice & Theory, 46, 118-134. Graves, S. C. (1982). Using Langrangean techniques to solve hierarchical production planning problems. Management Science, 28, 260-275. Hackman, S. T. (2008). Production Economics. Berlin: Springer. Hung, Y. F., & Hou, M. C. (2001). A Production Planning Approach Based on Iterations of Linear
Industrial Engineers, 18, 55-67.
CR IP T
Programming Optimization and Flow Time Prediction. Journal of the Chinese Institute of
Hung, Y. F., & Leachman, R. C. (1996). A production planning methodology for semiconductor manufacturing based on iterative simulation and linear programming calculations. IEEE Transactions on Semiconductor Manufacturing, 9, 257-269.
Irdem, D. F., Kacar, N. B., & Uzsoy, R. (2008). An experimental study of an iterative simulation-
AN US
optimization algorithm for production planning. 2008 Winter Simulation Conference. IEEE, Austin. 2176-2184.
Irdem, D. F., Kacar, N. B., & Uzsoy, R. (2010). An exploratory analysis of two iterative linear programming - simulation approaches for production planning. IEEE Transactions on Semiconductor Manufacturing, 23, 442-455.
M
Karmarkar, U. S. (1989). Capacity loading and release planning with work-in-progress (WIP) and lead-times. Journal of Manufacturing and Operations Management, 2, 105-123.
ED
Kefeli, A., Uzsoy, R., Fathi, Y., & Kay, M. (2011). Using a mathematical programming model to examine the marginal price of capacitated resources. International Journal of Production Economics, 131, 383-391.
PT
Kim, B., & Kim, S. (2001). Extended model for a hybrid production planning approach. International Journal of Production Economics, 73, 165-173.
CE
Kim, S. H., & Lee, Y. H. (2016). Synchronized production planning and scheduling in semiconductor fabrication. Computers and Industrial Engineering, 96, 72-85.
AC
Kim, S., & Uzsoy, R. (2008). Integrated planning of production and engineering process improvement. IEEE Transactions on Semiconductor Manufacturing, 21, 390-398.
Leachman, R. C., & Carmon, T. F. (1992). On capacity modeling for production planning with alternative machine types. IIE Transactions, 24, 62-72.
Luenberger, D. G. (1989). Linear and nonlinear programming. Massachusetts: Addison-Wesley Reading. Mesarovic, M. D., Macko, D., & Takahara, Y. (1970). Theory of Hierarchical, Multilevel, Systems: New York: Academic Press.
35
ACCEPTED MANUSCRIPT
Milne, R. J., Mahapatra, S., & Wang, C.-T. (2015). Optimizing Planned Lead Times for Enhancing Performance of MRP Systems. International Journal of Production Economics, 167, 220-231. Missbauer, H. (2002). Aggregate order release planning for time-varying demand. International Journal of Production Research, 40, 688-718. Missbauer, H. (2009). Models of the transient behaviour of production units to optimize the aggregate material flow. International Journal of Production Economics, 118, 387-397.
CR IP T
Missbauer, H., & Uzsoy, R. (2011). Optimization models of production planning problems. In K. G. Kempf, P. Keskinocak, & R. Uzsoy (Eds.), Planning Production and Inventories in the Extended Enterprise: A State of the Art Handbook. Berlin: Springer, 437-508.
Modigliani, F., & Hohn, F. E. (1955). Production planning over time and the nature of the expectation and planning horizon. Econometrica, 23, 46-66.
of Technology, Eindhoven, Netherlands.
AN US
Negenman, E. G. (2000). Material Coordination under Capacity Constraints. Eindhoven University
O'Regan, D., Meehan, M., & Agarwal, R. P. (2001). Contractions. In Fixed Point Theory and Applications (pp. 1-11). Cambridge: Cambridge University Press.
Puergstaller, P., & Missbauer H. (2012). Rule-based vs. optimization-based order release in workload
Economics, 140, 670-680.
M
control: A simulation study of a MTO manufacturer. International Journal of Production
Riaño, G. (2003). Transient Behavior of Stochastic Networks: Application to Production Planning
ED
with Load-Dependent Lead Times. Ph.D. dissertation, Georgia Institute of Technology, School of Industrial and Systems Engineering, Atlanta, GA. Schneeweiss, C. (2003). Distributed Decision Making. Berlin: Springer.
PT
Shanthikumar, J. G., & Sargent, R. G. (1983). A unifying view of hybrid simulation/analytic models and modeling. Operations Research, 31, 1030-1052.
CE
Voss, S., & Woodruff, D. L. (2006). Introduction to Computational Optimization Models for Production Planning in a Supply Chain. (2nd ed.). Berlin: Springer.
AC
Wolosewicz, C., Dauzère-Pérès, S., & Aggoune, R. (2015). A Lagrangian heuristic for an integrated lot-sizing and fixed scheduling problem. European Journal of Operational Research, 244, 3-12.
Zäpfel, G. (1984). Systemanalytische Konzeption der Produktionsplanung und –steuerung für Betriebe der Fertigungsindustrie. In C. Zink (Ed.), Sozio-technologische Systemgestaltung als Zukunftsaufgabe. München (Munich): Carl Hanser Verlag (in German). Zhao, S., Wang, L., & Zheng, Y. (2014). Integrating production planning and maintenance: an iterative method. Industrial Management & Data Systems, 114, 162-182.
36