Convex Polar Second-Order Taylor Approximation of AC Power Flows: A Unit Commitment Study

Modern mixed-integer quadratic solvers generally handle binary variables more efficiently than nonlinear mixed-integer solvers. This is relevant to the power system operation models as the unit commitment formulations typically contain a large number of binary variables. This paper investigates how to achieve the accuracy level close to the one of the exact nonlinear models, but by utilising convex models and solvers. The presented unit commitment model is based on a Taylor-series expansion where both the voltage magnitude and angle are quadratically constrained. To achieve high accuracy, the model takes advantage of the meshed transmission network structure that enables replacement of the quadratic inequality constraints that cause constraint relaxation errors with the linear equality constraints. Quadratic constraints to be replaced as well as the operating point parameters are determined based on the presolve. The first presented case study validates the model's accuracy and the convergence of the iterative algorithm, while the second is a non-iterative full unit commitment problem. Unit commitment results show superior accuracy and similar computation times to the existing quadratic formulations on one hand and faster computation times than the exact nonlinear polar formulation on the other.


Sets and Indices N
Set of buses, indexed by i and j. N P Tuple set of paired buses aligned with branch E orientation, indexed by (i, j). R Set of reference buses, indexed by i. E, E R Tuple set of branches, forward and reverse orientation, indexed by (e, i, j). Branch π-section susceptances. τ e , σ e Branch tap magnitude and shift angle. P g k , P g k Generator minimum and maximum active power production.
Second-order Taylor series voltage magnitude term approximation.

Binary variables x t,k
Generator activity state indicator. y t,k , z t,k Generator start-up and shut-down indicators.

I. INTRODUCTION
A. Motivation U NIT commitment is an optimization problem that determines the least-cost production of the generators to satisfy the demand while considering the generators' physical limitations. Besides generator capacity limits, it typically encompasses output ramp limits, minimum up and down times and start-up costs, whose modeling requires computationally expensive binary variables. Computational burden of the unit commitment problem can be assessed from two perspectives: the binary formulation tractability and the grid formulation tractability.
Novel binary formulations are typically studied without the network constraints to better demonstrate and isolate the source of numerical difficulties. To this end, there is even a unit commitment benchmark generally accepted by the scientific community [1], but without the network constraints. However, the ISO-type electricity markets consider transmission constraints already at the market-clearing phase. Hence, we focus on the network-constrained unit commitment problem and deliver a new transmission system power flow formulation that retains an accuracy level close to the exact nonlinear AC models, but allows for the use of more specific and performant mixed-integer quadratically constrained quadratic program (MIQCQP) solvers as opposed to the mixed-integer nonlinear (MINLP) ones. It is a common practice in transmission system modeling to reduce a MINLP to a mixed-integer linear (MILP), e.g. by applying DC network approximation. However, linear models are inaccurate when it comes to modeling reactive power flows, voltage magnitudes and losses. On the other hand, the existing convex quadratic approximations or relaxations of AC power flows may achieve good or even perfect accuracy when there are no quadratic constraint relaxation errors. These are a consequence of the convexification process that requires the quadratic equality constraints to be relaxed into quadratic inequalities. However, when such errors do occur, they are very large. Since unit commitment is a multi-period optimization problem thus simulating the grid under various conditions, including generating units nonconvexities, the chances of having relaxation errors in at least one of the simulated periods are relatively high, rendering the relaxation model inaccurate. A motivation to develop a model based on the Taylor expansion comes from the thought that constraint relaxation errors, which in this case can occur due to convexification of the second-order parts of the expansion, can be avoided without otherwise gross changes to the power flow equations by upfront neglecting these terms on per-branch basis as determined by the presolve. Constraint relaxation errors due to convexification are avoided since when the second-order terms are neglected, the resulting constraints are linear equalities instead of quadratic inequalities. Additionally, meshed transmission network structure acts favourably as it provides a sufficient number of quadratic constraints to preserve the accuracy and iterative convergence of the algorithm that reruns the model around an updated Taylor operating point despite some neglected second-order terms. A parallel can be drawn to the Newton's power flow calculation method that converges faster when using higher-order modifications [2]. Fast convergence is important as it allows us to achieve high warm-start single iteration accuracy using the approximate warm-start operating point parameters obtained by first solving the problem in its continuous version where binary variables are relaxed into continuous in the range from 0 to 1. Utilising the described features and the proposed transmission system power flow formulation, we solve the network-constrained unit commitment problem while cooptimizing the real and reactive powers to exploit the unused monetary value in the traditional separate optimization, as demonstrated in [3], but by using an implementation without any loop statements.

