Research Article | Open Access
Identification of Multiple-Mode Linear Models Based on Particle Swarm Optimizer with Cyclic Network Mechanism
This paper studies the metaheuristic optimizer-based direct identification of a multiple-mode system consisting of a finite set of linear regression representations of subsystems. To this end, the concept of a multiple-mode linear regression model is first introduced, and its identification issues are established. A method for reducing the identification problem for multiple-mode models to an optimization problem is also described in detail. Then, to overcome the difficulties that arise because the formulated optimization problem is inherently ill-conditioned and nonconvex, the cyclic-network-topology-based constrained particle swarm optimizer (CNT-CPSO) is introduced, and a concrete procedure for the CNT-CPSO-based identification methodology is developed. This scheme requires no prior knowledge of the mode transitions between subsystems and, unlike some conventional methods, can handle a large amount of data without difficulty during the identification process. This is one of the distinguishing features of the proposed method. The paper also considers an extension of the CNT-CPSO-based identification scheme that makes it possible to simultaneously obtain both the optimal parameters of the multiple submodels and a certain decision parameter involved in the mode transition criteria. Finally, an experimental setup using a DC motor system is established to demonstrate the practical usability of the proposed metaheuristic optimizer-based identification scheme for developing a multiple-mode linear regression model.
The derivation of reasonable mathematical models is the most important part of designing and analyzing control systems; thus, many theoretical and applied studies have been devoted to this research subject [1–6]. Since the early 1940s, linear time-invariant models have garnered considerable attention  and are widely used, because they are considered mathematically easy to analyze and characterize, compared to nonlinear or time-varying systems. However, it is obvious that many real-world mechanical systems involve several types of nonlinearities, such as friction  and backlash , which affect the behavior of systems. Further, there may exist another type of nonlinearity that results from changing the operating point of the system. Therefore, advanced modeling schemes that address the practical difficulties faced by engineers in capturing the complex behavior of systems have received significant research interest.
One attractive approach is to introduce switched linear and piecewise affine (PWA) models (see [10, 11] and the references therein). The development of such models is of considerable importance; thus, many modeling schemes have been widely studied over the last decade in the control system literature. A variety of research trends and methodologies for this topic can be found in the studies by Bako  and Ohlsson and Ljung . This increasing interest is mainly due to the universal approximation properties of PWA maps for some classes of nonlinear systems [14, 15]. The essential feature of a switched linear system is that it has multiple modes corresponding to the response resulting from switching among a finite set of linear subsystems, which enables simple representation of a wide range of practical systems (see  and the references therein). However, at the same time, key difficulties in the system identification processes arise from multiple-mode nature of a switched linear system.
Inferring a multiple-mode model from a set of finite input-output measurements is a highly complex process that requires simultaneous estimation of both the mode transitions and the linear subsystems. The underlying identification problem is inherently nonconvex and admits multiple local solutions, so the formulated optimization problem is not well posed . As mentioned in Ohlsson and Ljung , most conventional schemes can be said to be based on local searches, in which inferring an accurate model depends strongly on good initialization of the optimization. Therefore, convex reformulation of the identification problem to indirectly avoid the drawback related to local searching has been extensively studied in recent years (see, e.g., [12, 13, 17]). This approach usually introduces a relaxation of the nonconvex identification problem; however, the equivalence of optimal solutions to the original nonconvex and transformed convex problems is guaranteed only under certain conditions on the data. Note that the usual difficulty of verifying or obtaining those conditions may limit the practicality of such an approach . Conversely, it is worth mentioning that many existing approaches are often based on prior knowledge of the mode transition. Against this backdrop, efforts have been made for many years to develop efficient modeling methods that do not require such impractical prior knowledge and guarantee that a global optimum is found. One of the very few methods available for achieving these modeling objectives is the mixed integer quadratic programming approach proposed by Roll et al. . However, this program may result in an NP-hard problem, which is quite difficult to solve and tends to be computationally very expensive. This approach is therefore practically applicable only to very small problems with relatively few data .
This paper studies the metaheuristic optimization-based identification of a multiple-mode system, which consists of a finite set of linear regression representations of subsystems, from a collection of input-output data. Note that the multiple-mode linear regression model is simple to describe and captures essential properties of multiple-mode models. In addition, its identification problem is identical to parameter estimation of the subsystems included in PWA systems, without any prior knowledge of their mode transition. Therefore, the investigation of the multiple-mode linear regression model can play an essential role in a variety of control engineering problems. This paper first shows how to reduce the identification problem for a multiple-mode model into an optimization problem. Then, the cyclic-network-topology-based constrained particle swarm optimizer (CNT-CPSO) is applied to solve the formulated nonconvex optimization problem, with no prior knowledge regarding the mode transitions between subsystems. Note that compared to the standard PSO, the CNT-CPSO scheme exhibits improved performance when searching for the global optimum , especially when the optimization problem is ill-conditioned. Then, an experimental setup using a DC motor system is established to verify the performance of the proposed CNT-CPSO-based identification for multiple-mode linear regression models. Note that DC motor systems, which are common components in many practical control systems, usually have nonlinear friction, the modeling of which has always been a difficult and challenging problem (see, e.g., [8, 21–23]). The experimental results show that a finite set of linear regression models, identified in a simple manner via the proposed PSO-based scheme, reconstructs the input-output map with acceptable accuracy; this demonstration consequently verifies the powerful and practical applicability of the proposed multiple-mode modeling scheme. The above facts were examined preliminarily in Maruta et al.  by some of the authors of this paper, and the superior performance of the CNT-CPSO-based modeling method compared to that of the conventional PSO was clearly demonstrated. Our first study in Maruta et al.  aimed at exhaustively optimizing all the unknown parameters of the subsystems via CNT-CPSO with ease; however, its optimization results did not directly provide any information about the mode switching mechanism. Note that when an input-output map must be predicted online using the identified multiple-mode model, knowledge of this mechanism is essential in determining which one among the subsystems is active. Therefore, the original modeling approach is extended to provide not only the optimal parameters of multiple subsystems, but also the decision parameters involved in the mode transition criteria. Its effectiveness is examined using a set of finite input-output measurements obtained from the experimental DC motor system mentioned above. It is worth mentioning that, thanks to the flexibility of the PSO, such an extended identification problem can be solved without any difficulty by directly applying the CNT-CPSO algorithm with no modification.
The remainder of this paper is organized as follows. In Section 2, the multiple-mode linear regression model and its identification issue are considered, and then the reduction of the multiple-mode model identification to an optimization problem is presented in detail. Section 3 briefly reviews the CNT-CPSO algorithm and summarizes a concrete procedure for CNT-CPSO-based identification of a multiple-mode linear regression model. Experimental results for a DC motor system are presented and discussed in Section 4. This section also examines an extended identification methodology for simultaneously finding both the optimal parameters of the subsystems and the decision parameters involved in the mode transition criteria suitable for a DC motor system. Finally, concluding remarks are provided in Section 5.
2. Multiple-Mode Linear Regression Model and Its Identification Problem
This section first introduces the concept of a multiple-mode linear regression model and establishes the identification issues. Then, a method for reducing the identification problem for the multiple-mode model into an optimization problem is described.
2.1. Model Description and Problem Statement
The multiple-mode linear regression model formulated ascan be regarded as describing the linear relationship between a measured signal , designated as the output of the considered system, and the regressor vector, . Here, , where is the parameter vector of the th submodel. Assume that a collection of data sets is measurable, and the number of linear submodels is given a priori.
Let denote the estimate of for . Then, our objective is stated as follows. Given the discrete data for of the switched system (1), find the estimate of the set of coefficient vectors . Figure 1 is an illustrative example, in which multiple circles denote discrete pairs of the generated data, . Figures 1(a) and 1(b) show two linear regression models corresponding to and , respectively, which are explained as follows.(i)If , which is the case in Figure 1(a), the model is equivalent to the standard linear regression model. Thus, its identification is just to estimate one approximating line for a given discrete data set.(ii)Figure 1(b) captures the key feature of a multiple-mode linear regression model with . In this case, the identification of the two-mode model that best matches the given discrete data is equivalent to finding the two lines shown in the figure.
(a) Standard linear regression model
(b) Two-mode linear regression model
2.2. Formulation of Optimization Problem for Multiple-Mode Regression Model Identification
Now suppose that a collection of data sets is given. The identification problem stated in Section 2.1 is then reduced to the following form of optimization problem:where the objective function is defined as the squared sum of errors between the measured discrete data points and the closest linear model chosen in the estimated linear model set. Figure 2 illustrates the calculation of the objective function for and , where data points are denoted as circles and estimated linear models are denoted as lines. The figure shows that, in this situation, the objective function value is the sum of the squares of the length of the dashed lines. Therefore, an optimal solution of the problem (2) can be found by searching , the estimation of , which is obtained from a set of optimal coefficient vectors that best matches the given discrete data. Note that this procedure is computationally expensive; thus, a special algorithm may be required.
For the above optimization problem for multiple-mode linear regression model identification, the following theorem presents the condition under which the estimated set of optimal coefficient vectors becomes uniquely identical to the true set of coefficient vectors.
Theorem 1. Assume that a collection of data sets and the total number of modes, , are given a priori. Denote the total number of data samples in the case that corresponds to the th mode (i.e., ) as . Thus, . For and , let be the time instant satisfying . Then, if the given data set satisfiesthe optimal solution of the optimization problem (2) satisfies , which means that the obtained estimate of the coefficient vector set coincides with the true one.
Its proof, which is explained as follows, is self-evident. The optimal solution mentioned above lets the objective function in (2) be zero. The objective function is clearly nonnegative, and condition (3) ensures the uniqueness of the solution. Note that Theorem 1 ensures the validity of the reduction of the considered identification problem to the optimization problem in the form of (2). However, the formulated optimization problem is inherently ill-conditioned and nonconvex; thus, its global optimum cannot be obtained just by applying conventional deterministic optimization tools or a conventional PSO technique, as mentioned in Section 1. Therefore, in the following section, a recently proposed PSO mechanism, the CNT-CPSO [20, 24], is introduced to avoid entrapment in the suboptimal solutions.
3. The CNT-CPSO-Based Identification Scheme for Multiple-Mode Linear Regression Models
In this section, the CNT-CPSO-based direct identification scheme for multiple-mode linear regression models is described in detail. The design parameter vector is defined as . Let denote the position vector of the th particle of the swarm at the th iteration in the PSO; thus, . Here, the following notation is used for concise representation.Then, the proposed identification scheme for the multiple-mode linear regression model consists of the following steps.
Step 0. Set the initial iteration step to and initialize particles with randomly chosen , where . Then, and are set aswhere the even-numbered (≤np) denotes the number of neighbors the th particle has and denotes the objective function, which is defined in this study as follows:Step 1. If the termination criterion is satisfied, the optimization algorithm is terminated with the following optimal solution:
Otherwise, go to Step .
Step 2. Apply the following evolutionary update law to all the particles:where is equal to the zero vector , andThe inertia factor , cognitive scaling factor , and social scaling factor in (8) are given by the designer. The random numbers and are uniformly distributed in and represent the stochastic behaviors of the PSO. Next, set and then determine and as follows:
Next, go to Step .
Following the above identification procedure using the CNT-CPSO tool, the formulated optimization problem (2) can be solved directly without any difficulty, and eventually, our objective of deriving a multiple-mode linear regression model is achieved. It is important to note that the proposed CNT-CPSO-based scheme handles the identification of PWA systems without any prior knowledge of their mode transitions and can handle a large number of data samples without trouble. These are the most notable features of our CNT-CPSO-based identification methodology. The above fact is thoroughly demonstrated via experimental studies in the following section.
4. Experimental Validation of CNT-CPSO-Based Identification: Multiple-Mode Linear Regression Model for a DC Motor System
Experiments were conducted with a measured set of input-output data from the DC motor system shown in Figure 3 to examine the performance of the CNT-CPSO-based identification scheme for the development of a multiple-mode linear regression model. To this end, our methodology described in Section 3 was applied to the target system, and the reconstruction of the input-output map based on the obtained multiple submodels was studied.
4.1. Identification Problem Formulation
The experimental DC motor system was excited with a sinusoidal electric current (Figure 4(a)), and its angular velocity (Figure 4(b)) and derivative (Figure 4(c)) were measured. From these input-output data, a collection of data sets used for the identification process was defined as , where , and .
The multiple-mode linear regression model for is then constructed as follows:where , , and are the mode-dependent coefficients corresponding to the torque, viscous friction, and Coulomb friction, respectively. Let , , and be defined as , , and , respectively. Then, the model equation (11) is equivalently modified aswhich indicates that its identification problem can be handled by the scheme described in Sections 2 and 3. Here, it is assumed that the number of linear submodels is given a priori as (i.e., ). Under the above model formulation, the proposed CNT-CPSO-based identification procedure was applied to the optimization problem in the form of (2). The optimization results and some discussion are presented in detail in the following subsection.
4.2. Experimental Results and Discussion
To find the estimate of , the identification problem mentioned above was optimized using , , , , and (the maximum PSO iteration number). Then, the following set of optimal coefficient vectors was obtained:Three types of model output corresponding to , , and (i.e., , ) are plotted in Figure 5.
From the parameter vectors ( ()) and the corresponding model outputs , one can infer that the estimate of in (12), which takes one of at each time instant, can be derived as the final outcome using the following formula:The time histories of the obtained , , and are presented in Figure 6. It can easily be confirmed from Figure 6 and the plot in Figure 4 that the three modes are associated with the normal rotation state (), reverse rotation state (), and stall state () of the target DC motor. This classification of the set of identified parameter vectors is probably quite natural. The estimate of the angular acceleration, (=), can then be reconstructed on the basis of , , and in Figure 6 and is presented in Figure 7. This figure shows that the estimate obtained via the three-mode linear model better matches the measured discrete data of than that obtained using the ordinary linear regression model. This result verifies the effectiveness of the proposed identification for a multiple-mode linear regression model.
Some remarks on the characteristics observed from the set of identified parameter vectors are as follows. The Coulomb friction , shown in Figure 6(c), acts in the direction opposite to , which agrees with the well-known behavior of a DC motor, in which this friction depends only on the direction of the angular velocity and not on the magnitude of the velocity. However, the behavior of the identified viscous friction coefficient , shown in Figure 6(b), may disagree with the actual phenomenon in which the viscous friction is proportional to the velocity and reaches zero at zero angular velocity. Note that the viscous friction coefficient of (the parameter vector corresponding to the stall state of the target DC motor) takes a larger value, , than (=−2.53347) and (=−2.60725). This large coefficient value causes the viscous friction term, , in (11) to take a large value when the DC motor is nearly in its stall state. This phenomenon surely lowers the identification accuracy, and Figure 7 confirms that a small but nonnegligible behavioral discrepancy between and occurs when the DC motor is in the stall state. Conversely, once the set of optimal coefficient vectors, , is obtained via the identification scheme proposed in Section 3, the time variation should be found by applying (14). This knowledge of the time variation of makes offline reconstruction of the input-output map possible, as shown in Figure 7. However, the procedure described above does not directly provide any information about the mode switching mechanism. Note that if the input-output map must be predicted online using the identified multiple-mode model, knowledge of this mechanism is essential to determine which one among all the submodels is active. To overcome these shortcomings, the following subsection considers an extended identification scheme for the development of a multiple-mode linear regression model.
4.3. Extended Identification Scheme and Its Experimental Validation
This section considers an extension of the CNT-CPSO-based identification for a multiple-mode linear regression model that enables us to find not only the optimal parameters of the submodels, but also a certain decision parameter involved in the mode transition criteria. Its effectiveness is examined using a set of finite input-output measurements obtained from the experimental DC motor system described above.
The experimental results presented in the previous subsection show that mode switching, , is associated with transitions among the normal rotation state, reverse rotation state, and stall state of the target DC motor. This fact definitely implies that the occurrence of mode switching can be identified by observing the online measured angular velocity, , of the DC motor at each time instant . Therefore, the mode transition criteria depending on are formulated.
The multiple-mode linear regression model (11) is rewritten by applying the forward difference formula to the angular acceleration term as follows:where (=0.01 (s)) denotes the sampling time interval, , , and . Then, the set of the parameter vectors of the three subsystems, , and the mode transition criteria are defined as follows:where (>0) denotes the threshold velocity that enables us to determine which one among the three submodels is active; its value should be optimized along with those of the other parameters in . From (16), one can easily see that , , and correspond to the normal rotation state, reverse rotation state, and stall state of the DC motor, respectively.
Under the above problem formulation, our CNT-CPSO-based identification procedure described in Section 3 is applied to find the optimal values of , , (), and . Let the th particle of the swarm be defined as . Then, the objective function can be modified as follows:where , with , switches at each time instant as It is worth mentioning that, because of the flexibility of the PSO, the above extended identification problem can be solved without any difficulty by directly applying the CNT-CPSO algorithm with no modification. The CNT-CPSO was then run with a PSO setup identical to that described in Section 4.2. The experimental results are summarized below.
The identified parameter vector set and threshold velocity are as follows:where and Note that this extended scheme gives the standard for judging which mode is active as an output of the identification procedure in terms of the observed angular velocity . Therefore, the task of finding a suitable parameter vector at each time instant as presented in (14), which should also be performed after the identification procedure, is not required. Note that the formulation in (14) only allows us to find which one among minimizes the output prediction error (i.e., ) at each time instant . Further, this knowledge of the mode switching criteria may make it possible to predict the system output online, unlike the identification scheme in Section 4.2, which makes only offline reconstruction of the input-output map possible. The predicted angular velocity , which is obtained on the basis of the identified three-mode linear regression model with and the mode transition criteria in (19), is compared with the measured angular velocity in Figure 8. The two angular velocities are perfectly identical, and this fact demonstrates the validity of the extended identification methodology. The time variations of , , and are shown in Figure 9. This figure shows that and exhibit similar trends to those of and in Figure 6; however, the viscous friction coefficient shows entirely different behavior from that of in Figure 6(b). In fact, agrees well with the viscous friction phenomena, because becomes zero (i.e., ) when the DC motor is nearly in its stall state. Therefore, the viscous friction term also approaches zero. This superior identification result clearly leads to better prediction of the angular acceleration . Figure 10 shows the predicted obtained using the model parameters given in (20). Compared to the angular acceleration estimated using (13) and shown in Figure 7, the prediction of was remarkably improved, especially when the DC motor was in the stall state. The above observations clearly verify the effectiveness of the extended identification methodology.
In this paper, the CNT-CPSO-based direct identification of a multiple-mode system was studied, and two new strategies were introduced. The first strategy was to introduce a multiple-mode system consisting of a finite set of linear regression representations of subsystems and then to reduce the identification problem for such a multiple-mode model into an optimization problem. The target systems include the subsystems of PWA systems, and the introduced method is applicable regardless of the mode transition mechanism of the target system. The second strategy was to adopt a metaheuristic optimizer, the CNT-CPSO algorithm, which was developed relatively recently by some of the authors of this paper. This tool plays a key role in addressing some complex difficulties arising due to the inherent ill-conditioned and nonconvex nature of the formulated optimization problem. Then, a concrete procedure for applying the CNT-CPSO-based identification methodology to develop a multiple-mode linear regression model was described in detail. This scheme requires no prior knowledge of mode transitions between subsystems and, unlike some conventional methods, can handle a large amount of data without difficulty during the identification process, which is another distinguishing feature of the proposed method. Finally, an experiment was conducted on a DC motor system to evaluate and demonstrate the practical usability of the proposed metaheuristic optimizer-based identification scheme. Using the same experimental setup, an extension of the CNT-CPSO-based identification scheme was also proposed and examined. Comparison of the experimental results verified that this extension makes it possible to simultaneously obtain both the optimal parameters of multiple submodels and the designated decision parameter involved in the introduced mode transition criteria.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2016R1D1A1B03935288) and the Chung-Ang University Excellent Student Scholarship in 2016.
- L. Ljung, System Identification: Theory for the User, Prentice-Hall, 2nd edition, 1999.
- H. Garnier and L. Wang, “Identification of continuous-time models from sampled data,” in Advances in Industrial Control, Springer, Berlin, Germany, 2008.
- R. Isermann and M. Münchhof, Identification of Dynamic Systems: An Introduction with Applications, Springer, Berlin, Germany, 2011.
- Q. Zhang, Q. Wang, and G. Li, “Switched system identification based on the constrained multi-objective optimization problem with application to the servo turntable,” International Journal of Control, Automation and Systems, vol. 14, no. 5, pp. 1153–1159, 2016.
- J. Yao, Y. Huang, and Z. Ji, “Parameter identification of wiener models by multi-innovation algorithms,” International Journal of Control, Automation, and Systems, vol. 11, no. 6, pp. 1170–1176, 2013.
- D.-L. Zhang, Y.-G. Tang, J.-H. Ma, and X.-P. Guan, “Identification of wiener model with discontinuous nonlinearities using differential evolution,” International Journal of Control, Automation and Systems, vol. 11, no. 3, pp. 511–518, 2013.
- R. Tóth, Modeling and Identification of Linear Parameter-Varying Systems, Springer, Berlin, Germany, 2010.
- B. Armstrong-Hélouvry, P. Dupont, and C. C. De Wit, “A survey of models, analysis tools and compensation methods for the control of machines with friction,” Automatica, vol. 30, no. 7, pp. 1083–1138, 1994.
- M. Nordin and P.-O. Gutman, “Controlling mechanical systems with backlash—a survey,” Automatica, vol. 38, no. 10, pp. 1633–1649, 2002.
- S. Paoletti, A. L. Juloski, G. Ferrari-Trecate, and R. Vidal, “Identification of hybrid systems: a tutorial,” European Journal of Control, vol. 13, no. 2-3, pp. 242–260, 2007.
- A. Garulli, S. Paoletti, and A. Vicino, “A survey on switched and piecewise affine system identification,” in Proceedings of the 16th IFAC Symposium on System Identification, pp. 344–355, October 2012.
- L. Bako, “Identification of switched linear systems via sparse optimization,” Automatica, vol. 47, no. 4, pp. 668–677, 2011.
- H. Ohlsson and L. Ljung, “Identification of switched linear regression models using sum-of-norms regularization,” Automatica, vol. 49, no. 4, pp. 1045–1050, 2013.
- L. Breiman, “Hinging hyperplanes for regression, classification, and function approximation,” IEEE Transactions on Information Theory, vol. 39, no. 3, pp. 999–1013, 1993.
- J.-N. Lin and R. Unbehauen, “Canonical piecewise-linear approximations,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 39, no. 8, pp. 697–699, 1992.
- M. J. Feiler and K. S. Narendra, “Simultaneous identification and control of time-varying systems,” in Proceedings of the 45th IEEE Conference on Decision and Control (CDC '06), pp. 1093–1098, December 2006.
- N. Ozay, M. Sznaier, C. M. Lagoa, and O. I. Camps, “A sparsification approach to set membership identification of switched affine systems,” IEEE Transactions on Automatic Control, vol. 57, no. 3, pp. 634–648, 2012.
- T. P. Dinh, H. M. Le, H. A. L. Thi, and F. Lauer, “A difference of convex functions algorithm for switched linear regression,” IEEE Transactions on Automatic Control, vol. 59, no. 8, pp. 2277–2282, 2014.
- J. Roll, A. Bemporad, and L. Ljung, “Identification of piecewise affine systems via mixed-integer programming,” Automatica, vol. 40, no. 1, pp. 37–50, 2004.
- I. Maruta, T.-H. Kim, D. Song, and T. Sugie, “Synthesis of fixed-structure robust controllers using a constrained particle swarm optimizer with cyclic neighborhood topology,” Expert Systems with Applications, vol. 40, no. 9, pp. 3595–3605, 2013.
- H. D. Taghirad and P. R. Bélanger, “Modeling and parameter identification of harmonic drive systems,” Journal of Dynamic Systems, Measurement and Control, Transactions of the ASME, vol. 120, no. 4, pp. 439–444, 1998.
- E. R.-W. Bai, “Parametrization and adaptive compensation of friction forces,” International Journal of Adaptive Control and Signal Processing, vol. 11, no. 1, pp. 21–31, 1997.
- C. Canudas De Wit and P. Lischinsky, “Adaptive friction compensation with partially known dynamic friction model,” International Journal of Adaptive Control and Signal Processing, vol. 11, no. 1, pp. 65–80, 1997.
- I. Maruta, T. Sugie, and T.-H. Kim, “Identification of multiple mode models via distributed particle swarm optimization,” in Proceedings of the 18th IFAC World Congress, pp. 7743–7748, 2011.
Copyright © 2017 Tae-Hyoung Kim et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.