Imposing fair penalty to the harmonic sources based on the measurement data

Improving network power quality through harmonic reduction requires recognition of Harmonic Sources (HSs) to drive them to compensate their harmonics. This paper pro-poses a new method for equitable distribution of the Harmonic Compensation Cost of the network among the HSs based on the Harmonic Contribution Matrix. Each element of the Harmonic Contribution Matrix is the harmonic contribution of a speciﬁc source to the harmonic voltage of a speciﬁc bus. The output of the proposed method is a penalty curve for each HS over time. The amount of the ﬁne estimated for each individual HS per hour is a function of not only the contribution of that HS to the harmonic voltage of different buses, but also the contribution of the Harmonic Compensation Cost from the perspective of each bus, nominal voltage of each bus, and the sensitivity of each bus to the harmonic voltage. The proposed algorithm is evaluated on the IEEE 14-bus network and Esfahan regional electrical power network in Iran. The simulation results demonstrate the capability of the proposed method to allocate the hourly penalty curve to the HSs.

the case if each of these harmonics were present alone based on "second summation law" in the standard IEC 61000-3-6 [7].
In a restructured power grid with proper information and data, monitoring and control platforms, the grid voltage quality can be improved by the implementation of an economical mitigation solution such as the appropriate algorithm for the pricing of the harmonic contamination [8,9]. Due to the widespread use of non-linear loads, the use of real-time methods in data analysis has become very important for analysing the harmonic voltage data to real-time monitor and control the system [10,11]. For many countries, subscribers with low power factors are charged [12], so non-linear loads can receive a penalty determined according to their contribution to the harmonic voltage [13,14] or other indices such as harmonic adjusted power factor or service quality index [15]. The harmonic sources have to pay the penalty until the appropriate compensator is installed on their internal networks.
The voltage quality is concerned with deviations of the voltage from the ideal situation. The ideal voltage is a single-and constant-frequency sine wave and constant magnitude. Maintaining the voltage quality in the standard limitations is a network operator's responsibility. In the liberalized power network, the network operator spends the resources obtained from penalties on the HSs to purchase and install a suitable compensator or to pay to buy and install various Active Power Filter (APF) resources participating in the harmonic mitigation [1,16,17]. By applying these methods, the network managers would be able to improve the voltage quality and minimize the resulting damages [2,8,9,18]. On the other hand, when a network operator imposes a penalty on the harmonic sources, the subscriber compares the penalties with the cost of filter installation to decide among three options including installation of a compensator, paying the penalty, or a combination of them. Thus, determining the fair amount of the fine is considered as an important step. Up to now, various methods have been proposed, to determine the amount of the penalty. Some of these methods are based on the estimation of the cost of harmonic damages to penalize harmonic sources [1,14,[18][19][20]. Another set of methods seeks pricing through the determination of the HCC [8,21], which have better computability and are more reasonable for the HSs compared to the first group methods. The second advantage is that the specified penalty is comparable to the cost of harmonic power filter installation. The third group uses a combination of the first and second methods [22]. However, it seems unfair to use the combined method. Based on the combined method, the harmonic sources have to pay costs associated with both harmonic compensation and damages caused by voltage harmonics, including increasing the energy loss, decreasing the instrument lifetime, and relay malfunction, while with a full compensation, damages are minimized and it is unfair to consider them for harmonic pricing. Consequently, the determined surcharge imposed on the HSs, using the third method will be high because this method is associated with the possible errors of the first method and the HCC, which is the basis of the second method. These drawbacks make the third method inappropriate.
Previous models for harmonic penalty calculations presented in literature, have directly used the current measured at the PCC of customer to determine the penalty [18,19]. However, for example when a linear load is fed from the polluted bus, the current measured at PCC of that load also contains harmonic currents which in turn, causes an important error in calculating the penalty of that load. Hence, in order to measure the effect of any customer on the network and penalising them fairly, we have to consider the harmonic contribution of measured customer (harmonic sources) to the harmonic voltage. In our previous study [8], we proposed an approach, based on the results of harmonic contribution calculation at the PCC for pricing the harmonic pollution. This method includes some important advantages, for example, containing proper combination with an appropriate contribution determination method, and providing more fair compensation fee from the penalties paid by the HSs compared to the other previous references. In addition, this method acts according to the estimation of the harmonic contribution of each source to the harmonic voltage of Point of Com- Some subscribers may be more sensitive and pay more for higher voltage quality. If the quality of voltage provided to these subscribers fails to meet the required level, the network administrator will have to pay them compensetion [15,23]. The power network administrator should be able to serve the revenue received from harmonic sources, not only for harmonic compensation, but also for paying the penalties to the sensitive subscribers; until the appropriate compensation is made. Thus, when a HS is near to the sensitive load, this source must pay more penalty to compensate costs. In the proposed method a coefficient is introduced as the "Importance Factor (IF)" over time for different buses.
From the power network manager's point of view, the importance of the harmonic compensation varies across the different network buses. For example, the buses that feed the harmonic-sensitive loads should be given priority in compensation process. The different harmonic sources have a different contribution to the harmonic voltage of a certain bus over time. As a result, if the appropriate method is used to calculate the continuous harmonic contribution matrix of the network, it can be used to make a better penalty scheme and prioritize compensation at different network buses [11].
To address these challenges, we designed a modular method combining HCM, Norton equivalent, and harmonic compensation cost calculation as well as fairly distributing HCC to the harmonic sources, for the first time ( Figure 1). This study focuses on proposing a new algorithm for the fourth module. Previously, a method for harmonic contribution determination was successfully developed by the authors [11]. Based on the results, obtained from the contribution of each source to the voltage harmonic of each network bus at every time interval (T), the HCM was formed, which is used as the first module in our new four-module method.
In general, to the best of our knowledge, the main contributions of this paper can be summarized as follows: • To the best of our knowledge, it is for the first time that the multi-point harmonic contribution calculation method is used to determine the penalty of harmonic sources. This method made us able to consider the effects of each harmonic source to the harmonic voltage of various buses of the power network in addition to the PCC of the harmonic sources. • To the best of our knowledge, it is for the first time that a method is proposed to calculate the hourly penalty for each HS which makes subscribers able to optimally schedule their harmonic loads. • The proposed method requires only harmonic voltage and current magnitudes collected continuously by existing PQ monitors (it does not require the synchronized phase angles), which makes this method more practical and cost effective.
The proposed method is modular where any other methods can be used in modules of HCM, Norton equivalent and harmonic compensation cost calculation with no need to change the fourth module for fairly imposing the penalty on the harmonic sources.
Accordingly, the flowchart describes the proposed method in Section 2. Secondly, various steps of the algorithm are presented in detail, in Section 3. In section 4, the proposed method is implemented on the IEEE 14-bus network. The results are discussed in Section 5. In Section 6, an experimental study on Isfahan regional electric power network in Iran is presented.

