# Parameter Estimation of the Lotka–Volterra Model with Fractional Order Based on the Modulation Function and Its Application

Authors: Ying Hao; Mingshun Guo
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The Lotka–Volterra model is widely applied in various fields, and parameter estimation is important in its application. In this study, the Lotka–Volterra model with universal applicability is established by introducing the fractional order. Modulation function is multiplied by both sides of the Lotka–Volterra model, and the model is converted into linear equations with parameters to be estimated by the fractional integration method. The parameters are obtained by solving the equations. The state of the system is estimated by shifted Chebyshev polynomial. Last, the implementation program of the model is compiled. The concrete implementation method of the improved model is proposed by an example in this study.


---

## Body

## 1. Introduction

In the 1940s, Lotka (1925) and Volterra (1926) jointly laid the theoretical foundation for the Lotka–Volterra model. As an extension of the logistic model, this model was used to study the competition and cooperation between different populations in the process of biological evolution. The Lotka–Volterra model is widely applied in the financial industry, resource development, Internet of Things, currency transactions, national economy, food, and other fields [1–8].To make the Lotka–Volterra model better adapted to applications, a stochastic Lotka–Volterra model disturbed by G-Brownian motion in the framework of nonlinear expectation was proposed by He et al. [9], which considered the uncertainty of the noise variance. Vadillo [10] studied this problem using three different stochastic models built from a single Lotka–Volterra deterministic model. More concretely, they study their mean-extinction time which satisfies the backward Kolmogorov differential equation, a linear second-order partial differential equation with variable coefficients. An extended Lotka–Volterra model was proposed by Stefan et al. [11], and they studied the microbial dynamics of two interactions in the chemostat. Tahara et al. [12] demonstrated that the Lotka–Volterra system can be stabilized by inducing small immigration of the prey or predator population, and they also checked the effect of the inclusion of the nonlinear interaction term to the stability of the Lotka–Volterra system. Stomov et al. [13] extended the existing N-species impulsive competitive models to an uncertainty case, and the main results are obtained by using Lyapunov-type functions and a comparison principle. Gavina et al. [14] built a Lotka–Volterra competition equation with a nonlinear crowding effect.Estimating the unknown parameters of the Lotka–Volterra model is a critical issue. A Lotka–Volterra competitive model parameter estimation method based on the grey direct modeling method was proposed by Whu and Wanghai-Jun [15]. Using the least squares method, parameter estimation of a stochastic Lotka–Volterra model driven by small Levy noises from discrete observations was proposed by Wei [16]. A polynomial approximation method for the solutions of continuous-time delayed population models was presented by Yuzbasi and Karasayir [17]. They transformed nonlinear algebraic equations for the solution and proposed a technique of residual correction. A regularization method was applied to parameter inversion of the Lotka–Volterra cooperative system with the periodic boundary conditions by Xu et al. [18]. The unknown parameters of the Lotka–Volterra system can be estimated using the approximate Bayesian computation method by Skvortsov et al. [19]. Perturbed Lyapunov function methods were used to solve multispecies Lotka–Volterra models with regime switching modulated by a continuous-time Markov chain involving a small parameter in Wang et al.’s study [20]. The Lindstedt–Poincare method was applied to the Lotka–Volterra model by Amore and Fernandez [21]. They utilized an efficient systematic algorithm, there were an unprecedented number of perturbation corrections for the two dynamical variables, and the frequency was obtained.Empirical formulas for complex dynamical processes in some fields are often expressed in the form of power-law functions, and the standard “gradient” law is not satisfied by the corresponding mechanical constitutive relation. Mechanical and physical processes with historical memory and spatial global correlation can be described by fractional differential operators concisely and accurately. In this study, a universal applicability model is established by implanting a fractional-order model into the Lotka–Volterra model. The antinoise property of the algorithm is improved by introducing integrals into the algebraic formula of the variable of the modulation function. The algebraic expression of the algorithm makes it nonasymptotically convergent. Unknown parameter estimation in the Lotka–Volterra model is accurately realized.


## 2. Preliminary Knowledge

### 2.1. Fractional Calculus

