Models and Algorithms of Inventory Control in Case of Uncertainty

Proceedings of the 13th IFAC Symposium on Information Control Problems in Manufacturing Moscow, Russia, June 3-5, 2009

Models and Algorithms of Inventory Control in Case of Uncertainty Alexander Mandel V.A. Trapeznikov Institute of Control Sciences of the RAS, 65, Profsoyuznaya str., 117997, Moscow Russia (Tel. 495-334-89-69, e-mail: [email protected]) Abstract: The paper discusses a set of models and algorithms for inventory control with uncertainty following the methodology of adaptive control theory and the theory of expert-statistical data processing. ∞

1. INTRODUCTION

Cn* ( x) = min{A × 1(u) + cu + g ( x, u) + α ∫ Cn*−1( x + u − ξ )dF (ξ )} (1)

Since the mid 1980s in the Institute of Control Sciences RAS (ICS) research has been under way (Lototsky, Mandel), 1987; Lototsky, Mandel, 1991; Belyakov, Lapin, Mandel, 2005; Mandel, Semenov, 2008) aimed at developing methods of inventory control in case of uncertainty of both, baseline and current, data describing the status of supply systems using adaptive and robust algorithms. In the mid 1990s in the Institute of Control Sciences RAS (ICS) a tool for integration of heterogeneous information within the same control system referred to as expertstatistical methods (ESM) for data processing was proposed (Borzenko, Lototsky, Mandel, 1990; Mandel’, 1996, 1997). Main applications of the new tool included different types of social and economic systems and inventory control systems, in particular. The authors review briefly inventory control methods in case of uncertainty and for different assumptions of the features of the stochastic demand formation process. The models and methods discussed are based on general methods of the adaptation theory (Tsypkin, 1968; Tsypkin, 1995) and expertstatistical approach in the form of the analog method (Belyakov, Mandel, 2002a; Belyakov, Mandel, 2002b; Mandel, 2004; Belyakov, Mandel, Semenov, 2008).

u ≥0

0

where C n* ( x ) is an optimal value of total average expenditures (accounting for discount) of the n-step process which starts at an assumed stock level x (i.e. the stock level accounting for the back-ordered demand and “stock in transit”, see (Hadley, Whitin, 1969)) and is implemented using an optimal inventory control strategy, 1 (.) is function of the unit step (Heaviside function), A is fixed expenditure for one shipment, c is the purchasing price for a unit of sale, u is value of replenishment order (decision variable) for moment of time n steps before the planning horizon, ξ is one step demand value, F(ξ) is unknown demand distribution function for ξ, dF(ξ) is Stieltjes differential, α is the discount coefficient, 0 ≤ α ≤ 1, while g ( x, u ) is an average value of one-step expenditures x +u

∞

0

x +u

g ( x, u ) = h ∫ ( x + u − z )dF ( z ) + d

(2)

where h is the unit cost coefficient for stock storage, d is the unit shortage cost, while u ≥ 0 is an order size for stock replenishment. 2.2. Adaptive algorithms for evaluation of “myopic” strategy parameters

2. INVENTORY CONTROL SYSTEMS IN THE CONTEXT OF STATIONARY DEMAND 2.1. Background In multi-step problems within the inventory control theory one of the most used criteria for selection of optimal control strategies is that of the minimum of total average expenditures (Hadley, Whitin, 1969). Most adequate to this separable criterion is a mathematical tool designed for selection and analysis of features of the optimal strategies which is the dynamic programming technique. The corresponding equations of the discrete dynamic programming when the demand is described in probabilistic terms can be presented as follows (Barladyan et al., 2006):