General description of the proposed method
As depicted in Figure 2, the new four-module method is run in five steps where module 4 is implemented in steps 4 and 5.
Step 1: Calculation of the HCM In this paper, HCM is calculated based on the method presented in [11]. In this method after measuring the magnitude of the harmonics of the voltage of the target buses-as well as the current and harmonic voltage magnitudes of the harmonic sources, the HCM is determined. This stage provides the basics for the sharing of the HCC among the harmonic sources.
Step 2: Calculation of the Norton equivalent circuit At this point, according to the measurement data from the PCC of loads and sources, the Norton equivalent circuit is formed from the perspective of each target bus based on the method presented and verified in [11]. The buses feeding loads or connected to sources are considered as target buses.
Step 3: The HCC estimation from the perspective of each target bus At this stage, the cost of the design, construction, operation and maintenance of a suitable power harmonic filter to compensate the harmonics of the current of each individual target bus is calculated. This cost is considered as the HCC from the perspective of each individual target bus. Using a system of penalization, the HCC from the perspective of each individual target bus can be provided within Yn years. Thus, assuming that the load curve is approximately constant over this specified time period, the amount of the penalty required from the perspective of each bus for a time interval of t (min) is predicted.
Step 4: Calculation of the penalty correction factor The total harmonic voltage in the grid comes from the several harmonic sources. If the amount of compensation cost in the PCC of each harmonic source is calculated and multiplied by the harmonic contribution of that source to the harmonic voltage of that bus, the actual compensation cost will be estimated by summing up these values. In other words, if all the harmonic sources at the PCC are fully compensated, the harmonic contamination in the network will be totally eliminated. Using step 3, the HCC from the perspective of each individual target bus is calculated and its sum will be greater than the actual cost of compensation for all the harmonic sources. As a result, a penalty correction factor can be defined by dividing the sum of compensation cost fraction from the perspective of each individual target bus by the actual cost of compensation for all the HSs to calculate the real amount of penalty from the perspective of each individual target bus.
Step 5: Imposing the penalty to the harmonic sources The real compensation cost from the perspective of each target bus, obtained in step 4, is distributed among the HSs based on the contribution of each source to the harmonic voltage of each target bus. Since the importance of harmonic contamination differs from the perspective of different buses, a coefficient is multiplied by the real cost of compensation to consider this difference, resulting in the determination of the amount of the penalty that should be imposed on each harmonic source.

