Power System Oscillation Mode Prediction Based on the Lasso Method

This paper utilizes modern statistical and machine learning methodology to predict the oscillation mode of interest in complex power engineering systems. The damping ratio of the electromechanical oscillation mode is formulated as a function of the power of the generators and loads as well as bus voltage magnitudes in the entire power system. The celebrated Lasso algorithm is implemented to solve this high-dimension modeling problem. By the nature of the $L_{1}$ design, the Lasso algorithm can automatically render a sparse solution, and by eliminating redundant features, it provides desirable prediction power. The resultant model processes a simple structure, and it is easily interpretable. The precision of our sparse modeling framework is demonstrated in the context of an IEEE 50-Generator 145-Bus power network and an online learning framework for the power system oscillation mode prediction is also provided.


I. INTRODUCTION
Inter-area oscillations are one of the main concerns of small-signal stability in power systems [1]- [4]. These oscillations are inherent to power systems and can limit the power exchanged between areas if they are insufficiently damped [5]. Conventionally, the oscillation mode information can be calculated with the power system model linearized on a specific operation point (equilibrium point) [6]. However, this method is not practical in further online applications because the system changes on different operating points, and the correlation between the damping of an electromechanical oscillation mode and operating scenario is not explicit. Our goal in this paper is to bridge this gap by investigating how real-time changes in operating scenarios impact the damping of the oscillation mode -the prediction of mode damping on different operating scenarios.
Power system oscillation mode prediction is a helpful tool for day-ahead planning activities and online generator re-dispatch. Using the correlations between the damping of an electromechanical oscillation mode and operating scenarios, The associate editor coordinating the review of this manuscript and approving it for publication was Hui Liu . transmission system operators (TSOs) can optimize the 24h ahead load flow operating points with the desired power exchanged between areas and sufficient damping. On the other hand, by using this correlation and the operation data from the Phasor Measurement Units (PMUs), TSOs can monitor the electromechanical oscillation mode in real-time, and adjust system operation (generator re-dispatch) in an emergency for satisfied damping.
The aforementioned problem becomes increasingly attractive as a result of the increasing quantity of available measurements, such as Wide Area Monitoring Systems (WAMS). Various works have attempted to explore the possibility of establishing a bridge between the mode damping and the system operating condition. Moreover, these works split into two categories: 1) Find the direct relationship f (·) between the mode damping ζ and the system operation point P 0 , that is ζ = f (P 0 ); 2) Find the sensitivity g(·) of the mode damping ζ on the system operation point P 0 , that is ζ = g( P) P 0 , which is called power system online sensitivity identification (OSI).
Between these two categories, the former approach dominates the other. For instance, linear regression was implemented in [7], [8] to identify the damping sensitivity at an operation point, while weighted linear regression was employed in [9] to find critical variables in the identification procedure. Apart from the regression-type methods, in [10] the problem was formulated as a pattern classification problem, and decision trees were used to predict well or poorly damped oscillations using power flow data. In [11], neural networks were also applied to find the relationship between system variables and a dominant mode.
Some works concentrating on the power system online sensitivity identification were also brought to our attentions. A locally weighted linear regression (LWLR) method was applied in [12] to predict damping ratio of a dominant mode online. A k-nearest neighbors locally weighted ridge regression was used in [13] to predict electromechanical oscillation mode. Furthermore, the locally weighted ridge regression integrated with norm-2 Tikhonov-Phillips regularization was examined in [14] to overcome the collinearity problem encountered in power system electromechanical oscillation mode prediction.
However, most of these papers were dealing with relatively small-size power systems. Actually, when the power system of interest is of hundreds of buses or even more, the aforementioned methods such as multiple linear regression or ridge regression can hardly work, due to the well-known ''curse of dimensionality'' phenomenon [15] known to data scientists and statisticians. In order to overcome this problem, in this paper, we apply a sparse machine learning algorithm called the Lasso (least absolute shrinkage and selection operator), which was widely applied to power system community, such as power system transient stability problems [16], [17], total transfer capability online estimation [18], and voltage stability margin online monitoring [19]. This paper examines this approach in assessing, understanding, and predicting the small-signal behavior of large interconnected power systems and estimates the direct relationship between the mode damping and the system operation point. The validity of our approach is demonstrated by an IEEE 50-Generator 145-Bus power network, and the prediction accuracy is compared with two competing methods recently applied in the field. Namely, they are the multiple linear regression, and the ridge regression methods. This paper is organized as follows: Section II presents a brief overview of the relationship between electromechanical oscillation mode damping and operating scenario. System identification methods using the Lasso are introduced in Section III. Test cases and rules for generating multiple random power flow study cases for extensive eigenvalue analysis and eventually for system identification are discussed in Section IV. Results are presented and discussed in the context of two dominant electromechanical modes in Section V. Section V also shows the extension to apply our methodology in the online power system oscillation mode prediction. Section VI concludes the method on this paper.

II. THE RELATIONSHIP BETWEEN ELECTROMECHANICAL OSCILLATION MODE DAMPING AND OPERATING SCENARIO
To explain the relationship between electromechanical oscillation mode damping and operating scenario clearly, we use the classical electromechanical model [6], [20]- [22] with loads modeled as constant impedances. Consider an n-Machine, N -Bus power system. Machine i is modeled as a constant voltage E i behind a transient reactance x di . The motion of the i-th machine rotor angle δ i is modeled as P mi is the input mechanical power, P ei is the output electrical power, V j_re and V j_im are the real and imaginary parts of the bus voltage phasor at bus j which is the terminal bus of Machine i,V j = V j θ j = V j_re +jV j_im is the load bus voltage phasor, δ i is the i-th machine angle, ω i is the i-th machine speed, m i = 2H i where H i is the inertia of the i-th machine most commonly used in power system stability simulation programs, and D i is the i-th machine's damping coefficient. Note that j denotes the imaginary sign if it is not used as an index.
For each load bus i, the power flow equation is: where R L_ik , X L_ik , and B L_ik are the resistance, reactance, and line charging, respectively, of the line connecting buses i and k. Moreover, P ei and Q ei are the i-th generator's active and reactive electrical output power, respectively, if bus i is a generator bus. G i and B i are the load conductance and susceptance at bus i. V is the N -vector of the load bus voltage phasors. Equations (1), (2), and (5) can be combined into a vector form where M is the diagonal machine inertia matrix, f is a vector of acceleration torques, and g is the power flow equation of the transmission network. We linearize (8) on the operating scenario P 0 (δ 0 , V 0 ) (steady-state/equilibrium point) to obtain VOLUME 8, 2020 the linear model where δ is an n-vector of the machine angle deviations from δ 0 , ω is an n-vector of the machine speed deviations from ω 0 , and V is an N -vector of the load bus voltage phasor deviations from V 0 . D is an n × n diagonal matrix with elements being the damping coefficients. C is an n × n diagonal matrix with elements being the constant 2π f 0 . The n×n matrix K 1 , n×N matrix K 2 , and N ×n matrix K 3 consist of the partial derivatives of the power transfer between the machines and their terminal buses. The N × N matrix K 4 is the admittance matrix of the network. K 1 , K 2 K 3 , and K 4 are related to the equilibrium point P 0 (δ 0 , V 0 ). Considering that K 4 is nonsingular, we can substitute (11) into (9) to eliminate V to obtain a linearized electromechanical model reduced to the machine internal nodes: where We denote A 2n×2n as system state matrix in power system small-signal analysis. The damping ratio ζ i of each mode λ i (for i = 1, 2, . . . , n) is determined by the state matrix A 2n×2n , which is formulated by physical laws described above. As long as the grid topology is steady, the value of A 2n×2n is determined by component parameters (such as matrix M −1 D) and system operating conditions (matrix M −1 K). Therefore, in the case of a developed power grid whose grid topology and physical parameters are rarely changed, it is reasonable that the following link holds: where P 0 is the operating scenario. In other words, given an operating scenario P 0 , (K 1 , can be computed numerically step by step. And the corresponding ζ i is unique to the operating scenario P 0 . Therefore, the link in (14) can be reformulated as where f (·) is a function with input P 0 and output ζ i . The eigenvalue calculation of a matrix can refer to root calculation of this matrix's characteristic polynomial. Galois' theory indicates that there is no solution in radicals to general polynomial equations of degree five or higher with arbitrary coefficients [23]. This means that the eigenvalue of a n order matrix, for n ≥ 5, cannot be solved algebraically. Therefore, examining the inverse of matrix K 4 and the eigenvalue calculation of matrix A, it's impossible to get the explicit symbolic expression of the function f (·) with input P 0 and output ζ i .
Thus the relationship between equilibrium point P 0 and i-th mode's damping ratio ζ i can extend to a data-driven damping ratio prediction problem. The machine learning algorithms we use in modeling ζ i = f (P 0 ) will be described in Section.III.

III. SYSTEM IDENTIFICATION METHODS USING THE LASSO
In system identification, the training data of the form It is worth noting that the input variable is also sometimes referred to as the ''observational variable'', and the output variable is called the ''response variable''. Then system identification tries to learn the mapping between X i and Y i , i.e., f : R p → R 1 . In the power system oscillation mode prediction problem, for the training case indexed by i, the input variable X i = [X i 1 , · · · , X i p ] corresponds to the voltage magnitudes together with the active and reactive power of the generators and loads in the entire power system, whereas the output variable Y i corresponds to the damping ratio of a certain oscillation mode of interest. The aim of using artificial intelligence methods is to learn the relationship between them, so that with new observations, one can use the model to predict the responses accordingly, and the prediction error would be as small as possible. In order to achieve this task, one may use a large class of possible models. It is, however, practical to start with a simple yet very useful model based on the linear mapping where ε i is the noise process. Here we refer to the discussions in the Subsection.V-D for more detailed justification of this simple model in the application of ultra-high-dimensional learning problems. For the model in (16), one can write The most straightforward approach to solve (16) is through minimizing the mean square error and the solution is the well known least square (LS) solution [24] expressed bŷ The least-square formula provides a solution to the linear model, yet it naturally has the drawback of numerical instability, i.e., a small change in the data, such as when new observations come to be available, may result in a large variation of the estimateβ. Also, the matrix X T X might not be full rank, creating a problematic scenario in the matrix inversion process.
In order to overcome the defects in the classical least square approach, the so-called ridge regression was developed in the 1970s [25], [26]. It is to minimize the least square criterion together with a penalty term on the weights. Namely, where λ is a regularization parameter that needs to be specified. It controls the amount of shrinkage in the minimization criterion. In practice, one can select λ according to some re-sampling methods, i.e., cross validation. It is worth to mention that the ridge regression also has close-form solution where I p is a p × p identity matrix. The ridge regression method processes the numerical stability in the solution and has been widely applied in engineering and applied science. It was applied in the small signal stability problem in [12].
In the late 1990s, the so-called ''least absolute shrinkage and selection operator'' (Lasso) was developed [27]. It is also based on minimization of the least square criterion together with a penalty term. However, the penalty term is the L 1 norm of the weights, rather than the L 2 norm case as in the ridge regression scenario. Specifically, the Lasso has the following formβ where λ is the regularization term that needs to be specified. It is worth mentioning that according to the dual optimization theory, the ridge regression and the Lasso also have the following equivalent primal forms, respectively, β Lasso,primal = arg min β,||β|| 1 <s where || · || 2 and || · || 1 denote the L 2 norm and the L 1 norm. Furthermore, s is a term that can be expressed alternatively in terms of λ.
The difference between Lasso and ridge regression can be shown in the following illustration in Fig. 1 [28].
In Fig. 1, it shows an example of a bivariate case. The sum of squares is shown, and the shaded area corresponds to the weight constraints. The left panel corresponds to the L 2 case, and the shaded area is a circular shape. In higher dimensional cases, it would be a hyper-ball. The right panel corresponds to the L 1 case, and the shaded area is a rectangular shape. It can be seen that sometimes the optimal solution can happen at the corner points in the Lasso case, suggesting that part of all the features can be set to zero weight. Such a scenario can never happen in the ridge case, meaning that ridge regression can never render a sparse solution. In high-dimension problems, it is known that the proper sparse solution can greatly enhance the model prediction accuracy due to the well known ''curse of dimensionality'' [15] in data science. As a consequence, it is expected that in our damping ratio prediction problem, the implementation of the Lasso algorithm can greatly increase the model prediction performance.
Unfortunately, the solution of Lasso does not have closedform formula as the least square case (19) or the ridge case (21). In practice, after specification of regularization parameter λ, people use coordinate optimization method to calculate the solution for Lasso. This algorithm is shown in Algorithm 1, and it is also called the ''shooting algorithm'' for the Lasso.
Note that in Algorithm 1, (c) + equals to c if c > 0, and equals to 0 otherwise. Besides the β j , j = 1, · · · , p, β 0 can be estimated by taking the empirical mean of the sample {Y 1 , · · · , Y n }. It is also worth mentioning, in order for each feature to have equal strength, the user should standardize the data so that each covariate of X (except for the first column) should have zero mean and unit variance. This can be done by pre-processing data when the user subtracts the mean and then divides by the standard deviation for each feature. It is worth mentioning that it is standard practice for statistical researchers to pre-process data in such a way, before applying the machine learning / system identification methods.
In practice, the specification of the regularization parameter λ plays a critical role. It can be seen that larger λ will For j = 1, · · · , p 5: 6: where −j is the same as β [m] except the j-th component is set to zero vector, = n −1 X X, andˆ jj is the j-th diagonal component ofˆ , 7: until numerical convergence shrink more weights to zero. In practice, one can use crossvalidation methods to select λ. Namely, on a grid of candidate λ values, estimate the model according to Algorithm 1, find the cross-validation errors, and specify in the final modeling the value of λ according to the one leads to the smallest crossvalidation error.

IV. DATA GENERATION A. GENERATING MULTIPLE POWERFLOW CASES FOR EIGENVALUE ANALYSIS
The system base case powerflow is randomly adjusted to simulate the real power system. The training data were generated in the following way to randomly generate powerflow cases whilst not exceeding the system limits that maintain realistic operation of the components. First, every load is randomly adjusted within ±25% variation from the base case. For the kth steady-state point, active and reactive power consumption on the i-th load bus are given by where P  Secondly, the active power injection and voltage at the i-th generator bus are given by the result of the optimal power flow (OPF) with the loads' new active and reactive power. The generator cost data are provided by Matpower software package [29]. If the optimal power flow converges, the OPF solution will be the k-th steady-state point.
After the k-th steady-state point is carried out, the initial value of each component can be calculated to form the system state matrix. So we can get the eigenvalue and eigenvector matrix of the system state matrix.

B. MODE TRACKING
Let row vector W 0 i and column vector V 0 i denote the normalized left and right eigenvector of the specific mode λ i at the base case P 0 , respectively. After each adjustment of the loads, the optimal power flow is solved, and an eigenvalue analysis is performed. Then we can denote column vector V k as the normalized right eigenvector of mode λ k (k = 1, 2, · · · , n) at the new steady-state point P 0 .
For one specified electromechanical mode λ i , the following method has been implemented to make sure the mode has been tracked correctly after each steady-state point P 0 was produced: 1) Exclude those modes whose highest participation factors are not the angle and speed states; 2) Use a reasonable range of values to bound the frequency of the electromechanical mode λ i (for instance, we can use 0.5 ± 0.2 Hz on the 0.5-Hz mode), and exclude those modes which fall outside this range; 3) Calculatek = arg min k=1,2,··· ,m ||W 0 i V k − 1||, where m is the number of mode set after excluding on Step 1 and 2, arg min is argument of the minimum. λk is the specific electromechanical mode we want to track. For the aforementioned Step 3, in an equilibrium point, the normalized left and right eigenvectors corresponding to different eigenvalues are orthogonal, W i V j = 0 for i = j. In the cases that eigenvectors correspond to the same eigenvalue, we have W i V i = 1.
Physically speaking, the right eigenvectors V 0 i and V i of the same electromechanical mode λ i at the power base P 0 and the new steady-state point P 0 respectively, must be similar, i.e., V 0 i ≈ V i . Thus, according to the aforementioned expression W 0 is the number of the mode λ i we want to track at the new steady-state point P 0 .

C. DATE RECORD
For each adjustment of the loads, the optimal power flow is solved, and an eigenvalue analysis is performed. The active and reactive power generation, as well as the bus voltage magnitudes (the bus voltage angles are not recommended as features) are recorded. The frequency and damping ratio of specific electromechanical mode are also recorded.

V. SIMULATION STUDIES A. IEEE 50-GENERATOR 145-BUS SYSTEM TEST CASE
The performance of our proposed modeling methodology is evaluated using an IEEE 50-Generator 145-Bus system test case. Its network topology structure is shown in Fig. 2. This test system consists of 145 buses, 50 generating units (7 machines use sub-transient reactance model and 43 machines use classical model), 64 impedance loads, 453 transmission lines. We refer to [29] for a detailed description of this power system. In this 145 bus system, we have 373 features: 1) Features #1 to #145 correspond to magnitudes of bus voltage. 2) Features #146 to #195, and Features #196 to #245 correspond to active power and reactive power of generator units sorted by bus number, respectively. 3) Features #246 to #309, and Features #310 to #373 correspond to active power and reactive power of loads sorted by bus number, respectively. Four hundred and ninety-nine (499) different operating scenarios are generated by the method proposed in Section.IV. Two different modes following are adapted to demonstrate this method. These two modes are: 1) The intra-area 1.3-Hz mode related to Generators within Area-1.
2) The inter-area 0.33-Hz mode in which generators in Area-2 oscillates against those in Area-3. We denote these modes by Mode 1 and Mode 2. It is worth noting that Mode 1 is a simple oscillation mode corresponding to an intra-area oscillation, while Mode 2 is a more complexed oscillation mode corresponding to an inter-area oscillation. In the following parts, we refer to Mode 1 as the simple mode, and Mode 2 as the hard mode. We use the first 250 cases as the training set, and the rest 249 cases as the testing set.

B. THE INTRA-AREA MODE
From the small-signal analysis, one can see that the 1.3-Hz mode is an intra-area mode in which generators on bus-98, 100 and 103 oscillate against generator on bus-89 within Area-1. The top 10 participation factor states of this mode are shown in Table 1, and the mode shape is shown in Fig. 3.
After the proper selection of regularization parameter λ described in Section.III, the final model consists of only 8 features. The strength of these features are visually demonstrated in Fig. 4. We refer to the beginning of Section V-A for the description of these features. Besides, we also show these features in terms of their physical meaning in Table 2. Furthermore, Fig. 5 plots the intra-area 1.3-Hz mode's damping ratio versus the top 6 dominant features using the sample data.  By applying the model to the testing data set, we find the root mean sqaure (RMSE) testing error is 0.00007359, or 0.007359%. Alternatively, if we evaluate the testing error by the mean absolute error (MAE), it is 0.00005518, or 0.005518%. The comparison between the first 100 predicted values and their true values are shown in Fig. 6. We can see that there is almost no noticeable difference between the two sets of values.
From Fig. 5, the correlation between intra-area 1.3-Hz mode damping and the top 6 dominant features look consistent with the weight of corresponding features acquired VOLUME 8, 2020   by Lasso in Table 2. Moreover, the results from Lasso indicate that the active power of generators on bus-89 and bus-108 (close to bus-98) are negatively correlated with the intra-area 1.3-Hz mode damping, which are compatible with the premise derived from small-signal analysis that this mode is an intra-area mode in which generators on bus-98, 100 and 103 oscillate against generator on bus-89 within Area-1.
The intra-area modes are typically determined by only a few number of machines. This explains why we observe very accurate prediction results. Our modeling also identify those most important features automatically.

C. THE INTER-AREA MODE
From the small-signal analysis, one can see that the 0.33-Hz mode is an inter-area mode in which the generators in Area-2 (generators on bus-128, 104, 105, 121, etc) oscillate against those in Area-3 (generators on bus-135, 136, 137 and 139). The top 10 participation factor states of this mode are shown in Table 3 and the mode shape is shown in Fig. 7.
After the proper selection of regularization parameter λ described in Section.III, the final model consists of only 24 features. The strength of these features are visually   demonstrated in Fig. 8. We refer to the beginning of Section V-A for the description of these features. Besides, we also show the top 10 most prominent features in terms of their physical meaning in Table 4. Furthermore, Fig. 9 plots the intra-area 0.33-Hz mode's damping ratio versus the top 6 dominant features using the sample data. By applying the model to the testing data set, we find the root mean sqaure (RMSE) testing error is 0.0134, or 1.34%. Alternatively, if we evaluate the testing error by the mean absolute error (MAE), it is 0.0083, or 0.83%. The comparison between the first 100 predicted values and their true values are illustrated in Fig. 10.
From Fig. 9, the correlation between inter-area 0.33-Hz mode damping and the top 6 dominant features look consistent with the weight of corresponding features acquired by Lasso in Table 4. And results from Lasso indicate that active and reactive power of generators/loads within Area-2 and Area-3 are correlated with the inter-area 0.33-Hz mode damping, which are compatible with the premise derived from small-signal analysis that this mode is an inter-area mode in which generators in Area-2 oscillate against those in Area-3.

D. DISCUSSIONS AND COMPARISON WITH OTHER ALGORITHMS
In Section V-C, we have applied the Lasso algorithms in the damping ratio prediction problem of an inter-area oscillation mode (a hard mode). In comparison, we also apply the multiple linear regression and the ridge regression using all the 373 candidate features. The results are shown in Table 5. From Table 5, we can see that multiple linear regression can hardly do the work in our problem, in which the dimensionality equals to 373. This confirms the discussions on its drawbacks in Section III. Actually, when this method is applied in [7], the researchers only apply the modeling using a subset of all the features based on their subjective choice originated from the industry experience. Comparing to this method, the Lasso can automatically decide which features should exist in the final model, thus circumvent the subjectivity.
Our method also outperforms the ridge regression approach. It is conjectured that in larger power networks (e.g., thousands of buses, or more), the superiority of our method will be better demonstrated.
So far we have applied high-dimensioned linear modeling in the power system damping mode prediction problem. VOLUME 8, 2020 One may come to the natural question how good is it to represent an unknown nonlinear system by a linear model. Generally speaking, such representation will incur an irreducible error called ''approximation error'', and [30] shows that under a wide range of conditions, the approximation error is quite affordable. This is reflected by the fact that our final modeling result. A 1.34% RMSE for the Interarea Mode and a negligible prediction error for the Intraarea Mode have achieved the modeling accuracy confirming the acceptable level of the approximation error due to the linear model. As matter of fact, the ''curse of dimensionality'' phenomenon [15] suggests that in ultra-high-dimensionality problem, the dimension-reduction is always much more crucial compared to the adoption of linear structures to represent the general nonlinear system. This reflects the merits of our paper.

E. IMPLEMENTATION IN ONLINE LEARNING
So far we have introduced the Lasso approach in the power system oscillation mode prediction problem. Generally speaking, the classical Lasso approach is an offline methodology. Nonetheless, Combing with the so-called Kalman filter [31], this approach can be used for online learning scenarios.
The Kalman filter is a type of techniques that allow the model to be constantly updated, as new training data come to availability. Due to the recursive nature of the Lasso shotting algorithm, the Kalman filter for Lasso might be hard to derive. On the other hand, since the Lasso technique quite often reduces the dimensionality dramatically, one can always use Lasso as the first step, and on the reduced feature space, apply Kalman filter version of classical linear regression instead.
With a bit of abuse of notations, let where {(1),(2),. . . ,(q)} are a subset of index, and have been selected optimally by the Lasso method. Letβ be the model based on dataset D n = {(X 1 , Y 1 ), . . . , (X n , Y 1 )}. Now new data comes to the availability, then Y and X are the updated data sets. Let β denote the updated model parameters. Then following Equation (19): Then we can easily derive (X X + X n+1 X n+1 )β = X Y + X n+1 Y n+1 .

VI. CONCLUSION
In this paper, we have applied the Lasso method to model the power system oscillation mode damping prediction problem. The sparse nature of our methodology enables the modeling to be applied in large-scale power systems, and automatic feature selection is achieved simultaneously with the regression process. The performance of the method is evaluated using the IEEE 50-Generator 145-Bus testing case. Our method shows a stronger prediction accuracy comparing to two competing methods recently applied in the field. The correlations between the electromechanical oscillation mode damping and operating scenarios can enable the power engineering users to monitor and control some power system quantities to ensure the desirable level of damping, even if she/he processes little background understanding about the system. Importantly, the automatic feature selection property of the proposed approach can largely reduce the range of these power system quantities affecting the specific oscillation mode in a large power network. The framework of industry online application using the Lasso algorithm was shown in Fig. 11. Instead of using Prony analysis at the first step in Fig. 11, some novel power system oscillation mode estimation methods [32]- [34] can also calculate the power oscillation mode of interest from power system sampled signal in real-time. In a practical system, noises exist in measurement data. Therefore, the sensitivity analysis of the training data can be the further work of this study. HAOYONG CHEN (Senior Member, IEEE) received the B.S., M.S.E., and Ph.D. degrees in electrical engineering from Xi'an Jiaotong University, in 1995, 1997, and 2000, respectively. He was an Influential Professor of electrical engineering (in both academia and industry). He has 14 years' research and teaching experience as a University Faculty Member of electrical engineering. He is currently the Director of the Institute of Power Economics and Electricity Markets, South China University of Technology. He has been the Principal of many key national projects of China and has close connections with the Chinese Power Industry. His current research interests include excellent research record on modeling/optimization/control of power systems and electricity markets.