TRANSIENT BASED METHOD FOR IDENTIFYING FAULTS IN A HIGH / MEDIUM VOLTAGE ELECTRIC POWER TRANSMISSION SYSTEM

Abstract

The invention relates to a transient based method for identifying faults in an electric power transmission and/or distribution system (100) having at least one current transmission line (L.sub.12, L.sub.13, L.sub.23) comprising the following steps: —generation of a physical model of the at least one current transmission line (L.sub.12, L.sub.13, L.sub.23), the physical model depending on the fault parameters and describing the behavior of voltage and/or current transients due to the fault in the at least one current transmission line, fault parameters comprising a fault location parameter on said current transmission line (L.sub.12, L.sub.13, L.sub.23) and a fault impedance parameter, —measurement of voltage and/or currents evolution at least at one specific location in the said power system (100), —iterative simulation of the voltage and/or current evolution by the physical model at the measurement point with a set of fault parameters where at each step of iteration, simulated and measured voltage and/or current evolutions are compared and the set of fault parameters is adapted according to a convergence criterion, —identification of a fault with its fault parameters when convergence of the measured voltage and/or current evolutions and simulated voltage and/or current evolutions is reached in a limited number of iterations.

Claims

1. A transient based method for identifying faults in an electric power transmission and/or distribution system comprising a current transmission line, said method comprising generating a physical model of the current transmission line, the physical model depending on fault parameters and describing the behavior of voltage and/or current transients due to a fault in the current transmission line, said fault parameters comprising a fault location parameter on said current transmission line and a fault impedance parameter, measuring voltage and/or current evolution at a specific location in said power system, iteratively simulating the voltage and/or current evolution by the physical model at a measurement point with a set of fault parameters, wherein each step of iteration comprises comparing simulated and measured voltage and/or current evolutions and adapting the set of fault parameters according to a convergence criterion, and identifying a fault and fault parameters thereof upon having reached convergence of the measured voltage and/or current evolution and the simulated voltage and/or current evolution.

2. The method according to claim 1, wherein the measurements of voltage and/or current evolution are carried out at several specific locations in the said power system, each corresponding to a measurement point and wherein the steps of iteratively simulating of the voltage and/or current evolution and identifying a fault with its fault parameters are carried out for each measurement point independently from the other measurement points.

3. The method according to claim 1, wherein the physical model is based on a telegraph equation for describing the fault parameters and the behavior of travelling waves due to the fault in said electrical power transmission and/or distribution system.

4. The method according to claim 3, wherein only traveling waves are taken into account when iteratively simulating the voltage and/or current evolution that can reach the measurement point before a predetermined maximum time t.sub.max between one and ten milliseconds.

5. The method according claim 1, further comprising extending the physical model to take into account soil resistivity.

6. The method according to claim 5, wherein extending the physical model for taking into account soil resistivity comprises implementing a filtering function.

7. The method according to claim 1, wherein said electrical power transmission and/or distributor system is a high or medium voltage DC grid comprising overhead current transmission lines and/or underground and/or submarine current transmission lines.

8. The method according to claim 1, wherein said electrical power transmission and/or distributor system is a high or medium voltage AC grid.

9. The method according to claim 1, wherein said physical model takes into account modal voltage and current.

10. The method according to claim 1, wherein the convergence criterion comprises the minimization of a cost function of the difference of simulated and measured voltage and/or current evolutions.

11. The method according to claim 1, wherein identifying identification of a fault comprises verifying that the obtained fault parameters belong to a predefined validity domain.

12. The method according to claim 1, wherein each iteration step comprises completing voltage and/or current evolutions at specific locations in said electrical power transmission and/or distribution system with ongoing measurements at said measurement point.

13. The method according to claim 1, wherein the specific location is at one of an end of said current transmission line (L.sub.12, L.sub.13, L.sub.23), a relay, a breaker module, and a station.

14. The method according to claim 1, further comprising detecting unusual behavior with respect to a nominal functioning of the electrical power transmission and/or distribution system and thereafter iteratively simulating the voltage and/or current evolution by the physical model.

15. The method according to claim 14, wherein detecting said unusual behavior comprises detecting at least one of a voltage drop, a current increase, and a crossing of a threshold.

16. The method according to claim 1, further comprising, after having identified said fault, generating an output in order to activate protection circuits in said electrical power transmission and/or distribution system.

17. An apparatus comprising a fault-identification module for identifying faults in an electrical power transmission and/or distribution system having a current transmission line (L.sub.12, L.sub.13, L.sub.23) the fault-identification module being configured for: memorization of a physical model of the electrical power transmission and/or distribution system with said current transmission line, the physical model describing fault parameters and behavior of travelling waves due to a fault in the electrical power transmission and/or distribution system and, in case of a fault, fault parameters comprising a fault-location parameter on said current transmission line (L.sub.12, L.sub.13, L.sub.23) and a fault-impedance parameter, measurement of voltage and/or current evolution at a specific location in said power system, iterative simulation of the voltage and/or current evolution by the physical model at the measurement point with a set of fault parameters, wherein, at each iteration of said iterative simulation, simulated and measured voltage and/or current evolutions are compared and the set of fault parameters is adapted according to a convergence criterion, identification of a fault and fault parameters thereof when convergence of the measured voltage and/or current evolutions and simulated voltage and/or current evolutions has been reached.

18. The apparatus according to claim 17, wherein the physical model is based on a telegraph equation for describing the fault parameters and the behavior of travelling waves due to the fault in said electrical power transmission and/or distribution system.

19. The apparatus according to claim 17, wherein the physical model is extended to take into account soil resistivity.

20. The apparatus of claim 17, further comprising a high or medium voltage electric power transmission and/or distribution system that comprises said fault identification module for controlling at least one protection circuit in a high/medium voltage electric power transmission and/or distribution system in case of a fault identification.

Description

BRIEF DESCRIPTIONS OF THE DRAWINGS

[0050] Other advantages and characteristics will appear with the reading of the description of the following figures:

[0051] FIG. 1 shows in a synoptic scheme of an example of a three nodes MTDC network,

[0052] FIG. 2 shows a Bewley lattice diagram showing the position and direction of motion of every generated, reflected, and transmitted wave caused by a fault on a transmission line,

[0053] FIG. 3 shows four graphs, first and second graphs from top to bottom show the frequency dependency of the line characteristic impedance Z.sub.c and third and fourth graphs from top to bottom show the frequency dependency of the line propagation constant γ,

[0054] FIG. 4 shows two graphs modelling voltage and current in function of time for soil resistivity effect modeling using ARX model,

[0055] FIG. 5 shows a graph filter coefficients in function of fault distance,

[0056] FIG. 6 shows two graphs modelling voltage and current in function of time for soil resistivity effect modeling using ARX11 model,

[0057] FIG. 7A shows a synoptic presentation of a first example of an extended physical model taking into account soil resistivity effect,

