# Payload Parameter Identification of a Flexible Space Manipulator System via Complex Eigenvalue Estimation

Authors: Zhiyu Ni; Shunan Wu; Yewei Zhang; Zhigang Wu
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

Manipulator systems are widely used in payload capture and movement in the ground/space operation due to their dexterous manipulation capability. In this study, a method for identifying the payload parameters of a flexible space manipulator using the estimated system of complex eigenvalue matrix is proposed. The original nonlinear dynamic model of the manipulator is linearized at a selected working point. Subsequently, the system state-space model and corresponding complex eigenvalue parameters are determined by the observer/Kalman filter identification algorithm using the torque input signal of the motor and the vibration output signals of the link. Therefore, the inertia parameters of the payload, that is, the mass and the moment of inertia, can be derived from the identified complex eigenvalue system and mode shapes by solving a least-squares problem. In numerical simulations, the proposed parameter identification method is implemented and compared with the classical recursive least-squares and affine projection sign algorithms. Numerical results demonstrate that the proposed method can effectively estimate the payload parameters with satisfactory accuracy.


---

## Body

## 1. Introduction

Manipulator systems are an important component of the spacecraft structure. These systems are currently widely applied to the on-orbit servicing missions, such as space object capture and structure assembly [1–3]. Manipulator systems in space are generally designed as long span and lightweight to increase workspace and reduce launch mass. Thus, the flexibility of links is evident. When the manipulator carries a large endpoint payload or captures a space object, the vibration of the system is obvious, even with only a slow movement. For example, the operators of the shuttle remote manipulator system (SRMS) spend 20 to 40 s waiting for oscillation decay after maneuvering the arm to avoid the vibration problem in an on-orbit operation [4]. In the operation of the flexible space manipulator system, guaranteeing an accurate positioning in pick-and-place tasks is a crucial aspect that must be addressed. One of the major obstacles to overcome is designing a control method capable of canceling the vibrations when the system dynamic model is affected by the changes in the payload parameter. If these changes are disregarded in the controller design, then the algorithm may lose accuracy and effectiveness in the vibration suppression, and the system may become unstable in some cases.Therefore, the payload parameters are identified to obtain the new system characteristic after the operation. The identified parameters can be used to adjust the adaptive control law and provide an important reference for improving the overall performance of the system [5–7]. In particular, if the manipulator is used to capture a payload with an unknown property (e.g., a noncooperative object), then accurately identifying object parameters can be implemented to recalculate uncertain parameters periodically and update the certainty equivalence controller with unexpected parameter variations conveniently [8]. Payload object capture and corresponding parameter identification experiments using the space manipulator have been successfully implemented in certain on-orbit spacecraft, such as the Engineering Test Satellite-VII (ETS-VII) [9], Japanese Experiment Module Remote Manipulator System (JEMRMS) [10], and Orbital Express [11]. However, these experiments primarily focus on cooperative target capture. Thus, the payload parameters are identified, and all manipulators involved are considered rigid structures. Therefore, certain problems regarding structural flexibility and unknown payload parameters in the manipulator system require further study.Various identification methods have been proposed to estimate the inertia parameters (i.e., mass, mass center, and moment of inertia) of the manipulator and the unknown payload. Most of these methods are least-squares (LS) techniques, which are validated by JEMRMS experiments using space manipulators [9, 12–14]. Based on the planning for optimal maneuvers, the inertia parameters of the ETS-VII spacecraft and the manipulator payload are obtained by solving a nonlinear LS problem [9]. However, alternative methods, such as the linear and angular momentum conservation [15, 16], Kalman filtering [17], and the algebraic identification approach [18], are available. Nguyen-Huynh and Sharf developed a new inertia parameter identification method based on the momentum conservation equation and recursive least-squares (RLS) estimation after the space manipulator grasped an unknown tumbling target [15]. Bruggemann et al. also proposed a robotic-based identification method of unknown spacecraft inertia properties using the momentum conservation [19]. Liu et al. introduced a recursive differential evolution algorithm to identify the inertia parameters of an unknown target and revise the friction parameters of the space manipulator joints simultaneously [20]. San-Millan et al. developed a technique for real-time identification of the payload mass parameter for a single-link flexible manipulator based on the algebraic identification algorithm [21]. However, several of these studies only considered the manipulator to be a completely rigid structure and neglected the influence of flexible vibration [9, 15, 20]. Furthermore, many inertia parameter identification approaches may be unsuitable for the space structure because the available system excitation and measurement signals are finite in the space environment. Therefore, developing a general on-orbit identification algorithm is necessary.Notably, the system model and its corresponding eigenvalue parameters may be changed due to variations of the manipulator configuration and payload mass to determine the payload parameters of the flexible space manipulator. Thus, the corresponding system identification issue becomes sophisticated. In this study, only the situation for a flexible single-link space manipulator in a selected working point is considered to simplify the analysis process. Then, the original nonlinear dynamic manipulator model is degraded as a linear model by local linearization technology. Consequently, if the system eigenvalue parameters can be obtained using certain identification algorithms, then the corresponding relationship between the payload and the system eigenvalue parameters can be established from the linearized dynamic equation, and the payload parameter can be determined from the identified structural model parameters.In the existing time-domain identification approaches for the model parameters, the eigensystem realization algorithm (ERA) by Juang and Pappa [22] has been successfully applied for parameter identification experiments and operational modal analysis repeatedly [23, 24]. The ERA is a typical system realization method. By constructing the system Hankel matrix in the discrete-time domain, the singular value decomposition (SVD) is used to determine the system order and separate the noise signals. Then, the modal parameters are obtained from the identified state-space model using system impulse responses. However, the slow decaying response for lightly damped systems may produce a large Hankel matrix and long computation time. Therefore, a type of extended ERA called the observer/Kalman filter identification (OKID) method was developed to solve the aforementioned problem using system input and output (I–O) signals simultaneously [25]. The OKID method is conducted by adding an observer into the system and placing the desired eigenvalues. Thus, the observer Markov parameters are deadbeat. The ERA–OKID series algorithms were proven to be reliable identification algorithms with a certain degree of noise immunity due to SVD technology. These methods have already been used for the identification of the state-space model and modal parameters of the ground structure and spacecraft [22, 26–28]. However, studies involving the identification of manipulator parameters using the ERA series algorithm are limited.The current study mainly is aimed at investigating the payload parameter identification of a flexible manipulator system. A novel method is proposed to determine the unknown payload inertia parameter using the complex eigenvalue estimation of the system. Different from the frequently used inertia parameter identification methods that regard the manipulator as a rigid body [9, 15, 20], the manipulator vibration due to link flexibility is studied, and the nonlinear dynamic equation is established. Subsequently, the nonlinear dynamic equation is linearized at a selected working point, and then the corresponding state-space model is obtained. Moreover, using the designed I–O signal system, the complex eigenvalue parameters of the system are determined by the OKID algorithm. Thus, the payload parameters are derived using the LS method based on the identified eigenvalues and system complex modal analysis. Through comparison with the classical RLS method and the affine projection sign algorithm (APSA), the numerical simulation results demonstrate that the proposed method can be used for payload parameter identification of flexible space manipulator with satisfactory identification accuracy. In practical applications, the proposed method can provide a reference for the capturing process of the unknown object, and the manipulator structural modal and the endpoint payload parameters at the working point can be identified simultaneously.The remainder contents of this paper are organized as follows. The dynamic model of the flexible single-link manipulator is established and further be linearized at the selected working point in Section2. In Section 3, the OKID algorithm is briefly reviewed, and the identification of the payload parameter is studied based on complex eigenvalue estimation. Section 4 is dedicated to the analysis of the numerical simulations, and the identification results of payload parameters using the proposed method are compared with those of the classical RLS and APSA algorithm. Some conclusions are presented in Section 5.


## 2. Model Description of Space Manipulator

In this section, the nonlinear dynamic equation of a flexible single-link manipulator system will be established using the Lagrange method, and then the linearized dynamic model is obtained by the local linearization technology.


### 2.1. Dynamic Modeling Based on Lagrange Method

