Phase identification of power distribution systems

12607662 ยท 2026-04-21

Assignee

Inventors

Cpc classification

International classification

Abstract

Method and system for phase detecting of hybrid wye-connected and delta-connected loads in an unbalanced power distribution system. A data-driven event-triggered phase identification algorithm is presented in which the exceeding current fluctuation events are used to build the relationship between the feeder head and the undetermined phase load and mutual information index is used to assess the dependence level of the feeder head on the undetermined phase load. The connected and non-connected phase of wye-connected and delta-connected single-phase load are determined as the phase having highest and lowest dependences of the phases of feeder head on the load, respectively.

Claims

1. A phase detecting system for detecting phases of loads in a power distribution system, comprising: a signal interface configured to receive measurements from sensors connected to feeder heads and loads on the power distribution system, wherein the measurements include a real power, a reactive power, and a voltage magnitude for each phase of the feeder heads and each phase of the loads, wherein the loads include a set of loads with undetermined phases; at least one processor; and a non-transitory computer-readable medium having stored thereon a set of instructions for detecting phases of the loads, which when performed by the at least one processor, causes the at least one processor to performs steps of: storing, by selecting from the measurements, for a set of time intervals, current magnitude timeseries for each phase of the feeder heads and each phase of the loads with undetermined phases in a memory, wherein the current magnitude timeseries include current magnitudes at measured time intervals, wherein the current magnitudes are computed from the real powers, the reactive powers, and the voltage magnitudes, wherein each of the real powers, each of the reactive powers, and each of the voltage magnitudes corresponds to an interval among the set of measuring intervals; selecting a load indicating a greatest current magnitude from the set of loads with undetermined phases as a study load by a statistic detecting algorithm; determining exceeding current fluctuation events from the current magnitude timeseries at each phase of the feeder heads and the study load based on a threshold current range; computing a first marginal probability distribution for the exceeding current fluctuation events on each phase of the feeder heads, a second marginal probability distribution for the exceeding current fluctuation events at the study load, and a joint probability distribution for joint exceeding current fluctuation events for each combination of phases of the feeder heads and the study load; determining a phase for the study load based on dependence of each phase of the feeder heads on the study load according to mutual information index determined by the first and second marginal probability distributions of the exceeding current fluctuation events and the joint probability distribution of the joint exceeding current fluctuation events; updating the current magnitude timeseries at the feeder heads by a current disaggregation algorithm for equivalent current of the study load based on measured real powers, reactive powers and voltage magnitudes at the study load, and removing the study load from the set of loads with undetermined phases stored in the memory; and determining if there is a load indicating an undetermined phase in the set of loads with undetermined phases in the memory, then continuing the steps of the selecting through the updating until the set of loads with undetermined phase becomes empty, otherwise generating and outputting a dataset of the determined phases for loads.

2. The phase detecting system of claim 1, wherein the receiving is initiated by transmitting a command signal to the sensors using the signal interface.

3. The phase detecting system of claim 1, wherein the set of time intervals is defined as a collection of measuring intervals of smart meters for predefined days and a predefined range of hours of a day based on sampling rates of the sensors.

4. The phase detecting system of claim 1, wherein the power distribution system includes single-phase loads, two-phase loads, three-phase loads, or a combination thereof.

5. The phase detecting system of claim 1, the set of loads with undetermined phases include single phase loads, equivalent single phase loads converted from two-phase loads, and equivalent single phase loads converted from three-phase loads.

6. A phase detecting method for detecting phases of loads in a power distribution system, comprising steps of: receiving measurements from sensors via a signal interface, wherein the sensors are connected to feeder heads and loads on the power distribution system, wherein the measurements include a real power, a reactive power, and a voltage magnitude for each phase of the feeder heads and each phase of the loads, wherein the loads include a set of loads with undetermined phases; storing, by use of the measurements, for a set of time intervals, current magnitude timeseries for each phase of the feeder heads and each phase of the loads with undetermined phases in a memory, wherein the current magnitude timeseries include current magnitudes at measured time intervals, wherein the current magnitudes are computed from the real powers, the reactive powers, and the voltage magnitudes, wherein each of the real powers, each of the reactive powers, and each of the voltage magnitudes corresponds to an interval among the set of measuring intervals; selecting a load indicating a greatest current magnitude from the set of loads with undetermined phases as a study load by a statistic detecting algorithm; determining exceeding current fluctuation events from the current magnitude timeseries at each phase of the feeder heads and the study load based on a threshold current range; computing a first marginal probability distribution for the exceeding current fluctuation events on each phase of the feeder heads, a second marginal probability distribution for the exceeding current fluctuation events at the study load, and a joint probability distribution for joint exceeding current fluctuation events for each combination of phases of the feeder heads and the study load; determining a phase for the study load based on dependence of each phase of the feeder heads on the study load according to mutual information index determined by the first and second marginal probability distributions of the exceeding current fluctuation events and the joint probability distribution of the joint exceeding current fluctuation events; updating the current magnitude timeseries at the feeder heads by a current disaggregation algorithm for equivalent current of the study load based on measured real powers, reactive powers and voltage magnitudes at the study load, and removing the study load from the set of loads with undetermined phases stored in the memory; and determining if there is a load indicating an undetermined phase in the set of loads with undetermined phases in the memory, then continuing the steps of the selecting through the updating until the set of loads with undetermined phase becomes empty, otherwise generating and outputting a dataset of the determined phases for loads.

7. The phase detecting method of claim 6, wherein the receiving is initiated by transmitting a command signal to the sensors using the signal interface.

8. The phase detecting method of claim 6, wherein the set of time intervals is defined as a collection of measuring intervals of smart meters for predefined days and a predefined range of hours of a day based on sampling rates of the sensors.

9. The phase detecting method of claim 6, wherein the power distribution system includes single-phase loads, two-phase loads, three-phase loads, or a combination thereof.

