# Stochastic Pre-Event Preparation for Enhancing Resilience of Distribution
Systems with High DER Penetration

Qianzhi Zhang<EMAIL_ADDRESS>Zhaoyu Wang  Shanshan Ma  Anmar Arif
Department of Electrical and Computer Engineering, Iowa State University,
Ames, IA 50011, USA  School of Electrical, Computer and Energy Engineering,
Arizona State University, Tempe, AZ 85287, USA  Department of Electrical
Engineering, King Saud University, Riyadh 11451, Saudi Arabia

###### Abstract

This paper proposes a stochastic optimal preparation and resource allocation
method for upcoming extreme weather events in distribution systems, which can
assist utilities to achieve faster and more efficient post-event restoration.
With the objective of maximizing served load and minimizing operation cost,
this paper develops a two-stage stochastic mixed-integer linear programming
(SMILP) model. The first-stage determines the optimal positions and numbers of
mobile resources, fuel resources, and labor resources. The second-stage
considers network operational constraints and repair crew scheduling
constraints. The proposed stochastic pre-event preparation model is solved by
a scenario decomposition method, Progressive Hedging (PH), to ease the
computational complexity introduced by a large number of scenarios.
Furthermore, to show the impact of solar photovoltaic (PV) generation on
system resilience, we consider three types of PV systems during power outage
and compare the resilience improvements with different PV penetration levels.
Numerical results from simulations on a large-scale (more than 10,000 nodes)
distribution feeder have been used to validate the scalability and
effectiveness of the proposed method.

###### keywords:

Pre-event preparation , progressive hedging , PV systems , resource allocation
, two-stage stochastic model

††journal: Renewable and Sustainable Energy Reviews

##  1 Introduction

Extreme weather events have brought significant damage to power grid
infrastructure and caused 50%-60% of power outages in the U.S. [1]. Among
those outages, around 90% of them were due to failures in distribution systems
[2]. After severe weather events, the major challenge for utilities is the
shortage of various resources to repair damage and restore power supply. Pre-
event resource allocation is one of the most effective ways to mitigate
extreme events’ impacts on power distribution system. It can allocate
appropriate amounts of flexible resources to optimal positions before the
extreme events. These flexible resources include emergency power supply
resources, equipment resources and labor resources. Therefore, pre-event
preparation enables faster and more efficient post-event restoration of the
power gird.

There are exist studies that have investigated resource allocation problems
for the resilience enhancement of electric distribution systems. In [3, 4, 5],
proactive resource management in microgrids and proactive operation strategies
in distribution systems are considered to enhance system resilience during
extreme events. In [6], the number and location of depots are determined at
the pre-disturbance stage to manage the available resources. In [7], repair
crews are pre-allocated to depots and integrated with restoration process to
enhance the resilience of electric distribution systems. In [8], a two-stage
stochastic model is developed to select staging locations and allocate repair
crews for disaster preparation, while considering distribution system
operation and crew routing constraints. In [9], the authors developed a
stochastic model for optimizing proactive operation actions. The study
optimized the topology of the network and position of crews for upcoming
disturbances. In [10] and [11], a two-stage framework is developed to position
mobile emergency generators (MEGs) for pre- and post-disasters. Mobile energy
storage devices (MESs) are investigated in [12] and [13] for resilience
enhancement of power distribution systems. However, there remain limitations
in the above studies on pre-event preparation and resource allocation. These
limitations are described in the following:

(1) Pre-allocation of various flexible resources: In practice, pre-event
preparation includes allocation of various flexible resources, such as MEGs,
MESs, fuel resources for diesel generators, and repair crews. The optimal
allocation of those flexible resources can help utilities to achieve faster
and more efficient post-event power restoration. However, previous studies
mainly focused on allocating specific flexible resources, rather than
formulating a complete optimization problem to pre-allocate various flexible
resources together.

(2) Impacts of solar PV power on system resilience: Due to intermittent
characteristic of traditional distributed energy resources (DERs), such as
solar power, PV systems are not considered as a reliable resilient solution.
However, the distributed nature of PV power can contribute to a more resilient
power system. In practice, PV systems can be coupled with energy storage
technology, to enable continues operation during outages [14]. However,
different types of PV systems are ignored in most existing research.

(3) Scalability of the solution algorithm: On one side, the stochastic pre-
event preparation model may suffer from computational inefficiency due to a
large number of scenarios; on the other side, a limited number of scenarios
may influence the stability and quality of the solutions. Therefore, the
trade-off between computation time and solution accuracy needs to be studied
for stochastic pre-event preparation methods. In addition, a large-scale
system is needed to verify the scalability of solution algorithms.

To address these challenges, we propose a two-stage stochastic mixed integer
linear program (SMILP) for pre-event preparation with pre-allocation of mobile
resources, fuel resources and labour resources. Furthermore, the proposed pre-
event preparation model considers different types of PV systems and
facilitates the benefits of leveraging high PV penetration for improving the
resilience of distribution grids. In this paper, resilience improvement is
quantified by the increased served load and reduced outage duration. To deal
with the massive computation burden, the proposed two-stage stochastic pre-
event preparation problem is solved by a scenario decomposition method,
Progressive Hedging (PH) [15], while maintaining the accuracy and stability of
the solution [16]. Also, the quality of the solution is validated by a
multiple replication procedure (MRP). The main contribution of this paper is
three-folded:

  * 1.

We propose a two-stage SMILP model for pre-event preparation, where the first-
stage allocates MEGs, MESs, fuel, and repair crews, while the second-stage
considers distribution system operation and repair crew scheduling
constraints.

  * 2.

The proposed model considers three types of PV systems. We also demonstrate
the improvement of resilience and the reduction of outage duration with
different PV penetration levels.

  * 3.

The proposed solution algorithm is tested through a solution validation method
to show its quality. In addition, a large-scale system, consisting of more
than 10,000 nodes, is used to verify the scalability of the proposed pre-event
preparation model.

The remainder of the paper is organized as follows: Section 2 describes the
proposed two-stage SMILP for pre-event preparation and resource allocation.
Section 3 presents the PH solution algorithm, convergence analysis and
solution validation. Simulation results and conclusions are given in Section 4
and Section 5, respectively.

##  2 Two-stage Stochastic Pre-event Preparation Model

The general framework of the proposed two-stage stochastic pre-event
preparation model is shown in Figure 1.

Figure 1: The proposed two-stage stochastic pre-event model.

Damage scenarios for extreme weather events are generated based on: (1)
identification of extreme weather events, such as flood, hurricane and winter
storm; (2) extreme weather event data and metric; (3) fragility model of test
systems, which describes the behavior of components under extreme weather
events; (4) damage status of components in test system subject to specific
extreme weather events. To approximate the impact of extreme weather events to
grid infrastructures, damage scenarios can be generated by mapping the weather
data set to failure probability of grid infrastructures. Adopted from [17],
the failure probability of an overhead line being damaged by hurricane can be
expressed as follows:

$p_{l,ij}(w(t))=1-\prod_{k=1}^{m}\Big{(}1-p_{l_{k}}(w(t))\Big{)}\prod_{k=1}^{n}\Big{(}1-p_{fc,k}(w(t))\Big{)}$
(1)

where $p_{l,ij}(w(t))$ is the failure probability of the overhead line $ij$
with wind speed $w(t)$. $p_{l_{k}}(w(t))$ is defined as the conditional
failure probability of pole $k$ at line $ij$ as a log-normal cumulative
distribution function (CDF) of the wind speed $w(t)$, which is expressed in
equation (2). $m$ and $n$ are the number of distribution poles supporting line
$ij$ and the number of conductor wires between two adjacent poles at line
$ij$, respectively. In equation (3), $p_{fc,k}(w(t))$ represents the failure
probability of conductor $k$ between two poles.

$p_{l_{k}}(w(t))=\Phi\Big{[}\ln(\frac{w(t)/m_{R}}{\xi_{R}})\Big{]}$ (2)
$p_{fc,k}(w(t))=(1-p_{u})\max\Big{(}p_{fw,k}(w(t)),\alpha
p_{ftr,k}(w(t))\Big{)}$ (3)

where $m_{R}$ and $\xi_{R}$ are the median capacity and the logarithmic
standard deviation of intensity measurement, respectively. $p_{fw,k}(w(t))$
represents the direct wind-induced failure probability of conductor $k$ and
$p_{ftr,k}(w(t))$ represents the fallen tree-induced failure probability of
conductor $k$. $p_{u}$ is the probability that conductor $k$ is underground,
which is more invulnerable to extreme weather events. $\alpha$ represents the
average tree-induced damage probability of overhead conductors. More details
of weather forecasting methodologies, line fragility models and scenario
generation can be found in [18].

As shown in Figure 1, the proposed SMILP pre-event preparation model has two
stages: (i) Flexible resources are allocated in the first-stage, including the
optimal number and position of MEGs and MESs, allocation of available fuel to
generators, and pre-position of repair crews to depots. (ii) The second-stage
optimizes the operation of the distribution system and assign crews to the
damaged components. Constraints in the second-stage includes unbalanced
optimal power flow constraints, network reconfiguration and isolation
constraints, and repair crew scheduling constraints.

###  2.1 Model Objective Function

The objective function (2.1) is set to minimize operation costs and maximize
load served. There are three cost related terms in the objective, cost of fuel
$C^{\rm F}$, cost of switching operation $C^{\rm SW}$, and cost of load
shedding $C_{i}^{\rm D}$. The objective is formulated as follows:

$\displaystyle\min\sum_{\forall s}Pr(s)$ $\displaystyle\Big{(}C^{\rm F}r^{\rm
F}\sum_{\forall t}\sum_{\forall\phi}\sum_{\forall i}P^{\rm
G}_{i,\phi,t,s}+C^{\rm SW}\sum_{\forall t}\sum_{\forall k\in\Omega_{\rm
SW}}\gamma_{k,t,s}$ $\displaystyle+\sum_{\forall
t}\sum_{\forall\phi}\sum_{\forall i}C_{i}^{\rm
D}(1-y_{i,t,s})d_{i,\phi,t}^{\rm p}\Big{)}$ (4)

where $Pr(s)$ is the probability of occurrence for scenario $s$, $r^{\rm F}$
is the rate of fuel consumption of a generator, and $P_{i,\phi,t,s}^{\rm G}$
is the active power output for fuel-based generator at bus $i$, phase $\phi$,
time $t$, and scenario $s$. Binary variable $\gamma_{k,t,s}$ represents the
status of each switch, if switch $k$ is operated at time $t$ on scenario $s$,
then $\gamma_{k,t,s}=1$. The binary variable $y_{i,t,s}$ represents the status
of load at bus $i$, time $t$, and scenario $s$. If the demand
($d_{i,\phi,t}^{\rm p}$) is served, then $y_{i,t,s}=1$.

###  2.2 First-Stage Constraints

The first-stage constraints revolve around pre-allocating four critical
resources that will be utilized after an extreme event: (i) MEGs, (ii) MESs,
(iii) fuel and (iv) repair crews.

####  2.2.1 Mobile Resources Allocation Constraints

Mobile resources can be used to restore energy for isolated areas that are not
damaged, and to restore critical customers. In addition, fuel management is
critical after an extreme event to operate emergency generators. Distributing
fuel after an extreme event maybe difficult due to road conditions. As for
repair crews, pre-assigning them to different locations provides a faster and
more organized response. The constraints for allocating the mobile resources
are modeled as follows:

$\sum_{\forall i\in\Omega_{\rm CN}}n_{i}^{\rm MEG}=N^{\rm MEG}$ (5)
$\sum_{\forall i\in\Omega_{\rm CN}}n_{i}^{\rm MES}=N^{\rm MES}$ (6)
$n_{i}^{\rm MEG}+n_{i}^{\rm MES}\leq N_{i}^{\rm MU},\forall i\in\Omega_{\rm
CN}$ (7)

where binary variables $n_{i}^{\rm MEG}$ and $n_{i}^{\rm MES}$ equal 1 if an
MEG and MES are allocated to bus $i$, respectively. The set $\Omega_{\rm CN}$
represents the set of candidate buses for MEGs and MESs. Constraints (5) and
(6) indicates that the number of installed MEGs and MESs are equal to the
number of available devices ($N^{\rm MEG}$ and $N^{\rm MES}$). We assume that
each bus can only have a limited number of mobile units $N_{i}^{\rm MU}$,
which is enforced by (7).

####  2.2.2 Fuel Resources Allocation Constraints

Define the set $\Omega_{\rm G}=\Omega_{\rm EG}\cup\Omega_{\rm CN}$, where
$\Omega_{\rm EG}$ is the set of buses that have fuel-based emergency
generators. The fuel allocated to $\Omega_{\rm G}$ must be limited to the
available amount of fuel. We model the fuel allocation constraints as follows:

$\sum_{\forall i\in\Omega_{G}}n^{\rm Fuel}_{i}\leq N^{\rm Fuel}$ (8)
$F_{i}^{\rm G}\leq n^{\rm Fuel}_{i}\leq F_{i}^{\rm max},\forall
i\in\Omega_{\rm G}$ (9)

