# Simultaneous Identification of Model Structure and the Associated Parameters for Linear Systems Based on Particle Swarm Optimization

Authors: Semin Chun; Tae-Hyoung Kim
Journal: Complexity
 (2018)

---

## Abstract

In this study, a novel easy-to-use meta-heuristic method for simultaneous identification of model structure and the associated parameters for linear systems is developed. This is achieved via a constrained multidimensional particle swarm optimization (PSO) mechanism developed by hybridizing two main methodologies: one for negating the limit for fixing the particle’s dimensions within the PSO process and another for enhancing the exploration ability of the particles by adopting a cyclic neighborhood topology of the swarm. This optimizer consecutively searches the dimensional optimum of particles and then the positional optimum in the search space, whose dimension is specified by the explored optimal dimension. The dimensional optimum provides the optimal model structure, while the positional optimum provides the optimal model parameters. Typical numerical examples are considered for evaluation purposes, which clearly demonstrate that the proposed PSO scheme provides novel and powerful impetus with remarkable reliability toward simultaneous identification of model structure and unknown model parameters. Furthermore, identification experiments are conducted on a magnetic levitation system and a robotic manipulator with joint flexibility to demonstrate the effectiveness of the proposed strategy in practical applications.


---

## Body

## 1. Introduction

One of the major research directions in modern control theory has been directed at developing simple yet reliable control strategies to tackle the increasing complexity of practical control problems. From this viewpoint, model-based control has been widely accepted as a crucial tool for achieving the target system’s precise control performance. A major prerequisite in developing an advanced model-based controller is therefore to obtain a highly reliable dynamic model upon which the controller design procedure will be based. Identification is a powerful technique for building a precise dynamic mathematical model of the target system from input-output data corrupted by noises. It usually involves the following three interrelated steps: design of an experiment, construction of a candidate dynamic model structure, and estimation of the model parameters from the measurements [1].The dynamic model is usually subject to uncertainties in the parameters involved, and the performance of the developed model-based controller is sensitive to the associated parameters’ values. Therefore, some well-known parameter estimation schemes, which depend on the type of dynamic model one is dealing with for system identification, have gained widespread attention for developing model-based controllers [2–6]. Although canonical estimation methodology has been applied extensively and successfully, the typical nonconvex nature of the criterion function with many local minima, on which the parameter estimation algorithm depends, may cause significant numerical search problems [7]. Therefore, recently, an alternative research trend that uses computing techniques including swarm intelligence for parameter estimation has emerged [8]. For instance, [9] recently proposed the hybrid particle swarm optimization (PSO) and gravitational search algorithm for parameter estimation of infinite impulse response (IIR) systems. Kumar and Rawat [10] introduced the cuckoo search algorithm to estimate the optimal coefficients of finite impulse response-fractional order differentiator models. Mostajabi et al. [11] conducted a comparative study demonstrating that two well-known population-based optimization techniques, genetic algorithm and PSO, can perform better than the recursive least square algorithm in parameter estimation problems for IIR models. Pal et al. [12] examined the parameter estimation of a linear unstable system by using the craziness-based PSO algorithm. Peng and Wang [13] proposed a hybrid method combining the tissue P system and artificial bee colony for IIR system identification. Alsmadi et al. [14] examined some population-based solution approaches—invasive weed optimization, firefly algorithm, GA, and PSO—and their applicability to model order reduction problems. Dhaliwal and Dhillon [15] employed the cat swarm optimization algorithm with the incorporation of oppositional learning strategy for designing an optimal and stable digital IIR band-pass filter. These studies verified that swarm-based optimization techniques can achieve good performance in the target system’s parameter estimation problems. Note, however, that all aforementioned methods examine only the performance of parameter estimation in case the model order and structure are assumed to be known a priori.In the system identification process, the construction of a candidate dynamic model structure/order is clearly the most important and most difficult step [4]. The model structure selection may compromise the optimality properties [7]. The extensively used natural choice is to base the model on physical insight or interpretation, which leads to the so-called white-box model structure, which relates the laws of physics with the corresponding physical parameters. If some of these physical parameters are subject to uncertainties or not well known, the so-called gray-box modeling approach is usually adopted to preserve the physical meaning of the uncertain model parameters. On the other hand, in control applications, it usually suffices to adopt linear models that do not necessarily refer to the underlying theoretical physical knowledge of the target system. Such a model is usually called the black-box model [4]. This viewpoint clearly indicates that the aforementioned swarm-based methods [9–15] can be categorized as the gray-box modeling technique. Note that most other previously proposed swarm-based optimization methods mainly focus on the estimation of unknown parameters in case the model structure/order is fixed a priori, which implies that these conventional methods are not directly applicable to the black-box modeling problems. The development of an efficient black-box modeling methodology, in which swarm-based optimization is incorporated, is thus necessary; however, until recently, very little research has been conducted regarding the same. Mohd Azmi et al. [16–18] and Ibrahim et al. [19] proposed the identification methods using multiswarm-based optimization techniques for simultaneously estimating the model order and uncertain parameters involved. Their strategy is described briefly as follows. Let the maximum order of the appropriate mathematical models (transfer function models) be N, which means that one best model is chosen among the models of orders (i.e., the order of denominators) from 1 to N. Note that the nth-order transfer model consists of n types of numerator structures. Therefore, NN+1/2 types of models are examined for finding one that best matches the response data. Then, their core strategy for black-box identification is one-to-one correspondence between multiswarm optimizers consisting of N-independent swarms and model sets classified by orders ranging from 1 to N. This simply means that each swarm optimizer is assigned to examine models that have an identical particular order, and then multiple swarm optimizers run independently in parallel to achieve their respective identification tasks. The above scheme, however, may not be a refined one from the viewpoint that its mechanism is equivalent to the methodology such that a single swarm-based optimizer performs the identification tasks consecutively from the 1st-order model to the Nth-order model. To the best of our knowledge, no existing study, except for the above ones [16–19], has yet adequately addressed a complete method based on swarm-based optimization techniques for black-box identification, including optimal model structure estimation directly from input-output data.System identification is the procedure of deriving a mathematical model from input-output measurement data, and its first step is model estimation, which is the procedure of fitting a model with a specific model structure. The variation in the orders of numerator and denominator polynomials, which denotes the transfer function model structure variation, results in variation in the total number of uncertain model parameters to be identified. Since the particles’ dimension in the PSO algorithm equals the total number of model parameters, it should be variable according to a model structure to be evaluated in model estimation problems. However, many conventional PSO algorithms have a critical drawback, in that the particle dimensions remain fixed throughout the iterations. Such a drawback impedes the application of many PSO variants to model structure estimation studies in the field of system identification. To deal with the above difficulty, this study develops a novel constrained multidimensional PSO (MD-PSO) with a cyclic network topology-based neighborhood structure. The proposed optimization scheme specialized for system identification hybridizes two methodologies: one for negating the limit for fixing the particle dimensions within the PSO process and another for enhancing the exploration ability of the particles by adopting a cyclic network topology-based neighborhood structure. This novel scheme modifies the original MD-PSO algorithm of [20] so as to be compatible with the cyclic neighborhood topology [21] and then to be ultimately applicable to the simultaneous identification of the model structure and the associated parameters for linear systems.The key features of our hybridized PSO algorithm are as follows. The introduced cyclic network topology-based neighborhood structure of particles can effectively regulate the particles’ exploration abilities and enable them to locate a promising region in the hyperdimensional search space, which was originally proposed by one of the authors of the present paper [21]. Such an improved exploration ability of the particles is essential to cope with the typical nonconvex nature of the criterion function with many local minima in the process of system identification [7]. On the other hand, our PSO algorithm searches both (i) the dimensional optimum of particles within the preassigned dimension range and (ii) the positional optimum in the hyperdimensional search space, whose dimension is specified by the aforementioned optimal dimension of particles. Note that the dimensional optimum provides the optimal orders of numerator and denominator polynomials in the transfer function, while the positional optimum provides the optimal values of such numerator/denominator polynomials’ unknown parameters. The overall PSO procedure for finding an optimal model is briefly summarized as follows. First, the dimensional optimum of particles is searched via the optimal model structure-searching PSO algorithm. To apply this PSO mechanism to the model structure estimation problem, a candidate model structure is labeled with a unique model identifier number in advance. This enables us to transform the problem of searching the dimensional optimum of particles into that for exploring the optimal model identifier number (Table 1). Then, the optimal model structure-searching PSO process moves a particle staying at a position providing a certain model identifier number information to another position for renewing the previous model identifier number. The speed and position update formula in this PSO process follow the particle swarm optimizer with cyclic neighborhood topology, where each particle shares information through a fixed nearby-neighbor interaction structure with a series of successively numbered particles. The unique orders of the numerator and denominator polynomials of a candidate model are extracted from the model identifier number that equals to the renewed position of a particle used in the optimal model structure-searching PSO algorithm. Such information about the polynomials’ orders defines the total number of unknown parameters to be identified, which enables one to set the particle dimensions for exploring the positional optimum (i.e., optimal model parameters). Once the particle’s dimension is determined, its positional optimum is searched using the optimal model parameter-searching PSO algorithm. This optimizer also adopts the cyclic neighborhood topology of the swarm to effectively regulate the exploration abilities of the particles in the hyperdimensional search space, whose dimension is determined from the particle’s position renewed via the optimal model structure-searching PSO algorithm as a result. Next, the optimality of the particle’s dimension (i.e., the estimated orders of numerator and denominator polynomials) and position (i.e., the estimated model parameters) obtained is evaluated using a criterion function derived from the system identification viewpoint. Finally, the two aforementioned consecutive optimizers, the optimal model structure-searching and optimal model parameter-searching PSO algorithms, are run iteratively until the stopping criterion is reached. The proposed methodology provides a novel tool for simultaneous model structure and its parameter identification and enables one to easily overcome difficulties such as over- and underestimation of a model. The superior identification capability of the proposed method is verified through simulation examples. In addition, two experimental setups for the magnetic levitation system and the robotic manipulator with joint flexibility are established for system identification. The experimental results for a simultaneous identification of both the model structure and unknown model parameters demonstrate the superior reliability and validity of the proposed method.Table 1
Model structure of then^th-order transfer function (Dmin=1 case)
Then^th-order transfer function model
Model structure
b0sn^+a1sn^−1+⋯+an^
b0s+b1sn^+a1sn^−1+⋯+an^
…
b0sn^−2+b1sn^−3+⋯+bn^−2sn^+a1sn^−1+⋯+an^
b0sn^−1+b1sn^−2+⋯+bn^−1sn^+a1sn^−1+⋯+an^
Order of the numerator polynomial
0
1
…
n^−2
n^−1
Number of the unknown parameters
n^+1
n^+2
…
2n^−1
2n^
Model identifier number
n^n^−12+1
n^n^−12+2
…
n^n^−12+n^−1
n^n^−12+n^


