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

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.


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, modelbased 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 modelbased 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][3][4][5][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 responsefractional 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 socalled 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][10][11][12][13][14][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][17][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, N N + 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][17][18][19], has yet adequately addressed a complete method based on swarmbased 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 topologybased 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 2 Complexity 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.

Problem Formulation
Assume that a finite-dimensional, linear time-invariant, continuous-time plant to be identified is described in the compact transfer function form y t = G * p u t with p derivative operator (i.e., p α z t = d α z t /dt α ), where y t and u t are, respectively, the noise-free response and the excitation.Let G * s be a rational form in s, defined as where s is the Laplace transform variable, and A * s and B * s are assumed to be coprime and Let θ * denote the parameter vector, which includes the dynamic plant model's unknown parameters stacked column-wise as 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 Ĝ 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 Ĝ 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, A set of candidate models (i.e., model structure) constructed to be used in the identification procedure is described as follows.Consider the nth-order transfer function model with D min ≤ n ≤ D max , where D min and D max are assumed to be given in advance.Since the model order is equal to the order of the denominator polynomial, the nth-order model has a total of n types of numerator polynomials, as shown in Table 1.Let the model identifier number be denoted by where n, 0 ≤ m < n, is the order of the numerator polynomial.
The lowest and highest identifier numbers are, respectively, IdNo D min , 0 and IdNo D max , D max − 1 , and the total number of candidate models to be evaluated equals to IdNo D max , D max − 1 − IdNo D min , 0 + 1.Therefore, the aforementioned model structure estimation problem (ID-1) becomes the problem of finding an optimal IdNo n, m .Then, the nth-order transfer function model with order m of the numerator polynomial is denoted as Ĝ s, θ n, m , where is a vector of unknown parameters in the numerator and denominator polynomials.

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 a stable target plant, the input signal u t , t 0 ≤ t ≤ t N is applied to the plant G * p , which gives rise to an output signal y t , t 0 ≤ t ≤ t N .Assuming that N + 1 sampled measurements of the input-output signal u t k and y t k ≔ y t k + w t k , where k = 0, 1, … , N and w t s is a stochastic noise sequence with zero mean, are available.
For a stable plant G * p , the criterion function introduced for the open-loop time-domain estimation of both the model structure Ĝ p and its uncertain parameter vector θ n, m from the sampled measurement data is then formulated as where y OL−plant t ≔ G * p u t + w t is the measured output, y OL−model t ≔ Ĝ p, θ n, m u t 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 inherently unstable 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 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: 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 a stable plant, the criterion function introduced for the open-loop frequency-domain estimation of the model structure/parameters becomes 5 Complexity subject to the constraint condition (8).In (12), G * jω denotes the frequency response function of the plant, Ĝ jω, θ n, m is the estimated frequency response function, U jω and Y jω 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 for frequency-domain closed-loop identification of an unstable plant is formulated as where the constraint condition is identical to (11).

Constrained MD-PSO with a Cyclic Neighborhood Topology for System Identification
Let the notations P STR i and P PAR i be defined, respectively, as follows: (i) P STR i : the ith particle in the swarm used for searching the optimal model structure (i.e., an optimal model identifier number IdNo n, m ) and (ii) P PAR i : 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 IdNo n, m .
where x IdNo n,m ,i k denotes the parameter vector of a model specified by the identifier number IdNo n, m as follows:  are uniformly distributed in 0, 1 and represent the stochastic behaviors of PSO.Note that the maximum and minimum limits v str max and v str min of the velocity serve as constraints to avoid explosion.The above two operators are applied in the following two ways: where the notation • denotes the floor operator that is introduced because x str,i k must be an integer.In ( 17), x Pbest str,i k ∈ ℝ denotes the personal best position of P STR i (i.e., the best model identifier number that the individual particle P STR i has achieved so far), where "best" means that the transfer function model with the parameter vector x IdNo n,m ,i k in (15) specified by the model identifier number IdNo n, m = x str,i k minimizes the given fitness function and is defined as follows:

22
where L x IdNo n,m ,i • denotes one of the criterion functions previously defined as f OL−time x IdNo n,m ,i • in (7), (12), and f CL−f req x IdNo n,m ,i • in (13).An illustrative representation of the procedure for the determination of x Pbest str,i k is shown in Figure 1 and is explained as follows: (i) identify IdNo n, m that is identical to x str,i ℓ , where 0 ≤ ℓ ≤ k, (ii) calculate the objective function L x IdNo n,m ,i ℓ | IdNo n,m ←x str,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 x Pbest str,i k as (22).On the other hand, x Sbest str,i k is determined as

