**Optimal Control of Integrated Production – Forecasting System**

R. Hedjar, L. Tadj and C. Abid

Additional information is available at the end of the chapter

**1. Introduction**

[71] Wognum, P.M., Faber, C., 2003, Infrastructures for collaboration in virtual organiza‐ tions, Int. Journal of Computer Applications in Technology, Geneva, Vol. 18, Issue

[72] Zweger, A., Tolle, M., Vesterager, J., 2002, VERAM: Virtual Enterprise Reference Ar‐ chitecture and Methodology, Global Engineering and Manufacturing in Enterprise Networks GLOBEMEN, eds. I. Karvonen et. al., Julkaisija Utgivare Publisher, Fin‐

1-4.

116 Decision Support Systems

land.

The production planning problem has received much attention, and many sophisticated models and procedures have been developed to deal with this problem. Many other compo‐ nents of production systems have also been taken into account by researchers in so called integrated systems, in order to achieve a more effective control over the system.

In this work, optimal control theory is used to derive the optimal production rate in a manu‐ facturing system presenting the following features: the demand rate during a certain period depends on the demand rate of the previous period (dependent demand), the demand rate depends on the inventory level, items in inventory are subject to deterioration, and the firm can adopt a periodic or a continuous review policy. Also, we are using the fact that the cur‐ rent demand is related to the previous demand in order to integrate the forecasting compo‐ nent into the production planning problem. The forecast of future demand for the products being produced is needed to plan future activities. Forecasting information is an important input in several areas of manufacturing activity. This problem has been considered in the literature. The proposed approach is different from that of other authors which is mainly based on time-series. In [1], the authors deal with the interaction between forecasting and stock control in the case of non-stationary demand. In [2], the authors assume a distribution for the unknown demand, estimate its parameters and replace the unknown demand pa‐ rameters by these estimates in the theoretically correct model. In [3], the authors propose an approach to evaluate the impact of interaction between demand forecasting method and stock control policy on the inventory system performances. In [4], the authors present a sup‐ ply chain management framework based on model predictive control (MPC) and time series forecasting. In [5], the authors consider a data-driven forecasting technique with integrated inventory control for seasonal data and compare it to the traditional Holt-Winters algorithm

© 2012 Hedjar et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Hedjar et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

for random demand with a seasonal trend. In [6], the authors assess the empirical stock con‐ trol performance of intermittent demand estimation procedures. In [7], the authors study two modifications of the normal distribution, both taking non-negative values only.

The goal of this chapter is to study the same problem in periodic and continuous review pol‐ icy context, knowing that the inventory can be reviewed continuously or periodically. In a continuous-review model, the inventory is monitored continually and production/order can be started at any time. In contrast, in periodic-review models, there is a fixed time when the

Optimal Control of Integrated Production – Forecasting System 119

We assume that the firm has set an inventory goal level, a demand goal rate and a produc‐ tion goal rate, to build the objective function of our model. The inventory goal level is a safe‐ ty stock that the company wants to keep on hand. The demand goal rate is the amount that the company wishes to sell per unit of time. The production goal rate is the most efficient rate desired by the firm. The objective is to determine the optimal production rate that will keep the inventory level, the demand rate, and the production rate as close as possible to the inventory goal level, the demand goal rate, and the production goal rate, respectively.

Therefore, we deal with a dynamic problem and the solution sought, the optimal production rate, is a function of time. The problem is then represented as an optimal control problem with two state variables, the inventory level and the demand rate, and one control variable,

The rest of this chapter is organized as follows. In section 2, the notation used is introduced and the dynamics of the system are described for both periodic and continuous review sys‐ tems. In section 3, the optimal solution is computed for each case. Simulations are conducted in section 4 to verify the results obtained theoretically in section 3. Section 5 summarizes the

Consider a manufacturing firm producing units of an item over some time interval 0, *T* , where *T* >0. Let *I*(*t*), *D*(*t*), and *P*(*t*) represent the inventory level, the demand rate, and the

goals at time *t*. Also, let *h* , *K*, *r* represent the penalties for each variable to deviate from its

^(*t*) 2}*dt* <sup>+</sup> <sup>1</sup>

(*<sup>T</sup>* ) <sup>2</sup> <sup>+</sup> *KT <sup>D</sup>*(*<sup>T</sup>* ) - *<sup>D</sup>*

^ (*t*), and *<sup>P</sup>*

^

<sup>2</sup> {*h <sup>T</sup> I*(*T* ) - *I*

^ (*<sup>T</sup>* ) <sup>2</sup>

^

^

(*t*) represent the corresponding

} gives the salvage value of the

(*t*), the objective function is

^ (*<sup>T</sup>* ) <sup>2</sup>

} (1)

(*<sup>T</sup>* ) <sup>2</sup> <sup>+</sup> *KT <sup>D</sup>*(*<sup>T</sup>* ) - *<sup>D</sup>*

^ (*t*), *D*

inventory is reviewed and a decision is made whether to produce/order or not.

the rate of manufacturing.

**2. Model formulation**

In (1), the expression <sup>1</sup>

expressed as

