OLTC includes two parts, “transformer” and “tap changer.” Unlike traditional transformers, an OLTC sets the tap position as an unknown variable during modelling. By controlling the amplitude of the secondary side voltage to meet a certain dead band range, the OLTC adjusts the tap position to maintain the voltage level at the load center within a certain error range, thereby enhancing the power quality for the user. This paper proposes a droop control model for the OLTC, where the amplitude of the secondary side voltage and the tap turns ratio of the OLTC align with the droop control curve.

Define the column name of the nodal association matrix to represent the node i − a , i − b , i − c , j − a , j − b , j − c ,the line name corresponds to the branch connected to i − a ,the branch connected to j − a ,the branch connected to i − b ,the branch connected to j − b ,the branch connected to i − c ,the branch connected to j − c ,where i 、 j represent node; a 、 b 、 c represented by three phases of each node. The nodal correlation matrix is as shown in Equation 19: n o d e C = b r a n c h 1 − 1 1 1 − 1 1 − 1 1 1 (19)

Where represents the association between a branch and a node, with its direction flowing out of the node; and represents the association between a branch and a node, with its direction flowing into the node.

Taking phase A as an example, according to the law of energy conservation on the primary and secondary sides of the transformer, the voltage and current on the winding satisfy the following relationship: I ˙ i − a I ˙ j − a = − 1 t a U ˙ i − a U ˙ j − a − I ˙ j − a / y a = t a (20)

Here, I ˙ i − a and I ˙ j − a represent the branch currents flowing through the branches connected by i − a and j − a , respectively. y a = 1 / R a + j X a represents the equivalent internal impedance of phase A of the transformer, t a represents the turns ratio of phase A.

The relationship between branch current and node voltage is expressed as shown in Equation 21: I ˙ i − a = − 1 t a I ˙ j − a = − 1 t a U ˙ j − a − U ˙ i − a t a y a = y a t a 2 U ˙ i − a − y a t a U ˙ j − a I ˙ j − a = U ˙ j − a − U ˙ i − a t a y a = − y a t a U ˙ i − a + y a U ˙ j − a (21)

Equation 20 is expressed as a matrix as shown in Equation 22: I ˙ i − a I ˙ j − a = y a t a 2 − y a t a − y a t a y a U ˙ i − a U ˙ j − a (22)

Similarly, the relationship between the branch current and node voltage for phases B and C is derived, and the relationship between the three-phase branch current and voltage is ultimately obtained as follows: I ˙ = Y pr U ˙ Y pr = y a t a 2 − y a t a − y a t a y a y b t b 2 − y b t b − y b t b y b y c t c 2 − y c t c − y c t c y c (23)

In Equation 23: I ˙ = I ˙ i − a I ˙ j − a I ˙ i − b I ˙ j − b I ˙ i − c I ˙ j − c T represents the branch current matrix of the OLTC, U ˙ = U ˙ i − a U ˙ j − a U ˙ i − b U ˙ j − b U ˙ i − c U ˙ j − c T represents the node voltage matrix of the OLTC; Y pr is defined as the original admittance matrix; t b , t c respectively represent the B and C phase transformation turns ratios.

By combining the original admittance matrix with the node-branch incidence matrix, the node admittance matrix is obtained from the following equation: Y b u s = C T Y pr C

The branch power connected to each node is calculated using the node injection power expression as shown in Equations 24 and 25: P i − φ Reg = U i − φ ∑ β = A , B , C U j − β G i j − φ β Reg ⁡ cos ⁡ θ i j − φ β + B i j − φ β Reg ⁡ sin ⁡ θ i j − φ β (24) Q i − φ Reg = U i − φ ∑ j = A , B , C U j − β − B i j − φ β Reg ⁡ cos ⁡ θ i j − φ β + G i j − φ β Reg ⁡ sin ⁡ θ i j − φ β (25) G i j − φ β Reg and B i j − φ β Reg represent the real and imaginary parts of the corresponding elements in the node admittance matrix Y b u s of the on-load tap changer, G i j − φ β Reg represents the mutual conductance between node i − φ and node j − β , B i j − φ β Reg represents the mutual susceptance between node i − φ and node j − β ; The row name and column name of Y b u s are expressed as: node i − a , i − b , i − c , j − a , j − b , j − c .

The following describes the modelling of the regulator part in Figure 2.

The tap position of OLTC are regulated through a control circuit. When voltage control, direct control cannot be carried out on the high-voltage circuit. Therefore, voltage and current transformers are used to construct a simulated circuit - the control circuit. The control circuit is a proportional model of the actual line. For example, if the actual line transformer has a secondary side rated at 2.4 kV, the control circuit operates at a voltage level of 120 V. By monitoring the voltage of the control circuit, the voltage at the load center of the actual line can be determined. If the voltage is below the normal level, indicating that the voltage at the loading center is below the position of the normal level, the tap changer is adjusted to raise the voltage at the loading center. For three-phase supply voltages below 10 kV, the allowable deviation is 7% of the rated voltage. For single-phase supply voltages of 220 V, the allowable deviation is +7% and −10% of the rated voltage. R k j , e q and X k j , e q are the proportional impedances of the actual line in the control circuit and are typically known quantities.

