# Optimal Control Strategies in an Alcoholism Model

Authors: Xun-Yang Wang; Hai-Feng Huo; Qing-Kai Kong; Wei-Xuan Shi
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper presents a deterministic SATQ-type mathematical model (including susceptible, alcoholism, treating, and quitting compartments) for the spread of alcoholism with two control strategies to gain insights into this increasingly concerned about health and social phenomenon. Some properties of the solutions to the model including positivity, existence and stability are analyzed. The optimal control strategies are derived by proposing an objective functional and using Pontryagin’s Maximum Principle. Numerical simulations are also conducted in the analytic results.


---

## Body

## 1. Introduction

Alcoholism, also known as alcohol dependence, is a disease that includes the desire for alcohol and continuing to drink it despite its negative effect on individual’s health, relationships, and social status [1]. Similar to all other drug addictions, alcoholism can be regarded as a treatable disease. The World Health Organization estimates that about 140 million people throughout the world suffer from alcohol dependence with related problems, such as being sick, losing a job, among a host of other things [2]. Particularly, young people’s alcoholism problem is a major concern to public health. US surveys indicate that approximately 90% of college students have consumed alcohol at least once [3], and more than 40% of college students have engaged in binge drinking [4, 5]. Unfortunately, the biological mechanisms underpinning alcoholism are not known; however, risk factors include social environment, stress, mental health, genetic sensitivity, age, ethnic group, and sex [6, 7]. Long-term alcohol abuse will produce negative changes in the brain such as tolerance and physical dependence. The subtle changes make it difficult for the alcoholics to stop drinking and result in alcohol withdrawal symptoms upon discontinuation of alcohol consumption. Alcohol damages almost all parts of the body and contribute to a number of human diseases including but not limited to liver cirrhosis, pancreatitis, heart disease, and sexual dysfunction and can eventually be deadly [8]. Damage to the central and peripheral nervous systems can take place from sustained alcohol consumption [9–13].Although alcoholism is becoming more and more dangerous and serious as well as a widespread social phenomenon, only much less work has been done in the mathematical modelling of alcoholism as a growing health problem, including a few studies which offered some mathematical approaches to understand the growing burden of alcoholism [10, 14–19]. In [10], a SIR-type model was proposed; the authors used standard contact rate between susceptibles and alcoholism, getting alcoholism reproductive number and discussing the existence and stability of two equilibria. In [14], a framework where drinking was modeled as a socially contagious process in low- and high-risk connected environments was introduced; they found that high levels of social interaction between light and moderate drinkers in low-risk environments can diminish the importance of the distribution of relative drinking times on the prevalence of heavy drinking. In [15], neurophysiological examinations of 100 long-term alcohol dependent patients, who were having neuropsychiatric treatment, showed symptoms of polytopic damage of the peripheral and central nervous system. The results showed that for recognition of the damage an extensive diagnostic programme must be used. In [16], the authors considered a kind of binge drinking model with two equal infectivity drunk states; mathematical analyses established that the global dynamics of the model were determined by the basic reproduction number. In [17], the authors modified the model from [16]; that is, they considered different infectivity of two drunk states, and a SEIR-type model of alcoholism was thus presented, in which two alcohol related states were involved, namely, no alcohol dependent consumers D(t) and alcohol dependent consumers A(t). In [18], the authors formulated a deterministic model for evaluating the impact of heavy alcohol drinking on the reemerging gonorrhea epidemic, and both analytical and numerical results were provided to ascertain whether heavy alcohol drinking had an impact on the transmission dynamics of gonorrhea. The approach of the literature [18] was very meaningful, since it provided a new direction of thinking when the cross-infection between alcoholism and other pathological diseases occurs. In recent monograph [19], the authors also proposed a SIR-type model to investigate alcohol abuse phenomenon and generated some useful insights; for example, the basic reproductive number was not always the key to controlling drinking within the population. For other papers that study the model of giving up smoking or quitting drinking, please see [20, 21] and references cited therein.As living standard and health awareness get improved, more and more people who fall into binge drinking state are actively seeking the quitting alcoholism measures and treatment methods [1, 11, 22]. In [22], treatment strategy was introduced into a simple SIR-type alcoholics quitting model, in which the authors used the bilinear incidence to depict the “infection" between the occasional drinkers S and problem drinkers D. Motivated by some aforementioned documents [10, 19, 22], in this paper, we will formulate a more reasonable alcoholics quitting model. The fact that our model is more reasonable is embodied from the following three aspects.(1)
Taking into account that alcoholism is a widespread social phenomenon, so the standard incidence is superior to bilinear incidence when we portray the relationship between the alcoholism and the susceptibles during the course of infection. While in [22], the authors adopted bilinear incidence, we will adopt standard incidence in this paper.(2)
Since alcohol is harmful to health, moreover, as we all know, alcoholism is treatable if we can take approximate measure in time, for example, artificial isolation from alcoholisms, medications, persuasion, and education programing on alcoholism. So it is necessary to take effective measures to avoid alcohol or to treat after alcoholism. Documents [10, 19] have not considered these aspects.(3)
Since there is effective prevention and treatment in describing the phenomenon of alcoholism, there are some people who will never drink due to successful prevention or some people who no longer drink after successful treatment. Therefore, when we formulate the model in this paper, it's reasonable to introduce a new compartmentQ, the people in which will never drink for ever. Obviously, the models of [10, 19, 22] are not involving the quitting compartment Q.Based on the above considerations, we will premeditate two treating methods, namely, prevention of susceptibles from alcoholism and treatment on alcoholism as control variables; hence, we will derive a SATQ-type model. We note the fact that many authors are interested in solving optimal control problems, such as cost minimization and optimal control of various disease, especially with biological background and various mathematical models [22–24]. In this paper, we will propose an objective functional which considers not only alcohol quitting effects but also the cost of controlling alcohol. Then, we consider a range of issues related to the optimal control with the method of Pontryagin’s Maximum Principle, including optimal control existence, uniqueness, and characterization.The organization of this paper is as follows. In the next section, the alcoholism model with prevention for the susceptibles and treatment for alcoholism is formulated. In Section3, the basic reproduction number and the existence of equilibria are investigated. The stability of the disease free and endemic equilibria is proved in Section 4. Optimal control strategies by the classic method of PMP (Pontryagin’s Maximum Principle) are discussed in Section 5. In Section 6, we give some numerical simulations. We give some discussions and conclusions in the last section.


## 2. The Model Formulation and Some Fundamental Properties

In this section, we introduce a mathematical model with prevention and treatment for the alcoholism and then study some important properties such as the boundness and positivity of its solutions.


### 2.1. Model Formulation and Parameter Explanation

The total population is partitioned into four compartments: the susceptible compartmentS which refers to the persons who never drink or drink moderately without affecting the physical health, the alcoholism compartment A which refers to the persons who binge drink and affect the physical health seriously, the treatment compartment T which refers to the persons who have been receiving treatments by taking pills or other medical interventions after alcoholism, and the quitting compartment Q which refers to the persons who recover from alcoholism after treatment and stay off alcohol hereafter. In this paper, we focus on a closed environment, such as a community, a university, or a village. So the total number of population to be considered is a constant; we denote it as N. The population flow among those compartments is shown in the following diagram (Figure 1).Figure 1
Transfer diagram for the dynamics of alcoholism model.The schematic diagram leads to the following system of ordinary differential equations:(1)S′=μN-(1-u1)βSAN-μS,A′=(1-u1)βSAN+ξT-(u2+μ)A,T′=u2A-(μ+ξ+δ)T,Q′=δT-μQ.
Here, μN is the birth number of the population; μ is the natural death rate of the population; u1 is the fraction of the susceptible individuals who successfully avoid to stay off the alcoholism; u2 is the fraction of the alcoholics who take part in treatments; here, 0≤ui≤1, i=1,2, and they will be considered as two control variables in Section 5; β is the transmission coefficient of the “infection” for the susceptible individuals from the alcoholic individuals; ξ is the rate coefficient of the person who fail to be treated and return to the alcoholism compartment mostly due to their own weak will; δ is the rate coefficient of the person who have received effective treatment and recovered from alcoholism forever.


### 2.2. Boundedness of Solutions to System and Positively Invariant Region

It is important to show positivity and boundedness for the system (1) as they represent populations. Firstly, we present the positivity of the solutions. System (1) can be put into the matrix form
(2)X′=G(X),
where X=(S,A,T,Q)T∈R4 and G(X) is given by
(3)G(X)=(G1(X)G2(X)G3(X)G4(X))=(μN-(1-u1)βSAN-μS(1-u1)βSAN+ξT-(u2+μ)Au2A-(μ+ξ+δ)TδT-μQ).
It is easy to check that
(4)Gi(X)|Xi(t)=0,Xt∈C+≥0,i=1,2,3,4.
Due to Lemma  2 in [25], any solution of (1) is  X(t)∈R+4 for all t≥0.We denoteN(t)=S(t)+A(t)+T(t)+Q(t); summing equations in (1) yields
(5)dN(t)dt=0,
so N(t)=S(t)+A(t)+T(t)+Q(t)=constant (denoted as N), and the set
(6)Ω={(S,A,T,Q)∈R+4:S+A+T+Q≤N}
is a positively invariant region for (1). Therefore, we will consider the global stability of (1) on the set Ω.


## 2.1. Model Formulation and Parameter Explanation

The total population is partitioned into four compartments: the susceptible compartmentS which refers to the persons who never drink or drink moderately without affecting the physical health, the alcoholism compartment A which refers to the persons who binge drink and affect the physical health seriously, the treatment compartment T which refers to the persons who have been receiving treatments by taking pills or other medical interventions after alcoholism, and the quitting compartment Q which refers to the persons who recover from alcoholism after treatment and stay off alcohol hereafter. In this paper, we focus on a closed environment, such as a community, a university, or a village. So the total number of population to be considered is a constant; we denote it as N. The population flow among those compartments is shown in the following diagram (Figure 1).Figure 1
Transfer diagram for the dynamics of alcoholism model.The schematic diagram leads to the following system of ordinary differential equations:(1)S′=μN-(1-u1)βSAN-μS,A′=(1-u1)βSAN+ξT-(u2+μ)A,T′=u2A-(μ+ξ+δ)T,Q′=δT-μQ.
Here, μN is the birth number of the population; μ is the natural death rate of the population; u1 is the fraction of the susceptible individuals who successfully avoid to stay off the alcoholism; u2 is the fraction of the alcoholics who take part in treatments; here, 0≤ui≤1, i=1,2, and they will be considered as two control variables in Section 5; β is the transmission coefficient of the “infection” for the susceptible individuals from the alcoholic individuals; ξ is the rate coefficient of the person who fail to be treated and return to the alcoholism compartment mostly due to their own weak will; δ is the rate coefficient of the person who have received effective treatment and recovered from alcoholism forever.