B. Literature Review
Unit commitment research started with the development of the branch-and-bound algorithm [4], which is the basis of modern mixed-integer solvers. The MILP approach is still considered as the state-of-the-art due to its computational tractability. The early works grasped the unit commitment problem in its simple form without the network constraints, thus focusing on tightening the binary formulation, which can be expressed using three such variables per generator, as in [5], [6] and [7], or using a single binary variable as in [8]. Because of reliance on the branch-and-bound algorithm to solve mixed-integer problems, less variables does not necessarily imply better computational tractability. Size of the problem can be decreased by clustering similar generators [9], but this requires simplifications that reduce accuracy and applicability.
The subsequent unit commitment research branches out in multiple directions, mainly focusing on uncertainties [10] security constraints [11], and network constraints [12]. Security constraints add an additional contingency scenario dimension to the unit commitment problem drastically increasing the problem size. Our work and this literature review are focused on network constraints, whose inclusion in the unit commitment model also has a detrimental effect on the computational time. The inclusion of both security and network constraints forms an even more demanding problem, however in this work we focus on the network constraints to better isolate their difficulties and features. To reduce the problem complexity, DC optimal power flow is widely used. This transmission grid approximation results in good accuracy of the active power flows [13]. There are various attempts to generalize the DC model to include losses and reactive power, e.g., by expanding the first-order Taylor series around the operating point and adding only the voltageangle-dependent nonconvex piece-wise linear losses [14], which require integer variables, or by linear loss estimation [15], which needs a penalty factor to prevent negative losses. Quadratic approximations [16] are much more accurate, but also more computationally demanding than the linear models. Work [17] proposes a Taylor-based piece-wise linearization with no integer variables, i.e. linear programming approximation of AC power flows (LPAC) of the initially quadratic approximation, that considers angle-only dependent losses. In a post-publication in PowerModels package [18], the author implemented an enhanced variant of the model with better accuracy by using the quadratic losses constraint instead of the piece-wise linear one.
Relaxations, on the other hand, have less persistent accuracy that is highly dependent on the test case. Performance of Jabr's (Second-order Cone Programming -SOCP) [19], quadraticconvex (QC) [20] and Shor's (semi-definite programming -SDP) [21] relaxations were analysed in [22], which showed that Jabr's relaxation is dominated in terms of tightness by the both remaining formulations. Finally, there are nonconvex exact rectangular current-voltage [23] and voltage-based rectangular [24], [25] and polar [24], [26] formulations that can be directly utilised for unit commitment, but with the highest computational burden. An overview of the described optimal power flow (OPF) techniques is provided in Table I. Our work builds upon the quadratic implementation of the LPAC model by introducing the voltage magnitude-dependent losses and by expanding the Taylor series around a general operating point.

C. Paper Contribution and Structure
Contribution of the paper consists of the following: r We develop new transmission network AC equations based on the Taylor's expansion that approximates the secondorder voltage terms. The approximation consists of distributing power losses to both branch ends based on the forward-and reverse-orientation power flows.
r We develop a presolve technique for deciding whether to use the quadratic or the linear form of power flow constraints to avoid constraint relaxation errors due to convexification.
r The resulting MIQCQP solution is obtained much quicker than the MINLP solution without sacrificing accuracy.
Rest of the paper is structured as follows. Section II mathematically derives and states the proposed model. It is divided in three subsections: SubSection II-A presents the Taylor expansion analysis, SubSection II-B introduces the presolve technique and subSection II-C presents the model components. Case study Section III consists of the three main parts: the description and set-up III-A; and two case studies. In the first case study III-B we solve optimal power flow problem on a number of networks to demonstrate convergence and accuracy of the model. The second case study III-C solves unit commitment problems to demonstrate computational tractability and accuracy of the model. Section IV provides relevant conclusions and guidelines for future work.

