Control of Supply Chains

Abstract

This chapter aims to apply control principles in the discrete-time control of Supply Chains. The primary objective of the control is to keep the inventory levels (state variables) steady at their predetermined values and reduce any deviations to zero in the shortest possible time. The disturbances are induced by demand deviations from the planned/anticipated levels. The replenishment flows are the control variables. Thus, the control action is very similar to a “Linear Regulator with zero set-point”. A novel development in this chapter is the use of direct Operator methods to solve the system Difference Equations, thereby obviating the need for Z-Transforms, block diagrams and transfer functions of classical control theory. This chapter provides a novel application of control theory as well as an easier method of solution.

1. Introduction

In the field of business, one of the most important constituents of a business system, and one which can give a business a cutting-edge over the others, is the supply chain (SC) of the business, and its effective operation and management.

A fundamental strategy associated with supply chains is that of a ‘Responsive Chain’. A responsive chain focuses on its ability to respond quickly to demand changes and meet them within the shortest possible time. A responsive chain, by its very definition, has necessarily to be able to respond swiftly to unanticipated changes in demand, and this is determined largely by the dynamics of the system. Under normal circumstances, good responsiveness can be achieved through the maintenance of adequate pipeline inventories throughout the system, and one of the factors impacting this decision in a major way, is again the dynamics of the system. Hence, given the impact of the dynamics of the SC system on its design and operation as mentioned above, it is imperative to include the tenets of dynamic analysis in its design and operation. Thus system-dynamic methods are of significant value and utility in both SC design and SC operation. The design aspects pertain to the setting of inventory levels and the operational aspects to choosing the appropriate system controls for good operational performance.

Thus, we are led naturally to the use of system-dynamic methods and concepts of control theory for effective control of a supply chain. In the following sections we explore the dynamic modeling and analysis concepts in effectively controlling a supply chain that would yield good performance characteristics.

2. Supply chain dynamics and notation

Supply chain dynamics essentially deals with the dynamic behavior and temporal variation of the inventories and flows in the system over time when subjected to demand disturbances.

We now look at a simple three stage serial supply chain, a schematic diagram of which is given below in Figure 1. Each stage represents an inventory storage facility in the chain. Stage 1 represents the raw material input storage to a manufacturing plant, and stage 2 the finished goods inventory storage at the manufacturing plant. Stage 3 represents the finished goods warehouse (W/H) at the downstream end of the chain which meets the demand for the finished product, and from which material is shipped out to customers. The inventory levels at each stage are the state variables of the system, and the material flows between the stages the control variables.

Figure 1.

A three-stage serial supply chain.

In constructing a dynamic model of a supply chain system, we adopt the standard schematic diagram (Figure 1), system notation, and model variables and equations below, as commonly found in control theory literature [1, 2, 3, 4, 5, 6, 7, 8, 9].

We use a discrete time representation of the SC system, as this is more in keeping with the prevailing practices in the SC industry, wherein the inventory levels (state variables) are recorded at the end of each period or day. In most cases the state variable records are updated at the end of each day. Hence, we take our period to be a single day. However, we need to also emphasize that this need not always be so and would be dependent entirely on the convenience of the SC practitioners. And hence the state variables are recorded at epochs corresponding to the end of each period. The flow variables however are aggregated and taken to occur within and up to the end of each period.

The inventory and flow variables are subscripted to indicate the stages in the chain. We now give the detailed representation below:

$y_i\left(k\right)$: inventory at stage i of the chain at time epoch k

with: i = 1, representing (the raw material) the upstream end of the production facility.

i = 2 representing (the finished goods) the downstream end of the production facility.

i = 3 representing (the finished goods) the warehouse.

$q_i^{\prime }\left(k\right)$: material flow in period (k-1,k] into stage i of the chain