Constraint (8) limits the total amount of allocated fuel to the amount of fuel
available ($N^{\rm{Fuel}}$), where $n^{\rm Fuel}_{i}$ is the amount of fuel
allocated to the generator at bus $i$. Constraint (9) limits the amount of
fuel on each site, where $F_{i}^{\rm G}$ is the amount of fuel already present
for the generator at bus $i$, and $F_{i}^{\rm max}$ represents the maximum
capacity of fuel at bus $i$.

####  2.2.3 Repair Crew Allocation Constraints

In order to allocate the repair crews, we divide the network into different
regions $\Omega_{\rm R}$. Each region will be assigned with different crews,
who will conduct the repairs in that region. The repair crews are pre-
positioned to the regions using constraints (10) and (11), as follows:

$\sum_{\forall r\in\Omega_{\rm R}}n^{\rm Crew}_{r}=N^{\rm Crew}$ (10)
$N_{r}^{\rm Crew,min}\leq n^{\rm Crew}_{r}\leq N_{r}^{{\rm Crew,max}},\forall
r\in\Omega_{\rm R}$ (11)

where $n^{\rm Crew}_{r}$ is the number of repair crews in region $r$ and
$N^{\rm Crew}$ is the total number of crews. The number of repair crews is
limited in each region, using $N_{r}^{\rm Crew,min}$ and $N_{r}^{{\rm
Crew,max}}$, depending on the size and capacity of the staging locations.

###  2.3 Second-Stage Constraints

In the second-stage of the proposed pre-event preparation model, the
constraints of PV systems and repair crew dispatch are mainly discussed. The
model also considers unbalanced power flow constraints, voltage constraints
and reconfiguration constraints [19, 20].

####  2.3.1 PV System Constraints

To fully investigate the impact of PV systems on system resilience, three
types of PV systems [20] are considered in the second-stage, $\Omega_{\rm
PV}=\Omega^{\rm G}_{\rm PV}\cup\Omega^{\rm H}_{\rm PV}\cup\Omega^{\rm C}_{\rm
PV}$: (i) Type 1: on-grid (grid-following) PV ($\Omega^{\rm G}_{\rm PV}$),
where during an outage, the PV is switched off. (ii) Type 2: hybrid on-
grid/off-grid PV + energy storage system (ESS) ($\Omega^{\rm H}_{\rm PV}$),
where the PV system operates on-grid in normal conditions, and off-grid during
an outage. (iii) Type 3: grid-forming PV + ESS with grid-forming capability
($\Omega^{\rm C}_{\rm PV}$), this system can restore part of the network that
is not damaged if the fault is isolated. The output power of the PV systems is
determined using the following equations:

$0\leq P^{\rm PV}_{i,\phi,t,s}\leq\frac{I_{r_{i,t,s}}}{1000W/m^{2}}P_{i}^{\rm
rate},\forall i\in\Omega_{\rm PV}/\Omega^{\rm G}_{\rm PV},\phi,t,s$ (12)
$0\leq P^{\rm
PV}_{i,\phi,t,s}\leq\chi_{i,t,s}\frac{I_{r_{i,t,s}}}{1000W/m^{2}}P_{i}^{\rm
rate},\forall i\in\Omega^{\rm G}_{\rm PV},\phi,t,s$ (13) $(P^{\rm
PV}_{i,\phi,t,s})^{2}+(Q^{\rm PV}_{i,\phi,t,s})^{2}\leq(S^{\rm
PV}_{i})^{2},\forall i\in\Omega_{\rm PV}/\Omega^{\rm G}_{\rm PV},\phi,t,s$
(14) $(P^{\rm PV}_{i,\phi,t,s})^{2}+(Q^{\rm
PV}_{i,\phi,t,s})^{2}\leq\chi_{i,t,s}(S^{\rm PV}_{i})^{2},\forall
i\in\Omega^{\rm G}_{\rm PV},\phi,t,s$ (15)

The active power output $P^{\rm PV}_{i,\phi,t,s}$ of a PV depends on the
rating of the solar cell $P_{i}^{\rm rate}$ and the solar irradiance
$I_{r_{i,t,s}}$ [21]. The generated output power from the PV can be determined
in constraints (12) and (13), respectively. The binary variable $\chi_{i,t,s}$
equals 1 if bus $i$ is energized at time $t$ and scenario $s$. Using advanced
PV smart inverters [22], the PVs can provide reactive power support $Q^{\rm
PV}_{i,\phi,t,s}$, which is constrained by the capacity $S^{\rm PV}_{i}$ in
(14) and (15). During an outage, on-grid PVs are disconnected and the on-site
load is not served by the PVs, therefore, constraints (13) and (15) are
multiplied by $\chi_{i,t,s}$. PV systems of types $\Omega^{\rm C}_{\rm PV}$
and $\Omega^{\rm H}_{\rm PV}$ can disconnect from the grid and serve the on-
site load.

An example network with damaged line is given in Figure 2, where the network
is divided into three islands due to the damaged line. In this work, we assume
that the network can be restored using the grid-forming sources in
$\Omega^{\rm C}_{\rm PV}\cup\Omega_{\rm G}$. While PV system in types
$\Omega^{\rm G}_{\rm PV}$ or $\Omega^{\rm H}_{\rm PV}$ can connect to the grid
only after the PV bus is energized. Island A has a grid-forming generator,
therefore, a microgrid is created and the PV system can participate. Island B
must be isolated because of the damaged line. Island C does not have any grid-
forming generators; hence, it will not be active and the grid-tied PV will be
disconnected.

Figure 2: A single line diagram of an example network with one damaged line.

To determine the connection status of the PV systems, we design a virtual
network in parallel to the distribution network. The example network shown in
Figure 2 is transformed to a virtual network shown in Figure 3.

Figure 3: A virtual network created for the example network in Figure 2.

A virtual network with virtual sources, loads, and flow is built to identify
if an island can be energized by grid-forming generators. Each grid-forming
generator is replaced by a virtual source with infinite capacity. Other power
sources without grid-forming capability (e.g., grid-tied PVs) are removed. The
actual loads are replaced by virtual loads with magnitude of 1. The virtual
network scheme is modeled using constraints (16)-(20).

$\sum_{\forall j\in\Omega^{\rm C}_{\rm PV}\cup\Omega_{\rm G}}v_{j,t,s}^{\rm
S}+\sum_{\forall k\in\Omega_{\rm K}(.,i)}v_{k,t,s}^{\rm
f}=\chi_{i,t,s}+\sum_{\forall k\in\Omega_{\rm K}(i,.)}v_{k,t,s}^{\rm
f},\forall i,t,s$ (16) $-(u_{k,t,s})M\leq v_{k,t,s}^{\rm
f}\leq(u_{k,t,s})M,\forall k\in\Omega_{\rm K},t,s$ (17) $0\leq v_{k,t,s}^{\rm
S}\leq(n_{i}^{\rm MEG}+n_{i}^{\rm MES})M,\forall i\in\Omega_{\rm CN},t,s$ (18)
$\chi_{i,t,s}\geq y_{i,t,s},\forall i\in\Omega_{\rm N}/\\{\Omega^{\rm C}_{\rm
PV}\cup\Omega_{\rm PV}^{\rm H}\cup\Omega_{\rm G}\\},t,s$ (19)
$\chi_{i,t,s}+n_{i}^{\rm MEG}+n_{i}^{\rm MES}\geq y_{i,t,s},\forall
i\in\Omega_{\rm CN},t,s$ (20)

A power-balance equation is added for each virtual bus, which means that if
the virtual load at a bus is served, then that bus is energized. Therefore,
for islands without grid-forming generators, all buses will be de-energized as
the virtual loads in the island cannot be served. Constraint (16) is the node
balance constraint for the virtual network. Virtual source $v^{\rm S}$ is
connected to buses with power sources that have the capability to restore the
system. The variable $v_{k}^{\rm f}$ represents the virtual flow on line $k$
and each bus is given a load of 1 that is multiplied by $\chi_{i}$. Therefore,
$\chi_{i}=1$ (bus $i$ is energized) if the virtual load can be served by a
virtual source and 0 (bus $i$ is de-energized) otherwise. The virtual flow is
limited by (17). The limits are multiplied by the status of the line
($u_{k,t,s}$) so that the virtual flow is 0 if a line is disconnected. The
virtual source can be used only if a generator is installed, as enforced by
(18). Define $\Omega_{\rm N}$ as the set of all buses. If bus $i$ is de-
energized, then the load must be shed (19), unless bus $i$ has a local power
source with disconnect switch. Constraint (20) is similar to (19) but with the
presence of mobile sources.

####  2.3.2 Repair Crews Constraints

In the second-stage, repair crews are assigned to damaged components that are
in the area at which they are positioned. Note that the travel time is
neglected in this study, as the travel distances between components in the
same area is assumed to be small. An example for crew assignment is given in
Figure 4, where two working areas are assigned for the crews. In this example,
four damaged lines in Area 1 will be repaired by crews 1-3, while crews 4 and
5 are responsible for the two damaged lines in Area 2. The repair crews
constraints are formulated as follows:

$\sum_{\forall k\in\Omega_{\rm DL(s)}}z_{k,t,s}\leq n^{\rm Crew}_{r},\forall
r,t,s$ (21) $\sum_{\forall t}z_{k,t,s}\leq T^{r}_{k,s},\forall k\in\Omega_{\rm
DL(s)},s$ (22)
$\frac{1}{T^{r}_{k,s}}\sum_{\tau=1}^{t-1}z_{k,\tau,s}-1+\epsilon\leq
u_{k,t,s}\leq\frac{1}{T^{r}_{k,s}}\sum_{\tau=1}^{t-1}z_{k,\tau,s},\forall
k\in\Omega_{\rm DL(s)},t,s$ (23)

Figure 4: A crew assignment example with 2 depots and 5 crews.

Define $z_{k,t,s}$ as a binary variable that equals 1 if line $k$ is being
repaired at time $t$ on scenario $s$, and $\Omega_{\rm DL(s)}$ as the set of
damaged lines on scenario $s$. Constraint (21) limits the number of repairs
being conducted in each area according to the number of crews $n^{\rm
Crew}_{r}$ available. Constraint (22) defines the repair time for each damaged
line. The line status $u_{k,t,s}$ equals 0 until the repair process is
conducted for $T_{k,s}^{r}$ time periods. Based on constraint (23), let
$T_{k,s}^{r}=3$, $z_{k,t,s}=\\{0,0,1,1,1,0,0\\}$, then
$u_{k,t,s}=\\{0,0,0,0,0,1,1\\}$. For example, when $t=6$ and $\epsilon=0.001$,
then constraint (23) becomes $0.668\leq u_{k,6,s}\leq 1$, therefore,
$u_{k,6,s}=1$.

####  2.3.3 Network Operational Constraints

The next set of constraints are related to the operation of distribution
systems. We consider unbalanced power flow equations, radiality constraints,
fuel consumption, and energy storage constraints. The unbalanced distribution
system constraints are given below:

$\displaystyle\sum_{b\in\Omega_{\rm K}(i,.)}P_{b,\phi,t,s}^{\rm
K}-\sum_{k\in\Omega_{\rm K}(.,i)}$ $\displaystyle P_{k,\phi,t,s}^{\rm
K}=P^{\rm G}_{i,\phi,t,s}+P^{\rm PV}_{i,\phi,t,s}$ $\displaystyle+(P^{\rm
Ch}_{i,\phi,t,s}-P^{\rm Dis}_{i,\phi,t,s})-y_{i,t,s}d^{P}_{i,\phi,t},\forall
i,\phi,t,s$ (24) $\displaystyle\sum_{b\in\Omega_{\rm
K}(i,.)}Q_{b,\phi,t,s}^{\rm K}-\sum_{k\in\Omega_{\rm
K}(.,i)}Q_{k,\phi,t,s}^{\rm K}$ $\displaystyle=Q^{\rm G}_{i,\phi,t,s}+Q^{\rm
PV}_{i,\phi,t,s}$ $\displaystyle+Q^{\rm
ESS}_{i,\phi,t,s}-y_{i,t,s}d^{Q}_{i,\phi,t},\forall i,\phi,t,s$ (25)
$-u_{k,t,s}P_{k}^{\rm K,max}\leq P_{k,\phi,t,s}^{\rm K}\leq
u_{k,t,s}P_{k}^{\rm K,max},\forall k\in\Omega_{\rm K},\phi,t,s$ (26)
$-u_{k,t,s}Q_{k}^{\rm K,max}\leq Q_{k,\phi,t,s}^{\rm K}\leq
u_{k,t,s}Q_{k}^{\rm K,max},\forall k\in\Omega_{\rm K},\phi,t,s$ (27) $0\leq
P_{i,\phi,t,s}^{\rm G}\leq P_{i}^{\rm G,max},\forall i\in\Omega_{\rm
EG},\phi,t,s$ (28) $0\leq Q_{i,\phi,t,s}^{\rm G}\leq Q_{i}^{\rm G,max},\forall
i\in\Omega_{\rm EG},\phi,t,s$ (29) $0\leq P_{i,\phi,t,s}^{\rm G}\leq
n_{i}^{\rm{MEG}}P_{i}^{\rm G,max},\forall i\in\Omega_{\rm CN},\phi,t,s$ (30)
$0\leq Q_{i,\phi,t,s}^{\rm G}\leq n_{i}^{\rm{MEG}}Q_{i}^{\rm G,max},\forall
i\in\Omega_{\rm CN},\phi,t,s$ (31)
$\begin{split}U_{i,\phi,t,s}-U_{j,\phi,t,s}&\geq
2(\hat{R}_{ij}P_{ij,\phi,t,s}^{\rm K}+\hat{X}_{ij}Q_{ij,\phi,t,s}^{\rm K})\\\
&+(u_{k,t,s}+p_{ij,\phi}-2)M,\forall k,ij\in\Omega_{\rm
K},\phi,t,s\end{split}$ (32) $\begin{split}U_{i,\phi,t,s}-U_{j,\phi,t,s}&\leq
2(\hat{R}_{ij}P_{ij,\phi,t,s}^{\rm K}+\hat{X}_{ij}Q_{ij,\phi,t,s}^{\rm K})\\\
&+(2-u_{k,t,s}-p_{ij,\phi})M,\forall k,ij\in\Omega_{\rm
K},\phi,t,s\end{split}$ (33) $\chi_{i,t,s}U^{\rm min}_{i}\leq
U_{i,\phi,t,s}\leq\chi_{i,t,s}U^{\rm max}_{i},\forall i,\phi,t,s$ (34)
$\sum_{k\in\in\Omega_{\rm B(l)}}u_{k,t,s}\leq|\Omega_{\rm B(l)}|-1,\forall
l\in\Omega_{\rm loop},t,s$ (35)

