# Dynamical Behavior of a Predator-Prey System Incorporating a Prey Refuge with Impulse Effect

Authors: Zuxiong Li; Xue Yang; Shengnan Fu
Journal: Complexity
 (2022)

---

## Abstract

A predator-prey model with Holling II functional response incorporating a prey refuge with impulse effect is considered in this paper. With the help of the Floquet theory of impulsive differential equations, local stability and global attractivity of the boundary periodic solution of the system are derived, and then sufficient conditions for global asymptotic stability of the boundary periodic solution are obtained. Next, the permanence of the system is proved by constructing a Lyapunov function. Further, by applying the bifurcation theory of impulsive differential equations, conditions under which the system has a positive periodic solution are obtained. Finally, numerical simulations are presented to illustrate the analytical results.


---

## Body

## 1. Introduction

The study of population dynamic system has always been the focus of scholars. Population dynamics behavior is always affected by many factors. Smith [1] pointed out that scholars often study the effects on the behaviors such as territorial behavior and migration. The challenge of reducing the negative environmental impacts of land use and maintaining economic and social benefits is considered [2]. Wildlife species are also affected by agricultural activities, such as farming, intercropping, drainage, rotation, grazing, and the widespread use of pesticides and fertilizers [3]. One of the more relevant behavioral traits that affect the dynamics of predator-prey systems is the use of spatial refuges by the prey. Some fractions of the prey population are partially protected against predators by using refuges. González-Olivares et al. [4] proposed a predator-prey system with Holling type II functional response incorporating a constant prey refuge:(1)dxtdt=αxt1−xtK−βxt−myt1+axt−m,dytdt=−dyt+cβxt−myt1+axt−m.They investigated the local stability of equilibria and the existence of limit cycle of the system. Chen et al. [5] discussed the instability and global stability properties of the equilibria and the existence and uniqueness of limit cycle of the above system. The predator-prey systems with prey refuge have been focused by several authors. Tripathi et al. studied the model with Beddington-DeAngelis type function response incorporating a prey refuge [6, 7], the model with time delay and prey refuge [8], and the model for prey with variable rates in protected areas [9]. Mondal and Samanta [10] investigated the model with prey refuge dependent on both species and constant harvest in predator [10] and the model with nonlinear prey refuge [11, 12]. Many models with prey refuge have been studied; one can also refer to [13–21] and the references cited therein.It has been noticed that many dynamical systems may be disturbed by impulsive factors. With profound understanding of human nature, the theory of impulsive differential equations [22–24] becomes more perfect. The impulsive differential equation has become a widely concerned subject in recent years, and it is more appropriate to apply to biological systems for the actual reality, for example [25–35].Shea and Amarasekare [33] pointed out that “conservation, harvesting, and pest control are three aspects of the same general problem: population management.” People seek to maintain exploited populations at productive levels by harvesting.The results show that the population model with refuge is practical. Many periodic factors, such as seasonal periodic change, migration, and harvest, have a great impact on population dynamics. The model with refuge and impulsive effect has also been studied by many scholars, one can refer to [13, 34, 35]. Driven by the above reasons, we consider system (1) with prey subject to periodic impulse harvesting, the model is as follows:(2)dxtdt=αxt1−xtK−βxt−myt1+axt−m,dytdt=−dyt+cβxt−myt1+axt−m,t≠nT,△xt=−pxt,△yt=0,t=nT.where xt and yt are the population density of prey and predator, respectively, α represents the intrinsic growth rate, K is the environment carrying capacity, β is the capture rate of predator, d is the death rate of the predator, c is the conversion factor denoting the number of newly born predators for each captured prey, m is a constant number of prey using refuges, which protect prey from predation, and 0<p<1 is impulse capture rate of prey. We assume that the population density of prey is more than the constant mm<xt, t, d, K, α, β, c, a, and m are all positive constants, △xt=xt+−xt, △yt=yt+−yt, and T is the period of the impulse effect.The rest of this paper is organized as follows. Some definitions and preliminary lemmas are introduced in Section2. Existence and global attractivity of boundary periodic solution of system (2) are discussed in Section 3. Section 4 contributes to permanence of system (2), and Section 5 aims at discussion of positive periodic solution. Section 6 contains some numerical simulations and illustrates our analytical results. Some simple discussions are given in Section 7.


