**Open Access Library Journal**

Vol.05 No.07(2018), Article ID:85987,15 pages

10.4236/oalib.1104716

Studying the Changes of an Optimal Trajectory

Ahmed Mohamed Mohamed El-Sayed^{ }

Al-Obour High Institute for Management & Informatics, Department of Basic Science, Obour, Egypt

Copyright © 2018 by author and Open Access Library Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: June 9, 2018; Accepted: July 13, 2018; Published: July 16, 2018

ABSTRACT

This paper provides an extension to an optimal control problem using the negative logarithm of deterioration and spoilage function as total cost. This function must be minimized at the end of planning period depending on the alternative quadratic exponential form. The co-state variable ${\lambda}_{0}\left(t\right)$ has negative values along the optimal trajectory according to the Pontryagin Minimum Principle (PMP). The different values of this co-state variable are investigated using initial values for the optimal control rates, separately. The controlled system according to each value is presented. Studying the behavior of optimal inventory levels, the optimal production rates, and the optimal spoilage function, it is our optimal solution along the optimal trajectory. The effectiveness of increasing and decreasing the co-state values on the optimal trajectory especially at the end of planning period is investigated. Also, the sensitivity analysis that reflects the effect of changes of different parameters (the deterioration and spoilage parameters values, and the initial values of inventory levels and production rates) on the optimal solution is explained with many different cases. Finally, we compared, numerically, the results for using these different co-state values with the results obtained when this value is negative.

**Subject Areas:**

Applied Statistical Mathematics

**Keywords:**

Sensitivity, Simulation, Co-State Variables, Multi-Item Inventory System, Demand Rates, Deterioration, Spoilage Function

1. Introduction

An optimal control problem of multi-item inventory model has a wide importance in practice. El-Sayed [1] has studied the effect of different types of demand rates on the total cost function, which refers to minimize the negative value of logarithm of deterioration and spoilage function at the end of planning period, using the Pontryagin Minimum Principle (PMP) for −1 value for the co-state variable ${\lambda}_{0}\left(t\right)$ . Now we need to know whether the changing of co-state values of this variable has effects on the optimal trajectory or not. So, in this paper, we will extend this study using different co-state values, whether less or more than negative one, for the co-state variable ${\lambda}_{0}\left(t\right)$ , using initial values for the optimal control rates. The purpose of this paper is obtaining the optimal inventory levels, the optimal production rates, and the minimum value of deterioration and spoilage function as the total cost at the end of planning period, and also studying the behavior of an optimal trajectory for each value, studying the sensitivity analysis of changing the parameters values of the model. As it is expected, the optimal inventory levels of the two items could be affected by these changes as we shall see later.

So, we concentrated on the sensitivity analysis of changing the spoilage parameters on the optimal solution especially at the end of planning period. Finally, to complete this study, we must compare the obtained results for these different co-state values with the results which are obtained when the co-state variable equals −1.

Zhao and Prentice [2] have presented the quadratic exponential form (QEF) for the two correlated variables ${X}_{1},{X}_{2}$ as:

$f\left({x}_{1},{x}_{2}\right)=\frac{1}{{\displaystyle \underset{{x}_{1},{x}_{2}}{\sum}\mathrm{exp}\left\{{\theta}_{1}{x}_{1}+{\theta}_{2}{x}_{2}+{\theta}_{12}{x}_{1}{x}_{2}\right\}}}\mathrm{exp}\left\{{\theta}_{1}{x}_{1}+{\theta}_{2}{x}_{2}+{\theta}_{12}{x}_{1}{x}_{2}\right\}$ . (1)

El-Sayed [1] has supposed that ${\theta}_{1}={\psi}_{1}{u}_{1}$ , ${\theta}_{2}={\psi}_{2}{u}_{2}$ and ${\theta}_{12}={\psi}_{12}{u}_{1}{u}_{2}$ , where the spoilage parameters ${\psi}_{1}$ , ${\psi}_{2}$ and ${\psi}_{12}$ depend on the control variables ${U}_{1}$ and ${U}_{2}$ , ${\psi}^{\prime}s>0$ .

Since ${\theta}_{1}$ , ${\theta}_{2}$ and ${\theta}_{12}$ are the deterioration parameters, ${\theta}^{\prime}s>0$ .