Before establishing the dynamic model of the flexible single-link manipulator, some assumptions are introduced as follows:(1)
The manipulator moves in a plane, and the gravity influence is ignored
(2)
The flexible link complies with the Euler-Bernoulli beam theory, and the axial deformations as well as nonlinear geometric effect due to bending are negligible
(3)
The link’s cross-sectional area remains constant along the link and the material is homogeneous, so the linear density and Young’s modulus of the link are constant
(4)
In the paper, the payload is simplified as a mass point. Therefore, only consider the identification problem of the mass and moment of inertia parameters of the payload
(5)
The flexible of the joint and the system damping are negligible in this paperBased on these assumptions, the manipulator coordinate system and corresponding parameters are defined as follows: the originO of the inertial coordinate frame O‐XYZ and the origin O1 of the link coordinate frames O1‐x1y1z1 are both selected at the joint. The rotation of the link coordinate frame with respect to the inertial coordinate frame O‐XYZ is defined as θ. In addition, the link’s length, mass per unit length, Young modules, and second moment of area are represented by L, ρ, E, and I, respectively. The input torque of joint is τ and the lateral displacement of link is defined as wx,t. The mass and the moment of inertia of the endpoint payload are represented by me and Je, respectively. A simplified model of the flexible single-link manipulator is shown in Figure 1.Figure 1
Single-link manipulator with endpoint payload.If the manipulator is considered as a cantilever beam structure and the first two modal shapes of the link are selected, then by the assumed modes method, the lateral displacementwx,t can be expressed as follows:
(1)wx,t=ϕ1xη1t+ϕ2xη2t,where the two mode shape functions are selected as follows:
(2)ϕix=chβix−cosβix−chβiL+cosβiLshβiL+sinβiLshβix−sinβix,i=1,2.For any point on the link, the position vectorr and corresponding velocity vectors ṙ in inertial coordinate system can be represented, respectively, by:
(3)r=xcosθ−wsinθxsinθ+wcosθ,ṙ=−θ̇xsinθ−ẇsinθ−θ̇wcosθθ̇xcosθ+ẇcosθ−θ̇wsinθ.The system generalized coordinateq is selected as q=θη1η2T, and the kinetic energy TL of the link can be expressed as follows:
(4)TL=12ρ∫0LṙTṙdx=12q̇TρL33+ρLη12+ρLη220.5688ρL20.0908ρL20.5688ρL2ρL00.0908ρL20ρLq̇=12q̇TM1q̇,and the kinetic energy Te of endpoint payload can be denoted as follows:
(5)Te=12meṙeTṙe+12Jeθ̇+w′L,t2=12q̇TmeL2+me2η1−2η22+Je2meL+2.7528Je−4meL−9.5614Je2meL+2.7528Je4me+7.5779Je−4me−26.3206Je−2meL−9.5614Je−4me−26.3206Je4me+91.4204Jeq̇=12q̇TM2q̇,where ṙe is the end position vector of the link. Therefore, the system total kinetic energy T can be expressed as follows:
(6)T=TL+Te=12q̇TM1+M2q̇=12q̇TMq̇,where M is system generalized mass matrix.The system total potential energyU only includes the elastic potential energy for the link and thus it can be written as follows:
(7)U=12EI∫0L∂2w∂x22dx=12qTEI∫0L0000ϕ1″2ϕ1″ϕ2″0ϕ1″ϕ2″ϕ1″2dxq=12qTKq,where K is system generalized stiffness matrix. Consequently, based on the Lagrange method, the dynamic equation of the manipulator can be denoted as follows:
(8)Mq¨+Ṁq̇−∂T∂q+Kq=u,where Ṁq̇ and ∂T/∂q are the Coriolis and Centripetal forces, respectively. u is the generalized force and u=τ00T.Substituting the generalized mass matrix and stiffness matrix into Equation (8), the system dynamic equation can be rewritten as follows:
(9)Mqq¨+Eq,q̇q̇+Kq=u,and the detailed elements in n×n matrices Mq, Eq,q̇, and K are given in Appendix A.


### 2.2. Linearization of Nonlinear Dynamic Model

This paper mainly focuses on the structural vibration and corresponding payload parameter identification problem, so the global motion of the manipulator in the inertial coordinate frameO‐XYZ will be ignored. In this case, based on the local linearization theory, the nonlinear dynamics Equation (9) can be linearized at the selected working point as follows [29]:
(10)Mq0δq¨t+CEδq̇t+K+KM+KEδqt=Lut,where the notation “⋅0” denotes the value at the working point. Therefore, the linearization state vector δqt=qt−q0. L is the input influence matrix and the input signal ut=τt. The linearized matrices CE, KM, and KE in the equation can be computed, respectively, by:
(11)CE=∂E∂q̇q0,q̇0=∂E∂θ̇,∂E∂η̇1,∂E∂η̇2q0,q̇0,KM=∂M∂qq0q¨0=∂M∂θq¨0,∂M∂η1q¨0,∂M∂η2q¨0q0,KE=∂E∂qq0,q̇0=∂E∂θ,∂E∂η1,∂E∂η2q0,q̇0.To simplify the equation derivation, here it is assumed that the modal displacement, velocity, and acceleration of the link at the working point are very small, that is,ηi0=η̇i0=η¨i0≈0 and i=1,2, so the terms CE, KM, and KE in Equation (10) will be zero. Then, the linearized dynamics Equation (10) can be further simplified and the corresponding measurement equation can be expressed, respectively, as follows:
(12)Mq0δq¨t+Kδqt=Lut,(13)yt=Caδq¨t+Cvδq̇t+Cdδqt,where
(14)Mq0=ρL33+meL2+Je0.5688ρL2+2meL+2.7528Je0.0908ρL2−2meL−9.5614Je0.5688ρL2+2meL+2.7528JeρL+4me+7.5779Je−4me−26.3206Je0.0908ρL2−2meL−9.5614Je−4me−26.3206JeρL+4me+91.4204Je,and yt is m×1 output vector, and Cd, Cv, and Ca are the m×n output influence matrices for displacement, velocity, and acceleration, respectively.Now, the linearization dynamic equation of the manipulator with payload at the selected working point is obtained. In the following section, the identification of the payload parametersme,Je in Equation (12) will be studied by using the system modal matrix estimation. It is worth noting that the linearization Equation (12) can only apply to deal with the small vibration problem at the working point, and it does not suitable to use for the trajectory tracking problem when the manipulator is moving in a wide range.


## 2.1. Dynamic Modeling Based on Lagrange Method

Before establishing the dynamic model of the flexible single-link manipulator, some assumptions are introduced as follows:(1)
The manipulator moves in a plane, and the gravity influence is ignored
(2)
The flexible link complies with the Euler-Bernoulli beam theory, and the axial deformations as well as nonlinear geometric effect due to bending are negligible
(3)
The link’s cross-sectional area remains constant along the link and the material is homogeneous, so the linear density and Young’s modulus of the link are constant
(4)
In the paper, the payload is simplified as a mass point. Therefore, only consider the identification problem of the mass and moment of inertia parameters of the payload
(5)
The flexible of the joint and the system damping are negligible in this paperBased on these assumptions, the manipulator coordinate system and corresponding parameters are defined as follows: the originO of the inertial coordinate frame O‐XYZ and the origin O1 of the link coordinate frames O1‐x1y1z1 are both selected at the joint. The rotation of the link coordinate frame with respect to the inertial coordinate frame O‐XYZ is defined as θ. In addition, the link’s length, mass per unit length, Young modules, and second moment of area are represented by L, ρ, E, and I, respectively. The input torque of joint is τ and the lateral displacement of link is defined as wx,t. The mass and the moment of inertia of the endpoint payload are represented by me and Je, respectively. A simplified model of the flexible single-link manipulator is shown in Figure 1.Figure 1
Single-link manipulator with endpoint payload.If the manipulator is considered as a cantilever beam structure and the first two modal shapes of the link are selected, then by the assumed modes method, the lateral displacementwx,t can be expressed as follows:
(1)wx,t=ϕ1xη1t+ϕ2xη2t,where the two mode shape functions are selected as follows:
(2)ϕix=chβix−cosβix−chβiL+cosβiLshβiL+sinβiLshβix−sinβix,i=1,2.For any point on the link, the position vectorr and corresponding velocity vectors ṙ in inertial coordinate system can be represented, respectively, by:
(3)r=xcosθ−wsinθxsinθ+wcosθ,ṙ=−θ̇xsinθ−ẇsinθ−θ̇wcosθθ̇xcosθ+ẇcosθ−θ̇wsinθ.The system generalized coordinateq is selected as q=θη1η2T, and the kinetic energy TL of the link can be expressed as follows:
(4)TL=12ρ∫0LṙTṙdx=12q̇TρL33+ρLη12+ρLη220.5688ρL20.0908ρL20.5688ρL2ρL00.0908ρL20ρLq̇=12q̇TM1q̇,and the kinetic energy Te of endpoint payload can be denoted as follows:
(5)Te=12meṙeTṙe+12Jeθ̇+w′L,t2=12q̇TmeL2+me2η1−2η22+Je2meL+2.7528Je−4meL−9.5614Je2meL+2.7528Je4me+7.5779Je−4me−26.3206Je−2meL−9.5614Je−4me−26.3206Je4me+91.4204Jeq̇=12q̇TM2q̇,where ṙe is the end position vector of the link. Therefore, the system total kinetic energy T can be expressed as follows:
(6)T=TL+Te=12q̇TM1+M2q̇=12q̇TMq̇,where M is system generalized mass matrix.The system total potential energyU only includes the elastic potential energy for the link and thus it can be written as follows:
(7)U=12EI∫0L∂2w∂x22dx=12qTEI∫0L0000ϕ1″2ϕ1″ϕ2″0ϕ1″ϕ2″ϕ1″2dxq=12qTKq,where K is system generalized stiffness matrix. Consequently, based on the Lagrange method, the dynamic equation of the manipulator can be denoted as follows:
(8)Mq¨+Ṁq̇−∂T∂q+Kq=u,where Ṁq̇ and ∂T/∂q are the Coriolis and Centripetal forces, respectively. u is the generalized force and u=τ00T.Substituting the generalized mass matrix and stiffness matrix into Equation (8), the system dynamic equation can be rewritten as follows:
(9)Mqq¨+Eq,q̇q̇+Kq=u,and the detailed elements in n×n matrices Mq, Eq,q̇, and K are given in Appendix A.


