Table of Contents Author Guidelines Submit a Manuscript
Computational and Mathematical Methods in Medicine
VolumeΒ 2011, Article IDΒ 163834, 9 pages
Research Article

Optimal Campaign in the Smoking Dynamics

Centre for Advanced Mathematics and Physics, National University of Sciences and Technology, H-12, Islamabad, 64000, Pakistan

Received 16 June 2010; Accepted 16 December 2010

Academic Editor: SivabalΒ Sivaloganathan

Copyright Β© 2011 Gul Zaman. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


We present the optimal campaigns in the smoking dynamics. Assuming that the giving up smoking model is described by the simplified 𝑃𝐿𝑆𝑄 (potential-light-smoker-quit smoker) model, we consider two possible control variables in the form of education and treatment campaigns oriented to decrease the attitude towards smoking. In order to do this we minimize the number of light (occasional) and persistent smokers and maximize the number of quit smokers in a community. We first show the existence of an optimal control for the control problem and then derive the optimality system by using the Pontryagin maximum principle. Finally numerical results of real epidemic are presented to show the applicability and efficiency of this approach.

1. Introduction

Application of control theory to epidemics is a very large field, and the study of epidemic models is strongly related to the possibility of evaluation of different control strategies: screening and educational campaigns [1], the vaccination campaign [2], and resource allocation [3]. A comprehensive survey of control theory applied to epidemiology was performed by Wickwire [4]. Many different models with different objective functions have been proposed (see [5–9]). A major difficulty in applying control theoretic methods to practical epidemiology problems is the commonly made assumption that one has total knowledge of the state of the epidemics [10].

In 2000, Castillo-Garsow et al. [11] for the first time proposed a simple mathematical model for giving up smoking. They consider a system with a total constant population which is divided into three classes: potential smokers, that is, people who do not smoke yet but might become smokers in the future (𝑃), smokers (𝑆), and people (former smokers) who have quit smoking permanently (𝑄). Sharomi and Gumel developed mathematical models by introducing mild and chain classes [12]. In their work they presented the development and public health impact of smoking-related illnesses. Zaman [13] extended the work of Castillo-Garsow et al. [11] and developed a model taking into account the occasional smokers compartment in the giving up smoking model and presented its qualitative behavior.

In this study we consider the model introduced by Zaman [13] and extend that once a smoker quits smoking he/she may become a potential smoker again and first discuss the dynamical behavior and then use optimal control theory for our optimal problem to minimize the number of light and persistent smokers and maximize the number of quit smokers in a community. In our control problem we first show the existence of an optimal control and then we derive the optimality system. The technical tool used in this work to determine the optimal strategy is the Pontryagin maximum principle. We emphasize that we have not set up a smoking model that particularly fit our optimization scheme. In order to do this we use optimal control strategies associated with two types of control, an education and a treatment campaigns to minimize the number of light and persistent smokers and maximize the number of quit smokers in a community. That is, we derive the optimality system consisting of the state and adjoint equations and solve numerically the system by using an iterative method. We also give an example of the real epidemic model by real data so that we illustrate how the optimal control theory can be applied in real situation.

The structure of this paper is organized as follows. First, we presented the qualitative behavior of the proposed model and then control system for the optimality and its existence, and the optimal control pairs are derived in Section 2. In Section 3, we solve numerically the optimality system by using real data, which consists of the original state system, the adjoint system, and their boundary conditions. Finally, we conclude by discussing results of the numerical simulation for our giving up smoking model.

2. Optimal Control Problem