So, we can use the normalizing term, $\underset{{x}_{1},{x}_{2}}{\sum}\mathrm{exp}\left\{{\theta}_{1}{x}_{1}+{\theta}_{2}{x}_{2}+{\theta}_{12}{x}_{1}{x}_{2}\right\}$ , in the function (1) to be rewritten in the exponential form, El-Sayed et al. [3] , as shown below:

$f\left({x}_{1},{x}_{2}\right)=\mathrm{exp}\left\{{\theta}_{1}{x}_{1}+{\theta}_{2}{x}_{2}+{\theta}_{12}{x}_{1}{x}_{2}-\mathrm{log}\left(1+{\psi}_{1}{u}_{1}+{\psi}_{2}{u}_{2}+{\psi}_{12}{u}_{1}{u}_{2}\right)\right\}$ . (2)

The minimizing of an integral of negative logarithm of this function can be used as the total cost. This cost reflects the levels of deterioration and spoilage of items at the end of planning period.

This paper can be organized as follows: Section 2 presents the mathematical model for the optimal control problem and the coresponding controlled systems. Section 3 presents the numerical solution for the controlled systems for different rates and different co-state values. Section 4 presents the sensitivity analysis of the model parameters and co-state variable ${\lambda}_{0}\left(t\right)$ . Finally Section 5 gives some conclusions.

2. The Model

Let us define the following parameters, which are used in the mathematical optimal control model:

${X}_{i}\left(t\right)$ : Inventory levels at time t.

${U}_{i}\left(t\right)$ : Production rates at time t.

T: Length of planning period.

${x}_{i0}$ : Initial inventory levels.

${a}_{ii}$ : Deteriration coefficient due to self-contact of item ${x}_{i}$ .

${a}_{ij}$ : Deteriration coefficient of ${x}_{i}$ due to presence a unit of ${x}_{j}$ , $i\ne j=1,2$ .

${D}_{i}\left({x}_{1},{x}_{2},t\right)$ : Demand rates of $\left({x}_{1},{x}_{2}\right)$ .

${\psi}_{1}$ : Spoilage rate of of ${x}_{1}$ , ${\psi}_{1}>0$ .

${\psi}_{2}$ : Spoilage rate of of ${x}_{2}$ , ${\psi}_{2}>0$ .

${\psi}_{12}$ : Spoilage rate of $\left({x}_{1},{x}_{2}\right)$ , jointly, ${\psi}_{12}>0$ .

${\theta}_{1}$ : Natural deterioration rate of
${x}_{1}$ ^{,
${\theta}_{1}>0$ }.

${\theta}_{2}$ : Natural deterioration rate of ${x}_{2}$ , ${\theta}_{2}>0$ .

${\theta}_{12}$ : Natural deterioration rate of $\left({x}_{1},{x}_{2}\right)$ , jointly, ${\theta}_{12}>0$ .

As we mentioned before, the integral of negative logarithm of the function (2), which represents the deterioration and spoilage function, is used as a cost function:

$\begin{array}{c}{x}_{0}\left(T\right)={\displaystyle \underset{0}{\overset{T}{\int}}-\mathrm{ln}f\left({x}_{1},{x}_{2}\right)\text{d}t}\\ ={\displaystyle \underset{0}{\overset{T}{\int}}\left[-{\theta}_{1}{x}_{1}-{\theta}_{2}{x}_{2}-{\theta}_{12}{x}_{1}{x}_{2}+\mathrm{log}\left(1+{\psi}_{1}{u}_{1}+{\psi}_{2}{u}_{2}+{\psi}_{12}{u}_{1}{u}_{2}\right)\right]\text{d}t}.\end{array}$ (3)

So, the problem can be formulated as:

$\text{Minimize}\left\{{x}_{0}\left(T\right)={\displaystyle \underset{0}{\overset{T}{\int}}\left[-{\theta}_{1}{x}_{1}-{\theta}_{2}{x}_{2}-{\theta}_{12}{x}_{1}{x}_{2}+\mathrm{log}\left(1+{\psi}_{1}{u}_{1}+{\psi}_{2}{u}_{2}+{\psi}_{12}{u}_{1}{u}_{2}\right)\right]\text{d}t}\right\},$ (4)