## 2. Problem Formulation

Assume that a finite-dimensional, linear time-invariant, continuous-time plant to be identified is described in the compact transfer function formyt=G∗put with p derivative operator (i.e., pαzt=dαzt/dtα), where yt and ut are, respectively, the noise-free response and the excitation. Let G∗s be a rational form in s, defined as
(1)G∗s=B∗sA∗s,where s is the Laplace transform variable, and A∗s and B∗s are assumed to be coprime and
(2)A∗s=sn∗+a1∗sn∗−1+a2∗sn∗−2+⋯+an∗∗,(3)B∗s=b0∗sm∗+b1∗sm∗−1+b2∗sm∗−2+⋯+bm∗∗,n∗>m∗.Letθ∗ denote the parameter vector, which includes the dynamic plant model’s unknown parameters stacked column-wise as
(4)θ∗≔a1∗,a2∗,…,an∗∗,b0∗,b1∗,…,bm∗∗T∈ℝn∗+m∗+1.A model structure is defined as a parameterized collection of models that describe the relations between the input and output signals of the plant [7]. In several practical situations, no prior information of the plant model structure is given, and thus, orders n∗ and m∗ for the numerator and denominator polynomials of G∗s are assumed to be unknown. From this viewpoint, our system identification for deriving a model from the measured input and output data is composed of two procedures: (ID-1) model structure estimation, which is the procedure of fitting a model form G^s with a specific structure G∗s, and (ID-2) model parameter estimation, which is the procedure through which the parameters of a candidate model G^s are updated by matching the output data calculated from the model predictions to the actual measurement data. Note that the model structure estimation results in the problem of finding an optimal order estimation n^,m^ of n∗,m∗. The model parameter estimation results in the problem of finding an optimal estimation θ^, whose dimension depends on n^,m^, that is,
(5)θ^n^,m^≔a^1,a^2,…,a^n^,b^0,b^1,…,b^m^∈ℝn^+m^+1.A set of candidate models (i.e., model structure) constructed to be used in the identification procedure is described as follows. Consider then^th-order transfer function model with Dmin≤n^≤Dmax, where Dmin and Dmax are assumed to be given in advance. Since the model order is equal to the order of the denominator polynomial, the n^th-order model has a total of n^ types of numerator polynomials, as shown in Table 1. Let the model identifier number be denoted by
(6)IdNon^,m^≔n^n^−12+m^+1−DminDmin−12,where n^, 0≤m^<n^, is the order of the numerator polynomial. The lowest and highest identifier numbers are, respectively, IdNoDmin,0 and IdNoDmax,Dmax−1, and the total number of candidate models to be evaluated equals to IdNoDmax,Dmax−1−IdNoDmin,0+1. Therefore, the aforementioned model structure estimation problem (ID-1) becomes the problem of finding an optimal IdNon^,m^. Then, the n^th-order transfer function model with order m^ of the numerator polynomial is denoted as G^s,θ^n^,m^, where θ^n^,m^≔a^1,a^2,…,a^n^,b^0,b^1,…,b^m^T∈ℝn^+m^+1 is a vector of unknown parameters in the numerator and denominator polynomials.


### 2.1. Criterion Functions for Fitting the Time-Domain Data