In this section, we begin our explorations with a 𝑃𝐿𝑆𝑄 model for a population as defined by the following system consisting of differential equations:𝑑𝑃(𝑑)𝑑𝑑=𝑏𝑁(𝑑)βˆ’π›½1𝑑𝐿(𝑑)𝑃(𝑑)βˆ’1ξ€Έ+πœ‡π‘ƒ(𝑑)+πœŒπ‘„(𝑑),𝑑𝐿(𝑑)𝑑𝑑=𝛽1𝐿(𝑑)𝑃(𝑑)βˆ’π›½2𝑑𝐿(𝑑)𝑆(𝑑)βˆ’2ξ€Έ+πœ‡πΏ(𝑑),𝑑𝑆(𝑑)𝑑𝑑=𝛽2𝐿(𝑑)𝑆(𝑑)βˆ’π›Ύ+𝑑3ξ€Έ+πœ‡π‘†(𝑑),𝑑𝑄(𝑑)𝑑𝑑𝑑=𝛾𝑆(𝑑)βˆ’4ξ€Έ+πœ‡+πœŒπ‘„(𝑑).(1) As described in [13] the variables 𝑃, 𝐿, 𝑆, and 𝑄 represent the potential, light (occasional), persistent, and quit smoker, respectively. The parameter 𝑏 represents the birth rate, πœ‡ is the natural death rate, and 𝛾 is the quit rate from smoking. The parameters 𝛽1 and 𝛽2 approximate the average number of contacts with smoker individuals needed to make light and persistent smoker, respectively, in each unit of time. 𝜌 is the rate of quit smoker becoming potential smoker again, and 𝑑1,  𝑑2,  𝑑3, and 𝑑4 represent the death rate of potential smoker, occasional smoker, persistent smoker, and quit smoker, respectively. In this model, we also use the classic mass action hypothesis for both positive transmission coefficients 𝛽1 and 𝛽2. Potential smokers acquire infection at per capita rate 𝛽1𝐿(𝑑). The total population size at time 𝑑 is denoted by 𝑁(𝑑) with 𝑁(𝑑)=𝑃(𝑑)+𝐿(𝑑)+𝑆(𝑑)+𝑄(𝑑). The dynamics of the total population are governed by the following differential equation:𝑑𝑁𝑑𝑑𝑑=(π‘βˆ’πœ‡)𝑁(𝑑)βˆ’1𝑃(𝑑)+𝑑2𝐿(𝑑)+𝑑3𝑆(𝑑)+𝑑4ξ€Έ.𝑄(𝑑)(2) If all individuals die at the same death rates, then we have 𝑑1=𝑑2=𝑑3=𝑑4=𝑑, and dynamics of the total population become𝑑𝑁𝑑𝑑=(π‘βˆ’π‘‘βˆ’πœ‡)𝑁(𝑑).(3) When the birth rate is equal to the disease and natural death rate then the total population remains constant, that is, 𝑏=𝑑+πœ‡. The total constant population models were studied by several authors; see, for example, [8, 14]. In our proposed model it is not biologically feasible if we consider 𝑏<𝑑+πœ‡, so we assume that 𝑏⩾𝑑+πœ‡. The unique positive epidemic equilibrium of the system (1) is given by𝑃⋆=𝛽2𝑆⋆+𝑑2ξ€Έ+πœ‡π›½1,𝐿⋆=𝛾+𝑑3+πœ‡π›½2,𝑆⋆=𝑑2𝑑+πœ‡ξ€Έξ€·4𝛽+πœ‡+πœŒξ€Έξ€·1𝛾+𝑑3ξ€Έ+πœ‡+𝛽2𝑑1+πœ‡ξ€Έξ€Έπ›½2𝛽1πœŒπ›Ύβˆ’π›½1𝛾+𝑑3ξ€Έ+πœŒβˆ’π›½2𝑑1𝑑+πœ‡ξ€Έξ€·4𝑄+πœ‡+πœŒξ€Έξ€ΈΓ—(πΎβˆ’1),⋆=𝑑2𝛽+πœ‡ξ€Έξ€·1𝛾+𝑑3ξ€Έ+πœ‡+𝛽2𝑑1+πœ‡ξ€Έξ€Έπ›½2𝛽1πœŒπ›Ύβˆ’π›½1𝛾+𝑑3ξ€Έ+πœŒβˆ’π›½2𝑑1𝑑+πœ‡ξ€Έξ€·4+πœ‡+πœŒξ€Έξ€ΈΓ—(πΎβˆ’1),(4) where 𝐾=𝑏𝛽1𝛽2𝑁⋆/((𝑑2+πœ‡)(𝛽1(𝛾+𝑑3+𝜌)+𝛽2(𝑑1+πœ‡))) represents the smoking generation number (basic reproductive number). It measures the average number of new smokers generated by single smoker in a population of potential smokers. The reproduction number is a concept in the epidemiology of infectious diseases [12]. It is a measure of how infectious a disease is and is required if you wish to calculate how many people you need to vaccinate if you are to achieve herd immunity. When somebody starts smoking, they may pass it on to nobody else or they may assist 1, 2, or more other people (who become secondary cases). The reproduction number, 𝐾, is the average (mean) number of secondary cases caused by each case of an infectious disease, during the infectious period. The number 𝐾 will, of course, depend on a large number of factors depending on the states of the diseases. In our proposed model when 𝐾>1 there exist persistent smokers in the community while for 𝐾<1 there are no smokers.

In order to understand the qualitative behavior of the giving up smoking model we find the different equilibrium position of the model.

Remark 1. The Jacobian matrix around the trivial equilibrium 𝐸𝑑=(0,0,0,0) is 𝐽𝑑=βŽ‘βŽ’βŽ’βŽ’βŽ£βˆ’π‘‘1βˆ’πœ‡00𝜌0βˆ’πœ‡βˆ’π‘‘20000βˆ’π‘‘3βˆ’π›Ύβˆ’πœ‡000π›Ύβˆ’π‘‘4⎀βŽ₯βŽ₯βŽ₯βŽ¦βˆ’πœ‡βˆ’πœŒ.(5) The eigenvalues of the Jacobian matrix around the trivial equilibrium 𝐸0=(0,0,0,0) are βˆ’π‘‘1βˆ’πœ‡, βˆ’π‘‘2βˆ’πœ‡, βˆ’π›Ύβˆ’π‘‘3βˆ’πœ‡, βˆ’π‘‘4βˆ’πœ‡βˆ’πœŒ. Thus all the roots have negative real part, which shows that the trivial equilibrium is locally stable.

Theorem 1. When 𝐾>1, the giving up smoking model (1) has unique positive epidemic equilibrium 𝐸⋆=(𝑃⋆,𝐿⋆,𝑆⋆,𝑄⋆), where 𝑃⋆,𝐿⋆,𝑆⋆, and 𝑄⋆ are defined in (4).

Theorem 2. The giving up smoking model (1) has 𝐸𝑙=(1,0,0,0) as a locally stable smoking-free equilibrium if and only if 𝛽1<𝑑2+πœ‡. Otherwise 𝐸𝑙 is an unstable smoking-free equilibrium.

Proof. We can examine by linearizing our proposed giving up smoking model (1) around 𝐸𝑙=(1,0,0,0) to get local stability. Hence, the equilibrium point around 𝐸𝑙 gives us the Jacobian matrix: 𝐽𝑙=βŽ‘βŽ’βŽ’βŽ’βŽ£βˆ’π‘‘1βˆ’πœ‡βˆ’π›½10𝜌0𝛽1βˆ’πœ‡βˆ’π‘‘20000βˆ’π‘‘3βˆ’π›Ύβˆ’πœ‡000π›Ύβˆ’π‘‘4⎀βŽ₯βŽ₯βŽ₯βŽ¦βˆ’πœ‡βˆ’πœŒ.(6) The eigenvalues of the Jacobian matrix 𝐽𝑙 around smoking-free equilibrium 𝐸𝑙=(1,0,0,0) are βˆ’π‘‘1βˆ’πœ‡, 𝛽1βˆ’π‘‘2βˆ’πœ‡, βˆ’π›Ύβˆ’π‘‘3βˆ’πœ‡, βˆ’π‘‘4βˆ’πœ‡βˆ’πœŒ. Thus, we deduce that all the roots have negative real part when 𝛽1<𝑑2+πœ‡ which shows that the smoking-free equilibrium is locally asymptotically stable.