subject to:

${\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{D}_{1}+{u}_{1},$ (5)

${\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{D}_{2}+{u}_{2},$ (6)

and

${x}_{1}\left(t\right)\ge \mathrm{0,}\text{\hspace{1em}}{x}_{2}\left(t\right)\ge \mathrm{0,}\text{\hspace{1em}}{u}_{1}\left(t\right)\ge \mathrm{0,}\text{\hspace{1em}}{u}_{2}\left(t\right)\ge 0,$ (7)

where,

$t\in T,\text{\hspace{0.17em}}D\left({x}_{1},{x}_{2},t\right)\ge \mathrm{0,}\text{\hspace{1em}}{\theta}_{1},{\theta}_{2},{\theta}_{12}>\mathrm{0,}\text{\hspace{1em}}{\psi}_{1},{\psi}_{2},{\psi}_{12}>0.$

Using the Pontryagin Minimum Principle (PMP), let us define ${\stackrel{\dot{}}{x}}_{0}=\frac{\partial {x}_{0}\left(T\right)}{\partial t}$ ,

and introduce the co-state variables ${\lambda}_{0}$ , ${\lambda}_{1}$ and ${\lambda}_{2}$ corresponding to the state variables ${X}_{0}$ , ${X}_{1}$ and ${X}_{2}$ respectively. From (4), (5) and (6), we can write the Hamiltonian function as follows:

$H={\lambda}_{0}{\stackrel{\dot{}}{x}}_{0}+{\lambda}_{1}{\stackrel{\dot{}}{x}}_{1}+{\lambda}_{2}{\stackrel{\dot{}}{x}}_{2},$ (8)

Moreover, to obtain the co-state equations and the Lagrange multipliers associated with the constraints (5) and (6), we formulate the Lagrangian function as follows:

$L=H+{\mu}_{1}{x}_{1}+{\mu}_{2}{x}_{2}+{\mu}_{3}{u}_{1}+{\mu}_{4}{u}_{2}$ , (9)

where, ${\mu}_{1}\left(t\right)\mathrm{,}{\mu}_{2}\left(t\right),{\mu}_{3}\left(t\right)\mathrm{,}{\mu}_{4}\left(t\right)$ are known as Lagrange multipliers. These Lagrange multipliers satisfy the conditions:

${\mu}_{1}\left(t\right)\ge 0,\text{\hspace{0.17em}}{\mu}_{2}\left(t\right)\ge 0,\text{\hspace{0.17em}}{\mu}_{3}\left(t\right)\ge 0,\text{\hspace{0.17em}}{\mu}_{4}\left(t\right)\ge 0,\text{\hspace{0.17em}}{\mu}_{1}{x}_{1}=0,\text{\hspace{0.17em}}{\mu}_{2}{x}_{2}=0,\text{\hspace{0.17em}}{\mu}_{3}{u}_{1}=0,\text{\hspace{0.17em}}{\mu}_{4}{u}_{2}=0.$ (10)

From (9), we can easily obtain the co-state equations

${\stackrel{\dot{}}{\lambda}}_{i}\left(t\right)=-\frac{\partial L}{\partial {x}_{i}},\text{\hspace{1em}}i=0,1,2,$ (11)

then,

${\stackrel{\dot{}}{\lambda}}_{0}\left(t\right)=-\frac{\partial L}{\partial {x}_{0}}=0,\text{\hspace{1em}}{\stackrel{\dot{}}{\lambda}}_{1}\left(t\right)=-\frac{\partial L}{\partial {x}_{1}},\text{\hspace{1em}}{\stackrel{\dot{}}{\lambda}}_{2}\left(t\right)=-\frac{\partial L}{\partial {x}_{2}},$ (12)

The first equation of (12) shows that the co-state variable ${\lambda}_{0}\left(t\right)$ remains constant along the optimal trajectory, and the Pontryagin principle requires that this constant should be a negative value, Sethi and Thompson [4] .

In this paper, we will use different values for this co-state variable ${\lambda}_{0}\left(t\right)$ .

${\lambda}_{0}=-10,{\lambda}_{0}=-2$ or ${\lambda}_{0}=-0.1,$ (13)

Substituting from (4), (5), (6), (8) and (13) in (9), we can write the Hamiltonian function, L, in the following form, first: when ${\lambda}_{0}\left(t\right)=-10$ :

$\begin{array}{c}L=10\left[{\theta}_{1}{x}_{1}+{\theta}_{2}{x}_{2}+{\theta}_{12}{x}_{1}{x}_{2}-\mathrm{log}\left(1+{\psi}_{1}{u}_{1}+{\psi}_{2}{u}_{2}+{\psi}_{12}{u}_{1}{u}_{2}\right)\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\lambda}_{1}\left[-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{D}_{1}+{u}_{1}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\lambda}_{2}\left[-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{D}_{2}+{u}_{2}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\mu}_{1}{x}_{1}+{\mu}_{2}{x}_{2}+{\mu}_{3}{u}_{1}+{\mu}_{4}{u}_{2}.\end{array}$ (14)

From conditions (7) and (10), we get

${\mu}_{1}\left(t\right)={\mu}_{2}\left(t\right)={\mu}_{3}\left(t\right)={\mu}_{4}\left(t\right)=0$ . (15)

Substituting from (13) and (14) into (12) we get

${\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{1}}+10{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{1}}+{a}_{21}{x}_{2}\right)-{\theta}_{1}-10{\theta}_{12}{x}_{2},$ (16)

${\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{2}}+10{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{2}}+{a}_{12}{x}_{1}\right)-{\theta}_{2}-10{\theta}_{12}{x}_{1},$ (17)

