# Transmission Investment Coordination using MILP Lagrange Dual Decomposition
and Auxiliary Problem Principle

Sambuddha Chakrabarti, Hosna Khajeh, Thomas R Nudell, Mohammad Reza
Hesamzadeh, and Ross Baldick S. Chakrabarti is currently a staff researcher at
Princeton University, Princeton, NJ, USA. He was previously a postdoctoral
researcher at KTH Royal Institute of Technology, Stockholm, Sweden. H. Khajeh
is a doctoral researcher at the School of Technology and Innovations,
University of Vaasa, Vaasa, Finland. T. Nudell is the cofounder and CEO of Piq
Energy Corp, Union City, CA, USA. M. R. Hesamzadeh is an associate professor
at KTH Royal Institute of Technology. R. Baldick is a professor emeritus with
the Department of Electrical & Computer Engineering, The University of Texas
at Austin, Austin, TX, USA. This work was supported in parts by Svenska
Kraftnät & by National Science Foundation (grant number:ECCS-1406894)

###### Abstract

This paper considers the investment coordination problem for the long term
transmission capacity expansion in a situation where there are multiple
regional Transmission Planners (TPs), each acting in order to maximize the
utility in only its own region. In such a setting, any particular TP does not
normally have any incentive to cooperate with the neighboring TP(s), although
the optimal investment decision of each TP is contingent upon those of the
neighboring TPs. A game-theoretic interaction among the TPs does not
necessarily lead to this overall social optimum. We, therefore, introduce a
social planner and call it the Transmission Planning Coordinator (TPC) whose
goal is to attain the optimal possible social welfare for the bigger
geographical region. In order to achieve this goal, this paper introduces a
new incentive mechanism, based on distributed optimization theory. This
incentive mechanism can be viewed as a set of rules of the transmission
expansion investment coordination game, set by the social planner TPC, such
that, even if the individual TPs act selfishly, it will still lead to the
TPC’s goal of attaining overall social optimum. Finally, the effectiveness of
our approach is demonstrated through several simulation studies.

###### Index Terms:

Distributed Optimization Theory, Incentive Mechanism Design, Investment
Coordination.

##  I Notations and Conventions

  * •

Sets

    * $\mathcal{N}$: Set of nodes.

    * $\mathcal{N}_{z}$: Set of nodes in region $z$.

    * $\mathcal{N}_{z^{\prime}}$: Set of nodes in region $z^{\prime}$ (i.e. neighbors of $z$).

    * $G_{z}$: Set of generators in region $z$.

    * $G_{zn}$: Set of generators in $z$ connected to node $n$.

    * $H_{z}$: Set of existing transmission lines in region $z$.

    * $H_{zn}$: Set of existing lines in $z$ connected to node $n$.

    * $\tilde{H}$: Set of existing lines shared between multiple regions

    * $K_{z}$: Set of candidate lines in region $z$.

    * $D_{z}$: Set of loads in region $z$.

    * $D_{zn}$: Set of loads in region $z$ connected to node $n$.

    * $Z$: Set of regions or sub-networks.

    * $S$: Set of scenarios.

    * $\dagger$ used to denote the transpose of a vector or matrix

  * •

Indices

    * $n$: Index for the elements of $\mathcal{N}$.

    * $g$: Index for the elements of $G$.

    * $h$: Index for the elements of $H$ and $\tilde{H}$.

    * $k$: Index for the elements of $K$.

    * $d$: Index for the elements of $D$.

    * $z$: Index for the elements of $Z$.

    * $s$: Index for the elements of $S$.

  * •

Parameters

    * $x_{h}$: Reactance of the transmission line $h$.

    * $\hat{x}_{k}$: Reactance of the candidate transmission line $k$.

    * $\overline{L}_{h}$: Capacity of the transmission line $h$.

    * $\hat{L}_{k}$: Maximum expansion capacity for the candidate line $k$.

    * $\mathbf{A}_{(n,h)}$: $(n,h)^{th}$ element of the node-to-branch incidence matrix $A$. The element is 1 when node $n$ is the sending end of the line $h$, -1 when it is the receiving end, and 0 otherwise.

    * $w_{s}$: Weight of scenario, $s$.

    * $T_{k}$: Lifetime of the candidate line $k$.

    * $r$: Interest rate.

    * $P_{sd}$: Demand of the $d^{th}$ load in the scenario $s$.

    * $\underline{P}_{sg}$, $\overline{P}_{sg}$: Lower and upper generating limits of generator $g$ in scenario $s$.

    * $i_{k}$, $j_{k}$: Sending and receiving ends of candidate line $k$ respectively.

    * $\eta,\gamma,\delta_{APP}$: First and second penalty parameters and the step-length parameter of the APP iterations, respectively.

  * •

Variables & Functions

    * $P_{sg}$: Power generation of generator $g$ in scenario $s$.

    * $\theta_{sn}$: Node voltage phase angle of node $n$ in scenario $s$.

    * $u_{zk}$: Binary decision variable with $0$ for no capacity expansion and $1$ for capacity expansion to the maximum limit for the candidate line $k$ in region $z$.

    * $\hat{P}_{sk}$: Power flow on candidate line $k$ in scenario $s$.

    * $P_{sh}$: Power flow on existing line $h$ in scenario $s$.

    * $C_{g}(P_{sg})$: Cost function of the generator $g$ (typically convex quadratic or convex cubic), reflecting the fuel cost and heat rate, for producing power, $P_{sg}$.

    * $C_{k}(\hat{L}_{k})$: Cost of building a new line $k$ of capacity $\hat{L}_{k}$.

    * $\lambda^{\sigma}_{zz^{\prime}i(h)},\lambda^{\sigma}_{zz^{\prime}i(h)}$: Dual variable in the Auxiliary Problem Principle (APP) iterations for consensus on flow value on the line $h$ shared between regions $z$ and $z^{\prime}$ at the iteration count $\sigma$.

##  II Introduction

The contemporary electric power network is a complex, large, and dynamically
engineered system. It is old, with the average age of transmission assets
nearly forty years old. Yet it is constantly evolving to accommodate growing
demand and also to incorporate the new technologies like renewable generation,
modern communication, and demand-response equipment. The consideration of
competitive market setting for the expansion planning problems makes it even
more complex [5]. As such, the expansions of both the generation, as well as
the transmission infrastructures are of paramount importance in the long-term
planning of power systems across the world. In this paper, we focus on the
problem of transmission investment, particularly in the situation where there
are multiple transmission planners who do not necessarily have any incentive
to cooperate with each other.

###  II-A Literature Review

Previous work has addressed the generation and transmission expansion planning
problem in a centralized manner. These studies mainly concentrated on how to
develop generation capacities over a specific time horizon. For example, [11]
assessed how the electric power industry can change and develop its mix of
generating capacities over time according to the future demand, electricity
and fuel prices, future regulations, and technological costs. In another
similar work introduced by [12], U.S. Environmental Protection Agency (EPA)
tried to provide optimal options for the generation capacity expansion
problem. The options include utilizing the demand-side participation,
renewable resources, and traditional generation capacities in the future.
Reference [38] also considered the effects of customer choices and end-user
services on the energy demand and therefore on the generation capacity
expansion model. In another work introduced by [13] and [21], a capacity
expansion model is solved through three different modules. The first module is
named supply module, which aims to obtain the cost-efficient investment and
operation of electrical sectors. The second module is a demand module trying
to maximize the levels of investment on end-user devices, while the third
module analyzes the optimal values and key parameters of renewable generators.
A few modules written in Python such as “Switch” also provide a thorough
generation expansion planning method with regard to the existing dispatch and
investment rules for different technologies. These technologies can include
storage and demand response resources [17].

In a few studies, the authors focused on the problem of transmission and
generation expansion planning, simultaneously. For instance, PyPSA has worked
on transmission expansion based on integer decision variables and big-M
relaxations [3]. Moreover, GenX model defines a generation expansion planning
model considering both investment on centralized and distributed generation
capacities, storage-based and demand-side technologies [15] and [16].
Additionally, the model takes into account the transmission network expansion
decisions as well, but on a regional aggregated or commercially significant
constraint basis. In all these studies, the problem of transmission and
generation expansion planning has been defined in a centralized fashion and
prescriptive mode.

In the context of real-world transmission networks, the problem of generation
and transmission expansion is seldom solved in a centralized manner since the
fundamental nature of the problem is that of decentralized decision-making. In
reality, there are multiple regions and each region has its separate
transmission and generation planner. Each transmission planner is in charge of
solving an expansion planning problem in its own region. However, the previous
literature did not address the issue of “how to realize the prescriptive
expansion recommendations or what kind of incentivizing schemes to design in
order to convert their results to reality.” There have been problems of
similar nature and mathematical structure to that of our present problem,
which has been solved for different application domains. For instance, in the
area of the Internet of Things (IoT), IoT devices cooperatively execute tasks
in a decentralized way [14], or [39] put forward an incentive mechanism in
which vehicle users deploy a non-cooperative game to create service requests,
and the system accordingly guides the best available offloading strategy.
Reference [37] suggested the decentralized routing way for Internet Energy
which maintains a certain utilization rate and facilitates fast learning from
the mechanism.

To resolve the multi-regional issue, some research discussed the planning
problem with multiple transmission planners (TPs) in different regions using
game-based approaches. However, the main issue related to these approaches is
that the problem may not lead to optimal points for all of the TPs. In this
context, authors of [32] have solved the multi-region transmission expansion
planning problem by applying the non-cooperative game-theoretic model in a bi-
level optimization formulation, through the use of Nash Equilibrium (NE)
solution concept. As illustrated in that paper, the outcome of such an
approach typically happens to be sub-optimal compared to the centralized
solution where it is assumed that all the different transmission planners
(TPs) merge into one entity. The authors also alluded to different
compensation mechanisms in that reference which can improve the results. In
[30], the authors have mentioned the interesting effect of “free-riding,” in
the context of such non-cooperative, multi-regional transmission planning,
which is the presence of entities that make benefit out of expansion by
others, without themselves contributing anything. The authors have studied the
impact of this free-riding on congestion revenues. In the reference [31], the
authors have tried to improve the social welfare as presented in [32] by
adopting different solution approaches, whereas in [33], the multi-regional
transmission planing has been explored in the presence of wind generation and
its associated uncertainty. References [36], [34], and [35] have explored the
transmission expansion planning in the presence of generation expansion
planning considering pro-active and reactive coordination methods. Reference
[29] is a good comprehensive collection of the work that has been performed so
far in this direction.

It is important to observe that consideration of contingencies and pre- and
post-contingency line flows are important metrics for the solution of the
transmission capacity expansion problem. For obvious reasons, the inclusion of
contingencies increases the computational complexity of the problem.
References [25], [22], [23], and [24] explored the issue of inclusion of
contingencies in the context of both the deterministic and stochastic settings
as well as prescribed ways to eliminate umbrella contingency scenarios and
expedite the solution through decomposition methods in their papers.
Closely aligned to this problem is the problem of distributed optimization in
networked cyber-physical system, in an environment of very limited information
and bandwidth, which is presented in the work, [20]. Although the previous
literature mostly results in sub-optimal coordination outcome, this paper
introduces a Transmission Planning Coordinator (TPC) that guides Transmission
Planners (TP) with a view towards a near global social optimum, as much as
possible . It is to be noted, however though, that achieving global optimality
is not the primal aim of this work (this is so because, this problem has
integer decision variables and is solved in a decentralized manner. Therefore,
in principle, achieving global optimality is presumably impossible); however,
this paper tries to:

  * •

Simulate how transmission expansion works in real-world multi-region settings.

  * •

Come up with a market mechanism that presumably works better in comparison to
non-cooperative game-theoretic interactions among the regions (in the absence
of a social planner)

Hence, in the present case, we will assume the existence of a central rule-
making agent or market overseeing authority, who sets the rules of the game,
and the individual transmission and generation asset owners, who independently
solve their own optimization sub-problem.This idea has also been discussed
recently by the Federal Energy Regulatory Commission (FERC) at their staff-led
workshop 111See:https://www.ferc.gov/news-events/events/staff-led-workshop-
establishing-interregional-transfer-capability-transmission.We apply the ideas
from the world of distributed optimization and message-passing frameworks to
motivate the setting up of such a mechanism design scheme. Introducing the TPC
also unlocks a number of new mechanisms to improve social outcomes, including
enabling truth-telling behavior among participating regions based on Vickrey-
Clarke-Groves auction and the generalization to Vickrey-Clarke-Groves
mechanism. This paper does not elaborate on Vikrey-Clarke-Groves auction
and/or Vickrey-Clarke-Groves mechanism but takes the necessary step of
introducing the TPC to enable further research. In our work, we have adopted
the simplest 50-50 cost sharing of the new shared transmission lines to be
constructed between two regions. Adopting the VCG mechanism for inducing
truth-telling behavior implies modification of this cost-sharing and cost
allocation to each region in a way to capture the perceived utilities of all
the other regions, by each region. This is proven to induce truth-telling
behavior among participating regions when it comes to revealing their real
preference as to whether or not to build the shared transmission lines. Our
present work opens up avenues for such future work.