In the case of time-domain identification, the criterion functions (i.e., objective functions) used in the optimization procedure for continuous-time model structure/parameter estimation directly from the sampled data are formulated as follows. For astable target plant, the input signal ut,t0≤t≤tN is applied to the plant G∗p, which gives rise to an output signal yt,t0≤t≤tN. Assuming that N+1 sampled measurements of the input-output signal utk and y~tk≔ytk+wtk, where k=0,1,…,N and wts is a stochastic noise sequence with zero mean, are available.For astable plant G∗p, the criterion function introduced for the open-loop time-domain estimation of both the model structure G^p and its uncertain parameter vector θ^n^,m^ from the sampled measurement data is then formulated as
(7)fOL−timeθ^n^,m^≔∑k=0Ny~OL−planttk−yOL−modeltk21/2,subject to
(8)ReλmaxG^p,θ^n^,m^<0,where y~OL−plantt≔G∗put+wt is the measured output, yOL−modelt≔G^p,θ^n^,m^ut is the simulated output, and λmax· denotes the pole with the greatest real part. Note that (8) is the constraint condition for guaranteeing the model stability.On the other hand, if the target plant has inherentlyunstable nature, the usual open-loop identification methodology cannot be adopted. In such situations, experimental data can only be obtained under closed-loop conditions, and then, the criterion function for time-domain identification becomes
(9)fCL−timeθ^n^,m^≔∑k=0Ny~CL−planttk−yCL−modeltk21/2.Here,yCL−modelt and y~CL−plantt are, respectively, defined as yCL−modelt≔M^p,θ^n^,m^ut and y~CL−plantt≔M∗put+wt, where the closed-loop systems M^·,· and M∗p stabilized by an identical known controller Cp are defined as
(10)M^p,θ^n^,m^≔G^p,θ^n^,m^Cp1+G^p,θ^n^,m^Cp,M∗p≔G∗pCp1+G∗pCp.Similarly to the case of an open-loop time-domain identification problem, the following constraint condition is introduced to guarantee the stability of the closed loop:(11)ReλmaxM^p,θ^n^,m^<0.


### 2.2. Criterion Functions for Fitting the Frequency-Domain Data

Frequency-domain identification is another methodology for developing a dynamic model, which has the following key advantages [22]. First, the output measurement noise does not bias the frequency response estimates, given that the noise is uncorrelated with the control input. Second, the frequency range used to fit each input-output pair can be selected individually. Thus, only the most accurate measurement data are involved in the identification process. Third, selection of a specific frequency range can be used as an effective tool for separating data that are relevant for system identification from irrelevant data such as noise and disturbances. Therefore, the frequency-domain technique enables the modeling of a continuous-time system from the sampled data in a straightforward fashion for a certain class of band-limited excitation signals [23].For astable plant, the criterion function introduced for the open-loop frequency-domain estimation of the model structure/parameters becomes
(12)fOL−freqθ^n^,m^≔∑k=0NG∗jωkUjωk−YOL−modeljωk21/2,YOL−modeljωUjω=G^jω,θ^n^,m^,subject to the constraint condition (8). In (12), G∗jω denotes the frequency response function of the plant, G^jω,θ^n^,m^ is the estimated frequency response function, Ujω and Yjω are the Fourier transforms of the input and output signals, and ωk (k=0,1,…,N) are the given discrete frequencies.On the other hand, the criterion function forfrequency-domain closed-loop identification of an unstable plant is formulated as
(13)fCL−freqθ^n^,m^≔∑k=0NM∗jωkUjωk−YCL−modeljω21/2,YCL−modeljωUjω=M^jω,θ^n^,m^,where the constraint condition is identical to (11).


## 2.1. Criterion Functions for Fitting the Time-Domain Data

In the case of time-domain identification, the criterion functions (i.e., objective functions) used in the optimization procedure for continuous-time model structure/parameter estimation directly from the sampled data are formulated as follows. For astable target plant, the input signal ut,t0≤t≤tN is applied to the plant G∗p, which gives rise to an output signal yt,t0≤t≤tN. Assuming that N+1 sampled measurements of the input-output signal utk and y~tk≔ytk+wtk, where k=0,1,…,N and wts is a stochastic noise sequence with zero mean, are available.For astable plant G∗p, the criterion function introduced for the open-loop time-domain estimation of both the model structure G^p and its uncertain parameter vector θ^n^,m^ from the sampled measurement data is then formulated as
(7)fOL−timeθ^n^,m^≔∑k=0Ny~OL−planttk−yOL−modeltk21/2,subject to
(8)ReλmaxG^p,θ^n^,m^<0,where y~OL−plantt≔G∗put+wt is the measured output, yOL−modelt≔G^p,θ^n^,m^ut is the simulated output, and λmax· denotes the pole with the greatest real part. Note that (8) is the constraint condition for guaranteeing the model stability.On the other hand, if the target plant has inherentlyunstable nature, the usual open-loop identification methodology cannot be adopted. In such situations, experimental data can only be obtained under closed-loop conditions, and then, the criterion function for time-domain identification becomes
(9)fCL−timeθ^n^,m^≔∑k=0Ny~CL−planttk−yCL−modeltk21/2.Here,yCL−modelt and y~CL−plantt are, respectively, defined as yCL−modelt≔M^p,θ^n^,m^ut and y~CL−plantt≔M∗put+wt, where the closed-loop systems M^·,· and M∗p stabilized by an identical known controller Cp are defined as
(10)M^p,θ^n^,m^≔G^p,θ^n^,m^Cp1+G^p,θ^n^,m^Cp,M∗p≔G∗pCp1+G∗pCp.Similarly to the case of an open-loop time-domain identification problem, the following constraint condition is introduced to guarantee the stability of the closed loop:(11)ReλmaxM^p,θ^n^,m^<0.


## 2.2. Criterion Functions for Fitting the Frequency-Domain Data

Frequency-domain identification is another methodology for developing a dynamic model, which has the following key advantages [22]. First, the output measurement noise does not bias the frequency response estimates, given that the noise is uncorrelated with the control input. Second, the frequency range used to fit each input-output pair can be selected individually. Thus, only the most accurate measurement data are involved in the identification process. Third, selection of a specific frequency range can be used as an effective tool for separating data that are relevant for system identification from irrelevant data such as noise and disturbances. Therefore, the frequency-domain technique enables the modeling of a continuous-time system from the sampled data in a straightforward fashion for a certain class of band-limited excitation signals [23].For astable plant, the criterion function introduced for the open-loop frequency-domain estimation of the model structure/parameters becomes
(12)fOL−freqθ^n^,m^≔∑k=0NG∗jωkUjωk−YOL−modeljωk21/2,YOL−modeljωUjω=G^jω,θ^n^,m^,subject to the constraint condition (8). In (12), G∗jω denotes the frequency response function of the plant, G^jω,θ^n^,m^ is the estimated frequency response function, Ujω and Yjω are the Fourier transforms of the input and output signals, and ωk (k=0,1,…,N) are the given discrete frequencies.On the other hand, the criterion function forfrequency-domain closed-loop identification of an unstable plant is formulated as
(13)fCL−freqθ^n^,m^≔∑k=0NM∗jωkUjωk−YCL−modeljω21/2,YCL−modeljωUjω=M^jω,θ^n^,m^,where the constraint condition is identical to (11).


## 3. Constrained MD-PSO with a Cyclic Neighborhood Topology for System Identification

Let the notationsPiSTR and PiPAR be defined, respectively, as follows: (i) PiSTR: the ith particle in the swarm used for searching the optimal model structure (i.e., an optimal model identifier number IdNon^,m^) and (ii) PiPAR: the ith particle in the swarm used for identifying the optimal parameters (i.e., an optimal estimation θ^∈ℝn^+m^+1) of the model specified by IdNon^,m^.


### 3.1. PositionxiPARk and Velocity viPARk of PiPAR

The position ofPiPAR is denoted by xiPARk as
(14)xiPARk≔xIdNoDmin,0,iTk,xIdNoDmin,1,iTk,…,xIdNoDmin,Dmin−1,iTkT,…,xIdNon,0,iTk,xIdNon,1,iTk,…,xIdNon,n−1,iTkT,…,xIdNoDmax,0,iTk,xIdNoDmax,1,iTk,…,xIdNoDmax,Dmax−1,iTkT,where xIdNon,m,ik denotes the parameter vector of a model specified by the identifier number IdNon,m as follows:
(15)xIdNon,m,ik≔a^1k,a^2k,…,a^nk,b^0k,b^1k,…,b^mkT∈ℝn+m+1,where a^ℓk and b^ℓk denote, respectively, the estimates of aℓ∗ and bℓ∗ in (2)-(3) at the kth iteration. Next, the velocity of PiPAR is denoted by viPARk, which has the same dimension as that of xiPARk, and the velocity corresponding to xIdNon,m,ik is denoted by vIdNon,m,ik∈ℝn+m+1.