## 2.2. Boundedness of Solutions to System and Positively Invariant Region

It is important to show positivity and boundedness for the system (1) as they represent populations. Firstly, we present the positivity of the solutions. System (1) can be put into the matrix form
(2)X′=G(X),
where X=(S,A,T,Q)T∈R4 and G(X) is given by
(3)G(X)=(G1(X)G2(X)G3(X)G4(X))=(μN-(1-u1)βSAN-μS(1-u1)βSAN+ξT-(u2+μ)Au2A-(μ+ξ+δ)TδT-μQ).
It is easy to check that
(4)Gi(X)|Xi(t)=0,Xt∈C+≥0,i=1,2,3,4.
Due to Lemma  2 in [25], any solution of (1) is  X(t)∈R+4 for all t≥0.We denoteN(t)=S(t)+A(t)+T(t)+Q(t); summing equations in (1) yields
(5)dN(t)dt=0,
so N(t)=S(t)+A(t)+T(t)+Q(t)=constant (denoted as N), and the set
(6)Ω={(S,A,T,Q)∈R+4:S+A+T+Q≤N}
is a positively invariant region for (1). Therefore, we will consider the global stability of (1) on the set Ω.


## 3. The Basic Reproduction Number and Existence of Alcoholism Equilibria

### 3.1. The Basic Reproduction NumberR0

In epidemiology, the basic reproduction number (sometimes called basic reproductive rate or basic reproductive ratio) of an infection is the number of infectious cases that one infectious case generates on average over the course of its infectious period. While in this context, it means the number of persons that an alcoholic will “infect” during his “infectious” period in the pure susceptible environment so that the infected persons will enter the alcoholism compartment. It is easy to see that the model has an alcohol free equilibriumE0=(S0,0,0,0)=(N,0,0,0). In the following, the basic reproduction number of system (1) will be obtained by the next generation matrix method formulated in [26].Letx=(A,T,Q,S)T, then system (1) can be written as
(7)dxdt=ℱ(x)-𝒱(x),
where
(8)ℱ(x)=((1-u1)βSAN000),𝒱(x)=((u2+μ)A-ξT(μ+ξ+δ)T-u2AμQ-δT(1-u1)βSAN+μS-μN).
The Jacobian matrices of ℱ(x) and 𝒱(x) at the alcohol free equilibrium E0 are, respectively,
(9)Dℱ(E0)=(F3×3000),D𝒱(E0)=(V3×30(1-u1)β00μ),
where
(10)F=((1-u1)β00000000),V=(u2+μ-ξ0-u2μ+ξ+δ00-δμ).
The basic reproduction number, denoted by R0, is thus given by
(11)R0=ρ(FV-1)=β(1-u1)(μ+ξ+δ)u2(μ+δ)+μ(μ+ξ+δ).
It is easy to see that both of the control parameters contributed to reducing the alcoholism. From this point, the control measures are meaningful.


### 3.2. Existence of Alcoholism Equilibrium

The endemic equilibriumE*(S*,A*,T*,Q*) of system (1) is determined by equations
(12)μN-(1-u1)βSAN-μS=0,(1-u1)βSAN+ξT-(u2+μ)A=0,u2A-(μ+ξ+δ)T=0,δT-μQ=0.
The third equation in (12) leads to
(13)A=μ+ξ+δu2T.
From the last equation in (12), we have
(14)Q=δTμ.
From the first equation of (12), and together with (13), we can get
(15)S=μNμ+(((1-u1)β)/N)A=μN2u2μNu2+(1-u1)βT(μ+ξ+δ).
Substituting (13)–(15) into the second equation of (12) gives
(16)μN-μ2N2u2μNu2+(1-u1)βT(μ+ξ+δ)+ξT-(μ+u2)μ+ξ+δu2T=0.
By simplifying (16), we can get
(17)T{[(μ+ξ+δ)2u2ξ(1-u1)β(μ+ξ+δ)-(1-u1)β(μ+ξ+δ)2(μ+u2)]T+σ}=0,
where
(18)σ=u2μN(1-u1)β(μ+ξ+δ)+(u2)2μNξ-μNu2(μ+u2)(μ+ξ+δ).Hence, we get two explicit solutions to (17); one is T0=0, which is corresponding to the alcohol free equilibria, and the other is
(19)T*=(σ)((1-u1)β(μ+ξ+δ)2(μ+u2)-u2ξ(1-u1)β(μ+ξ+δ)(μ+ξ+δ)2)-1=σ(1-u1)β(μ+ξ+δ)[μξ+(δ+μ)(μ+u2)],
which should be corresponding to the alcoholism equilibria on condition that T*>0; otherwise, the alcoholism equilibria are nonexistent. It is enough to show the positivity of σ to make sure the existence of alcoholism equilibria on the condition R0≥1. By some simple calculations, we simplify the expression of σ to be
(20)σ=μNu2{(μ+ξ+δ)(1-u1)βllllllllllllll-[(μ+u2)(δ+μ)+μξ]}.
Since R0>1 is equivalent to
(21)β(1-u1)(μ+ξ+δ)>u2(μ+δ)+μ(μ+ξ+δ),
the right side of this inequality is exactly equal to (μ+u2)(δ+μ)+μξ. Hence, we have proved the existence of T*>0, so are the alcoholism equilibria. We summarize this result in Theorem 1.Theorem 1.
For system (1), there is always an alcohol free equilibrium E0=(N,0,0,0). When R0>1, besides alcohol free equilibrium E0, system (1) also has a unique alcoholism equilibrium E*(S*,A*,T*,Q*), where
(22)S*=μN2u2μNu2+(1-u1)βT*(μ+ξ+δ),A*=μ+ξ+δu2T*,Q*=δT*μ,T*=σ(1-u1)β(μ+ξ+δ)[μξ+(δ+μ)(μ+u2)].


## 3.1. The Basic Reproduction NumberR0

In epidemiology, the basic reproduction number (sometimes called basic reproductive rate or basic reproductive ratio) of an infection is the number of infectious cases that one infectious case generates on average over the course of its infectious period. While in this context, it means the number of persons that an alcoholic will “infect” during his “infectious” period in the pure susceptible environment so that the infected persons will enter the alcoholism compartment. It is easy to see that the model has an alcohol free equilibriumE0=(S0,0,0,0)=(N,0,0,0). In the following, the basic reproduction number of system (1) will be obtained by the next generation matrix method formulated in [26].Letx=(A,T,Q,S)T, then system (1) can be written as
(7)dxdt=ℱ(x)-𝒱(x),
where
(8)ℱ(x)=((1-u1)βSAN000),𝒱(x)=((u2+μ)A-ξT(μ+ξ+δ)T-u2AμQ-δT(1-u1)βSAN+μS-μN).
The Jacobian matrices of ℱ(x) and 𝒱(x) at the alcohol free equilibrium E0 are, respectively,
(9)Dℱ(E0)=(F3×3000),D𝒱(E0)=(V3×30(1-u1)β00μ),
where
(10)F=((1-u1)β00000000),V=(u2+μ-ξ0-u2μ+ξ+δ00-δμ).
The basic reproduction number, denoted by R0, is thus given by
(11)R0=ρ(FV-1)=β(1-u1)(μ+ξ+δ)u2(μ+δ)+μ(μ+ξ+δ).
It is easy to see that both of the control parameters contributed to reducing the alcoholism. From this point, the control measures are meaningful.


## 3.2. Existence of Alcoholism Equilibrium

The endemic equilibriumE*(S*,A*,T*,Q*) of system (1) is determined by equations
(12)μN-(1-u1)βSAN-μS=0,(1-u1)βSAN+ξT-(u2+μ)A=0,u2A-(μ+ξ+δ)T=0,δT-μQ=0.
The third equation in (12) leads to
(13)A=μ+ξ+δu2T.
From the last equation in (12), we have
(14)Q=δTμ.
From the first equation of (12), and together with (13), we can get
(15)S=μNμ+(((1-u1)β)/N)A=μN2u2μNu2+(1-u1)βT(μ+ξ+δ).
Substituting (13)–(15) into the second equation of (12) gives
(16)μN-μ2N2u2μNu2+(1-u1)βT(μ+ξ+δ)+ξT-(μ+u2)μ+ξ+δu2T=0.
By simplifying (16), we can get
(17)T{[(μ+ξ+δ)2u2ξ(1-u1)β(μ+ξ+δ)-(1-u1)β(μ+ξ+δ)2(μ+u2)]T+σ}=0,
where
(18)σ=u2μN(1-u1)β(μ+ξ+δ)+(u2)2μNξ-μNu2(μ+u2)(μ+ξ+δ).Hence, we get two explicit solutions to (17); one is T0=0, which is corresponding to the alcohol free equilibria, and the other is
(19)T*=(σ)((1-u1)β(μ+ξ+δ)2(μ+u2)-u2ξ(1-u1)β(μ+ξ+δ)(μ+ξ+δ)2)-1=σ(1-u1)β(μ+ξ+δ)[μξ+(δ+μ)(μ+u2)],
which should be corresponding to the alcoholism equilibria on condition that T*>0; otherwise, the alcoholism equilibria are nonexistent. It is enough to show the positivity of σ to make sure the existence of alcoholism equilibria on the condition R0≥1. By some simple calculations, we simplify the expression of σ to be
(20)σ=μNu2{(μ+ξ+δ)(1-u1)βllllllllllllll-[(μ+u2)(δ+μ)+μξ]}.
Since R0>1 is equivalent to
(21)β(1-u1)(μ+ξ+δ)>u2(μ+δ)+μ(μ+ξ+δ),
the right side of this inequality is exactly equal to (μ+u2)(δ+μ)+μξ. Hence, we have proved the existence of T*>0, so are the alcoholism equilibria. We summarize this result in Theorem 1.Theorem 1.
For system (1), there is always an alcohol free equilibrium E0=(N,0,0,0). When R0>1, besides alcohol free equilibrium E0, system (1) also has a unique alcoholism equilibrium E*(S*,A*,T*,Q*), where
(22)S*=μN2u2μNu2+(1-u1)βT*(μ+ξ+δ),A*=μ+ξ+δu2T*,Q*=δT*μ,T*=σ(1-u1)β(μ+ξ+δ)[μξ+(δ+μ)(μ+u2)].


