Bogdan Marinescu
email: bogdan.marinescu@rte-france.com

David Petesch
email: david.petesch@rte-france.com

Three-Level Coordination in Power System Stabilization

Keywords: PSS tuning, robustness, coordination, speed governors

This paper presents a new methodology for the synthesis of Power System Stabilizers (PSSs) and speed governors in order to satisfy some recent objectives and constraints imposed by the evolution of large-scale interconnected power systems. As a consequence of the increasing size of these systems, the frequency of some inter-area modes diminish so much so that these oscillatory phenomena interact with the dynamics of the speed governors. These low-frequency inter-area modes are also more spread (i.e., they involve a larger number of generators) than the other inter-area modes of the system. This two particularities require several supplementary levels of coordination when tuning the PSSs and the governors of the machines involved in these oscillatory phenomena. To ensure this, the synthesis of the aforementioned regulators is done using a new control model which takes into account the key interactions of the dynamics of the power system. It thus allows the coordination of control actions among several generators as well as between the electrical and mechanical paths. This control framework ensures also an optimal trade-off performance/robustness. Validation simulations have been carried out with the Eurostag software on a realistic large-scale model of the interconnected European power system.

Introduction

Power system stability analysis is generally done in a local way, i.e., study of the local oscillatory phenomena or the so-called local modes which frequency is usually around 1Hz. The power system stabilizers (PSS), which are additional voltage regulation loops for the generators, are used to damp these local modes. In power systems there are also global oscillatory phenomena, i.e., the so-called inter-area modes (see, e.g., [START_REF] Rogers | Power System Oscillations[END_REF]). For the European system, the frequencies of these modes actually range between 0.4 and 0.17Hz. The inter-area modes have also to be damped using the same PSS loops used for the local modes which must thus be tuned (at least) for these two objectives. Notice that the frequency ranges of the two phenomena are slightly different. This is called in the sequel the first level of coordination.

Contrary to the local modes, the inter-area modes present two properties. First, they involve a large number of distant machines. Indeed, the damping cannot be achieved by tuning only one PSS, coordination between the tuning actions for different machines must also be envisaged. This is a second level of coordination which is required.

To ensure the two above mentioned levels of coordination, information from the whole system should be taken into account when tuning the parameters of the PSS. One way to do this is to use distant measures which can be provided by wide-area measurement systems (WAMS). An alternative is to use a reduced-order model of the whole system, called a control model. This second path is followed here as it does not need extra material for the implementation. In [START_REF] Gupta | Design of Decentralized Power System Stabilizers for Multimachine Power System using Model Reduction and Fast Output Sampling Techniques[END_REF] and related references, the control model is obtained by reducing a time-domain model of the overall power system using standard balancing techniques. Other approaches use diagonal modal decomposition for the same scope. None of these approaches can be used for largescale systems. To overcome this difficulty, we developed here a reduced-order model directly focused on the modes of interest, i.e., the modes to be damped. This model is refined by frequency identification for which only a frequency response (Bode plots) is necessary. The latter can be easily obtained even for large-scale models as indicated in Section 2. This kind of control model is of very low order and can be considered as an extension of the one usually used for the tuning of PSS and which consists of the sensitivities of the modes of interest with respect to the gains of the PSSs (see, e.g., [START_REF] Rouco | Coordinated Design of Multiple Controllers for Damping Power System Oscillations[END_REF]). It can be alternatively obtained by time-domain identification as in [START_REF] Grondin | An Approach to PSS Design for Transient Stability Improvement through Supplementary Damping of the Common Low-Frequency[END_REF], [START_REF] Kamwa | Low-Order Black-Box Models for Control System Design in Large Power Systems[END_REF], [START_REF] Kamwa | Small Signal Analysis of Hydro-Turbine Governors in Large Interconnected Power Plants[END_REF] and related references.

Next, the frequency of the inter-area modes depend on the size of the power system. The connection of the Turkey zone to the ENTSO-E (European Network of Trans-mission System Operators for Electricity) synchronous zone led to a new inter-area mode at 0.17Hz. Notice that before this extension, the slowest European interarea mode was at 0.22Hz. Moreover, the latest ENTSO-E interconnection feasibility study [START_REF] Luther | Feasibility Aspects of a Synchronous Coupling of the IPS/UPS with the UCTE[END_REF] has shown that, with this latter potential extension, an inter-area mode at 0.07Hz may exist. Low frequency inter-area modes have already been identified on the Hydro-Québec network [START_REF] Grondin | An Approach to PSS Design for Transient Stability Improvement through Supplementary Damping of the Common Low-Frequency[END_REF] and across the interconnection between the Northwest and Southwest American power systems [START_REF] Schleif | Damping of System Oscillations with a Hydrogenerating Unit[END_REF]. When the frequencies of the inter-area modes diminish so much, these oscillatory phenomena interact with the dynamics of the speed governors as shown in [START_REF] Rogers | Power System Oscillations[END_REF]. These interactions have been taken into account in the work presented here by coordinating the synthesis of the PSSs with the one of the speed governors. This corresponds to a third level of coordination. It has already been shown in [START_REF] Schleif | Damping of System Oscillations with a Hydrogenerating Unit[END_REF] that the governors can be used to damp system oscillations. In [START_REF] Kiss | Realization of a "Band-Stop" Device to Damp Inter-Area Oscillations with Intervention into the Turbine Governor[END_REF] these oscillations are "filtered" at the governor loop level. In [START_REF] Okou | Decentralized multivariable voltage and speed regulator for large-scale power systems with guarantee of stability and transient performance[END_REF], modern control methods are used to re-design at a glance the AVR and the governor functions of a given system which leads to a new controller given in a general (i.e., state-space) form. In the present work, standard PSSs (i.e., generally IEEE lead-lag forms [START_REF]IEEE Recommended Practice for Excitation System Models for Power System Stability Studies[END_REF] used in practice) are tuned in a coordinated manner with supplementary control actions in standard governors. For this, the control model mentioned above is extended in order to take into account the governor's dynamics.

Robustness is an important topic in modern control, especially in the case of power systems where, due to the large scale of the problem, the models used for the dynamic studies are a priori simplified. A lot of work has been dedicated to robust control techniques. The more direct approach consists in tuning the parameters of given standard PSSs' structures in order to achieve the desired performances and level of robustness. This leads to a problem of optimization under constraints [START_REF] Kamwa | Optimization-Based Tuning and Coordination of Flexible Damping Controllers for Bulk Power Systems[END_REF], [START_REF] Kamwa | Robust Design and Coordination of Multiple Damping Controllers Using Nonlinear Constraints Optimization[END_REF]. This framework has been used in the work reported here for the tuning of fixed standard structures of PSSs and speed governors for several machines in order to ensure the multilevel coordinations specified above.

The paper is organized as follows: in Section 2, a control model is built for the first two levels of coordination. It is further used in Section 3 for the control synthesis. In Section 4 it is explained why and how the synthesis of the PSSs and governors should be coordinated. The control model introduced in Section 2 as well as the control method in Section 3 are extended in Section 5 to the general case in which all the three levels of coordination can be achieved. Section 6 is devoted to conclusions and remarks.

All the simulations presented in this paper are carried out with the Eurostag transient stability simulator software [START_REF] Meyer | Eurostag, a Single Tool for Power System Simulation[END_REF].

Multi-frequency control model

Degrees of detail of the modeling

To synthesize a regulator, a reduced order dynamic model is needed for the system to be controlled. It is called a control model and its particularity is to capture only the features which are relevant to the control objectives. In general, the control model most used for the machines of the power systems consists in a machine connected to an infinite bus through a line of variable reactance. The value of the reactance of this line is used to model the short-circuit power at the grid bus to which the machine is connected. However, doing so, it is difficult to capture both local and inter-area oscillatory phenomena especially when they are of quite different frequencies. To move towards a coordination of tuning of the PSSs for several local and inter-area modes in the new context mentioned in the Introduction, sensitivities between the damping of the modes of interest and the gains of the PSSs selected for tuning were used as a control model in [START_REF] Rouco | Coordinated Design of Multiple Controllers for Damping Power System Oscillations[END_REF] and related references as follows : the sensitivity of a closed-loop eigenvalue λ i with respect to the gain K j of the PSS transfer function Γ j (s) is given by [START_REF] Taranto | Combined use of Analytic and Genetic Algorithms for Robust Coordinated Tuning of Power System Damping Controllers[END_REF]

∂λ i ∂K j = r jj i (1)
where r jj i is the residue of λ i in the open-loop transfer function H jj (s) from V s j to ω j in Fig. 1. V s j is the stabilizing signal of the AVR of machine j, while ω j is the speed of the same machine. Notice also that the transfer function mentioned above is computed before the installation of the PSS Γ j (s). This residue can be computed from the right (v i ) and the left (w i ) eigenvectors of λ i and the input (B j ), respectively the output (C j ) matrices of a minimal state-space realization of H jj (s):

r jj i = C j v i w T i B j . ( 2
)
The residues in (2) are usual outputs of softwares for small-signal analysis like, e.g., [START_REF] Rouco | A computer Package for Analysis of Small Signal Stability in Large Electric Power Systems[END_REF]. From (1) one deduces the impact of small variations of the gain K j of the modes of interest:

λ i = λ 0 i + r i K j ( 3
)
where λ 0 i is the mode computed on the open-loop, i.e., the situation in Fig. 1 without Γ j (s).

However, this static characterization can be further enriched to take into account the dynamics of interest of the overall system. Overall dynamic models are usually available for the interconnected systems. They consist in a detailed model of the machines of the system of significant installed power (usually the ones of more than 100 MW) along with their regulations and the 400 kV/220 kV transmission grid (see [START_REF] Breulmann | Analysis and Damping of Inter-area Oscillations in the UCTE-CENTREL Power-System[END_REF] and [START_REF] Luther | Feasibility Aspects of a Synchronous Coupling of the IPS/UPS with the UCTE[END_REF] for the case of the European system). Although simplified with respect to reality, this kind of model, called simulation model is too complex (about 20000 state variables in the case of the European system) for a control model. It is used for the full nonlinear simulation of the behavior of the power system, in particular for a posteriori validation of an already synthesized controller. However, it can be used to extract a suitable control model as shown below.

Choice of the control model structure

Let Λ = {λ 1 , ..., λ l } be the oscillatory modes (local and inter-area) to be damped and M = {M 1 , ..., M m } the machines for which PSS loops have been chosen to be installed or adapted to perform the damping task. The latter are among the machines with the greatest participation factors in the modes in the set Λ. The way to choose these machines have been addressed in several preceding works (see, e.g., [START_REF] Rogers | Power System Oscillations[END_REF], [START_REF] Rouco | Coordinated Design of Multiple Controllers for Damping Power System Oscillations[END_REF]) and the set M is considered here as input data.

The control model concerns the power system seen from the PSSs (see Fig. 1), i.e., the transfer matrix H(s) between the stabilizing signals V s i of the PSS loops of the machines in the set M and their speeds ω j :

ω(s) = H(s)V s (s) ( 4 )
where