In order to motivate the search for an algorithmic framework, based on which
we design our incentive mechanism, we observe the following characteristics of
the present optimization problem:

  * •

The centralized optimization problem is a mixed integer programming problem
(MIP) since it has discrete binary decision variables as well as continuous
ones. The binary variables decide whether to build a particular line or not
along a candidate line, and continuous variables pertaining the generator
dispatch and voltage phase angles.

  * •

The bigger problem is split into several regions and the candidate pathways.
Moreover, some existing lines are shared between two different regions.

  * •

The decision variables corresponding to the shared lines are also shared among
each region and there needs to be consensus among the values decided by the
different regions.

In order to satisfy the above requirements, we propose a two staged incentive
mechanism design, which consists of the following stages or layers:

  * •

The first stage, based on Distributed Stochastic Optimization, solves for the
integer variables and approximately solves for the continuous variables.

  * •

The second stage, based on Auxiliary Problem Principle (APP), solves for the
continuous variables more accurately, once the discrete variable values are
determined.

For the first stage, we use the algorithm introduced in [1] for solving
Stochastic Unit Commitment problem in a distributed manner. It is an
asynchronous distributed algorithm for maximizing the Langrangian dual problem
of the classic Unit Commitment (UC) formulation. The synchronous version of
the algorithm was presented in the earlier works [27], [28], and [26]. In
those papers, the authors have determined the commitment decisions, which are
binary discrete variables, of slow responding generators across different
scenarios, through the use of “non-anticipativity constraints.” Instead of a
scenario decomposition, in our work, we will be performing a region
decomposition. For the second stage of the mechanism design, we will be making
use of the Auxiliary Problem Principle (APP), which was introduced by Cohen
_et al._ in the seminal works [9] and [8]. It was subsequently used in power
flow problems for coarse decomposition, both for a region decomposition in
[2], [19], [18], and [10] and for decomposition across different dispatch
intervals and contingency scenarios in [4], [6], [7].
Reflecting on what was stated, the novelty of our paper can be summarized as
follows:

  * •

Previous approaches have taken a non-cooperative game theoretic approach. This
results in sub-optimal coordination outcomes. By introducing a TPC, individual
TPs are incentivized to converge towards the near global social optimum while
still following their locally selfish objectives. This approach has not, to
our knowledge, been proposed in the context of inter-regional transmission
planning.

  * •

Although this paper has been inspired by the BCD method, it deals with a
totally different problem (transmission expansion planning problem) and the
algorithm has been completely modified to be compliant with the new
application.

  * •

We propose the novel application of APP in which TPs reach a consensus
regarding the flow and voltage angles of inter-regional shared lines between
them.

The rest of the current paper is organized as follows: In Section III we
present the centralized transmission expansion planning for reference and for
the sake of completeness. Section IV presents the mathematical formulation of
both the stages of our distributed algorithmic incentive mechanism design and
the steps thereof. In Subsection IV-A, we will introduce the first stage of
the mechanism with its mathematical model. In Subsection IV-B, we mention the
steps of the first stage incentive mechanism design. Subsection IV-C details
the second stage mechanism design model, which is based on the APP algorithm,
while Subsection IV-D walks us through the steps of the second stage mechanism
design. We state the results of numerical simulations in Section V and in
Section VI, we draw some concluding remarks, while pointing to the future
research directions.

##  III Centralized Coordination (Benchmark)

The centralized coordination of transmission investment is modeled in (1) to
create a benchmark for the distributed coordination algorithm. In the
centralized version, the TPC minimizes the total operations and investment
costs in all regions. The model presented here is a slight modification of the
one presented in [32]:

$\displaystyle\min_{P_{sg},u_{zk}}\sum_{z\in Z}\sum_{s\in
S}\sum_{g\in{G}_{z}}w_{s}C_{g}(P_{sg})$ $\displaystyle+\sum_{z\in Z}\sum_{k\in
K_{z}}u_{zk}C_{k}({\hat{L}}_{k})\frac{r(1+r)^{T_{k}}}{(1+r)^{T_{k}}-1}$ (1a)
$\displaystyle\mbox{Subject\;to:\>}\sum_{g\in G_{zn}}P_{sg}-\sum_{d\in
D_{zn}}P_{sd}=$ $\displaystyle\sum_{h\in
H_{zn}}\mathbf{A}_{(n,h)}P_{sh}+\sum_{k\in{K_{z}},i_{k}=n}\hat{P}_{sk}-\sum_{k\in{K_{z}},j_{k}=n}\hat{P}_{sk}$
$\displaystyle\forall{z\in Z},\forall{s\in S},\forall{{n}\in\mathcal{N}_{z}}$
(1b) $\displaystyle
P_{sh}=\frac{1}{x_{h}}\sum_{n\in\mathcal{N}}\mathbf{A}_{(n,h)}\theta_{sn},\forall{z\in
Z},\forall{s\in S},\forall{{h}\in H_{z}}$ (1c)
$\displaystyle\hat{P}_{sk}=u_{zk}\frac{1}{\hat{x}_{k}}(\theta_{si}-\theta_{sj}),k\rightarrow(i,j)$
$\displaystyle\forall{z\in Z},\forall{s\in S},\forall{{k}\in K_{z}}$ (1d)
$\displaystyle-{\overline{L}}_{h}\leq
P_{sh}\leq{{\overline{L}}_{h}},\;\forall{z\in Z},\forall{s\in
S},\forall{{h}\in H_{z}}$
$\displaystyle-{u_{zk}{\hat{L}}_{k}}\leq\hat{P}_{sk}\leq{u_{zk}{\hat{L}}_{k}},\;\forall{z\in
Z},\forall{s\in S},\forall{{k}\in K_{z}}$
$\displaystyle{\underline{P}}_{sg}\leq
P_{sg}\leq{{\overline{P}}_{sg}},\;\forall{z\in Z},\forall{s\in
S},\forall{g\in{G_{z}}}$ (1g) $\displaystyle u_{zk}\in\\{0,1\\},\;\forall{z\in
Z},\forall{{k}\in K_{z}},\;$ (1h)

In the objective function, the first term is the expected operational cost and
the second term is the total investment amount taking into account the
interest of the transmission investment cost. Constraint (1b) represents the
power balance for each node. Constraints (1c) and (1d) define line flows for
the existing lines and the candidate lines, respectively, while
(LABEL:third_Step2:General_one) and (LABEL:third_Step2II:General) define the
line flow limits for existing and candidate lines, respectively. The last two
constraints are for generation limits and for permissible values that the
binary variables for constructing new lines can assume.

##  IV Distributed Coordination Mechanism

In this section, we present the design of a novel two-stage incentive
mechanism. The first stage is based on, and inspired by the Distributed
Stochastic Mixed Integer Optimization and the concept of scenario
decomposition. We extend the idea of a “scenario decomposition” to “area
decomposition.” The outcome of the first stage determines the investment
decisions as to whether or not and where to invest in building new
transmission lines. It also gives us the approximate values of tie-line flows
and generation values.

###  IV-A Mathematical Model of Stage I Mechanism Design

For the first stage, we present an algorithmic incentive mechanism design
based on the asynchronous distributed algorithm presented by the authors in
[1]. In constraint (2d), $M$ is a large number. The purpose of this constraint
is to linearize the constraint (1d) in Section III, which is non-linear due to
the presence of product terms.

$\displaystyle\min_{P^{a}_{g_{q}}}\sum_{z\in Z}\Bigg{(}\sum_{s\in
S}\sum_{g\in{G_{z}}}w_{s}C_{g}(P_{sg})$ $\displaystyle+\sum_{k\in
K_{z}}\Big{\\{}u_{zk}C_{k}({\hat{L}}_{k})\frac{r(1+r)^{T_{k}}}{(1+r)^{T_{k}}-1}\Big{\\}}\Bigg{)}$
(2a) $\displaystyle\mbox{Subject\;to:\>}\sum_{g\in G_{zn}}P_{sg}-\sum_{d\in
D_{zn}}P_{sd}=$ $\displaystyle\sum_{h\in
H_{zn}}\mathbf{A}_{(n,h)}P_{sh}+\sum_{k\in{K_{z}},i_{k}=n}\hat{P}_{sk}-\sum_{k\in{K_{z}},j_{k}=n}\hat{P}_{sk}$
$\displaystyle\forall{z\in Z},\forall{s\in S},\forall{{n}\in\mathcal{N}_{z}}$
(2b) $\displaystyle
P_{sh}=\frac{1}{x_{h}}\sum_{n\in\mathcal{N}}\mathbf{A}_{(n,h)}\theta_{sn},\forall{z\in
Z},\forall{s\in S},\forall{{h}\in H_{z}}$ (2c)
$\displaystyle-M(1-u_{zk})\leq\hat{P}_{sk}-\frac{1}{\hat{x}_{k}}(\theta_{si}-\theta_{sj})\leq
M(1-u_{zk})$ $\displaystyle k\leftarrow(i,j),\forall{z\in Z},\forall{s\in
S},\forall{{k}\in K_{z}}$ (2d) $\displaystyle-{\overline{L}}_{h}\leq
P_{sh}\leq{{\overline{L}}_{h}},\;\forall{z\in Z},\forall{s\in
S},\forall{{h}\in H_{z}}$
$\displaystyle-{u_{zk}{\hat{L}}_{k}}\leq\hat{P}_{sk}\leq{u_{zk}{\hat{L}}_{k}},\;\forall{z\in
Z},\forall{s\in S},\forall{{k}\in K_{z}}$
$\displaystyle{\underline{P}}_{sg}\leq
P_{sg}\leq{{\overline{P}}_{sg}},\;\forall{z\in Z},\forall{s\in
S},\forall{g\in{G_{z}}}$ (2g) $\displaystyle u_{zk}\in\\{0,1\\},\;\forall{z\in
Z},\forall{{k}\in K_{z}},\;$ (2h) Non-Anticipativity Constraints
$\displaystyle u_{zk}=u_{k}\leftrightarrow\pi_{zk},\;\forall k\in K$ (2i)
$\displaystyle|\hat{P}_{zsk}|=\hat{\mathbb{P}}_{sk}\leftrightarrow\mu_{zsk},\;\forall
k\in K$ (2j)
$\displaystyle|{P}_{zsh}|=\mathbb{P}_{sh}\leftrightarrow\mu_{zsh},\;\forall
h\in\tilde{H}$ (2k)
$\displaystyle\theta_{zsi(k)}=\phi_{si(k)}\leftrightarrow\xi_{zsi(k)},\;\forall
k\in K$ (2l)
$\displaystyle\theta_{zsj(k)}=\phi_{sj(k)}\leftrightarrow\xi_{zsj(k)},\;\forall
k\in K$ (2m)
$\displaystyle\theta_{zsi(h)}=\phi_{si(h)}\leftrightarrow\xi_{zsi(h)},\;\forall
h\in\tilde{H}$ (2n)
$\displaystyle\theta_{zsj(h)}=\phi_{sj(h)}\leftrightarrow\xi_{zsj(h)},\;\forall
h\in\tilde{H}$ (2o)
$\displaystyle\mathbb{P}_{sh}=\frac{1}{x_{h}}(\phi_{si(h)}-\phi_{sj(h)}),\forall{s\in
S},\forall{{h}\in\tilde{H}}$ (2p)
$\displaystyle-M(1-u_{k})\leq\hat{\mathbb{P}}_{sk}-\frac{1}{\hat{x}_{k}}(\phi_{si(k)}-\phi_{sj(k)})\leq
M(1-u_{k})$ $\displaystyle u_{k}\in\\{0,1\\},k\leftarrow(i,j),\forall{s\in
S},\forall{{k}\in K}$ (2q)

The region-decomposed Lagrangian dual of the above optimization problem is
stated below in mathematical expressions (3a)–(3c). This set of problems is
solved by the different TPs while being coordinated by the TPC. The penalty
cost appearing to the right of $\leftrightarrow$ plays the pivotal role of
motivating the TPs to perform such calculations and report the updated
variable values to the TPC.