The following notations and definitions are performed before the parameter and system state estimation proof.Mark:Rn represents the n-dimensional Euclidean space, and Rm×n represents the m×n set of real matrices. Given the column vector x=x1,x2,…,xnT∈Rn, the upper corner “T” represents the transpose. Given that A∈Rn×n, 1n=1,1,…,1T∈Rn. Besides, diagd1,d2,…,dn∈Rn×n represents an n-order diagonal matrix, and detA represents the determinant of matrix A. Given B∈Rm×n and C=Rp×q, B⊗C∈Rmp×nq represents the Kronecker product of matrix B and matrix C. Interval I=0,h∈R, α∈R+. Moreover, l=⌊α⌋+1, ⌊α⌋ and ⌈α⌉ represent the largest integer which is less than or equal to α and the smallest integer which is more than or equal to α, respectively. The d square integrable space of interval 0,ξ is ℓ20,ξ.Definition 1.
The left-sided Riemann–Liouville fractional integral of a functiong⋅ is defined by(1)Dt0ft≔ft,Dtαft≔1Γα∫0tfγt−γα−1dγ,where Γ⋅ is the gamma function.
It should be noted thatDtαft=Dt−αf is also used to refer to the −α-order Riemann–Liouville fractional derivative of f.Definition 2.
Theα-order left-sided Riemann–Liouville fractional derivative of a function f⋅ is defined by(2)Dtαft≔dldtlDtα−lft.Definition 3.
Theα-order left-sided Caputo fractional derivative of a function is defined by(3)CDtαft≔Dtα−lflt.Definition 4.
The right-sided Riemann–Liouville fractional integral of a functionf⋅ is defined by(4)Dt,h0ft≔ft,Dt,hαft≔−1lΓ1−α∫0hflγγ−tα+1−ldγ.Definition 5.
Theα-order right-sided Riemann–Liouville fractional derivative of a function f⋅ is defined by(5)Dt,hαft≔−1ldldtlDt,hα−lft.Definition 6.
Theα-order right-sided Caputo fractional derivative of a function g⋅ is defined by(6)DCtαft≔Dtα−lflt.Lemma 1.
For any interval0,h∈R, the following formula holds [22]:(7)∫0hgtDtαftdt=∫0hDCt,hαgtftdt+∑k=0α−1−1kgktDtα−1−kftt=0t=h.


### 2.2. Modulation Function

Definition 7.
Givena,b∈R, L∈N∗, if k=0,1,…,L−1, then function g⋅ meets the criteria C1g∈CLa,b, C2gka=0, C3gkb=0, and the function g⋅ is called an L-order modulation function on the interval a,b.


### 2.3. Shifted Chebyshev Polynomial

Shifted Chebyshev polynomial is defined by the following recursive formula:(8)Tξ,0t=1,Tξ,1t=t,Tξ,i+1t=22tξ−1Tξ,it−Tξ,i−1t,i=1,2,….Equation (8) is equivalent to the following form:(9)Tξ,i=i∑k=0i−1i−ki+k−1!22ki−k!2k!ξkxk.Here,Tξ,i=−1i and Tξ,iξ=1.Anyxt∈ξ20,ξ could be represented as a shifted Chebyshev polynomial form as follows:(10)xt=∑j=0∞cjTξ,jt,where cj=1/hj∫0ξutTξ,jtωξtdt, j=1,2,⋯. Especially, the preceding N+1− paragraph is truncated, and the following can be obtained:(11)x˜t=∑k=0NckTξ,kt.


## 2.1. Fractional Calculus