[0058] FIG. 7B shows a synoptic presentation of a second example of an extended physical model taking into account soil resistivity effect,

[0059] FIG. 8 shows an example of a DC protection relay

[0060] FIG. 9A shows an example of a simplified flowchart of the method and computer program according to the invention,

[0061] FIG. 9B shows a synthetic diagram of inputs and outputs of operations according to the simplified flowchart of FIG. 9A,

[0062] FIG. 10 shows a scheme of a meshed grid with single-conductor over-headlines and 3 MMC stations,

[0063] FIG. 11 shows the result of a first example of simulation in a graph the fault resistance in function of the fault distance with cost function iso-levels and the estimation trajectory,

[0064] FIG. 12 relates to FIG. 11 and shows in two graphs current and voltage transient behavior in function of time for comparing simulated measurements with the method according to the invention,

[0065] FIG. 13 shows the result of a second example of simulation in a graph the fault resistance in function of the fault distance with cost function iso-levels and the estimation trajectory,

[0066] FIG. 14 is related to FIG. 13 and shows in two graphs current and voltage transient behavior in function of time for comparing simulated measurements with the method according to the invention,

DETAILED DESCRIPTION

[0067] The embodiment(s) in the following description are only to be considered as examples. Although the description may refer to one or several embodiments, this does not mean inevitably that every reference concerns the same embodiment, or that the characteristics apply only to a single embodiment. Simple characteristics of various embodiments can be also combined to new embodiments that are not explicitly described.

[0068] In general, the invention applies to high and medium AC and DC voltage electric power transmission and/or distribution systems or grids.

[0069] By way of example, a specific embodiment related to high or medium voltage direct current network will be described hereafter.

[0070] FIG. 1 is a synoptic scheme of an example of a three nodes MTDC (Multi-Terminal high voltage Direct Current) network 100. The grid 100 comprises three DC current transmission lines L.sub.12, L.sub.13, L.sub.23. At each end of the transmission lines L.sub.12, L.sub.13, L.sub.23 are disposed respectively one DC breaking module R.sub.i,j (i, j=1, 2, 3, i≠j).

[0071] The ends of two neighboring DC breaking modules R.sub.i,j (i, j=1, 2, 3, i≠j), for example R.sub.31 and R.sub.32, are connected together to further DC breaking modules R.sub.i,j (i=1, 2, 3) which are further connected to DC/AC converters 113, 213, 313.

[0072] The breaking modules R.sub.i,j (i, j=1, 2, 3) comprise protection circuits that can be activated to protect the high voltage DC grid and in particular the converters 113, 213 and 313.

[0073] DC breaker module R.sub.11 and converter 113 can be disposed in one station. Similar DC breaker module R.sub.22 and the converter 213 can also be disposed in one station and DC breaker module R.sub.33 and the converter 313 can be disposed in one station either.

[0074] Physical Modeling

[0075] Faults occurring on DC transmission lines L.sub.12, L.sub.13, L.sub.23, can be described in the first milliseconds after the fault using the framework of traveling waves. Along a DC transmission line L.sub.12, L.sub.13, L.sub.23, current and voltage verifies the Telegraphist's equations that can be written in the Laplace domain as:

[00001]$\begin{array}{cc}\frac{{\partial }^{2}V}{\partial x^2}=Z\left(s\right)Y\left(s\right)V& \left(1\right)\\ \frac{{\partial }^{2}I}{\partial x^2}=Y\left(s\right)Z\left(s\right)I& \left(2\right)\end{array}$

[0076] where Z(s)=R+sL and Y (s)=G+sC are the transfer functions of the distributed series impedance and shunt admittance, respectively. When a fault occurs, voltage and current waves starting at fault location travel along the line according to

[00002]$\begin{array}{cc}V\left(s;d\right)=H\left(s;d\right)V_{\mathrm{init}}\mathrm{and}I\left(s;d\right)=Z_c^{-1}V\left(s;d\right)& \left(3\right)\end{array}$

[0077] where V.sub.init is the initial surge at fault location and d the distance to the fault.

[0078] Hence the characteristic impedance Zc and the propagation function H characterizes the behavior of the line such that,

[00003]$\begin{array}{cc}H\left(s,d\right)=\exp \left(-\sqrt{YZ}d\right)\mathrm{and}{Z}_c\left(s\right)=\sqrt{Z\text{/}Y}& \left(4\right)\end{array}$

[0079] When a change of propagation medium occurs, the forward wave V.sub.f gives rise to transmitted and reflected waves, V.sub.t and V.sub.r respectively. The associated voltage at the station then is V=V.sub.t=V.sub.f+V.sub.r=(1+K)V.sub.f=TV.sub.f, where the transmission and reflection coefficients T and K depend on the media characteristic impedance in particular of a related current transmission line L.sub.12, L.sub.13 or L.sub.23, i.e., for a wave from media of characteristic impedance Z1 to media of characteristic impedance Z2,

[00004]$\begin{array}{cc}{K}_{1.fwdarw.2}=\frac{Z_2-Z_1}{Z_2+Z_1}\mathrm{and}{T}_{1.fwdarw.2}=\frac{2Z_2}{Z_2+Z_1}.& \left(5\right)\end{array}$

[0080] Hence using those coefficients and the modeling of the line as described before one can track the traveling waves using a Bewley lattice diagram as shown in FIG. 2.

[0081] A Bewley lattice diagram shows the position and direction of motion of every generated, reflected, and transmitted wave caused by a fault on a transmission line.

[0082] As can be seen in FIG. 2, a fault has occurred, for example in current transmission line L.sub.13, at t.sub.fault at a distance d.sub.f,x from a station x (for example DC breaking module R.sub.31, or a station comprising DC breaker module R.sub.11 and the converter 113) and d.sub.f,y from a station y (for example DC breaking module R.sub.31, or a station comprising DC breaker module R.sub.33 and the converter 313). In this diagram station x and y can for example be respective DC breaking modules R.sub.13, R.sub.31; R.sub.12, R.sub.21; R.sub.13, R.sub.31 at the end of the current transmission lines L.sub.12, L.sub.13 or L.sub.3 or the stations described above. On this FIG. 2, one can see the propagation of two wave fronts W.sub.x and W.sub.y.

[0083] However, to derive explicit expressions for the different waves, one has to precise the different terms in the model. Regarding the transmission line, a single conductor overhead line is considered. The effect of soil resistivity in the ground return path is neglected at this point, so that the characteristic impedance Z and admittance Y depends only on the overhead line characteristics. Moreover a lossless approximation, where series resistance R and shunt conductance G are neglected, allows for explicit expressions of Zc and H as functions of the per unit length series inductance L and shunt capacitance C

[00005]$\begin{array}{cc}H\left(s\right)=\exp \left(-s\sqrt{LC}d\right)\mathrm{and}{Z}_c=\sqrt{L\text{/}C}.& \left(6\right)\end{array}$

