Kautz Function Based Continuous-Time Model Predictive Controller for Load Frequency Control in a Multi-Area Power System

A continuous-time Model Predictive Controller was proposed using Kautz function in order to improve the performance of Load Frequency Control (LFC). A dynamic model of an interconnected power system was used for Model Predictive Controller (MPC) design. MPC predicts the future trajectory of the dynamic model by calculating the optimal closed loop feedback gain matrix. In this paper, the optimal closed loop feedback gain matrix was calculated using Kautz function. Being an Orthonormal Basis Function (OBF), Kautz function has an advantage of solving complex pole-based nonlinear system. Genetic Algorithm (GA) was applied to optimally tune the Kautz function-based MPC. A constraint based on phase plane analysis was implemented with the cost function in order to improve the robustness of the Kautz function-based MPC. The proposed method was simulated with three area interconnected power system and the efficiency of the proposed method was measured and exhibited by comparing with conventional Proportional and Integral (PI) controller and Linear Quadratic Regulation (LQR).


Introduction
In large electrical power systems, the load changes cannot be determined due to which the deviations occur in real and reactive power towards an unstable condition. So there is a need for continuous regulation of electric power generation for an efficient and reliable power system operation. The change in real power causes a change in frequency whereas the change in reactive power causes a change in voltage magnitude. Load Frequency Control (LFC) is a control strategy to minimize the tie line power and frequency deviations while at the same time maintaining the frequency within accepted limits [Saadat (1999); Kothari (2003)]. Utilizing fixed controllers is observed as conventional method to solve LFC problem. Fixed controllers such as Integral controller (I), Proportional-Integral (PI) and Proportional-Integral-Derivative (PID) controller are widely used in process control.
algorithm. Kautz functions in MPC controller to improve control and stability of the autonomous car vehicle trajectory [Yakub and Mori (2014)]. Discrete and continuoustime MPC with both Laguerre and Kautz function was designed and explained with simple examples in the literature [Wang (2009)]. This paper presents continuous-time MPC using Kautz function to solve LFC (Load Frequency Control) problem in an interconnected power system. A phase plane constraint was included in the optimization problem to ensure the stability of the system during uncertain disturbances. The proposed method was simulated in three area hydro-thermal interconnected power system. The results were compared with a conventional PI controller and LQR controller to prove the effectiveness of the proposed method.
2 Modeling of multi-area power system The electrical power system was divided into a number of different areas. Each area was interconnected with its neighboring area through transmission lines called 'tie lines'. In an interconnected power system, the load deviation leads to frequency deviation and tie line power deviation. LFC was used to maintain the frequency and tie line power within the desired limits in each area. Thus, in order for the LFC to be efficiently controlled, the design of dynamic LFC model was formulated considering both frequency as well as tie line power. In the current research work, thermal, hydro and wind power plant were considered for the study purpose. The basic block diagrams of thermal and hydropower plants in an interconnected system are shown in the Fig. 1 and Fig. 2 [Kumar (2017); Shiroei (2014)]. The power system parameters considered for this study were listed in Appendix A. The mathematical equations for the dynamic model of both thermal and hydro plants are as follows [Kumar and Suhag (2017); Shiroei and Ranjbar (2014); Bangal (2009)]. Governor block (Thermal) is The mathematical equations for frequency deviation and tie-line power flow are same for all the interconnected systems which are represented as follows: Frequency deviation is given by the equation Tie line power flow is given by the equation To maintain area frequency and tie-line power exchange with neighboring areas within the preset values, a supplementary control with a linear combination of frequency and the net interchange power for each area 'i' known as the Area Control Error (ACE) is to be reduced.