The following notations and definitions are performed before the parameter and system state estimation proof.Mark:Rn represents the n-dimensional Euclidean space, and Rm×n represents the m×n set of real matrices. Given the column vector x=x1,x2,…,xnT∈Rn, the upper corner “T” represents the transpose. Given that A∈Rn×n, 1n=1,1,…,1T∈Rn. Besides, diagd1,d2,…,dn∈Rn×n represents an n-order diagonal matrix, and detA represents the determinant of matrix A. Given B∈Rm×n and C=Rp×q, B⊗C∈Rmp×nq represents the Kronecker product of matrix B and matrix C. Interval I=0,h∈R, α∈R+. Moreover, l=⌊α⌋+1, ⌊α⌋ and ⌈α⌉ represent the largest integer which is less than or equal to α and the smallest integer which is more than or equal to α, respectively. The d square integrable space of interval 0,ξ is ℓ20,ξ.Definition 1.
The left-sided Riemann–Liouville fractional integral of a functiong⋅ is defined by(1)Dt0ft≔ft,Dtαft≔1Γα∫0tfγt−γα−1dγ,where Γ⋅ is the gamma function.
It should be noted thatDtαft=Dt−αf is also used to refer to the −α-order Riemann–Liouville fractional derivative of f.Definition 2.
Theα-order left-sided Riemann–Liouville fractional derivative of a function f⋅ is defined by(2)Dtαft≔dldtlDtα−lft.Definition 3.
Theα-order left-sided Caputo fractional derivative of a function is defined by(3)CDtαft≔Dtα−lflt.Definition 4.
The right-sided Riemann–Liouville fractional integral of a functionf⋅ is defined by(4)Dt,h0ft≔ft,Dt,hαft≔−1lΓ1−α∫0hflγγ−tα+1−ldγ.Definition 5.
Theα-order right-sided Riemann–Liouville fractional derivative of a function f⋅ is defined by(5)Dt,hαft≔−1ldldtlDt,hα−lft.Definition 6.
Theα-order right-sided Caputo fractional derivative of a function g⋅ is defined by(6)DCtαft≔Dtα−lflt.Lemma 1.
For any interval0,h∈R, the following formula holds [22]:(7)∫0hgtDtαftdt=∫0hDCt,hαgtftdt+∑k=0α−1−1kgktDtα−1−kftt=0t=h.


## 2.2. Modulation Function

Definition 7.
Givena,b∈R, L∈N∗, if k=0,1,…,L−1, then function g⋅ meets the criteria C1g∈CLa,b, C2gka=0, C3gkb=0, and the function g⋅ is called an L-order modulation function on the interval a,b.


## 2.3. Shifted Chebyshev Polynomial

Shifted Chebyshev polynomial is defined by the following recursive formula:(8)Tξ,0t=1,Tξ,1t=t,Tξ,i+1t=22tξ−1Tξ,it−Tξ,i−1t,i=1,2,….Equation (8) is equivalent to the following form:(9)Tξ,i=i∑k=0i−1i−ki+k−1!22ki−k!2k!ξkxk.Here,Tξ,i=−1i and Tξ,iξ=1.Anyxt∈ξ20,ξ could be represented as a shifted Chebyshev polynomial form as follows:(10)xt=∑j=0∞cjTξ,jt,where cj=1/hj∫0ξutTξ,jtωξtdt, j=1,2,⋯. Especially, the preceding N+1− paragraph is truncated, and the following can be obtained:(11)x˜t=∑k=0NckTξ,kt.


## 3. Problem Description and Main Conclusions

### 3.1. Problem Description

Assume that multiple species in the same natural environment must affect the survival and development of each other. The fractional-order model is introduced to establish the general applicability Lotka–Volterra model in this study, which is as follows:(12)Dtαxit=xitbi+∑j=1naijxjt,(13)xit=θit,t∈0,h,where t∈I⊂R+U0, α is the fractional order, 0α1, xijt is the system state, θ· is the known function, θ refers to the initialization function of the state x, and aij and bi are the system parameters to be estimated.Furthermore, the numerical estimates ofxit in t∈h,+∞ need to be given, i=1,2,…,n.


### 3.2. Parameter and System State Estimation

As a robust and fast estimation method, the modulating functions method has been widely used in the fields of signal processing and automatic control for parameter identification [23], derivative estimation [24, 25], and integral estimation [26]. Recently, it has been extended to state estimation problems for fractional-order systems [27–30]. However, the models considered in previous works are linear. In this study, the modulating functions method is applied for both parameter and state estimation in the nonlinear fractional-order Lotka–Volterra model. First, according to the known initial value function, the unknown parameters aij and bi of system (12) are estimated, where i,j=1,2,…,n. On this basis, the approximate estimation in t∈h,+∞ is given.


#### 3.2.1. Parameter Estimation

