Frequency estimation

Abstract

Determining a property of an electrical signal at a node instant at the current instant (INSTC); determining a previous function indicative of a property of the signal at the node at a previous instant; determining an error function in a previous frequency estimation at the node at a previous instant; determining a first current function indicative of the property of the signal at the node at INSTC in accordance with the determined property, the previous function, and the error function; receiving, from another node, a second current function indicative of a property of an electrical signal at the at least one other node at INSTC; combining the first current function and the second current function to produce a current combined function (CCF) indicative of a property of the signal at the node at INSTC; and estimating a current frequency of the signal at the node in accordance with the CCF.

Claims

1. A method for estimating a frequency of an electrical signal at a node in an electrical system having a plurality of nodes, the method comprising: determining, at a node of an electrical system having a plurality of nodes, at least one property of an electrical signal at the node at a current instant in time; determining a previous function indicative of one or more properties of the electrical signal at the node at a previous instant in time; determining a function of an error in a previous frequency estimation at the node at a previous instant in time; determining a current function indicative of the one or more properties of the electrical signal at the node at the current instant in time in accordance with the at least one determined property of the electrical signal at the node at the current instant in time, the previous function indicative of the one or more properties of the electrical signal at the node at a previous instant in time, and the function of the error in the previous frequency estimation; receiving, from at least one other node of the plurality of nodes, a current function indicative of one or more properties of an electrical signal at the at least one other node at the current instant in time; combining the current function indicative of one or more properties of the electrical signal at the node and the current function indicative of one or more properties of the electrical signal at the at least one other node to produce a current combined function indicative of one or more properties of the electrical signal at the node at the current instant in time; and estimating a current frequency of the electrical signal at the node in accordance with the current combined function.

2. The method according to claim 1, wherein the one or more properties of the electrical signal that the previous and current functions are indicative of includes a frequency of the electrical signal.

3. The method according to claim 2, wherein the one or more properties of the electrical signal that the previous and current functions are indicative of includes a phase associated with the electrical signal.

4. The method according to claim 1, wherein the determined property of the electrical signal is a voltage, the method further comprising converting the three-phase electrical signal into a transform-domain, wherein the current function indicative of one or more properties of the electrical signal in the transform-domain.

5. The method according to claim 1, wherein the function of an error in the previous estimated frequency is determined by multiplication of a calculated error in a previous estimated frequency by a weighting, wherein the weighting is indicative of a level of confidence in the calculated error.

6. The method according to claim 1, wherein the function of an error is a function of a phase-only error.

7. The method according to claim 1, wherein the combining is a linear combination process in the original or a transform domain.

8. The method according to claim 1, further comprising: after receiving the current function indicative of one or more properties of the electrical signal at the at least one other node, applying a node weighting factor to the received current function before performing the step of combining, wherein the node weighting factor is indicative of a confidence in the accuracy of the received current function.

9. The method according to 1, wherein the current estimated frequency {circumflex over (f)}.sub.i(k) at the node i, is determined by: ${\hat{f}}_i\left(k\right)=\frac{1}{2T}{\sin }^{-1}\left(\mathrm{funct}\left(k\right)\right)$ wherein T is a sampling interval, and funct(k) is the current combined function at the current instant in time, k.

10. The method according to claim 1, wherein the previous function indicative of one or more properties of the electrical signal at the node at a previous instant in time, the current function indicative of one or more properties of the electrical signal at the node at the current instant in time, and the function received from the at least one other node indicative of one or more properties of the electrical signal at the at least one other node at the current instant in time, each comprise: a first function and a second function, wherein the first function utilizes the at least one determined property of the electrical signal, and the second function utilizes the complex conjugate of the at least one determined property of the signal; and the step of combining comprises: combining, at the node, the first function associated with the node with the first function associated with the at least one other node to produce a current combined first function; and combining, at the node, the second function associated with the node with the second function associated with the at least one other node to produce a current combined second function; and the step of estimating comprises estimating a current frequency of the at least one electrical characteristic at the node in accordance with the current combined first function and the current combined second function.

11. The method according to claim 10, wherein the first function and the second function also utilize a function of a calculated error.

12. The method according to claim 10, wherein the first function and the second function also utilize a function of the conjugate of a complex transform of the at least one determined property of the electrical signal at the node at the previous instant in time.

13. The method according to claim 10, wherein the first function, h, at the current instant in time, k, is determined by:
h(k)=h(km)+p(e(km))q(v*(km)) wherein m can take any integer value, h(km) is a previous first function, p(e(km)) is a function of a calculated error, is a weighting factor applied to the calculated error, and q(v*(km)) is a function of the conjugate of a complex transform of the at least one determined property of the electrical signal at the node at the previous instant in time.

14. The method according to claim 10, wherein the second function, g, at the current instant in time, k, is determined by:
g(k)=g(km)+p(e(km))q(v(km)) wherein m can take any integer value, g(km) is a previous second function, p(e(km)) is a function of the calculated error, is a weighting factor applied to the calculated error at the previous instant in time, and q(v(km)) is a function of a complex transform of the at least one determined property of the electrical signal at the node at a previous instant in time.