ω(s) = ⎡ ⎢ ⎣ ω 1 (s) . . . ω m (s) ⎤ ⎥ ⎦ and V s (s) = ⎡ ⎢ ⎣ V s 1 (s) . . . V sm (s) ⎤ ⎥ ⎦ . ( 5
)
This transfer matrix can be made available in practice using a mixed nonlinear/linear analysis framework for power systems like, e.g., [START_REF] Marinescu | A Unified Framework for Nonlinear Dynamic Simulation and Modal Analysis for Control of Large-Scale Power Systems[END_REF]. More precisely, the nonlinear system is linearized around a given operation point and, next, the frequency response from all the modes can be made available by numeric computation of the frequency response even for the large-scale systems. Bode plots of two entries of H(s) for the case of the ENTSO-E system are given in Fig. 3. A control model can be obtained from H(s) if each of its entries H ij (s) is written as a limited order development plus

a correction. Indeed, H ij (s) = n k=1 r ij k s-p k ,
where r ij k denotes the residue of the pole p k of H ij . Obviously, n, the total number of poles of H, is huge since it equals the order of the full simulation model (between 8000 and 20000 for the European system). However, for the control model, the concerned dynamics are defined within the frequency band of the modes in the set Λ. We propose as control model the following approximation of H ij (s), i, j ∈ {1, ..., m}:

Hij (s) = A(s) B(s) + P (s) Q(s) ; A(s) B(s) = l k=1 [ r ij k s -λ k + rij k s -λk ] ( 6 )
where r ij k , rij k , k = 1, ..., l are known (by rij k we denote the complex conjugate of r ij k ) and the polynomials P (s) and Q(s) are computed such that Hij (s) fits H ij (s) in the frequency working band mentioned above.

Test system

The techniques for the synthesis of the control model and the PSS loops investigated in this paper are tested and illustrated on a realistic large-scale representation of the interconnected European. More precisely, this is a representation of the European power system before the interconnection with zone 2 (Romania and Bulgaria) and Turkey. It consists of about 2000 buses, 2400 lines and 810 transformers. The generators of which power is greater than 100MW (about 400 machines) are represented by detailed dynamic models along with the detailed models of their regulations. The rest of the generation is considered as static injection at the load-flow stage. A winter peak load scenario is considered. The resulting linear model is described by about 8000 state variables. It is well-known that this system exhibits a low damped interarea oscillation around 0.22 Hz in which the generators of the eastern part of the grid are oscillating against the generators of the western part [START_REF] Breulmann | Analysis and Damping of Inter-area Oscillations in the UCTE-CENTREL Power-System[END_REF]. This phenomenon is represented by the first two modes of the linearized full model of which dampings are given in the first line of Table 1.

They are studied in this paper along with the one in the third column of the same table which is of different nature; it is an inter-area mode of the Spanish system at a slightly higher frequency (0.9 Hz) than the first two ones. It is chosen here along with the slow ones ( 1 and 2) in order to be in the situation which requires the "second level of coordination" mentioned in the Introduction. Thus, one has (6) with l = 3 and λ k , k = 1, ..., 3 given in Table 1.

Frequency identification of the parameters of the control model

The degrees of the polynomials P and Q in (3) are chosen such that Hij is strictly proper:

Q(s) = s v + q v-1 s v-1 + ... + q 0 , P(s) = p v-1 s v-1 + ... + p 0 (7)
To obtain a control model of order 2l + v, 2v parameters have to be identified as coefficients of P and Q. This is done using the Bode plots of the H ij (s) as in Fig. 3. They contain information in the wide frequency range of the full simulation model (i.e., about 8000 poles). However, for the stabilization problem we want to solve here, only the frequency range [ω - Λ ω + Λ ] which covers the set of modes Λ is of interest. For the chosen test system, this range is defined around the two major resonance peaks in Fig. 3 which correspond to λ 1 and λ 3 in Table 1 (the frequency of λ 2 is in between), i.e., ω - Λ = 0.15 Hz, ω + Λ = 2 Hz. An adequate control model should fit the Bode plots of the full simulation model in the frequency range [ω - Λ ω + Λ ]. The coefficients p i , q i in (7) of each transfer Hij (s) in [START_REF] Kamwa | Small Signal Analysis of Hydro-Turbine Governors in Large Interconnected Power Plants[END_REF] are thus computed via a frequency identification procedure based on a least squares objective function of the form

J ident = ω - Λ ≤ω k ≤ω + Λ [α k (A k -| Hij (iω k )|) 2 + β k (ϕ k -arctg( Hij (iω k ))) 2 ] ( 8
)
where A k and ϕ k are the values of the magnitude and respectively the phase of

H ij (iω k ) and i 2 = -1. (A k , ω k
) and (ϕ k , ω k ) are points of the Bode plots of the transfer of the full simulation model and are thus input data for the frequency identification problem

{p i , q i } i∈{0,...,v-1} = argmin{J ident }. ( 9
)
The weights α k , β k are used to manage the trade-off between magnitude and phase fitting and, eventually, to give priority to the fitting at specific frequencies.

A stability constraint on the roots of Q could be used in [START_REF] Kiss | Realization of a "Band-Stop" Device to Damp Inter-Area Oscillations with Intervention into the Turbine Governor[END_REF] to ensure a stable control model Hij (s).

This identification is repeated as a trial and error iterative procedure with increasing degrees for the polynomials P and Q above till an acceptable frequency response fitting is achieved.

Such a situation is shown in Fig. 3a for a transfer function of the diagonal of H(s). The responses of the full model are in solid lines and the reduced model is obtained with deg(P ) = 1 and deg(Q) = 2, i.e., for a control model of order 8.

For the transfers of the extra diagonal entries of H the fitting is more difficult since interactions among different machines should be captured. Fig. 3b shows the result obtained with the same order (v = 2) for two Spanish machines (i = Almaraz, j = Cofrentes). The resulting numerical values are given in Table 2. They were obtained with equal weights ( α k = β k = 1 for all k) and using standard Matlab routines (fminsearch/lsqnonlin). If necessary, accuracy of the identification can be further increased by increasing v. However, the order 2 has been retained for the example treated in this paper since the tuning provided with this model was satisfactory. Indeed, the main objective here is to tune with the simplest control model.

Optimal tuning of PSSs

The control model ( 6) is used to simultaneously tune the PSSs of the machines in the set M for desired target values ζ ref i , i = 1, ..., l of the damping of the modes in Λ. To ensure a standard IEEE structure for the PSSs, like, e.g.,

V s (s) ω(s) = K 1 + T 1 s 1 + T 2 s 1 + T 3 s 1 + T 4 s T 5 s 1 + T 6 s , (10)
(see [START_REF]IEEE Recommended Practice for Excitation System Models for Power System Stability Studies[END_REF]) the problem is formulated as an optimization one which gives an optimal set of gains and time constants for each machine (K and T i in [START_REF] Okou | Decentralized multivariable voltage and speed regulator for large-scale power systems with guarantee of stability and transient performance[END_REF]) in the set M.

The objective function

The objective function to be minimized should capture the dynamic performance specified, basically the damping of the oscillatory response of the system. In [START_REF] Simo | Validation of a new modal performance measure for flexible controllers design[END_REF] and [START_REF] Kamwa | Robust Design and Coordination of Multiple Damping Controllers Using Nonlinear Constraints Optimization[END_REF] a modal objective function is proposed:

J contr = l k=1 [(a k Ψ(σ k , ω k ) m i=1 m j=1 r ij k rij k )] 1/2 Ψ(σ k , ω k ) = σ 2 k +ω 2 k 2σ 3 k (e 2σ 3 k T σ 2 k +ω 2 k -1) (11)
where r ij k is the residue of the mode λ k = σ k ± jω k in the transfer H ij (s), a k a weighting function (a k = 1 if λ k is real and a k = 2 if λ k is complex) and T > 0 is the time horizon over which J contr is evaluated.

Roughly speaking, [START_REF]IEEE Recommended Practice for Excitation System Models for Power System Stability Studies[END_REF] contains the integral of the surface under the modal response envelope through a given time horizon T . First, the computation of this function is not straightforward in practice and, next, minimizing J contr given by ( 11) leads to a response with maximum damping which is not necessarily needed. Indeed, only the level ζ ref i , i = 1, ..., l of damping is required and not going below would allow one to settle a better trade-off robustness/performance (see, e.g., [START_REF] Skogestad | Multivariable Feedback Control[END_REF]). This led us to test a simpler and more direct index of performance:

J contr = l i=1 (ζ i -ζ ref i ) 2 . ( 12
)
Relation (3) can be exploited to obtain in a simplified manner the influence of the adjustment of the gains of the PSSs on the modes in Λ. Moreover, to reduce the size of the problem, it is considered that T 3 = T 1 and T 4 = T 2 in [START_REF] Okou | Decentralized multivariable voltage and speed regulator for large-scale power systems with guarantee of stability and transient performance[END_REF]. The optimal PSS parameters are thus:

{K * i , T * 1 , T * 2 } = argmin{J contr } ( 13
)
where J contr is now given by [START_REF] Kamwa | Optimization-Based Tuning and Coordination of Flexible Damping Controllers for Bulk Power Systems[END_REF]. The above objective functions have been compared on the Single Machine Infinite Bus (SMIB) system in Fig. 5 of which parameters are given in Appendix 1. Without the use of a PSS, the inter-area mode of this system is -0.0227 ± j0.4951, i.e., of damping ζ 0 = 4.58%. Problem (13) solved with [START_REF] Kamwa | Optimization-Based Tuning and Coordination of Flexible Damping Controllers for Bulk Power Systems[END_REF] for ζ ref = 10% leads to K = 46pu, T 1 = 2.17s and T 1 = 2.24s which improve in closed-loop the damping of the above inter-area mode to ζ = 10%. The same problem solved with (11) instead of ( 12) leads to K = 49.01pu, T 1 = 4.01s and T 1 = 4.01s which provide a damping of ζ = 0.1075. In this latter case the gain is larger and the damping is (a little) higher than the target. This is normal because in this case the damping target is not directly taken into account into the optimization as with [START_REF] Kamwa | Optimization-Based Tuning and Coordination of Flexible Damping Controllers for Bulk Power Systems[END_REF]. The trade-off between the gain magnitude and the obtained damping can de managed by adding constraints to problem [START_REF] Kamwa | Robust Design and Coordination of Multiple Damping Controllers Using Nonlinear Constraints Optimization[END_REF] as explained in the section below.

The constraints

Several types of constraints may be considered with ( 13) or [START_REF]IEEE Recommended Practice for Excitation System Models for Power System Stability Studies[END_REF]. First, physical bound constraints