Before the conclusion is given,a¯i=ai1,ai2,…,ain,biT,i=1,2,…,n. Since system (12) is known at t∈0,h, estimation of parameters aij and bi can be transformed into a linear system of equations about these parameters by the modulation function method. Then, a¯ii=1,2,…,n can be obtained by selecting the appropriate modulation function.Proposition 1.
Given the known state of system (12), xit, t∈0,h, and if there is a modulation function gi,mt, then(14)det∫0hGitxitXtdt≠0.
The system parametersaij and bi can be given by the following formula:(15)a¯i=∫0hGitxitXtdt−1∫0hDCt,tαgiTxitdt,where xt=x1t,x2t,…,xntT, x¯t=x1t,x2t,⋯,xnt,1T, Xt=1n+1⊗x¯Tt, git=gi,1t,gi,2t,…,gi,nt,gi,n+1tT, Git=diaggi,1t,gi,2t,…,gi,nt,gi,n+1t, i,j=1,2,…,n, and n,m=1,2,⋯,n+1.Proof.
Multiplygi,mt in both sides of equation (12) and the integrals from 0 to h. Obtain i=1,2,… and n,m=1,2,…,n+1.(16)∫0hgi,mtDtαxitdt=∫0hgi,mtxitbi+∑j=1naijxjtdt.
According to Lemma1, the left side of equation (16) becomes(17)∫0hgi,mtDt,hαgi,mtxtdt+∑k=0α−1−1kgi,mktDtα−1−kxitt=0t=h.
According to Definition7, from formula (17), the following equation can be obtained:(18)∫0hgi,mtDtαxitdt=∫0hDt,hαCgi,mtxitdt,where I=1,2,… and n,m=1,2,…,n+1. Therefore, according to formula (18),(19)∫0hGitxitXtdta¯i=∫0hCDt,hαgi,mTtxitdt.
Thus, Proposition1 can be proved.


#### 3.2.2. State Estimation

Suppose(20)C=c10c20⋯cn0c11c21⋯cn1⋮⋮⋱⋮c1Nc2N⋯cnN.According to the Chebyshev polynomial,x˜it=∑k=0NcikTξ,kt can be used to estimate the corresponding system state xit,i=1,2,… in t∈h,+∞ approximately, that is, estimates of the system state xit, t∈h,+∞, can be obtained by solving the algebraic equation C (cik can be obtained based on the known system state in t∈0,h).Proposition 2.
If there areN+1 time points, then t0,t1,⋯,tN∈0,h becomes(21)detTξ≠0,where Tξ=Tξ,0t0Tξ,1t0⋯Tξ,Nt0Tξ,0t1Tξ,1t1⋯Tξ,Nt1⋮⋮⋱⋮Tξ,0tNTξ,1tN⋯Tξ,NtN.
Then, the state of system (12) xit, t∈h,+∞ can be approximately given by x˜it=∑k=0NcikTξ,kt, and ciki=1,2,…,n,k=1,2,…,N can be given by the following formula:(22)C=Tξ−1θ,where θ=θ1t0θ2t0⋯θnt0θ1t1θ2t1⋯θnt1⋮⋮⋱⋮θ1tNθ2tN⋯θntN.Proof.
According to formula (11), the numerical solution of the system state xit, t∈h,+∞ can be obtained only by solving the corresponding coefficient ciki=1,2,…,n,k=1,2,…,N.
According to the known conditions and formula (13), the following can be obtained:(23)xitq=θitq,q=0,1,…,N.
The firstN+1 term x˜it is used to approximate xit, and then, system (23) is(24)∑k=0NcikTξ,ktq=θitq,q=0,1,…,N,that is,(25)TξC=θ.
Thus, according to formula (21), the proposition is proved.


## 3.1. Problem Description

Assume that multiple species in the same natural environment must affect the survival and development of each other. The fractional-order model is introduced to establish the general applicability Lotka–Volterra model in this study, which is as follows:(12)Dtαxit=xitbi+∑j=1naijxjt,(13)xit=θit,t∈0,h,where t∈I⊂R+U0, α is the fractional order, 0α1, xijt is the system state, θ· is the known function, θ refers to the initialization function of the state x, and aij and bi are the system parameters to be estimated.Furthermore, the numerical estimates ofxit in t∈h,+∞ need to be given, i=1,2,…,n.


## 3.2. Parameter and System State Estimation