15. The method according to claim 10, wherein the current function indicative of the one or more properties of the electrical signal at the node at the current instant in time is determined in accordance with the following state space model: State equation: $\left[\begin{array}{c}{h}_i\left(k\right)\\ g_i\left(k\right)\\ h_i^{*}\left(k\right)\\ g_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{h}_i\left(k-m\right)\\ g_i\left(k-m\right)\\ h_i^{*}\left(k-m\right)\\ g_i^{*}\left(k-m\right)\end{array}\right]+u_i\left(k\right)$ Observation equation for the node, i: $\left[\begin{array}{c}{v}_i\left(k\right)\\ v_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{v}_i\left(k-m\right)h_i\left(k\right)+v_i^{*}\left(k-m\right)g_i\left(k\right)\\ v_i^{*}\left(k-m\right)h_i^{*}\left(k\right)+v_i\left(k-m\right)g_i^{*}\left(k\right)\end{array}\right]+n_i\left(k\right)$ Wherein h.sub.i(km) is the first function indicative of the one or more properties of the electrical signal at the node, i, at a previous instant in time, km, g.sub.i(km) is the second function indicative of the one or more properties of the electrical signal at the node at the previous instant in time, km, u.sub.i(k) and n.sub.i(k) are functions indicative of an error in the previous frequency estimation, and v.sub.i(km) is the determined property of the electrical signal at the previous instant in time, km, at the node, i.

16. The method according to claim 10, wherein the current function indicative of the one or more properties of the electrical signal at the node, i, at the current instant in time is determined in accordance with the following state space model: $\left[\begin{array}{c}{h}_i\left(k\right)\\ g_i\left(k\right)\\ z_i\left(k\right)\\ h_i^{*}\left(k\right)\\ g_i^{*}\left(k\right)\\ z_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{h}_i\left(k-m\right)\\ g_i\left(k-m\right)\\ z_i\left(k-m\right)h_i\left(k-m\right)+z_i^{*}\left(k-m\right)g_i\left(k-m\right)\\ h_i^{*}\left(k-m\right)\\ g_i^{*}\left(k-m\right)\\ z_i^{*}\left(k-m\right)h_i^{*}\left(k-m\right)+z_i\left(k-m\right)g_i^{*}\left(k-m\right)\end{array}\right]+u_i\left(k\right)$ Observation equation for a node i: $\left[\begin{array}{c}{v}_i\left(k\right)\\ v_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{z}_i\left(k\right)\\ z_i^{*}\left(k\right)\end{array}\right]+n_i\left(k\right)$ Wherein h.sub.i(km) is the first function indicative of the one or more properties of the electrical signal at the node i at a previous instant in time, km, g.sub.i(km) is the second function indicative of the one or more properties of the electrical signal at the node at the previous instant in time, z.sub.i(km) is an estimation of the one or more properties of the electrical signal at the previous instant in time, km, u.sub.i(k) and n.sub.i(k) are functions indicative of an error in the previous frequency estimation, and v.sub.i(km) is the determined property of the electrical signal at the previous instant in time, km, at the node, i.

17. The method according to claim 1, wherein the current function indicative of the one or more properties of the electrical signal at the node is also determined in accordance with a conjugate of the at least one determined property of the electrical signal.

18. An apparatus for estimating a frequency, the apparatus comprising: a processing unit arranged to perform the method of claim 1.

19. A system comprising: a plurality of nodes, each node comprising a processing unit arranged to perform the method of claim 1, wherein each node of the plurality of nodes is coupled to at least one other node of the plurality of nodes.

20. A non-transitory computer readable medium comprising computer readable code operable, in use, to instruct a computer to perform the method of claim 1.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Exemplary embodiments of the invention shall now be described with reference to the drawings in which:

(2) FIG. 1 is a simplified diagram of the transmissions and distribution channels of a power grid system 100;

(3) FIG. 2 is a simplified diagram of a integrated frequency determination network; and

(4) FIG. 3 shows a structure of a node.

(5) Throughout the description and the drawings, like reference numerals refer to like parts.

SPECIFIC DESCRIPTION

(6) FIG. 1 is a simplified diagram of the transmission and distribution channels of a power grid system 100.

(7) Power is transmitted from a main power station (not shown) at a subtransmission voltage via a high voltage power line 101 and the subtransmission voltage is then received at a substation 102. The substation 102 steps the voltage down from the voltage of the high voltage power line 101 used for transmitting the power generated at the power station long distances to the substation and other substations, to a local distribution voltage delivered to an area local to the substation via a distribution voltage line 103. Dividing line 104 therefore marks a boundary between the components of the high-voltage transmission sector 105 of the system 100 and the local distribution sector 106 of the power grid system 100.

(8) The distribution voltage line 103 splits into multiple lines delivering power to different parts of the local area. In FIG. 1, the distribution voltage line splits into three sub-lines delivering power to three separate local grids 107, 108,109. The power delivered by the distribution voltage line 103 is unevenly split between the three local grids 107, 108, 109 by 20%, 50% and 30% respectively due to the demand of each local grid. For reasons of simplicity only one local grid 109 is depicted in detail. However, all of the local grids provide one or more loads that take power from the system and one or more generators that put power back into the system. Each load may be a particular building, such as a house, apartment block, office or factory, and each generator may be a local wind farm, household solar PV cell or any other suitable small-scale energy harvesting mechanism.

(9) Local grid 109 comprises a local pole transformer 110 which provides a final main power voltage step down to the standard mains power to be delivered to the various end points, which in the UK is 240V AC at 50 Hz. The local grid 109 provides three local supply lines 111, 112, 113, which each connect to various loads and/or generators. In this case local supply line 111 supplies a load 114 and a generator 115, local supply line 112 supplies three loads 116, 117, 118 and local supply line 113 supplies local supply line 119, 120. Each load L will demand different amounts of power and each generator will deliver different amounts of power and the delivery of such power will be time varying due to fluctuations in power generation due to, for example in the case of solar PV cells, the variation in sunlight intensity, or for wind turbines, the variation in wind velocity.

(10) Since the grid system of FIG. 1 has various inputs of power (i.e. generators) and outputs of power (i.e. loads) distributed throughout the system, the supply-demand power balance throughout the system will constantly be varying. The electrical signal at each node (i.e. generator/load) can be considered a separate electrical signal dependent on electrical signals at nearby nodes, or can be considered a single electrical signal that has differing characteristics at different nodes.