min *P*(*t*) *<sup>J</sup>* <sup>=</sup> <sup>1</sup> 2 *∫* 0 *T* {*h I*(*t*) - *I* ^

chapter and outlines some future research directions.

**2.1. Continuous review integrated production model**

goal. Then, the objective function *J*to minimize is given by

<sup>2</sup> {*h <sup>T</sup> I*(*T* ) - *I*

ending state. Using the shift operator defined by Δ *f* (*t*)= *f* (*t*) - *f*

^ (*t*) <sup>2</sup> <sup>+</sup> *<sup>r</sup> <sup>P</sup>*(*t*) - *<sup>P</sup>*

^

production rate at time *t*, respectively. Let *I*

(*t*) <sup>2</sup> + *K D*(*t*) - *D*

Many researchers have investigated the situation where the demand rate is dependent on the level of the on-hand inventory. In [8], the authors consider an inventory model under inflation for deteriorating items with stock-dependent consumption rate and partial back‐ logging shortages. In [9], the authors examine an inventory model for deteriorating items under stock-dependent demand and two-level trade credit. The reference [10] deals with a supply chain model for deteriorating items with stock-dependent consumption rate and shortages under inflation and permissible delay in payment. In [11], the authors deal with the optimal replenishment policies for non-instantaneous deteriorating items with stock-de‐ pendent demand. In [12], the authors investigate an inventory model with stock–dependent demand rate and dual storage facility. In [13], the authors develop a two warehouse inven‐ tory model for single vendor multiple retailers with price and stock dependent demand. In [14], the authors asses an integrated vendor-buyer model with stock-dependent demand. In [15], the authors study an EOQ model for perishable items with stock and price dependent demand rate. In [16], the authors develop the optimal replenishment policy for perishable items with stock-dependent selling rate and capacity constraint. In [17], the authors consider an inventory model for Weibull deteriorating items with price dependent demand and timevarying holding cost. In [18], the authors study fuzzy EOQ models for deteriorating items with stock dependent demand and nonlinear holding costs. In [19], the authors approach an extended two-warehouse inventory model for a deteriorating product where the demand rate has been assumed to be a function of the on-hand inventory. In [20], the authors investi‐ gate a channel who sells a perishable item that is subject to effects of continuous decay and fixed shelf lifetime, facing a price and stock-level dependent demand rate. In [21], the au‐ thors develop a mathematical model to formulate optimal ordering policies for retailer when demand is practically constant and partially dependent on the stock, and the supplier offers progressive credit periods to settle the account. The literature on stock-dependent de‐ mand rate is abundant. We have reported some of it here but only a comprehensive survey can summarize and classify it efficiently.

In [22], the authors review the most recent literature on deteriorating inventory models, clas‐ sifying them on the basis of shelf-life characteristic and demand variations. In [23], the au‐ thors introduce an order-level inventory model for a deteriorating item, taking the demand to be dependent on the sale price of the item to determine its optimal selling price and net profit. In [18], the authors formulate an inventory model with imprecise inventory costs for deteriorating items under inflation. Shortages are allowed and the demand rate is taken as a ramp type function of time as well. In [10], the authors model the retailer's cost minimiza‐ tion retail strategy when he confronts with the supplier trade promotion offer of a credit policy under inflationary conditions and inflation-induced demand. In [24], the authors de‐ velop two deterministic economic production quantity (EPQ) models for Weibull-distribut‐ ed deteriorating items with demand rate as a ramp type function of time.

The goal of this chapter is to study the same problem in periodic and continuous review pol‐ icy context, knowing that the inventory can be reviewed continuously or periodically. In a continuous-review model, the inventory is monitored continually and production/order can be started at any time. In contrast, in periodic-review models, there is a fixed time when the inventory is reviewed and a decision is made whether to produce/order or not.

We assume that the firm has set an inventory goal level, a demand goal rate and a produc‐ tion goal rate, to build the objective function of our model. The inventory goal level is a safe‐ ty stock that the company wants to keep on hand. The demand goal rate is the amount that the company wishes to sell per unit of time. The production goal rate is the most efficient rate desired by the firm. The objective is to determine the optimal production rate that will keep the inventory level, the demand rate, and the production rate as close as possible to the inventory goal level, the demand goal rate, and the production goal rate, respectively.

Therefore, we deal with a dynamic problem and the solution sought, the optimal production rate, is a function of time. The problem is then represented as an optimal control problem with two state variables, the inventory level and the demand rate, and one control variable, the rate of manufacturing.

The rest of this chapter is organized as follows. In section 2, the notation used is introduced and the dynamics of the system are described for both periodic and continuous review sys‐ tems. In section 3, the optimal solution is computed for each case. Simulations are conducted in section 4 to verify the results obtained theoretically in section 3. Section 5 summarizes the chapter and outlines some future research directions.

### **2. Model formulation**

for random demand with a seasonal trend. In [6], the authors assess the empirical stock con‐ trol performance of intermittent demand estimation procedures. In [7], the authors study

Many researchers have investigated the situation where the demand rate is dependent on the level of the on-hand inventory. In [8], the authors consider an inventory model under inflation for deteriorating items with stock-dependent consumption rate and partial back‐ logging shortages. In [9], the authors examine an inventory model for deteriorating items under stock-dependent demand and two-level trade credit. The reference [10] deals with a supply chain model for deteriorating items with stock-dependent consumption rate and shortages under inflation and permissible delay in payment. In [11], the authors deal with the optimal replenishment policies for non-instantaneous deteriorating items with stock-de‐ pendent demand. In [12], the authors investigate an inventory model with stock–dependent demand rate and dual storage facility. In [13], the authors develop a two warehouse inven‐ tory model for single vendor multiple retailers with price and stock dependent demand. In [14], the authors asses an integrated vendor-buyer model with stock-dependent demand. In [15], the authors study an EOQ model for perishable items with stock and price dependent demand rate. In [16], the authors develop the optimal replenishment policy for perishable items with stock-dependent selling rate and capacity constraint. In [17], the authors consider an inventory model for Weibull deteriorating items with price dependent demand and timevarying holding cost. In [18], the authors study fuzzy EOQ models for deteriorating items with stock dependent demand and nonlinear holding costs. In [19], the authors approach an extended two-warehouse inventory model for a deteriorating product where the demand rate has been assumed to be a function of the on-hand inventory. In [20], the authors investi‐ gate a channel who sells a perishable item that is subject to effects of continuous decay and fixed shelf lifetime, facing a price and stock-level dependent demand rate. In [21], the au‐ thors develop a mathematical model to formulate optimal ordering policies for retailer when demand is practically constant and partially dependent on the stock, and the supplier offers progressive credit periods to settle the account. The literature on stock-dependent de‐ mand rate is abundant. We have reported some of it here but only a comprehensive survey

In [22], the authors review the most recent literature on deteriorating inventory models, clas‐ sifying them on the basis of shelf-life characteristic and demand variations. In [23], the au‐ thors introduce an order-level inventory model for a deteriorating item, taking the demand to be dependent on the sale price of the item to determine its optimal selling price and net profit. In [18], the authors formulate an inventory model with imprecise inventory costs for deteriorating items under inflation. Shortages are allowed and the demand rate is taken as a ramp type function of time as well. In [10], the authors model the retailer's cost minimiza‐ tion retail strategy when he confronts with the supplier trade promotion offer of a credit policy under inflationary conditions and inflation-induced demand. In [24], the authors de‐ velop two deterministic economic production quantity (EPQ) models for Weibull-distribut‐

ed deteriorating items with demand rate as a ramp type function of time.

two modifications of the normal distribution, both taking non-negative values only.

can summarize and classify it efficiently.

118 Decision Support Systems

#### **2.1. Continuous review integrated production model**

Consider a manufacturing firm producing units of an item over some time interval 0, *T* , where *T* >0. Let *I*(*t*), *D*(*t*), and *P*(*t*) represent the inventory level, the demand rate, and the production rate at time *t*, respectively. Let *I* ^ (*t*), *D* ^ (*t*), and *<sup>P</sup>* ^ (*t*) represent the corresponding goals at time *t*. Also, let *h* , *K*, *r* represent the penalties for each variable to deviate from its goal. Then, the objective function *J*to minimize is given by

$$\min\_{\mathbf{P}(t)} \text{I} = \frac{1}{2} \Big[ \big[ h\big[ \mathbf{I}(t) \cdot \big[ \mathbf{I}(t) \big] \big] + \mathbf{K}\big[ \mathbf{D}(t) \cdot \big[ \mathbf{D}(t) \big] \big] + r\big[ \mathbf{P}(t) \cdot \big[ \mathbf{I}(t) \big] \big] \big] \big] \text{I} \\ + \frac{1}{2} \Big[ h\big[ \mathbf{I}(T) \big] \big] \big[ \mathbf{I}(T) \big] \big] + \mathbf{K}\_{\Gamma} \big[ \mathbf{D}(T) \big] \big \big[ \mathbf{I}(T) \big] \big] \big] \tag{1}$$

In (1), the expression <sup>1</sup> <sup>2</sup> {*h <sup>T</sup> I*(*T* ) - *I* ^ (*<sup>T</sup>* ) <sup>2</sup> <sup>+</sup> *KT <sup>D</sup>*(*<sup>T</sup>* ) - *<sup>D</sup>* ^ (*<sup>T</sup>* ) <sup>2</sup> } gives the salvage value of the ending state. Using the shift operator defined by Δ *f* (*t*)= *f* (*t*) - *f* ^ (*t*), the objective function is expressed as

$$\min\_{\boldsymbol{P}(t)} \mathcal{J} = \frac{1}{2} \left[ \mathbb{I} h \,\Delta^2 I(t) + K \,\Delta^2 \mathcal{D}(t) + r \,\Delta^2 \mathcal{P}(t) \right] dt + \frac{1}{2} \left[ h\_T \,\Delta^2 I(T) + K\_T \,\Delta^2 \mathcal{D}(T) \right] \tag{2}$$

To use a matrix notation, which is more convenient, let *<sup>X</sup>* (*t*)= <sup>Δ</sup>*I*(*t*) <sup>Δ</sup>*D*(*t*) and let *X <sup>A</sup>* <sup>2</sup> = *X <sup>T</sup> AX* . Then, the objective function (2) can be further rewritten as

