Optimal Control of Integrated Production – Forecasting System

1. Introduction

The production planning problem has received much attention, and many sophisticated models and procedures have been developed to deal with this problem. Many other components 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 manufacturing 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 current demand is related to the previous demand in order to integrate the forecasting component 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 parameters 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 supply 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 for random demand with a seasonal trend. In [6], the authors assess the empirical stock control performance of intermittent demand estimation procedures. In [7], the authors study two modifications of the normal distribution, both taking non-negative values only.

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 backlogging 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-dependent 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 inventory 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 time-varying 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 investigate 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 authors 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 demand 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, classifying them on the basis of shelf-life characteristic and demand variations. In [23], the authors 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 minimization 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 develop two deterministic economic production quantity (EPQ) models for Weibull-distributed 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 policy 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 production goal rate, to build the objective function of our model. The inventory goal level is a safety 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 systems. 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

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 $\widehat{I}(t)$, $\widehat{D}(t)$, and $\widehat{P}(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

$$\underset{P\left(t\right)}{\min }J=\frac{1}{2}\underset{0}{\overset{T}{\int }}\left\{h{\left[I\left(t\right)-\widehat{I}\left(t\right)\right]}^2+K{\left[D\left(t\right)-\widehat{D}\left(t\right)\right]}^2+r{\left[P\left(t\right)-\widehat{P}\left(t\right)\right]}^2\right\}dt +\frac{1}{2}\left\{h_T{\left[I\left(T\right)-\widehat{I}\left(T\right)\right]}^2+K_T{\left[D\left(T\right)-\widehat{D}\left(T\right)\right]}^2\right\}$$E1

In (1), the expression $\frac{1}{2}\left\{h_T{\left[I\left(T\right)-\widehat{I}\left(T\right)\right]}^2+K_T{\left[D\left(T\right)-\widehat{D}\left(T\right)\right]}^2\right\}$ gives the salvage value of the ending state. Using the shift operator defined by $\mathrm{\Delta }f\left(t\right)=f(t)-\widehat{f}(t)$, the objective function is expressed as

$$\underset{P(t)}{\min }J=\frac{1}{2}\underset{0}{\overset{T}{\int }}\left[h{\mathrm{\Delta }}^{2}I\left(t\right)+K{\mathrm{\Delta }}^{2}D\left(t\right)+r{\mathrm{\Delta }}^{2}P(t)\right]dt+\frac{1}{2}\left\{h_T{\mathrm{\Delta }}^{2}I\left(T\right)+K_T{\mathrm{\Delta }}^{2}D\left(T\right)\right\}$$E2

To use a matrix notation, which is more convenient, let $X(t)=\left[\begin{array}{c}\mathrm{\Delta }I(t)\\ \mathrm{\Delta }D(t)\end{array}\right]$ and let ${‖X‖}_A^2=X^{T}AX$. Then, the objective function (2) can be further rewritten as

$$\underset{P(t)}{\min }J=\frac{1}{2}\underset{0}{\overset{T}{\int }}\left[{‖X(t)‖}_H^2+r{\mathrm{\Delta }}^{2}P(t)\right]dt+\frac{1}{2}{‖X(T)‖}_{H_T}^2$$E3

where $H=\left[\begin{array}{cc}h& 0\\ 0& K\end{array}\right]$ and $H_T=\left[\begin{array}{cc}{h}_T& 0\\ 0& K_T\end{array}\right]$.

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\left(t\right)=P\left(t\right)-D\left(t\right)-\theta I(t)$$E4

with known initial inventory level $I(0)=I_0$ and

$$\frac{d}{dt}D\left(t\right)=aD\left(t\right)+bI(t)$$E5

with known initial demand rate $D(0)=D_0$ and $a<0$ for a stable demand. Since the goals $\widehat{I}(t)$, $\widehat{D}(t)$, and $\widehat{P}(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}\mathrm{\Delta }I\left(t\right)=\mathrm{\Delta }P\left(t\right)-\mathrm{\Delta }D\left(t\right)-\theta \mathrm{\Delta }I(t)$$E6

and

$$\frac{d}{dt}\mathrm{\Delta }D\left(t\right)=a\mathrm{\Delta }D\left(t\right)+b\mathrm{\Delta }I(t)$$E7

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\mathrm{\Delta }P(t)$$E8

where $A=\left[\begin{array}{cc}-\theta & -1\\ b& a\end{array}\right]$, $B=\left[\begin{array}{c}1\\ 0\end{array}\right]$, with initial condition $X\left(0\right)=X_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 controllable 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 $\theta$, the dynamics of the first state variable, the inventory level, are governed by the following difference equation:

$$I\left(k+1\right)= \alpha I\left(k\right)+P\left(k\right)-D(k)$$E9

where $\alpha =1-\theta$.

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\left(k+1\right)=a D\left(k\right)+b I(k)$$E10

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 $\widehat{P}(k)$, the inventory goal level $\widehat{I}(k)$, and the demand goal rate $\widehat{D}(k)$. If penalties $q_I$, $q_D$, and r are incurred for a variable to deviate from its respective goal, then the objective function to minimize is given by:

$$J\left(P,I,D\right)=\frac{1}{2}\sum _{k=0}^N\left[q_I{∆}^{2}I\left(k\right)+q_D{∆}^{2}D\left(k\right)+r {∆}^{2}P(k)\right]$$E11

where the shift operator Δ is defined by $∆f\left(k\right)=f\left(k\right)-\widehat{f}(k)$.

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

$$∆I\left(k+1\right)=\alpha ∆I\left(k\right)+∆P\left(k\right)-∆D(k)$$E12

$$∆D\left(k+1\right)=a ∆D\left(k\right)+b∆I(k)$$E13

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

$Z\left(k\right)=\left[\begin{array}{c}I(k)\\ D(k)\end{array}\right]$, $\widehat{Z}\left(k\right)=\left[\begin{array}{c}\widehat{I}(k)\\ \widehat{D}(k)\end{array}\right]$, $Z_0=\left[\begin{array}{c}{I}_0\\ D_0\end{array}\right]$, $Q=\left[\begin{array}{cc}{q}_I& 0\\ 0& q_D\end{array}\right]$.

The objective function (11) becomes

$$J\left(P,I,D\right)=\frac{1}{2}\sum _{k=0}^N\left[{‖∆Z(k)‖}_Q^2+r {∆}^{2}P(k)\right]$$E14

where ${‖X‖}_A^2=X^{T}A X$ while the dynamics (12)-(13) become

$$∆Z\left(k+1\right)=A ∆Z\left(k\right)+B ∆P(k)$$E15

where $A=\left[\begin{array}{cc}\alpha & -1\\ b& a\end{array}\right]$ and $B=\left[\begin{array}{c}1\\ 0\end{array}\right]$.

Thus we need to determine the production rates $P(k)$ at each sample that minimize the objective 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):