[0084] The relevance of this approximation is studied by plotting the evolution of the characteristic impedance Z.sub.c(s)=√{square root over (Z/Y)} and the propagation constant γ(s)=√{square root over (YZ)} as a function of the frequency. This shown in FIG. 3. where first and second graphs from top to bottom show the frequency dependency of the characteristic impedance Z.sup.c and third and fourth graphs from top to bottom show the frequency dependency of the propagation constant γ.

[0085] Since the transient behavior has frequency content above 1 kHz, the lossless approximation is close to the general lossy case. Moreover, the frequency dependency of the line parameters themselves (L; C) was investigated using frequency dependent line model calculation function from EMTP-RV, showing that the line parameters do not depend on frequency.

[0086] Regarding the power converters 113, 213 or 313 at each station, Modular Multilevel Converters (MMC) are considered. With the rise of current, power electronics within the MMC will protect by auto blocking if the current reaches 2 p.u. which corresponds to currents ranging from 3 to 9 kA, depending on the characteristics of the grid.

[0087] Since only modeling of the 2-3 first waves after the fault is necessary to identify and localize the fault, it is reasonable to consider those components are not blocked at a first time.

[0088] Under that assumption, one can use a RLC equivalent for the MMC, i.e., Z.sub.mmc(s)=R.sub.mmc+sL.sub.mmc+1/(sC1.sub.mmc). Regarding the initial surge at fault location V.sub.init, it can be modeled as a step voltage in series with the fault impedance. The step voltage is connected at the fault instant t.sub.f and has an amplitude opposite to the one just before fault occurrence at fault location,

[00006]$\begin{array}{cc}{V}_{\mathrm{init}}\left(s,t_f,Z_f\right)=-{\upsilon }_{\mathrm{prior}\mathrm{fault}}\frac{Z_c}{Z_c+2Z_f}\exp \left(-s{t}_f\right).& \left(7\right)\end{array}$

[0089] Note the prior fault voltage at the fault location ν.sub.prior fault, is not measured. However, the prior fault behavior can be described using the Telegraphist's equations in steady state,

[00007]$\frac{{\partial }^{2}V}{{\partial }^2{x}^2}=\mathrm{RG}V$

where R and G are the distributed series resistance and shunt conductance, respectively. Typical values of those parameters are small, hence the voltage do not vary significantly along the distance for lines of several hundred km long. Thus, v.sub.prior fault can be replaced by the steady state voltage measured at the relay just before the fault is detected. Finally, the fault impedance accounts for different effects such as the electric arc, the tower grounding impedance, and resistance of additional unknown objects in the current path. It is generally considered those different elements can be accounted for using a single unknown fault resistance R.sub.f. Hence knowing the characteristics of the different lines, stations and the topology of the gridone can derive expressions in the frequency domain for any traveling waves.

[0090] As an example, consider a fault occurred on the line e.sub.x,y between stations x and y, at a distance d.sub.f,x from the station x (FIG. 2). The nth forward and return waves at the station due to reflection at fault location is then expressed as a function of the fault parameters p in the Laplace domain as