### 3.2. Positionxstr,ik and Velocity vstr,ik of PiSTR

The position ofPiSTR is denoted by xstr,ik∈ℝ, and its velocity is denoted by vstr,ik∈ℝ. Note that, for given Dmin and Dmax, the boundary condition of xstr,ik∈ℝ,∀kis derived as
(16)1=IdNoDmin,0≤xstr,ik≤IdNoDmax,Dmax−1=DmaxDmax+1−DminDmin−12.


### 3.3. Positional Update Law forPiSTR

The positional update ofPiSTR (i.e., the update of model identifier number) is performed for each particle i∈1,npstr, where npstr denotes the total number of particles, as follows:
(17)vstr,ik+1=c0strvstr,ik+c1strr1,istrxstr,iPbestk−xstr,ik+c2strr2,istrxstr,iSbestk−xstr,ik,(18)xstr,ik+1=xstr,ik+Cvelstrvstr,ik+1,vminstr,vmaxstr,(19)xstr,ik+1←Cposstrxstr,ik+1,IdNoDmin,0,IdNoDmax,Dmax−1,where Cvelstr·,· and Cposstr·,· are the clamping operators applied over the velocity component vstr,i· and positional component xstr,i·, respectively. The inertia factor c0str, cognitive scaling factor c1str, and social scaling factor c2str in (17) are specified by the designer. The random numbers r1,istr and r2,istr are uniformly distributed in 0,1 and represent the stochastic behaviors of PSO. Note that the maximum and minimum limits vmaxstr and vminstr of the velocity serve as constraints to avoid explosion. The above two operators are applied in the following two ways:
(20)Cvelstrvstr,ik,vminstr,vmaxstr=vstr,ik,ifvminstr≤vstr,ik≤vmaxstr,vminstr,elseifvstr,ik<vminstr,vmaxstr,elsevstr,ik>vmaxstr,Cposstrxstr,ik,IdNoDmin,0,IdNoDmax,Dmax−1=xstr,ik,ifIdNoDmin,0≤xstr,ik≤IdNoDmax,Dmax−1,IdNoDmin,0,elseifxstr,ik<IdNoDmin,0,IdNoDmax,Dmax−1,elsexstr,ik>IdNoDmax,Dmax−1,where the notation · denotes the floor operator that is introduced because xstr,ik must be an integer.In (17), xstr,iPbestk∈ℝ denotes the personal best position of PiSTR (i.e., the best model identifier number that the individual particle PiSTR has achieved so far), where “best” means that the transfer function model with the parameter vector xIdNon,m,ik in (15) specified by the model identifier number IdNon,m=xstr,ik minimizes the given fitness function and is defined as follows:
(21)nPbest,ik,mPbest,ik,ℓPbest,ik=argminIdNon,m←xstr,iℓ,1≤ℓ≤kLxIdNon,m,iℓ,(22)xstr,iPbestk←nPbest,iknPbest,ik−12+mPbest,ik+1−DminDmin−12,where LxIdNon,m,i· denotes one of the criterion functions previously defined as fOL−timexIdNon,m,i· in (7), fCL−timexIdNon,m,i· in (9), fOL−freqxIdNon,m,i· in (12), and fCL−freqxIdNon,m,i· in (13). An illustrative representation of the procedure for the determination of xstr,iPbestk is shown in Figure 1 and is explained as follows: (i) identify IdNon,m that is identical to xstr,iℓ, where 0≤ℓ≤k, (ii) calculate the objective function LxIdNon,m,iℓ∣IdNon,m←xstr,iℓ, (iii) evaluate the objective function values to find the pair of n,m that achieves the minimum of the objective function value, and (iv) set xstr,iPbestk as (22).Figure 1
Illustrative representation of the procedure for determiningxstr,iPbestk.On the other hand,xstr,iSbestk is determined as
(23)nSbest,ik,mSbest,ik,ℓSbest,ik,jSbest,ik=argminj∈i−nsstr/2,…,i,…,i+nsstr/2minDmin≤n≤Dmax,0≤m≤n−1,1≤ℓ≤kLxIdNon,m,jℓ,(24)xstr,iSbestk←nSbest,iknSbest,ik−12+mSbest,ik+1−DminDmin−12,where even-numbered nsstr≤npstr is the number of neighbors of PiSTR and xIdNon,m,jℓ≔xIdNon,m,j−1modns,str+1ℓ for j<1 or nsstr+1≤j [24]. An illustrative representation of the procedure for determining xstr,iSbestk, where nsstr=2 (i.e., Pi−1STR and Pi+1STR are the neighboring particles of PiSTR), is shown in Figure 2 and explained as follows: (i) calculate the objective function values of all entries of xiPARk in (14), (ii) find the minimum entry of LxiPARk (e.g., LxIdNon2,m2,ik in Figure 2), (iii) repeat the above processes (i) and (ii) for the neighboring particles Pi−1STR and Pi+1STR (e.g., LxIdNon1,m1,i−1k and LxIdNon3,m3,i+1k are, respectively, the minimum entries of xi−1PARk and xi+1PARk in Figure 2), (iv) evaluate three minimum objective function values (i.e., LxIdNon2,m2,ik, LxIdNon1,m1,i−1k, and LxIdNon3,m3,i+1k) and LxIdNon¯,m¯,j¯ℓ¯ with n¯,m¯,j¯,ℓ¯←nSbest,ik−1,mSbest,ik−1,jSbest,ik−1,ℓSbest,ik−1 to find a n,m pair that achieves the minimum among the four objective function values, and (v) set xstr,iSbestk as (24).Figure 2
Illustrative representation of the procedure for determiningxstr,iSbestk.


### 3.4. Positional Update Law forPiPAR