DESCRIPTION OF THE DETAILS OF THE PROPOSED ALGORITHM
Original as mentioned above, the proposed algorithm consists of five steps. The details of each step are described below.

Calculation of the HCM
The first step is to determine the contribution of the different HSs to the harmonic voltage of the targeted buses in the grid. The HCM calculation should be done in such a way that the contribution curve of each harmonic source to the harmonic voltage of each individual target bus can be plotted over time. The HCM is calculated by the method presented in our previous work [11]. If the number of harmonic sources in the considered network is supposed to be equal to NS and the number of buses studied would be NB, the size of HCM will be NS × NB for each time interval. Thus, the sum of the elements in each column of this matrix is equal to the sum of the contribution of different HSs to the harmonic voltage of a bus that equals to one.

Calculation of Norton's equivalent circuit from the perspective of each target bus
The harmonic voltage magnitude generated by a low-current harmonic source connected to a network with a large equivalent impedance and a large harmonic source connected to a network with a small impedance can be equivalent to each other. Compensation for these two similar amounts of voltage harmonics requires completely different filters. Therefore, it is necessary to calculate the harmonic current corresponding to the harmonic voltage measured in the target buses.
The characteristic impedance of each network bus can be calculated for different frequency values by active and passive methods [24]. Non-invasive methods are used to harmonic impedance calculation using the harmonic voltage and current data measured at a PCC without harmful impact on power system [25][26][27]. Aiming at this goal, the equivalent impedance is calculated from the point of view of each target bus using the method presented in [11] and converting Thévenin equivalent to the Norton and finally using Equation (1). Knowing the value of the magnitude of the harmonics of the voltage in the target bus, the equivalent harmonic current can be calculated as shown in Figure 3.
This harmonic of current can generate a voltage as high as V PCC in the target bus in the absence of any other harmonic sources in the grid. This current is calculated using Equation (2). (All values are in the h th harmonic order.) Also I eq is used as an input for HCC calculation in next step.
where Z u is the Norton equivalent impedance of the utility, Z c is the Norton equivalent impedance of the customer and Z eq and I eq are the Norton equivalent impedance and current of the PCC, respectively. The limitation of the proposed method is that for calculating the Norton equivalent, the target bus should be a PCC of a load or several loads. Otherwise, it is necessary to measure the harmonic current of all connections of the bus.