## 2. Definitions and Preliminary Lemmas

LetR+=0,∞ and R+2=‖z∈R2z≥0. Denote the map f=f1,f2 defined by the right hand of system (2). Let V: R+×R+2⟶R+, then V is said to belong to class V0, if(1)
V is continuous in nT,n+1T and limt,u⟶nT,z,t>nTVt,u=VnT+,z(2)
V is locally Lipschitzian in xDefinition 1.
IfV∈V0, then for t,z∈nT,n+1T×R+2, the upper right derivative of Vt,x of system (2) is defined as D+Vt,z=limh⟶0sup1/hVt+h,z+hft,z−Vt,z.
It follows from the second equation of (2) that yt≥0 for y0≥0, so(3)dxtdt≥αxt1−xtK−βxtyt1+axt−m.
Therefore, the following Lemma1 is obvious.Lemma 1.
Ifxt is a solution of system (2) with x0+≥0, then xt≥0 for all t≥0.Definition 2.
System (2) is said to be permanent if there are constants m, M>0 (independent of initial value) and a finite time T0, such that m≤xt≤M and m≤yt≤M when t≥T0 for all solutions xt,yt, with all initial values x0+>0 and y0+>0.Lemma 2 (see [22]).
LetV:R+×Rn⟶R+ and V∈V0. Assume that(4)D+Vt,z≤ht,Vt,z,t≠nT,Vt,zt+≤ψnVt,z,t=nT,where h:R+×R+⟶R satisfies ψn:R+⟶R+ is nondecreasing; let rt be the maximal solution of the following scalar impulsive differential equation:(5)dutdt=ht,ut,t≠nT,ut+=ψnut,t=nT,u0+=u0≥0,existing on t0,+∞. Then, Vt0+,x0≤u0 implies that(6)Vt,zt≤rt,t≥t0,where zt=zt,t0,z0 is any solution of system (2) existing on t0,+∞.Now we first consider some basic properties of the following subsystem of system (2):(7)dxtdt=αxt1−xtK,t≠nT,△xt=−pxt,t=nT..Lemma 3 (see [36]).
System (7) has a positive periodic solution x∗t. Any solution xt of system (7) satisfies xt⟶x∗t as t⟶∞. It is easy to see that(8)x∗t=K1−p−exp−αTpexp−αt−nT+1−p−exp−αT,t∈nT,n+1T.Therefore,x∗t is the positive periodic solution of system (7).


## 3. Global Attractivity of Boundary Periodic Solution