## 4. Stability Analysis of Equilibria

For the convenience of subsequent proof, we denote a vectorX=(A,T,Q,S)T and
(23)f(X)=((1-u1)βSAN+ξT-(u2+μ)Au2A-(μ+ξ+δ)TδT-μQμN-(1-u1)βSAN-μS).
So the Jacobian matrix of f(x) about vector X is as the following:(24)J=∂f(X)∂X=((1-u1)βSN-(μ+u2)ξ0(1-u1)βANu2-(μ+ξ+δ)000δ-μ0-(1-u1)βSN00-μ-(1-u1)βAN).Theorem 2.
For system (1), the alcohol free equilibrium E0 is locally asymptotically stable if R0<1.Proof.
Since(25)J(E0)=((1-u1)β-(μ+u2)ξ00u2-(μ+ξ+δ)000δ-μ0-(1-u1)β00-μ),
we can easily get that two of the eigenvalues are λ1=λ2=-μ<0, while λ3,λ4 satisfy
(26)λ2+[2μ+ξ+δ+u2-(1-u1)β]λ+(μ+ξ+δ)(μ+u2-(1-u1)β)-u2ξ=0.
Thus,
(27)λ3+λ4=(1-u1)β-(μ+ξ+δ)-(μ+u2).
Since R0<1 is equivalent to
(28)β(1-u1)(μ+ξ+δ)<u2(μ+δ)+μ(μ+ξ+δ)<(μ+u2)(μ+ξ+δ),
so
(29)β(1-u1)<μ+u2,
and then
(30)λ3+λ4<-(μ+ξ+δ)<0,
while
(31)λ3λ4=(μ+ξ+δ)(μ+u2-(1-u1)β)-u2ξ.
Similarly from R0<1, we can derive the inequality
(32)-β(1-u1)(μ+ξ+δ)>-u2(μ+δ)-μ(μ+ξ+δ),
so
(33)λ3λ4>(μ+u2)(μ+ξ+δ)-u2ξ-u2(μ+δ)-μ(μ+ξ+δ).
It reduces to
(34)λ3λ4>0.
Hence, Reλ3<0, Reλ4<0. The proof is complete.Next, we will turn to investigate the global stability ofE0.Theorem 3.
For system (1), the alcohol free equilibrium E0 is globally asymptotically stable if R0<1.Proof.
Consider the subsystem of (1) as follows:
(35)A′=(1-u1)βSAN+ξT-(u2+μ)A,T′=u2A-(μ+ξ+δ)T,Q′=δT-μQ.
Equation (35) can be rewritten as
(36)(A˙T˙Q˙)=(F-V)(ATQ)-(1-SN)×(β(1-u1)00000000)(ATQ).
Since S≤N and 0≤u1≤1, then for all t>0, we can get
(37)(A˙T˙Q˙)≤(F-V)(ATQ).
According to Lemma 1 in [26], all the eigenvalues of matrix F-V have negative real parts, so the solutions of this subsystem are stable whenever R0<1. So (A(t),T(t),Q(t))→(0,0,0) as t→∞. By the comparison theorem [27], and based on the fact that the total population is constant N, it follows that (A(t),T(t),Q(t))→(0,0,0) and S(t)→N as t→∞. So the alcohol free equilibrium E0 is globally asymptotically stable; the proof is complete.Theorem 4.
For system (1), the alcoholism equilibrium E*(S*,A*,T*,Q*) is globally asymptotically stable if R0>1.Proof.
Since the total population in model (1) is a constant number N, in order to prove the global stability of system (1), it is sufficed to prove the corresponding stability of subsystem (38):
(38)S′=μN-(1-u1)βSAN-μS,A′=(1-u1)βSAN+ξT-(u2+μ)A,T′=u2A-(μ+ξ+δ)T.
We make normalization transform and still use the same symbols S,A,T to denote the variables; then (38) can be transformed into
(39)S′=μ-(1-u1)βSA-μS,sA′=(1-u1)βSA+ξT-(u2+μ)A,T′=u2A-(μ+ξ+δ)T.
From (39), we can easily know that the equilibria (S*,A*,T*) satisfy the following three equalities to be used later:
(40)μ=(1-u1)βS*A*-μS*,(1-u1)βS*A*+ξT*=(u2+μ)A*,u2A*=(μ+ξ+δ)T*.
LetV=x1(S-S*-S*ln⁡⁡(S/S*))+x2(A-A*-A*ln⁡⁡(A/A*))+x3(T-T*-T*ln⁡⁡(T/T*)); then
(41)V′|(39)=x1[S*Sμ-(1-u1)βSA-μSlllllllllll-S*Sμ+(1-u1)βS*A+μS*]+x2[A*A(1-u1)βSA+ξT-(u2+μ)Allllllllllllll-(1-u1)βSA*-A*AξT+(u2+μ)A*]+x3[T*Tu2A-(μ+ξ+δ)Tlllllllllllll-T*Tu2A+(μ+ξ+δ)T*]=x1[S*S(1-u1)βS*A*+μS*-(1-u1)βSA-μSllllllllll-S*S((1-u1)βS*A*+μS*)llllllllll+(1-u1)βS*A+μS*S*S]+x2[A*A(1-u1)βSA+ξT-(u2+μ)Alllllllllllll-(1-u1)βSA*-A*AξT+(u2+μ)A*]+x3[T*Tu2A-(μ+ξ+δ)Tlllllllllllll-T*Tu2A+(μ+ξ+δ)T*]=x1μS*(2-SS*-S*S)+[x1(1-u1)βS*A*+x2(u2+μ)A*llllllllll+x3(μ+ξ+δ)T*]-[x1(S*)2(1-u1)βA*S+x2(1-u1)βSA*llllllllll+x2A*ξTA+x3T*Tu2A(S*)2(1-u1)βA*S]+SA[-x1(1-u1)β+x2(1-u1)β]+A[-(u2+μ)x2+x1(1-u1)βS*+u2x3]+T[ξx2-(μ+ξ+δ)x3].
To eliminate the cross-term SA and two single-variable terms A and T, we let
(42)-x1(1-u1)β+x2(1-u1)β=0,-(u2+μ)x2+x1(1-u1)βS*+u2x3=0,ξx2-(μ+ξ+δ)x3=0.
By solving them, we can get
(43)x1=1,x2=1,x3=ξμ+ξ+δ=u2+μ-(1-u1)βS*u2.
Next, we let
(44)V1′=x1(1-u1)βS*A*+x2(u2+μ)A*+x3(μ+ξ+δ)T*,V2′=-[x1(S*)2(1-u1)βA*Slllllllllll+x2(1-u1)βSA*+x2A*ξTA+x3T*Tu2A(S*)2(1-u1)βA*S],
and then
(45)V′=x1μS*(2-SS*-S*S)+V1′+V2′.
Due to
(46)(u2+μ)A*=(1-u1)βS*A*+ξT*,x3(μ+ξ+δ)T*=ξT*x2=ξT*,
so
(47)V1′=2(1-u1)βS*A*+2ξT*,V2′=-[(S*)2(1-u1)βA*S+(1-u1)βSA*+A*ξTA+x3T*Tu2A(S*)2(1-u1)βA*S]≤-2[(S*)2(1-u1)βA*S·(1-u1)βSA*]1/2-2[A*ξTA·x3T*Tu2A]=-2(1-u1)βS*A*-2(x3A*ξu2T*)1/2=-2(1-u1)βS*A*-2[ξT*(x3u2A*)]1/2=-2(1-u1)βS*A*-2{ξT*[(u2+μ)A*-(1-u1)βS*A*]}1/2=-2(1-u1)βS*A*-2(ξT*ξT*)1/2=-2(1-u1)βS*A*-2ξT*.
Hence,
(48)V′=μS*(2-SS*-S*S)+V1′+V2′≤μS*(2-SS*-S*S)+2(1-u1)βS*A*+2ξT*-2(1-u1)βS*A*-2ξT*=μS*(2-SS*-S*S)≤0,V′=0 if and only if S=S*, A=A*, T=T*. According to LaSalle’s invariance principle [28], we can derive the conclusion that the alcoholism equilibria E*(S*,A*,T*,Q*) are globally asymptotically stable; the proof is complete.


## 5. Optimal Control Problem

### 5.1. The Existence of Optimal Control

