A Model for Demand Planning in Supply Chains with Congestion Effects

: This paper is concerned with demand planning for internal supply chains consisting of workstations, production facilities, warehouses, and transportation links. We address the issue of how to help a supplier ﬁrmly accept orders and subsequently plan to fulﬁll demand. We ﬁrst formulate a linear aggregate planning model for demand management that incorporates elements of order promising, recipe run constraints, and capacity limitations. Using several scenarios, we discuss the use of the model in demand planning and capacity planning to help a supplier ﬁrmly respond to requests for quotations. We extend the model to incorporate congestion effects at assembly and blending nodes using clearing functions; the resulting model is nonlinear. We develop and test two algorithms to solve the nonlinear model: one based on inner approximation and the other on outer approximation.


Introduction and Literature Review
The tremendous growth of electronic commerce has focused the attention of supplier firms on the need for both internal and external supply chain integration. To be competitive, it is now necessary not only to synchronize the activities of one's own internal supply chain, but also to develop and manage a partnership with external business partners and customers. Beamon [1] notes that there has been increasing attention paid to the performance, design, and analysis of the supply chain as a whole. Keskinocak and Tayur [2] discuss the role of quantitative models in the electronic commerce context. Customers have more power and more complex requirements in such a context, which makes flexibility and customer satisfaction the two major concerns in production-distribution system planning. Thus, supply chain firms are moving toward customer-oriented planning and decisionmaking systems (Askin and Goldberg [3]).
Electronic commerce also makes external supply chain integration possible. External supply chain integration includes negotiation-based due dates and price setting. Highspeed internet connections allow instant information exchange between supply chain entities, and consequently reduce transaction costs. A firm can work on the due dates, configuration, and price of products together with its business partners for mutual benefit. An internal integration plan is frequently updated as a result of negotiations, while at the same time providing feedback that impacts the negotiation process itself. Taken to the extreme, the electronic commerce context can allow a firm to establish a seamless interface with its external supply chain. This paper draws on several streams in the literature. The first stream is concerned with multi-plant coordination and material requirements planning (MRP) systems modelled at the aggregate level. Billington et al. [4] develop mathematical formulations for a time is captured by fitting a concave curve to the results of the simulation. Pahl et al. [21] summarize the literature on production planning with load-dependent lead times. In the context of a supply chain and operations planning system, Selcuk et al. [22] show that the shape of the clearing function plays an important role in the completion time of orders. Asmundsson et al. [23] develop a production planning model for a singlestage multi-product system that captures the nonlinear relationship between resource workload and lead times. They use outer linearization to linearize the clearing function and extend it to multi-stage systems. Missbauer and Uzsoy [24] examine the use of clearing functions in multi-stage production planning systems, and discuss how the shape of the clearing function at a particular stage could be dependent on the decisions made at earlier upstream stages. Kacar et al. [25] develop a linearized clearing-function-based production planning model and apply it to a large-scale wafer fabrication facility using two products and hundreds of operations. They use simulation to show that the clearing function model yields higher profit than the tradition LP-based production planning models. Charnsirisaksul et al. [26] look at integrated order selection and scheduling, where the manufacturer has the flexibility to choose lead times and uses that to maximize profit. More recently, Kefeli and Uzsoy [27] identify bottlenecks in production systems using dual prices in the presence of congestion at resource nodes. They compare the dual price information in the congestion-based production planning model in Asmundsson et al. [23] with the fixed lead time model that does not capture queuing behaviour (as used in MRP-based systems). An interesting conclusion in this paper is that improvements at a workcentre with non-zero dual prices, even though it is not the principal bottleneck, can lead to improvements in the overall system.
In bringing these three streams together (MRP/ERP, demand management, and congestion in production planning), this paper makes an important conceptual contribution. For example, the papers on congestion in production systems, such as Asmundsson et al. [23] and Srinivasan et al. [17], do not deal with material availability issues. On the other hand, the literature on MRP/ERP and demand management does not deal with congestion. From a practical standpoint, this is relevant to global network production. From a technical standpoint, this work could be seen as an important step to incorporate load-based lead times in MRP/ERP. The remainder of this paper is organized as follows. We present a linear optimization model for demand planning in Section 2 and illustrate the use of the model in Section 3. In Section 4, we develop an extension to the basic demand planning model that incorporates congestion effects at assembly and blending nodes through the use of recipes. This is the main contribution of the paper, since all previous work on congestion modelling in production planning has been developed for simple products or WIP. The resulting optimization model is nonlinear, and we develop linear programming-based algorithms to solve the model and report on the computational performance. Section 5 concludes the paper by summarizing the work and pointing out future research directions.

Background
The setting considered for our model is as follows. A supplier firm receives requests for quotations (RFQs) from its customers. Some of the quotations are accepted, whereupon they become committed demands, i.e., product and time commitments made by the firm to its customers. The problem of interest is how a supplier firm should respond to RFQs in a capacitated multi-period multi-commodity network, while at the same time planning and allocating its internal flows.
The demand-planning model described in this section involves a multi-commodity flow allocation problem; the model's objective is to minimize costs. When the problem is solved, the supplier firm is able to know what quotations to make in response to the RFQs (product and date). The model also gives the firm an estimate of the price to be quoted. The following representation is used for the supply network of the supplying firm in question:

•
Nodes in the network represent external nodes, such as suppliers or customers, and internal nodes, such as workstations, plants, and warehouses. Production nodes involve product transformation resulting from manufacturing, assembly, or blending based on recipes (see below). • Arcs in the network represent transportation between nodes. In order to simplify the notation in our model, activities and capacity restrictions are imposed on nodes instead of arcs. • A recipe is defined as a process step with input and output products. There is a subtle difference between the recipe concept and the bill-of-material concept; a recipe models inputs and outputs at the process level, while a bill of material is always defined at the product level. For example, a node could represent an assembly step. In this case, a step in an assembly bill of material could be regarded as a recipe with several discrete inputs and one output. However, if the node represents a plant, a recipe looks at inputs and outputs at the plant level. The intermediate steps in the bills of materials are collapsed to reflect that process level. A recipe could also represent blending with by-products and is general enough to model any discrete process.