First of all, we give the extinction of predator.Theorem 1.
Letxt,yt be any solution of system (2), if(9)T>−ln1−pα,−d−cβmqT+cβKαq2−Kaqlnqp+expαTqR+KRaqp+qR+KRa<0,then the stability of boundary periodic solution x∗t,0 is globally asymptotically stable, where q=1−am and R=1−p−exp−αT>0.Proof.
Firstly, we will prove the local stability of boundary periodic solution. Letxt=ut+x∗t and yt=vt, then the solution of variation equation of system (2) is(10)utvt=Φtu0v0,where Φt satisfies(11)dΦtdt=α−2αx∗tK−βx∗t−m1+ax∗t−m0−d+cβx∗t−m1+ax∗t−mΦt,and Φ0=I, in which I is the identity matrix. Hence, the fundamental solution matrix is(12)ΦT=exp∫0Tα−2αx∗tKdt∗0∫0T−d+cβx∗t−m1+ax∗t−mdt.
Then, the third equation and the forth equation of system (2) are linearized to be(13)unT+vnT+=1−p001unTvnT.
Stability of boundary periodic solution depends on the eigenvalue ofJ by(14)J=1−p001ΦT,and the eigenvalues of J are(15)μ1=1−pexp∫0Tα−2αx∗tKdt,μ2=exp∫0T−d+cβx∗t−m1+ax∗t−mdt.
According to the Floquet theory of impulsive equation, the boundary periodic solutionx∗t,0 is locally asymptotically stable if μ1<1 and μ2<1; that is,(16)T>−ln1−pα,−d−cβmqT+cβKαq2−Kaqlnqp+expαTqR+KRaqp+qR+KRa<0.
Secondly, we will prove the global attractivity of boundary periodic solutionx∗t,0. Let ϵ>0, so as to δ=∫0T−d+cβx∗t+ϵ−m/1+ax∗t+ϵ−mdt; that is.(17)−d−cβm−ϵ1−am−ϵT+cβKαq2−Kaqlnqp+expαTqR+KRaqp+qR+KRa<0,where q=1−am,R=1−p−exp−αT.
From the first equation of system (2), it is easy to obtain(18)dxtdt<αxt1−xtK.
Consider the following comparison system:(19)dztdt=αzt1−ztK,t≠nT,△zt=−pz,t=nT,z0+=x0+>0.
Whent is large enough, by Lemmas 2 and 3, we can obtain that(20)xt≤zt<x∗t+ϵ.
Fort≥0, we assume that xt≤zt<x∗t+ϵ, and then we have(21)dytdt≤−dy+cβx∗t−m+ϵy1+ax∗t−m+ϵ.
Integrating and solving the above system onnT,n+1T, we can derive that(22)yn+1T≤ynTexp∫nTn+T−d+cβx∗t−m+ϵ1+ax∗t−m+ϵdt.
Therefore,ynT≤y0+1−pnδ, then ynT⟶0 as n⟶∞. When t∈nT,n+1T, we can easily know that 0<ynT≤ynTH, so yt⟶0 as t⟶∞, where H=supu∈0,Texp∫0u−dcβx∗t−m+ϵ/1+ax∗t−m+ϵdt.
Next, we provext⟶x∗t as t⟶∞. By dxt/dt<αxt1−xt/K, we can obtain xt<x∗t. By Lemma 3, we can get(23)x∗t=K1−p−exp−αTpexp−αt−nT+1−p−exp−αT≤K1−p−exp−αT1−p1−exp−αT,t∈nT,n+1T.
Letm2=K1−p−exp−αT/1−p1−exp−αT, we have proved that xt<x∗t+ϵ for any ϵ>0 when t is large enough. Suppose that xt≤m2 as t>0. There exist T0>0 and m1>0 for t≥T0 such that yt≤m1. Assume that yt≤m1 as t>0, it follows by system (2)(24)dxtdt≥αxt1−xtK−βm1m2−m1+am2−m.
Consider the comparison system(25)dwtdt=αwt1−wtK−βm1m2−m1+am2−m,t≠nT,△wt=−pwt,t=nT,w0+=x0+>0.
By equation (25), we can obtain(26)w∗t=expMt−nT1−pQ11−expMT1−1−p1−expMT−Q1M+Q1M,where M=α1−m2/K=αpexp−αT/1−p1−exp−αT and Q1=βm1m2−m/1+am2−m.
By Lemmas2 and 3, one can obtain that wt⟶x∗t and wt⟶w∗t, as t⟶∞. So, for any ϵ2, there exists a T1, such that x∗t−ϵ2<xt<x∗t−ϵ2 as t>T1. Let ϵ2⟶0, then xt⟶x∗t as t⟶∞. The theorem is completely proved.


## 4. Permanence of System (2)