HCC estimation from the perspective of each bus
Compensation of the harmonics in the power grid is performed via various methods. Passive filters, especially at the high voltage levels, are one of the most commonly used harmonic power filters [28]. The inductor and capacitor values of these filters are adjusted to compensate a definite order of harmonics using the resonance phenomenon. Harmonic filters are designed to take into account a variety of conditions such as reactive power compensation at the main frequency, maintaining voltage within the permissible range at the main frequency, minimizing costs with regard to the operating capabilities etc. [29]. The capacitors and inductors required in passive filters are the major determinants of the cost of operating the filters [29,30]. At the high voltage levels, the capacitor cost is often the main determinant of the total cost of the passive filter. Capacitor price is influenced by two main factors, that is, capacitor rated voltage and capacitance. In passive filters with inductor in series with the capacitor, the capacitor has to be capable of withstanding a higher voltage than the nominal voltage of the installation bus.
Various methods are presented in literature for designing optimum filter based on the minimizing cost. In [29], cost function is calculated based on harmonic current and voltage magnitudes. In [31], the cost function was calculated based on harmonic current magnitude and unit costs of the capacitor and the reactor. Equation (3) is proposed to estimate the cost of a passive filter based on the harmonic current magnitude and the nominal voltage of the target bus. W1 to W5 can be calculated based on the regression method. In regression analysis, curve fitting is the process of specifying the model that provides the best fit to the specific curves in a certain dataset. The price list of the filters provided by manufacturers of the harmonic filters is the dataset used in regression analysis. The price of the filter in Equation (3) is determined by the two main factors, including the nominal voltage of the filter at connection point and the current passing through the filter. The goal of regression is minimizing root mean square error (RMSE) where CostF i is the cost of a suitable filter from the perspective of i th bus to compensate the harmonic current I hi at the voltage V ni . As mentioned, the method of this module (HCC calculation depicted in Figure 1) can be replaced by another method to estimate the price of the filter, which can estimate the active, passive or hybrid filter costs and does not interfere with the proposed algorithm structure.
The equivalent harmonic current, which should be compensated by the harmonic filter from the perspective of each target bus, can be calculated by determining the Norton equivalent impedance from the perspective of the same bus at any given time, based on its harmonic voltage, using Equation (2). Substituting the result of Equation (2) into Equation (3) independently, gives the HCC from the perspective of each target bus.

FIGURE 4 4-Bus example network
It is worth mentioning that the Norton equivalent impedance and the harmonic voltage of each bus are changing at any moment. However, these values are assumed to be constant over Y n years. During this period, the capital required for purchasing and installing a suitable harmonic filter system (CostF i ) has to be covered by the penalties. Then, the fraction of the penalty required from the perspective of each individual bus for a time interval of M (minutes), that is,P it , is calculated using Equation (4).
where H y is the number of hours in a year.

Calculation of the penalty correction factor
The HCC is independently calculated from the perspective of each target bus in step 3 (CostF i ). It is assumed that P it calculated within M = 60Y n H y minutes as an initial estimation of the required penalty from the perspective of i th target bus can afford the required capital for investing in application of a suitable harmonic filter system from the perspective of i th target bus. As summing up the required penalties from the perspective of all target buses results in a greater value than the actual cost of the harmonic compensation, a correction factor should be defined.
To simplify the correction factor calculation, a 4-bus network (Figure 4) is considered as an example. In this network the X and Y loads are non-linear; however, other loads are linear. In step one, the contribution of these two sources (X and Y) to the harmonic voltages of these four buses are calculated and shown as the blue numbers for each bus. In the second step, the Norton equivalent impedance is calculated from the perspective of each bus. Then, the HCC for a desired time interval from the perspective of each individual bus is calculated using Equations (3) and (4), which are shown as A-D $ for buses 1-4, respectively. If the compensation is carried out for the harmonic sources X and Y at buses 1 and 4, the harmonics in the network are almost completely compensated. Accordingly, the correction factor is calculated as Equation (5). By multiplying the correction factor by the initial estimation of the HCC from the perspective of each bus, the real amount of the penalty from the perspective of each individual bus is determined. In step 5, the real estimated penalties calculated in step 4 are distributed among the harmonic sources in accordance with the HCM calculated in step one to determine the final amount of penalty, which should be imposed on each HS. Thus, the amount of penalty imposed to the source X will be as Equation (6).