One of the most frequently discussed special cases of criterion (1) is a so-called “shortsighted”, or “myopic’, criterion of optimality for α = 0, when all current average expenditures are accounted for only at a given step of decision making. In this case the adaptive algorithms for recalculation of estimated parameters of the optimal twolevel inventory control system will be as follows (Barladyan et al,, 2006) ) Rˆ n( myop = Rˆ n( myop) − γ n [1( x n − z n ) − (d − с) /(h + d )] +1

223

(3)

) rˆn(+myop = 1

= rˆn( myop) − γ n′ [(d − h) /(d + h)rˆn( myop) − rˆn 1(rˆn( myop) − z n ) + + A /( h + d ) − z 1( Rˆ ( myop ) − z )1( z − rˆ ( myop ) )] , (4) n

978-3-902661-43-2/09/$20.00 © 2009 IFAC

∫ ( z − x − u)dF ( z ) ,

n

n

n

n

10.3182/20090603-3-RU-2001.0269

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

Dn =

where z n is the demand size at the n-th step, while coefficients γ n and γ n′ meet the conditions ∞

∞

∞

∞

n =1

n =1

n =1

n =1

∑ γ n = ∞, ∑ γ n2 < ∞ and ∑ γ ' n = ∞, ∑ γ ' n2 < ∞ .

(

)

rˆn +1 = rˆn + γ n'' [

3.1. Background

(6)

)

⎧ ⎫ (7) ⎞ 2 ⎛c ˆ ⎪ A + (c + h) Rn − crˆn zn − ⎜ + h ⎟ z n + hDn + ⎪ 2 1 ⎪ ⎪ ⎠ ⎝ + ⎨ ⎬], ˆ 2 − rˆ 2 ⎪ ( Rˆ n − rˆn ) 2 ⎪ R n n ˆ ˆ ˆ ˆ ⎪⎩+ (h + d )η2 ( Rn , rn ; ξ n ) + d ( Rn − rn ) zn − d ⎪⎭ 2

where coefficients γ n and γ n′ still satisfy: (5), function η 2 ( R, r; ξ ) is described as

η 2 ( R, r ; ξ ) =

Rч

R

r0

r

∫ ∫ 1( z − ξ )dzdx = ∫ max{x − ξ ,0}dx

⎧ R2 − r 2 − ξ ( R − r ), if ξ ≤ r , ⎪ ⎪ 2 2 2 ⎪⎪ R − ξ =⎨ − ξ ( R − ξ ), if r ≤ ξ < R, ⎪ 2 0, if R ≤ ξ , ⎪ ⎪ ⎪⎩

,

(8)

while the recurrent estimates of the average demand size z n , second moment z 2 n and dispersion D n are described by the following formulas zn =

1 n −1 z n −1 + ξ n , n n

z2 n =

1 n −1 2 z n −1 + ξ n2 , n n

(9) (10)

Considerable part of the sales revenue of many manufacturers comes from selling commodities of higher demand or “hits”. Obviously, formation of a line of such goods is directly related to the expert solution of active and passive marketing problems. In this case passive marketing implies the analytical study of the demand on the market for the subsequent generation of proposed new types of “superfashionable” products, whereas active marketing is interpreted here as the organization and running of productive advertising campaigns. After a portfolio of the “hits” is formed and a sales strategy is chosen including the advertising campaign timeline and sales schedule, the next high priority becomes the volume and schedule of production (if the manufacturing company has its own distribution network) or purchase (in case of a procurement company) of the respective goods.

1 ( A + Rˆ n (c + h) − crˆn + (c − h) z n ) ˆ Rn − rˆn

(

(11)

3. MODELS OF INVENTORY CONTROL IN CASE OF SPECULATIVE DEMAND

Adaptive algorithms for recalculation of estimated parameters of optimal two-level inventory control strategies for a restricted stock-out case are discussed in (Mandel, Belyakov, Semenov, 2008; Mandel, Semenov, 2008)

c ⎧ ⎫ A + (c + h) Rˆ n − crˆn z n − z 2 n + ⎪ ⎪⎪ 2 1 ⎪ ], ⎨ 2 2 ˆ − rˆ ⎬ R ( Rˆ n − rˆn ) 2 ⎪ n n ⎪ + (h + d )η 2 ( Rˆ n , rˆn ; ξ n ) − d 2 ⎩⎪ ⎭⎪

n −1 1 ≅ Dn −1 + (ξ n − zn ) 2 . n n

(5)

2.3. Adaptive algorithms for evaluation of optimal strategy parameters in a general case of the criterion (1)

c+h Rˆ n +1 = Rˆ n − γ n' [ zn − ˆ R n − rˆn

1 n 1 n (ξ k − zn ) 2 ≅ ∑ (ξ k − zn −1 ) 2 ≅ ∑ n k =1 n k =1