## 2.2. Linearization of Nonlinear Dynamic Model

This paper mainly focuses on the structural vibration and corresponding payload parameter identification problem, so the global motion of the manipulator in the inertial coordinate frameO‐XYZ will be ignored. In this case, based on the local linearization theory, the nonlinear dynamics Equation (9) can be linearized at the selected working point as follows [29]:
(10)Mq0δq¨t+CEδq̇t+K+KM+KEδqt=Lut,where the notation “⋅0” denotes the value at the working point. Therefore, the linearization state vector δqt=qt−q0. L is the input influence matrix and the input signal ut=τt. The linearized matrices CE, KM, and KE in the equation can be computed, respectively, by:
(11)CE=∂E∂q̇q0,q̇0=∂E∂θ̇,∂E∂η̇1,∂E∂η̇2q0,q̇0,KM=∂M∂qq0q¨0=∂M∂θq¨0,∂M∂η1q¨0,∂M∂η2q¨0q0,KE=∂E∂qq0,q̇0=∂E∂θ,∂E∂η1,∂E∂η2q0,q̇0.To simplify the equation derivation, here it is assumed that the modal displacement, velocity, and acceleration of the link at the working point are very small, that is,ηi0=η̇i0=η¨i0≈0 and i=1,2, so the terms CE, KM, and KE in Equation (10) will be zero. Then, the linearized dynamics Equation (10) can be further simplified and the corresponding measurement equation can be expressed, respectively, as follows:
(12)Mq0δq¨t+Kδqt=Lut,(13)yt=Caδq¨t+Cvδq̇t+Cdδqt,where
(14)Mq0=ρL33+meL2+Je0.5688ρL2+2meL+2.7528Je0.0908ρL2−2meL−9.5614Je0.5688ρL2+2meL+2.7528JeρL+4me+7.5779Je−4me−26.3206Je0.0908ρL2−2meL−9.5614Je−4me−26.3206JeρL+4me+91.4204Je,and yt is m×1 output vector, and Cd, Cv, and Ca are the m×n output influence matrices for displacement, velocity, and acceleration, respectively.Now, the linearization dynamic equation of the manipulator with payload at the selected working point is obtained. In the following section, the identification of the payload parametersme,Je in Equation (12) will be studied by using the system modal matrix estimation. It is worth noting that the linearization Equation (12) can only apply to deal with the small vibration problem at the working point, and it does not suitable to use for the trajectory tracking problem when the manipulator is moving in a wide range.


## 3. Payload Parameter Identification Based on the Modal Matrix Estimation

In this section, Equations (12) and (13) are transformed into the state-space form, and the system state-space model and corresponding eigenvalue matrix are identified by the OKID algorithm. Then, the system mass, stiffness, and damping matrices are computed on the basis of the identified eigenvalues. Finally, the inertial parameters of the endpoint payload in the linearized Equation (12) are obtained using the LS technology.


### 3.1. OKID Algorithm

Define a new state vectorxt=δηTtδη̇TtT and then Equations (12) and (13) can be rewritten as the following state-space form:
(15)ẋt=Acxt+Bcut,(16)yt=Cxt+Dut,where Ac, Bc, C, and D are the 2n×2n system, 2n×r input, m×2n output, and m×r direct transmission matrices in the continuous system, respectively:
(17)Ac=0I−Mq0−1K0,Bc=0Mq0−1L,C=Cd−CaM−1K0,D=CaM−1L,where I is the unit matrix. Since the following computations are implemented under discrete systems, Equations (15) and (16) of the continuous system can be further transformed into the discrete state-space equations of the following form as:
(18)xk+1=Axk+Buk,(19)yk=Cxk+Duk,where A and B are the discretized system and input matrices, respectively. Further, the observer form of Equation (18) can be expressed as follows:
(20)xk+1=A¯xk+B¯vk,where
(21)vk=ukyk,where G is an 2n×m arbitrary matrix chosen to make the matrix A¯ as stable as desired. Then, for the time step k+1,k+2,⋯,k+p, the measurement Equation (19) can be expressed, respectively, as follows:
(22)yk+1=CA¯xk+CB¯vk+Duk+1,yk+2=CA¯2xk+CA¯B¯vk+B¯vk+1+Duk+2,⋮yk+p=CA¯pxk+CA¯p−1B¯vk+⋯+CB¯vk+p−1+Duk+p,and the set of these equations for a sequence of k=0,1,⋯,l−1 can be written as follows:
(23)y¯=CA¯px+Y¯V¯,where
(24)y¯=yp,yp+1,⋯,yl−1,Y¯=D,CB¯,CA¯B¯,⋯,CA¯p−1B¯,V¯=upup+1⋯ul−1vp−1vp⋯vl−2⋮⋮⋱⋮v0v1⋯vl−p−1,where D,CB¯,CA¯B¯,L,CA¯p−1B¯ are called observer Markov parameters. The first term in Equation (23) represents the effect of the preceding p−1 time steps. It can select the proper matrix G to ensure that the observer model A¯p is asymptotically stable and all the states in the vector x are bounded, so if the time step is long enough, A¯p≈0 and thus Equation (23) can be approximated by neglecting the first term on the right-hand side as follows:
(25)y¯m×l=Y¯m×pm+r+rV¯pm+r+r×l.Then, the least-squares solution of the observer Markov parameter matrixY¯ in Equation (25) is as follows:
(26)Y¯=y¯V¯†,where the notation “†” is the Moore-Penrose inverse.If define the observer Markov parametersD,CB¯,CA¯B¯,⋯,CA¯p−1B¯ in matrix Y¯ as the following form:
(27)Y¯=D,CB¯,CA¯B¯,⋯,CA¯p−1B¯=Y¯0,Y¯1,Y¯2,⋯,Y¯p,then the observer Markov parameters Y¯i for each time step k=0,1,2,⋯ can be computed by:
(28)Y¯0=D,Y¯k=CA¯k−1B¯=CA+GCk−1B+GD,−CA+GCk−1G≜Y¯k1,−Y¯k2,k=1,2,⋯Based on the definition of Markov parameters, the system Markov parametersYik=1,2,⋯ can be recovered from the observer Markov parameters Y¯i by using the following relationship:
(29)Y0=Y¯0=D,Yk=Yk1−∑i=1kY¯i2Yk−i,fork=1,2,⋯,p,Yk=−∑i=1pY¯i2Yk−i,fork=p+1,p+2,⋯,∞.Using the system Markov parametersYk, the system Hankel matrix H can be written as follows:
(30)Hk−1=YkYk+1⋯Yk+β−1Yk+1Yk+2⋯Yk+β⋮⋮⋱⋮Yk+α−1Yk+α⋯Yk+α+β−2mα×rβ,where α and β are arbitrary positive integers, and the rank of the matrix Hk−1 should be larger than the system order 2n of the state-space model [22].Then, the SVD is implemented for the Hankel matrixH0:
(31)H0=RΣST,where the columns of the matrices R and S are orthonormal, and Σ is a rectangular matrix as follows:
(32)∑=∑2n000,with ∑2n is a diagonal matrix which includes the system singular values σ1,σ2,⋯,σ2n.Therefore, the desired minimum realization matricesA,B,C,D of this state-space system can be obtained by [30]:
(33)A^=∑2n−1/2R2nTH1S2n∑2n−1/2,(34)B^=∑2n1/2S2nTEr,(35)C^=EmTR2n∑2n1/2,where the superscript “^” in Equations (33)–(35) denotes the identified values of relevant matrices that are distinguished from the original values A,B,C, and the matrices R2n, S2n, Er, and Em in Equations (33)–(35) are represented by:
(36)R2n=R:,1:2n,S2n=S:,1:2n,Er=Ir×r,0r×rβ−1T,Em=Im×m,0m×mα−1T,and the expression “:,1:2n” denotes that the first 2n columns of the matrix are selected. The brief computation flow chart of system state-space model by using the OKID algorithm is shown in Figure 2.Figure 2
Computation procedure of the OKID algorithm.


