I Introduction
Piezoelectric actuators (PEAs) are widely used in the applications of precision positioning [1] and microvibration isolation [2, 3], thanks to their excellent characteristics including high resolution, fast response, and light weight. Unfortunately, the inherent hysteresis nonlinearity of the PEAs often significantly limits their operation accuracy.
To alleviate the hysteresis effects of the PEAs, many compensation methods have been developed [4, 5, 6, 7, 8]. Among the existing methods, feedforward control is commonly used because of its lowcost and no external sensors required [9]. Generally speaking, its principle consists in establishing the hysteresis model as accurately as possible and then employing the inverse of the model as a compensator to cascade with the controlled system. Motivated by this, there are numerous works in the literature concentrating on hysteresis modeling, such as the Preisach model [10, 11], the PrandtlIshlinskii model [12, 13], and the BoucWen (BW) model [14]. Compared with the first two models, the BW model has the advantage of computational simplicity, as it only requires one differential equation with a few parameters to describe the hysteresis behavior [15]. However, the classical BW (CBW) model is only efficient for the symmetric and rateindependent hysteresis description [16, 17]. Consequently, severe modeling errors would occur when the CBW model is utilized to represent the piezoelectric hysteresis nonlinearity which exhibits asymmetric and ratedependent characteristics. In this regard, various efforts have been recently performed to modify the CBW model. For example, Wang [18] and Zhu [19]
introduced a nonodd input function and an asymmetric formula into the CBW model to achieve the asymmetric hysteresis, respectively. For further describing the ratedependent behavior, Li
[16] developed a Hammerstein structure cascading the asymmetric BW model with the linear dynamics, and Zhu [17] proposed a generalized BW model with a frequency factor, which both make it possible to simultaneously characterize the asymmetric and ratedependent hysteresis effect.Although exciting results have been obtained in [16, 17], they are still not completely satisfactory. For instance, both modified BW hysteresis models in [16, 17] do not take into account the problem that there exists the inherent redundancy of parameters in the CBW model [20], which makes it more difficult to identify the model parameters. At the same time, redundant parameters also easily lead to the divergence of the numerical solution of the differential equation in the identification process. Moreover, most of the existing modified BW models are limited to an integerorder differential equation. It is known that fractional calculus extends the order of classical calculus from integer domain to complex domain, whose unique nonlocal memory effect provides an excellent potential for the application of hysteresis modeling [21, 22, 23, 24]. Specifically for the BW model, the fractional calculus could be a good choice for solving the ratedependent hysteresis problem. However, to the best of authors’ knowledge, there have been few works to apply the fractional calculus to the BW model for the hysteresis modeling of PEAs.
In this paper, a novel fractionalorder normalized BW (FONBW) model is proposed and investigated to improve the modeling accuracy of the hysteresis nonlinearity in a PEA system. The main contributions include: 1) the redundancy of parameters in the CBW model is considered and eliminated via the normalization processing, such that the model has fewer parameters; 2) a generalized input function is proposed to describe the hysteresis nonlinearity with the asymmetric behavior; 3) two fractional operators are introduced into the BW model for the first time to characterize the ratedependent hysteresis effect. The hysteresis parameters of the proposed FONBW model are identified by the selfadaptive differential evolution algorithm on a PEA system. Comparative experiments show that the developed model matches better the system responses than the CBW model.
The remainder of this paper is organized as follows. Section II presents the hysteresis modeling of the PEA system, and then the characteristics of the proposed FONBW model are discussed in Section III. Afterwards, the model identification experiments are conducted to verify the effectiveness of the proposed model in Section IV. Finally, Section V concludes this paper.
Ii Hysteresis Modeling of PEA system
Iia Review of the CBW Model
The CBW model was firstly proposed by Bouc and further modified by Wen for modeling the hysteresis in vibrational mechanics [15]. Owing to the capability of describing many categories of hysteresis and the benefit of simplicity in computing, the CBW model has been extensively applied in piezoelectric hysteresis modeling. In this model, the relationship between the output displacement and input voltage of the PEA system can be expressed in the form of
(1)  
(2) 
where is the input voltage to the PEA; represents the hysteresis output displacement, composed of an elastic term and a purely hysteretic term with the parameters , , and ; is an auxiliary hysteresis variable which is the solution of the nonlinear firstorder differential equation (2), and its shape and magnitude are determined by parameters , , and . Through proper choice of these parameter values, a wide range of hysteresis loops can be described. A detailed description of the relationship between the parameters and the hysteresis loop can be found in [25].
Fig. 1 shows the hysteresis curves generated by the CBW model with the different input frequencies, if the parameter values are chosen as and . It can be found that the hysteresis loop is not only symmetric around its center point, but also does not change with the input frequency. But actually, for the piezoelectric materials, they exhibit the asymmetric and ratedependent behaviors [26], which are unavailable for the CBW model. Therefore, to further improve the modeling accuracy, it is necessary to make a modification.
IiB Normalization Processing
Although several works [18, 19, 16, 17] have been devoted to modifying the CBW model to guarantee the asymmetric or ratedependent characteristic, they generally still suffer from the redundancy of parameters, which originally exists in the CBW model [20]. To facilitate the description of parameter redundancy, without loss of generality, taking two different sets of parameters for example: , , , , , , , where is a positive constant, and the initial condition . Thus, the CBW model in Eqs. (1) and (2) can be rewritten in the following two forms:
(3) 
and
(4) 
According to the given parameters, Eq. (4) can be transferred into
(5) 
Letting , Eq. (5) is then expressed as
(6) 
It is obviously found from Eqs. (3) and (6) that both models exactly deliver the same hysteresis loop for any input signal . That means the inputoutput behavior of the CBW model is not determined by a unique set of parameters . For this reason, it is necessary to normalize the CBW model, making the parameters defined in a unique way.
To this end, define the constants , , , , and the parameter variable . Substituting them into Eqs. (1) and (2) yields
(7)  
(8) 
which are the socalled normalized BW (NBW) model [20]. Note that the normalized form of the CBW model is exactly equivalent to its standard form (1) and (2), if the initial condition . Moreover, it also has the advantage of having only five parameters instead of the seven parameters for the standard form, which make it easier to be identified. Based on the normalized form (7) and (8), next, a modified FONBW model will be developed to describe the hysteresis effect in PEAs with asymmetric and ratedependent characteristics.
IiC Proposed FONBW Model
Inspired by the work in [27], to describe the hysteresis in PEAs which exhibits asymmetric characteristic, a generalized thorder polynomial input function with the memoryless and locally Lipschitz continuous properties is utilized by modifying Eq. (7) as follows:
(9) 
with
(10) 
where is the coefficient of the polynomial function . Combining Eqs. (9) and (10) with (8), the asymmetric NBW model is obtained. To validate the effectiveness of the normalization and asymmetric processing, the output of the asymmetric NBW model under the sinusoidal input signal is depicted in Fig. 2 by simulation. The parameters of the asymmetric NBW model are chosen as , and the other parameters are consistent with those of the CBW model in Fig. 1. As a comparison, the outputs of the CBW and NBW models are also shown in Fig. 2. It can be found that, compared with the CBW model, the NBW model has the same hysteresis loop with less parameters, and the asymmetric NBW model enables to describe the asymmetric behavior.
Furthermore, since the fractional calculus [21, 22, 23, 24] is generally recognized as an effective choice to describe the hysteresis loops with ratedependent characteristic, the idea is also introduced here to develop the FONBW model. To this end, the nonlinear hysteresis effect (8) is extended from the integer order to the fractional order as
(11) 
where the operator denotes the fractional derivative with the order . Here, if taking into account the Caputo’s definition ^{1}^{1}1The thorder Caputo fractional derivative of function with respect to time is defined as [28]: , where , , , is the Gamma function. for the fractionalorder derivatives, one can know that the calculation of the hysteresis variable in Eq. (11) not only depends on the current input voltage, but also requires all of its history states due to the existence of integral operation. That is to say, the nonlocal memory effect of the fractional calculus is essentially beneficial for the description of the inherent ratedependent property of the piezoelectric hysteresis. To the best of authors’ knowledge, this fact has not been recognized in any other modified BW models. Thus, the proposed FONBW model is finally synthesized by Eqs. (9)(11).
It is worth mentioning that the purpose of normalization processing is to make the hysteresis shapes of the CBW model determined by a unique parameter set, reducing the identification complexity. Although the number of parameters in the proposed FONBW model is ultimately bigger than that in the CBW model, it is beneficial for the description of asymmetric and ratedependent hysteresis effect. Actually, compared with the traditional modified BW models [16, 17, 18, 19], the proposed model indeed has an advantage in terms of the number of parameters, which will be seen in Table II.
Iii Discussion on Proposed FONBW Model
As stated in the previous section, the proposed FONBW model not only inherits the advantages of the CBW model, but also makes it possible to simultaneously achieve the asymmetric and ratedependent characteristics with the relative simple mathematic expression. In this section, these notable characteristics will be further discussed in combination with the model parameters.
Iiia Asymmetric Characteristic
It has been disclosed that the asymmetric characteristic of the proposed FONBW model results from the selection of polynomial function . In order to intuitively describe the influence of the order of on the asymmetric characteristic, without loss of generality, taking the order and the coefficients in Eq. (10) for example, the corresponding curve is demonstrated in Fig. 3, where the other curves are shown for the purpose of comparison.
As can be seen, when the input voltage is symmetrical about the origin (shown in Fig. 3(a)), nonodd term needs to exist in to produce the asymmetric behavior, while for the positive input voltage (shown in Fig. 3(b)), any second or higherorder polynomial function can result in the asymmetric characteristic. Moreover, the higher order of , the more generalized class of hysteresis shapes, but the more identification complexity. Therefore, the order of can be chosen by weighing the shapes of the hysteresis and number of parameters.
It should be noted that, the polynomial functions and shown in Fig. 3 are exactly the types adopted in the modified BW models [16, 18] to characterize the asymmetric hysteresis, respectively. Nevertheless, the polynomial input function in this work makes the proposed FONBW model accommodate a more generalized class of hysteresis shapes with the asymmetric behavior.
IiiB RateDependent Characteristic
Thanks to the unique nonlocal memory effect of fractional calculus, the proposed FONBW model enables to describe the ratedependent hysteresis in PEAs. Consequently, to more clearly seen the nature of the hysteresis model, Eq. (11) is broken into following four fractional differential equations:
(12)  
(13)  
(14)  
(15) 
It can be found that each equation represents a part of the hysteresis loop, i.e., Eq. (12) for the positive ascending branch, Eq. (13) for the positive descending branch, Eq. (14) for the negative descending branch and Eq. (15) for the negative ascending branch. They can be easily solved with the Simulink software by choosing proper parameter values, and the solutions are depicted in Fig. 4 with different simulation conditions.
Note that the parameters and also affect the hysteresis shape of the FONBW model, but they are not the cause resulting in the ratedependent characteristic. The influences of those three parameters have been discussed from a mathematical point of view in [20]. Thus, only the fractional orders and are considered here. In the whole simulations for Fig. 4, the parameters and are all set to be 1, but are taken into account the following three cases: 1) ; 2) ; 3) . From Fig. 4, it is seen that a large number of hysteresis shapes can be described by varying the parameters and . Moreover, with the input frequency increased, the hysteresis shape also changes regularly, which validates the ratedependent characteristic of the proposed FONBW model.
IiiC Inverse Multiplicative Structure
Besides the aforementioned characteristics, the proposed FONBW model also exhibits an attractive potential in the hysteresis compensation due to its simple inverse multiplicative structure.
Inspired by the idea of utilizing the CBW model to compensate static hysteresis directly without solving the model inversion [14], the inverse multiplicative structure of the proposed FONBW model is obtained from Eqs. (9)(11) by extracting the input voltage as
(16) 
with
(17) 
where is the given desired displacement. From the control point of view, Eq. (16) can be directly employed as a compensator for the hysteresis nonlinearity with the asymmetric and ratedependent behaviors. The corresponding block diagram is depicted in Fig. 5. Since and are already known during the modeling and identification, no additional calculation for the inverse model is required, which finally makes the compensator simple to be implemented.
Iv Model Identification and Verification
In this section, the parameters of the proposed FONBW model are identified on a real PEA system, and its effectiveness is also verified by comparative experiments.
Iva Experimental Setup
The experimental setup is shown in Fig. 6. The adopted PEA (model Pst120/7/20VS12 with maximal voltage of 120 V, CoreTomorrow Co., China) is actuated through a voltage amplifier (model E00.D6 with 6 channels from the CoreTomorrow Co.). Its output displacement is realtime measured by a laser displacement sensor (model LKH022, Keyence Co., Japan). The sensor output voltage signal is passed through a signal conditioner (model LKG5001 from the Keyence Co.), and then acquired by the A/D channel of a data acquisition card (model PCI6229 with 16bit A/D and D/A converters from NI Co.). The voltage control signal is produced by the D/A channel and then amplified 12 times via the voltage amplifier to drive the PEA. Programs are developed with Matlab/Simulink software on a host PC and downloaded through the TCP/IP mode to a target PC. In the process of experiments, the sampling frequency is set to be 10 kHz.
Moreover, the toolbox FOMCOM [29] is used here to implement the fractionalorder terms in the proposed FONBW model, which are approximated by a fifthorder refined Oustaloup filter and the frequency range is set as 0.011000 rad/s. Through a tradeoff between the model accuracy and the identification complexity, the order of polynomial function is chosen as .
IvB Parameter Identification
Due to the strong coupling effect of the hysteresis nonlinearity and dynamics behaviors in the PEA system, it is not easy to effectively identify the parameters of coupled dynamics through conventional two decoupled steps, i.e., the parameters for the linear and hysteresis components are individually identified. Fortunately, in this work, the linear dynamics is not necessary to be identified thanks to the inherent ratedependent property of the FONBW model.
Nevertheless, the identification of hysteresis models is still a challenging task, and many algorithms including the least mean square, genetic algorithm, and particle swarm optimization have been developed to solve this problem
[26]. In this work, a selfadaptive differential evolution (DE) algorithm [18] featuring simple evolutionary scheme and few tunable parameters is adopted as an illustration. The model identification process is carried out offline and summarized in Fig. 7, where the rootmeansquare (RMS) error is employed as the objective function for reflecting the modeling deviation(18) 
where
is the parameter vector to be identified,
is the number of sampling points, is the th measured displacement, and is the th output displacement generated by the FONBW model. Besides, a variableamplitude variablefrequency input voltage signal is also chosen as an illustration for the identification(19) 
FONBW model  CBW model  