Theorem 3. When 𝐾>1, the unique positive epidemic equilibrium 𝐸⋆=(𝑃⋆,𝐿⋆,𝑆⋆,𝑄⋆) of the giving up smoking model (1) is locally asymptotically stable.

Proof. For the proof of this theorem see Appendix A

In order to investigate an effective campaign to control smoking in a community which satisfies that the maximum number of occasional (light) and persistent smoker individuals is not larger than that of potential smoker individuals and more individuals are recovered after quitting the use of tobacco, we consider the system (1) and use two control variables to control both the occasional (light) smoker and persistent smoker population. In the system (1) we have four state variables 𝑃(𝑑), 𝐿(𝑑), 𝑆(𝑑), and 𝑄(𝑑). For the optimal control problem we consider the control variable 𝑒(𝑑)=(𝑒1(𝑑),𝑒2(𝑑))βˆˆπ‘ˆ relative to the state variables (𝑃(𝑑),𝐿(𝑑),𝑆(𝑑),𝑄(𝑑)), where π‘ˆ={(𝑒1,𝑒2)βˆ£π‘’π‘–(𝑑) is measurable and 0≀𝑒𝑖(𝑑)β‰€πœ”π‘–<∞,π‘‘βˆˆ[0,𝑑end], for 𝑖=1,2,}, is an admissible control set.

Next, we describe the role played by the first control 𝑒1(𝑑) and how it is incorporated in the system (1). It the beginning people start smoking occasionally for a variety of different reasons. Some think it looks cool. Others start because their family members or friends smoke. Statistics show that about 9 out of 10 tobacco users start occasionally and then gradually become persistent smoker. Furthermore there are strong evidence that the attitude towards smoking is starting from high (junior) school time. Thus, a campaign like movie particularly in a class or group to show them that smoking is more likely to develop cancer of the mouth and throat, bladder, kidney, liver, stomach, and pancreas which significantly raises risk of heart attack. Therefore control function 𝑒1(𝑑) represents education campaign level used to control smoking in a community. The control 𝑒2(𝑑) represents the level of treatment in the form of stop smoking injection (vaccine) that contains drugs that block nicotine receptors of the brain, which in turn helps in reducing the desire to smoke. The physical meaning of the control variable in this problem is that low levels of the number of occasional (light) and persistent smokers build. In case of no campaign or no treatment the number of persistent smoker increases while the number of former smoker decreases. Better and prefect time of campaign brings the number of occasional (light) and persistent smokers to a small level, potential smokers begin to build again and more individuals are recovered by quitting the use of tobacco. The effects of tobacco on both occasional (light) and persistent smokers are negative for quit smokers around them, so we wish to minimize both. Also small amount of control variables (campaign and treatment) are acceptable; therefore, we wish to penalize for amount too large, so quadratic terms for both control variables will be analyzed. Hence, our optimal control problem to minimize the objective functional is given by𝐽𝑒1,𝑒2ξ€Έ=ξ€œπ‘‡0𝛼1𝐿(𝑑)+𝛼21𝑆(𝑑)+2ξ€·πœ‰1𝑒21(𝑑)+πœ‰2𝑒22ξ€Έξ‚„(𝑑)𝑑𝑑(7) subject to 𝑑𝑃(𝑑)𝑑𝑑=𝑏𝑁(𝑑)βˆ’π›½1𝑑𝐿(𝑑)𝑃(𝑑)βˆ’1ξ€Έ+πœ‡π‘ƒ(𝑑)+πœŒπ‘„(𝑑),𝑑𝐿(𝑑)𝑑𝑑=𝛽1𝐿(𝑑)𝑃(𝑑)βˆ’π›½2𝑑𝐿(𝑑)𝑆(𝑑)βˆ’2+πœ‡+𝑒1ξ€Έ(𝑑)𝐿(𝑑),𝑑𝑆(𝑑)𝑑𝑑=𝛽2𝐿(𝑑)𝑆(𝑑)βˆ’π›Ύ+𝑑3+πœ‡+𝑒2ξ€Έ(𝑑)𝑆(𝑑),𝑑𝑄(𝑑)=𝑑𝑑𝛾+𝑒2𝑑(𝑑)𝑆(𝑑)βˆ’4ξ€Έ+πœ‡+πœŒπ‘„(𝑑)+𝑒1(𝑑)𝐿(𝑑),(8) with initial conditions𝑃(0)=𝑃0,𝐿(0)=𝐿0,𝑆(0)=𝑆0,𝑄(0)=𝑄0.(9) Here 𝛼𝑖 and πœ‰π‘– for 𝑖=1,2 are weight factors (positive constants) that represent a smoker level of acceptance of the control campaign. The goal is to minimize the occasional (light) and persistent smokers and the cost of implementing the control. Occasional smokers induce an optimal control 𝑒1(𝑑) before they become persistent smokers, and an optimal control 𝑒2(𝑑) should be provided to persistent smokers population.