CF =
Sum of the HCCs at PCC of HSs considering HCM Sum of the HCCs at all target buses Using this method, the higher penalty is imposed on source Y, because this source has a greater role in generating harmonic voltage in the grid. The higher contribution of Y in buses 2 and 3 may be due to the lower distance of these buses from bus 4 rather than bus 1 or the higher harmonic current of source Y. The correction factor of time interval T (CFT) is generally calculated as: where N bs is the number of buses to which the harmonic sources are connected, NB is the number of target buses and HCM bs,bs is the harmonic contribution of source bs to the harmonic voltage of its PCC.

Penalty determination for each HS
It is notable that the harmonic sources in the grid are responsible for the harmonic voltage of different buses; however, according to step 1, contribution of a certain HS to the harmonic voltage of different buses is different and varies over time. Hence considering HCM and initial estimation of required penalties from the perspective of each target bus and the correction factor obtained from steps 1-4, the amount of the penalty that should be imposed on each HS for a certain time interval, can be calculated by Equation (8).
Where the HCM is the harmonic contribution matrix of the harmonic voltage of the different buses in the time interval T and PS is the amount of penalty that should be imposed on each HS during the same time interval.
It is mentionable that the amount of penalty imposed on each HS is not just influenced by the HCC from the perspective of each target bus, but by the sensitivity of each bus to the harmonic voltage. Consequently, an Importance Factor (IF) for each bus should be defined, expressing its sensitivity to the harmonic voltage. If the IF factor is optimally determined by the network manager, the compensation costs and cost of subscribers affected by the off-contract or standard harmonic voltage will be equal to the total fines imposed on the harmonic sources. As a consequence, the final amount of the fine imposed on each bus (CPS) is calculated using Equation (9).
The corrected amount of the penalty imposed on each bus, can be calculated once every T m minutes over one hour. As a result, the average of 60/T m values of the final amount of the penalty calculated for each HS represents the amount of the penalty per hour-which should be imposed on each HS. This process is repeated for each hour over a billing cycle, for example, 1 month. So, the penalty, which should be imposed on each HS in a billing cycle (APs); is calculated using Equation (10).
where • refers to the Hadamard product and T is the number of interval, NT is the number of T m minutes interval in an hour, NB is the number of target buses, P b,T is the HCC from target bus b point of view, HC M s,b,T is the harmonic contribution of source s to the harmonic voltage of bus b at the T th interval. Moreover, I F b,T is the importance factor of bus b at the T interval. In this equation C F T is the correction factor of T th interval. In addition, a billing cycle of 30 days (720 h) has here been eventually considered.

CASE STUDY
The proposed algorithm is, moreover, implemented on the IEEE 14-bus power network to evaluate the capability of the method in order to determine the appropriate penalties,

Classification of consumers
Consumers in this sample network are classified into three categories, including urban, industrial type-one and industrial type-two. The loads, connected to buses 2 and 9, are industrial type-1 and load-connected to the bus 4 is industrial type-2. The industrial consumers without night shift are considered as industrial type-1, while consumers of industrial type-2 have a more uniform load curve and significant energy consumption at night. Other network loads fall into the urban category.

Load curves
For each type of consumers, described in Section 4.1, a normalized load curve was assumed to be multiplied by the original load value of each bus in 14-bus network to calculate their hourly load values. Table 1 presents these normalized load curves in detail.

Non-linear loads
Loads connected to buses 2, 4 and 9 are considered as nonlinear loads. The harmonic current of these sources over time is assumed to be a constant percentage relative to the main component current as shown in Table 2. Varying percentages can also be used in the proposed model.

Sensitivity of each bus to the harmonic voltage (importance factor)
The coefficient of importance for each bus expressing its sensitivity to the harmonic voltage varies not only over time, but also it depends on the type of the loads that each bus feeds. Table 1 shows the coefficients for the 14-bus network buses presented in Figure 5 over time.