$\begin{array}{cc}(\mathcal{P})& \left\{\begin{array}{c}\underset{P(t)}{\min }J=\frac{1}{2}\underset{0}{\overset{T}{\int }}\left[{‖X(t)‖}_H^2+r{\mathrm{\Delta }}^{2}P(t)\right]dt+\frac{1}{2}{‖X(T)‖}_{H_T}^2\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{ }\mathrm{t}\mathrm{o}\\ \frac{d}{dt}X\left(t\right)=A\mathrm{\Delta }X\left(t\right)+B\mathrm{\Delta }P\left(t\right), X\left(0\right)=X_0\end{array}\right.\end{array}$

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

$$\mathcal{H}\left(\mathrm{\Delta }P,X,\mathrm{\Lambda },t\right)=\frac{1}{2}\left[{‖X(t)‖}_H^2+r{\mathrm{\Delta }}^{2}P(t)\right]+{\mathrm{\Lambda }}^{\mathrm{T}}(t)\left[AX\left(t\right)+B\mathrm{\Delta }P(t)\right]$$E16

where $\mathrm{\Lambda }\left(t\right)=\left[\begin{array}{c}{\lambda }_1(t)\\ {\lambda }_2(t)\end{array}\right]$ is the adjoint function associated with the constraint (8). An optimal solution to the control problem $(\mathcal{P})$ satisfies several conditions. The first condition is the control equation ${\nabla }_{\mathrm{\Delta }P(t)}\mathcal{H}\left(\mathrm{\Delta }P,X,\mathrm{\Lambda },t\right)=0$ which is equivalent to

$$r\mathrm{\Delta }P(t)+B^T\mathrm{\Lambda }(\mathrm{t})=0$$E17

The second condition is the adjoint equation ${\nabla }_{\mathrm{\Delta }X(t)}\mathcal{H}\left(\mathrm{\Delta }P,X,\mathrm{\Lambda },t\right)=-\frac{d}{dt}\mathrm{\Lambda }(t)$ which is equivalent to

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

The next condition is the state equation ${\nabla }_{\mathrm{\Lambda }(t)}\mathcal{H}\left(\mathrm{\Delta }P,X,\mathrm{\Lambda },t\right)=\frac{d}{dt}X (t)$ which is equivalent to (8). Finally, the last condition is given by the initial and terminal conditions $X\left(0\right)=X_0$ and $\mathrm{\Lambda }\left(T\right)=H_{T}X\left(T\right)$.

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

$$\frac{d}{dt}X\left(t\right)=AX\left(t\right)-r^{-1}B{B}^T\mathrm{\Lambda }(t)$$E19

Model Solution

We propose the following two equivalent solution approaches to solve the optimal control problem $(\mathcal{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 backward sweep method of Bryson and Ho [26], we let

$$\mathrm{\Lambda }\left(t\right)=S\left(t\right)X(t)$$E20

where $S(t)=\left[\begin{array}{cc}{s}_1(t)& 0\\ 0& s_2(t)\end{array}\right]$.

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

$$\frac{d}{dt}\mathrm{\Lambda }\left(t\right)=\left[\frac{d}{dt}S\left(t\right)+S\left(t\right)A-r^{-1}S\left(t\right)B{B}^{T}S\left(t\right)\right]X(t)$$E21

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

$$\frac{d}{dt}\mathrm{\Lambda }\left(t\right)=\left[-H-A^{T}S\left(t\right)\right]X(t)$$E22

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

$$\frac{d}{dt}S\left(t\right)=-H-S\left(t\right)A-A^{T}S\left(t\right)+r^{-1}S\left(t\right)B{B}^{T}S\left(t\right)$$E23

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)$$E24

The Riccati equation (23) becomes

$$\frac{d}{dt}E\left(t\right)F^{-1}\left(t\right)-E\left(t\right)F^{-1}\left(t\right)\frac{d}{dt}F\left(t\right)F^{-1}\left(t\right)=-H-E\left(t\right)F^{-1}\left(t\right)A-A^{T}E\left(t\right)F^{-1}\left(t\right)+r^{-1}E(t)F^{-1}(t)B{B}^{T}E(t)F^{-1}(t)$$E25

Multiplying this expression from the right by $F$ yields

$$\frac{d}{dt}E\left(t\right)-E\left(t\right)F^{-1}\left(t\right)\frac{d}{dt}F\left(t\right)=-HF-E\left(t\right)F^{-1}\left(t\right)AF-A^{T}E\left(t\right)+r^{-1}E(t)F^{-1}(t)B{B}^{T}E(t)$$E26

Now, set

$$\frac{d}{dt}E\left(t\right)=-HF(t)-A^{T}E(t)$$E27

Then, we have

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

Multiplying this expression from the left by ${\left(E{F}^{-1}\right)}^{-1}$ yields

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

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

$$\left[\begin{array}{c}\frac{d}{dt}E\left(t\right)\\ \cdots \cdots \cdots \\ \frac{d}{dt}F\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}{-A}^T& ⋮& -H\\ \cdots \cdots \cdots & ⋮& \cdots \cdots \cdots \\ {-r}^{-1}B{B}^T& ⋮& A\end{array}\right]\left[\begin{array}{c}E(t)\\ \cdots \cdots \cdots \\ F(t)\end{array}\right]$$E30

Call $G\left(t\right)=\left[\begin{array}{c}E(t)\\ F(t)\end{array}\right]$. The differential equations (30) become of the form

$$\frac{d}{dt}G\left(t\right)=\mathrm{\Gamma }G\left(t\right),G\left(t_0\right)\mathrm{ }\mathrm{g}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}, \mathrm{\Gamma }\mathrm{ }\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}$$E31