Approximation and realization
Some of the systems or network might be difficult to solve due to its nonlinear structure. Such systems need to be approximated for obtaining desired response [Kautz (1954); Magni and Scattolini (2004)]. In approximation, the input function was approximated with some intermediate functions and then was optimized. Once the function was approximated to the desired form, it was converted back to a physically-realizable function. The basic steps of approximation are as follows: Step 1: Simplification of function In approximation theory, the first step is to express in detail about the operation or response of a function (which are difficult to analyze) in the form of simpler functions. If an arbitrary response ( ) is considered for a system which can be approximated by a Fourier series expansion over the interval (0, ∞) as where ( ) denotes a set of exponentially-damped sinusoids and orthogonal The coefficients of the expansion Then the arbitrary response of the first 'n' terms is * ( ) = ∑ ( )

=1
(10) Step 2: Determination of optimum pole locations The location of the poles and zeros gives an approximate understanding of the system's response characteristics. The roots of the characteristic equation, also known as eigenvalues of the system is equal to the poles of the system response. The proper selection of pole locations guides to the simplified form of the coefficient .
Step 3: Orthonormal set The orthonormal condition for the function with unity weight is expressed as Step 4: Formulation of The coefficient is formulated and evaluated in terms of orthonormal set ( ).
Step 5: Minimization of error The mean-squared error between ( ) and * ( ) has to be minimized to find the optimal pole location and to approximate the given arbitrary function identical to the physical system. If the ( ) is assumed as a piecewise continuous function, the integral squared error in an orthonormal expansion is given by the following equation When expanding the Eq. (12) and assuming the error as '0' then, The Eq. (13) can be written as Thus the equation (14) becomes the condition to be satisfied for proper convergence of the desired response.

Kautz function
The orthonormal basis functions are mostly used for system approximation and identification. A dynamic model with OBF gives approximate knowledge about the dominant dynamics such as time constants and damping factors of the system which can be included in the system identification process [Akcay and Ninness (1999)]. Some of the orthonormal basis functions that are widely used are Laguerre base function and Kautz base function. Kautz functions are an extended generalization of Laguerre functions. Unlike Laguerre functions which only deal with identical real poles, Kautz functions possess the ability to approximate the system function with both real and complex poles. Thus Kautz function-based dynamic model can approximate the system in resonance.

Real pole
If all the poles are real, non-identical and all poles >0 whose pole locations lie at -1 , -2 … -…, -, then the Kautz network can be represented as where =2 … and ' ' is the maximum number of poles. If all the pole locations are identical, then this is called a Laguerre function.

Complex pole
Let the non-identical complex poles are denoted as − ± for = 1,2, … , where − > 0 for all . Then Kautz networks with complex poles are represented as and ' ' is the maximum number of terms

State space form
Kautz functions can be expressed in a state space form as follows where ( ), ( ) and ( ) are the state vector, the unit impulse function and the Kautz network terms respectively. , and are matrices determined by the locations of the poles.

Model predictive controller (MPC)
Model Predictive Controller (MPC) has the ability to solve multivariable problem and handle complex system constraints. These advantages make MPC a well suited option for power system and to overcome LFC problem [Shiroei, Toulabi and Ranjbar (2013)]. MPC uses the explicit model to predict the future response of the system. It calculates the controller input by solving objective function with constraints through optimization. Once the optimal trajectory of the future response of the system was calculated, only the first input sample of the control signal would be taken for system processes based on Receding Horizon Control (RHC) whereas the remaining input samples are ignored. This optimization is performed again for next time instant with a new set of system information. State space model is the most common type of MPC structure used to represent the dynamical system. The mathematical model of a linear time-invariant multivariable system with -input, -output, -disturbances and -states can be expressed in a state space form as follows: where ( ), ( ), ( ) and ( )are state vector, output vector, control vector, and disturbance vector respectively. The matrix coefficients , , and are system matrix, the input matrix and disturbance matrix respectively. The number of output should be less than the number of input .

Kautz function-based MPC (K-MPC)
The modeling and implementation of the Kautz function-based MPC controller for continuous-time application are described here as per the literature [ Wang (2009)]. The K-MPC method was designed by describing system model in continuous-time while computations of the optimization problem were performed in discrete time (in off-line) [Minh and Rashid (2012); Magni and Scattolini (2004)]. The optimization problem for control trajectory included manipulated variables whereas its optimal tuning using GA is explained below.