### 3.2. Similarity Transformation of State-Space Model

In the previous section, the OKID algorithm was used to obtain a set of identified discrete model parametersA^,B^,C^. The two sets of state-space models, namely, A^,B^,C^ and A,B,C, have the same system I–O relationship, but the detailed element values in the two sets of models, such as the values of matrices A^ and A, are different. On this basis, a system includes an infinite set of the state-space model Ai,Bi,Ci,i=1,2,⋯,∞. Therefore, when a set of the state-space model is identified by the OKID algorithm, if the detailed elements for the original state-space model are estimated, then a similar transformation should be first performed between the original and identified state-space models by using a transformation matrix to obtain the original model parameters A,B,C [31]. In this case, the identified and original state-space models have the following relationships:
(37)A^=TAT−1,B^=TB,C^=CT−1,where T is the transformation matrix. In practical situation, the original input matrix C is usually priori known. For example, Equation (17) indicates that the matrix C was constructed by the output influence matrices Cd, Cv, and Ca, which were obtained from certain prior knowledge, such as finite element analysis. Therefore, T=CC∧†, and the original state-space model A,B,C was determined by Equation (37).


### 3.3. Identification of System Modal Parameters

The system state-space model Equations (18) and (19) can be decoupled into n pairs of complex conjugate modes by solving the eigenproblem of A. Therefore, the system’s modal parameters (frequencies, damping ratios, and mode shape matrices) can be computed from the system state-space model matrices A,C. Firstly, the eigenvalue decomposition of system matrix A is as follows [32]:
(38)A=Ψ∑Ψ−1,where
(39)∑=ΛdΛd∗∈ℝ2n×2nis a diagonal matrix with complex eigenvalue pairs λj,λj∗,j=1,2,⋯,n, where the notation “∗” denotes the conjugate. Ψ∈ℝ2n×2n is the corresponding complex eigenvector matrix.Then, the complex eigenvaluesλcj,j=1,2,⋯,n for the continuous-time state-space model (15) and (16) can be expressed as follows:
(40)λcj=1Δtlnλj,where Δt is the sampling time. Consequently, the jth normal undamped circular frequency ωj, damping ratio ξj, and experimental mode shapes matrix Φ at the sensor locations for the dynamic system (12) and (13) can be obtained, respectively, as follows:
(41)ωj=Reλcj2+Imλcj2,ξj=ReλcjReλcj2+Imλcj2,Φ=CΨ,where Reλcj and Imλcj denote the real part and imaginary part of the complex eigenvalue λcj, respectively. Therefore, for the continuous-time system model (15) and (16), the complex eigenvalue matrix Λc∈ℝn×n and corresponding complex eigenvector matrix Ψc can be expressed as follows:
(42)Λc=diagλcj,Λc∗=diagλcj∗,Ψc=Ψ=ψψ∗ψΛcψ∗Λc∗,where Ψc is a 2n×2n matrix containing the corresponding complex mode shape ψ∈ℝn×n, which can be expressed as follows:
(43)ψ=ψ1ψ2⋯ψn=ψ11ψ12⋯ψ1nψ21ψ22⋯ψ2n⋮⋯⋱⋮ψn1ψn2⋯ψnn.Then, the complex eigenvalue matrixΛc and complex eigenvector Ψc for the dynamics Equations (12) and (13) are determined by using the OKID algorithm. In the next section, the parameters Λc,Ψc will be used to estimate the inertia parameters me and Je of the endpoint payload object based on complex modal analysis [33, 34].


### 3.4. Inertia Parameter Estimation of Endpoint Payload Object

After the identification of complex eigenvalueΛc and corresponding eigenvector matrices Ψc of the dynamic system (12) and (13), Equation (12) can be written using a canonical form by
(44)0Mq0Mq00δq̇tδq¨t=−K00Mq0δqtδq̇t+L0ut.Furthermore, the generalized eigenproblem can be written in a symmetric form as follows:(45)−K00Mq0ψψΛc=0Mq0Mq00ψψΛcΛc.By simple algebraic manipulations, the orthogonality conditions in Equation (45) can be represented as the following equivalent eigenvalue expression:
(46)Mq0ψΛc2+KψΛc=0.The symmetric formulation in Equation (46) leads to a general solution for the inverse damped vibration problem. Then, the following n complex eigenequations can be obtained:
(47)λcj2Mq0+Kψj=0n×1,j=1,2,⋯,n.For ann-order system as Equation (47), if the generalized form of the elements in the matrices Mq0 and K is defined as follows, then
(48)Mq0=m11m12⋯m1nm21m22⋯m2n⋮⋮⋱⋮mn1mn2⋯mnn,K=k11k22⋱knn,and the complex mode shape ψj∈ℝn×1,j=1,2,⋯,n in Equation (42) can be expressed as follows:
(49)ψj=ψ1jψ2j⋮ψnj,j=1,2,⋯,n.Therefore, the corresponding generalized form of Equation (47) can be further expressed as follows:
(50)λcj2Mq0+Kψj=λcj2m11m12⋯m1nm21m22⋯m2n⋮⋮⋱⋮mn1mn2⋯mnn+k11k22⋱knnψ1jψ2j⋮ψnj=λcj2m11ψ1j+λcj2m12ψ2j+⋯+λcj2m1nψnj+k11ψ1jλcj2m21ψ1j+λcj2m22ψ2j+⋯+λcj2m2nψnj+k22ψ2j⋮λcj2mn1ψ1j+λcj2mn2ψ2j+⋯+λcj2mnnψnj+knnψnj=0n×1,j=1,2,⋯,n.The present study shows the system ordern=3 in Equation (47), and three eigenvalues, namely, λc1, λc2, and λc3, are available. Considering that the element k11 in the stiffness matrix K is equal to zero (the relevant element values in Appendix A), the first eigenvalue λc1 is zero, and the two other eigenvalues λc2 and λc3 are nonzero. Therefore, based on the description elements of mass matrix Mq0 in Equation (12), Equation (50) for the arbitrary nonzero eigenvalues (e.g., λc2) can be denoted as follows:
(51)λc22Mq0+Kψ2=λc22m11ψ12+λc22m12ψ22+λc22m13ψ32+k11ψ12λc22m21ψ12+λc22m22ψ22+λc22m23ψ32+k22ψ22λc22m31ψ12+λc22m32ψ22+λc22m33ψ32+k33ψ32=03×1,where m11, m12,…, m33 can refer to Equation (12) which contains the unknown parameters me and Je that must be identified.In Equation (51), the physical parameters ρ and L of the link and the stiffness k11,k22,k33, respectively, in matrix K are constant. Thus, these parameters can be regarded as known prior knowledge, and the eigenvalue λc2 and complex mode shape ψ2=ψ12ψ22ψ32T can be identified by the OKID method. Consequently, only two unknown payload parameters, me and Je, must be identified. Equation (51) can be further written to the following generalized LS form:
(52)Ξw=δ,where
(53)Ξ≜a11a12a21a22a31a32,w=meJe,δ=b1b2b3.The detailed element values in matrixΞ and vector δ are presented in Appendix B. By using the following LS method, me and Je can then be estimated as follows:
(54)w^=m^eJ^e=Ξ†δ,where the notation “^” denotes the identified value. Therefore, on the basis of Equation (54), the identification of the payload parameters me and Je is reduced to an LS problem with a unique solution.In Equation (51), only two parameters me and Je must be identified in this study and system order n=3 indicates that three modal equations can be used to estimate the two unknown parameters. If certain prior knowledge in the mass matrix Mq0, such as the length L of the link and the coefficients of me and Je, can be obtained by system modeling, then the identified values m^e and J^e can be theoretically calculated only using arbitrary one system eigenvalue (e.g., λc2 or λc3). However, in practical application, the influence of modeling error and the model order truncation is observed, and then the coefficients of me and Je in the mass matrix Mq0 may not be the same as the theory values. Therefore, all the nonzero eigenvalues and corresponding complex mode shapes of the system can also be used to obtain multiple modal equations, as in Equation (51). Then, the LS estimation in Equation (52) can be implemented. This estimation contains all the system modal information to overcome the inaccuracy problem of the matrix coefficients as much as possible.