$\displaystyle\max_{\mathbf{\pi},\mathbf{\mu},\mathbf{\xi}}f_{0}(\mathbf{\pi},\mathbf{\mu},\mathbf{\xi})+\sum_{z\in
Z}f_{z}(\mathbf{\pi},\mathbf{\mu},\mathbf{\xi}),\;\text{where}$ (3a)
$\displaystyle f_{0}(\mathbf{\pi},\mathbf{\mu},\mathbf{\xi})$
$\displaystyle=\inf_{u_{k},\hat{\mathbb{P}}_{sk},\mathbb{P}_{sh},\phi_{k},\phi_{h}}\Bigg{(}-\sum_{z\in
Z}\sum_{k\in K}\\{u_{k}\pi_{zk}+\sum_{s\in S}(\hat{\mathbb{P}}_{sk}\mu_{zsk}$
$\displaystyle+\phi_{si(k)}\xi_{zsi(k)}+\phi_{sj(k)}\xi_{zsj(k)})\\}-\sum_{z\in
Z}\sum_{h\in\tilde{H}}\sum_{s\in S}\\{\mathbb{P}_{sh}\mu_{zsh}$
$\displaystyle+\phi_{si(h)}\xi_{zsi(h)}+\phi_{sj(h)}\xi_{zsj(h)}\\}:\text{Subject
to: }(\ref{third_Step9:General})-(\ref{third_Step10:General})\Bigg{)}$ (3b)
$\displaystyle
f_{z}(\mathbf{\pi},\mathbf{\mu},\mathbf{\xi})=\inf_{P_{g},u_{zk}}\Bigg{(}\sum_{s\in
S}\sum_{g\in{G_{z}}}w_{s}C_{g}(P_{sg})+$ $\displaystyle\sum_{k\in
K_{z}}\Big{\\{}u_{zk}C_{k}({\hat{L}}_{k})\frac{r(1+r)^{T_{k}}}{(1+r)^{T_{k}}-1}+u_{zk}\pi_{zk}$
$\displaystyle+\sum_{s\in
S}(|\hat{P}_{zsk}|\mu_{zsk}+\theta_{zsi(k)}\xi_{zsi(k)}+\theta_{zsj(k)}\xi_{zsj(k)})\Big{\\}}$
$\displaystyle+\sum_{s\in
S}\sum_{h\in\tilde{H}}\Bigg{\\{}|{P}_{zsh}|\mu_{zsh}+\theta_{zsi(h)}\xi_{zsi(h)}+\theta_{zsj(h)}\xi_{zsj(h)}\Bigg{\\}}$
$\displaystyle:\text{Subject to:
}(\ref{second_Step2:General})-(\ref{sixth_Step2II:General})\Bigg{)}$ (3c)

Figure 1: Illustration of the incentive mechanism as applied to transmission
expansion planning

If all the generator cost curves are considered as linear or piece-wise
linear, then both (3b) and (3c) are simply Mixed Integer Linear Programming
(MILP) problems, for which there exist standard algorithms for solving.

###  IV-B Steps of Stage I Mechanism Design

We will now state the steps of our algorithmic incentive mechanism. These are
very similar to the ones that appear in [1], with the only significant
differences being instead of scenario decomposition, we do a region or area
decomposition here. In our mechanism formulation, we will apply the Block
Coordinate Descent (BCD) method to reach the consensus regarding the decision
to build new transmission lines. The following are the steps of the
algorithmic incentive mechanism design:

  * •

Each TP of region $z$ evaluates its own
$f_{z}(\pi^{\nu(z)},\mu_{s}^{\nu(z)},\xi_{s}^{\nu(z)})$ based on its current
updated rewards/penalties for consensus that it receives from the TPC.

  * •

Each TP passes to the TPC, the calculated optimizers from its own optimization
subproblem. The TPC evaluates
$f_{0}(\tilde{\pi}_{z},\tilde{\mu}_{zs},\tilde{\xi}_{zs})$,
$u^{*}_{k},\hat{\mathbb{P}}^{*}_{sk},\mathbb{P}^{*}_{sh},\phi^{*}_{si(k)},\phi^{*}_{sj(k)},\phi^{*}_{si(h)},\phi^{*}_{sj(h)}$
(which are the minimizers of
$f_{0}(\tilde{\pi}_{z},\tilde{\mu}_{zs},\tilde{\xi}_{zs})$),
$f_{z}(\pi^{\nu(z)},\mu_{s}^{\nu(z)},\xi_{s}^{\nu(z)})$, $P^{*}_{sg}$,
$u^{*}_{zk}$.

  * •

Based on the optimal decision variable values from the last step, the TPC
updates the rewards/penalties according to:

    * $\pi^{\nu(z)+1}=\pi^{\nu(z)}+\frac{\tilde{\alpha}^{\nu}}{\tilde{\beta}_{z}}(u^{*}_{zk}-u^{*}_{k})$.

    * $\mu_{s}^{\nu(z)+1}=\mu_{s}^{\nu(z)}+\frac{\tilde{\alpha}^{\nu}}{\tilde{\beta}_{z}}(|\hat{P}^{*}_{zsk}|-\hat{\mathbb{P}}^{*}_{sk})$.

    * $\mu_{s}^{\nu(z)+1}=\mu_{s}^{\nu(z)}+\frac{\tilde{\alpha}^{\nu}}{\tilde{\beta}_{z}}(|P^{*}_{zsh}|-\mathbb{P}^{*}_{sh})$.

    * $\xi_{s}^{\nu(z)+1}=\xi_{s}^{\nu(z)}+\frac{\tilde{\alpha}^{\nu}}{\tilde{\beta}_{z}}(\theta^{*}_{zsi(k)}-\phi^{*}_{si(k)})$.

    * $\xi_{s}^{\nu(z)+1}=\xi_{s}^{\nu(z)}+\frac{\tilde{\alpha}^{\nu}}{\tilde{\beta}_{z}}(\theta^{*}_{zsj(k)}-\phi^{*}_{sj(k)})$.

    * $\xi_{s}^{\nu(z)+1}=\xi_{s}^{\nu(z)}+\frac{\tilde{\alpha}^{\nu}}{\tilde{\beta}_{z}}(\theta^{*}_{zsi(h)}-\phi^{*}_{si(h)})$.

    * $\xi_{s}^{\nu(z)+1}=\xi_{s}^{\nu(z)}+\frac{\tilde{\alpha}^{\nu}}{\tilde{\beta}_{z}}(\theta^{*}_{zsj(h)}-\phi^{*}_{sj(h)})$.

  * •

Each TP computes a lower bound on its objective function by relaxing the
integer variables and passes on this information to the TPC, which also
computes its lower bound. It then adds up all the regional lower bounds for
the current iteration as follows:
$LB^{new}=LB_{0}(\tilde{\pi}_{z},\tilde{\mu}_{zs},\tilde{\xi}_{zs})+\sum_{z\in
Z}LB_{z}(\pi^{\nu(z)},\mu_{s}^{\nu(z)},\xi_{s}^{(z)})$.

  * •

Each TP updates the above and passes them to the TPC, in response to which,
the TPC does the following:

    * Step 1: updates values of rewards/penalties (dual variables) $\pi^{\nu(z)+1}$, $\pi^{\nu(z)}$, $\mu_{s}^{\nu(z)+1}$, $\mu_{s}^{\nu(z)}$, $\xi_{s}^{\nu(z)+1}$, $\xi_{s}^{\nu(z)}$.

    * Step 2: updates iteration count $\nu(z):=\nu(z)+1$

    * Step 3: calculates new lower bound of objective.

    * Step 4: calculates $f_{0}(\pi^{\nu(z)},\mu_{s}^{\nu(z)},\xi_{s}^{\nu(z)})$

  * •

This step is for the primal recovery. In this step, each TP fixes the right
hand side of the non-anticipativity constraints to the optimal values of the
left-side, like $u_{k}=u^{*}_{zk}$ from the last iteration and evaluates the
optimization problem for only the TPC and its own region. Let’s call the
optimum, for each region as upper bound and indicate it as $UB_{z}$ and for
the TPC, as $UB^{TPC}_{z}$. It then conveys this information to the TPC.

  * •

The TPC determines the global $UB$ as
$UB=\text{min}\\{UB^{TPC}_{z}\\}+\sum_{z\in Z}UB_{z}$

  * •

The TPC terminates the algorithm when $1-(LB/UB)\leq\epsilon$, where
$\epsilon$ is a predecided tolerance and/or, when consensus is reached.

This algorithm is an asynchronous counterpart of a previously presented
algorithm in [26], [28], and [27]. In all these previous references, the
decomposition has been carried out across different scenarios among the shared
variables (which happened to be the decision variables for switching states of
slowly responding generators in the context of solving Stochastic Unit
Commitment problem). Eventually an attempt is made to reach a consensus
regarding the optimal values. Block Coordinate Descent (BCD) is an approximate
method, provided that the error in the subgradient is bounded, and the
stepsize in the algorithm is diminishing [1]. In this situation, the algorithm
converges to an approximate solution of its dual problem if the original
problem is convex (and if not, the dual gives a lower bound). Thus, in the
proposed application, it is expected that when the algorithm converges and the
error is bounded, it leads to the approximate optimum.

###  IV-C Mathematical Model of Stage II Mechanism Design

The second stage of the mechanism is based on “Auxiliary Problem Principle
(APP),” which decides the generation values and the tie-line flows more
accurately (if the accuracy in the first stage for these is not high enough).
Let us first focus on the problem instances (3c) which are the area decomposed
versions of (2) with some modifications. In order to consider the inter-region
transmission expansion, or the horizontal coordination as well as the existing
lines shared between two different regions, it is necessary to calculate the
flows on these existing and potential shared lines. In the previous section
the algorithm presented helps us determine which candidate lines will actually
be considered for constructing new transmission lines as well as determining
the first-order approximation of flows on both the existing lines as well as
the “to-be-constructed lines.” Once we have reached a consensus regarding the
optimal values of the discrete decision variables in the second stage, we
explicitly solve for the continuous variables. Our goal here is to come up
with an incentive mechanism design such that the individual TPs are
incentivized to act in a way that reaches the global optimum of the above-
mentioned problem. In order to achieve this end it is important that each TP
reach consensus about the flows on the shared lines. This is achieved with the
_Auxiliary Problem Principle (APP)_ layer of the incentive mechanism. We
assume each of the shared lines belonging to each region is replaced by a
fictitious generator as shown in Fig. 2. In Fig. 2 we have shown part of a
power system network which is split into two regions on the two sides of the
solid and dashed lines. There are two shared transmission lines, one existing
(the solid line) and one candidate line (the dashed one) shared between the
two regions. These two shared lines are notionally replaced by four
generators—two for each line. The dashed generators represent the candidate
transmission line and the solid generators are the ones for the existing
transmission line. When the original centralized optimization problem is
decomposed across the different regions, the flows on each line are treated
equivalently to the outputs of these fictitious generators. However, each TP
has to have some belief about the value of the shared variables and also about
what the other TP’s beliefs are about those respective variables. The APP
iterations attempt to reach a consensus on those beliefs.

Figure 2: Illustration of the APP layer of the incentive mechanism:
centralized problem

Let us consider a line $h\in\tilde{H}$ shared between two regions, $z_{1}$ and
$z_{2}$. The beliefs regarding the values of the node voltage angles at the
two ends of the lines are as follows:

  * •

$\theta_{z_{1},i(h)}$: Belief of $z_{1}$ about the voltage phase angle of the
node belonging to itself

  * •

$\theta_{z_{1},j(h)}$: Belief of $z_{1}$ about the voltage phase angle of the
node belonging to $z_{2}$

  * •

$\theta_{z_{2},i(h)}$: Belief of $z_{2}$ about the voltage phase angle of the
node belonging to itself

  * •

$\theta_{z_{2},j(h)}$: Belief of $z_{2}$ about the voltage phase angle of the
node belonging to $z_{1}$

Ideally $\theta_{z_{1},i(h)}=\theta_{z_{2},j(h)}$ and
$\theta_{z_{2},i(h)}=\theta_{z_{1},j(h)}$ must be satisfied $\forall
h\in\tilde{H}$. However, this need not necessarily be true at the end of the
algorithmic incentive mechanism computations introduced in Section IV. The
optimization problem formulations stated below describes the _APP_ iterations
to achieve the consensus between the beliefs regarding the voltage phase
angles and consequently the flows on the shared lines.