10. The phase detecting method of claim 9, wherein the set of loads with undetermined phases are single-phase loads connected between a phase and a neutral when a wye-connection is applied, or connected between two phases when a delta-connection is applied.

11. The phase detecting method of claim 6, the set of loads with undetermined phases include single phase loads, equivalent single phase loads converted from two-phase loads, and equivalent single phase loads converted from three-phase loads.

12. The phase detecting method of claim 6, wherein the statistic detecting algorithm selects the load indicating a greatest current magnitude based on sample means of current magnitudes from the set of loads with undetermined phases.

13. The phase detecting method of claim 6, wherein the current disaggregation algorithm calculates each of the current magnitude timeseries of the feeder heads at phase interval t, F D R ( t ) as absolute value of a corresponding complex current at phase interval t, . FDR ( t ) , FDR ( t ) = .Math. "\[LeftBracketingBar]" . FDR ( t ) .Math. "\[RightBracketingBar]" ; wherein complex current, . F D R ( t ) is calculated by dividing a conjugate of complex power at phase interval t by a conjugate complex voltage at phase interval t, . FDR ( t ) = P FDR ( t ) - jQ FDR ( t ) V FDR ( t ) e - j V ; wherein complex power is defined as P F D R ( t ) + j Q F D R ( t ) , P F D R ( t ) and Q F D R ( t ) are measured real power and reactive power at phase interval t; wherein the complex voltage is defined by a magnitude V F D R ( t ) and voltage phase angle on phase , .sub.V.sub.; wherein V F D R ( t ) is a measured voltage magnitude at phase interval t; wherein .sub.V.sub., is set as 0 when =A, 120, when =B and 120 when =C; wherein A, B and C are candidate wye-connected phase for the power distribution system.

14. The phase detecting method of claim 13, wherein the current disaggregation algorithm updates the current magnitude of the feeder heads at phase interval t, F D R ( t ) by subtracting the contribution of current of the study load k at phase interval t when load k is wye-connected and its phase is determined as ; wherein the updated current magnitude of phase interval t, F D R ( t ) is calculated as absolute value of a corresponding updated complex current at phase interval t, F D R ( t ) , F D R ( t ) = .Math. "\[LeftBracketingBar]" F D R ( t ) .Math. "\[RightBracketingBar]" ; wherein the updated complex current is calculated as deducting equivalent current at load k from the latest complex current of feeder head on phase at interval . FDR ( t ) . FDR ( t ) - P k LD ( t ) - jQ k LD ( t ) V k LD ( t ) e - j V ; wherein P k L D ( t ) , Q k L D ( t ) and V k L D ( t ) are measured real power, reactive power and voltage magnitude at load k and interval t.

15. The phase detecting method of claim 13, wherein the current disaggregation algorithm updates the current magnitudes of the feeder heads at phase .sub.1 and phase .sub.2 of interval t, 1 F D R ( t ) and 2 F D R ( t ) by subtracting the contributions of currents of the study load k when load k is delta-connected and its phases are determined as .sub.1 and .sub.2; wherein the updated current magnitude of phase .sub.1 and phase .sub.2 of interval t, 1 F D R ( t ) and 2 F D R ( t ) is calculated as absolute values of corresponding updated complex currents at phase .sub.1 and phase .sub.2 of interval t, . 1 F D R ( t ) and . 2 F D R ( t ) , 1 F D R ( t ) = .Math. "\[LeftBracketingBar]" . 1 F D R ( t ) .Math. "\[RightBracketingBar]" , 2 F D R ( t ) = .Math. "\[LeftBracketingBar]" . 2 F D R ( t ) .Math. "\[RightBracketingBar]" ; wherein the updated complex currents are calculated as deducting the contribution of equivalent currents at load k from the latest complex currents of feeder head on at phase .sub.1 and phase .sub.2 of interval t, . 1 F D R ( t ) . 1 F D R ( t ) - P k LD ( t ) - jQ k LD ( t ) V k LD ( t ) e - j V 1 2 , . 2 F D R ( t ) . 2 F D R ( t ) - P k LD ( t ) - jQ k LD ( t ) V k LD ( t ) e - j V 1 2 ; wherein P 1 F D R ( t ) , Q 1 F D R ( t ) and V 1 F D R ( t ) , and P 2 F D R ( t ) , Q 2 F D R ( t ) and V 2 F D R ( t ) are the measured real powers, reactive powers and voltage magnitudes at feeder head phases .sub.1 and .sub.2 at interval t, respectively; wherein .sub.V.sub.12 is a phase angle for voltage between phase .sub.1 and phase .sub.2, and set as 30 when .sub.1.sub.2=AB, 90 when .sub.1.sub.2=BC, and 150 when .sub.1.sub.2=CA; wherein AB, BC and CA are candidate delta-connected phases for the power distribution system.

16. The phase detecting method of claim 6, wherein the connected phase of load k, .sub.ke is determined as the phase with highest dependence of the phase of the feeder head on load k among three phases ={A, B, C} if load k is a wye-connected single-phase load: ke = arg ( max ) MI ( ; k ) wherein MI(; k) is an index to assess the dependence of phase of feeder head on the load k.

17. The phase detecting method of claim 6, wherein the connected phases of load k, .sub.ke are determined as the remaining phases by removing the phase with lowest dependence of feeder head on load k from the full set of three phases ={A, B, C} if the load is a delta-connected single-phase load: ke = - arg ( min ) MI ( ; k ) wherein MI(; k) is an index to assess the dependence of phase of feeder head on the load k.