Various deterministic optimal control models and their existence are investigated by several authors [8, 9, 15]. First we shall show the existence for the control system (8). Let 𝑃(𝑑), 𝐿(𝑑), 𝑆(𝑑), and 𝑄(𝑑) be state variables with control variable 𝑒(𝑑)=(𝑒1(𝑑),𝑒2(𝑑))βˆˆπ‘ˆ. For existence we consider the control system (8) with initial conditions in (9). Then we can rewrite the system (8) in the following form:𝑉𝑑=𝑀𝑉+πœ’(𝑉),(10) where βŽ‘βŽ’βŽ’βŽ’βŽ£π‘„βŽ€βŽ₯βŽ₯βŽ₯βŽ¦βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ£π‘‰=𝑃(𝑑)𝐿(𝑑)𝑆(𝑑)(𝑑),𝑀=βˆ’π‘‘1ξ€·π‘‘βˆ’πœ‡00𝜌0βˆ’2+πœ‡+𝑒1ξ€Έξ€·(𝑑)0000βˆ’π›Ύ+𝑑3+πœ‡+𝑒2ξ€Έ0(𝑑)0𝑒1(𝑑)𝛾+𝑒2𝑑(𝑑)βˆ’4ξ€ΈβŽ€βŽ₯βŽ₯βŽ₯βŽ₯⎦,⎑⎒⎒⎒⎣+πœ‡+πœŒπœ’(𝑉)=𝑏𝑁(𝑑)βˆ’π›½1𝛽𝑃(𝑑)𝐿(𝑑)1𝑃(𝑑)𝐿(𝑑)+𝛽2𝛽𝐿(𝑑)𝑆(𝑑)2𝐿0⎀βŽ₯βŽ₯βŽ₯⎦,(𝑑)𝑆(𝑑)(11)and 𝑉𝑑 denotes the derivative of 𝑉 with respect to the time 𝑑. The system (10) is a nonlinear system with a bounded coefficient. We set𝐺(𝑉)=𝑀𝑉+πœ’(𝑉).(12) The second term on the right-hand side of (12) satisfies ||πœ’ξ€·π‘‰1ξ€Έξ€·π‘‰βˆ’πœ’2ξ€Έ||≀Λ1||𝑃1(𝑑)βˆ’π‘ƒ2||(𝑑)+Ξ›2||𝐿1(𝑑)βˆ’πΏ2||(𝑑)+Ξ›3||𝑆1(𝑑)βˆ’π‘†2||,(𝑑)(13) where the positive constants Ξ›1, Ξ›2, and Ξ›3 are independent of state variables 𝑃(𝑑), 𝐿(𝑑), and 𝑆(𝑑)≀𝑁(𝑑), respectively. Also we get ||𝐺𝑉1ξ€Έξ€·π‘‰βˆ’πΊ2ξ€Έ||||𝑉≀Λ1βˆ’π‘‰2||,(14) where Ξ›=Ξ›1+Ξ›2+Ξ›3+‖𝑀‖<∞. Thus it follows that the function 𝐺 is uniformly Lipschitz continuous. From the definition of π‘ˆ and the restriction on 𝑃(𝑑), 𝐿(𝑑), 𝑆(𝑑), and 𝑄(𝑑)β‰₯0 we can see that a solution of the system (10) exists (see [16]).

To illustrate how to solve a control problem actually, first we should find Hamiltonian for the optimal control problem (7)–(9). To do this, we define the Hamiltonian 𝐻 for the control problem as follows:𝐻𝑃,𝐿,𝑆,𝑄,𝑒1,𝑒2,πœ†1,πœ†2,πœ†3,πœ†4ξ€Έ=𝐿,𝑆,𝑒1,𝑒2ξ€Έ+πœ†1(𝑑)𝑔1(𝑑)+πœ†2(𝑑)𝑔2(𝑑)+πœ†3(𝑑)𝑔3(𝑑),(15) where (𝐿,𝑆,𝑒1,𝑒2)=𝛼1𝐿(𝑑)+𝛼2𝑆(𝑑)+(1/2)(πœ‰1𝑒21(𝑑)+πœ‰2𝑒22(𝑑)) is the Lagrangian and 𝑔𝑖 for 𝑖=1,2,3 is the right-hand side of the differential equations of the state system (8), respectively.

Theorem 4. There exists an optimal control pair π‘’βˆ—=(π‘’βˆ—1,π‘’βˆ—2)βˆˆπ‘ˆ such that π½ξ€·π‘’βˆ—1,π‘’βˆ—2ξ€Έ=min𝑒1,𝑒2ξ€Έβˆˆπ‘ˆπ½ξ€·π‘’1,𝑒2ξ€Έ,(16) subject to the control system (8) with initial conditions (9).

Proof. To prove the existence of an optimal control, we have to show the following. (1) The control and state variables are nonnegative values.(2) The control π‘ˆ set is convex and closed.(3) The RHS of the state system is bounded by linear function in the state and control variables.(4) The integrand of the objective functional is concave on π‘ˆ.(5) There exist constants such that the integrand in (7) of the objective functional is satisfied.
In order to verify these conditions, we use a result by Fister [5]. We note that the solutions are bounded. The set of all the control variables 𝑒(𝑑)=(𝑒1(𝑑),𝑒2(𝑑))βˆˆπ‘ˆ, is also convex and closed by definition. The optimal system is bounded which determines the compactness needed for the existence of the optimal control. In addition, the integrand in functional (7) is given by 𝛼1𝐿(𝑑)+𝛼2𝑆(𝑑)+(1/2)(πœ‰1𝑒21(𝑑)+πœ‰2𝑒22(𝑑)) and is convex on the control set π‘ˆ. Also we can easily see that there exist a constant 𝜎>1 and positive numbers πœ‚1 and πœ‚2 such that 𝐽𝑒1,𝑒2ξ€Έβ‰₯πœ‚2+πœ‚1ξ‚€||𝑒1||2+||𝑒2||2ξ‚πœŽ/2,(17) which completes the existence of an optimal control.

Pontryagin maximum principal (PMP) states necessary conditions [17] that must hold on an optimal trajectory. These conditions help us in calculation of the state, adjoint, and control variables in the optimal control problems. PMP can be used as both a computational technique and analytic technique. The technical tools used to determine the optimal strategy is given in the following form.

If (π‘₯βˆ—(𝑑),π‘’βˆ—(𝑑)) is an optimal solution of an optimal control problem, then there exists a nontrivial vector function πœ†(𝑑)=(πœ†1(𝑑),πœ†2(𝑑),…,πœ†π‘›(𝑑)) satisfying the following inequalities: π‘₯ξ…žξ€·(𝑑)=πœ•π»π‘‘,π‘₯βˆ—(𝑑),π‘’βˆ—ξ€Έ(𝑑),πœ†(𝑑),ξ€·πœ•πœ†0=πœ•H𝑑,π‘₯βˆ—(𝑑),π‘’βˆ—ξ€Έ(𝑑),πœ†(𝑑),πœ†πœ•π‘’ξ…žξ€·(𝑑)=βˆ’πœ•π»π‘‘,π‘₯βˆ—(𝑑),π‘’βˆ—ξ€Έ(𝑑),πœ†(𝑑),πœ•π‘₯(18) where β€œβ€²β€ denotes the derivative with respect to time 𝑑. Now we apply the necessary conditions to the Hamiltonian 𝐻 in (15).