In order to investigate an effective campaign to control alcoholism in a community which pursue the goals of the minimized alcoholisms and more recovered individuals, we will reconsider the system (1) and use two control variables to reduce the numbers of alcoholics. The difference is that we will change the parameters u1,u2 into control variable u1(t),u2(t). Their aforementioned definitions allow us to do so. u1(t) is used to limit the proportion of the susceptible individual to contact with alcoholism, usually by propaganda and education, so that the susceptible individual can stay off alcoholism consciously and be free of “infection,” we can understand the effect of u1(t) is to prevent the the susceptible from contacting with the alcoholism. The control variable u2(t) is used to control the alcoholism to take appropriate treatment measures, such as taking pills or seeking other medical help. However, just as a coin has two sides, there will be a lot of costs generated during the control process. So it is advisable to balance between the costs and the alcohol effects. In view of this, our optimal control problem to minimize the objective functional is given by
(49)J(u1,u2)=∫0tf‍[A(t)+c12u12(t)+c22u22(t)]dt,
which subjects to system
(50)S′=μN-(1-u1(t))βSAN-μS,A′=(1-u1(t))βSAN+ξT-(u2(t)+μ)A,T′=u2(t)A-(μ+ξ+δ)T,Q′=δT-μQ,
with initial conditions
(51)S(0)=S0,A(0)=A0,T(0)=T0,Q(0)=Q0.
Here, ui(t)∈U≜{(u1,u2)∣ui(t)ismeasurableand0≤ui(t)≤1,forallt∈[0,tf]}, tf is the end time to be controlled, U is an admissible control set, ci, and i=1,2, are weight factors (positive constants) that adjust the intensity of two different control measures.Next, we will investigate the existence of the optimal control of the above-mentioned problem.Theorem 5.
There exists an optimal control pairu*=(u1*,u2*)∈U such that
(52)J(u1*,u2*)=min⁡J(u1,u2),u1(t),u2(t)∈U
subjects to the control system (1) with initial conditions (50).Proof.
To prove the existence of an optimal control, according to the classic literature [29], we have to show the following.(1)
The control and state variables are nonnegative values.(2)
The control setU is convex and closed.(3)
The right side of the state system is bounded by linear function in the state and control variables.(4)
The integrand of the objective functional is concave onU.(5)
There exist constantsd1,d2>0 and α>1 such that the integrand L(t;u1;u2)≜A(t)+(c1/2)u12(t)+(c2/2)u22(t) of the objective functional satisfies
(53)L(t;u1;u2)≥d1(|u1|2+|u2|2)α/2-d2statements (1), (2) and (3) are obvious satisfied, we only need to test and verify the latter two ones. Since the four state variables have been all proved to be up bounded by N, we will get the following equalities:
(54)S′≤μN,A′≤(1-u1(t))βS+ξT,T′≤u2(t)A,Q′≤δT,
so the fourth condition is set up. As for the last condition,
(55)L(t;u1;u2)≥d1(|u1|2+|u2|2)α/2-d2
is also true, when we choose d1=min⁡{c1/2,c2/2}, and for all d2∈R+,α=2. The proof is complete.We next come to the core of this section.


### 5.2. The Characterization of the Optimal Control

With the existence of the optimal control pairs established, we now present the optimality system and use a result from [30]; we can easily know the existence of the solutions to the optimality system (71) which will be gotten later. Firstly, we come to discuss the theorem that relates to the characterization of the optimal control. The optimality system can be used to compute candidates for optimal control pairs. To do this, we begin by defining an augmented Hamiltonian H with penalty terms for the control constraints as follows:
(56)H=A(t)+c12u12(t)+c22u22(t)+λ1[μN-(1-u1(t))βSAN-μS]+λ2[(1-u1(t))βSAN+ξT-(u2(t)+μ)A]+λ3[u2(t)A-(μ+ξ+δ)T]+λ4(δT-μQ)-w11u1(t)-w12(1-u1(t))-w21u2(t)-w22(1-u2(t)),
where wij(t)≥0 are the penalty multipliers satisfying
(57)w11(t)u1(t)=w12(t)(1-u1(t))=0atoptimalcontrolu1*,w21(t)u2(t)=w22(t)(1-u2(t))=0atoptimalcontrolu2*.Theorem 6.
Given optimal control pairs(u1*,u2*) and solutions S(t),A(t),T(t),Q(t) of the corresponding state system (50), there exist adjoint variables λi, i=1,2,3,4, satisfying
(58)λ1′=λ1(1-u1(t))βAN+μλ1-λ2(1-u1(t))βAN,λ2′=-1+λ1(1-u1(t))βSN-λ2(1-u1(t))βSN+λ2(μ+u2(t))-λ3u2(t),λ3′=-λ2ξ+λ3(μ+ξ+δ)-λ4δ,λ4′=μλ4,
with the terminal conditions
(59)λi(tf)=0,i=1,2,3,4.
Furthermore, (u1*,u2*) are represented by
(60)u1*=min⁡(1,max⁡(0,βSA(λ2-λ1)c1N)),u2*=min⁡(1,max⁡(0,A(λ2-λ3)c2)).Proof.
According to Pontryagin Maximum Principle [29–31], we first differentiate the Hamiltonian operator H, with respect to states. Then the adjoint system can be written as
(61)λ1′=-∂H∂S,λ2′=-∂H∂A,λ3′=-∂H∂T,λ4′=-∂H∂Q.
The terminal condition (56) of adjoint equations is given by λi(tf)=0, i=1,2,3,4.
To obtain the necessary conditions of optimality (59), we also differentiate the Hamiltonian operator H, with respect to U=(u1,u2) and set them equal to zero; then
(62)∂H∂u1=c1u1(t)+λ1βSAN-λ2βSAN-w11+w12=0,∂H∂u2=c2u2(t)-λ2A+λ3A-w21+w22=0.
By solving the optimal control, we obtain(63)u1*=1c1[(λ2-λ1)βSAN+w11-w12].
To determine an explicit expression for the optimal control withoutw11 and w12, a standard optimality technique is utilized [29]. We consider the following three cases.(i)
On the set{t∣0<u1*(t)<1}, we have w11(t)=w12(t)=0. Hence, the optimal control is
(64)u1*=1c1[(λ2-λ1)βSAN].(ii)
On the set{t∣u1*(t)=1}, we have w11(t)=0. Hence,
(65)1=u1*(t)=1c1[(λ2-λ1)βSAN-w12].This implies that(66)1c1(λ2-λ1)βSAN≥1sincew12(t)≥0.(iii)
On the set{t∣u1*(t)=0}, we have w12(t)=0. Hence,
(67)0=u1*(t)=1c1[(λ2-λ1)βSAN+w11].This implies that(68)1c1(λ2-λ1)βSAN≤0sincew11(t)≥0.
Combining these results, the optimal controlu1*(t) is characterized as
(69)u1*=min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}.
Using the similar arguments, we can also obtain the other optimal control function
(70)u2*=min⁡{1,max⁡{0,(λ2-λ3)Ac2}}.
The proof is complete.We point out that the optimality system consists of the state system (50) with the initial conditions S(0),A(0),T(0),Q(0), the adjoint (or costate) system (58) with the terminal conditions (59), and the optimality condition (60). Any optimal control pairs must satisfy this optimality system. For the convenience of subsequent numerical simulation in Section 6, we give the optimality system as follows:
(71)S′=μN-(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βSAN-μS,A′=(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))βSAN+ξT-(min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)+μ)A,T′=min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)A-(μ+ξ+δ)T,Q′=δT-μQ,λ1′=λ1(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βAN+μλ1-λ2(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βAN,λ2′=-1+λ1(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βSN-λ3·min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)-λ2(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}})×βSN+λ2(μ+min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)),λ3′=-λ2ξ+λ3(μ+ξ+δ)-λ4δ,λ4′=μλ4,S(0)=S0,A(0)=A0,T(0)=T0,Q(0)=Q0,λi(tf)=0,i=1,2,3,4.


### 5.3. The Uniqueness of Optimal Control