### 3.5. Summary of Identification Procedures for Payload Parameters

The parameter identification procedures of the endpoint payload object of the manipulator system can be summarized as follows:Step 1.
According to the input and output signals, the system model parametersA,B,C are estimated using the OKID method and matrix similarity transformation by Equations (33), (34), (35), and (37).Step 2.
The complex eigenvaluesλc and complex mode shape ψ of the dynamics Equations (12) and (13) are determined from the identified system model parameters A,B,C by Equations (38), (40), (42), and (43).Step 3.
Using the complex eigenvaluesλc and complex mode shape ψ, the complex eigenequations of the manipulator system operates with payload can be established by Equations (47) and (51).Step 4.
The identified massm^e and moment of inertia J^e of unknown endpoint payload can be computed based on the least-squares estimation of Equation (54).


## 3.1. OKID Algorithm

Define a new state vectorxt=δηTtδη̇TtT and then Equations (12) and (13) can be rewritten as the following state-space form:
(15)ẋt=Acxt+Bcut,(16)yt=Cxt+Dut,where Ac, Bc, C, and D are the 2n×2n system, 2n×r input, m×2n output, and m×r direct transmission matrices in the continuous system, respectively:
(17)Ac=0I−Mq0−1K0,Bc=0Mq0−1L,C=Cd−CaM−1K0,D=CaM−1L,where I is the unit matrix. Since the following computations are implemented under discrete systems, Equations (15) and (16) of the continuous system can be further transformed into the discrete state-space equations of the following form as:
(18)xk+1=Axk+Buk,(19)yk=Cxk+Duk,where A and B are the discretized system and input matrices, respectively. Further, the observer form of Equation (18) can be expressed as follows:
(20)xk+1=A¯xk+B¯vk,where
(21)vk=ukyk,where G is an 2n×m arbitrary matrix chosen to make the matrix A¯ as stable as desired. Then, for the time step k+1,k+2,⋯,k+p, the measurement Equation (19) can be expressed, respectively, as follows:
(22)yk+1=CA¯xk+CB¯vk+Duk+1,yk+2=CA¯2xk+CA¯B¯vk+B¯vk+1+Duk+2,⋮yk+p=CA¯pxk+CA¯p−1B¯vk+⋯+CB¯vk+p−1+Duk+p,and the set of these equations for a sequence of k=0,1,⋯,l−1 can be written as follows:
(23)y¯=CA¯px+Y¯V¯,where
(24)y¯=yp,yp+1,⋯,yl−1,Y¯=D,CB¯,CA¯B¯,⋯,CA¯p−1B¯,V¯=upup+1⋯ul−1vp−1vp⋯vl−2⋮⋮⋱⋮v0v1⋯vl−p−1,where D,CB¯,CA¯B¯,L,CA¯p−1B¯ are called observer Markov parameters. The first term in Equation (23) represents the effect of the preceding p−1 time steps. It can select the proper matrix G to ensure that the observer model A¯p is asymptotically stable and all the states in the vector x are bounded, so if the time step is long enough, A¯p≈0 and thus Equation (23) can be approximated by neglecting the first term on the right-hand side as follows:
(25)y¯m×l=Y¯m×pm+r+rV¯pm+r+r×l.Then, the least-squares solution of the observer Markov parameter matrixY¯ in Equation (25) is as follows:
(26)Y¯=y¯V¯†,where the notation “†” is the Moore-Penrose inverse.If define the observer Markov parametersD,CB¯,CA¯B¯,⋯,CA¯p−1B¯ in matrix Y¯ as the following form:
(27)Y¯=D,CB¯,CA¯B¯,⋯,CA¯p−1B¯=Y¯0,Y¯1,Y¯2,⋯,Y¯p,then the observer Markov parameters Y¯i for each time step k=0,1,2,⋯ can be computed by:
(28)Y¯0=D,Y¯k=CA¯k−1B¯=CA+GCk−1B+GD,−CA+GCk−1G≜Y¯k1,−Y¯k2,k=1,2,⋯Based on the definition of Markov parameters, the system Markov parametersYik=1,2,⋯ can be recovered from the observer Markov parameters Y¯i by using the following relationship:
(29)Y0=Y¯0=D,Yk=Yk1−∑i=1kY¯i2Yk−i,fork=1,2,⋯,p,Yk=−∑i=1pY¯i2Yk−i,fork=p+1,p+2,⋯,∞.Using the system Markov parametersYk, the system Hankel matrix H can be written as follows:
(30)Hk−1=YkYk+1⋯Yk+β−1Yk+1Yk+2⋯Yk+β⋮⋮⋱⋮Yk+α−1Yk+α⋯Yk+α+β−2mα×rβ,where α and β are arbitrary positive integers, and the rank of the matrix Hk−1 should be larger than the system order 2n of the state-space model [22].Then, the SVD is implemented for the Hankel matrixH0:
(31)H0=RΣST,where the columns of the matrices R and S are orthonormal, and Σ is a rectangular matrix as follows:
(32)∑=∑2n000,with ∑2n is a diagonal matrix which includes the system singular values σ1,σ2,⋯,σ2n.Therefore, the desired minimum realization matricesA,B,C,D of this state-space system can be obtained by [30]:
(33)A^=∑2n−1/2R2nTH1S2n∑2n−1/2,(34)B^=∑2n1/2S2nTEr,(35)C^=EmTR2n∑2n1/2,where the superscript “^” in Equations (33)–(35) denotes the identified values of relevant matrices that are distinguished from the original values A,B,C, and the matrices R2n, S2n, Er, and Em in Equations (33)–(35) are represented by:
(36)R2n=R:,1:2n,S2n=S:,1:2n,Er=Ir×r,0r×rβ−1T,Em=Im×m,0m×mα−1T,and the expression “:,1:2n” denotes that the first 2n columns of the matrix are selected. The brief computation flow chart of system state-space model by using the OKID algorithm is shown in Figure 2.Figure 2
Computation procedure of the OKID algorithm.


## 3.2. Similarity Transformation of State-Space Model

In the previous section, the OKID algorithm was used to obtain a set of identified discrete model parametersA^,B^,C^. The two sets of state-space models, namely, A^,B^,C^ and A,B,C, have the same system I–O relationship, but the detailed element values in the two sets of models, such as the values of matrices A^ and A, are different. On this basis, a system includes an infinite set of the state-space model Ai,Bi,Ci,i=1,2,⋯,∞. Therefore, when a set of the state-space model is identified by the OKID algorithm, if the detailed elements for the original state-space model are estimated, then a similar transformation should be first performed between the original and identified state-space models by using a transformation matrix to obtain the original model parameters A,B,C [31]. In this case, the identified and original state-space models have the following relationships:
(37)A^=TAT−1,B^=TB,C^=CT−1,where T is the transformation matrix. In practical situation, the original input matrix C is usually priori known. For example, Equation (17) indicates that the matrix C was constructed by the output influence matrices Cd, Cv, and Ca, which were obtained from certain prior knowledge, such as finite element analysis. Therefore, T=CC∧†, and the original state-space model A,B,C was determined by Equation (37).


## 3.3. Identification of System Modal Parameters

The system state-space model Equations (18) and (19) can be decoupled into n pairs of complex conjugate modes by solving the eigenproblem of A. Therefore, the system’s modal parameters (frequencies, damping ratios, and mode shape matrices) can be computed from the system state-space model matrices A,C. Firstly, the eigenvalue decomposition of system matrix A is as follows [32]:
(38)A=Ψ∑Ψ−1,where
(39)∑=ΛdΛd∗∈ℝ2n×2nis a diagonal matrix with complex eigenvalue pairs λj,λj∗,j=1,2,⋯,n, where the notation “∗” denotes the conjugate. Ψ∈ℝ2n×2n is the corresponding complex eigenvector matrix.Then, the complex eigenvaluesλcj,j=1,2,⋯,n for the continuous-time state-space model (15) and (16) can be expressed as follows:
(40)λcj=1Δtlnλj,where Δt is the sampling time. Consequently, the jth normal undamped circular frequency ωj, damping ratio ξj, and experimental mode shapes matrix Φ at the sensor locations for the dynamic system (12) and (13) can be obtained, respectively, as follows:
(41)ωj=Reλcj2+Imλcj2,ξj=ReλcjReλcj2+Imλcj2,Φ=CΨ,where Reλcj and Imλcj denote the real part and imaginary part of the complex eigenvalue λcj, respectively. Therefore, for the continuous-time system model (15) and (16), the complex eigenvalue matrix Λc∈ℝn×n and corresponding complex eigenvector matrix Ψc can be expressed as follows:
(42)Λc=diagλcj,Λc∗=diagλcj∗,Ψc=Ψ=ψψ∗ψΛcψ∗Λc∗,where Ψc is a 2n×2n matrix containing the corresponding complex mode shape ψ∈ℝn×n, which can be expressed as follows:
(43)ψ=ψ1ψ2⋯ψn=ψ11ψ12⋯ψ1nψ21ψ22⋯ψ2n⋮⋯⋱⋮ψn1ψn2⋯ψnn.Then, the complex eigenvalue matrixΛc and complex eigenvector Ψc for the dynamics Equations (12) and (13) are determined by using the OKID algorithm. In the next section, the parameters Λc,Ψc will be used to estimate the inertia parameters me and Je of the endpoint payload object based on complex modal analysis [33, 34].