[00008]$\begin{array}{cc}\{\begin{array}{c}{V}_{f,n}\left(s,p\right)={\left(K_x{K}_f\right)}^{n-1}{H\left(s,d_{f,x}\right\}}^{2n-1}{V}_{\mathrm{init}}\left(s,t_f\right)\\ V_{r,n}\left(s,p\right)=K_x{V}_{f,n}\left(s,p\right)\end{array}& \left(8\right)\end{array}$

where K.sub.x and K.sub.f are the reflection coefficients from the line to the station x and to the fault respectively, as defined in equation (5).

[0091] Using the previous methodology, one can then get back to temporal domain expressions using inverse Laplace transform. This is due to the different hypotheses made before, such as the loss-less approximation and the linear modeling of the MMC stations. Examples of such formulas for the first wave in equation (8) are given below

[00009]${\upsilon }_{f,1}\left(t\right)=-\frac{V_f{R}_f}{Z_c+2R_f}u\left(t-t_d\right)$$\mathrm{and}$${\upsilon }_{r,1}\left(t\right)=\frac{2\sqrt{C_{\mathrm{mmc}}}{V}_f{Z_c}^2}{\left(Z_c+2R_f\right)A}\sinh \left(\frac{\left(t-t_d\right)A}{4\sqrt{C_{\mathrm{mmc}}}{L}_{\mathrm{mmc}}}\right)\times \exp \left(-\frac{\left(2R_{\mathrm{mmc}}+Z_c\right)\left(t-t_d\right)}{4L_{\mathrm{mmc}}}\right)u\left(t-t_d\right)$

[0092] where A=√{square root over (C.sub.mmc(2R.sub.mmc+Z.sub.c).sup.2−16L.sub.mmc)} is a scalar coefficient depending solely on the line surge impedance and the MMC parameters. The time t.sub.d is the detection time corresponding to the arrival time of the first wave. It can be expressed as a function of the fault time t.sub.f, the fault distance d.sub.f and the wave speed c.sub.w=1/√{square root over (LC)}; t.sub.d=t.sub.f+d.sub.f/c.sub.w. Similar expressions can be obtained for the other waves due to reflection and the fault and for waves due to reflection at other stations. The total voltage at the station is then the sum of the measured prior fault voltage ν(t=t.sub.d and the different waves arriving at the station as derived in the Bewley diagram in FIG. 2.

[00010]$\begin{array}{cc}{\upsilon }_m\left(t\right)=\upsilon \left(t_d^{-}\right)+{\upsilon }_{f,1}\left(t\right)+{\upsilon }_{r,1}\left(t\right)+{\upsilon }_{f,2}\left(t\right)+{\upsilon }_{r,2}\left(t\right)+& \left(9\right)\end{array}$

[0093] The number of reflections taken into account depends on a trade off between the complexity and desired time validity range of the model.

[0094] In most of the cases, in order to detect and locate faults rapidly for allowing to shut down for example a current transmission line where a fault occurred, the first wave coming from the fault along with the first reflected wave due to the reflection at the fault or at the remote station are needed. In many cases, the very first wave coming directly after fault occurrence is sufficient. Depending on the requested precision needed, more waves (2 or 3) even up to 4 waves can be taken into account to identify and locate the fault.

[0095] Since the characteristic impedance is a scalar, the current derivation can be deduced from the voltage.

[00011]$\begin{array}{cc}{i}_m\left(t\right)=i\left(t_d^{-}\right)+Z_c^{-1}({\upsilon }_{f,1}\left(t\right)-{\upsilon }_{r,1}\left(t\right)+{\upsilon }_{f,2}\left(t\right)\mathrm{.Math.})& \left(10\right)\end{array}$

[0096] The obtained formulas depend on the known characteristics of the grid, and on measured quantities such as the detection time and the prior fault conditions and also on the fault parameters p=(e.sub.x,y; d.sub.f,x; d.sub.f,y; R.sub.f; t.sub.f).sup.T that are unknown and need to be identified.

[0097] In conclusion, an analytical parametric model (physical model) based on wave propagation and allowing to predict the current and voltage evolution in case of fault occurrence, depending on a set of physical fault parameters has been developed here above. This physical model can be inverted in a way that when a fault is suspected, one can determine the fault parameters from current and voltage measurements available at dedicated points in the grids for example at a relay or a DC breaker circuit.

[0098] It is important to notice that only one physical model is necessary and that by the way of inversion, fault parameters can be determined in using only the current and voltage measurements realized some milliseconds upon appearance of the faults.

[0099] Soil Resistivity Modeling

[0100] As mentioned, the physical model above neglects the effect of the soil resistivity, noted ρ. In many cases, in particular if the value of the soil resistivity is very low, for instance in case of return through see water, soil resistivity may be neglected. Besides, the effect of soil resistivity increases with the fault distance. Hence if the fault is suspected to be quite close to the relay (up to 5 km for instance) the soil resistivity can also be neglected. However, in a more developed approach soil resistivity may be taken into account.

[0101] While different ways to model physically the behavior of ground return exist, they do not lead to an analytic and explicit solution of the current and voltage transient evolution.

[0102] In order to complete soil resistivity effects into the previously obtained model a behavioral approach is proposed taking the previously developed physical model neglecting earth influence (ν.sub.m (p; t); i.sub.m (p; t)) as input and produces an new model (ν.sub.m,g (p; t; ρ); i.sub.m,g (p; t; ρ)) that still depends on the fault parameters.

[0103] To do this, an extended model (ARX model—Autoregressive model with eXternal inputs) is proposed that is a filtering of the previous physical model:

[00012]$\begin{array}{cc}{y}_{m,g}\left(p,\rho ,t_k\right)=G\left(z^{-1}\right)y_m\left(p,t_k\right)& \left(11\right)\end{array}$

[0104] where G(z.sup.−1)=(B(z.sup.−1)/A(z.sup.−1)) is a function to be determined. The polynomial A(z.sup.−1)=1−a.sub.1z.sup.−1− . . . −a.sub.naZ.sup.na is an autoregressive part that represents the inductive effect of the soil resistivity while B(z.sup.−1)=z.sup.nk(b.sub.0+b.sub.1z.sup.−1+ . . . +b.sub.naz.sup.nb) is a exogenous part that takes the previous model as external input.

[0105] Coefficients n.sub.a and n.sub.b are the orders of the polynomials A and B respectively, n.sub.k is a dead-time that represents a delay in the system response.

[0106] Note the models γ.sub.m=(ν.sub.m; i.sub.m).sup.T and γ.sub.m,g=(ν.sub.m,g; i.sub.m,g).sup.T include both current and voltage. In order to identify coefficients a.sub.i; b.sub.i the measurements model is to be considered:

[00013]$\begin{array}{cc}A\left(z^{-1}\right)y\left(t_k\right)=B\left(z^{-1}\right)y_m\left(p,t_k\right)+\mathrm{.Math.}\left(k\right)& \left(12\right)\end{array}$

[0107] where γ=(ν; i).sup.T are the measurements and ε(k) is a measurement error term. Hence the parameters a.sub.i; b.sub.i can be identified using least squares by minimizing ε.sup.Tε. One then must select appropriate orders by comparing the goodness of fit for different values of (n.sub.d; n.sub.b; n.sub.k).

[0108] FIG. 4 shows two graphs modelling voltage and current in function of time for soil resistivity effect modeling using an ARX 11 model (meaning the orders of the polynomials are n.sub.a=1, n.sub.b=0, i.e., one coefficient a.sub.1 and one coefficient b.sub.0 are to be identified) as explained above.

[0109] The measurements data are obtained using Electro-Magnetic Transient (EMT) software EMTP-RV. One can see that the filtering stage applied to the approximate physical model improves the accuracy of the model) with respect to the simulated measurement data (in dots EMTP simulation).

[0110] As suggested in equation (11), the transfer function G that models the effect of soil resistivity depends on the value of ρ. Moreover, the inductive effect of the ground is expected to increase with the distance to the fault, so that G also depends on the fault parameters p through d.sub.f. Hence, for a given set of fault parameters p and a given value of soil resistivity ρ, one can find coefficients a.sub.i(p; ρ); b.sub.i(p; ρ) that defines a filter G(p; ρ). The evolution of the filter parameters as function of the fault distance can be made explicit using some interpolation function such as

[00014]$a_1\left(d_f\right)={\alpha }_{1,0}+{\alpha }_{1,1}{d}_f^{-1}+{\alpha }_{1,2}{d}_f^{-2}{b}_0\left(d_f\right)={\beta }_{0,0}+{\beta }_{0,1}{d}_f^{-1}+{\beta }_{0,2}{d}_f^{-2}.$

[0111] An example of the evolution of such coefficients with the fault distance and the corresponding fitting function for an ARX 11 model is plotted on FIG. 5 showing a graph of filter coefficients in function of fault distance d.sub.f. Note the inductive effect, represented by parameter a.sub.1 is more important as distance increases whereas the approximate physical model, through b.sub.0, is predominant for small d.sub.f.

[0112] Hence the extended model is able to predict the effect of soil resistivity on the first wave using a filter G(d.sub.f; ρ) that depends explicitly on the fault parameters through d.sub.f. To adapt this approach to the following waves due to reflection at the fault or at the remote station, it is proposed to aggregate the effect of soil resistivity using a single filter on each wave. Those filters are the one obtained previously through the modeling of the first wave but the fault distance is adjusted to account for the total distance traveled by a given wave. For instance the first wave due to reflection at the fault location travels thrice the fault distance and is thus obtained by filtering the model with soil effect γ.sub.m,2 through the filter given by G(3d.sub.f; ρ). Similarly the first wave coming from the remote station uses a filter G(2D-d.sub.f; ρ). Such an approach makes the hypothesis that the effect of soil resistivity can be lumped in a single filter for each wave but it allows to capture most of the effect of the soil using an single adaptable filter. FIG. 6 shows two graphs modelling voltage and current in function of time for soil resistivity effect modeling using ARX11 model with a fault at 30 km from a relay station. As can be seen on the FIG. 6, the extended physical model correctly fits with the simulated measurement data (in dots EMTP simulation)

[0113] FIG. 7A shows a synoptic presentation of a first example of the extended physical model taking into account soil resistivity effect as presented here above.

[0114] The obtained model is a combination of physical and behavioral approaches and depends explicitly on the fault parameters.

[0115] Thus, in measuring voltage and current at dedicated points in the grid, by inverting the extended physical model when a fault is suspected, one can determine the fault parameters from those current and voltage measurements available at dedicated points in the grids for example at a relay or a DC breaker circuit, and this even in the case soil resistivity cannot be neglected.

[0116] FIG. 7B shows a synoptic presentation of a second example of the extended physical model taking into account soil resistivity effect as presented here above. With respect to FIG. 7A, the physical model may be aggregated with the behavioral model with a part that does not pass into the behavioral model. As mentioned above, it may be sufficient to use the physical model if soil resistivity effects are negligible (if the fault distance is small or the value of the soil resistivity is small). Otherwise, it is preferable to use the extended model with the behavioral approach.

[0117] Systematic Modeling of Traveling Waves

[0118] Some electric power transmission and/or distribution systems may include several junction points where waves are reflected and transmitted. This is for example the case when current transmission lines are hybrid, i,e., composed of overhead current transmission lines and/or underground and/or submarine current transmission lines.

[0119] In such a network a systematic approach is needed to account for any possible wave traveling through a grid and to filter out non-pertinent signals. When a fault occurs on a hybrid transmission line, it is important to identify which section is faulty because faults on overhead line are generally temporary. This allows reclosing the line after some time typically 200-500 ms whereas faults on underground or submarine cables are generally permanent and reclosing should be avoided.

[0120] A protected network can be described by a graph g containing a set of vertexes ν.sub.f representing nodes between two or more current transmission line segments, like stations or junctions between overhead current transmission lines and underground current transmission lines. A vertex is denoted q.sub.i. The graph contains a set of edges ε.sub.f representing the network current transmission lines segments, which can be overhead current transmission lines or underground current transmission lines. The edge between the nodes q.sub.i, q.sub.j is noted e.sub.i,j. A fault location corresponds to a node q.sub.f and the considered measurement point is the node q.sub.s.

[0121] All possible traveling waves in a grid due to a fault can be described by every possible paths between the fault location and the measurement point, the node q.sub.s. All the possible paths can therefore be described by the set

[00015]${\mathrm{��}}_{q_f.fwdarw.q_s}=\left\{\pi =\left(q^1,\mathrm{.Math.},q^n\right)|q^1=q_f,q^n=q_s,q^i\in v_f,\left(q^i,q^{i+1}\right)\in {\mathrm{.Math.}}_f\right\}$

[0122] In order to compute only finite number of traveling waves, one only considers the waves that can reach the node q.sub.s before some maximum time t.sub.max. For instance one possible value is t.sub.max=2d.sub.L/c, where d.sub.L is the total length of the protected line and c is the speed of light. t.sub.max is typically in a range between 1-10 ms, more specifically between 1-5 ms.

[0123] Consider a particular voltage wave traveling through a path π. As the loss-less Telegraph equations correspond to a pure delay, the output of the physical model can be written in the Laplace domain as

[00016]$V_{\pi }^{\phi }\left(s,p\right)=V_{\pi ,i}\left(s,R_f\right)\exp \left(-{\tau }_{\pi }\left(d_f\right)s\right)$

[0124] where τ.sub.π accounts for the total delay due to the propagation of the wave through the path π, which depends on the fault distance. ν.sub.π,i accounts for all the reflections and transmissions that occur at each junction. For a given path π, one can gather the coefficients involved at each junction

[00017]$\underset{i=1}{\overset{n-1}{.Math.}}{J}_{\left(q^{i-1},q^i\right).fwdarw.\left(q^i,q^{i+1}\right)}\left(s,r_f\right)$ [0125] Where the coefficients can be reflection or transmission coefficients

[00018]$J_{\left(q^{i-1},q^i\right).fwdarw.\left(q^i,q^{i+1}\right)}=\{\begin{array}{c}{T}_{\left(q^{i-1},q^i\right).fwdarw.\left(q^i,q^{i+1}\right)},q^{i-1}\ne q^{i+1}\\ K_{\left(q^{i-1},q^i\right).fwdarw.\left(q^i,q^{i+1}\right)},q^{i-1}=q^{i+1}\end{array}$ [0126] By convention the first coefficient corresponds to the initial surge at the fault location. The interface part ν.sub.π,i of the physical model gathers all those coefficients, in addition to the initial step occurring at the fault instant

[00019]$V_{\pi ,i}=\underset{i=1}{\overset{n-1}{.Math.}}{J}_{\left(q^{i-1},q^i\right).fwdarw.\left(q^i,q^{i+1}\right)}\left(s,r_f\right)\times \exp \left(-{\mathrm{st}}_f\right)$

[0127] The physical model computes all possible traveling waves on the network that can reach the observation point before some maximum time t.sub.max with respect to the fault inception time at node q.sub.s. Those different traveling waves can be displayed in the form of a Bewley Lattice diagram.

[0128] Behavioral Modeling

[0129] A wave traveling through an edge e either corresponding to a portion of over-head current transmission line or an underground current transmission line is both distorted and attenuated. Those effects are in particular due to the frequency dependent behavior of the current transmission lines and the influence of the soil resistivity. To account for such effects, a behavioral approach is applied using the step response of the transmission lines. The step response of a given edge e can be obtained using EMT simulation software such as EMTP, knowing the geometry of the transmission lines and the characteristics of the conductors. In particular, the step response depends on the length of the considered segment d.sub.e and the value of the soil resistivity ρ.sub.e Considering a set of known step responses b.sub.d,ρ(k) for various edge lengths (d.sub.1; d.sub.2, . . . , d.sub.n) and various soil resistivity (ρ.sub.1; ρ.sub.2; . . . ; ρ.sub.n). In order to get a step response for any length and any soil resistivity, the following interpolation for the step response b.sub.d,ρ(k) is proposed, where ρ.sub.1<ρ<ρ.sub.2 and d.sub.1<d<d.sub.2.

[00020]$b_{d,\rho }\left(k\right)=\frac{b_{d,{\rho }_2}\left(k\right)-b_{d,{\rho }_1}\left(k\right)}{\log \left({\rho }_2\right)-\log \left({\rho }_1\right)}\left(\log \left(\rho \right)-\log \left({\rho }_1\right)\right)+b_{d,{\rho }_1}\left(k\right)$$b_{d,{\rho }_1}\left(k\right)=\frac{b_{d_2,{\rho }_1}\left(k\right)-b_{d_1,{\rho }_1}\left(k\right)}{\left(d_2-d_1\right)}\left(d-d_1\right)+b_{d_1,{\rho }_1}\left(k\right)$

[0130] The step response of a given edge e is differentiated to obtain the impulse response h.sub.e(k). The total distortion and attenuation of a wave is obtained by combining the impulse responses of all the traveled edges.

[0131] The principle of cascading the different impulse response filters is shown in Figure XX. Hence, the total output of the model of a given wave is then obtained as.

[00021]$v_{\pi }^m=v_{\pi }^{\phi }*\underset{e\in \pi }{.Math.}{h}_e$

[0132] Where the inversion of Laplace domain expressions can be done using the inverse Fourier transform. The total voltage model is composed of all the waves travelling through a path π from the fault location to the observation point in less than a time t.sub.max.

[00022]$v^m\left(p,t_k\right)=\underset{\pi \in {\mathrm{��}}_{q_f.fwdarw.q_s,t_{m\mathrm{ax}}}}{.Math.}{v}_{\pi }^m\left(p,t_k\right)$

[0133] Fault Identification and Localization

[0134] Taking into account the above exposures, in only monitoring current i and voltage v at dedicated locations of the grid, in particular at the terminals of the currents transmission lines like L.sub.12, L.sub.13 or L.sub.23 for example, one can invert the calculations in order to determine the relevant fault parameters, in particular localization of the fault. It is important to notice that each dedicated location can deliver the information of the relevant fault parameters, in particular identification and localization of a fault without taking into account the measurements realized at another location.

[0135] FIG. 8 shows an example of relay R.sub.i,j (i, j=1, 2, 3) equipped with a fault identification module 20. The part of DC breaking module R.sub.i,j (i, j=1, 2, 3) dedicated to break the circuit is known to the man skilled in the art and not described in more detail.

[0136] Such fault identification module 20 can also be installed in other units of the grid and in other locations without departing from the invention. In particular one may envisage that such fault identification module 20 is installed in a converter station or more generally in a relay.

[0137] For improving protection of the grid 100, it seems appropriate that one fault identification module 20 is installed at each end of a current transmission line L.sub.12, L.sub.13 or L.sub.23 allowing to identify rapidly a fault and control and command for example the two DC breaking module R.sub.13 and R.sub.31 if the fault has occurred on current transmission line L.sub.13, R.sub.23 and R.sub.32 if the fault has occurred on current transmission line L.sub.23, or R.sub.12 and R.sub.21 if the fault has occurred on current transmission line L.sub.12. As will be set out hereafter, the fault identification modules 20 that are integrated into relays R.sub.i,j (i, j=1, 2, 3) function independently from each other, at their own location.

[0138] The fault identification module 20 comprises a current measurement input 22 and a voltage measurement input 24. Therefore current and voltage measurement are realized at one location, here for example within relays R.sub.i,j (i, j=1, 2, 3). Each input is for example connected to dedicated sensors that are known in the art in particular for HVDC grids. The sensors may comprise optical, electrical or other sensors. They may function based on diverse physical effects, like for example Hall effect, Faraday effect, Zeeman or Stark effect in order to cite only some sensors suitable to furnish relevant current and voltage measurements in a high voltage DC environment.

[0139] The fault identification module 20 can be considered as an autonomous local unit and comprises for example one or several memory 26 for storing a computer program configured to calculate on the basis of preregistered parameters, like line parameters, station parameters and parameters relating to the grid topology and in particular the physical parameters of the current transmission lines L.sub.12, L.sub.13 or L.sub.23 (conductor diameter, resistivity of the line per unit length, shielding parameters etc) and the current and voltage measurements received from the input the fault parameters. The memory(-ies) can also memorize the preregistered parameters. The preregistered parameters can be determined by calculation, specific measurements, for example before entering in service of the grid, or by calibration, for example in running the grid 100 according to a specific calibration protocol or according to a nominal service.

[0140] The fault identification module 20 comprises in addition for example one or several processors 28 for executing the computer program for calculation of the fault parameters.

[0141] In case the fault identification module 20 establishes the presence of a fault for example on the current transmission line L.sub.12, L.sub.13 or L.sub.23 to which it is dedicated, the identification module 20 can be configured to output a control signal for example on output 30, in particular to activate a specific breaking module, for example DC breaking module R.sub.i,j (i, j=1, 2, 3).

[0142] It is important to notice that the fault identification module only needs for reliable functioning the current and voltage measurement at one location, like for example at one end of the current transmission lines L.sub.12, L.sub.13 or L.sub.23, where other known methods rely on measurements that are conducted on both ends of the transmission line. This is quite an advantage, in particular for allowing much faster detection as there is no need for consolidation of measurements realized at different locations.

[0143] FIG. 9A shows an example of a simplified flowchart of the method and computer program according to the invention that may be executed by the processor of the fault identification module 20 and FIG. 9B shows a synthetic diagram of inputs and outputs of operations according to the simplified flowchart of FIG. 9A.

[0144] In step 600, the method starts. This start can be triggered for example when an unusual behavior with respect to the nominal behavior of the grid 100 or of the transmission lines L.sub.12, L.sub.13 or L.sub.23 is detected for example via sensors like current sensors i (i, j=1, 2, 3, i≠j) or voltage sensors ν.sub.i (i=1, 2, 3) which are symbolized in FIG. 9B or 10.

[0145] Such unusual behavior may be detected using a threshold on the voltage drop, on the rate of change of the voltage, or on the increase of the current, or a combination of both.

[0146] In one embodiment, one may monitor in a sliding window voltage ν and current i, for determining such unusual behavior.

[0147] In another embodiment, one may consider thresholds that are for example exceeded.

[0148] Then in step 602, the fault detection module 20 is fed with the current and voltage transients (current and voltage evolutions) that are measured for example via sensors like current sensors i.sub.i,j (i, j=1, 2, 3, i≠j) or voltage sensors ν.sub.i (i=1, 2, 3).

[0149] Typical voltage and current transients are shown in FIG. 6. It is clear that the sampling of the current and voltage measurements must be fast enough to accurate represent the evolution of current and voltage. As can be seen on FIG. 6, a transient may last 600 μs and the first voltage drop occurs in a window of less than 100 μs. Thus sampling should be for example as fast 10 μs. Due to the fast dynamics of the traveling waves, a sampling frequency of 500 to 1000 kHz is required.

[0150] In step 604, in using the physical model described above or the extended physical model taking into account soil resistivity (FIG. 7A or 7B), calculations are carried out by the fault detection module 20 in order to determine the unknown fault parameters for example d.sub.f, R.sub.f. Indeed, Given a value of fault parameters p, one can compare the outputs from the analytical model embedded in the relay, and the available measurements, as shown in FIG. 9B.

[0151] To do this, the method and the computer program executed by the processor 28 may for example use an appropriate cost function that evaluates the mismatch between the model outputs {circumflex over (ν)}.sub.x,e(p; t); {circumflex over (.Math.)}.sub.x,e(p; t) and the measurements (ν.sub.x,e; i.sub.x,e).

[0152] The relay x (for example DC breaking module R.sub.13) and line e (for example current transmission line L.sub.12) being fixed, they are omitted in the subscript in the following section. One possible choice of cost function uses the least squares:

[00023]$C^n\left(p\right)=\lambda \underset{k=1}{\overset{n}{.Math.}}{\left(v\left(t_k\right)-\hat{v}\left(p,t_k\right)\right)}^2+\mu \underset{k=1}{\overset{n}{.Math.}}{\left(i\left(t_k\right)-\widehat{\mathrm{.Math.}}\left(p,t_k\right)\right)}^2$$C^n\left(p\right)={F^{\left(n\right)}\left(p\right)}^T{F}^{\left(n\right)}\left(p\right)$

[0153] where n is the number of available data samples, and λ,μ are weighting factors between current and voltage. Either λ or μ can be set to 0 if one wishes to use only current or voltage measurement. In order to estimate the fault parameters, one aims at minimizing the cost function, i.e., find {circumflex over (p)}.sup.(n) such that

[00024]$\begin{array}{cc}{\hat{p}}^{\left(n\right)}=\arg \underset{p}{\min }{C}^{\left(n\right)}\left(p\right).& \left(15\right)\end{array}$

[0154] Note the number of unknown parameters in p that are to be identified can be reduced using prior knowledge on the line total length D.sub.x,y=d.sub.f,x+d.sub.f,y. For a given number of observations n, iterative optimization techniques such as Gauss-Newton or Levenberg-Marquardt may be used to find the optimum. Starting from an initial guess p.sub.0 the algorithm proceeds step by step using a correction term δ.sub.k, i.e.:

[00025]$p_{k+1}^{\left(n\right)}=p_k^{\left(n\right)}+{\delta }_k$

[0155] Moreover, measurements are added to the cost function after each step/iteration to improve the estimate. Hence the estimated parameters from the last step are used as initial guess for the next step.

[00026]$p_0^{\left(n+\Delta n\right)}=p_1^{\left(n\right)}.$

[0156] In particular, the method can begin at the very first detection of the fault and the iteration is enriched with new incoming measurements in real time. Δn is the number of additional data points at each estimation step and is considered constant. A smaller Δn means new data are added more frequently, increasing the precision of the algorithm, but the total number of iterations will also increase, making the algorithm slower. Δn should then be tuned to meet both speed and precision requirements. The algorithm proceeds until a stopping condition is reached, which is detailed hereafter.

[0157] In step 606, the method stops as soon as it has successfully converged to a satisfying estimate of the fault parameters, meaning the fault is identified to be on the protected line.

[0158] In case of successful detection of the fault on the dedicated current transmission line L.sub.12, L.sub.13 or L.sub.23, a control signal may be outputted in step 608 by the fault identification module 20 in order to activate for example the DC breaker R.sub.13.

[0159] One should emphasize that each fault detection module 20 is autonomous and relies on its own current and voltage measurements at its own location. Therefore, fault identification can be carried out with current and voltage measurements at a single location which allows quite faster detection of a fault and allows to take the necessary protection steps for example by activating a specific DC breaker for protection purposes. This is also a feature distinguishing the present method to those of the prior art that may allow to locate a fault, but that are not fast enough to activate protection devices in an electric power transmission and/or distribution system.

[0160] The method also monitors the number of iterations for example and if no convergence criterion is met (step 610), the method may stop in step 612 after reaching a maximum number of iterations, corresponding to a maximum measurement window, or after receiving communications from neighboring relays indicating the fault has been identified on an adjacent line.

[0161] Two indicators may for example be used to evaluate the convergence and consistency of the results with the hypothesis that the monitored line is faulty. This is done in step 606 after each iteration of the estimation algorithm so that the fault can be identified as fast as possible.

[0162] First the estimated parameters must belong to a certain validity domain that represents plausible values for the fault parameters such as

[00027]$D_p=\left\{\left(d_f,R_f\right)|d_{min}\le d_f\le d_{m\mathrm{ax}},R_{min}\le R_f<R_{\mathrm{ma}x}\right\}$

[0163] where (d.sub.min; d.sub.max) defines the portion of the line monitored by the relay and (R.sub.min; R.sub.max) the range of fault resistance that requires fast decision, since high value of Rf indicates the fault is not critical and more time is available to take action. Typical boundary values are d.sub.min=0 km; d.sub.max=90% D.sub.tot and R.sub.max=0Ω; R.sub.max=200Ω.

[0164] Second, a consistency test is used to evaluate the confidence one has in the results of the estimation. One possible approach used here may be the Cramér-Rao Lower Bound (CRLB) to benefit from prior knowledge on sensors and which gives a reasonable evaluation of the efficiency of the estimation process. Under the “regularity” conditions and considering the estimator is unbiased the CRLB states that the following matrix inequality holds

[00028]$\begin{array}{cc}{P}_{\mathrm{ls}}\left(p^{*}\right)\ge I^{-1}\left(p^{*}\right)& \left(16\right)\end{array}$

[0165] where P.sub.ls is the covariance of the estimator and I the Fisher information, both with respect to the true vector of the parameters p*, in practice replaced by the estimated vector {circumflex over (p)}. An estimate of the covariance of the maximum likelihood estimator is

[00029]$P_{\mathrm{ls}}={{\sigma }^2\left({\mathrm{��}}_F^T\left(p^{*}\right){\mathrm{��}}_F\left(p^{*}\right)\right)}^{-1}$

where JF is the Jacobian of the F vector as introduced in equation (14) with respect to the parameters p. The noise variance σ.sup.2 can be replaced by an unbiased estimate

[00030]${\hat{\sigma }}^2=\frac{1}{n_t-n_p}{F^{\left(n\right)}\left(p^{*}\right)}^T{F}^{\left(n\right)}\left(p^{*}\right).$

[0166] The Fisher information is given by

[00031]$I={\mathrm{��}}_F^T\left(p^{*}\right){.Math.}^{-1}{\mathrm{��}}_F\left(p^{*}\right)$

[0167] where Σ is the measurement noise covariance matrix, supposed known using prior information on sensors. As the estimation algorithm proceeds, the estimated parameters will get closer to the true value of the parameters and the estimator covariance closer to the inverse of the Fisher information. Hence we test the consistency of the estimation process using a threshold on the CRLB equation (16), for instance

[00032]$S\ge .Math.P_{\mathrm{ls}}\left(\hat{p}\right)-I^{-1}\left(\hat{p}\right).Math..$

[0168] Such threshold must be tuned to ensure both selectivity and sensitivity, i.e., to avoid false positives and true negatives. Tuning is done heuristically considering in particular limit cases where faults occurring on an adjacent line are harder to distinguish from fault occurring within the validity domain D.sub.p.

[0169] Simulation Results

[0170] In order to test the proposed algorithm in open-loop, measurements were simulated using EMT software EMTP-RV. Consider a meshed grid with single-conductor over-headlines and 3 MMC stations as shown in FIG. 10). Numerical values for the line and station physical characteristics as well as the parameters used in the model are presented in Tables I, II and IV, V. MMC are simulated with switching function of arms and a wide-band model is used for the transmission lines. Models for sensors are also included to reflect measurements errors (Table III). Numerical values used for tuning the algorithm are listed in Table VI.

TABLE-US-00001 TABLE I CHARACTERISTICS OF THE MMC STATIONS Rated power (MW) 1000 DC rated voltage 320 Arm inductance (p.u.) 0.15 Transformer resistance (p.u.) 0.001 Capacitor energy in each submodule (kJ/MVA) 40 Conduction losses of each IGBT/diode 0.001 Number of sub-modules per arm 400

TABLE-US-00002 TABLE II OVERHEAD-LINE CHARACTERISTICS L.sub.13 L.sub.12 L.sub.23 Length (km) 200 100 150 DC resistance (mΩ/km) 24 Outside diameter (cm) 4.775 Horizontal distance (m) 5 Vertical height at tower (m) 30 Vertical height at mid-span (m) 10 Soil resistivity (Ωm) 100

TABLE-US-00003 TABLE III SENSOR SPECIFICATIONS Current Voltage Accuracy class 1% 1% Bandwidth (kHz) 300 300 Sampling frequency (MHz)  1  1 Resolution (bits)  16  16

TABLE-US-00004 TABLE IV EQUIVALENT PARAMETERS OF THE MMC STATIONS Equivalent inductance (mH) 8.1 Equivalent resistance (Ω) 0.4 Equivalent capacitance (μF) 391

TABLE-US-00005 TABLE V OVERHEAD-LINE DISTRIBUTED PARAMETERS L.sub.13 L.sub.12 L.sub.23 Series resistance R (mΩ/km) 24 Series inductance L (mH/km) 1.45 Shunt capacitance C (nF/krn) 7.68 Shunt conductance G (nS/km) 0.2

TABLE-US-00006 TABLE VI ALGORITHM TUNING PARAMETERS. Parameter Value Initial point, p.sub.0 R.sub.init = 1 Ω, d.sub.init = 6 km Added data points per iteration, Δn 10 Consistency threshold, S  1 Maximum measurement window 3D/c.sub.w

[0171] To illustrate the operation of the protection algorithm, two typical use cases are studied where a pole to ground fault occurs between station 3 and 1, in the grid of FIG. 10. Consider in a first example a fault located at d*.sub.f=30 km from station 313 with an impedance of R*.sub.f=20Ω. Once the fault is detected at relay R.sub.31 (=DC breaking module with integrated fault identification module 20), the behavior of the estimation algorithm can be represented by plotting the evolution of the estimated parameters and the cost function levels, which also changes after each iteration due to new measurements being added. Hence the cost function and the estimated parameters are plotted in FIG. 11 after 5 iterations, when the method stops.

[0172] The algorithm stops and identifies the faults after 0.056 ms measurement window, after satisfying both feasibility and confidence criteria. Final estimated parameters are =25 km and =22Ω. Predicted voltage and current are plotted against simulated measurements in FIG. 12. Here the first wave seems to be sufficient to estimate the fault parameters and identify the fault. This allows ultra-fast operation but limits the precision of the estimated parameters. A longer measurement window of 0.250 ms gives the following estimated parameters =31 km and =18Ω.

[0173] To test for the remote fault case, consider a in a second example fault at d*.sub.f=140 km from station 313 on the same line with an impedance of R*.sub.f=40Ω. In that case, the relay R.sub.31 (DC breaking module with integrated fault identification module 20) identifies the fault after 0.116 ms measurement window with estimated parameters =92 km and =64Ω.

[0174] Evolution of the estimation algorithm is displayed on FIG. 13, analogue to FIG. 11 and starting from the same initial point. The corresponding transient voltage and current, analogue to FIG. 12, are displayed FIG. 14. The precision of the estimated parameters decreased compared to the previous first example considered as “close fault” case. This is largely due to the decreasing precision of soil resistivity effects with larger fault distance. As for the previous example, a longer measurement window would increase the precision. For instance a measurement window of 0.670 ms leads to estimated parameters =131 km and =39Ω.

[0175] Thus a the method according to the invention may carry out the measurements during a very short measurement window, for example less than 0.250 ms, in particular less than 150 ms and specifically less than 100 ms.

[0176] The inventors carried out extensive simulations and validated the behavior of the method on a wide range of fault characteristics.

[0177] Taking full benefit of the information contained in the first waves after fault occurrence, identification of the fault can be performed in less than 1 ms, which is sufficient to trigger for example associated DC breakers to protect the grid 100. Contrarily to existing single ended algorithm, DC reactors are not required to achieve selectivity. Extensive simulations in EMTP software showed the method is applicable on a wide range of fault cases regarding fault distance and resistance.

[0178] One understands that the method according the invention allows a very fast fault identification and localization which allows to protect the grid in a more efficient and reliable way.

[0179] The invention has been described in relationship with Delta high voltage DC grids (FIGS. 1 and 10), but the invention also applies to other high voltage DC grids, even more complex.

[0180] A DC transmission line in bipolar architecture or symmetrical monopolar architecture is composed by a multi-conductor system (e.g. the positive and negative poles) and hence there are mutual impedances and mutual admittances between the conductors. The measured phase voltages and currents can be transformed in modal circuit that are not coupled. The mathematical basis for the transformation are well known from the literature and resides in the theory of eigenvalues and eigenvectors of the impedance and admittance matrices (see for example E. W. Kimbark, “Transient overvoltages caused by monopolar ground fault on bipolar dc line: theory and simulation,” IEEE Trans. Power App. Syst., vol. PAS-89, no. 4, pp. 584-592, April 1970. The proposed algorithm can be applied taking into account the modal voltage and currents. The cost function can compare simulated modal voltage and current with measured voltage and current transformed in modal domain.

[0181] The proposed algorithm can also be applied to identify faults on a power transmission line in an AC electrical power transmission and/or distribution system.

[0182] After a short circuit fault on an AC transmission line, voltage and current transients can be measured at the line end as this has been described above with respect to the DC grid.

[0183] Two cases may be considered:

[0184] In a first case, the fault generates voltage and current travelling waves that propagate along the transmission lines. Such travelling waves can be described using Telegraphist's equations, similarly to what has been presented for a DC system here above.

[0185] It can be noted that AC line model employed to describe travelling wave propagation is the same as the one presented for the DC line. In order to use the single phase Telegraphist's equations as explained hereabove, it is necessary to decouple the n-phases (normally 3 phases) of the AC system using a transformation matrix as known from literature, see for example M. Bollen, “On travelling-wave-based protection of high voltage networks,” Ph.D. dissertation, TU Eindhoven, Eindhoven, The Netherlands, 1989.

[0186] High frequency model (1-10 kHz) of the AC substation equipment are also well know from literature. As an example a transformer connected to the AC substation could be represented by a simple inductance in series with a resistance.

[0187] Pre-fault voltage is measured at the line end. In case of very long line, correction of the phase voltage along the line can be taken into account. This can allow a better calculation of the pre-fault voltage at the fault point.

[0188] According to a second case, concerns the pre-fault voltage at fault location is close to zero the generated travelling waves may not be detectable by the instrument transformers. In this particular case the fault identification algorithm could not be based on Telegraphist's equations. Nevertheless a physical model capable of calculate the post fault current and voltage transients could still be identified, for example considering simple equations based on concentrated model of the line.