Constraints (2.3.3) and (2.3.3) are the active and reactive nodal power
balance constraints, where $P_{ij,\phi,t,s}^{\rm K}$ and $Q_{ij,\phi,t,s}^{\rm
K}$ are the active and reactive line flows, and $P^{\rm G}_{i,\phi,t,s}$ and
$Q^{\rm G}_{i,\phi,t,s}$ are the power outputs of the generators. The active
charging/discharging and reactive power outputs of energy storage systems are
denoted by $P^{\rm Ch}_{i,\phi,t,s}$, $P^{\rm Dis}_{i,\phi,t,s}$ and $Q^{\rm
ESS}_{i,\phi,t,s}$. Constraints (26)-(27) represent the active and reactive
power limits of the lines, where the limits ($P_{k}^{\rm K,max}$ and
$Q_{k}^{\rm K,max}$) are multiplied by the line status binary variable
$u_{k,t,s}$. Therefore, if a line is disconnected or damaged, power cannot
flow through it. Constraints (28)-(29) limit the output of the generators to
$P_{i}^{\rm G,max}$ and $Q_{i}^{\rm G,max}$. Similarly, we limit the output of
the MEGs in (30)-(31) if an MEG is installed ($n_{i}^{\rm MEG}=1$).

Constraints (32) and (33) calculate the voltage difference along line $k$
between bus $i$ and bus $j$, where $U_{i,\phi,t,s}$ is the square of voltage
magnitude of bus $i$. We use the big-M method to relax constraints (32) and
(33), if lines are damaged or disconnected. $\hat{R}_{ij}$ and $\hat{X}_{ij}$
are the unbalanced three-phase resistance matrix and reactance matrix of line
$ij$, which can be referred to [22]. The vector $p_{ij,\phi}$ represents the
phases of line $ij$. Constraint (34) guarantees that the voltage is limited
within a specified region ($U^{\rm min}_{i}$ and $U^{\rm max}_{i}$), and is
set to 0 if the bus is in an outage area. Constraint (35) can guarantee the
radiality network during the network reconfiguration. In this paper, we assume
that all the possible loops can be identified by depth-first search method.
The set of loops are given by $\Omega_{\rm loop}$, and the set of switches in
loop $l$ is given by $\Omega_{\rm B(l)}$. For each fuel-based generator, the
total fuel consumption $F_{i,s}$ is limited by the available fuel resources
$n^{\rm Fuel}_{i}$ in constraint (36), as follows:

$F_{i,s}=r^{\rm f}\sum_{\forall t}\sum_{\forall\phi}P^{\rm G}_{i,\phi,t,s}\leq
n^{\rm Fuel}_{i},\forall i\in\Omega_{\rm G},\phi,t,s$ (36)

Next, we model the operation constraints for ESSs and MESs. The constraints
include the change in state of charge (SOC), charging and discharging limits,
and reactive power limits. Let $\Omega_{\rm ES}$ be the set of buses with
ESSs, and $\Omega_{\rm ESC}=\Omega_{\rm ES}\cup\Omega_{\rm CN}$. We can then
define the energy storage constraints as follows:

$\begin{split}E^{\rm SOC}_{i,t,s}=&E^{\rm SOC}_{i,t-1,s}+\\\ &\Delta
t\frac{(\sum_{\forall\phi}P^{\rm Ch}_{i,\phi,t,s}\eta_{\rm
Ch}-\sum_{\forall\phi}P^{\rm Dis}_{i,\phi,t,s}/\eta_{\rm Dis})}{E^{\rm
Cap}_{i}},\forall i\in\Omega_{\rm ESC},\phi,t,s\end{split}$ (37) $E^{\rm
SOC,min}_{i}\leq E^{\rm SOC}_{i,t,s}\leq E^{\rm SOC,max}_{i},\forall
i\in\Omega_{\rm ESC},t,s$ (38) $0\leq P^{\rm Ch}_{i,\phi,t,s}\leq
h_{i,t,s}P^{\rm Ch,max}_{i},\forall i\in\Omega_{\rm ESC},\phi,t,s$ (39) $0\leq
P^{\rm Dis}_{i,\phi,t,s}\leq(1-h_{i,t,s})P^{\rm Dis,max}_{i},\forall
i\in\Omega_{\rm ESC},\phi,t,s$ (40) $-Q_{i}^{\rm ESS,max}\leq
Q_{i,\phi,t,s}^{\rm ESS}\leq Q_{i}^{\rm ESS,max},\forall i\in\Omega_{\rm
ES},\phi,t,s$ (41) $0\leq P^{\rm Ch}_{i,\phi,t,s}\leq n_{i}^{\rm MES}P^{\rm
Ch,max}_{i},\forall i\in\Omega_{\rm CN},\phi,t,s$ (42) $0\leq P^{\rm
Dis}_{i,\phi,t,s}\leq n_{i}^{\rm MES}P^{\rm Dis,max}_{i},\forall
i\in\Omega_{\rm CN},\phi,t,s$ (43) $-n_{i}^{\rm MES}Q_{i}^{\rm ESS,max}\leq
Q_{i,\phi,t,s}^{\rm ESS}\leq n_{i}^{\rm MES}Q_{i}^{\rm ESS,max},\forall
i\in\Omega_{\rm CN},\phi,t,s$ (44)

Constraint (37) determines the state of charge of ESSs ($E^{\rm
SOC}_{i,t,s}$). $E^{\rm Cap}_{i}$ denotes the maximum capacity of the storage
system. To ensure safe ESS operation, the SOC and charging ($P^{\rm
Ch}_{i,\phi,t,s}$) and discharging ($P^{\rm Dis}_{i,\phi,t,s}$) power of ESSs
are constrained as shown in (38)-(40). Here, $E^{\rm SOC,min}_{i}$, $E^{\rm
SOC,max}_{i}$, $P^{\rm Ch,max}_{i}$ and $P^{\rm Dis,max}_{i}$ define the
permissible range of SOC, and maximum charging and discharging power,
respectively. In constraints (39)-(40), the binary variable $h_{i,t,s}$
indicates that ESSs cannot charge and discharge at the same time instant. The
ESS charging/discharging efficiency are represented by $\eta_{\rm
Ch}$/$\eta_{\rm Dis}$. The reactive power of ESS, $Q_{i,\phi,t,s}^{\rm ESS}$,
is kept within maximum limit, $Q_{i}^{\rm ESS,max}$, through constraint (41).
For MES units, we add constraint (42)-(43) so that if $n_{i}^{\rm MES}=0$, the
output power is 0 at bus $i$. The same method is applied for the reactive
power in (44).

##  3 Solution Algorithm

When the number of scenarios is finite, a two-stage stochastic problem can be
modeled as a single-stage large linear programming model, where each
constraint in the problem is duplicated for each realization of the random
data. For problems where the number of realization is too large or infinite,
the Monte Carlo sampling technique can be used to generate a manageable number
of scenarios. In this work, we use the scenario decomposing method PH to solve
the proposed two-stage stochastic pre-event preparation problem.

###  3.1 Two-stage Progressive Hedging Algorithm

The proposed two-stage stochastic pre-event preparation problem (2.1)-(44) can
be compactly reformulated with an extensive form (EF) as follows:

$\xi=\min_{x,y_{s}}a^{T}x+\sum_{\forall s}Pr(s)b^{T}_{s}y_{s}$ (45)
$\text{s.t.}\hskip 5.69054pt(x,y_{s})\in Q_{s},\forall s$ (46)

where $a$ and $b_{s}$ are vectors containing the coefficients associated with
the compact first-stage variable $x$ and compact second-stage variable $y_{s}$
in the objective (45), respectively. The constraint (46) represents the
subproblem constraints that ensure a feasible solution. The PH algorithm
decomposes the extensive form into scenario-based subproblems, by relaxing the
non-anticipativity of the first-stage variables. Hence, with the total number
$S$ of scenarios, the proposed stochastic pre-event preparation problem is
decomposed into $S$ subproblems. The proposed two-stage PH algorithm is
presented in Algorithm 1. Define $\tau$ as iteration number, $\rho$ as a
penalty factor and $\epsilon$ as a termination threshold. The PH algorithm
starts by solving the subproblems with individual scenarios. Note that for an
individual scenario, the two-stage model is reformulated to a single-level
problem. In Step 4, the first-stage solution obtained from Step 2 is
aggregated to obtain the expected value $\bar{x}$. Step 5 calculates the value
of the multiplier $\eta_{s}$. In Step 8, the subproblems are solved, where
each subproblem is augmented with a linear term proportional to the multiplier
$\eta^{\tau-1}_{s}$ and a squared two norm term penalizing the difference of
$x$ from $\bar{x}^{\tau-1}$. Steps 9-10 are similar as Steps 4-5. The
algorithm terminates once all first-stage decisions $x_{s}$ converge to a
common $\bar{x}$.

Algorithm 1 The Two-Stage PH Algorithm

1:Initialization: Let $\tau:=0$.

2:For all $s\in S$, compute.

3:$x_{s}^{(\tau)}:=\operatorname*{arg\,min}_{x}\\{a^{T}x+b^{T}_{s}y_{s}:(x,y_{s})\in
Q_{s}\\}$.

4:$\bar{x}^{(\tau)}:=\sum_{\forall s\in S}Pr(s)x^{(\tau)}_{s}$.

5:$\eta^{(\tau)}_{s}:=\rho(x_{s}^{(\tau)}-\bar{x}^{(\tau)})$.

6:$\tau:=\tau+1$.

7:For all $s\in S$, compute.

8:$x_{s}^{(\tau)}:=\operatorname*{arg\,min}_{x}\\{a^{T}x+b^{T}_{s}y_{s}+\eta^{(\tau-1)}_{s}x+\frac{\rho}{2}\|x_{s}^{(\tau)}-\bar{x}^{(\tau)}\|^{2}:(x,y_{s})\in
Q_{s}\\}$.

9:$\bar{x}^{(\tau)}:=\sum_{\forall s\in S}Pr(S)x^{(\tau)}_{s}$.

10:$\eta^{(\tau)}_{s}:=\eta^{(\tau-1)}_{s}+\rho(x_{s}^{(\tau)}-\bar{x}^{(\tau)})$.

11:if $\sum_{\forall s\in
S}Pr(s)\|x_{s}^{(\tau)}-\bar{x}^{(\tau)}\|\leq\varepsilon$ then

12: Go to Step 5.

13:else

14: terminate.

15:end if

###  3.2 Convergence and Solution Validation

As shown in Algorithm 1, the convergence metric $g^{\tau}$ of progressive
hedging algorithm at each iteration $\tau$ is expressed as the deviation from
the mean summed across all first-stage variables $x_{s}(\tau)$ and the average
value of the first-stage variable $\bar{x}^{\tau}$ as follows:

$g^{\tau}=\sum_{s\in S}Pr(s)\|x_{s}(\tau)-\bar{x}^{\tau}\|$ (47)

Numerical results for convergence analysis are given in case study section. In
order to test the solution quality based on the limited generated damage
scenarios, we follow the suggestion from [20] and apply MRP to test the
stability and quality of the candidate solutions, as shown in Algorithm 2. MRP
is to repeat the procedure of generating $S$ scenarios and solving the
proposed model for $S$ times and construct the confidence interval (CI) for
the optimality gap. The detailed steps in MRP is shown in Algorithm 2, where
$\bar{G_{n}}(n_{g})$ is the gap estimate and $s^{2}_{G}(n_{g})$ is the sample
variance.