(11) Due to the variations in the power balance throughout the system of FIG. 1, the system runs the risk of various problems occurring such as imbalances between the generation and load, single and dual phase faults, dual character of load-supply, harmonics, and transient stability, as previously discussed. However, the grid system 100 of FIG. 1 is specifically designed to minimise the system failures that occur within the grid system 100, as discussed below.

(12) Each node is arranged to measure power system parameters such as voltages and/or currents and to estimate the frequency of the AC voltage at the node in accordance with those measurements. This frequency estimation can then be used to detect faults, for example by means of complex circularity analysis of the voltage.

(13) The system is then further designed to improve the speed and accuracy of frequency estimation by enabling each node to share information regarding its own frequency estimation with some or all of the neighbouring nodes, or remote nodes via a communication link. Each node then utilises the information received from neighbouring nodes in its own frequency estimation, as explained in detail below with reference to FIG. 2.

(14) FIG. 2 is a simplified diagram of an integrated frequency estimation network. FIG. 2 relates to specifically to local grid 109 of FIG. 1. However, it will be appreciated that the network could also include other local grids, but for ease of explanation only local grid 109 is considered in this document. The integrated frequency determination network of FIG. 2 comprises a plurality of interconnected nodes, wherein each node is any position on the local grid 109 of FIG. 1 where a measurement of the electrical characteristics at that position is taken. In general, each node will correspond to the position of a load or power source within the grid, but this is not always necessary.

(15) In general, the system of FIG. 2 works as follows. Each node is arranged to perform a frequency estimation method to determine the instantaneous frequency of the electrical signal at that node. However, in order to improve the performance and accuracy of the frequency estimation, each node shares information relating to its frequency determination with neighbouring nodes. Each node then uses the information collected directly from its own node along with information collected from neighbouring nodes regarding previous frequency estimations in order to assist in the next frequency estimation at that node. Hence, the integrated frequency determination network is a multinodal network utilising a multi-input frequency determination scheme, based on the past measurements (e.g. regression or state space approach). Such information sharing improves the accuracy and speed at which the frequency is determined and robustness against noise, harmonics, and faults in communication links.

(16) In FIG. 2, the nodes are communicatively coupled via communications links provided via the power lines, or otherwise. For example, systems such as broadband over power lines (BPL) may be utilised. However, it will be appreciated that it is envisaged that in alternative embodiments of the invention nodes within such a network may be communicatively connected in ways other than by using BPL. For example, all nodes within a network may be wirelessly connected with one another so that all nodes share information with one another. Alternatively, all nodes may be connected to each other via a direct wired network or via a central server, which routes the information between all nodes. Hence, an internet-based connection system may be utilised.

(17) Each node in the network receives information from different nodes. For example, in FIG. 2, node 114 receives information from nodes 119, 116 and 115, node 115 only receives information from node 114, node 119 receives information from node 114, 116, and 120, node 116 receives information from node 119, 114 and 117, node 117 receives information from node 116 and node 118, while node 118 receives information from node 117. Some nodes therefore receive information from more nodes than others. The nodes that receive more information are those nodes that are connected to a greater number of different nodes. However, the receipt of larger amounts of information should assist such nodes in identifying system errors more quickly. Furthermore, the local and distributed nature of the system means that each interconnected node is in fact passing information on to nodes that the node itself is not directly connected to. Hence, the network and in particular the passing of information through the network can be complex.

(18) The frequency determination performed at each node shall now be described in detail.

(19) At each node in the network a three phase voltage measurement is taken at a particular instance in time, k. Hence, at each node, three time-domain voltage measurements are obtained v.sub.a(k), v.sub.b(k) and v.sub.c(k). Each of these phase voltages can be represented mathematically in discrete time form as shown by equation 1:

(20) $\begin{array}{cc}{v}_a\left(k\right)=V_a\left(k\right)\cos \left(kT+\right)v_b\left(k\right)=V_b\left(k\right)\cos \left(kT+-\frac{2}{3}\right)v_c\left(k\right)=V_c\left(k\right)\cos \left(kT++\frac{2}{3}\right)& \mathrm{Equation}1\end{array}$
where V.sub.a(k), V.sub.b(k), and V.sub.c(k) are the voltage amplitudes of the first, second and third phase voltage signals respectively, is the angular frequency, k is the instant in time of the voltage measurement, T is the sampling period, and is the phase.

(21) It can be seen from equation 1 that each of the three voltages of the three phase voltage signal has an associated system frequency, and their phasors rotate in the complex plane according to that frequency. However, as discussed previously, determining the system frequency from any one of the three voltages does not provide an accurate and robust measure of the overall system frequency at that node at the instance in time, k. Instead, it is necessary to determine the frequencies associated with all three voltages.

(22) In order to simplify the process of determining the overall frequency of a voltage at a node at an instance in time, k, Clarke's transform is utilised, which employs the orthogonal 0 transformation to map the time-dependent three-phase voltage into a zero-sequence v.sub.0 and the direct- and quadrature-axis components v.sub. and v.sub..

(23) $\begin{array}{cc}\left[\begin{array}{c}{v}_0\left(k\right)\\ v_{}\left(k\right)\\ v_{}\left(k\right)\end{array}\right]=\sqrt{\frac{2}{3}}\left[\begin{array}{ccc}\frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\\ 1& -\frac{1}{2}& -\frac{1}{2}\\ 0& \frac{\sqrt{3}}{2}& -\frac{\sqrt{3}}{2}\end{array}\right]\left[\begin{array}{c}{v}_a\left(k\right)\\ v_b\left(k\right)\\ v_c\left(k\right)\end{array}\right]& \mathrm{Equation}2\end{array}$