⎧ ⎨ ⎩ T 1 i ≥ 0 T 2 i ≥ 0 K i ≥ 0 , i ∈ {1, ..., m} (14)
contain the physical bounds for the gains and the time constants of the PSSs. An additional upper bound K i ≤ K max i can be added to ( 14) to ensure implementable solutions, i.e., to avoid mathematical solutions with unrealistic high gain. Next, several nonlinear constraints are added. The main one consists in a stability condition of the resulting closed-loop (H(s) feedback connected with Γ(s) in Fig. 1. The eigenvalues of state matrix of the closed-loop should be in the left-half plane:

Max{Re(λ cl k )} < 0, k ∈ {1, ..., 2l + v + 3m}. ( 15
)
Notice that λ cl k are functions of K i , T 1 i and T 2 i of the PSSs Γ(s) i of type [START_REF] Okou | Decentralized multivariable voltage and speed regulator for large-scale power systems with guarantee of stability and transient performance[END_REF]. In [START_REF] Kamwa | Robust Design and Coordination of Multiple Damping Controllers Using Nonlinear Constraints Optimization[END_REF] and [START_REF] Marinescu | Robust Coordinated tuning of Parameteres of Standard Power Systems Stabilizers for Local and Global Grid Objectives[END_REF], it was shown how constraints given by H ∞ -norm of adequate transfer functions may be used to ensure robustness.

Validation on the European power system

For the test system presented in Section 2.3, the Spanish machines Almaraz, Cofrentes and PGR were chosen to damp the modes selected in Table 1 since they have high participation in these modes and they are not already equipped with PSSs. The damping target is ζ ref = 10% for each of the three modes. This objective is interpreted as follows: mode #1 is poorly damped and ζ 1 should thus be increased, but the damping actions should be chosen in order to not degrade the damping of the other two modes, one directly concerned by the east-west oscillation and the other one local to the Spanish system.

Linear validation

The results obtained minimizing (13) (first position in Table 3) are compared with the ones provided by the tuning proposed in [START_REF] Rouco | Coordinated Design of Multiple Controllers for Damping Power System Oscillations[END_REF] where only a static control model is used (second position in the same Table 3). First, the gains obtained with the procedure presented above are lower than the ones in [START_REF] Rouco | Coordinated Design of Multiple Controllers for Damping Power System Oscillations[END_REF] and thus the robustness is improved in general. This difference can also be noticed in the Bode plots of the open-loop transfer functions of the control system, i.e., of Γ i (s)H ii (s), given in Fig. 2 for i = Almaraz. The damping factors obtained when the optimal parameters are used for PSSs are given in the second line of Table 1. Also, notice that only two machines (Almaraz and Cofrentes) are needed to reach the damping objectives with the tuning procedure presented here since the gain obtained for PGR is K * = 0 as shown in Table 3. Notice that for this example no robustness constraints have been taken into account. The optimization problem (13) has been solved with standard Matlab routines (fmincon).

Full nonlinear validation

It is now shown that the damping level achieved in the preceding section using linear models is still valid in nonlinear simulation. For the latter validation, the system is simulated in detail using the Eurostag software for transient stability analysis [START_REF] Meyer | Eurostag, a Single Tool for Power System Simulation[END_REF]. Fig. 4 gives the speed responses of machines Almaraz to a short-circuit of 200ms at the grid connection point of the same machine. It can be seen that the damping of the response is improved (the more damped curve corresponds to the situation when the 2 PSSs described in the first entries of Table 3 are used), while the less damped response is obtained on the initial situation when no PSSs are used on the machines in the set M).

Remark

Notice that, as shown in Table 1, that damping target ζ ref = 10% is achieved for mode 3 and (almost for) mode 2 but not for mode 1 although the damping of the latter mode has been significantly improved. This is normal since it is the result of a constrained optimization. From the physical point of view, this means that the damping means are not sufficient. To further improve the damping in these situations, the class M of machines on which PSS are installed/tuned should be enlarged.

Coordinated tuning of PSSs and governors

In this section it is shown in which cases interactions between the PSS and governor dynamics may occur and how they can be overcome.

Interactions PSSs/governors

Consider first the simplest situation of the SMIB system in Fig. 5. This system has an inter-area mode which frequency depends (mainly) on the inertia H of the generator and the reactance X of the line. The latter parameters (given in Appendix 1 along with all the data of the test system) are chosen such that the frequency of this mode is f = 0.07Hz. This simple example reproduces a type of inter-are mode which exists on a large-scale power system with weak interconnections as the one treated in Section 5.3.

When a standard AVR is used, this mode is poorly damped (ζ = 0.93%) as shown by the solid line response in Fig. 6a to a 100ms short-circuit applied at the generator bus.

If a PSS loop is added, the oscillations are damped as shown by the response in dotted line in the same figure. The PSS structure [START_REF] Okou | Decentralized multivariable voltage and speed regulator for large-scale power systems with guarantee of stability and transient performance[END_REF] with T 1 = T 3 and T 2 = T 4 was used and its parameters have been tuned for the target ζ * = 7.5% with the classical technique based on the residues recalled in Section 2.1. This led to

K P SS = 49.7pu, T 1P SS = 2.3s, T 2P SS = 2.13s. ( 16
)
The latter two simulations have been carried out only with voltage regulations for the machine, i.e., its mechanical torque has been considered constant. This is the usual way of tuning the PSS since the voltage and mechanical dynamics are sufficiently decoupled (in frequency) to be considered independent. In other words, the adjustment of the voltage regulators have no, or little influence, on the turbine regulation and vice versa. As a consequence, the two classes of regulators are synthesized independently and, afterwards, used together on the machine. This is no longer valid in the situation considered in this paper. Indeed, if the frequency of the modes to be damped by the PSS is too low, the voltage correction actions may act also on the turbine and vice versa. Consider the structure in Fig. 7 which corresponds to a standard governor if the lead-lag bloc in dotted line is not used (i.e., T 1GOV = T 2GOV = 0). In this case, a phase lag action is induced by the PI bloc in Fig. 7. If the PSS is tuned for a frequency close to the governors dynamics and if it should provide a phase lead (as it is the case for the parameters ( 16)), this action is thus diminished by the PI bloc of the governor and the PSS is thus less efficient when used in conjunction with the governor. This is the case of the example treated here. If the governor presented above with R = 20pu, K P = 0.5pu, T I = 10s is used in conjunction with the PSS tuned before, the damping target is no longer achieved as shown by the response in dotted line in Fig. 6a.

Mixed control synthesis

To overcome the problem formulated in the section above, the structure of the governor has been enriched with the lead-lag filter in dashed line in Fig. 7. This allows one to shape the frequency profile of the governor without modifying the parameters obtained for the initial design and which ensure the turbine performances, i.e., the droop R, the gain K P and the time constant T I of the PI bloc. The tuning of the parameters T 1GOV and T 2GOV of this new lead-lag filter is done simultaneously and in a coordinated way with the tuning of the parameters of the PSS. They are the result of an constrained optimization of the costs function [START_REF] Kamwa | Optimization-Based Tuning and Coordination of Flexible Damping Controllers for Bulk Power Systems[END_REF] as in Sections 3.1 and 3.2. Notice, that in this process, the parameters R, K P and T I are considered as constants. For the SMIB example presented in Section 4.1, the target of this optimization is ζ * = 10% and it is reached with K P SS = 48.9pu, T 1P SS = 0.48s, T 2P SS = 0.62s, T 1GOV = 1.17s and T 2GOV = 0.26s. This can also be checked with the response of the nonlinear simulation in solid line in Fig. 6b obtained with the latter parameters.

The same situation is reproduced on the 4-machines test system in Fig. 8 of which parameters are given in Appendix 1. Indeed, the lines are NHVC1-NHVA1, NHVC1-NHVD1 and NHVC2-NHVB1 are longer than the other branches of the grid and, as a consequence, the inter-area mode related to the swing of generator GENC against the three machines in the south of the system is of low frequency (0.1Hz) and poorly damped. This mode can be observed in solid line in Fig. 9a which shows the responses to a short-circuit at bus NHVCEQ. This mode is damped when a standard PSS is added to machine GENC (dotted-line response in Fig. 9a) but this improvement is lost when a standard governor is also used for the same machine. In the latter case, the system is unstable (dashed-line response in Fig. 9a). As in the SMIB case, this problem is overcome with the coordinated synthesis of the PSS-governor for which satisfactory responses are obtained as shown in Fig. 9b (solid line).

General multi-level tuning procedure

The results of sections 2, 3 and 4.2 are now aggregated into a general control procedure in order to satisfy the three levels of coordination described in the Introduction.

General control model

The control model introduced in Section 2 is extended in order to take into account also the dynamics of the turbines of the m generators selected for the control. It consists of the transfer matrix H(s) in Fig. 10 defined by

P ω = H(s) P m V ( 17
)
where P = P 1 . . . P m T is the vector which contains the electrical powers of the involved machines, ω = ω 1 . . . ω m T the vector of the rotor speeds, P m = P m 1 . . . P mm T the vector of the mechanical powers and Using the methodology given in Section 2, a reduced order representation is constructed for the transfer matrix H(s) defined above in the range of frequency of the modes to be damped (in the set Λ).

V = V ref -V S with V ref = V ref 1 . . . V refm T the

Tuning procedure

The control model above is used to tune the parameters of the PSSs Γ i , i = 1, ..., m and of the lead-lag filters of the governors Γ i , i = 1, ..., m. The same procedure based on the constraint optimization technique given in Section 3 can be used. The main control objective is the damping of the modes directly taken into account in the cost function [START_REF] Kamwa | Optimization-Based Tuning and Coordination of Flexible Damping Controllers for Bulk Power Systems[END_REF]. Stability and robustness specifications can be taken into account by adding supplementary constraints to the optimization problem

{K * i , T * 1 i , T * 2 i , T * 1GOV i , T * 2GOV i , } = argmin{J contr } (18)
as shown in Section 3, [START_REF] Kamwa | Robust Design and Coordination of Multiple Damping Controllers Using Nonlinear Constraints Optimization[END_REF] and [START_REF] Marinescu | Robust Coordinated tuning of Parameteres of Standard Power Systems Stabilizers for Local and Global Grid Objectives[END_REF].

Large-scale test

The general methodology summarized in the preceding section has been applied to a large-scale power system. A benchmark has been constructed from the ENTSO-E model presented in Section 2.3 as follows: the dynamic model has been triplicated and the three samples have been interconnected as shown in Fig. 11. The resulting benchmark reproduces in an exploitable manner the interaction phenomena described in the Introduction. More precisely:

• the resulting power system model has a realistic large-scale. It consists of about 1200 generators which is the actual size of the European interconnected power system (which includes Romania, Bulgaria and Turkey zones).

• the series interconnection in Fig. 11 generates a low-frequency inter-area mode.

The most involved machines are among the nuclear plants of the French zone. The shape of the oscillation is zone ENTSO-E 2 against ENTSO-E 1 and 3 in Fig. 11 and the frequency of the mode is 0.09Hz. Indeed, Fig. 12 provides in dashed and dashed-dotted lines the speed responses of machines BUGEY of zones ENTSO-E 2 and 3 to a 150ms short-circuit at the connection bus of machine BUGEY of zone ENTSO-E 2. It can be noticed the phase opposition of the two oscillatory responses as well as a common period which correspond to the frequency of the mode. If the specific four-loops voltage regulators (see, e.g., [START_REF] Bénéjean | A New Voltage Regulator for Large French Alternators[END_REF] and related references) which control the French nuclear plants is replaced by a standard IEEE AVR (type IEEE X4 -Type 4 excitation system) [START_REF]IEEE Recommended Practice for Excitation System Models for Power System Stability Studies[END_REF] without PSS for the three most participant machines in the oscillation (namely BUGEY, GRAVELINE and PALUEL in zone ENTSO-E 1), the damping of the aforementioned mode drops to ζ = 3.48% (see again the two responses aforementioned). The context is thus similar to the one of [START_REF] Luther | Feasibility Aspects of a Synchronous Coupling of the IPS/UPS with the UCTE[END_REF] in which a very spread and poorly damped inter-area mode has been put into evidence at 0.07Hz.