Once the positionxstr,ik+1 of PiSTR is updated as explained in the above section, the positional update of PiPAR (i.e., the update of xiPARk in (14)) is performed for each particle, i∈1,nppar, as follows. First, the model identifier number xstr,ik+1 in (19) designates one of the entries in xiPARk in (14) to be updated via IdNon,m←xstr,ik+1. The parameter vector, xIdNon,m,ik, of the model indexed by IdNon,m is then updated as follows:
(25)vIdNon,m,ik+1=c0parvIdNon,m,ik+c1parr1,iparxIdNon,m,iPbestk−xIdNon,m,ik+c2parr2,iparxIdNon,m,iSbestk−xIdNon,m,ik,(26)xIdNon,m,ik+1=xIdNon,m,ik+vIdNon,m,ik+1,where three factors c0par, c1par, and c2par are specified by the designer, and r1,ipar and r2,ipar denote random numbers. In (25), xIdNon,m,iPbestk and xIdNon,m,iSbestk are updated at each iteration as
(27)xIdNon,m,iPbestk=argmin1≤ℓ≤kLxIdNon,m,iℓ,xIdNon,m,iSbestk=argargminj∈i−nspar/2,…,i,…,i+nspar/2min1≤ℓ≤kLxIdNon,m,jℓ.Note that forIdNo·,·≠xstr,ik+1, vIdNo·,·,ik+1←vIdNo·,·,ik and xIdNo·,·,ik+1←xIdNo·,·,ik.The above two consecutive updates forPiSTR and PiPAR are iteratively performed until the stopping criterion is reached. Once the optimization process is terminated, the optimum of xstr,i· provides the information IdNon,m about the optimal orders of numerator and denominator polynomials in the transfer function, and the optimum of xIdNon,m,i· provides the optimal values of such numerator/denominator polynomials’ unknown parameters.Remark 1.
The fitness of the candidate model structure and its parameters explored by particles is evaluated based on the criterion function defined in Section2. For the fitness evaluation, it necessarily introduces a technique that can handle the constraint condition such as (8) or (11), which are generally treated as hard constraints in system identification problems. It is therefore crucial to account for the given constraint conditions in any form within the PSO framework in order to calculate the fitness value of each particle and, ultimately, find the optimal model structure/parameters. The constraint-handling strategy utilized in this study is adopted from Chun et al. [25]. This technique enables the original criterion function subject to constraints to be converted into a pseudofunction that ensures all particles move toward a feasible region of the search space at the beginning stage of evolution and then explore the feasible region to find an optimal solution. Note that, because of the flexibility of PSO, such a redefinition of the criterion function can be used without any problem arising, regardless of the objective and constraint functions. Therefore, if the above technique is applied to the system identification problems considered in this study, a global optimal solution (or, at least, a feasible solution) that guarantees the given constraints can be efficiently found, which is experimentally demonstrated in the following sections.Remark 2.
Comparing to the conventional PSO-based white- or gray-box modeling techniques, the additional computational complexity mainly results from the optimization process ofxstr,ik∈ℝ, since the dimension of xIdNon,m,ik varies at each iteration according to xstr,ik updated at the same iteration. However, because xstr,ik is only a one-dimensional parameter, the associated computational burden is not too high and further the convergence speed of xstr,ik can be improved when Dmin and Dmax are reasonably adopted. Therefore, the computational complexity of the proposed method is comparable to that of the existing PSO-based system identification techniques.


## 3.1. PositionxiPARk and Velocity viPARk of PiPAR

The position ofPiPAR is denoted by xiPARk as
(14)xiPARk≔xIdNoDmin,0,iTk,xIdNoDmin,1,iTk,…,xIdNoDmin,Dmin−1,iTkT,…,xIdNon,0,iTk,xIdNon,1,iTk,…,xIdNon,n−1,iTkT,…,xIdNoDmax,0,iTk,xIdNoDmax,1,iTk,…,xIdNoDmax,Dmax−1,iTkT,where xIdNon,m,ik denotes the parameter vector of a model specified by the identifier number IdNon,m as follows:
(15)xIdNon,m,ik≔a^1k,a^2k,…,a^nk,b^0k,b^1k,…,b^mkT∈ℝn+m+1,where a^ℓk and b^ℓk denote, respectively, the estimates of aℓ∗ and bℓ∗ in (2)-(3) at the kth iteration. Next, the velocity of PiPAR is denoted by viPARk, which has the same dimension as that of xiPARk, and the velocity corresponding to xIdNon,m,ik is denoted by vIdNon,m,ik∈ℝn+m+1.


## 3.2. Positionxstr,ik and Velocity vstr,ik of PiSTR

The position ofPiSTR is denoted by xstr,ik∈ℝ, and its velocity is denoted by vstr,ik∈ℝ. Note that, for given Dmin and Dmax, the boundary condition of xstr,ik∈ℝ,∀kis derived as
(16)1=IdNoDmin,0≤xstr,ik≤IdNoDmax,Dmax−1=DmaxDmax+1−DminDmin−12.


## 3.3. Positional Update Law forPiSTR

The positional update ofPiSTR (i.e., the update of model identifier number) is performed for each particle i∈1,npstr, where npstr denotes the total number of particles, as follows:
(17)vstr,ik+1=c0strvstr,ik+c1strr1,istrxstr,iPbestk−xstr,ik+c2strr2,istrxstr,iSbestk−xstr,ik,(18)xstr,ik+1=xstr,ik+Cvelstrvstr,ik+1,vminstr,vmaxstr,(19)xstr,ik+1←Cposstrxstr,ik+1,IdNoDmin,0,IdNoDmax,Dmax−1,where Cvelstr·,· and Cposstr·,· are the clamping operators applied over the velocity component vstr,i· and positional component xstr,i·, respectively. The inertia factor c0str, cognitive scaling factor c1str, and social scaling factor c2str in (17) are specified by the designer. The random numbers r1,istr and r2,istr are uniformly distributed in 0,1 and represent the stochastic behaviors of PSO. Note that the maximum and minimum limits vmaxstr and vminstr of the velocity serve as constraints to avoid explosion. The above two operators are applied in the following two ways:
(20)Cvelstrvstr,ik,vminstr,vmaxstr=vstr,ik,ifvminstr≤vstr,ik≤vmaxstr,vminstr,elseifvstr,ik<vminstr,vmaxstr,elsevstr,ik>vmaxstr,Cposstrxstr,ik,IdNoDmin,0,IdNoDmax,Dmax−1=xstr,ik,ifIdNoDmin,0≤xstr,ik≤IdNoDmax,Dmax−1,IdNoDmin,0,elseifxstr,ik<IdNoDmin,0,IdNoDmax,Dmax−1,elsexstr,ik>IdNoDmax,Dmax−1,where the notation · denotes the floor operator that is introduced because xstr,ik must be an integer.In (17), xstr,iPbestk∈ℝ denotes the personal best position of PiSTR (i.e., the best model identifier number that the individual particle PiSTR has achieved so far), where “best” means that the transfer function model with the parameter vector xIdNon,m,ik in (15) specified by the model identifier number IdNon,m=xstr,ik minimizes the given fitness function and is defined as follows:
(21)nPbest,ik,mPbest,ik,ℓPbest,ik=argminIdNon,m←xstr,iℓ,1≤ℓ≤kLxIdNon,m,iℓ,(22)xstr,iPbestk←nPbest,iknPbest,ik−12+mPbest,ik+1−DminDmin−12,where LxIdNon,m,i· denotes one of the criterion functions previously defined as fOL−timexIdNon,m,i· in (7), fCL−timexIdNon,m,i· in (9), fOL−freqxIdNon,m,i· in (12), and fCL−freqxIdNon,m,i· in (13). An illustrative representation of the procedure for the determination of xstr,iPbestk is shown in Figure 1 and is explained as follows: (i) identify IdNon,m that is identical to xstr,iℓ, where 0≤ℓ≤k, (ii) calculate the objective function LxIdNon,m,iℓ∣IdNon,m←xstr,iℓ, (iii) evaluate the objective function values to find the pair of n,m that achieves the minimum of the objective function value, and (iv) set xstr,iPbestk as (22).Figure 1
Illustrative representation of the procedure for determiningxstr,iPbestk.On the other hand,xstr,iSbestk is determined as
(23)nSbest,ik,mSbest,ik,ℓSbest,ik,jSbest,ik=argminj∈i−nsstr/2,…,i,…,i+nsstr/2minDmin≤n≤Dmax,0≤m≤n−1,1≤ℓ≤kLxIdNon,m,jℓ,(24)xstr,iSbestk←nSbest,iknSbest,ik−12+mSbest,ik+1−DminDmin−12,where even-numbered nsstr≤npstr is the number of neighbors of PiSTR and xIdNon,m,jℓ≔xIdNon,m,j−1modns,str+1ℓ for j<1 or nsstr+1≤j [24]. An illustrative representation of the procedure for determining xstr,iSbestk, where nsstr=2 (i.e., Pi−1STR and Pi+1STR are the neighboring particles of PiSTR), is shown in Figure 2 and explained as follows: (i) calculate the objective function values of all entries of xiPARk in (14), (ii) find the minimum entry of LxiPARk (e.g., LxIdNon2,m2,ik in Figure 2), (iii) repeat the above processes (i) and (ii) for the neighboring particles Pi−1STR and Pi+1STR (e.g., LxIdNon1,m1,i−1k and LxIdNon3,m3,i+1k are, respectively, the minimum entries of xi−1PARk and xi+1PARk in Figure 2), (iv) evaluate three minimum objective function values (i.e., LxIdNon2,m2,ik, LxIdNon1,m1,i−1k, and LxIdNon3,m3,i+1k) and LxIdNon¯,m¯,j¯ℓ¯ with n¯,m¯,j¯,ℓ¯←nSbest,ik−1,mSbest,ik−1,jSbest,ik−1,ℓSbest,ik−1 to find a n,m pair that achieves the minimum among the four objective function values, and (v) set xstr,iSbestk as (24).Figure 2
Illustrative representation of the procedure for determiningxstr,iSbestk.