(24) For a balanced system, V.sub.a(k)=V.sub.b(k)=V.sub.c(k) and thus v.sub.0(k)=0, v.sub.(k)=A cos(kT+), and v.sub.(k)=A cos(kT++/2), where v.sub.(k) and v.sub.(k) are the orthogonal coordinates of a point whose position is time variant at a rate proportional to the system frequency. In practice, only v.sub. and v.sub. are used and the resulting complex voltage signal v(k) is given by:
v(k)=v.sub.(k)+jv.sub.(k) Equation 3
where j={square root over (1)}.

(25) There is no loss in information in using this representation, and this voltage also serves as the desired signal in adaptive frequency estimation and can be calculated recuresively from:
v(k+1)=Ae.sup.j((k+1)T+)=v(k)e.sup.jT Equation 4

(26) Where the instantaneous system frequency f is represented by the phasor e.sup.jT(f=/2).

(27) In normal operating conditions, samples of v(k) are located on a circle in the complex plane with a constant radius A. For a constant sampling frequency, the probability density function of v(k) is rotation invariant since v and ve.sup.j have the same distribution for any real . This in turn means that v(k) is second order circular, and in this case, the frequency estimation can be performed adequately by a standard linear adaptive filter, based on the strictly linear model, on which techniques such as the Complex Least Mean Square (CLMS) algorithm are built, and can be summarised as:
{circumflex over (v)}(k+1)=v(k)w(k)
e(k)=v(k+1){circumflex over (v)}(k+1)
w(k+1)=w(k)+p(e(k))q(v*(k)) Equation 5

(28) Where w(k) is the weight coefficient at time instant k, {circumflex over (v)}(k+1) is the estimate of the desired signal v(k+1), e(k) is the estimation error, and is the step-size. For the CLMS, the functions above take the form p(e(k))=e(k) and q(v(k))=v(k).

(29) Comparing Equation 5 and the strictly linear model in equation 4, the system frequency can be estimated from:

(30) 0$\begin{array}{cc}\hat{f}\left(k\right)=\frac{1}{2T}{\sin }^{-1}\left(\mathrm{Im}\left(w\left(k\right)\right)\right)& \mathrm{Equation}6\end{array}$
where Im() is the imaginary part operator.

(31) While frequency can be estimated from equation 6, such an equation does not produce optimal results because it only takes into consideration the second order circularity, and is thus only suitable for balanced systems, for which their Clarke's voltage would be depicted by a circle in a circularity diagram.

(32) As such, in alternative embodiments of the invention alternative frequency estimation methods are utilised, which produce optimal or near optimal results, as discussed below.

(33) When the three-phase power system deviates from its normal condition, such as when the three channel voltages exhibit different levels of dips or transients, or the phase difference between two channels deviates from the nominal

(34) $\frac{2}{3},$
v.sub.a(k), v.sub.b(k) and v.sub.c(k) are not identical, and samples of v(k) are not allocated on a circle with a constant radius, the distribution of v(k) is rotation dependent (noncircular), and the signal is accurately expressed only by using a widely linear model, that is a model utilising both second-order circular and non-circular information (covariance and pseudocovariance), as defined by:

(35) $\begin{array}{cc}v\left(k\right)=A\left(k\right)e^{j\left(kT+\right)}+B\left(k\right)e^{-j\left(kT+\right)}\mathrm{Where}A\left(k\right)=\frac{\sqrt{6}\left(V_a\left(k\right)+V_b\left(k\right)+V_c\left(k\right)\right)}{6}\mathrm{and}B\left(k\right)=\frac{\sqrt{6}\left(2V_a\left(k\right)-V_b\left(k\right)-V_c\left(k\right)\right)}{12}-j\frac{\sqrt{2}\left(V_b\left(k\right)-V_c\left(k\right)\right)}{4}& \mathrm{Equation}7\end{array}$

(36) In other words, when V.sub.a(k), V.sub.b(k) and V.sub.c(k) are not identical, B(k) is not equal to 0, introducing a certain degree of non-circularity. Since the widely linear model is a second order optimal estimator for noncircular data, both v(k) and its complex conjugate, v*(k), should be considered in the frequency estimation in unbalanced cases. In order to improve performance, the second order non-circularity (improperness) is therefore also taken into account, i.e.:
{circumflex over (v)}(k+1)=v(k)h(k)+v*(k)g(k) Equation 8

(37) Where h(k) and g(k) are the filter weight coefficients corresponding to the standard, h, and second, g, updates at time instant k, respectively, wherein the estimation error e(k) is defined as e(k)=v(k+1){circumflex over (v)}(k+1).

(38) The weight update of the so-called Augmented Complex Least Mean Square (ACLMS) algorithm, designed for training the widely linear adaptive filters is given by:
h(k+1)=h(k)+p(e(k))q(v*(k))
g(k+1)=g(k)+q(e(k))q(v(k)) Equation 9
where p(e(k))=e(k) and q(v(k))=v(k).

(39) It is noted that there are several stochastic gradient type algorithms to train the widely linear adaptive filters, such as the normalised complex least mean square (NCLMS), affine projection (AP), least mean phase (LMP), and their sign and leaky versions, depending on the choice of function used for the functions p() and q(). The weight update of all the algorithms use Equation 9 and the specific functions p() and q() is summarised as:

(40) $\begin{array}{cc}\mathrm{NCLMS}\mathrm{algorithm}\text{:}p\left(e\left(k\right)\right)=e\left(k\right),q\left(v\left(k\right)\right)=\frac{v\left(k\right)}{{.Math.v\left(k\right).Math.}^2}\mathrm{AP}\mathrm{algorithm}\text{:}p\left(e\left(k\right)\right)=e\left(k\right),q\left(V\left(k\right)\right)=V\left(k\right){\left(V^T\left(k\right)V^{*}\left(k\right)+I\right)}^{-1}\mathrm{LPM}\mathrm{algorithm}\text{:}p\left(e\left(k\right)\right)=<v\left(k+1\right)-<\hat{v}\left(k+1\right),q\left(v\left(k\right)\right)=\frac{jv\left(k\right)}{{\hat{v}}^{*}\left(k+1\right)}& \mathrm{Equation}10\end{array}$

(41) Where {circumflex over (v)}(k+1) is the at least one property of an electrical signal determined by estimation, e(k) is the function of an error in the frequency estimation where defined as e(k)=v(k+1){circumflex over (v)}(k+1), V(k) is the input matrix, defined as:

(42) $V\left(k\right)=\left[\begin{array}{c}v\left(k\right),v\left(k-1\right),\mathrm{.Math.},v\left(k-m\right)\\ v\left(k-1\right),v\left(k-2\right),\mathrm{.Math.},v\left(k-m-1\right)\\ \mathrm{.Math.}\\ v\left(k-n\right),v\left(k-n-1\right),\mathrm{.Math.},v\left(k-m-n\right)\end{array}\right]$
is a regularisation term, which is kept to be a small positive constant, I is an identity matrix, is an angle of a complex variable, and p() and q() are functions.

(43) To introduce the corresponding ACLMS-based frequency estimation method for the three-phase unbalanced power system, by substituting Equation 7 into 8, the estimate {circumflex over (v)}(k+1) can be expressed as:

(44) $\begin{array}{cc}\begin{array}{c}\hat{v}\left(k+1\right)=A\left(k\right)h\left(k\right)e^{j\left(kT+\right)}+B\left(k\right)h\left(k\right)e^{-j\left(kT+\right)}+\\ A^{*}\left(k\right)g\left(k\right)e^{-j\left(kT+\right)}+B^{*}\left(k\right)g\left(k\right)e^{j\left(kT+\right)}\\ =\left(A\left(k\right)h\left(k\right)+B^{*}\left(k\right)g\left(k\right)\right)e^{j\left(kT+\right)}+\\ \left(A^{*}\left(k\right)g\left(k\right)+B\left(k\right)h\left(k\right)\right)e^{-j\left(kT+\right)}\end{array}& \mathrm{Equation}11\end{array}$

(45) While from Equation 7, the expression for v(k+1) can be rewritten as:
v(k+1)=A(k+1)e.sup.jTe.sup.j(kT+)+B(k+1)e.sup.jTe.sup.j(kT+) Equation 12

(46) Therefore, at the steady state, the first term on the right hand side (RHS) of Equation 12 can be estimated approximately by its counterpart in Equation 11; hence, the term e.sup.jT containing the frequency information can be estimated as:

(47) $\begin{array}{cc}{e}^{j\hat{}T}=\frac{A\left(k\right)h\left(k\right)+B^{*}\left(k\right)g\left(k\right)}{A\left(k+1\right)}& \mathrm{Equation}13\end{array}$

(48) Comparing the second term on the RHS of Equation 11 and 12, the evolution of the term e.sup.jT can be expressed as:

(49) $\begin{array}{cc}{e}^{-j\hat{}T}=\frac{A^{*}\left(k\right)g\left(k\right)+B\left(k\right)h\left(k\right)}{B\left(k+1\right)}& \mathrm{Equation}14\end{array}$

(50) Upon taking the complex conjugate, we obtain:

(51) $\begin{array}{cc}{e}^{j\hat{}T}=\frac{A\left(k\right)g^{*}\left(k\right)+B^{*}\left(k\right)h^{*}\left(k\right)}{B^{*}\left(k+1\right)}& \mathrm{Equation}15\end{array}$

(52) We here adopt the assumptions held implicitly in frequency estimation by adaptive filtering algorithms that A(k+1)A(k) and B(k+1)B(k), and thus, Equation 13 and Equation 15 can be simplified as:

(53) $\begin{array}{cc}{e}^{j\hat{}T}=h\left(k\right)+\frac{B^{*}\left(k\right)}{A\left(k\right)}g\left(k\right)& \mathrm{Equation}16\\ e^{j\hat{}T}=h^{*}\left(k\right)+\frac{A\left(k\right)}{B^{*}\left(k\right)}{g}^{*}\left(k\right)& \mathrm{Equation}17\end{array}$

(54) Since the coefficient A(k) is real-valued, whereas B(k) is complex-valued, and thus, B*(k)/A(k)=(B(k)/A(k))*. Since Equation 16 should be equal to Equation 17, using a(k)=(B(k)/A(k))*, we can find the form of a(k) by solving the following quadratic equation with complex-valued coefficients:
g(k)a.sup.2(k)+(h(k)h*(k))a(k)g*(k)=0 Equation 18

(55) Where the discriminant of this quadratic equation is given by:

(56) 0$\begin{array}{cc}\begin{array}{c}=\sqrt{{\left(h\left(k\right)-h^{*}\left(k\right)\right)}^2+4{.Math.g\left(k\right).Math.}^2}\\ =2\sqrt{-{\mathrm{Im}}^2\left(h\left(k\right)\right)+{.Math.g\left(k\right).Math.}^2}\end{array}& \mathrm{Equation}19\end{array}$

(57) Since a(k) is complex-valued, the discriminant is negative, and the two roots can be found as:

(58) $\begin{array}{cc}{a}_1\left(k\right)=\frac{-j\mathrm{��}\left(h\left(k\right)\right)+j\sqrt{{\mathrm{Im}}^2\left(h\left(k\right)\right)-{.Math.g\left(k\right).Math.}^2}}{g\left(k\right)}{a}_2\left(k\right)=\frac{-j\mathrm{��}\left(h\left(k\right)\right)-j\sqrt{{\mathrm{Im}}^2\left(h\left(k\right)\right)-{.Math.g\left(k\right).Math.}^2}}{g\left(k\right)}& \mathrm{Equation}20\end{array}$

(59) From Equation 16, the phasor e.sup.j{circumflex over ()}T is estimated either by using h(k)+a.sub.1(k)g(k) or h(k)+a.sub.2(k)g(k). Since the system frequency is far smaller than the sampling frequency, the imaginary part of e.sup.j{circumflex over ()}T is always positive, thus excluding the solution based on a.sub.2(k). The system frequency {circumflex over (f)}(k) is therefore estimated in the form:

(60) $\begin{array}{cc}\hat{f}\left(k\right)=\frac{1}{2T}{\sin }^{-1}\left(\mathrm{Im}\left(h\left(k\right)+a_1\left(k\right)g\left(k\right)\right)\right)& \mathrm{Equation}21\end{array}$

(61) It is noted that, in addition to equation 21, there are several functions to obtain the estimated system frequency, and our aim here is to claim any functions involving h(km) and g(km), where m can take any integer value, that can be used to calculate the frequency estimate.

(62) The above described systems for frequency calculation relate to frequency calculation at a single node, without interaction between the nodes. The improved global frequency determination shall therefore now be discussed in detail.

(63) Whether using the widely linear or state space models, at each node it is possible to derive the functions h and g from the previous error, complex voltage samples, sampling time interval T and the phase . In particular, the h and g parameters of each node i are determined from Equation 9 for the widely linear model or from any one of Equations 22, 26, 28, 30 and 32 for the state space models.
h.sub.i(k+1)=h.sub.i(k)+p(e.sub.i(k))q(v.sub.i*(k))
g.sub.i(k+1)=g.sub.i(k)+p(e.sub.i(k))q(v.sub.i(k)) Equation 22

(64) Where h.sub.i(k) and g.sub.i(k) are the h and g estimate at node i at time k, e.sub.i(k) is an error calculated using previous error estimates of the parameters h and g together with the signals observed or measured at node i, namely v.sub.i(km) where the value of m can be any integer, and functions p and q are as in equations 9 and 10. Different methods for estimating h.sub.i(k) and g.sub.i(k) can be achieved, such as the normalised least mean square (NLMS) and least mean phase (LMP), depending on the choice of function used for the functions p() and q(). The variables h.sub.i(k) and g.sub.i(k) are the diffused or global h and g estimates at node i at time instant k, and are computed as:

(65) $\begin{array}{cc}{\underset{\_}{h}}_i\left(k\right)=\underset{l{N}_i}{.Math.}c\left(i,l\right)h_l\left(k\right){\underset{\_}{g}}_i\left(k\right)=\underset{l{N}_i}{.Math.}c\left(i,l\right)g_l\left(k\right)& \mathrm{Equation}23\end{array}$

(66) Where all the h estimates received by node i, as well as the h estimate from node i itself are weighted using the node dependent weighting factors or coefficients c(i,l), while the expression N.sub.i represents the neighbourhood of node i, which consists of node i and all the other nodes that can communicate or are connected with node i. It is appreciated that the weighing factors c(i,l) used for computing the diffused variable h.sub.i(k) can either be the same as or different from those used for g.sub.i(k). Moreover, the diffused variables h.sub.i(k) and g.sub.i(k) can be alternatively computed from the h and g estimates through nonlinear mappings such that:
h.sub.i(k)=r(h.sub.N.sub.i) Equation 24
g.sub.i(k)=s(g.sub.N.sub.i) Equation 25

(67) Wherein r() and s() are linear or nonlinear functions, and h.sub.N.sub.i is a vector consisting of all the h estimates from the neighbourhood of node i, similarly, g.sub.N.sub.i is a vector consisting of all the g estimates from the neighbourhood of node i.

(68) In other words, when each node has determined its local h and g value it then shares this information with each of the connected nodes in the network. Each node then receives h and g values from any nodes with which the respective node is connected. Each node then determines a global h and g, based on each of the local h and g values it has received and the h and g value that it has determined.

(69) In practice, each node combines the received h values and received g values independently. A set of weighting factors, one associated with each local node, including itself, are applied to each h and g value in the summation process. The weighting factors or coefficients c(i,l) are indicative of the perceived accuracy of the information associated with each node. Hence, some nodes which provide accurate frequency estimation will be heavily weighted in the summation, while others, which do not provide accurate frequency estimation, are less heavily weighted in the summation. A weighting factor of zero could be utilised to discard information from a particular node. The node dependent weightings may be continually updated based upon error determinations.

(70) More specifically, the weighing factors or coefficients c (or r and s) from each node i in the multi-node network can be chosen from a number of different possible schemes. These schemes include, but are not limited to: Methods for which the weighing factors are chosen according to the distributed network topology, such as the nearest neighbourhood method, which involve setting the weighing factors according to the number of nodes each node can communicate with; Methods involving setting weighing factors according to errors, be they functions of the estimated or calculated instantaneous errors, or predetermined errors; or Methods involving using predetermined weighing factors without any estimation or calculation.

(71) In an alternative embodiment of the invention the variables, h.sub.i and g.sub.i are estimated using state space modelling techniques, including Kalman and particle filters. Within such frameworks, the estimation problem is described by state space models, which consists of state and observation equations, after which the actual algorithm is implemented using the state space model. The state equation describes how the variables to be estimated (h.sub.i and g.sub.i) evolve with time, while the observation equation describes how the measurements (v.sub.i) are related to the state variables (h.sub.i and g.sub.i).