Due to the a priori boundedness of the state, adjoint functions, and the resulting Lipschitz structure of the ODEs, we can obtain the uniqueness of the optimal control.Lemma 7 (see [23]).
The functionu*(s)=min⁡(b,max⁡(s,a)) is Lipschitz continuous in s, where a<b are some fixed positive constants.Theorem 8.
For allt∈[0,tf], the solution to the optimality system (71) is unique.Proof.
Suppose(S,A,T,Q,λ1,λ2,λ3,λ4) and (S¯,A¯,T¯,Q¯,λ¯1,λ¯2,λ¯3,λ¯4) are two different solutions of our optimality system (71). Let
(72)S=eλtm,A=eλtn,T=eλtp,Q=eλtq,λ1=e-λtr,λ2=e-λts,λ3=e-λtw,λ4=e-λtv,S¯=eλtm¯,A¯=eλtn¯,T¯=eλtp¯,Q¯=eλtq¯,λ¯1=e-λtr¯,λ¯2=e-λts¯,λ¯3=e-λtw¯,λ¯4=e-λtv¯,
where λ>0 is to be chosen.
Accordingly, we have(73)u1*(t)=min⁡{1,max⁡{0,βmn(s-r)eλtc1N}},u2*(t)=min⁡{1,max⁡{0,(s-w)nc2}},u¯1*(t)=min⁡{1,max⁡{0,βm¯n¯(s¯-r¯)eλtc1N}},u¯2*(t)=min⁡{1,max⁡{0,(s¯-w¯)n¯c2}}.
Now we substitute S=eλtm, A=eλtn, T=eλtp, Q=eλtq, λ1=e-λtr,  λ2=e-λts,  λ3=e-λtw,  λ4=e-λtv and S¯=eλtm¯,  A¯=eλtn¯,  T¯=eλtp¯,  Q¯=eλtq¯,  λ1¯=e-λtr¯,  λ2¯=e-λts¯,  λ3¯=e-λtw¯,  λ4¯=e-λtv¯  into the first ODE of (71), respectively; then we can obtain
(74)m˙+λm=μNe-λt-(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})×βmneλtN-μm,
for m and m¯, respectively. Similarly, we can derive
(75)n˙+λn=(1-min⁡{1,max⁡{1,βmn(s-r)eλtc1N}})βmneλtN+ξp-(min⁡{1,max⁡{0,n(s-w)c2}}+μ)n,
for n and n¯, respectively;
(76)p˙+λp=min⁡{1,max⁡{0,n(s-w)c2}}n-(μ+ξ+δ)p,
for p and p¯, respectively;
(77)q˙+λq=δp-μq,
for q and q¯, respectively;
(78)r˙-λr=r(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})βneλtN+μr-s(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})×βneλtN,
for r and r¯, respectively;
(79)s˙-λs=-eλt+r(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})βmeλtN-s(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})×βmeλtN+s(min⁡{1,max⁡{0,n(s-w)c2}}+μ)-w(min⁡{1,max⁡{0,n(s-w)c2}}),
for s and s¯, respectively;
(80)w˙-λw=-sξ+(μ+ξ+δ)w-δv,
for w and w¯, respectively;
(81)v˙-λv=μv,
for v and v¯, respectively.
By Lemma7, we can obtain
(82)|u1*(t)-u¯1*(t)|≤βeλtc1N|mn(s-r)-mn¯(s¯-r¯)|,|u2*(t)-u¯2*(t)|≤1c2|n(s-w)-n¯(s¯-w¯)|.
The equations form,n,p,q,r,s,w,v and the equations for m¯,n¯,p¯,q¯,r¯,s¯,w¯,v¯ are subtracted, respectively; then we multiply each equation by appropriate difference of functions and integrate from 0 to tf. Next, we add all eight integral equations and some inequality techniques to obtain uniqueness. The following calculation is similar; for the sake of simplicity, we only take m and m¯ for an example:(83)m˙-m¯˙+(μ+λ)(m-m¯)=βeλtN[-(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})mn+(1-min⁡{1,max⁡{0,βmn¯(s¯-r¯)eλtc1N}})mn¯].
Multiplying both sides of (83) by (m-m¯) and integrating from 0 to tf gives
(84)12(m-m¯)2(tf)+(μ+λ)∫0tf‍(m-m¯)2dt=∫0tf‍(m-m¯)βeλtNlllllllllll×[-(1-u1*)mn+(1-u¯1*)mn¯]dt=∫0tf‍(m-m¯)βeλtN[(u¯1*-1)(mn-mn¯+mn¯-m¯n¯)lllllllllllllllllllllllllllllllllll+mn(u1*-u¯1*)]dt≤βeλtfN∫0tf‍(m-m¯)llllllllllllllllllll×[u¯1*-1|(m|n-n¯|+|m-m¯|n¯)+mnβeλtc1N|llllllllllllllllllllllllll×mn(s-r)-mn¯(s¯-r¯)|(m|n-n¯|+|m-m¯|n¯)+mnβeλtc1N|]dt≤βeλtfN∫0tf‍(m-m¯)llllllllllllllllllll×[βeλtc1N|u¯1*-1|(m|n-n¯|+|m-m¯|n¯)+mnβeλtc1Nlllllllllllllllllllllllllll|mn(s-s¯)+(mn-mn¯)s¯llllllllllllllllllllllllllll-(mn(r-r¯)+(mn-mn¯)r¯)|βeλtc1N]dt≤βeλtfN∫0tf‍[|m-m¯|2|u¯1*-1||m||m-m¯||n-n¯|llllllllllllllllllllllll+|m-m¯|2|n¯|+|B||s¯||n¯||m-m¯|2llllllllllllllllllllllll+|B||mn||m-m¯||r-r¯|llllllllllllllllllllllll+|B||mn||m-m¯||s-s¯|llllllllllllllllllllllll+|B||s¯||m||m-m¯||n-n¯|llllllllllllllllllllllll+|B||r¯||m||m-m¯||n-n¯|llllllllllllllllllllllll+|B||r¯||n¯||m-m¯|2]dt≤βeλtfN∫0tf‍[|u¯*-1|2|m|((m-m¯)2+(n-n¯)2)llllllllllllllllllllllll+|n¯|(m-m¯)2llllllllllllllllllllllll+|B||mn|2((m-m¯)2+(s-s¯)2)llllllllllllllllllllllll+|B||s¯||m|2((m-m¯)2+(n-n¯)2)llllllllllllllllllllllll+|B||s¯||n¯|(m-m¯)2llllllllllllllllllllllll+|B||mn|2((m-m¯)2+(r-r¯)2)llllllllllllllllllllllll+|B||m||r¯|2((m-m¯)2+(n-n¯)2)llllllllllllllllllllllll+|B||r¯||n¯|(m-m¯)2|u¯*-1|2]dt.
In the above derivation, we use many scaling techniques for inequality or absolute inequality. Particularly, what should be noted is that to get the first inequality of above derivation, we use the estimation of|u1*(t)-u¯1*(t)| which has been given before; besides, for the sake of convenience, we note B=mn(βeλtf/N). Furthermore, we notice that the coefficients of all the eight terms in the last formula: (m-m¯)2+(n-n¯)2, (m-m¯)2, (m-m¯)2+(s-s¯)2, (m-m¯)2+(n-n¯)2, (m-m¯)2, (m-m¯)2+(r-r¯)2, (m-m¯)2+(n-n¯)2, (m-m¯)2, namely, (|u¯1*-1|/2)|m|, |n¯|, (|B||mn|)/2, (|B||s¯||m|)/2, |B||s¯|n¯, (|B||mn|)/2, (|B||m||r¯|)/2, |B||r¯||n¯| are nonnegative and bounded. So there exists a positive constant c2 such that
(85)12(m-m¯)2(tf)+(μ+λ)∫0tf‍(m-m¯)2dt≤c2βeλtfN∫0tf‍((m-m¯)2+(n-n¯)2lllllllllllllllllllllllllllll+(s-s¯)2+(r-r¯)2)dt.
Combining eight of these inequalities gives(86)12(m-m¯)(tf)+12(n-n¯)(tf)+12(p-p¯)(tf)+12(q-q¯)(tf)+12(r-r¯)(0)+12(s-s¯)(0)+12(w-w¯)(0)+12(v-v¯)(0)+(μ+λ)×∫0tf‍{(m-m¯)2+(n-n¯)2+(p-p¯)2lllllllllllllllllll+(q-q¯)2+(s-s¯)2llllllllllllllllll+(r-r¯)2+(w-w¯)2+(v-v¯)2)}dt≤B¯∫0tf{(m-m¯)2+(n-n¯)2+(p-p¯)2‍llllllllllllllllll+(q-q¯)2+(r-r¯)2+(s-s¯)2lllllllllllllllll+(w-w¯)2+(v-v¯)2}dt.
Thus, from the above inequality we can conclude that
(87)(μ+λ-B¯)∫0tf‍{(m-m¯)2+(n-n¯)2llllllllllllllllllllllllllll+(p-p¯)2+(q-q¯)2+(r-r¯)2llllllllllllllllllllllllllll+(s-s¯)2+(w-w¯)2+(v-v¯)2}dt≤0,
where B¯ depends on the coefficients and the bounds depend on m,n,p,q,r,s,w,v. If we choose λ such that μ+λ>B¯, then m=m¯, n=n¯, p=p¯, q=q¯, r=r¯, s=s¯, w=w¯, and v=v¯. Hence, the solution to the optimality system is unique. The proof is complete.


## 5.1. The Existence of Optimal Control

In order to investigate an effective campaign to control alcoholism in a community which pursue the goals of the minimized alcoholisms and more recovered individuals, we will reconsider the system (1) and use two control variables to reduce the numbers of alcoholics. The difference is that we will change the parameters u1,u2 into control variable u1(t),u2(t). Their aforementioned definitions allow us to do so. u1(t) is used to limit the proportion of the susceptible individual to contact with alcoholism, usually by propaganda and education, so that the susceptible individual can stay off alcoholism consciously and be free of “infection,” we can understand the effect of u1(t) is to prevent the the susceptible from contacting with the alcoholism. The control variable u2(t) is used to control the alcoholism to take appropriate treatment measures, such as taking pills or seeking other medical help. However, just as a coin has two sides, there will be a lot of costs generated during the control process. So it is advisable to balance between the costs and the alcohol effects. In view of this, our optimal control problem to minimize the objective functional is given by
(49)J(u1,u2)=∫0tf‍[A(t)+c12u12(t)+c22u22(t)]dt,
which subjects to system
(50)S′=μN-(1-u1(t))βSAN-μS,A′=(1-u1(t))βSAN+ξT-(u2(t)+μ)A,T′=u2(t)A-(μ+ξ+δ)T,Q′=δT-μQ,
with initial conditions
(51)S(0)=S0,A(0)=A0,T(0)=T0,Q(0)=Q0.
Here, ui(t)∈U≜{(u1,u2)∣ui(t)ismeasurableand0≤ui(t)≤1,forallt∈[0,tf]}, tf is the end time to be controlled, U is an admissible control set, ci, and i=1,2, are weight factors (positive constants) that adjust the intensity of two different control measures.Next, we will investigate the existence of the optimal control of the above-mentioned problem.Theorem 5.
There exists an optimal control pairu*=(u1*,u2*)∈U such that
(52)J(u1*,u2*)=min⁡J(u1,u2),u1(t),u2(t)∈U
subjects to the control system (1) with initial conditions (50).Proof.
To prove the existence of an optimal control, according to the classic literature [29], we have to show the following.(1)
The control and state variables are nonnegative values.(2)
The control setU is convex and closed.(3)
The right side of the state system is bounded by linear function in the state and control variables.(4)
The integrand of the objective functional is concave onU.(5)
There exist constantsd1,d2>0 and α>1 such that the integrand L(t;u1;u2)≜A(t)+(c1/2)u12(t)+(c2/2)u22(t) of the objective functional satisfies
(53)L(t;u1;u2)≥d1(|u1|2+|u2|2)α/2-d2statements (1), (2) and (3) are obvious satisfied, we only need to test and verify the latter two ones. Since the four state variables have been all proved to be up bounded by N, we will get the following equalities:
(54)S′≤μN,A′≤(1-u1(t))βS+ξT,T′≤u2(t)A,Q′≤δT,
so the fourth condition is set up. As for the last condition,
(55)L(t;u1;u2)≥d1(|u1|2+|u2|2)α/2-d2
is also true, when we choose d1=min⁡{c1/2,c2/2}, and for all d2∈R+,α=2. The proof is complete.We next come to the core of this section.


## 5.2. The Characterization of the Optimal Control