$\displaystyle\min_{P_{sg},\theta}\Bigg{(}\sum_{s\in
S}\sum_{g\in{G_{z}}}w_{s}C_{g}(P_{sg})$
$\displaystyle+\sum_{z^{\prime}}\sum_{h\in\tilde{H}_{z}}\Big{\\{}\eta\Big{(}\theta_{z,i(h)}(\theta^{\sigma}_{z,i(h)}-\theta^{\sigma}_{z^{\prime},j(h)})$
$\displaystyle+\theta_{z,j(h)}(\theta^{\sigma}_{z,j(h)}-\theta^{\sigma}_{z^{\prime},i(h)})\Big{)}+\frac{\gamma}{2}\Big{(}||\theta_{z,i(h)}-\theta^{\sigma}_{z,i(h)}||^{2}_{2}$
$\displaystyle+||\theta_{z,j(h)}-\theta^{\sigma}_{z,j(h)}||^{2}_{2}\Big{)}+\lambda^{\sigma}_{zz^{\prime}i(h)}\theta_{z,i(h)}+\lambda^{\sigma}_{zz^{\prime}j(h)}\theta_{z,j(h)}\Big{\\}}\Bigg{)}$
(4a) $\displaystyle\mbox{Subject\;to:\>}\sum_{g\in G_{zn}}P_{sg}-\sum_{d\in
D_{zn}}P_{sd}=$ $\displaystyle\sum_{h\in
H_{zn}}\mathbf{A}_{(n,h)}P_{sh}+\sum_{h\in{\tilde{H}_{zn}}}P_{sh}$
$\displaystyle\forall{z\in Z},\forall{s\in S},\forall{{n}\in\mathcal{N}_{z}}$
(4b) $\displaystyle
P_{sh}=\frac{1}{x_{h}}\sum_{n\in\mathcal{N}}\mathbf{A}_{(n,h)}\theta_{sn},\forall{z\in
Z},\forall{s\in S},\forall{{h}\in H_{z}}$ (4c) $\displaystyle
P_{sh}=\frac{1}{x_{h}}(\theta_{szi(h)}-\theta_{szj(h)})$ $\displaystyle
h\leftarrow(i,j),\forall{z\in Z},\forall{s\in S},\forall{{h}\in\tilde{H}_{z}}$
(4d) $\displaystyle-{\overline{L}}_{h}\leq
P_{sh}\leq{{\overline{L}}_{h}},\;\forall{z\in Z},\forall{s\in
S},\forall{{h}\in H_{z}}$
$\displaystyle-{\overline{L}}_{h}\leq{P}_{sh}\leq{\overline{L}}_{h},\;\forall{z\in
Z},\forall{s\in S},\forall{{h}\in\tilde{H}_{z}}$
$\displaystyle{\underline{P}}_{sg}\leq
P_{sg}\leq{{\overline{P}}_{sg}},\;\forall{z\in Z},\forall{s\in
S},\forall{g\in{G_{z}}}$ (4g) Dual Variable Updates
$\displaystyle\lambda^{\sigma+1}_{zz^{\prime}i(h)}=\lambda^{\sigma}_{zz^{\prime}i(h)}+\delta_{APP}(\theta^{\sigma+1}_{z,i(h)}-\theta^{\sigma+1}_{z^{\prime},j(h)})\;\forall
h\in\tilde{H}$ (4h)
$\displaystyle\lambda^{\sigma+1}_{zz^{\prime}j(h)}=\lambda^{\sigma}_{zz^{\prime}j(h)}+\delta_{APP}(\theta^{\sigma+1}_{z,j(h)}-\theta^{\sigma+1}_{z^{\prime},i(h)})\;\forall
h\in\tilde{H}$ (4i)

###  IV-D Steps of Stage II Mechanism Design

In the incentive mechanism design we assume the presence of a TPC as before.
For example, in Europe the TPC is ENTSO-E. In North America, the closest
entity to a TPC is FERC. This entity sets the penalty and designs the market
rules, in order to incentivize the individual owners of transmission sub-
network, to act in a way that solves the problem (4). In the following section
we will explicitly state the steps or the rules of the second stage of the
incentive mechanism design. In this stage each TP solves for the continuous
variables more precisely. The steps of this stage of the mechanism can be
described as follows:

  * •

Each region $z\in Z$ receives the values of the previous beliefs of voltage
phase angles of the other node of the shared lines belonging to the adjoining
regions, from the TPC. Regions $z$ also receive the previous updates of the
rewards/penalties, which are the Lagrange Multiplier or dual variable values
($\theta^{\sigma}_{z^{\prime},i(h)},\theta^{\sigma}_{z^{\prime},j(h)}$ in
(4)).

  * •

Each region independently solves the optimization problem (4a)–(4g) and
broadcasts the optimal decision variable values of its own voltage phase
angles of the nodes of the shared lines to the TPC.

  * •

The TPC calculates the rewards/penalties (Lagrange multipliers) by solving
constraints (4h)–(4i) and sends the updated values to the respective regions.

  * •

The TPC stops the process when
$(\theta^{\sigma+1}_{z,i(h)}-\theta^{\sigma+1}_{z^{\prime},j(h)})^{2}+(\theta^{\sigma+1}_{z,j(h)}-\theta^{\sigma+1}_{z^{\prime},i(h)})^{2}\leq\epsilon_{APP}$
for all the pairs of regions and all the shared lines, where $\epsilon_{APP}$
is the predetermined APP tolerance.

###  IV-E The Proposed Organizational Set-up

In summary, this paper suggests the following set-up:

  * •

Step 1: the TPC announces the incentive mechanism for stage I and asks the TPs
to submit their candidate lines for inter- and intra-regional planning.

  * •

Step 2: the TPC solves its optimization problem in stage I and sends the
penalty and rewards to TPs.

  * •

Step 3: TPs correct their decisions based on the penalty and rewards they have
received. The new decisions are sent to TPC

  * •

Step 4: the TPC iterates between Steps 2 and 3 to reach a consensus.

  * •

Step 5: Once the consensus is reached the inter- and intra-regional planning
results will be published. The TPC announces the incentive mechanism for stage
II and the same steps as stage I are run for the second stage. Finally, with
the final results, TPs can start constructing the planned lines.

##  V Numerical Results

###  V-A 2-region System

In this sub-section, an illustrative two-region power system is considered to
present our novel mechanism design compared to a game-theoretic design. For
the system shown in Figure 3, Tables I and II provide the system details.

Figure 3: Two-region power system for comparison between game-theoretic and mechanism design results TABLE I: The data of the two-region system Parameter | Region-1 | Region-2
---|---|---
Gen. Range (MW) | 0-3000 | 0-3000
Gen. Cost ($/MWh) | 50 for 0-1800 MW | 10
| 200 for 1800-300 MW |
Load (MW) | 2000 | 500
First-case Investment ($) | 1000 | 1000
second-case Investment ($) | 20000 | 20000
TABLE II: Shared Lines of two-Region System Parameter | Existing | Candidate
---|---|---
Capacity (MW) | 150 | 1350
Reactance (pu) | 9X | X

The line along the candidate line is constructed iff both TP-1 and TP-2 think
that the line should be constructed. We have considered here two cases. In the
first case, the total investment cost is $2000 and is shared equally between
TP-1 and TP-2 if the line is constructed.

Figure 4: Game matrix for two-region power system

The two-player cost matrix for the game of strategic interaction between TP-1
and TP-2 in this case is shown at the top of Fig. 4. The numbers shown are the
costs faced by the individual regions. The numbers 1 and 0 represent the
decisions to build and to not build the transmission line by each of the
players, respectively. The social cost for the (1,1) entry is $47000, whereas
for each of the other three boxes is $378500. Hence, the box that is in orange
(corresponding to the transmission line being built) is clearly the social
optimum (SO, or Pareto optimum). The game-theoretic Nash Equilibrium (NE) are
the boxes in green, and we can see that there are two of them.

In the second case, we consider the situation where the cost for building the
new line is $40000 which is shared equally among the two players. If one
player decides to build the line when another one does not, then the player
who decides to build the line bribes the other one in the hope of changing
their mind in the future. The cost matrix for such a game is shown at the
bottom half of Fig. 4. Here, the SO is still the box (1,1), whereas now there
is only one NE in box (0,0). In both the cases examined, the SO and NE are not
the same strategy combination. This is a serious drawback of the game-
theoretic interaction, which our mechanism attempts to resolve.

In order to show how our mechanism can help resolve the diffrence between SO
and NE, we now refer to the modified model of the system shown in the bottom
half of Fig. 3. In this diagram, we have introduced the TPC with the green
lines representing the message exchange pertaining to the rewards/penalties
and tentative decisions for building the lines. The shared lines have been
split up and represented as virtual generators, with the line capacities for
forward and reverse flows acting as generating limits and the investment cost
can be interpreted as the start-up cost of these fictitious generators. The
mechanism is implemented in different steps or iterations. Before the start of
each round or iteration of the mechanism, the TPC broadcasts the most recent
estimate of the rewards or penalties ($\pi,\mu_{k},\mu_{h}$ respectively refer
to the rewards/penalties associated with the non-anticipativity constraints
for consensus for decision to build transmission line, flow on the candidate
line, and flow on the existing shared line). Each region then solves its own
optimization problem (calculating $u^{*},P^{*}_{k},P^{*}_{h},P^{*}_{g}$, which
respectively refer to the integer variable for decision to build transmission
line, flow on the candidate line, flow on the existing shared line, and the
generator output) and relays its tentative decision to the TPC, following
which the TPC updates the rewards/penalties. This goes on until consensus is
reached. At the beginning of the first iteration, all the rewards/penalties
are set to zero. The details of each step are shown in Table III:

TABLE III: Iterates of the incentive mechanism Iterates | Iter-1 | Iter-2 | Iter-3 | Iter-4
---|---|---|---|---
$u^{TP1*}$ | 1 | 1 | 1 | 1
$P^{TP1*}_{k}$ | 1350 | 180 | 180 | 1350
$P^{TP1*}_{h}$ | 150 | 20 | 20 | 150
$P^{TP1*}_{g}$ | 500 | 1800 | 1800 | 500
$u^{TP2*}$ | 0 | 0 | 1 | 1
$P^{TP2*}_{k}$ | 0 | 0 | -1350 | -1350
$P^{TP2*}_{h}$ | 0 | 0 | -150 | -150
$P^{TP2*}_{g}$ | 500 | 500 | 2000 | 2000
$u^{MO*}$ | 0 | 1 | 1 | 1
$P^{MO*}_{k}$ | 0 | 1350 | 1350 | 1350
$P^{MO*}_{h}$ | 0 | 150 | 150 | 150
$\pi^{TP1*}$ | 1 | 1 | 1 | 1
$\mu^{TP1*}_{k}$ | 1350 | 180 | -990 | -990
$\mu^{TP1*}_{h}$ | 150 | 20 | -110 | -110
$\pi^{TP2*}$ | 0 | -1 | -1 | -1
$\mu^{TP2*}_{k}$ | 0 | 1350 | 1350 | 1350
$\mu^{TP2*}_{h}$ | 0 | 150 | 150 | 150
$(1-UB^{*}/LB^{*})$ | -0.01 | 15.78 | 0.0002 | 0

As we can see from Table III, indeed our proposed mechanism reaches the social
optimum and the finally settled values of the rewards/penalties are precisely
the incentivizing factors needed to force the behavior of each TP away from
the NE and towards the SO.

###  V-B Generalized/Multi-region System

In Fig. 5, we have considered a bigger system, which consists of three
regions. Regions 1, 2, and 3, respectively are the IEEE 14-, 30-, and 5-node
systems. The red lines are the inter-regional shared candidate lines, the blue
lines represent the intra-regional candidate lines, and the black ones are the
existing shared lines.

Figure 5: The multi-regional system considered for the simulation Figure 6:
Variation of convergence criteria for the first stage according to iterations

We apply the stage-I incentive mechanism design algorithm for 1000 iterations
in order to illustrate the trend and variations of the convergence curve. We
get the convergence at the 85th iteration as illustrated in Fig. 6. The
convergence curve indicates the difference between the low and high bounds of
the distributed problem. Fig. 6 explains that the convergence curve gradually
reaches a stable value after 500 iterations. The big spikes in the convergence
curve, for example, the ones happening at iterations 123, 195, 335, and 456
indicate the discrepancy between the binary decision of building the fourth
line. They state that the decision of the TPC does not match those of TPs. The
smaller oscillations refer to a minor disagreement in the operational
decisions. They can be disregarded since their domain is relatively small and
the aim of the first stage is to determine binary decision variables related
to the candidate lines’ constructions.