where

$\mathrm{\Gamma }=\left[\begin{array}{ccc}{-A}^T& ⋮& -H\\ \cdots \cdots \cdots & ⋮& \cdots \cdots \cdots \\ {-r}^{-1}B{B}^T& ⋮& A\end{array}\right]=\left[\begin{array}{cccc}\theta & -b& -h& 0\\ 1& -a& 0& -K\\ {-r}^{-1}& 0& -\theta & -1\\ 0& 0& b& a\end{array}\right]$

The boundary conditions $S(N)=H$ are equivalent to $E(N)F^{-1}(N)=H$ or $E(N)=H$ and $F\left(N\right)=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\left(t\right)=e^{\mathrm{\Gamma }(t-t_0)} G\left(t_0\right)$$E32

We recall that for a given a constant matrix $\mathrm{\Gamma }$, the matrix exponential $e^{\mathrm{\Gamma }t}$ is found as $e^{\mathrm{\Gamma }t}=P{e}^{Dt}{P}^{-1}$, where $D$ is the diagonal matrix whose elements are the eigenvalues of $\mathrm{\Gamma }$ and $P$ is the matrix whose columns are the corresponding eigenvectors. Thus,

$$G\left(t\right)=P{e}^{D(t-t_0)}{P}^{-1}G\left(t_0\right)$$E33