•
A new demand or RFQ is defined as a demand under negotiation. In this case, lateness is minimized through "artificial" penalties, which are penalties set by user discretion.

•
As negotiation evolves, a new demand turns into a committed demand when a quotation (which is a response to an RFQ) is accepted by a customer.
The following are assumed: • Flows arrive or leave at trans-shipment and demand nodes exactly midway through a time period. They stay for a minimum of one period at a node. • At production nodes, it is assumed that flows arrive at the beginning of a time period; a fraction β of flow may be used for production during the time period (0 ≤ β ≤ 1).

Notation
The notation used is consistent with Venkatadri et al. [13]. The indices used in the formulation are: The directed arcs Γ(t,j,t + l(j,k,m),k) represent allowable routings from the output of resource j to the input of resource k for item m in time period t, where l(j,k,m) is the lead time between nodes j and k for item m. Supplier nodes do not have predecessors and demand nodes do not have successors.
The decision variables are as follows: Logistics 2021, 5, 3 5 of 24 f m t,j,t + l(j,k,m),k flow of item m on arc (t,j,t + l(j,k,m),k) x m t,j initial inventory of item m at node j ∈ (I∪D) in time period t N t , j,r actual number of production runs made using recipe r ∈ R j at node j ∈ IP in time period t upper bound on the quantity of item m to order from node j ∈ S in time period t C t,j available capacity of node j ∈ I in time period t. In the case of trans-shipment nodes, C t,j is in storage units because space availability is the primary concern in trans-shipment resources.
β a value between 0 and 1 indicating the percentage of flow coming into an IP node during a time period that can be used for production Parameters relating to assembly or blending are as follows: OPT j,r output item set at node j ∈ IP, using recipe r ∈ R j IPT j,r input item set at node j ∈ IP, using recipe r ∈ R j P j recipe set at node j ∈ IP RO m j,r number of units of output item m∈ OPT j,r produced when recipe r ∈ R j is run at node j RI m j,r number of units of input item m∈ IPT j,r consumed when recipe r ∈ R j is run at node j V j,r capacity utilized by one production run of recipe r ∈ R j at node j ∈ IP. V j,r and C t,j in IP nodes are usually in available machine time. IPT j,r and OPT j,r come from the bill of materials and describe recipe r.
The parameters relating to demand planning are as follows: The ordering cost s m t,j"k,t+l(j,k,m) changes depending on the supplier selected and may change for the same supplier at different time periods. For a node k ∈ (I∪D), the lead time may be different depending on the connected supplier. Parameter c j is only dependent on resource j, which can be interpreted as labour cost or facility maintenance cost. Parameter c m j is dependent on both item m and resource j, which can be interpreted as production or handling cost. All other unit node costs can also be incorporated into either c j or c m j .

The Demand Planning Model
Without loss of generality, it is assumed that lateness to committed demand is allowed with a penalty.
Our objective (1) is to minimize overall costs: The first term in (1) is the flow-related total cost FC(f t,j,k,m ), defined as: According to the assumptions, flows stay for at least one time period when going through a node. Thus, the ordering, node-independent, and node-dependent costs, as seen in the first three terms of FC(f t,j,k,m ), are flow dependent. The fourth term in the flow cost function FC(f t,j,k,m ) represents the arc costs, which usually include transportation and holding costs.
The second term in (1) is the inventory cost for all nodes in {I} and {D}. If the inventory at {D} is not the firm's responsibility, it can be omitted. The third and fourth terms in the objective function are the lateness costs for committed demand and new demand, respectively. The lateness penalty for committed demand should reflect the value of contractual penalty and loss of goodwill. The lateness penalty for new demand can be assigned a very small number to give incentive for the model to satisfy new demand when applicable. It is also a goodwill cost, but in the sense that it should reflect the goodwill incurred by trying to satisfy an RFQ. The last two terms in the objective function are revenues for satisfying committed and new demand, respectively. They provide further incentives for order completion.
Constraint (2) is the inventory conservation equation for each item at trans-shipment nodes. Inventory is dependent on the inflow and outflow at the node.
Constraint (3) is the inventory conservation equation for each input item at production nodes. Inventory is dependent on the inflow and quantity consumed in assembly at the node. Constraint (4) is inventory conservation equation for each output item at production nodes. Inventory is dependent on the quantity transferred in assembly and outflow at the node. ∑ r∈R j N t,j,r RI m j,r and ∑ r∈R j N t,j,r RO m j,r represent the number of input components consumed for assembly over all recipes used and the number of output components produced, respectively.
x m t+1,j = x m t,j + ∑ k∈B m j f m t−l(k,j,m),k,t,j − ∑ r∈R j N t,j,r RI m j,r ∀j ∈ IP, ∀t, ∀m ∈ IPT j,r x m t+1,j = x m t,j − ∑ k∈A m j f m t,j,t+l(j,k,m),k + ∑ r∈R j N t,j,r RO m j,r ∀j ∈ IP, ∀t, ∀m ∈ OPT j,r Constraint (5) is the inventory conservation equation for each input item at demand nodes. Inventory is dependent on the inflow and satisfied demand.
Constraint (6) limits the outflow at each trans-shipment node.
Constraint (7) imposes capacity limitations on production nodes.
Constraints (8) and (9) are imposed to ensure that the outflow does not exceed the initial inventory.
∑ k∈A m j f m t,j,t+l(j,k,m),k ≤ x m t,j ∀j ∈ IT, ∀t, ∀m (8) ∑ k∈A m j f m t,j,t+l(j,k,m),k ≤ x m t,j ∀j ∈ IP, ∀t, ∀m ∈ OPT j,r Similarly, Constraint (10) ensures that the total input item m consumed in time period t is less than or equal to the sum of its initial inventory and a proportion of the inflow in time period t decided by the parameter β. ∑ ∀r∈R j N t,j,r RI m j,r ≤ x m t,j + β ∑ k∈B m j f m t−l(k,j,m),k,t,j ∀j ∈ IP, ∀t, ∀m ∈ IPT j,r Equalities (11) and (12) specify upper and lower bounds on the quantity of raw materials that can be bought from all suppliers.
∑ k∈A m j f m t,j,t+l(j,k,m),k ≤ UQ m t,j ∀j ∈ S, ∀t, ∀m (11) LQ m t,j ≤ ∑ k∈A m j f m t,j,t+l(j,k,m),k ∀j ∈ S, ∀t, ∀m (12) Equalities (13) and (14) define the decision variables sd m t,j (satisfied committed demand) and sn m t,j (satisfied new demand). In addition, sd m t,j and sn m t,j take into account delayed demand accumulated from previous time periods.
Constraint sets (15) and (16) are bounds on inventory at nodes and non-negativity requirements, respectively. If lateness is not desired or allowed, the model can be easily simplified.
LX m j ≤ x m t,j ≤ UX m j ∀j ∈ I, ∀t, ∀m All The optimization model then is:

Using the Model
DPP is a versatile model that can be used for order promising, internal flow planning, capacity planning, and bottleneck elimination. The dual price on the capacity bundle Constraints (6) and (7) can be interpreted as the improvement in the objective function if an additional resource j unit in time period t is made available. The additional resource unit means extra storage units or production time in trans-shipment nodes or production nodes, respectively. Bottlenecks can be identified by looking at these dual prices.

Bottleneck Elimination for Capacity Planning
The direct way to eliminate an identified bottleneck is to increase capacity. An updated model with increased capacity could be solved to find and alleviate the next bottleneck resource. Another option is to modify the capacity constraint for capacity planning. Instead of treating capacity as just one parameter, it may be broken up into several levels. For example, there may be two capacity categories at a node j, regular time capacity and overtime capacity, with overtime capacity having a higher cost. In general, there could be n capacity categories. The decision variable C n,t,j represents the capacity in category n at node j in time period t. The capacity parameter UC n,t,j is the upper bound of C n,t,j . The unit capacity cost in category n at node j in time period t is CC n,t,j . By replacing Constraints (6) and (7) by (17) and (18), respectively, and adding Constraint (19) and Bound (20), we get the new model DPP .
C n,t,j ≤ UC n,t,j ∀j ∈ IP, ∀n, ∀t (19) C n,t,j ≥ 0 ∀j ∈ I, ∀n, ∀t (20) The DPP objective function (1) is modified as follows by adding the term C n,t,j : Model DPP is now solved as: Minimize (21) Subject to: (2) to (5), (8) to (20) Assuming that CC m,t,j (for example, regular time cost) < CC m+1,t,j (for example, overtime cost) DPP will always choose less expensive capacity modes before choosing more expensive capacity. Thus, explicit information about the bottleneck could be obtained, and the bottleneck could then be removed by adding new capacity tiers.
Several other extensions to DPP are possible; however, they are not discussed due to lack of space. One such extension is explicitly modelling customer orders and the lateness of those orders. Details can be found in Wang [28]. The DPP model is quite close to the models in Chen et al. [29] and Chen et al. [30], who pioneered the research literature in order promising. Ball et al. [31] present a generic framework for order promising and discuss applications at Dell and Toshiba. There are other models in the literature of this nature. As already mentioned, Arntzen et al. [8] formulated a large-scale model for global production and distribution (GSCM) that incorporates multi-product bill of materials for a large electronic supply chain. Shapiro et al. [32] developed a strategic production and distribution planning application for a large consumer products company. Degbotse et al. [33] presented Logistics 2021, 5, 3 9 of 24 a semiconductor supply chain planning strategy implemented at IBM. They decomposed the problem by dividing the bills of materials product structure horizontally and vertically into complex and simple portions for the stages of semiconductor manufacturing; the complex portion was solved with an MIP and the simple portion with heuristics containing embedded LPs. Fordyce et al. [34] discussed planning and scheduling in semiconductorbased packaged goods companies.

Illustrative Example
The business case for such systems is seen in the electronics and aviation industries, where queuing effects are experienced due to variations in processing times. For example, in the electronic industry, semiconductor wafers undergo reentrant process steps through several workstations, effectively creating a queuing network. In the aviation industry, machining can be quite intricate and varies every time a component is made, making the processing time non-deterministic. In addition, in the defense aviation industry, assembly time is non-deterministic because it is not only complex, but also involves extensive testing.
Consider a simple internal supply chain network with three assembly plants, two warehouses, one upstream supplier, and three downstream customers ( Figure 1). Material flow starts from suppliers, undergoes transformation and finally reaches customers. NC, TC, and LT are abbreviations for the node cost, transportation cost, and lead time, respectively. Note that direct shipping is allowed between assembly node 4 in stage 3 and customer node 7 in stage 5, with a higher transportation cost. Node 7 is supplied by nodes 4, 5, and 6, and node 8 is supplied by nodes 5 and 6. Node 9 is supplied only by node 6. We assume that the supplier's supplier charges a constant price for raw materials. The business case for such systems is seen in the electronics and aviation industries, where queuing effects are experienced due to variations in processing times. For example, in the electronic industry, semiconductor wafers undergo reentrant process steps through several workstations, effectively creating a queuing network. In the aviation industry, machining can be quite intricate and varies every time a component is made, making the processing time non-deterministic. In addition, in the defense aviation industry, assembly time is non-deterministic because it is not only complex, but also involves extensive testing.
Consider a simple internal supply chain network with three assembly plants, two warehouses, one upstream supplier, and three downstream customers ( Figure 1). Material flow starts from suppliers, undergoes transformation and finally reaches customers. NC, TC, and LT are abbreviations for the node cost, transportation cost, and lead time, respectively. Note that direct shipping is allowed between assembly node 4 in stage 3 and customer node 7 in stage 5, with a higher transportation cost. Node 7 is supplied by nodes 4, 5, and 6, and node 8 is supplied by nodes 5 and 6. Node 9 is supplied only by node 6. We assume that the supplier's supplier charges a constant price for raw materials. The bill of materials is shown in Figure 2. Tables 1-4 show the parameter settings. Raw material (items 1 to 6) are purchased by the firm from a supplier in stage 1. They then enter assembly nodes 2 or 3, where two recipes are used. Sub-final items 3 and 7 leave stage 2 and enter assembly node 4 for final production. End product 8 goes through either warehouse 5 or 6 in stage 4 to reach the customers in the final stage. We also set an upper bound on space requirements at assembly nodes of 600 and 300 units for input items and for output items, respectively. The length of each time period is one week. The bill of materials is shown in Figure 2. Tables 1-4 show the parameter settings. Raw material (items 1 to 6) are purchased by the firm from a supplier in stage 1. They then enter assembly nodes 2 or 3, where two recipes are used. Sub-final items 3 and 7 leave stage 2 and enter assembly node 4 for final production. End product 8 goes through either warehouse 5 or 6 in stage 4 to reach the customers in the final stage. We also set an upper bound on space requirements at assembly nodes of 600 and 300 units for input items and for output items, respectively. The length of each time period is one week.   Table 3. Capacity consumed per time period by a production run in production nodes. 4  1  1  4  2  2  5  1  1  5  2  2  6 3 3    Table 3. Capacity consumed per time period by a production run in production nodes. 4  1  1  4  2  2  5  1  1  5  2  2  6 3 3 The scenarios below show how the model can be used for demand planning. In all the scenarios below, the value of β in DPP was assumed to be 0. All problems below were run on a Pentium III Personal Computer with a Celeron CPU running at 850 MHZ using MINOS 5.5. The problems ran in less than a few seconds. The objective function values below are in the same units (dollars, Euros, yen, etc.) as the unit, unit holding, and transportation costs.