For the control circuit, the current rating of the primary winding of the current transformer is set as C T P , and the current rating of the secondary winding is set as C T S (usually taken as 5). The voltage transformer turns ratio is set as N P T . First, according to Ohm’s law, the equivalent impedance of the three-phase line is calculated as shown in Equation 26: R l i n e + j X l i n e = U ˙ k − U ˙ j I ˙ k j (26)

It should be noted that the equivalent impedance of the three-phase line is not the actual line impedance; it is the turns ratio of the actual measured voltage on the secondary side to the load center voltage difference to the load side current under the rated turns ratio.

The equivalent impedance of the line compensator is calculated from the equivalent impedance of the three-phase line as R k j , e q 、 X k j , e q using the formula as shown in Equation 27: R e q , k j + j X e q , k j = R l i n e + j X l i n e · C T P N P T · C T S (27)

The current in the compensator branch is obtained from the actual line current I ˙ k j and the current transformer I ˙ k j , e q . I ˙ k j , e q = I ˙ k j · C T S C T P (28)

In Equation 28: I ˙ k j = I ˙ k j − a I ˙ k j − b I ˙ k j − c T represents the three-phase current in the line.

The voltage difference U ˙ d r o p of the compensator branch is obtained as shown in Equation 29: U ˙ d r o p = R k j , e q + j X k j , e q · I ˙ k j , e q (29)

Finally, the control voltage of the compensator branch is as shown in Equation 30: U ˙ r = U ˙ 3 N P T − U ˙ d r o p (30)

Applying the droop control to the control voltage and tap changer turns ratio, taking phase A as an example, the principle of droop control for the OLTC is introduced. Define Δ t a , Δ t b , Δ t c as the adjustment amounts of the tap changer turns ratio for phases A, B, and C relative to their respective initial turns ratios, and add to the objective function as shown in Equation 31: min   Δ t a 2 + Δ t b 2 + Δ t c 2 (31)

The physical meaning of the above equation is that when the control voltage is within the dead zone, the OLTC tap position remains unchanged.

The Δ t a − U r − a droop control function of the variation in the variable turns ratio variation and the control voltage is shown in Equation 31: Δ t a = Δ t max − Δ t 0   U min ≤ U r − a < U l k d r 1 U r − a − U d b l   U l ≤ U r − a ≤ U d b l 0   U d b l < U r − a < U d b h k d r 2 U r − a − U d b h   U d b h ≤ U r − a ≤ U h t min − Δ t 0   U h < U r − a ≤ U max (32)

The relationship between Δ t a and the control voltage U r − a satisfies the following curve in Figure 3:

Where, Δ t max and Δ t min respectively represent the upper and lower limits of the variation in turns ratio Δ t ; Δ t 0 represents the change in tap position corresponding to the initial turns ratio; U d b h and U d b l respectively represent the upper and lower bounds of the voltage dead zone; U h and U l respectively represent the upper and lower bounds of the voltage droop control curve; U max and U min respectively represent the upper and lower bounds of the bus voltage magnitude; k d r 1 and k d r 2 respectively represent the droop control coefficients, and k d r 1 < 0 , k d r 2 < 0 . When the voltage U r − a is within the dead zone range, the OLTC does not actuate, and the variation in turns ratio is 0.

The relationship between the droop control coefficient and the variation in the turns ratio and the voltage is as shown in Equation 33: k d r 1 = Δ t max U l − U d b l k d r 2 = − Δ t min U d b h − U h (33)

From Figure 3, it can be seen that the characteristic curve represents a piecewise function. When the system operates near the turning point, the derivative is discontinuous, the algorithm search direction is uncertain, and it is difficult to smoothly switch operating curves, so the piecewise droop control function has strong non-smooth characteristics and is difficult to solve. Usually, mixed integer nonlinear programming is used to describe piecewise function control characteristics, but this method has low computational efficiency. In this paper, a fitting function is used to approximate the piecewise function, and the fitted droop control function is shown in Equation 34: Δ t a = Δ t max + k d r 1 α ln 1 + e α U r − a − U l − ln 1 + e α U r − a − U d b l + k d r 2 α ln 1 + e α U r − a − U d b h − ln 1 + e α U r − a − U h (34)

Where, α is the fitting coefficient, set α = 500 in this paper.

Approximately fitting the piecewise droop control function can transform it into a smooth function that is continuously differentiable, as shown in Figure 4. The smoothing function can avoid sudden changes in derivative order during algorithm iteration calculations, thus improving the convergence of the algorithm.