A. Taylor Expansion Analysis
Our analysis of the optimal power flow starts from the polar formulation for the branch power flow in equations (1.1) and (1.2). For clarity of the analysis, a general branch, which encompasses lines and transformers according to the PowerModels 0.13 standard [18], is simplified to a line without the shunt sections, i.e. tap is 1∠0 • and g fr e , g to e , b fr e , b to e = 0. Also, to shorten the expressions, two substitutions are made: Otherwise, the final presented model (2.1)-(2.20) and all the benchmarked models include a general branch without any simplifications or substitutions.
Full second-order Taylor series for the branch active and reactive power flow is structurally written in expressions (1.3) and (1.4). The first line contains the zeroth-order part (in blue), the next three lines the first-order part (in green) and the last six lines the second-order part (in red) of the Taylor series expanded over variables V t,i , V t,j and θ t,i,j around the operating point parameters Variables representing a change from the operating point (delta variables) are at the beginning of the row, while the corresponding coefficients are within the square brackets.
Expressions (1.3) and (1.4) are nonconvex because i) they are equalities, and no quadratic equality is convex; ii) (V Δ t,j ) 2 term is multiplied by zero. Convexification of the polar AC OPF by Taylor series was studied in [17], where the model was obtained by ignoring all but one of the second-order terms, (θ Δ t,i,j ) 2 , and by taking θ op t,i,j = 0, which leaves the model voltage-wise loosely constrained, i.e. without the voltage second-order terms.
We performed extensive tests to determine the importance of different second-order terms by removing them one-by-one from the series in expressions (1.3) and (1.4) and analysing accuracy and iterative algorithm convergence which reruns the model around an operating point obtained from a previous solve. Since the Taylor series is nonconvex (quadratic equality), nonlinear IPOPT solver was used. As a result, all combinations without the (V Δ t,i ) 2 term failed to algorithmically fully converge and the version with only (V Δ t,i ) 2 and (θ Δ t,i,j ) 2 terms converged extremely slowly, requiring over 100 iterations. However, the versions with terms had little effect. Also, reactive power voltage second-order terms, which in the tests were shown to be better ignored than approximated, have small impact on convergence and overall accuracy due to an absence of reactive power in the objective function. Based on these results, we include the quadratic angle term (θ Δ t,i,j ) 2 convexified by a relaxation in constraint (1.5), as well as the voltage (V Δ t,i ) 2 and V Δ t,i · V Δ t,j terms for active power flow by convex approximation in (1.6). The relaxation in constraints (1.5) and (1.6) refers to the swap of the equality sign with the inequality. However, since the constraints are obtained by the Taylor expansion, which is an approximation, even with inequality sign they are still an approximation, albeit convex.
Constraints (1.5) and (1.6) are embedded in the main model as constraints (2.9.1) and (2.8.1), respectively, without simplifications and substitutions introduced at the beginning of this section. To retain similarity with the existing formulations in the literature, the right-hand side of (1.5) is a second-order Taylor series of a cosine function and, once multiplied with a parameter in the power flow constraints (2.4)-(2.7), it forms a part of the power flow zero-order Taylor term and (θ Δ t,i,j ) 2 term. The approximation in expression (1.6) is obtained by Taylor terms from both the forward and the reverse orientations of the branch active power flow, effectively representing losses evenly distributed between both branch ends in active power constraints (2.4) and (2.5). Second-order voltage approximation is always at least marginally convex, assuming g e > 0.
While the second-order voltage approximation from (1.6) could mathematically be applied to reactive power as well (by simply extracting and transferring factor g e from the approximation into the active power flow constraints and −b e for the reactive power flow constraints), it was determined that this approximation is in some cases inadequate for reactive power, as it leads more commonly to an infeasible model using the flat start operating point assumptions, i.e. V op t,i = 1 p.u. and θ op t,i = 0 rad. Reactive power flows are generally about an order of magnitude lower than active power flows, while the dominant factor for the approximation in reactive power flows, susceptance b e , is an order of magnitude higher than the dominant factor, conductance g e , in active power flows. Applying approximation only to active power flows does not limit its purpose to better constrain voltages since active power has much stronger effect on the objective function and thus on constraining the voltages. Also, no accuracy drawbacks were observed when applying the approximation only on active power flows in conditions where b e /g e was close to 1 and in lightly loaded networks. Constraint (2.8.1), the subsequent of constraint (1.6), is applied in the model as a constraint for every branch in the forward orientation E rather than bus-pairs N P , so it can account for parallel conductances (g fr e and g to e ) and differing taps in parallel branches.