## 3.4. Inertia Parameter Estimation of Endpoint Payload Object

After the identification of complex eigenvalueΛc and corresponding eigenvector matrices Ψc of the dynamic system (12) and (13), Equation (12) can be written using a canonical form by
(44)0Mq0Mq00δq̇tδq¨t=−K00Mq0δqtδq̇t+L0ut.Furthermore, the generalized eigenproblem can be written in a symmetric form as follows:(45)−K00Mq0ψψΛc=0Mq0Mq00ψψΛcΛc.By simple algebraic manipulations, the orthogonality conditions in Equation (45) can be represented as the following equivalent eigenvalue expression:
(46)Mq0ψΛc2+KψΛc=0.The symmetric formulation in Equation (46) leads to a general solution for the inverse damped vibration problem. Then, the following n complex eigenequations can be obtained:
(47)λcj2Mq0+Kψj=0n×1,j=1,2,⋯,n.For ann-order system as Equation (47), if the generalized form of the elements in the matrices Mq0 and K is defined as follows, then
(48)Mq0=m11m12⋯m1nm21m22⋯m2n⋮⋮⋱⋮mn1mn2⋯mnn,K=k11k22⋱knn,and the complex mode shape ψj∈ℝn×1,j=1,2,⋯,n in Equation (42) can be expressed as follows:
(49)ψj=ψ1jψ2j⋮ψnj,j=1,2,⋯,n.Therefore, the corresponding generalized form of Equation (47) can be further expressed as follows:
(50)λcj2Mq0+Kψj=λcj2m11m12⋯m1nm21m22⋯m2n⋮⋮⋱⋮mn1mn2⋯mnn+k11k22⋱knnψ1jψ2j⋮ψnj=λcj2m11ψ1j+λcj2m12ψ2j+⋯+λcj2m1nψnj+k11ψ1jλcj2m21ψ1j+λcj2m22ψ2j+⋯+λcj2m2nψnj+k22ψ2j⋮λcj2mn1ψ1j+λcj2mn2ψ2j+⋯+λcj2mnnψnj+knnψnj=0n×1,j=1,2,⋯,n.The present study shows the system ordern=3 in Equation (47), and three eigenvalues, namely, λc1, λc2, and λc3, are available. Considering that the element k11 in the stiffness matrix K is equal to zero (the relevant element values in Appendix A), the first eigenvalue λc1 is zero, and the two other eigenvalues λc2 and λc3 are nonzero. Therefore, based on the description elements of mass matrix Mq0 in Equation (12), Equation (50) for the arbitrary nonzero eigenvalues (e.g., λc2) can be denoted as follows:
(51)λc22Mq0+Kψ2=λc22m11ψ12+λc22m12ψ22+λc22m13ψ32+k11ψ12λc22m21ψ12+λc22m22ψ22+λc22m23ψ32+k22ψ22λc22m31ψ12+λc22m32ψ22+λc22m33ψ32+k33ψ32=03×1,where m11, m12,…, m33 can refer to Equation (12) which contains the unknown parameters me and Je that must be identified.In Equation (51), the physical parameters ρ and L of the link and the stiffness k11,k22,k33, respectively, in matrix K are constant. Thus, these parameters can be regarded as known prior knowledge, and the eigenvalue λc2 and complex mode shape ψ2=ψ12ψ22ψ32T can be identified by the OKID method. Consequently, only two unknown payload parameters, me and Je, must be identified. Equation (51) can be further written to the following generalized LS form:
(52)Ξw=δ,where
(53)Ξ≜a11a12a21a22a31a32,w=meJe,δ=b1b2b3.The detailed element values in matrixΞ and vector δ are presented in Appendix B. By using the following LS method, me and Je can then be estimated as follows:
(54)w^=m^eJ^e=Ξ†δ,where the notation “^” denotes the identified value. Therefore, on the basis of Equation (54), the identification of the payload parameters me and Je is reduced to an LS problem with a unique solution.In Equation (51), only two parameters me and Je must be identified in this study and system order n=3 indicates that three modal equations can be used to estimate the two unknown parameters. If certain prior knowledge in the mass matrix Mq0, such as the length L of the link and the coefficients of me and Je, can be obtained by system modeling, then the identified values m^e and J^e can be theoretically calculated only using arbitrary one system eigenvalue (e.g., λc2 or λc3). However, in practical application, the influence of modeling error and the model order truncation is observed, and then the coefficients of me and Je in the mass matrix Mq0 may not be the same as the theory values. Therefore, all the nonzero eigenvalues and corresponding complex mode shapes of the system can also be used to obtain multiple modal equations, as in Equation (51). Then, the LS estimation in Equation (52) can be implemented. This estimation contains all the system modal information to overcome the inaccuracy problem of the matrix coefficients as much as possible.


## 3.5. Summary of Identification Procedures for Payload Parameters

The parameter identification procedures of the endpoint payload object of the manipulator system can be summarized as follows:Step 1.
According to the input and output signals, the system model parametersA,B,C are estimated using the OKID method and matrix similarity transformation by Equations (33), (34), (35), and (37).Step 2.
The complex eigenvaluesλc and complex mode shape ψ of the dynamics Equations (12) and (13) are determined from the identified system model parameters A,B,C by Equations (38), (40), (42), and (43).Step 3.
Using the complex eigenvaluesλc and complex mode shape ψ, the complex eigenequations of the manipulator system operates with payload can be established by Equations (47) and (51).Step 4.
The identified massm^e and moment of inertia J^e of unknown endpoint payload can be computed based on the least-squares estimation of Equation (54).


## 4. Numerical Simulations

In simulations, a flexible single-link manipulator model is first established, and the designed motor torques and the computed system responses are used as the I–O signals for identification. Subsequently, the system state-space model and corresponding eigenvalue matrix are determined using the OKID algorithm. Finally, the identification of payload parameters is performed using the presented approach, the RLS method, and the APSA algorithm.


### 4.1. Simulation Parameters

The structural parameters of the manipulator are listed in Table1, where the original values of the payload mass me and moment of inertia Je are 100 kg and 200 kg m2, respectively. The system working point q0=θ0,η10,η20 is selected as θ0=30∘ and η10=η20=0.Table 1
Physical parameters of the flexible single-link manipulator for simulations.
Component
Symbol
Value
Unit
Link length
L
1.00
(m)
Length density of link
ρ
19.5
(kg m-1)
Elastic modulus of link
E
20.00
(GPa)
Inertia moment of link
I
5.20×10−7
(m4)
Payload mass
me
100.00
(kg)
Payload moment of inertia
Je
200.00
(kg m2)The input torqueτ of the motor is designed as a sinusoidal signal τ=50sin3.5t0≤t≤1.8s as shown in Figure 3 to ensure that the manipulator motion is near the working point and the original nonlinear model can be linearized. Identification under disturbance is a disadvantage of the OKID method. The OKID method is an LS method that aims at minimizing the error between real and identified outputs. Therefore, if the SNR is low, then the noise usually has high energy. Moreover, distinguishing the noise from the true system signals using the ERA–OKID series method is difficult. Thus, only 5% of Gaussian random measurement noise is added in the system.Figure 3
Designed input torque signalτ=50sin3.5t (0≤t≤1.8s).For the working point {θ0=30∘,η10=η20=0}, the corresponding system responses without and with payload are shown in Figures 4 and 5, respectively. The original nonlinear dynamic equation can be linearized in this I–O condition due to the small system outputs (the response signal of the rotation angle δθ is less than 5°).Figure 4
Response signals of rotation angleδθ and angular velocity δθ̇ at the working point.
(a)
(b)Figure 5
Response signals of modal displacementsδηi and modal velocities δη̇i at the working point i=1,2.
(a)
(b)
(c)
(d)