Changes in the network connections
It is anticipated that there are two changes in network connections every 24 h. The first change is the outage of the transformer connected between buses 4 and 9 (first outage). It starts at 10:00 and ends at 15:00. The second change is the outage of the line between buses 2 and 4 (second outage) that starts at 20:00 and ends at 24:00. These outages are shown by arrows in Figure 5.

RESULTS AND DISCUSSIONS
The proposed algorithm was implemented on the IEEE 14-bus network in DIgSILENT software with the settings described in Section 4 for the 5th and 7th harmonic order, and the results are discussed step by step. The whole simulation time by CPU i7 2.7GHZ and 12 GB RAM is 8.3 s. The softwares used are Python and DPL (DIgSILENT Programming Language).

Calculation of the HCM
Assuming that measurements are taken every 6 min, there will be 10 records per hour and 240 records per day. Using every The harmonic contribution curve of the load 9 to the 5th harmonic voltage of different buses The harmonic contribution curve of the load 4 to the 7th harmonic voltage of different buses three consecutive data sets, the contribution of each HS to the harmonic voltage of different buses is calculated [16]. Figure 6 shows the contribution curve of load 9 to the harmonic voltage of various buses at the 5th harmonic order. Figure 7 shows the contribution of load 4 to the harmonic voltage of the various buses at the 7th harmonic order. For a better observation, only several buses are displayed in each figure.
As shown in Figure 6, when the first outage takes place (records no. 100-150), the contribution of load 9 to the harmonic voltage of buses at 132 kV (i.e. buses 2 and 4) at 5th harmonic order decreases, while its contribution to 33 kV buses (i.e. buses 6, 9, and 13) increases due to the increased equivalent impedance. A reverse effect can be seen in Figure 7. As shown in Figure 7, during the first outage, the contribution of load 4 to the harmonic voltage of buses at 132 kV (i.e. buses 2 and 4) at 7th harmonic order increases, while its contribution to 33 kV buses (i.e. buses 6, 9, and 13) decreases. Hence, the HCM can be updated in response to any change occurring in the operating condition of the network.

Calculation of Norton's equivalent circuit from the perspective of each bus
In this step, the Norton equivalent impedance is determined from the perspective of each bus based on the measurement data. Figure 8 shows a view of the Norton equivalent impedance curve from the perspective of buses 5 and 14 at the 5th and 7th harmonic orders. As can be seen, the impedance magnitude from the perspective of bus 5 is greater than that of the perspective of bus 14. On the other hand, during the first outage (records no. 100-150), there is a drastic increase in impedance magnitude from the perspective of bus 5 at the 7th harmonic order. Then, the equivalent harmonic current that should be compensated by the filter from the perspective of each target bus is calculated using Equation (2).

HCC estimation from the perspective of each bus
After determining the equivalent harmonic current that has to be compensated by the filter from the perspective of each individual target bus, the HCC from the perspective of that target bus is estimated using Equation (3). The coefficients used in Equation (3) are presented in Table 3. These weights are calculated based on the price list of harmonic filters in Iran.
The fraction of the HCC for a time interval of 6 min, which is an initial estimation of the penalty from the perspective of each target bus, is calculated using Equation (4). This calculation is repeated every 6 min over 24 h for all buses at the 5th and 7th harmonic order; the results are shown in Figures 9 and 10, respectively.

FIGURE 9
The HCC for each time interval between records at 5th harmonic order

FIGURE 10
The HCC for each time interval between records at 7th harmonic order As shown in Figure 9, the highest and lowest costs at 5th harmonic order are related to bus 4 and bus 12, respectively. The lowest cost associated with bus 12 is due to its high impedance magnitude and the low level of the harmonic voltage of the 5th order in this bus. As can be seen from Figure 10, for the 7th harmonic order, the highest cost is obtained from the perspective of bus 4. Because of limitation explained in Section 3.2, buses #1, #7 and #8 are excluded from the analysis and results.
The load curve data summarized in Table 1 show that the maximum amount of industrial loads often takes place from