Base-Case Scenario: Quoting Due Dates to Customers
It is assumed that initial inventories are 200 units in IT nodes, and 100 units and 50 units in IP nodes for input and output products, respectively. Table 5 shows the RFQ data at customer nodes 7, 8, and 9 for time periods 6 to 12. Let us assume that none of this demand is committed, that is, all demand values are treated as new demand or RFQ with assumed lateness cost on new demand. Solving DPP, the optimal objective function value obtained is 227,204.74, showing that all deliveries can be made within 12 time periods. However, this does not mean every delivery can be made in time. Scenario 1 (base case) in Table 6 shows the total amount of lateness to customers 7, 8, and 9. The lateness values in the table should be used by the firm to promise orders. Assume that customer 8 increases demand from time period 10 from 55 units to 90 units even before the firm has had a chance to respond to the RFQ in Table 6. By the running the model again, it is seen that this change causes increase in lateness for customers 7 and 9, and the optimal objective function value is 278,427.98 (see Scenario 2 in Table 6). The increase in objective function value includes two types of cost: the cost to push more product through the system and the cost of late deliveries at customers 7 and 9, beyond the lateness they agreed to. Clearly, the firm should charge customer 8 at least the difference in cost (278,427.98 − 227,204.74 = 51,223.24) in order to accept the increase.

New Recipe
Assume that the customer at node 8 comes in with an RFQ for a second item that has a new bill of material ( Figure 3) in addition to the older RFQ for item 8, all before the firm has responded to any RFQ. The new final item 12 is assumed to have higher RFQ lateness cost than item 8. The customer asks for an additional 25 and 15 units of item 12 at time periods 11 and 12, respectively.
Objective Function Value = 159,620.03

Order Change: Increased Demand for Customer 8
Assume that customer 8 increases demand from time period 10 from 55 units to 90 units even before the firm has had a chance to respond to the RFQ in Table 6. By the running the model again, it is seen that this change causes increase in lateness for customers 7 and 9, and the optimal objective function value is 278,427.98 (see Scenario 2 in Table 6). The increase in objective function value includes two types of cost: the cost to push more product through the system and the cost of late deliveries at customers 7 and 9, beyond the lateness they agreed to. Clearly, the firm should charge customer 8 at least the difference in cost (278,427.98 − 227,204.74 = 51,223.24) in order to accept the increase.

New Recipe
Assume that the customer at node 8 comes in with an RFQ for a second item that has a new bill of material (Figure 3) in addition to the older RFQ for item 8, all before the firm has responded to any RFQ. The new final item 12 is assumed to have higher RFQ lateness cost than item 8. The customer asks for an additional 25 and 15 units of item 12 at time periods 11 and 12, respectively.
After solving DPP, the result shows that this scenario causes even more lateness for the other two customers. The optimal objective function value is 257674.85 (Scenario 3 in Table 6). Thus, the firm should charge the customer at least the difference in cost (257,674.85 − 227,204.74 = 30,470.11) to accept this change to the RFQ.

Order Cancellation
Going back to the basic scenario, it is assumed that the RFQs are accepted. Subsequently, customer at node 8 cancels his or her order in time period 10. DPP is solved again to evaluate the consequence of this cancellation. The optimal objective function value drops to 159,620.03 (Scenario 4 in Table 6), and lateness for the other two customers is greatly improved. Only customers 7 and 9 will have late deliveries; in each case, 6.67 units of demand in time period 11 will be delayed by 1 period.

Demand Planning and Capacity Bottleneck Alleviation
Another feature of this model is its ability to identify the resource bottlenecks of each time period dynamically. The goal is minimizing total cost, instead of lateness. Four decision factors determine how to alleviate a bottleneck: dual price, sensitivity, cost of adding one unit of capacity, and the ceiling on how much capacity can be added. At any given point in time, the resource with the highest positive value after subtracting the corresponding cost of adding a unit of capacity from its dual price is regarded as the bottleneck at that time period. The amount of capacity to be added is decided by the sensitivity range After solving DPP, the result shows that this scenario causes even more lateness for the other two customers. The optimal objective function value is 257674.85 (Scenario 3 in Table 6). Thus, the firm should charge the customer at least the difference in cost (257,674.85 − 227,204.74 = 30,470.11) to accept this change to the RFQ.