### 4.2. Identification of System State-Space Model and Frequencies

Using the designed I–O signals in Section4.1, the state-space model of the manipulator with payload is identified using the OKID algorithm, in which the parameters of the OKID algorithm are α=β=5, and the sampling time is Δt=0.001s. Then, the prediction errors of the first three system output signals y1,y2,y3 using the OKID algorithm are provided in Figure 6(a), and corresponding variance is shown in Figure 6(b). Moreover, the singular value magnitude of the system is shown as Figure 7. In Figure 7, it is easy to determine the system order n is n=6 from the truncation of the singular value curve.Figure 6
Prediction error for each step and computation variance of the OKID method.
(a)
(b)Figure 7
Singular values magnitude of the system by the OKID method.According to the similarity transformation mentioned in Section3.2, the different sets of system state-space models satisfy the same I–O relationship. Therefore, the same test inputs are applied to the original and identified state-space models with the zero initial state condition to verify the accuracy of the identified model. The corresponding system responses for the two sets of state-space models at selected working points are shown in Figure 8. The test responses of the estimated state-space models are consistent with those of the original system. The results verify that the identified system model parameters are reasonably accurate.Figure 8
Test response signals of rotation angleδθ, angular velocity δθ̇, 1st modal displacement δη1, and 1st modal velocity δη̇1 at the working point.
(a)
(b)
(c)
(d)Using the identified state-space model, the values of the 2nd and 3rd eigenvalues of the system are computed and presented in Table2, and the corresponding complex mode shapes are provided in Table 3. The values of 1st-order eigenvalue and mode shapes are zero because they denote the rigid rotation motion of the manipulator. The results in Tables 2 and 3, respectively, illustrate that the algorithm can effectively identify the eigenvalues and mode shapes of the manipulator system at the working point.Table 2
Average values of system complex eigenvaluesλc between the original and identified system.
System order
Original values
Identified values
Real
Imag
Real
Imag
2nd
2.5887e−13
2.9452e+01
1.9474e−04
2.9452e+01
3rd
2.8334e−06
9.8878e+02
−2.2964e−07
9.8878e+02Table 3
Average values of system complex mode shapesψ between the original and identified system matrix Ac.
Original values
Identified values
Real
Imag
Real
Imag
Mode 2
−1.9915e−11
−3.3954e−02
−2.4143e−06
−3.2966e−02
4.7128e−12
8.0349e−03
1.6449e−08
7.8011e−03
−1.1913e−12
−2.0313e−03
−3.0973e−09
−1.9722e−03
Mode 3
−2.8823e−12
1.0113e−03
−3.7741e−11
8.8203e−04
1.6052e−12
5.6450e−04
2.9017e−11
4.9232e−04
1.6084e−13
−5.6848e−05
1.3828e−11
4.9579e−05


### 4.3. Identification of Payload Parameters

The proposed method in Section3.4 is implemented once the state-space model and complex eigenvalue parameters are obtained using the OKID algorithm to determine the inertia parameters me,Je of payload. In addition to the proposed algorithm, the RLS and APSA methods are also employed in this simulation to identify the payload parameters.For Equation (52), the standard regression form of the RLS algorithm is as follows:
(55)w^t+1=w^t+Wtδt+1−Ξt+1w^t,Wt+1=ZtΞTt+1γI+Ξt+1ZtΞTt+1−1,Zt+1=1γZt−Wt+1Ξt+1Zt,where I is a unit matrix, and γ is a forgetting factor. Similarly, the APSA method for Equation (52) can also be formulated as follows [30, 35]:
(56)w^t+1=w^t+μXTtsgnΔtsgnΔTtXtXTtsgnΔt+χ,Xt=Ξt;Ξt−1;⋯;Ξt−M+1,Yt=δt;δt−1;⋯;δt−M+1,Δt=Yt−Xtw^t,where μ is a forgetting factor and χ represents the regularization parameter, which should be a positive number.The identification result for each computation may vary due to the influence of the random measurement noise of the system. Thus, a total ofN=10 simulations are conducted for each algorithm. Moreover, the arithmetic mean value (AMV), absolute error (AE), mean absolute deviation (MAD), standard deviation (SD), and relative standard deviation (RSD) of the payload mass parameter me are, respectively, defined as follows:
(57)AMV≜m^¯e=1N∑i=1Nm^ei,AE=m^ei−me,i=1,2,⋯,N,MAD=1N∑i=1Nm^ei−me,SD=1N∑i=1Nm^ei−m^¯e2,RSD=SDAMV×100%,where N is the simulation times, and the corresponding AMV, AE, MAD, SD, and RSD for the moment of inertia Je are similarly defined. The identified payload parameters m^e,J^e and relevant AMV using the three approaches are shown in Table 4. The AEs of the payload parameters using the three methods (10 simulations) are also shown in Figure 9.Table 4
Identification results of the payload parameters (10 simulations).
Simulation no.
Massm^e (kg)
Moment of inertiaJ^e (kg m2)
Proposed method
RLS algorithm
APSA method
Proposed method
RLS algorithm
APSA method
1
99.6703
99.2730
100.1103
201.1777
199.8258
199.4947
2
101.5996
99.2729
101.2102
198.3134
197.8104
198.8921
3
101.6014
98.7180
98.9655
199.2484
197.9906
195.2928
4
102.5959
99.2093
101.0607
202.4915
199.5080
198.9590
5
99.4926
97.8406
99.3710
199.5160
197.7941
199.4952
6
102.4911
95.3401
98.2425
199.4401
195.5997
198.7861
7
101.3892
95.4243
99.8859
198.4393
196.6997
198.0171
8
100.4916
98.0171
97.3875
200.4166
196.9188
199.8103
9
101.5955
98.9264
102.8673
199.3512
198.0419
198.8404
10
98.9814
99.0898
101.0377
198.4253
199.3663
199.9468
AMV
100.9909
98.1111
100.0139
199.6820
197.9555
198.7535Figure 9
Absolute error (AE) of the payload parametersm^e and J^e (10 simulations for each algorithm).The AMV results in Table4 demonstrate that the proposed method can estimate the payload parameters more accurately than the other approaches. In Figure 9, the circular lines for the three methods are used to distinguish the maximum error distribution of the payload parameters around the original value (namely, the center of the figure). Similarly, Figure 9 also indicates that the proposed method better approximates the original value of the system than the two other classical algorithms.The MAD, SD, and RSD of the payload parameters for the three methods are shown in Figures10 and 11. It can be seen that the proposed method has satisfactory identification accuracy and the corresponding RSD is less than 5%. Figures 9–11 demonstrate that and the computational errors by the proposed method are smaller than those by the classical RLS and APSA algorithms.Figure 10
Computation results of the payload parametersm^e and J^e for the three methods: (a) mean absolute deviation and (b) standard deviation (10 simulations).
(a)
(b)Figure 11
Relative standard deviation result of the payload parametersm^e and J^e for the three methods (10 simulations).


## 4.1. Simulation Parameters

The structural parameters of the manipulator are listed in Table1, where the original values of the payload mass me and moment of inertia Je are 100 kg and 200 kg m2, respectively. The system working point q0=θ0,η10,η20 is selected as θ0=30∘ and η10=η20=0.Table 1
Physical parameters of the flexible single-link manipulator for simulations.
Component
Symbol
Value
Unit
Link length
L
1.00
(m)
Length density of link
ρ
19.5
(kg m-1)
Elastic modulus of link
E
20.00
(GPa)
Inertia moment of link
I
5.20×10−7
(m4)
Payload mass
me
100.00
(kg)
Payload moment of inertia
Je
200.00
(kg m2)The input torqueτ of the motor is designed as a sinusoidal signal τ=50sin3.5t0≤t≤1.8s as shown in Figure 3 to ensure that the manipulator motion is near the working point and the original nonlinear model can be linearized. Identification under disturbance is a disadvantage of the OKID method. The OKID method is an LS method that aims at minimizing the error between real and identified outputs. Therefore, if the SNR is low, then the noise usually has high energy. Moreover, distinguishing the noise from the true system signals using the ERA–OKID series method is difficult. Thus, only 5% of Gaussian random measurement noise is added in the system.Figure 3
Designed input torque signalτ=50sin3.5t (0≤t≤1.8s).For the working point {θ0=30∘,η10=η20=0}, the corresponding system responses without and with payload are shown in Figures 4 and 5, respectively. The original nonlinear dynamic equation can be linearized in this I–O condition due to the small system outputs (the response signal of the rotation angle δθ is less than 5°).Figure 4
Response signals of rotation angleδθ and angular velocity δθ̇ at the working point.
(a)
(b)Figure 5
Response signals of modal displacementsδηi and modal velocities δη̇i at the working point i=1,2.
(a)
(b)
(c)
(d)