0.1811  0.1547  
0.5552  
0.0364  
0.0272  
2.0006  1.0003  
0.9557      
0.6220     
The identified model parameters are shown in Table I. With the obtained model, Fig. 8 illustrates the model output and experimental results. As a contrast, the parameters of the CBW model are also identified in the same way. It should be mentioned that two auxiliary constants are respectively introduced in place of the terms and to make the identification more easier. It can be seen from Fig. 8 that, compared with the CBW model, the predicted result of the FONBW model matches better the experimental response with the RMS error about 1.46% relative to the motion range of the PEA system, which also validates the accuracy of the identified FONBW model.
To further intuitively demonstrate the advantage of fewer parameters in the proposed FONBW model, Table II lists a quantified comparison with existing modified BW models [16, 17, 19, 18]. Note that although the modified BW model in [18] has fewer parameters than the proposed FONBW model, but it only describes the asymmetric hysteresis, which does not meet the requirements of high modeling accuracy in a wide frequency range.
IvC Experimental Verification
Frequency  1 Hz  2 Hz  5 Hz  10 Hz  15 Hz  20 Hz 

FONBW  1.44%  2.23%  2.50%  3.17%  3.63%  3.93% 
CBW  3.35%  4.12%  8.07%  15.16%  21.36%  26.66% 
To evaluate the effectiveness and feasibility of the identified FONBW model, the sinusoidal input signals with different frequency are applied to the PEA system.
Fig. 9 shows the comparisons of hysteresis curves of the experimental results and the identified models under different input frequencies. The corresponding predicted errors are also demonstrated in Fig. 10. For the quantified comparison, Table III lists the RMS errors of the proposed FONBW model and CBW model with respect to the motion range of the PEA system.
It can be observed that the CBW model cannot exactly describe the complicated hysteresis loops with asymmetric and ratedependent behaviors. Large errors exist between the model output and experimental results. Specifically, with the increase of the input frequencies, the RMS errors deteriorate sharply from 3.35% to 26.66%. In contrast, the FONBW model have a better performance on the asymmetric and ratedependent hysteresis description. The corresponding RMS errors increase slowly from 1.44% to 3.93% in a relative wide frequency band between 1 Hz and 20 Hz. Therefore, the proposed FONBW model can well solve the asymmetric and ratedependent behaviors of the piezoelectric hysteresis nonlinearity with fewer parameters, which obviously validates its effectiveness and feasibility.
V Conclusion
In this paper, a new FONBW model is proposed to describe the asymmetric and ratedependent piezoelectric hysteresis nonlinearity for improving the modeling accuracy. Considering the redundancy of parameters in the CBW model, the normalization processing is firstly conducted such that the model parameters are determined in a unique way. Based on the normalized BW model, a generalized input function is then utilized to represent the asymmetric hysteresis behavior. Furthermore, two fractionalorder operators are introduced into the BW model in place of the integerorder differentials, resulting in the ratedependent characteristic. Model parameters are identified by the selfadaptive DE algorithm. Comparative experimental results on a PEA system validate the effectiveness and superiority of the developed model on the asymmetric and ratedependent hysteresis description. Thanks to the attractive advantages, the proposed modeling approach provides a wide range of possibilities for the modelbased control techniques. Future work will focus on the hysteresis compensation by adopting the proposed inversemultiplicativestructurebased feedforward control approach to improve the tracking performance of the PEAs.
References
 [1] G.Y. Gu, L.M. Zhu, C.Y. Su, H. Ding, and S. Fatikow, “Modeling and control of piezoactuated nanopositioning stages: A survey,” IEEE Transactions on Automation Science and Engineering, vol. 13, no. 1, pp. 313–332, 2016.
 [2] S. Kang, H. Wu, X. Yang, Y. Li, and Y. Wang, “Fractionalorder robust model reference adaptive control of piezoactuated active vibration isolation systems using output feedback and multiobjective optimization algorithm,” Journal of Vibration and Control, p. 1077546319875260, 2019.
 [3] X. Yang, H. Wu, B. Chen, S. Kang, and S. Cheng, “Dynamic modeling and decoupled control of a flexible stewart platform for vibration isolation,” Journal of Sound and Vibration, vol. 439, pp. 398–412, 2019.
 [4] Y. Qin, Y. Tian, D. Zhang, B. Shirinzadeh, and S. Fatikow, “A novel direct inverse modeling approach for hysteresis compensation of piezoelectric actuator in feedforward applications,” IEEE/ASME Transactions on Mechatronics, vol. 18, no. 3, pp. 981–989, 2012.