Theorem 5. Let π‘ƒβˆ—(𝑑), πΏβˆ—(𝑑), π‘†βˆ—(𝑑), and π‘„βˆ—(𝑑) be optimal state solutions with associated optimal control variables π‘’βˆ—1(𝑑) and π‘’βˆ—2(𝑑) for the optimal control problem (7)–(9). Then there exist adjoint variables πœ†1(𝑑), πœ†2(𝑑), πœ†3(𝑑), and πœ†4(𝑑) satisfing πœ†ξ…ž1(𝑑)=𝛽1ξ€·πœ†1(𝑑)βˆ’πœ†2𝐿(𝑑)βˆ—(𝑑)+πœ†1𝑑(𝑑)1ξ€Έ,πœ†+πœ‡ξ…ž2(𝑑)=𝛽1ξ€·πœ†1(𝑑)βˆ’πœ†2(𝑃𝑑)βˆ—(𝑑)+𝛽2ξ€·πœ†2(𝑑)βˆ’πœ†3(𝑆𝑑)βˆ—(𝑑)+πœ†2𝑑(𝑑)2ξ€Έ+ξ€·πœ†+πœ‡2(𝑑)βˆ’πœ†4𝑒(𝑑)1(𝑑)βˆ’π›Ό1,πœ†ξ…ž3(𝑑)=𝛽2ξ€·πœ†2(𝑑)βˆ’πœ†3(𝐿𝑑)βˆ—(𝑑)+πœ†3(𝑑)𝛾+𝑑3ξ€Έ+ξ€·πœ†+πœ‡3(𝑑)βˆ’πœ†4𝑒(𝑑)2(𝑑)βˆ’π›Ό2,πœ†β€²4(ξ€·πœ†π‘‘)=4(𝑑)βˆ’πœ†1(𝑑)𝜌+πœ†4(𝑑𝑑)4ξ€Έ,+πœ‡(19) with transversality conditions (or boundary conditions) πœ†π‘–(𝑇)=0,𝑖=1,2,3.(20) Furthermore, optimal control pairs are given as follows: π‘’βˆ—1ξƒ―ξƒ―ξ€·πœ†(𝑑)=maxmin1(𝑑)βˆ’πœ†3ξ€Έ(𝑑)𝐿(𝑑)πœ‰1,πœ”1𝑒,0,(21)βˆ—2ξƒ―ξƒ―ξ€·πœ†(𝑑)=maxmin2(𝑑)βˆ’πœ†3𝑆(𝑑)(𝑑)πœ‰2,πœ”2ξƒ°ξƒ°,0.(22)

Proof. For the proof of this theorem see Appendix B.

Here we call formulas (21) and (22) for π‘’βˆ— the characterization of the optimal control. The optimal control and the state are found by solving the optimality system, which consists of the state system (8), the adjoint system (20), boundary conditions (9) and (20), and the characterization of the optimal control (21) and (22). To solve the optimality system we use the initial and transversality conditions together with the characterization of the optimal control pair (π‘’βˆ—1(𝑑),π‘’βˆ—2(𝑑)) given by (21) and (22). By substituting the values of π‘’βˆ—1(𝑑) and π‘’βˆ—2(𝑑) in the control system (8) we get the following system:π‘‘π‘ƒβˆ—(𝑑)𝑑𝑑=π‘π‘βˆ—(𝑑)βˆ’π›½1πΏβˆ—(𝑑)π‘ƒβˆ—ξ€·π‘‘(𝑑)βˆ’1𝑃+πœ‡βˆ—(𝑑)+πœŒπ‘„βˆ—(𝑑),π‘‘πΏβˆ—(𝑑)𝑑𝑑=𝛽1πΏβˆ—(𝑑)π‘ƒβˆ—(𝑑)βˆ’π›½2πΏβˆ—(𝑑)π‘†βˆ—βˆ’ξƒ©π‘‘(𝑑)2ξƒ―ξƒ―ξ€·πœ†+πœ‡+maxmin2(𝑑)βˆ’πœ†4𝐿(𝑑)βˆ—(𝑑)πœ‰1,πœ”1ξƒ°,0ξƒ°ξƒͺΓ—πΏβˆ—(𝑑),π‘‘π‘†βˆ—(𝑑)𝑑𝑑=πœ‰2πΏβˆ—(𝑑)π‘†βˆ—βˆ’ξƒ©(𝑑)𝛾+𝑑3ξƒ―ξƒ―ξ€·πœ†+πœ‡+maxmin3(𝑑)βˆ’πœ†4𝑆(𝑑)βˆ—(𝑑)πœ‰2,πœ”2ξƒ°,0ξƒ°ξƒͺΓ—π‘†βˆ—(𝑑),π‘‘π‘„βˆ—(𝑑)=ξƒ©ξƒ―ξƒ―ξ€·πœ†π‘‘π‘‘π›Ύ+maxmin3(𝑑)βˆ’πœ†4(𝑆𝑑)βˆ—(𝑑)πœ‰2,πœ”2𝑆,0ξƒ°ξƒͺβˆ—βˆ’ξ€·π‘‘(𝑑)4𝑄+πœ‡+πœŒβˆ—ξƒ―ξƒ―ξ€·πœ†(𝑑)+maxmin2(𝑑)βˆ’πœ†4𝐿(𝑑)βˆ—(𝑑)πœ‰1,πœ”1𝐿,0βˆ—(𝑑).(23) With the Hamiltonian π»βˆ— at (π‘ƒβˆ—,πΏβˆ—, π‘†βˆ—, π‘„βˆ—, π‘’βˆ—1, π‘’βˆ—2, πœ†1, πœ†2, πœ†3, πœ†4),π»βˆ—=𝛼1πΏβˆ—(𝑑)+𝛼2π‘†βˆ—+1(𝑑)2βŽ‘βŽ’βŽ’βŽ£πœ‰1ξƒ©ξƒ―ξƒ―ξ€·πœ†maxmin2(𝑑)βˆ’πœ†4𝐿(𝑑)βˆ—(𝑑)πœ‰1,πœ”1ξƒ°,0ξƒ°ξƒͺ2+πœ‰2ξƒ©ξƒ―ξƒ―ξ€·πœ†maxmin3(𝑑)βˆ’πœ†4𝑆(𝑑)βˆ—(𝑑)πœ‰2,πœ”2ξƒ°,0ξƒ°ξƒͺ2⎀βŽ₯βŽ₯⎦+πœ†1(𝑑)π‘”βˆ—1(𝑑)+πœ†2(𝑑)π‘”βˆ—2(𝑑)+πœ†3(𝑑)π‘”βˆ—3(𝑑),(24) where π‘”βˆ—π‘–(𝑑) for 𝑖=1,2,3 is the right-hand side of the differential equations of the state system (23), respectively. To find out the optimal control and state we will numerically solve the system (23) and Hamiltonian (24).

