# The Optimal Control and MLE of Parameters of a Stochastic Single-Species System

Authors: Huili Xiang; Zhijun Liu
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

This paper investigates the optimal control and MLE (maximum likelihood estimation) for a single-species system subject to random perturbation. With the help of the techniques of stochastic analysis and mathematical statistics, sufficient conditions for the optimal control threshold value, the optimal control moment, and the maximum likelihood estimation of parameters are established, respectively. An example is presented to illustrate the feasibility of our theoretical results.


---

## Body

## 1. Introduction

The Malthus model is usually expressed as(1.1)dx(t)dt=rx(t),
where x(0)=x0>0, x(t) stands for the density of species x at t moment, and r is the intrinsic growth rate. As everyone knows, model (1.1) has epoch-making significance in mathematics and ecology and later, many deterministic mathematical models have been widely studied (see [1–5]). In fact, a population system is inevitably affected by the environmental noise in the real world. As a consequence, it is reasonable to study a corresponding stochastic model. Notice that some recent results, especially on optimal control, for the following stochastic model
(1.2)dx(t)x(t)=rdt+σdw(t),
have been obtained (see [6–9]), where  w(t)  stands for the standard Brownian motion. However, for some pest populations, their generations are nonoverlapping (e.g., poplar and willow weevil, osier weevil and paranthrene tabaniformis) and the discrete models are more appropriate than the continuous ones. Compared with the continuous ones, the study on discrete mathematical models is more challenging. Inspired by [1–12], in this paper we will consider the following discrete model of system (1.2)
(1.3)x(n)=x(n-1)exp(r+σεn),
where x(0)>0,  εi~N(0,1),i=1,2,…,n,  and any two of them are independent.  σ stands for the noise intensity. We will focus on the optimal control threshold value, the optimal control moment, and the maximum likelihood estimation of parameters. To the best of our knowledge, no work has been done for system (1.3).The rest of this paper is organized as follows. In Section2, some preliminaries are introduced. In Section 3, we give three results of this paper. As applications of our main results; an example is presented to illustrate the feasibility of our theoretical results in Section 4.


## 2. Preliminaries