Figure 7: Variation of convergence criteria for the second stage according to
iterations

According to the given system, Region 2 has a high load, and expensive
generation whereas Regions 1 and 3 have low and moderate loads with cheap and
moderately cheap generations, respectively. We have observed that in this
situation, the optimal decision is to build candidate line 4 which is shared
between Regions 1 and 2.
In addition, Fig. 7 shows the variation of the convergence curve for the
second stage. We have initialized variables to zeros. Hence the convergence
starts with a small value. After almost 1000 iterations, all TPs reach an
agreement on the operational variables. The obtained convergence criterion is
equal to 0.012.
Finally, Table IV compares the intra-regional candidate decisions considering
centralized and the proposed distributed models. Correspondingly, Table V
compares the inter-regional candidate decisions when the results are
determined by the centralized and the proposed distributed methods.
We should note that the centralized framework is just a mathematical benchmark
in our paper and it does not exist in the real-world situations. In real-life
situations, some level of coordination always exists. In our paper, we suggest
an optimization-based coordination for the investment coordination problem. We
have shown that our mechanism can reach a close-to-global solution as compared
to our theoretical benchmark model. At our close-to-optimal solution, each TP
receives a share of total social welfare which is optimal for the whole inter-
connected system. Hence, it will achieve a solution close to the original
optimal solution. The error is mainly due to the non-convex nature of the
whole optimization problem.
We have run this simulation on a desktop computer using Julia programming
language with JuMP optimization interface and with several different solvers
(for the first stage we have used MILP solvers such as Gurobi, HiGHS, GLPK,
and for the second stage which is a Quadratic Program (QP) we have used
solvers such as Gurobi and IPOPT). Codes are available to be run and verified
at http://github.com/sambuddhac/Horizontal_Proper.

TABLE IV: Intra-regional Candidate Line Construction Decisions Mechanism | Region-1 | Region-2 | Region-3
---|---|---|---
Centralized | 10 | 1,10,11,14,16 | –
Dist. Mech. | – | – | –
TABLE V: Inter-regional Candidate Line Construction Decisions Mechanism | Region-1 | Region-2 | Region-3 | TPC | Decision
---|---|---|---|---|---
Centralized | 4 | 4,5 | 5 | 4,5 | 4,5
Dist. Mech. | 4 | 4 | – | 4 | 4

##  VI Conclusion

In this paper, we have addressed the important problem of transmission
investment when there are multiple regions and new lines built between them.
We have proposed a two-staged incentive mechanism design consisting of a
region-decomposed distributed Lagrangian method. In our numerical study, we
have analyzed an illustrative two-region network and a larger and more general
three-region system. The two-region network results intuitively analyzed and
validated our incentive mechanism approach, which includes a TPC, compared to
the game-theoretic Nash equilibrium approach that does not have a TPC [31, 32,
36]. We showed that our proposed incentive mechanism indeed reached the Pareto
optimal solution which is also the socially optimal solution. For the three-
region case, we implemented our proposed mechanism design in a computational
framework to simulate the mechanism performance. We presented the convergence
characteristics of both the stages of our mechanism and discussed the results.
Future research steps include:

  * •

integrating our model with generation expansion and consider policies such as
capacity reserve margin, renewable portfolio standard and emissions margin.

  * •

including $(N-1)$ contingency criteria within our transmission expansion
planning problem.

  * •

including temporal variation, inter-temporal constraints, and time-domain
reduction.

As a limitation of implementing the proposed algorithms, one challenge would
be related to defining the appropriate stepsize in the first and second
algorithms in order to help the convergence of the algorithm and avoid
unbounded error as well as a compromise between the computational costs and
accuracy of the final solutions. This issue can be further analyzed and
resolved in future research.

## References

  * [1] S I A Aravena, A Papavasiliou and A Papalexopoulos  “A Distributed Computing Architecture for the Large-Scale Integration of Renewable Energy and Distributed Resources in Smart Grids”, 2017
  * [2] R Baldick, B H Kim, C Chase and Y Luo  “A Fast Distributed Implementation of Optimal Power Flow”  In _IEEE Transactions on Power Systems_ 14.3 IEEE, 1999, pp. 858–864
  * [3] Tom Brown, Jonas Hörsch and David Schlachtberger  “PyPSA: Python for Power System Analysis”  In _arXiv preprint arXiv:1707.09913_ , 2017
  * [4] S Chakrabarti  “Post-Contingency States Representation & Redispatch for Restoration in Power Systems Operation”, 2017
  * [5] S Chakrabarti  “Study of UPLAN based Resources Planning & Analysis by Power Generation Utilities in the Deregulated Electricity Market”, 2010
  * [6] Sambuddha Chakrabarti and Ross Baldick  “Look-Ahead SCOPF (LASCOPF) for Tracking Demand Variation via Auxiliary Proximal Message Passing (APMP) Algorithm”  In _arXiv_ 1 arXiv, 2019
  * [7] Sambuddha Chakrabarti and Ross Baldick  “Look-Ahead SCOPF (LASCOPF) for Tracking Demand Variation via Auxiliary Proximal Message Passing (APMP) Algorithm”  In _International Journal of Electrical Power & Energy Systems_ 116 Elsevier, 2020, pp. 105533
  * [8] G Cohen  “Auxiliary Problem Principle and Decomposition of Optimization Problems”  In _Journal of optimization Theory and Applications_ 32.3 Springer, 1980, pp. 277–305
  * [9] G Cohen  “Optimization by Decomposition and Coordination: A Unified Approach”  In _IEEE Transactions on automatic control_ 23.2 IEEE, 1978, pp. 222–232
  * [10] R Ebrahimian and R Baldick  “State Estimation Distributed Processing [for power systems]”  In _IEEE Transactions on Power Systems_ 15.4 IEEE, 2000, pp. 1240–1246
  * [11] US EIA  “The Electricity Market Module of the National Energy Modeling System: Model Documentation 2020”  In _US Energy Information Administration, Washington, DC, Tech. Rep_ , 2020
  * [12] EPA  “Documentation for EPA’s Power Sector Modeling Platform v6-November 2018 Reference Case”  US Environmental Protection Agency Washington, DC, 2018
  * [13] Jonathan Ho et al.  “Regional Energy Deployment System (ReEDS) Model Documentation: Version 2020”, 2021
  * [14] Wenjing Hou et al.  “Incentive-driven task allocation for collaborative edge computing in industrial internet of things”  In _IEEE Internet of Things Journal_ 9.1 IEEE, 2021, pp. 706–718
  * [15] Jesse D Jenkins and Nestor A Sepulveda  “Enhanced Decision Support for a Changing Electricity Landscape: The GenX Configurable Electricity Resource Capacity Expansion Model”  MIT Energy Initiative, 2017
  * [16] Jesse David Jenkins  “Electricity System Planning with Distributed Energy Resources: New Methods and Insights for Economics, Regulation, and Policy”, 2018
  * [17] Josiah Johnston, Rodrigo Henriquez-Auba, Benjamín Maluenda and Matthias Fripp  “Switch 2.0: A Modern Platform for Planning High-Renewable Power Systems”  In _SoftwareX_ 10 Elsevier, 2019, pp. 100251
  * [18] B H Kim and R Baldick  “A Comparison of Distributed Optimal Power Flow Algorithms”  In _IEEE Transactions on Power Systems_ 15.2 IEEE, 2000, pp. 599–604
  * [19] B H Kim and R Baldick  “Coarse-Grained Distributed Optimal Power Flow”  In _IEEE Transactions on Power Systems_ 12.2 IEEE, 1997, pp. 932–939
  * [20] S Magnússon  “Bandwidth Limited Distributed Optimization with Applications to Networked Cyberphysical Systems”, 2017
  * [21] Trieu Mai et al.  “Resource Planning Model: An Integrated Resource Planning and Dispatch Tool for Regional Electric Systems”, 2013
  * [22] Mohammad Majidi-Qadikolai and R. Baldick  “A Generalized Decomposition Framework for Large-Scale Transmission Expansion Planning”  In _IEEE Transactions on Power Systems_ PP, 2017, pp. 1–1  DOI: 10.1109/TPWRS.2017.2724554
  * [23] Mohammad Majidi-Qadikolai and R. Baldick  “Integration of $N-1$ Contingency Analysis With Systematic Transmission Capacity Expansion Planning: ERCOT Case Study”  In _IEEE Transactions on Power Systems_ 31, 2015, pp. 1–12  DOI: 10.1109/TPWRS.2015.2443101
  * [24] Mohammad Majidi-Qadikolai and R. Baldick  “Reducing the number of candidate lines for high level Transmission Capacity Expansion Planning under uncertainties”, 2015, pp. 1–6  DOI: 10.1109/NAPS.2015.7335117
  * [25] Mohammad Majidi-Qadikolai and R. Baldick  “Stochastic Transmission Capacity Expansion Planning With Special Scenario Selection for Integrating N-1 Contingency Analysis”  In _IEEE Transactions on Power Systems_ 31, 2016, pp. 1–12  DOI: 10.1109/TPWRS.2016.2523998
  * [26] A Papavasiliou and S S Oren  “Multiarea Stochastic Unit Commitment for High Wind Penetration in a Transmission Constrained Network”  In _Operations Research_ 61.3 INFORMS, 2013, pp. 578–592
  * [27] A Papavasiliou, S S Oren and R P O’Neill  “Reserve Requirements for Wind Power Integration: A Scenario-based Stochastic Programming Framework”  In _IEEE Transactions on Power Systems_ 26.4 IEEE, 2011, pp. 2197–2206
  * [28] A Papavasiliou, S S Oren and B Rountree  “Applying High Performance Computing to Transmission-Constrained Stochastic Unit Commitment for Renewable Energy Integration”  In _IEEE Transactions on Power Systems_ 30.3 IEEE, 2015, pp. 1109–1120
  * [29] Y Tohidi  “Optimal Long-Term Generation-Transmission Planning in the Context of Multiple TSOs”, 2016
  * [30] Y Tohidi and M R Hesamzadeh  “Free Riding Effect in Multi-National Transmission Expansion Planning”  In _Innovative Smart Grid Technologies Europe (ISGT EUROPE), 2013 4th IEEE/PES_ , 2013, pp. 1–5  IEEE
  * [31] Y Tohidi and M R Hesamzadeh  “Multi-National Transmission Planning Using Joint and Disjoint Solutions”  In _European Energy Market (EEM), 2013 10th International Conference on the_ , 2013, pp. 1–6  IEEE
  * [32] Y Tohidi and M R Hesamzadeh  “Multi-Regional Transmission Planning as a Non-Cooperative Decision-Making”  In _IEEE Transactions on Power Systems_ 29.6 IEEE, 2014, pp. 2662–2671
  * [33] Y Tohidi and M R Hesamzadeh  “Multi-Regional Transmission Planning Under Interdependent Wind Uncertainty”  In _Energy Conference (ENERGYCON), 2014 IEEE International_ , 2014, pp. 1474–1479  IEEE
  * [34] Y Tohidi, M R Hesamzadeh and K Ostman  “Reactive Coordination of Transmission-Generation Investment Planning”  In _European Energy Market (EEM), 2015 12th International Conference on the_ , 2015, pp. 1–5  IEEE
  * [35] Y Tohidi, M R Hesamzadeh and F Regairaz  “Sequential Coordination of Transmission Expansion Planning with Strategic Generation Investments”  In _IEEE Transactions on Power Systems_ 32.4 IEEE, 2017, pp. 2521–2534
  * [36] Y Tohidi, L Olmos, M Rivier and M R Hesamzadeh  “Coordination of Generation and Transmission Development Through Generation Transmission Charges—A Game Theoretical Approach”  In _IEEE Transactions on Power Systems_ 32.2 IEEE, 2017, pp. 1103–1114
  * [37] Kun Wang et al.  “A survey on energy internet: Architecture, approach, and emerging technologies”  In _IEEE systems journal_ 12.3 IEEE, 2017, pp. 2403–2416
  * [38] D Young et al.  “US-REGEN Model Documentation”, 2018
  * [39] Feng Zeng, Yaojia Chen, Lan Yao and Jinsong Wu  “A novel reputation incentive mechanism and game theory analysis for service caching in software-defined vehicle edge computing”  In _Peer-to-Peer Networking and Applications_ 14 Springer, 2021, pp. 467–481