To calculate the above production or procurement schedules one needs an adequate model of demand for the types of goods newly launched on the market. All ultrafashionable goods have one feature in common, i.e. the lion’s share of their sales is limited in time which normally is several weeks (usually from 7-8 to 14-16). After the consumer boom is over the goods are either replaced on the market by the new “hits” or slip to the category of “classics”, i.e. goods which enjoy stable but dozen-fold lower demand. It is noteworthy that the revenue received from the sales of ultrafashionable goods frequently accounts for half of the total revenue of trading and manufacturing companies. Information concerning the demand for these goods during the total period of their active sales is based on assumptions only, hence any mistakes in the assumptions can result in gross losses – from the lost consumer demand, if underestimated, to formation of the absurdly large stock of goods the demand for which had been overestimated. 3.2. Expert-Statistical Procedure of Demand Prediction by the Analog Method

The method rests on the assumption that when forecasting short time series the experts in the absence of essential a priori information on the demand for the newly launched goods do their forecasting with whatever information they might have on the objects or processes whose background is well known to them. It is further assumed that the number of such objects or processes is large enough and that the space of features of the objects which actually is the contents of the experience that the experts have yields itself to clear-cut or, which is frequently the case, rather vague classification

224

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

(Bauman, 1988; Bauman, Dorofeuyk, 1982, Bauman, Dorofeuyk, 1999). Formally this implies that for a given functional of the metric type ⎛ M Ai D = ∑ p Ai χ⎜⎜ A i i ⎝ p

⎞ ⎟, ⎟ ⎠

(12)

where χ is convex function, Ai are classes of points (places), p Ai are a priori probabilities of classes Ai , while M Ai are

first non-normalized moments of classes Ai , we can build algorithms (recurrent, for instance) of search for its extreme points. The attribution to a class (if it is well-defined), is established by a characteristic function, taking the value 1 on the object belonging to a fixed class (otherwise, it is 0), while if the class is not distinct (fuzzy), then the function of attribution hi ( x) : 0 ≤ hi ( x) ≤ 1 with ∑ (hi ( x)) 2 = 1 . i

As shown in (Belyakov, Mandel, 2002a; Mandel, 2004; and Mandel, Belyakov, Semenov, 2008) interaction of the expert and the expert-statistical prediction system (ESPS) provides the set Z of analogs of the predictable object (PO) under consideration. For this set, the ESPS database contains information about “complete”, that is, represented by much longer time series, realizations of operation of the analog objects. This information is represented by the set {x k ( n ), k ∈ Z , n = 0, 1, 2,... , N 1 } where N 1 >> N . Additionally, sets of values of the similarity, {lk , k ∈ Z } , and scale, {s k , k ∈ Z } , coefficients are given. To predict the values of the PO time series at time n, n > N, the following formula can now be used: yˆ ( n ) = L−1 ∑ α k l k s k x ( n ) ,

(13)

k∈Z

where L =

∑ lk

. For N > 0, the values of the coefficients

k∈Z

Consider a problem of inventory control in a warehouse for a finite time interval with statistically unspecified demand characteristics in the presence of a priori information generated through the use of the expert-statistical approach. The problem is to choose an optimal control strategy, i.e. the sequence of requests for stock replenishment which minimizes total average expenditures. Let us use the adaptive dynamic programming technique for the solution of the problem. In the system under consideration the control is done in the finite time interval T . Following regular time intervals of duration τ (it is assumed that T = Nτ , where N is a positive integer) the system state (stock level) is monitored. We will count time from the end of the period and denote as t n the time moment n steps prior to the end of interval T (so called “backward” time). In every sub-interval [t n ; t n −1 ) a decision is made on initiating (or not initiating) a request for replenishment of stock u n . We assume that all requests are fully met by the supplier. Let F (x, θ ) be the function of the demand distribution in the interval [t n ; t n −1 ) , which is a function of the unknown vector of parameters θ ; G0 (θ ) is a priori distribution of vector θ (at time t N ); G k (θ ) is a posteriori distribution of vector θ (at time t N − k , i.e. at the k-th step of the “forward” time, n + k = N); C n* ( x, G k ) is as low as minimum value of total T provided that at moment t n the stock was x; whereas a posteriori distribution of the vector of parameters θ at the kth step in the forward time (k + n = N) is G k (θ ) .