Order Cancellation
Going back to the basic scenario, it is assumed that the RFQs are accepted. Subsequently, customer at node 8 cancels his or her order in time period 10. DPP is solved again to evaluate the consequence of this cancellation. The optimal objective function value drops to 159,620.03 (Scenario 4 in Table 6), and lateness for the other two customers is greatly improved. Only customers 7 and 9 will have late deliveries; in each case, 6.67 units of demand in time period 11 will be delayed by 1 period.

Demand Planning and Capacity Bottleneck Alleviation
Another feature of this model is its ability to identify the resource bottlenecks of each time period dynamically. The goal is minimizing total cost, instead of lateness. Four decision factors determine how to alleviate a bottleneck: dual price, sensitivity, cost of adding one unit of capacity, and the ceiling on how much capacity can be added. At any given point in time, the resource with the highest positive value after subtracting the corresponding cost of adding a unit of capacity from its dual price is regarded as the bottleneck at that time period. The amount of capacity to be added is decided by the sensitivity range and the ceiling. Demand managers can find and remove core bottlenecks using the algorithm described below:

Bottleneck Alleviation Algorithm
Step 1: Solve the model DPP and obtain resource dual prices and sensitivity ranges.
Step 2: Find the resource with the highest value after subtracting the corresponding cost of adding capacity from its dual price. If the value is positive, go to step 3; otherwise, stop.
Step 3: Add capacity up to the ceiling or upper bound determined by the sensitivity range of the bottleneck.
Step 4: Go to step 1. The procedure thus alleviates bottlenecks by solving the model repeatedly after capacity modifications. Consider the base scenario and the effects of adding overtime capacity, which cannot exceed 15% of the regular capacity for all production. Assume that the cost of adding one unit of capacity is 100 for node 4 and 130 for node 5. After DPP is solved, the dual prices on the capacity constraint are as specified in Run 1 of Table 7. It is clear that node 5 at time period 2 is the bottleneck. Twenty-three extra capacity units are added to the bottleneck (this is less than 15% of the basic capacity), based on the sensitivity range (823.33-800), and the model is run again. The optimal objective function value is 225,907.56 (see Run 2 of Table 7). The total lateness for all nodes drops to 178.52, and node 4 at time period 2 becomes the new bottleneck. In run 3, the problem is solved again with 615 units of capacity in node 4, time period 2, because the sensitivity range indicates that the capacity can be increased up to 615.56 units. After this run, the optimal objective function value is 221,929.13, and the total lateness across all nodes is 171. Node 4 at time period 2 is still the bottleneck (run 3 of Table 7). Although increasing capacity up to 708.67 is the best solution, only 75 more capacity units can be added to the bottleneck because of the 15% ceiling rule. In the last run shown in Table 7, the optimal objective function value decreases to 212,500.85 (decreased by 4.24% from the third run), and the total lateness for all nodes is 133.51 (decreased by 21.92% from the third run). The bottleneck moves back to node 5 at time period 2. The process of bottleneck elimination can be continued in this fashion until it is not economical to alleviate the bottleneck.

Modelling Congestion Effects
In this section, we extend the model developed in Section 2 to incorporate congestion effects.

Clearing Functions
The throughput of a network node (output in a time period) depends on both the nominal capacity of the node and the WIP of input product available, the former being the upper bound on throughput. Without any congestion in the system, the clearing function is linear (function 1 in Figure 4). Note that function 1 is implicit in the DPP model described in Section 2 and is defined by Constraints (7) and (10). However, when there is congestion in the system, the WIP versus throughput relationship is more like function 2 shown in Figure 4. Karmarkar [18] and Srinivasan et al. [17] introduced the concept of the concave nonlinear capacity function to link WIP and lead time in capacity planning for discrete time period models by developing clearing functions based on queuing models. The concavity assumption is appropriate for most production facilities, where production rate increases asymptotically to a limit as the WIP level increases. Karmarkar [18] developed clearing function formulations for a single commodity model in an M/M/1 queuing system, and a multi-commodity model in an M/G/1 queuing system, using continuous time periods. An analogous clearing function for a single-product discrete time period model based on input/output control was also discussed. Conway et al. [35] found that the shape of the capacity versus WIP curve was sharply concave for serial production lines. Bhatnagar et al. [36] also observed the concavity of the clearing function in assembly systems using a simulation approach.
indicates that the capacity can be increased up to 615.56 units. After this run, the optimal objective function value is 221,929.13, and the total lateness across all nodes is 171. Node 4 at time period 2 is still the bottleneck (run 3 of Table 7). Although increasing capacity up to 708.67 is the best solution, only 75 more capacity units can be added to the bottleneck because of the 15% ceiling rule. In the last run shown in Table 7, the optimal objective function value decreases to 212,500.85 (decreased by 4.24% from the third run), and the total lateness for all nodes is 133.51 (decreased by 21.92% from the third run). The bottleneck moves back to node 5 at time period 2. The process of bottleneck elimination can be continued in this fashion until it is not economical to alleviate the bottleneck.

Modelling Congestion Effects
In this section, we extend the model developed in Section 2 to incorporate congestion effects.

Clearing Functions
The throughput of a network node (output in a time period) depends on both the nominal capacity of the node and the WIP of input product available, the former being the upper bound on throughput. Without any congestion in the system, the clearing function is linear (function 1 in Figure 4). Note that function 1 is implicit in the DPP model described in Section 2 and is defined by Constraints (7) and (10). However, when there is congestion in the system, the WIP versus throughput relationship is more like function 2 shown in Figure 4. Karmarkar [18] and Srinivasan et al. [17] introduced the concept of the concave nonlinear capacity function to link WIP and lead time in capacity planning for discrete time period models by developing clearing functions based on queuing models. The concavity assumption is appropriate for most production facilities, where production rate increases asymptotically to a limit as the WIP level increases. Karmarkar [18] developed clearing function formulations for a single commodity model in an M/M/1 queuing system, and a multi-commodity model in an M/G/1 queuing system, using continuous time periods. An analogous clearing function for a single-product discrete time period model based on input/output control was also discussed. Conway et al. [35] found that the shape of the capacity versus WIP curve was sharply concave for serial production lines. Bhatnagar et al. [36] also observed the concavity of the clearing function in assembly systems using a simulation approach.