Algorithm 2 Multiple Replication Procedure

1:Input: Value $\alpha\in(0,1)$ (e.g., $\alpha=0.05$), sample size $n$,
replication size $n_{g}$ and a candidate solution $\hat{x}\in X$.

2:Output: Approximate $(1-\alpha)$ as the level confidence interval on
$\mu_{\hat{x}}$.

3:For $k=1,2,...,n_{g}$.

4:Sample i.i.d. observations $\zeta^{k_{1}},\zeta^{k_{2}},...,\zeta^{k_{n}}$
from the distribution of $\zeta$.

5:Solve ($\rm SP_{n}$) using $\zeta^{k_{1}},\zeta^{k_{2}},...,\zeta^{k_{n}}$
to obtain $x^{k*}_{n}$.

6:$G^{k}_{n}(\hat{x}):=n^{-1}\sum_{j=1}^{n}(f(\hat{x},\zeta^{kj})-f(x^{k*}_{n},\zeta^{kj}))$.

7:$\bar{G_{n}}(n_{g}):=\frac{1}{n_{g}}\sum_{k=1}^{n_{g}}G^{k}_{n}(\hat{x})$.

8:$s^{2}_{G}(n_{g}):=\frac{1}{n_{g}-1}\sum_{k=1}^{n_{g}}(G^{k}_{n}(\hat{x})-\bar{G_{n}}(n_{g}))^{2}$.

9:$\epsilon:=t_{n_{g}-1,\alpha}S_{G}(n_{g})/\sqrt{n_{g}}$.

10:Obtain one-sided CI on $[0,\bar{G_{n}}(n_{g})+\epsilon_{g}]$.

##  4 Case Study

In this section, a large-scale system is used as a test case to verify the
scalability and effectiveness of the two-stage stochastic pre-event resource
allocation model. This large-scale system consists of 3 existing test systems,
EPRI ckt5, ckt7 systems [23], and IEEE 8500 bus system [24], Following the
suggestions from [25], the cost parameters in the simulation are $C^{\rm
D}=14\$/kWh$, $C^{\rm SW}=8\$$, $C^{\rm F}=1\$/L$ and $r^{\rm F}=0.3L/kWh$.
The stochastic models and algorithms are implemented using the PySP package in
Pyomo [26]. IBM’s CPLEX 12.6 mixed-integer solver is used to solve all
subproblems. The experiments were performed on Iowa State University’s Condo
cluster, whose individual blades consist of two 2.6 GHz 8-Core Intel E5-2640
v3 processors and 128 GB of RAM.

###  4.1 Pre-Event Preparation Results

In this case, we have included 9 depots that are hosting a total of 27 crews,
9 dispatchable DGs, 8 MEGs, 3 MESs, 123 switches, 5 small PVs, 15 large PVs,
and 12 ESSs. The 9 DGs are rated at 300 kW and 250 kVAr. The 5 small PVs are
rated at 11kW$\sim$22kW. The 15 large PVs are rated at 500 kW. The 12 ESSs are
rated at 500 kW/ 3500 kWh. The pre-event preparation model of the large-scale
system is solved in 10.2 hours with 10 damage scenarios. The first-stage
decision variables (locations of MEGs, MESs and crews) are shown in Figure 5.
27 crews are allocated to 9 different depots. The value inside the crew depot
in Figure 5 represents the number of crews assigned to that depot. Areas with
large number of crews indicates that the lines in the area have high damage
probabilities.

Figure 5: Resource allocation of large-system with the proposed model.

As discussed in Section 3.2, we use the convergence metric to evaluate the
convergence speed of the proposed model. At the same time, we also compare the
computational speed with and without a soft-start solution. Soft-start
solution means that the previous computed solution in other instance will be
used as the starting point. The comparison result is shown in Figure 6. If the
convergence metric reaches the convergence threshold 0.01, the algorithm will
stop and obtain the optimal solution. The instance with soft-start solution
converges at 57 iteration and takes 10.2 hours. The case without soft-start
solution converges after 100 iteration and takes 24.3 hours. To test the
solution quality with MRP, based on the limited generated damage scenarios,
the one-sided CI of the obtained solution is $[0,12.48\%]$. This small gap
indicates that our solution is stable and of high quality.

Figure 6: The convergence metric comparison with and without soft-start
solutions.

To evaluate the performance of the developed pre-event preparation model, the
model is compared to a base model. The base case is generated by the following
steps: (i) one MEG are prepositioned at the substations. (ii) Extra MEG are
prepositioned at high-priority loads. (iii) PV and ESS are not considered.
(iv) Fuel is allocated to the MEGs such that they can operate for at least 24
hours. (v) Crews are allocated evenly between depots. In this work, we
calculate average outage duration by dividing the sum of outage durations for
the loads with the number of loads. To compare the performance of the proposed
model and the base model, we generate a random scenario and test the response
of the system. The generated scenario has 103 damaged lines and they were
aggregated to 34 damaged areas in Figure 7. Each circle represents the repair
time needed for the specific damaged area considering all the aggregated
damaged lines.

Figure 7: Aggregated damaged areas.

The comparison between the base model and the proposed method is shown in
Figure. 8. In the base model, the total restored energy is 231,422.38 kWh and
the average outage duration is 14.69 hours. In the proposed method, the total
restored energy is 291,727.48 kWh and the average outage duration is 11.28
hours. Therefore, approximately 20.67% more loads are served by the proposed
method and the outage duration decreased by 30.22%.

Figure 8: Comparison between base model and proposed method.

###  4.2 Impacts of Solar PV on System Resilience

To show the advantages of the PV systems, we test the response of the system
with the proposed method and different PV penetration levels. As discussed in
Section 2.3.1, we consider three types of PV: (i) Type I PV, which represents
residential PV panels and the rated capacity is assumed to be 6 kW; (ii) Type
II PV, which represents mid-size PV systems and the rated capacity is assumed
to be 48 kW; (iii) Type III PV, which represents large utility PV farm and the
rated capacity is assumed to be 2000 kW. Based on the number of different
types of PVs, we define 6 PV penetration levels as 9%, 27%, 45%, 63%, 81%, and
99%. The number of Type I, II and III PVs for each PV penetration levels is
summarized in Table 1. To better collaborate the setting of PV penetration,
the number of dispatchable DGs has been changed to 10 and the positions of
those DGs have been changed accordingly. The rest of case settings keep the
same.

Table 1: Number of different types of PV | PV Penetration
---
Percentage
| Type I
---
PV
| Type II
---
PV
| Type III
---
PV
9% | 8 | 1 | 1 |
27% | 24 | 4 | 3 |
45% | 40 | 7 | 5 |
63% | 63 | 9 | 7 |
81% | 72 | 12 | 9 |
99% | 88 | 15 | 11 |

Based on the results of Figure. 9, it can be observed that different PV
penetration levels have different allocation results of the flexible
resources, including the positions of MEGs, MESs and number of repair crews.

Figure 9: Pre-event resource allocation results with different PV penetration
levels.

Figure. 10 shows the percentage of power served during the event, and after
the repair process starts. Table 2 and Table 3 compare the amount of load
served and average outage duration with different levels of PV penetration.

Figure 10: Load served percentage comparison of proposed model with various PV
penetration level and base model.

Based on the results from Figure. 10, Table 2 and Table 3, it can be seen that
the penetration of PV contributes to enhancing system resilience.
Approximately 31.13% more loads are served compared to the base model when the
proposed method with 99% PV penetration is used. Also, the average outage
duration decreased by 31.12%. However, compared with the case of 81% PV
penetration level, the proposed method with 99% PV penetration does not have
significant improvement.

Table 2: The amount of load served and resilience improvement with different level of PV penetration | PV Penetration
---
Percentage
| Load Served
---
(kWh)
| Resilience Improvement
---
Percentage(%)
0 | 251,210.72 | -
9% | 318,668.37 | 26.85
27% | 335,525.77 | 33.56
45% | 336,710.74 | 34.04
63% | 344,588.22 | 37.17
81% | 360,668.04 | 43.57
99% | 364,785.93 | 45.21
Table 3: The amount of average outage duration and outage decreased percentage with different level of PV penetration | PV Penetration
---
Percentage
| Average Outage
---
Duration (hour)
| Outage Decreased
---
Percentage(%)
0 | 14.69 | -
9% | 12.33 | 16.07
27% | 11.72 | 20.22
45% | 11.65 | 20.69
63% | 11.21 | 23.69
81% | 10.45 | 28.86
99% | 10.12 | 31.11

##  5 Conclusion

In this paper, we develop a two-stage stochastic pre-event resource allocation
method for upcoming extreme events, which enables faster and more efficient
post-event restoration. The proposed pre-event method leverages the pre-
allocation of mobile resources, fuel resources and labor resources. It also
facilitates the benefits of distributed PV systems in resilience improvement
of distribution grids. According to the case studies, we have the following
observations: (i) Compared to the base model without pre-event resource
allocation, the proposed pre-event preparation model can serve more loads and
reduce the outage duration. (ii) Based on the response of the system with
different PV penetration levels, it can be observed that the proposed pre-
event preparation model with high PV penetration can further improve system
resilience and reduce the outage duration. Therefore, PV systems can play a
critical role in improving distribution grid resilience and further promote
the renewable energy deployment. (iii) By considering the trade-off between
solution accuracy and computation efficiency, the result of MRP indicates that
the proposed model’s solutions with a limited number of scenarios can be very
stable and of high quality. The scalability of the proposed pre-event
preparation model is verified with a large-scale system.

## Acknowledgement

This work was supported by the U.S. Department of Energy Wind Energy
Technologies Office under Grant DE-EE0008956.

## References

  * Salman et al. [2015] A. M. Salman, Y. Li, M. G. Stewart,  Evaluating system reliability and targeted hardening strategies of power distribution systems subjected to hurricanes,  Reliability Engineering & System Safety 144 (2015) 319–333.
  * of the President [2020] E. O. of the President, Economic benefits of increasing electric grid resilience to weather outages, Technical Report, White House Tech. Rep., 2020\.
  * Gholami et al. [2019] A. Gholami, T. Shekari, S. Grijalva,  Proactive management of microgrids for resiliency enhancement: An adaptive robust approach,  IEEE Trans. Sustain. Energy 10 (2019) 470–480.
  * Wang et al. [2017] C. Wang, Y. Hou, F. Qiu, S. Lei, K. Liu,  Resilience enhancement with sequentially proactive operation strategies,  IEEE Trans. Power Syst. 32 (2017) 2847–2857.
  * Panteli et al. [2017] M. Panteli, P. Mancarella, D. N. Trakas, E. Kyriakides, N. D. Hatziargyriou,  Metrics and quantification of operational and infrastructure resilience in power systems,  IEEE Trans. Power Syst. 32 (2017) 4732–4742.
  * Wang et al. [2004] S. Wang, B. R. Sarker, L. Mann, E. Triantaphyllou,  Resource planning and a depot location for electric power restoration,  Euro. J. Oper. Res. 155 (2004) 22–43.
  * Arif et al. [2018] A. Arif, Z. Wang, J. Wang, C. Chen,  Power distribution system outage management with co-optimization of repairs, reconfiguration, and DG dispatch,  IEEE Trans. Smart Grid 9 (2018) 4109–4118.
  * Arif et al. [2020] A. Arif, Z. Wang, C. Chen, B. Chen,  A stochastic multi-commodity logistic model for disaster preparation in distribution systems,  IEEE Trans. Smart Grid 11 (2020) 565–576.
  * Taheri et al. [2019] B. Taheri, A. Safdarian, M. Moeini-Aghtaie, M. Lehtonen,  Enhancing resilience level of power distribution systems using proactive operational actions,  IEEE Access 7 (2019) 137378–137389.
  * Lei et al. [2018] S. Lei, J. Wang, C. Chen, Y. Hou,  Mobile emergency generator pre-positioning and real-time allocation for resilient response to natural disasters,  IEEE Trans. Smart Grid 9 (2018) 2030–2041.
  * Lei et al. [2019] S. Lei, C. Chen, H. Zhou, Y. Hou,  Routing and scheduling of mobile power sources for distribution system resilience enhancement,  IEEE Trans. Smart Grid 10 (2019) 5650–5662.
  * Kim and Dvorkin [2019] J. Kim, Y. Dvorkin,  Enhancing distribution system resilience with mobile energy storage and microgrids,  IEEE Trans. Smart Grid 10 (2019) 4996–5006.
  * Samara et al. [2020] S. Samara, M. F. Shaaban, A. H. Osman,  Optimal management of mobile energy generation and storage systems,  IEEE Access 8 (2020) 203890–203900.
  * Belding et al. [2020] S. Belding, A. Walker, A. Watson, Will solar panels help when the power goes out?, Technical Report, National Renewable Energy Tech. Rep., 2020.
  * Rockafellar and Wets [1991] R. T. Rockafellar, R. J.-B. Wets,  Scenarios and policy aggregation in optimization under uncertainty,  Mathematics of operations research 16 (1991) 119–147.
  * Watson and Woodruff [2011] J.-P. Watson, D. L. Woodruff,  Progressive hedging innovations for a class of stochastic mixed-integer resource allocation problems,  Computational Management Science 8 (2011) 355–370.
  * Ma et al. [2018] S. Ma, B. Chen, Z. Wang,  Resilience enhancement strategy for distribution systems under extreme weather events,  IEEE Trans. Smart Grid 9 (2018) 1442–1451.
  * Ma et al. [2019] S. Ma, S. Li, Z. Wang, F. Qiu,  Resilience-oriented design of distribution systems,  IEEE Trans. Power Syst. 34 (2019) 2880–2891.
  * Chen et al. [2018] B. Chen, C. Chen, J. Wang, K. L. Butler-Purry,  Sequential service restoration for unbalanced distribution systems and microgrids,  IEEE Trans. Power Syst. 33 (2018) 1507–1520.
  * Arif et al. [2020] A. Arif, Z. Wang, C. Chen, J. Wang,  Repair and resource scheduling in unbalanced distribution systems using neighborhood search,  IEEE Trans. Smart Grid 11 (2020) 673–685.
  * Melhem et al. [2018] F. Y. Melhem, O. Grunder, Z. Hammoudan, N. Moubayed,  Energy management in electrical smart grid environment using robust optimization algorithm,  IEEE Trans. Industry Applications 54 (2018) 2714–2726.
  * Zhang et al. [2019] Q. Zhang, K. Dehghanpour, Z. Wang,  Distributed CVR in unbalanced distribution systems with PV penetration,  IEEE Trans. Smart Grid 10 (2019) 5308–5319.
  * EPRI [2019] EPRI, OPENDSS test circuits, Apr.2019. URL: https://sourceforge.net/p/electricdss/discussion/beginners.html.
  * Arritt and Dugan [2010] R. F. Arritt, R. C. Dugan,  The IEEE 8500-node test feeder,  in: IEEE PES T&D Conference, 2010, pp. 1–6.
  * Ma [2020] S. Ma,  Resilience-oriented design and proactive preparedness of electrical distribution system,  PhD Thesis (2020).
  * Hart et al. [2017] W. E. Hart, C. D. Laird, J.-P. Watson, D. L. Woodruff, G. A. Hackebeil, B. L. Nicholson, J. D. Siirola, Pyomo-optimization modeling in python, volume 67, Springer, 2017\.