2

⎛ ⎞ min L ∑n =1 ⎜ ∑α k l k s k x( n ) ⎟ . {α , k∈Z } ⎝ k∈Z ⎠ N

3.4.Inventory control: adaptive approach

average expenditures from time t n to the end of sales period

α k , k ∈ Z , in (13) are established from:

−2

− the expert for some reasons is unwilling or finds it difficult to reveal the interval or point estimates of the future values of time series; − there is expert information about the PO that allows one to classify (identify) it with one or another similarity class; − there exists a representative set of statistical information about a substantial number of objects from the given knowledge domain.

(14)

k

If N = 0, that is, if there is no data sample on PO at all, then all α k , k ∈ Z , are assumed to be equal to 1. 3.3. Analog method: field of application and practical recommendations

It is recommended to use the expert-statistical prediction procedures based on the analog method if:

The problem is to determine such a sequence of requests for un , stock replenishment (the control strategy) n = N , N - 1, …, 1 , which minimizes total average expenditures for the whole period T, i.e. the strategy which determines the functional C N* ( x, G0 ) .

The following equations of adaptive dynamic programming hold

− there is no statistical information about the PO or prediction can be based only on the subjective information;

225

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

where ξ k is white noise with zero expectation and it is assumed that ξ k < bk . It is further assumed that the effective

x +u

Cn* ( x, Gk ) = min[( A + cu ) + + h ∫ [ ∫ ( x + u − ξ ) dF (ξ ,θ ) + u ≥0

Θ

0

∞

demand measurement is done with random errors according to

(15)

∫ (ξ − x − u )dF (ξ ,θ ) +

+d

ς k = zk + ηk ,

x +u ∞

+ α ∫ Cn*−1 ( x + u − ξ , Gk +1 ) f (ξ ,θ )dξ ]dGk (θ ),

(19)

where η k is also discrete white noise with zero mean.

0

where Θ is a set of possible values of parameter θ.

We use as the criterion either total average expenditures (1), or the mean square deviation minimum of the stock level R.

Calculation of a posterior distributions Gk (θ ) using a priori at the k-th step (forward time) distribution Gk −1 (θ ) is done as

The above inventory control problem can be solved by means of the Kalman filtration.

∫ f (x

(1)

,θ ) f ( x

( 2)

,θ )... f ( x

(k )

,θ ) dG0 (θ )

Θ

=

f (x

(k )

∫ f (x

4.2. Inventory control algorithms using the Kalman filtration

f ( x (1) ,θ ) f ( x ( 2) ,θ )... f ( x ( k ) ,θ )dG0 (θ )

dGk (θ ) =

,θ ) dGk −1 (θ )

(k )

,θ ) dGk −1 (θ )

Using any of the criteria discussed in the previous section the inventory control algorithm is given here as

(16)

,

u k +1 = (R − z k + zˆ k +1 )+ ,

Θ

where f ( x (i ) ,θ ) is the density of probabilistic distribution F ( x, θ ) in the observed at the i-th step of the forward

time point x (i ) , whereas G0 (θ ) is a priori distribution of parameter θ .

Hence, the best possible algorithm of search for an optimal inventory control strategy in the form of a totality of stock replenishment solutions {u n , n = 1,2,..., N } can be written as C0* ( x, GN ) ≡ 0 for ∀GN ,

(20)

where u k +1 is the order size at time k + 1, R is either the upper level in the optimal two-level hierarchical inventory control strategy in terms of criterion (1), or the specified value of the stock level with respect to which the supply system is to be stabilized, and zˆk +1 в is the best, in terms of the least squares, forecast of the demand at the (k + 1)-th step. So, if the regression coefficients a k and bk in (18) are known, then the solution of the inventory control problem is reduced to the formation of a demand forecast using the following algorithm of the Kalman filtration

x +u

Cn* ( x, Gk ) = min[( A + cu) + + h ∫ [ ∫ ( x + u − ξ)dF(ξ, θ) + u≥0

Θ

zˆ k +1 = ak z k + bk + γ k (ς k − z k ),

0

∞

+ d ∫ (ξ − x − u)dF(ξ, θ) +

γk = (17).

Dk2+1

x +u ∞

+ β ∫ Cn*−1 ( x + u − ξ, Gk +1 ) f (ξ, θ)dξ]dGk (θ),

=

/

(

Dk2

Dk2

2

)

+ ση , 2

+ σ ξ − a k2 Dk4

/

(

Dk2

z 0 = 0, 2

)

+σξ ,

(21) (22)

D02

= 0.

(23)

5. EXAMPLE OF ANALOG METHOD APPLICATION

0

where n = 1,2,..., N and Gk +1 is defined by formula (16).

4. INVENTORY CONTROL MODELS IN THE CONTEXT OF THE AUTOREGRESSIVE DEMAND CHARACTERIZATION 4.1. Background Let the demand formation process zk, k – forward time, be described by the first-order autoregressive process (Lototsky, Mandel, 1992)

z k +1 = ak z k + bk + ξ k ,

Dk2

(18)

As described in 3.1, a period of the speculative demand for most of the ultrafashionable goods normally lasts for several weeks. Let us refer to the respective time interval T = Nτ (where τ is a time unit, in our case it is one week) as the sales period. In other words, a sales period is the time interval during which the demand essentially exceeding the average level is registered. After the sales period is over the demand drops dramatically, i.e. by dozen fold. The present author has a database with the sales history for several thousands of item types of a major Russian commercial chain. Many of the items fall in the category of the “hits”, i.e. goods enjoying the speculative demand. The analysis showed three most frequently encountered demand curves in the sales interval for the goods of the speculative demand: 1)

226

“one-humped” demand curves (Fig. 1);

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

2)