In this section, we summarize several definitions, assumptions, and lemmas which are useful for the later sections.Definition 2.1.
Only when the quantity of pest population reachesU one starts to control the pest population, and the real number U is called to be a control threshold value.Definition 2.2.
Until theN0th generation, the total quantity of pest population first reaches the control threshold value, then one says that  N0 is the first reaching time.Two main goals of this paper are to seek the optimal control threshold value and the optimal control moment from the point of view of the lowest control cost. Considering that the practical control to some pest population must be in the limited time range, we give the first assumption:(H1)
n≤n0, where n is the number of generation of pest population in a control period and n0 is a positive integer.DenoteM=max{x(0),x(1),…,x(n0)}. Usually, at the beginning, the number of pest population is very small, so we give the second assumption:(H2)
The first reaching timeN0>0.Let the life period of pest populationx be τ, we should annihilate pest at N0τ moment from the point of view of the lowest control cost. We further give the third assumption.(H3)
The number of pest populationx will not reach the extent which can cause damage again after being annihilated.By(H3), we have
(2.1)P{M<U}=P{N0>n0}orP{M≥U}=P{N0≤n0}.
So we can give the expression of the total loss caused by pest and expending for annihilating pest, respectively. It is obvious that the loss caused by pest population comes from the quantity of population and damaging time. We need to the fourth assumption(H4)  The generation of pest population is nonoverlapping.On one hand, the loss caused by pest can be expressed as(2.2)S(N0)=k1∑n=0N0E[x(n)]τ,
where k1 stands for the loss caused by unit number pest in one generation, E[x(n)] is the mean function of x(n). On the other hand, the expending for annihilating pest can be expressed as
(2.3)k2H(M-U),
where H(x) is defined by
(2.4)H(x)={1,x>0,0,x≤0,
that is,
(2.5)H(M-U)={1,M>U,0,M≤U,
where  k2 stands for the expending for annihilating  pest once. Since  S(N0) is dependent on random variable N0 and k2H(M-U) is dependent on random variable M and threshold value U, the total cost is a random variable, which can be expressed as
(2.6)J(U)=E[S(N0)]+E[k2H(M-U)]=E[S(N0)]+k2P{M>U}.
Thus, we need to search for U* such that J(U*) is minimum and consequently, we can give the optimal control moment.Next, we will give some lemmas which are very important to the proofs of three theorems in the following section.Lemma 2.3.
The solution of system (1.3) can be expressed as
(2.7)x(n)=x(0)exp(nr+σ∑i=1nεi).Proof.
By (1.3), we have
(2.8)x(n)x(n-1)=exp(r+σεn),(n=1,2…).
Thus, one has
(2.9)x(1)x(0)=exp(r+σε1),x(2)x(1)=exp(r+σε2),⋮x(n)x(n-1)=exp(r+σεn).
By a simplification, we obtain
(2.10)x(n)x(0)=exp(nr+σ∑i=1nεi),
that is,
(2.11)x(n)=x(0)exp(nr+σ∑i=1nεi).Lemma 2.4.
IfE(x(n)) is the mean-value function of the solution of system (1.3), then one has
(2.12)E[x(n)]=x(0)exp[n(r+σ22)].Proof.
One has(2.13)E[x(n)]=E[x(0)exp(nr+∑i=1nεi)]=E[x(0)exp(nr)exp(σ∑i=1nεi)]=x(0)exp(nr)E[exp(σ∑i=1nεi)].
Let  Y=∑i=1nεi. Since  εi~N(0,1),i=1,2,…,n, we have
(2.14)Y~N(0,n),
and the probability density function of random variable Y is
(2.15)f(y)=12πnexp(-y22n).
It follows from (2.13) and (2.15) that we have
(2.16)E[x(n)]=x(0)exp(nr)E[exp(σY)]=x(0)exp(nr)∫-∞+∞12πnexp(-y22n)dy=x(0)exp(nr)12πn∫-∞+∞exp(σy--y22n)dy=x(0)exp(nr)12πn∫-∞+∞‍exp[(y-σn)2-(σn)22n]dy=x(0)exp(nr)exp(nσ2/2)2πn∫-∞+∞exp[(y-σn)22n]dy=x(0)exp(nr)exp(nσ2/2)2π∫-∞+∞exp[((y-σn)/n)22]d(y-σnn)=x(0)exp[n(r+σ22)].Lemma 2.5.
Let the life period of pest population beτ, let  k1 be the loss caused by unit number pest in one generation, and let  k2 be the expending for annihilating pests once time. The loss caused by pest can be expressed as
(2.17)S(N0)=k1x(0)1-exp(r+(σ2/2))[1-exp(N0(r+σ22))].Proof.
Consider(2.18)S(N0)=k1τ∑n=0N0x(n)=k1τ{x(0)+x(0)exp(r+σ22)+x(0)exp2(r+σ22)+⋯+x(0)expN0(r+σ22)}=k1τx(0)1-exp(N0(r+(σ2/2)))1-exp(r+(σ2/2))=k1τx(0)1-exp(r+(σ2/2))[1-exp(N0(r+σ22))].Lemma 2.6.
LetP{N0=k}=Pk,  One has
(2.19)Pk=[1-Φ(ln(U/x(0))-krσk)]∏n=1k-1Φ(ln(U/x(0))-nrσn),k=1,2,…,n0,
where
(2.20)Φ(x)=∫-∞x‍12πexp(-t22)dt.Proof.
By the definition ofN0, we have
(2.21)x(N0)≥U,x(n)<U,n=0,1,…,N0-1.
Then
(2.22)P{N0=k}=P{x(k)≥U,x(n)<U,(n=0,1,…,k-1)}=P{x(k)≥U}∏n=1k-1P{x(n)<U}.
By(H2), we have
(2.23)P{x(0)<U}=1.
Furthermore,
(2.24)P{x(k)≥U}=1-P{x(k)<U}.
By (2.24), we have
(2.25)P{N0=k}=[1-P{x(k)<U}]∏n=1k-1P{x(n)<U}.
Moreover, one has
(2.26)P{x(n)<U}=P{x(0)exp(nr+σ∑i=1nεi)<U}=P{exp(nr+σ∑i=1nεi)<Ux(0)}=P{(nr+σ∑i=1nεi)<lnUx(0)}=P{∑i=1nεi<ln(U/x(0))-nrσ}=P{∑i=1nεin<ln(U/x(0))-nrσn}=Φ(ln(U/x(0))-nrσn),
and then we obtain
(2.27)Pk=[1-Φ(ln(U/x(0))-krσk)]∏n=1k-1Φ(ln(U/x(0))-nrσn),k=1,2,…,n0.Lemma 2.7.
The mean-value function of the loss caused by pest population is(2.28)E[S(N0)]=∑k=1n0{[1-exp(k(r+σ22))][1-Φ(ln(U/x(0))-krσk)]×∏n=1k-1‍Φ(ln(U/x(0))-nrσn)}.Proof.
By the definition of mean value function, we have(2.29)E[S(N0)]=∑k=1n0{[1-exp(k(r+σ22))]P{N0=k}},
then by Lemma 2.6, we obtain
(2.30)E[S(N0)]=∑k=1n0{[1-exp(k(r+σ22))][1-Φ(ln(U/x(0))-krσk)]×∏n=1k-1‍Φ(ln(U/x(0))-nrσn)}.Lemma 2.8.
The following equality holds(2.31)P{M>U}=∑k=1n0-1{[1-Φ(ln(U/x(0))-krσk)]∏n=1k-1Φ(ln(U/x(0))-nrσn)}.Proof.
By the definitions ofM and U, we have
(2.32)P{M>U}=P{N0<n0}=∑k=1n0-1P{N0=k}=∑k=1n0-1{[1-Φ(ln(U/x(0))-krσk)]∏n=1k-1Φ(ln(U/x(0))-nrσn)}.


## 3. Main Results

In this section, we give three main results. We first give the optimal control threshold value.Theorem 3.1.
If the assumptions(H1)–(H4) are satisfied, then the optimal control threshold value of system (1.3) is the minimal nonnegative solution of the following equation about U(3.1)k1x(0)∑k=1n0p+k2(1-exp(r+(1/2)σ2))∑k=1n0-1pexp(r+(1/2)σ2)-1=0,
where
(3.2)p1=2-k∏n=1k-1[1+f(22nσln(U-nrx(0)x(0)))],p2=exp(-12kσ2)(ln(U-krx(0)x(0)))2,p3=kπσ(krx(0)-U),p=-2p1p2p3,f(x)=2π∫0x‍e-t2dt.Proof.
By Lemmas2.3–2.8, we obtain that the total loss can be expressed as
(3.3)J(U)=E[S(N0)]+k2P(M>U)=∑k=1n0{[1-exp(k(r+σ22))][1-Φ(ln(U/x(0))-krσk)]∏n=1k-1Φ(ln(U/x(0))-nrσn)}+k2∑k=1n0-1{[1-Φ(ln(U/x(0))-krσk)]∏n=1k-1Φ(ln(U/x(0))-nrσn)}.
A calculation leads to
(3.4)J′(U)=k1x(0)∑k=1n0p+k2(1-exp(r+(1/2)σ2))∑k=1n0-1pexp(r+(1/2)σ2)-1.
Denote U* is the minimal nonnegative solution of the above equation, it follows from (3.4) that J′(U*)=0 and U* is the optimal control threshold value of system (1.3). The proof of Theorem 3.1 is complete.In the following, we give the optimal control moment.Theorem 3.2.
If the assumptions(H1)–(H4) hold, then the optimal control moment of system (1.3) can be expressed as
(3.5)T0=τ∑k=1n0{k[1-Φ(ln(U*/x(0))-krσk)]∏n=1k-1Φ(ln(U*/x(0))-nrσn)}.
where U* is defined in Theorem 3.1 and τ is the life period of the pest population.Proof.
By the definition ofN0, we have T0=E[N0τ]. Furthermore, it follows from Lemma 2.6 that
(3.6)T0=τ∑k=1n0kP{N0=k}=τ∑k=1n0{k[1-Φ(ln(U*/x(0))-krσk)]∏n=1k-1Φ(ln(U*/x(0))-nrσn)}.
The proof of Theorem 3.2 is complete.Finally, we give the estimate of the maximum likelihood estimations of the parametersr and σ of system (1.3).Theorem 3.3.
Letr^ and σ^ be the maximum likelihood estimations of the parameters r and σ, one has
(3.7)r^=y(n)-y(0)n,σ^=[∑i=1n[y(i)-y(i-1)-r^]2n]1/2,
where y(n)=lnx(n).Proof.
From system (1.3), we have
(3.8)lnx(n)-lnx(n-1)=r+σεn,
let y(n)=lnx(n), then we obtain
(3.9)y(n)-y(n-1)=r+σεn.
Since εn i.i.d N(0,1), we have [y(n)-y(n-1)] i.i.d N(r,σ2). Let x(n) be the quantity of the nth generation pest population, we can obtain corresponding values y(0),y(1),…,y(n), then the likelihood function of parameters  r and  σ is
(3.10)L(r,σ)=∏i=1n12πσexp(-[y(i)-y(i-1)-r]22σ2)=(12πσ)nexp(-12σ2)∑i=1n[y(i)-y(i-1)-r]2.
Further, we have
(3.11)lnL(r,σ)=-nln2πσ-12σ2∑i=1n[y(i)-y(i-1)-r]2.
From (3.11), we obtain the following likelihood equation
(3.12)∂lnL(r,σ)∂σ=-nσ+∑i=1n[y(i)-y(i-1)-r]2σ3=0,∂lnL(r,σ)∂r=y(n)-y(0)-nrσ2=0,
and the maximum likelihood estimations of r and σ are
(3.13)r^=y(n)-y(0)n,σ^=[∑i=1n[y(i)-y(i-1)-r^]2n]1/2.
The proof of Theorem 3.3 is complete.


## 4. An Example

In this section, to illustrate the feasibility of our theoretical results, we will give the following example.Example 4.1.
Consider the following system(4.1)x(n)=x(n-1)exp(0.3+0.08εn).
The choose the loss caused by the unit number pest k1=0.8, the expending for annihilating pest once k2=0.2, and initial value x(0)=0.01, n0=5. By Theorems 3.1 and 3.2, we can obtain the approximates of the optimal control threshold value U*=0.1168326 and the optimal control moment T0=3.3217.Next, we give the MLE of the parametersr and σ to compare the true value with estimation. In Table 1, for the given true value of parameters r=0.3 and σ=0.1, the number of the sample “size n” increases from 500 to 2000, the data of the columns r-MLE and σ-MLE are obtained by the average of 10 MLEs from the data coming from system (1.3). The columns of AE shows the absolute error of MLE. Table 1 shows that, with the augment of the number of the sample, the absolute error of MLE of r and σ will decrease, which implies that it is reasonable to estimate the parameters of system (1.3) by MLE.Table 1
The average value and absolute error of MLE of parameters with different number of sample.
True
Size
Aver
AE
(
r
,
σ
)
n
r-MLE
σ-MLE
r
σ
500
0.3049901
0.1000706
0.0049901
0.0000706
(
0.3
,
0.1
)
1000
0.3048290
0.0999396
0.0048290
0.0000604
2000
0.3040423
0.1000054
0.0040423
0.0000054



---

*Source: 676871-2012-10-17.xml*