18. The phase detecting method of claim 6, wherein the dependence of feeder head on load k for a given phase is defined using a mutual information index, MI(; k) calculated by: MI ( ; k ) = .Math. u { 1.2 } .Math. v { 1.2 } p , k ( u , v ) log p , k ( u , v ) p ( u ) p k ( v ) , wherein u represents an exceeding current fluctuation event at feeder head on phase , v represents an exceeding current fluctuation event at load k; wherein 1 stands for an upper-bound exceeding current fluctuation event, and 2 stands for an lower-bound exceeding current fluctuation event; p.sub.(u) is the marginal probability of exceeding current fluctuation events u occurred at the feeder head on phase , p.sub.k(v) is the marginal probability of exceeding current fluctuation events v occurred at load k, p.sub.,k(u, v) is the joint probability of exceeding current fluctuation events occurred at both phase of feeder head and load k.

19. The phase detecting method of claim 6, wherein the exceeding current fluctuation event is occurred at a time interval t at phase of feeder head if the differential current at feeder head F D R ( t ) is beyond a range of thresholds, [ L k n , L k p ] defined for a given load k; wherein an upper-bound exceeding current fluctuation event is occurred if F D R ( t ) L k p ; wherein a lower-bound exceeding current fluctuation event is occurred if F D R ( t ) L k n ; wherein differential current at time interval t, F D R ( t ) is the deviation of current magnitudes at phase of the feeder head between interval t and interval (t1).

20. The phase detecting method of claim 6, wherein the exceeding current fluctuation event is occurred at a time interval t at load k if the differential current at load k, k L D ( t ) is beyond a range of thresholds, [ L k n , L k p ] defined for load k; wherein an upper-bound exceeding current fluctuation event is occurred if k L D ( t ) L k p ; wherein a lower-bound exceeding current fluctuation event is occurred if k L D ( t ) L k n ; wherein differential current at time interval t, k L D ( t ) is the deviation of current magnitudes at load k between interval t and interval (t1).

21. The phase detecting method of claim 6, wherein the exceeding current fluctuation event is determined based on a range of thresholds, [ L k n , L k p ] ; wherein L k p = k L D _ + ( k L D ) , L k n = k L D _ - ( k L D ) ; wherein k L D _ and ( k L D ) are the mean and the standard deviation of differential current magnitude for the time series k L D over the set of measuring intervals.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The presently disclosed embodiments will be further explained with reference to the attached drawings. The drawings shown are not necessarily to scale, with emphasis instead generally being placed upon illustrating the principles of the presently disclosed embodiments.

(2) While the following identified drawings set forth presently disclosed embodiments, other embodiments are also contemplated, as noted in the discussion. This disclosure presents illustrative embodiments by way of representation and not limitation. Numerous other modifications and embodiments can be devised by those skilled in the art which fall within the scope and spirit of the principles of the presently disclosed embodiments.

(3) FIG. 1A is a flowchart for the steps of phase identification method for a power distribution system, according to embodiments of the present disclosure;

(4) FIG. 1B is a schematic of phase identification system of a power distribution system, according to embodiments of the present disclosure;

(5) FIG. 2A is a schematic illustrating IEEE 34 bus test system;

(6) FIG. 2B is a schematic illustrating single-phase loads and three-phase loads with wye and delta connections;

(7) FIG. 3 is a schematic illustrating a typical solar irradiance variation within 24 hours of a typical day;

(8) FIG. 4A is a schematic illustrating meter readings for a bus 822 of the IEEE 34 bus test system with measurement rates as 60 readings per hour, according to embodiments of the present disclosure;

(9) FIG. 4B is a schematic illustrating meter readings for a bus 822 of the IEEE 34 bus test system with measurement rates as 12 readings per hour, according to embodiments of the present disclosure; and

(10) FIG. 5 is a schematic illustrating a comparative results obtaining by using a kNN algorithm with different number of neighbors, according to embodiments of the present disclosure.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

(11) The present invention relates generally to power distribution systems, and more particularly to resilient distribution system infrastructure planning. The following description provides exemplary embodiments only, and is not intended to limit the scope, applicability, or configuration of the disclosure. Rather, the following description of the exemplary embodiments will provide those skilled in the art with an enabling description for implementing one or more exemplary embodiments. Contemplated are various changes that may be made in the function and arrangement of elements without departing from the spirit and scope of the subject matter disclosed as set forth in the appended claims.

(12) Specific details are given in the following description to provide a thorough understanding of the embodiments. However, understood by one of ordinary skill in the art can be that the embodiments may be practiced without these specific details. For example, systems, processes, and other elements in the subject matter disclosed may be shown as components in block diagram form in order not to obscure the embodiments in unnecessary detail. In other instances, well-known processes, structures, and techniques may be shown without unnecessary detail in order to avoid obscuring the embodiments. Further, like reference numbers and designations in the various drawings indicated like elements.

(13) FIG. 1A is a flowchart for the steps of phase identification method for a power distribution system, according to embodiments of the present disclosure.

(14) The method 100A comprises of several steps as described below. Step 120 identify the a set of loads without phase labels, and collects readings of smart meters at the feeder head, and customer loads with unknown phase labels; Step 122 chooses the load with heaviest current among all loads with undetermined phases as the study load; Step 124 computes marginal probabilities and joint probability of current exceeding fluctuations for each phase between the study load and the feeder head based on corresponding differential current data; Step 126 determines the phase for the study load based on mutual information index determined based marginal and joint probabilities between the feeder head and the study load; Step 128 updates feeder head currents by deducting the contribution of currents of the study load from feeder head currents, and remove the study load from the set of loads with undetermined phases; Step 130 checks if the set of loads with undetermined phases is empty. If yes, go to step 132. Otherwise, go to step 122. Step 132 outputs the results of phase identification to distribution management system, or distribution control system.

(15) FIG. 1B is a schematic of phase identification system of a power distribution system, according to embodiments of the present disclosure.