With the existence of the optimal control pairs established, we now present the optimality system and use a result from [30]; we can easily know the existence of the solutions to the optimality system (71) which will be gotten later. Firstly, we come to discuss the theorem that relates to the characterization of the optimal control. The optimality system can be used to compute candidates for optimal control pairs. To do this, we begin by defining an augmented Hamiltonian H with penalty terms for the control constraints as follows:
(56)H=A(t)+c12u12(t)+c22u22(t)+λ1[μN-(1-u1(t))βSAN-μS]+λ2[(1-u1(t))βSAN+ξT-(u2(t)+μ)A]+λ3[u2(t)A-(μ+ξ+δ)T]+λ4(δT-μQ)-w11u1(t)-w12(1-u1(t))-w21u2(t)-w22(1-u2(t)),
where wij(t)≥0 are the penalty multipliers satisfying
(57)w11(t)u1(t)=w12(t)(1-u1(t))=0atoptimalcontrolu1*,w21(t)u2(t)=w22(t)(1-u2(t))=0atoptimalcontrolu2*.Theorem 6.
Given optimal control pairs(u1*,u2*) and solutions S(t),A(t),T(t),Q(t) of the corresponding state system (50), there exist adjoint variables λi, i=1,2,3,4, satisfying
(58)λ1′=λ1(1-u1(t))βAN+μλ1-λ2(1-u1(t))βAN,λ2′=-1+λ1(1-u1(t))βSN-λ2(1-u1(t))βSN+λ2(μ+u2(t))-λ3u2(t),λ3′=-λ2ξ+λ3(μ+ξ+δ)-λ4δ,λ4′=μλ4,
with the terminal conditions
(59)λi(tf)=0,i=1,2,3,4.
Furthermore, (u1*,u2*) are represented by
(60)u1*=min⁡(1,max⁡(0,βSA(λ2-λ1)c1N)),u2*=min⁡(1,max⁡(0,A(λ2-λ3)c2)).Proof.
According to Pontryagin Maximum Principle [29–31], we first differentiate the Hamiltonian operator H, with respect to states. Then the adjoint system can be written as
(61)λ1′=-∂H∂S,λ2′=-∂H∂A,λ3′=-∂H∂T,λ4′=-∂H∂Q.
The terminal condition (56) of adjoint equations is given by λi(tf)=0, i=1,2,3,4.
To obtain the necessary conditions of optimality (59), we also differentiate the Hamiltonian operator H, with respect to U=(u1,u2) and set them equal to zero; then
(62)∂H∂u1=c1u1(t)+λ1βSAN-λ2βSAN-w11+w12=0,∂H∂u2=c2u2(t)-λ2A+λ3A-w21+w22=0.
By solving the optimal control, we obtain(63)u1*=1c1[(λ2-λ1)βSAN+w11-w12].
To determine an explicit expression for the optimal control withoutw11 and w12, a standard optimality technique is utilized [29]. We consider the following three cases.(i)
On the set{t∣0<u1*(t)<1}, we have w11(t)=w12(t)=0. Hence, the optimal control is
(64)u1*=1c1[(λ2-λ1)βSAN].(ii)
On the set{t∣u1*(t)=1}, we have w11(t)=0. Hence,
(65)1=u1*(t)=1c1[(λ2-λ1)βSAN-w12].This implies that(66)1c1(λ2-λ1)βSAN≥1sincew12(t)≥0.(iii)
On the set{t∣u1*(t)=0}, we have w12(t)=0. Hence,
(67)0=u1*(t)=1c1[(λ2-λ1)βSAN+w11].This implies that(68)1c1(λ2-λ1)βSAN≤0sincew11(t)≥0.
Combining these results, the optimal controlu1*(t) is characterized as
(69)u1*=min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}.
Using the similar arguments, we can also obtain the other optimal control function
(70)u2*=min⁡{1,max⁡{0,(λ2-λ3)Ac2}}.
The proof is complete.We point out that the optimality system consists of the state system (50) with the initial conditions S(0),A(0),T(0),Q(0), the adjoint (or costate) system (58) with the terminal conditions (59), and the optimality condition (60). Any optimal control pairs must satisfy this optimality system. For the convenience of subsequent numerical simulation in Section 6, we give the optimality system as follows:
(71)S′=μN-(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βSAN-μS,A′=(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))βSAN+ξT-(min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)+μ)A,T′=min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)A-(μ+ξ+δ)T,Q′=δT-μQ,λ1′=λ1(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βAN+μλ1-λ2(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βAN,λ2′=-1+λ1(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}}(t))×βSN-λ3·min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)-λ2(1-min⁡{1,max⁡{0,1c1[(λ2-λ1)βSAN]}})×βSN+λ2(μ+min⁡{1,max⁡{0,(λ2-λ3)Ac2}}(t)),λ3′=-λ2ξ+λ3(μ+ξ+δ)-λ4δ,λ4′=μλ4,S(0)=S0,A(0)=A0,T(0)=T0,Q(0)=Q0,λi(tf)=0,i=1,2,3,4.


## 5.3. The Uniqueness of Optimal Control

Due to the a priori boundedness of the state, adjoint functions, and the resulting Lipschitz structure of the ODEs, we can obtain the uniqueness of the optimal control.Lemma 7 (see [23]).
The functionu*(s)=min⁡(b,max⁡(s,a)) is Lipschitz continuous in s, where a<b are some fixed positive constants.Theorem 8.
For allt∈[0,tf], the solution to the optimality system (71) is unique.Proof.
Suppose(S,A,T,Q,λ1,λ2,λ3,λ4) and (S¯,A¯,T¯,Q¯,λ¯1,λ¯2,λ¯3,λ¯4) are two different solutions of our optimality system (71). Let
(72)S=eλtm,A=eλtn,T=eλtp,Q=eλtq,λ1=e-λtr,λ2=e-λts,λ3=e-λtw,λ4=e-λtv,S¯=eλtm¯,A¯=eλtn¯,T¯=eλtp¯,Q¯=eλtq¯,λ¯1=e-λtr¯,λ¯2=e-λts¯,λ¯3=e-λtw¯,λ¯4=e-λtv¯,
where λ>0 is to be chosen.
Accordingly, we have(73)u1*(t)=min⁡{1,max⁡{0,βmn(s-r)eλtc1N}},u2*(t)=min⁡{1,max⁡{0,(s-w)nc2}},u¯1*(t)=min⁡{1,max⁡{0,βm¯n¯(s¯-r¯)eλtc1N}},u¯2*(t)=min⁡{1,max⁡{0,(s¯-w¯)n¯c2}}.
Now we substitute S=eλtm, A=eλtn, T=eλtp, Q=eλtq, λ1=e-λtr,  λ2=e-λts,  λ3=e-λtw,  λ4=e-λtv and S¯=eλtm¯,  A¯=eλtn¯,  T¯=eλtp¯,  Q¯=eλtq¯,  λ1¯=e-λtr¯,  λ2¯=e-λts¯,  λ3¯=e-λtw¯,  λ4¯=e-λtv¯  into the first ODE of (71), respectively; then we can obtain
(74)m˙+λm=μNe-λt-(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})×βmneλtN-μm,
for m and m¯, respectively. Similarly, we can derive
(75)n˙+λn=(1-min⁡{1,max⁡{1,βmn(s-r)eλtc1N}})βmneλtN+ξp-(min⁡{1,max⁡{0,n(s-w)c2}}+μ)n,
for n and n¯, respectively;
(76)p˙+λp=min⁡{1,max⁡{0,n(s-w)c2}}n-(μ+ξ+δ)p,
for p and p¯, respectively;
(77)q˙+λq=δp-μq,
for q and q¯, respectively;
(78)r˙-λr=r(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})βneλtN+μr-s(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})×βneλtN,
for r and r¯, respectively;
(79)s˙-λs=-eλt+r(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})βmeλtN-s(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})×βmeλtN+s(min⁡{1,max⁡{0,n(s-w)c2}}+μ)-w(min⁡{1,max⁡{0,n(s-w)c2}}),
for s and s¯, respectively;
(80)w˙-λw=-sξ+(μ+ξ+δ)w-δv,
for w and w¯, respectively;
(81)v˙-λv=μv,
for v and v¯, respectively.
By Lemma7, we can obtain
(82)|u1*(t)-u¯1*(t)|≤βeλtc1N|mn(s-r)-mn¯(s¯-r¯)|,|u2*(t)-u¯2*(t)|≤1c2|n(s-w)-n¯(s¯-w¯)|.
The equations form,n,p,q,r,s,w,v and the equations for m¯,n¯,p¯,q¯,r¯,s¯,w¯,v¯ are subtracted, respectively; then we multiply each equation by appropriate difference of functions and integrate from 0 to tf. Next, we add all eight integral equations and some inequality techniques to obtain uniqueness. The following calculation is similar; for the sake of simplicity, we only take m and m¯ for an example:(83)m˙-m¯˙+(μ+λ)(m-m¯)=βeλtN[-(1-min⁡{1,max⁡{0,βmn(s-r)eλtc1N}})mn+(1-min⁡{1,max⁡{0,βmn¯(s¯-r¯)eλtc1N}})mn¯].
Multiplying both sides of (83) by (m-m¯) and integrating from 0 to tf gives
(84)12(m-m¯)2(tf)+(μ+λ)∫0tf‍(m-m¯)2dt=∫0tf‍(m-m¯)βeλtNlllllllllll×[-(1-u1*)mn+(1-u¯1*)mn¯]dt=∫0tf‍(m-m¯)βeλtN[(u¯1*-1)(mn-mn¯+mn¯-m¯n¯)lllllllllllllllllllllllllllllllllll+mn(u1*-u¯1*)]dt≤βeλtfN∫0tf‍(m-m¯)llllllllllllllllllll×[u¯1*-1|(m|n-n¯|+|m-m¯|n¯)+mnβeλtc1N|llllllllllllllllllllllllll×mn(s-r)-mn¯(s¯-r¯)|(m|n-n¯|+|m-m¯|n¯)+mnβeλtc1N|]dt≤βeλtfN∫0tf‍(m-m¯)llllllllllllllllllll×[βeλtc1N|u¯1*-1|(m|n-n¯|+|m-m¯|n¯)+mnβeλtc1Nlllllllllllllllllllllllllll|mn(s-s¯)+(mn-mn¯)s¯llllllllllllllllllllllllllll-(mn(r-r¯)+(mn-mn¯)r¯)|βeλtc1N]dt≤βeλtfN∫0tf‍[|m-m¯|2|u¯1*-1||m||m-m¯||n-n¯|llllllllllllllllllllllll+|m-m¯|2|n¯|+|B||s¯||n¯||m-m¯|2llllllllllllllllllllllll+|B||mn||m-m¯||r-r¯|llllllllllllllllllllllll+|B||mn||m-m¯||s-s¯|llllllllllllllllllllllll+|B||s¯||m||m-m¯||n-n¯|llllllllllllllllllllllll+|B||r¯||m||m-m¯||n-n¯|llllllllllllllllllllllll+|B||r¯||n¯||m-m¯|2]dt≤βeλtfN∫0tf‍[|u¯*-1|2|m|((m-m¯)2+(n-n¯)2)llllllllllllllllllllllll+|n¯|(m-m¯)2llllllllllllllllllllllll+|B||mn|2((m-m¯)2+(s-s¯)2)llllllllllllllllllllllll+|B||s¯||m|2((m-m¯)2+(n-n¯)2)llllllllllllllllllllllll+|B||s¯||n¯|(m-m¯)2llllllllllllllllllllllll+|B||mn|2((m-m¯)2+(r-r¯)2)llllllllllllllllllllllll+|B||m||r¯|2((m-m¯)2+(n-n¯)2)llllllllllllllllllllllll+|B||r¯||n¯|(m-m¯)2|u¯*-1|2]dt.
In the above derivation, we use many scaling techniques for inequality or absolute inequality. Particularly, what should be noted is that to get the first inequality of above derivation, we use the estimation of|u1*(t)-u¯1*(t)| which has been given before; besides, for the sake of convenience, we note B=mn(βeλtf/N). Furthermore, we notice that the coefficients of all the eight terms in the last formula: (m-m¯)2+(n-n¯)2, (m-m¯)2, (m-m¯)2+(s-s¯)2, (m-m¯)2+(n-n¯)2, (m-m¯)2, (m-m¯)2+(r-r¯)2, (m-m¯)2+(n-n¯)2, (m-m¯)2, namely, (|u¯1*-1|/2)|m|, |n¯|, (|B||mn|)/2, (|B||s¯||m|)/2, |B||s¯|n¯, (|B||mn|)/2, (|B||m||r¯|)/2, |B||r¯||n¯| are nonnegative and bounded. So there exists a positive constant c2 such that
(85)12(m-m¯)2(tf)+(μ+λ)∫0tf‍(m-m¯)2dt≤c2βeλtfN∫0tf‍((m-m¯)2+(n-n¯)2lllllllllllllllllllllllllllll+(s-s¯)2+(r-r¯)2)dt.
Combining eight of these inequalities gives(86)12(m-m¯)(tf)+12(n-n¯)(tf)+12(p-p¯)(tf)+12(q-q¯)(tf)+12(r-r¯)(0)+12(s-s¯)(0)+12(w-w¯)(0)+12(v-v¯)(0)+(μ+λ)×∫0tf‍{(m-m¯)2+(n-n¯)2+(p-p¯)2lllllllllllllllllll+(q-q¯)2+(s-s¯)2llllllllllllllllll+(r-r¯)2+(w-w¯)2+(v-v¯)2)}dt≤B¯∫0tf{(m-m¯)2+(n-n¯)2+(p-p¯)2‍llllllllllllllllll+(q-q¯)2+(r-r¯)2+(s-s¯)2lllllllllllllllll+(w-w¯)2+(v-v¯)2}dt.
Thus, from the above inequality we can conclude that
(87)(μ+λ-B¯)∫0tf‍{(m-m¯)2+(n-n¯)2llllllllllllllllllllllllllll+(p-p¯)2+(q-q¯)2+(r-r¯)2llllllllllllllllllllllllllll+(s-s¯)2+(w-w¯)2+(v-v¯)2}dt≤0,
where B¯ depends on the coefficients and the bounds depend on m,n,p,q,r,s,w,v. If we choose λ such that μ+λ>B¯, then m=m¯, n=n¯, p=p¯, q=q¯, r=r¯, s=s¯, w=w¯, and v=v¯. Hence, the solution to the optimality system is unique. The proof is complete.