In the next solution approach, which also leads to a set of homogeneous equations, we will show how the constant term $G\left(t_0\right)$ 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 $(\mathcal{P})$ are:

$I\left(t\right)=\widehat{I}\left(t\right)+∆I\left(t\right)$, $D\left(t\right)=\widehat{D}\left(t\right)+∆D\left(t\right)$, $P\left(t\right)=\widehat{P}\left(t\right)+∆P\left(t\right)$,

where $∆I\left(t\right)$ and $∆I\left(t\right)$ are solutions of the linear equation:

$\frac{d}{dt}\left[\begin{array}{c}∆I(t)\\ ∆D(t)\end{array}\right]=\left(A-r^{-1}{B}^{T}E\left(t\right){F(t)}^{-1}\right)\left[\begin{array}{c}∆I(t)\\ ∆D(t)\end{array}\right]$,

with initial condition $X\left(0\right)=\left[\begin{array}{c}∆I(0)\\ ∆D(0)\end{array}\right]$

and

$∆P\left(t\right)=-r^{-1}{B}^{T}E\left(t\right){F\left(t\right)}^{-1}\left[\begin{array}{c}∆I(t)\\ ∆D(t)\end{array}\right].$

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

$$\left[\begin{array}{c}\frac{d}{dt}X\left(t\right)\\ \cdots \cdots \cdots \\ \frac{d}{dt}\mathrm{\Lambda }\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}A& ⋮& {-r}^{-1}B{B}^T\\ \cdots \cdots \cdots & ⋮& \cdots \cdots \cdots \\ -H& ⋮& -A^T\end{array}\right]\left[\begin{array}{c}X\left(t\right)\\ \\ \mathrm{\Lambda }\left(t\right)\end{array}\right]$$E34

Let $Z\left(t\right)=\left[\begin{array}{c}X\left(t\right)\\ \mathrm{\Lambda }\left(t\right)\end{array}\right]$. Then, the vector-matrix state equation (32) can be rewritten as

$$\frac{d}{dt}Z\left(t\right)=\mathrm{\Phi }Z(t)$$E35

where

$\mathrm{\Phi }=\left[\begin{array}{ccc}A& ⋮& -r^{-1}B{B}^T\\ \cdots \cdots \cdots & ⋮& \cdots \cdots \cdots \\ -H& ⋮& -A^T\end{array}\right]=\left[\begin{array}{cccc}-\theta & -1& -r^{-1}& 0\\ b& a& 0& 0\\ -h& 0& \theta & -b\\ 0& -K& 1& -a\end{array}\right]$

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

The matrix $\mathrm{\Phi }$ has four distinct eigenvalues $m_i, i=\mathrm{1,2},\mathrm{3,4}$. The explicit expressions of these eigenvalues are easily obtained using some mathematical software with symbolic computation capabilities such as MATHCAD, MAPLE, or MATLAB. These expressions are lengthy and thus are not reproduced here. The explicit expressions of the corresponding eigenvectors are also obtained using the same software. Then, the solution to (33) is given by

$$Z\left(t\right)=\mathrm{\phi }(t)Z\left(0\right)$$E36

Note that the first two components of $Z\left(t\right)$ thus computed form the state vector $X\left(t\right)$ whose components are $\mathrm{\Delta }I(t)$ and $\mathrm{\Delta }D(t)$, while the last two components form the co-state vector $\mathrm{\Lambda }(t)$ whose components are ${\lambda }_1\left(t\right)$ and ${\lambda }_2\left(t\right)$. In what follows, we show how $\mathrm{\phi }(t)$ and $Z\left(0\right)$ are determined using the two initial conditions $I\left(0\right)=I_0, D\left(0\right)=D_0$ and the terminal conditions ${\lambda }_1\left(T\right)=h_T\mathrm{\Delta }I\left(T\right), {\lambda }_2\left(T\right)=K_T\mathrm{\Delta }D(T)$.