FIGURE 11
The magnitude of harmonic voltage at 7th harmonic order 10 a.m. to 3 p.m. Hence, the harmonic value of voltage in different buses approaches its peak value during these hours. If the first outage did not occur, the highest harmonic voltage for all buses would occur during the mentioned hours. As shown in Figure 11, due to the first outage, the effect of loads 2 and 4 (especially load 4), on the harmonic voltage of the 33 kV section of the network has been reduced. Therefore, the HCC from the perspective of 33 kV buses decrease significantly during the first outage ( Figure 10). The results presented in this section confirm the capability of the proposed algorithm to respond favourably to any change occurring in a network.

Calculation of the penalty correction factor
The initial estimated penalties from the perspective of each target bus are calculated and summed up, every 6 min for 24 h. The results are shown at the 5th and 7th harmonic orders in Figure 12. The area under the penalty curve indicates the amount of penalty per day. Assuming that all days throughout Y n years have identical penalty curves and Y n = 2 years, the total penalty per day at the 5th and 7th harmonic order are equal to $40,001 and $51,581, respectively.
This assumption is only used for calculating the hourly share of the total compensation cost. In fact, the capital required for purchasing and installing a suitable harmonic filter system should be covered by the fines received within Y n years. It is necessary to find out how much fine should be received per hour within Y n years. This assumption is repeated for each hour. In fact, the situation of the network is always changes, thus, hour by hour, different voltage harmonics and the harmonic contribution matrix are applied to the subscribers' penalties.
However, the compensation is, indeed, done in harmonic sources (i.e. buses 2, 4 and 9). The penalty from the perspective of each bus connected to a harmonic source is multiplied The CF at 5th and 7th harmonic order by the contribution of that particular HS to the harmonic voltage of that particular bus. The sum of these values represents the actual compensation cost.
As summarized in Table 4, the actual compensation costs at the 5th and 7th harmonic orders are equal to $5932 and $9281, respectively, which are smaller than the corresponding values of the total penalty per day obtained from summing up the penalties from the perspective of all buses. Table 4 also shows the initial estimation of penalties from the perspective of each individual harmonic source per day at the 5th and 7th harmonic orders.
Then, the correction factor is calculated every 6 min for 24 h according to Equation (8) to be used for calculation of the real amount of penalty from the perspective of each target bus. Figure 13 shows the curve of correction factor calculated for the 5th and 7th harmonic order.

Determination of the penalties to the HSs
The penalizing system should impose fine on the harmonic sources; hence, in this step, the HCM is used to convert real

FIGURE 14
Imposed penalty on the HSs at 5th harmonic order-without considering IF

FIGURE 15
Imposed penalty on the HSs at 7th harmonic order-without considering IF amount of penalty from the perspective of each target bus to amount of penalty that should be imposed on each harmonic source. Figures 14 and 15 show the amount of penalty that should be imposed on each HS according to Equation (8), every 60 min for 24 h at the 5th and 7th harmonic orders, respectively.
To understand the effect of penalty correction factor on determination of the real amount of penalties that should be imposed on each HS, Equation (8) is used with and without considering CFT. The results are summarized in Tables 5 and 6, respectively.
Comparing Tables 4 and 6 demonstrates how considering the contribution of each HS to harmonic voltage of each bus can  result in a more equitable distribution of the actual compensation costs among various harmonic sources. As mentioned in Sections 3-5, to calculate the final amount of penalty that should be imposed on each HS, the sensitivity of each bus to the harmonic voltage should be taken into account according to Equation (9). As a result, a harmonic source generating harmonic voltage in a more sensitive bus has to pay a higher amount of penalty compared to other HSs, which will excite that HS to reduce its harmonic pollution and consequently to decrease its harmonic bills. Applying the coefficient of importance presented in Table 1 according to Equation 9 resulted in Figures 16 and 17, which illustrate the final amount of penalty imposed on each HS every hour for 24 h at the 5th and 7th harmonic orders, respectively.
The final amounts of penalties imposed on each HS per day at the 5th and 7th harmonic orders are summarized in Table 7.
Comparing the results presented in Table 7 with those presented in Table 6 shows the effect of sensitivity of each bus to the harmonic voltage. In addition, comparing the last column of these two tables shows that considering the sensitivity leads to an increase in the total amount of penalty obtained from summing up the amounts of penalty for each HS. Consequently, considering the IF defined by the network administrator not only provides the required budget for compensation, but also pays the additional penalty to the sensitive customers who might be damaged by the harmonic voltage. Given the average energy price of $0.15 per kWh, the monthly energy bill is calculated as shown in Table 8.
As it can be seen in Table 8, the penalty factor defined as the ratio of harmonic penalty (for 5th and 7th harmonic orders) to energy bill ranged from 6.8% to 8.8%. The load 4 makes larger contribution to the harmonic voltage of the network (especially to the harmonic voltage of sensitive buses), so it has a higher penalty factor ratio than that of the other two HSs.