with boundary conditions

${\lambda}_{i}\left(T\right)\ne 0,\text{\hspace{1em}}i=1,2$ (18)

where T is the length of planning period which can be suggested.

To obtain the optimal production rates (control variables) ${U}_{i},i=1,2$ , we differentiate the Lagrange function (14) with respect to ${u}_{1},{u}_{2}$ respectively and putting it equal to zero, we get

$\frac{\partial L}{\partial {u}_{1}}=-\frac{10\left({\psi}_{1}+{\psi}_{12}{u}_{2}\right)}{1+{\psi}_{1}{u}_{1}+{\psi}_{2}{u}_{2}+{\psi}_{12}{u}_{1}{u}_{2}}+{\lambda}_{1}=0.$

$\frac{\partial L}{\partial {u}_{2}}=-\frac{10\left({\psi}_{2}+{\psi}_{1}{u}_{1}\right)}{1+{\psi}_{1}{u}_{1}+{\psi}_{2}{u}_{2}+{\psi}_{12}{u}_{1}{u}_{2}}+{\lambda}_{2}=0.$

Then,

${U}_{1}^{*}\left(t\right)=\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{2}}{10\left({\psi}_{1}+{\psi}_{12}{u}_{2}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{1}\ne 0$ (19)

${U}_{2}^{*}\left(t\right)=\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{1}}{10\left({\psi}_{2}+{\psi}_{12}{u}_{1}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{2}\ne 0$ (20)

Since, ${u}_{1}$ and ${u}_{2}$ are start values of the production rates. Then using the Equations (5), (6), (16) and (17), we get the controlled system of non-linear ordinary differential equations:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{D}_{1}+{u}_{1}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{D}_{2}+{u}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{1}}+10{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{1}}+{a}_{21}{x}_{2}\right)-{\theta}_{1}-10{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{2}}+10{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{2}}+{a}_{12}{x}_{1}\right)-{\theta}_{2}-10{\theta}_{12}{x}_{1}\end{array}\}$ (21)

We can construct this system when ${\lambda}_{0}=-2$ :

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{D}_{1}+{u}_{1}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{D}_{2}+{u}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{1}}+2{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{1}}+{a}_{21}{x}_{2}\right)-{\theta}_{1}-2{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{2}}+2{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{2}}+{a}_{12}{x}_{1}\right)-{\theta}_{2}-2{\theta}_{12}{x}_{1}\end{array}\}$ (22)