To determine $\mathrm{\phi }(\mathrm{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,

$$\begin{array}{c}\mathrm{\phi }\left(t\right)=Y{e}^{Dt}{Y}^{-1}=\sum _{i=1}^{4}Y\left(:,i\right)Y^{-1}\left(i,:\right)e^{m_{i}t}\end{array}$$E37

where$Y\left(:,i\right)$ is the $i$th column of $Y$ and $Y^{-1}\left(i,:\right)$ is the $i$th row of $Y^{-1}$.

To determine $Z\left(0\right)=\left[\begin{array}{c}X\left(0\right)\\ \mathrm{\Lambda }\left(0\right)\end{array}\right]$, we recall that $X\left(0\right)=\left[\begin{array}{c}\mathrm{\Delta }I(0)\\ \mathrm{\Delta }D(0)\end{array}\right]$ is known while $\mathrm{\Lambda }\left(0\right)=\left[\begin{array}{c}{\mathrm{\lambda }}_1(0)\\ {\mathrm{\lambda }}_2(0)\end{array}\right]$ is not. However, using the final value $\mathrm{\Lambda }\left(T\right)$, we can find $\mathrm{\Lambda }\left(0\right)$ as follows. From (36), we have at $t=T$,

$Z\left(T\right)=\mathrm{\phi }(T)Z\left(0\right)$

which can be rewritten as

$\left[\begin{array}{c}X\left(T\right)\\ \mathrm{\Lambda }\left(T\right)\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\phi }}_1(T)& {\mathrm{\phi }}_2(T)\\ {\mathrm{\phi }}_3(T)& {\mathrm{\phi }}_4(T)\end{array}\right]\left[\begin{array}{c}X\left(0\right)\\ \mathrm{\Lambda }\left(0\right)\end{array}\right]$

Using the terminal condition

$\mathrm{\Lambda }\left(T\right)=H_{T}X\left(T\right)$

we have

$\left[\begin{array}{c}X\left(T\right)\\ H_{T}X\left(T\right)\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\phi }}_1(T)& {\mathrm{\phi }}_2(T)\\ {\mathrm{\phi }}_3(T)& {\mathrm{\phi }}_4(T)\end{array}\right]\left[\begin{array}{c}X\left(0\right)\\ \mathrm{\Lambda }\left(0\right)\end{array}\right]$

from which, it follows

$$\mathrm{\Lambda }\left(0\right)={\left(H_T{\mathrm{\phi }}_2\left(T\right)-{\mathrm{\phi }}_4\left(T\right)\right)}^{-1}({\mathrm{\phi }}_3\left(T\right)-H_T{\mathrm{\phi }}_1(T))X\left(0\right)$$E38

Finally, the optimal solutions of the problem $(\mathcal{P})$ using the second method are:

$I\left(t\right)=\widehat{I}\left(t\right)+∆I\left(t\right)$, $D\left(t\right)=\widehat{D}\left(t\right)+∆D\left(t\right)$, $P\left(t\right)=\widehat{P}\left(t\right)+∆P\left(t\right)$,

where $∆I\left(t\right)$ and $∆I\left(t\right)$ are solutions of the linear equation:

$\left[\begin{array}{c}∆I(t)\\ ∆D(t)\end{array}\right]={\phi }_1\left(t\right)\left[\begin{array}{c}∆I(0)\\ ∆D(0)\end{array}\right]+{\phi }_2\left(t\right)\mathrm{\Lambda }\left(0\right)$

$\mathrm{\Lambda }\left(t\right)={\phi }_3\left(t\right)\left[\begin{array}{c}∆I(0)\\ ∆D(0)\end{array}\right]+{\phi }_4\left(t\right)\mathrm{\Lambda }\left(0\right)$

and

$∆P\left(t\right)=-r^{-1}{B}^T\mathrm{\Lambda }\left(t\right) .$

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: $\mathrm{\Lambda }\left(k\right)=\left[\begin{array}{c}{\lambda }_I(k)\\ {\lambda }_D(k)\end{array}\right]$

Then, the Lagrangian function is given by

$$\mathcal{L}\left(P,Z,\mathrm{\Lambda }\right)=\sum _{k=0}^N\left\{\frac{1}{2}\left[{‖∆Z\left(k\right)‖}_Q^2+r {∆P\left(k\right)}^2\right]+{\mathrm{\Lambda }\left(k+1\right)}^T\left[-\mathrm{\Delta }Z\left(k+1\right)+A \mathrm{\Delta }Z\left(k\right)+B \mathrm{\Delta }P(k)\right]\right\}$$E39