24
where even-numbered n str s ≤n str p is the number of neighbors of P STR i and x IdNo n,m ,j ℓ ≔ x IdNo n,m , j−1 mod n s,str +1 ℓ for j < 1 or n str s + 1 ≤ j [24].An illustrative representation of the procedure for determining x Sbest str,i k , where n str s = 2 (i.e., P STR i−1 and P STR i+1 are the neighboring particles of P STR i ), is shown in Figure 2 and explained as follows: (i) calculate the objective function values of all entries of x PAR i k in ( 14), (ii) find the minimum entry of L x PAR i k (e.g., L x IdNo n 2 ,m 2 ,i k in Figure 2), (iii) repeat the above processes (i) and (ii) for the neighboring particles P STR i−1 and P STR i+1 (e.g., L x IdNo n 1 ,m 1 ,i−1 k and L x IdNo n 3 ,m 3 ,i+1 k are, respectively, the minimum entries of x PAR i−1 k and x PAR i+1 k in Figure 2), (iv) evaluate three minimum objective function values (i.e., L to find a n, m pair that achieves the minimum among the four objective function values, and (v) set x Sbest str,i k as (24).

i
. Once the position x str,i k + 1 of P STR i is updated as explained in the above section, the positional update of P PAR i (i.e., the update of x PAR i k in ( 14)) is performed for each particle, i ∈ 1, n par p , as follows.First, the model identifier number x str,i k + 1 in (19) designates one of the entries in x PAR i k in (14) to be updated via IdNo n, m ← x str,i k + 1 .The parameter vector, x IdNo n,m ,i k , of the model indexed by IdNo n, m is then updated as follows: The above two consecutive updates for P STR i and P PAR i are iteratively performed until the stopping criterion is reached.Once the optimization process is terminated, the optimum of x str,i • provides the information IdNo n, m about the optimal orders of numerator and denominator polynomials in the transfer function, and the optimum of x IdNo n,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 Section 2. 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 IdNo (n1,m1−1), i (l) x IdNo (n1,m1), i (l) x IdNo (n1,m1+1), i (l) x Pbest (k) 8 Complexity 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 of x str,i k ∈ ℝ, since the dimension of x IdNo n,m ,i k varies at each iteration according to x str,i k updated at the same iteration.However, because x str,i k is only a one-dimensional parameter, the associated computational burden is not too high and further the convergence speed of x str,i k can be improved when D min and D max are reasonably adopted.Therefore, the computational complexity of the proposed method is comparable to that of the existing PSO-based system identification techniques.

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][27][28][29]) is considered.
The above continuous-time model is simulated using a sampling time T s = 0 001s and a pseudorandom binary sequence (PRBS) input shown in Figure 3(a).It is assumed that only the contaminated outputs y t k ≔ y t k + w t k (t k = 0, T s , 2T s , … , 5000T s ), where w t s 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 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 IdNo n, m = 13 with n = 4 and m = 1. Figure 4(a) shows the convergence characteristics of IdNo n, m for  The convergence characteristics of the criterion function values for all runs are shown in Figure 4(b).The comparative validation performance is shown graphically in Figure 5. Figure 5(a) shows a simulation run with the estimated model Ĝ s in (29) and comparison with measurement y t k .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).
Ĝ s MATLAB,NSR=10% = −15 54s + 2 178 s 2 + 0 8901s + 4 261 30 To show the real performance, the true error y t k − Ĝ p | NSR=10% u t k (rather than the measured error y t k − Ĝ p | NSR=10% u t k ) is illustrated in Figure 5(b).The plotted error series reveal the advantage of the estimated model Ĝ 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.On the other hand, our approach and the MATLAB toolbox-based method are applied for developing the optimal model when NSR = 50% (Figure 3 The behaviors of the outputs of Ĝ s | NSR=50% and Ĝ s | MATLAB,NSR=50% are shown in Figure 6(a), where the simulated output of Ĝ s | NSR=50% matches the measured data y t k 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 Ĝ s | NSR=50% is different from that of the target system G * s .However, the Hankel norm of the error system, G * s − Ĝ s | NSR=50% H , is 0.2157.The above facts again confirm that the estimated model Ĝ s | NSR=50% is quite accurate.

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 Figure 7 to validate the performance of the proposed constrained MD-PSObased 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 openloop 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  11 Complexity point with the PID controller C s = −0 1s 2 − 2s − 2/s, where the gains were manually tuned.Furthermore, a dynamic model was developed by applying the frequencydomain 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 T s = 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.
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: 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: Ĝconv2 s = b 0 s 3 + a 1 s 2 + a 2 s + a 3 34 After fixing the model structure as (34) (i.e., D min = D max = n = 3 and m = 0), only the parameter identification methodology of Section 3 was applied for estimating uncertain θ n, m ≔ a 1 , a 2 , a 3 , b 0 ∈ ℝ 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: a 1 = 98 45, a 2 = −552 9, a 3 = −9 279 × 10 4 , and b 0 = −1 348 × 10 5 .However, the poor agreement of the simulated frequency response of Ĝconv2 s 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:  (34) and Ĝmaglev s in (35), with that experimentally measured using the actual maglev system.This figure shows that the response of our fourth-order model Ĝmaglev s closely matches the experimental data with an error level smaller than that in the Ĝconv2 s case.Note that the simple second-order model Ĝconv1 s 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., Ĝconv2 s ) 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 singlelink flexible joint module).The experimental setup, shown in Figure 9(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 12 Complexity 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 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 statespace form was derived: α, θ, α T and A, B, C, D, and E are given as 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 = 1 1 0 0 x, 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.
To identify a model of acceptable accuracy from the measured output data, the proposed MD-PSO-based identification with 1 ≤ n ≤ 10 and various PSO parameters in Table 2 was performed.The identified optimal model structure and parameters are as follows: The improvement achieved in the model quality is illustrated in Figure 11(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 4thorder 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 The Lagrangian equations of motion are as follows:

40
where T output is the motor torque that serves as input to the system and B link is the equivalent viscous damping coefficient of the link.The equation of motion is then obtained from (39) and (40) as Next, the electrical equation of the Quanser SRV02 plant is derived using Kirchhoffs voltage law as follows: where V m = u is the motor input voltage, I m is the motor current, and L m is the inductance.The torque T output in ( 40) and (41) on the load from the motor is defined as Finally, the full linear model of the considered flexible joint manipulator system can be formulated from (41), (42), and (43) as follows: for x = θ, α, θ, α, I m T , x 2 = x 4 , B eq J eq x 3 + B link J eq x 4 + η g η m K t K g J eq x 5 , x 4 = − K stif f J eq J arm J eq + J arm x 2 + B eq J eq x 3

−
B link J eq J arm J eq + J arm x 4 − η g η m K t K g J eq x 5 ,

44
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 L m in (42) is much less than the resistance R m , 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 B link 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.

Conclusion
In this study, we investigated a novel, easy-to-use metaheuristic 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 topologybased 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 PSObased 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.

3. 1 .
Position x PAR i k and Velocity v PAR i k of P PAR i .The position of P PAR i is denoted by x PAR i k as

Figure 1 :Figure 2 :
Figure 1: Illustrative representation of the procedure for determining x Pbest str,i k (b) and 50% for Figure 3(c).The upper and lower bounds of order of the denominator polynomial n are set to 1 = D min ≤ n ≤ 10 = D max .This implies that 55 = IdNo D max , D max − 1 − Id No D min , 0 + 1 candidate models are evaluated for optimal model structure identification.The values of various parameters for our PSO algorithm are shown in Table

Figure 5 :Figure 4 :
Figure 5: Validation of the model identification results for NSR = 10% case.

Figure 6 :
Figure 6: Validation of the model identification results for NSR=50% case.

Ĝmaglev s = − 1 35 Figure 8
Figure8compares the frequency responses produced by the two models, Ĝconv2 s in(34) and Ĝmaglev s in(35), with that experimentally measured using the actual maglev system.This figure shows that the response of our fourth-order model Ĝmaglev s closely matches the experimental data with an error level smaller than that in the Ĝconv2 s case.Note that the simple second-order model Ĝconv1 s 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., Ĝconv2 s ) 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 singlelink 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

Figure 8 :
Figure 8: Measurement of the maglev system and simulated frequency responses of the identified models.
(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 :Figure 10 :
Figure 9: Description of the flexible joint manipulator system.

Figure 11 :
Figure 11: Measured and simulated outputs of the flexible joint manipulator system.

3 Complexity Table 1 :
Model structure of the nth-order

9
Here, y CL−model t and y CL−plant t are, respectively, defined as y CL−model t ≔ M p, θ n, m u t and y CL−plant t ≔ M * p u t + w t , where the closed-loop systems M •, • and M * p stabilized by an identical known controller C p are defined as IdNo n,m ,i k is denoted by v IdNo n,m ,i k ∈ ℝ n+m+1 .3.2.Position x str,i k and Velocity v str,i k of P STR is denoted by x str,i k ∈ ℝ, and its velocity is denoted by v str,i k ∈ ℝ.Note that, for given D min and D max , the boundary condition of x str,i k ∈ ℝ, ∀kis derived as str pos x str,i k , IdNo D min , 0 , IdNo D max , D max − 1 = x str,i k , if IdNo D min , 0 ≤ x str,i k ≤ IdNo D max , D max − 1 , IdNo D min , 0 , else if x str,i k < IdNo D min , 0 , IdNo D max , D max − 1 , else x str,i k > IdNo D max , D max − 1 , par 2,i denote random numbers.In (25), x Pbest IdNo n,m ,i k and x Sbest IdNo n,m ,i k are updated at each iteration as 1≤ℓ≤k L x IdNo n,m ,j ℓ 27 Note that for IdNo •, • ≠ x str

Table 2 :
Parameters for PSO algorithms.

Table 3 :
System parameters of flexible joint manipulator.