Modeling trajectory of the future control signal
The moving time window of the optimization process is assumed to be between ≤ ≤ + . ' ' is considered as the current time and ' ' as the maximum duration of the moving time window. The Eq. (14) stated that, for a stable closed loop system, the control trajectory ( ) should converge exponentially to zero when there is a transient state response due to load disturbance. If the disturbance is continuous, then the control trajectory ( ) converges to a specific constant value, instead of zero. So, it is wise to optimize derivative of the control trajectory ̇( ). The derivative of the control trajectory ̇( ) is expressed as follows where ( )is set of OBF's and is the vector of coefficients.

Augmented state space model
As the Eq. (19) concludes, there is a need for integral action occurs once the derivative of the control trajectory was obtained. So the design of system state space model need to be modified to an augmented state space form. With an assumption that ( ) =̇ is the augmented state space form of (18a) and (18b) with ̇( ) as control input, this can be defined as follows. [̇( Where an identity matrix A zero matrix ̇( ) A derivative of the control trajectory ( ) Continuous-time integrated zero-mean white noise ( ) A constant vector of a reference signal The dimension of the augmented state space is 2 = + .

Selection of poles
Kautz function requires the optimal selection of poles. The poles can be calculated from the priori information of Linear Quadratic Regulator (LQR) using the augmented state space model of the system to be solved. The closed-loop eigenvalues of the system are calculated by solving the characteristic equation where and are matrix coefficients of an augmented state space model, is the eigenvalues of system matrix( − ), is an identity matrix and is a gain matrix obtained by minimizing the cost function of LQR over a finite time horizon. Assuming = , the quadratic cost function of the LQR is expressed as: where = × , is a matrix coefficient of augmented state space model of the system, is a constant multiplier and × is an identity matrix.

Computation of Kautz function
The optimal poles (eigenvalues), obtained by solving the Eq. (21), was used for the computation of Kautz network. Then this Kautz network terms again were used for calculating the state space parameters of (17a) and (17b). Kautz function, using state space parameters, is defined as where = 1,2, ⋯ , When substituting (23) in (19), it gives the derivative of the control trajectory ̇( ) in terms of Kautz function as where ( )is a Kautz function and is the vector of coefficients.

Design of optimal controller
The main objective of the K-MPC is to design an optimal control trajectory for MPC by minimizing the quadratic cost function. The optimal control trajectory ̇( ) is designed using a controller gain matrix . Assuming that there is a time ' ' which lies within the interval 0 ≤ ≤ and = , the optimal cost function of a system is expressed as: where is a diagonal matrix and each diagonal element of is a real matrix with dimension ' 2 × 2 '.
As the first term of (25) is negligible and it is the only term in the cost function dependent on , the optimal is = −Ω −1 Ψ ( ) (28)

Computation of Coefficient vector
The cost function Eq. (25) clearly states that the coefficient vector plays a vital role in the design of the control trajectory. The vector of coefficients is calculated from the Eq.

Computation of matrices and
The matrices and can be evaluated using the Eqs. (26) and (27) in discrete form as they do not rely on current time . By dividing time in discrete as = 0, ℏ, ⋯ , ℏ, such that = ℏ, the matrices and can be evaluated as follows: where ℏ is a constant step size.

Convolution integral ̂( )
The matrices Ω and Ψ are the functions of ( ). So ( ) has to be calculated first followed by ( ) can be calculated by computing ̂( ) . The -th input of the convolution integral ̂( ) is obtained by solving the linear algebraic equation Above equation is solved in discrete form to obtain optimal ̂( ) . After ̂( ) is calculated, ( ) can be calculated by substituting ̂( ) in ( ) =̂( ) (32)

Implementation of control signal
According to Receding Horizon Control (RHC), the first signal of ̇( ) alone is considered for the control input of the system. Assuming a random time ' ', the derivative of the optimal control is ̇( )can be again defined by splitting into state and output component as, where and are the gain matrices for state variables and system output.
is a diagonal matrix with each diagonal values, a gain value of the integral controller.
Finally, the control trajectory was obtained by integrating the derivative of the control trajectory. Therefore, Substituting the Eq. (35) in (36), the control trajectory in terms of state variables is given by where ( ) is a reference signal.

Evaluation of system response
To evaluate the response of a closed loop system, the root of the characteristic equation of the closed-loop system is required. To find the response of the system, the derivative of the control trajectory (33) The Eq. (39) is the system output response of the closed loop feedback system in ' ' domain. By taking the inverse Laplace transform of the Eq. (39) results in the differential equation of the system is obtained in time domain. This differential equation is used as the input for the objective function for optimization in the following section.

Constrained optimization problem
In this paper, to obtain an optimal controlled response, one of the well-known performance criteria was used as an objective function. An Integral Absolute Error (IAE) is a performance measure which integrates absolute of error over time. Let be the objective function to be solved. The objective function is given by Subject to, > 0 and < 0 (41) where = Determinant of ( − )

= Eigenvalues of ( − )
To ensure the stability of the system, the Eq. (41) was implemented into the system cost function as a penalty factor. Stability constraints were applied as a penalty factor in the cost function of MPC problem [Magni and Scattolini (2004)]. Inserting stability constraints in optimization problem for uncertainty is described in the literature as well [Bemporad and Morari (1999)].

Phase-plane analysis of eigenvalues
In this paper, a new technique was proposed by implementing the constraints using penalty factor in the cost function, by phase plane analysis, in order to determine an optimal solution in MPC such that the system remains stable. Phase plane analysis is an analysis method to observe the features of a dynamic system in a graphical manner. To observe the system behavior, the solution of the system at equilibrium is to be determined which is called a 'critical point' and the vector plane that represents the behavior of the system is known as 'phase portrait'. The path of the solution in phase portrait is viewed as a moving particle in a curve or line. Eigenvalues and eigenvectors are the most convenient ways of representing the differential equation for phase plane analysis. Depending on whether the eigenvalues are real or complex and positive or negative, the system behavior can be determined as stable or not. This kind of analysis is well suitable for oscillatory systems.

Optimal tuning of
The eigenvalues of the system matrix ( − ) are the poles of the given system. The gain matrix has to be selected optimally for a stable operation of the system. The optimal ensures the derivative of the control trajectory to exponentially decay to zero. For optimal , the parameters such as , ℏ, , and have to be fine tuned.
There are several optimization algorithms available for finding optimal solution for the given optimization problem. If the optimization problem is a multimodal type, then conventional iterative type algorithm cannot guarantee an optimal solution. As the multiarea power system is a complex and nonlinear problem, which is difficult to optimize with conventional methods, a meta-heuristic algorithm such as the Genetic Algorithm was applied to the proposed methodology for optimal tuning.

Genetic algorithm
Genetic Algorithm (GA) is an Evolutionary Algorithm (EA) based on Charles Darwin theory of natural evolutionary process in biological systems. GA is a prominent optimization algorithm which can hold large search space and has a greater chance of determining the global best solution for an optimization problem [McCall (2005)]. GA also has the capability to handle multiple variables which makes it suitable for solving both constraint and unconstraint optimization problems. The standard procedure of GA consists of biologically-inspired operators known as selection, crossover, and mutation [Rao, Rao and Dattaguru (2004)]. GA starts with the initialization of the population and the population consists of individuals whom are randomly generated. Each individual of the population is represented as 'real value' or 'strings'. At each generation, the potential of every chromosome is evaluated with a value known as 'fitness' by solving the objective function. Once the fitness of each chromosome is evaluated, GA chooses healthier chromosomes among all the chromosomes based on fitness value. Crossover is a recombination of chromosomes from the selection process to form a new set of chromosomes. In mutation process, every gene position, in an individual chromosome of a newly-produced chromosome, is interchanged with randomly generated numbers to form a new population for next generation. The whole process is repeated again with the new population produced after mutation until the best solution is reached or a stopping criterion is reached.

Algorithm steps for the proposed methodology
The proposed methodology to tune with GA algorithm is summarized in the following steps: Step 1 Initialize GA parameters such as search space , number of chromosomes ℂ, number of generations , crossover rate ℛ and mutation rate ℳℛ.
Step 2 Initialize populations ℙ for search parameters , ℏ, , and with random values and set generation counter as 0.
Step 3 Calculate eigenvalues using LQR with parameter .
Step 4 Calculate Kautz state space matrices , and with eigenvalues.
Step 7 Calculate matrices and with parameters ℏ, and .
Step 8 Calculate the feedback gain matrix with parameter .
Step 9 Evaluate the system response using ( − ).
Step 10 Evaluate the fitness value for each chromosome in the current generation.
Step 11 For each population, the parent chromosomes are selected based on fitness values.
Step 12 Parent chromosomes are assigned with random numbers ' '. If ℛ > , crossover takes place and new children are produced. Else no children are produced.
Step 13 The child chromosomes are assigned with random numbers ' '. If ℳℛ > , mutation takes place and the new population is produced.
Step 14 Increment the counter and repeat the steps 3 to 12 till the stopping criterion is reached.
Step 15 The final solution is the parameters to form optimal .

Simulation results and discussion
The three area power system examples from the literature [Bangal (2009)] were used for the simulation. Area 1 and Area 2 were two identical thermal non-reheat power plants whereas Area 3 was a hydropower plant. The dynamic system for MPC was designed assuming that all the areas would be provided with load change. The experiment was conducted using two case studies. In Case 1, a step load was applied to area 1 and area 2 alone whereas in the Case 2, step load change was applied to all areas. In both the cases, step load change with magnitude 0.1 pu was applied at a time instant of 10 s. The system responses were compared with conventional PI controller, LQR and Kautz function-based MPC (K-MPC).
Case 1: The output responses for the case study 1 are shown in the Fig. 3. In case 1, all three controller performances were found to be almost identical. As area 3 was not applied with step load change, the oscillations in both frequency deviation and tie line power of area 3 alone were much lesser when compared to other two area frequency deviation and tie line power. Oscillations in area 2 responses were more. All controller outputs were stable even if the system was modeled assuming all areas would have the same load. But none of the controller performances seem to be optimal for all the three areas. Area 1 system output for LQR was better than both K-MPC and PI. Area 2 system output for K-MPC was better than both LQR and PI. PI controller showed much lesser performance than both K-MPC and LQR. Case 2: In this case, all the areas were provided with step load change. The output responses for case study 2 are shown in the Fig. 4. It is evident from the comparison of output responses that the proposed method exhibited optimal performance when compared to other two controllers. The responses of the PI controller were better than LQR unlike as in Case 1. LQR struggled to maintain the stability of system response.

Conclusion
A Kautz function-based MPC (K-MPC) with phase plane stability was proposed for the optimal control of LFC problem. For simplicity, the design of the MPC used in this problem is of a centralized type. GA was presented to optimize the parameters of K-MPC.
To prove the efficiency and robustness of the proposed method, it was examined on a multi-area interconnected power system problem that consisted of two thermal non-reheat power plants and a hydropower plant through simulation. Simulations were carried out using MATLAB\Simulink. The experiment was conducted as two case studies by including and excluding the load in the hydropower plant. The K-MPC is basically a PI controller. So, the proposed method was compared with both LQR and the conventional PI controller. The simulation results inferred that when the stability constraint is included in the cost function, the proposed method is able to provide reliable and stable closed loop output, even in case of any uncertainties present in the system. If the proposed method is implemented in the system without uncertainties, it is predicted to have better performance than other methods.