An Analysis of the Systematic Error of a Remote Method for a Wattmeter Adjustment Gain Estimation in Smart Grids

Smart energy meters supporting bidirectional data communication enable novel remote error monitoring applications. This research targets characterization of the systematic worst-case error of the previously published remote watthour meter’s gain estimation method based on the comparison of synchronous measurements by the reference and meter under test. To achieve the research aim a methodology based on global maximization of the systematic error objective function assuming the typical low voltage electrical distribution network operation parameters ranges as defined by the standard recommendations for network design. To cross verify the reliability of the assessed solutions the suggested error analysis methodology was implemented utilizing two stochastic global extremum search techniques (genetic algorithms, pattern search) and the third one utilizing nonlinear programming solver. It was determined that the wattmeter adjustment gain worst-case error does not exceed 0.5% if the remote wattmeter monitored load power factor is larger than 0.1 and a network is designed according to the recommendation of the acceptable voltage drop less than 5%. For a load exhibiting power factor larger than cos φ = 0.9 the worst-case error was found to be less than 0.1%. It is concluded therefore that considering the systematic worst-case error the previously suggested remote wattmeter adjustment gain estimation method is suitable for remote error monitoring of Class 2 and Class 1 wattmeters.


Introduction
The increasing amount of distributed electrical energy generation, smart revenue and subaccounting meters, interconnected via the advanced measurement infrastructure (bidirectional data communication channels between meters and data collection modules) are the key features of the evolving modern smart electrical grids. The demand the metrological reliability of the readings (trusted metering) of measurement devices cannot be underestimated. Not only it is crucial for maintaining the correct energy consumption billing (revenue metering) but also provides the electricity consumer with their energy consumption budget (submetering, home automation, etc.) which on its way enables changing the energy consumption behavior towards better energy efficiency of buildings, logistics, manufacturing processes, etc. The revenue metering (Level 2 in Figure 1) is a traditional A new method for remote estimation of energy meter's adjustment factor was proposed in [3]. This method is aimed to determine the multiplicative correction factor for a separate energy meter using a remote reference energy meter. This method can be applied for verification of the selected meter despite the presence of any other power consumption occurrences by the same grid loads. The meter under inspection (later in the paper called remote wattmeter RemW) and the reference meter are installed in the electrical grid in such a way that the reference meter (later in the paper called reference wattmeter RefW) is able to measure energy consumption of the part of electrical grid including the electrical load monitored by RemW ( Figure 2). The data delivery and measurement synchronization are mandatory between RefW and RemW, and a reasonably small (around 10 W) switchable verification load has to be built into the RemM or connected as an external add-on module.
During the derivation of the RemW adjustment gain, certain assumptions influencing the systematic error of the coefficient estimation are made. These assumptions are related to the parameters of the electrical grid interconnecting RemW and RefW. The systematic error of the estimation was not yet investigated for wide range of electrical distribution networks. Usually, it is of interest to explore the worst-case error (WCE) possible in the real world electrical grid for any A new method for remote estimation of energy meter's adjustment factor was proposed in [3]. This method is aimed to determine the multiplicative correction factor for a separate energy meter using a remote reference energy meter. This method can be applied for verification of the selected meter despite the presence of any other power consumption occurrences by the same grid loads. The meter under inspection (later in the paper called remote wattmeter RemW) and the reference meter are installed in the electrical grid in such a way that the reference meter (later in the paper called reference wattmeter RefW) is able to measure energy consumption of the part of electrical grid including the electrical load monitored by RemW ( Figure 2). The data delivery and measurement synchronization are mandatory between RefW and RemW, and a reasonably small (around 10 W) switchable verification load has to be built into the RemM or connected as an external add-on module. conditions of the network and connected loads always vary within some bounds. Derivation of analytical expression of the error as a function of grid and load parameters is a complex task because of various possible grid topologies and load characteristics including active and reactive components, harmonic distortion, etc. Even though, it might be possible to derive an analytical expression of the error but searching for the maximum WCE value and the corresponding set of grid and load parameters might be a challenging nonlinear constrained optimization task. The goal of the research presented hereain is to determine the systematic WCE of an energy meter's remote adjustment gain estimation procedure according to the method described in [3]. The other constituent of the total error is a random error, which was already partially explored in [3,4]. The idea of a remote watthour meter's online adjustment gain estimation method [3] is quite new and dependence of its systematic error on circuit parameters is not investigated yet. Field calibrationbased methods [5][6][7] utilize portable calibrators directly connectable to units under calibration and the calibration procedure is executed offline, i.e., with a disconnected user loads. Impedance between calibration unit and portable calibrator is small and is neglected. Specifics of the investigated remote online adjustment gain estimation method is that some irrefutable impedance separates a reference unit RefW and unit under adjustment RemW. Furthermore, if an adjustment gain evaluation procedure is carried out online, the error of the procedure is also a function of the connected loads' size and type. There are some contributions to solve similar problems presented by other researches. G. B. Samson et al. proposed an online auto-calibration technique of system for energy consumption measurement based on low price but low accuracy Hall effect sensors (HES) [2]. Sensors for monitoring all loads connected to the household grid are installed at outputs of each breaker. The sensors are calibrated according to measurements of a single high precision current transformer installed at the main power line. To improve accuracy of the entire HES network a least mean square algorithm was applied to adjust the gain of each HES sensor in the network. The authors present comparative analysis of HES sensors network error with and without calibration for particular loads connected, i.e., three common residential appliances. Results of MATLAB simulation for the case with a higher number of sensors were also presented but detailed analysis of error dependence on system load types and magnitudes were omitted. Furthermore, the procedure proposed [2] differs in nature from the procedure introduced by authors in [3] and later explored in our work presented in the publication.
Since an evaluation of all possible combinations of external loads (types, magnitudes and power factors) in our research is costly to implement using physical models, a numerical simulation was applied. Whereas a problem is to find an extremum value for the defined bounds of parameters, the problem could be treated as a multi-parameter optimization task with non-linear constrains and parameter bounds. The genetic algorithms (GA) are often used to solve this type of During the derivation of the RemW adjustment gain, certain assumptions influencing the systematic error of the coefficient estimation are made. These assumptions are related to the parameters of the electrical grid interconnecting RemW and RefW. The systematic error of the estimation was not yet investigated for wide range of electrical distribution networks. Usually, it is of interest to explore the worst-case error (WCE) possible in the real world electrical grid for any newly proposed measurement method. By specifying the WCE despite the electrical grid configuration and load parameters it is possible to identify the applicability of the method for remote verification of the energy meters. The real electrical grid is designed according to various recommendations and regulations. Therefore, the set of parameters characterizing real operating conditions of the network and connected loads always vary within some bounds. Derivation of analytical expression of the error as a function of grid and load parameters is a complex task because of various possible grid topologies and load characteristics including active and reactive components, harmonic distortion, etc. Even though, it might be possible to derive an analytical expression of the error but searching for the maximum WCE value and the corresponding set of grid and load parameters might be a challenging nonlinear constrained optimization task.
The goal of the research presented hereain is to determine the systematic WCE of an energy meter's remote adjustment gain estimation procedure according to the method described in [3]. The other constituent of the total error is a random error, which was already partially explored in [3,4]. The idea of a remote watthour meter's online adjustment gain estimation method [3] is quite new and dependence of its systematic error on circuit parameters is not investigated yet. Field calibration-based methods [5][6][7] utilize portable calibrators directly connectable to units under calibration and the calibration procedure is executed offline, i.e., with a disconnected user loads. Impedance between calibration unit and portable calibrator is small and is neglected. Specifics of the investigated remote online adjustment gain estimation method is that some irrefutable impedance separates a reference unit RefW and unit under adjustment RemW. Furthermore, if an adjustment gain evaluation procedure is carried out online, the error of the procedure is also a function of the connected loads' size and type. There are some contributions to solve similar problems presented by other researches. G. B. Samson et al. proposed an online auto-calibration technique of system for energy consumption measurement based on low price but low accuracy Hall effect sensors (HES) [2]. Sensors for monitoring all loads connected to the household grid are installed at outputs of each breaker. The sensors are calibrated according to measurements of a single high precision current transformer installed at the main power line. To improve accuracy of the entire HES network a least mean square algorithm was applied to adjust the gain of each HES sensor in the network. The authors present comparative analysis of HES sensors network error with and without calibration for particular loads connected, i.e., three common residential appliances. Results of MATLAB simulation for the case with a higher number of sensors were also presented but detailed analysis of error dependence on system load types and magnitudes were omitted. Furthermore, the procedure proposed [2] differs in nature from the procedure introduced by authors in [3] and later explored in our work presented in the publication.
Since an evaluation of all possible combinations of external loads (types, magnitudes and power factors) in our research is costly to implement using physical models, a numerical simulation was applied. Whereas a problem is to find an extremum value for the defined bounds of parameters, the problem could be treated as a multi-parameter optimization task with non-linear constrains and parameter bounds. The genetic algorithms (GA) are often used to solve this type of optimization problems, including those coming from electric power system field. The Nonintrusive Load Monitoring system (NILM) approach based on genetic algorithms [8] can be considered as an example. The NILM [9][10][11][12] method is meant to disaggregate the consumption of the entire active power measured by standard digital watt-hour meter into its major electrical loads. GAs together with the support vector machines are the main tools of the proposed framework [13], allowing to detect any abnormal consumer behavior in the power distribution sector and identify probable power grid frauds. Other groups of power systems problems to solve using the GA are the energy consumption optimization tasks [14,15] and optimizations for forecasting purposes [16,17]. The nature of the problem under consideration is very similar to the optimal power flow [18], optimal automation devices [19] or distributed generation sources [20] allocation, optimal operating and scheduling of micro grids [21] and other problems [22,23] but with different objective functions. These are the entire tasks based on network operation mode calculations with the pre-defined non-linear constraints to find global extremum, i.e., the minimal losses, the minimal or the maximal flows, the minimal investments, or the worst-case error for the case under consideration.
The main contribution by authors in this publication is the estimation of WCE of remote gain adjustment estimation method [3] by means of solving the global nonlinear optimization problem using three different solvers for cross verification of the results reliability. The objective function for the error maximization was also composed twofold and each optimization problem separately solved by three different in nature optimization techniques, namely genetic algorithm, pattern search and nonlinear optimizer FMINCON. First, the analytical expression of the objective function was derived and used to find the WCE and electrical grid (wiring losses) and load parameters (active power, power factor). Second, the objective function was implemented using a Simulink model which was called by the used optimizers in each iteration step. The second approach is more universal, because any electrical grid configuration or load models can be applied in the WCE analysis. However, it imposes much higher requirements upon computational resources compared to analytically expressed objective function as reported in the comparative analysis contributed by the authors.
This work is structured as follows: Section 2 presents the description of the researched remote energy meter adjustment gain estimation method and its systematic error analysis; Section 3 presents the proposed methodology; Section 4 presents the method's systematic error estimation results obtained by several methods described in the Section 3.