The variation of the ratio of the voltage regulator and the tap variable satisfies the relationship is shown in Equation 34: Δ t a = x d − a · Δ d (35)

In the equation, x d − a is the phase A tap-changer variable; Δ d is the turns ratio change corresponding to a 1-tap adjustment, which is typically taken as Δ d = 0.00625   p . u . in distribution networks.

In this paper, the discrete tap-changer variable is first treated as a continuous variable for state estimation calculations to obtain a continuous optimal solution for state estimation. Subsequently, a positive curvature quadratic penalty function is introduced in the objective function, and the state estimation is continued using the continuous optimal solution as the initial value to obtain an integer solution for the tap-changer setting.

Taking phase A as an example, the positive curvature quadratic penalty function is shown in Figure 5. x d 0 、 x d 1 、 x d 2 are any three continuous discrete component values. Define the neighborhood R x d − a as shown in Equation 36: R x d − a = x d − a | x d 1 − 1 2 Δ d ≤ x d − a ≤ x d 1 + 1 2 Δ d (36)

Where, x d 1 is the neighborhood center, which is the closest discrete grading value determined according to the continuous optimal solution.

In the optimization process, when the value of x d − a is in the neighborhood of the above definition, the penalty function ε x d − a is introduced: ε x d − a = 1 2 μ d x d − a − x d 1 2 (37)

In the equation, μ d is the penalty factor, a known quantity; the penalty function will force x d − a within the neighborhood to approach the neighborhood center. From Equation 37, we can obtain the first and second derivatives of the quadratic penalty function within the neighborhood. ∂ ε x d − a ∂ x d = μ d x d − a − x d 1 ∂ 2 ε x d − a ∂ x d 2 = μ d (38)

From Equation 38, it can be observed that for the Newton method, the introduction of a penalty function in the objective function will result in the incremental inclusion of the first and second derivative terms in the Jacobian matrix and the Hessian matrix during the optimization iteration. Through this procedure, not only can rapid convergence of OLTC tap positions be achieved, but also the positive definiteness of the iteration matrix can be strengthened, thereby improving algorithm convergence. From a perspective in the field of power systems, this approach enables efficient adjustment of OLTC tap positions and enhances algorithm convergence by modifying derivative terms in iterative matrices.

Nonlinear constraints in the model of OLTC mainly fall into two categories: the first category is branch power constraints; the second category is droop control nonlinear constraints. Due to the introduction of the unknown variable t of the voltage regulator in the branch power constraint and the logarithmic function in the droop control nonlinear constraints, the aforementioned physical linearization methods are no longer applicable. In this paper, it is proposed to utilize a first-order Taylor expansion for physical linearization, followed by using Partial Least Squares Regression (PLSR) to calculate compensation errors, resulting in the development of a hybrid data-physical-driven Linearization model.

Firstly, introduce the principle of linearization of power constraints. For ease of description, define the relationship between the branch power of the transformer and the variables as shown in Equation 39: P i − φ Reg = f i − φ P X Reg Q i − φ Reg = f i − φ Q X Reg X Reg = θ i − ABC U i − ABC θ j − ABC U j − ABC t ABC T (39)

Where X Reg represents the variables of the on-load tap-changing transformer, composed of the phase angle and voltage magnitude of the primary and secondary sides, and the three-phase turns ratio; f i − φ P and f i − φ Q respectively represent the functional relationships between the active and reactive power of phase φ branch and the variables of the on-load tap-changing transformer.

The first-order Taylor expansion is performed at the initial point X Reg , 0 of the variable, as shown in Equation 40: f i − φ P X Reg ≈ f i − φ P X Reg , 0 + f i − φ P ′ X Reg , 0 · X Reg − X Reg , 0 + Δ P i − φ Reg f i − φ Q X Reg ≈ f i − φ Q X Reg , 0 + f i − φ Q ′ X Reg , 0 · X Reg − X Reg , 0 + Δ Q i − φ Reg (40)

In the equations, f i − φ P ′ and f i − φ Q ′ represent the partial derivatives of the active and reactive power functions of the phase φ branch with respect to each variable in X Reg ; Δ P i − φ Reg and Δ Q i − φ Reg respectively represent the compensating errors of the active and reactive power of the phase φ branch, and the compensating errors are also calculated using the PLSR method.

Next, the linearization principle of the droop control function of the OLTC is introduced. Similar to the linearization principle of branch power, the functional relationship between the turns ratio change and the control voltage is defined in Equation 41: Δ t = f Δ t U r (41)

Where, f Δ t represents the functional relationship between the variable ratio variation and the secondary side control voltage.

The first-order Taylor series expansion is performed at the initial value point U r , 0 of the variable, as shown in Equation 42: f Δ t U r ≈ f Δ t U r , 0 + f Δ t ′ U r , 0 · U r − U r , 0 + Δ t r (42)