• in opposition to study [START_REF] Luther | Feasibility Aspects of a Synchronous Coupling of the IPS/UPS with the UCTE[END_REF], this benchmark is exclusively constructed with European data.

When the target damping of the inter-area mode mentioned above is fixed at ζ ref = 10%, the control synthesis methodology given in Section 5.2 leads to the optimal parameters given in Table 4 for the three French machines chosen for the control in zone ENTSO-E 1. The obtained parameters have been implemented in the structures given by [START_REF] Okou | Decentralized multivariable voltage and speed regulator for large-scale power systems with guarantee of stability and transient performance[END_REF] and Fig. 7 on the selected machines in the nonlinear model of the interconnected benchmark in Fig. 11. The speed response of machine BUGEY of zone ENTSO-E 2 given in solid line in Fig. 12 shows that the tuned regulators are efficient, i.e., the target value of 10% of the damping of the inter-area mode is reached.

Remark

Notice that in the methodology of the mixed-control proposed in Sections 4.2 and 5.2 the governor structure is only modified by adding a lead-lag filter (Fig. 7) which is tuned on the frequency of the inter-area modes which have to be damped. It acts thus only when the latter modes are excited and its impact in these situations on the mechanical torque is limited as shown in Fig. 13a. This figure shows the response of the mechanical torque of one of the machines used for the synthesis presented for the large-scale test to the short-circuit at its grid connection point which excites the target low-frequency (0.09Hz) inter-area mode. Also, the initial function of the governor is not deteriorated by the added lead-lag filter as shown by the response in Fig. 13b to a 3% step increase of the power reference of the same machine.

Conclusions

This paper shows how several levels of coordination can be simultaneously achieved for the synthesis of PSSs. These coordinations, especially the one between the tuning of the PSSs and of the governors of the chosen generators, opens the way to stabilize very large interconnected power systems in which the frequencies of the slowest interarea modes come into the turbine frequency range. The size and the structure of the actual European interconnected power system are close to the conditions mentioned above. The next potential extensions of this system [START_REF]Feasibility Study on the Interconnection Variants for the Integration of the Baltic States to the EU Internal Electricity Market[END_REF], [START_REF] Luther | Feasibility Aspects of a Synchronous Coupling of the IPS/UPS with the UCTE[END_REF] might cause such interactions and the methodology presented here is a way to overcome these difficulties.

The principles of coordination presented here can be applied to other types of PSSs like, e.g., the multi-band PSS [START_REF]IEEE Recommended Practice for Excitation System Models for Power System Stability Studies[END_REF], [START_REF] Grondin | An Approach to PSS Design for Transient Stability Improvement through Supplementary Damping of the Common Low-Frequency[END_REF] or to coordinate the damping actions of the generators with the ones of other grid devices as, for example, some FACTS (see [START_REF] Cai | Simultaneous Coordinated Tuning of PSS and FACTS Damping Controllers in Large Power Systems[END_REF] and related references).

The key point of the proposed methodology is the use of an efficient control model for the synthesis of the control laws which allows, in particular, to deal with the large-scale systems. This model is built to take into account the interactions of the main dynamics of the system.

The parameters of the regulators are computed via a constraint optimization in order to keep standard structures for the PSSs and governors. First, the constraints of this problem can be adapted in order to ensure different aspects of performances and robustness (especially by taken into account several operating points of the system) for the resulting closed-loop as shown, for example, in [START_REF] Jabr | A Sequential Conic Programming Approach for the Coordinated and Robust Design of Power System Stabilizers[END_REF] and related references.

Next, algorithms with better performances may be used to solve the optimization problem [START_REF] Marinescu | A Unified Framework for Nonlinear Dynamic Simulation and Modal Analysis for Control of Large-Scale Power Systems[END_REF]. Indeed, the optimization methodologies have been the subject of ongoing enhancements over the years. A review of the most used in power systems is given in [START_REF] Eslami | Power System Stabilizer Design Based on Optimization Techniques[END_REF]. Finally, the optimization itself can be replaced by a synthesis technique with better performances if new structures (like, e.g., state-space forms) are allowed for the regulators instead of the actual standard ones. From this point of view, the methodology presented here is an open framework.

  vector of references of voltage control loops, V S = V S 1 . . . V Sm T the vector of the PSS outputs, Γ(s) = diag{Γ 1 (s), ..., Γ m (s)}, Γ(s) = diag{Γ 1 (s), ..., Γ m (s)} where Γ i (s) is the transfer function of the PSS of machine i, Γ i (s) is the transfer function of the governor of machine i and i ∈ {1, ..., m}.

Figure . 1 :Figure . 3 :Figure . 4 :Figure . 5 :Figure . 7 :FigureFigure . 9 :Figure . 10 :

 13457910 Figure .1: Control model

Table . 1

 . : Damping ζ[%] of the modes in ΛTable .3: Comparison of the stabilizer parameters achieved with coordinated tuning and with robust coordinated tuning Table .4: Parameters tuned for PSSs and governors of the large-scale test system Table .5: Participations and right eigenvectors' phase for the 0.07Hz mode of the large-scale system

					K * [pu]	T * 1 = T * 3 [s]	T * 2 = T * 4 [s]
		Almaraz 0.59 / 2.43 0.22 / 0.21 0.05 / 0.02
		Cofrentes 4.81 / 5.74 0.59 / 0.58 0.15 / 0.05
		PGR		0 / 0.54	0.71 / 0.61 0.15 / 0.06
					BUGEY GRAVELINE PALUEL
			K P SS [pu]	28.5	27.6		28.5
			T 1P SS [s]	2.1	2.6		2.1
			mode 1 0.23Hz mode 2 0.24Hz mode 3 0.91Hz T 2P SS [s] 2.8 2.4 2.7
	without PSSs T 1GOV [s]	3.87 2.1	11.7 2.1	2.1	6.25
		with PSSs T 2GOV [s]	8.43 2.8	9.97 2.8	2.8	11.53
		Table .2: Coefficients of the control model used in Fig. 3a
			s 6	s 5	s 4	s 3	s 2	s 1	s 0
		P(s)						0.001 -0.11
		Q(s)					1	1.55 44.20
		A(s) Zone ENTSO-2 -0.001 -0.01 -0.10 -0.30 -0.33 -0.52 Zones ENTSO-1 and 3
	Machine	B(s) 1 Participation [%] Phase [ • ] Machine 1.23 38.59 22.53 169.97 46.80 193.85 Participation [%] Phase [ • ]
	2 BUGEYSF1 70.9		1.2	3 BUGEYSF1 20.6	-175.7
	2 GRAVSF1	88.2		1.0	1 PENLYSF1 29	-165.1
	2 PALUESF1 100		-0.4	1 PALUESF1 32.1	174.9

Acknowledgment

The authors would like to thank Dr. B. Mallem for preliminary investigations on the topic of governors control at the end of his Ph.D stay at RTE.