11institutetext:  LESIA, Observatoire de Paris, Université PSL, CNRS, 5 Place
Jules Janssen, 92190 Meudon, France 22institutetext: Univ. Grenoble Alpes,
CNRS, IPAG, 38000 Grenoble, France 33institutetext: Pixyl S.A. La Tronche,
France

# Updated characterization of long-period single companion by combining radial
velocity, relative astrometry, and absolute astrometry

F. Philipot, Please send any requests to<EMAIL_ADDRESS>A.-M.
Lagrange 1122 P. Rubini 33 F. Kiefer 11 A. Chomez 11

(Received XXX / Accepted XXX)

###### Abstract

Context. Thanks to more than 20 years of monitoring, the radial velocity (RV)
method has detected long-period companions (P ¿ 10yr) around several dozens of
stars. Yet, the true nature of these companions remains unclear because of the
uncertainty as to the inclination of the companion orbital plane.

Aims. We wish to constrain the orbital inclination and the true mass of long-
period single companions.

Methods. We used a Markov Chain Monte Carlo (MCMC) fitting algorithm to
combine RV measurements with absolute astrometry and, when available, relative
astrometry data.

Results. We have lifted the sin(i) indetermination for seven long-period
companions. We find true masses in the planetary mass range for the candidate
planets detected in the following systems: Epsilon Indi A, HD 13931, HD
115954, and HD 222155. The mass of HD 219077 b is close to the deuterium-
burning limit and its nature is uncertain because of the imprecise mass of the
host star. Using additional RV measurements, we refine the orbital parameters
of HIP 70849 b and find a mass in the planetary range. By combining RV data
with absolute and relative astrometry, we significantly improve the
characterization of HD 211847 B and properly determine its mass, which appears
to be in the low-mass star range. This work illustrates how Gaia and Hipparcos
allow for the orbital properties and masses of long-period RV companions to be
further constrained.

###### Key Words.:

Techniques: radial velocities – Techniques: high angular resolution – Proper
motions – Stars: planetary systems – Stars: brown dwarfs – Stars: low-mass

##  1 Introduction

In the last decade, several long-period giant planets have been detected using
the radial velocity (hereafter RV) method thanks to the increasing temporal
baselines of different surveys (Mayor et al. (2011), Wittenmyer et al. (2020),
Rosenthal et al. (2021)). Yet, a precise determination of the orbital
parameters and mass of the planets is very difficult when the orbital period
is much larger than the RV time baseline. As a consequence, the radial
distribution of planets beyond 8-10 au – such as those found by Fernandes et
al. (2019) and Fulton et al. (2021) based on the results of the two long RV
surveys of Mayor et al. (2011) and Rosenthal et al. (2021), respectively – are
questionable. This unfortunately prevents an accurate comparison with
formation model outputs from being made.

Combining RV data with other methods such as relative or absolute astrometry
can, in principle, improve the orbital characterization of these companions.
Furthermore, it can also remove the uncertainty of the orbital inclination and
then allow us to determine the true mass of the planets.

Coupling RV data with relative astrometry from direct imaging (hereafter DI)
or interferometry has been, however, limited to very few cases since high-
contrast imaging (hereafter HCI) or interferometry observations favor young
systems to minimize the flux contrast between the star and its companion while
RV observations favor old and inactive stars which produce low RV jitters.
However, when possible, such a coupling is very efficient. An illustration is
the HD 7449 system for which the outer companion was first reported as a
planet candidate using only RV data (Mayor et al. (2011), Wittenmyer et al.
(2019)), and it was then identified as a low-mass star by combining RV data
with HCI observations (Rodigas et al. 2016).

In the 2000s, the combination of RV data and absolute astrometry, thanks to
the Fine-Guidance-Sensor onboard the Hubble Space Telescope, also allowed for
the inclination of a few stellar systems to be constrained and a few candidate
planets to be confirmed (Benedict et al. (2002), Benedict et al. (2006)),
while others were finally identified as brown dwarfs or low-mass stars (Bean
et al. (2007); Benedict et al. (2010)). Today, the position and proper motion
measurements obtained with the telescopes Hipparcos (Perryman et al. (1997),
van Leeuwen (2007)) and Gaia (Gaia Collaboration 2020) allow us to combine the
RV data and more precise absolute astrometry for a large number of systems.
Since the publication of the first Gaia data release (DR1), a few studies have
proven the efficiency of combining RV data with absolute and/or relative
astrometry to improve the constraints on the orbital parameters and mass of a
companion (Grandjean et al. (2019), Brandt et al. (2019), Damasso et al.
(2020), Lagrange et al. (2020), Nielsen et al. (2020), Venner et al. (2021),
Brandt et al. (2021a), Brandt et al. (2021b), Kiefer et al. (2021), Li et al.
(2021), Feng et al. (2022)).

In this paper, we focus on seven long-period single companions detected by the
RV method, and combine the available RV data with Hipparcos and Gaia early
data release 3 (hereafter EDR3) absolute astrometry and, when available,
relative astrometry, to improve the orbital parameters and determine the true
mass of these companions. In Section 2, we describe our target selection
method and present the RV, HCI, and astrometric data used in our study.
Section 3 presents the method used to perform the orbital fitting and, there,
we provide the new orbital parameters and mass found for each target. Finally,
we discuss the results in Section 4.

##  2 Target selection and data

###  2.1 Target selection

We first selected the planetary systems in the exoplanet.eu catalog (Schneider
et al. 2011) for which a single companion has been reported with a semi-major
axis greater than 5 au using the RV method. Twenty-five companions were found
with such criteria. For nine of them (HD 13724 B, HD 25015 b, HD 181234 b, and
HD 219828 B (Feng et al. 2022); HD 92987 B (Venner et al. 2021); HIP 36985 B
(Biller et al. 2022); and HD 98649 b, HD 196067 b, and HD 221420 B (Li et al.
2021)), the orbital parameters and the true mass have already been properly
determined in previous studies. We, therefore, do not consider them in the
present study.

We first discarded the HD 95872 system because no Hipparcos data were
available. We then discarded four systems for which the available RV time
series did not cover both extrema of the RV variations and the orbital period
could not be properly determined (HD 26161, HD 120066, HD 150706, and HD
213472). In those four cases, the combination of RV and absolute astrometry
did not allow us to constrain the orbital parameters, the orbital inclination,
or the true mass of the companion. Finally, we discarded four systems for
which the orbital period was well covered by the RV data, but the coupling
with absolute astrometry did not allow us to constrain the orbital inclination
(HD 136925, HD 190984, HD 220773, and HD 238914). Indeed, the variations in
position and acceleration of the proper motion of these stars were too small
to constrain the orbital inclination of the companion due to the low mass of
the companion (¡ $\sim$2 MJup) and/or the distance of the system.

Thus, we were left with seven systems for which the addition of absolute
astrometry and/or new RV measurements and/or relative astrometry measurements
allowed us to determine the exact nature of the companion: Epsilon Indi Ab, HD
13931, HD 1159554, HD 211847, HD 219077, HD 222155, and HIP 70849. For three
of these companions (Epsilon Indi Ab, HD 211847 B, and HD 219077 b), a first
estimation of their orbital inclination and true mass has been obtained by
combining RV data and absolute astrometry. Yet, thanks to additional data or
more precise astrometric measurements, we obtained more precise and
significantly different results from those reported in the previous studies
for six of these companions. In the case of HD 219077 b, the differences were
mainly found for the mass of the companion. They are mainly due to the
uncertainties as to the host star’s mass.

###  2.2 RV data

The RV data used in this study were obtained with different spectrographs
between 1997 and 2021. The HARPS (Mayor et al. 2003) data were taken from the
ESO archives; the ELODIE (Baranne et al. 1996) and SOPHIE (Perruchot et al.
2008) data were retrieved from the OHP archives; and the CORALIE (Queloz et
al. 1999), the HIRES (Vogt et al. 1994), the UVES (Dekker et al. 2000), the
AAT (Diego et al. 1990), the CES (Enard 1982) long camera (LC), and the CES
very long camera (VLC) data were taken from the literature.

As instrument upgrades can lead to new RV offsets, the same instrument before
and after a major upgrade is considered as two different instruments.
Consequently, HARPS data obtained before and after the optical fiber upgrade
in 2015 (Lo Curto et al. 2015) are referred to as H03 and H15, respectively.
The SOPHIE data obtained before and after the spectrograph upgrade in 2011
(Bouchy et al. 2013) are referred to as SOPHIE and SOPHIE+, respectively. The
HIRES data obtained before and after the upgrade of the spectrograph in 2004
(Tal-Or et al. 2019) are referred to as Hir94 and Hir04, respectively.
Finally, the CORALIE spectrograph had two major upgrades in 2007 (Ségransan et
al. 2010) and in 2014. The data obtained before 2007 and after 2014 are
referred to as C98 and C14, respectively, and the data obtained between 2007
and 2014 are referred to as C07.

###  2.3 Direct imaging data

In three cases, HCI data are available in the SPHERE archive and can provide
relative astrometry. The three targets were observed in angular (and spectral)
differential imaging (A(S)DI, Marois et al. (2006)) using the telescope in
pupil tracking mode. The standard observing mode of SPHERE was used, with
IRDIS (Dohlen et al. 2008) dual band images at H2 and H3 (K1 and K2,
respectively) and IFS (Claudi et al. 2008) data covering the YJ (YJH,
respectively) bands. The observing log is given in Table.1. Whenever possible,
the robust PACO A(S)DI algorithm (Flasseur et al. (2020a), Flasseur et al.
(2020b)) was used. The processing step takes advantage of the developments
made to the COBREX data center pipeline (the prereduction improvement as well
as the improvement of the detection capability of PACO). If a dataset did not
sufficiently cover the field-of-view rotation to apply ADI-based algorithms,
the SPECAL (Galicher et al. 2018) No-ADI algorithm was used.

In those three systems, only one companion was detected (HD 211847 B). The
detected companion is characterized in Table.2. No detection above 5$\sigma$
was found around HD 219077. Six sources were detected around HIP 70849 but,
given their position in a color-magnitude diagram and their separations, they
are likely background sources.

Table 1: Observing logs.

STAR DATE OBS FILTER DIT(s)$\times$Nframea $\Delta$ PA ($\degree$)a Seeing
(”)b Airmassb $\tau_{0}$ (ms)a,b Program ID HIP 70849 2015-05-05 DB_H23 64x64
40.4 1.12 1.08 0.0012 095.C-0298(A) HD 211847 2015-06-10 DB_K12 64x8 5.8 1.31
1.05 0.0025 095.C-0476(A) HD 219077 2015-06-09 DB_K12 8x64 3.3 1.25 1.31
0.0026 095.C-0476(A)

