# Actionable Three-Phase Infeasibility Optimization with Varying Slack Sources
††thanks:  * denotes the corresponding author

Elizabeth Foster1∗, Timothy McNamara1, Amritanshu Pandey1,2 and Larry Pileggi1
{emfoster, tmcnama2, amritanp<EMAIL_ADDRESS>1 Electrical &
Computer Engineering, Carnegie Mellon University, Pittsburgh, PA  2 Electrical
& Computer Engineering, University of Vermont, Burlington, VT

###### Abstract

Modern distribution grids that include numerous distributed energy resources
(DERs) and battery electric vehicles (BEVs) will require simulation and
optimization methods that can capture behavior under infeasible operating
scenarios to assess reliability. A three-phase infeasibility analysis is based
on localizing and identifying power deficient areas in distribution feeders
via a non-convex optimization that injects and subsequently minimizes slack
sources, subject to AC network constraints. In our power flow optimization, we
introduce and compare three different representations for slack resources:
current, power, and admittance. Operational bounds are incorporated to ensure
realistic solutions for all cases. We find that the slack power formulation
converges fastest across our testcases of up to 1599 nodes, and both slack
power and slack admittance generally return smaller net power requirements for
network feasibility as compared with slack current models. Formulations with
slack power models and slack admittance models are demonstrated to effectively
represent reactive power compensation on a taxonomy feeder.

###### Index Terms:

distribution grid optimization, three-phase infeasibility analysis, voltage
bounds, optimal power flow

##  I Introduction

Modern electric grids are seeing an influx of distributed energy resources
(DERs) and smart technologies like battery electric vehicles (BEVs) at the
distribution level [1] [2] [3]. These new resources and technologies present
opportunities to achieve climate goals, amongst other benefits; however, they
also present new challenges as system planners attempt to build out the
infrastructure required to incorporate these assets while maintaining grid
reliability and resiliency [2] [3].

System planners use tools that are generally limited to traditional methods
like three-phase power flow, which outputs nodal voltages and line/transformer
flows under different feasible operating conditions [3]. Three-phase power
flow will struggle to adequately assess the challenges presented by BEV and
DER growth because scenarios requiring upgrades might fail to converge due to
an infeasible network with unsatisfied AC network constraints. For example,
simulation alone cannot identify optimal locations for needed infrastructure
upgrades that will increase electric load, like fast charging stations, or
areas of weakness due to insufficient reactive power compensation in the
network. Tools for distribution grid power flow offer limited, if any,
optimization capabilities [3]. New distribution planning tools and
capabilities are needed that are capable of running optimizations while
satisfying three-phase AC unbalanced network constraints. These new tools will
help system planners make optimal decisions on updating and installing new
infrastructure when combined with operational, financial, and policy
constraints [4].

While availability of such complete three-phase power flow tools are limited,
three-phase optimization is an active and broad research area [5]. Several
optimization approaches are targeted to three-phase optimal power flow (OPF),
which typically minimizes economic dispatch subject to network constraints.
These OPF methods have often been formulated using power balance network
equations as seen in [6], which applies OPF to active distribution networks.
The authors in [7] extend traditional power-balance constrained OPF to include
slack variables that indicate constraint violations when there is no feasible
solution. Other three-phase optimizations include technical objectives, like
those that minimize losses [8], while still being subject to network equality
constraints.

Formulations using power-balance equations [6] [7] [8] are highly nonlinear
and therefore do not scale well in general. This behavior is shown in [9],
which developed a four-wire OPF with both current-voltage (I-V) and power
balance formulations and found that power-balance did not converge for 35% of
testcases. Not only was I-V faster across the moderately sized testcases where
power balance did converge, but the I-V formulation also returned more
practical solutions. In the power balance formulation for 4-wire systems, some
of the nodes were non-physically grounded, thereby resulting in erroneous
voltage solutions [9]. Furthermore, many of these three-phase network-
constrained optimizations will fail when there is no feasible power flow
solution if slack variables or constraint violations are not considered.

A three-phase infeasibility analysis (TPIA) is another type of three-phase
optimization formulated subject to current-based network constraints for use
in distribution grid planning. TPIA localizes and identifies deficient areas,
in addition to the amount of the deficiency, in situations where power flow
fails [10]. At a high level, it is a non-convex optimization study that
identifies scenarios where AC three-phase power flow constraints are
infeasible by using a set of variable slack sources with a physical meaning to
compensate for mismatch in the nodal equations. In TPIA, solutions with
nonzero slack sources occur when power flow would otherwise fail.

Initially, the concept of quantifying infeasibility used power balance network
constraints for positive sequence systems to find a solvable boundary for
otherwise insolvable power flow simulations [11]. We formulated an
infeasibility analysis using slack current sources and KCL network constraints
for transmission networks in [12] and three-phase networks in [10], where we
additionally developed an L1-norm objective formulation that effectively
localized the sources of infeasibility. Our work in [10] and [12] focused only
on slack current models that can be directly added into the KCL network
constraints. Neither work included any analysis on whether the algorithmic
performance would change if slack current sources were replaced with either
slack power or slack admittance sources, both of which can represent physical
grid assets. In addition, operational constraints, like nodal voltages bounds,
were not included in the prior formulation.

In this paper, we evaluate slack model choices by analyzing both problem
complexity and real-world applications. We ensure the solutions from our TPIA
formulations are feasible and actionable by using inequality constraints to
incorporate bus voltage magnitude restrictions. We give insight into the
mathematical complexity of each formulation in Section III and empirically
verify computational complexity using convergence properties from infeasible
networks in Section V. To evaluate the practicality of each slack source, we
consider the relationship between physical assets that address real-world
challenges, like capacitor banks, and the information provided from each TPIA
formulation. We find that slack power or slack admittance best represent
reactive power compensation because real power or conductance slack variables
can be set to zero, which allows only reactive power slack injections.

##  II Background

The background for this work is primarily the concept of infeasibility
analysis and the circuit heuristics that are utilized for three-phase
optimization.

###  II-A Infeasibility Analyses

Infeasibility analysis was first introduced for transmission networks in [11],
and later in terms of an equivalent circuit for power grids in [12]. [12]
introduced slack current sources at each node in the positive sequence circuit
model of the power grid. It minimized the norm of the slack current sources to
identify when positive sequence power flow was infeasible and returned the
amount of missing current at each node in the system. [13] extended [12] to
localize the sources of infeasibility by using an L1 regularization term in a
least-squares objective. In [10], we introduced a scalable TPIA framework
using three-phase current-based network constraints with two different
objective functions; one using a least-squares objective and the other an
L1-norm objective, which localized the solution to a subset of nodes.

However, the equivalent circuit approaches in [10] [12] [13] cannot directly
relate slack source values with realizable physical assets on the grid.
Furthermore, they do not enforce physical bounds on the grid voltages, thereby
producing results with non-actionable infeasible solutions. Additionally,
minimizing current produces solutions where the net power injection might be
impractical given the varying base voltages across the distribution grid
(i.e., providing a 1 amp current injection at a 13.2 kV bus requires much more
power than supplying 1 amp at a 240 V bus). In this paper, we address these
limitations by: (a) formulating TPIA with alternative slack sources (slack
admittance and power); (b) restructuring the objective function to minimize
net complex power for all slack formulations; and (c) enforcing operational
bounds using inequality constraints.

###  II-B Circuit Heuristics for Three-phase Optimization

We leverage and adapt some of techniques that have been used for circuit
simulation to improve the performance of our three-phase power grid
optimization. One such technique is homotopy, which we have adapted for power
systems to attain convergence for difficult power flow and optimization
problems. To mitigate the convergence challenges due to a lack of good initial
guess, Tx-stepping (a power system specific homotopy method) begins by
relaxing the network constraints by adding large admittances in parallel to
all network branches and transformers, which effectively shorts most of the
buses in the network. An analysis of the relaxed network is first solved and
then the synthetic large admittance is gradually removed in small steps. The
analysis is solved repeatedly using the prior solution as an initial guess for
subsequent steps until the original complex problem is solved [14]. We will
apply Tx-stepping to solve TPIA in this work. To further ensure robust
convergence of TPIA when there are inequality constraints from operational
limits, we further leverage a limiting technique for Newton-Raphson (NR) that
is based on our understanding of individual device model characteristics.
First developed for diodes, this heuristic dampens the update steps of primal
and dual variables associated with inequality constraints so that their
respective values never enter infeasible regions [15].

##  III Methodology

Next we describe three versions of the TPIA formulation using different slack
sources: TPIA-I with slack current sources; TPIA-PQ with slack power sources;
and TPIA-GB with slack admittance sources. Our prior work in [10] covered a
version of TPIA-I, but for completeness, we reproduce it here. To ensure that
various TPIA formulations provide meaningful solutions, we introduce
operational limits in Section III-F.

In TPIA, if the network is feasible, the slack sources converge to zero.
Otherwise, nonzero slack sources occur and can be used to identify and
localize grid infeasibility for planning insight or corrections during
operation. In this paper, we consider various slack formulations to compare
their solution complexity and practical applicability in grid operation and
planning. All three model formulations provide flexibility to impose upper and
lower limits on the slack resources (e.g., constraining a susceptance $B$ per
phase $\phi$ by its upper rating: $B^{\phi}_{s}\leq\overline{B^{\phi}_{s}}$
[16]).

###  III-A Three-Phase Power Flow Constraints

We use nonlinear KCL network equations to represent the steady-state grid
physics in the TPIA formulations. Due to the modular circuit-based approach,
the TPIA framework can model any three-wire or four-wire distribution feeders
with various transformer configurations and distribution assets like inverter
based distributed generators, voltage regulators, triplex loads, capacitors,
fuses, shunts, and switches. The models for these devices are detailed in
[17].

We describe the KCL-based three-phase AC network constraints in (1) and (2),
where $f^{\phi}_{r}(x)$ and $f^{\phi}_{i}(x)$, are the respective rectangular
real and imaginary nodal current equations. $x$ is a vector of power flow
variables that includes the real and imaginary voltages, $V^{\phi}_{R}$ and
$V^{\phi}_{I}$, per phase $\phi$. $G$ and $B$ are matrices of network
conductances and susceptances. $I^{NL}_{R}$ and $I^{NL}_{I}$ are nonlinear
currents from either PV or PQ buses, but for brevity, we represent them as
though they are from PQ buses and omit the reactive power variable and voltage
setpoint equation for PV buses. Given a solution $x$ for a feasible network,
$f^{\phi}_{r}(x)$ and $f^{\phi}_{i}(x)$ are equal to 0, but in infeasible
networks, there will be a nonzero mismatch at some of the nodal equations.

$\displaystyle
f_{r}^{\phi}(x)=GV^{\phi}_{R}-BV^{\phi}_{I}+I^{NL}_{R}(V^{\phi}_{R},$
$\displaystyle V^{\phi}_{I})=0$ (1) $\displaystyle
f_{i}^{\phi}(x)=BV^{\phi}_{R}+GV^{\phi}_{I}+I^{NL}_{I}(V^{\phi}_{R},$
$\displaystyle V^{\phi}_{I})=0$ (2)

The vectorized nonlinear current equations $I^{NL}_{R}$ (3) and $I^{NL}_{I}$
(4) are functions of power injections. We later define slack power sources
with a modified version of (3) and (4).

$\displaystyle I^{NL}_{R}(V^{\phi}_{R},V^{\phi}_{I})=(P^{\phi}\odot$
$\displaystyle V_{R}^{\phi}+Q^{\phi}\odot V_{I}^{\phi})$
$\displaystyle\oslash((V_{R}^{\phi})^{2}+(V_{I}^{\phi})^{2})$ (3)
$\displaystyle I^{NL}_{I}(V^{\phi}_{R},V^{\phi}_{I})=(P^{\phi}\odot$
$\displaystyle V_{I}^{\phi}-Q^{\phi}\odot V_{R}^{\phi})$
$\displaystyle\oslash((V_{R}^{\phi})^{2}+(V_{I}^{\phi})^{2})$ (4)

###  III-B General TPIA Format

The general format for TPIA stays the same across the different slack variable
formulations. The objective function in 5a represents the sum of squares of
the complex power injection, $g_{m}^{\phi}(s)$, by the set of slack variables
$s$ for all $n$ buses and all connected phases $\phi$. We select this
objective function to address some of the challenges discussed in Section
II-A. The variable $\alpha$ represents a scaling parameter necessary to ensure
robust convergence in the presence of inequality constraints. Each formulation
is subject to the network constraints modified by $h^{\phi}_{r}(s)$ and
$h^{\phi}_{i}(s)$, which represent the real (5b) and imaginary (5c) current
injections from each slack source.

$\displaystyle\min\quad$
$\displaystyle\sum_{m=1}^{n}\sum_{\phi\in\\{A,B,C\\}}\alpha g^{\phi}_{m}(s)$
(5a) s.t. $\displaystyle f_{r}^{\phi}(x)-h_{r}^{\phi}(s)=0$ (5b)
$\displaystyle f_{i}^{\phi}(x)-h_{i}^{\phi}(s)=0$ (5c)

The functions for $g_{m}^{\phi}(s),h^{\phi}_{r}(s)$, and $h^{\phi}_{i}(s)$
will change depending on the choice of slack source. This is to ensure that we
compare equivalent quantities in the objective.

Although it can be shown that the subsequent three formulations for TPIA-I,
TPIA-PQ, and TPIA-GB are mathematically equivalent, the reformulations can
affect the nonlinearities in the objective function and the constraints, and
with different initial conditions, the choice of formulation can result in a
different solution trajectory.

###  III-C TPIA-I Formulation

In TPIA-I, we use real and imaginary slack current sources, $i_{fR}^{\phi}$
and $i_{fI}^{\phi}$, to develop insight into how much current (or post-
processed power) is needed at a given location. The objective function sums
$g^{\phi}_{m}(s)$, which is the square of complex power at node $m$ and phase
$\phi$ (6a). Complex power is calculated by the product of complex voltage and
conjugate complex current. The current injections in (6b) and (6c) are simply
defined by slack current sources themselves.

$\displaystyle g_{m}^{\phi}(s)$
$\displaystyle=\left((i_{fR,m}^{\phi})^{2}+(i_{fI,m}^{\phi})^{2}\right)\left((V_{R,m}^{\phi})^{2}+(V_{I,m}^{\phi})^{2}\right)$
(6a) $\displaystyle h_{r}^{\phi}(x)$ $\displaystyle=i_{fR}^{\phi}$ (6b)
$\displaystyle h_{i}^{\phi}(x)$ $\displaystyle=i_{fI}^{\phi}$ (6c)

Overall, the slack current sources add the least complexity to the network
constraints by introducing only linear terms but greater nonlinear complexity
to the objective because of the product of squared variables. It is important
to note, however, that TPIA-I cannot restrict injections to only reactive
power without introducing additional nonlinear constraints that set real power
to zero.

###  III-D TPIA-PQ Formulation

TPIA-PQ can provide the same insight as TPIA-I but has expanded use-cases when
considering real-world assets that provide only reactive power compensation
like static synchronous compensators (STATCOMs). The objective function sums
$g_{m}^{\phi}(s)$, which is the individual squares of the real and reactive
slack power sources, $P_{s}$ and $Q_{s}$, at bus $m$ and phase $\phi$ (7a).
The KCL current injections from the slack power sources, $I^{s}_{R}$ (7b) and
$I^{s}_{I}$ (7c), are equivalent to (3) and (4) except the original constants
$P$ and $Q$ are now variables $P_{s}$ and $Q_{s}$.

$\displaystyle g_{m}^{\phi}(s)$
$\displaystyle=(P_{s,m}^{\phi})^{2}+(Q_{s,m}^{\phi})^{2}$ (7a) $\displaystyle
h_{r}^{\phi}(x)$
$\displaystyle=I^{s}_{R}(V^{\phi}_{R},V^{\phi}_{I},P_{s}^{\phi},Q_{s}^{\phi})$
(7b) $\displaystyle h_{i}^{\phi}(x)$
$\displaystyle=I^{s}_{I}(V^{\phi}_{R},V^{\phi}_{I},P_{s}^{\phi},Q_{s}^{\phi})$
(7c)

Slack power sources introduce the most complexity to the KCL equations since
they are rational functions of voltages with a sum of squares in the
denominator. However, the objective function is only quadratic, which makes it
the least complex of the three formulations. TPIA-PQ can model only reactive
power injections by simply setting $P_{s,m}=0$ at each bus $m$. This
formulation, which we call TPIA-Q, requires no additional constraints.

###  III-E TPIA-GB Formulation

Since system planners might want to explore installing capacitor banks for
reactive power compensation, we develop an admittance-based TPIA formulation.
In TPIA-GB, we use slack conductance $G_{s}$ and slack susceptance $B_{s}$
variables to represent the slack admittance sources. The objective function
sums $g_{m}^{\phi}(s)$, which is the product of admittance and squared
voltages at bus $m$ and phase $\phi$ (8a). The current injections from the
slack admittance sources are found by multiplying voltage and admittance in
(8b) and (8c).

$\displaystyle g_{m}^{\phi}(s)$
$\displaystyle=\left((G_{s,m}^{\phi})^{2}+(B_{s,m}^{\phi})^{2}\right)\left((V^{\phi}_{R,m})^{2}+(V^{\phi}_{I,m})^{2}\right)^{2}$
(8a) $\displaystyle h_{r}^{\phi}(x)$ $\displaystyle=-G_{s}^{\phi}\odot
V^{\phi}_{R}+B_{s}^{\phi}\odot V^{\phi}_{I}$ (8b) $\displaystyle
h_{i}^{\phi}(x)$ $\displaystyle=-G_{s}^{\phi}\odot
V^{\phi}_{I}-B_{s}^{\phi}\odot V^{\phi}_{R}$ (8c)

The terms added to the network constraints by slack admittance sources are
bilinear, which is less complex than terms added by TPIA-PQ, but more complex
than TPIA-I. However, (8c) is the most nonlinear objective of the three
formulations. While TPIA-GB includes both $G_{s}$ and $B_{s}$, we recognize
that a formulation including only $B_{s}$ is more realistic since it
represents capacitive behavior. To allow only reactive power compensation, we
develop TPIA-B, which requires no additional constraints and sets $G_{s,m}=0$
on all buses.

###  III-F Operational Limits

Formulations in Sections III-C through III-E may help identify infeasible
networks, but their output may not be usable by grid planners. In some
particularly stressed networks, the bus voltage magnitudes in a TPIA solution
may lie outside of the acceptable range, rendering any information from TPIA
frameworks non-actionable. Therefore, to ensure that bus voltage magnitudes
are within specific ranges, we include inequality constraints on the bus
voltage magnitudes for all $n$ buses and phases $\phi$. In (9),
$\overline{V^{\phi}_{i}}$ and $\underline{V^{\phi}_{i}}$ are the constant-
valued lower and upper voltage magnitude bounds:

$\displaystyle(\overline{V^{\phi}_{i}})^{2}\leq(V_{R,i}^{\phi})^{2}+(V_{I,i}^{\phi})^{2}\leq(\underline{V^{\phi}_{i}})^{2},\quad\forall
i\in n$ (9)

(9) introduces additional nonlinearities to TPIA, so heuristics from Section
II-B are necessary for robust convergence.

###  III-G Optimization Solution Quality

For each of the above formulations, we solve the optimization problem using
the primal-dual interior point (PDIP) approach [18] aided by the circuit-
simulation heuristics from Section II-B. We seek a local minimizer by solving
for a set of primal and dual variables that satisfy the perturbed first-order
KKT conditions of the Lagrangian. The underlying problem is non-convex due to
nonlinear network constraints, so no standard approach exists to obtain the
global optima. However, our prior works in [17] and [19] have proposed
heuristics that can be used to improve convergence to practical and meaningful
solutions and to avoid saddle points.

##  IV Experimental Set-up

To evaluate the computational efficiency of the TPIA formulations, we
incorporated the new models and techniques into our three-phase power grid
implementation in [20]. Specifically, we implemented TPIA-PQ, TPIA-GB, and
TPIA-I formulations with a common objective and enforced voltage bounds. We
tested our approaches on the PNNL taxonomy feeders, where modifications are
made to some of them to increase the load factor to stress the networks and
induce infeasibility [21]. The testcases are available at [22] and described
in [10]. They represent a small urban center on the U.S. West Coast with 76
nodes (R1-12.47-3_OV); an urban area in the U.S. Northeast with 800 nodes
(R2-25.00-1_OV); and an urban area with rural spur in the U.S. Southeast with
1599 nodes (R4-12.47-1), which is the original version of the testcase and
exhibits solutions with voltages outside of operational ranges.

For consistency, we include slack variables at all nodes in the system.
However, in TPIA, slack variables can be included at only a subset of system
nodes based on the planner’s or operator’s knowledge. To demonstrate our
operational bounds, we used $0.9$ per unit (p.u.) of the nominal bus voltage
for $\overline{V^{\phi}_{i}}$ and $1.1$ p.u for $\underline{V^{\phi}_{i}}$ for
all results except in the reactive power use case where voltage bounds are
more strictly defined. Diode limiting is used for all voltage-bounded
formulations. Tx-stepping is only used when explicitly mentioned.

##  V Results & Discussion

Here we show the results and discuss performance and practical considerations
for the TPIA formulations.

###  V-A Impact of Voltage Bounds

Fig. 1 demonstrates the impact of voltage bounding on the solution for
R2-25.00-1_OV. In this case, 13 buses have at least one phase with a voltage
magnitude below 10% of its nominal value, as shown by the blue dots. When
voltage limits are included, all bus voltage magnitudes (red dots) are kept
within $\pm$10% by introducing infeasibility injections of 18% and 35% of the
total P and Q consumed by the feeder, respectively.

Figure 1: The impact of TPIA-PQ injections on satisfying bus voltage magnitude
limits for R2-25.00-1_OV.

Since omitting voltage limits may lead to non-actionable solutions, the
following results include voltage limits for all cases. Fig. 2 shows the
number of iterations it took for TPIA-I, TPIA-PQ, and TPIA-GB to converge for
the three infeasible networks. TPIA-PQ took the fewest iterations.

Figure 2: NR iterations in our PDIP solver [20]for convergence of each TPIA
formulation and testcase. Hatched filling implies that Tx-stepping was needed.

The formulations also yielded varying performances regarding the amount of
power required to make the network feasible. Table I reports the sum of
absolute values of the complex power from infeasibility sources dispatched
across all cases and formulations. While the optimization sought to minimize
the sum of squares of complex power injections, we report the simple sum
because those values are more intuitive for grid planning. For these cases,
TPIA-PQ and TPIA-GB mostly generate solutions with less slack power dispatched
than TPIA-I.

TABLE I: Total complex power injections | TPIA-I | TPIA-PQ | TPIA-GB
---|---|---|---
Case |  $\Sigma S_{f}$ (kVA) |  $\Sigma S_{f}$ (kVA) |  $\Sigma S_{f}$ (kVA)
R1-12.47-3_OV | 6.78e3+j4.65e3 | 2.03e2+j2.39e2 | 2.03e2+j4.34e2
R2-25.00-1_OV | 1.52e4+j1.03e4 | 5.81e3+j1.07e4 | 5.69e3+j1.26e4
R4-12.47-1 | 8.90e2+j8.40e2 | 1.37e3+j1.81e2 | 1.16e2+j1.93e2

###  V-B Reactive Power Compensation Use-Case

TPIA can help grid planners with optimal reactive power compensation. To
demonstrate its utility, we used the testcase R4-12.47-1, which, without
reactive power compensation, has 12 system buses (out of 1599 nodes) outside
of 0.95 and 1.05 p.u. voltage range, a range commonly enforced by real-world
utilities. To enforce the voltages with reactive power compensation, we ran
TPIA-B and found that the network needed 5 MVAR of reactive power to achieve
acceptable voltages. We then modified the testcase to include 40 capacitors
based on values from the TPIA-B solution and placed them at their respective
locations. We subsequently ran the modified testcase in TPIA-B without voltage
bounds, where it converged in 4 iterations and had all bus voltages between
0.95 and 1.05 p.u. without needing any infeasibility injections (which is
indicative of a power flow solution).

##  VI Conclusion

Three-phase analysis tools used for distribution grid planning must evolve to
account for the growing complexity at the grid edge. We developed a three-
phase infeasibility analysis with current, power, and admittance slack models
that can provide insight on where infeasibilities occur and the extent of the
infeasibility. To produce realistic results, we construct voltage bounds for
each formulation. We find that the formulations with slack reactive power or
admittance can be modified to represent specific reactive power assets like
STATCOMs and capacitors, respectively. Across the testcases that we studied,
the voltage-bounded slack power formulation converged the fastest and both
voltage-bounded slack power and admittance formulations generally returned
smaller net complex power injections than the voltage-bounded slack current
formulation.

## References

  * [1] International Energy Agency, “Global EV Outlook 2021”, April 2021. Available: https://www.iea.org/reports/global-ev-outlook-2021.
  * [2] National Academy of Sciences, Engineering, and Medicine. The Future of Electric Power in the United States. Washington, D.C: the National Academies Press, 2021.
  * [3] J. S. Homer, et. al. “Electric Distribution System Planning with DERs- High-level Assessment of Tools and Methods”, Pacific Northwest National Lab, Richland, WA, March 2020.
  * [4] “Reimagining the Grid”. Southern California Edison, Rosemead, CA, USA. Dec. 2020 [Online]. Available: https://www.edison.com/home/our-perspective/reimagining-the-grid.html
  * [5] A. Keane et al., “State-of-the-Art Techniques and Challenges Ahead for Distributed Generation Planning and Optimization,” IEEE Trans. on Power Sys., vol. 28, no. 2, May 2013.
  * [6] S. Gill, I. Kockar, and G. Ault, “Dynamic Optimal Power flow for Active Distribution Networks,” IEEE Trans. on Power Sys., vol. 29, no. 1, Jan. 2014.
  * [7] A. S. Zamzam, N. D. Sidiropoulos, and E. Dall’Anese, “Beyond Relaxation and Newton-Raphson: Solving AC-OPF for Multi-Phase Systems with Renewables,” IEEE Trans. on Smart Grid, vol. 9, no. 5, Sept. 2016.
  * [8] L. Ochoa and G. Harrison, “Minimizing Energy Losses: Optimal Accommodation and Smart Operation of Renewable Distributed Generation,” IEEE Trans.on Power Sys., vol. 26, no. 1, Feb. 2011.
  * [9] S. Claeys, F. Geth, and G. Deconinick, “Optimal Power Flow in Four-Wire Distribution Networks: Formulation and Benchmarking,” Electric Power Sys. Res., vol. 213, Dec. 2022.
  * [10] E.Foster, A. Pandey, and L. Pileggi, “Three-Phase Infeasibility Analysis for Distribution Grid Studies,” Electric Power Sys. Res., vol. 212, Nov. 2022.
  * [11] T. J. Overbye,“A power flow measure for unsolvable cases,” IEEE Trans. on Power Systems, vol. 9, no. 3, Aug. 1994.
  * [12] M. Jereminov, D. M. Bromberg, A. Pandey, M. R. Wagner and L. Pileggi, “Evaluating Feasibility Within Power Flow,” IEEE Trans.on Smart Grid, vol. 11, no. 4, July 2020.
  * [13] S. Li, A. Pandey, A. Agarwal, M. Jereminov and L. Pileggi, ”A LASSO-Inspired Approach for Localizing Power System Infeasibility,” presented at 2020 IEEE PES General Meeting, Montreal, QC, Canada, Aug. 2-6, 2020.
  * [14] A. Pandey, M. Jereminov, M. R. Wagner, G. Hug and L. Pileggi, ”Robust Convergence of Power Flow Using TX Stepping Method with Equivalent Circuit Formulation,” presented at the Power Sys. Comp. Conf. (PSCC), Dublin, Ireland, June 11-15, 2018.
  * [15] M. Jerminov, A. Pandey, and L. Pileggi, “Equivalent Circuit Formulation for Solving AC Optimal Power Flow,” IEEE Trans. on Power Sys., vol. 34, no. 3, May 2019.
  * [16] IEEE Standard for Shunt Power Capacitors, IEEE Std 18-2012 (Revision of IEEE Std 18-2002), Feb. 15, 2013.
  * [17] A. Pandey, et al., “Robust Power Flow and Three-Phase Power Flow Analyses,” in IEEE Trans. on Power Sys., vol. 34, no. 1, Jan. 2019.
  * [18] S. Boyd, S. P. Boyd, and L. Vandenberge, Convex Optimization, Cambridge University Press, 2004.
  * [19] T. McNamara, A. Pandey, A. Agarwal, and L. Pileggi, “Two-stage homotopy method to incorporate discrete control variables into AC-OPF,” Electric Power Sys. Res., vol. 212, Nov. 2022.
  * [20] A. Pandey, N. Turner-Bandele, E. Foster, and T. McNamara, SUGAR-D: A Distribution Systems Analysis and Optimization Tool v2.02, CMU ECE Pileggi Lab Group, Pittsburgh, PA, 2022.
  * [21] K.P. Schneider, et al., “Modern Grid Initiative Distribution Taxonomy Final Report,” Pacific Northwest National Lab, Richland, WA, Nov. 2008.
  * [22] E. Foster, TPIA Experiment Testcases, 2022, Github Repository. https://github.com/emfoster/TPIA-Experiment-Testcases


# ENHANCING THE ACCURACY OF THE TAYLOR POLYNOMIAL BY DETERMINING THE REMAINDER
TERM

J.S.C Prentice<EMAIL_ADDRESS>
( )

###### Abstract

We determine the Lagrange function in Taylor polynomial approximation by
solving an appropriate initial-value problem. Hence, we determine the
remainder term which we then approximate by means of a natural cubic spline.
This results in a significant improvement in the quality of the Taylor
approximation. We observe improvements in the accuracy of the approximation of
many orders of magnitude, including a case when the independent variable $x$
lies beyond the relevant radius of convergence.

Key words: Taylor, approximation, remainder, Lagrange function, cubic spline

##  1 Introduction

In a recent paper [1], we determined the remainder term in Lagrange
interpolation by solving suitable initial-value problem. In this paper, we
apply the same notion to Taylor approximation. We derive a differential
equation that will allow us to compute the reminder term (in Lagrange form) of
a first-order Taylor expansion. We then consider how a cubic spline can be
used to approximate the remainder term, leading to a significant improvement
in the accuracy of the approximation.

##  2 Relevant Concepts

Let $y\left(x\right)$ be a real-valued function, and assume that
$y\left(x\right)$ is as differentiable as is required in this paper (we will
make comments regarding this point when necessary). Taylor’s theorem [2]
provides the following result:

$y\left(x\right)=y\left(x_{0}\right)+y^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)+\frac{y^{\prime\prime}\left(\xi_{x}\right)}{2}\left(x-x_{0}\right)^{2}$
(1)

where $x_{0}<\xi_{x}<x.$ The third term on the RHS is the remainder term,
presented here in Lagrange form [3]. There are other representations of the
remainder term, such as the Cauchy form [3]

$y^{\prime\prime}\left(\xi_{C}\right)\left(x-\xi_{C}\right)\left(x-x_{0}\right),\text{
\ }x_{0}<\xi_{C}<x$

and the more general form

$y^{\prime\prime}\left(\xi_{C}\right)\left(x-\xi_{g}\right)\left(x-x_{0}\right)\left(\frac{g\left(x\right)-g\left(x_{0}\right)}{g^{\prime}\left(\xi_{g}\right)}\right),\text{
\ }x_{0}<\xi_{g}<x$ (2)

of which the Lagrange and Cauchy forms are particular cases. In (2), the
function $g$ is continuous on $\left[x_{0},x\right]$ and differentiable with a
non-vanishing derivative on $\left(x_{0},x\right).$ However, it is the
Lagrange form that will most easily serve our purposes here.

If we define

$T_{1}\left(x\right)\equiv
y\left(x_{0}\right)+y^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)$

to be the so-called Taylor polynomial of first degree, then we have

$y\left(x\right)-T_{1}\left(x\right)=\frac{y^{\prime\prime}\left(\xi_{x}\right)}{2}\left(x-x_{0}\right)^{2}.$
(3)

##  3 The Lagrange Function

By differentiating (3) with respect to $x,$we find

$y^{\prime}\left(x\right)-T_{1}^{\prime}\left(x\right)=y^{\prime\prime}\left(\xi_{x}\right)\left(x-x_{0}\right)+\frac{y^{\prime\prime\prime}\left(\xi_{x}\right)\left(x-x_{0}\right)^{2}}{2}\frac{d\xi_{x}}{dx}.$
(4)

Since the LHS of this expression exists, we must assume that the RHS also
exists. Hence, we must assume that $\frac{d\xi_{x}}{dx}$ exists. This implies
that $\xi_{x}$ is a function of $x,$ and we will sometimes write
$\xi\left(x\right)$ or simply $\xi$ in place of $\xi_{x}.$ We will refer to
$\xi\left(x\right)$ as the Lagrange function. In (4), we have used the
notation

$\displaystyle y^{\prime\prime}\left(\xi\right)$ $\displaystyle\equiv$
$\displaystyle y^{\prime\prime}\left(\xi\left(x\right)\right)$ $\displaystyle
y^{\prime\prime\prime}\left(\xi\right)$ $\displaystyle\equiv$
$\displaystyle\frac{dy^{\prime\prime}\left(\xi\left(x\right)\right)}{d\xi\left(x\right)}=\left.\frac{dy^{\prime\prime}\left(x\right)}{dx}\right|_{x=\xi}.$

Rearranging (4) and using
$T_{1}^{\prime}\left(x\right)=y^{\prime}\left(x_{0}\right)$ gives

$\frac{d\xi}{dx}=\frac{2\left(y^{\prime}\left(x\right)-y^{\prime}\left(x_{0}\right)-y^{\prime\prime}\left(\xi\right)\left(x-x_{0}\right)\right)}{y^{\prime\prime\prime}\left(\xi\right)\left(x-x_{0}\right)^{2}}.$
(5)

This is an initial-value problem that can, in principle, be solved to yield
the Lagrange function $\xi\left(x\right)$ for a suitable initial value. Once
$\xi\left(x\right)$ is known, the remainder term is easily computed. Hence, if
$T_{1}\left(x\right)$ is taken as an approximation to $y\left(x\right),$ the
resulting approximation error is known. We could also obtain differential
equations for $\xi_{C}$ and $\xi_{g}$ in a similar fashion, but we believe
these would be more complex than (5), and so we work only with $\xi$ in this
paper - an application of Occam’s Razor, we suppose.

There is a subtle point to be made: in (1), the remainder is derived using the
Mean Value Theorem [4], hence the constraint $x_{0}<\xi_{x}<x.$ Our assumption
above that $\xi_{x}=\xi\left(x\right)$ may be at odds with this constraint.
The Lagrange form of the remainder term does not explicitly require that
$\xi_{x}$ be a smooth function of $x.$ This is not a concern for us. Our task
here is simply to find values for $\xi$ that satisfy (1), by solving (5),
whether or not they satisfy the constraint. However, we will demonstrate later
how the solution(s) of (5) can be made consistent with the constraint.

Lastly, note that our analysis up to this point requires that
$y\left(x\right)$ is three times differentiable, at least.

##  4 Numerical Examples

###  4.1 First example

For our first example, we consider

$y\left(x\right)=e^{x/5}\sin x$

for $x\in\left[1,10\right].$ So we have $x_{0}=1.$ To find a suitable initial
value for $\xi,$ we choose a point very close to $x_{0}\ $(which we denote
$x_{z}$), say $x_{z}=1.0005$. We then solve

$y\left(x_{z}\right)=y\left(x_{0}\right)+y^{\prime}\left(x_{0}\right)\left(x_{z}-x_{0}\right)+\frac{y^{\prime\prime}\left(\xi_{z}\right)}{2}\left(x_{z}-x_{0}\right)^{2}$
(6)

numerically to find $\xi_{z}=1.000167.$ Next, we use the initial value
$\left(x_{z},\xi_{z}\right)$ to solve (5) using a seventh-order Runge-Kutta
(RK7) method [5] to find $\xi\left(x\right).$ This then allows us to find the
remainder

$R_{\xi}\left(x\right)\equiv\frac{y^{\prime\prime}\left(\xi\left(x\right)\right)}{2}\left(x-x_{0}\right)^{2},$

which can be compared with the actual remainder

$R_{act}\left(x\right)\equiv
y\left(x\right)-T_{1}\left(x\right)=y\left(x\right)-y\left(x_{0}\right)-y^{\prime}\left(x_{0}\right)\left(x-x_{0}\right).$

We can measure the quality of $R_{\xi}\left(x\right)$ through the device of

$\Delta R\left(x\right)\equiv R_{act}\left(x\right)-R_{\xi}\left(x\right).$

For this example, we find

$\max_{\left[1,10\right]}\left|\Delta R\left(x\right)\right|=5.4\times
10^{-13},$

suggesting that our calculation of $R_{\xi}\left(x\right)$ is very accurate,
indeed.

It transpires that (6) has another solution, leading to $\xi_{z}=3.157781.$ We
use this initial value in the same way as before, eventually finding

$\max_{\left[1,10\right]}\left|\Delta R\left(x\right)\right|=3.0\times
10^{-13},$

which shows accuracy similar to the first case.

In Figure 1, we show $\xi\left(x\right)$ for both initial values. On this plot
we show the constant $x_{0}=1$ and the line $y=x.$ These are the bounds on the
constraint $x_{0}<\xi_{x}<x,x\in\left[1,10\right].$ For
$\left(x_{z},\xi_{z}\right)=\left(1.0005,1.000167\right),$ we see that
$\xi\left(x\right)$ lies between the two bounds up to $x\approx 5.4,$ and lies
below $x_{0}=1$ thereafter. By contrast, for
$\left(x_{z},\xi_{z}\right)=\left(1.0005,1.000167\right),$
$\xi\left(x\right)>x$ up to $x\approx 2.8$ and thereafter lies between the
bounds. Neither of these Lagrange functions satisfies the constraint on the
entire interval $\left[1,10\right].$ However, the Lagrange form of the
remainder term, derived using the Mean Value Theorem, does not require that
$\xi_{x}$ must be smooth or even continuous. It simply states that, for each
$x\in\left[1,10\right],$ there does exist an $\xi_{x}$ such that
$x_{0}<\xi_{x}<x$ and $\xi_{x}$ satisfies (6). If we choose values of
$\xi\left(x\right)$ from the first Lagrange function for
$x\in\left[1,4\right],$ and values of $\xi\left(x\right)$ from the second
Lagrange function for $x\in\left(4,10\right],$ we necessarily satisfy Taylor’s
theorem for this example.

In Figure 2 we show $R_{\xi}\left(x\right)$ and $R_{act}\left(x\right),$ and
in Figure 3 we show $\Delta R\left(x\right).$

###  4.2 Second example

For our second example, we consider

$y\left(x\right)=\ln\left(1+x\right)$

for $x\in\left[0,10\right].$ We have $x_{0}=0.$ Again, we choose
$x_{z}=1.0005$. We then solve (6) to find $\xi_{z}=1.67\times 10^{-4}.$
Application of RK7 then gives the Lagrange function $\xi\left(x\right),$ shown
in Figure 4, along with the bounds of the constraint $x_{0}=0$ and the line
$y=x.$ We see that $\xi\left(x\right)$ lies between the bounds for all
$x\in\left[0,10\right].$ In Figure 5 we show $R_{\xi}\left(x\right)$ and
$R_{act}\left(x\right),$ and in Figure 6 we show $\Delta R\left(x\right).$ We
determine

$\max_{\left[0,10\right]}\left|\Delta R\left(x\right)\right|=1.3\times
10^{-13}.$

##  5 Polynomial Approximation of the Remainder

Now that we have determined the Lagrange function $\xi\left(x\right)$ and,
hence, the remainder term $R_{\xi}\left(x\right)$ for each example, it seems
reasonable to attempt to approximate the remainder term by means of a
polynomial. This polynomial can be added to the original Taylor polynomial
$T_{1}\left(x\right)$ to yield a more accurate approximation than
$T_{1}\left(x\right).$ If we denote the polynomial approximation to
$R_{\xi}\left(x\right)$ by $P_{R}\left(x\right),$ we would then present

$T_{1}\left(x\right)+P_{R}\left(x\right)$

as an approximation to $y\left(x\right)$ on the given interval.

We choose to use a natural cubic spline [6] to form $P_{R}\left(x\right),$ for
several reasons: we have the RK nodes at our disposal (we used $10000$ nodes
in the RK computations); a cubic polynomial combines with the factor
$\left(x-x_{0}\right)^{2}$ in $R_{\xi}\left(x\right)$ to yield a polynomial of
degree five, at most; and such splines can be generated very efficiently on
our computational platform [7]. Also, it is possible to estimate a bound on
the accuracy of a natural spline although, as will be seen, the bound is not
tight.

Defining

$\displaystyle\Delta_{T}$
$\displaystyle\equiv\max_{I}\left|y\left(x\right)-T_{5}\left(x\right)\right|$
$\displaystyle\Delta_{CS}$
$\displaystyle\equiv\max_{I}\left|y\left(x\right)-\left(T_{1}\left(x\right)+P_{R}\left(x\right)\right)\right|$

where $T_{5}\left(x\right)$ is the Taylor polynomial of fifth degree, and $I$
is the relevant interval of approximation, we show results in Table 1.

$y\left(x\right)$ | $I$ | $\Delta_{T}$ | $\Delta_{CS}$ | $B_{U}$
---|---|---|---|---
$e^{x/5}\sin x$ | $\left[1,10\right]$ | $5.8\times 10^{2}$ | $5.4\times 10^{-13}$ | $5.1\times 10^{-10}$
$\ln\left(1+x\right)$ | $\left[0,10\right]$ | $1.8\times 10^{4}$ | $1.4\times 10^{-13}$ | $8.6\times 10^{-9}$

Table 1: Results for cubic spline and Taylor approximations.

In Table 1, $B_{U}$ is an upper bound on the error in the cubic spline
approximation, estimated using [8]

$B_{U}=72h^{4}\max_{I}\left|y^{\left(6\right)}\left(x\right)\right|$ (7)

where $h$ is the uniform RK stepsize (see Appendix). The bounds are clearly
generous and, for each example, the actual error is less than the bound. It is
abundantly clear that $T_{1}\left(x\right)+P_{R}\left(x\right)$ is a far
better approximant than $T_{5}\left(x\right).$ This is particularly true for
the second example, whose Taylor series has a radius of convergence of
$\left|x\right|<1$ (hence the very large value of $\Delta_{T}).$

##  6 Concluding Comments

We have shown how the Lagrange function in Taylor polynomial approximation can
be determined by solving an appropriate initial-value problem. This allows the
remainder term to be determined. The remainder term can then be approximated
by means of a polynomial, and this can result in a significant improvement in
the quality of the Taylor approximation overall. We have demonstrated this
effect using a cubic spline, and we note improvements in the accuracy of the
approximation of many orders of magnitude, including the case when the
independent variable $x$ lies beyond the radius of convergence. This speaks to
the potential value of the idea presented here, and in [1].

This paper is intended as a demonstration of an idea and, as such, our
analysis has not been exhaustive. Further studies should consider the effect
of error control in the RK solution, and how this will affect the quality of
the polynomial approximation of the remainder term. It is also feasible to
approximate the remainder term by means of a least-squares fit, which might be
worth considering. Also, the multidimensional case must be investigated, which
would require the solution of a system of differential equations, but all
these aspects will be reserved for future work.

## References

  * [1] J.S.C. Prentice, Determining the Rolle function in Lagrange interpolatory approximation, arXiv.org, Cornell University Library, 2018. [arXiv: 1810.09601].
  * [2] M. Spivak, Calculus 3rd ed., Publish or Perish, Houston, 1994.
  * [3] T. Apostol, Calculus, Wiley, 1967.
  * [4] E. Boman and R. Rogers, Lagrange’s Form of the Remainder, Chapter 5 in Real Analysis, LibreTexts${}^{\text{TM}}$ Mathematics (https://math.libretexts.org).
  * [5] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, Chichester, 2003.
  * [6] R.L. Burden & J.D. Faires, Numerical Analysis 9th ed., Brooks/Cole, Boston, 2011.
  * [7] Intel i7 9700K CPU, Strix Z390-F mainboard, 48GB DDR4 RAM, Windows 10 64-bit, Matlab R2022b.
  * [8] G. Birkhoff and C. de Boor, Error Bounds for Spline Interpolation, Journal of Mathematics and Mechanics, 13, 5: 827-835, 1964.

##  7 Appendix

Here we derive the bound in (7). We refer to [8]. Following [8] (see their eqn
(27)), but using our own notation, we have

$\displaystyle e\left(x\right)$ $\displaystyle\equiv
y^{\prime\prime}\left(x\right)-P_{R}\left(x\right)$
$\displaystyle\left|e^{\prime\prime\prime}\right|$
$\displaystyle=\left|y^{\left(5\right)}\left(x\right)-P_{R}^{\prime\prime\prime}\left(x\right)\right|\leqslant
3M\left(1+M\right)^{2}\max_{I}\left|\frac{d^{4}y^{\prime\prime}}{dx^{4}}\right|\max_{I}\left|h_{i}\right|$

where $h_{i}$ denotes the spacing between the RK nodes (which is not
necessarily uniform) and

$M=\frac{\max_{I}\left|h_{i}\right|}{\min_{I}\left|h_{i}\right|}.$

In our calculations, the RK nodes are uniformly spaced (denoted $h$) so that
$M=1.$ Hence,

$\left|y^{\left(5\right)}\left(x\right)-P_{R}^{\prime\prime\prime}\left(x\right)\right|\leqslant
12\max_{I}\left|y^{\left(6\right)}\right|h.$

Eqn (29) in [8] gives the recursion

$\displaystyle\left|y^{\left(4\right)}\left(x\right)-P_{R}^{\prime\prime}\left(x\right)\right|$
$\displaystyle\leqslant\left(3\right)\left|y^{\left(5\right)}\left(x\right)-P_{R}^{\prime\prime\prime}\left(x\right)\right|h$
$\displaystyle\left|y^{\left(3\right)}\left(x\right)-P_{R}^{\prime}\left(x\right)\right|$
$\displaystyle\leqslant\left(2\right)\left|y^{\left(4\right)}\left(x\right)-P_{R}^{\prime\prime}\left(x\right)\right|h$
$\displaystyle\left|y^{\prime\prime}\left(x\right)-P_{R}\left(x\right)\right|$
$\displaystyle\leqslant\left(1\right)\left|y^{\left(3\right)}\left(x\right)-P_{R}^{\prime}\left(x\right)\right|h$

so that

$\displaystyle\left|y^{\prime\prime}\left(x\right)-P_{R}\left(x\right)\right|$
$\displaystyle\leqslant
6\left|y^{\left(5\right)}\left(x\right)-P_{R}^{\prime\prime\prime}\left(x\right)\right|h^{3}$
$\displaystyle\leqslant 72\max_{I}\left|y^{\left(6\right)}\right|h^{4}.$

Note that the above analysis requires that $y\left(x\right)$ is six times
differentiable, at least.