3. Numerical Results

In this section, we shall solve numerically an optimal control problem for the 𝑃𝐿𝑆𝑄 model. Here we obtain the optimality system from the state and adjoint equations. The optimal control problem strategy is obtained by solving the optimal system, which consists of eight ordinary differential equations and boundary conditions. The optimality system can be solved by using an iterative method by Runge-Kutta fourth scheme [18]. Using an initial guess for the control variables, 𝑒1(𝑑) and 𝑒2(𝑑), the state variables, 𝑃, 𝐿, 𝑆, and 𝑄 are solved forward in time, and then the adjoint variables, πœ†π‘– for 𝑖=1,2,3,4, are solved backwards in time. If the new values of the state and adjoint variables differ from the previous values, the new values are used to update 𝑒1(𝑑) and 𝑒2(𝑑), and the process is repeated until the system converges. For the 𝑃𝐿𝑆𝑄 model presented in this work, the state system is given by (8) with the initial conditions given by (9). The adjoint system is given by (20) with the final time conditions given by (21) and the characterization of the optimal control by (21) and (22).

To determine the numerical simulation of the optimality system we use a small time step size Δ𝑑=0.005. We consider that the optimal campaign continues for 30 day, and use a real estimate of parameter value represented in Table 1. For reasonable values of the parameter we restate the idea presented in [8, 13, 19]. When a person first becomes a smoker it is not likely that she/he quits for several years since tobacco contains nicotine, which is shown to be an addictive drug. We assume 1/𝛾 to be a value between 15 and 25 years, that is, 20. Hence, the value of 1/𝛾 is set to 7300 days. The death rate of each individual is different from others and depends on the real life situation. Here, we assume that 𝑑4≀𝑑1≀𝑑2≀𝑑3. The natural death rate πœ‡ is per 1000 per year (currently 8 in the US). There is strong evidence that the attitude towards smoking is starting from high (junior) school time. Therefore 1/𝑏 is set to be 1095 days. Death from lung cancer was the leading cause by smoking, with a rate of 37 per 100,000 individuals [20]. The reasonable value of parameters 𝛽1 and 𝛽2 which represent the average number of contacts with smoker individuals needed to make light and persistent smoker, respectively, in each unit of time can be obtained by using the same techniques presented in [13]. We consider the real data used in [8] for the 𝑃(0)=153,  𝑆(0)=79, and assume that 𝐿(0)=40. and 𝑄(0)=9.

Table 1: Reasonable values of the parameter.

We presented in Figure 1 the potential smoker in two systems (1) and (8). The plain line denotes the population of potential smokers in the system (1) without control while the dotted line denotes the population of potential smokers in the system (8) with control. The potential smokers in the two systems (1) and (8) are almost the same in around days 1–10. After 10 days a slight increase appears in potential smokers population. Figure 2 represents the population of occasional smokers in two systems (1) without control and (8) with control. The population of occasional smokers sharply decreases from the first day of infection in systems both with control and without control and reaches its minimum number of occasional smokers 𝐿=39.3 and πΏβˆ—=39 at the end of the optimal campaign.

Figure 1: The plot shows the population of potential smokers 𝑃 both with control and without control.
Figure 2: The plot shows occasional smoker 𝐿 both with control and without control.

Figure 3 represents the persistent smokers in the two systems (1) and (8). The persistent smokers population in both systems (1) and (8) decreases more sharply than the light (occasional) smokers population and about 𝑆=19 and π‘†βˆ—=16.8 persistent smokers population remains at the end of the optimal campaign. Figure 4 represents the numerical results of both systems (1) and (8) of the quit smokers population. After the optimal campaign the light (occasional) smokers and persistent smokers populations decreases while the quit smoker population increases.

Figure 3: The plot shows persistent smoker 𝑆 both with control and without control.
Figure 4: The plot shows quit smoker 𝑄 both with control and without control.

In this work we use one set of parameters for both dynamical systems (1) and (8). Simulations with different sets of parameter values can be used in the future to obtain a sampling of possible behaviors of a dynamical system.

4. Conclusion

In this paper we considered the model introduced by Zaman [13] and extend that once a smoker quits he/she may becomes a potential smoker again, and then we used optimal control theory to minimize the number of light and persistent smokers and maximize the number of quit smokers in a community. The control plots obtained from the numerical simulation in this paper represent that the numbers of light and persistent smokers decreases and the number of quit smokers increases in the optimality system. We also showed that for certain values of control rate there exists its corresponding optimal solution. We considered a special disease (smoking) in a specific community as a realistic model, and we hope that the approach introduced in this paper will be applicable in other epidemic models beyond the giving up smoking model.


A. Proof of Theorem 3