$$\min\_{P(t)} J = \frac{1}{2} \Big[ \Big\| X \left( t \right) \Big\|\_{H}^{2} + r \Delta^{2} P(t) \Big\| dt + \frac{1}{2} \Big\| X \left( T \right) \Big\|\_{H\_{T}}^{2} \tag{3}$$

where *<sup>H</sup>* <sup>=</sup> *<sup>h</sup>* <sup>0</sup> <sup>0</sup> *<sup>K</sup>* and *HT* <sup>=</sup> *h <sup>T</sup>* 0 0 *KT* .

Two state equations are used to describe the dynamics of our system. The variations of the inventory level and demand rate are governed by the following state equations

$$\frac{d}{dt}I(t) = P(t) \cdot D(t) \cdot \Theta I(t) \tag{4}$$

**2.2. Periodic review integrated production model**

demand rate *D*(*k*).

where *α* =1 - *θ*.

equation:

goal rate *P*

^

shift operator Δ to get:

*<sup>D</sup>*(*k*) , *<sup>Z</sup>*

^ (*k*)= *<sup>I</sup>* ^ (*k*) *D* ^ (*k*) , *Z*<sup>0</sup> <sup>=</sup>

*<sup>Z</sup>*(*k*)= *<sup>I</sup>*(*k*)

function to minimize is given by:

lowing difference equation:

In the periodic review model, the time interval 0, *T* is divided into *N* subintervals of equal length. During period *k*, the plant manufactures units of some product at the control‐ lable rate *P*(*k*), and the state of the system is represented by the inventory level *I*(*k*) and the

Assuming that the initial inventory level is *I*(0)= *I*0 and that the units in stock deteriorate at a rate *θ*, the dynamics of the first state variable, the inventory level, are governed by the fol‐

Also, and as mentioned in the introduction and previous paragraph, assuming a dependent demand rate and a stock-dependent demand rate with initial value *D*(0)=*D*0, the dynamics of the second state variable, the demand rate, are governed by the following difference

It is assumed that the firm has set for each period *k* the following targets: the production

*qD*, and *r* are incurred for a variable to deviate from its respective goal, then the objective

Since the target variables satisfy the dynamics (9) and (10), these can be rewritten using the

It is more convenient to write the model using a matrix notation. To this end, let

*qI* 0 0 *qD* .

*I*0 *D*0 , *Q* = ^

where *a* and *b* are positive constants, with 0<*a* <1, for a stable demand.

(*k*), the inventory goal level *I*

*<sup>J</sup>*(*P*, *<sup>I</sup>*, *<sup>D</sup>*)= <sup>1</sup>

where the shift operator Δ is defined by ∆ *f* (*k*)= *f* (*k*) - *f*

<sup>2</sup> ∑ *k*=0 *N*

*I*(*k* + 1)= *α I*(*k*) + *P*(*k*) - *D*(*k*) (9)

Optimal Control of Integrated Production – Forecasting System 121

*D*(*k* + 1)=*a D*(*k*) + *b I*(*k*) (10)

*qI*∆<sup>2</sup> *<sup>I</sup>*(*k*) <sup>+</sup> *qD*∆<sup>2</sup> *<sup>D</sup>*(*k*) <sup>+</sup> *<sup>r</sup>* <sup>∆</sup><sup>2</sup> *<sup>P</sup>*(*k*) (11)

^ (*k*). If penalties *qI* ,

(*k*), and the demand goal rate *D*

^ (*k*).

∆*I*(*k* + 1)=*α* ∆ *I*(*k*) + ∆*P*(*k*) - ∆*D*(*k*) (12)

∆*D*(*k* + 1)=*a* ∆*D*(*k*) + *b*∆*I*(*k*) (13)

with known initial inventory level *I*(0)= *I*0 and

$$\frac{d}{dt}D(t) = aD(t) + bI(t) \tag{5}$$

with known initial demand rate *D*(0)=*D*0 and *a* <0 for a stable demand. Since the goals *I* ^ (*t*), *D* ^ (*t*), and *<sup>P</sup>* ^ (*t*) also follow the dynamics (4)-(5), we can use the shift operator defined above to express the state equations (4) and (5) as

$$\frac{d}{dt}\Delta I(t) = \Delta P(t) \text{ - } \Delta D(t) \text{ - } \theta \Delta I(t) \tag{6}$$

and

$$\frac{d}{dt}\Delta D(t) = a\Delta D(t) + b\Delta I(t) \tag{7}$$

The state equations (6)-(7) can also be written in matrix form as

$$\frac{d}{dt}X\left(t\right) = AX\left(t\right) + B\Delta P(t) \tag{8}$$

where *A*<sup>=</sup> -*<sup>θ</sup>* -1 *<sup>b</sup> <sup>a</sup>* , *<sup>B</sup>* <sup>=</sup> <sup>1</sup> <sup>0</sup> , with initial condition *<sup>X</sup>* (0)= *<sup>X</sup>*0.

#### **2.2. Periodic review integrated production model**

In the periodic review model, the time interval 0, *T* is divided into *N* subintervals of equal length. During period *k*, the plant manufactures units of some product at the control‐ lable rate *P*(*k*), and the state of the system is represented by the inventory level *I*(*k*) and the demand rate *D*(*k*).

Assuming that the initial inventory level is *I*(0)= *I*0 and that the units in stock deteriorate at a rate *θ*, the dynamics of the first state variable, the inventory level, are governed by the fol‐ lowing difference equation:

$$I(k+1) = \alpha \ I(k) + P(k) \text{ - } D(k) \tag{9}$$

where *α* =1 - *θ*.

min *P*(*t*)

120 Decision Support Systems

*X <sup>A</sup>*

*D*

and

^ (*t*), and *<sup>P</sup>*

^

where *A*<sup>=</sup> -*<sup>θ</sup>* -1

*<sup>b</sup> <sup>a</sup>* , *<sup>B</sup>* <sup>=</sup> <sup>1</sup>

where *<sup>H</sup>* <sup>=</sup> *<sup>h</sup>* <sup>0</sup>

*<sup>J</sup>* <sup>=</sup> <sup>1</sup> 2 *∫* 0 *T h* Δ<sup>2</sup>

*I*(*t*) + *K*Δ<sup>2</sup>

min *P*(*t*)

<sup>0</sup> *<sup>K</sup>* and *HT* <sup>=</sup>

with known initial inventory level *I*(0)= *I*0 and

to express the state equations (4) and (5) as

*d*

*d*

The state equations (6)-(7) can also be written in matrix form as

*d*

<sup>0</sup> , with initial condition *<sup>X</sup>* (0)= *<sup>X</sup>*0.

*<sup>J</sup>* <sup>=</sup> <sup>1</sup> 2 *∫* 0 *T* *D*(*t*) + *r*Δ<sup>2</sup>

To use a matrix notation, which is more convenient, let *<sup>X</sup>* (*t*)= <sup>Δ</sup>*I*(*t*)

<sup>2</sup> + *r*Δ<sup>2</sup>

<sup>2</sup> = *X <sup>T</sup> AX* . Then, the objective function (2) can be further rewritten as

*X* (*t*) *<sup>H</sup>*

.

inventory level and demand rate are governed by the following state equations

*h <sup>T</sup>* 0 0 *KT*

*d*

*d*

*P*(*t*) *dt* +

*P*(*t*) *dt* +

Two state equations are used to describe the dynamics of our system. The variations of the

with known initial demand rate *D*(0)=*D*0 and *a* <0 for a stable demand. Since the goals *I*

(*t*) also follow the dynamics (4)-(5), we can use the shift operator defined above

1

<sup>2</sup> *X* (*T* ) *HT*

*dt I*(*t*)=*P*(*t*) - *D*(*t*) - *θI*(*t*) (4)

*dt D*(*t*)=*aD*(*t*) + *bI*(*t*) (5)

*dt* Δ*I*(*t*)=Δ*P*(*t*) - Δ*D*(*t*) - *θ*Δ*I*(*t*) (6)

*dt* Δ*D*(*t*)=*a*Δ*D*(*t*) + *b*Δ*I*(*t*) (7)

*dt X* (*t*)= *AX* (*t*) + *B*Δ*P*(*t*) (8)

1 <sup>2</sup> {*<sup>h</sup> <sup>T</sup>*Δ<sup>2</sup>

*<sup>I</sup>*(*<sup>T</sup>* ) <sup>+</sup> *KT*Δ<sup>2</sup>

*D*(*T* )} (2)

<sup>Δ</sup>*D*(*t*) and let

^ (*t*),

<sup>2</sup> (3)

Also, and as mentioned in the introduction and previous paragraph, assuming a dependent demand rate and a stock-dependent demand rate with initial value *D*(0)=*D*0, the dynamics of the second state variable, the demand rate, are governed by the following difference equation:

$$D\begin{pmatrix} k+1 \end{pmatrix} = a \ D(k) + b \ I(k) \tag{10}$$

where *a* and *b* are positive constants, with 0<*a* <1, for a stable demand.

It is assumed that the firm has set for each period *k* the following targets: the production goal rate *P* ^ (*k*), the inventory goal level *I* ^ (*k*), and the demand goal rate *D* ^ (*k*). If penalties *qI* , *qD*, and *r* are incurred for a variable to deviate from its respective goal, then the objective function to minimize is given by:

$$J(P, I, D) = \frac{1}{2} \sum\_{k=0}^{N} \left[ q\_I \Delta^2 I(k) + q\_D \Delta^2 D(k) + r \Delta^2 P(k) \right] \tag{11}$$

where the shift operator Δ is defined by ∆ *f* (*k*)= *f* (*k*) - *f* ^ (*k*).

Since the target variables satisfy the dynamics (9) and (10), these can be rewritten using the shift operator Δ to get:

$$
\Delta I(k+1) = \alpha \,\,\Delta I(k) + \Delta P(k) \,\,\text{-}\,\Delta D(k) \tag{12}
$$

$$
\Delta D(k+1) = a \,\,\Delta D(k) + b \,\Delta I(k) \tag{13}
$$

It is more convenient to write the model using a matrix notation. To this end, let

*<sup>Z</sup>*(*k*)= *<sup>I</sup>*(*k*) *<sup>D</sup>*(*k*) , *<sup>Z</sup>* ^ (*k*)= *<sup>I</sup>* ^ (*k*) *D* ^ (*k*) , *Z*<sup>0</sup> <sup>=</sup> *I*0 *D*0 , *Q* = *qI* 0 0 *qD* . The objective function (11) becomes

$$f(P, I, D) = \frac{1}{2} \sum\_{k=0}^{N} \left[ \|\Delta Z(k)\|\|\_{Q}^{2} + r \,\,\Delta^{2} \, P(k) \right] \tag{14}$$


The next condition is the state equation ∇Λ(*t*)H(Δ*P*, *<sup>X</sup>* , <sup>Λ</sup>, *<sup>t</sup>*)= *<sup>d</sup>*

*d*

ward sweep method of Bryson and Ho [26], we let

Note that, using the control equation (17), the state equation (8) becomes

*dt <sup>X</sup>* (*t*)= *AX* (*t*) - *<sup>r</sup>* -1

Λ(*T* )=*HT X* (*T* ).

*Model Solution*

problem (P).

where *S*(*t*)=

*3.1.1. First solution approach*

*s*1(*t*) 0 <sup>0</sup> *<sup>s</sup>*2(*t*) .

the change of variable (20) yields

*d dt* <sup>Λ</sup>(*t*)= *<sup>d</sup>*

*d*

matrix equations. Let

(8). Finally, the last condition is given by the initial and terminal conditions *X* (0)= *X*<sup>0</sup> and

We propose the following two equivalent solution approaches to solve the optimal control

In this approach, we need to solve a Riccati equation as we will see below. To use the back‐

Differentiating (20) with respect to *t* and then using successively the state equation (19) and

To solve Riccati equation (23), we use a change of variable to reduce it to a pair of linear

*dt <sup>S</sup>*(*t*) <sup>+</sup> *<sup>S</sup>*(*t*)*<sup>A</sup>* - *<sup>r</sup>* -1

Also, using the change of variable (20), the adjoint equation (18) becomes

Equating expressions (21) and (22), we obtain the following Riccati equation

*dt <sup>S</sup>*(*t*)= - *<sup>H</sup>* - *<sup>S</sup>*(*t*)*<sup>A</sup>* - *<sup>A</sup><sup>T</sup> <sup>S</sup>*(*t*) <sup>+</sup> *<sup>r</sup>* -1

*d*

*dt* <sup>Λ</sup>(*t*)=*HX* (*t*) <sup>+</sup> *<sup>A</sup><sup>T</sup>* <sup>Λ</sup>(t) (18)

*dt X* (*t*) which is equivalent to

*BB <sup>T</sup>* Λ(*t*) (19)

Optimal Control of Integrated Production – Forecasting System 123

Λ(*t*)=*S*(*t*)*X* (*t*) (20)

*dt* <sup>Λ</sup>(*t*)= -*<sup>H</sup>* - *<sup>A</sup><sup>T</sup> <sup>S</sup>*(*t*) *<sup>X</sup>* (*t*) (22)

*S*(*t*)*BB <sup>T</sup> S*(*t*) *X* (*t*) (21)

*S*(*t*)*BB <sup>T</sup> S*(*t*) (23)

where *X <sup>A</sup>* <sup>2</sup> = *X <sup>T</sup> A X* while the dynamics (12)-(13) become

$$
\Delta Z(k+1) = A\_1 \Delta Z(k) + B\_1 \Delta P(k) \tag{15}
$$

where *A*<sup>=</sup> *<sup>α</sup>* -1 *b a* and *<sup>B</sup>* <sup>=</sup> <sup>1</sup> 0 .

Thus we need to determine the production rates *P*(*k*) at each sample that minimize the ob‐ jective function (14), subject to the state equation (15).

#### **3. Optimal control**

#### **3.1. Optimal control of continuous review model**

Given the preceding definitions, the optimal control problem is to minimize the objective function (3) subject to the state equation (8):

$$\text{(P)} \quad \begin{cases} \min\_{P(t)} J = \frac{1}{2} \int\_{0}^{T} \|\mathbb{I}\left[X\left(t\right)\right]\|\_{H}^{2} + r\Delta^{2}P(t)\Big\|dt + \frac{1}{2} \|\mathbb{I}\left[X\left(T\right)\right]\|\_{H\_{T}}^{2} \\\\ \text{subject to} \\\\ \frac{d}{dt}X\left(t\right) = A\Delta X\left(t\right) + B\Delta P(t), \ \text{ $X\left(0\right) = X\_{0}$ } \end{cases}$$

To use Pontryagin principle, see for example the reference [25], we introduce the Hamiltoni‐ an

$$\mathsf{H}(\Delta P, \ X, \ \Lambda, \ t) = \frac{1}{2} \mathsf{I} \left[ \|X(t)\|\_{H}^{2} + r\Lambda^{2}P(t) \right] + \Lambda^{T}(t) \mathsf{I} \left[ AX(t) + B\Lambda P(t) \right] \tag{16}$$

where Λ(*t*)= *λ*1(*t*) *<sup>λ</sup>*2(*t*) is the adjoint function associated with the constraint (8). An optimal sol‐ ution to the control problem (P) satisfies several conditions. The first condition is the control equation ∇Δ*<sup>P</sup>*(*t*)H(Δ*P*, *X* , Λ, *t*)=0 which is equivalent to

$$\left(r\Delta P(t) + B^{\top}\Lambda(t) = 0\right) \tag{17}$$

The second condition is the adjoint equation ∇Δ*<sup>X</sup>* (*t*)H(Δ*P*, *<sup>X</sup>* , <sup>Λ</sup>, *<sup>t</sup>*)= - *<sup>d</sup> dt* Λ(*t*) which is equivalent to

$$-\frac{d}{dt}\Lambda(t) = HX\left(t\right) + A\left(^{T}\Lambda(t)\right)\tag{18}$$

The next condition is the state equation ∇Λ(*t*)H(Δ*P*, *<sup>X</sup>* , <sup>Λ</sup>, *<sup>t</sup>*)= *<sup>d</sup> dt X* (*t*) which is equivalent to (8). Finally, the last condition is given by the initial and terminal conditions *X* (0)= *X*<sup>0</sup> and Λ(*T* )=*HT X* (*T* ).

Note that, using the control equation (17), the state equation (8) becomes

$$\frac{d}{dt}X\left(t\right) = AX\left(t\right) \cdot r^{-1}BB^{\;T}\Lambda(t) \tag{19}$$

#### *Model Solution*

The objective function (11) becomes

where *X <sup>A</sup>*

122 Decision Support Systems

where *A*<sup>=</sup> *<sup>α</sup>* -1

*b a*

**3. Optimal control**

*<sup>J</sup>* <sup>=</sup> <sup>1</sup> 2 *∫* 0 *T*

*d*

(P) {

an

min *P*(*t*)

where Λ(*t*)=

equivalent to

and *<sup>B</sup>* <sup>=</sup> <sup>1</sup>

0 .

jective function (14), subject to the state equation (15).

**3.1. Optimal control of continuous review model**

<sup>2</sup> + *r*Δ<sup>2</sup>

subject to

*dt X* (*t*)= *A*Δ*X* (*t*) + *B*Δ*P*(*t*), *X* (0)= *X*<sup>0</sup>

equation ∇Δ*<sup>P</sup>*(*t*)H(Δ*P*, *X* , Λ, *t*)=0 which is equivalent to

*P*(*t*) *dt* +

<sup>2</sup> *X* (*t*) *<sup>H</sup>*

1

<sup>2</sup> + *r*Δ<sup>2</sup>

The second condition is the adjoint equation ∇Δ*<sup>X</sup>* (*t*)H(Δ*P*, *<sup>X</sup>* , <sup>Λ</sup>, *<sup>t</sup>*)= - *<sup>d</sup>*

function (3) subject to the state equation (8):

*X* (*t*) *<sup>H</sup>*

<sup>H</sup>(Δ*P*, *<sup>X</sup>* , <sup>Λ</sup>, *<sup>t</sup>*)= <sup>1</sup>

*λ*1(*t*)

*<sup>J</sup>*(*P*, *<sup>I</sup>*, *<sup>D</sup>*)= <sup>1</sup>

<sup>2</sup> ∑ *k*=0 *N*

<sup>2</sup> = *X <sup>T</sup> A X* while the dynamics (12)-(13) become

∆*Z*(*k*) *<sup>Q</sup>*

Thus we need to determine the production rates *P*(*k*) at each sample that minimize the ob‐

Given the preceding definitions, the optimal control problem is to minimize the objective

<sup>2</sup> *X* (*T* ) *HT* 2

To use Pontryagin principle, see for example the reference [25], we introduce the Hamiltoni‐

ution to the control problem (P) satisfies several conditions. The first condition is the control

*<sup>λ</sup>*2(*t*) is the adjoint function associated with the constraint (8). An optimal sol‐

<sup>2</sup> + *r* ∆<sup>2</sup> *P*(*k*) (14)

*P*(*t*) + ΛT(*t*) *AX* (*t*) + *B*Δ*P*(*t*) (16)

*dt* Λ(*t*) which is

*r*Δ*P*(*t*) + *B <sup>T</sup>* Λ(t)=0 (17)

∆*Z*(*k* + 1)= *A* ∆*Z*(*k*) + *B* ∆*P*(*k*) (15)

We propose the following two equivalent solution approaches to solve the optimal control problem (P).

#### *3.1.1. First solution approach*

In this approach, we need to solve a Riccati equation as we will see below. To use the back‐ ward sweep method of Bryson and Ho [26], we let

$$
\Lambda(t) = S(t)X(t) \tag{20}
$$

where *S*(*t*)= *s*1(*t*) 0 <sup>0</sup> *<sup>s</sup>*2(*t*) .

Differentiating (20) with respect to *t* and then using successively the state equation (19) and the change of variable (20) yields

$$\frac{d}{dt}\Delta(t) = \left[\frac{d}{dt}S(t) + S(t)A - r^{-1}S(t)BB^{\top}S(t)\right]X(t)\tag{21}$$

Also, using the change of variable (20), the adjoint equation (18) becomes

$$\frac{d}{dt}\Lambda(t) = \begin{bmatrix} \cdot\\ \cdot \end{bmatrix} \begin{bmatrix} \cdot\\ \cdot \end{bmatrix} \begin{bmatrix} S \ \end{bmatrix} \begin{bmatrix} \end{bmatrix} \begin{bmatrix} \end{bmatrix} \tag{22}$$

Equating expressions (21) and (22), we obtain the following Riccati equation

$$\frac{d}{dt}S(t) = -\boldsymbol{H} - S(t)\boldsymbol{A} - \boldsymbol{A}^T S(t) + \boldsymbol{r}^{-1} S(t)\boldsymbol{B}\boldsymbol{B}^T S(t) \tag{23}$$

To solve Riccati equation (23), we use a change of variable to reduce it to a pair of linear matrix equations. Let

$$S(t) = E(t)F^{-1}(t)\tag{24}$$

Γ=


*e* <sup>Γ</sup>*<sup>t</sup>* =*Pe Dt*

*I*(*t*)= *I* ^

*d dt*

and

∆*I*(*t*)

∆*P*(*t*)= - *r* -1

*P* -1

of the problem (P) are:

<sup>∆</sup>*D*(*t*) =(*<sup>A</sup>* - *<sup>r</sup>* -1

(*t*) + ∆*I*(*t*), *D*(*t*)=*D*

with initial condition *<sup>X</sup>* (0)= <sup>∆</sup>*I*(0)

*3.1.2. Second solution approach*


*BB <sup>T</sup>* ⋮ *A*

=

*θ* -*b* -*h* 0 1 -*a* 0 -*K* -*r*-1 0 -*θ* -1 0 0 *b a*

*G*(*t*)=*e*

is the matrix whose columns are the corresponding eigenvectors. Thus,

*G*(*t*)= *Pe*

^ (*t*) <sup>+</sup> <sup>∆</sup>*D*(*t*), *P*(*t*)=*<sup>P</sup>*


∆*D*(0)

state equation (19) are equivalent to the vector-matrix state equation

<sup>∆</sup>*D*(*t*) .

where ∆*I*(*t*) and ∆*I*(*t*) are solutions of the linear equation:

*B <sup>T</sup> E*(*t*)*F* (*t*)

*<sup>B</sup> <sup>T</sup> <sup>E</sup>*(*t*)*<sup>F</sup>* (*t*)-1 <sup>∆</sup>*I*(*t*)

*F* (*N* )= *I*4 where *I*4 denotes the identity matrix of order 4. The linear equations (32) can be solved in terms of a matrix exponential. The homogeneous set of equations has the solution

We recall that for a given a constant matrix Γ, the matrix exponential *e* <sup>Γ</sup>*<sup>t</sup>* is found as

In the next solution approach, which also leads to a set of homogeneous equations, we will show how the constant term *G*(*t*0) is obtained. For this approach, after finding *E*(*t*)and *F* (*t*), the desired result *S*(*t*) is obtained by using the change of variable (24). The optimal solutions

This approach avoids the Riccati equation as we will see. The adjoint equation (18) and the

^(*t*) <sup>+</sup> <sup>∆</sup>*P*(*t*),

, where *D* is the diagonal matrix whose elements are the eigenvalues of Γ and *P*

Γ(*t*-*t*0)

*D*(*t*-*t*0) *P* -1 (*N* )=*H* or *E*(*N* )=*H* and

*G*(*t*0) (32)

Optimal Control of Integrated Production – Forecasting System 125

*G*(*t*0) (33)

The boundary conditions *S*(*N* )=*H* are equivalent to *E*(*N* )*F* -1

The Riccati equation (23) becomes

$$\frac{d}{dt}\mathbb{E}\{t\}\mathbb{F}^{-1}(t) - \mathbb{E}\{t\}\mathbb{F}^{-1}(t)\frac{d}{dt}\mathbb{F}(t)\mathbb{F}^{-1}(t) = -H - \mathbb{E}(t)\mathbb{F}^{-1}(t)A - A^\top\mathbb{E}(t)\mathbb{F}^{-1}(t) + r^{-1}\mathbb{E}(t)\mathbb{F}^{-1}(t)\mathbb{B}\mathbb{B}^\top\mathbb{E}(t)\mathbb{F}^{-1}(t) \tag{25}$$

Multiplying this expression from the right by *F* yields

$$\frac{d}{dt}E\{t\} \text{ - } E\{t\}F^{-1}\{t\} \frac{d}{dt}F\{t\} = \text{ - } HF \text{ - } E\{t\}F^{-1}\{t\}AF \text{ - } A^T E\{t\} + r^{-1}E\{t\}F^{-1}\{t\}BB^T E\{t\} \tag{26}$$

Now, set

$$\frac{d}{dt}E\,\mathrm{(t)} = -\,HF\,\mathrm{(t)} - A\,^TE\,\mathrm{(t)}\tag{27}$$

Then, we have

$$E\{t\}F^{-1}(t)\frac{d}{dt}F(t) = E\{t\}F^{-1}(t)AF(t) \text{ - }r^{-1}E\{t\}F^{-1}(t)BB^{\top}E(t) \tag{28}$$

Multiplying this expression from the left by (*EF* -1)-1 yields

$$\frac{d}{dt}F\left(t\right) = AF\left(t\right) - r^{-1}BB^{-T}E\left(t\right)\tag{29}$$

Equations (27) and (29) now give two sets of linear equations

$$
\begin{bmatrix}
\frac{d}{dt}E(t) \\
\cdots \\
\frac{d}{dt}F(t)
\end{bmatrix} = \begin{bmatrix}
\cdot A^T & \vdots & \cdot H \\
\cdots & \cdots & \vdots & \cdots & \cdots \\
\cdot \cdot \,^{-1}B B \,^T & \vdots & A
\end{bmatrix} \begin{bmatrix}
E(t) \\
\cdots \\
\cdot \end{bmatrix} \tag{30}
$$

Call *G*(*t*)= *<sup>E</sup>*(*t*) *<sup>F</sup>* (*t*) . The differential equations (30) become of the form

$$\frac{d}{dt}G(t) = \Gamma G(t),\ G(t\_0) \text{ given}, \quad \Gamma \text{ constant} \tag{31}$$

where

$$
\Gamma = \begin{bmatrix}
\cdot \cdot A^T & \vdots & \cdot H \\
\cdot \cdot \cdot \cdot \cdots & \cdot \cdot & \cdot \cdot \cdots \cdots \\
\cdot \cdot r^{-1} B B^T & \vdots & A \\
\end{bmatrix} = \begin{bmatrix}
\theta & \cdot b & \cdot h & 0 \\
1 & \cdot a & 0 & \cdot K \\
\cdot r^{-1} & 0 & \cdot \theta & \cdot 1 \\
0 & 0 & b & a
\end{bmatrix}
$$

The boundary conditions *S*(*N* )=*H* are equivalent to *E*(*N* )*F* -1 (*N* )=*H* or *E*(*N* )=*H* and *F* (*N* )= *I*4 where *I*4 denotes the identity matrix of order 4. The linear equations (32) can be solved in terms of a matrix exponential. The homogeneous set of equations has the solution

$$G(t) = e^{\Gamma(t - t\_0)} \ G(t\_0) \tag{32}$$

We recall that for a given a constant matrix Γ, the matrix exponential *e* <sup>Γ</sup>*<sup>t</sup>* is found as *e* <sup>Γ</sup>*<sup>t</sup>* =*Pe Dt P* -1 , where *D* is the diagonal matrix whose elements are the eigenvalues of Γ and *P* is the matrix whose columns are the corresponding eigenvectors. Thus,

$$\mathbf{G}(t) = \mathbf{P}e^{D(t \cdot t\_0)} \mathbf{P}^{\cdot 1} \mathbf{G}(t\_0) \tag{33}$$

In the next solution approach, which also leads to a set of homogeneous equations, we will show how the constant term *G*(*t*0) is obtained. For this approach, after finding *E*(*t*)and *F* (*t*), the desired result *S*(*t*) is obtained by using the change of variable (24). The optimal solutions of the problem (P) are:

$$I(t) = \stackrel{\wedge}{I}(t) + \Delta I(t), \; D(t) = \stackrel{\wedge}{D}(t) + \Delta D(t), \; P(t) = \stackrel{\wedge}{P}(t) + \Delta P(t), \; t$$

where ∆*I*(*t*) and ∆*I*(*t*) are solutions of the linear equation:

$$
\frac{d}{dt} \begin{bmatrix} \Delta I(t) \\ \Delta D(t) \end{bmatrix} = \begin{pmatrix} A \ \cdot \ r^{-1} \boldsymbol{B}^{\top} E(t) F(t)^{-1} \end{pmatrix} \begin{bmatrix} \Delta I(t) \\ \Delta D(t) \end{bmatrix}
$$

with initial condition *<sup>X</sup>* (0)= <sup>∆</sup>*I*(0) ∆*D*(0)

and

*S*(*t*)=*E*(*t*)*F* -1

*dt <sup>F</sup>* (*t*)*<sup>F</sup>* -1(*t*)= - *<sup>H</sup>* - *<sup>E</sup>*(*t*)*<sup>F</sup>* -1(*t*)*<sup>A</sup>* - *<sup>A</sup><sup>T</sup> <sup>E</sup>*(*t*)*<sup>F</sup>* -1(*t*) <sup>+</sup> *<sup>r</sup>* -1

*dt <sup>F</sup>* (*t*)= - *HF* - *<sup>E</sup>*(*t*)*<sup>F</sup>* -1(*t*)*AF* - *<sup>A</sup><sup>T</sup> <sup>E</sup>*(*t*) <sup>+</sup> *<sup>r</sup>* -1

*dt <sup>F</sup>* (*t*)=*E*(*t*)*<sup>F</sup>* -1(*t*)*AF* (*t*) - *<sup>r</sup>* -1

*dt <sup>F</sup>*(*t*)= *AF* (*t*) - *<sup>r</sup>* -1


*BB <sup>T</sup>* ⋮ *A*

*<sup>F</sup>* (*t*) . The differential equations (30) become of the form

The Riccati equation (23) becomes

Multiplying this expression from the right by *F* yields

*d*

*dt <sup>E</sup>*(*t*)*<sup>F</sup>* -1(*t*) - *<sup>E</sup>*(*t*)*<sup>F</sup>* -1(*t*) *<sup>d</sup>*

*dt <sup>E</sup>*(*t*) - *<sup>E</sup>*(*t*)*<sup>F</sup>* -1(*t*) *<sup>d</sup>*

*<sup>E</sup>*(*t*)*<sup>F</sup>* -1(*t*) *<sup>d</sup>*

Multiplying this expression from the left by (*EF* -1)-1

*d dt E*(*t*) ⋯ ⋯ ⋯ *d dt F* (*t*)

*d*

*d*

Equations (27) and (29) now give two sets of linear equations

=


*d*

124 Decision Support Systems

*d*

Now, set

Then, we have

Call *G*(*t*)= *<sup>E</sup>*(*t*)

where

(*t*) (24)

(*t*)*BB <sup>T</sup> E*(*t*)*F* -1

(*t*) (25)

(30)

(*t*)*BB <sup>T</sup> E*(*t*) (26)

(*t*)*BB <sup>T</sup> E*(*t*) (28)

*BB <sup>T</sup> E*(*t*) (29)

*E*(*t*)*F* -1

*E*(*t*)*F* -1

*dt <sup>E</sup>*(*t*)= - *HF* (*t*) - *<sup>A</sup><sup>T</sup> <sup>E</sup>*(*t*) (27)

*E*(*t*)*F* -1

yields

*E*(*t*) ⋯ ⋯ ⋯ *F*(*t*)

*dt G*(*t*)=Γ*G*(*t*), *G*(*t*0) given, Γ constant (31)

$$
\Delta P(t) = -r^{-1} \mathcal{B}^{\ \ \ \Gamma} E(t) F(t)^{-1} \begin{bmatrix} \Delta I(t) \\ \Delta D(t) \end{bmatrix}.
$$

#### *3.1.2. Second solution approach*

This approach avoids the Riccati equation as we will see. The adjoint equation (18) and the state equation (19) are equivalent to the vector-matrix state equation

$$
\begin{bmatrix}
\frac{d}{dt}X\{t\} \\
\cdots \\
\frac{d}{dt}\Lambda(t)
\end{bmatrix} = \begin{bmatrix}
A & \vdots & \text{\textquotedbl{}}\,^{-1}BB^{\text{T}} \\
\cdots & \cdots & \vdots & \cdots & \cdots \\
\cdot\,^{d}H & \vdots & \text{\textquotedbl{}}\,^{-1}\Lambda^{\text{T}}
\end{bmatrix} \begin{bmatrix}
X\{t\} \\
\Lambda(t)
\end{bmatrix} \tag{34}
$$

Let *Z*(*t*)= *<sup>X</sup>* (*t*) <sup>Λ</sup>(*t*) . Then, the vector-matrix state equation (32) can be rewritten as

$$\frac{d}{dt}Z\left(t\right) = \Phi Z\left(t\right) \tag{35}$$

To determine *Z*(0)= *<sup>X</sup>* (0)

which can be rewritten as

φ1(*T* ) φ2(*T* ) φ3(*T* ) φ4(*T* )

> φ1(*T* ) φ2(*T* ) φ3(*T* ) φ4(*T* )

Using the terminal condition

*X* (0) Λ(0)

> *X* (0) Λ(0)

> > (*T* ) - φ<sup>4</sup>

^ (*t*) <sup>+</sup> <sup>∆</sup>*D*(*t*), *P*(*t*)=*<sup>P</sup>*

where ∆*I*(*t*) and ∆*I*(*t*) are solutions of the linear equation:

(*t*)Λ(0)

(*T* ))-1(φ<sup>3</sup>

^(*t*) <sup>+</sup> <sup>∆</sup>*P*(*t*),

Here also we assume that the system state is available during each period *k*. To use the standard Lagrangian technique, we introduce the discrete Lagrange multiplier vector:

Finally, the optimal solutions of the problem (P) using the second method are:

Λ(0)=(*HT* φ<sup>2</sup>

+ *φ*<sup>2</sup>

(*t*)Λ(0)

*t* =*T* ,

*X* (*T* ) <sup>Λ</sup>(*<sup>T</sup>* ) <sup>=</sup>

we have

*I*(*t*)= *I* ^

> ∆*I*(*t*) <sup>∆</sup>*D*(*t*) <sup>=</sup>*φ*<sup>1</sup>

Λ(*t*)=*φ*<sup>3</sup>

∆*P*(*t*)= - *r* -1

and

Λ(*k*)=

(*t*)

*λ<sup>I</sup>* (*k*) *λD*(*k*)

*X* (*T* ) *HT <sup>X</sup>* (*<sup>T</sup>* ) <sup>=</sup>

*Z*(*T* )=φ(*T* )*Z*(0)

Λ(*T* )=*HT X* (*T* )

from which, it follows

(*t*) + ∆*I*(*t*), *D*(*t*)=*D*

∆*I*(0) ∆*D*(0)

+ *φ*<sup>4</sup>

**3.2. Optimal control of periodic review model**

Then, the Lagrangian function is given by

(*t*)

∆*I*(0) ∆*D*(0)

*B <sup>T</sup>* Λ(*t*) .

<sup>Λ</sup>(0) , we recall that *<sup>X</sup>* (0)= <sup>Δ</sup>*I*(0)

However, using the final value Λ(*T* ), we can find Λ(0) as follows. From (36), we have at

<sup>Δ</sup>*D*(0) is known while Λ(0)=

Optimal Control of Integrated Production – Forecasting System 127

(*T* ) - *HT* φ1(*T* ))*X* (0) (38)

λ1(0)

<sup>λ</sup>2(0) is not.

where

$$\mathbf{Q} = \begin{bmatrix} \mathbf{A} & \vdots & \mathbf{-}r^{-1}BB^{\top} \\ \cdots & \cdots & \cdots \\ \cdot H & \vdots & \cdots \end{bmatrix} \begin{bmatrix} \mathbf{-}\boldsymbol{\Theta} & \mathbf{-1} & \mathbf{-}r^{-1} & \mathbf{0} \\ \mathbf{b} & \mathbf{a} & \mathbf{0} & \mathbf{0} \\ \cdot \boldsymbol{h} & \mathbf{0} & \boldsymbol{\Theta} & \mathbf{-}b \\ \mathbf{0} & \cdot\mathbf{K} & \mathbf{1} & \cdots \end{bmatrix}$$

Expression (35) is a set of 4 first-order homogeneous differential equations with constant co‐ efficients. It is similar to (31) and it has a solution similar to (33). We give here the explicit solution.

The matrix Φ has four distinct eigenvalues *mi* , *i* =1,2, 3,4. The explicit expressions of these eigenvalues are easily obtained using some mathematical software with symbolic computa‐ tion capabilities such as MATHCAD, MAPLE, or MATLAB. These expressions are lengthy and thus are not reproduced here. The explicit expressions of the corresponding eigenvec‐ tors are also obtained using the same software. Then, the solution to (33) is given by

$$
\nabla Z(t) = \!\!\!\varphi(t) Z(0) \tag{36}
$$

Note that the first two components of *Z*(*t*) thus computed form the state vector *X* (*t*) whose components are Δ*I*(*t*) and Δ*D*(*t*), while the last two components form the co-state vector Λ(*t*) whose components are *λ*<sup>1</sup> (*t*) and *λ*<sup>2</sup> (*t*). In what follows, we show how φ(*t*) and *Z*(0) are determined using the two initial conditions *I*(0)= *I*0, *D*(0)=*D*<sup>0</sup> and the terminal conditions *λ*1 (*T* )=*h <sup>T</sup>*Δ*I*(*T* ), *λ*<sup>2</sup> (*T* )= *KT*Δ*D*(*T* ).

To determine φ(t), introduce the diagonal matrix *M* =*diag*(*m*1, *m*2, *m*3, *m*4) and denote by *Y* the matrix whose columns are the corresponding eigenvectors. Then,

φ(*t*)=*Y e Dt Y* -1 = ∑ *i*=1 4 *Y* (:, *i*)*Y* -1(*i*, :)*e mi <sup>t</sup>* (37)

where*Y* (:, *i*) is the *i*th column of *Y* and *Y* -1(*i*, :) is the *i*th row of *Y* -1 . To determine *Z*(0)= *<sup>X</sup>* (0) <sup>Λ</sup>(0) , we recall that *<sup>X</sup>* (0)= <sup>Δ</sup>*I*(0) <sup>Δ</sup>*D*(0) is known while Λ(0)= λ1(0) <sup>λ</sup>2(0) is not. However, using the final value Λ(*T* ), we can find Λ(0) as follows. From (36), we have at

*t* =*T* ,

*Z*(*T* )=φ(*T* )*Z*(0)

*d dt X* (*t*) ⋯ ⋯ ⋯ *d dt* Λ(*t*)

Let *Z*(*t*)= *<sup>X</sup>* (*t*)

126 Decision Support Systems

*A* ⋮ -*r* -1

⋯ ⋯ ⋯ ⋮ ⋯ ⋯ ⋯ -*H* ⋮ -*A<sup>T</sup>*

Λ(*t*) whose components are *λ*<sup>1</sup>

(*T* )=*h <sup>T</sup>*Δ*I*(*T* ), *λ*<sup>2</sup>

*BB <sup>T</sup>*

The matrix Φ has four distinct eigenvalues *mi*

=

where

Φ=

solution.

*λ*1

=

*d*


*A* ⋮ -*r*-1

⋯ ⋯ ⋯ ⋮ ⋯ ⋯ ⋯ -*H* ⋮ -*A<sup>T</sup>*

<sup>Λ</sup>(*t*) . Then, the vector-matrix state equation (32) can be rewritten as

Expression (35) is a set of 4 first-order homogeneous differential equations with constant co‐ efficients. It is similar to (31) and it has a solution similar to (33). We give here the explicit

eigenvalues are easily obtained using some mathematical software with symbolic computa‐ tion capabilities such as MATHCAD, MAPLE, or MATLAB. These expressions are lengthy and thus are not reproduced here. The explicit expressions of the corresponding eigenvec‐

Note that the first two components of *Z*(*t*) thus computed form the state vector *X* (*t*) whose components are Δ*I*(*t*) and Δ*D*(*t*), while the last two components form the co-state vector

determined using the two initial conditions *I*(0)= *I*0, *D*(0)=*D*<sup>0</sup> and the terminal conditions

To determine φ(t), introduce the diagonal matrix *M* =*diag*(*m*1, *m*2, *m*3, *m*4) and denote by *Y*

*Y* (:, *i*)*Y* -1(*i*, :)*e*

*mi*

tors are also obtained using the same software. Then, the solution to (33) is given by

(*t*) and *λ*<sup>2</sup>

the matrix whose columns are the corresponding eigenvectors. Then,

where*Y* (:, *i*) is the *i*th column of *Y* and *Y* -1(*i*, :) is the *i*th row of *Y* -1

*Y* -1 = ∑ *i*=1 4

(*T* )= *KT*Δ*D*(*T* ).

φ(*t*)=*Y e Dt*

*BB <sup>T</sup>*

*X* (*t*)

*dt Z*(*t*)=Φ*Z*(*t*) (35)

, *i* =1,2, 3,4. The explicit expressions of these

*Z*(*t*)=φ(*t*)*Z*(0) (36)

(*t*). In what follows, we show how φ(*t*) and *Z*(0) are

*<sup>t</sup>* (37)

.

<sup>Λ</sup>(*t*) (34)

which can be rewritten as

$$
\begin{bmatrix} X \begin{pmatrix} T \end{pmatrix} \\ \Lambda \begin{pmatrix} T \end{pmatrix} \end{bmatrix} = \begin{bmatrix} \varphi\_1(T) & \varphi\_2(T) \\ \varphi\_3(T) & \varphi\_4(T) \end{bmatrix} \begin{bmatrix} X \begin{pmatrix} 0 \end{pmatrix} \\ \Lambda \begin{pmatrix} 0 \end{pmatrix} \end{bmatrix}
$$

Using the terminal condition

Λ(*T* )=*HT X* (*T* )

we have

$$
\begin{bmatrix} X \begin{pmatrix} T \end{pmatrix} \\ H\_T X \begin{pmatrix} T \end{pmatrix} \end{bmatrix} = \begin{bmatrix} \varphi\_1(T) & \varphi\_2(T) \\ \varphi\_3(T) & \varphi\_4(T) \end{bmatrix} \begin{bmatrix} X \begin{pmatrix} 0 \end{pmatrix} \\ \Lambda(0) \end{bmatrix}
$$

from which, it follows

$$\Lambda(0) = \left(H\_T \varphi\_2(T) \text{ - } \varphi\_4(T)\right) \text{-} \left(\varphi\_3(T) \text{ - } H\_T \varphi\_1(T)\right) X \text{ (0)}\tag{38}$$

Finally, the optimal solutions of the problem (P) using the second method are:

$$I(t) = \stackrel{\wedge}{I}(t) + \Delta I(t), \\ D(t) = \stackrel{\wedge}{D}(t) + \Delta D(t), \\ P(t) = \stackrel{\wedge}{P}(t) + \Delta P(t),$$

where ∆*I*(*t*) and ∆*I*(*t*) are solutions of the linear equation:

$$
\begin{bmatrix}
\Delta I(t) \\
\Delta D(t)
\end{bmatrix} = \varphi\_1(t) \begin{bmatrix}
\Delta I(0) \\
\Delta D(0)
\end{bmatrix} + \varphi\_2(t) \Lambda(0)
$$
 
$$
\Lambda(t) = \varphi\_3(t) \begin{bmatrix}
\Delta I(0) \\
\Delta D(0)
\end{bmatrix} + \varphi\_4(t) \Lambda(0)
$$

and

$$
\Delta P(t) = -r^{-1} \mathcal{B}^{\,T} \Lambda(t) \, \dots
$$

#### **3.2. Optimal control of periodic review model**

Here also we assume that the system state is available during each period *k*. To use the standard Lagrangian technique, we introduce the discrete Lagrange multiplier vector:

Λ(*k*)= *λ<sup>I</sup>* (*k*) *λD*(*k*)

Then, the Lagrangian function is given by

$$\Delta(P, Z\_{\prime\prime}, \Lambda) = \sum\_{k=0}^{N} \left\{ \frac{1}{2} \left[ \left\| \Delta Z(k) \right\|\_{Q}^{2} + r\_{\prime} \Delta P(k)^{2} \right] + \Lambda(k+1)^{T} \left[ -\Delta Z(k+1) + A\_{\prime} \Delta Z(k) + B\_{\prime} \Delta P(k) \right] \right\} \tag{39}$$

The necessary optimality conditions are the control equation

$$\frac{\partial \mathcal{L}}{\partial \Delta P(k)} = 0 \iff r \,\,\Delta P(k) + B^T \,\,\Lambda(k+1) = 0\tag{40}$$

=Q ΔZ(k) + A*<sup>T</sup> S*(*k* + 1) *A* Δ*Z*(*k*) - *B V* (*k* + 1)Δ*Z*(*k*) , ={Q + A*<sup>T</sup> S*(*k* + 1)*A* - *A<sup>T</sup> S*(*k* + 1)*B V* (*k* + 1)}Δ*Z*(*k*).

*A<sup>T</sup> S*(*k* + 1)*B* 1 + *r* -1

duction rates, we use the change of variable (42) to get from the adjoint equation (41),

*B <sup>T</sup> S*(*k* + 1)*B* -1

Parameters Values Nonmonetary parameters Length of planning horizon *T* =2.5

Monetary parameters Penalty for production rate deviation *r* =0.1

where the optimal state vector Δ*Z*(*k*) is found in expression (49) above and *S*(*k*) is found in

Also, from the dynamics (15), we have the optimal production rates

^(*k*) - *<sup>r</sup>* -1 <sup>1</sup> <sup>+</sup> *<sup>r</sup>* -1

(DTRE) given by the recursive relation

*S*(*k*)=*Q* + A*<sup>T</sup> S*(*k* + 1)*A* - *r* -1

*P*(*k*)=*P*

**Table 1.** Data for continuous-review model

so that

expression (47).

Finally, the matrices *S* can be computed from the backward discrete time Ricatti equation

The boundary condition *S*(*N* )=*Q* follows from Δ*P*(*N* )=0. Now, to obtain the optimal pro‐

*B <sup>T</sup> S*(*k* + 1)*B* -1

*AS*(*k* + 1)Δ*Z*(*k* + 1)=*S*(*k*)Δ*Z*(*k*) - *A*Δ*Z*(*k*) (48)

Δ*Z*(*k* + 1)= *AS*(*k* + 1) -1 *S*(*k*) - *Q* Δ*Z*(*k*) (49)

Initial inventory level *I*<sup>0</sup> =5 Initial demand rate *D*<sup>0</sup> =2 Demand rate coefficient *a*=0.8 Inventory level coefficient *b*=1 Deterioration rate θ =0.1

Penalty for inventory level deviation *h* =4 Penalty for demand rate deviation *K* =10 Inventory salvage value *hT* =100 Demand salvage value *KT* =100

*B <sup>T</sup> S*(*k* + 1)*A* (47)

Optimal Control of Integrated Production – Forecasting System 129

*B <sup>T</sup> S*(*k* + 1)*A* Δ*Z*(*k*) (50)

and the adjoint equation

$$\frac{\partial \mathcal{L}}{\partial \Delta Z(k)} = 0 \iff Q \Delta Z(k) \text{ - } \Lambda(k) + A^T \Lambda(k+1) = 0. \tag{41}$$

In order to solve these equations, we use the backward sweep method of Bryson and Ho (1975), who treat extensively in their book the problem of optimal control and estimation. They detail two methods for solving the Riccati equation arising in linear optimal control problem, the first one being the transition matrix method and the second being the back‐ ward sweep method. Let

$$
\Lambda(k) = S\left(k\right) \,\Delta Z\left(k\right),
\tag{42}
$$

The control equation (40) becomes

$$\begin{split} \Delta P(k) &= -\boldsymbol{r}\cdot\boldsymbol{^{1}B}^{T}\boldsymbol{\Lambda}(k+1), = -\boldsymbol{r}\cdot\boldsymbol{^{1}B}^{T}\boldsymbol{S}(k+1)\boldsymbol{\Delta}Z(k+1), \\ &= -\boldsymbol{r}\cdot\boldsymbol{^{1}B}^{T}\boldsymbol{S}(k+1)[\boldsymbol{\Lambda}\,\boldsymbol{\Delta}Z(k) + \boldsymbol{B}\,\boldsymbol{\Delta}P(k)]. \end{split} \tag{43}$$

Solving for Δ*P*(*k*), we get

$$
\Delta P(k) = -V(k+1)\,\Delta Z(k)\tag{44}
$$

where

$$V(k+1) = r^{-1} \left[ 1 + r^{-1} B^{\ \ T} \ S(k+1) \ B \right]^{-1} B^{\ \ T} \ S(k+1) A^{\ \tau}$$

Now the adjoint equation (41) becomes

$$
\Delta\Lambda(k) = Q \,\Delta Z(k) + A^T \,\Lambda(k+1) \tag{45}
$$

so that

$$\begin{aligned} \mathbf{S}(k) \,\Delta Z(k) &= \mathbf{Q} \,\Delta Z(k) + \mathbf{A}^T \mathbf{S}(k+1) \boldsymbol{\Delta} Z(k+1), \\ \mathbf{S} &= \mathbf{Q} \,\Delta Z(k) + \mathbf{A}^T \mathbf{S}(k+1) \mathbf{I} \,\Delta Z(k) + \mathbf{B} \,\Delta P(k) \mathbf{J}, \end{aligned} \tag{46}$$

$$\begin{aligned} &= \mathbf{Q} \cdot \boldsymbol{\Delta} \mathbf{Z}(\mathbf{k}) + \mathbf{A}^{\top} \boldsymbol{\mathcal{S}}(\boldsymbol{k}+1) \mathbf{I} \boldsymbol{A} \cdot \boldsymbol{\Delta} \mathbf{Z}(\mathbf{k}) \mathbf{-} \boldsymbol{B} \, V(\mathbf{k}+1) \boldsymbol{\Delta} \mathbf{Z}(\mathbf{k}) \mathbf{J}, \\ &= \left\{ \mathbf{Q} + \mathbf{A}^{T} \boldsymbol{\mathcal{S}}(\boldsymbol{k}+1) \mathbf{A} \, -\boldsymbol{A}^{T} \boldsymbol{\mathcal{S}}(\boldsymbol{k}+1) \mathbf{B} \, V(\mathbf{k}+1) \right\} \boldsymbol{\Delta} \mathbf{Z}(\mathbf{k}). \end{aligned}$$

Finally, the matrices *S* can be computed from the backward discrete time Ricatti equation (DTRE) given by the recursive relation

$$S(k) = Q + A^T S(k+1)A - r^{-1}A^T \ S(k+1)B \left[1 + r^{-1}B^T S(k+1)B\right]^{-1} B^T \ S(k+1)A \tag{47}$$

The boundary condition *S*(*N* )=*Q* follows from Δ*P*(*N* )=0. Now, to obtain the optimal pro‐ duction rates, we use the change of variable (42) to get from the adjoint equation (41),

$$AS\,(k+1)\Delta Z\,(k+1) = S\,(k)\Delta Z\,(k) - A\Delta Z\,(k)\tag{48}$$

so that

L(*P*, *Z*, Λ)= ∑

128 Decision Support Systems

and the adjoint equation

ward sweep method. Let

Solving for Δ*P*(*k*), we get

*V* (*k* + 1)=*r* -1 1 + *r* -1

where

so that

The control equation (40) becomes

Δ*P*(*k*)= - *r* -1

= - *r* -1

*B <sup>T</sup> S*(*k* + 1) *B* -1

Now the adjoint equation (41) becomes

*k*=0 *N* { 1

<sup>2</sup> ∆*Z*(*k*) *<sup>Q</sup>*

The necessary optimality conditions are the control equation

∂ L

∂ L

<sup>2</sup> + *r* ∆*P*(*k*)<sup>2</sup> + Λ(*k* + 1)*<sup>T</sup>* -Δ*Z*(*k* + 1) + *A* Δ*Z*(*k*) + *B* Δ*P*(*k*) } (39)

<sup>∂</sup> <sup>∆</sup> *<sup>P</sup>*(*<sup>k</sup>* ) =0 <sup>⇔</sup>*<sup>r</sup>* <sup>∆</sup>*P*(*k*) <sup>+</sup> *<sup>B</sup> <sup>T</sup>* <sup>Λ</sup>(*<sup>k</sup>* <sup>+</sup> 1)=0 (40)

<sup>∂</sup> <sup>∆</sup> *<sup>Z</sup>* (*<sup>k</sup>* ) =0 <sup>⇔</sup>*<sup>Q</sup>* <sup>∆</sup>*Z*(*k*) - <sup>Λ</sup>(*k*) <sup>+</sup> *<sup>A</sup><sup>T</sup>* <sup>Λ</sup>(*<sup>k</sup>* <sup>+</sup> 1)=0. (41)

Λ(*k*)=*S*(*k*) Δ*Z*(*k*), (42)

*B <sup>T</sup> S*(*k* + 1)Δ*Z*(*k* + 1),

*<sup>B</sup> <sup>T</sup> <sup>S</sup>*(*<sup>k</sup>* <sup>+</sup> 1) <sup>A</sup> <sup>Δ</sup>*Z*(*k*) <sup>+</sup> *<sup>B</sup>* <sup>Δ</sup>*P*(*k*) . (43)

Δ*P*(*k*)= - *V* (*k* + 1) Δ*Z*(*k*) (44)

Λ(*k*)=*Q* Δ*Z*(*k*) + *A<sup>T</sup>* Λ(*k* + 1) (45)

=Q <sup>Δ</sup>Z(k) <sup>+</sup> <sup>A</sup>*<sup>T</sup> <sup>S</sup>*(*<sup>k</sup>* <sup>+</sup> 1) *<sup>A</sup>* <sup>Δ</sup>*Z*(*k*) <sup>+</sup> *<sup>B</sup>* <sup>Δ</sup>*P*(*k*) , (46)

In order to solve these equations, we use the backward sweep method of Bryson and Ho (1975), who treat extensively in their book the problem of optimal control and estimation. They detail two methods for solving the Riccati equation arising in linear optimal control problem, the first one being the transition matrix method and the second being the back‐

*B <sup>T</sup>* Λ(*k* + 1), = - *r* -1

*B <sup>T</sup> S*(*k* + 1)*A*

*S*(*k*) ΔZ(k)=Q ΔZ(k) + A*<sup>T</sup> S*(*k* + 1)Δ*Z*(*k* + 1),

$$
\Delta Z\,(k+1) = \mathsf{f}\,AS\,(k+1)\mathsf{I}^{\mathsf{I}}\mathsf{I}\,S\,(k) \cdot \mathsf{Q}\,\mathsf{I}\,\Delta Z\,(k)\tag{49}
$$

Also, from the dynamics (15), we have the optimal production rates

$$P(k) = \stackrel{\wedge}{P}(k) \cdot r^{-1} \stackrel{\wedge}{1} + r^{-1}B^{\;T} \; S(k+1)B \stackrel{\bullet}{J}^{-1}B^{\;T}S(k+1)A \; \Delta Z(k) \tag{50}$$

where the optimal state vector Δ*Z*(*k*) is found in expression (49) above and *S*(*k*) is found in expression (47).


**Table 1.** Data for continuous-review model