“two-humped” demand curves (Fig. 2);

3)

“disordered” demand curves (Fig. 3).

The demand curves referred to as “disordered”, present as the alteration of more than two peaks (maxima) and lows (minima) of the demand in the sales interval. Consequently, the first-, second-, and third-type of the demand were classified as three classes, and their data (within each class) have been normalized so that to reduce all curves to the same maximum amplitude and the same sales period. Specifically, each of the demand curves ξ i (t ) in the sales period [0, Ti] was replaced with curves ϑi (t ′) , obtained from ξ i (t ) through the substitution of the variables: t ′ = t / Ti , ϑ i (t ′) = ξ i (t ) max ξ i (t ) .

Fig. 2. Second type demand curve.

t∈[ 0,Ti ]

Subsequently, from each of the assemblies of the normalized curves referred to a specific class, “similar” curves had been selected (for instance, such curves for the class of “onehumped” whose time t ′ , of reaching max ξ i (t ) , coincided). t∈[ 0,Ti ]

Next, for the identified sub-assemblies the estimated pairwise correlation coefficients were calculated for nearest neighbor, time-wise, points (separated one from the other by one time step which, prior to normalization, was equal to τ). We will not present the kind of the respective relationships between the estimated correlation coefficients and time t ′ , just mention several remarkable facts which will be extremely relevant in this discussion: 1)

for curves of class 1 demand (“one-humped” curves) the pair-wise correlation coefficient almost the whole sales period was above 0,85;

2)

for curves of class 2 demand (“two-humped” curves) the pair-wise correlation coefficient almost the whole sales period was above 0,6;

3)

for curves of class 3 demand (“disordered” curves) the pair-wise correlation coefficient almost the whole sales period rarely exceeded the threshold of 0,1.

Hence, one can identify such sub-assemblies (synonyms: subtypes, sub-classes) of the 1st and 2nd type curves (especially in case of “one-humped” demand curves), for which in the normalized version almost all the deterministic dependencies of demand on time are realized (almost the entire randomness is “concentrated” in Ti only) and max ξ i (t ) . t∈[ 0,Ti ]

Fig. 1. First type demand curve.