B. Presolve Technique
An important problem feature is a meshed structure of the transmission networks. It makes all the delta variables V Δ t,i and θ Δ t,i quadratically constrained despite removing some quadratic inequality constraints and swapping them with the equality constraints since the delta variables appear in multiple instances of (2.8.1) and (2.9.1) constraints. Swapping the quadratic inequality constraints with linear equality constraints avoids constraint relaxation errors due to convexifications, i.e. errors due to a deviance from the inequality boundary, without a significant loss in a single iteration accuracy or multiple iterations convergence. Constraint (2.8.1) is a second-order voltage approximation without the simplifications introduced in the analysis, while (2.8.2) is its linear alternative. Similarly, (2.9.1) and (2.9.2) are quadratic and linear representations of the Taylor series cosine. Quadratic forms of the constraints are used only if the respective Boolean parameter Λ t,e or Γ t,i,j is true and if conductance g e is positive in the case of voltage approximation constraint (2.8.1).
Our proposition is to determine the value of the Boolean parameters, i.e. decide whether to use linear or quadratic forms of the voltage and cosine approximation constraints, in the presolve computation step. The AC unit commitment is a difficult problem mostly due to binary variables needed for generators in combination with network constraints, which slow down the solution process and limit the selection of solvers that can be used. However, if binary variables are fixed, i.e. replaced with a parameter, the problem is relatively easy to solve even if the network constraints are nonconvex or nonlinear. Furthermore, since the constraint marginal represents the sensitivity of the objective function on adding a small positive constant to the right-hand side of the constraint, sign of the equality constraint marginal, which is computed by default by, e.g., IPOPT and Knitro solvers, indicates if this constraint would be binding if it was relaxed into an inequality constraint. For constraint (2.8.1) to be binding, due to its greater-or-equal sign, q V t,e should have a tendency to be as small as possible. Adding a positive constant to the right-hand side of (2.8.1) would in this binding scenario increase q V t,e and thus worsen the objective function, i.e. the marginal would be positive. Oppositely, for constraint (2.9.1) to be binding, since its less-or-equal sign, its marginal needs to be negative. As a result, a nonconvex presolve with only quadratic approximation constraints (2.8.1) and (2.9.1) in the equality form and fixed binary variables is proposed. Values of the fixed variables are determined simultaneously with an approximate operating point as described in SubSection III-C. This way, the main unit commitment solve will have only the convex quadratic approximation constraints that the presolve flagged as binding. All other quadratic constraints are replaced with linear ones to avoid constraint relaxation errors.

C. Optimization Model
This section presents the whole network-constrained unit commitment model. The objective function (2.1) is a variable generation and start-up cost minimization. While the ISO markets typically consider piecewise linear cost curves, we tend to avoid further approximations and thus use the quadratic cost curves. Constraints (2.2) and (2.3) are the bus balance constraints, (2.4)-(2.7) are power flow equations which also contain second-order term approximation variables q V t,e and cos t,i,j from (2.8.1)-(2.9.2). Constraints (2.10) and (2.11) are generator production constraints that disable production when a generator is inactive (x t,k = 0). (2.12) is the branch apparent power constraint for both orientations, (2.13) is the reference bus angle constraint, (2.14) and (2.15) are voltage magnitude and bus-pair angle constraints, and (2.16) is the generator rampup and -down constraint. Constraints (2.17) and (2.18) model the interaction between the generator activity binary variable and the start-up and shut-down binary variables. Constraints (2.19) and (2.20) ensure generator minimum up and down time requirements. Such start-up and shut-down formulation using three binary variables was first presented in [6] and was proven efficient in [27]. The presented problem is of MIQCQP [28] class with convex objective function and constraints, except for binary variables.

