Loading [MathJax]/jax/output/SVG/jax.js
Research article

On entire solutions of certain type of nonlinear differential equations

  • Received: 06 June 2020 Accepted: 15 July 2020 Published: 28 July 2020
  • MSC : 34M10, 30D05, 30D35

  • In this paper, we shall extend some results regarding the growth estimate of entire solutions of certain type of linear differential equations to that of nonlinear differential equations. Moreover, our results will include several known results for linear differential equations obtained earlier as special cases.

    Citation: Fengrong Zhang, Linlin Wu, Jing Yang, Weiran Lü. On entire solutions of certain type of nonlinear differential equations[J]. AIMS Mathematics, 2020, 5(6): 6124-6134. doi: 10.3934/math.2020393

    Related Papers:

    [1] Álvaro Rodríguez Abella, Melvin Leok . Discrete Dirac reduction of implicit Lagrangian systems with abelian symmetry groups. Journal of Geometric Mechanics, 2023, 15(1): 319-356. doi: 10.3934/jgm.2023013
    [2] Elias Maciel, Inocencio Ortiz, Christian E. Schaerer . A Herglotz-based integrator for nonholonomic mechanical systems. Journal of Geometric Mechanics, 2023, 15(1): 287-318. doi: 10.3934/jgm.2023012
    [3] Valentin Duruisseaux, Melvin Leok . Time-adaptive Lagrangian variational integrators for accelerated optimization. Journal of Geometric Mechanics, 2023, 15(1): 224-255. doi: 10.3934/jgm.2023010
    [4] Xavier Rivas, Daniel Torres . Lagrangian–Hamiltonian formalism for cocontact systems. Journal of Geometric Mechanics, 2023, 15(1): 1-26. doi: 10.3934/jgm.2023001
    [5] Marcin Zając . The dressing field method in gauge theories - geometric approach. Journal of Geometric Mechanics, 2023, 15(1): 128-146. doi: 10.3934/jgm.2023007
    [6] Robert I McLachlan, Christian Offen . Backward error analysis for conjugate symplectic methods. Journal of Geometric Mechanics, 2023, 15(1): 98-115. doi: 10.3934/jgm.2023005
    [7] Markus Schlarb . A multi-parameter family of metrics on stiefel manifolds and applications. Journal of Geometric Mechanics, 2023, 15(1): 147-187. doi: 10.3934/jgm.2023008
    [8] Jacob R. Goodman . Local minimizers for variational obstacle avoidance on Riemannian manifolds. Journal of Geometric Mechanics, 2023, 15(1): 59-72. doi: 10.3934/jgm.2023003
    [9] Jordi Gaset, Arnau Mas . A variational derivation of the field equations of an action-dependent Einstein-Hilbert Lagrangian. Journal of Geometric Mechanics, 2023, 15(1): 357-374. doi: 10.3934/jgm.2023014
    [10] Alexander Pigazzini, Cenap Özel, Saeid Jafari, Richard Pincak, Andrew DeBenedictis . A family of special case sequential warped-product manifolds. Journal of Geometric Mechanics, 2023, 15(1): 116-127. doi: 10.3934/jgm.2023006
  • In this paper, we shall extend some results regarding the growth estimate of entire solutions of certain type of linear differential equations to that of nonlinear differential equations. Moreover, our results will include several known results for linear differential equations obtained earlier as special cases.


    Spacecraft, underwater vehicles, aerial vehicles, and mobile robots require accurate knowledge of their orientation with respect to a known inertial frame. Typically, attitude estimators rely on the measurements of angular velocity and known inertial vectors in the body-fixed frame. Therefore, the rigid body attitude estimation problem using angular velocity and inertial vectors measurements in a body-fixed frame has been widely studied in past research. In practice, angular velocity is measured at a higher rate than inertial vectors. In this work*, we address the attitude estimation problem given multi-rate measurements of inertially fixed vectors and angular velocity by minimizing the "energy" stored in the state estimation errors.

    *A preliminary arxiv version of this article is available at [1]

    One of the earliest solutions to attitude determination from vector measurements is the TRIAD algorithm, used to determine the rotation matrix from two linearly independent inertial vector measurements [2]. In [3], Wahba presented the attitude determination problem as an optimization problem using three or more vector measurements where the cost function is the sum of the squared norms of vector errors. Various methods have been proposed in the literature to solve the Wahba's problem. Davenport was the first to reduce the Wahba's problem to finding the largest eigenvalue and the corresponding eigenvector of the so-called Davenport's K-matrix [4]. In a similar approach, Mortari presented the EStimator of the Optimal Quaternion (ESOQ) algorithm in[5], which provides the closed-form expressions of a 4×4 matrix's eigenvalues and then computes the eigenvector associated with the greatest of them, representing the optimal quaternion. The QUEST algorithm of [6] determines the attitude that achieves the best-weighted overlap of an arbitrary number of reference vectors. A Singular Value Decomposition (SVD) based method of solving the Wahba's problem was proposed by [7]. Markley also devised a Fast Optimal Matrix Algorithm (FOMA) to solve Wahba's problem in [8]. A coordinate-free framework of geometric mechanics was used in [9] to obtain a solution to Wahba's problem with robustness to measurement noise. [10] provides a numerical solution to Wahba's problem.

    Attitude estimation methods based on minimizing "energy" stored in state estimation errors can be found in [11,12,13,14]. Prior research that has designed attitude estimation schemes based on geometric mechanics includes [15,16,17,18]. Comprehensive surveys of various attitude estimation methods are available in [19,20]. However, most of the aforementioned schemes for attitude estimation work only in continuous time or measurement rich environments and neglect the sparsity in the inertial vector measurements. Inertial vector measurements are usually obtained with the help of Sun (and star) sensors or magnetometers which are accurate but suffer from lower sampling rates. [21] provides one of the first solutions to rigid body attitude estimation with multi-rate measurements using uncertainty ellipsoids. A recursive method based on the cascade combination of an output predictor and an attitude observer can be found in [22]. An attitude estimation scheme on the Special Orthogonal group using intermittent body-frame vector measurements was presented in [23]. In [13], a filtering scheme in continuous-time is proposed by applying the Lagrange-d'Alembert principle on suitably formulated artificial kinetic and potential energy functions where the authors formulate filter equations assuming that inertial vector measurements and angular velocity measurements are available synchronously and continuously. A discrete-time estimation scheme in the presence of multi-rate measurements is proposed in [14]. This work presents the Lyapunov analysis but does not contain a variational interpretation.

    The previous works in attitude estimations are presented in a continuous time domain and then discretized for numerical implementation. This voids the theoretical guarantee of asymptotic stability provided by the continuous-time estimator. Furthermore, most of the previous work in attitude estimation do not consider the realistic scenario when the inertial vectors and angular velocities are measured at a different rate. In the current work, we try to overcome these shortcoming by focusing on developing an asymptotically stable optimal geometric discrete-time attitude estimator based on the minimization of "energy" stored in the errors of the state estimators in the presence of multi-rate measurements. The measurements can be corrupted by noise and we do not assume any specific statistics (like normal distribution) on the measurement noise. However, the noise is assumed to be bounded. We represent the attitude as a rotation matrix which precludes potential singularity issues due to local coordinates (such as Euler angles). The multi-rate discrete-time filtering scheme presented here is obtained by applying the discrete-time Lagrange-D'Alembert principle [24] on a discrete-time lagrangian followed by a discrete-time Lyapunov analysis using a Lyapunov candidate that depends on the state estimation errors. The filtering scheme provided is asymptotically stable with an almost global region of convergence.

    This paper is organized as follows. In the section 2, the attitude estimation problem is formulated as Wahba's optimization problem and then some important properties of the Wahba's cost function are presented. In the section 3, continuous-time rigid body attitude kinematics has been discretized and the propagation model for the measurements in the multi-rate measurement case is presented. Section 4 contains the application of variational mechanics to obtain a filter equation for attitude estimation. The filter equations obtained in the 4 are proven to be asymptotically stable with an almost global domain of convergence by deriving an appropriate dissipation torque using the discrete-time Lyapunov method in the 5. Filter equations are numerically verified with realistic measurements (corrupted by bounded noise) in the 6. Finally, 7 presents the concluding remarks with contributions and future work.

    We define the trace inner product on Rm×n as

    A1,A2:=trace(AT1A2). (2.1)

    The group of orthogonal frame transformations on R3 is defined by O(3):={QR3×3|det(Q)=±1}. The Special orthogonal group on R3 is denoted as SO(3) defined as SO(3):={RR3×3|RTR=RRT=I3}. The corresponding Lie algebra is denoted as so(3):={MR3×3|M+MT=0}. Let ()×:R3so(3)R3×3 be the skew-symmetric matrix cross-product operator and denotes the vector space isomorphism between R3 and so(3):

    v×=[v1v2v3]×=[0v3v2v30v1v2v10]. (2.2)

    Further, let vex():so(3)R3 be the inverse of ()×. exp():so(3)SO(3) is the map defined as

    exp(M)=i=01k!Mk. (2.3)

    We define Ad:SO(3)×so(3)so(3) as

    AdRΩ×=RΩ×RT=(RΩ)×. (2.4)

    In the rest of the article, the text "Consider the time interval [t0,T]R+", indicates that the estimation process will be carried out for the time interval [t0,T] and the time interval is divided into N equal sub-intervals [ti,ti+1] for i=0,1,,N with tN=T with time step size is denoted as, h:=ti+1ti; unless specified otherwise. Furthermore, throughout this article we use ()m to denote measured quantities. For example, if Ω is used to denote angular velocity, Ωm denotes measured value of the angular velocity.

    For the attitude estimation, we consider kN known and linearly independent inertial vectors in the body-fixed frame. Let's denote these vectors in the body-fixed frame by umjR3 for j=1,,k, where k2. Note, that k2 is necessary for determining the attitude uniquely. When k=2, the cross product of the two measured vectors is used as the third independent measurement. Let ejR3 be the corresponding known inertial vectors. We denote the true vectors in the body-fixed frame by uj:=RTej, where R is the rotation matrix of the body-fixed frame with respect to the inertial frame. This rotation matrix provides a coordinate-free global and unique description of the attitude of the rigid body. Define the matrix populated by all k measured vectors expressed in the body-fixed frame as column vectors as

    Um=[um1um2um1×um2]R3×3whenk=2and,Um=[um1um2umk]R3×kwhenk>2, (2.5)

    and the corresponding inertial frame vectors as

    E=[e1e2e1×e2]R3×3whenk=2and,E=[e1e2ek]R3×kwhenk>2. (2.6)

    The true body vector matrix is as below:

    U=RTE=[u1u2u1×u2]R3×3whenk=2and,U=RTE=[u1u2uk]R3×kwhenk>2. (2.7)

    The optimal attitude determination problem using a set of vector measurements is finding an estimated rotation matrix ˆRSO(3), such that the weighted sum of squared norms of the vector errors

    sj=ejˆRumj (2.8)

    is minimized. This attitude determination problem is known as Wahba's problem and consists of minimizing

    U(ˆR,Um)=12kj=1wj(ejˆRumj)T(ejˆRumj) (2.9)

    with the respect to ˆRSO(3), and the weights wj>0 for all j{1,2,,k}. we can express (2.9) as

    U(ˆR,Um)=12EˆRUm,(EˆRUm)W, (2.10)

    where Um is given by (2.5), E is given by (2.6), and W=diag(wi) is the positive definite diagonal matrix of the weight factors for the measured directions. W in (2.10) can be generalized to be any positive definite matrix. It is to be noted that, having different weights for different vector measurements is helpful as in general we can have some sensors that provide more accurate measurements than others. Therefore, it makes sense to provide more weight to those accurate measurements during attitude estimation.

    This is an observer design that, like any other observer, can filter measurement noise from realistic measurements. As observer designs are based on perfect (deterministic) measurements for which their stability properties apply, it is useful to confirm these properties through numerical simulations. We have Um=U=RTE in the absence of measurement errors or noise. Let Q=RˆRTSO(3) denote the attitude estimation error. The following lemmas from [13] stated here without proof give the structure and characterization of critical points of the Wahba's cost function.

    Lemma 2.1. Let rank(E) = 3 and the singular value decomposition of E be given by

    E:=UEΣEVTEwhereUEO(3),VESO(m).ΣEDiag+(3,m), (2.11)

    and Diag+(n1,n2) is the vector space of n1×n2 matrices with positive entries along the main diagonal and all the other components zero. Let σ1,σ2,σ3 denote the main diagonal entries of ΣE. Further, Let W from (2.10) be given by

    W=VEW0VTEwhereW0Diag+(m,m) (2.12)

    and the first three diagonal entries of W0 are given by

    w1=d1σ21,w2=d2σ22,w3=d3σ23whered1,d2,d3>0. (2.13)

    Then, K=EWET is positive definite and

    K=UEΔUTEwhereΔ=diag(d1,d2,d3) (2.14)

    is its eigen-decomposition. Moreover, if didj for ij and i,j{1,2,3} then IQ,K is a Morse function whose set of critical points given as the solution of SK(Q):=vex(KQQTK)=0:

    CQ:={I,Q1,Q2,Q3}whereQi=2UEaiaTiUTEI (2.15)

    and ai is the ith column vector of the identity matrix ISO(3).

    Lemma 2.2. Let K=EWET have the properties given by lemma 2.1. Then the map SO(3)QIQ,KR with critical points given by (2.15) has a global minimum at the identity ISO(3), a global maximum and two hyperbolic saddle points whose indices depend on the distinct eigenvalues d1,d2, and d3 of K.

    Consider the time interval [t0,T]R+. Let the true angular velocity in the body-fixed frame be denoted by ΩR3. The true and measured angular velocities at the time instant ti will be denoted by Ωi and Ωmi respectively. Further, let Ui and Umi denote the matrix formed by true and measured inertial vectors in the body-fixed frame at the time instant ti respectively. The assumption is that angular velocity measurements and inertial vectors measurements in the body-fixed frame are coming at a different but constant rate. In general coarse rate gyros have much higher sampling rate than that of a coarse attitude sensor. Therefore, in a realistic scenario, angular velocities are measured at a higher rate than the inertial vector measurements in the body-fixed frame. Therefore, we assume that the measurements of angular velocity (Ωm) are available after each time interval h say, Ωm0,Ωm1,,ΩmN while, inertial vector measurements in the body-fixed frame are available after time interval nh,nN say, Um0,Umn,Um2n,.

    We have, U=RTE. Therefore, at time instants ti and ti+1, the following relations will hold true respectively; Ui=RTiEi,Ui+1=RTi+1Ei+1. Here, Ri and Ri+1 are the rotation matrices from body-fixed frame to inertial frame at time instants ti and ti+1 respectively. Ei=Ei+1=E are the corresponding known vectors expressed in the inertial frame. Note that the vectors are fixed in the inertial frame and do not change with the time. The continuous time attitude kinematics are,

    ˙R=RΩ×. (3.1)

    We discretize the kinematics in (3.1) as follows,

    Ri+1=Riexp(h2(Ωi+1+Ωi)×), (3.2)

    where, exp():so(3)SO(3) is the map defined as,

    exp(M)=k=01k!Mk. (3.3)

    Using (2.7) and the discretization from (3.2) we get,

    Ui+1=exp(h2(Ωi+1+Ωi)×)RTiEi=exp(h2(Ωi+1+Ωi)×)Ui. (3.4)

    For the instants of time when inertial vector measurements in the body-fixed frame are not available we will use (3.4) to obtain the missing values of Umi. This implies that for the time instants (n1)h<ti<nh,nN, by employing the propagation scheme in (3.4), we propagate direction vector measurements between the instants at which they are measured, using the angular velocity measurements that are obtained at a faster rate. We now formalize the aforementioned inertial vector measurement model as below,

    ˜Umi:={Umi,ifimodn=0exp(h2(Ωmi1+Ωmi)×)˜Umi1,otherwise. (3.5)

    Note that in the absence of measurements errors, we have Ωmi=Ωi,i{0,1,,N}. Also, Umi=Ui for the time instants when inertial vector measurements are available. Now, at time instant t0, we have ˜Um0=Um0=U0 and Ωm0=Ω0. Using (3.5) at time instant t1, noting that Ωm1=Ω1, we get ˜Um1=exp(h2(Ω0+Ω1)×)U0. Comparing it with (3.4), we have ˜Um1=U1. Using the relation from (2.7) we have ˜Um1=RT1E1. Similarly, combining (3.4), and (3.5), and using the relation in (2.7) we get the following relation for all i{0,1,,N} in the absence of measurement errors:

    ˜Umi=RTiEi. (3.6)

    Consider the time interval [t0,T]R+. Let the true angular velocity in the body-fixed frame be denoted by ΩR3. The true and measured angular velocities at the time instant ti will be denoted by Ωi and Ωmi respectively. It is assumed that a measurement Um0 is available at time instant t0. Now the next measurement will be available after a time interval nh where n is an integer variable taking values between n1 and n2 which are user provided. Note that since n is varying, that sampling rate of inertial vector sensor is not fixed but varying over time. Considering the discrete kinematics of R and U from (3.2) and (3.4) we propose the following discrete-time vector measurement model,

    ˜Umi:={Umi,ifmeasurement availableexp(h2(Ωmi1+Ωmi)×)˜Umi1,otherwise. (3.7)

    Note that now we have obtained values of Ωmi and ˜Umi for each time instant ti and similar to the previous arguments, it can be proven that in the absence of measurement errors, we have Ωmi=Ωi and ˜Umi=RTiEi.

    Again consider the time interval [t0,T]R+. Let the angular velocity measurements ΩmR3 be available after each time interval of h. Now suppose the inertial vector measurements are available after a time interval rh where r=n1/n2 is a rational number with n1,n2N. Let n be the greatest common factor of n1 and n2. Let h=nh/n2. Now, consider the time interval [t0,T]R+ divided into N equal sub-intervals [ti,ti+1] for i=0,1,,N with tN=T and let ti+1ti=h. Note that a measurement Ωm will be available after a time interval of n2h/n and a measurement Um will be available after time interval rh=n1n2×n2nh=n1h/n.

    We define the measurement model for Ωm as

    ˜Ωmi:={Ωmi,ifimodn2/n=0Ωmi1,otherwise. (3.8)

    for all i=0,1,,N. Now we can define the vector measurement model to be,

    ˜Umi:={Umi,ifimodn1/n=0exp(h2(Ωmi1+Ωmi)×)˜Umi1,otherwise. (3.9)

    Note that the division of the time interval [t0,T] with the time step size of h is finer than that of with the time step size of h with the relation as h=n2nh. Therefore for each i{0,1,,N} we have j=n2ni{0,1,,n} such that ti=tj. Redefining Ωmi:=˜Ωmn2ni, and ˜Umi=˜Umn2ni for all i{0,1,,N}. Note that now we have obtained values of Ωmi and ˜Umi for each time instant ti and similar to previous arguments, it can be proven that in the absence of measurement errors, we have Ωmi=Ωi and ˜Umi=RTiEi.

    Observe that in any of the aforementioned three cases, the measurement models provide values Ωmi and ˜Umi for all the discrete-time instants ti,i{0,1,,N} in the time interval [t0,T].

    The value of the Wahba's cost function at each instant encapsulates the error in the attitude estimation. We can consider the Wahba's cost function as an artificial potential energy-like term. Therefore using (2.10) we have

    Ui=U(ˆRi,˜Umi)=12EiˆRi˜Umi,(EiˆRi˜Umi)Wi, (4.1)

    where ˜Umi is according to the inertial vector propagation model presented in section 3. The term encapsulating the "energy" in the angular velocity estimation error is denoted by the map Tv defined as

    Tvi:=Tv(ˆΩi,Ωmi,ˆΩi+1,Ωmi+1):=m2(Ωmi+Ωmi+1ˆΩiˆΩi+1)T(Ωmi+Ωmi+1ˆΩiˆΩi+1), (4.2)

    where m>0 is a scalar and ˆΩi stands for the value of estimated angular velocity at the discrete time instant ti. We can write (4.2) in terms of the angular velocity estimation error ωi:=ΩmiˆΩi as

    Tv(ωi,ωi+1)=m2(ωi+ωi+1)T(ωi+ωi+1). (4.3)

    The discrete time Lagrangian can be written as:

    L(ˆRi,˜Umi,ωi,ωi+1)=Tv(ωi,ωi+1)U(ˆRi,˜Umi)=m2(ωi+ωi+1)T(ωi+ωi+1)12EiˆRi˜Umi,(EiˆRi˜Umi)Wi. (4.4)

    We use variational mechanics [25,26] approach to obtain an optimal estimation scheme for the constructed Lagrangian. In the absence of a dissipative term, the variational approach would result in (a generalization to the Lie group of) the Euler-Lagrange equations that we obtain in the context of optimal control. The generalization would be an Euler-Poincare equation on SO(3) [27]. Therefore, if the estimation process is started at time t0, then the discrete-time action functional corresponding to the discrete-time Lagrangian (4.4) over the time interval [t0,T] can be expressed as

    sd(L(ˆRi,˜Umi,ωi,ωi+1))=hNi=0(L(ˆRi,˜Umi,ωi,ωi+1))=hNi=0{m2(ωi+ωi+1)T(ωi+ωi+1)12EiˆRi˜Umi,(EiˆRi˜Umi)Wi}. (4.5)

    Consider attitude state estimation in discrete-time in the presence of multirate measurements with noise and initial state estimate errors. Applying the discrete-time Lagrange-d'Alembert principle [24] from variational mechanics to the action functional sd(L(ˆRi,˜Umi,ωi,ωi+1)) given by (4.5), in the presence of a dissipation term on ωi=ΩmiˆΩi, leads to the following attitude and angular velocity filtering scheme.

    Proposition 1. Consider the time-interval [t0,T]. Consider the multi-rate measurement models from Section 3 such that we have values Ωmi and ˜Umi for all the discrete-time instants ti,i{0,1,,N} in the time interval [t0,T]. Let the Wi be chosen such that Ki=EiWiETi satisfies eigen-decomposition condition (2.14) of lemma 2.1. Also, let τDiR3 denote the value of the dissipation torque at the time instant ti. A discrete-time optimal filter obtained by applying the discrete-time Lagrange-d'Alembert principle would be as follows:

    {ˆRi+1=ˆRiexp(h2(ˆΩi+1+ˆΩi)×)m(ωi+2+ωi+1)=exp(h2(ˆΩi+2+ˆΩi+1)×){m(ωi+1+ωi)+h2SLi+1(ˆRi+1)h2τDi+1}ˆΩi=Ωmiωi, (4.6)

    where SLi(ˆRi)=vex(LTiˆRiˆRTiLi)R3, Li=EiWi(˜Umi)T, and (ˆR0,ˆΩ0)SO(3)×R3×3 are initial estimated states.

    Proof. Consider a first variation in the discrete attitude estimate as

    δˆRi=ˆRiΣ×i, (4.7)

    where ΣiR3 represents a variation for the discrete attitude estimate. For fixed end-point variations, we have Σ0=ΣN=0. A first order approximation is to assume that ˆΩ× and δˆΩ× commute. Taking the first variation of the discrete-time attitude kinematics according to the first equation of (4.6) and comparing with (4.7) we get

    δˆRi+1=δˆRiexp(h2(ˆΩi+1+ˆΩi)×)+h2ˆRiexp(h2(ˆΩi+1+ˆΩi)×)δ(ˆΩi+1+ˆΩi)×=ˆRi+1Σ×i+1. (4.8)

    (4.8) can be rearranged to

    ˆRi+1h2δ(ˆΩi+1+ˆΩi)×=ˆRi+1Σ×i+1ˆRi+1Adexp(h2(ˆΩi+1+ˆΩi)×)Σ×ih2δ(ˆΩi+1+ˆΩi)×=Σ×i+1Adexp(h2(ˆΩi+1+ˆΩi)×)Σ×i. (4.9)

    (4.9) can be equivalently written as an equation in R3 as follows:

    h2δ(ˆΩi+1+ˆΩi)=Σi+1exp(h2(ˆΩi+1+ˆΩi)×)Σi. (4.10)

    Note that, ωi=ΩmiˆΩi gives us

    δ(ωi+1+ωi)=δ(ˆΩi+1+ˆΩi). (4.11)

    Consider artificial potential energy term as expressed in (4.1). Taking its first variation with the respect to estimated attitude ˆR, we get

    δUi=12{δˆRi˜Umi,(EiˆRi˜Umi)Wi+EiˆRi˜Umi,(δˆRi˜Umi)Wi}=δˆRi˜Umi,(EiˆRi˜Umi)Wi=ˆRiΣ×i,(EiˆRi˜Umi)Wi=trace((˜Umi)TΣ×iˆRTi(EiˆRi˜Umi)Wi)=trace((Σ×i)T˜UmiWiETiˆRi)=Σ×i,˜UmiWiETiˆRi=12Σ×,˜UmiWiETiˆRiˆRTiEiWi(˜Umi)T=12Σ×i,LTiˆRiˆRTiLi=STLi(ˆRi)Σi. (4.12)

    Consider the first variation in the artificial kinetic energy Tv(ωi,ωi+1) as in (4.3) with the respect to the angular velocity estimation error,

    δTvi=m(ωi+ωi+1)Tδ(ωi+ωi+1). (4.13)

    Using the results in (4.10) and (4.11) we get

    δTvi=m(ωi+ωi+1)Tδ(Ωi+Ωi+1)=2hm(ωi+ωi+1)T(exp(h2(ˆΩi+1+ˆΩi)×)ΣiΣi+1). (4.14)

    The first variation of the discrete-time action sum in (4.5) using (4.12) and (4.14) can be written as

    δsd=hNi=0{δTviδUi}=Ni=0{2m(ωi+ωi+1)Texp(h2(ˆΩi+1+ˆΩi)×)Σi2m(ωi+ωi+1)TΣi+1hSTLi(ˆRi)Σi}. (4.15)

    Applying the discrete-time Lagrange-d'Alembert principle to the attitude motion, we obtain

    δsd+hN1i=0τTDiΣi=0N1i=0{2m(ωi+ωi+1)T[exp(h2(ˆΩi+1+ˆΩi)×)ΣiΣi+1]hSTLi(ˆRi)Σi+hτTDiΣi}=0. (4.16)

    For 0i<N, (4.16) leads to

    2m(ωi+2+ωi+1)Texp(h2(ˆΩi+2+ˆΩi+1)×)2m(ωi+1+ωi)ThSTLi+1(ˆRi+1)+hτTDi+1=02mexp(h2(ˆΩi+2+ˆΩi+1)×)(ωi+2+ωi+1)=2m(ωi+1+ωi)+hSLi+1(ˆRi+1)hτDi+1, (4.17)

    which in turn leads to the second filter equation.

    Remark 1. Please note that, we introduced a virtual dissipation torque in the Proposition 1 and the performance of the estimator will depend on the value of dissipation torque. In the next section, we determine the value of dissipation torque via discrete-Lyapunov analysis so that the resulting estimator is almost globally asymptotically stable.

    For the Lyapunov stability of the filter equations, we need to construct a suitable Lyapunov candidate function. We use the Wahba's cost function expressed in (4.1) as the artificial potential energy which encapsulates the error in the estimation of attitude. A new term encapsulating the "energy" in the angular velocity estimation error can be constructed as the map Tl:R3×R3R defined as

    Tli:=Tl(ˆΩi,Ωmi):=m2(ΩmiˆΩi)T(ΩmiˆΩi), (5.1)

    where m>0 is a scalar same as before. Further, (5.1) can be written in terms of angular velocity estimation error, ωi:=ΩmiˆΩi as follows:

    Tl(ωi)=m2(ωi)T(ωi). (5.2)

    In the absence of measurement errors, we have ˜Umi=RTiEi. Therefore we can we can write (4.1) in terms of state estimation error Qi=RiˆRTi as

    U(ˆRi,˜Umi)=12EiˆRiRTiEi,(EiˆRiRTiEi)Wi=IRiˆRTi,EiWiETiUi=U(Qi)=IQi,KiwhereKi=EiWiETi. (5.3)

    The weights Wi's are chosen such that Ki is always positive definite with distinct eigenvalues according to lemma 2.1.

    Theorem 5.1. Consider the time-interval [t0,T]. Consider the multi-rate measurement models from Section 3 such that we have values Ωmi and ˜Umi for all the discrete-time instants ti,i{0,1,,N} in the time interval [t0,T]. Then the estimation scheme in Proposition 1 with the following value of the dissipation torque:

    τDi+1=1h{2m(ωi+1+ωi)+hSLi+1(ˆRi+1)2mm+lexp(h2(ˆΩi+2+ˆΩi+1)×)[2mωi+1+kphSLi+1(ˆRi+1)]}, (5.4)

    leads to the estimation scheme

    {ωi+1=1m+l[(ml)ωi+kphSLi(ˆRi)]ˆΩi=ΩmiωiˆRi+1=ˆRiexp(h2(ˆΩi+1+ˆΩi)×), (5.5)

    where SLi(ˆRi)=vex(LTiˆRiˆRTiLi)R3, Li=EiWi(˜Umi)T, l>0, lm and kp>0, which is asymptotically stable at the estimation error state (Q,ω):=(I,0) (Qi=RiˆRTi) in the absence of measurement errors. Further, the domain of attraction of (I,0) is a dense open subset of SO(3)×R3.

    Proof. Using the third equation from (5.5) we have

    Qi+1=Ri+1ˆRTi+1=QiˆRiexp(h2(ωi+1+ωi)×)ˆRTi. (5.6)

    We choose the following discrete-time Lyapunov candidate:

    Vi:=V(Qi,ωi):=kpUi+Tli, (5.7)

    where kp>0 is a constant.

    The stability of the attitude and angular velocity estimation error can be shown by analyzing ΔVi=kpΔUi+ΔTli.

    Assuming Ki to be constant and letting K=Ki=Ki+1 we obtain

    ΔUi=Ui+1Ui=IQi+1,KIQi,KΔUi=QiQi+1,K=ΔQi,K, (5.8)

    where ΔQi=Qi+1Qi. Therefore,

    ΔQi=Qi+1Qi=Qi[ˆRiexp(h2(ωi+1+ωi)×)ˆRTiI]. (5.9)

    Considering the first order expansion of exp(h2(ˆωi+1+ˆωi)×) as

    exp(h2(ωi+1+ωi)×)I+h2(ωi+1+ωi)×, (5.10)

    we have

    ΔQi=Qi[ˆRi(I+h2(ωi+1+ωi)×)ˆRTiI]=h2Qi(ˆRi(ωi+1+ωi)׈RTi)=h2Qi(ˆRi(ωi+1+ωi))×. (5.11)

    It has to be noted that approximation in (5.10) is accurate for small values of h and may affect the stability results for very high values of h. In the absence of measurement errors, we have ˜Umi=RTiEi. Therefore, it follows that

    ΔUi=h2Qi(ˆRi(ωi+1+ωi))×,K=h2Ri(ωi+1+ωi)׈RTi,EiWiETi=h2(ωi+1+ωi)׈RTi,RTiEiWiETi=h2(ωi+1+ωi)׈RTi,˜UmiWiETi, (5.12)

    and noting that Li=EiWi(˜Umi)T, we get

    ΔUi=h2(ωi+1+ωi)×,LTiˆRi=h4(ωi+1+ωi)×,LTiˆRiˆRTiLi=h2(ωi+1+ωi)TSLi(ˆRi). (5.13)

    Similarly we can compute the change in the kinetic energy as follows:

    ΔTli=Tl(ωi+1)Tl(ωi)=(ωi+1+ωi)Tm2(ωi+1ωi)ΔTli=(ωi+1+ωi)Tm2(ωi+1ωi). (5.14)

    The change in the value of the candidate Lyapunov function can be computed as,

    ΔVi=Vi+1Vi=ΔTi+kpΔUi=12(ωi+1+ωi)T(m(ωi+1ωi)kphSLi(ˆRi)). (5.15)

    Similarly, we obtain

    ΔVi+1=12(ωi+2+ωi+1)T(m(ωi+2ωi+1)kphSLi+1(ˆRi+1)). (5.16)

    Substituting the value of ωi+2 from the filtering scheme presented in Proposition 1 we get

    ΔVi+1=12(ωi+2+ωi+1)T{exp(h2(ˆΩi+2+ˆΩi+1)×){m(ωi+1+ωi)+hSLi+1(ˆRi+1)hτDi+1}2mωi+1kphSLi+1(ˆRi+1)}. (5.17)

    Now, for ΔV to be negative definite for all i we require

    exp(h2(ˆΩi+2+ˆΩi+1)×){m(ωi+1+ωi)h2τDi+1+h2SLi+1(ˆRi+1)}2mωi+1kphSLi+1(ˆRi+1)=l(ωi+2+ωi+1), (5.18)

    where l>0,lm, and ΔVi+1 simplifies to

    ΔVi+1=l2(ωi+2+ωi+1)T(ωi+2+ωi+1). (5.19)

    Substituting ωi+2 from the third equation presented in Proposition 1 into (5.18),

    exp(h2(ˆΩi+2+ˆΩi+1)×){m(ωi+1+ωi)h2τDi+1+h2SLi+1(ˆRi+1)}2mωi+1kphSLi+1(ˆRi+1)=lmexp(h2(ˆΩi+1+ˆΩi)×){m(ωi+1+ωi)h2τDi+1+h2SLi+1(ˆRi+1)}, (5.20)

    which further simplifies to,

    m+lmexp(h2(ˆΩi+2+ˆΩi+1)×){m(ωi+1+ωi)h2τDi+1+h2SLi+1(ˆRi+1)}=2mωi+1+kphSLi+1(ˆRi+1), (5.21)

    which upon simple manipulations yields (5.4). It can be seen that after substituting (5.4) into Proposition 1, we obtain

    ωi+2=1m+l[(ml)ωi+1+kphSLi+1(ˆRi+1)]. (5.22)

    (5.22) can also be rewritten as

    ωi+1=1m+l[(ml)ωi+kphSLi(ˆRi)], (5.23)

    in terms of ωi,ωi+1 and SLi(ˆRi). From (5.19), ΔVi can be written as

    ΔVi=l2(ωi+1+ωi)T(ωi+1+ωi). (5.24)

    We employ the discrete-time La-Salle invariance principle from [28] considering our domain (SO(3)×R3) to be a subset of R12. We use Theorem 6.3 and Theorem 7.9 from Chapter-1 of [28]. For this we first compute E:={(Qi,ωi)SO(3)×R3|ΔVi(Qi,ωi)=0}={(Qi,ωi)SO(3)×R3|ωi+1+ωi=0}. From (5.6), ωi+1+ωi=0 implies that

    Qi+1=Qi. (5.25)

    Also, from (5.13) we have ΔU=0 whenever ωi+1+ωi=0. This implies that the potential function, which is a Morse function according to lemma 2.1, is not changing and therefore has converged to one of its stationary points. Stationary points of the Morse function IQ,K are characterized as the solutions of

    SK(Qi)=0vex(KQiQTiK)=0KQi=QTiK. (5.26)

    Now, Li=EiWi(˜Umi)T=EiWi(RTiEi)T=(EiWiETi)Ri=KRi, which further gives us

    (SLi(ˆRi))×=LTiˆRiˆRTiLi=RTiKˆRiˆRTiKRi. (5.27)

    Multiplying (5.27) from the right hand side by ˆRTi and from the left hand side by ˆRi,

    ˆRi(SLi(ˆRi))׈RTi=ˆRiRTiKKRiˆRTi=QTiKKQi. (5.28)

    At the critical points from (5.26), the right side of the above expression vanishes. Therefore, as ˆRi is an orthogonal matrix, the following holds true at the critical points:

    (SLi(ˆRi))×=0SLi(ˆRi)=0. (5.29)

    Substituting this information in (5.23) yields,

    ωi+1=1m+l(ml)(ωi). (5.30)

    Now if, ωi+1+ωi=0, we have

    2mm+lωi=0ωi=ωi+1=0. (5.31)

    This leads to the conclusion that the set of estimation errors, E={(Qi,ωi)SO(3)×R3|QiCQ,ωi=0}, is the largest invariant set for the estimation error dynamics, and we obtain M=E={(Qi,ωi)SO(3)×R3|QiCQ,ωi=0}. Therefore, we obtain the positive limit set as the set,

    I:=MV1i(0)={(Q,ω)SO(3)×R3|QCQ,ω=0}. (5.32)

    In the absence of measurement errors, all the solutions of this filter converge asymptotically to the set I. More specifically, the attitude estimation error converges to the set of critical points of IQ,K. The unique global minimum of this function is at (Q,ω)=(I,0) from lemma 2.2. Therefore, (Q,ω)=(I,0) is locally asymptotically stable. The remainder of this proof is similar to the last part of the proof of stability of the variational attitude estimator in [13].

    Consider the set,

    C=I(I,0) (5.33)

    which consists of all the stationary states that the estimation errors may converge to, besides the desired estimation error state (I,0). Note that all states in the stable manifold of a stationary state in C will converge to this stationary state. From the properties of the critical points QiCQ(I) of Φ(K,IQ) given in lemma 2.2. we see that the stationary points in I(I,0)={(Qi,0):QiCQ(I)} have stable manifolds whose dimensions depend on the index of Qi. Since the angular velocity estimate error ω converges globally to the zero vector, the dimension of the stable manifold MSi of (Qi,0)SO(3)×R3 is

    dim(MSi)=3+(3index of Qi)=6index of Qi. (5.34)

    therefore, the stable manifolds of (Q,ω)=(Qi,0) are three-dimensional, four dimensional, or five-dimensional, depending on the index of QiCQ(I) according to (5.34). Moreover, the value of the Lyapunov function V(Qi,ωi) is non decreasing (increasing when (Qi,ωi)I) for trajectories on these manifolds when going backwards in time. This implies that the metric distance between error states (Q,ω) along these trajectories on the stable manifolds MSi grows with the time separation between these states, and this property does not depend on the choice of the metric on SO(3)×R3. Therefore, these stable manifolds are embedded (closed) sub-manifolds of SO(3)×R3 and so is their union. Clearly, all states starting in the complement of this union, converge to the stable equilibrium (Q,ω)=(I,0); therefore the domain of attraction of this equilibrium is,

    DOA(I,0)=SO(3)×R3{3i=1MSi} (5.35)

    which is a dense open subset of SO(3)×R3.

    This section presents numerical simulation results of the discrete-time estimator presented in the section 5. The rigid body is assumed to have an initial attitude and angular velocity given by,

    R0=expmSO(3)((π4×[47,27,57]T)×),
    andΩ0=π60×[1.2,2.1,1.9]Trad/s.

    W is selected based on the measured set of inertial vectors E at each instant such that it satisfies lemma 2.1. Initially estimated states have the following initial estimation errors:

    Q0=expmSO(3)((π2.5×[47,27,57]T)×),
    andω0=[0.001,0.002,0.003]Trad/s.

    It has been assumed that there are at most 9 inertially known directions that are being measured by the sensors attached to the rigid body. The number of observed direction can vary randomly between 2 to 9 at each time instant. In the case where the number of observed directions is 2, the cross product of the two measurements is used as the third measurement. The standard rigid body dynamics are used to produce true states of the rigid body by applying sinusoidal forces. The observed directions in the body-fixed frame are simulated with the help of the aforementioned true states. The true quantities are disturbed by bounded, random noise with zero mean to simulate realistic measurements. Based on coarse attitude sensors like sun sensors and magnetometers, a random noise bounded in magnitude by 2.4 is added to the matrix U=RTE to generate measured Um. Similarly, a random noise bounded in magnitude by 0.97/s, which is close to real noise levels of coarse rate gyros, is added Ω to generate measured Ωm.

    The estimator is simulated over a time interval of T = 60s, with a step-size of h=0.01s. The inertial scalar gain is m=1.5 and the dissipation term is chosen to be l=0.3. The difference of sampling rate between measurements of angular velocity and measurements inertial vectors in body-fixed frame is taken to be n=10. Furthermore, the value of gain kp is chosen to be kp=1. The principle angle ϕ of the rigid body's attitude estimation error Q is shown in the Figure 1. Components of estimation error ω in the rigid body's angular velocity are shown in Figure 2.

    Figure 1.  Case-1: Principle angle of the attitude estimation error.
    Figure 2.  Case-1: Angular velocity estimation error.

    The estimator is simulated over a time interval of T = 60s, with a step-size of h=0.01s. The inertial scalar gain is m=1.5 and the dissipation term is chosen to be l=0.3. The difference of sampling rate between measurements of angular velocity and measurements inertial vectors in body-fixed frame is to varying randomly but bounded between 10 and 30. Furthermore, the value of gain kp is chosen to be kp=1. The principle angle ϕ of the rigid body's attitude estimation error Q is shown in the Figure 3. Components of estimation error ω in the rigid body's angular velocity are shown in Figure 4.

    Figure 3.  Case-2: Principle angle of the attitude estimation error.
    Figure 4.  Case-2: Angular velocity estimation error.

    The estimator is simulated over a time interval of T = 60s, with a step-size of h=0.008s. The inertial scalar gain is m=2.5 and the dissipation term is chosen to be l=0.5. The assumption is that the angular velocity measurements are available after each time interval of 0.008s and inertial vector measurements are available after each time interval of 0.05s. Furthermore, the value of gain kp is chosen to be kp=10. The principle angle ϕ of the rigid body's attitude estimation error Q is shown in the Figure 5. Components of estimation error ω in the rigid body's angular velocity are shown in Figure 6.

    Figure 5.  Case-3: Principle angle of the attitude estimation error.
    Figure 6.  Case-3: Angular velocity estimation error.

    We have discussed the performance of the estimator in various scenarios but with the fixed values of initial estimation errors. It remains to see whether the estimator perform good in case of large initial estimation errors. For this part, we restrict our attention to the Case-2 i.e., when the difference of sampling rate between measurement of angular velocity and measurement of inertial vectors is time varying and bounded between 10 and 30. We check the performance of the estimator for three different values of initial estimation errors in angular velocity and attitude. The results of the same are plotted in Figure 7.

    Figure 7.  Performance of the estimator for different values of initial estimation errors.

    On comparing Figures 7a and 7b, when the initial estimation error in principle angle of attitude estimation is increased to about 2.5 times and estimation error in angular velocity is increased to 10 times, convergence takes about 10s more time. However, on comparing Figure 7a and 7c, it can be see that if the initial attitude estimation error is fixed then change in initial angular velocity estimation error(100 times the original values) has no effect on the convergence time. Therefore, it can be concluded that increasing the initial attitude estimation error can increase the convergence time but increase in initial angular velocity estimation error has no effect on the convergence time.

    In this section, we discuss how convergence time can be decreased by increasing the value of gain. Same as before, we restrict our attention to the Case-2 i.e., when the difference of sampling rate between measurement of angular velocity and measurement of inertial vectors is time varying and bounded between 10 and 30. We consider very high errors in initial estimation errors and plot the performance of the estimator for three different values of gain kp in Figure 8. As it can be seen, increasing the value of gain significantly improves the convergence time.

    Figure 8.  Performance of the estimator for different values of gain kp.

    We develop a geometric attitude and angular velocity estimation scheme using the discrete-time Lagrange-D'Alembert principle followed by discrete-time Lyapunov stability analysis in the presence of multi-rate measurements. The attitude determination problem from two or more vector measurements in the body-fixed frame is formulated as Wahba's optimization problem. To overcome the multi-rate challenge, a discrete-time model for attitude kinematics is used to propagate the inertial vector measurements forward in time. The filtering scheme is obtained with the aid of appropriate discrete-time Lagrangian and Lyapunov functions consisting of Wahba's cost function as an artificial potential term and a kinetic energy-like term that is quadratic in the angular velocity estimation error. As it can be observed, the Lyapunov function is not constructed from the same artificial potential and kinetic energy terms that are used for constructing the Lagrangian. There are mainly two reasons behind this; 1) the filtering scheme obtained by applying the discrete Lagrange-d'Alembert principle is implicit in nature and therefore it can increase computational load and runtime making it difficult to use for real-time applications. Therefore, an explicit filtering scheme is more desirable, 2) we also need the filtering scheme to be asymptotically stable. A Lyapunov function, different from system energy constructed appropriately, helps us meet both the requirements. The explicit filtering scheme obtained after the Lyapunov analysis was proven to be asymptotically stable in the absence of measurement noise and the domain of convergence is proven to be almost global. Numerical simulations were carried out with realistic measurement data corrupted by bounded noise. Numerical simulations verified that the estimated states converge to a bounded neighborhood of (I,0). Furthermore, the rate of convergence of the estimated states to the real state can be controlled by choosing appropriate gains. Future endeavors are towards obtaining a discrete-time optimal attitude estimator in the presence of multi-rate measurements when there a constant or slowly time-varying bias in the measurements of angular velocity while also obtaining a bound on the state estimation errors when there is measurement noise in the inertial vector measurements and the angular velocity measurements.

    All authors hereby declare that there are no conflicts of interest in this paper.



    [1] S. Bank, Some results on analytic and meromorphic solutions of algebraic differential equations, Adv. Math., 15 (1975), 41-62. doi: 10.1016/0001-8708(75)90124-3
    [2] S. Bank, J. Langley, Oscillation theory for higher order linear differential equations with entire coeffients, Complex Var. Elliptic Eq., 16 (1991), 163-175.
    [3] Z. X. Chen, The growth of solutions of f+ezf+Q(z)f=0 where the order (Q)=1, Sci. China Ser. A, 45 (2002), 290-300.
    [4] Z. X. Chen, K. H. Shon, The hyper order of solution of second order differential equations and subnomal solutions of periodic equations, Taiwanese J. Math., 14 (2010), 611-628. doi: 10.11650/twjm/1500405809
    [5] Z. X. Chen, C. C. Yang, Some further results on the zeros and growths of entire solutions of second order linear differential equations, Kodai Math. J., 22 (1999), 273-285. doi: 10.2996/kmj/1138044047
    [6] A. E. Eremenko, Meromorphic solutions of algebraic differential equations, Russ. Math. Surv., 37 (1982), 61-95.
    [7] M. Frei, Über die subnormalen Löungen der Differentialgleichung ω+ezω+(konst.)ω=0, Comment. Math. Helv. 36 (1962), 1-8.
    [8] F. Gackstatter, I. Laine, Zur Theorie der gewAhnlichen Differentialgleichungen im Komplexen, Ann. Polo. Math., 38 (1980), 259-287. doi: 10.4064/ap-38-3-259-287
    [9] A. A. Gol'dberg, On the single valued integrals of differential equations of the first order, Ukrainian Math. J., 8 (1956), 254-261 (Russian).
    [10] G. Gundersen, On the question of whether f+ezf+B(z)f=0 can admit a solution f0 of finite order, Proc. Roy. Soc. Edinburgh Sect. A, 102 (1986), 9-17. doi: 10.1017/S0308210500014451
    [11] W. K. Hayman, Meromorphic functions, Clarendon Press, Oxford, 1964.
    [12] S. Hellerstein, J. Miles, J. Rossi, On the growth of solutions of certain linear differential equations, Ann. Acad. Sci. Fenn., Ser. A 1 Math., 17 (1992), 343-365. doi: 10.5186/aasfm.1992.1723
    [13] I. Laine, Nevanlinna theory and complex differential equations, Walter de Gruyter, Berlin/New York, 1993.
    [14] P. Li, Entire solutions of certain type of differential equations II, J. Math. Anal. Appl., 375 (2011), 310-319. doi: 10.1016/j.jmaa.2010.09.026
    [15] P. Li, C. C. Yang, On the nonexistence of entire solutions of certain type of nonlinear differential equations, J. Math. Anal. Appl., 320 (2006), 827-835. doi: 10.1016/j.jmaa.2005.07.066
    [16] L. W. Liao, C. C. Yang, J. J. Zhang, On meromorphic solutions of certain type of non-linear differential equations, Ann. Acad. Sci. Fenn. Math., 38 (2013), 581-593. doi: 10.5186/aasfm.2013.3840
    [17] M. Ozawa, On a solution of ω+ezω+(az+b)ω=0, Kodai Math. J., 3 (1980), 295-309.
    [18] N. Steinmetz, Meromorphic solutions of second-order algebraic differential equations, Complex Var. Elliptic Eq., 13 (1989), 75-83.
    [19] N. Toda, On algebroid solutions of algebraic differential equations in the complex plane II, J. Math. Soc. Japan, 45 (1993), 705-717. doi: 10.2969/jmsj/04540705
    [20] J. Wang, W. R. Lü, The fixed points and hyper-order of solutions of second order linear differential equations with meromorphic coefficients, Acta Math. Appl. Sinica, 27 (2004), 72-80 (Chinese).
    [21] E. T. Whittaker, G. N. Watson, A course of modern analysis, Cambridge University Press, Cambridge, 1927.
    [22] H. Wittich, Neuere Untersuchungen über eindeutige analytische Funktionen, Springer, Berlin Heidelberg, 1968.
    [23] C. C. Yang, On deficiencies of differential polynomials, Math. Z., 116 (1970), 197-204. doi: 10.1007/BF01110073
    [24] C. C. Yang, H. X. Yi, Uniqueness theory of meromorphic functions, Science Press, Beijing/New York, 2003.
    [25] K. Yosida, A generalization of Malmquist's theorem, Japan. J. Math., 9 (1933), 253-256.
    [26] J. Zhang, L. W. Liao, On entire solutions of a certain type of nonlinear differential and difference equations, Taiwanese J. Math., 15 (2011), 2145-2157. doi: 10.11650/twjm/1500406427
    [27] J. J. Zhang, X. Q. Lu, L. W. Liao, On exact transcendental meromorphic solutions of nonlinear complex differential equations, Houston J. Math., 45 (2019), 439-453.
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4132) PDF downloads(246) Cited by(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog