A FAST MULTI-STEP PREDICTION AND ROLLING OPTIMIZATION EXCITATION CONTROL METHOD FOR MULTI-MACHINE POWER SYSTEM

Original scientific paper A fast excitation predictive control method for multi-machine power system is presented. The multi-step prediction technique is realized via system dynamic model. Some inequality constraints on states, inputs and outputs are considered in rolling optimization. The Gramian balanced reduction technique and the improved optimization algorithm are used in order to save the time of open-loop optimization in model predictive control. A 50machine power system is used to verify the effectiveness of this approach. Compared with simulated results under different controllers, this method can greatly reduce the calculating-time. The voltages of generator terminals are maintained within the set points. The stability of power system is improved.


Introduction
The synchronous generator excitation control system can maintain the voltage of generator terminals and pivot points at the reasonable range.It is one of the most efficient methods to improve power system stability.Lots of excitation control approaches have been applied to power system generators [1][2][3][4][5], such as robust control, optimal control and model predictive control (MPC).
MPC is an open-loop optimization and close-loop control method based on a system dynamic model [5][6][7][8][9][10][11][12][13][14][15][16][17].It can use multi-step prediction technique to predict the dynamic change trend of the system.The constraints on states, inputs and outputs are considered in rolling optimization.Compared with other methods, MPC can greatly improve the closed-loop control performance.The control target, even dynamic model constraints are easily adjusted in rolling optimization.
With the development of smart grid, power system will be autonomous, predictable, self-healing.In order to get a good control results, power system needs to make a beforehand response to the uncertain factors.So many researchers have been applying MPC to power system recently [5][6][7][8][9][10][11][12][13][14][15].Reference [6] demonstrates that MPC for energy storage systems can improve the dispatchability of wind power plants.The application of MPC can ensure the transient stability of a power system after a contingency in reference [10].A method of decentralized nonlinear excitation predictive controllers for multimachine power system to enhance their transient stability is presented in reference [13].However, the main obstacle for application of MPC is to solve the optimization problem quickly within one sampling period for rolling optimization, especially when the states' dimension is growing rapidly with the increase of scale and complexity of the system model.
For MPC of generator excitation (MPCGE) in power system, most researchers directly derive the analytical and/or single-step prediction optimization control law based on power system model [13][14][15].These methods need not solve dynamic optimization problems.They cannot realize multi-step prediction.The constraints on system states, inputs and outputs are also not considered.So they are lacking in flexibility of the generic MPC.
In this paper, a fast MPCGE method is presented.In this method multi-step prediction and rolling optimization are realized.Some inequality constraints on system states, inputs and outputs are considered.The Gramian balanced realization technique and the improved interior-point method are used in order to speed up numerical calculation.A fast algorithm, i.e. low-rank Cholesky factor-alternation direction implicit (LRCF-ADI), is used to solve the large-scale Lyapunov equations in the Gramian balanced reduction.

Multi-machine power system excitation model predictive control
Fig. 1 is the simulation model of excitation multi-step predictive control for multi-machine power system.The model is composed of the nonlinear multi-machine power system model and its linear excitation multi-step predictive controller.
In Section 2.1, it accounts for a nonlinear multimachine power system model -detailed generators, transmission network and loads.Each generator is described by a 6-order utility dynamic model.Loads are composed of the constant impedance load and the nonconforming loads.The non-conforming loads contain the constant current load and the constant power load, according to [18][19][20].The process of load recovery is considered in Section 2.1.The linear model of nonlinear multi-machine power system model is used to design the linear excitation multi-step predictive controller in Section 2.2.The principles of Section 2.1 and 2.2 are given in the following.

Nonlinear Multi-machine Power System Model
A generator can be described by the following 6order dynamic model according to [18].

(
) The detailed expressions of d-axis and q-axis voltages are described as, , , , , , , where i = 1, 2, …, h, and h is the number of generators.The equations of transmission network are expressed by The connection among generators, transmission network and loads is explained in the following.
1) I xyG and V xyG are transformed to I dq and V dq respectively via the following relations I xyG =SCI dq and V xyG =CV dq .I dq can be described as a function of V dq in Eq. (2).C is a xy-dq coordinates transformation matrix.S is a capacity conversion matrix.
2) The constant impedance loads are added in the Y.
3) The non-conforming loads are connected to the transmission network Eq.(3) via the following Eq.( 4), which are functions of the voltages of load buses [19,20].

(
) where P L and Q L are functions of the voltage vector V xyL .The nonlinear interface equations can be described as a function of voltage vectors of generators and nonconforming loads as follows, ) Eq. ( 1) and Eq. ( 5) can be simplified as where, The active power vector P e and the terminal voltage vector V t of generators are chosen as outputs of nonlinear multi-machine power system model and they are described as follows, e d d q q P = V I +V I (7a) Eq. ( 7) can be simplified as follows, ( ) = y h x (8) where, , T m e t t h 6) and Eq. ( 8) are used as a controlled objection to simulate actual multi-machine power system.The linear form of Eq. ( 6) and Eq. ( 8) is used to research excitation multi-step predictive controller.

Linear excitation multi-step predictive controller
The linear model of Eq. ( 6) and Eq. ( 8), which includes generators, transmission network and loads, can be expressed as follows [21].
Linear excitation multi-step predictive controller chooses the least-square residual of system states and inputs as an objective function.Its equality constraints are the linear dynamic equations of power system Eq.(9).Its inequality constraints are the limits of states, inputs and outputs.The p-step MPCGE optimization problem at the instant t is described as follows.The optimization control u * (t) is acted as the actual input in the close-loop control.
where, τ = 0, 1, …, p−1, Q and R correspond to the weighting matrices; p is the number of predictive step, x(τ), y(τ) and u(τ) correspond to the predictive values of states, outputs and inputs at the instant (t+τ), respectively.
x min , x max , y min , y max , u min and u max are limits of states, outputs and inputs respectively.Eq. ( 10b) is the discrete form for Eq. ( 9).The predictive horizon T p in MPC can be calculated by T p =s×p.s is the size of predictive step in Eq. (10b), and its value is equal to the sampling period T s .The control horizon T c is equal to T p .
Generally, the matrices A, B and C are calculated by a stable operating point.When the operating condition of power system is changed, the stable operating point will be changed, and the matrices should be calculated again via the new one in order to get a more accurate predictive model.

The challenge and novelty relative to previous work
The following problems can be solved by the MPCGE Eq. ( 10) for multi-machine power system: 1) For the traditional certain excitation control methods i.e., proportion integral differential (PID) controller in [22] and linear optimal (LO) controller, they cannot predict and use the dynamic information of power system in future.PID controller can only implement constrains using saturations about Ef , and LO controller cannot consider constraints on states, inputs and outputs.
2) For other MPCGE methods (mentioned in Section 1), they cannot realize multi-step prediction and rolling optimization.At the same time, it is difficult to consider inequality constraints on system states, inputs and outputs.
3) If one of the control targets, constraints and dynamic model of the system needs to be adjusted in traditional methods, the control strategy has to be redesigned.
Although the proposed method is well used to handle the problems above, how to solve the dynamic optimization problem Eq. ( 10) quickly within a sampling period is still a main challenge to this method, especially when the dimension of power system dynamic model is high.

The Gramian balanced reduction technique and improved 3.1 Linear MPCGE based on Gramian balanced reduction model
For linear model order reduction problems, as considered in systems theory and control of ordinary differential or partial differential equations, balanced truncation and related methods are the methods of choice since they have some desirable properties: it preserves the stability of the system, and it provides a global computable error bound between the transfer function of the original system and that of the reduced-order system [23,24].
As first computational step, the computation of the controllability Gramian matrix W C and observability Gramian matrix W O is required.The Gramians are equivalently given by the solutions of two large-scale Lyapunov equations [25].Solving these equations is the challenge of balanced truncation method.So a fast algorithm, i.e.LRCF-ADI is used [26][27][28].The basic idea of LRCF-ADI is taking low-rank matrices to approximate the Gramians.In this algorithm, some ADI parameters ζ i are used to reconstruct Lyapunov equations to compute W C and W O .
Generally W C and W O are unequal, so the influences of a state variable on inputs and outputs of system are not consistent, i.e., a certain state variable has a great influence on the inputs, but it has little influence on the outputs, or vice verse.Moore proposed a balanced realization reduction approach in [29] in order to compromise this influence.This method is used to convert the original system to a balanced one where the controllable and observable Gramian matrices are equal, through the balanced transformation matrix T.
Let x = Tx  , the balanced system of Eq. ( 9) can be obtained.

= + x Ax Bu
where is the vector of states in balanced system,

C = CT 
Referring to [30], Hankel singular values σ i of balanced system and Galerkin projection matrix P are used to truncate those states which have little influence on the inputs and outputs of Eq. ( 11), according to the threshold ε in Eq. (12).P is comprised by identity matrix and zero matrix.
where,  is the order of reduction system, so the balanced reduction model of linear dynamic model Eq. ( 9) can be calculated.
is the vector of retained principal states of the balanced system Eq.(11) In the balanced reduction system Eq.( 13), x   has no specific physical meaning, and it cannot be measured directly.However, the balanced reduction system is usable.x   can be calculated by balanced transformation matrix T and Galerkin projection matrix P, thus, , is a linear combination regarding x.Even though there is only one state in balanced reduction system, the state also contains the main dynamic information of all states in the original system.
According to Eq. ( 13) and the transformation,  , the linear MPCGE Eq. ( 10) is changed into the MPCGE based on reduced model (MPCGERM), Subject to ( 1) ( ) ( ) ,min ,max ( ) where, Q  represents the weighting matrix of the balanced reduction system states, and ,min x   , ,max x   , are limits.The relationship of relative parameters between Eq. ( 10) and Eq. ( 14

Improved interior-point method
Improved interior-point algorithm (IIPA) is used to solve the optimization problem Eq. ( 14) in order to further decrease the numerical calculating time.Eq. ( 14) can be simplified as where,  (16)   where k>0 is a obstacle factor, and v is a dual variable.
According to KKT conditions and Taylor expansion, Newton method is used to solve Eq. ( 17).The solution (Δz,Δv) is used to correct (z,v).The (z,v) has to meet the inequality constraints in iterations in order to find the optimization solution.
The warm start technique is used to reduce iteration times.At the instant t, it takes the previously computed trajectory of the instant t − 1 as a starting point.Roughly speaking, in MPC the current t * z can be computed by working out starting points int z for p-step prediction.This technique takes the previously computed

) t p t p t p
the constraints, int t z will also meet them, so the iteration times are decreased.The iteration of solving Eq. ( 17) is not stopped until the norm of the residual becomes small enough, or the iteration times reached the limit N max (e.g., N max =50).Much more calculating time will be spent in this case, and this paper gives a small value to N max based on the following reasons: on the one hand, the start warm technique decreases the iteration times; on the other hand, in close-loop control MPC only uses the current control u * (t), and the multi-step prediction can ensure that u * (t) does not have a bad influence on the future dynamic behavior of the system.Therefore, a small value is reasonable.The limit of iteration times N max can be chosen in 5-10 by extensive trials.The controlled result in this case has a little deviation compared with the one where N max is 50.
The warm start and the setting of iteration times' limit can overcome the main obstacle of MPC to a certain degree.Fig. 2 is a block diagram of the control analyzed from Section 3.1 to 3.2.

Simulations and analysis
The proposed method has been applied to a 50machine power system which includes 49 generators, a reference machine, 145 nodes and 453 lines in [35].In view of the constant impedance and non-conforming loads, the dimension of multi-machine power system linear model Eq.( 9) can reach 294-order, and the number of input variables is 49.
Fault Scenario: An instantaneous three-phase short circuit fault occurred at t=1.00s in an arbitrary line, it was cleared after holding 0.10 s, then a three-phase recloser operated at t=2. 10s and made the system return to the normal operation condition.IIPA is used in this section.
The method of how to choose the predictive step p and the sampling period T s is explained as the following.For stationary T s , with the increase of p, the control performance and the time in numerical calculation are also improved steadily.The choice of p itself is a trade off.As far as T s , if it is too long, the control result will be bad owing to lack of sampling information, as there will be a big pressure on calculating time.The choice of T s is also a tradeoff.Therefore a number of simulated trials are used to get a compromise.
It is difficult to provide and analyze all generators' dynamic behaviors for 50-machine power system in this paper, and some of them are chosen discretionarily as representative examples.

Analysis of control performances for different excitation controllers
The control performances of PID controller, LO controller and MPCGE based on full-order model (MPCGEFM) controller with two different kinds of constrains and control targets (MPCGEFM1 and MPCGEFM2) are discussed in the following.In simulation, the predictive step and the sampling period are 3 and 20 ms, respectively.The predictive horizon is 60ms.The maximum iteration times N max is 5. Generator 43 is taken as an example.The simulated results for it under different controllers are shown in Fig. 3 including curves of active power, terminal voltage and angular velocity.
As shown in Fig. 3, the maximum deviations under the control of LO and MPCGEFM are smaller than those of PID, so are oscillation times.The reason is that PID is not an optimization controller, thus is, there is no optimization function.The prediction and the inequality constraints on states and outputs cannot be considered in PID.The active power curves in Fig. 3 (a) are taken as an example to further illustrate the advantages of MPCGEFM.
i) Multi-step prediction and inequality constraints can be considered in MPCGEFM.The maximum deviation under the control of LO is larger than the one of MPCGEFM as same as oscillation amplitude.The reason is that dynamic change trend of the system in future and inequality constraints are not considered in LO where the optimal target is considered only.
ii) The optimization target and constraints in MPCGEFM can be adjusted easily.In Fig. 3 (a), the curves under the controls of MPCGEFM1 and MPCGEFM2 are different.The reason is that the weighting matrix Q and the limits about outputs y referring to Eq. (10a) and Eq.(10e) in MPCGEFM1 are different from the ones in MPCGEFM2.LO failed in this regard.The reason is that LO is only a determinate optimization control law (u * =−Kx).If the control target is adjusted, the law has to be re-designed.Therefore, to a certain degree, the control performance and flexibility of MPCGEFM is better than the ones of PID, LO.

Analysis of excitation predictive control based balanced reduction model
The comparison is made in Section 4.1 to illustrate the advantages of MPCGEFM.In this section, the model order reduction of 50-machine power system is analyzed.The reduced model forms MPCGERM to decrease the numerical calculating time.The distribution of the Hankel singular values in the balanced system is shown in Fig. 4. When the threshold ε is 1.0, the 294-order linear dynamic model of 50-machine power system can be projected onto Gramian balanced model.The balanced model can form MPCGE based on full-order balanced model (MPCGEFBM) controller.As shown in Fig. 4, the Hankel singular values after the σ 75 are very close to zero.The first 75 states contain most energy of full-order system.It seems feasible that the 294-order balanced system can be reduced to 75-order one, even lower.The smaller the ε is, the lower the order is, and also the bigger the distortion of the reduced system dynamic behavior is.
The active power, terminal voltage and angular frequency curves of reduced (23, 30 and 39 orders) and unreduced (full-order) systems with the same PIDC are shown in Fig. 5.The threshold ε is 1.0, 0.9, 0.8 and 0.7, and the corresponding order is 294, 39, 30 and 23 respectively.Generator 6 is nearby the location of threephase short circuit fault and taken as an example.
As shown in Fig. 5, the deviations between reduced system and unreduced one are increasing with the decreasing of the order of reduced system.Compared with other-order systems, there are large deviations in dynamic behaviors for the 23-order system.The bigger the distortion is, the worse the control performance of low-order MPCGERM may be.The 30-order and 39order models are chosen to form corresponding MPCGERMs respectively.In the following, a comparison is made among MPCGEFM, MPCGEFBM, 39-order MPCGERM and 30-order MPCGERM.
Generator 40 is taken as an example.The simulated results for it under these controllers are shown in Fig. 6.In simulation, the predictive step and the sampling period are 2 and 20 ms, respectively.The predictive horizon is 40 ms.The maximum iteration times N max is 5.As can be seen from Fig. 6, the control performance of MPCGEFBM is the same as that of MPCGEFM.The 30-order MPCGERM has some deviations compared with MPCGEFM.The 39-order MPCGERM is a better one.The simulated results show that MPCGERM is feasible, and the Gramian balanced reduction technique can be applied to a large-scale power system.

Statistics analysis of computational burden, calculating time, and memory storage
In many excitation control strategies, such as PIDC, LOC and some simple excitation predictive control laws mentioned in Section 1, the related parameters even control laws have been designed offline.There is no rolling optimization in these controllers.So these controllers can save more numerical calculating time compared with MPCGEFM.However, as shown in Section 4.1 and 4.2, the control performance and flexibility of them are worse than MPCGEFM.The proposed method coordinates the numerical calculating time and the control performance of MPCGEFM.The Gramian balanced reduction technique and improved the interior-point method are introduced to decrease numerical calculating time.In the following, some statistics are provided to analyze the numerical calculating time of the proposed method.10) with the same predictive step are approximate.As shown in Tab.1, the time of numerical calculation is gradually increasing with the lengthening of p for T 1 , T 2 and T 3 .IIPA can cut down the calculating time most effectively from the μ 13 and μ 23 .and CN 4 represent their computational burden to a certain degree respectively.A simplified efficiency ρ ij =(CN i −CN j )/CN i ×100% is defined, which can reflect the decreased degree of calculating burden in MPCGERM.
As can be seen from Tab. 2 and Fig. 6, the 39-order MPCGERM optimization problem with p=2 can be solved within a sampling period (T s =20 ms) and the control performance is good.Analyzing lots of simulated data, the numerical calculating time of MPCGERM is gradually lengthening with the increase of p.The Gramian balanced reduction technique can effectively lessen the calculating time, and sometimes it makes optimization problem solved within a sampling period.Some configurations would appear to be actually applicable.
Statistics of the memory storage (MS) occupied by PID, LO, MPCGEFM and 39-order MPCGERM in the processes of simulation are shown in Tab. 3.Where MS F , MS R , MS P and MS L are the memory storage occupied by MPCGEFM, 39-order MPCGERM, PID and LO, respectively.A storage-efficiency η=(MS F −MS R )/MS F ×100% is defined to show the effectiveness of the Gramian balanced reduction technique in term of lessening memory storage.
As shown in Tab. 3, MS F and MS R are gradually lengthening with the increase of p.The Gramian balanced reduction technique can decrease the memory storage occupied by MPCGEFM.MS R is larger than MS P and MS L .The reason is that PID and LO are determinate control laws, and MPCGERM needs rolling optimization.The Fortran, C, C++, and other programming languages can be used for speeding up this method and reducing memory storage further.
These simulated results in Section 4 show that MPCGERM can greatly reduce the calculating-time.At the same time, the stability of power system is improved.

Conclusion
A fast linear excitation multi-step prediction and rolling optimization control method is presented in this paper.The numerical calculating time is greatly shortened by the Gramian balanced reduction technique and improved interior-point algorithm.Compared with PIDC, LOC and other excitation predictive control laws, the multi-step prediction and rolling optimization are realized in the presented approach.At the same time, the method provides more flexibility of the generic MPC.Last but not least, the control performance of MPCGERM is similar to the one of MPCGEFM.The voltages of generator terminals are maintained within the set points.The stability of power system is improved.

Figure 1
Figure 1Multi-machine power system excitation predictive control Gramian balanced reduction can greatly decrease the number of constraints in Eq.(10).

Figure 2
Figure 2 Block diagram of the control analysed

Figure 4
Figure 4 Hankel singular values of the balanced system

Figure 5
Simulated results for the unreduced and reduced systems

Figure 6
Simulated results for generator 40 under different controllers In Tab. 1, a comparison is made among the optimization algorithm in MATLAB Toolbox (Fuction f mincon, FC), interior-point algorithm (IPA) and IIPA with different predictive steps.The IPA and IIPA are programmed in MATLAB.The numerical solutions of optimization problem Eq. (

Table 1
Numerical calculating time of three algorithms with different predictive steps

Table 2
Numerical calculating time of MPCGEFM and 39-order MPCGERM with different predictive steps 3is the mean numerical calculating time of 39order MPCGERM.CN 3 is the number of constrains (10bf) in MPCGEFM.CN 4 is the number of constrains (13b-f) in 39-order MPCGERM.The optimization problems of MPCGEFM and MPCGERM are solved by IIPA.CN3