Appendix 1: Data of the test systems SMIB test system (Fig. 5

, T e =0.1s, K e =1pu; GOVERNOR: constant mechanical power the line: R=0pu, X=0.7pu (S base =100MVA and V base =400kV) load-flow : P=100MW, Q=35MVar, V=24kV, V inf =400kV Intermediate-scale test system (Fig. 8): common generators' parameters: the same as for SMIB generators' inertias [MW s/MVA] : GENA1 10, GENB1 5, GENB2 5, GENC 10 block transformers: 24kV/380kV, rate=1.1pu, X=0.00769pu for GENA1, X=0.02251pu for GENB1, X=0.09pu for GENB2, X=0.00769pu for GENC (base 24kV) regulators: AVR of each machine: the same as for SMIB with K e =1pu; K=40pu; T e =0.1s for GENA1 and GENB1 and K=30pu; T e =1s for GENB2 and GENC GOVERNOR for GENA1, GENB1, GENB2:


Leyla Arioua
email: leyla.arioua@ens-cachan.fr

Bogdan Marinescu
email: bogdan.marinescu@rte-france.com

Large-Scale Implementation and Validation of Multivariable Control With Grid Objectives of an HVDC Link

The HVDC links are increasingly used not only to interconnect asynchronous AC systems but are also embedded into a same meshed AC power system. Thanks to its speed and flexibility, the HVDC technology is able to provide transmission system advantages as transfer capacity enhancement and power flow control. In addition, studies have shown that the way of controlling the HVDC converters impacts the stability of the AC system. This can be particularly exploited to enhance the dynamic power system performances during transients. In this paper a robust multivariable control design for HVDC link converters is proposed. It is based on the coordination of the control actions of the HVDC converters and the use of a control model which takes into account the dynamics that mostly impact stability of the neighbor zone of the HVDC link. This new methodology was used to synthesize the controller for an actual grid 1000 MW HVDC link reinforcement project called "Midi-Provence" in the southern part of the French grid. The synthesis, implementation and validation processes are presented in detail. The new controller is tested in comparison with the standard vector control. A large-scale dynamic model of the whole European power system, currently used and updated by the European TSO's for the interconnection studies has been used with Eurostag simulation software.

I. INTRODUCTION

The HVDC link is a mean of transmission of electric power based on high power electronics. Thanks to its speed and flexibility, the HVDC technology is able to provide the transmission system advantages as transfer capacity enhancement and power flow control [START_REF] Hingorani | Understanding FACTS: concepts and technology of flexible AC transmission systems[END_REF]. Initially, it was used in power systems to interconnect asynchronous AC systems. The ends of the HVDC link are electrically independent one from each other and this avoids the propagation of perturbations between the two AC grids as it is the case of the England-France interconnection [START_REF] Goodrich | The 2000 MW HVDC link between England and France[END_REF]. Nowadays, HVDC links are increasingly embedded into a same meshed AC system in order to enhance the grid's transmission capability and flexibility of the power system. In this context, the HVDC link co-exists with other AC system elements as for instance AC lines and generators. Several projects of insertion of this type of device in a meshed AC system are underway in Europe, as for instance France-Spain HVDC project [3], France-Italy interconnection and Midi-Provence HVDC project [4]. This latter is our case of study. In fact, a 1000 MW HVDC link will be inserted in France between the two areas Fos and Gaudière, in order to enhance capacity of power transmission and the stability of the Midi-Provence region. The need of such HVDC link is due to the fact that only one sigle 400 kV AC line interconnects the two areas (see Fig. 1).

It was shown that the strategy used to control the HVDC converters can impact the stability of the system in which the link is embedded [START_REF] Hauer | Robustness issues in stability control of large electric power systems[END_REF]- [START_REF] Henry | Influence of Embedded HVDC Transmission on System Security and AC Network Performance[END_REF]. This conducted us in a previous paper [START_REF] Arioua | Control of HVDC links with Overall Large-Scale Grid Objectives[END_REF], to take into account this context of use of the HVDC during the synthesis of converters controllers. In fact, we have developed a methodology which takes into account, in a control model, the neighbor zone of the HVDC link and the dynamics that most impact transient stability. This control model captures the dynamics that are relevant for the transient stability of the power system. The controller synthesis is done on the base of this control model. As a consequence, power system performances were enhanced in addition to local performances. Also in other papers [START_REF] Hu | Improvement of transient stability in AC system by HVDC Light[END_REF], [START_REF] Mao | Application of an optimal coordinated control strategy to VSC HVDC[END_REF], the transient stability is considered as an a priori objective for the synthesis of coordinated controllers for HVDC. However, the effectiveness of these methods had not been proved in the case of large-scale power systems. The present paper shows that the strategy developed in [START_REF] Arioua | Control of HVDC links with Overall Large-Scale Grid Objectives[END_REF] is applicable to a realistic large-scale case which is the full model of the European power system. Moreover, the control model is enriched with a family of perturbed models, which allowed us to take some critical fault situations into consideration during the synthesis. This lead us to use a more efficient control design (H-infinity theory). The proposed control is tested in comparison with the standard vector control via simulations performed with Eurostag [START_REF] Meyer | EUROSTAG, a single tool for power system simulation[END_REF]. The rest of the paper is organized as follows: in Section II, the context of this work is described. Then in Section III, the method for developing the control model for the power system is presented. This control model is used in Section IV to synthesize a robust multivariable outputfeedback controller for the HVDC converters. Section V deals with the computation and implementation of the controller. Finally, in Section VI, simulation results are presented to illustrate the performances and the robustness of the new controller.

II. CONTEXT OF THE STUDY

Midi-Provence project aims to strengthen the link between the two dynamic regions, Provence-Alpes-Côte d'Azur and Languedoc-Roussillon by inserting an HVDC link of 1000 MW capacity between them. This two regions are now directly connected by a single 400 kV AC line going from Tavel (Est of Gard) to La Gaudière (in the Aude) (see Fig. 1). In case of breakdown of this link, the electric power security of much of the south of France would be weakened. In some worst scenarios, this situation could lead to a major power cut in the south of France, and even in-beyond Europe, to balance supply and demand. This HVDC project responds thus to a need of enhancement of both stability and transport capacity in this region [4].

III. CONTROL MODEL SYNTHESIS

The methodology developed in [START_REF] Arioua | Control of HVDC links with Overall Large-Scale Grid Objectives[END_REF] is briefly recalled here and enriched in order to take into account not only the nominal model of the power system, but also a family of perturbed models which correspond to the critical faults. The resulting control model must preserve the dynamics that most impact the stability of the region neighbouring the HVDC link. For this reason, it was necessary to start from a full model of a sufficient large zone arround the HVDC link. In this study, we considered the whole European power system.

A. Description of the European power system

We use for the European power electric system, the reference model which is maintained and generally used for interconnection studies by the European TSO's. It is a detailed non-linear model including generators along with their regulations (AVRs, PSS and Governors) and where the high voltage network (225, 400kV ) is modeled. It consists of 1121 generators, 7625 nodes, 10404 lines, 2550 transformers and 458 GVA global apparent power. The Midi-Provence HVDC link is of 230 km length and 1000 MW power capacity. It is placed as shown in Fig. 1. Note that we have considered a 2020 forecast winter situation.

B. Structure of the control model

The control model is composed of two parts: the AC power system (grid and machines) and the HVDC link. In the method developed, we retain in the control model the most critical machines, i.e., those which cause the system loss of synchronism after a severe perturbation. A transient stability margin of the power system is estimated by the Critical Clearing Time (CCT) which is defined as the maximal fault duration for which the system remains transiently stable [START_REF] Kundur | Power stability and control[END_REF]. The instability is then manifested by the loss of synchronism of a group of machines. The most critical of them are retained in the control model. The irrelevant topology of the electric network is also reduced.

C. Construction of the control model

1) The study area selection: The study area is the zone which is impacted by the HVDC link and corresponds to the neighbor zone of the HVDC interconnection of which transient stability depends on the HVDC dynamics (see Fig. 1). The extent of this zone is detected by standard stability studies usually run by TSOs.

2) Critical machines selection: The strategy of machines selection is based on the computation of CCT at well-chosen points of the study area. As the study area, these points are determined by the same type of a priori stability studies run by TSO's. The method is more detailed in our previous work [START_REF] Arioua | Control of HVDC links with Overall Large-Scale Grid Objectives[END_REF]. We have retained in our case 14 machines in the study area defined above.

3) Reduction of the topology: Reduction of the topology of the power system is also performed. The buses to which the retained machines are connected are kept along with the neighbor zone buses. The rest of the buses and branches are replaced by a WARD-PV method [START_REF] Baldwin | Dynamic Ward equivalents for transient stability analysis[END_REF], i.e., by equivalent impedances and injectors.

4) Perturbed models: The model described above can be considered as nominal in the sense that it corresponds to the no incident situation of the power system. A family of perturbed models was also generated. Each member of this family corresponds to a model under a critical fault situation. The electric system is then considered as a set of perturbed systems. The most severe fault generates the the worst-case situation (see Section IV-C).

D. Modeling of the VSC-HVDC

A voltage-source converters (VSC) based HVDC is considered as in Fig. 2a. The converters are VSC employing IGBT power semiconductors, one operating as a rectifier and the other as an inverter. The two converters are connected through a DC cable. These converters have the ability to rapidly control the transmitted active power, and also to independently exchange reactive power with the AC system at each end. The high frequency switching operation of the power electronics is neglected and each converter can be considered as an ideal sinusoidal voltage source whose magnitude U c and phase angle θ can be controlled (see [START_REF] Latorre | Active and reactive power control of a VSC-HVDC[END_REF]). The VSC-HVDC is thus modelled as two sources, each one in series with the converter transformer (see Fig. 2b where Ūc1 = U c1i + jU c1j , Ūc2 = U c2i + jU c2j and θ 1 , θ 2 are, respectively, the magnitudes and the phase angles of the voltage sources, P 1 , P 2 , and, Q 1 , Q 2 are, respectively, active and reactive powers exchanged with the power system and X T 1 , X T 2 are the transformer reactances). In this model of HVDC, it is assumed that the DC voltage is kept close to its rated voltage, because its dynamic is much faster than the one of currents. Therefore, the losses of the converters are assumed constant, regardless of the current through the converters. They are represented as a constant active load. However, the losses of the DC cables are neglected [START_REF] Latorre | Active and reactive power control of a VSC-HVDC[END_REF]. This leads to

P 1 = (U 1 (sin θ 1 U c1i -cos θ 1 U c1j ))/X T 1 Q 1 = (U 2 1 -U 1 (cos θ 1 U c1i -sin θ 1 U c1j ))/X T 1 P 2 = (U 2 (sin θ 2 U c2i -cos θ 2 U c2j ))/X T 2 Q 2 = (U 2 2 -U 2 (cos θ 2 U c2i -sin θ 2 U c2j ))/X T 2 .
(1) In addition, notice that P 2 = -P 1 . The control variables can be expressed as follows:

U c1i = U c1i0 + ∆U c1i , U c1j = U c1j0 + ∆U c1j , U c2i = U c2i0 + ∆U c2i , U c2j = U c2j0 + ∆U c2j (2)

E. Validation of the control model

The resulting control model is validated by non-linear simulations comparing responses of the reduced and the full model when a fault occurred near the HVDC link. Theses simulations were done using Eurostag software. Fig. 3a and Fig. 3b for instance, illustrate the speed responses of two generators retained in the control model G 1 and G 2 . We can observe that the two responses of the control model and the full model are closed to each other for this two generators.

IV. SYNTHESIS OF THE CONTROL LAW

In our previous work [START_REF] Arioua | Control of HVDC links with Overall Large-Scale Grid Objectives[END_REF], an output-feedback controller for HVDC converters has been developed using an LQ pole placement [START_REF] Arioua | Control of HVDC links with Overall Large-Scale Grid Objectives[END_REF]. In this present work an H-infinity control theory is used to take into account not only the nominal model of the system, but also the family of disturbed models presented in Section III-C4. This allowed us to improve both performances and robustness which allows a more robust control design.

A. Design Requirements

The synthesis of the HVDC controller is done such that the performances are ensured for several cases of fault (i.e, the enhancement of transient stability) along with local ones which are tracking of references for active power, reactive power. For these purpose, a robust multivariable control is applied. To avoid the use of remote variables (e.g., machines speed, angles, ...), an output feedback structure described in Fig. 4 is used. Only variables available at converter stations are used and a coordinated controller is developed for both converters. A linear approximation of the control model around an equilibrium point is used for the synthesis.

B. Linear approximation of the control model

The control model includes the machines retained in the control model and their regulations (AVRs and governors) be represented by the non linear differential and algebraic equations