The necessary optimality conditions are the control equation

$$\frac{\partial \mathcal{L}}{\partial ∆P(k)}=0 \iff r ∆P\left(k\right)+B^T\mathrm{\Lambda }\left(k+1\right)=0$$E40

and the adjoint equation

$$\frac{\partial \mathcal{L}}{\partial ∆Z(k)}=0 \iff Q∆Z\left(k\right)-\mathrm{\Lambda }\left(k\right)+A^T\mathrm{\Lambda }\left(k+1\right)=0.$$E41

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 backward sweep method. Let

$$\mathrm{\Lambda }\left(k\right)=S\left(k\right)\mathrm{ }\mathrm{\Delta }Z\left(k\right),$$E42

The control equation (40) becomes

$$\mathrm{\Delta }P\left(k\right)=-r^{-1}{B}^T\mathrm{\Lambda }\left(k+1\right)$$, $$= -r^{-1}{B}^{T}S\left(k+1\right)\mathrm{\Delta }Z(k+1)$$, $$= -r^{-1}{B}^{T}S\left(k+1\right)\left[\mathrm{A}\mathrm{ }\mathrm{\Delta }Z\left(k\right)+B \mathrm{\Delta }P(k)\right]$$. E43

Solving for $\mathrm{\Delta }P(k)$, we get

$$\mathrm{\Delta }P\left(k\right)=-V\left(k+1\right)\mathrm{ }\mathrm{\Delta }Z(k)$$E44

where

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

Now the adjoint equation (41) becomes

$$\mathrm{\Lambda }\left(k\right)=Q \mathrm{\Delta }Z\left(k\right)+A^T\mathrm{\Lambda }\left(k+1\right)$$E45

so that

$$S\left(k\right)\mathrm{ }\mathrm{\Delta }\mathrm{Z}\left(\mathrm{k}\right)=\mathrm{Q}\mathrm{ }\mathrm{\Delta }\mathrm{Z}\left(\mathrm{k}\right)+{\mathrm{A}}^{T}S\left(k+1\right)\mathrm{\Delta }Z(k+1)$$, $$=\mathrm{Q}\mathrm{ }\mathrm{\Delta }\mathrm{Z}\left(\mathrm{k}\right)+{\mathrm{A}}^{T}S\left(k+1\right)\left[A \mathrm{\Delta }Z\left(k\right)+B \mathrm{\Delta }P(k)\right]$$, $$=\mathrm{Q}\mathrm{ }\mathrm{\Delta }\mathrm{Z}\left(\mathrm{k}\right)+{\mathrm{A}}^{T}S\left(k+1\right)\left[A \mathrm{\Delta }Z\left(k\right)-B V\left(k+1\right)\mathrm{\Delta }Z\left(k\right)\right]$$, $$=\left\{\mathrm{Q}\mathrm{ }+{\mathrm{A}}^{T}S\left(k+1\right)A-A^{T}S\left(k+1\right)B V(k+1)\right\}\mathrm{\Delta }Z(k)$$. E46

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

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

The boundary condition $S(N)=Q$ follows from $\mathrm{\Delta }P\left(N\right)=0$. Now, to obtain the optimal production rates, we use the change of variable (42) to get from the adjoint equation (41),

$$AS\left(k+1\right)\mathrm{\Delta }Z\left(k+1\right)=S\left(k\right)\mathrm{\Delta }Z\left(k\right)-A\mathrm{\Delta }Z(k)$$E48

so that

$$\mathrm{\Delta }Z\left(k+1\right)={\left[AS(k+1)\right]}^{-1}\left[S\left(k\right)-Q\right]\mathrm{\Delta }Z(k)$$E49

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

$$P\left(k\right)=\widehat{P}\left(k\right)-r^{-1}{\left[1+r^{-1}{B}^T S\left(k+1\right)B\right]}^{-1}{B}^{T}S\left(k+1\right)A \mathrm{\Delta }Z(k)$$E50

where the optimal state vector $\mathrm{\Delta }Z(k)$ is found in expression (49) above and S(k) is found in expression (47).

Parameters		Values
Nonmonetary parameters	Length of planning horizon	$T=2.5$
	Initial inventory level	$I_0=5$
	Initial demand rate	$D_0=2$
	Demand rate coefficient	$a=0.8$
	Inventory level coefficient	$b=1$
	Deterioration rate	$\theta =0.1$