## 3.4. Positional Update Law forPiPAR

Once the positionxstr,ik+1 of PiSTR is updated as explained in the above section, the positional update of PiPAR (i.e., the update of xiPARk in (14)) is performed for each particle, i∈1,nppar, as follows. First, the model identifier number xstr,ik+1 in (19) designates one of the entries in xiPARk in (14) to be updated via IdNon,m←xstr,ik+1. The parameter vector, xIdNon,m,ik, of the model indexed by IdNon,m is then updated as follows:
(25)vIdNon,m,ik+1=c0parvIdNon,m,ik+c1parr1,iparxIdNon,m,iPbestk−xIdNon,m,ik+c2parr2,iparxIdNon,m,iSbestk−xIdNon,m,ik,(26)xIdNon,m,ik+1=xIdNon,m,ik+vIdNon,m,ik+1,where three factors c0par, c1par, and c2par are specified by the designer, and r1,ipar and r2,ipar denote random numbers. In (25), xIdNon,m,iPbestk and xIdNon,m,iSbestk are updated at each iteration as
(27)xIdNon,m,iPbestk=argmin1≤ℓ≤kLxIdNon,m,iℓ,xIdNon,m,iSbestk=argargminj∈i−nspar/2,…,i,…,i+nspar/2min1≤ℓ≤kLxIdNon,m,jℓ.Note that forIdNo·,·≠xstr,ik+1, vIdNo·,·,ik+1←vIdNo·,·,ik and xIdNo·,·,ik+1←xIdNo·,·,ik.The above two consecutive updates forPiSTR and PiPAR are iteratively performed until the stopping criterion is reached. Once the optimization process is terminated, the optimum of xstr,i· provides the information IdNon,m about the optimal orders of numerator and denominator polynomials in the transfer function, and the optimum of xIdNon,m,i· provides the optimal values of such numerator/denominator polynomials’ unknown parameters.Remark 1.
The fitness of the candidate model structure and its parameters explored by particles is evaluated based on the criterion function defined in Section2. For the fitness evaluation, it necessarily introduces a technique that can handle the constraint condition such as (8) or (11), which are generally treated as hard constraints in system identification problems. It is therefore crucial to account for the given constraint conditions in any form within the PSO framework in order to calculate the fitness value of each particle and, ultimately, find the optimal model structure/parameters. The constraint-handling strategy utilized in this study is adopted from Chun et al. [25]. This technique enables the original criterion function subject to constraints to be converted into a pseudofunction that ensures all particles move toward a feasible region of the search space at the beginning stage of evolution and then explore the feasible region to find an optimal solution. Note that, because of the flexibility of PSO, such a redefinition of the criterion function can be used without any problem arising, regardless of the objective and constraint functions. Therefore, if the above technique is applied to the system identification problems considered in this study, a global optimal solution (or, at least, a feasible solution) that guarantees the given constraints can be efficiently found, which is experimentally demonstrated in the following sections.Remark 2.
Comparing to the conventional PSO-based white- or gray-box modeling techniques, the additional computational complexity mainly results from the optimization process ofxstr,ik∈ℝ, since the dimension of xIdNon,m,ik varies at each iteration according to xstr,ik updated at the same iteration. However, because xstr,ik is only a one-dimensional parameter, the associated computational burden is not too high and further the convergence speed of xstr,ik can be improved when Dmin and Dmax are reasonably adopted. Therefore, the computational complexity of the proposed method is comparable to that of the existing PSO-based system identification techniques.


## 4. Numerical Simulation and Discussion

To appraise the performance of the proposed constrained MD-PSO approach, a typical 4th-order nonminimal phase transfer function model appearing often in the literature (e.g., [26–29]) is considered.
(28)G∗s=B∗sA∗s=−6400s+1600s4+5s3+408s2+416s+1600.The above continuous-time model is simulated using a sampling timeTs=0.001s and a pseudorandom binary sequence (PRBS) input shown in Figure 3(a). It is assumed that only the contaminated outputs y~tk≔ytk+wtk (tk=0,Ts,2Ts,…,5000Ts), where wts is a stochastic noise sequence with zero mean, are measurable. In what follows, we adopt the measurable noises so that the noise-to-signal ratio (NSR) becomes around 10% for Figure 3(b) and 50% for Figure 3(c). The upper and lower bounds of order of the denominator polynomial n^ are set to 1=Dmin≤n^≤10=Dmax. This implies that 55=IdNoDmax,Dmax−1−IdNoDmin,0+1 candidate models are evaluated for optimal model structure identification. The values of various parameters for our PSO algorithm are shown in Table 2. The model structure and its associated unknown parameters are identified by minimizing the criterion function (7) and (8) in Section 2.1, using our constrained MD-PSO algorithm in Section 3. For the case of NSR = 10%, the optimal model identifier number in (6) found via the algorithm in Section 3.3 is IdNon^,m^=13 with n^=4 and m^=1. Figure 4(a) shows the convergence characteristics of IdNon^,m^ for kmaxstr=2000 iterations of 30 independent runs. Then, the identified optimal model is as follows:
(29)G^sNSR=10%=−6430.441s+1611.574s4+5.003s3+407.941s2+416.009s+1593.404.Figure 3
PRBS input signal (a) and measured behavior ofy~tk when NSR = 10% (b) and NSR = 50% (c).
(a)
(b)
(c)Table 2
Parameters for PSO algorithms.
Parameters
Simulation studies
Experimental studies
Ball levitation
Flexible joint manipulator
c0str
1
1
1
c1str, c2str
1.5389
1.5389
1.5389
vminstr
−5
−5
−15
vmaxstr
5
5
15
Dmin
1
1
1
Dmax
10
10
10
c0par
0.7543
0.7543
0.7543
c1par, c2par
1.5389
1.5389
1.5389
npstr, nppar
200
200
1000
nsstr, nspar
20
20
5
kmaxstr, kmaxpar
2000
4000
6000Figure 4
Optimization results: convergence properties for NSR = 10% case.
(a)
(b)The convergence characteristics of the criterion function values for all runs are shown in Figure4(b). The comparative validation performance is shown graphically in Figure 5. Figure 5(a) shows a simulation run with the estimated model G^s in (29) and comparison with measurement y~tk. Note that a MATLAB toolbox for system identification offers some tools to assist in the task of model order selection. However, the model structure and parameter identification result in (30) leads to a poor identification performance, which can be verified from Figure 5(a).
(30)G^sMATLAB,NSR=10%=−15.54s+2.178s2+0.8901s+4.261.Figure 5
Validation of the model identification results for NSR = 10% case.
(a)
(b)To show the real performance, the true errorytk−G^p∣NSR=10%utk(rather than the measured error y~tk−G^p∣NSR=10%utk) is illustrated in Figure 5(b). The plotted error series reveal the advantage of the estimated model G^s, that is, our identification approach is more robust against unknown and unpredictable measurement noises than the method using a MATLAB toolbox for system identification. Note that the Hankel norm is generally known to be a sensible measure of the I/O characteristics of the system. The Hankel norms of error systems, G∗s−G^s∣NSR=10%H and G∗s−G^s∣MATLAB,NSR=10%H, are, respectively, 0.0867 and 2.0583. The calculated Hankel norms and Figure 5 show that the model (29) developed via the proposed modeling methodology shows a very good model fidelity.On the other hand, our approach and the MATLAB toolbox-based method are applied for developing the optimal model when NSR =50% (Figure 3(c)), and the following results are obtained:
(31)G^sNSR=50%=1.763s2−6422.455s+1598.403s4+5.007s3+408.034s2+420.246s+1616.024,(32)G^sMATLAB,NSR=50%=−17.13s+3.432.The behaviors of the outputs ofG^s∣NSR=50% and G^s∣MATLAB,NSR=50% are shown in Figure 6(a), where the simulated output of G^s∣NSR=50% matches the measured data y~tk very well even in the presence of heavy measurement noise. This fact can also be confirmed from Figure 6(b). Note that the order of the numerator polynomial of the estimated model G^s∣NSR=50% is different from that of the target system G∗s. However, the Hankel norm of the error system, G∗s−G^s∣NSR=50%H, is 0.2157. The above facts again confirm that the estimated model G^s∣NSR=50% is quite accurate.Figure 6
Validation of the model identification results for NSR=50% case.
(a)
(b)