Method Description and Assumptions
A method for remote wattmeter active power adjustment gain estimation was introduced in [3].
The key assumptions made for the method equations derivation are:

‚
There is a reference meter installed at the input of low voltage distribution grid similarly as shown in Figure 1 or  The aim of the method is to determine an estimate kp of active power adjustment gain of the RemW. The definition of the power adjustment gain is: where P is the actual power and P˚is the indication of active power by the RemW. The methodology of the method development is based on some later stated assumptions followed by the derivation of expression for the adjustment gain. Afterwards, the adequacy of the made assumptions is verified by estimating systematic error. If the error is found to fall in to the acceptable range, then it is confirmed that the method is valid in the specified range of electrical network and load parameters. In this research authors focused on the analysis of the systematic WCE due to the assumptions made in the derivation of the remote adjustment gain expression.
A network model including wattmeters RefW and RemW, input energy source U 1 , consumer load Load 3, network loads Load 1 and Load 2, Wiring between RefW and RemW, and stimulus load P S are shown in Figure 3. Complex powers of each load and losses are denoted by S, active power by P, and reactive power by Q. where is the actual power and * is the indication of active power by the RemW. The methodology of the method development is based on some later stated assumptions followed by the derivation of expression for the adjustment gain. Afterwards, the adequacy of the made assumptions is verified by estimating systematic error. If the error is found to fall in to the acceptable range, then it is confirmed that the method is valid in the specified range of electrical network and load parameters. In this research authors focused on the analysis of the systematic WCE due to the assumptions made in the derivation of the remote adjustment gain expression.
A network model including wattmeters RefW and RemW, input energy source U1, consumer load Load 3, network loads Load 1 and Load 2, Wiring between RefW and RemW, and stimulus load PS are shown in Figure 3. Complex powers of each load and losses are denoted by S, active power by P, and reactive power by Q.
SW position 2: 2 = 1(2) + 2(2) + 3(2) + + (2) where 1 , 2 , 3 are powers at the point of RefW corresponding to switch SW position. The active powers 1 , 2 , 3 are the powers of the corresponding loads (Load 1, Load 2, Load 3). In Equations (2)-(3) used notation ( ) k is used to indicate load location according to the schematics in Figure 3, while m is used to indicate position of the switch SW. is a grid wiring power loss when the load is disconnected. ∆ ( ) is the increase of network wiring active power losses due to connection of the : , are the differences of corresponding powers caused by the voltage change 2 due to connection and disconnection of : Active power at the connection point of the reference wattmeter (RefW) could be decomposed to the sum of network load powers and network wiring loss power for each SW ( Figure 3) position: SW position 2 : SW position 3 : P r3 " P 1p3q`P2p3q`P3p3q`PWp3q`∆ P 2 p∆U 2 q`∆P 3 p∆U 2 q, where P r1 , P r2 , P r3 are powers at the point of RefW corresponding to switch SW position. The active powers P 1 , P 2 , P 3 are the powers of the corresponding loads (Load 1, Load 2, Load 3). In Equations (2)-(3) used notation P kpmq k is used to indicate load location according to the schematics in Figure 3, while m is used to indicate position of the switch SW. P W is a grid wiring power loss when the load P S is disconnected. ∆P W pP S q is the increase of network wiring active power losses due to connection of the P S : Energies 2019, 12, 37 6 of 26 ∆P 2 p∆U 2 q, ∆P 3 p∆U 2 q, are the differences of corresponding powers caused by the voltage change U 2 due to connection and disconnection of P S : where the voltage difference ∆U 2 is caused by voltage losses in network wiring due injection of load power P S . Powers that should be measured by the RemW for each switch SW position can be expressed: SW position 2 : P m2 " P 3p2q`PS , SW position 3 : where P m1 , P m2 , P m3 are the powers at the connection point of RemW for corresponding switch SW position ( Figure 3). Having unadjusted indications of the RemW wattmeter Pm 1 , Pm 2 , Pm 3 and considering Equation (1) the following expressions can be written: Substituting Equations (11)- (13) to Equations (2)-(4) it is possible to rewrite: SW position 2 : P r2 " P 1p2q`P2p2q`kp¨Pm2`PWp2q`∆ P Wp2q pP S q, SW position 3 : P r3 " P 1p3q`P2p3q`kp¨Pm3`PWp3q`∆ P 2 p∆U 2 q.
The main assumption of the method described in [3] is that P 1 , P 2 , P 3 and P W are constant during the remote adjustment gain estimation procedure. Therefore, these assumptions are described as: Based on the experimental data collected from typical buildings (office, multi-apartment block, private house), it was shown in [4] that this assumption is valid with a high probability, for instance, it is 90% probability that active power consumption is fluctuations are within 1% from its nominal value during the period of 3 s.
The second assumption of the method [3] is expresed by the approximation: The assumption (18) is equivalent to the following approximations: The assumption (19) introduces a systematic error of the adjustement gain estimate kp according to the proposed method. Its impact on the systematic error is analyzed in detail aiming to verify if the assumption (19) forced errors to exceed the target level. To say in advance, it is shown at the end of the paper that even with the assumption (19) the adjustment gain estimation systematic errors stay within the level suitable to remotely monitor errors of Class 1 and Class 2 meters.
According to schematics in Figure 3 and neglecting the influence of RemW internal power losses the following equality is valid: By subtracting both equality sides of (15) from both sides of equality (14) and taking into consideration (17), (18) and (20) it can be derived: Power indications of RemW Pm 1 and Pm 2 are related to actual powers P m1 and P m2 by gain coefficient k p : Because the method of remote adjustment gain estimation relies on the accuracy of RefW, its power indications Pr 1 and Pr 3 are assumed to be equal to the real (measured) power values, i.e.,: To calculate kp according to Equation (21) only unadjusted wattmeter indications are available. Therefore, from substituting (22)- (25) to (21) it can be rewritten: The relative error of the wattmeter adjustment gain is defined by: Substituting Equation (26) to (27) it is possible to reduce the latter to:

Analytical Expression of Relative Error
In this section we seek to derive an analytical expression relating δk p to network wiring and wattmeter load parameters. In [3] it was shown that losses in the internal resistance of a wattmeter RemW due to the current sensing shunt resistance (a typical value is 0.5 mΩ) are negligible because the equivalent resistance of Load 3 is many times larger. For example, 10 kW active power load in 230 V AC grid corresponds to around 5 Ω and increases with lower active powers. Regarding the fact that the internal impedance of a wattmeter could be neglected, the loads Load 2 and Load 3 could be considered as connected to the same potential point of the circuit. Therefore, both loads could be replaced by equivalent load: Energies 2019, 12, 37 8 of 26 The replacement (29) does not influence operation mode of the remaining part of the circuit, i.e., it does not influence I W , P W , U 2 , ∆P W pP S q, ∆U 2 , P s , and ∆P W pP S q. Let us denote the sum of Load 2 and Load 3 active powers: P l " P 2`P3 , (30) and the correspondingly the change of it due to the change of voltage ∆U 2 by: The results of numerator and the denominator of (28) are indifferent to the connection point of the Load 2. Indeed, considering (2), (4), (8), (9), (17), (30), and (31) when the Load 2 is connected in front of the wattmeter: and when the Load 2 is connected behind of the wattmeter: Since the replacement of loads Load 2 and Load 3 with an equivalent load S l " S 2`S3 does not affect the error expression δk p (28) and Load 1 does not influence any of the component in equations (32)-(35), the equivalent circuit shown in Figure 4 is later used for the derivation of analytical expression of δk p . The load S l is modeled by a serial connection of resistive, inductive and capacitive components, while wiring impedance is modeled by a serial connection of resistance and inductance.
Energies 2018, 11, x FOR PEER REVIEW 8 of 25 Since the replacement of loads Load 2 and Load 3 with an equivalent load = 2 + 3 does not affect the error expression (28) and Load 1 does not influence any of the component in equations (32)-(35), the equivalent circuit shown in Figure 4 is later used for the derivation of analytical expression of . The load is modeled by a serial connection of resistive, inductive and capacitive components, while wiring impedance is modeled by a serial connection of resistance and inductance. Equations (2)- (4) and (8)-(10) for the circuit shown in Figure 4 could be rewritten as: SW position 2: 2 = (2) + + (2) (17) to (28) the analytical expression for the relative error can be derived: .
(42) Equations (2)- (4) and (8)-(10) for the circuit shown in Figure 4 could be rewritten as: SW position 2 : P r2 " P lp2q`PS`PWp2q`∆ P W pP S q, SW position 3 : P r3 " P lp3q`PWp3q`∆ P l p∆U 2 q, SW position 1 : SW position 2 : P m2 " P lp2q`PS , SW position 3 : P m3 " P lp3q`∆ P l p∆U 2 q.  (17) to (28) the analytical expression for the relative error δk p can be derived: By denoting: Equation (42) can be rewritten: Denoting the wiring loss in the case of SW position 1 or 2 (P S is connected) by P W pP S q, the following expression is written: Then Equation (43) can be rewritten as: The admittance of the overall circuit for the switch SW in position 3, could be expressed as: The admittance of the overall circuit for the switch SW in position 1 or 2 could be expressed as: The equivalent conductance of load and switching resistor R S is: The equivalent susceptances of load and switching resistor R S is: Recalculating G Sl and B Sl to series connected resistance and reactance yields: The current through a network wiring in the case of SW position 3 according to Ohms law is: and in the case of SW position 1 or 2 it is: The difference of power losses in wiring with a connected and disconnected R S is: The second error component in (45) caused by a difference ∆P l∆U 2 " ∆P ∆U2`∆ P ∆U3 with a connected and disconnected R S could be expressed by: Voltage at the input of RemW when SW is at the position 3 (R S is disconnected) could be expressed using the vector diagram method: Voltage U 2S on the consumer side when R S is connected could be expressed as: where ϕ is the phase angle between voltage U 1 and current I W : and ϕ S is a phase angle between voltage U 1 and current I Ws : The power of load without connected P S and disconnected P S are correspondingly: Then the difference of consumer power ∆P l p∆U 2 q with a connected and disconnected P S can be expressed as: The power difference ∆P l p∆U 2 q expressed as a function of circuit parameters and voltage U 1 : Power P S consumed by R S in the denominator of the (42)-(44) is: Substitution of (56), (65) and (66) to the Equation (42) provides the expression of the error as a function of circuit parameters: The analytical expression (67) can be used to find the worst-case error and the corresponding combination of network parameters. However, deriving the analytical expression for the maximum of δk p is difficult due to the group of influencing parameters cos ϕ, Moreover, it is necessary to describe constraints of network operation mode, restricting parameters like voltage in every point of the network, current ratings in the wiring, etc. The complexity of the problem related to the multi-parameter optimization in circuits with multiple operation mode constraints and parameter bounds yielded to the choice of its solution applying numerical optimization methods.