## 6. Numerical Simulation

### 6.1. The Simulation of State System (1) without Control Parameters

For the sake of simplicity but without loss of generality, we will perform the numerical simulation of state system (1) with parameters u1=0, u2=0. Before illustrating the analytic properties of the alcoholism model (1), we will target the populations in the environment of a community or a university, for example, the school of material science and engineering in our university, that is, Lan zhou University of Technology (LUT for short), owing to the accurate and available information we can obtain. Referring to the information provided by the admissions office of LUT, this school will enroll almost 1200 undergraduates and almost 300 various postgraduates at the beginning of fall semester; at the same time, there will be almost 1500 various students graduated and left this school, so the scale of students in school remained almost 6000; we can take the total population N=6000. In this simulation, we will take September as the initial time and units in one week, period in one year. According to the investigations of the student union implemented in September every year, we can take initial values as S(0)=4500, A(0)=1000, T(0)=300, and Q(0)=200. It seems that the alcoholism is a little bit more, but it is rather natural because many freshmen feel confused when they are faced with the new environment and a new lifestyle; many of them have no better choice but gather together to drink in small groups to mediate the anxiety and get to know each other; over time, some of them develop the habit of drinking. To a certain extent, for example, the frequent drinking badly affects their study; we can classify them into the alcoholism compartment. Other initial values seem more reasonable, so we need no more explanation. As we know, alcoholism death is seldom happen within one year, so we omit mortality from alcoholism; then how to understand the recruitment rate as well as natural death rate μ? We can treat the freshmen admission as the recruitment population and graduation students as the natural “death” parts. So we can take μ=1500/6000=0.25, which is exactly consistent with the value in [16]. As for the infection rate β and recovery rate δ, we will let them be variables, since the drinking behaviors are related to many factors such as the season and the pressure.According to the data we get from the student union, we chooseξ=0.4. To summarize, we list the values of the parameters in Table 1. Using the values of parameters in Table 1, we can plot Figures 2 and 3 which are on condition R0<1 and R0≥1, respectively. From Figure 2, we easily know when R0<1 holds; the solution of system (1) tends to the alcohol free equilibrium E0 and verifies the global stability of E0. While seen from Figure 3, we also know that if R0≥1 holds, the solution of system (1) tends to the alcoholism equilibrium E* and verifies the global stability of E*.Table 1
The parameters description of model (1).
Parameter
Description
Values
μ
Natural birth rate or death rate
0.25
β
Transmission coefficient between alcoholism and susceptibles
Variables
N
The total populations to be considered
6000
δ
The rate of populations quitting from alcoholism permanently after treatment
Variables
ξ
The rate of populations failed in treatment and returned to be alcoholic
0.4Figure 2
WhenR0<1, the alcohol free equilibrium E0 is globally asymptotically stable (β=0.2,δ=0.3, and R0=0.71698).Figure 3
WhenR0≥1, the alcoholism equilibrium E*=(4466,476,32,26) corresponding to the given parameters is globally asymptotically stable (β=0.3,δ=0.2, and R0=1.08511).


### 6.2. The Sensitivity ofR0 about Two Control Parameters

Although from the expression of the model reproduction numberR0, we can easily find out the fact that the two control variables, that is, u1 and u2, attribute to reducing the severity of alcoholism; we will still depict the graph between R0 and the two control variables to see more intuitiveness see Figure 4. It seems from the figure that R0 is a monotonically decreasing function about two control parameters, so it is advisable to take two approaches simultaneously to control the alcoholics.Figure 4
The relationship betweenR0 and two control variables u1 and u2.


### 6.3. The Simulation of Optimality System (71)

In this subsection, we will investigate numerically the optimal solution to optimality system (71) by numerical method from [32]; the optimality system is solved with a fourth-order Runge-Kutta scheme. Beginning with a guess for the control variables, the state system is solved forward in time and then those values of state are used to solve the adjoint equations backward in time. The controls are updated at the end of each iteration using the values of optimal controls obtained lastly. The iterations continue until convergence takes place.In the simulations, we choose the available variable values as Table1 shows; besides, β=0.3,δ=0.2. The initial value of model (1) is assumed to be S(0)=4500, A(0)=1000, T(0)=300, and Q(0)=200 as before.The ideal weights in objective functional are very difficult to obtain in reality; it needs much work on data mining and fitting. Hence, the acquisition of appropriate practical weights is still a difficult problem and remains for further investigations. The cost associated withA(t) and u1(t) mainly includes the cost of dangerous behaviours during the alcoholism time and educating the public, while the cost associated with u2(t) mainly comes from health professional and the medical resource including medicines and nursing care. In view of this and taking the expressions of u1,u2 into account, after many numerical simulations, we finally give weighting coefficients as c1=102;c2=104. It should be pointed out that the weights here are of only theoretical interest to reveal the control strategies proposed in this paper. Another point to note is that the maximum control is very difficult to achieve in reality, so we will omit the situation of the maximum control during the series of simulations.Next, we will make some necessary instructions and explanations to the above simulation graphs. Figures5, 6, 7, and 8 depict the number of four compartments under different control levels when we choose the weight coefficients in objective function to be c1=102;c2=104. From the four simulation graphs, we can observe the following simple facts, in reducing the total number of alcoholisms and increasing the number of susceptibles; the effectiveness of various control measures is as follows: optimal control is evidently better than middle control, and middle control is better than single control u2, single control u2 is better than single control u1, while single control u1 is much better than no control.Figure 5
Number of the susceptibles when we choose the weights in objective function arec1=102;c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 6
Number of the alcoholics when we choose the weights in objective function arec1=102; c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 7
Number of the persons in treatment when we choose the weights in objective function arec1=102; c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 8
Number of people in quitting compartment when we choose the weights in objective function arec1=102; c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figures9 and 10 depict the optimal control law of u1,u2, respectively. In the beginning of the simulation, the control effort of u1 should be decreased from 0.65 to 0.5 within the first month, and over the next week, it should be increased to the maximum control until 50 weeks, then rapidly decreased to 0 at the end of the simulation. As for the control u2, it should start from around 0.5 due to the initial alcoholics then increase to 0.55 within one week since the rapid infection and next decrease to almost 0 since the effectiveness of treatment in three weeks, but with the infection going on, the control effort of u2 should gradually increase to the maximum and maintain this level until the tenth week for the purpose of consolidation therapy and preventing rebound; hereafter, it should be gradually decreased to the level of almost 0.43 until the fifty weeks then quickly decreased to 0 in the end.Figure 9
Figures of the optimal controlu1 when we choose the weight in objective function are c1=102; c2=104.Figure 10
Figures of the optimal controlu2 when we choose the weight in objective function are c1=102; c2=104.In order to investigate the influence of different weight coefficients in the objective functional on the effect of controlling, at the same time, for a better comparison, we will change the weight coefficients in objective function intoc2=104; c1=102, and we will list the corresponding numerical simulation results as Figures 11–16 show.Figure 11
Number of the susceptibles when we choose the weights in objective function arec1=104, c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 12
Number of the alcoholics when we choose the weights in objective function arec1=104, c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 13
Number of the persons in treatment when we choose the weights in objective function arec1=104, c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 14
Number of people in quitting compartment when we choose the weights in objective function arec1=104,  c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 15
Figures of the optimal controlu1 when we choose the weight in objective function are c1=104, c2=102.Figure 16
Figures of the optimal controlu2 when we choose the weight in objective function are c1=104, c2=102.When we change the weight coefficients in objective function intoc1=104, c2=102, we find that the results of simulations derived from the graphs are very similar to the ones before. We speculate that the most likely reasons of this result are due to three respects; one is that the weight coefficients are not too sensitive in the numerical simulation, and another possible reason is that both of the two controls are important, in some sense, equivalently important. The last but not the most unlikely reason is that we have not found the most appropriate weight coefficients in the simulation, which is very difficult to find as previously mentioned.


