Optimal Campaign in the Smoking Dynamics

Abstract

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) 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 𝛽₁ and 𝛽₂ 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 𝑑₁,  𝑑₂,  𝑑₃, and 𝑑₄ 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 𝛽₁ and 𝛽₂. Potential smokers acquire infection at per capita rate 𝛽₁𝐿(𝑡). The total population size at time 𝑡 is denoted by 𝑁(𝑡) with 𝑁(𝑡)=𝑃(𝑡)+𝐿(𝑡)+𝑆(𝑡)+𝑄(𝑡). The dynamics of the total population are governed by the following differential equation:𝑑𝑁𝑑𝑡=(𝑏−𝜇)𝑁(𝑡)−𝑑¹𝑃(𝑡)+𝑑₂𝐿(𝑡)+𝑑₃𝑆(𝑡)+𝑑₄𝑄(𝑡).(2) If all individuals die at the same death rates, then we have 𝑑₁=𝑑₂=𝑑₃=𝑑₄=𝑑, 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𝑃^⋆=𝛽²𝑆^⋆+𝑑₂+𝜇𝛽₁𝐿_⋆,=𝛾+𝑑³+𝜇𝛽₂𝑆_⋆,=𝑑²+𝜇𝑑₄+𝜇+𝜌𝛽₁𝛾+𝑑₃+𝜇+𝛽₂𝑑₁+𝜇𝑄×𝛽₂(𝐾𝛽−₁𝜌1𝛾)−,_⋆𝛽₁𝛾+𝑑₃+𝜌−𝛽₂𝑑₁+𝜇𝑑₄+𝜇+𝜌=𝑑²+𝜇𝛽₁𝛾+𝑑₃+𝜇+𝛽₂𝑑₁+𝜇×𝛽₂(𝐾𝛽−₁𝜌1𝛾)−,𝛽₁𝛾+𝑑₃+𝜌−𝛽₂𝑑₁+𝜇𝑑₄+𝜇+𝜌(4) where 𝐾=𝑏𝛽₁𝛽₂𝑁^⋆/((𝑑₂+𝜇)(𝛽₁(𝛾+𝑑₃+𝜌)+𝛽₂(𝑑₁+𝜇))) 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 𝐽_𝑡=⎡⎢⎢⎢⎣−𝑑¹0−𝜇−𝜇0−𝑑₂0𝜌00−𝑑₃0000−𝛾𝛾−𝜇−𝑑₄0−𝜇−𝜌⎤⎥⎥⎥⎦.(5) The eigenvalues of the Jacobian matrix around the trivial equilibrium 𝐸₀=(0,0,0,0) are −𝑑₁−𝜇, −𝑑₂−𝜇, −𝛾−𝑑₃−𝜇, −𝑑₄−𝜇−𝜌. 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 𝛽₁<𝑑₂+𝜇. 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: 𝐽_𝑙=⎡⎢⎢⎢⎣−𝑑¹0−𝜇𝛽₁−𝛽₁0𝜌0−𝜇0−𝑑₂−𝑑₃0000−𝛾𝛾−𝜇−𝑑₄0−𝜇−𝜌⎤⎥⎥⎥⎦.(6) The eigenvalues of the Jacobian matrix 𝐽_𝑙 around smoking-free equilibrium 𝐸_𝑙=(1,0,0,0) are −𝑑₁−𝜇, 𝛽₁−𝑑₂−𝜇, −𝛾−𝑑₃−𝜇, −𝑑₄−𝜇−𝜌. Thus, we deduce that all the roots have negative real part when 𝛽₁<𝑑₂+𝜇 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 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡))∈𝑈 relative to the state variables (𝑃(𝑡),𝐿(𝑡),𝑆(𝑡),𝑄(𝑡)), where 𝑈={(𝑢₁,𝑢₂)∣𝑢_𝑖(𝑡) is measurable and 0≤𝑢_𝑖(𝑡)≤𝜔_𝑖<∞,𝑡∈[0,𝑡_end], for 𝑖=1,2,}, is an admissible control set.