Monetary parameters	Penalty for production rate deviation	$r=0.1$
	Penalty for inventory level deviation	$h=4$
	Penalty for demand rate deviation	$K=10$
	Inventory salvage value	$h_T=100$
	Demand salvage value	$K_T=100$

Table 1.

Data for continuous-review model

4. Simulation results

4.1. Simulation of continuous review model

To illustrate numerically the results obtained, firstly we present some simulations for optimal control of the continuous-review integrated production-forecasting system with stock-dependent demand and deteriorating items. The data used in this simulation is presented in Table 1.

Using the MATLAB software, we implemented the results of the previous section and obtained the graphs below. Figures 1, 2 and 3 show the variations of $\mathrm{\Delta }I(t)$, $\mathrm{\Delta }D(t)$, and $\mathrm{\Delta }P(t)$. We observe that they all converge toward zero, as desired.

Using equation (3), the optimal cost is found to be J = 3047.93. A sensitivity analysis is performed in order to assess the effect of some of the system parameters on the optimal cost. The analysis is conducted by keeping the values of the parameters at the base values shown in Table 1 and varying successively one parameter at a time. We were interested in the effect on the value of the optimal objective function $J$ of the parameters $a$, $b$, and $\theta$, that we varied from 0.1 to 0.9. Table 2 summarizes the results of the sensitivity analysis.

Figure 1.

Variation of the optimal inventory level

Figure 2.

Variations of the optimal demand rate

Figure 3.

Variations of the optimal production rate

$a/b/\theta$	$J(a)$	$J(b)$	$J(\theta )$
0.1	5822.40	23426.21	3047.93
0.2	489.77	7648.73	463.79
0.3	485.62	5560.00	461.53
0.4	481.57	740.12	459.30
0.5	477.58	659.34	457.09
0.6	473.67	600.85	454.91
0.7	469.84	556.02	452.76
0.8	466.08	520.24	450.63
0.9	462.38	490.83	448.52

Table 2.

Sensitivity analysis

As can be seen, the objective function decreases as any of the three parameters increases.

4.2. Simulation of periodic review model

In this second part of the simulation, we illustrate the results obtained on the optimal control of the periodic review integrated production-forecasting system with stock-dependent demand and deteriorating items. Thus, consider the production planning problem for a firm for the next $T$ units of time. Divide this interval into $N$ subintervals of equal length. Assume the product in stock deteriorates at the rate $\theta$. Assume also the variations of the demand rate occur according to the dynamics (10). The firm has set the following targets. For $k=1,\cdots ,N$, the goal inventory level and goal demand rate are assumed to be as follows:

$\widehat{I}\left(k\right)=5+1.5 sign\left(sin\left(\frac{2\pi k}{40}\right)\right)$ and $\widehat{D}\left(k\right)=2+0.5 sign\left(sin\left(\frac{2\pi k}{15}\right)\right)$

where the sign function of a real number $x$ is defined by