Theorem 2.
There exists a constantS>0 such that xt<S and yt<S for each solution xt,yt of system (2) when t is large enough.Proof.
Suppose thatxt,yt is any solution of system (2). Let V=cxt+yt and 0<λ<d; we calculate the upper right derivative of Vt along system (2) and get the following inequality:(27)D+Vt+λVt=cαxt1−xtK−dyt+cλxt+λyt=−cαx2tK+cα+λxt+λ−dyt≤Kcα+λ24α,t≠nT.
Therefore, it is bounded. We choose anM0 such that(28)D+Vt≤−λVt+M0.
It is easy to know that(29)VnT+=ynT+1−pxnT≤VnT.
According to Lemma2, we have(30)Vt≤V0+−M0λexp−λt+M0λ,t∈nT,n+1T.
Therefore,Vt is ultimately bounded. There exists a constant S>0 such that any solution xt,yt of system (2) satisfies xt≤S and yt≤S for t large enough.Theorem 3.
System (2) is permanent if(31)−d−cβm1−amT+cβKαq2−Kaqlnqp+expαTqR+KRaqp+qR+KRa>0,where q=1−am and R=1−p−exp−αT.Proof.
Supposext,yt is any solution of system (2); we have proved there exists an M such that xt≤M and yt≤M. Let h=K1−p−exp−αT/1−exp−αT−ϵ; we can know xt>h when t is large enough. We will find a constant m4 such that yt≥m4 when t is large enough, where m4=m3exp1+n2+n3Tδ2 and m3 is a constant. Next, we consider the following two cases. Firstly, there exists a sufficient small ϵ1 such that(32)δ1=∫0T−d+cv∗t−m−ϵ11+ax∗t−m−ϵ1dt.
We can proveyt<m3 is not satisfied for any t; otherwise,(33)dxtdt>αxt1−m2K−βm2−mm31+am2−m.
Consider the comparison system(34)dvtdt=αvt1−vtK−βm3m2−m1+am2−m,t≠nT,△vt=−pvt,t=nT,v0+=x0+>0.
It is not hard to getv∗t=expMt−nT1−pQ21−expMT/1−1−p1−expMT−Q2/M+Q2/M,where(35)M=α1−m2K=αpexp−αT1−p1−exp−αT,Q2=βm3m2−m1+am2−m.
Therefore, there exists at2>0 such that v∗t−ϵ1<v∗t≤xt for t>t2. Let n1∈ℕ and n1T≥t1 as n>n1. It is easy to know that(36)dytdt≥yt−d+cβv∗t−ϵ1−m1+av∗t−ϵ1−m.
Integrating and solving the above system onnT,n+1T, we can derive that(37)yn+1T≥ynTexp∫nTn+1T−d+cβv∗t−ϵ1−m1+av∗t−ϵ1−mdt=ynTexpδ1.
Then,j⟶∞,yn+jT≥ynTexpjδ1⟶∞; it contradicts with yt<m3. Thus, there is a t0>0 such that yt0≥m3.
Secondly, ift>t0, then yt≥m3. Otherwise, let t∗=inft>t0t|yt<m3, we can easily know yt≥m3 when t∈t1,t∗, choose n2,n3∈ℕ such that expn2+qδ2Texpn3δ1T>1, where δ2=−d+cβh−m/1+ah−m. Let T¯=n2+n3T, then there is a t2∈n1T,n1+1T+T¯ such that yt2≥m3. Otherwise, yt<m3 as t∈n1T,n1+1T+T¯. From system (25), we can obtain(38)yn1+1+n2+n−3T≥yn1+1+n2expn3δ1.as dyt/dt=−dyt+cβxt−myt/1+axt−m, so dy/dt≥−dyt+cβh−myt/1+ah−m.
Integrating and solving the above system ont∗,n1+n2+1T, we can derive that yn1+n2+1T≥m3expn2+1δ2T, so yn1+1+n2+n−3T≥m3expn2+1δ2Texpn3δ1T>m3, it contradicts with the assumption yt<m3. Let t¯=inft|yt≥m3,t>t∗, we have xt=m3exp1+n2+n3Tδ2≜m4. For t>t¯, as xt¯≥m3, we discuss it in the same way, at last we obtain yt≥m4, for any t≥t1.


## 5. Existence of Positive Periodic Solution