As a robust and fast estimation method, the modulating functions method has been widely used in the fields of signal processing and automatic control for parameter identification [23], derivative estimation [24, 25], and integral estimation [26]. Recently, it has been extended to state estimation problems for fractional-order systems [27–30]. However, the models considered in previous works are linear. In this study, the modulating functions method is applied for both parameter and state estimation in the nonlinear fractional-order Lotka–Volterra model. First, according to the known initial value function, the unknown parameters aij and bi of system (12) are estimated, where i,j=1,2,…,n. On this basis, the approximate estimation in t∈h,+∞ is given.


### 3.2.1. Parameter Estimation

Before the conclusion is given,a¯i=ai1,ai2,…,ain,biT,i=1,2,…,n. Since system (12) is known at t∈0,h, estimation of parameters aij and bi can be transformed into a linear system of equations about these parameters by the modulation function method. Then, a¯ii=1,2,…,n can be obtained by selecting the appropriate modulation function.Proposition 1.
Given the known state of system (12), xit, t∈0,h, and if there is a modulation function gi,mt, then(14)det∫0hGitxitXtdt≠0.
The system parametersaij and bi can be given by the following formula:(15)a¯i=∫0hGitxitXtdt−1∫0hDCt,tαgiTxitdt,where xt=x1t,x2t,…,xntT, x¯t=x1t,x2t,⋯,xnt,1T, Xt=1n+1⊗x¯Tt, git=gi,1t,gi,2t,…,gi,nt,gi,n+1tT, Git=diaggi,1t,gi,2t,…,gi,nt,gi,n+1t, i,j=1,2,…,n, and n,m=1,2,⋯,n+1.Proof.
Multiplygi,mt in both sides of equation (12) and the integrals from 0 to h. Obtain i=1,2,… and n,m=1,2,…,n+1.(16)∫0hgi,mtDtαxitdt=∫0hgi,mtxitbi+∑j=1naijxjtdt.
According to Lemma1, the left side of equation (16) becomes(17)∫0hgi,mtDt,hαgi,mtxtdt+∑k=0α−1−1kgi,mktDtα−1−kxitt=0t=h.
According to Definition7, from formula (17), the following equation can be obtained:(18)∫0hgi,mtDtαxitdt=∫0hDt,hαCgi,mtxitdt,where I=1,2,… and n,m=1,2,…,n+1. Therefore, according to formula (18),(19)∫0hGitxitXtdta¯i=∫0hCDt,hαgi,mTtxitdt.
Thus, Proposition1 can be proved.


### 3.2.2. State Estimation

Suppose(20)C=c10c20⋯cn0c11c21⋯cn1⋮⋮⋱⋮c1Nc2N⋯cnN.According to the Chebyshev polynomial,x˜it=∑k=0NcikTξ,kt can be used to estimate the corresponding system state xit,i=1,2,… in t∈h,+∞ approximately, that is, estimates of the system state xit, t∈h,+∞, can be obtained by solving the algebraic equation C (cik can be obtained based on the known system state in t∈0,h).Proposition 2.
If there areN+1 time points, then t0,t1,⋯,tN∈0,h becomes(21)detTξ≠0,where Tξ=Tξ,0t0Tξ,1t0⋯Tξ,Nt0Tξ,0t1Tξ,1t1⋯Tξ,Nt1⋮⋮⋱⋮Tξ,0tNTξ,1tN⋯Tξ,NtN.
Then, the state of system (12) xit, t∈h,+∞ can be approximately given by x˜it=∑k=0NcikTξ,kt, and ciki=1,2,…,n,k=1,2,…,N can be given by the following formula:(22)C=Tξ−1θ,where θ=θ1t0θ2t0⋯θnt0θ1t1θ2t1⋯θnt1⋮⋮⋱⋮θ1tNθ2tN⋯θntN.Proof.
According to formula (11), the numerical solution of the system state xit, t∈h,+∞ can be obtained only by solving the corresponding coefficient ciki=1,2,…,n,k=1,2,…,N.
According to the known conditions and formula (13), the following can be obtained:(23)xitq=θitq,q=0,1,…,N.
The firstN+1 term x˜it is used to approximate xit, and then, system (23) is(24)∑k=0NcikTξ,ktq=θitq,q=0,1,…,N,that is,(25)TξC=θ.
Thus, according to formula (21), the proposition is proved.


## 3.2.1. Parameter Estimation