(72) Any one of the following state space models could be utilised for the determination of the variables h.sub.i and g.sub.i.

(73) State Space Model 1:

(74) State equation:

(75) $\begin{array}{cc}\left[\begin{array}{c}{h}_i\left(k\right)\\ g_i\left(k\right)\\ h_i^{*}\left(k\right)\\ g_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{h}_i\left(k-m\right)\\ g_i\left(k-m\right)\\ h_i^{*}\left(k-m\right)\\ g_i^{*}\left(k-m\right)\end{array}\right]+u_i\left(k\right)& \mathrm{Equation}26\end{array}$

(76) Observation equation for a node i in a network:

(77) $\begin{array}{cc}\left[\begin{array}{c}{v}_i\left(k\right)\\ v_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{v}_i\left(k-m\right)h_i\left(k\right)+v_i^{*}\left(k-m\right)g_i\left(k\right)\\ v_i^{*}\left(k-m\right)h_i^{*}\left(k\right)+v_i\left(k-m\right)g_i^{*}\left(k\right)\end{array}\right]+n_i\left(k\right)& \mathrm{Equation}27\end{array}$

(78) Where v.sub.i(k) is Clarke's transformed voltage at node i at time instant k, the variables u.sub.i(k) and n.sub.i(k) are noise models with possibly varying statistics, which can be used to model inaccuracies in the system, as well as to set the convergence rates of algorithms (including Kalman filters), while the symbol (.)* is the complex conjugate operator. The time index variable m can take multiple values but is typically set to 1.

(79) State Space Model 2:

(80) State equation:

(81) $\begin{array}{cc}\left[\begin{array}{c}{h}_i\left(k\right)\\ g_i\left(k\right)\\ z_i\left(k\right)\\ h_i^{*}\left(k\right)\\ g_i^{*}\left(k\right)\\ z_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{h}_i\left(k-m\right)\\ g_i\left(k-m\right)\\ z_i\left(k-m\right)h_i\left(k-m\right)+z_i^{*}\left(k-m\right)g_i\left(k-m\right)\\ h_i^{*}\left(k-m\right)\\ g_i^{*}\left(k-m\right)\\ z_i^{*}\left(k-m\right)h_i^{*}\left(k-m\right)+z_i\left(k-m\right)g_i^{*}\left(k-m\right)\end{array}\right]+u_i\left(k\right)& \mathrm{Equation}28\end{array}$

(82) Observation equation for a node i in a network:

(83) $\begin{array}{cc}\left[\begin{array}{c}{v}_i\left(k\right)\\ v_i^{*}\left(k\right)\end{array}\right]=\left[\begin{array}{c}{z}_i\left(k\right)\\ z_i^{*}\left(k\right)\end{array}\right]+n_i\left(k\right)& \mathrm{Equation}29\end{array}$

(84) The state and observation noise variables, that is variables u.sub.i(k) and n.sub.i(k), in State Space Model 2 are not necessarily the same as those in State Space Model 1, due to the differing natures of the models. State Space Model 2 contains the variable z.sub.i(k), which is used to estimate Clarke's transformed voltage v.sub.i(k). Consequently, this makes State Space Model 2 more robust to systems with increased observation noise.

(85) State Space Model 3:

(86) State equation:

(87) $\begin{array}{cc}\left[\begin{array}{c}{h}_i\left(k\right)\\ g_i\left(k\right)\end{array}\right]=\left[\begin{array}{c}{h}_i\left(k-m\right)\\ g_i\left(k-m\right)\end{array}\right]+u_i\left(k\right)& \mathrm{Equation}30\end{array}$

(88) Observation equation for a node i in a network:
v.sub.i(k)=v.sub.i(km)h.sub.i(k)+v.sub.i*(km)g.sub.i(k)+n.sub.i(k) Equation 31

(89) State Space Model 3 is essentially the strictly linear version of State Space Model 1, and is more suited to cases where the state and observation noises have circular distributions.

(90) State Space Model 4:

(91) State equation:

(92) $\begin{array}{cc}\left[\begin{array}{c}{h}_i\left(k\right)\\ g_i\left(k\right)\\ z_i\left(k\right)\end{array}\right]=\left[\begin{array}{c}{h}_i\left(k-m\right)\\ g_i\left(k-m\right)\\ z_i\left(k-m\right)h_i\left(k-m\right)+z_i^{*}\left(k-m\right)g_i\left(k-m\right)\end{array}\right]+u_i\left(k\right)& \mathrm{Equation}32\end{array}$

(93) Observation equation for a node i in a network:
v.sub.i(k)=z.sub.i(k)+n.sub.i(k) Equation 33

(94) State Space Model 4 is the strictly linear version of State Space Model 2, and is more suited to systems with circular noise distributions.

(95) The statistics of the noise models, u.sub.i(k) and n.sub.i(k), in the state space models above can be different for each state space models, and vary with time. These statistics do not necessarily need to conform to the true underlying statistics, and can be chosen, computed or estimated using a number of different methods. For example, the covariance matrix for n.sub.i(k) can be set to a multiple of the identity matrix (e.g. 100I), while simultaneously the covariance matrix for u.sub.i(k) is set to a zero matrix or a multiple of the identity matrix (e.g. 0.1I). In fact, the statistics used for these noise models, as well as that for the state estimate itself, can be varied according to calculated or estimated errors.