A. Description and Set-Up
Although the main case study is elaborated in SubSection III-C, where we demonstrate the effectiveness of the proposed model, we start by applying the model on a number of OPF problems to show convergence speed and accuracy of the model. In all cases we compare the obtained results, both in terms of accuracy and computational effort, to the exact nonlinear polar model [24], i.e. based on constraints (1.1) and (1.2), which serves as a reference. Additionally, in the unit commitment case study, our results are compared with the LPAC approximation in the warm-start and the quadratic implementation [17], [18] and with the QC relaxation [20]. Our implementations of the existing Run nonconvex model 5: Select constraints by evaluating marginals 6: Run convex model results in Table II  7: Update operating point 8: until | * |gap < 0.005% gap to exact polar models were verified in a side-by-side comparison to match the PowerModels [18] implementation. In SubSection III-B all models were solved using IPOPT solver that has proven to be numerically highly robust and was set up to run with the HSL linear MA27 and scaling MC19 modules, option to always apply scaling and at most 500 solver iterations. In SubSection III-C the convex MIQCQP/MISOCP models were solved using Gurobi 9.0.2, while all the other models using Knitro 11.1. Both Gurobi and Knitro were run under the default settings on Intel i5 7600 CPU on 4 threads on a machine with 16 GB of RAM, while IPOPT was run on only one thread since it showed no performance scaling using additional threads. The solvers were instructed to compute to full optimality, unless the time limit of 1 h is reached. GAMS 31.1.1 was used as a modelling language and all continuous variables were declared as free variables. In the SubSection III-B, variable start values used to generate the model's Jacobian and Hessian matrices for the nonlinear solvers were reset to flat start assumptions after every computation. The specifics on bounds and variable starting values are important since they affect nonlinear solvers. All variable and parameter units are in p.u. or dimensionless.

B. Convergence and Accuracy Demonstration
This section demonstrates the accuracy of the model using OPF under assumption that the constraint type selection is known, i.e. determined by steps 4 and 5 in Algorithm 1. In an iterative procedure, the model is first run in its nonconvex form in which constraints (2.8.1) and (2.9.1) are equalities. Signs of their marginal values determine which constraints would deviate from the inequality boundary if they were relaxed. Subsequently, the model is rerun but in its convex form with selected quadratic or linear constraints to avoid constraint relaxation errors due to convexifications. The convex model updates the operating point voltage and angle parameters initially set to 1 p.u. and 0 rd. The procedure continues until fully converged in comparison to the exact nonlinear polar model [24], as described in Algorithm 1. The purpose of this case study is purely to explore the model's behavior on a large number of cases and not to solve an OPF as there are better existing options for problems with no integer variables, e.g. using the exact nonlinear models.
A well-established power grid model benchmark PGLib-OPF v19.05 [29] is used to perform a single time step optimization. Thus, the unit commitment and ramping constraints (2.16)-(2.20) as well as the unit commitment binary variables are  [29] removed. The resulting convex and nonconvex models minimize (2.1) subject to (2.2)-(2.15) over variables The nonconvex model has only constraints (2.8.1) and (2.9.1) in the equality form, while the convex one has both linear and relaxed quadratic forms of constraints determined by the nonconvex presolve.
Results of the described iterative procedure for the convex step 6 in Algorithm 1 are shown in Table II. The iteration gap is defined as a percentage deviation from the exact nonlinear polar solution. Other metrics used are the computation time, the number of linear voltage constraints (2.8.2) with positive g e , the number of linear cosine constraints (2.9.2) and the number of quadratic constraints (2.8.1) and (2.9.1) that deviated from the inequality boundary. Total number of voltage linear constraints is equal to the number of branches with nonpositive g e plus those decided to be replaced by the iteration's presolve. For a better overview of the grids' sizes, there are also columns with number of branches and bus-pairs, while the grid names contain the number of buses.
The model converges fully on all grids within three iterations from the first feasible iteration step. Full convergence results in 0.00% approximation error since the Taylor expansion is exact at the expansion point and there are no relaxation errors related to constraints (2.8.1) and (2.9.1), as displayed in Table II. This is because the presolve removes those relaxed quadratic constraints that would deviate from the inequality boundary with their linear equality alternatives. However, the flat start operating point assumption at the first iteration is not good enough for the model to provide a feasible solution in all cases.
The solver usually does not recognise infeasible models, but reaches iteration limit marked by "it/inf" in the gap column. Despite being infeasible, a solution is returned by the solver and used to update the operating point for the next iteration. The only two grids that did not converge within four iterations are 6468_rte and 6495_rte. Very similar grids, 6470_rte and 6515_rte, did converge because they got updated with more favorable operating points quicker than the former two grids. It took seven iterations for 6468_rte and 6495_rte grids to fully converge as their first feasible solution was reached in the fifth iteration. Fast convergence indicates that very high accuracy can be achieved if the model is warm-started with a reasonably good operating point and constraints to be linearized selected by the presolve, which is pivotal for SubSection III-C.
The model's working principle can be well described on the 3_lmbd grid, which is small but congested. The second-order voltage approximation variable seldom deviates from the inequality boundary in (2.8.1) since such deviation would only increase active power losses at both branch ends. Cosine variable deviation from its bound in quadratic constraint (2.9.1) causes much higher reactive power losses than the active power losses due to branch susceptance being much higher than conductance. In 3_lmbd grid the ratio b e /g e is −30 for the branch with (2.9.2) applied. That branch in the optimal solution is congested and also producing reactive power. Thus, reactive power losses are favorable as they enable greater transfer of active power. To prevent false losses by the cosine approximation constraint relaxation, the presolve determined that the linear equality constraint should be applied. However, congestion is not a necessary condition for potential cosine That occurs more commonly with high negative b e /g e ratios.
Iterations that were infeasible or reached the iteration limit have an unrealistic number of linear constraints due to equivocal information from the proposed nonconvex presolve that is also infeasible. It is important that the network data contains entries for branch conductances as otherwise the model would contain only linear voltage constraints (2.8.2), making it insufficiently quadratically constrained for convergence and warm-start accuracy. This case study demonstrates, however, that even with partial conductance data, the model operates well as grid 179_goc has 27% of all branches with null conductances.