In this section, we will study the existence of positive periodic solution by bifurcation theory of impulsive differential equations. In order to be consistent with [37], we denote xt, yt as x1t, x2t, respectively. System (2) is(39)dx1tdt=αx1t1−x1tK−βx1t−mx2t1+ax1t−m≜F1x1t,x2t,dx2tdt=−dx2t+cβx1t−mx2t1+ax1t−m≜F2x1t,x2t,t≠nT,x1nT+=1−px1nT≜θ1x1nT,x1nT,x2nT+=x2nT≜θ2x1nT,x1nT,t=nT.By formally deriving the equation(40)ddtΦt,x0=FΦt,x0,which characterized the dynamics of the unperturbed flow associated to the first two equations in (2), we obtain that(41)ddtDxΦt,x0=DxFΦt,x0DxΦt,x0.This relation will be integrated in what follows in order to compute the components ofDxΦt,x0 explicitly. Firstly, it is clear that(42)Φt,x0=ϕ1t,x0,ϕ2t,x0.Then, we derive that(43)ddt∂ϕ1∂x1∂ϕ1∂x2∂ϕ2∂x1∂ϕ2∂x2=α−2αx∗tK−βx∗t−m1+ax∗t−m0−d+cβx∗t−m1+ax∗t−m∂ϕ1∂x1∂ϕ1∂x2∂ϕ2∂x1∂ϕ2∂x2t,x0.The initial condition att=0 being DxΦ0,x0=I. Here I is the identity matrix. It follows that(44)∂ϕ2t,x0∂x1=exp∫0t−d+cβx∗s−m1+ax∗s−mds∂ϕ20,x0∂x1.AsDxΦ0,x0=I, we can get ∂ϕ2t,x0/∂x1=0 for t≥0. We can obtain from (43) that(45)ddt∂ϕ1t,x0∂x1=α−2αx∗tK∂ϕ1t,x0∂x1,ddt∂ϕ1t,x0∂x2=α−2αx∗tK∂ϕ1t,x0∂x2ddt∂ϕ1t,x0∂x2−βx∗t−m1+ax∗t−m∂ϕ2t,x0∂x2,ddt∂ϕ2t,x0∂x2=−d+βx∗t−m1+ax∗t−m∂ϕ2t,x0∂x2.According to the initial condition, we obtain that(46)∂ϕ1t,x0∂x1=exp∫0tα−2αx∗sKds,∂ϕ1t,x0∂x2=∫0texp∫utα−2αx∗rKdr−βx∗u−m1+ax∗u−mexp∫0u−d+cβx∗r−m1+ax∗r−mdrdu,∂ϕ2t,x0∂x2=exp∫0t−d+βx∗r−m1+aβx∗r−mdr.From (4), we obtain that DxN0,0,0=I−DxΦT0,x0,which implies(47)DxN0,0,0=a0′b0′0d0′,where a0′, b0′, and d0′ are(48)a0′=1−1−pexpαT0−2ln1−p−exp−αTexpαT0+p1−exp−αT>0,b0′=p−1∫0T0exp∫uT0α−2αx∗rKdr−βx∗u−m1+ax∗u−mexp∫0u∫0T0−d+cβx∗r−m1+ax∗r−mdrdu>0,d0′=1−exp∫0T0−d+cβx∗r−m1+ax∗r−mdr,where T0 is the root of d0′=0.Firstly, we make sure the sign ofC, where(49)C=2∂θ2∂x2b0′a0′∂2ϕ2T0,x0∂x1∂x2−∂θ2∂x2∂2ϕ2T0,x0∂x22.For determiningC, we must calculate ∂2ϕ2T0,x0/∂x22 and ∂2ϕ2T0,x0/∂x1∂x2.We have(50)∂2ϕ2T0,x0∂x22=∫0T0exp∫uT0−d+cβx∗r−m1+ax∗r−mdrcβ1+ax∗u−m2exp∫uT0α−2αx∗rKdr−βx∗u−m1+a∗u−mexp∫0p∫uT0−d+cβx∗r−m1+ax∗r−mdrdpdu<0.Then, we have(51)∂2ϕ2T0,x0∂x1∂x2=∫0T0cβ1+ax∗u−m2duexp∫0T0−d+cβx∗r−m1+ax∗r−mdr>0.Therefore,C>0.Secondly, we make sure the sign ofB, where(52)B=−∂2θ2∂x1∂x2∂ϕ1T0,x0∂T¯+∂ϕ1T0,x0∂x11a0′∂ϕ1T0,x0∂T¯∂ϕ2T0,x0∂x2−∂θ2∂x2∂2ϕ2T0,x0∂T¯∂x2+∂2ϕ2T0,x0∂x1∂x21a0′θ1x1∂ϕ1T0,x0∂T¯.Therefore, we have to calculate∂2ϕ2T0,x0/∂T¯∂x2,(53)∂2ϕ2T0,x0∂T¯∂x2=−d+cβx∗T0−m1+ax∗T0−mexp∫0T0−d+cβx∗r−m1+ax∗r−mdr.Then, we have(54)B=−−d+cβx∗T0−m1+ax∗T0−mexp∫0T0−d+cβx∗r−m1+ax∗r−mdr.In order to make sure the sign ofB, let(55)ft=−d+cβx∗t−m1+ax∗t−m.Then,(56)f′t=cβ1+ax∗u−m2K1−p−exp−αTpαexp−αT1−p−exp−αT+pexp−αt2>0.Asd0′=0, we can obtain ∫0T0ftdt=0. Because ft is an increasing function, so fT0>0. According to equation, we can know B<0, thus BC<0. According to [37], the following conclusion can be derived.Theorem 4.
System (2) has a critical bifurcation; there exists a positive periodic solution of system (2) if t>T0, near T0.