Notes : a: DIT is the detector integration time per frame. $\Delta$ PA is the
amplitude of the parallactic rotation. $\tau_{0}$ corresponds to the coherence
time. b are values extracted from the updated DIMM information, averaged over
the sequence.

Table 2: Relative astrometry for HD 211847 companion.

Sources JD - 2400000 IRDIS filter SEP (mas) PA ($\deg$) HD 211847 B 57183.39
K12 220$\pm$4.73 194.5$\pm$2.23

Notes : The errors displayed are 1$\sigma$. The relative astrometry combines
the astrometry measured in the dual bands.

###  2.4 Absolute astrometry

We used measurements from Hipparcos obtained around epoch 1991.25 and from
Gaia EDR3 (Gaia Collaboration et al. (2016), Gaia Collaboration et al. (2021))
obtained around epoch 2016.0. For each target, the stellar acceleration was
determined from the proper motion and the position values were measured by
Hipparcos and Gaia with an interval of about 25 years. We considered the
proper motion values published by Brandt (2021) in the Hipparcos-Gaia Catalog
of Accelerations (HGCA). Moreover, a more accurate tangential proper motion
($\mu_{Hip-EDR3}$) was estimated by the difference between the position
measurements obtained by Hipparcos and Gaia divided by the time interval
between the two measurements ($\sim$25 years). The proper motion values used
for each star are given in table 4.

##  3 Updated orbital parameters and mass

###  3.1 Orbit fitting

Orbits were fitted using a custom MCMC tool, based on the emcee 3.0 python
package (Foreman-Mackey et al. 2013). It uses a mixture of move functions
(such as the differential evolution move function) to alleviate potential
multimodality issues. The Hipparcos/Gaia data processing uses the HTOF package
(Brandt et al. 2021) and borrows large sections of the orvara code (Brandt et
al. 2021) for the likelihood computation. The HTOF package (Brandt et al.
2021) was used to fit the intermediate astrometric data (IAD) from Hipparcos,
based on the 1997 (Esa 1997) and 2007 (van Leeuwen 2007) reductions and from
Gaia, thanks to the Gaia Observation Forecast Tool (GOST) which allowed us to
obtain the estimated Gaia observations and scan angles for each target, in
order to reproduce proper motion and position of each observation. Using the
Hipparcos and Gaia positions and the temporal baseline, the algorithm derived
a tangential proper motion value that allowed us to better constrain the
orbital fit when combined with RV data.

We considered ten free parameters for each system: the semi-major axis (a),
the eccentricity, the orbital inclination (i), the host star mass, the
companion mass, the longitude of ascending node ($\Omega$), the argument of
periastron ($\omega$), the phase, a stellar jitter, and the distance of the
system. In addition, to combine data from different instruments, we added an
instrumental offset for each instrument as a free parameter of the model (see
above). Finally, we considered uniform priors for all fitting parameters,
except for the host star mass and the distance of the system for which we
considered Gaussian priors.

###  3.2 Results

####  3.2.1 Epsilon Indi A

Epsilon Indi is a triple system with a $0.76\pm 0.04$ M☉, K2V star (Epsilon
Indi A) and a binary composed of a $75.0\pm 0.8$ MJup, T1.5 brown dwarf
(Epsilon Indi B) and a $70.1\pm 0.7$ MJup, T6 brown dwarf (Epsilon Indi C)
separated by about 2.6 au (Dieterich et al. 2018). The projected separation
between the binary brown dwarfs and the star is about 1460 au. Combining RV
data and absolute astrometry based on Hipparcos and the Gaia data release 2
(DR2) measurements, Feng et al. (2019) reported a giant planet with a semi-
major axis of $11.55^{+0.98}_{-0.86}$ au, a mass of $3.25^{+0.39}_{-0.65}$
MJup, an inclination of $64.25^{+13.80}_{-6.09}$°, and an eccentricity of
$0.26^{+0.07}_{-0.03}$. Yet, the Gaia EDR3 proper motion and position
measurements are significantly more precise compared to the Gaia DR2
measurement and they significantly improve the characterization of Epsilon
Indi Ab.

We used 4278 RV measurements obtained with the HARPS spectrograph between 2003
and 2016111The 3636 RV data obtained between Julian days 2455790 and 2455805
were obtained to study high-frequency oscillations of the star. These data
were measured with high cadence, which led to a significantly lower signal-to-
noise ratio compared to the other data., 163 RV measurements obtained with the
UVES spectrograph between 1996 and 2017, 72 RV measurements obtained with the
LC spectrograph between 1992 and 1997, and 53 RV measurements obtained with
the VLC spectrograph between 2000 and 2006. We also combined these RV data
with absolute astrometry based on Hipparcos and the Gaia EDR3 measurements
(Fig.1). We found significantly different orbital parameters with a semi-major
axis of $8.8^{+0.2}_{-0.1}$ au, a mass of $3.0\pm 0.1$ MJup, an inclination of
$91^{+4}_{-5}$°, and an eccentricity of $0.48\pm 0.01$. It is important to
note that if we consider only the 539 HARPS RV data with a signal-to-noise
ratio greater than 110 and thus remove the high cadence observations made in
August 2011, we find similar solutions.

Figure 1: Orbital fits for Epsilon Indi Ab. Top: Fit of the Epsilon Indi A RV
data corrected from the instrumental offset (V0). Bottom: Fit of the Epsilon
Indi A astrometric acceleration in right ascension (left) and declination
(right). The black points correspond to the measurements obtained with
Hipparcos (1991.25) and Gaia EDR3 (2016.0). In each plot, the black curve
shows the best fit. The color bar indicates the log likelihood of the
different fits plotted.

####  3.2.2 HD 13931

HD 13931 is a $1.02\pm 0.05$ M☉(Rosenthal et al. 2021), G0V star. Based on 66
RV measurements obtained with the HIRES spectrograph between 1998 and 2019,
Rosenthal et al. (2021) reported a giant planet with a semi-major axis of
$5.323\pm 0.091$ au, a minimum mass of $1.911^{+0.077}_{-0.076}$ MJup , and an
eccentricity of $0.02^{+0.021}_{-0.014}$.

We combined these RV data with absolute astrometry (Fig.2). As the RV baseline
is much larger than the orbital period, the orbital parameters are well-
constrained. As expected, we found a semi-major axis and an eccentricity very
close to those reported by Rosenthal et al. (2021) with a = $5.33\pm 0.09$ au
and e ¡ 0.04. Using, in addition, the absolute astrometry, we found an orbital
inclination of either $39^{+13}_{-8}$° or $141^{+9}_{-18}$° and a true mass of
$3.1^{+0.8}_{-0.7}$ MJup.

Figure 2: Orbital fits for HD 13931 b. Top: Fit of the HD 13931 RV data
corrected from the instrumental offset (V0). Bottom: Fit of the HD 13931
astrometric acceleration in right ascension (left) and declination (right).
The black points correspond to the measurements obtained with Hipparcos
(1991.25) and Gaia (2016.0). In each plot, the black curve shows the best fit.
The color bar indicates the log likelihood of the different fits plotted.

####  3.2.3 HD 115954

HD 115954 is a $1.18\pm 0.06$ M☉, G0V star (Demangeon et al. 2021). Based on
four RV measurements obtained with the ELODIE spectrograph between 2004 and
2005 and 45 RV measurements obtained with the SOPHIE spectrograph between 2009
and 2018, Demangeon et al. (2021) reported a giant planet with a semi-major
axis of $5.00^{+1.3}_{-0.36}$ au, a minimum mass of $8.29^{+0.75}_{-0.58}$
MJup , and an eccentricity of $0.487^{+0.095}_{-0.041}$.

We combined these RV data with the absolute astrometry data (Fig.3). We found
a semi-major axis compatible with Demangeon et al. (2021) with a =
$4.5^{+0.2}_{-0.1}$ au and an eccentricity of $0.46\pm 0.03$. Using, in
addition, the absolute astrometry, we found an orbital inclination of
$92^{+17}_{-16}$° and a true mass of $8.5^{+0.6}_{-0.4}$ MJup.

Figure 3: Orbital fits for HD 115954 b. Top: Fit of the HD 115954 RV data
corrected from the instrumental offset (V0). Bottom: Fit of the HD 115954
astrometric acceleration in right ascension (left) and declination (right).
The black points correspond to the measurements obtained with Hipparcos
(1991.25) and Gaia (2016.0). In each plot, the black curve shows the best fit.
The color bar indicates the log likelihood of the different fits plotted.

####  3.2.4 HD 211847

HD 211847 is a $0.94\pm 0.04$ M☉, G5V star (Sahlmann et al. 2011). Sahlmann et
al. (2011) reported a brown dwarf candidate orbiting around HD 211847 based on
31 RV measurements obtained with the CORALIE spectrograph between 2002 and
2009. However, only one minimum of the HD 211847 B RV curve was covered by the
dataset. Thus, the orbital parameters and minimum mass reported in this study
are poorly constrained. Using the Levenberg-Marquardt method, they found
ranges corresponding to a $3\sigma$ confidence interval for the semi-major
axis, the eccentricity, and the minimum mass of 4.6–42 au, 0.48–0.95, and
16.3–24.3 MJup , respectively. Moutou et al. (2017) obtained one HCI detection
with SPHERE of HD 211847 B for a projected separation of 11.3 au. Using the
BT-Settl models (Allard 2014), they fit the HD 211847 B spectrum and found a
low stellar mass of $155\pm 9$ MJup assuming an age of 3 Gyr for the host
star. Based on the result of Sahlmann et al. (2011), Moutou et al. (2017)
estimated the inclination of the companion orbit to be around seven°.
Recently, combining the CORALIE RV measurement and the absolute astrometry,
Feng et al. (2022) reported HD 211847 B as a brown dwarf with a semi-major
axis of $4.514^{+0.458}_{-0.287}$ au, a mass of $55.32^{+1.335}_{-18.48}$
MJup, an inclination of $163.649^{+36.239}_{-5.017}$°, and an eccentricity of
$0.419^{+0.035}_{-0.064}$.

We combined the RV dataset used by Sahlmann et al. (2011), the relative
astrometry observation obtained with SPHERE in June 2015, and the absolute
astrometry (Fig.4). Adding one relative astrometry observation allowed us to
properly constrain the orbital parameters and the mass of HD211847 B with
results significantly different from those reported by Feng et al. (2022). We
found a semi-major axis of $6.78\pm 0.08$ au and an eccentricity of
$0.59^{+0.01}_{-0.02}$. Using, in addition, the absolute astrometry, we found
an orbital inclination of $172.3^{+0.05}_{-0.04}$° and a true mass of $148\pm
5$ MJup. We note that by taking only RV data and absolute astrometry into
account, we found very poorly constrained solutions with large uncertainties
as to the semi-major axis (16 - 30 au) and mass (80 - 140 MJup). Moreover, the
solutions found are not in agreement with those reported by Feng et al. (2022)
or with the solutions found when adding the HCI data.

Figure 4: Orbital fits for HD 211847 B. Top left: Fit of the HD 211847 RV data
corrected from the instrumental offset (V0). Top right: Fit of HD 211847
relative astrometry data. The red cross corresponds to the measurement
obtained with SPHERE. Bottom: Fit of the HD 211847 astrometric acceleration in
right ascension (left) and declination (right). The black points correspond to
the measurements obtained with Hipparcos (1991.25) and Gaia (2016.0). In each
plot, the black curve shows the best fit. The color bar indicates the log
likelihood of the different fits plotted.

####  3.2.5 HD 219077

Based on 63 CORALIE RV measurements obtained between 1999 and 2012 and 30
HARPS RV measurements obtained between 2003 and 2012, Marmier et al. (2013)
reported a very eccentric giant planet with a semi-major axis of $6.22\pm
0.09$ au and a minimum mass of $10.39\pm 0.09$ MJup. It is important to note
that the RV data used by Marmier et al. (2013) are not publicly available.
Based on 72 pieces of RV data obtained with the AAT spectrograph between 1998
and 2015, Kane et al. (2019) reported slightly different properties for HD
219077 with a semi-major axis of $7.03^{+0.20}_{-0.21}$ au and a minimum mass
of $13.40^{+0.76}_{-0.78}$ MJup. These differences are probably due to the
different assumptions on the mass of the star. Indeed, Marmier et al. (2013)
used the values reported by Hipparcos (M⋆= $1.05\pm 0.02$ M☉) while Kane et
al. (2019) used the values reported in Valenti & Fischer (2005) (M⋆= $1.51\pm
0.13$ M☉). Recently, Feng et al. (2022) combined the AAT RV measurements used
by Kane et al. (2019) and 33 HARPS RV measurements obtained between 2003 and
2012 with the absolute astrometry based on Hipparcos and the Gaia EDR3
measurements and found a semi-major axis close to Marmier et al. (2013), an
orbital inclination of $90.178^{+9.527}_{-9.462}$°, and a true mass of
$9.620^{+1.001}_{-0.733}$ MJup. The prior on the mass of the star is not
given.