Next, we describe the role played by the first control 𝑢₁(𝑡) 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 𝑢₁(𝑡) represents education campaign level used to control smoking in a community. The control 𝑢₂(𝑡) 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𝐽𝑢₁,𝑢₂=^𝑇0𝛼¹𝐿(𝑡)+𝛼₂𝑆(𝑡)+12𝜉¹𝑢²¹(𝑡)+𝜉₂𝑢²²(𝑡)𝑑𝑡(7) subject to 𝑑𝑃(𝑡)𝑑𝑡=𝑏𝑁(𝑡)−𝛽¹𝑑𝐿(𝑡)𝐿(𝑡)𝑃(𝑡)−𝑑₁+𝜇𝑃(𝑡)+𝜌𝑄(𝑡),𝑑𝑡=𝛽¹𝑑𝑆(𝑡)𝐿(𝑡)𝑃(𝑡)−𝛽₂𝐿(𝑡)𝑆(𝑡)−𝑑₂+𝜇+𝑢₁(𝑡)𝐿(𝑡),𝑑𝑡=𝛽²𝑑𝑄(𝑡)𝐿(𝑡)𝑆(𝑡)−𝛾+𝑑₃+𝜇+𝑢₂(𝑡)𝑆(𝑡),𝑑𝑡=𝛾+𝑢²(𝑡)𝑆(𝑡)−𝑑₄+𝜇+𝜌𝑄(𝑡)+𝑢₁(𝑡)𝐿(𝑡),(8) with initial conditions𝑃(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 𝑢₁(𝑡) before they become persistent smokers, and an optimal control 𝑢₂(𝑡) 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 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡))∈𝑈. 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 𝑉=⎡⎢⎢⎢⎣𝑄𝑃𝑆𝐿((((𝑡𝑡𝑡𝑡))))⎤⎥⎥⎥⎦,𝑀=⎡⎢⎢⎢⎢⎣−𝑑¹0−𝜇−𝑑₂00𝜌0+𝜇0+𝑢₁(𝑡)−𝛾+𝑑₃000𝑢₁+𝜇+𝑢₂(𝑡)0𝜒(𝑉)=⎡⎢⎢⎢⎣𝑏𝑁(𝑡)(𝑡)−𝛽₁𝛾+𝑢₂(𝑡)−𝑑₄+𝜇+𝜌⎤⎥⎥⎥⎥⎦,𝛽𝑃(𝑡)𝐿(𝑡)₁𝑃(𝑡)𝛽𝐿(𝑡)+𝛽₂𝐿(𝑡)𝑆(𝑡)₂𝐿(𝑡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 ||𝜒𝑉₁−𝜒𝑉₂||≤+ΛΛ₁||₃𝑃₁(𝑡)−𝑃₂(𝑡)||+Λ₂||𝐿₁(𝑡)−𝐿₂(𝑡)||||𝑆₁(𝑡)−𝑆₂(𝑡)||,(13) where the positive constants Λ₁, Λ₂, and Λ₃ are independent of state variables 𝑃(𝑡), 𝐿(𝑡), and 𝑆(𝑡)≤𝑁(𝑡), respectively. Also we get ||𝐺𝑉₁−𝐺𝑉₂||≤Λ||𝑉₁−𝑉₂||,(14) where Λ=Λ₁+Λ₂+Λ₃+‖𝑀‖<∞. 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:𝐻𝑃=,𝐿,𝐿𝑆,,𝑆𝑄,,𝑢𝑢₁₁,𝑢₂,𝜆₁,𝜆₂,𝜆₃,𝜆₄,𝑢₂+𝜆₁(𝑡)𝑔₁(𝑡)+𝜆₂(𝑡)𝑔₂(𝑡)+𝜆₃(𝑡)𝑔₃(𝑡),(15) where (𝐿,𝑆,𝑢₁,𝑢₂)=𝛼₁𝐿(𝑡)+𝛼₂𝑆(𝑡)+(1/2)(𝜉₁𝑢²¹(𝑡)+𝜉₂𝑢²²(𝑡)) 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∈𝑈}𝐽𝑢₁,𝑢₂,(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 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡))∈𝑈, 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)(𝜉₁𝑢²¹(𝑡)+𝜉₂𝑢²²(𝑡)) and is convex on the control set 𝑈. Also we can easily see that there exist a constant 𝜎>1 and positive numbers 𝜂₁ and 𝜂₂ such that 𝐽𝑢₁,𝑢₂≥𝜂₂+𝜂₁||𝑢¹||²+||𝑢₂||²^𝜎/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 𝜆(𝑡)=(𝜆₁(𝑡),𝜆₂(𝑡),…,𝜆_𝑛(𝑡)) 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 𝜆₁(𝑡), 𝜆₂(𝑡), 𝜆₃(𝑡), and 𝜆₄(𝑡) satisfing 𝜆𝜆^1(𝑡)=𝛽₁𝜆₁(𝑡)−𝜆₂(𝑡)𝐿^∗(𝑡)+𝜆₁(𝑡)𝑑₁+𝜇_2,(𝑡)=+𝛽₁𝜆₂𝜆₁(𝑡)−𝜆₂(𝑡)𝑃^∗(𝑡)+𝛽₂𝜆₂(𝑡)−𝜆₃(𝑡)𝑆^∗(𝑡)𝜆_3(𝑡)𝑑₂+𝜇+𝜆₂(𝑡)−𝜆₄(𝑡)𝑢₁(𝑡)−𝛼₁,(𝑡)=+𝛽₂𝜆𝜆₂(𝑡)−𝜆₃(𝑡)𝐿^∗(𝑡)+𝜆₃(𝑡)𝛾+𝑑₃+𝜇₃𝜆′₄(𝑡)−𝜆₄(𝑡)𝑢₂(𝑡)−𝛼₂,(𝑡)=𝜆₄(𝑡)−𝜆₁(𝑡)𝜌+𝜆₄(𝑡)𝑑₄+𝜇,(19) with transversality conditions (or boundary conditions) 𝜆_𝑖(𝑇)=0,𝑖=1,2,3.(20) Furthermore, optimal control pairs are given as follows: 𝑢^∗1(𝑡)=maxmin𝜆¹(𝑡)−𝜆₃(𝑡)𝐿(𝑡)𝜉₁,𝜔₁𝑢,0,(21)_∗2(𝑡)=maxmin𝜆²(𝑡)−𝜆₃(𝑡)𝑆(𝑡)𝜉₂,𝜔₂,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:𝑑𝑃^∗(𝑡)=𝑑𝑏𝑡𝑁_∗𝑑𝐿_∗(𝑡)−𝛽₁𝐿^∗(𝑡)𝑃^∗(𝑡)−𝑑₁+𝜇𝑃^∗(𝑡)+𝜌𝑄^∗(𝑡),(𝑡)=𝑑𝛽𝑡₁−𝐿𝑑^∗(𝑡)𝑃^∗(𝑡)−𝛽₂𝐿^∗(𝑡)𝑆^∗(𝑡)₂+𝜇+maxmin𝜆²(𝑡)−𝜆₄(𝑡)𝐿^∗(𝑡)𝜉₁,𝜔₁×𝐿_∗,0𝑑𝑆_∗(𝑡),(𝑡)=𝑑𝜉𝑡₂−𝐿𝛾^∗(+𝑡)𝑑𝑆^∗₃(𝑡)+𝜇+maxmin𝜆³(𝑡)−𝜆₄(𝑡)𝑆^∗(𝑡)𝜉₂,𝜔₂×𝑆_∗,0𝑑𝑄_∗(𝑡),(𝑡)=𝑑𝑡𝛾+maxmin𝜆₃(𝑡)−𝜆₄(𝑡)𝑆^∗(𝑡)𝜉₂,𝜔₂,0𝑆^∗−𝑑(𝑡)₄+max+𝜇m+in𝜌𝑄^∗𝜆(𝑡₂)(𝑡)−𝜆₄(𝑡)𝐿^∗(𝑡)𝜉₁,𝜔₁,0𝐿^∗(𝑡).(23) With the Hamiltonian 𝐻^∗ at (𝑃^∗,𝐿^∗, 𝑆^∗, 𝑄^∗, 𝑢^∗1, 𝑢^∗2, 𝜆₁, 𝜆₂, 𝜆₃, 𝜆₄),𝐻^∗=+𝛼₁1𝐿^∗(𝑡)+𝛼₂𝑆^∗(𝑡)2⎡⎢⎢⎣𝜉¹maxmin𝜆²(𝑡)−𝜆₄(𝑡)𝐿^∗(𝑡)𝜉₁,𝜔₁,0²+𝜉₂maxmin𝜆³(𝑡)−𝜆₄(𝑡)𝑆^∗(𝑡)𝜉₂,𝜔₂,0²+𝜆₁⎤⎥⎥⎦(𝑡)𝑔^∗1(𝑡)+𝜆₂(𝑡)𝑔^∗2(𝑡)+𝜆₃(𝑡)𝑔^∗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, 𝑢₁(𝑡) and 𝑢₂(𝑡), 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 𝑢₁(𝑡) and 𝑢₂(𝑡), 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 𝑑₄≤𝑑₁≤𝑑₂≤𝑑₃. 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 𝛽₁ and 𝛽₂ 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.