## 5. Experimental Validations of Constrained MD-PSO-Based System Identification Method

The first experiment was conducted using the magnetic levitation (maglev) system (model number 33-210 from Feedback Instruments Ltd.) shown in Figure7 to validate the performance of the proposed constrained MD-PSO-based system identification method. Note that this maglev setup serves as a simple model of practical devices whose applications can be found in maglev trains and magnetic bearings and which have achieved immense popularity recently [30]. The overall maglev system consists of connection interface panels along with a mechanical-electrical part. A coil is mounted on the mechanical part, and an infrared sensor is also attached to this part. The sensor measures the ball’s position which is converted to a corresponding sensor voltage. The analog-to-digital (A/D) interface converts the measured voltage to digital signal that is fed to the controller. The control signal produced by the implemented controller is then fed to digital-to-analog (D/A) interface which generates analog control signal related to an equivalent coil current. The flow of current in the coil finally produces the force generated by the magnetic field that counterbalances the gravitational pull of the ball [30]. The detailed experimental setup has also been described in detail in Debdoot et al. [31] and Pati et al. [32]. For the maglev system, the open-loop identification method cannot be adopted because of its inherently unstable nature. Therefore, the model of the maglev system was identified in closed-loop operation after the experimental system was stabilized about the operating point with the PID controller Cs=−0.1s2−2s−2/s, where the gains were manually tuned. Furthermore, a dynamic model was developed by applying the frequency-domain identification method in this experimental study. To achieve this aim, a set of frequency-domain measurements was obtained by injecting a phase-shifted multisine signal with frequency components of 0.01−30 (rad/s) and a sampling frequency of 0.01 (rad/s) into the closed-loop maglev system. The maglev system is manufactured for sampling time Ts=0.001s. The experimental frequency response shown in Figure 8 was obtained via discrete Fourier transform analysis of the measured time-domain response of the closed-loop maglev system.Figure 7
Magnetic levitation experimental setup.Figure 8
Measurement of the maglev system and simulated frequency responses of the identified models.In Debdoot et al. [31], Swain et al. [30], and Pati et al. [32], the linearized model structure in terms of ball position and electromagnetic coil current, which was derived based on theoretical physical knowledge of the maglev system, was derived as follows:
(33)G^conv1s=−3518.85s2−2180.On the other hand, Vinodh Kumar and Jerome [33] introduced the following form of a linearized dynamic model for a maglev system, whose nonlinear mathematical model was derived from the first principle as follows:
(34)G^conv2s=b0s3+a1s2+a2s+a3.After fixing the model structure as (34) (i.e., Dmin=Dmax=n^=3 and m^=0), only the parameter identification methodology of Section 3 was applied for estimating uncertain θ^n^,m^≔a1,a2,a3,b0∈ℝ4. The criterion function given in (13) was introduced for the frequency-domain closed-loop identification of the unstable maglev system. The parameter estimation results were as follows: a1=98.45, a2=−552.9, a3=−9.279×104, and b0=−1.348×105. However, the poor agreement of the simulated frequency response of G^conv2s to the experimental measurement data can be confirmed from Figure 8. Thus, to obtain a more accurate model, our MD-PSO-based identification with 1≤n^≤10 and various PSO parameters in Table 2 was performed. The identified optimal model structure and parameters were as follows:
(35)G^maglevs=−1.001×105s−4.465×106s4+90.255s3+5430.614s2−2.215×105s−2.703×106.Figure8 compares the frequency responses produced by the two models, G^conv2s in (34) and G^maglevs in (35), with that experimentally measured using the actual maglev system. This figure shows that the response of our fourth-order model G^maglevs closely matches the experimental data with an error level smaller than that in the G^conv2s case. Note that the simple second-order model G^conv1s in (33) cannot reflect the important properties of the frequency response of the maglev system.On the other hand, the above fact may seem natural, since increasing the order of the identified model can reduce the modeling error due to unmodeled dynamics at high frequencies. This means that if the modeling error is not ignorable, one should decide whether to increase the order of the model, by considering the working frequency span of the target system [34]. However, the usual modeling methodology using the information about the underlying dynamics (e.g., G^conv2s) provides no knowledge concerning how much the orders of numerator and denominator polynomials should be increased. In such a case, the selection of a suitable model structure can totally rely on experts’ personal experiences, which may not ensure the optimality concerning the model’s reliability and validity. Our identification method, by contrast, provides a straightforward and systematic approach for selecting an optimal model structure and the associated parameters simultaneously to fit arbitrary frequency response data to, without relying on an expert’s knowledge. From the above viewpoint, the proposed identification scheme can play a powerful role as a simple and practical tool for developing a model that optimally captures some key dynamic characteristics of the target system.The second experiment was conducted using a robotic manipulator with joint flexibility (Quanser’s rotary single-link flexible joint module). The experimental setup, shown in Figure9(a), comprises a Quanser QPIDe data acquisition and control board, a Quanser VoltPAQ-X2 amp, a Quanser SRV02 plant, a Quanser rotary flexible joint module, and a PC equipped with necessary software, including the Quanser QUARC 2.5. This rotary single-link flexible joint module consists of a rigid beam mounted on a flexible joint that rotates via a dc motor, and the joint deflection is measured using the implemented sensors. More details can be found in Quanser’s rotary flexible joint workbook. In this experiment, a PRBS, whose amplitude switches between ±1, was used to excite the system. The input u and experimental output y signals were sampled at a uniform sampling interval of 0.002 s, over a period of 10 s (N=5001 data points). The input in Figure 10(a) is the voltage of the servomotor and the output γ (tip angular position) in Figure 10(b) is the summation of θ (motor rotation angle) and α (deflection angle of the flexible link), as shown in Figure 9(b).Figure 9
Description of the flexible joint manipulator system.
(a)
(b)Figure 10
PRBS input signal (a) and measured output behavior (b) of flexible joint manipulator.
(a)
(b)Ahmad et al. [36], Talole et al. [37], and Kandroodi et al. [35] described the modeling of the above flexible joint manipulator system, where Euler-Lagrange formulation was considered in characterizing the dynamic behavior of the system. Then, the following linear model represented in a state-space form was derived:
(36)ẋ1=x3,ẋ2=x4,ẋ3=Ax2+Bx3+Cu,ẋ4=Dx2+Ex3−Cu,where x=θ,α,θ̇,α̇T and A, B, C, D, and E are given as
(37)A=KstiffJeq,B=−ηmηgktkmKg2+BeqRmJeqRm,C=ηmηgKtKgJeqRm,D=−KstiffJeq+JarmJeqJarm,E=ηmηgktkmKg2+BeqRmJeqRm.In their research, the values of the parameters in (37) were defined as those in Table 3. For the model derived in (36), the simulated output data, y=1100x, corresponding to the PRBS input signal of Figure 10(a) is plotted in Figure 11(a). This figure shows that the simulated output of the model in (36) cannot match the measured output data, which is also verified from the output error plot in Figure 11(b). On the other hand, Quanser’s rotary flexible joint workbook introduced an identical model structure given in (36), but the system parameters were set differently to those of Ahmad et al. [36], Talole et al. [37], and Kandroodi et al. [35], as shown in Table 3. However, such a model of Quanser’s workbook also cannot guarantee a sufficient model quality and accuracy, as shown in Figure 11.Table 3
System parameters of flexible joint manipulator.
Symbol
Quantity
Value
Ahmad et al. [36] and Kandroodi et al. [35]
Quanser’s workbook
Rm
Armature resistance (Ω)
2.6
2.6
Km
Motor back-EMF constant (V rad/s)
0.00767
0.00768
Kt
Motor torque constant (N M/A)
0.00767
0.00768
Jarm
Total arm inertia (kg m2)
0.0035
0.003616
Jeq
Equivalent inertia (kg m2)
0.0026
0.00208
Kg
High gear ratio
14 : 5
14 : 5
Kstiff
Joint stiffness
1.2485
1.4467
Beq
Equivalent viscous damping of motor (N.M.S/rad)
0.004
0.015
ηg
Gearbox efficiency
0.9
0.9
ηm
Motor efficiency
0.69
0.69Figure 11
Measured and simulated outputs of the flexible joint manipulator system.
(a)
(b)To identify a model of acceptable accuracy from the measured output data, the proposed MD-PSO-based identification with1≤n^≤10 and various PSO parameters in Table 2 was performed. The identified optimal model structure and parameters are as follows:
(38)G^manipulators=2.771×107s2+8.665×108s+1.888×108s5+3.777×106s4+5.417×107s3+5.321×108s2+1.842×108s+4.372×104.The improvement achieved in the model quality is illustrated in Figure11(a), which compares the simulated time history response for the model (38) with the measured data. The remarkable model accuracy became even more obvious when the time history errors between the outputs of the three models and actual measurement were compared, as demonstrated in Figure 11(b).Finally, the theoretical validity of our 5th-order model structure in (38) is explained as follows. The state-space model (36) introduced in Ahmad et al. [36], Talole et al. [37], and Kandroodi et al. [35] was converted into the 4th-order transfer function model. To obtain a dynamic model of the underlying system using Lagrangian function, the energy expression, which is the difference between the total kinetic energy (T) and the total potential energy (V), is derived as
(39)L=T−V=12Jeqθ̇2+12Jarmθ̇+α̇2−12Kstiffα2.The Lagrangian equations of motion are as follows:(40)∂∂t∂L∂θ̇−∂L∂θ=Toutput−Beqθ̇,∂∂t∂L∂α̇−∂L∂α=−Blinkα̇,where Toutput is the motor torque that serves as input to the system and Blink is the equivalent viscous damping coefficient of the link. The equation of motion is then obtained from (39) and (40) as
(41)Jeqθ¨+Jarmθ¨+α¨=Toutput−Beqθ̇,Jarmθ¨+α¨+Kstiffα=−Blinkα̇.Next, the electrical equation of the Quanser SRV02 plant is derived using Kirchhoffs voltage law as follows:(42)Vm−RmIm−Lmİm−Kmθ̇=0,where Vm=u is the motor input voltage, Im is the motor current, and Lm is the inductance. The torque Toutput in (40) and (41) on the load from the motor is defined as
(43)Toutput=ηgηmKtKgIm.Finally, the full linear model of the considered flexible joint manipulator system can be formulated from (41), (42), and (43) as follows: for x=θ,α,θ̇,α̇,ImT,
(44)ẋ1=x3,ẋ2=x4,ẋ3=KstiffJeqx2−BeqJeqx3+BlinkJeqx4+ηgηmKtKgJeqx5,ẋ4=−KstiffJeqJarmJeq+Jarmx2+BeqJeqx3−BlinkJeqJarmJeq+Jarmx4−ηgηmKtKgJeqx5,ẋ5=−KmLmx3−RmLmx5+1Lmu,which is obviously transformed to the 5th-order transfer function model. The above fact provides a concrete foundation for the feasibility of our 5th-order model in (38) identified successfully with no aforementioned theoretical knowledge. Note that most studies including Ahmad et al. [36], Talole et al. [37], and Kandroodi et al. [35] assumed that the motor inductance Lm in (42) is much less than the resistance Rm, and so, it can be ignored, which results in the conventional model given in (36) and (37). In addition, the viscous damping coefficient of the link Blink is assumed to be negligible in most studies. However, the experimental results plotted in Figure 11 clearly demonstrate the insufficient numerator- and denominator-order identification of the conventional 4th-order transfer function model for capturing the dynamic behavior of the flexible joint manipulator system.