The Nonlinear Clearing Function Model
A new continuous variable, Zt,j,r ( 0 ≥ ), is introduced to model capacity congestion: Zt,j,r maximum number of production runs made using recipe r ∈ Rj at node j ∈ IP in time period t.

The Nonlinear Clearing Function Model
A new continuous variable, Z t,j,r (≥ 0), is introduced to model capacity congestion: Z t,j,r maximum number of production runs made using recipe r ∈ R j at node j ∈ IP in time period t.
The clearing function illustrations in Figure 4 involve only one product. When several products are involved, the input item that constrains the maximum number of production runs in time period t for a recipe represents the critical WIP for that time period. Therefore, the following constraint set is added to the formulation: RI m j,r ∀j ∈ IP, ∀t, ∀r ∈ R j , ∀m ∈ IPT j,r Z t,j,r ≥ 0 ∀j ∈ IP, ∀t, ∀r ∈ R j The ratio on the right-hand side of (22) represents the number of runs of recipe r at node (t,j) that can be made based on initial inventory of input product m. Since a recipe has many input products, the available (maximum) number of production runs is constrained by the product for which this ratio is the smallest.
The clearing functions described in this paper so far are based on the single-product queue. In order to operationalize these clearing functions, we have to ensure that only one product is produced at a resource at a given time. We first examine the case when only one recipe can be run at a node in a given time period. In this case, the actual number of production runs that can be made at a production node depends on the nominal capacity of the node C t,j , the capacity usage V j,r , and the clearing function, defined as g j,r (Z t,j,r ), and is modelled by Constraint (24). Constraint (7) of the formulation is now redundant (since we will show in the next section the value of the clearing function g j,r (Z t,j,r ) approaches 1 as Z t,j,r approaches infinity), ensures that the sum of the production usage of all recipes at a node does not exceed its nominal capacity. In order to ensure that only one recipe can be run at a time during a given time period, we introduce a new zero/one variable, Y t,j,r , which takes a value of 1 if recipe r is run at node (t,j). In Constraint (25), M is a large number, and therefore, the binary variable Y t,j,r is set to 1 when N t,j,r is non-zero. Constraint (26) enforces the one-recipe at a node-time rule. Constraint (27) simply defines the Y t,j,r 's as binary variables.
The non-linear extension of DPP may now be written as: DPPNL: Minimize (1) Subject to: (2) to (6), (8) to (16) and (22) to (27) It must be noted that both g j,r (Z t,j,r ) and N t , j,r are defined as continuous variables, while in reality, they should both be integers. This may not be an unreasonable assumption in the rolling-horizon planning context, where only the decisions in period one are implemented and all the decisions in subsequent periods are re-planned. Therefore, a simple rounding down of these variable values, while not optimal for discrete parts assembly processes, would be feasible and practical.
As mentioned before, a big difference between the approach in the literature and this work is that instead of working directly with products, the recipes are the virtual entities queueing at the nodes. This is an important contribution of this paper, because we model congestion in assembly or blending through recipes.
In order to allow several recipes to be run at a node in the same time period, we can modify Constraint (24) in a manner analogous to the capacity constraints in Asmundsson et al. [23]. The capacity for each recipe r is partitioned through Constraint (28), separable by recipes as follows: This method uses a partitioning variable, α tjr , for each recipe r as a multiplier in the right-hand side of Constraint (28), with the additional condition that the sum of the α tjr 's is 1, as shown in Constraint (29). Constraints (25) to (27) may be dropped for this case. Also, the variable V j is used instead of variable V j,r , because the average capacity used by one production run of recipe r is now defined in terms of the node and is an average across all recipes at the node.
In what follows, three different clearing functions are adapted from the literature to our model. We will show that they are concave and continuously differentiable, making them computationally tractable (Asmundsson et al. [23] assume that the clearing function for the partitioned capacity case in Constraint (28) is concave). The assumption behind the clearing functions above, which are all separable by recipe r, is that during a period t, one recipe is run repeatedly for steady-state queues to form. In other words, it is assumed that different recipe runs are not mixed together because the queuing behaviour of a multi-class node is different from the queuing behaviour of a single-class node.

Input/Output Clearing Function
Karmarkar [18] employs the following function based on the M/M/1 queue: where W t-1 is the initial WIP in time period t, R t the constant release rate in time period t, X t the actual production rate in time period t, and P the maximum production possible in time period t. The constant k (>0) is used to generate a family of clearing functions. The second argument of the above function is that production cannot also exceed the total of WIP in the previous period and the current release. In our model, Constraint (22) ensures this. Also, is the maximum production rate per time period at node (j,t), using recipe r, with Z t,j,r representing the available (maximum) number of runs, and N t,j,r being the actual number of runs. Hence, the analogous form of the constraint corresponding to the input/output clearing function of Karmarkar [18] is as follows: N t,j,r ≤ C t,j V j,r Z t,j,r Z t,j,r + k ∀j ∈ IP, ∀t, ∀r ∈ R j Clearly, lim Z t,j,r →∞ Z t,j,r Z t,j,r + k = 1 Strictly speaking, this constraint models the clearing function for a single server station. However, since DPPNL is a supply chain network model, constraint may be appropriate where a bottleneck M/M/1 workstation is identified. Alternatively, individual M/M/1 workstations could be modelled in the supply chain network (the model size will then be large).