(16) The phase identification system 100B includes a human machine interface (HMI) 167 connectable with a keyboard 111 and a pointing device/medium 112, a processor 155, a storage device 154, a memory 137, a network interface controller 163 (NIC) connectable with a network 151 including local area networks and internet network, a display interface 161 connected to a display device 165, an input interface 139 connectable with an input device 135, a printer interface 133 connectable with a printing device 131. Further, the NIC 163 includes a signal interface configured to receive measurements 195 from smart sensors 105 arranged in the power distribution system 115 through the communication 180. The memory 137 is configured to load the phase identification program 159 by associating with the storage device 154 when executing the method implemented in 100A. In some cases, the memory 137 and the storage device 154 may be referred to as a memory.

(17) The phase identification system can receive electric signals 195 indicating timeseries measurements of currents of loads arranged in a power distribution system 115 via the network 151 connected to the NIC 163. The network 151 is connected to an outside system(s) 101 that can provide control signals to the power distribution system 115 and associated smart meters 105 installed at the loads 107 in the power distribution systems 115 for performing measurement collecting of the feeder head 109 and loads 107. The feeder head may be multiple header heads 109 if multiple feeders are existing. Further, the phase identification 100B can provide the outside system 101 phase identification data (signals) via the network 151 so that the outside system 101 can trigger data collection associated with smart meters arranged in the power distribution system 115. Further, the phase identification system 100B can be controlled from the outside system 101 by receiving data (signals) of the phase identification system 100B via the network 151.

(18) The storage device 154 includes parameters 158 for phase identification of power distribution system with respect to the power distribution system 115 and a phase identification program module 159. The input device/medium 139 may include modules that read programs stored on a computer readable recording medium (not shown).

(19) For identifying connected phases of loads in the power distribution system 115, the phase identification system 100B may receive the measuring data of the loads 107 and feeder head 109 from the smart meters 105 arranged in the customer sides and the feeder head of the power distribution system 115.

(20) When the phase identification system 100B receives the command from the outside system 101 (such as distribution management system, or distribution control system) to identify all missing phase labels for all loads in the power distribution system 115. It first loads parameters 158 of phase identification from the memory storage unit 154 into the memory 137, then execute the phase identification steps described in mutual information based phase identification program 159. Before implementing the phase identification modules, the phase identification system 100B first identify a set of loads without phase labels and collecting readings of smart meters at the feeder head, and loads with unknown phase labels; then iteratively execute following steps: 1). choosing the load with heaviest current among all loads with undetermined phases as the study load, 2). computing marginal probabilities and joint probability of current exceeding fluctuations for each phase between the study load and the feeder head based on corresponding differential current data, 3). determining the phase for the study load based on mutual information index determined based marginal and joint probability between feeder head and the study load, 4). removing the study load from the set of loads with undetermined phases, and updating feeder head currents by deducting the contribution of currents of the study load from the feeder head currents, and 5). checking if the set of loads with undetermined phases is empty. If no, repeat the steps 1)-5) for next heaviest load. If yes, then the iteration process is completed, and results for phase identification are transmitted through display interface 161 to display on display device 165. Meanwhile, the phase identification results are also submitting to the outside system 101 to update its power distribution system topology model or execute other tasks based on those labeled load phase data.

Power Distribution System and Load Connections

(21) FIG. 2A is a schematic illustrating IEEE 34 bus test system. This one-line diagram shows a feeder fed from a substation 210 and then delivered powers to customer loads 220, 230, 240 through distribution lines. It contains both three-phase loads 220, and single-phase loads 230, 240. An advanced sensor or a smart meter 234 is installed at the feeder head. Each load is equipped with a smart meter 235.

(22) Three phase load 220 is a load connected to all three phases A, B, and C. Single-phase load is a load connected only one phase 230, 240. For example, load 230 is connected to single phase A, and load 240 is connected to phase B.

(23) Dependent on power distribution system configuration, the load can connect with the system using wye connection or delta connection.

(24) FIG. 2B is a schematic illustrating single-phase loads and three-phase loads with wye and delta connections. The power is fed from a substation bus 215 to customer loads 250, 260, 270 and 280. The smart meters 234 and 235 are installed at the feeder head and loads of the power distribution system 115. The left side of the transformer 265 is a three-phase three-wire system, including wires for phases A, B and C. The right side of the transformer 265 is a three-phase four-wire system, including wires for phases A, B, C and the neutral. In the figure, feeder head has installed 3 smart meters, and one for each phase. Each single-phase load also installed with a smart meter.

(25) Load 250 is a three-phase delta-connected load in which each phase load is connected between two phases. Load 260 is a single-phase delta-connected load which connected between two phases.

(26) Load 270 is a three-phase wye-connected load in which each phase load is connected between one phase and a neutral. Load 280 is a single-phase wye-connected load which connected between a phase and a neutral.

(27) Data Pre-Processing

(28) Consider a power distribution system, where smart meters are installed at a subset of the load buses. The load buses are assumed to be a mixture of wye and delta connections. Each smart meter can provide the following three measurements: real power, reactive power, and the magnitude of voltage. Accordingly, the magnitude of current of load k can be determined at each smart meter, as:

(29) I k = P k 2 + Q k 2 V k , k K ( 1 )
where I.sub.k is the magnitude of current of load k, P.sub.k and Q.sub.k are the real power and the reactive power of load k, V.sub.k is the magnitude of voltage of load k, and all values are given in per unit. In addition, all of the data are measured in time series by the smart meters on the customer primary voltage sides with a predefined sampling rate.

(30) Let K be the set of indexes for all loads without phase labels. It is worthy to mention that all two-phase or three-phase load without phase labels are converted into equivalent single-load loads to be considered. A set custom character.sup.LD is used to combine the time series currents in the given set of measuring intervals T for all the loads as

(31) L D = { I k L D ( t ) , t T , k K } | T | | K | .
|T| and |K| are the total numbers of measuring intervals and loads, respectively.

(32) Another advanced sensor on the feeder head can be used to measure the currents in three phases at the feeder head. Note that the three phases are in set ={A, B, C}. Here we also use the magnitude of currents at the feeder head and note it as