ẋ = f (x, U) y = g(x, U) (3
U=[U c1i , U c1j , U c2i , U c2j ] y=[P 1 , Q 1 , Q 2 ]. y ref =[P ref 1 , Q ref 1 , Q ref 2 ].
The linearization of the system (3) around a given equilibrium The model of the European power system is of a very large scale (17020 sate variables and 43940 total variables). Thank's to the control model developed in Section III, we can reduce the order of the model considered for the control to 310 state variables with a total of 1117 variables. By performing a minimal realization using Robust control toolbox of Matlab, a 30th-order control model is obtained. Notice that dedicated tools (like, e.g, Eurostag [START_REF] Meyer | EUROSTAG, a single tool for power system simulation[END_REF]), provide the linear model (4) for large-scale power systems.

C. Mixed sensitivity H-infinity control design

Let P (s) be the nominal design model described in (4). Let P ∆ (s) be the family of perturbed linearized models as described in section III-C4. This set can be modeled as an uncertain model P ∆ (s) = (1 + ∆(s)W t (s))P (s) as shown in Fig. 5 where W t (s) is the uncertainty weight and ∆(s) the normalized model deviation with respect to the nominal model. Fig. 6 shows the singular values (frequency responses for multivariable systems) of the deviation of all the perturbed models with respect to the nominal one and the weight W t (s) which is used to capture how the relative uncertainty varies with frequency. This weight function allowed us to take into account the worst-case situation in the control model.

In the control theory, it exists methods which permit to consider the worst-situation described above, during the synthesis of control law. Among theses methods, the Mixed sensitivity H-infinity control design. The sensitivity S(s) = (I + P (s)K(s)) -1 is the closed loop transfer function matrix from the reference inputs to the error. The complementary sensitivity T (s) = P (s)K(s)(I + P (s)K(s)) -1 is the closed-loop transfer function matrix from the references to the measured outputs where P is the nominal model, K the controller, and I identity matrix. The principle of this control design is that if it exists a controller K(s) given by its transfer matrix which ensures W t (s)T (s) ∞ ≤ 1 for all frequencies, the closed-loop system is ensured stable for all plants in P ∆ (see [START_REF] Skogestad | Multivariable feedback control[END_REF]). The desired performances are achieved if W s (s)S(s) ∞ ≤ 1. W s (s) and W t (s) are the weights used to shape the sensitivity Block diagram of the generalized plant in the mixed H-infinity synthesis [START_REF] Skogestad | Multivariable feedback control[END_REF] and the complementary sensitivity respectively. The weight W t (s) is used to ensure robustness against model uncertainties. The weight W s (s) is based on shaping the sensitivity for desired tracking performances. The state space model is augmented with the matrix weights W t (s) and W s (s) as shown in Fig. 7.

The objective of the mixed sensitivity problem is to find a stabilizing controller K which minimizes.

W s (s)S(s) W t (s)T (s) ∞ (5)
In this work W t (s) is the uncertainty weight and W s (s) is chosen in the form (see [START_REF] Zhou | Essentials of robust control[END_REF])

W s (s) = s/ k √ M s + ω b s + ω b k √ ε k (6)
where ω b is the bandwidth, M s the peak sensitivity , ε the steady state error desired with respect to a step input, and k is the order of the weighting function. Notice that making ε small enforce the integral action (in our case ε = 10 -4 )

V. COMPUTATION AND IMPLEMENTATION OF THE CONTROLLER The multivariable robust controller K is calculated by solving the Ricatti equation (see [START_REF] Zhou | Essentials of robust control[END_REF]) via a function of Robust control toolbox of Matlab. The final controller is given in its state representation (A K , B K , C K , D K ) and it is of the same order as the plant model augmented by the weights W s and W t . In our case, the model obtained after reduction and minimal realization is of 30th-order. The controller obtained is of 40th-order which is too high for the implementation. For this reason, it is a common practice to reduce the order of the controller. A reduction method based on the singular values [START_REF] Zhou | Essentials of robust control[END_REF] helped us to generate a controller of 5th-order. The truncations has been performed with some existing functions in Robust control design toolbox of Matlab. This multivariable output controller can be written in its matrix transfer function form K(s) = C K (sI -A K ) -1 B K where A K is a state matrix of 5th order.

A coordinated structure

The controller transfer matrix obtained above can be detailed

∆U 1 ∆U 2 = K 11 (s) K 12 (s) K 21 (s) K 22 (s) K(s) ∆y 1 ∆y 2 (7)
∆U 1 = ∆U c1i ∆U c1j , ∆U 2 = ∆U c2i ∆U c2j (8)
and

∆y 1 = ∆P 1 ∆Q 1 , ∆y 2 = ( ∆Q 2 ) . (9)
The transfer matrices K 12 (s) and K 21 (s) are the cross-terms of the controller and ensure the coordination between the control of the two HVDC converters. In the standard vector control, HVDC links are controlled with a decentralized structure, which is equivalent to consider a priori that the two HVDC converters don't interact (K 12 (s) = 0 and K 21 (s) = 0). This hypothesis is not valid when the HVDC is used in the new grid context mentioned in the introduction. Each term of the K ij transfer matrix is in the form of :

k ij = α ij s 4 n=1 1 + τ n s 1 + γ n s (10)
VI. SIMULATION TESTS This section illustrates the effectiveness of the H-infinity controller proposed. Simulation tests using Eurostag software [START_REF] Meyer | EUROSTAG, a single tool for power system simulation[END_REF] were performed. The new synthesized controller is tested on the full model introduced in Section III-A. An export of 1000 MW is considered from Fos to Gaudière via the DC link. The robust H-infinity controller is compared with the standard vector controller (see for example [START_REF] Li | Control of HVDC Light System Using Conventional and Direct Current Vector Control Approaches[END_REF]) which has a cascade structure with an inner loop more rapid than the outer one. The proportional and integral gains of this control are synthesized using standard criteria for electrical drives. The two controllers are tuned to satisfy almost the same performance specifications (the usual time setting for HVDC voltage and power control).

A. Transient stability

In a first scenario, a symmetrical double fault at Tavel-Tamareau, cleared after 100 ms by opening the lines at their both ends is simulated. Fig. 9a and Fig. 9b show the responses of active and reactive powers exported. It can be observed that the dynamic responses for this case of fault are better with the new controller especially for the reactive power. Indeed, more reactive power is needed with the vector control. Moreover, a limitation is reached in this case during the transient. In a second scenario, a symmetrical short-circuit is applied at one of the HVDC terminals (in the Fos area). As for the previous case the transient responses are enhanced in comparison to the standard vector control. This is confirmed by analysis of the transient stability of the zone. In fact, Table I shows that the critical clearing time are augmented when the new control is used. We have considered the situation of a 450 MW import of active power from Gaudière area to Fos. As for the previous tests, a 100 ms symmetrical fault is applied at one of the HVDC terminals (Gaudière area). The active and reactive powers imported are represented in Fig. 10a and Fig. 10b and both are comparable with the ones obtained when the export scenario is considered (Fig. 9a and Fig. 9b). This confirm the good robustness of the performances of the proposed controller against variation of operating conditions.

VII. CONCLUSION

A large-scale implementation and validation of a multivariable controller of an HVDC link has been presented. In fact, a new methodology for developing the control law for HVDC link converters has been used. It is based on the use of a control model which takes into account the region of the The main advantage of this new methodology is that the main dynamics of the full power system are taken into account via a control model. It has been shown that this strategy can be applied to a realistic large scale power system.

Moreover, an H-infinity control design has been used. It allowed us to consider at the synthesis stage some faulted situations and this contributed to enhance robustness along with performances of the HVDC link controller and a better fault recovery.

The control model proposed in this methodology is a way to connect the system expertise to the control one. More specially, in more cases, the manufactures do not possess the knowledge and the data to construct such a model and, as a consequence, the control of the DC link is synthesized with a very simplified model of the surrounding system. Most of the time, infinite buses or large power generators are considered connected to the two ends. The control model proposed here is rich enough to improve the performances of the controllers and simple enough to be built and provided by TSO's to the manufactures.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Midi-Provence interconnected power system

Fig. 4 .

 4 Fig. 3. Speed responses to a 100 ms short-circuit in the two cases of full and reduced models (p.u)

Fig. 5 .

 5 Fig. 5. Plant with multiplicative output uncertainty

Fig. 6 .

 6 Fig. 6. Relative gaps vs. Magnitude of Wt(s)

Fig. 8 .

 8 Fig. 8. Responses to a 100 ms Tavel-Tamareau short-circuit

Fig. 9 .

 9 Fig. 9. Responses to a 100 ms short-circuit near Fos area in case of 1000 MW power export

Fig. 10 .

 10 Fig. 10. Responses to a 100 ms short-circuit near Gaudière area in case of 450 MW power import

TABLE I .

 I CRITICAL CLEARING TIME VALIDATION OF THE TWO

		CONTROLLERS
	Fault	Standard	New
	simulated	controller	controller
		CCT [ms]	CCT [ms]
	f ault1	216	224
	f ault2	192	204
	f ault3	183	195
	f ault4	181	192
	f ault5	194	205

th Power Systems Computation Conference Wroclaw, Poland -August 18-22, 2014 th Power Systems Computation Conference Wroclaw, Poland -August 18-22, 2014 th Power Systems Computation Conference Wroclaw, Poland -August 18-22, 2014 th Power Systems Computation Conference Wroclaw, Poland -August 18-22, 2014 th Power Systems Computation Conference Wroclaw, Poland -August 18-22, 2014 th Power Systems Computation Conference Wroclaw, Poland -August 18-22, 2014 ACKNOWLEDGMENT The authors would like to thank Prof. Eric Monmasson, University of Cergy-Pontoise and Dr. Alexandre Parisot, head of Integration of New Technology division (INT) of RTE for their helpful suggestions and remarks.


M Belhocine
email: mohamed.belhocine@ens-cachan.fr

Member B Marinescu
email: bogdan.marinescu@rte-france.com

A mix balanced-modal truncations for power systems model reduction

In this paper, a new model reduction technique mixing the balanced and modal truncations is proposed for power systems. Usually, only the power transmission lines are approximated such that, when used with the rest of the power system, the phenomena of interest are reproduced. Our approach is structural since it is based on the reduction techniques of the dynamical systems and consider the full power system rather than the lines alone. It is shown that, to obtain a good approximation with the balanced or modal truncations, a large number of dynamics are needed to be kept in the reduced model. By mixing the balanced and modal truncations, a much lower order reduced model is obtained and the behaviour of the full system is faithfully reproduced. This result is useful in practice since it opens the way to improve the numerical simulations and behaviour analysis of power systems by making easy the digital implementations. This new technique is checked by simulation tests in Matlab and EMTP.

I. INTRODUCTION

A power transmission system is a set of different electric components interconnected by electric transmission lines. Indeed, the full dynamic model of such a system is usually constructed by putting together the dynamic models of each system component with those of transmission lines. Generally, the electric components like machines and loads are finite dimensional systems described by ordinary differential equations (ODE) unlike the transmission lines which are distributed parameters (DP) systems described by partial differential equations (PDE). Therefore, the full model (i.e., all interconnected models) is an infinite dimensional system which can reproduce responses to all kinds of perturbations, particularly electromagnetic transients or wave propagation through the lines.

In practice, this full and detailed model is difficult to be used, and it is usually replaced by one which uses a priori simplified representation of the lines like, for example, the π-model [START_REF] Kundur | Power System Stability and Control[END_REF]. However, since the approximation of the line models is a priori imposed and not issued from the analysis of the full model with respect to the purposes of the modelization process (i.e., the phenomena which are supposed to be reproduced with the resulting simplified model), there is no guarantee on the adequacy of the resulting simplified power system model. For this reason, a structural truncation (i.e., based on the analysis of the structure of the whole system) of the full model by using techniques like balanced [START_REF] Glover | Realisation and Approximation of Linear Infinite-Dimensional Systems with Error Bounds[END_REF], [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF] or modal truncations [START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF] is a good manner to preserve the phenomena of interest by keeping only the most important dynamics in the reduced model.

The infinite dimensional dynamic models of the transmission power systems have some particularities and the direct application of the balanced or modal truncations give not satisfactory results. This is due to the infinite bandwidth of such systems (resonance peaks infinitely repeated). Physically, this is related directly to the wave propagation phenomena and means that the signals of all frequencies can propagate through the lines. As a consequence, the Hankel singular values involved in balanced truncation decay slowly and there is no significant gap between a previous and next value. This means that all modes are important from an energetic point of view and the balanced truncation give a good result only if a large number of modes are retained in the reduced model. The same conclusion comes out also for modal truncation.

To overcome this difficulty, in [START_REF] Lefebvre | Reduced order transmission line models for power system analysis[END_REF] it is proposed a balanced truncation but only for a specific, a priori fixed, range of frequency. This is done by weighting the initial dynamic model which represents the physical power system by artificial lowpass blocs which leads to a gap in the Hankel singular values. However, this way of doing is not compliant with our objectives since, as mentioned before, the resulting reduced model must be used for the analysis of dynamic phenomena of which connection to the frequency band restrict is not obvious.

In this paper, we propose an alternative method to obtain much lower order reduced models for power systems by mixing the balanced and modal truncations. More precisely, only some both strongly controllable and observable modes (i.e., more important form the energetic point of view) are kept in the modal truncation.

In addition to the low order of the resulting reduced system, this method is directly related to the structure of the whole system and not only to the one of the transmission lines. As a consequence, the resulting reduced models are better suited to reproduce structural phenomena like interarea and hypo-synchronous oscillations. Also, the adequacy of the truncation to the reduction objectives is quantitatively measured a priori (i.e., before the truncation in terms of magnitude of Hankel singular values), and not a posteriori like for the existing methods for which a simplified model is first proposed by expert intuition for the line model. Next, the phenomena of the resulting reduced model of the whole power system are checked a posteriori which leads to an uncertain trial and error process. The power system in Fig. 1 is composed by three subsystems: A generator assumed to be a voltage source V 0 (t) behind an impedance Z 0 (which is the usual simplified model used for voltages studies (see, e.g., [START_REF] Kundur | Power System Stability and Control[END_REF])), a transmission line of length for which the transverse conductance is neglected, see, e.g., [START_REF] Miano | Transmission Lines and Lumped Circuits[END_REF] and a resistive load Z L .

II. POWER SYSTEM MODEL DESCRIPTION

The dynamic model of the full system is thus an interconnection between the dynamic models of the elements mentioned above.

A. Distributed parameters model for power systems

First, the distributed parameters model of a transmission line can be written as ∂v(x, t) ∂x

= -L ∂i(x, t) ∂t -Ri(x, t), ∂i(x, t) ∂x = -C ∂v(x, t) ∂t , (1
) for 0 ≤ x ≤ , where v(x, t) and i(x, t) are respectively the voltage and current along the line, and the positive constant parameters (i.e., independent of frequency) R, L, C are, respectively, the resistance, the inductance and the capacitance per unit length. Next, the interconnection between the transmission line, the load and the generator in Fig. 1 is performed by using the Kirchhoff's laws and the boundary variables of (1) as follows

v (x, t) | x=0 = V 0 (t) -Z 0 i (x, t) | x=0 v (x, t) | x= = Zi (x, t) | x= (2)
Thus, equations ( 1)-( 2) form the full model of the power system (in Fig. 1) which can be written also under the form ż1 ż2

ż(t) = 0 -1 C ∂ ∂x -1 L ∂ ∂x -R L A z 1 z 2 z(t) Bz (t) = V 0 (t) , Cz (t) = y (t) , z (0, x) = z 0 (x) , (3)
where z 1 (t) = v (x, t), z 2 (t) = i (x, t) and with

B [z 1 (t) z 2 (t)] T = v (0, t) -Z 0 i (0, t) = V 0 (t) the input and C [z 1 (t) z 2 (t)] T = v ( , t) = y (t) , the output.
Equation ( 3) is a general form of an abstract boundary control system, where A, B, C are operators defined on a given Hilbert space Z (see, e.g., [START_REF] Curtain | An Introduction to Infinite Dimensional Linear Systems Theory[END_REF]). It can be written under the infinite dimensional general state representation as

ż (t) = Ãz (t) + BV 0 (t) , y (t) = Cz (t) , z (0, x) = z 0 (x) (4
) where à generates an analytic C 0 -semigroup T (t) = e Ãt on the space Z and the operator B is deduced from (3) (see, e.g., [START_REF] Curtain | An Introduction to Infinite Dimensional Linear Systems Theory[END_REF] or [START_REF] Tucsnak | Observation and Control for Operator Semigroups[END_REF] for more details and examples).

B. Ad-hoc finite dimensional model of the power system Fig. 2. the π-model of a transmission line As previously mentioned, in practice, before aggregating the model (3) of the full power system, the transmission line model ( 1) is replaced by a finite dimensional equivalent like the π-model in Fig. 2. Hence, the dynamic model of the full system is given by

ẋ1 ẋ2 = 0 1 -2R +2Z L LC 2 Z L -2L+RC Z L LC Z L x1 x2 + 0 2 LC 2 Ve (t) (5)
with x 1 = V s (t) and x 2 = ẋ1 , instead of (3).

Obviously, the model ( 5) cannot be used instead of (3) to reproduce all physical phenomena of interest in power system analysis. It is generally used only to perform loadflow computations, to analyse transient stability and to design voltage controllers.

III. STRUCTURAL TRUNCATION OF THE POWER SYSTEM

MODEL

A structural approximation of the power system described by the distributed parameters model [START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF], is a way to find a finite dimensional reduced model where only the important dynamics (i.e., which have a major impact on the behaviour of the full power system) are kept. In the balanced truncation, e.g., [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF] these important dynamics are selected from an energetic point of view based on the Hankel singular values of the input/output system. In the modal truncation, (e.g., [START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF]) they are selected following their impact on the infinite norm of the input/output system. Both techniques are recalled below and their performances are compared when applied to the simple test system in Fig. 1.

A. Modal truncation

First, consider the model (4) for the particular case of an open line and ideal voltage source (i.e., the load Z L = ∞ and Z 0 = 0 in Fig. 1). Obviously, this case does not occur in practice but it is useful to give a first view on the dynamic of the system. As the system becomes simpler in this case, analytic expressions of the transfer function as well as of the trajectories of the system are available.

More precisely, V ( , s)

V0 (s) = 1 ∞ n=1 4 2 LC (2n-1) 2 π 2 s 2 + ( R L )4 2 LC (2n-1) 2 π 2 s + 1 (6)
is the transfer function between the input voltage of sending end V 0 (t) and the output voltage at distance from the sending end v ( , t).

At this stage, a way to obtain a finite dimensional approximation of ( 4) is the modal truncation of the transfer function [START_REF] Glover | Realisation and Approximation of Linear Infinite-Dimensional Systems with Error Bounds[END_REF]. For n = 1, the poles are

p 1,2 = - 1 2   R L ± R L 2 - π 2 8 8 2 LC   . (7)
Also, the transfer function of the π-model (5) for Z L = ∞ and Z 0 = 0, has the poles

λ 1,2 = - 1 2   R L ± R L 2 - 8 LC 2   . ( 8
)
The only difference between ( 7) and ( 8) is factor π 2 8 . As the latter factor is close to unity, the first order modal truncation of ( 6) is close to the π-model (5) from a structural point of view (see [START_REF] Belhocine | On the connection between used models of power transmission lines[END_REF] for details). For trajectories, the Fig. 3 shows a comparison of the step responses of the full power system in Fig. 1 for which a π, respectively, a DP model of line is used. 1 The global shapes of the two trajectories in Fig. 3 are similar. The differences come from the transient dynamics due to the interactions between the waves during the propagation and they are observed in the upper and lower parts of the trajectories provided by the DP model. Thus, the πmodel (first order modal truncation) captures an important component in the behaviour of the DP one but not the dynamics related to the waves propagation. In the latter case, this truncation at the first order is not sufficient and a large number of fractions is needed to be kept in [START_REF] Glover | Realisation and Approximation of Linear Infinite-Dimensional Systems with Error Bounds[END_REF].

1) General approach: Starting from a general state representation (4) involving operators, the partial fractions expansion