and when ${\lambda}_{0}=-0.1$

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{D}_{1}+{u}_{1}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{D}_{2}+{u}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{1}}+0.1{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{1}}+{a}_{21}{x}_{2}\right)-{\theta}_{1}-0.1{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(\frac{\partial D{}_{2}}{\partial {x}_{2}}+0.1{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{\lambda}_{1}\left(\frac{\partial D{}_{1}}{\partial {x}_{2}}+{a}_{12}{x}_{1}\right)-{\theta}_{2}-0.1{\theta}_{12}{x}_{1}\end{array}\}$ (23)

The optimal control variables can be constructed when ${\lambda}_{0}=-2$ :

${U}_{1}^{*}\left(t\right)=\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{2}}{2\left({\psi}_{1}+{\psi}_{12}{u}_{2}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{1}\ne 0,$ (24)

${U}_{2}^{*}\left(t\right)=\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{1}}{2\left({\psi}_{2}+{\psi}_{12}{u}_{1}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{2}\ne 0,$ (25)

and when ${\lambda}_{0}\left(t\right)=-0.1$ :

${U}_{1}^{*}\left(t\right)=\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{2}}{0.1\left({\psi}_{1}+{\psi}_{12}{u}_{2}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{1}\ne 0,$ (26)

${U}_{2}^{*}\left(t\right)=\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{1}}{0.1\left({\psi}_{2}+{\psi}_{12}{u}_{1}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{2}\ne 0.$ (27)

This system can be used to describe the time evolution of inventory levels and production rates. The analytical solution of this system is very difficult and then we can solve it numerically.

3. Numerical Solution

The solution of optimal control problem of this model will be carried out using Pontryagin Minimum Principle (PMP). The numerical solution is to be necessary when the analytical solution is absence for the non-linear systems (21, 22 and 23). In this solution we solve the non-linear ordinary differential equations using Runge-Kutta method, using the initial and boundary values for ${x}_{1}\left(t\right),{x}_{2}\left(t\right),{\lambda}_{1}\left(t\right)$ and ${\lambda}_{2}\left(t\right)$ . For simplicity we supposed that the initial values ${u}_{10},{u}_{20}$ can be used alternative to ${u}_{1},{u}_{2}$ in the Equations (19), (20), (24), (25), (26) and (27) to obtain the optimal production rates ${U}_{1}^{*}\left(T\right)$ and ${U}_{2}^{*}\left(T\right)$ respectively. Also, we will use these initial values to obtain the optimal total cost ${x}_{0}^{*}\left(T\right)$ as it is in the Equation (4).

The numerical solution can be explained by different types of demand as:

1) The demand rates are constant:

$D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}.$

2) The demand rates are linear functions of inventory levels and time:

$D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}+{w}_{i}{x}_{i}.$

3) The demand rates are logistic functions of inventory levels and time:

$D\left({x}_{1},{x}_{2},t\right)=2{x}_{i}\left({\kappa}_{i}-{x}_{i}\right).$

4) The demand rates are periodic functions of time:

$D\left({x}_{1},{x}_{2},t\right)=1-{b}_{i}\mathrm{cos}\left(t\right).$

where ${\gamma}_{i}\mathrm{,}{w}_{i}\mathrm{,}{\kappa}_{i}$ and ${b}_{i}\left(i=1,2\right)$ are positive constants.

Table 1 presents the values of system parameters and the initial states which are used in the numerical examples for four cases of demand rate functions as follows:

The next subsections explain the controlled system for each case of the demand rates functions with different co-state value ${\lambda}_{0}\left(t\right)$ as shown below.

3.1. Co-State Value ${\lambda}_{0}=-10$

In this subsection we will use different demand rates with ${\lambda}_{0}=-10$ .

3.1.1. Constant Rates

We will present the model with demand function as constant rates, $D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}$ .