Before the conclusion is given,a¯i=ai1,ai2,…,ain,biT,i=1,2,…,n. Since system (12) is known at t∈0,h, estimation of parameters aij and bi can be transformed into a linear system of equations about these parameters by the modulation function method. Then, a¯ii=1,2,…,n can be obtained by selecting the appropriate modulation function.Proposition 1.
Given the known state of system (12), xit, t∈0,h, and if there is a modulation function gi,mt, then(14)det∫0hGitxitXtdt≠0.
The system parametersaij and bi can be given by the following formula:(15)a¯i=∫0hGitxitXtdt−1∫0hDCt,tαgiTxitdt,where xt=x1t,x2t,…,xntT, x¯t=x1t,x2t,⋯,xnt,1T, Xt=1n+1⊗x¯Tt, git=gi,1t,gi,2t,…,gi,nt,gi,n+1tT, Git=diaggi,1t,gi,2t,…,gi,nt,gi,n+1t, i,j=1,2,…,n, and n,m=1,2,⋯,n+1.Proof.
Multiplygi,mt in both sides of equation (12) and the integrals from 0 to h. Obtain i=1,2,… and n,m=1,2,…,n+1.(16)∫0hgi,mtDtαxitdt=∫0hgi,mtxitbi+∑j=1naijxjtdt.
According to Lemma1, the left side of equation (16) becomes(17)∫0hgi,mtDt,hαgi,mtxtdt+∑k=0α−1−1kgi,mktDtα−1−kxitt=0t=h.
According to Definition7, from formula (17), the following equation can be obtained:(18)∫0hgi,mtDtαxitdt=∫0hDt,hαCgi,mtxitdt,where I=1,2,… and n,m=1,2,…,n+1. Therefore, according to formula (18),(19)∫0hGitxitXtdta¯i=∫0hCDt,hαgi,mTtxitdt.
Thus, Proposition1 can be proved.


## 3.2.2. State Estimation

Suppose(20)C=c10c20⋯cn0c11c21⋯cn1⋮⋮⋱⋮c1Nc2N⋯cnN.According to the Chebyshev polynomial,x˜it=∑k=0NcikTξ,kt can be used to estimate the corresponding system state xit,i=1,2,… in t∈h,+∞ approximately, that is, estimates of the system state xit, t∈h,+∞, can be obtained by solving the algebraic equation C (cik can be obtained based on the known system state in t∈0,h).Proposition 2.
If there areN+1 time points, then t0,t1,⋯,tN∈0,h becomes(21)detTξ≠0,where Tξ=Tξ,0t0Tξ,1t0⋯Tξ,Nt0Tξ,0t1Tξ,1t1⋯Tξ,Nt1⋮⋮⋱⋮Tξ,0tNTξ,1tN⋯Tξ,NtN.
Then, the state of system (12) xit, t∈h,+∞ can be approximately given by x˜it=∑k=0NcikTξ,kt, and ciki=1,2,…,n,k=1,2,…,N can be given by the following formula:(22)C=Tξ−1θ,where θ=θ1t0θ2t0⋯θnt0θ1t1θ2t1⋯θnt1⋮⋮⋱⋮θ1tNθ2tN⋯θntN.Proof.
According to formula (11), the numerical solution of the system state xit, t∈h,+∞ can be obtained only by solving the corresponding coefficient ciki=1,2,…,n,k=1,2,…,N.
According to the known conditions and formula (13), the following can be obtained:(23)xitq=θitq,q=0,1,…,N.
The firstN+1 term x˜it is used to approximate xit, and then, system (23) is(24)∑k=0NcikTξ,ktq=θitq,q=0,1,…,N,that is,(25)TξC=θ.
Thus, according to formula (21), the proposition is proved.


## 4. Example Analysis