Appendices

A. Proof of Theorem 3

The Jacobian matrix of the giving up smoking model (1) around 𝐸_⋆=(𝑃^⋆,𝐿^⋆,𝑆^⋆,𝑄^⋆) is 𝐽^⋆=⎡⎢⎢⎢⎢⎣𝐴^𝑃𝐵𝐴_𝐿0𝜌_𝑃𝑜𝐵𝐶_𝐿𝐿𝐵_𝑆000𝐶𝛾_𝑆𝐷0_𝑄⎤⎥⎥⎥⎥⎦,(A.1) where 𝐵_𝑃𝐴_𝑃=−𝛽₁𝐿^⋆−𝑑₁+𝜇,𝐴_𝐿=−𝛽₁𝑃^⋆,=𝛽₁𝐿𝐵^⋆,𝐵_𝑃=𝛽₁𝑃^⋆−𝛽₂𝑆^⋆−𝑑₂+𝜇_𝑆,𝐶_𝑆=−𝛽₂𝐿^⋆,𝐶_𝐿=𝛽₂𝑆^⋆,=𝛽₂𝐿^⋆−𝛾+𝑑₃+𝜇,𝐷_𝑄=−𝑑₄+𝜇+𝜌.(A.2) The characteristic equation of 𝐽^⋆ is 𝜆⁴+𝑎₁𝜆³+𝑎₂𝜆²+𝑎₃𝜆+𝑎₄=0,(A.3) where 𝑎𝑎₁₂=−𝐴_𝑃+𝐵_𝐿+𝐶_𝑆+𝐷_𝑄,𝑎₃=𝐶_𝑆𝐷_𝑄+𝐴_𝑃𝐵_𝐶+𝐶_𝑆+𝐷_𝑄𝐴_𝑃+𝐵_𝐿−𝐴_𝐿𝐵_𝑃,𝑎=𝐴_𝐿𝐵_𝑃𝐶_𝑆+𝐷_𝑄−𝐶_𝑆𝐷_𝑄𝐴_𝑃+𝐵_𝐿+𝐴_𝑃𝐵_𝐶𝐶_𝑆+𝐷_𝑄₄,=𝐴_𝑃𝐵_𝐿𝐶_𝑆𝐷_𝑄−𝐴_𝐿𝐵_𝑃𝐶_𝑆𝐷_𝑄−𝜌𝛾𝐵_𝑃𝐶_𝐿.(A.4) By Routh-Hurwitz conditions the equilibrium state is asymptotically stable if 𝑎₁>0, 𝑎₃>0, 𝑎₄>0 and 𝑎₁𝑎₂𝑎₃>𝑎²³+𝑎²¹𝑎₄. 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(𝑡)=−𝜕𝐻𝜕𝑃=𝛽¹𝜆₁𝜆(𝑡)−𝜆₂(𝑡)𝐿(𝑡)+𝜆₁(𝑡)𝑑₁+𝜇_2,(𝑡)=−𝜕𝐻=𝛽₁𝜕𝐿+𝛽₂𝜆₁(𝑡)−𝜆₂(𝑡)𝑃(𝑡)+𝜆𝜆₂(𝑡)−𝜆₃(𝑡)𝑆(𝑡)+𝜆₂(𝑡)𝑑₂+𝜇₂𝜆_3(𝑡)−𝜆₄(𝑡)𝑢₁(𝑡)−𝛼₁,(𝑡)=−𝜕𝐻=𝛽₂𝜕𝑆+𝜆𝜆₂(𝑡)−𝜆₃(𝑡)𝐿(𝑡)+𝜆₃(𝑡)𝛾+𝑑₃+𝜇₃𝜆_4(𝑡)−𝜆₄(𝑡)𝑢₂(𝑡)−𝛼₂,(𝑡)=−𝜕𝐻𝜕𝑄=𝜆⁴(𝑡)−𝜆₁(𝑡)𝜌+𝜆₄(𝑡)𝑑₄+𝜇.(B.1) Using the optimality conditions and the property of the control space 𝑈 for the control variable 𝑢₁, we obtain 𝜕𝐻𝜕𝑢₁=0⟹𝜉₁𝑢₁(𝑡)+−𝜆₂(𝑡)+𝜆₄(𝑡)𝐿(𝑡)=0.(B.2) By using the property of the control space, we obtain 𝑢^∗1(𝑡)=⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩0if𝜆²(𝑡)−𝜆₄(𝑡)𝐿(𝑡)𝜉₁𝜆≤0,₂(𝑡)−𝜆₄(𝑡)𝐿(𝑡)𝜉₁if0<𝜆²(𝑡)−𝜆₄(𝑡)𝐿(𝑡)𝜉₁,<𝜔₁𝜔₁if𝜆²(𝑡)−𝜆₄(𝑡)𝐿(𝑡)𝜉₁≥𝜔₁.(B.3) This can be rewritten in compact notation𝑢^∗1(𝑡)=maxmin𝜆²(𝑡)−𝜆₄(𝑡)𝐿(𝑡)𝜉₁,𝜔₁,0.(B.4) Similar to the control variable 𝑢₂, we get 𝑢^∗2(𝑡)=⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩0if𝜆³(𝑡)−𝜆₄(𝑡)𝑆(𝑡)𝜉₂𝜆≤0,₃(𝑡)−𝜆₄(𝑡)𝑆(𝑡)𝜉₂if0<𝜆³(𝑡)−𝜆₄(𝑡)𝑆(𝑡)𝜉₂𝜔₂<𝜔₂,if𝜆³(𝑡)−𝜆₄(𝑡)𝑆(𝑡)𝜉₂≥𝜔₂.(B.5) This can also be rewritten similar to the above compact form as𝑢^∗2(𝑡)=maxmin𝜆³(𝑡)−𝜆₄(𝑡)𝑆(𝑡)𝜉₂,𝜔₂,0.(B.6) This completes the proof.

Acknowledgment

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