For this study, as the data used by Marmier et al. (2013) are not available,
we considered the HARPS and AAT RV measurements used on Feng et al. (2022) and
the 65 CORALIE RV measurements available on the DACE
archive222https://dace.unige.ch obtained between 1999 and 2012333The CORALIE
RV data available on DACE and the HARPS RV data available on the ESO archive
cover the same time base as those used by Marmier et al. (2013).. For the mass
of the star, we considered the value given by Kervella et al. (2022) based on
the Gaia DR3 results (M⋆= $1.15\pm 0.06$ M☉). We combined the RV data with the
absolute astrometry (Fig.5). We found a semi-major axis and an eccentricity
close to those reported in the previous studies with a = $6.4\pm 0.1$ au and e
= $0.769\pm 0.002$ and an orbital inclination close to that of Feng et al.
(2022) with either i = $83\pm 3$° or i = $97\pm 3$°. Considering the star mass
found by Kervella et al. (2022), we found a planetary mass at $11.3\pm 0.4$
MJup. However, considering that the mass used in Valenti & Fischer (2005)
would lead to a mass close to the deuterium-burning limit, $M_{b}$ = $13.6\pm
0.5$ MJup. Due to the uncertainties on the mass of the host star, it is not
possible to determine the exact nature of HD 219077 b.

Figure 5: Orbital fits for HD 219077 b. Top: Fit of the HD 219077 RV data
corrected from the instrumental offset (V0). Bottom: Fit of the HD 219077
astrometric acceleration in right ascension (left) and declination (right).
The black points correspond to the measurements obtained with Hipparcos
(1991.25) and Gaia (2016.0). In each plot, the black curve shows the best fit.
The color bar indicates the log likelihood of the different fits plotted.

####  3.2.6 HD 222155

HD 222155 is a $1.13\pm 0.11$ M☉, G2V star (Boisse et al. 2012). Based on 44
RV measurements obtained with the ELODIE spectrograph between 1997 and 2005
and 67 RV measurements obtained with the SOPHIE spectrograph between 2007 and
2011, Boisse et al. (2012) reported a giant planet with a semi-major axis of
$5.1^{+0.6}_{-0.7}$ au, a minimum mass of $1.90^{+0.67}_{-0.53}$ MJup , and an
eccentricity of $0.38^{+0.28}_{-0.32}$.

We considered 31 additional pieces of SOPHIE RV data obtained between 2011 and
2016. We combined the RV data with the absolute astrometry (Fig.6). We found
orbital parameters within the error bars associated with the values found by
Boisse et al. (2012) with a = $4.7\pm 0.1$ au and e = $0.34\pm 0.09$. As the
RV baseline is now much larger than the orbital period, the orbital parameters
are better constrained. Using, in addition, the absolute astrometry, we found
an orbital inclination of either $66^{+14}_{-11}$° or $115^{+13}_{-16}$° and a
true mass of $2.1^{+0.3}_{-0.2}$ MJup.

Figure 6: Orbital fits for HD 222155 b. Top: Fit of the HD 222155 RV data
corrected from the instrumental offset (V0). Bottom: Fit of the HD 222155
astrometric acceleration in right ascension (left) and declination (right).
The black points correspond to the measurements obtained with Hipparcos
(1991.25) and Gaia (2016.0). In each plot, the black curve shows the best fit.
The color bar indicates the log likelihood of the different fits plotted.

####  3.2.7 HIP 70849

HIP 70848 is a $0.63\pm 0.03$ M☉, K7V star (Ségransan et al. 2011). Ségransan
et al. (2011) reported the first detection of HIP 70849 b based on 18 RV
measurements obtained with the HARPS spectrograph between 2006 and 2010.
However, only one minimum of the HIP 70849 b RV curve was covered by the
dataset. The observations carried out by Ségransan et al. (2011) led to poorly
constrained orbital parameters and minimum mass. Using a genetic algorithm
followed by MCMC simulations, they reported a semi-major axis between 4.5 and
36 au, a minimum mass between 3 and 15 MJup , and an eccentricity between 0.47
and 0.96 with ranges corresponding to a $3\sigma$ confidence interval.

We considered 39 additional pieces of HARPS RV data obtained between 2011 and
2021. We combined the RV data with the absolute astrometry (Fig.7). With these
additional observations, the dataset then covered two minimum and one maximum
of the RV curve of HIP 70849 b and this allowed us to properly constrain the
properties of the companion. We found a semi-major axis of
$3.99^{+0.06}_{-0.07}$ au and an eccentricity of $0.65^{+0.02}_{-0.01}$.
Using, in addition, the absolute astrometry, we found an orbital inclination
of $96\pm 16$° and a true mass of $4.5^{+0.4}_{-0.3}$ MJup.

Figure 7: Orbital fits for HIP 70849 b. Top: Fit of the HIP 70849 RV data
corrected from the instrumental offset (V0). Bottom: Fit of the HIP 70849
astrometric acceleration in right ascension (left) and declination (right).
The black points correspond to the measurements obtained with Hipparcos
(1991.25) and Gaia (2016.0). In each plot, the black curve shows the best fit.
The color bar indicates the log likelihood of the different fits plotted.
Table 3: Summary of posteriors obtained with our MCMC algorithm.

Parameter Eps Ind A HD 13931 HD 115954 HD 211847 HD 219077 HD 222155 HIP 70849
a (au) ${8.8}_{-0.1}^{+0.2}$ $5.33\pm 0.09$ ${4.5}_{-0.1}^{+0.2}$ $6.78\pm
0.08$ $6.4\pm 0.1$ $4.7\pm 0.1$ ${3.99}_{-0.07}^{+0.06}$ Period (days)
${10932}_{-228}^{+266}$ ${4442}_{-46}^{+49}$ ${3258}_{-190}^{+179}$
${6199}_{-46}^{+52}$ ${5514}_{-39}^{+44}$ ${3470}_{-106}^{+102}$ $3649\pm 18$
Eccentricity $0.48\pm 0.01$ ¡ 0.04 $0.46\pm 0.03$ $0.769\pm 0.002$
${0.59}_{-0.02}^{+0.01}$ $0.34\pm 0.09$ ${0.65}_{-0.01}^{+0.02}$ Inclination
(°) $91^{+4}_{-5}$ $39^{+13}_{-8}$ or $141^{+9}_{-18}$ ${59}_{-4}^{+5}$ or
$127\pm 4$ ${172.3}_{-0.4}^{+0.5}$ $83\pm 3$ or $97\pm 3$ ${66}_{-11}^{+14}$
or ${115}_{-16}^{+13}$ $96\pm 16$ Mass (MJup ) $3.0\pm 0.1$
$3.1^{+0.8}_{-0.7}$ ${8.5}_{-0.4}^{+0.6}$ $148\pm 5$ $11.3\pm 0.4$
${2.1}^{+0.3}_{-0.2}$ ${4.5}_{-0.3}^{+0.4}$ $\Omega$ (°) $58\pm 5$
${343}_{-19}^{+17}$ or ${110}_{-24}^{+19}$ ${211}_{-28}^{+25}$ $184\pm 5$
${135}_{-21}^{+38}$ or ${347}_{-30}^{+28}$ ${264}_{-33}^{+34}$ or
${180}_{-35}^{+34}$ $35\pm 6$ $\omega$ (°) $85\pm 3$ 74 - 227
${173}_{-8}^{+7}$ ${168}_{-4}^{+5}$ $56.2\pm 0.4$ 153 - 217 $182\pm 1$ Phase
$0.37\pm 0.01$ 0.38 - 0.86 $0.40_{-0.06}^{+0.09}$ $0.420\pm 0.004$
$0.354_{-0.003}^{+0.002}$ $0.98_{-0.13}^{+0.04}$ $0.745_{-0.006}^{+0.007}$
Jitter (m/s) $3.37\pm 0.04$ ${2.9}\pm 0.3$ $6.6_{-1.1}^{+1.3}$
${11.0}_{-2.0}^{+2.8}$ ${4.3}_{-0.2}^{+0.3}$ ${12.7}_{-0.8}^{+0.9}$
${6.7}_{-0.7}^{+0.8}$ H03 = $-39972\pm 1$ Hir94 = $-13\pm 1$ ELODIE =
${-14757}_{-41}^{+52}$ C98 = $6907\pm 17$ C98 = ${-30867}_{-1}^{+2}$ ELODIE =
$-3999\pm 3$ H03 = $53\pm 1$ Instrumental H15 = ${-39952}_{-1}^{+2}$ Hir04 =
$-8.0\pm 0.5$ SOPHIE = ${-14768}_{-8}^{+9}$ C07 = $6849\pm 22$ C07 =
${-30864}_{-9}^{+8}$ SOPHIE = ${-3950}_{-2}^{+3}$ H15 = ${65}_{-1}^{+2}$
offset (m/s) UVES = $-5\pm 1$ SOPHIE+ = ${-14743}_{-3}^{+4}$ H03 =
${-30830}_{-2}^{+1}$ SOPHIE+ = ${-3923}_{-3}^{+4}$ LC = ${-39978}_{-1}^{+2}$
AAT = $-68\pm 1$ VLC = ${-39976}_{-1}^{+2}$

Notes: The results were obtained by combining RV, absolute astrometry, and,
when available, relative astrometry. We provide $68\%$ confidence intervals
for each parameter and the median is only given when the probability
distribution has a profile close to a Gaussian distribution.

##  4 Summary and concluding remarks

Combining RV measurements from various spectrographs with absolute astrometry
based on Hipparcos and Gaia EDR3 data and, when available, relative
astrometry, we determined the orbital parameters and, in particular, the
orbital inclination and the true mass of seven long-period single companions
detected by the RV method. Figure 8 summarizes the true mass and the semi-
major axis of the companions and compares them with the previous estimations.
Clearly, Gaia EDR3 data allow for a better determination of these companions’
orbital parameters and mass. All of these companions have true masses of 2
MJup or more and orbit between 3.9 - 9 au from their stars. Absolute
astrometry would probably help to determine the true mass of planets with a
period larger than the duration of Gaia DR3 observations (P ¿ $\sim$1000 d)
down to 1 MJup , provided the RV variations are well covered and the
variations in position and acceleration of the proper motion of the star are
large enough. In practice, in most cases, when the period is not well
constrained by the RV data, the impact of the coupling of RV data with
absolute astrometry is more limited. An illustration of this is the case of HD
211847 B for which, by combining RV data that cover only a minimum of the RV
variations with absolute astrometry, Feng et al. (2022) reported a mass of
about 55 MJup, corresponding to a brown dwarf. Yet, HCI revealed a companion
and the fit of the RV and the relative and absolute astrometry leads to a mass
of about 150 MJup instead.

We conclude that Gaia/Hipparcos can help to further constrain the orbital
parameters of long-period RV planets, providing good coverage of the RV
variations is available. Otherwise, additional information is needed, such as
relative astrometry, provided by DI or interferometry.

Figure 8: Update of the orbital parameters and masses of the seven analyzed
systems thanks to the combination of absolute astrometric and RV data and,
when available, absolute astrometry data. For each system, a dotted line
between two solutions was drawn to allow for the different solutions obtained
to be compared.

###### Acknowledgements.

This study was funded by a grant from PSL/OCAV. This project has also received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (COBREX; grant agreement n°
885593). This work presents results from the European Space Agency (ESA) space
mission Gaia. Gaia data are being processed by the Gaia Data Processing and
Analysis Consortium (DPAC). Funding for the DPAC is provided by national
institutions, in particular the institutions participating in the Gaia
MultiLateral Agreement (MLA). The Gaia mission website is
https://www.cosmos.esa.int/gaia. The Gaia archive website is
https://archives.esac.esa.int/gaia. This publications makes use of the The
Data & Analysis Center for Exoplanets (DACE), which is a facility based at the
University of Geneva (CH) dedicated to extrasolar planets data visualisation,
exchange and analysis. DACE is a platform of the Swiss National Centre of
Competence in Research (NCCR) PlanetS, federating the Swiss expertise in
Exoplanet research. The DACE platform is available at https://dace.unige.ch.
This research has made use of the SIMBAD database and VizieR catalogue access
tool, operated at CDS, Strasbourg, France. Based on data retrieved from the
SOPHIE archive at Observatoire de Haute-Provence (OHP), available at
atlas.obs-hp.fr/sophie. Based on spectral data retrieved from the ELODIE
archive at Observatoire de Haute-Provence (OHP). Based on observations
collected at the European Southern Observatory under ESO programme(s)
183.C-0972(A), 108.22KV.001, 108.22KV.002, 072.C-0488(E), 108.222V.001,
090.C-0421(A), 085.C-0019(A), 087.C-0831(A), 095.C-0551(A), 0100.C-0097(A),
0103.C-0432(A), 0101.C-0379(A), 096.C-0460(A), 0102.C-0558(A), 093.C-0409(A),
098.C-0366(A), 099.C-0458(A), 089.C-0732(A), 091.C-0034(A), 196.C-1006(A),
0102.C-0558(A), 106.21R4.001, 183.C-0972(A), 075.C-0332(A), 192.C-0852(A),
190.C-0027(A), 091.C-0936(A), 085.C-0063(A), 086.C-0284(A), 076.C-0155(A),
077.C-0101(A), 082.C-0212(A).