## 6. Numerical Simulations

In this section, we investigate the effects of impulsive perturbations on system (2) by using numerical method to illustrate our theoretical results. For convenience, we assume keeping some parametric values of system (2) as α=0.5, K=4, c=0.5, β=0.4, m=2, a=0.2, and d=0.2. By simple calculation, we can get T0≈2.4023 for(57)d0′=1−exp∫0T0−d+cβx∗r−m1+ax∗r−mdr=0.Now, we choose different parameterT to illustrate our main results by numerical simulation. Let T=1.5, x0=3.3, and y0=0.05; bounded periodic solution x∗t,0 is globally asymptotically stable (see Figure 1); By Theorem 2, consider the comparison system, let T=3, x0=3.3, and y0=0.14, and one can get system (2) is permanent (see Figure 2). By Theorem 4, let T=2.42 nearly to T0, x0=3.3, and y0=0.005; system (2) has a critical bifurcation (Figure 3). Let T=2.38, x0=3.3, and y0=0.005; the number of predators population yt of system (2) decreases significantly (Figure 4). Further, take the smaller y0=0.00001, T=2.38, and x0=3.3, and we can obtain that the predator will be extinct (Figure 5).Figure 1
Thext time series and yt time series are simulated by numerical integration of system (2) with T=1.5, x0=3.3, and y0=0.05. Bounded periodic solution x∗t,0 is globally asymptotically stable.Figure 2
Thext time series and yt time series are simulated by numerical integration of system (2) with T=3 and x0=3.3y0=0.14. System (2) is permanent.Figure 3
Thext time series and yt time series are simulated by numerical integration of system (2) with T=2.42, nearly to T0, x0=3.3, and y0=0.005. System (2) has a critical bifurcation.Figure 4
Thext time series and yt time series of system (2) with T=2.38, x0=3.3, and y0=0.005.Figure 5
Thext time series and yt time series of system (2) with T=2.38, x0=3.3, and y0=0.00001.We simulate the influence of the refuge constantm on the system (see Figures 6 and 7) and the influence of pulse time T on the system (see Figure 8).Figure 6
Bifurcation diagram of system (2) affected by prey refuge constant number m with T=2.
(a)(b)Figure 7
Bifurcation diagram of system (2) affected by prey refuge constant number m with T=4.
(a)(b)Figure 8
Bifurcation diagram of system (2) affected by impulsive periodic T.
(a)(b)


## 7. Discussion

We have studied a predator-prey model with Holling II functional response incorporating a prey refuge and impulsive effect. Bounded periodic solutionx∗t,0 is globally asymptotically stable if T<T0. System (2) is permanent if T>T0. By Figures 1 and 2, we show the facts above. System (2) has a critical bifurcation; there exists a positive periodic solution of system (2) if t>T0, nearby T0. The facts above are shown in Figure 3. By Theorem 1, one can know that the predator is extinct when T is little smaller than T0, as shown in Figure 5. It is easily known that T0 is a threshold. From the view point of biology, if the harvest period T of the prey is too short T<T0, that is, the harvest frequency is relatively frequent, the number of the prey is bound to be small, which cannot guarantee the food supply of the predator, resulting in the final extinction of the predator. On the contrary, if the harvest period exceeds the threshold T>T0, the food supply of the predator is sufficient, and the system is persistent. This provides a good theoretical basis for how to select appropriate parameters to maintain ecological balance.The influence of refuge constantm on system (2) can be seen in Figures 6 and 7. This shows that the predator becomes extinct and the number of prey remains unchanged when the prey refuge constant reaches a certain value. This shows that when the number of prey entering the shelter reaches a certain amount, the prey that the predator can catch is less and insufficient, and the supply of prey is in short supply, resulting in the extinction of the predator. It is consistent with the reality. The impact of the change of pulse time T on system (2) is shown in Figure 8, this shows that the selection of pulse period T also has an obvious impact on the system. From b of Figure 8, we can see that the predator yt is persistent when T>T0≈2.4023.The effects of different pulse forms on the dynamic behavior of the system are also very different, so it is necessary for us to study the effects of these different pulse forms on the dynamic behavior of the system. The state feedback control of the system is also close to the real problem, we will study it in the following work.



---

*Source: 2422923-2022-07-14.xml*