G (s) = ∞ k=0 r k s -β k (9)
1 Notice that, to get the trajectories of the DP model, the finite difference method over a finite grid with a step of 10 -1 unit length in [0 ] was used.

of the transfer function can be obtained ( [START_REF] Curtain | Transfer functions of distributed parameter systems: A tutorial[END_REF], [START_REF] Curtain | An Introduction to Infinite Dimensional Linear Systems Theory[END_REF], [START_REF] Tucsnak | Observation and Control for Operator Semigroups[END_REF]), where β k are the poles and r k their residues. The modal truncation consists in keeping only a finite number of these fractions, i.e., the N th order modal truncation of G (s) is

G N (s) = N -1 k=0 r k s -β k .
One has the following bound for the H ∞ truncation error:

G (s) -G N (s) ∞ ≤ ∞ k=N |r k | |Re (β k ) | . ( 10
)
The latter bound shows that the partial fractions in (9) of large H ∞ norm are important to be kept in the reduced model in order to lower the truncation error.

2) Numerical approximation: For the general case of the system in Fig. 1 (i.e, with a load) and for power systems of large scale, it is not possible to analytically obtain the partialfraction expansion of the transfer function [START_REF] Miano | Transmission Lines and Lumped Circuits[END_REF]. In this case, these quantities are numerically approximated by finding an approximation (A n , B n , C n ) of the operators Ã, B, C in (4) where n is the order of the numerical approximation scheme. Now, consider the system (4) with the flowing parameters R = 0.5Ω/Km, L = 2.10 It can also be seen from Table I, that only one pole (i.e., -5.62.10 2 ) of the π-model is close to a pole (-5.48.10 2 ) of the DP model. Its rank (k = 591) is also very high in the partial fractions expansion of the DP model. Therefore, the reduction link explained above between the π and the DP model for open lines, is not true in this case. Thus, from a systemic point of view, when only the transmission lines are replaced by the π-model of power systems in Fig. 1, the resulting simplified model of the power system cannot be viewed as a structural truncation of the full one. However, the mode of null frequency (i.e., real) kept by the π-model plays an important role in the behaviour of the DP model as shown in Fig. 6 and it is also kept in the reduced model (of order 981) issued from the modal truncation. In addition, the modes judged important in the DP model according to the H ∞ norm are in high frequencies (i.e., pole with large imaginary part) as shown in Table I. Thus the modes of low frequencies (small imaginary parts) are of high rank, and to keep them (like the pole -5.48.10 2 ) in the reduced model, a large number of modes should be kept in the end.

It is obvious from the results exhibited in this section that the modal truncation does not give satisfactory results from a practical point of view, since to obtain a good approximation with this method, a large number of modes should be kept which leads to a high order reduced model which cannot be exploited in practice.

B. Balanced truncation

The balanced truncation is an energy based truncation for which the Hankel singular values are used to keep into the reduced model the highest energy states. As the latter correspond to the most controllable and observable modes, the balanced truncation is performed in two steps. First, get the balanced realization ( 11)- [START_REF] Singler | Balanced POD for model reduction of linear PDE systems: convergence theory[END_REF]. Next, truncate the system S b Ãb , Bb , C b at order r by inspecting the Hankel singular values, i.e., r such that σ r σ r+1 . The balanced realisation of the model ( 4) is given by (see, e.g., [START_REF] Curtain | An Introduction to Infinite Dimensional Linear Systems Theory[END_REF], [START_REF] Singler | Model Reduction of Linear PDE Systems: A Continuous Time Eigensystem Realization Algorithm[END_REF] or [START_REF] Glover | Realisation and Approximation of Linear Infinite-Dimensional Systems with Error Bounds[END_REF])