M/D/1 Clearing Function
This function is adapted from Karmarkar et al. [15], who use an M/G/1 model to study lead time and WIP as a function of batching policy for multiple products. Consider the set IP as an open queuing network where each node is a workstation. Also, assume that there is only one recipe per production node. Since each run represents the processing of a group of components (from the input components set), and assuming that the assembly will not proceed until there is at least one group of components available, ρ j,r = I sever utilization per time period at node j using recipe r. T q = average time spent in queue. Note that batch sizes are considered to be one, and variations in processing time are ignored (since the M/D/1 assumes constant service time).
Using the Pollaczek-Khintchine formula for the M/D/1 queue: Applying Little's law (Hopp and Spearman [37]): Solving for N t,j,r gives: To use the same functional form as in Constraint (24), the above equation is rewritten as: Again, This clearing function representation is applicable for an M/D/1 bottleneck station at an aggregate node in the supply chain network. Alternatively, the workstations in the queuing network can be represented as individual nodes in the network.

General Clearing Function
This function is adapted from Srinivasan et al. [17], and is modelled as follows: where µ can be assigned any value to generate a family of clearing functions. Setting µ= V j,r C t,j results in the value of the function being always below the 45-degree line. To use the same functional form as in Constraint (24), the congestion constraint is written as: As before, lim This clearing function is generic and useful when the queuing behaviour at a node is complex and the underlying process is not understood. However, there may empirical data to fit the clearing function using regression to choose the value of the parameter µ.
The function g j,r (Z t,j,r ) is concave in all three cases. A question that arises is which of these functions should be used. The assumptions behind the first two clearing functions (Sections 4.2.1 and 4.2.2) should be looked at carefully before applying either. In the case of the general clearing function (Section 4.2.3), it can be used almost universally. In any case, the parameters chosen should reflect the throughput characteristics of the production system. DPPNL has a linear objective function with linear constraints and a non-linear constraint, either Constraint (30), Constraint (31) or Constraint (32). Therefore, DPPNL can be solved numerically using a nonlinear programming package such as MINOS, which uses the projected augmented Lagrangian method of Robinson [38]. The next section shows how two linear programming-based algorithms can be used to solve the model.

Algorithms for DPPNL
Two different algorithms are proposed for DPPNL: inner approximation using piecewise linear programming and outer approximation using the Kelley's cutting plane method.

Inner Approximation
The clearing function g j,r (Z t,j,r ) can be represented as a summation of continuous single-variable piecewise linear functions to give an inner approximation. One of the most important issues in piecewise linear programming is whether the adjacency criterion will hold. Fortunately, the adjacency criterion is automatically satisfied by the optimal solution in a convex separable program under which our problem is a subset (Simmons [39]). Thus, no binary variables are necessary to explicitly model adjacencies.
To implement the piecewise linear approximation, I pieces (I + 1 points) are used for each Z t,j,r . Herein, ν i j,r is defined as the coordinate of each piece variable, and λ i t,j,r as the new decision variable in the piecewise linearized function. The concave function is evaluated at each point g j,r (ν i j,r ).
Since g j,r (Z t,j,r ) = I ∑ i=0 g j,r (ν i j,r )λ i t,j,r ∀j ∈ IP, ∀t, ∀r ∈ R j , g j,r (Z t,j,r ) is replaced with the piecewise linearized summation I ∑ i=0 g j,r (ν i j,r )λ i t,j,r and Constraints (33), (34), and (35) are introduced, as shown below: Z t,j,r = I ∑ i=0 ν i j,r λ i t,j,r ∀j ∈ IP, ∀t, ∀r ∈ R j (34) 0 ≤ λ i t,j,r ≤ 1 ∀j ∈ IP, ∀t, ∀r ∈ R j The optimal solution to DPPNL can be obtained to any degree of accuracy by increasing the number of pieces. However, the formulation is entirely linear and can be solved efficiently.

Outer Approximation
Unlike the inner approximation method that solves one large DPPNL problem, the Kelley's cutting plane algorithm (KCP) may be used as an outer approximation technique. Outer approximation requires several iterations, as the problem size grows slowly with each iteration. A concave polyhedron is formed, iteration by iteration, to cover the original smooth concave function. The problem is first solved with the nonlinear constraint replaced by linear constraints, which are above the corresponding concave constraints and thus provide a feasible region larger than what is allowable. Constraints are generated by adding additional tangential lines to the concave function, based on the previous solution, at the point where the solution is infeasible for the original problem. The problem is solved again until the stopping condition is met; the solution should be feasible or the decrease in optimal objective function value must be less than the scalar ε. The cutting plane method thus develops an outer approximation. An optimal or near optimal solution can be reached by choosing an appropriate scalar, and the number of iterations can increase if more accuracy is desired.
Let N i t,j,r , Z i t,j,r denote the solution to iteration i. Let O i be the optimal objective function value to iteration i. Is g j,r (Z i t,j,r ) − (a i t,j,r + b i t,j,r Z i t,j,r ) < 0 and O i−1 − O i > ε ? If no: stop; the optimal solution has been found. If yes: I = i + 1 Write new constraint: Go to step 2 If no: stop. Convergence is guaranteed when a feasible solution exists and the nonlinear constraints are convex and continuously differentiable (Zangwill [40], Luenberger and Ye [41], pp. 463-465). Proposition 1. The I/O clearing function in Section 4.2.1 is concave and continuously differentiable. g j,r (Z t,j,r ) = C t,j V j,r Z t,j,r Z t,j,r + κ Taking the first derivative, g j,r (Z t,j,r ) = C t,j V j,r κ (Z t,j,r +κ) 2 , it can be seen that g j,r (Z t,j,r ) is continuously differentiable in the domain Z t,j,r ≥ 0 (since k > 0). 3 , is always negative. Therefore, g j,r (Z t,j,r ) is concave.
The M/D/1 clearing function in Section 4.2.2 is concave and continuously differentiable. g j,r (Z t,j,r ) = C t,j V j,r (Z t,j,r + 1 − Z t,j,r 2 + 1) Taking the first derivative, g j,r (Z t,j,r ) = C t,j V j,r (1 − 1/ Z t,j,r 2 + 1), it can be seen that g j,r (Z t,j,r ) is continuously differentiable in the domain Z t,j,r ≥ 0.
The second derivative, g j,r (Z t,j,r ) = −C t,j V j,r /( Z t,j,r 2 + 1) 3/2 , is always negative. Therefore, g j,r (Z t,j,r ) is concave.  g j,r (Z t,j,r ) = 1 − e −µZ t,j,r Taking the first derivative, g j,r (Z t,j,r ) = µe −µZ t,j,r , it can be seen that g j,r (Z t,j,r ) is continuously differentiable in the domain Z t,j,r ≥ 0. The second derivative, g j,r (Z t,j,r ) = −µ 2 e −µZ t,j,r , is always negative. Therefore, g j,r (Z t,j,r ) is concave. Note that (24) may be rewritten as follows: N t,j,r − C tj V j,r g j,r (Z t,j,r ) ≤ 0 ∀j ∈ IP, ∀t, ∀r ∈ R j From Propositions 1, 2, and 3, g j,r (Z t,j,r ) is concave and continuously differentiable in each of the three cases for the clearing function. Hence, the −g j,r (Z t,j,r )'s are differentiable convex functions, and the entire constraint set in DPPNL is convex. Luenberger and Ye [41] (pp. 463-465) show that the Kelley's cutting plane method converges to a limit solution for such a problem.