(33) I FDR ,
for different phases and put it into a set custom character.sup.FDR. Given the set of measuring intervals T, we have

(34) F D R = { I F D R ( t ) , t T , } | T | 3 .

(35) Since the differential time series data contain potentially useful information, we create the differential set custom character.sup.LD(t)=custom character.sup.LD(t)custom character.sup.LD(t1) and note it as custom character.sup.LDcustom character.sup.(|T|1)|K|. A similar differential process is also used for custom character.sup.FDR and note as custom character.sup.FDRcustom character.sup.(|T|1)3.

(36) custom character.sup.FDR can be a set of actually current measurements over time at the feeder head, or a set of equivalent currents at feeder head contributed by a specified set of loads, such as all loads except a specified load. The equivalent feeder head currents can be computed by deducting the equivalent currents contributed by the specified load from the measured feeder head currents. The equivalent currents contributed by the specified load are computed by assuming voltages at the bus that the load connected to are balanced, that is any two of three phases are 120 degree apart.

(37) The set of measuring intervals, T can strategically be selected according to the penetrations of PVs and smart meters. For example, we can use the time-series data from 1 AM to 6 AM local time to avoid noise caused by PVs.

(38) The next step is to figure out the useful information from both the customer load measurements and the feeder head measurements. There is no doubt that any customer load would have some fluctuations due to its load variation over time (usually represented by a load profile). The fluctuations also show on the feeder head on the same phase.

(39) For example, suppose a customer use an appliance for two minutes from 6:00 to 6:02 AM, the current on the customer side would show an increase fluctuation between 5:59 and 6:00 AM when the customer turns on the appliance. The current would show a decrease fluctuation between 6:02 to 6:03 when the customer turns off the appliance. In the meantime, the sensor on the source side (i.e., at the feeder head) would also capture the same fluctuations if both the customer load and the feeder head meters are synchronized on the same phase.

(40) The question is, how significant are the fluctuations in correct identification of the phase? To answer this question, we first define an event as a fluctuation beyond the range

(41) [ L k n , L k p ]
for any given k in the load set K, where:

(42) a ) L k p = k LD _ + ( k LD ) ( 2 a ) b ) L k n = k LD _ - ( k LD ) ( 2 b )

(43) Note that

(44) k LD _
is the mean and

(45) ( k LD )
is the standard deviation of signal

(46) k L D .
From statistics, the probability of a normal distribution is 66.67% in the range of mean valuestandard deviation. Since we are more interested in significant fluctuations, we choose only fluctuations that fall beyond the range

(47) 0 [ L k n , L k p ]
as events to be studied in the present invention.

(48) It should be noticed that the selection of events is crucial because not every point in

(49) k L D
is useful for phase identification. The same two minutes example of turning on an appliance can be used to explain the reason. If the customer does not turn on other appliances, the PV installed at the customer side has no power supply before sunrise, and the only appliance is working in a fixed power, then no fluctuation will be seen in both the customer side and the feeder head side during the two minutes. In such a case, no event signature can be captured even we compare the same phase from the customer load and the feeder head. Therefore, we can only analyze the load phase if the event signatures are clear in both the customer side and the feeder head side.

(50) Here, we use a relatively simple way to select events, as it was explained in the previous paragraph. Other complex methods, can be used for the purpose of selecting the events, such as by using moving average or least square error. However, from the numerical test in this invention, there does not seem to be a need to make the process of selecting the events unnecessarily complicated.

Mutual Information Among Joint Events

(51) Suppose a joint event at both customer load k and feeder head is beyond the range of

(52) [ L k n , L k p ] ,
it does not mean that the load k and feeder head are on the same phase. For example, if an event at load k shows fluctuation greater than

(53) L k p
while the same event at phase of feeder head shows fluctuation less than

(54) L k n ,
it is highly possible that customer k does not belong to phase . However, if an event at both customer load k and phase of the feeder head is greater than

(55) L k p ,
, it does not mean customer load k belongs to phase because it is possible that different customers have similar load profiles. Therefore, we use cross entropy from the information theory to solve this issue.

(56) Suppose m joint events are chosen between a pair of customer load k and feeder head phase in , we can obtain the probability table as follows:

(57) p , k ( u , v ) , ( 3 ) Where u = { 1 if m : FDR ( m ) L k p 2 if m : FDR ( m ) L k n ( 4 ) v = { 1 if m : k LD ( m ) L k p 2 if m : k LD ( m ) L k n ( 5 )

(58) Here p.sub.,k(1,1) shows that the probability of fluctuations in both feeder head phase and customer load k are greater than

(59) L k p ;
; p.sub.,k(1,2) shows that the probability of fluctuations in the feeder head at phase is greater than threshold

(60) L k p
while the fluctuations in the load k is less than the threshold

(61) L k n ;
p.sub.,k(2,1) shows the probability of fluctuations in feeder head phase is less than threshold

(62) 0 L k n
while the fluctuations in load k is greater than threshold

(63) L k p ;
and p.sub.,k(2,2) is the probability of fluctuation at less than

(64) L k n
in both feeder head at phase and customer load k.

(65) We can transform the data from the joint events m into a probability relationship, such as p.sub.,k(u, v) in equation (3), for different relationships between feeder head phase and the load k. Since equation (3) is based on probability, it is possible to have bias, which can be avoided by using enough events. In other words, the larger amount of joint events used, the more accurate the results will be.

(66) For example, we can only analyze the load phase when number of joint events m is not less than 12. Note that the predetermined period in the testing is 7 days as predefined days. Since we use the differential time series data, we have 60 minutes6 hours minus 1 for a day. The total 2,513 time series differential data is obtained by multiplying the 359 minutes by 7 days. Therefore, the selection of joint events m12 is close to 0.5% of the time series data in the predetermined period, which is big enough due to numerical verification in this invention. Suppose we would like a more reliable result, we can set a bigger limitation for m but it also means that some phase of the loads may not be recognized. Note that equation (3) can also be extended into:

(67) p ( u ) , and p k ( v ) , ( 6 ) Where p ( u ) = { p , k ( 1 , 1 ) + p , k ( 1 , 2 ) if u = 1 p , k ( 2 , 1 ) + p , k ( 2 , 2 ) if u = 2 ( 7 a ) p k ( v ) = { p , k ( 1 , 1 ) + p , k ( 2 , 1 ) if v = 1 p , k ( 2 , 1 ) + p , k ( 2 , 2 ) if v = 2 ( 7 b )

(68) In fact, p.sub.(u) and p.sub.k(v) in equation (6) are the marginal probability of equation (3). p.sub.(u=1) is the probability that events m at feeder head phase are greater or equal to

(69) L k p
and p.sub.(u=2) is the probability that events m at feeder head phase are less or equal to

(70) L k n
Similarly, p.sub.k(v=1) and p.sub.k(v=2) is the probability that events m at load k are greater or equal to

(71) L k p
and less or equal to

(72) L k n ,
respectively.

(73) One can follow equation (3) and (6) to build the probability table between the phase and load k. After that, an index is needed to figure out which feeder head phase in is the most related to load k when the load belongs to a wye connection. The mutual information index we use here is summarized as follows:

(74) MI ( ; k ) = .Math. u { 1.2 } .Math. v { 1.2 } p , k ( u , v ) log p , k ( u , v ) p ( u ) p k ( v ) ( 8 ) Where k , e = arg max MI ( ; k ) ( 9 )

(75) The .sub.ke is the estimated phase of load k in wye connection. It should be noted that the present invention can also be used to recognize the delta connected load by:

(76) ke = arg min MI ( ; k ) and ( 10 ) ke = - ke ( 11 )

(77) As shown above that regarding wye and delta connections, the highest correlated phase is considered as the estimation result for wye and the lowest correlated phase is considered as not the phase for a delta load, respectively.

Additional Details of the Present Invention

(78) It shall be pointed out three additional details with regards to the present invention.

(79) First, it is possible that two cases with the same mutual information index in equation (8), even if both cases are in totally opposite relationships. For example, suppose a probability table U is p(1,1)=0.1, p(1,2)=0.4, p(2,1)=0.4, p(2,2)=0.1, the result will be the same as another table V which is p(1,1)=0.4, p(1,2)=0.1, p(2,1)=0.1, p(2,2)=0.4. Although the resulting values of inputting U and V in equation (8) are the same, they are different in the sense that in V, and k are positively correlated, while in U they are not. Therefore, we will select results with higher value of p(1, 1) and p(2, 2) if two cases have the same results of (8). The reason is that p(1, 1) and p(2, 2) imply that phase and k show the same pattern, which also imply that phase and k can be on the same phase if the load k is in wye connection.

(80) Second, since the present invention focuses on the signal from the feeder head and customer load side, the value of the signal itself can impact the results. In other words, it may be easier to recognize the phase of customers with heavier loads as opposed to customers that consume less power. The reason is because fluctuations during a heavier load can be recognized more easily than those with light loads. To avoid the setback, the present invention identifies the phase from the heaviest load first. After recognizing the phase, we subtract the customer-side current from the corresponding feeder head phase if the customer load is wye connected. If the customer load is delta connected, then we will subtract two corresponding phases from the feeder head side. Then, we calculate the new probability table for the second largest load, and so on. Note that this is a valid assumption because the public utilities should have a power consumption plan when they supply the power to customer loads in the power distribution systems. The heaviest load is selected based a statistic detecting algorithm based on sample means of current magnitude time series of loads in which the load with largest sample mean of current magnitude is chosen.

(81) The feeder head current update process is implemented using a current disaggregation algorithm by assuming three phase angles are balanced, and 120 apart, and three phase angles for voltages with wye-connection are set as .sub.V.sub.A=0, .sub.V.sub.B=120, .sub.V.sub.C=120. Three phase angles for voltages with delta-connection are set as .sub.V.sub.AB=30, .sub.V.sub.BC=90, .sub.V.sub.CA=150

(82) The feeder head current before removing any loads are calculated as:

(83) 0 . FDR ( t ) = P FDR ( t ) - jQ FDR ( t ) V FDR ( t ) e - j V , ( 12 ) FDR ( t ) = .Math. "\[LeftBracketingBar]" . FDR ( t ) .Math. "\[RightBracketingBar]" { A , B , C } , t T

(84) Where

(85) P F D R ( t ) , Q F D R ( t ) and V F D R ( t )
are the measured real power, reactive power and voltage magnitude at feeder head phase at interval t, and all values in per unit.

(86) . F D R ( t )
is the complex current at feeder head phase at interval t, and the feeder head current magnitude used for phase identification,

(87) F D R ( t )
is calculated as the absolute value of complex current

(88) . F D R ( t )
as shown in (12).

(89) The feeder head currents are updated each time when the connected phases of any load originally missing phase label has been determined.

(90) If load k is wye-connected and its phase is determined as , the magnitude of current at feeder head on phase at interval t,

(91) F D R ( t )
can be updated using (13) by subtracting the contribution of current of load k from the latest feeder head current on phase at interval t:

(92) . FDR ( t ) . FDR ( t ) - P k LD ( t ) - jQ k LD ( t ) V LD ( t ) e - j V , ( 13 ) FDR ( t ) = .Math. "\[LeftBracketingBar]" . FDR ( t ) .Math. "\[RightBracketingBar]" { A , B , C } , t T

(93) Where all values in per unit,

(94) P k L D ( t ) , Q k L D ( t ) and .Math. "\[LeftBracketingBar]" V k L D .Math. "\[RightBracketingBar]" ( t )
are the measured real power, reactive power and voltage magnitude at load k and interval t

(95) If load k is delta-connected and its phases are determined as .sub.1 and .sub.2, the currents at feeder head on phases as .sub.1 and .sub.2 at interval t,

(96) 1 F D R ( t )
and

(97) 2 F D R ( t )
can be updated using (14) by subtracting the equivalent contributions of current of load k from the measured currents at feeder head phases .sub.1 and .sub.2:

(98) 0 . 1 FDR ( t ) . 1 FDR ( t ) - P k LD ( t ) - jQ k LD ( t ) V k LD ( t ) e - j V 1 2 , ( 14 ) . 2 FDR ( t ) . 2 FDR ( t ) - P k LD ( t ) - jQ k LD ( t ) V k LD ( t ) e - j V 1 2 , 1 FDR ( t ) = .Math. "\[LeftBracketingBar]" . 1 FDR ( t ) .Math. "\[RightBracketingBar]" , 2 FDR ( t ) = .Math. "\[LeftBracketingBar]" . 2 FDR ( t ) .Math. "\[RightBracketingBar]" , 1 2 { AB , BC , CA } , t T

(99) Where all values in per unit,

(100) P 1 F D R ( t ) , Q 1 F D R ( t ) and V 1 F D R ( t ) ,
and

(101) P 2 F D R ( t ) , Q 2 F D R ( t ) and V 2 F D R ( t )
are the measured real powers, reactive powers and voltage magnitudes at feeder head phases .sub.1 and .sub.2 at interval t, respectively.

(102) Third, the use of cross entropy in the present invention can perform better than the traditional use of cross correlation. Correlation is a linear relationship. However, the relationship between the load k and its corresponding feeder head phase could be positively related but not linearly related. Thus, it is better to use cross entropy to solve the problem as we proposed here.

Test Results of an Exemplar System

(103) The IEEE 34 bus test system as shown in FIG. 2A is used as an exemplar system to demonstrate the usage of the present invention. This system includes thirty loads, half in wye connection and half in delta connection. A set of sample load profiles with size of 100 have been used to simulate the load variations for the loads. Since the IEEE 34 test system only has thirty loads while the sample load profiles have one hundred load profiles, we shuffle the load profiles and assign random ones to the loads in the IEEE 34 test system in the same day with the power factor set at 0.95.

(104) The IEEE 34 test system did not include PV inverters. Therefore, we add PVs at each bus with the penetration rate from 0% to 100%. Note that each PV rated power is 250 kW. It takes irradiation and temperature as inputs.

(105) FIG. 3 is a schematic illustrating a typical solar irradiance variation within 24 hours of a typical day. The horizontal axis 310 represent the hour of the day, and the vertical axis 320 represents the irradiance measured in watt per squared meter. The curve 330 displays the variation of solar irradiance over the hour of the day.

(106) Considering the sizes of geographic area for a power distribution system should be limited, the same weather conditions are used for all PVs and the PVs are behind-the-meter (i.e. we cannot see the power and current generated by the PVs). All the simulations are conducted by using on Simulink.

(107) We will use advanced smart meters that can record not only real power but also reactive power and magnitude of voltage per minute to identify phases of loads. The above-described mutual-information based approach is applied to recognize the phase of a wye or a delta load under different PV penetration rates when a power distribution system is mixture with wye and delta loads. Different scenarios that may impact the accuracy of the present invention are analyzed, including the penetration rate of smart meters, change of the measurement rate, and the measurement error of the smart meter.

Different PV Penetration Rates with Different Time Period

(108) Table I shows the results of the present invention applied to Different PV Penetration Rates with Different Time Period. The PV penetration rate is 0%, 33.3%, 66.7%, and 100%, corresponding to PVs are 0, 10, 20, and 30 out of all 30 loads in the IEEE 34 bus test system.

(109) The accuracy is 100% for all the PV penetration rates when the present invention works with the time series differential data from 1 AM to 6 AM. The main reason is that PVs have less impacts to the event signature due to less irradiance as shown in FIG. 3 while the present invention can still recognize the event signatures between the feeder head and customer-side. The results of a different period, 11 AM to 16 PM, are also listed in the same table to show the impact of PVs. Although the present invention can still work with 100% accuracy when no PV installed, the accuracy of the present invention can decrease to 60% with PVs. The main reason is that the present invention depends on the event signature and the contribution of behind-the-meter PVs can bring challenges to identify the phase. One may noticed that the accuracy decrease to 60% when the PV penetration rate is 66.7% but increase to 70% when the PV penetration rate is 100%. A possible reason is that the simulation in this study assumes all the PVs in the IEEE 34 bus test system have the same irradiation and temperature input. This assumption may help the present invention separate loads profile easily when the PV penetration rate is 100%. For the rest of the analysis in different factors that have different challenges, we use the time series differential data from 1 AM to 6 AM for analysis.

(110) TABLE-US-00001 TABLE I Accuracy of the present invention on different pv penetration rate in different period PV Penetration Data Period Rate (%) 1 AM to 6 AM 11 AM to 16 PM 0 100.0 100.0 33.3 100.0 86.7 66.7 100.0 60.0 100.0 100.0 70.0

Different Smart Meters' Penetration Rates

(111) TABLE-US-00002 TABLE II Accuracy with different smart meter(SM) penetration rate SM Penetration rate (%) 96.7 67.7 33.3 3.3 Average Accuracy (%) 99.7 99.1 98.2 98.0

(112) Since in the present invention, we remove the corresponding current on the feeder head once we identify the phase on the customer load side, it may be interesting to know the accuracy if the smart meter penetration rate is not 100%. We want to find out how the present invention will perform if we cannot remove the current from the corresponding phase. Here, 100 tests are run with different penetration rates of smart meters (SM). Note that in the 100 tests, the deployment of smart meters are randomly picked. When SM penetration rate is 96.7%, it means only one load in the IEEE 34 bus test system is not using smart meter. When SM penetration rate is 3.3%, it means that only one load in the IEEE 34 bus test system has smart meter. The results are listed in Table II. In the worst case where SM penetration rate is 3.3%, the average accuracy is 98%, which shows the present invention can work even if only one customer uses the smart meter.

Results with Different Measurement Rate in Smart Meters

(113) One limitation of the present invention is its reliance on the high measurement rate of smart meters (i.e. per minute measurement or 60 readings per hour) which capture more fluctuations on both the customer load and the feeder head. It is estimated that the present invention will yield less satisfying results as the measurement rate decrease. Here we use the same IEEE 34 bus system with different measurement rates to show the impact. Two lower measurement rates are selected, which are 20 and 12 readings per hour, corresponding to readings of every three and five minutes. Note that these are accumulated readings and therefore, the fluctuations between the current reading and the previous reading will be smaller and less frequent in every three and five minutes as opposed to in every minute. The results are shown in Table III.

(114) TABLE-US-00003 TABLE III accuracy with different smart meter measurement rates Measurement Rate Per Hour 60 20 12 Accuracy (%) 100.0 90.0 76.7

(115) FIG. 4A is a schematic illustrating meter readings for a bus 822 of the IEEE 34 bus test system with measurement rates as 60 readings per hour, according to embodiments of the present disclosure. The horizontal axis 410 represent the hour of the day, and the vertical axis 420 represents the current magnitude of bus 822. The circles 430 represent the measurements collected at 60 readings per hour.

(116) FIG. 4B is a schematic illustrating meter readings for a bus 822 of the IEEE 34 bus test system with measurement rates as 12 readings per hour, according to embodiments of the present disclosure. The horizontal axis 415 represent the hour of the day, and the vertical axis 425 represents the current magnitude of bus 822. The circles 435 represent the measurements collected at 12 readings per hour.

(117) The accuracy decreases from 100% to 90.0% and 76.7%. An example of bus 822 in FIGS. 4A and 4B convince that

(118) 8 2 2 L D
will be less events beyond

(119) [ L k p , L k n ]
when the measurement rate is lower. When the event m is less, the present invention can be less accurate. It is due to the fact that the fluctuations are not significant enough for the present invention to identify the phase.

Measurement Error at Smart Meters

(120) Measurement error is possible for all kind of sensors. It is also important to know the impacts of measurement error on the present invention. Here we assume the measurement error in the smart meter is a normal distribution with zero mean and the standard deviations are 0.01, 0.05, and 0.1. With 100 Monte Carlo tests, the average accuracy of the present invention is summarized as Table IV below. The average accuracy will drop to 97.2%, 90.4%, and 68.9%, corresponding to different standard deviation of error measurement.

(121) TABLE-US-00004 TABLE IV Accuracy with different smart meter measurement error Measurement Error (S.D.) 0.01 0.05 0.1 Average Accuracy (%) 97.2 90.4 68.9 Lowest Accuracy (%) 96.7 80.0 60.0

Comparison with Cross Correlation and kNN

(122) Next, we compare the present invention with correlation and kNN (k nearest neighbor, note that we add a bar to avoid confusion to the customer load k). Since the load can be either a wye or delta connection, we should explain how we use correlation and kNN to solve the problem.

(123) For correlation, we use the same process, including using the same differential time series data and removing the corresponding current after the phase is recognized. Regarding wye and delta connections, the highest correlated phase is considered as the estimation result for wye and the lowest correlated phase is considered as not the phase for a delta load, respectively.

(124) For kNN method, we build the data set using the same differential time series data. By considering the combination of different feeder head phases and load k, we label the correct combination as 1 and incorrect one as 0. The dimension of the data set for the kNN is R.sup.502790. The dimension of 5,027 is twice of the time series differential data we used in the present invention (i.e. 2,513) add 1 due to the label. The dimension of 90 comes from the combination of 30 loads3 phases. The data set is divided into 80% in a training set and 20% in a testing set.

(125) FIG. 5 is a schematic illustrating a comparative results obtaining by using a kNN algorithm with different number of neighbors, according to embodiments of the present disclosure. The horizontal axis 510 represent the number of neighbors, and the vertical axis 520 represents the accuracy of phase identification. The curve 530 represents the accuracy variation with the number of neighbors.

(126) A preliminary search of k=3 shows in FIG. 5. The data set is shuffled every time before run kNN.

(127) The results for all three methods are shown in Table V based on the averaging of the results across one hundred cases. As we can see, the present invention provides a significantly higher accuracy. Note that, in the kNN process, the current is not removed from the feeder head once we obtain the phase of load.

(128) TABLE-US-00005 TABLE V Accuracy of different methods Method Proposed Method Correlation kNN Average Accuracy (%) 100.0 83.3 59.7

(129) The above-described embodiments of the present invention can be implemented in any of numerous ways. For example, the embodiments may be implemented using hardware, software or a combination thereof. When implemented in software, the software code can be executed on any suitable processor or collection of processors, whether provided in a single computer or distributed among multiple computers. Such processors may be implemented as integrated circuits, with one or more processors in an integrated circuit component. Though, a processor may be implemented using circuitry in any suitable format.

(130) Also, the embodiments of the invention may be embodied as a method, of which an example has been provided. The acts performed as part of the method may be ordered in any suitable way. Accordingly, embodiments may be constructed in which acts are performed in an order different than illustrated, which may include performing some acts simultaneously, even though shown as sequential acts in illustrative embodiments.

(131) Use of ordinal terms such as first, second, in the claims to modify a claim element does not by itself connote any priority, precedence, or order of one claim element over another or the temporal order in which acts of a method are performed, but are used merely as labels to distinguish one claim element having a certain name from another element having a same name (but for use of the ordinal term) to distinguish the claim elements.

(132) Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.