## 6.1. The Simulation of State System (1) without Control Parameters

For the sake of simplicity but without loss of generality, we will perform the numerical simulation of state system (1) with parameters u1=0, u2=0. Before illustrating the analytic properties of the alcoholism model (1), we will target the populations in the environment of a community or a university, for example, the school of material science and engineering in our university, that is, Lan zhou University of Technology (LUT for short), owing to the accurate and available information we can obtain. Referring to the information provided by the admissions office of LUT, this school will enroll almost 1200 undergraduates and almost 300 various postgraduates at the beginning of fall semester; at the same time, there will be almost 1500 various students graduated and left this school, so the scale of students in school remained almost 6000; we can take the total population N=6000. In this simulation, we will take September as the initial time and units in one week, period in one year. According to the investigations of the student union implemented in September every year, we can take initial values as S(0)=4500, A(0)=1000, T(0)=300, and Q(0)=200. It seems that the alcoholism is a little bit more, but it is rather natural because many freshmen feel confused when they are faced with the new environment and a new lifestyle; many of them have no better choice but gather together to drink in small groups to mediate the anxiety and get to know each other; over time, some of them develop the habit of drinking. To a certain extent, for example, the frequent drinking badly affects their study; we can classify them into the alcoholism compartment. Other initial values seem more reasonable, so we need no more explanation. As we know, alcoholism death is seldom happen within one year, so we omit mortality from alcoholism; then how to understand the recruitment rate as well as natural death rate μ? We can treat the freshmen admission as the recruitment population and graduation students as the natural “death” parts. So we can take μ=1500/6000=0.25, which is exactly consistent with the value in [16]. As for the infection rate β and recovery rate δ, we will let them be variables, since the drinking behaviors are related to many factors such as the season and the pressure.According to the data we get from the student union, we chooseξ=0.4. To summarize, we list the values of the parameters in Table 1. Using the values of parameters in Table 1, we can plot Figures 2 and 3 which are on condition R0<1 and R0≥1, respectively. From Figure 2, we easily know when R0<1 holds; the solution of system (1) tends to the alcohol free equilibrium E0 and verifies the global stability of E0. While seen from Figure 3, we also know that if R0≥1 holds, the solution of system (1) tends to the alcoholism equilibrium E* and verifies the global stability of E*.Table 1
The parameters description of model (1).
Parameter
Description
Values
μ
Natural birth rate or death rate
0.25
β
Transmission coefficient between alcoholism and susceptibles
Variables
N
The total populations to be considered
6000
δ
The rate of populations quitting from alcoholism permanently after treatment
Variables
ξ
The rate of populations failed in treatment and returned to be alcoholic
0.4Figure 2
WhenR0<1, the alcohol free equilibrium E0 is globally asymptotically stable (β=0.2,δ=0.3, and R0=0.71698).Figure 3
WhenR0≥1, the alcoholism equilibrium E*=(4466,476,32,26) corresponding to the given parameters is globally asymptotically stable (β=0.3,δ=0.2, and R0=1.08511).


## 6.2. The Sensitivity ofR0 about Two Control Parameters

Although from the expression of the model reproduction numberR0, we can easily find out the fact that the two control variables, that is, u1 and u2, attribute to reducing the severity of alcoholism; we will still depict the graph between R0 and the two control variables to see more intuitiveness see Figure 4. It seems from the figure that R0 is a monotonically decreasing function about two control parameters, so it is advisable to take two approaches simultaneously to control the alcoholics.Figure 4
The relationship betweenR0 and two control variables u1 and u2.


## 6.3. The Simulation of Optimality System (71)

In this subsection, we will investigate numerically the optimal solution to optimality system (71) by numerical method from [32]; the optimality system is solved with a fourth-order Runge-Kutta scheme. Beginning with a guess for the control variables, the state system is solved forward in time and then those values of state are used to solve the adjoint equations backward in time. The controls are updated at the end of each iteration using the values of optimal controls obtained lastly. The iterations continue until convergence takes place.In the simulations, we choose the available variable values as Table1 shows; besides, β=0.3,δ=0.2. The initial value of model (1) is assumed to be S(0)=4500, A(0)=1000, T(0)=300, and Q(0)=200 as before.The ideal weights in objective functional are very difficult to obtain in reality; it needs much work on data mining and fitting. Hence, the acquisition of appropriate practical weights is still a difficult problem and remains for further investigations. The cost associated withA(t) and u1(t) mainly includes the cost of dangerous behaviours during the alcoholism time and educating the public, while the cost associated with u2(t) mainly comes from health professional and the medical resource including medicines and nursing care. In view of this and taking the expressions of u1,u2 into account, after many numerical simulations, we finally give weighting coefficients as c1=102;c2=104. It should be pointed out that the weights here are of only theoretical interest to reveal the control strategies proposed in this paper. Another point to note is that the maximum control is very difficult to achieve in reality, so we will omit the situation of the maximum control during the series of simulations.Next, we will make some necessary instructions and explanations to the above simulation graphs. Figures5, 6, 7, and 8 depict the number of four compartments under different control levels when we choose the weight coefficients in objective function to be c1=102;c2=104. From the four simulation graphs, we can observe the following simple facts, in reducing the total number of alcoholisms and increasing the number of susceptibles; the effectiveness of various control measures is as follows: optimal control is evidently better than middle control, and middle control is better than single control u2, single control u2 is better than single control u1, while single control u1 is much better than no control.Figure 5
Number of the susceptibles when we choose the weights in objective function arec1=102;c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 6
Number of the alcoholics when we choose the weights in objective function arec1=102; c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 7
Number of the persons in treatment when we choose the weights in objective function arec1=102; c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 8
Number of people in quitting compartment when we choose the weights in objective function arec1=102; c2=104 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figures9 and 10 depict the optimal control law of u1,u2, respectively. In the beginning of the simulation, the control effort of u1 should be decreased from 0.65 to 0.5 within the first month, and over the next week, it should be increased to the maximum control until 50 weeks, then rapidly decreased to 0 at the end of the simulation. As for the control u2, it should start from around 0.5 due to the initial alcoholics then increase to 0.55 within one week since the rapid infection and next decrease to almost 0 since the effectiveness of treatment in three weeks, but with the infection going on, the control effort of u2 should gradually increase to the maximum and maintain this level until the tenth week for the purpose of consolidation therapy and preventing rebound; hereafter, it should be gradually decreased to the level of almost 0.43 until the fifty weeks then quickly decreased to 0 in the end.Figure 9
Figures of the optimal controlu1 when we choose the weight in objective function are c1=102; c2=104.Figure 10
Figures of the optimal controlu2 when we choose the weight in objective function are c1=102; c2=104.In order to investigate the influence of different weight coefficients in the objective functional on the effect of controlling, at the same time, for a better comparison, we will change the weight coefficients in objective function intoc2=104; c1=102, and we will list the corresponding numerical simulation results as Figures 11–16 show.Figure 11
Number of the susceptibles when we choose the weights in objective function arec1=104, c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 12
Number of the alcoholics when we choose the weights in objective function arec1=104, c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 13
Number of the persons in treatment when we choose the weights in objective function arec1=104, c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 14
Number of people in quitting compartment when we choose the weights in objective function arec1=104,  c2=102 and under different optimal control strategies, that is, (1) without control: u1=u2=0; (2) single control: u1=0.5, u2=0; (3) single control: u1=0, u2=0.5; (4) middle control: u1=u2=0.5; (5) optimal control: u1=u1*, u2=u2*.Figure 15
Figures of the optimal controlu1 when we choose the weight in objective function are c1=104, c2=102.Figure 16
Figures of the optimal controlu2 when we choose the weight in objective function are c1=104, c2=102.When we change the weight coefficients in objective function intoc1=104, c2=102, we find that the results of simulations derived from the graphs are very similar to the ones before. We speculate that the most likely reasons of this result are due to three respects; one is that the weight coefficients are not too sensitive in the numerical simulation, and another possible reason is that both of the two controls are important, in some sense, equivalently important. The last but not the most unlikely reason is that we have not found the most appropriate weight coefficients in the simulation, which is very difficult to find as previously mentioned.


## 7. Conclusions

In this paper, we formulate an alcoholics quitting model and firstly investigate the variation discipline of various populations from the perspective of global stability; then we propose an objective functional to examine two different control measures (i.e., prevention and treatment) on the effect of alcohol. The basic reproduction number of the model was derived and the global stability of the two equilibria is given. From the expression of the basic reproduction numberR0 and related numerical simulation, we can easily see that the two control strategies are effective in the alcoholics process.Using Pontryagin’s Maximum Principle, we firstly determine the necessary conditions for existence of optimal control pairs. The uniqueness of the solution to the optimality system (71) is derived by the classical method of contradiction. Numerical simulations of the model suggest that the two different groups of weights in the objective function have much similar effects on the transmission of the alcoholism; from this point, the two control measures are almost equally important in controlling the alcoholism, although they will probably have great influences on the cost of the objective function. From the simulation figures, it seems that the effect of optimal control, which is measured by the reduction in the number of alcoholics and the increase in the number of susceptibles, is much better than other control strategies as noted earlier in the simulation section. According to the real-time curve of two optimal controls, we point out the specific implementation methods of optimal control which can be achieved in practice.



---

*Source: 954069-2014-03-02.xml*