Computational Performance
In this section, results obtained by implementing the solution algorithms for DPPNL are discussed. All problems were run on a Pentium III Personal Computer with a Celeron CPU running at 850 MHZ.

Comparing Inner and Outer Approximation
For the sample case discussed in the paper, the RFQs in Table 8 are used to solve the problem using inner approximation. The optimal solution converges to 48,323 when the number of pieces is 55. Figure 5 shows the results.

Computational Performance
In this section, results obtained by implementing the solution algorithms for DPPNL are discussed. All problems were run on a Pentium III Personal Computer with a Celeron CPU running at 850 MHZ.

Comparing Inner and Outer Approximation
For the sample case discussed in the paper, the RFQs in Table 8 are used to solve the problem using inner approximation. The optimal solution converges to 48,323 when the number of pieces is 55. Figure 5 shows the results. The same problem is solved using the outer approximation method. An objective function value of 48,323 is reached from the ninth iteration (cut) onwards, with no significant change. The results are shown in Figure 6.  The same problem is solved using the outer approximation method. An objective function value of 48,323 is reached from the ninth iteration (cut) onwards, with no significant change. The results are shown in Figure 6. As seen above, both methods converge to an objective function value of approximately 48,323. Outer approximation provides a more accurate result, but the gap between the starting value and the final value is quite large. On the other hand, inner approximation gives a less precise result, but the difference between the starting value and the final value is much smaller.
It is interesting to notice that outer approximation solves the fully relaxed problem at first and then tightens the constraint set as close to the nonlinear constraint as possible, iteration by iteration. In other words, the method always underestimates the real optimal objective function value. Inner approximation does exactly the opposite, i.e., it always overestimates the real optimal objective function value. Thus, the inner approximation result is the upper bound and the outer approximation result the lower bound on the optimal solution.

Model Solution Time for Different Problem Sizes
Different sizes of problems in the linear and nonlinear DPP without lateness models were tested to calibrate the model convergence and running time. A separate two-level bill of material with four items in the problems was added to get 12 commodities and two final items in the network. The production horizon for each problem class was 10 time periods. The general clearing function was used for the problem. Table 9 summarizes the run time results obtained using inner and outer approximation. Outer approximation provides better performance than inner approximation in terms of CPU time as problem size increases, as can be seen in Table 9. As seen above, both methods converge to an objective function value of approximately 48,323. Outer approximation provides a more accurate result, but the gap between the starting value and the final value is quite large. On the other hand, inner approximation gives a less precise result, but the difference between the starting value and the final value is much smaller.
It is interesting to notice that outer approximation solves the fully relaxed problem at first and then tightens the constraint set as close to the nonlinear constraint as possible, iteration by iteration. In other words, the method always underestimates the real optimal objective function value. Inner approximation does exactly the opposite, i.e., it always overestimates the real optimal objective function value. Thus, the inner approximation result is the upper bound and the outer approximation result the lower bound on the optimal solution.

Model Solution Time for Different Problem Sizes
Different sizes of problems in the linear and nonlinear DPP without lateness models were tested to calibrate the model convergence and running time. A separate two-level bill of material with four items in the problems was added to get 12 commodities and two final items in the network. The production horizon for each problem class was 10 time periods. The general clearing function was used for the problem. Table 9 summarizes the run time results obtained using inner and outer approximation.
Outer approximation provides better performance than inner approximation in terms of CPU time as problem size increases, as can be seen in Table 9.
The augmented Lagrangian method implemented in MINOS 5.5 was used to evaluate the quality of the inner and outer approximation solutions. It may be seen in Table 10 that the solutions obtained by the approximations were reasonably close to the optimal objective function value obtained by the augmented Lagrangian method.

Conclusions
In this paper, we formulated a congestion model for of demand planning in supply chains that is general enough for the discrete-product or process industry. We illustrated the usefulness of the basic model in demand management via a series of scenarios in which the firm responds to RFQs, changes in demand, engineering changes, and due-date changes. The model's use in capacity management was illustrated with an example. We extended the model to incorporate congestion effects using clearing functions that work at the recipe level, which is very general. The resulting model was nonlinear, and we developed and tested two linear programming-based algorithms to solve the nonlinear model. The performance of the two algorithms, one based on inner approximation and the other on outer approximation, was very good. This would allow for large-scale practical application of the model. This research can be enhanced in many ways. One limitation of the model is that the number of recipes is a continuous variable, which is somewhat limiting for assembly systems. The queuing models considered are basic and assume that one recipe is run for a fairly long period of time for the queue to be stable. The concept of dual prices in Srinivasan et al. [17] and Kefeli and Uzsoy [27] can be applied in identifying not only capacity but also material bottlenecks. These limitations can be addressed in future research. Institutional Review Board Statement: It's not applicable for the study not involving humans or animals.
Informed Consent Statement: It's not applicable for the study not involving humans or animals.

Data Availability Statement:
The data presented in this study are available within the article.

Conflicts of Interest:
The authors declare no conflict of interest.