## References

  * Allard (2014) Allard, F. 2014, in Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham, Vol. 299, 271–272
  * Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373
  * Bean et al. (2007) Bean, J. L., McArthur, B. E., Benedict, G. F., et al. 2007, AJ, 134, 749
  * Benedict et al. (2010) Benedict, G. F., McArthur, B. E., Bean, J. L., et al. 2010, AJ, 139, 1844
  * Benedict et al. (2002) Benedict, G. F., McArthur, B. E., Fredrick, L. W., et al. 2002, AJ, 124, 1695
  * Benedict et al. (2006) Benedict, G. F., McArthur, B. E., Gatewood, G., et al. 2006, AJ, 132, 2206
  * Biller et al. (2022) Biller, B. A., Grandjean, A., Messina, S., et al. 2022, A&A, 658, A145
  * Boisse et al. (2012) Boisse, I., Pepe, F., Perrier, C., et al. 2012, A&A, 545, A55
  * Bouchy et al. (2013) Bouchy, F., Díaz, R. F., Hébrard, G., et al. 2013, A&A, 549, A49
  * Brandt et al. (2021a) Brandt, G. M., Brandt, T. D., Dupuy, T. J., Li, Y., & Michalik, D. 2021a, AJ, 161, 179
  * Brandt et al. (2021b) Brandt, G. M., Dupuy, T. J., Li, Y., et al. 2021b, AJ, 162, 301
  * Brandt et al. (2021) Brandt, G. M., Michalik, D., Brandt, T. D., et al. 2021, The Astronomical Journal, 162, 230
  * Brandt et al. (2021) Brandt, G. M., Michalik, D., Brandt, T. D., et al. 2021, AJ, 162, 230
  * Brandt (2021) Brandt, T. D. 2021, ApJS, 254, 42
  * Brandt et al. (2019) Brandt, T. D., Dupuy, T. J., & Bowler, B. P. 2019, AJ, 158, 140
  * Brandt et al. (2021) Brandt, T. D., Dupuy, T. J., Li, Y., et al. 2021, The Astronomical Journal, 162, 186
  * Claudi et al. (2008) Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70143E
  * Damasso et al. (2020) Damasso, M., Sozzetti, A., Lovis, C., et al. 2020, A&A, 642, A31
  * Dekker et al. (2000) Dekker, H., D’Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 534–545
  * Demangeon et al. (2021) Demangeon, O. D. S., Dalal, S., Hébrard, G., et al. 2021, A&A, 653, A78
  * Diego et al. (1990) Diego, F., Charalambous, A., Fish, A. C., & Walker, D. D. 1990, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 1235, Instrumentation in Astronomy VII, ed. D. L. Crawford, 562–576
  * Dieterich et al. (2018) Dieterich, S. B., Weinberger, A. J., Boss, A. P., et al. 2018, ApJ, 865, 28
  * Dohlen et al. (2008) Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70143L
  * Enard (1982) Enard, D. 1982, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 331, Instrumentation in Astronomy IV, 232–242
  * Esa (1997) Esa, . 1997, VizieR Online Data Catalog, I/239
  * Feng et al. (2019) Feng, F., Anglada-Escudé, G., Tuomi, M., et al. 2019, MNRAS, 490, 5002
  * Feng et al. (2022) Feng, F., Butler, R. P., Vogt, S. S., et al. 2022, ApJS, 262, 21
  * Fernandes et al. (2019) Fernandes, R. B., Mulders, G. D., Pascucci, I., Mordasini, C., & Emsenhuber, A. 2019, ApJ, 874, 81
  * Flasseur et al. (2020a) Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020a, A&A, 637, A9
  * Flasseur et al. (2020b) Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020b, A&A, 634, A2
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Fulton et al. (2021) Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, ApJS, 255, 14
  * Gaia Collaboration (2020) Gaia Collaboration. 2020, VizieR Online Data Catalog, I/350
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Galicher et al. (2018) Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92
  * Grandjean et al. (2019) Grandjean, A., Lagrange, A. M., Beust, H., et al. 2019, A&A, 629, C1
  * Kane et al. (2019) Kane, S. R., Dalba, P. A., Li, Z., et al. 2019, AJ, 157, 252
  * Kervella et al. (2022) Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7
  * Kiefer et al. (2021) Kiefer, F., Hébrard, G., Lecavelier des Etangs, A., et al. 2021, A&A, 645, A7
  * Lagrange et al. (2020) Lagrange, A. M., Rubini, P., Nowak, M., et al. 2020, A&A, 642, A18
  * Li et al. (2021) Li, Y., Brandt, T. D., Brandt, G. M., et al. 2021, AJ, 162, 266
  * Lo Curto et al. (2015) Lo Curto, G., Pepe, F., Avila, G., et al. 2015, The Messenger, 162, 9
  * Marmier et al. (2013) Marmier, M., Ségransan, D., Udry, S., et al. 2013, A&A, 551, A90
  * Marois et al. (2006) Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556
  * Mayor et al. (2011) Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv e-prints, arXiv:1109.2497
  * Mayor et al. (2003) Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20
  * Moutou et al. (2017) Moutou, C., Vigan, A., Mesa, D., et al. 2017, A&A, 602, A87
  * Nielsen et al. (2020) Nielsen, E. L., De Rosa, R. J., Wang, J. J., et al. 2020, AJ, 159, 71
  * Perruchot et al. (2008) Perruchot, S., Kohler, D., Bouchy, F., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70140J
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49
  * Queloz et al. (1999) Queloz, D., Mayor, M., Weber, L., et al. 1999, arXiv e-prints, astro
  * Rodigas et al. (2016) Rodigas, T. J., Arriagada, P., Faherty, J., et al. 2016, ApJ, 818, 106
  * Rosenthal et al. (2021) Rosenthal, L. J., Fulton, B. J., Hirsch, L. A., et al. 2021, ApJS, 255, 8
  * Sahlmann et al. (2011) Sahlmann, J., Ségransan, D., Queloz, D., et al. 2011, A&A, 525, A95
  * Schneider et al. (2011) Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79
  * Ségransan et al. (2011) Ségransan, D., Mayor, M., Udry, S., et al. 2011, A&A, 535, A54
  * Ségransan et al. (2010) Ségransan, D., Udry, S., Mayor, M., et al. 2010, A&A, 511, A45
  * Tal-Or et al. (2019) Tal-Or, L., Trifonov, T., Zucker, S., Mazeh, T., & Zechmeister, M. 2019, MNRAS, 484, L8
  * Valenti & Fischer (2005) Valenti, J. A. & Fischer, D. A. 2005, VizieR Online Data Catalog, J/ApJS/159/141
  * van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
  * Venner et al. (2021) Venner, A., Vanderburg, A., & Pearce, L. A. 2021, AJ, 162, 12
  * Vogt et al. (1994) Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 2198, Instrumentation in Astronomy VIII, ed. D. L. Crawford & E. R. Craine, 362
  * Wittenmyer et al. (2019) Wittenmyer, R. A., Clark, J. T., Zhao, J., et al. 2019, MNRAS, 484, 5859
  * Wittenmyer et al. (2020) Wittenmyer, R. A., Wang, S., Horner, J., et al. 2020, MNRAS, 492, 377

##  Appendix A Proper motion values

Table 4: Proper motion values from HGCA.

Star Eps ind A HD 13931 HD 115954 HD 211847 HD 219077 HD 222155 HIP 70849
$\mu_{Hip}^{\alpha}$ (mas/yr) $3964.6\pm 0.4$ $98.8\pm 0.8$ $-74.6\pm 0.7$
$56.1\pm 1.1$ $477.5\pm 0.4$ $195.3\pm 0.5$ $-47.0\pm 2.1$
$\mu_{Hip}^{\delta}$ (mas/yr) $-2537.1\pm 0.4$ $-183.7\pm 0.6$ $20.9\pm 0.6$
$14.1\pm 0.8$ $-424.9\pm 0.4$ $-117.8\pm 0.6$ $-203.3\pm 1.9$
$\mu_{EDR3}^{\alpha}$ (mas/yr) $3966.7\pm 0.1$ $98.57\pm 0.04$ $-74.77\pm
0.02$ $44.43\pm 0.03$ $478.30\pm 0.03$ $195.31\pm 0.02$ $-44.05\pm 0.02$
$\mu_{EDR3}^{\delta}$ (mas/yr) $-2536.2\pm 0.1$ $-183.41\pm 0.04$ $21.49\pm
0.02$ $9.66\pm 0.04$ $-424.43\pm 0.04$ $-117.34\pm 0.02$ $-201.58\pm 0.3$
$\mu_{Hip-EDR3}^{\alpha}$ (mas/yr) $3965.02\pm 0.01$ $98.45\pm 0.03$
$-74.79\pm 0.02$ $47.90\pm 0.04$ $478.36\pm 0.01$ $195.25\pm 0.02$ $-44.43\pm
0.06$ $\mu_{Hip-EDR3}^{\delta}$ (mas/yr) $-2537.25\pm 0.01$ $-183.51\pm 0.02$
$-21.41\pm 0.02$ $-10.44\pm 0.03$ $-424.40\pm 0.01$ $-117.39\pm 0.02$
$-202.05\pm 0.04$

Notes: $\mu_{Hip}$ corresponds to the proper motion obtained by Hipparcos.
$\mu_{EDR3}$ corresponds to the proper motion obtained by Gaia EDR3.
$\mu_{Hip-EDR3}$ corresponds to the proper motion obtained by the Hipparcos-
Gaia EDR3 positional difference.

##  Appendix B MCMC priors

Table 5: Priors considered for each free parameter.

Parameter Eps ind A HD 13931 HD 115954 HD 211847 HD 219077 HD 222155 HIP 70849
a (au) [1,20] [1,10] [1,10] [1,100] [1,10] [1,10] [1,10] Eccentricity [0,0.95]
[0,0.95] [0,0.95] [0,0.95] [0,0.95] [0,0.95] [0,0.95] Inclination (°) [0,180]
[0,180] [0,180] [0,180] [0,180] [0,180] [0,180] Mass (MJup ) [1,20] [1,20]
[1,20] [1,500] [1,20] [1,20] [1,20] $\Omega$ (°) [0,360] [0,360] [0,360]
[0,360] [0,360] [0,360] [0,360] $\omega$ (°) [0,360] [0,360] [0,360] [0,360]
[0,360] [0,360] [0,360] Phase [0,1] [0,1] [0,1] [0,1] [0,1] [0,1] [0,1] Jitter
(m/s) [0,10] [0,10] [0,10] [0,20] [0,10] [0,20] [0,10] Star mass (M☉) $0.76\pm
0.04$ $1.02\pm 0.05$ $1.18\pm 0.06$ $0.94\pm 0.05$ $1.15\pm 0.06$ $1.13\pm
0.11$ $0.63\pm 0.03$ Distance (pc) $3.622\pm 0.004$ $44.2\pm 1.4$ $218\pm 2$
$50.6\pm 3.3$ $29.3\pm 0.2$ $49.1\pm 0.9$ $24.0\pm 0.7$ H03: [-41,-39] Hir94:
[-1,1] ELODIE: [-15,-13] C98: [5,7] C98: [-31,-29] ELODIE: [-5,-3] H03: [-1,1]
Instrumental H15: [-41,-39] Hir04: [-1,1] SOPHIE: [-15,-13] C07: [5,7] C07:
[-31,-29] SOPHIE: [-5,-3] H03: [-1,1] offset (km/s) UVES: [-1,1] SOPHIE+:
[-15,-13] H03: [-31,-29] SOPHIE+: [-5,-3] LC: [-41,-39] AAT: [-1,1] VLC:
[-41,-39]

##  Appendix C MCMC results

Figure 9: Corner plot of the posteriors’ fit of Epsilon Indi A combined RV and
absolute astrometry. An offset of 39.9 km/s was added to V0, V1, V3, and V4 to
improve readability.  Figure 10: Corner plot of the posteriors’ fit of HD
13931 combined RV and absolute astrometry.  Figure 11: Corner plot of the
posteriors’ fit of HD 115954 combined RV and absolute astrometry. An offset of
14 km/s was added to V0, V1, and V2 to improve readability.  Figure 12: Corner
plot of the posteriors’ fit of HD 211847 combined RV, relative astrometry, and
absolute astrometry. An offset of 6 km/s was subtracted to V0 and V1 to
improve readability.  Figure 13: Corner plot of the posteriors’ fit of HD
219077 combined RV and absolute astrometry. An offset of 30 km/s was added to
V0, V1, and V2 to improve readability.  Figure 14: Corner plot of the
posteriors’ fit of HD 222155 combined RV and absolute astrometry. An offset of
3.9 km/s was added to V0, V1, and V2 to improve readability.  Figure 15:
Corner plot of posteriors fit of HIP 70849 combined RV and absolute
astrometry.