## 4.2. Identification of System State-Space Model and Frequencies

Using the designed I–O signals in Section4.1, the state-space model of the manipulator with payload is identified using the OKID algorithm, in which the parameters of the OKID algorithm are α=β=5, and the sampling time is Δt=0.001s. Then, the prediction errors of the first three system output signals y1,y2,y3 using the OKID algorithm are provided in Figure 6(a), and corresponding variance is shown in Figure 6(b). Moreover, the singular value magnitude of the system is shown as Figure 7. In Figure 7, it is easy to determine the system order n is n=6 from the truncation of the singular value curve.Figure 6
Prediction error for each step and computation variance of the OKID method.
(a)
(b)Figure 7
Singular values magnitude of the system by the OKID method.According to the similarity transformation mentioned in Section3.2, the different sets of system state-space models satisfy the same I–O relationship. Therefore, the same test inputs are applied to the original and identified state-space models with the zero initial state condition to verify the accuracy of the identified model. The corresponding system responses for the two sets of state-space models at selected working points are shown in Figure 8. The test responses of the estimated state-space models are consistent with those of the original system. The results verify that the identified system model parameters are reasonably accurate.Figure 8
Test response signals of rotation angleδθ, angular velocity δθ̇, 1st modal displacement δη1, and 1st modal velocity δη̇1 at the working point.
(a)
(b)
(c)
(d)Using the identified state-space model, the values of the 2nd and 3rd eigenvalues of the system are computed and presented in Table2, and the corresponding complex mode shapes are provided in Table 3. The values of 1st-order eigenvalue and mode shapes are zero because they denote the rigid rotation motion of the manipulator. The results in Tables 2 and 3, respectively, illustrate that the algorithm can effectively identify the eigenvalues and mode shapes of the manipulator system at the working point.Table 2
Average values of system complex eigenvaluesλc between the original and identified system.
System order
Original values
Identified values
Real
Imag
Real
Imag
2nd
2.5887e−13
2.9452e+01
1.9474e−04
2.9452e+01
3rd
2.8334e−06
9.8878e+02
−2.2964e−07
9.8878e+02Table 3
Average values of system complex mode shapesψ between the original and identified system matrix Ac.
Original values
Identified values
Real
Imag
Real
Imag
Mode 2
−1.9915e−11
−3.3954e−02
−2.4143e−06
−3.2966e−02
4.7128e−12
8.0349e−03
1.6449e−08
7.8011e−03
−1.1913e−12
−2.0313e−03
−3.0973e−09
−1.9722e−03
Mode 3
−2.8823e−12
1.0113e−03
−3.7741e−11
8.8203e−04
1.6052e−12
5.6450e−04
2.9017e−11
4.9232e−04
1.6084e−13
−5.6848e−05
1.3828e−11
4.9579e−05


## 4.3. Identification of Payload Parameters

The proposed method in Section3.4 is implemented once the state-space model and complex eigenvalue parameters are obtained using the OKID algorithm to determine the inertia parameters me,Je of payload. In addition to the proposed algorithm, the RLS and APSA methods are also employed in this simulation to identify the payload parameters.For Equation (52), the standard regression form of the RLS algorithm is as follows:
(55)w^t+1=w^t+Wtδt+1−Ξt+1w^t,Wt+1=ZtΞTt+1γI+Ξt+1ZtΞTt+1−1,Zt+1=1γZt−Wt+1Ξt+1Zt,where I is a unit matrix, and γ is a forgetting factor. Similarly, the APSA method for Equation (52) can also be formulated as follows [30, 35]:
(56)w^t+1=w^t+μXTtsgnΔtsgnΔTtXtXTtsgnΔt+χ,Xt=Ξt;Ξt−1;⋯;Ξt−M+1,Yt=δt;δt−1;⋯;δt−M+1,Δt=Yt−Xtw^t,where μ is a forgetting factor and χ represents the regularization parameter, which should be a positive number.The identification result for each computation may vary due to the influence of the random measurement noise of the system. Thus, a total ofN=10 simulations are conducted for each algorithm. Moreover, the arithmetic mean value (AMV), absolute error (AE), mean absolute deviation (MAD), standard deviation (SD), and relative standard deviation (RSD) of the payload mass parameter me are, respectively, defined as follows:
(57)AMV≜m^¯e=1N∑i=1Nm^ei,AE=m^ei−me,i=1,2,⋯,N,MAD=1N∑i=1Nm^ei−me,SD=1N∑i=1Nm^ei−m^¯e2,RSD=SDAMV×100%,where N is the simulation times, and the corresponding AMV, AE, MAD, SD, and RSD for the moment of inertia Je are similarly defined. The identified payload parameters m^e,J^e and relevant AMV using the three approaches are shown in Table 4. The AEs of the payload parameters using the three methods (10 simulations) are also shown in Figure 9.Table 4
Identification results of the payload parameters (10 simulations).
Simulation no.
Massm^e (kg)
Moment of inertiaJ^e (kg m2)
Proposed method
RLS algorithm
APSA method
Proposed method
RLS algorithm
APSA method
1
99.6703
99.2730
100.1103
201.1777
199.8258
199.4947
2
101.5996
99.2729
101.2102
198.3134
197.8104
198.8921
3
101.6014
98.7180
98.9655
199.2484
197.9906
195.2928
4
102.5959
99.2093
101.0607
202.4915
199.5080
198.9590
5
99.4926
97.8406
99.3710
199.5160
197.7941
199.4952
6
102.4911
95.3401
98.2425
199.4401
195.5997
198.7861
7
101.3892
95.4243
99.8859
198.4393
196.6997
198.0171
8
100.4916
98.0171
97.3875
200.4166
196.9188
199.8103
9
101.5955
98.9264
102.8673
199.3512
198.0419
198.8404
10
98.9814
99.0898
101.0377
198.4253
199.3663
199.9468
AMV
100.9909
98.1111
100.0139
199.6820
197.9555
198.7535Figure 9
Absolute error (AE) of the payload parametersm^e and J^e (10 simulations for each algorithm).The AMV results in Table4 demonstrate that the proposed method can estimate the payload parameters more accurately than the other approaches. In Figure 9, the circular lines for the three methods are used to distinguish the maximum error distribution of the payload parameters around the original value (namely, the center of the figure). Similarly, Figure 9 also indicates that the proposed method better approximates the original value of the system than the two other classical algorithms.The MAD, SD, and RSD of the payload parameters for the three methods are shown in Figures10 and 11. It can be seen that the proposed method has satisfactory identification accuracy and the corresponding RSD is less than 5%. Figures 9–11 demonstrate that and the computational errors by the proposed method are smaller than those by the classical RLS and APSA algorithms.Figure 10
Computation results of the payload parametersm^e and J^e for the three methods: (a) mean absolute deviation and (b) standard deviation (10 simulations).
(a)
(b)Figure 11
Relative standard deviation result of the payload parametersm^e and J^e for the three methods (10 simulations).


## 5. Conclusions

A payload parameter identification method of a flexible space manipulator system is proposed in this study. By linearizing the nonlinear dynamic equation of the manipulator at selected working points, the system state-space model and the corresponding eigenvalue parameters are, respectively, determined using the OKID algorithm and the similarity transformation. Moreover, the payload parameters, namely, the mass and moment of inertia, can be derived by solving an LS problem based on the complex modal analysis of the system. The payload parameters obtained in the simulations by the proposed method are compared with those by the classical RLS and APSA algorithms. The results demonstrate that the proposed method can successfully be used to identify the manipulator payload parameters.Certain problems related to the proposed method should be improved. (1) Only the planar situation is considered in this study, and the endpoint payload is simplified as a mass point. Therefore, the identification of the mass center is ignored. (2) In the identification procedures, the prior knowledge of model parameters, such as the coefficients in the mass matrix, is required. Otherwise, the accuracy of the solution will be affected. (3) Model order truncation is conducted in the established manipulator system, and only the first two modal shapes of the link in Equation (1) are selected. However, if the number of selected vibration modes of the link is increased, then the establishment of the dynamic model would be considerably complicated. In future works, the extension of the proposed method to three-dimensional objects should be investigated, and a study will be conducted to reduce the effect of model truncation on the identification results.



---

*Source: 5142925-2020-02-05.xml*