|  Sambuddha Chakrabarti (S’14-M’22) is an energy-systems researcher with the
Andlinger Center for Energy + the Environment at Princeton University. From
2017 to 2020 he was a postdoctoral researcher at National Renewable Energy
Laboratory (NREL), KTH, Linkoping University, Arizona State University, and
The University of Texas at Austin. He received his MS and Ph.D. from the
Department of Electrical and Computer Engineering at The University of Texas
at Austin. His research focuses on application of Operations Research methods
to power network planning and operations problems
---|---
|  Hosna Khajeh (S’10) received the M.Sc. (Tech.) degree in electrical
engineering (power systems) from Semnan University, Semnan, Iran in 2016. She
is currently pursuing a doctoral degree at the University of Vaasa, Vaasa,
Finland. Her research interests include electricity and ancillary service
markets, energy flexibility management, forecasting, and trading structures.
---|---
|  Thomas R Nudell received the B.A. magna cum laude degree in physics from
Kalamazoo College in 2009, the M.S. and Ph.D. degrees from North Carolina
State University in 2011 and 2015, respectively, both in electrical
engineering. He was then a postdoctoral research associate with the Active
Adaptive Control Laboratory, Massachusetts Institute of Technology, where he
investigated resilient interconnected infrastructures and transactive control
for microgrid energy management. His professional goals include realizing the
rapid and widespread integration of renewable and energy efficient resources
in the electric power grid. Dr Nudell is the cofounder and CEO of Piq Energy
Corp, building autonomous grid planning and operation software.
---|---
|  Mohammad Reza Hesamzadeh (IEEE-SM’13) received his Docent from KTH Royal
Insitute of Technology, Sweden, and his PhD from Swinburne University of
Technology, Australia, in 2013 and 2010, respectively. He was a Post-Doctoral
Fellow at KTH in 2010-2011 where he is currently a faculty member. He is also
a Faculty Affiliate at Program on Energy and Sustainable development (PESD),
Stanford University and Research Affiliate at German Institute for Economic
Research (DIW Berlin). Dr Hesamzadeh is a senior memeber of IEEE, and a member
of informs, IAEE and Cigre. He has won several awards for his papers in
different professional events and conferences. He has published two books and
numerous papers on electricity market design and analysis. He also served as
editor of IEEE Transactions on Power Systems and Guest Editor of several other
journals. Dr. Hesamzadeh has been providing advice and consulting services on
different energy market issues to both private and government sections over
the last several years.
---|---
|  Ross Baldick (S’90-M’91-SM’04-F’07) recieved the B.Sc. degree in
Mathematics and Physics and the B.E. degree in electrical engineering from the
University of Sydney, Sydney Australia and the M.S. and Ph.D. degrees in
electrical engineering from the University of California, Berkeley, in 1988
and 1990 respectively. From 1991 to 1992, he was a Post-Doctoral fellow with
the Lawrence Berkeley Laboratory, Berkeley, CA. In 1992 and 1993, he was an
Assistant Professor with the Worcester Polytechnic Institute, Worcester, MA.
He is currently a Professor Emeritus with the Department of Electrical and
Computer Engineering, The University of Texas at Austin, Austin, TX.
---|---


# Transport signatures of phase fluctuations in superconducting qubits

M. Wisne  Department of Physics and Astronomy, Northwestern University, 2145
Sheridan Road, Evanston, IL 60208, USA  Y. Deng  Department of Physics and
Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208,
USA  H. Cansizoglu  Rigetti Computing, Berkeley, CA 94710, USA  C. Kopas
Rigetti Computing, Berkeley, CA 94710, USA  J. Y. Mutus  Rigetti Computing,
Berkeley, CA 94710, USA  V. Chandrasekhar  Department of Physics and
Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208,
USA

###### Abstract

Josephson junctions supply the nonlinear inductance element in superconducting
qubits. In the widely used transmon configuration, where the junction is
shunted by a large capacitor, the low charging energy minimizes the
sensitivity of the qubit to charge noise while maintaining the necessary
anharmonicity to qubit states. We report here low-frequency transport
measurements on small standalone junctions and identically fabricated
capacitively-shunted junctions that show two distinct features normally
attributed to small capacitance junctions near zero bias: reduced switching
currents and prominent finite resistance associated with phase diffusion in
the current-voltage characteristic. Our transport data reveals the existence
of phase fluctuations in transmons arising from intrinsic junction
capacitance.

††preprint: APS/123-QED

Josephson junctions (JJs) form the heart of superconducting qubits, arguably
the leading technology for quantum computing and quantum sensing. In its
simplest form, a JJ is a tunnel junction between two superconductors. Two
energy scales govern the physics of such a junction: the Josephson energy
$E_{J}$, which is proportional to the superconducting critical current $I_{c}$
of the junction, and the charging energy $E_{C}$, inversely proportional to
the intrinsic capacitance $C$ of the junction. JJs have been studied now for
more than four decades, first as essential elements of superconducting
electronics and sensors such as superconducting quantum interference devices
(SQUIDs) as well as a fundamental physics playground for macroscopic quantum
effects [1, 2, 3, 4, 5]. Adding to this rich history, JJs now form the
critical component of superconducting quantum circuits.

The first JJs were many tens of square microns in area, with correspondingly
large critical currents $I_{c}$ and therefore large $E_{J}$. These JJs had
relatively large intrinsic capacitances $C$ between electrodes and therefore
low $E_{C}$, and thus were in the regime $E_{C}/E_{J}\ll 1$. On the other
hand, JJs in current manifestations of superconducting qubits are submicron in
size, and $E_{C}$ is more comparable to $E_{J}$. Such junctions are the focus
of this paper and are referred to as small-$C$ JJs from this point forward. It
has been known for some years that small-$C$ junctions are particularly
sensitive to noise [6, 7, 8]. In order to decrease the susceptibility to
charge noise in superconducting qubits, most qubits are fabricated in the so-
called transmon configuration, where the JJ is shunted by an on-chip capacitor
$C_{S}\gg C$. Such transmon qubits show significant improvement in coherence
times over unshunted JJs [9, 10, 11].

Here we report dc and low frequency electrical transport measurements on
small-$C$ unshunted (“standalone”) JJs with $E_{C}/E_{J}\geq 1$ as well as
shunted (“transmon”) JJs where the charging energy is several times reduced.
Somewhat surprisingly, we find that all devices show two features associated
with small capacitance junctions subject to noise in their environment: a
finite zero-bias resistance in the nominally superconducting regime of the JJ
current-voltage characteristic (IVC), and a switching current $I_{sw}$ out of
the superconducting regime that is much lower than the expected critical
current $I_{c}$ based on the measured characteristics of the JJ. Most notably,
we find that there is no significant difference between the dc transport
characteristics of unshunted and shunted JJs with nominally identical junction
areas and fabrication methods, contrary to theoretical predictions. dc
transport characterization on superconducting qubits is not routinely
performed: our results suggest that such measurements may be useful in
optimizing qubit fabrication to increase coherence times.

##  I Background

The dynamics of a JJ can be described by the quantum Langevin equation for a
particle of mass $m=(\hbar/2e)^{2}C$ in a potential $U(\phi)$ subject to a
time- ($t$) dependent random force $F(t)$ [12]

$m\ddot{\phi}+\int_{0}^{\infty}\mu(t-t^{\prime})\phi(t^{\prime})dt^{\prime}+\frac{\partial
U}{\partial\phi}=F(t).$ (1)

where $\phi$ is the phase difference between the two superconductors forming
the JJ. $\mu(t)$ is the so-called memory function related to the
autocorrelation of the random force $F(t)$ [12]

$\displaystyle\frac{1}{2}$
$\displaystyle<F(t)F(t^{\prime})+F(t^{\prime})F(t)>$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{0}^{\infty}d\omega\;\Re[\tilde{\mu}(\omega)]\hbar\omega\coth(\hbar\omega/2k_{B}T)\cos\omega(t-t^{\prime}),$
(2)

where $\tilde{\mu}(\omega)$ is the Fourier transform of $\mu(t)$, $k_{B}$ is
Boltzmann’s constant and $T$ is the temperature. Equation 1 is identical to
the equation for the Brownian motion of a particle of mass $m$ in a potential
$U$ subject to the random force $F(t)$ [13]. This equation is a generalized
form of the more familiar resistively and capacitively shunted junction (RCSJ)
model for a JJ, where $U(\phi)$ is the tilted washboard potential
$U(\phi)=-E_{J}\cos(\phi)-I\phi$, $E_{J}=(\hbar/2e)I_{c}$ being the Josephson
energy and $I$ the current applied through the junction [14]. $I_{c}$ is given
by the Ambegaokar-Baratoff relation [15]

$I_{c}R_{n}=\frac{\pi\Delta(T)}{2e}\tanh(\Delta(T)/2k_{B}T)$ (3)

where $R_{n}$ is the normal state resistance of the junction and $\Delta(T)$
the temperature-dependent gap. The additional term $F(t)$ which is absent in
the conventional RCSJ model represents the noise sources influencing the JJ.
Classically, one expects a particle with smaller mass to be more affected by a
given random force $F(t)$. Indeed, since the mass of the fictitious phase
particle is proportional to the intrinsic capacitance $C$, junctions with
smaller $C$ tend to be more susceptible to noise.

Let us consider some simple limits of Eqn. 1 to gain a better understanding of
the dynamics of our phase particle and to define some parameters of the
problem. If we ignore $F(t)$ and its resulting dissipative term (the second
term with the integral in Eqn. 1) and consider the case with $I=0$, we can
expand the wasboard potential about its minimum at $\phi=0$ to obtain

$m\ddot{\phi}+\frac{1}{2}m\omega_{0}^{2}\phi^{2}=0$ (4)

which is the equation of a simple harmonic oscillator with the fundamental
oscillation frequency given by
$\hbar\omega_{0}=\sqrt{4e^{2}E_{J}/C}=\sqrt{8E_{C}E_{J}}$, where
$E_{C}=e^{2}/2C$ is the single electron charging energy. The quantum analog of
this equation relevant for superconducting qubits corresponds to a harmonic
oscillator with energy levels $E_{n}=(n+1/2)\hbar\omega_{0}$. Considering the
full anharmonic cosine potential results in the unequal level spacings
required for qubit operation. Note that to first order, $\omega_{0}$ is
independent of the size of the junction; $I_{c}$ and $C$ scale with the area
of the junction, so that the product $E_{J}E_{C}$ is independent of the area
of the junction, while $E_{C}/E_{J}$ increases as the junction size decreases.

To include the effects of the random force and dissipation, it is useful to
Fourier transform Eqn. 1. The response of the system can then be represented
by $\tilde{\phi}(\omega)=\alpha(\omega)\tilde{F}(\omega)$, where
$\tilde{\phi}$ and $\tilde{F}$ are the Fourier transforms of $\phi(t)$ and
$F(t)$, respectively, and the response function $\alpha(\omega)$ is given by
[12]

$\alpha(\omega)=\frac{1}{-m\omega^{2}+i\omega\tilde{\mu}(\omega)+m\omega_{0}^{2}}.$
(5)

The poles of $\alpha(\omega)$ give the excitation energies of the system. If
$\tilde{\mu}(\omega)=0$, we regain our earlier result for the excitation
energy. It would seem that including a finite $\tilde{\mu}(\omega)$ would
introduce a decay time for the oscillation, or a finite lifetime for the state
in the quantum limit. However, this depends on the specific form of
$\tilde{\mu}(\omega)$. In our case, $\tilde{\mu}(\omega)$ can be related to
the equivalent admittance of the environment coupled to the JJ,
$\tilde{\mu}(\omega)=(\hbar/2e)^{2}Y(\omega)$ [16], where $Y(\omega)$ is the
equivalent frequency dependent admittance considered connected in parallel to
the junction. Let us now consider the transmon configuration, where the JJ is
shunted by a large capacitor $C_{S}$, $Y(\omega)=i\omega C_{S}$. In this case,
it is easy to show that the excitation energy is then modified to
$\hbar\omega_{0}=\sqrt{8E_{J}E_{C}^{\prime}}$, where
$E_{C}^{\prime}=e^{2}/2(C+C_{S})$, as expected.

More generally, the environmental coupling to the JJ, which includes the
impedances of the measurement electronics, the wiring in the cryostat, etc.,
may include an ohmic component that will give rise to dissipation. In the
context of the RCSJ model, dissipation is embedded in a shunt resistance $R$
in parallel to the JJ, so that $Y(\omega)=1/R$. While the exact origin of $R$
seems to be a matter of debate, it is commonly thought that $R$ is the
resistance associated with quasiparticle transport across the junction [17]
combined with any ohmic contributions from the environment such as the source
impedance of the electronics connected to the JJ [6, 16, 18]. Associated with
$R$ are fluctuations in the current $\Delta I(t)$ attributed to Johnson noise
and are described by the fluctuating random force $F(t)$ in the Langevin
equation (Eqn. 1).

In principle, the Josephson phase particle, treated as a quantum variable, can
thermally activate or tunnel through the barrier between adjacent minima of
the tilted cosine potential. This is more likely for currents $I$ approaching
$I_{c}$, where the potential barrier is substantially reduced. In underdamped
junctions where motion of the phase particle does not dissipate much energy,
there is a high likelihood that once the particle thermally activates or
tunnels out of the well, it will continue to time evolve in the so-called
free-running state, leading to a finite voltage across the JJ [14]. At this
point, the junction is said to have entered the quasiparticle regime, where
quasiparticle transport across the junction results in ohmic dissipation.
Consequently, the switching current $I_{sw}$ at which the JJ transitions to
the finite voltage state is typically less than the critical current $I_{c}$
expected from the Ambegaokar-Baratoff relation Eqn. 3. Thus, one expects
small-$C$ junctions to be more susceptible to early switching, such that the
ratio $I_{sw}/I_{c}$ decreases with increasing $E_{C}/E_{J}$, a trend that has
been recently verified experimentally [19].

Relevant for the operation of qubits, at $I\sim 0$ there is a finite tunneling
probability between adjacent minima of the periodic washboard potential. If
there exists an available state in a neighboring well, the particle can
experience a phase jump of $2\pi$ without entering the full running state.
Over time, many successive phase jumpss results in phase diffusion associated
with a finite voltage for currents $I<I_{sw}$ in the nominally superconducting
regime of the JJ. The corresponding “zero-bias” resistance $R_{0}$ has been
observed for many years in small-$C$ JJs, beginning with Ref. [20].
Theoretically, one expects $R_{0}$ to increase super-exponentially with
increasing $E_{C}/E_{J}$, a result that has also been verified experimentally
[21, 22, 23, 16].

With the amount of attention JJs have garnered in the last half-century, one
may assume our understanding of single junction dynamics to be complete.
However, in anticipation of the following discussion, we reiterate two
important pieces of information. First, we point out that phase tunneling is
not exclusive to junctions where $E_{C}\sim E_{J}$. The extent to which these
phase fluctuations negatively influence qubit operation is unknown. Second, we
emphasize the incomplete nature of our understanding of how environmental
impedance influences such JJ phase dynamics, especially in the context of the
high-impedance environment of superconducting qubits. A goal of this paper is
to draw attention to features inherent to single Josephson junctions that
could pose challenges to superconducting qubit performance.

##  II Device fabrication and measurement details

The Al/AlOx/Al junctions intended for superconducting qubits included in this
study were fabricated using electron-beam lithography without a suspended
bridge [24]. Details of the fabrication process can be found in the
Supplementary Materials, but we point out key features relevant to our
discussion here. The areas of the junctions were defined by a fixed width of
240 nm and lengths varying from 260 to 1510 nm, yielding nominal overlap areas
between 0.06 and 0.36 $\mu\text{m}^{2}$. Two different sets of devices with
these junction areas were fabricated and measured: single “standalone”
junctions with Al bonding pads and single junctions that were connected to Nb
shunting capacitors in a transmon configuration. The intrinsic capacitance of
the JJs is estimated to be approximately 4-22 fF, depending on junction size,
while the additional shunt capacitance $C_{S}$ in the transmon configuration
is calculated to be $\sim$100 fF by electrostatic simulation prior to
fabrication.

To address concerns over intrinsic defects associated with this particular
fabrication process, we also fabricated and measured standalone junctions
using the familiar Dolan shadow evaporation technique [25] to compare with the
bridgeless junctions. A double angle of $\pm 15^{\circ}$ in PMMA/MMA bilayer
resist with 100 mTorr of oxygen introduced between Al depositions was used to
create Dolan junctions with similar areas to the bridgeless process.

Transport measurements were performed in an Oxford MX100 dilution refrigerator
at temperatures down to $\sim$20 mK using a four-terminal measurement
configuration. For the dc IV measurements, the dc current was applied using a
custom room temperature current source with a nominal output impedance of
$\sim 10^{13}\;\Omega$, while the voltage across the JJ was measured using a
custom room temperature voltage preamp with an input impedance of
$10^{9}\;\Omega$. For the differential resistance ($dV/dI$) measurements, a
low frequency (10s of Hz) ac current was superposed on the dc current, and the
resulting ac signal was demodulated using a PAR124 lock-in amplifier. In order
to capture the full extent of the features observed at low bias, ac currents
of order 25 pA rms were used. It is important to note that unlike some other
works on small capacitance JJs, no large on-chip or low-temperature resistors
were used to isolate the JJ. This means that no offset voltages due to the
Coulomb blockade of charge were expected in the IV characteristics, a fact
that was confirmed by applying a large magnetic field to suppress
superconductivity. With such a field applied, the IV characteristic appeared
perfectly linear and antisymmetric with respect current drive, with a slight
increase in slope near $I=0$ associated with the dynamical Coulomb blockade
[26] discernible only in the $dV/dI$ traces (see Supplementary Materials).

##  III Experimental results

Figure 1 shows a typical IVC for a standalone JJ in the regime
$E_{C}/E_{J}\sim 1/10$ with area $\sim 0.25\;\mu\text{m}^{2}$. Since several
features in this IVC are critical to our discussion of transmons, we explain
them in detail in the following text.

Figure 1: Simultaneously measured IVC and dV/dI on small standalone junctions
fabricated using a bridgeless evaporation technique (top) and traditional
Dolan shadow evaporation (bottom) at 20 mK. Inset: scanning electron
microscope images of junctions. Scale bar denotes $1\;\mu\text{m}$.

Before cooling superconducting qubits in a dilution refrigerator, $E_{J}$ is
often estimated from the measured junction resistance at room temperature and
the assumed gap $\Delta$ of aluminum. However, one can also directly extract
$E_{J}$ from the IVC at low temperature using Eqn. 3 where
$\Delta(T\rightarrow 0)$ and $R_{n}$ are measured quantities. Consider, for
example, the bridgeless junction IVC in Fig. 1. The dotted lines indicate the
slope at high bias ($R_{n}$) and the onset of the quasiparticle branch
($V=2\Delta/e$) used to calculate $E_{J}$ and $I_{c}$ from Eqn. 3. With
$R_{n}=11.42\;\Omega$ and $\Delta=188\;\mu{eV}$, the Josephson energy $E_{J}$
of the bridgeless junction is found to be $55\;\mu\text{eV}$. Since tunnel
junction resistances increases at low temperature and small differences in
aluminum disorder vary $\Delta$, we suspect that, when possible, $E_{J}$
extracted from the IVC provides a more accurate description for junction
energy compared to room temperature approximations.

While $E_{C}$ cannot be directly determined from the IVC, an estimate of $C$
given the dielectric constant of aluminum oxide and junction size yields an
approximate charging energy $E_{C}=4\;\mu eV$ [27]. For the bridgeless
standalone junction in Fig. 1, then, $E_{C}/E_{J}\sim 1/14$ and
$I_{C}=27\;\text{nA}$. Notably, however, the current at which the junction is
measured to switch occurs at a much smaller value, the average being
$I_{sw}\sim 3\;\text{nA}$ for this particular sample. The switching current is
also reduced in the Dolan junction, though the ratio $I_{c}/I_{sw}$ is larger,
agreeing with trends in the literature for $E_{J}$ and $I_{sw}$ [19].

Figure 2: Transport measurements of standalone JJ (left) and junction shunted
by a 100 fF on-chip capacitor (right). Standalone JJ IVC (a) with optical
image of entire device (top inset) and transmon JJ IVC (b) with schematic
showing circular shunt capacitor and protruding wirebond pads (top inset).
Arrows show direction of bias current sweep magnified near zero bias.
Standalone and transmon JJ $dV/dI$ (c) and (d), respectively, with nonlinear
$R_{0}$ around zero bias (insets) magnified to show experiment (red) and
theory (black).

One may observe switching currents $I_{sw}\ll I_{c}$ for a several reasons,
e.g. inadvertent electron heating due to the biasing circuit. At millikelvin
temperature, however, we expect quantum phase fluctuations to dictate the
switching behavior. Since $I_{sw}$ is the maximum current carried in the phase
diffusion regime of the IVC before spontaneously transitioning into a higher
voltage state, it gives the experimentalist an idea of the degree of
environmental noise coupled to the system. $I_{sw}$ is a stochastic variable
that depends weakly on dc slew rate, or how quickly the potential is tilted; a
slower tilt allows more time for probabilistic particle escape, which causes
early activation into the running state. Varying the slew rate from 1 aA/s to
1 $\mu$A/s extends the average $I_{sw}$ from 2 nA to 5 nA (see Supplementary
Material for further details). For the purposes of this experiment, we take
$I_{sw}$ to be the peak value of the switching distribution.

Though not critical to the aim of this paper, we point out the signature
feature of underdamped junctions is hysteresis in the IVC, as seen in Fig. 1.
The potential tilt required to relocalize an underdamped phase particle is
steeper than the tilt at $I=I_{sw}$. The so-called retrapping current $I_{r}$
is therefore a smaller fraction of $I_{c}$ in underdamped junctions. The
difference between $I_{sw}$ and $I_{r}$ appears as hysteretic behavior in both
directions of current sweep of the IVC.

The finite slope in Fig. 1 where $I<I_{sw}$ depicts phase diffusion of a
continuously tunneling Josephson particle described earlier in the text. Phase
diffusion in small JJs has been experimentally observed for decades [21, 22,
6, 23] with various theoretical explanations [7, 28, 18, 8]. In the present
case, the zero-bias resistances $R_{0}=dV/dI$ as $I\rightarrow 0$ agrees with
past experiments on underdamped JJs of similar sizes [19]. Moreover, we also
observe non-linear dependence of $R_{0}$ with current bias $I$ which, until
very recently, lacked a theoretical model [29] that we briefly summarize here.
When current across the junction is increased, there are two competing
phenomena describing the resistance near zero bias: exponentially increasing
tunneling probability from the diminishing height of the cosine potential and
decreasing ground state energy of the phase particle. These combine to form
the nonlinearity of $R_{0}(I)$ observed in Fig. 1.

It has been suspected that the measurement circuit influences $R_{0}$ much the
same way that $R$ in the RCSJ model depends on environmental impedance.
However, replacing the high-impedance current source with a pseudo-current
source setup utilizing a voltage-biased $100\;\text{k}\Omega$ source resistor
in series with the JJ had no effect on the phase diffusion slope of the IVC.
Additionally, Ref. [29] found that, compared to the ratio $E_{C}/E_{J}$,
environmental admittance $Y(\omega)$ has little effect on $R_{0}$.

As $I\rightarrow 0$, $R_{0}=3.63$ k$\Omega$ in the bridgeless junction and
$R_{0}=600$ k$\Omega$ in the Dolan junction. As previously discussed, zero-
bias resistance $R_{0}$ represents phase fluctuations or the continuous
diffusion of phase through the cosine potential. $R_{0}$ on the order of
$10^{3}\;\Omega$ at milliKelvin temperature is attributed predominately to
phase tunneling since the crossover temperature from tunneling to thermal
activation [30] greatly exceeds the temperatures used in this experiment. We
assume low electron temperature for several reasons. First, milliKelvin
temperature is substantiated in the sharpness of the quasiparticle tunneling
regime of the IVC. Deliberately increasing temperature results in Fermi-
function-like rounding at $V=2\Delta/e$ as $T\rightarrow\;$300 mK. The onset
of the quasiparticle branchs in Fig. 1 are considered substantially sharper.
Second, we note that the region of differential negative resistance observed
in the $dV/dI$ disappears in warm environments. This region is associated with
the time-correlated tunneling of single Cooper pairs, with the back-bending
region in the IVC defined as the so-called “Bloch nose” [31]. The Bloch nose
feature is washed out at higher temperature, similar to switching current.

##  IV Discussion