Substituting in the controlled system (21) by the constant demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{\gamma}_{1}+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{10\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{\gamma}_{2}+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{10\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(10{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-10{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(10{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-10{\theta}_{12}{x}_{1}\end{array}\}$ (28)

Solving the controlled system (28) numerically, we get some results as displayed in Table 2.

3.1.2. Linear Rates

Also, we will present the model with demand function as linear rates, $D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}+{w}_{i}{x}_{i}$ . Substituting in the controlled system (21) by the linear demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({w}_{1}+{\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{\gamma}_{1}+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{10\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({w}_{2}+{\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{\gamma}_{2}+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{10\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left({w}_{1}+10{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-10{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left({w}_{2}+10{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-10{\theta}_{12}{x}_{1}\end{array}\}$ (29)

Solving the controlled system (29) numerically, we get some results as displayed in Table 2.

3.1.3. Logistic Rates

We present the model with demand function as logistic rates, $D\left({x}_{1},{x}_{2},t\right)=2{x}_{i}\left({\kappa}_{i}-{x}_{i}\right)$ .

Table 1. Values and initial states of system parameters.

Table 2. The optimal solution when ${\lambda}_{0}=-10$ .

Substituting in the controlled system (21) by the logistic demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left(2\left({\kappa}_{1}-{x}_{1}\right)+{\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{10\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}(2({\kappa}_{2}-{x}_{2})+{\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2})+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{10\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(2\left({\kappa}_{1}-2{x}_{1}+{a}_{11}{x}_{1}\right)+10{\theta}_{1}+{a}_{12}{x}_{2}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-10{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(2\left({\kappa}_{2}-2{x}_{2}+{a}_{22}{x}_{2}\right)+10{\theta}_{2}+{a}_{21}{x}_{1}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-10{\theta}_{12}{x}_{1}\end{array}\}$ (30)

Solving the controlled system (30) numerically, we get some results as displayed in Table 2.

3.1.4. Periodic Rates

Finally, we will present the model with demand function as periodic rates, $D\left({x}_{1},{x}_{2},t\right)=1-{b}_{i}\mathrm{cos}\left(t\right)$ . Substituting in the controlled system (21) by the periodic demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-1+{b}_{1}\mathrm{cos}\left(t\right)+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{10\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-1+{b}_{2}\mathrm{cos}\left(t\right)+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{2}{u}_{10}}{10\left({\psi}_{1}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(10{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-10{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(10{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-10{\theta}_{12}{x}_{1}\end{array}\}$ (31)

Solving the controlled system (31) numerically, we get some results as displayed in Table 2.

As we see from Table 2, the otimal value of the cost function when ${\lambda}_{0}=-10$ is achieved in the logistic rate (0). The optimal production rates are simillar (0.87, 0.86) in all cases. The inventory levels are smaller (0.81, 0.87) in the linear case but higher (1.16, 1.14) in the periodic case.

3.2. Co-State ${\lambda}_{0}=-2$

Also, in this subsection we use the previous demand rates with ${\lambda}_{0}=-2$ .

3.2.1. Constant Rates

We will present the model with demand function as constant rates, $D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}$ . Substituting in the controlled system (22) by the constant demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}\mathrm{=}-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{\gamma}_{1}+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{2\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}\mathrm{=}-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{\gamma}_{2}+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{2\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}\mathrm{=}{\lambda}_{1}\left(2{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-2{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}\mathrm{=}{\lambda}_{2}\left(2{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-2{\theta}_{12}{x}_{1}\end{array}\}$ (32)

Solving the controlled system (32) numerically, we get some results which are displayed in Table 3.

3.2.2. Linear Rates

Also, we will present the model with demand function as linear rates, $D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}+{w}_{i}{x}_{i}$ . Substituting in the controlled system (22) by the linear demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({w}_{1}+{\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{\gamma}_{1}+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{2\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({w}_{2}+{\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{\gamma}_{2}+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{2\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left({w}_{1}+2{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-2{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left({w}_{2}+2{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-2{\theta}_{12}{x}_{1}\end{array}\}$ (33)

Solving the controlled system (33) numerically, we get the results displayed in Table 3.

3.2.3. Logistic Rates

We present the model with demand function as logistic rates, $D\left({x}_{1},{x}_{2},t\right)=2{x}_{i}\left({\kappa}_{i}-{x}_{i}\right)$ . Substituting in the controlled system (22) by the logistic demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left(2\left({\kappa}_{1}-{x}_{1}\right)+{\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{2\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left(2\left({\kappa}_{2}-{x}_{2}\right)+{\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{2\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(2\left({\kappa}_{1}-2{x}_{1}+{a}_{11}{x}_{1}\right)+2{\theta}_{1}+{a}_{12}{x}_{2}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-2{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(2\left({\kappa}_{2}-2{x}_{2}+{a}_{22}{x}_{2}\right)+2{\theta}_{2}+{a}_{21}{x}_{1}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-2{\theta}_{12}{x}_{1}\end{array}\}$ (34)

Solving the controlled system (30) numerically, we get the results are displayed in Table 3.

3.2.4. Periodic Rates

Finally, we will present the model with demand function as periodic rates, $D\left({x}_{1},{x}_{2},t\right)=1-{b}_{i}\mathrm{cos}\left(t\right)$ . Substituting in the controlled system (22) by the periodic demand rates, we have the controlled system:

Table 3. The optimal solution when ${\lambda}_{0}=-2$ .

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-1+{b}_{1}\mathrm{cos}\left(t\right)+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{2}{u}_{20}}{2\left({\psi}_{1}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-1+{b}_{2}\mathrm{cos}\left(t\right)+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{2\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(2{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-2{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(2{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-2{\theta}_{12}{x}_{1}\end{array}\}$ (35)

Solving the controlled system (35) numerically, we get the results can be displayed in Table 3.

We will use the parameters values as they are in Table 1 without changing in all cases.

As we see from Table 3, the otimal value of the cost function when ${\lambda}_{0}=-2$ is equal (0) in all cases. Also, the optimal production rates are simillar (0.36, 0.31) in all cases. The optimal inventory levels are smaller (0.65, 0.68) in the logistic case but higher (1.57, 1.69) in the periodic case.

3.3. Co-State Value ${\lambda}_{0}=-0.1$

Finally, we use the previous demand rates with ${\lambda}_{0}=-0.1$ .

3.3.1. Constant Rates

We will present the model with demand function as constant rates, $D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}$ .

Substituting in the controlled system (23) by the constant demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{\gamma}_{1}+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{1}{u}_{20}}{0.1\left({\psi}_{2}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{\gamma}_{2}+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{0.1\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left(0.1{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-0.1{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left(0.1{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-0.1{\theta}_{12}{x}_{1}\end{array}\}$ (36)

Solving the controlled system (36) numerically, we get the results displayed in Table 4.

3.3.2. Linear Rates

Also, we will present the model with demand function as linear rates, $D\left({x}_{1},{x}_{2},t\right)={\gamma}_{i}+{w}_{i}{x}_{i}$ . Substituting in the controlled system (23) by the linear demand rates, we have the controlled system:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}=-{x}_{1}\left({w}_{1}+{\theta}_{1}+{a}_{12}{x}_{2}+{a}_{11}{x}_{1}\right)-{\gamma}_{1}+\frac{1}{{\lambda}_{1}}-\frac{1+{\psi}_{1}{u}_{20}}{0.1\left({\psi}_{2}+{\psi}_{12}{u}_{20}\right)}\\ {\stackrel{\dot{}}{x}}_{2}=-{x}_{2}\left({w}_{2}+{\theta}_{2}+{a}_{21}{x}_{2}+{a}_{22}{x}_{2}\right)-{\gamma}_{2}+\frac{1}{{\lambda}_{2}}-\frac{1+{\psi}_{1}{u}_{10}}{0.1\left({\psi}_{2}+{\psi}_{12}{u}_{10}\right)}\\ {\stackrel{\dot{}}{\lambda}}_{1}={\lambda}_{1}\left({w}_{1}+0.1{\theta}_{1}+{a}_{12}{x}_{2}+2{a}_{11}{x}_{1}\right)+{a}_{21}{\lambda}_{2}{x}_{2}-{\theta}_{1}-0.1{\theta}_{12}{x}_{2}\\ {\stackrel{\dot{}}{\lambda}}_{2}={\lambda}_{2}\left({w}_{2}+0.1{\theta}_{2}+{a}_{21}{x}_{1}+2{a}_{22}{x}_{2}\right)+{a}_{12}{\lambda}_{1}{x}_{1}-{\theta}_{2}-0.1{\theta}_{12}{x}_{1}\end{array}\}$ (37)

Table 4. The optimal solution when ${\lambda}_{0}=-0.1$ .

Solving the controlled system (37) numerically, we get the results displayed in Table 4.

3.3.3. Logistic Rates

We present the model with demand function as logistic rates, .

Substituting in the controlled system (23) by the logistic demand rates, we have the controlled system:

(38)

Solving the controlled system (38) numerically, we get the results displayed in Table 4.

3.3.4. Periodic Rates

Finally, we will present the model with demand function as periodic rates,. Substituting in the controlled system (23) by the periodic demand rates, we have the controlled system:

(39)

Solving the controlled system (39) numerically, we get the results displayed in Table 4.

We will use the parameters values as they are in Table 1. But we changed the next values, since the models are more sensititvity for changing these paprameters:

Constant case

Linear case

Logistic case

Periodic case

As we see from Table 4, the minimum cost when is achieved in the constant case (0.67). Also, the optimal production rates are simillar (0, 0) in all cases. The optimal inventory levels are higher in the constant case but smaller in the logistic and periodic cases (0.13, 0.21).

4. Sensitivity Analysis

In this section will compare between the results when the co-state value, using the parameters values are used in Table 1. Then we can study the sensitivity analysis for cthese changes in all cases.

The optimal production rates, when, is become:

(40)

(41)

Also, we can use the initial values and alternatively to and respectively.

The controlled systems are become in each case as follow.

4.1. Constant Rates

The controlled system is become:

(42)

4.2. Linear Rates

The controlled system is become:

(43)

4.3. Logistic Rates

The controlled system is become:

(44)

4.4. Periodic Rates

The controlled system is become:

(45)

Constant case

Linear case

Logistic case

Periodic case

The optimal cost is become as in the function (4) after replacing, in the Equations (40) and (41) by and respectively. The optimal solution is presented when in Table 5 as shown below.

As we see from Table 5, the minimum cost when is achieved in the logistic case (0.02). Unlike all previous cases for the values of, the optimal production rates are not simillar. These rates are higher (0.85, 0.76) in the linear case but smaller (0.48, 0.26) in the constant case . The optimal inventory levels are higher (1.86, 1.95) in the periodic case but smaller (0.89, 0.90) in the logistic case.

So, we can conclude that the co-state value is more sensitive for changing the values of parameters and actually more effectiveness on the optimal solution.

5. Conclusions

In this study, we discussed the optimal control problem using the deterioration and spoilage function, as the total cost must be minimized at the end of planning period, depending on the alternative quadratic exponential form (AQEF). We have used different values for the co-state variable which has negative values along the optimal trajectory. Also, we studied the effectiveness of increasing and decreasing these values on the behavior of optimal trajectory and then the optimal solution (the inventory levels, the production rates and the total cost) at the end of planning period. Also, we explained the sensitivity analysis for the effect of changing the values of model parameters (especially the deterioration and spoilage parameters and the initial values for inventory levels and production rates) on the optimal solution when the variable. Also, we compared the results that obtained when the variable, with the results that obtained when and have conducted that the model is more sensitive for changing the values of. Finally, we concluded that the production rates are similar in all cases for the

Table 5. The optimal solution when.

co-state variable except when. Also, the optimal inventory levels of the two items are affected by these changes in all cases.

Cite this paper

El-Sayed, A.M.M. (2018) Studying the Changes of an Optimal Trajectory. Open Access Library Journal, 5: e4716. https://doi.org/10.4236/oalib.1104716

References

- 1. El-Sayed, A. (2017) Alternative Quadratic Exponential Form to Optimal Control Problem. Journal of Statistical Science, ISSR, Cairo University, 61, 1.
- 2. Zhao, L.P. and Prentice, R. L. (1990) Correlated Binary Regression Using a Quadratic Exponential Model. Biometrika, 77, 642-648.

https://doi.org/10.1093/biomet/77.3.642 - 3. El-Sayed, A., Islam, M. and Alzaid, A. (2013) Estimation and Test of Measures of Association for Correlated Binary Data. Bulletin of the Malaysian Mathematical Sciences Society, 36, 985-1008.
- 4. Sethi, S.P. and Thompson, G.L. (2000) Optimal Control Theory: Applications to Management Science and Economics. Kluwer Academic Publishers, Boston, Dordrecht, London.