The Jacobian matrix of the giving up smoking model (1) around 𝐸⋆=(𝑃⋆,𝐿⋆,𝑆⋆,𝑄⋆) is 𝐽⋆=βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ£π΄π‘ƒπ΄πΏπ΅0πœŒπ‘ƒπ΅πΏπ΅π‘†0π‘œπΆπΏπΆπ‘†000π›Ύπ·π‘„βŽ€βŽ₯βŽ₯βŽ₯βŽ₯⎦,(A.1) where 𝐴𝑃=βˆ’π›½1πΏβ‹†βˆ’ξ€·π‘‘1ξ€Έ+πœ‡,𝐴𝐿=βˆ’π›½1𝑃⋆,𝐡𝑃=𝛽1𝐿⋆,𝐡𝑃=𝛽1π‘ƒβ‹†βˆ’π›½2π‘†β‹†βˆ’ξ€·π‘‘2ξ€Έ,𝐡+πœ‡π‘†=βˆ’π›½2𝐿⋆,𝐢𝐿=𝛽2𝑆⋆,𝐢𝑆=𝛽2πΏβ‹†βˆ’ξ€·π›Ύ+𝑑3ξ€Έ+πœ‡,𝐷𝑄𝑑=βˆ’4ξ€Έ.+πœ‡+𝜌(A.2) The characteristic equation of 𝐽⋆ is πœ†4+π‘Ž1πœ†3+π‘Ž2πœ†2+π‘Ž3πœ†+π‘Ž4=0,(A.3) where π‘Ž1𝐴=βˆ’π‘ƒ+𝐡𝐿+𝐢𝑆+𝐷𝑄,π‘Ž2=𝐢𝑆𝐷𝑄+𝐴𝑃𝐡𝐢+𝐢𝑆+𝐷𝑄𝐴𝑃+π΅πΏξ€Έβˆ’π΄πΏπ΅π‘ƒ,π‘Ž3=𝐴𝐿𝐡𝑃𝐢𝑆+π·π‘„ξ€Έβˆ’ξ€·πΆπ‘†π·π‘„ξ€·π΄π‘ƒ+𝐡𝐿+𝐴𝑃𝐡𝐢𝐢𝑆+𝐷𝑄,π‘Žξ€Έξ€Έ4=π΄π‘ƒπ΅πΏπΆπ‘†π·π‘„βˆ’π΄πΏπ΅π‘ƒπΆπ‘†π·π‘„βˆ’πœŒπ›Ύπ΅π‘ƒπΆπΏ.(A.4) By Routh-Hurwitz conditions the equilibrium state is asymptotically stable if π‘Ž1>0, π‘Ž3>0, π‘Ž4>0 and π‘Ž1π‘Ž2π‘Ž3>π‘Ž23+π‘Ž21π‘Ž4. Thus it follows that the unique positive equilibrium 𝐸⋆ of the giving up smoking model, which exists if 𝐾>1, is always locally asymptotically stable. This completes the proof.

B. Proof of Theorem 5

To determine the adjoint equations and the transversality conditions, we use Hamiltonian (15). From setting 𝑃(𝑑)=π‘ƒβˆ—(𝑑), 𝐿(𝑑)=πΏβˆ—(𝑑), 𝑆(𝑑)=π‘†βˆ—(𝑑), and 𝑄(𝑑)=π‘„βˆ—(𝑑) and differentiating the Hamiltonian (15) with respect to 𝑃, 𝐿, 𝑆, and 𝑄, we obtain πœ†ξ…ž1(𝑑)=βˆ’πœ•π»πœ•π‘ƒ=𝛽1ξ€·πœ†1(𝑑)βˆ’πœ†2ξ€Έ(𝑑)𝐿(𝑑)+πœ†1𝑑(𝑑)1ξ€Έ,πœ†+πœ‡ξ…ž2(𝑑)=βˆ’πœ•π»πœ•πΏ=𝛽1ξ€·πœ†1(𝑑)βˆ’πœ†2𝑃(𝑑)(𝑑)+𝛽2ξ€·πœ†2(𝑑)βˆ’πœ†3(𝑑)𝑆(𝑑)+πœ†2(𝑑𝑑)2ξ€Έ+ξ€·πœ†+πœ‡2(𝑑)βˆ’πœ†4𝑒(𝑑)1(𝑑)βˆ’π›Ό1,πœ†ξ…ž3(𝑑)=βˆ’πœ•π»πœ•π‘†=𝛽2ξ€·πœ†2(𝑑)βˆ’πœ†3ξ€Έ(𝑑)𝐿(𝑑)+πœ†3ξ€·(𝑑)𝛾+𝑑3ξ€Έ+ξ€·πœ†+πœ‡3(𝑑)βˆ’πœ†4𝑒(𝑑)2(𝑑)βˆ’π›Ό2,πœ†ξ…ž4(𝑑)=βˆ’πœ•π»=ξ€·πœ†πœ•π‘„4(𝑑)βˆ’πœ†1(𝑑)𝜌+πœ†4(𝑑𝑑)4ξ€Έ.+πœ‡(B.1) Using the optimality conditions and the property of the control space π‘ˆ for the control variable 𝑒1, we obtain πœ•π»πœ•π‘’1=0βŸΉπœ‰1𝑒1ξ€·(𝑑)+βˆ’πœ†2(𝑑)+πœ†4ξ€Έ(𝑑)𝐿(𝑑)=0.(B.2) By using the property of the control space, we obtain π‘’βˆ—1⎧βŽͺβŽͺβŽͺβŽͺ⎨βŽͺβŽͺβŽͺβŽͺβŽ©ξ€·πœ†(𝑑)=0if2(𝑑)βˆ’πœ†4𝐿(𝑑)(𝑑)πœ‰1ξ€·πœ†β‰€0,2(𝑑)βˆ’πœ†4𝐿(𝑑)(𝑑)πœ‰1ξ€·πœ†if0<2(𝑑)βˆ’πœ†4𝐿(𝑑)(𝑑)πœ‰1,<πœ”1πœ”1ξ€·πœ†if2(𝑑)βˆ’πœ†4ξ€Έ(𝑑)𝐿(𝑑)πœ‰1β‰₯πœ”1.(B.3) This can be rewritten in compact notationπ‘’βˆ—1ξƒ―ξƒ―ξ€·πœ†(𝑑)=maxmin2(𝑑)βˆ’πœ†4ξ€Έ(𝑑)𝐿(𝑑)πœ‰1,πœ”1ξƒ°ξƒ°,0.(B.4) Similar to the control variable 𝑒2, we get π‘’βˆ—2⎧βŽͺβŽͺβŽͺβŽͺ⎨βŽͺβŽͺβŽͺβŽͺβŽ©ξ€·πœ†(𝑑)=0if3(𝑑)βˆ’πœ†4𝑆(𝑑)(𝑑)πœ‰2ξ€·πœ†β‰€0,3(𝑑)βˆ’πœ†4𝑆(𝑑)(𝑑)πœ‰2ξ€·πœ†if0<3(𝑑)βˆ’πœ†4𝑆(𝑑)(𝑑)πœ‰2<πœ”2,πœ”2ξ€·πœ†if3(𝑑)βˆ’πœ†4ξ€Έ(𝑑)𝑆(𝑑)πœ‰2β‰₯πœ”2.(B.5) This can also be rewritten similar to the above compact form asπ‘’βˆ—2ξƒ―ξƒ―ξ€·πœ†(𝑑)=maxmin3(𝑑)βˆ’πœ†4ξ€Έ(𝑑)𝑆(𝑑)πœ‰2,πœ”2ξƒ°ξƒ°,0.(B.6) This completes the proof.