Transport characterization was performed on over 20 capacitively shunted and
unshunted junctions fabricated using bridgeless evaporation and Dolan double
angle techniques. A main finding of this paper is the surprising similarity
between the transport characteristics of single standalone junctions and
identically fabricated junctions shunted with a large on-chip capacitor. As
previously mentioned, the Langevin model predicts that as the effective device
capacitance is increased with respect to constant $E_{J}$, one expects a
decreased susceptibility to the random force $F(t)$ that results in reduced
switching currents out of a local potential minimum. According to our
measurements, however, there is little difference in $I_{sw}$ for capacitively
shunted and unshunted junctions. Figure 2 shows transport measurements
performed on two junctions with nearly identical $E_{J}$ but different
$E_{C}$, one in a transmon configuration, the other standalone. Despite
$E_{C}/E_{J}$ ranging from $1/10$ to $1/50$, switching currents in both
devices agree very well with previous data for $E_{C}\sim E_{J}>k_{b}T$ [19].
This is an unexpected result. dc transport theory developed for JJs where
$E_{C}/E_{J}\ll 1$ – under which one may traditionally categorize transmons –
does not adequately explain the reduced switching current observed in these
junctions. It is suggested that delta-correlated Johnson noise of the biasing
electronics generates incoherent phase slips which drive the junction into the
voltage state $V=2\Delta/e$ [28]. However, the high-impedance current source
used in this study lie outside this regime and are not considered to influence
switching currents. This heavily implies that the intrinsic capacitance of the
junction, not the shunt, dictates phase dynamics. It appears that while
$E_{C}/E_{J}$ determines qubit resonance and anharmonicity, it is not quite as
straightforward a measure of the degree of phase fluctuations in capacitively
shunted junctions.

Figure 3: Simulated zero-bias resistance $R_{0}$ as a function of
$E_{C}/E_{J}$ using quantum Langevin equation from Ref. [29]. $10^{3}\;\Omega$
resistance corresponds to unity $E_{C}/E_{J}$ regardless of shunt resistance
$R$. The value of $E_{J}$ is fixed and corresponds to $I_{C}=30$ nA.

Additionally, transmon phase fluctuations mirror that of the standalone
junction in the zero-bias resistance $R_{0}$. The differential resistances for
the shunted and unshunted JJs in Fig. 2 closely resemble each other. Extracted
from the $dV/dI$, we find $R_{n}=11.42$ for the standalone junction and
$R_{n}=11.06\;\text{k}\Omega$ for the transmon, implying a similarity in
junction area between the two. As $I\rightarrow 0$, $R_{0}=$ 2.78 k$\Omega$
for the transmon and $R_{0}=$ 3.63 k$\Omega$ for the standalone JJ. These
numbers are typical for junctions of this size measured in this study. The
fact that $R_{0}$ does not vanish in the capacitively shunted junction is
somewhat surprising. According to Ref. [29], $R_{0}$ depends exponentially on
$E_{C}/E_{J}$ for junctions with $E_{J}\approx 50\;\mu\text{eV}$. We would
expect, then, a drastic decrease in $R_{0}$ with decreasing $E_{C}/E_{J}$.
However, equating the measured $R_{0}$ to the current-dependent tunneling rate
suggests the effective ratio $E_{C}/E_{J}$ is very nearly the same between the
shunted and unshunted junction. In fact, according to calculations using the
quantum Langevin equation and dissipative quantum tunneling for underdamped
junctions [18, 29], the magnitude of $R_{0}$ for the standalone and transmon
JJs correspond to $E_{C}/E_{J}\sim 1$. See Fig. 3. The energy ratio
$E_{C}/E_{J}$ calculated from the IVC and junction area is incommensurate with
the measured zero-bias resistance; one expects $R_{0}\ll 1$ k$\Omega$ for
junctions with $E_{C}/E_{J}<1$. In single junctions, $R_{0}$, which describes
quantum tunneling, has been to found to scale inversely to $E_{J}$ [32]. To
our knowledge, however, measurements on small-$C$ single junctions with
$R_{0}\sim 1$ k$\Omega$ and well-defined switching current have never been
reported in the literature. Puzzling still is the robust nature of $R_{0}$,
which varies little between capacitively shunted and unshunted junctions and
implies the external shunt capacitor has little effect on the presence of
phase fluctuations.

The authors note an absence of transport measurements performed on single
Josephson junctions destined for superconducting qubits in the literature.
Recent devices utilizing SQUIDs, JJ arrays, on-chip impedances and cold
resistors tend to focus on ultra-small junctions where $E_{C}/E_{J}\leq 1$.
With the data reported in this paper, we emphasize isolated, inherent
Josephson phase dynamics that lie at the heart of the qubit. We believe dc
transport is a powerful tool in the growing repertoire of qubit
characterization.

In summary, we observe two features characteristic of small-capacitance JJs –
reduced switching current and zero-bias resistance – in the transmon
configuration, suggesting that phase fluctuations, which signify intrinsic
noise, ignore the shunt capacitance and could fundamentally limit the
performance of superconducting qubits. Our findings suggest the intrinsic
capacitance of the junction matters more than previously assumed, and might
motivate work on qubits that incorporate the shunt into the junction itself,
as in the merged-element transmon design [33], which could potentially
mitigate phase fluctuations.

This work is supported by the U.S. Department of Energy, Office of Science,
National Quantum Information Science Research Centers, Superconducting Quantum
Materials and Systems Center (SQMS) under Contract No. DEAC02-07CH11359.

## References

  * Voss and Webb [1981] R. F. Voss and R. A. Webb, Macroscopic Quantum Tunneling in Nb Josephson Junctions, Physical Review Letters 47, 265 (1981).
  * Devoret  _et al._ [1985] M. H. Devoret, J. M. Martinis, and J. Clarke, Measurements of Macroscopic Quantum Tunneling out of the Zero-Voltage State of a Current-Biased Josephson Junction, Physical Review Letters 55, 1908 (1985).
  * Ambegaokar  _et al._ [1982] V. Ambegaokar, U. Eckern, and G. Schön, Quantum Dynamics of Tunneling between Superconductors, Physical Review Letters 48, 1745 (1982).
  * Fisher and Zwerger [1985] M. P. A. Fisher and W. Zwerger, Quantum Brownian motion in a periodic potential, Physical Review B 32, 6190 (1985).
  * Clarke and Schön [1986] J. Clarke and G. Schön, Macroscopic Quantum Phenomena in Josephson Elements, Europhysics News 17, 94 (1986).
  * Kautz and Martinis [1990] R. L. Kautz and J. M. Martinis, Noise-affected I \- V curves in small hysteretic Josephson junctions, Physical Review B 42, 9903 (1990).
  * Biswas and Jha [1970] A. C. Biswas and S. S. Jha, Effect of Thermal Noise on the dc Josephson Effect, Physical Review B 2, 2543 (1970).
  * Caldeira and Leggett [1981] A. O. Caldeira and A. J. Leggett, Influence of Dissipation on Quantum Tunneling in Macroscopic Systems, Physical Review Letters 46, 211 (1981).
  * Nguyen  _et al._ [2019] L. B. Nguyen, Y.-H. Lin, A. Somoroff, R. Mencia, N. Grabon, and V. E. Manucharyan, High-coherence fluxonium qubit, Phys. Rev. X 9, 041041 (2019).
  * Gyenis  _et al._ [2021] A. Gyenis, P. S. Mundada, A. Di Paolo, T. M. Hazard, X. You, D. I. Schuster, J. Koch, A. Blais, and A. A. Houck, Experimental realization of a protected superconducting circuit derived from the $0$–$\pi$ qubit, PRX Quantum 2, 010339 (2021).
  * Somoroff  _et al._ [2023] A. Somoroff, Q. Ficheux, R. A. Mencia, H. Xiong, R. Kuzmin, and V. E. Manucharyan, Millisecond coherence in a superconducting qubit, Phys. Rev. Lett. 130, 267001 (2023).
  * Ford  _et al._ [1988] G. Ford, J. Lewis, and R. O’Connell, Dissipative quantum tunneling: quantum langevin equation approach, Physics Letters A 128, 29 (1988).
  * Balakrishnan [2021] V. Balakrishnan,  _Elements of Nonequilibrium Statistical Mechanics_ (Springer, Cham, 2021).
  * Tinkham [2004] M. Tinkham,  _Introduction to Superconductivity_, 2nd ed. (Dover Publications, 2004).
  * Ambegaokar and Baratoff [1963] V. Ambegaokar and A. Baratoff, Tunneling Between Superconductors, Physical Review Letters 10, 486 (1963).
  * Joyez  _et al._ [1999] P. Joyez, D. Vion, M. Götz, M. H. Devoret, and D. Esteve, The Josephson Effect in Nanoscale Tunnel Junctions, Journal of Superconductivity 12, 757 (1999).
  * Orr  _et al._ [1986] B. G. Orr, J. R. Clem, H. M. Jaeger, and A. M. Goldman, Phase fluctuations in josephson junctions, Phys. Rev. B 34, 3491 (1986).
  * Hu and O’Connell [1992] G. Y. Hu and R. F. O’Connell, Low-voltage resistance in small Josephson junctions, Journal of Physics: Condensed Matter 4, 9635 (1992).
  * Lu  _et al._ [2023] W.-S. Lu, K. Kalashnikov, P. Kamenov, T. J. DiNapoli, and M. E. Gershenson, Phase Diffusion in Low-EJ Josephson Junctions at Milli-Kelvin Temperatures, Electronics 12, 416 (2023).
  * Martinis  _et al._ [1985] J. M. Martinis, M. H. Devoret, and J. Clarke, Energy-level quantization in the zero-voltage state of a current-biased josephson junction, Phys. Rev. Lett. 55, 1543 (1985).
  * Iansiti  _et al._ [1987] M. Iansiti, A. T. Johnson, W. F. Smith, H. Rogalla, C. J. Lobb, and M. Tinkham, Charging energy and phase delocalization in single very small Josephson tunnel junctions, Physical Review Letters 59, 489 (1987).
  * Iansiti  _et al._ [1988] M. Iansiti, A. T. Johnson, C. J. Lobb, and M. Tinkham, Crossover from Josephson Tunneling to the Coulomb Blockade in Small Tunnel Junctions, Physical Review Letters 60, 2414 (1988).
  * Martinis and Kautz [1989] J. M. Martinis and R. L. Kautz, Classical phase diffusion in small hysteretic Josephson junctions, Physical Review Letters 63, 1507 (1989).
  * Nersisyan  _et al._ [2019] A. Nersisyan, S. Poletto, N. Alidoust, R. Manenti, R. Renzas, C.-V. Bui, K. Vu, T. Whyland, Y. Mohan, E. A. Sete, S. Stanwyck, A. Bestwick, and M. Reagor, Manufacturing low dissipation superconducting quantum processors, in  _2019 IEEE International Electron Devices Meeting (IEDM)_ (2019) pp. 31.1.1–31.1.4.
  * Dolan [1977] G. J. Dolan, Offset masks for lift‐off photoprocessing, Applied Physics Letters 31, 337 (1977).
  * Grabert [1991] H. Grabert, Single charge tunneling: A brief introduction, Z. Phys. B 85, 319 (1991).
  * [27] In a parallel plate approximation, $\text{C}=\epsilon_{r}\epsilon_{0}\frac{A}{d}\approx 20\;\text{fF}$ using $d=1$ nm between electrodes.
  * Ivanchenko and Zil’Berman [1969] Y. M. Ivanchenko and A. Zil’Berman, The josephson effect in small tunnel contacts, Sov. Phys. JETP 28, 1272 (1969).
  * Chandrasekhar [2024] V. Chandrasekhar, Current dependence of the low bias resistance of small capacitance josephson junctions, Physics Letters A 507, 129493 (2024).
  * Blackburn  _et al._ [2016] J. A. Blackburn, M. Cirillo, and N. Grønbech-Jensen, A survey of classical and quantum interpretations of experiments on josephson junctions at very low temperatures, Physics Reports 611, 1 (2016).
  * Corlevi [2006] S. Corlevi,  _Quantum effects in nanoscale Josephson junction circuits_ , Ph.D. thesis, KTH Royal Institute of Technology (2006).
  * Watanabe and Haviland [2003] M. Watanabe and D. B. Haviland, Quantum effects in small-capacitance single josephson junctions, Phys. Rev. B 67, 094505 (2003).
  * Mamin  _et al._ [2021] H. Mamin, E. Huang, S. Carnevale, C. Rettner, N. Arellano, M. Sherwood, C. Kurter, B. Trimm, M. Sandberg, R. Shelby, M. Mueed, B. Madon, A. Pushp, M. Steffen, and D. Rugar, Merged-element transmons: Design and qubit performance, Phys. Rev. Appl. 16, 024023 (2021).