Objective Functions, Constraints and Bounds of Independent Parameters
In this research we seek to determine the systematic worst-case error (WCE)ˇˇδk pˇm ax of the method described in [3], considering a low voltage electrical distribution grid design and operation mode constraints. The value ofˇˇδk pˇm ax is important for the judgement about the precision and applicability of the method, while the parameters of grid and loads corresponding to the WCE are important for searching the improvement of the method precision. In [3], a wattmeter adjustment gain estimation uncertainty was preliminary explored by analyzing its random and systematic components. However, the systematic error was sampled only for some specific combinations of network and load parameters that were expected to cause the largest errors. In the current publication the search of the WCE is approached by utilizing global optimization techniques, that are perform extremum (WCE) location over the whole possible range of independent parameters. Because the method introduced in [3] is very new, its systematic WCE analysis is still absent.
The objective function used to find WCE ofˇˇδk pˇm ax is defined in (68), constraints function in (69) and parameter bounds in (70):ˇˇδ k pˇm ax " max where input vector X " pR W , R l , cos ϕq represents electrical parameters (R,L,C) of grid wiring and loads. Inequality (69) represents nonlinear constraints. In particular, H`X˘includes maximum rated current I Wmax through the wiring (Figure 4) and the maximum rated consumer power P lmax . Inequality (70) represents bounds of input parameters. X min is the vector representing the minimal possible values of circuit parameters. X max is the vector representing the largest acceptable values of circuit parameters. The implementation of the optimization procedure is presented by the flowchart in Figure 5. The procedure of the input vector X initialization is called randomPopulation(p_s) ( Figure 5).  In the following research the objective and constraint functions are implemented either using derived analytical error expressions (67) or by simulating a Simulink model shown in Figure 6. The second approach does not require analytical expression of the objective function and any grid can be modeled. The advantage of the second approach is that it can be easily adapted to find the worst-case error for a particular circuit. The drawback on the second approach is that the problem assigned is challenging in terms of computing requirements. Flow chart of the optimization by using genetic algorithm for both approaches presented in Figure 5. Both methods use the same logic except for calculation of the objective and constraint functions. Each time when the objective and constraint functions are called from the optimization algorithm, the first method calculates values of the objective and constraint functions applying the derived expressions. The second approach loads Simulink models, sets up network parameters of the models, solves the models, calculates and returns results of objective and constraint functions. To compute the result of the objective function, the Simulink model is solved twice. To reduce the computational load a single model corresponding SW1 and SW2 switch cases with two power measurement blocks for remote wattmeter powers P m1 and P m2 calculation is applied. To compute parameters of the constraint, function the Simulink model is solved once (SW1 case). Afterwards, the model with a circuit where the load R S is disconnected is called out (SW3 case). Implementation of the objective and constraint functions utilizing Simulink model are called objectiveN(P(t), LoadType) and constraintsN(P(t), LoadType). Correspondingly, the objective and constraint functions implementations utilizing analytical expressions are called objectiveA(P(t), LoadType) and constraintsA(P(t), LoadType).

Digital Optimization Techniques
The optimization problem stated in this paper was never solved before and a priori knowledge about the objective function local and global extremums is absent. Therefore, to cross verify the obtained solutions several strongly different optimization techniques should be employed. Also, it is obvious that the optimization problem is nonlinear, multivariable, and includes nonlinear constraints. genetic algorithm (GA), pattern search (PS) methods and nonlinear programing solver FMINCON [24] were selected to conduct the WCE search. Key features and their fit to the solution of the stated optimization problem are shortly summarized below.
The GA is a stochastic optimization method based on Darwin's theory of natural selection and the mechanism of population genetics. It initializes the first population randomly and through a process of evolution tends to find the most fit individual, i.e., combination of input parameter values that correspond to the maximum (minimum) value of the objective function. Each new generation produced by selection, crossover and mutation operators. GA method was successfully applied to solve nature related problems just for different goals [25]. The GA method is easy adaptable for a particular problem, does not require prioritization, scale, or weight objectives, and is an attractive alternative to other optimization methods because of its robustness in finding global optimum. Since the input values of the GA are taken recursively, the algorithm operates purposefully and, in many cases, would give the result faster in comparison to pure random search (PRS) or Monte Carlo simulation where inputs are chosen in purely random manner. Speed of an

Digital Optimization Techniques
The optimization problem stated in this paper was never solved before and a priori knowledge about the objective function local and global extremums is absent. Therefore, to cross verify the obtained solutions several strongly different optimization techniques should be employed. Also, it is obvious that the optimization problem is nonlinear, multivariable, and includes nonlinear constraints. genetic algorithm (GA), pattern search (PS) methods and nonlinear programing solver FMINCON [24] were selected to conduct the WCE search. Key features and their fit to the solution of the stated optimization problem are shortly summarized below.
The GA is a stochastic optimization method based on Darwin's theory of natural selection and the mechanism of population genetics. It initializes the first population randomly and through a process of evolution tends to find the most fit individual, i.e., combination of input parameter values that correspond to the maximum (minimum) value of the objective function. Each new generation produced by selection, crossover and mutation operators. GA method was successfully applied to solve nature related problems just for different goals [25]. The GA method is easy adaptable for a particular problem, does not require prioritization, scale, or weight objectives, and is an attractive alternative to other optimization methods because of its robustness in finding global optimum. Since the input values of the GA are taken recursively, the algorithm operates purposefully and, in many cases, would give the result faster in comparison to pure random search (PRS) or Monte Carlo simulation where inputs are chosen in purely random manner. Speed of an algorithm convergence is important for the problem because calculation of the objective function is time consuming especially for Simulink-based approach.
To verify the obtained solutions, the pattern search (PS) method is used in our research. PS is a derivative-free algorithm of a class of direct search (DS) methods for solving optimization problems [26]. The algorithm recursively searches a set of points around the current point, so called mesh, to find closer to optimum (lower/higher) value of objective function in each iteration. The mesh formed in the previous iteration by adding the current point to a scalar of vectors is called a pattern [26]. If the mesh improves the objective function of the current point, it becomes the new point, otherwise, the new mesh is formed, and calculations are repeated. The procedure is iterated until stopping criteria is met. A stand-alone PS algorithm tends to fall into local minima (maxima) [27]. According to the method application recommendations [27] to find the global extremum, multiple starting points were defined and implemented. The points are initiated in a random way from the defined parameter bounds. Half of the points are initiated for a set of parameters with a higher probability of an optimum point. Estimation of ranges of the higher optimum point probability are presented further in the section "Parameter Bounds". The DS methods as GA were successfully applied for some power system optimization tasks also [28,29].
Nonlinear programing solver FMINCON based on the interior point algorithm was the third optimization method implemented for problem solving. FMINCON is a solver used to find the minimum of the constrained multivariable function [24]. Contrary to the first two methods it is a gradient-based method that is designed to work on problems where the objective and constraint functions are both continuous and have continuous first order derivatives [24]. As well as PS, the FMINCON was utilized to search for the extremum point from set of different starting points. A fraction of them was chosen from region of higher probability of the optimum point.
Usually optimization solvers are designed to search for an objective function minimum. However, according to the definition of WCE in Equation (68) it is required to find the maximum of the objective function. Therefore, the objective function is transformed into the inverse function:ˇδ k pˇm ax " min´ˇˇ`δk p`X˘˘´1ˇ¯.
Implementation cost and complexity of GA, PS and FMINCON are not considered in this publication because MATLAB toolboxes were used in research. Objective and constraints functions were implemented in MATLAB script language according to Equations (28), (67) and (75) and related expressions.

Electrical Parameters Calculation
The input parameters of the objective function are: load resistance R l , wire resistance R W and the load power factor cos ϕ. Since the inductive and capacitive loads have the opposite signs of reactance, there is an infinite amount of possible combinations of inductance L l and capacitance C l that meet particular load power factor cos ϕ for the fixed R l value. Therefore, submission of the load inductance L l and capacitance C l to objective function as input parameters are fated to poor convergence of the solution. Therefore, the power factor cos ϕ was chosen as an input parameter to define the reactive load. While the linear loads could be RL or RC type, separate GA searches were established for both cases allowing just RL type load for one case and just RC for another. Parameters L l and C l are calculated inside of the objective function from the resistance R l and cos ϕ according to the following formulas: R l tanparccospcos ϕqq ω , C l " 10 10 F; i. e. , X Cl « 0 Ω. (72) where ω " 2π f is angular frequency. Expressions (72) were used for the case of RL type load, while expressions (73) for the case of RC type load. Inductance of the wiring cable L W is calculated inside of the objective function according to expression: where R W1m is a resistance of 1 meter cable length, L W1m is 1 meter inductance of the same cable. Expressions (72)-(74) were used to calculate load and wiring electrical parameters in the optimization algorithm explained in Figure 5.

Nonlinear Inequality Constraints
The used nonlinear inequality constraints denoted H`X˘in Equation (69) are: Load resistance R l , wire resistance R W and load power factor cos ϕ are the input parameters of the constraint function (75) as well as for the objective function. Parameters L l , C l and L W are calculated by applying formulas (72)-(74). When the analytical expression of the objective function is used, the value I Ws is calculated by applying the Equation (55). Calculation of value P ls is defined in Equation (63). When Simulink model ( Figure 6) is used, the I Ws and P ls " P m2 are estimated by solving the circuit presented in Figure 4.
The current I Wmax is set according to the calculated value for a particular defined maximum rated load power P lmax nominal voltage U N and power factor cos ϕ N : where P lmax is the maximal consumer active power for a particular case under consideration. A typical power factor for residential load cos ϕ N " 0.9 is used in this case. The nominal voltage U N " 230 V.