The author wishes to thank two anonymous referees for helpful suggestion to improve the quality of this paper.


  1. C. Castilho, β€œOptimal control of an epidemic through educational campaigns,” Electronic Journal of Differential Equations, vol. 2006, pp. 1–11, 2006. View at Google Scholar Β· View at Scopus
  2. N. G. Becker and D. N. Starczak, β€œOptimal vaccination strategies for a community of households,” Mathematical Biosciences, vol. 139, no. 2, pp. 117–132, 1997. View at Publisher Β· View at Google Scholar Β· View at Scopus
  3. M. L. Brandeau, G. S. Zaric, and A. Richter, β€œResource allocation for control of infectious diseases in multiple independent populations: beyond cost-effectiveness analysis,” Journal of Health Economics, vol. 22, no. 4, pp. 575–598, 2003. View at Publisher Β· View at Google Scholar Β· View at PubMed Β· View at Scopus
  4. K. Wickwire, β€œMathematical models for the control of pests and infectious diseases: a survey,” Theoretical Population Biology, vol. 11, no. 2, pp. 182–238, 1977. View at Google Scholar Β· View at Scopus
  5. K. R. Fister, S. Lenhart, and J. S. McNally, β€œOptimizing chemotherapy in an HIV model,” Electronic Journal of Differential Equations, vol. 1998, no. 32, pp. 1–12, 1998. View at Google Scholar
  6. D. Kirschner, S. Lenhart, and S. Serbin, β€œOptimal control of the chemotherapy of HIV,” Journal of Mathematical Biology, vol. 35, no. 7, pp. 775–792, 1997. View at Google Scholar
  7. H. R. Joshi, β€œOptimal control of an HIV immunology model,” Optimal Control Applications and Methods, vol. 23, no. 4, pp. 199–213, 2002. View at Publisher Β· View at Google Scholar
  8. G. Zaman, Y. H. Kang, and I. H. Jung, β€œStability analysis and optimal vaccination of an SIR epidemic model,” BioSystems, vol. 93, no. 3, pp. 240–249, 2008. View at Publisher Β· View at Google Scholar Β· View at PubMed
  9. Y. H. Kang, S. Lenhart, and V. Protopopescu, β€œOptimal control of parameters and input functions for nonlinear systems,” Houston Journal of Mathematics, vol. 33, no. 4, pp. 1231–1256, 2007. View at Google Scholar
  10. K. Dietz and D. Schenzle, β€œMathematical models for infectious disease statistics,” in A Celebrationf of Statistics, ISI Centenary, A. C. Atkinson and S. E. Fienberg, Eds., pp. 167–204, Springer, New York, NY, USA, 1985. View at Google Scholar
  11. C. Castillo-Garsow, G. Jordan-Salivia, and A. Rodriguez Herrera, β€œMathematical models for the dynamics of tobacco use, recovery, and relapse,” Technical Report Series BU-1505-M, Cornell University, Ithaca, NY, USA, 2000. View at Google Scholar
  12. O. Sharomi and A. B. Gumel, β€œCurtailing smoking dynamics: a mathematical modeling approach,” Applied Mathematics and Computation, vol. 195, no. 2, pp. 475–499, 2008. View at Publisher Β· View at Google Scholar
  13. G. Zaman, β€œQualitative behavior of giving up smoking model,” Bulletin of the Malaysian Mathematical Sciences Society. In press.
  14. A. J. Arenas, J. A. Moraño, and J. C. Cortés, β€œNon-standard numerical method for a mathematical model of RSV epidemiological transmission,” Computers and Mathematics with Applications, vol. 56, no. 3, pp. 670–678, 2008. View at Publisher Β· View at Google Scholar
  15. G. Zaman, Y. H. Kang, and I. H. Jung, β€œOptimal treatment of an SIR epidemic model with time delay,” BioSystems, vol. 98, no. 1, pp. 43–50, 2009. View at Publisher Β· View at Google Scholar Β· View at PubMed
  16. G. Birkhoff and G. C. Rota, Ordinary Differential Equations, John Wiley & Sons, New York, NY, USA, 4th edition, 1989.
  17. I. K. Morton and L. S. Nancy, Dynamics Optimization: The Calculus of Variations and Optimal Control in Economics and Management, Elsevier Science, New York, NY, USA, 2000.
  18. S. Lenhart and J. T. Workman, Optimal Control Applied to Biological Models, Mathematical and Computational Biology Series, Chapman & Hall/CRC, London, UK, 2007.
  19. American National Institute of Drug Abuse, β€œCigarettes and Other Nicotine Products,”
  20. National Cancer Information Center, β€œDeath statistics based on vital registration,” February 2006,