Fig. 3. Third type demand curve. To the contrary, it is much more natural to describe the 3rd type demand curves using the model of a sequence of mutually independent random values which might well have one common probability distribution for all demand values (the model of stationary discrete in time “white noise”). 6. AKNOWLEGMENT This work is supported by Russian Foundation for Basic Research, project 09-07-00195-а. 7. CONCLUSIONS A set of models and algorithms of inventory control with uncertainty is discussed for the stationary, speculative and autoregressive description of a stochastic process of demand formation. The algorithms discussed use general techniques of the adaptive system theory, Kalman filtration and methods of the expert-statistical data processing. REFERENCES Barladyan, I.I., N.I. Borzenko, A.V. Lapin, A.S. Mandel and A.B. Tokmakova. (2006). Comparative Analysis of Myopic and Long-sighted Inventory Control Policies under Uncertainty. In: Izbrannye Trudy 2 Mezhdunarodnoi Konferentsii po Problemam Upravleniya, IPU RAN, Moscow (in Russian). Bauman, E.V. (1988). Methods of Fuzzy Classification (Variational Approach). Automation and Remote Control, v. 49, no. 12. Bauman, E.V. and A.A. Dorofeyuk (1982). Recurrent Algorithms of Automatic Classification. Automation and Remote Control, v. 43, no. 3. Bauman, E.V. and A.A. Dorofeyuk (1999). Classification Data Analysis. In: Izbrannye Trudy Mezhdunarodnoi

227

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

Konferentsii po Problemam Upravleniya (Int. Conf. on Control Sciences, Selected Papers), v. 1. SINTEG, Moscow (in Russian). Bellman, R. and R. Kalaba. (1960). Dynamic programming and adaptive processes: mathematical foundation. IRE Trans. On AC, v. AC-5, no. 1. Belyakov, A.G. and A.S. Mandel’ (2002a). Prediction of Time Series on the Basis of Analog Method (Elements of the Theory of Expert-Statistical Systems). Preprint of Institute of Control Sciences, Moscow (in Russian). Belyakov, A.G. and A.S. Mandel’ (2002b). Analysis of Reliability of Conclusions Generated by Means of ExpertStatistical Systems. Preprint of Institute of Control Sciences, Moscow (in Russian). Belyakov, A.G., A.V. Lapin and A.S. Mandel' (2005). Inventory Control for Goods in Speculative Demand. Problemy Upravleniya, v. 6. (in Russian). Borzenko N.I., V.A. Lototsky and A.S. Mandel' (1990). Expert Statistical Systems for Demand Prediction and Inventory Control. In: Evaluation of Adaptive Control Strategies in Industrial Applications – IFAC Workshop Series, no. 7. Pergamon Press, Oxford. Hadley, G. and Т.M. Whitin (1969). Analysis of Inventory Systems. Prentice Hall, Inc. Englewood Cliffs, New Jersey. Kalaba, R. (1966). Mathematical Aspects of Adaptive Regulation. In: Modeli biologicheskih system. Mir publishers, Moscow (in Russian). Lototsky, V.A. and A.S. Mandel' (1987). Adaptive Inventory Control. In: Preprints of X IFAC World Congress (Munich), v. 10, subv. 6. DI/VDE, Dusseldorf. Lototsky, V.A. and A.S. Mandel' (1991). Models and Inventory Control Techniques. Nauka publishers, Moscow (in Russian). Mandel’, A.S. (1996). Expert-Statistical Systems in Control and Information Processing. Part I. Pribory&Sistemy Upravleniya, No. 12 (in Russian). Mandel’, A.S. (1997). Expert-Statistical Systems in Control and Information Processing. Part II. Pribory&Sistemy Upravleniya, no. 1 (in Russian). Mandel’, A.S. (2004). Method of Analogs in Prediction of Short Time Series: An Expert-statistical Approach. Automation and Remote Control, v. 65, no. 4. Mandel’, A.S., A.G. Belyakov and D.A. Semenov (2008). Expert-Statistical Processing of Data and the Method of Analogs in Solution of Applied Problems in Control Theory. In: Preprints of 17th IFAC World Congress. July 6-11, 2008, Seoul, Korea. Mandel’, A.S. and D.A. Semenov (2008). Adaptive Algorithms for Estimation of Optimal Inventory Control Strategy Parameters under Restricted Stock-Out. Automation and Remote Control, v. 69, no. 6. Tsypkin, Ya.Z. (1968). Adaptation and Learning in Automatic Systems. Nauka publishers, Moscow (in Russian). Tsypkin, Ya.Z. (1995). Identification Information Theory. Nauka publishers, Moscow (in Russian).

228