(96) The initialisation for all the algorithms based on the state space models and gradient decent techniques described above, are such that the h.sub.i and g.sub.i variables at each node i are set using the best estimates available for these variables at each node, which can be based on the expected nominal operating frequency at each node. However, it is appreciated that the h.sub.i and g.sub.i variables can be initialised using a number of different methods, including randomly initialisation, using predetermined values or estimated (or computed) prior to starting the algorithms.

(97) Once the global h and g functions are determined at each node, it is possible for each node to estimate the frequency of the voltage at the node. Firstly, the derivation of the frequency calculation utilised by each node will be explained.

(98) 0$\begin{array}{cc}{f}_i\left(k\right)=\frac{1}{2T}{\sin }^{-1}\left(\mathrm{Im}\left\{{\underset{\_}{h}}_i\left(k\right)+{\underset{\_}{a}}_i\left(k\right){\underset{\_}{g}}_i\left(k\right)\right\}\right)& \mathrm{Equation}34\end{array}$

(99) Where T is the sampling period, Im{.} is the imaginary part of a complex number, and

(100) ${\underset{\_}{a}}_i\left(k\right)=\frac{-\mathrm{jIm}\left\{{\underset{\_}{h}}_i\left(k\right)\right\}+j\sqrt{{\mathrm{Im}}^2\left\{{\underset{\_}{h}}_i\left(k\right)\right\}-|{\underset{\_}{g}}_i\left(k\right){|}^2}}{{\underset{\_}{g}}_i\left(k\right)}.$

(101) A global frequency estimator, consisting of the all the individual estimates by all the nodes in the network, can then be derived as:

(102) $\begin{array}{cc}{f}_{\mathrm{global}}\left(k\right)=\frac{1}{2T}{\sin }^{-1}\left(\mathrm{Im}\left\{{\underset{\_}{h}}_{\mathrm{global}}\left(k\right)+{\underset{\_}{}}_{\mathrm{global}}\left(k\right){\underset{\_}{g}}_{\mathrm{global}}\left(k\right)\right\}\right){\underset{\_}{}}_{\mathrm{global}}\left(k\right)=\frac{\begin{array}{c}-\mathrm{jIm}\left\{{\underset{\_}{h}}_{\mathrm{global}}\left(k\right)\right\}+\\ j\sqrt{{\mathrm{Im}}^2\left\{{\underset{\_}{h}}_{\mathrm{global}}\left(k\right)\right\}-|{\underset{\_}{g}}_{\mathrm{global}}\left(k\right){|}^2}\end{array}}{{\underset{\_}{g}}_{\mathrm{global}}\left(k\right)}\mathrm{Where}{\underset{\_}{h}}_{\mathrm{global}}\left(k\right)=\frac{1}{L}\underset{i=1}{\overset{L}{.Math.}}{\underset{\_}{h}}_i\left(k\right)\mathrm{and}{\underset{\_}{g}}_{\mathrm{global}}\left(k\right)=\frac{1}{L}\underset{i=1}{\overset{K}{.Math.}}{\underset{\_}{g}}_i\left(k\right)& \mathrm{Equation}35\end{array}$

(103) It is noted that, similarly to equation 21, in equations 34 and 35 there are several functions to obtain the estimated system frequency, and our aim here is to claim any functions involving h(km) and g(km), where m can take any integer value, that can be used to calculate the frequency estimate.

(104) It will be appreciated that the above-mentioned method for utilising resources from multiple nodes in order to provide an improved frequency estimation could equally utilise a strictly linear model, as discussed previously. When such a strictly linear model is used the weight coefficient received from each node is combined to obtain a combined function. The combined function may be determined in a similar way to the methods described above. Finally, the frequency can then be determined based on this combined function.

(105) The structure of the frequency estimation apparatus at example node 200 shall now be described with reference to FIG. 3.

(106) Node 200 receives power from an supply power cable 201 at power input 202, and the received power is then output to the local area, which may be a building or such like, via the power output 203 and the local power distribution line 204. Sensor 205 measures a voltage at the node 200 and supplies this node to processor 206, which performs the frequency estimation. Memory 207 is provided for use during the frequency estimation process. A communication unit 208 is then coupled to the processor and to the internal power connection 209 in order to transmit frequency estimation information at node 208 to neighbouring nodes, and receive frequency estimation information from neighbouring nodes.

(107) While the frequency estimation at node 200 has been depicted as a standalone unit, it will be appreciated that such functionality could be incorporated within a standard smart meter. Such smart meters include the processor, memory, sensor, and communication unit required by the frequency estimation apparatus. Consequently, the frequency estimation process could be incorporated in software downloaded onto a smart meter. The smart meter hardware could then implement the frequency estimation software. Hence, the various methods described above may be implemented by a computer program, as discussed below.

(108) The computer program may include computer code arranged to instruct a computer to perform the functions of one or more of the various methods described above. A smart meter may be considered a computer in such circumstances. The computer program and/or the code for performing such methods may be provided to an apparatus, such as a computer or smart meter, on a computer readable medium. The computer readable medium could be, for example, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, or a propagation medium for data transmission, for example for downloading the code over the Internet. Non-limiting examples of a physical computer readable medium include semiconductor or solid state memory, magnetic tape, a removable computer diskette, a random access memory (RAM), a read-only memory (ROM), a rigid magnetic disc, and an optical disk, such as a CD-ROM, CD-R/W or DVD.

(109) An apparatus such as a computer may be configured in accordance with such computer code to perform one or more processes in accordance with the various methods discussed above.

(110) It will be appreciated that any of the state space models or regressions models described above may be utilised for estimating the frequency of any signal and not just an electrical signal. Furthermore, it will be appreciated that any of the state space models and regression models based on equation (10) described above can be utilised for estimating frequency at a single node without receiving and utilising information from neighbouring nodes.