Parameter Bounds
Parameter bounds in (70) X min˜Xmax are vectors of a circuit and operation mode parameters defining a range of values for independent input parameters: The upper bound of a network wiring active resistance R Wmax is calculated according to the standard limits of allowable voltage drop in low voltage network. In many countries voltage drop in low voltage distribution network cannot exceed 5% [30]. Since for the fundamental frequency X W ! R W , the inductance of the cable could be neglected. To calculate upper bound of the wire resistance, the simplified expression is used: The lower bound of wire active resistance is accepted to be: Bounds for R l were chosen in connection to nominal network voltage and the maximum rated power measured by the RemW: where P lmin " 1 W in all of the following simulations (corresponds to R lmax " 52, 900 Ω), and P lmax is specified along with the presented results in the next chapter, for example P lmax " 10 kW with power factor cos ϕ min " 0.1, corresponds to R lmin " 0.0529 Ω.
The upper power factor bound of cos ϕ in Equation (77) corresponds to pure active load (cos ϕ " 1). The lower bound for the power factor cos ϕ min is in the range from 0.1 to 0.9 and is dependent on the simulation settings specified in the next section.

Initial Points Selection
To improve a convergence to a solution, some fraction of initial population (for GA) or a fraction of initial points (for PS, FMINCON) were forced to meet the parameter ranges with higher probability of the optimum point. To determine the parameter range, the objective function was calculated by changing one parameter by fixed step for a wide range while the rest of the parameters were kept constant. The results are presented in Figures 7-9.

Initial Points Selection
To improve a convergence to a solution, some fraction of initial population (for GA) or a fraction of initial points (for PS, FMINCON) were forced to meet the parameter ranges with higher probability of the optimum point. To determine the parameter range, the objective function was calculated by changing one parameter by fixed step for a wide range while the rest of the parameters were kept constant. The results are presented in Figures 7-9.          According to the results, bounds for a half of initial population defined as follows: R lmin ď R l ď 50 Ω, cos ϕ min ď cos ϕ ď 1.

(82)
All members of the initial population (for GA case) or initial points (for PS, FMINCON cases) were initialized randomly. During the initialization half of the initial population (for GA) or initial points (for PS, FMINCON) were forced to meet the bounds defined by Equations (77)-(81). The second half of the initial population (initial points) was forced to meet the bounds defined in Equation (82).

The Network Parameters Corresponding to the Worst-Case Error
Following the aim of the research to identify the WCE and the corresponding network and load parameters, a combination of typical loads, wire resistances, and power factor values were used to search forˇˇδk pˇm ax according to Equations (68)-(70). Though the input vector in (70) is X " pR W , R l , cos ϕq, the most often measured parameters by the metering instruments are load power P l and power factor cos ϕ. Also, a distribution network design recommendation include restrictions on power losses and transferred power ratio. The results presented in this chapter are aimed to reveaľˇδ k pˇm ax dependence on P l , cos ϕ and P w {P l . The complete set of the obtained results are shown in Tables 1 and 2, while some of dependencies are plotted graphically in Figures 10-13.

Method
Load Type  Table 2. The adjustment gain relative error dependence on the maximum rated load power 5 kW, 7 kW, and 10 kW. another vector of independent variables yielding the same maximum error can only be important for understanding the worst case network configuration but does not affect the characterization of the remote estimation method.
To visualize the influence of the wiring resistance on the adjustment gain error, the GA search was completed for the upper bounds 1 = 0.75 , 2 = 0.5 and 3 = 0.25 in addition to = 0.952 Ω. The obtained results are plotted in Figure 10. The results of presented research indicate that the WCE values were found for the highest wire impedance values = . The trend of the increasing error with each increment of the wiring impedance supports the assumption that the source of the error is the wiring impedance between RefW and RemW. The assumption could be made by analyzing the Equations (36), (42), (43) and (67). Zero wiring impedance ( = 0 Ω and = 0 Ω) results in zero values of both ∆ ( ) and ∆ (∆ 2 ), and according to the equation (67), consequently to = 0. On the other hand, it could be a future research direction to estimate wiring losses on-line and look for the modification of the adjustment gain estimation method by including information about the wiring resistance or losses and expecting to reduce its WCE.
According to Table 1, the WCE in case of the active load cos( ) = 1 corresponds to the smallest load (the upper bound of ⋆ ), but WCE corresponds to the largest load (the lower bound of ⋆ ) when cos( ) ≠ 1. To investigate the haziness a more detail exploration of WCE dependency on power factor is conducted in a range of load power factor 0.9 ≥ cos( ) ≥1 and plotted in Figure 11. The explanation to the form of the dependence in Figure 11 can be given by analyzing Equations (42)-(45). The results presented in Figure 11 show that by incrementing the minimal , the WCE is found corresponding to lowest load ( = ) and (∆ ) dominates over (∆ ). Despite this discontinuity of the dependence (cos( ) ), the modulus of as it can be seen from Table 1 and Figure 10, is always decreasing in response to increase of the load power factor cos( ) . Figure 11. Dependence of the wattmeter gain WCE on power factor in the range of a load power factor 0.9 ≥ cos( ) ≥1.  The simulation results corresponding to = 5, 7 and 10 kW are presented in Table 2. It can be seen from Table 2 that the load power factor reduction always degrades the error | | despite the largest allowed load power . The comparison of the results in Table 1 and Table 2 reveals that the error values are similar for the same load power factor in the whole explored the minimal power factor range (0.1 ≥ cos ( ) ≥ 0.9). It could be explained by the fact that the product • = Δ = 0.05 • is constant. Cable impedance must be less, while the maximal current is relatively higher for higher powers to meet the requirement of the largest allowable voltage drop (5%). Equal or lower than the allowable impedance of the wires must be assured by selecting proper diameter of wires during the electrical distribution grid design stage. The results also show that the increase of gives the same impact on results as proportional relative increase of for the power factor bounds resulting in negative error values. Summarizing up the results, it can be clearly stated that the WCE is obtained from combination of the highest wiring impedance ( , ), the highest allowable current , and the lowest load power factor cos( ) ≤ cos ( ) values. To visualize the dependence of the error modulus on load value , which is calculated according to Equation (62), the direct error value calculations are carried out using Equation (67) and Figure 12 is plotted. Figure 12 indicates that a real error | | in case of a particular can be significantly less than | | which corresponds to the maximum considered load. The real error value could be estimated by measuring and cos( ) and utilizing expression approximated from data shown in Figure 12. However, since | | is less that the target accuracy of 1%, in the current stage of the development, we did not attempt to derive the precise error estimation but only focused on its upper limit | | . In Figure 13 there are presented WCE and load power factor dependencies on power loss to load power ratio ⋆ / ⋆ . Power loss to load power ratio is often accepted as selection criteria for power distribution grid cables. It can be seen that decreasing the power factor of a load herewith also increasing loss to load power ratio, causes the increase of the wattmeter gain estimation error.

Requirements upon Computational Resources
The aim of this section is to compare the computational cost of all applied global search algorithms combined with the methods of objective function implementation. The analytical derivation of the expression of WCE is based on the mentioned approximations. In opposite, the Simulink model based objective function does not include any approximate assumptions. Therefore, it is of important to estimate if using the analytical expression is reasonable in a sense of speeding up calculations. For the case of analytically expressed objective function optimization it is of interest to discover which of algorithm (GA, PS, FMINCON) perform faster in finding the same solution. The objective function is unique because it is derived from new method equations introduced in [3]. Therefore, it was not solved by someone else and computational cost reported in the analysis below cannot be compared over data of any other references.
The obtained wattmeter adjustment gain WCE values along with calculation times and population size of GA or number of initial points of PS are presented in Table 3. All computations were performed on Dell PowerEdge R730 server with two Intel(R) Xeon(R) CPU E5-2620 v3 @ 2.40 GHz processors, together containing 12 physical cores and 24 logical cores with 32 GB operating memory. The power management profile of the computer operational system was set to performance level during the simulations. The computations were implemented and performed using MATLAB Parallel computing toolbox. The comparative speed analysis was carried out for a load power factor within the range 0.9 ≥ cos ≥ 1. The GA search based on objective and constraint functions implemented applying analytical expressions yielded the result much faster in comparison to Simulink based implementation of Simulation results in the case of the maximum load power P max " 2.5 pkWq and pure active load (R), active-inductive (RL) and active-capacitive (RC) type load are presented in Table 1. The set of selected parameters is considered as typical for the case when the reference meter is a revenue meter and the meter under adjustment is a sub-accounting meter like a smart socket in a residential apartment.
The differences of the error obtained by different in their nature search methods GA, PS and FMINCON are marginal. The differences of the error in cases of RL type and RC type load are also minor with a trend for RC type loads to produce slightly higher error values. The practical meaning of the differences is negligible, therefore only results of the error estimations for the more common RL type loads are presented further.
It is the value of WCE |δk P | max is the most important parameter characterizing remote adjustment gain estimation method accuracy. The same goal attained by three different search methods ensures that the global maximum was found. The possibility that there might exist yet another vector of independent variables yielding the same maximum error can only be important for understanding the worst case network configuration but does not affect the characterization of the remote estimation method.
To visualize the influence of the wiring resistance on the adjustment gain error, the GA search was completed for the upper bounds R Wmax1 " 0.75 R Wmax , R Wmax2 " 0.5 R Wmax and R Wmax3 " 0.25 R Wmax in addition to R Wmax " 0.952 Ω. The obtained results are plotted in Figure 10.
The results of presented research indicate that the WCE values were found for the highest wire impedance values R W " R Wmax . The trend of the increasing error with each increment of the wiring impedance supports the assumption that the source of the error is the wiring impedance between RefW and RemW. The assumption could be made by analyzing the Equations (36), (42), (43) and (67). Zero wiring impedance (R W " 0 Ω and X W " 0 Ω) results in zero values of both ∆P W pP S q and ∆P l p∆U 2 q, and according to the equation (67), consequently to δk p " 0. On the other hand, it could be a future research direction to estimate wiring losses on-line and look for the modification of the adjustment gain estimation method by including information about the wiring resistance or losses and expecting to reduce its WCE.
According to Table 1, the WCE in case of the active load cospϕq min " 1 corresponds to the smallest load (the upper bound of R ‹ l ), but WCE corresponds to the largest load (the lower bound of R ‹ l ) when cospϕq min ‰ 1. To investigate the haziness a more detail exploration of WCE dependency on power factor is conducted in a range of load power factor 0.9 ě cospϕq min ě1 and plotted in Figure 11. The explanation to the form of the dependence in Figure 11 can be given by analyzing Equations (42)-(45). The results presented in Figure 11 show that by incrementing the minimal bound of power factor leads to reduction of the error modulus, while δk Pmax value always remains negative. In this range the error component related to the load power change δk p p∆P l q dominates over the error component related to change of wiring losses δk p p∆P W q (43). By reaching some threshold power factor cospϕq minTHR , the sign of δk Pmax turns positive. In the range cospϕq min ą cospϕq minTHR , the WCE is found corresponding to lowest load (R l " R lmax ) and δk p p∆P W q dominates over δk p p∆P l q. Despite this discontinuity of the dependence δk pmax`c ospϕq min˘, the modulus of δk Pmax as it can be seen from Table 1 and Figure 10, is always decreasing in response to increase of the load power factor cospϕq min .
The simulation results corresponding to P max " 5, 7 and 10 kW are presented in Table 2. It can be seen from Table 2 that the load power factor reduction always degrades the error |δk P | max despite the largest allowed load power P max . The comparison of the results in Tables 1 and 2 reveals that the error values are similar for the same load power factor in the whole explored the minimal power factor range (0.1 ě cospϕq min ě 0.9). It could be explained by the fact that the product I Wmax¨RWmax " ∆U max " 0.05∆U N is constant. Cable impedance must be less, while the maximal current is relatively higher for higher powers to meet the requirement of the largest allowable voltage drop (5%). Equal or lower than the allowable impedance of the wires must be assured by selecting proper diameter of wires during the electrical distribution grid design stage. The results also show that the increase of R Wmax gives the same impact on results as proportional relative increase of I Wmax for the power factor bounds resulting in negative error values. Summarizing up the results, it can be clearly stated that the WCE is obtained from combination of the highest wiring impedance (R Wmax , L Wmax ), the highest allowable current I Wmax , and the lowest load power factor cospϕq ď cospϕq min values.
To visualize the dependence of the error modulus on load value P l , which is calculated according to Equation (62), the direct error value calculations are carried out using Equation (67) and Figure 12 is plotted. Figure 12 indicates that a real error |δk P | in case of a particular P l can be significantly less than |δk P | max which corresponds to the maximum considered load. The real error value could be estimated by measuring P l and cospϕq and utilizing expression approximated from data shown in Figure 12. However, since |δk P | max is less that the target accuracy of 1%, in the current stage of the development, we did not attempt to derive the precise error estimation but only focused on its upper limit |δk P | max .
In Figure 13 there are presented WCE and load power factor dependencies on power loss to load power ratio P ‹ W {P ‹ l . Power loss to load power ratio is often accepted as selection criteria for power distribution grid cables. It can be seen that decreasing the power factor of a load herewith also increasing loss to load power ratio, causes the increase of the wattmeter gain estimation error.
This dependence is nonlinear and if P ‹ W P ‹ l p%q ă 5% the error value is much below the 1% target. Even in the case of abnormal level of power losses in the network (Figure 4) wires, the method yields wattmeter adjustment gain estimation error less than 0.5%.

Requirements upon Computational Resources
The aim of this section is to compare the computational cost of all applied global search algorithms combined with the methods of objective function implementation. The analytical derivation of the expression of WCE is based on the mentioned approximations. In opposite, the Simulink model based objective function does not include any approximate assumptions. Therefore, it is of important to estimate if using the analytical expression is reasonable in a sense of speeding up calculations. For the case of analytically expressed objective function optimization it is of interest to discover which of algorithm (GA, PS, FMINCON) perform faster in finding the same solution. The objective function is unique because it is derived from new method equations introduced in [3]. Therefore, it was not solved by someone else and computational cost reported in the analysis below cannot be compared over data of any other references.
The obtained wattmeter adjustment gain WCE values along with calculation times and population size of GA or number of initial points of PS are presented in Table 3. All computations were performed on Dell PowerEdge R730 server with two Intel(R) Xeon(R) CPU E5-2620 v3 @ 2.40 GHz processors, together containing 12 physical cores and 24 logical cores with 32 GB operating memory. The power management profile of the computer operational system was set to performance level during the simulations. The computations were implemented and performed using MATLAB Parallel computing toolbox. The comparative speed analysis was carried out for a load power factor within the range 0.9 ě cos ϕ ě 1. The GA search based on objective and constraint functions implemented applying analytical expressions yielded the result much faster in comparison to Simulink based implementation of objective and constraints functions despite for the 100 times larger search population size. The higher population size caused slightly higher extremum value with much higher probability of the result being a global minimum. PS implementation applying analytical expressions gave a very close result in noticeably lower time in comparison to the same implementation of the GA method. FMINCON based search implementation of analytical expressions has shown a noticeable speed advantage in comparison to the rest methods and gave very close result as can be seen from Table 3. Therefore, the results in Table 3 support the statement, that using analytical expression does not lead to the loss of search precision but significantly reduce computational cost. Also, FMINCON seems to obtain the same solution of the optimization problem (66) several times faster than GA and PS methods.

Conclusions
The systematic worst-case error (WCE) of a previously introduced wattmeter adjustment gain remote estimation method was studied using the simplified electrical energy delivery model enabling to derive an analytical expression of the error (objective function) and constraint functions. To find the WCE, the error expression was maximized using two stochastic search techniques (genetic algorithm, pattern search) and the nonlinear programming solver. The adequacy of the error analytical expression and its maximization solution was verified by the second approach based on numerical simulation of electrical circuit model. Both solutions obtained produced very close results with practically insignificant differences. The second approach utilizing numerical electrical circuit model universal in a sense that it enables to find the WCE together with the corresponding electrical network configuration without analytical expression of the objective function. However, the requirements for computational cost is many times higher for the second approach. For the considered network model the computation time required to reach the solution was over 200 times longer to the search of WCE using the derived analytical expression. Therefore, it we the second approach is recommended for the verification of the correctness of problem solution, but analytical expression is reasonable for The research has revealed that the systematic error of the method is influenced by the wattmeter monitored load active power, its power factor and losses in the network wiring between reference and wattmeter under test. It was found that the WCE is below 0.5% if the load power factor is larger than 0.1. The WCE is less than 0.1% for a load exhibiting power factor larger than 0.9 in a distribution network designed following the recommendations to ensure the grid voltage magnitude drop not more than 5%. The WCE of the proposed method is a conservative characterization of error range. Therefore, the systematic error of the described remote adjustment gain estimation method does not exceed error limits of Class 2 and Class 1 watthour revenue meters.