$r_3^{\text{'}}\left(k\right)$: finished goods demand at warehouse in period (k-1,k]

The dynamic equations of the system are written using deviation variables, as under:

$$x_i\left(k\right)=y_i\left(k\right)-y_i^0\left(k\right)\text{for}\mathrm{i}=\mathrm{1,2,3}$$

$$q_i\left(k\right)=q_i^{\prime }\left(k\right)-q_i^0\left(k\right)\text{for}\mathrm{i}=\mathrm{1,2,3}$$

$$r_3\left(k\right)=r_3^{\prime }\left(k\right)-r_3^0\left(k\right)$$

where,

$y_i^0\left(k\right)$ is the desired or planned inventory level at time k

$q_i^0\left(k\right)$ is the desired or planned flow in period (k-1,k]

$r_3^0\left(k\right)$ is the anticipated or forecasted demand at the warehouse in period (k-1,k], based on which the desired inventory levels in the chain have been set.

$x_i\left(k\right)$ is the inventory deviation at time k.

$q_i\left(k\right)$ is the material flow deviation in period (k-1,k].

$r_3\left(k\right)$ is the deviation in demand observed at the warehouse in (k-1, k].

Under consistent units (in equivalent units of finished goods) for all inventory and flows in the system, the dynamic equations of the system are written as:

$$x_i\left(k+1\right)=x_i\left(k\right)+q_i\left(k+1\right)-q_{i+1}\left(k+1\right)$$E1

where for i = 3, $q_4\left(k\right)=r_3\left(k\right)$, the demand outflow from the warehouse. The system behavior is controlled by the replenishment policies followed in the system.

Now, since the demand is a stochastic variable, it can be split into two components, viz., a mean demand component, and a stochastic component which represents the random variations over and above the mean demand. The mean demand is what is predicted using forecasting techniques and is the planned offtake and the demand that the system is designed to meet in the normal course. Since prediction can never be exact, the residual variation is the stochastic component.

In such stochastic systems, the demand has an additional stochastic term, which is a white noise term, given by $\epsilon \left(k\right)\sim \mathit{WN}\left(0,{\sigma }^2\right)$ for all k. The sequence of stochastic disturbance components over each period is an independent and identically distributed sequence of $N\left(0,{\sigma }^2\right)$ random variables.

The standard initial conditions for the system are: $x_i\left(k\right)=0$ for all $k\le 0$ and $r_3\left(k\right)=0$ for all $k\le 1$, i.e. the system is at zero deviation at time k = 0, and the first deviation in demand is felt at the end of the first period, at k = 1.

The demand disturbance could be of the following Types:

1. a sudden shock demand increase, represented by a Dirac delta input function.

2. a sudden and sustained increase in demand, represented by a step input.

3. a demand with an increasing trend, represented by a ramp input function.

4. a demand with second order (quadratic), third order (cubic), or higher-order trends represented by higher-order polynomial input functions.

5. a demand with seasonality, represented by a sinusoidal input function.

6. A random component which is represented by a White Noise process.

The first five components pertain to and define the mean demand disturbance, while the last represents the residual random variations over and above the mean demand disturbance.

The demand can have any combination of, or even all of the above components in the mean demand, and, of course, the additional stochastic component represented by a White Noise process. Such demands are often seen in supply chain warehouses.

Thus, a demand disturbance at the downstream end (at the warehouse) provides the perturbation to the system. The system is controlled through regulation/control of the replenishment flows which are the control variables in the system.

Now from the above, we can see that the system response will have two components, the Mean Response, and a stochastic component of the response. It is the mean response that characterizes the system behavior, while the stochastic component is characterized by the inventory variance.

3. Dynamic supply chain performance metrics and control action triggers

3.1 Dynamic performance metrics

The performance metrics of a supply chain that would be of interest to us from a dynamic performance point of view and control system design are explained below.

Since the demand input to the system suddenly increases, we can expect the system response to lag the demand and to fluctuate, as the system scrambles to catch up with the increased demand. Accordingly, the key performance indicators that would be of importance and interest to us from the point of view of both design and operation are the following [1, 2, 3, 4, 5, 6, 7, 8, 9]:

1. Permanent depletion of the inventory level (the offset), if any.

2. The trough value or the lowest dip in the inventory levels (the undershoot).

3. The amplitude of fluctuations of the inventory.

4. The center-line about which fluctuations occur.

5. The fraction of time the inventory level stays depleted (in the negative region).

6. The limiting inventory variance.

The first indicator shows whether the system is able to catch up with the demand and whether it is ultimately restored to its original level, or not. The second indicator, the trough value, represents the lowest point or value that the inventory level is likely to touch, and impacts the base stock levels that would have to be carried to maintain uninterrupted material flows in the system, (adequate pipeline inventories). The third defines the magnitude of fluctuations of the inventory levels, and high values would naturally be undesirable; and we would like to damp them down to zero as quickly as possible. The fourth indicates the center-line about which fluctuations occur and is indicative of the average inventory levels, which we would like to keep at positive levels for comfortable operation. While the fifth can be taken to be a surrogate measure of the stock-out risk, and the higher the fraction of time in the negative region (with a depleted inventory level), the higher could be taken to be the implied stock-out risk. The last, the inventory variance is the variation that we could expect even after the system has been restored to its original operating levels and is commonly taken as a measure of the robustness of the system to random demand variations.

3.2 Control action triggers

The two most common triggers for initiation of replenishment control action in the system are [1, 2, 3, 4, 5, 6, 7, 8, 9]:

1. The inventory levels at the W/H, as found in logistic systems and warehouses,

2. The demand variations at the W/H, as in electronic-data-interchange systems.

Thus, the control flows into the W/H are set to a function of the latest available inventory deviations and latest available Demand deviations. Thus, we have,

$$q_3\left(k+1\right)=f\left[x_3\left(k-1\right),x_3\left(k-2\right),..,x_3\left(k-r\right);r_3\left(k-1\right),r_3\left(k-2\right),\dots r_3\left(k-p\right)\right]$$E2

We discuss these points in more detail for a single stage system first.

3.3 Single stage control notation

To differentiate between the standard abbreviations in conventional control theory, we adopt the following notation used to indicate the type of control used.

P, PI, PID would represent Proportional, Proportional-integral, and Proportional-integral-derivative control as usual. Additionally, MA denotes a ‘Moving Average’ type of control (explained in detail subsequently).

An ‘I’ within parenthesis would represent ‘inventory-triggered’ control, while ‘D’ within parenthesis would denote a ‘Demand-triggered’ control. Also, ‘ID’ within parenthesis would denote a control which would have both inventory-triggered and Demand-triggered components. We give examples below.

Thus P(D) (pronounced as “P of D”) denotes Demand-triggered Proportional control, PI(I) denotes Inventory-triggered PI control. Similarly, PID(ID) (pronounced as “PID of ID”) denotes PID control with inventory-triggered and Demand-triggered components.

We also could have cases of multiple inventory triggers and multiple-demand triggers. These will be denoted as follows:

$X\left(I_n{D}_m\right)$ control would indicate a control of type X (P, PI, PID, MA, Composite) with n-inventory trigger terms and m-demand trigger terms. Thus, the control $\mathit{PD}\left(I_5{D}_3\right)$ (pronounced as “PD of I5D3”) would be Proportional-derivative control with 5 inventory triggers and 3 demand triggers.

We now first look at a single stage system below which we take to be the warehouse end of a SC.

4. Single stage control

4.1 The two response components

Since the demand has both a deterministic component as well as a stochastic component, the response of the system can also be broken down into two components, viz., a deterministic part which is the mean response, and a random or stochastic part characterized by the inventory variance. It is the mean response that portrays the behavior of the system and yields the performance indicators of the system. Whereas, the stochastic part, represents the random fluctuations that have to be accounted for even after the system is brought under full control.

From the above discussion we can see the close parallel with conventional feedback control theory. Here the information about the warehouse inventory (the state variable) is fed back to the system for initiating replenishment flow control action. Additionally, we can also have controls wherein the disturbance is also directly fed back to the system for initiation of control action.

4.2 Single stage system controls and transportation lags

In close parallel with classical feedback control theory, the control flows are functions of the state variables and the demand or input perturbation to the system. The types of functions used also closely parallel classical control theory. And hence we have P, PI, and PID controls.

Additionally, we also have Moving Average (MA) controls, wherein the control flow is set to a weighted Moving Average of the latest available inventory deviations of up to ‘r’ periods back. The parameter ‘r’ is termed the ‘Order’ of the Moving Average.

We first look at Proportional Controls.

Now, for Proportional Controls of the P(I) type, the control flows into the warehouse in stage 3 of the chain would be given by:

$$q_3\left(k+1\right)=K_3{x}_3\left(k-1-l_3\right)$$E3

where $l_3$ is the transportation lag in the flows into the warehouse from the upstream unit of the system, i.e., the finished goods inventory at the manufacturing plant in our case. And $K_3$ is the constant of proportionality between the control flow and the latest available inventory deviation based on which the order is initiated (and hence ‘Proportional’ to the error in conventional control theory).

The development herein corresponds to the ‘Regulator Problem’ with ‘set point’ of zero. And hence control of a Responsive Chain parallels the regulator problem in conventional control theory.

In conventional modeling of SCs, the lag is taken as the number of periods strictly between the period of order initiation and the period of arrival of the consignment, not including the period of arrival of the consignment. Thus, the lag is taken to be zero if the replenishment consignment arrives in the period immediately succeeding the period of order initiation. Instantaneous replenishment is not envisaged and is very rare in SC contexts. Arrival of the consignment within the same period of order initiation is also not envisaged and is very rare in such contexts. The earliest arrival of an ordered consignment is taken to be the immediately succeeding period, for which we take the lag as zero.

This convention is based on what is normally followed in the industry and practice, as well as the literature on dynamic modeling of SCs.

Hence in our further development of dynamic models of SCs we take the lag as zero if an order initiated in period k i.e., in the interval (k – 1, k] arrives in period (k + 1), i.e., in the interval (k, k + 1].

For a transportation lag of ‘l’ periods, an order placed in the interval (k – 1, k] would arrive in the interval (k + l, k + l + 1].

Hence for the warehouse under zero lag the control flows for P(I) control would be set as:

$$q_3\left(k+1\right)=K_3{x}_3\left(k-1\right)$$E4

Now it is to be noted that the order is initiated in period k i.e., in the interval (k – 1, k] based on the latest fully observed inventory deviation which in this case would be that at the start of period k, which is the inventory recorded at the end of period (k-1) i.e., $x_3\left(k-1\right)$. We take it that since the inflows in period k would be the aggregate of flows in the interval (k – 1, k], and the order is initiated within the period k (and not at time point k), the latest available fully observed inventory level would be $x_3\left(k-1\right)$ and not $x_3\left(k\right)$ because $x_3\left(k\right)$ would not be available to the order initiator within period k. The value of $x_3\left(k\right)$ is recorded at the closure of period k, after the cessation of all activities including the action of ordering, of period k. Thus, replenishment orders need to be initiated within and before the closure of a period and not at the end of a period.

Thus, warehouse records would be updated at the closure of the day’s operations and would show the closing inventory and the replenishment orders placed during the day. The replenishment orders would have been placed based on the previous day’s closing stock and would arrive during the course of the next day if the lag is zero.

This is the convention that we will follow in the further development of the models.

We next look at the next type of control, which is the PI(I) control.

For the Proportional-Integral (PI(I) type) control case with zero lag, the control flows at the warehouse would be set as under:

$$q_3\left(k+1\right)=K_3{x}_3\left(k-1\right)+K_c\sum _{m=0}^{k-1}{x}_3\left(m\right)$$E5

where the second term is the integral term in our discrete-time system, and K_c is the proportionality constant (gain term) factor of the integral of the error.

Next we have the PID(I) control wherein the control flows into the warehouse under zero lag would be set to:

$$q_3\left(k+1\right)=K_3{x}_3\left(k-1\right)+K_c\sum _{m=0}^{k-1}{x}_3\left(m\right)+K_d\left(x_3\left(k-1\right)-x_3\left(k-2\right)\right)$$E6

where the last term represents the derivative term in our discrete-time system, and K_d is the proportionality constant factor (gain term) of the derivative component of the control.

Next we have the Moving Average (MA) type of control (MA(I) type), for which the control flows into the warehouse under zero lag would be set to:

$$q_3\left(k+1\right)=K_3^1{x}_3\left(k-1\right)+K_3^{2}x\left(k-2\right)+K_3^{3}x\left(k-3\right)+\dots \dots +K_3^{r}x\left(k-r\right)$$E7

where the r is the order of the moving average, the Ks are the control parameters (the weights) of the MA terms. Thus, the control flow is set to a weighted moving average of the latest available fully observed inventory levels up to r period back.

The above controls discussed above are the conventional inventory-triggered schemes. In all these types of controls, we could additionally have demand-triggered terms also, like the PI(ID), PID(ID), MA(ID) etc.

We discuss some of them below for single stage systems.

5. Solution of single stage systems

Firstly, we note herein that in solving for the response of a supply chain system, we will not use the Z-transform nor block diagrams and transfer functions as in conventional control theory. Rather we will work directly on the system difference equation in the time domain itself. And instead, we will use direct Operator methods to obtain the system solution and response (rather than the transformed equations and inverse transforms). This is one of the advantages of this modeling paradigm.

Another advantage of the type of discrete time modeling taken up here is that transportation lags do not result in differential-delay equations as in conventional continuous-time control theory. Rather, transportation lags would increase the order of the resulting system difference equation, which is expected to be easier to solve than differential-delay equations.

We now take up the simplest form of control which is the P(I) control and illustrate the formulation of the system equation and its solution method.

5.1 P(I) control under zero lag

The control is an inventory-triggered Proportional control. And we take the Replenishment Lag = 0 in the simplest case.

The flow balance equation for the warehouse is given by:

$$x_3\left(k+1\right)=x_3\left(k\right)+q_3\left(k+1\right)-r_3\left(k+1\right)$$E8

The control flows into the warehouse are given by:

$$q_3\left(k+1\right)=K_3{x}_3\left(k-1\right)$$E9

It is to be noted that the value of $K_3$ < 0, since the control flows and inventory deviation have to be of opposite sign, i.e., when the inventory deviation is (−)ve (inventory is at a depleted level) then the flow deviation has to be (+)ve to enable more/extra flow into the warehouse to make up the inventory shortfall. Likewise, if the inventory deviation is (+)ve (extra inventory in stock) then the flow deviation has to be (−)ve (i.e., reduced flow) to reduce the inflow into the warehouse to maintain inventory levels.

Thus, we can clearly see that the controls are of the ‘feedback’ type, and they seek to keep the inventory deviation at zero level, which is just the’ Linear Regulator with zero set-point’ in standard control theory.

Thus, substituting for the control flow into the flow balance eqn. Above, yields the system-dynamic eqn. For the warehouse as:

$$x_3\left(k+1\right)\equiv x_3\left(k\right)+K_3{x}_3\left(k-1\right)-r_3\left(k+1\right)$$E10

$$\text{Or equivalently},x_3\left(k+1\right)-x_3\left(k\right)-K_3{x}_3\left(k-1\right)\equiv -r_3\left(k+1\right)$$E11

The above is the deterministic part of the system equation. Since demand is a stochastic variable with a stochastic component, the complete system equation is given by:

$$x_3\left(k+1\right)-x_3\left(k\right)-K_{3}x\left(k-1\right)\equiv -r_3\left(k+1\right)-\epsilon \left(k+1\right)$$E12

which has both parts. To solve the system equation completely, we split it into its two components and solve for each of the components separately. We hence solve the following two equations, one each for the deterministic part and the stochastic part.

$$x_3^{\det }\left(k+1\right)-x_3^{\det }\left(k\right)-K_3{x}_3^{\det }\left(k-1\right)\equiv -r_3\left(k+1\right)\text{for the deterministic part},\text{and},$$E13

$$x_3^{\mathit{\text{stoc}}}\left(k+1\right)-x_3^{\mathit{\text{stoc}}}\left(k\right)-K_3{x}_3^{\mathit{\text{stoc}}}\left(k-1\right)\equiv -\epsilon \left(k+1\right)\text{for the stochastic part}.$$E14

The first is a deterministic Linear Difference Equation, and the second, a Stochastic Linear Difference Equation (SDE).

Both are second-order Linear Difference Equations (LDEs) in the state variable $x_3\left(k\right)$ and can be solved by direct Operator methods for any type of input or disturbance term on the RHS of the system LDE.

An excellent treatment of difference calculus and solution methods for LDEs is given in [10]. We follow the methods given therein.

In order to solve the LDE, we first introduce the Forward Shift Operator E as under: ${\mathit{Ex}}_3\left(k\right)\equiv x_3\left(k+1\right)$, with the important property that: $E\left[{\mathit{Ex}}_3\left(k\right)\right]=E^2{x}_3\left(k\right)=x_3\left(k+2\right)$. Before using the Operator, we first write the LDE in standard form as under: (the lowest time value is taken as k):

$$x_3\left(k+2\right)-x_3\left(k+1\right)-K_3{x}_3\left(k\right)\equiv -r_3\left(k+2\right)$$E15

Now using the forward Shift Operator E, we can write the LDE in Operator form as:

$$\left[E^2-E-K_3\right\}{x}_3^{\det }\left(k\right)\equiv -r_3\left(k+2\right)\text{for the deterministic part},\text{and},$$E16

$$\left[E^2-E-K_3\right\}{x}_3^{\mathit{\text{stoc}}}\left(k\right)\equiv -\epsilon \left(k+2\right)\text{for the stochastic part}.$$E17

5.1.1 The mean response: solution of the deterministic LDE

We first look at the deterministic part of the solution below, which will yield the mean response. For notational convenience, we drop the superscript on $x_3\left(k\right)$.

$$\left[{\mathrm{E}}^2-\mathrm{E}-{\mathrm{K}}_3\right]{\mathrm{x}}_3\left(\mathrm{k}\right)\equiv -{\mathrm{r}}_3\left(\mathrm{k}+2\right)$$

Now this is an LDE of order two (the order being the highest power of the Operator E).

We write the Characteristic Equation of the LDE as [10]:

$${\alpha }^2-\alpha -K_3=0$$E18

to determine the characteristic roots of the LHS Operator.

$$\text{From which we get the roots}\mathrm{as}:{\alpha }_{1/2}=\frac{1\pm \sqrt{1+4K_3}}{2}.$$E19

The stability of the system is entirely controlled by the roots of the LHS Operator. And elementary analysis leads to the following stability conditions:

$K_3\ge 0$: instability

$-1/4\le K_3<0$: stability with non-oscillatory response

$-1\le K_3<-1/4$: stability with oscillatory behavior

$K_3<-1$: instability with oscillatory behavior

Now we take up the solution of the system LDE for a unit step increase in demand, i.e.,

$$r_3\left(k+1\right)=1,\forall k\ge 0$$E20

Substituting for the demand disturbance in the system equation yields the LDE:

$$\left[E^2-E-K_3\right]x_3\left(k\right)\equiv -1$$E21

which we can call the “Original Non-Homogeneous Eqn.” (O-NHE).

We first look at the solution of the homogeneous LDE (i.e., with RHS = 0). The homogeneous LDE is:

$$\left[E^2-E-K_3\right]x_3\left(k\right)\equiv 0$$E22

which upon factoring the LHS Operator can be written as: (${\lambda }_1,{\lambda }_2$ being the distinct roots of the LHS Operator):

$$\left[\left(E-{\lambda }_1\right)\left(E-{\lambda }_2\right)\right]x_3\left(k\right)\equiv 0$$E23

which has the solution as: $x_3\left(k\right)\equiv C_1{\lambda }_1^k+C_2{\lambda }_2^k$ for the case of distinct roots.

For a repeated root, the solution is given by: $x_3\left(k\right)\equiv \left(C_0+C_{1}k\right){\lambda }^k$ for a repeated root of algebraic multiplicity two.

Now that we have the solution of the homogeneous LDE, we next look for a particular solution of the Original Non-Homogeneous LDE (O-NHE).

A standard method of solution of the Non-homogeneous eqn. is by the ‘Annihilator Method’ [10].

We look for the Operator that annihilates the RHS terms of the O-NHE, say A(E). Then operating by the Annihilator on both sides of the O-NHE yields:

$$A\left(E\right)\left[\left(E-{\lambda }_1\right)\left(E-{\lambda }_2\right)\right]x_3\left(k\right)\equiv A\left(E\right)r_3\left(k+2\right)\equiv 0$$E24

which is a homogeneous LDE albeit of a higher order, but which can be solved by factorizing the LHS Operator. As an example, in our case of a unit step disturbance, $r_3\left(k\right)\equiv -1,\forall k\ge 1$.

And the Annihilator is given by:

$$A\left(E\right)\equiv \left(E-1\right),\because \left(E-1\right)r_3\left(k\right)\equiv -\left(E-1\right)\left(1\right)=0$$

And hence the equivalent Homogeneous LDE is given by:

$$\left[\left(E-1\right)\left(E-{\lambda }_1\right)\left(E-{\lambda }_2\right)\right]x_3\left(k\right)\equiv 0$$

which has the solution:

$$x_3\left(k\right)\equiv D\left(1^k\right)+C_1{\lambda }_1^k+C_2{\lambda }_2^k\equiv D+C_1{\lambda }_1^k+C_2{\lambda }_2^k\text{for the case of distinct roots}.$$E25

$$x_3\left(k\right)\equiv D{\left(1\right)}^k+\left(C_0+C_{1}k\right){\lambda }^k\text{for the case of repeated roots}$$E26

Now we note that the O-NHE being of order two, will admit only two undetermined constants. The above solution, however, has three undetermined constants. The third was introduced by us due to the Annihilator. Hence to determine the extra constant D in the solution, we substitute the solution into the O-NHE to determine D. Thus, and noting that the terms involving the roots of the LHS operator of the O-NHE are precisely the homogeneous solution terms of the O-NHE, we only need substitute the extra terms introduced by the Annihilator into the O-NHE. Hence, we have:

$$\left[\left(E-{\lambda }_1\right)\left(E-{\lambda }_2\right)\right]D\equiv \left[E^2-E-K_3\right]D=-1$$

Noting that E(D) = D itself ($\because E\left(D\right)\equiv D\left(k+1\right)\equiv D,\forall k$), we obtain $D\left(k+2\right)-D\left(k+1\right)-K_{3}D\equiv D-D-K_{3}D=-1$ from which we readily obtain

$$D=1/K_3\text{which is}\left(-\right)\mathrm{ve}\mathrm{due}\text{to the}\left(-\right)\mathrm{ve}sign\text{of}{K}_3.$$E27

Thus, the full solution is given by:

$$x_3\left(k\right)\equiv 1/K_3+C_1{\lambda }_1^k+C_2{\lambda }_2^k\text{for the case of distinct roots},\text{and}$$E28

$$x_3\left(k\right)\equiv 1/K_3+\left(C_0+C_{1}k\right){\lambda }^k\text{for the case of repeated roots}.$$E29

Now, for stable solutions, we will have $\left|{\lambda }_1\right|<1,\left|{\lambda }_2\right|<1$, and hence the terms involving the roots of the LHS Operator decay to zero for large k, leaving the ‘Offset’ term as $1/K_3$ which is the Offset value. The Offset value is (−)ve because of the (−)ve sign of $K_3$.

The Damping rate, which is the rate at which the fluctuations decay to zero are given by the magnitudes of the roots of the LHS Operator $\left({\left|{\lambda }_1\right|}^k,{\left|{\lambda }_2\right|}^k\right)$, which can clearly be seen to decrease in magnitude with decrease in magnitude of $\left|K_3\right|$. However, the Offset value increases in magnitude with decrease of $K_3$.

We take the case of distinct roots first, say for a value of $K_3=-1/2$, in the stable region.

Hence, we have the system equation as:

$$\left[E^2-E-K_3\right]x_3\left(k\right)\equiv \left[E-E-\left(-1/2\right)\right]x_3\left(k\right)\equiv -1$$E30

The roots of the LHS Operator are readily obtained as: $\lambda =\left(1/2\right)\pm \left(1/2\right)j,\left|{\lambda }_1\right|<1,\left|{\lambda }_2\right|<1,j=\sqrt{-1}$, which can be written in the form.

$\lambda =\rho e^{\pm \mathit{j\theta }}\equiv \left(1/\sqrt{2}\right)e^{\pm \mathit{j\pi }/4}$, which hence yields the homogeneous solution as:

$\left\{A\mathrm{Cos}\left(\mathit{k\pi }/4\right)+B\mathrm{Sin}\left(\mathit{k\pi }/4\right)\right\}{\left(1/\sqrt{2}\right)}^k$, and hence the full solution as:

$$x_3\left(k\right)\equiv \left(1/K_3\right)+\left\{A\mathrm{Cos}\left(\mathit{k\pi }/4\right)+B\mathrm{Sin}\left(\mathit{k\pi }/4\right)\right\}{\left(1/\sqrt{2}\right)}^k$$E31

$$\equiv -2+\left\{A\mathrm{Cos}\left(\mathit{k\pi }/4\right)+B\mathrm{Sin}\left(\mathit{k\pi }/4\right)\right\}{\left(1/\sqrt{2}\right)}^k,$$E32

which shows a Damping Rate of the Order of $O{\left(1/\sqrt{2}\right)}^k=O{\left(0.707\right)}^k$, with an Offset of - 2 units. The response if plotted in Figure 2, and the ‘Undershoot’ is obtained as −2.5 units.

Figure 2.

P(I) control.

We next take a value of $K_3=-1/4$ for the case of repeated roots, again in the stable region. The roots of the LHS Operator are obtained as: ${\lambda }_1={\lambda }_2=\lambda =1/2$, and the full solution as:

$$x_3\left(k\right)\equiv \left(1/K_3\right)+\left\{C_0+C_{1}k\right\}{\left(1/2\right)}^k\text{which for}{K}_3=-1/4\text{yields}$$E33

$$x_3\left(k\right)\equiv -4+\left\{C_0+C_{1}k\right\}{\left(1/2\right)}^k,$$E34

which shows a Damping Rate of the Order of $O{\left(1/2\right)}^k=O{\left(0.5\right)}^k$, but with a high Offset value of - 4 units, and an Undershoot of - 4 units.

Thus, while we would like the oscillations to be damped out rapidly, this would compromise on the Offset value, which would impact the base stock requirements of the system.

Thus, for the practitioner, the implications are quite clear: there is a trade-off between stability and rapid damping on the one hand, and the base stock requirements of the system (for low stock-out risk) on the other. The higher the damping (stability) required, the higher would be the base-stock requirements to keep stock-out risk low; and alternatively, the higher the (stability) damping achieved, the higher would be the stock-out risk at fixed base-stock levels.

We next obtain the undetermined constants in the general solutions above, using the Initial Conditions (ICs) of the system, as under:

The standard ICs of the system are as: $\left\{x_3\left(k\right)\equiv 0,r_3\left(k\right)\equiv 0,\forall k\le 0\right\}$.

Now the system LDE: $\left[E^2-E-K_3\right]x_3\left(k\right)\equiv -r\left(k+2\right)\equiv -1,\forall k\ge 0$ is valid for $\forall k\ge 0$.

And hence the ICs for our warehouse system can be obtained from the system equation itself using the standard system ICs and yields: $\left\{x_3\left(0\right)=0,x_3\left(1\right)=-1\right\}$. Substituting these ICs into the solutions above yields the full solutions as:

$$\text{For}{K}_3=-1/2:x_3\left(k\right)\equiv -2+2\mathrm{Cos}\left(\mathit{k\pi }/4\right){\left(1/\sqrt{2}\right)}^k,k\ge 0$$E35

$$\text{For}{K}_3=-1/4:x_3\left(k\right)\equiv -4+\left(4+2k\right){\left(1/2\right)}^k,k\ge 0.$$E36

We additionally examine the case for $K_3=-1$ (the maximum possible magnitude for stability), the response is given by

$x_3\left(k\right)\equiv -1+\left[A\mathrm{Cos}\left(\mathit{k\pi }/3\right)+B\mathrm{Sin}\left(\mathit{k\pi }/3\right)\right]{\left(1\right)}^k,k\ge 0$, which using the LDE ICs yields:

$$x_3\left(k\right)\equiv -1+\mathrm{Cos}\left(\mathit{k\pi }/3\right)+\left(1/\sqrt{3}\right)\mathrm{Sin}\left(\mathit{k\pi }/3\right),k\ge 0,$$E37

which can be simplified to.

$$x_3\left(k\right)=\left\{-1-\left(2/\sqrt{3}\right)\mathrm{Sin}\left(\left(k-1\right)\pi /3\right)\right\}H\left(k-1\right),k\ge 1,x_3\left(0\right)=0.$$E38

where H(.) is the unit Heaviside step function and yields a sinusoidal pattern with a center-line of – 1 and constant amplitude of $2/\sqrt{3}$ and the maximum negative deviation in inventory, the undershoot, equal to −2.

This last case is that of marginal stability characterized by constant amplitude perpetual oscillations which never die down to zero.

The responses for the three cases above are plotted in Figure 2.

5.1.2 The limiting inventory variance: solution of the SDE

We next look at the determination of the limiting inventory variance, which is a measure of the variation that we could expect even after the system (mean inventory level) has been restored to its original value.

Also, for the behavior of the mean response as well as our inferences from it to be meaningful, it is necessary that the limiting inventory variance be bounded and finite. We can then expect the inventory levels to be within the band given by: $x_3{\left(k\right)}^{\det }\pm 3\sqrt{{lim}_{t\to \infty }var\left(x_3\left(k\right)\right)}$, where $x_3{\left(k\right)}^{\det }$ is the mean response that has been determined above from the deterministic system LDE.

In order to determine the limiting inventory variance, we make use of the stochastic component of the system equation and determine the stochastic part of the response. For our system the Stochastic LDE (or SDE) is as under:

$$\left[E^2-E-K_3\right]x_3\left(k\right)\equiv -\epsilon \left(k+2\right),$$E39

where the term on the RHS of the SDE is the random variation represented by a White Noise Process, with $\epsilon \left(k\right)\sim \mathit{WN}\left(0,{\sigma }^2\right)$.

We can note that the LHS Operator is again the same as in the deterministic part of the system LDE that we have solved for above. Hence the roots of the LHS Operator remain unaltered in the SDE also.

Now following the method used in [11], we can note that if unity is not a root of the LHS Operator, then the SDE admits as a solution, an infinite weighted Moving Average representation in terms of the White Noise disturbance terms as under:

$$x_3{\left(k\right)}^{\mathit{\text{stoc}}}\equiv \sum _{l=0}^{k-1}{\beta }_l\epsilon \left(k-l\right)$$E40

where the ${\beta }_l$ s are the weights.

And hence the stochastic part of the solution of the system eqn. can be written as an infinite linear combination of the white noise disturbance terms.

Now since the individual white noise terms are Uncorrelated and Normal, i.e., with $\epsilon \left(k\right)\sim N\left(0,{\sigma }^2\right)$, and, with $\mathit{Cov}\left(\epsilon \left(i\right),\epsilon \left(j\right)\right)={\delta }_{\mathit{ij}}{\sigma }^2,{\delta }_{\mathit{ij}}=\left[\begin{array}{c}0,i\ne j\\ 1,i=j\end{array}\right.$, the limiting variance of x(k) is given by:

$${lim}_{k\to \infty }var\left(x_3\left(k\right)\right)={lim}_{k\to \infty }var\left(\left\{\sum _{l=0}^{k-1}{\beta }_l\epsilon \left(k-l\right)\right\}\right)={lim}_{k\to \infty }\left\{\sum _{l=0}^k{\beta }_l^2\right\}{\sigma }^2$$E41

In order to solve for the weighting terms, the βs, we substitute the solution into the system SDE above, as under:

$$\left[E^2-E-K_3\right]x_3{\left(k\right)}^{\mathit{\text{stoc}}}\equiv -\epsilon \left(k+2\right),\text{valid for}\mathrm{all}k\ge 0,$$E42

which is

$$\left[E^2-E-K_3\right]\left\{\sum _{l=0}^{k-1}{\beta }_l\epsilon \left(k-l\right)\right\}\equiv -\epsilon \left(k+2\right),k\ge 0.$$E43

Another and more convenient way to write the SDE is to use the Backward Shift Operator L, defined by:

$$L\mathrm{x}\left(\mathrm{k}\right)=\mathrm{x}\left(\mathrm{k}-1\right),\mathrm{i}.\mathrm{e}.,L=E^{-1},E=L^{-1}.$$E44

The SDE for the βs becomes:

$$\left[1-L-K_3{L}^2\right]\left\{\sum _{l=0}^{k-1}{\beta }_l\epsilon \left(k-l\right)\right\}\equiv -\epsilon \left(k\right),k\ge 0,$$E45

which is

$$\left[1-L-K_3{L}^2\right]\left\{{\beta }_0\epsilon \left(k\right)+{\beta }_1\epsilon \left(k-1\right)+{\beta }_2\epsilon \left(k-2\right)+{\beta }_3\epsilon \left(k-3\right)\dots \right\}\equiv -\epsilon \left(k\right),k\ge 0$$E46

Now comparing coefficients of ε(k) for each k, yields the system of equations as below:

LDE term	Coefft of ε(k)	Coefft of ε(k-1)	Coefft of ε(k-2)	Coefft of ε(k-3)	Coefft of ε(k-4)	Coefft of ε(k-5)	….	Coefft of ε(k)
1	${\beta }_0$	${\beta }_1$	${\beta }_2$	${\beta }_3$	${\beta }_4$	${\beta }_5$	…..	${\beta }_k$
- L	0	$-{\beta }_0$	$-{\beta }_1$	$-{\beta }_2$	$-{\beta }_3$	$-{\beta }_4$	…..	$-{\beta }_{k-1}$
$-K_3{L}^2$	0	0	$-K_3{\beta }_0$	$-K_3{\beta }_1$	$-K_3{\beta }_2$	$-K_3{\beta }_3$	…..	$-K_3{\beta }_{k-2}$
RHS	-1	0	0	0	0	0	0	0

The above set of equations yields: ${\beta }_0=-1,{\beta }_1-{\beta }_0=0,\&,{\beta }_k-{\beta }_{k-1}-K_3{\beta }_{k-2}=0,\forall k\ge 2$, which is an LDE for the βs, as.

$$\left[E^2-E-K_3\right]{\beta }_k\equiv 0,\forall k\ge 0,\text{with}\mathrm{ICs}:\left\{{\beta }_0=-1={\beta }_1\right\}.$$E47

We can note that the LHS Operator is the same as for the system LDE, and hence has the same characteristics and form of solution.

We first illustrate the computation for the stable case $K_3=-1/2$, for which the solution form has already been obtained as: $\left(A\mathrm{Cos}\left(\mathit{k\pi }/4\right)+B\mathrm{Sin}\left(\mathit{k\pi }/4\right)\right){\left(1/\sqrt{2}\right)}^k$, and hence plugging in the ICs, we have the solution for βs as:

$${\beta }_k\equiv \left(-\mathrm{Cos}\left(\mathit{k\pi }/4\right)-\mathrm{Sin}\left(\mathit{k\pi }/4\right)\right){\left(1/\sqrt{2}\right)}^k,k\ge 0.$$

To obtain the limiting inventory variance, we firstly note that:

$${\beta }_k^2={\left(-\mathrm{Cos}\left(\mathit{k\pi }/4\right)-\mathrm{Sin}\left(\mathit{k\pi }/4\right)\right)}^2{\left(1/2\right)}^k=\left(1-2\mathrm{Sin}\mathit{k\pi }/2\right){\left(1/2\right)}^k\le 3{\left(1/2\right)}^k.$$E48

And hence,

$${\mathit{Lim}}_{k\to \infty }{\sum }_{l=0}^k{\beta }_l^2\le 1+1+3{\left(1/2\right)}^2\left\{1+\left(1/2\right)+{\left(1/2\right)}^2+{\left(1/2\right)}^3+\dots \dots \right\}=2+3/2=3.5.$$E49

Hence, we have: ${\mathit{Lim}}_{k\to \infty }var\left(x_3\left(k\right)\right)\le 3.5{\sigma }^2$, showing that the inventory variance is bounded and finite for this case.

We next illustrate the computation for the marginally stable case $K_3=-1$, for which the solution form has already been obtained as: $A\mathrm{Cos}\left(\mathit{k\pi }/3\right)+B\mathrm{Sin}\left(\mathit{k\pi }/3\right)$, and hence plugging in the ICs, we have the solution for the βs as:

$${\beta }_k\equiv -\mathrm{Cos}\left(\mathit{k\pi }/3\right)-\left(1/\sqrt{3}\right)\mathrm{Sin}\left(\mathit{k\pi }/3\right)\equiv \left(2/\sqrt{3}\right)\mathrm{Cos}\left(\mathit{k\pi }/3-\pi /6\right),k\ge 0$$E50

And we can see from the above that the βs oscillate in value, from −1 to 0 to +1 infinitely often, i.e., the sequence {…0, −1, −1, 0, +1,+1, 0…..} repeats infinitely often. And hence the series ${\sum }_{k=0}^{\infty }{\beta }_k^2$ diverges to infinity.

Hence in this case the limiting inventory variance is not bounded and not finite.

It can similarly be shown that the limiting inventory variance is not bounded for unstable solutions also.

Thus, the limiting inventory variance will be bounded only for stable cases.

5.1.3 The complete solution

We can also note from the above discussion and work up that we have also obtained the stochastic part of the solution, as: $x_3{\left(k\right)}^{\mathit{\text{stoc}}}\equiv {\sum }_{l=0}^{k-1}{\beta }_l\epsilon \left(k-l\right)$, where the βs are as has been obtained above. Thus, the complete solution can be written as:

$$x_3\left(k\right)\equiv x_3{\left(k\right)}^{\det }+x_3{\left(k\right)}^{\mathit{\text{stoc}}}\equiv x_3{\left(k\right)}^{\det }+\sum _{l=0}^{k-1}{\beta }_l\epsilon \left(k-l\right)$$E51

where $x_3{\left(k\right)}^{\det }$ is the mean response derived from the solution of the deterministic part of the LDE, and the βs in the stochastic part are as obtained above by solution of the stochastic part of the system equation, the SDE.

The above representation proves useful for simulation purposes.

We next take up the P(ID) control. We discuss only the deterministic system LDE hereafter, since the stochastic part of the solution and the limiting inventory can be obtained by methods similar to that discussed above in all cases to follow.

5.2 P(ID) control under zero lag

In this type of control an additional demand-triggered component is also added to the control thereby making it more proactive. The control initiates corrective replenishment action no sooner than a demand deviation is observed. It does not wait for an inventory deviation to take place before initiating replenishment action though it does have an inventory-triggered component also.

The replenishment control flow is given by:

$$q_3\left(k+1\right)\equiv K_3{x}_3\left(k-1\right)+K_0^3{r}_3\left(k-1\right)$$E52

where the first term is the inventory-triggered component and the second the demand-triggered component. Substituting for the control flow into the system equation yields:

$$x_3\left(k+1\right)\equiv x_3\left(k\right)+K_3{x}_3\left(k-1\right)+K_0^3{r}_3\left(k-1\right)-r_3\left(k+1\right)$$E53

which can be written as

$$x_3\left(k+1\right)-x_3\left(k\right)-K_3{x}_3\left(k-1\right)\equiv K_0^3{r}_3\left(k-1\right)-r_3\left(k+1\right)$$E54

We can note that the addition of the demand-triggered component has left the LHS of the LDE unaltered. Thus, the LHS Operator of the LDE remains the same and is unaffected by addition of demand-triggered components to the control. The system eqn. can hence be written in Operator form as:

$$\left[E^2-E-K_3\right]x_3\left(k\right)\equiv K_0^3{r}_3\left(k\right)-r_3\left(k+2\right)\text{valid in}k\ge 0$$E55

Since $r_3\left(k\right)\equiv 1,\forall k\ge 1$, the system LDE can be written as:

$$\left[E^2-E-K_3\right]x_3\left(k\right)\equiv K_0^3-1\text{valid in}k\ge 1,$$E56

$$\text{with the}\mathrm{ICs}\mathrm{now}\mathrm{as}:\left\{x_3\left(1\right)=-1,x_3\left(2\right)=-2\right\}.$$E57

The ICs for the LDE have been obtained from the system Eq. (55) above using the standard system ICs; $\left\{\left(x_3\left(k\right),r_3\left(k\right)\right)\equiv \left(0,0\right),\forall k\le 0\right\}$ applied to the Eq. (55) above.

Since the LHS Operator is the same as for the earlier P(I) control, the stability analysis remains the same as earlier, as also the roots of the LHS Operator for various values of the inventory-trigger parameter discussed earlier, i.e.,$\left\{K_3=-1/4,-1/2,-1\right\}$.

The solutions are the same as given earlier in Eqs. (25) and (26).

And substituting the solution back into the O-NHE yields the value of the extra constant D as

$$D=\frac{K_0^3-1}{K_3},$$E58

which hence yields the solutions for the two cases as:

$$x_3\left(k\right)\equiv \left(K_0^3-1\right)/K_3+C_1{\lambda }_1^k+C_2{\lambda }_2^k\text{for the case of distinct roots}$$E59

$$x_3\left(k\right)\equiv \left(K_0^3-1\right)/K_3+\left(C_0+C_{1}k\right){\lambda }^k\text{for the case of repeated roots}(K_3=-1/4)$$E60

where the offset term is now given by $\left(K_0^3-1\right)/K_3$ for both cases.

The important point to note in the above solution is that the offset can be made zero by choice of the demand-trigger parameter as $K_0^3=1$.

And hence we can observe the enhanced response of the P(ID) control over the earlier P(I) control, in that the Offset can now be controlled by us by choice of $K_0^3$.

In fact, we can also achieve a (+)ve value of the offset by choosing $K_0^3\ge 1$.

We can now obtain the full solutions for the three cases above, using the LDE ICs $\left\{x_3\left(1\right)=-1,x_3\left(2\right)=-2\right\}$. We take $K_0^3=1$ to be able to obtain zero offset. Hence, we have:

$$\text{For}{K}_3=-1/4\left(\text{the repeated roots case}\right):x_3\left(k\right)\equiv \left(4-6k\right){\left(1/2\right)}^k,k\ge 1$$E61

$$\text{For}{K}_3=-1/2:x_3\left(k\right)\equiv \left(2\sqrt{5}\right)\mathrm{Cos}\left(\mathit{k\pi }/4-\varphi \right){\left(1/\sqrt{2}\right)}^k,tan\varphi =2,k\ge 1$$E62

$$\text{For}{K}_3=-1\left(\text{the marginal stability case}\right):x_3\left(k\right)\equiv 2\mathrm{Cos}\left(\left(k+1\right)\pi /3\right),k\ge 1$$E63

The solution curves are plotted in Figure 3, from which we can see that the response in all cases has zero offset.

Figure 3.

P(ID) control.

We can similarly extend the modeling and analysis to PI(I), PID(I), and MA(ID) controls.

We could also have different types of input disturbances as indicated earlier in Section 3.1. Additionally, the third dimension of our analysis could be to have non-zero lags, i.e., lags of one, two periods, and so on. Cases with non-zero and higher lags will result in higher order system LDEs, and the LHS Operator would be of a higher order.

Further details of the above can be obtained in [1, 2, 3, 4, 5, 6, 7, 8, 9].

6. Controls for multi-stage supply chains

We look at the serial supply chain system as given in Figure 1.

We can see from Figure 1 that the immediately succeeding downstream stage in a supply chain will provide the “demand perturbation” for the immediately preceding stage. Thus, the demand perturbation at the warehouse at the downstream end will successively be felt up the chain. And the single-stage analysis described above can be used in turn for each stage of the chain.

For non-serial supply chains, the arguments are similar and single-stage analysis can be used as described above.

Details of some of these analyses can be found in [1, 2, 3, 4, 5, 6, 7, 8, 9].

7. Conclusion

This chapter has presented the application of control concepts to the control of supply chains. The state variables have been taken to be the inventory levels, while the control variables are the replenishment flows into the various stages of the system. The conventional P, PI, PID controls have been discussed, as also some newer forms of control which are especially applicable to supply chains and warehouses. The performance of P(I) and P(ID) controls have been derived in detail, and their performance analyzed.

A significant feature of this chapter is that the conventional block diagrams and transfer functions of conventional control theory have not been used. Rather direct Operator Methods have been used to good advantage to solve the system equations.