## 6. Conclusion

In this study, we investigated a novel, easy-to-use meta-heuristic optimization scheme for simultaneous identification of the model structure and the associated parameters for linear systems. The main tool used in this scheme was the PSO algorithm, which was developed by hybridizing two main methodologies: one for negating the limit for fixing the particle’s dimension within the PSO process and another for enhancing the exploration ability of the particles by adopting a cyclic neighborhood topology of the swarm. Their key features are as follows. First, the cyclic network topology-based neighborhood structure of particles can effectively regulate the particles’ exploration abilities and enable them to locate a promising region in the hyperdimensional search space. Such an improved exploration ability of the particles is essential to cope with the typical nonconvex nature of the criterion function with many local minima in the process of system identification. Second, in contrast to the existing PSO variants, our PSO algorithm can consecutively search both (i) the dimensional optimum of particles and (ii) the positional optimum in the hyperdimensional search space, whose dimension is specified by the aforementioned optimal dimension of particles. As a result, the dimensional optimum provides the optimal orders of numerator and denominator polynomials in the transfer function, while the positional optimum provides the optimal values of such numerator/denominator polynomials’ unknown parameters. This novel PSO scheme features superior reliability and high-level practicality and is easily implemented for an identification procedure of simultaneously fitting a model with a specific model structure and estimating unknown model parameters. It was clearly demonstrated through typical numerical examples that the proposed PSO scheme provides novel and powerful impetus with remarkable reliability toward simultaneous identification of the model structure and unknown model parameters. The experiments conducted on the magnetic levitation system and robotic manipulator with joint flexibility validated the effectiveness of the proposed strategy in practical applications. Future work will focus on verifying the versatility of the developed PSO-based scheme in the identification of the general class of systems that may include linear time-varying, linear parameter-varying, and linear systems with time delay.



---

*Source: 2713684-2018-09-24.xml*