EXPERIMENTAL ANALYSIS
Esfahan electrical power grid with more than 700 buses and a load with approximate capacity of 3.4 GW is used for experimental analysis. The overall scheme of this network, the location of the harmonic sources and the results of harmonic contribution of two harmonic sources are shown in [11]. As it can be seen, HSs are shown by numbers 1 and 2, and the location of the target buses are marked by numbers 1 to 4. The simulation time by CPU i7 2.7GHZ and 12GB RAM is 6.5s. The average HC of sources 1 and 2 to the harmonic voltage of all target points (7th harmonic order) calculated over a 2-day period of simultaneous measurement is illustrated in Table 9.
Based on the proposed method, the HCC from the perspective of all target buses is separately calculated using Equation (3) and presented in Figure 18. Weights of this equation are calculated based on the price list of harmonic filters in Iran that is presented in Table 3.
Nominal voltage of all target buses is 63KV. The equivalent impedance from the perspective of P4 is the lowest. Because of this fact, P it is highest for this bus at the study time. The hourly HCC from each target bus point of view is highly dependent on the hourly harmonic voltage magnitude of that bus. The harmonic voltage magnitude of 7 th harmonic order is given in [11]. The average equivalent impedance at 7 th harmonic order from the perspective of each target bus is presented in Table 10.

FIGURE 18
The HCC for each hour at 7th harmonic order from the perspective of each target bus The correction factor is calculated using Equation (7) for each hour in the study time and depicted in Figure 19. Therefore, the maximum values of CF are obtained at the times when the highest harmonic voltage magnitude has occurred.
According to [11], the magnitude of harmonic voltage at 7 th harmonic order on the target buses are under the standard threshold. On the other hand, as presented in Table 9, harmonic sources 1 and 2 have individually a small HC to the harmonic voltage of target buses. Therefore, the hourly penalty imposed to each harmonic source, shown in Figure 20, is less than $4. The annual harmonic penalties for HS1 and HS2 are calculated as $20,357 and $10,026, respectively (the Importance Factor (IF) is considered 1 in this case study).

FIGURE 19
The CF at 7th harmonic order FIGURE 20 Imposed penalty on the HSs at 7th harmonic order

CONCLUSION
Fair penalization system is a way for motivating harmonic sources to keep their harmonics within the defined limits and standards; therefore, by improving the network power quality, sensitive customers can achieve their desired voltage quality. To this aim, a modular method was proposed in this paper for the estimation of the amount of penalty that should be imposed on each HS to provide the required budget for harmonic compensation based on the contribution of each HSs to the harmonic voltage of different buses. In this method, the location of the harmonic sources in the power grid, changes of the network connections over time, contribution of different sources to the harmonic voltage of different buses, nominal voltage of the buses, sensitivity of the customers fed from each bus and the amount of harmonic voltage in different buses are taken into account. The output of the proposed method is a continuous penalty curve for each harmonic source-which can be used for the estimation of the annual penalties of harmonic sources and for the harmonics compensation planning. By adjusting the IF and estimating the HCC, the network manager can provide not only the necessary capital to compensate the harmonic voltage made by the harmonic sources, but also provide additional penalty to pay to the sensitive customers who may be damaged by the harmonic voltage.