In the formula, f Δ t ′ represents the derivative of the variable turns ratio variation function to the variable U r ; Δ t r represents the compensation error of the variable turns ratio variation.

Then, the PLSR algorithm is used for data-driven error compensation, the specific principle and mathematical expressions are as follows:

The compensation error Δ y of active power is defined in Equation 43: Δ y = ξ y ′ + η Δ y = Δ P P V − φ m Δ Q P V − φ m T y ′ = P l o a d P V − φ Q l o a d P V − φ T (43)

Among them, ξ is the coefficient matrix and η is the constant matrix, which are obtained by PLSR fitting is the independent variable matrix of load composition; P l o a d P V − φ and Q l o a d P V − φ represent the vectors composed of active and reactive loads of the grid-connected node φ phase of the system.

The independent variable set and the dependent variable set can be obtained from the results of the power flow calculation. After obtaining the set of independent variables R and the dependent variable set Z, it is standardized. R * = R − R ¯ S R − 1 Z * = Z − Z ¯ S Z − 1 (44)

In Equation 44: R * , Z * represent the standardized independent variables and dependent variable sets; R ¯ and S R respectively represent the mean and standard deviation of the independent variable set R.

For the standardized independent variable and dependent variable sets, the least squares regression analysis method is used to calculate the fitting coefficients, with the knowledge in the field of power systems: C = PLS R * , Z * (45) η = Z ¯ − R ¯ S R C ⊙ S Z (46)

In Equations 45 and 46: C represents the regression coefficient matrix obtained by PLSR. ξ and η represent the coefficient matrix and the constant matrix in the linear regression equation, respectively. Thus, the compensation error Δ y can be obtained. Finally, the correction equation of the measurement is shown in Equation 47: f P V − φ m − p X P V ≈ f P V − φ m − p X P V 0 + ∇ f P V − φ m − p X P V 0 × X P V − X P V 0 + Δ P P V − φ m f P V − φ m − Q X P V ≈ f P V − φ m − Q X P V 0 + f P V − φ m − Q X P V 0 × X P V − X P V 0 + Δ Q P V − φ m (47)

According to the two-stage bad data identification model proposed in this paper, 10 bad data are set as shown in Table 3. At the same time, the false alarm rate is set to P e = 0.0025 . By querying the standard normal distribution table, the normal range of normalized residuals can be obtained, and the detection threshold τ is set to 3 and the threshold D for the projection statistics is set to 1.

According to Table 4, among the 10 bad data, the projection statistical method effectively identified only 8 pcs bad data, while all other measurements with projection statistics values greater than 1 resulted in false alarms. Similarly, the maximum normalized residual method also failed to accurately identify the bad data. Therefore, both traditional statistical methods for bad data identification have certain limitations.

Due to the potential for misjudgments and missed detections by traditional statistical methods when multiple measurement bad data points are present in the system, these methods may not accurately identify measurement bad data. In contrast, the MILP bad data identification model used in the second stage can accurately identify all measurement bad data in the set identified by traditional methods when misjudgments or missed detections occur in the first stage. Therefore, in the first stage of this paper, all 265 lever measurements in the system were placed into a suspicious measurement set. The maximum normalized residual method was then used to identify the remaining 137 measurements, thereby validating the accuracy of the proposed bad data identification model and its efficiency compared to the single-stage model. Four sets of tests were conducted to compare the bad data identification method proposed in this paper with the traditional MINLP-based bad data identification method:

Test 1: Compression of the suspicious measurement set using the MILP-based two-stage bad data identification method proposed in this paper;

Test 2: Compression of the suspicious measurement set using the traditional MINLP-based two-stage bad data identification method;

Test 3: No compression of the suspicious measurement set, directly using the MILP-based single-stage bad data identification method;

Test 4: No compression of the suspicious measurement set, using the traditional MINLP-based single-stage bad data identification method.

Models were built for both methods in the GAMS optimization software, different solvers were called for calculation, and the results were compared and analyzed.

According to Table 5, for the modified IEEE-33 node test feeder case chosen in this paper, the two-stage model reduces the number of 0/1 integer variables and improves identification accuracy. However, for traditional MINLP-based bad data identification models, feasible solutions could not be obtained with most solvers, with correct results only achievable using the BONMIN solver. In contrast, the MILP-based bad data identification model achieved accurate identification results with most solvers. In terms of identification efficiency, even with the two-stage model, the MINLP identification model’s solving time reached the computational limit of 1,000 s, and the solver exited abnormally, indicating limited applicability. For the MILP-based bad data identification method, using the two-stage model, the SCIP solver improved solving efficiency by approximately 23 times, and the CPLEX solver achieved the highest efficiency, requiring only 0.546 s.