$sign\left(x\right)=\left\{\begin{array}{cc}-1& \mathrm{i}\mathrm{f} x<0,\\ 0& \mathrm{i}\mathrm{f} x=0,\\ 1& \mathrm{i}\mathrm{f} x>0.\end{array}\right.$

We have to note that the goal inventory level and the goal demand rate were constant in the continuous review case.

The goal production rate is then computed using

$\widehat{P}\left(k\right)=\widehat{D}\left(k\right)+\theta \widehat{I}(k)$

where we assume that the inventory goal level is constant over a certain range. The penalties for deviating from these targets are $q_I$ for the inventory level, $q_D$ for the demand rate, and $r$ for the production rate. The data are summarized in Table 3.

Parameters	Values
planning horizon length $T$	10
number of subintervals $N$	51
coefficient for demand dynamics $a$	0.1
coefficient for demand dynamics $b$	0.2
deterioration rate $\theta$	0.01
deviation cost for inventory level $q_I$	20
deviation cost for demand rate $q_D$	15
deviation cost for production rate $r$	0.01

Table 3.

Data for periodic-review model

For the periodic review case, the simulation results are shown in the graphs below. Figure 4 shows the variations of the optimal inventory level and the inventory goal level. We observe that except for the early transient periods, $I(k)$ follows $\widehat{I}(k)$very closely.

Figure 4.

Optimal and inventory goal levels

Figure 5 shows the variations of the optimal demand rate and the demand goal rate. We observe that except for the early transient periods, $D(k)$ follows $\widehat{D}(k)$ very closely.

Finally, Figure 6 shows the variations of the optimal production rate and the production goal rate. We again observe that except for the early transient periods, $P(k)$ follows $\widehat{P}(k)$ very closely.

The optimal cost is found to be $J=10.3081$. Here also a sensitivity analysis is performed in order to assess the effect of some of the system parameters on the optimal cost. The analysis is conducted by keeping the values of the parameters at the base values shown in Table 3 and varying successively one parameter at a time. We were interested in the effect on $J$ of the parameters $a$, $b$, and $\theta$, that we varied from 0.1 to 0.9. Table 4 summarizes the results of the sensitivity analysis.

Figure 5.

Optimal and demand goal levels

Figure 6.

Optimal and production goal rates

$\text{a/b/θ}$	$\text{J(a)}$	$\text{J(b)}$	$\text{J(θ)}$
0.1	10.3081	10.0807	10.3072
0.2	10.3173	10.3081	10.3064
0.3	10.3339	10.6868	10.3056
0.4	10.3601	11.2166	10.3050
0.5	10.4004	11.8973	10.3044
0.6	10.4630	12.7284	10.3036
0.7	10.5651	13.7098	10.3040
0.8	10.7444	14.8413	10.3033
0.9	11.0880	16.1226	10.3032

Table 4.

Effect of the parameters $a$, $b$ and $\theta$ on the optimal cost $J$

The second column of Table 4 shows that the optimal cost increases as $a$ increases. The third column of Table 4 shows that the optimal cost increases also as $b$ increases. The effect of $b$ on $J$ is however more significant than the effect of $a$. Finally, column 4 of Table 4 shows that the optimal cost decreases as $\theta$ increases. The effect of $\theta$ is however almost negligible.

5. Conclusion

In this chapter we have used optimal control theory to derive the optimal production rate in a manufacturing 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 adopts either a continuous or periodic review policy. In contrast to most of the existing research which uses time series forecasting models, we propose a new model, namely, the demand dynamics equation. This model approaches realistic problems by integrating the forecasting component into the production planning problem with deteriorating items and stock dependent demand under continuous-review policy. Simulations were conducted in order to show the performance of the obtained solution. The theoretical and the simulations results allow gaining insights into operational issues and demonstrating the scope for improving stock control systems.

Of course, as with any research work, this study is not without limitations. The main contribution of our model is equation (5) where we use the demand from the previous period to predict the demand in the current period. The main limitation of that equation is that it involves two coefficients. We have assumed in this chapter that the parameters $a$ and $b$ of the demand state equation are known. However, in real life, that may not be the case. We are currently further investigating this model to estimate these parameters in the case when they are unknown, using self-tuning optimal control.

Another research direction would be to use a predictive control strategy where, given the current inventory level, the optimal production rates to be implemented at the beginning of each of the following periods over the control horizon, are determined. Model predictive control (or receding-horizon control) strategies have gained wide-spread acceptance in industry. It is also well-known that these models are interesting alternatives for real-time control of industrial processes. In the case where the above parameters $a$ and $b$ are unknown, the self-tuning predictive control can be applied. The proposed control algorithm estimates online these coefficients and feeds the controller to take the optimal production decision.

Note that our state equations are linear and thus linear model predictive control (LMPC), which is widely used both in academic and industrial fields, can be used. Nonlinear model predictive control (NMPC) can be used in case one of the state equations is nonlinear, for example, if equation (5) were of the form

$$\frac{d}{dt}D\left(t\right)=aD\left(t\right)+\alpha {I(t)}^{\beta }$$E51

NMPC has gained significant interest over the past decade. Various NMPC strategies that lead to stability of the closed-loop have been developed in recent years and key questions such as the efficient solution of the occurring open-loop control problem have been extensively studied.

The case combining unknown coefficients and a nonlinear relationship between the demand rate and the on-hand inventory yields a very complex, highly nonlinear process for which there is no simple mathematical model. The use of fuzzy control seems particularly well appropriate. Fuzzy control is a technique that should be seen as an extension to existing control methods and not their replacement. It provides an extra set of tools which the control engineer has to learn how to use where it makes sense. Nonlinear and partially known systems that pose problems to conventional control techniques can be tackled using fuzzy control.

Acknowledgement

This work has been supported by the Research Center of College of Computer and Information Sciences, King Saud University.