Ãb ij = √ σiσj σ 2 i -σ 2 j σ j f i (0) T f j (0) -σ i g i (0) T g j (0) , i = j, Ãb ii = -1 2 f i (0) T f i (0) = -1 2 g i (0) T g i (0) , ( 11
)
Bb = [ √ σ 1 f 1 (0) , √ σ 2 f 2 (0) , . . .] T , ( 12
)
C b = [ √ σ 1 g 1 (0) , √ σ 2 g 2 (0) , . . .] . ( 13
)
where, σ 1,...,∞ are the Hankel singular values of the system, f 1,...,∞ (t), g 1,...,∞ (t), are the singular vectors of the Hankel operator [START_REF] Tucsnak | Observation and Control for Operator Semigroups[END_REF]. From ( 11)-( 13) it is clear that σ i , f i , g i are firstly needed to balance the system and finally truncate it. In the literature (see, e.g., [START_REF] Curtain | An Introduction to Infinite Dimensional Linear Systems Theory[END_REF]), these quantities are well defined and their definitions are recalled bellow.

Γ (•) = ∞ 0 CT (t) B (•) dt, (14)
2) singular vectors: The singular vectors f i , g i are called Schmidt pair of Γ and satisfy

Γf i (t) = σ i g i (t) , Γ * g i (t) = σ i f i (t) . ( 15
)
Notice that in the case of infinite dimensional systems f i (t) , g i (t) are functions and not vectors as in the finite dimensional case.

Since the σ i and the Schmidt pair (f i , g i ) can be computed in practice only for a few cases [START_REF] Opmeer | Model reduction for distributed parameter systems: A functional analytic view[END_REF], the algorithms proposed e.g., in [START_REF] Singler | Model Reduction of Linear PDE Systems: A Continuous Time Eigensystem Realization Algorithm[END_REF] and [START_REF] Saragih | Model reduction for infinite dimensional systems using reciprocal transformation[END_REF] give a numerical approximation of these quantities as in Section III-A.2 and the r th order truncation is performed when σ r σ r+1 . Notice also that other methods like [START_REF] Singler | Balanced POD for model reduction of linear PDE systems: convergence theory[END_REF] can be used for finite dimensional approximation.

Consider the system (4) with the same parameters as in Section III-A. Fig. 4 shows the Hankel singular values of (4) obtained by finite dimensional approximation as in Section III-A.2 with n = 1000. It is clear form this figure that the Hankel singular values decay slowly and there is no significant gap indicating the order of truncation r (i.e., like mentioned above) but just two small gaps (at state 428 the gap is of 2, 9.10 -3 between the previous and next state and of 1, 2.10 -3 at state 868) which means, that from an energetic point of view, all modes play an important role in the inputoutput transfer or equivalently all modes are both strongly controllable and observable which means that the system it is hard to truncate. The result of Fig. 4 can be explained also from the gain Bode diagram in Fig. 5. Indeed, the system has an infinite bandwidth and the resonance peaks caused by the infinite number of complex conjugate poles have approximatively the same magnitude and they are infinitely repeated. As a consequence, the system is sensitive to input in any frequency range, i.e., any mode when excited (by putting the frequency of input exactly the same with the frequency of the mode) produces the same input-output gain, which means that all modes are important to ensure the bandwidth of the system. Indeed, the response of the reduced model computed with r = 428 which corresponds to the first gap in Fig. 4, are given in Figs. 5 and6. The frequency response in Fig. 5 shows a bad approximation of the modes of the full model especially for the low frequencies. The step responses presented in Fig. 6 show that the reduced model approximates the full one just at the start and thereafter the trajectories become very different from the ones provided by the DP model and the model implemented in EMTP (Electromagnetic Transient Program) simulator [START_REF] Dommel | Digital Computer Solution of Electromagnetic Transients in Single and Multiphase Networks[END_REF].

When r = 868 which corresponds to the second gap in Fig. 4, the resulting reduced model is much better but contain a too large number of modes. Finally, as for the modal truncation, the number of dynamics to be kept in the reduced model is large and to overcome this difficulty, a new alternative technique is proposed in the next section.

IV. NEW TRUNCATION METHODOLOGY: MIX BALANCED AND MODAL TRUNCATION

As none of the balanced or modal truncations can give a low order reduced model for the power system, the idea is to take the advantages of each technique by mixing them. Indeed, the main advantage of the modal truncation is that the modes are exactly kept but their choice is not always adapted to capture the main dynamics in the system. For the balanced truncation the choice of the dynamics to be kept is useful since it is based on the energy criteria which are close to the physic but the truncation when there is not a significant gap between consecutive Hankel singular values leads to a bad approximation of the modes. The new proposed methodology is a modal truncation for which the modes are chosen from an energetic point of view, i.e., only the modes which participate actively in the states which correspond to large Hankel singular values are kept.

A key point of this approach is the link between the states and the modes of the system which is established using the participation factors [START_REF] Verghese | Selective modal analysis with applications to electric power systems, Part II: The dynamic stability problem[END_REF]. More precisely, let Thus, each row i of the matrix P contains the net participations of all modes in the i th state.

(A n , B n , C n ) be
Here, the participation factors are used to select the modes to be kept in the modal truncation, i.e., from a state for which the Hankel singular value is large, the participation factors are computed to identify the modes which have large participation in this state and finally keep them in the modal truncation.

Hence, the steps of the methodology we propose are: 1) Compute a balanced realisation as mentioned in Section III-B. 2) Select the states for which the Hankel singular values are large. 3) Compute the participation factors of each mode into the states selected at step 2. 4) Select the modes with major participation in the states selected in step 2 by using the participation factors. 5) Construct the reduced model by modal truncation to the modes found at step 4. As previously, the system (4) is considered with the parameters given in Section III-A. The reduced model obtained by using this new technique is of order 171, which is much lower than the reduced models obtained by modal (order 981) or balanced (order 868) truncations. Fig. 7 shows a comparison between the step responses of the new reduced order model, the DP model and the π-model. The new reduced model, reproduces faithfully the behaviour of the full model which means that all dynamics of interest are kept. Also, Table II shows that the most important (from an energetic point of view) modes in this new reduced model have low frequencies unlike to the modes of Table I given by the modal truncation. This can be seen more precisely from the Bode diagrams of Fig. 8 which shows that the V. CONCLUSION This paper presents a structural approximation of power systems by a new reduction methodology which mixes the balanced and modal truncations.

From the systemic point of view, this new methodology is advantageous since it allows one to capture the important dynamics of the system with a reduced model of much lower order than the ones provided by modal and balanced based truncations. This advantage is capital in the case where the system has many (infinite number) of important dynamics like the power system treated here.

From the physical point of view, the advantage of this approach is to be structural, i.e., to consider the entire system when the reduction question is addressed and not to substitute a priori simplified models for parts of the system like its transmission lines. The new resulting reduced models are more adequate to reproduce the phenomena of interest in the power system and need less trial and error iterations.

This systemic and modal view will be further exploited to unify and classify some models actually used in different simulators. In this manner, the Transmission Systems Operators (TSOs) could improve their simulation methodology and tools by using a reduced number of simulation models/softwares and better connecting these models to physical phenomena like inter-area and hypo-synchronous oscillations. This will also allow the TSOs to take into account the new electric components (like, e.g., energy sources, High Voltage Direct Current links (HVDC), electronic switches,..).

Fig. 1 .

 1 Fig. 1. Power system

Fig. 3 .

 3 Fig. 3. Step responses with of π and DP models for the line

  -2 H/Km, C = 10 -8 F/Km, = 50Km, Z L = 500Ω, Z 0 = 0, n = 1000. Table I presents the poles and the infinite norm of each partial fraction in the transfer function of the π and DP model. The infinite norm of the partial fractions of the DP model decreases very slowly and their is no important gap between two consecutive norms of partial fractions (i.e., G r ∞ G r+1 ∞ ) which might give an index for the modal truncation. For this reason, taking, e.g., = 0.1 as H ∞ error bound of (10), i.e, n k=N |r k | |Re(β k )| = , the resulting reduced model of the modal truncation is of order 981. As a large number of dynamics are kept in this reduced model, its trajectories are almost the same with the ones provided by the DP model shown in Fig.6. This means that this reduced model is adequate to reproduce all phenomena of interest in the power system.

1 )

 1 Hankel singular values: The Hankel singular values are σ i = λ i (W o W c ) where Γ is the Hankel operator

Fig. 4 .

 4 Fig. 4. Hankel singular values of the DP model.

Fig. 5 .

 5 Fig. 5. Bode diagrams of balanced truncation reduced model and DP model.

Fig. 6 .

 6 Fig. 6. Step responses comparisons between balanced truncation reduced model, π, DP and EMTP models.

  a numerical approximation of operators Ã, B, C like in Section III-A.2. The participation factors are the elements of the matrix P = p 1 , p 2 , . . . , p n with p i = φ 1i ψ i1 φ 2i ψ i2 • • • φ ni ψ in T where φ ki is the kth entry of the ith right eigenvector Φ i of A n and ψ ik is the kth entry of the ith left eigenvector Ψ i of A n . The following properties hold Ψ i Φ i = 1 and n i=1 p ki = n k=1 p ik = 1.

Fig. 7 .

 7 Fig. 7. Comparison between step responses of the new reduced model, π and DP models.TABLE II FIRST POLES OF THE NEW REDUCED MODEL states poles frequency (KHz) x1 λ1;2 = 5, 35.10 2 ± 8, 88.10 3 i 1, 41 λ3;4 = 5.35.10 2 ± 4.44.10 3 i 0, 70 x2 λ5;6 = 5, 35.10 2 ± 1, 33.10 3 i 2, 11 λ1;2 1, 41 x3 λ7;8 = 5, 35.10 2 ± 1, 77.10 3 i 2, 81 λ5;6 2, 11

Fig. 8 .

 8 Fig. 8. Bode diagram-comparison between the new reduced, the π and DP models.

TABLE I POLES

 I AND H∞ NORM OF EACH PARTIAL FRACTION • ∞= |r k | |Re(λ k )|

		full model		π-model
	k	poles (λ k,k+1 )	• ∞	poles	• ∞
	1; 2	-3, 94.10 2 ± 6, 85.10 5 i	1, 0000789	-5, 62.10 2	1, 03
	3; 4	-4, 14.10 2 ± 6, 34.10 5 i	1, 0000788	-74, 621.10 2	0, 07
	5; 6	-4, 12.10 2 ± 6, 38.10 5 i	1, 0000786
	7; 8	-3, 96.10 2 ± 6, 81.10 5 i	1, 0000782
	9; 10	-3, 97.10 2 ± 6, 77.10 5 i	1, 0000777
	11; 12	-4, 09.10 2 ± 6, 46.10 5 i	1, 0000771
	13; 14	-4, 11.10 2 ± 6, 42.10 5 i	1, 0000764
	15; 16	-4, 06.10 2 ± 6, 53.10 5 i	1, 0000757
	17; 18	-4, 00.10 2 ± 6, 69.10 5 i	1, 0000745
	19; 20	-3, 99.10 2 ± 6, 73.10 5 i	1, 0000737
	. . .
	591	-5, 48.10 2	0, 975424
	. . .

ACKNOWLEDGMENT

The authors would like to thank Dr. A. Parisot from RTE and Dr. H. Mounier from LSS for their support.