C. Unit Commitment
We use a 24-hour network-constrained unit commitment to evaluate accuracy and computational performance of the model. Due to availability of a limited amount of unit commitment network data, grids from the OPF benchmark [29] are adapted for unit commitment purposes using generic ramp limits, start-up costs, minimum up-and down-times, load curve and generator costs. Ramp limits are set so the mean of the allowed production range can be reached within a single time period, i.e. (|P g k | + |P g k |)/2, start-up costs are set to 1500 cost units for all generators, minimum up-and down-times are set to 2 hours for generators with ≤1 p.u. (100 MW) and 4 hours for larger ones, generator fixed costs c k are modified to 10% of the linear costṡ c k if their original value is zero. Loads use scaling factors for the winter weekday periods from IEEE RTS-96 [30].
The proposed quadratic approximation model is designed to operate well in the vicinity of the operating point parameters V op t,i and θ op t,i . Thus, the first step is to provide a good operating point by solving the continuous exact nonlinear model with relaxed binary variables. The next step is to run the nonconvex continuous quadratic presolve, i.e. with (2.8.1) and (2.9.1) as equalities, using fixed (replaced with a parameter) binary variables and around the operating point computed in the previous solve. The fixed binary variables inherit relaxed, i.e. continuous, values. Its solution is the same as in the first step, but it returns the constraints' marginal values whose sign determines the constraint type (linear or quadratic) for every branch and bus pair of the main unit commitment solve. The third step is to run the mixed-integer unit commitment. The last, fourth step is to determine the approximation error by running the exact continuous polar formulation [24], but with binary variables fixed to the solution of the mixed-integer unit commitment run from the previous step. Inexact models may return binary variable values for which no solution is possible. In that case the approximation error is marked with "inf" in Table III. Compared to the existing LPAC, the proposed model has an additional, but easy to solve, second step to select constraints to be linearized. The described procedure is itemized in Algorithm 2. Optimization source code is provided on GitHub [31].
The simulation results are provided in Table III. Convex nature of the presented model allows for the use of traditional solvers typically used to solve MILP problems as they handle binary variables more efficiently than MINLP solvers. Thus, the presented model outperforms the exact polar MINLP in all but one test case, with computation times up to 77 times faster for 39_epri grid. Test cases 24_ieee_rts and 73_ieee_rts are difficult due to a large number of generators resulting in a large number of binary variables, which severely impacts computational tractability. It was found that in these two test cases the other solvers outperformed Gurobi. For the proposed MIQCQP model applied to 24_ieee_rts, CPLEX 12.10 reached 0.09% MIP gap as compared to 1.26% achieved by Gurobi. The best performance on 73_ieee_rts was achieved with Xpress 8.8, which reached 2.78% MIP gap, while Gurobi did not find a feasible solution. Furthermore, this MIP gap is much lower than the 10.29% gap obtained using Knitro 11.1 with the exact polar MINLP model. Both grids at the solution point had 0.00% objective function approximation errors. Test cases 162_ieee_dtc and 179_goc were infeasible even for the QC relaxation and were eliminated from the study, as well as networks larger than 200_tamu as they are too difficult to solve.
The most significant result of the study is very low, almost nonexistent, approximation error. For 13 test cases, the highest objective function error is −0.02%. The presented model is very accurate around a broad vicinity of the operating point, which can be sufficiently well estimated by solving the model with relaxed binary variables. Furthermore, at the operating point the quadratic constraints that can cause large errors by deviating from the inequality bounds are replaced with the linear ones. Quadratic constraints that finally do deviate from the inequality bound in the unit commitment solve cause a very small relaxation error. If this error was larger, quadratic constraints would deviate at the operating point and thus be replaced with linear ones. Furthermore, we evaluate the approximation errors of the individual variables using average and maximum normalized distance of solution variables [32]. Normalization is carried out by dividing the absolute variable error value by the variable feasible range. For P t,e,i,j and Q t,e,i,j , the feasible range is 2S e , for V t,i it is V i − V i and for the branch angle difference θ t,i,j it is θ i,j − θ i,j . The results are displayed in Table IV. The highest average active power flow variable distance is 0.06%, reactive power flow distance 0.14%, voltage distance 1.42% and branch  angle distance 0.02%, which is significantly lower than what can be expected from the competing QC and even Shor's model that normally achieve overall average distances in the range from 5 to 10% (see Fig. 3 in [32]). The approximation errors can be even further reduced by iteratively running the unit commitment solve by updating the operating point and retesting the constraints for that new operating point. However, as an alternative, due to a low objective function approximation error, we recommend using the verified solution, i.e. Algorithm's 2 step 4, as the final solution to the unit commitment problem.
Achieving great accuracy in unit commitment is important for two reasons: a) active power losses are only about 2% of the total production, thus the accuracy needs to be far greater than 2% to account for them properly; b) inaccurate solutions can result in infeasible decisions for generator activity status, i.e. binary variables. The LPAC model in its more accurate quadratic variant, despite being warm-started similarly like the proposed model, has the highest approximation errors −0.72%, −0.39% and −0.33% and three out of 13 cases are infeasible. QC relaxation, despite proven to be tighter than the Jabr's relaxation, is still overly optimistic and thus results in 7 infeasible cases, while the greatest error of feasible models is −0.56%. These results highlight the importance of quadratically constraining the voltage in combination with the constraint type selection procedure. Computation-time-wise, the proposed MIQCQP model is in between the faster LPAC approximation and the slower QC approximation.

IV. CONCLUSION
The presented case studies demonstrate that the proposed MIQCQP model is moderately accurate assuming flat start operating point parameters, but almost perfectly accurate if the assumed operating point is in a broader vicinity of the optimal solution. Accuracy is achieved by having both the voltage magnitude and the angle quadratically constrained and by replacing some of the quadratic inequality constraints with linear equality constraints to avoid constraint relaxation errors due to convexification as determined by the presolve. Meshed transmission network structure acts favourably to further constrain the model so it can withstand the replacement of some quadratic constraints with linear ones and stay accurate and convergent.
Since the model's power-flow-related constraints are convex quadratic, the model supports solvers that handle binary variables more efficiently than the MINLP solvers. The proposed model displays average computational tractability for the MIQCQP/MISOCP problem class, which means it is slightly slower than the quadratic implementation of the LPAC approximation [17], [18] and quicker than the QC relaxation [20]. Accuracy of the existing quadratic models is insufficient for network-constrained unit commitment leading to frequent infeasible decisions for binary variables or significant approximation or relaxation errors.