The survival and development of different populations in the same natural environment must influence each other. The degree of interaction among multiple populations can be analyzed quantitatively by the Lotka–Volterra model under the premise of a certain population size (n). The Lotka–Volterra model with fractional-order multicompetition relations in this study has more general application. The corresponding examples are discussed for n = 3 and n = 2, respectively.Example 1.
The Lotka–Volterra competition model is a dynamic model, and the cooperative relationship between enterprises can be reflected by this model [31]. Three companies with more invention patents in the new energy vehicle industry are chosen by this study. The number of patent holdings in the future can be predicted by the Lotka–Volterra model with the fractional order. In order to meet the requirements of the algorithm for sample data, the patent dates of enterprises are fitted to make the data continuous. Equation (12) can be written as(26)Dtαxit=xitbi+∑j=13aijxjt,i=1,2,3,where xit is the patent number of enterprise i, xjt is the patent number of enterprise j, bi is the change rate of each enterprise’s own patents, and aij is the influence intensity of enterprise j on enterprise i.
The distribution of the patent holdings of the three enterprises in the past five years is given in Table1.
According to the parameter identification (as shown in Figure1), the number of patent granted to the Zhejiang Geely Holding Group Co., Ltd. is influenced mostly by BYD Company Limited, and the number of patent granted to Beijing New Energy Automobile Co., Ltd. is influenced least by BYD Company Limited. The number change of patent granted of the three enterprises is less affected by the enterprises themselves and more affected by other enterprises in the same system. In addition, a12<0, a13<0, a21>0, and a31<0; it shows that in terms of patent granted number, BYD Company Limited and Zhejiang Geely Holding Group Co., Ltd. are in the predator-prey system and BYD Company Limited and Beijing New Energy Automobile Co. Ltd. are in the mutualism system.
Figure2 shows that the growth trend of the patent holdings of BYD Company Limited has a significant slowing trend while patent holdings of Zhejiang Geely Holding Group Co., Ltd. and Beijing New Energy Automobile Co., Ltd. are in a period of rapid growth after 10 years. Patent holdings of BYD Company Limited are influenced by another two enterprises increasingly seriously. Moreover, the patent holdings of Beijing New Energy Automobile Co. Ltd. will surpass Zhejiang Geely Holding Group Co., Ltd. in five years. The gap between patent holdings of the two enterprises will widen as time goes on.Table 1
Distribution of patent holdings in recent years.
Enterprise20142015201620172018BYD Company Limited13202589419955176067Zhejiang Geely Holding Group Co., Ltd.6521085179623383080Beijing New Energy Automobile Co., Ltd.219701156823302801Figure 1
Results of parameter identification.Figure 2
Forecasting results of enterprise patent quantity for the next 30 years.Example 2.
Research on parameter estimation of patent application number of two competitive enterprises is in order to reflect the relationship between patent applied number and patent granted number more objectively. The insufficiency of decision condition in Example1 is supplemented. Equation (12) can be written as(27)Dtαxit=xitbi+∑j=12aijxjt,i=1,2,where xit is the patent applied number of enterprise i, xjt is the patent applied number of enterprise j, bi is the patent application change rate of each enterprise, and aij is the influence intensity of enterprise j on enterprise i.
Patents applied number of the two enterprises from 2014 to 2018 is given in Table2.
According to the parameter identification (as shown in Figure3), in terms of the patent application number, BYD Company Limited and Zhejiang Geely Holding Group Co., Ltd. are in the predator-prey system. The change of the patent applied number of Zhejiang Geely Holding Group Co., Ltd. is influenced mostly by BYD Company Limited, while the influence of it owns much less. It shows that the innovation ability of Zhejiang Geely Holding Group Co., Ltd. is slightly lower than BYD Company Limited, and BYD Company Limited has a leading edge in technology.Table 2
Distribution of patent applications in recent years.
Enterprise20142015201620172018BYD Company Limited18431923319732092597Zhejiang Geely Holding Group Co., Ltd.1452859155512762800Figure 3
Results of parameters estimation.


## 5. Conclusion

In this study, the modulating functions method was applied to estimate unknown parameters, and the state of the Lotka–Volterra model was defined by a fractional-order nonlinear system. First, by multiplying the modulating functions to the considered model and applying fractional-order integration by parts formulas, the differential equations can be transformed into a set of integral equations where the unknown parameters are linear concerning the integrals. Then, the unknowns are estimated by solving the linear system of the obtained integral equations. Second, by taking the truncated expansion of the state based on the shifted Chebyshev polynomials, the unknown coefficients are obtained by solving a linear system using the initialization function of the state. Then, these estimated coefficients are used to predict the future behavior of the state.



---

*Source: 6645059-2021-07-14.xml*