[5]
Y. Li and Q. Xu, “Adaptive sliding mode control with perturbation estimation and pid sliding surface for motion tracking of a piezodriven micromanipulator,”
IEEE Transactions on control systems technology, vol. 18, no. 4, pp. 798–810, 2009.  [6] Q. Xu, “Identification and compensation of piezoelectric hysteresis without modeling hysteresis inverse,” IEEE Transactions on Industrial Electronics, vol. 60, no. 9, pp. 3927–3937, 2012.
 [7] Y. Fan and U.X. Tan, “Design of a feedforwardfeedback controller for a piezoelectricdriven mechanism to achieve highfrequency nonperiodic motion tracking,” IEEE/ASME Transactions on Mechatronics, vol. 24, no. 2, pp. 853–862, 2019.
 [8] Z. Li, J. Shan, and U. Gabbert, “Inverse compensation of hysteresis using krasnoselskiipokrovskii model,” IEEE/ASME Transactions on Mechatronics, vol. 23, no. 2, pp. 966–971, 2018.
 [9] K. K. Leang, Q. Zou, and S. Devasia, “Feedforward control of piezoactuators in atomic force microscope systems,” IEEE Control Systems Magazine, vol. 29, no. 1, pp. 70–82, 2009.
 [10] H. Tang and Y. Li, “Feedforward nonlinear pid control of a novel micromanipulator using preisach hysteresis compensator,” Robotics and ComputerIntegrated Manufacturing, vol. 34, pp. 124–132, 2015.
 [11] G. Song, J. Zhao, X. Zhou, and J. A. De AbreuGarcía, “Tracking control of a piezoceramic actuator with hysteresis compensation using inverse preisach model,” IEEE/ASME transactions on mechatronics, vol. 10, no. 2, pp. 198–209, 2005.
 [12] G.Y. Gu, L.M. Zhu, and C.Y. Su, “Modeling and compensation of asymmetric hysteresis nonlinearity for piezoceramic actuators with a modified prandtl–ishlinskii model,” IEEE Transactions on Industrial Electronics, vol. 61, no. 3, pp. 1583–1595, 2013.
 [13] M. Al Janaideh, M. Rakotondrabe, and O. Aljanaideh, “Further results on hysteresis compensation of smart micropositioning systems with the inverse prandtl–ishlinskii compensator,” IEEE Transactions on Control Systems Technology, vol. 24, no. 2, pp. 428–439, 2016.
 [14] M. Rakotondrabe, “Bouc–wen modeling and inverse multiplicative structure to compensate hysteresis nonlinearity in piezoelectric actuators,” IEEE Transactions on Automation Science and Engineering, vol. 8, no. 2, pp. 428–431, 2010.
 [15] M. Ismail, F. Ikhouane, and J. Rodellar, “The hysteresis boucwen model, a survey,” Archives of Computational Methods in Engineering, vol. 16, no. 2, pp. 161–188, 2009.
 [16] C.X. Li, L.L. Li, G.Y. Gu, and L.M. Zhu, “Modeling of ratedependent hysteresis in piezoelectric actuators using a hammersteinlike structure with a modified boucwen model,” in International Conference on Intelligent Robotics and Applications. Springer, 2016, pp. 672–684.
 [17] W. Zhu and X.T. Rui, “Hysteresis modeling and displacement control of piezoelectric actuators with the frequencydependent behavior using a generalized bouc–wen model,” Precision Engineering, vol. 43, pp. 299–307, 2016.
 [18] G. Wang, G. Chen, and F. Bai, “Modeling and identification of asymmetric bouc–wen hysteresis for piezoelectric actuator via a novel differential evolution algorithm,” Sensors and Actuators A: Physical, vol. 235, pp. 105–118, 2015.
 [19] W. Zhu and D.h. Wang, “Nonsymmetrical bouc–wen model for piezoelectric ceramic actuators,” Sensors and Actuators A: Physical, vol. 181, pp. 51–60, 2012.
 [20] F. Ikhouane, J. E. Hurtado, and J. Rodellar, “Variation of the hysteresis loop with the bouc–wen model parameters,” Nonlinear Dynamics, vol. 48, no. 4, pp. 361–380, 2007.
 [21] Z. Zhu and X. Zhou, “A novel fractional order model for the dynamic hysteresis of piezoelectrically actuated fast tool servo,” Materials, vol. 5, no. 12, pp. 2465–2485, 2012.
 [22] Y. Liu, J. Shan, U. Gabbert, and N. Qi, “Hysteresis and creep modeling and compensation for a piezoelectric actuator using a fractionalorder maxwell resistive capacitor approach,” Smart Materials and Structures, vol. 22, no. 11, p. 115020, 2013.
 [23] Z. Zhu, S. To, Y. Li, W.L. Zhu, and L. Bian, “External force estimation of a piezoactuated compliant mechanism based on a fractional order hysteresis model,” Mechanical Systems and Signal Processing, vol. 110, pp. 296–306, 2018.
 [24] C. Ding, J. Cao, and Y. Chen, “Fractionalorder model and experimental verification for broadband hysteresis in piezoelectric actuators,” Nonlinear Dynamics, pp. 1–11, 2019.
 [25] T. T. Baber and Y.K. Wen, “Stochastic equivalent linearization for hysteretic, degrading, multistory strutctures,” University of Illinois Engineering Experiment Station. College of Engineering, Tech. Rep., 1980.
 [26] V. Hassani, T. Tjahjowidodo, and T. N. Do, “A survey on hysteresis modeling, identification and control,” Mechanical systems and signal processing, vol. 49, no. 12, pp. 209–233, 2014.
 [27] S. Bashash and N. Jalili, “A polynomialbased linear mapping strategy for feedforward compensation of hysteresis in piezoelectric actuators,” Journal of Dynamic Systems, Measurement, and Control, vol. 130, no. 3, p. 031008, 2008.
 [28] C. A. Monje, Y. Chen, B. M. Vinagre, D. Xue, and V. FeliuBatlle, Fractionalorder systems and controls: fundamentals and applications. Springer Science & Business Media, 2010.
 [29] A. Tepljakov, E. Petlenkov, and J. Belikov, “Fomcom: a matlab toolbox for fractionalorder system identification and control,” International Journal of Microelectronics and Computer Science, vol. 2, no. 2, pp. 51–62, 2011.
Comments
There are no comments yet.