# Learning-based State Estimation in Distribution Systems with Limited Real-
Time Measurements

J. G. De la Varga, , S. Pineda, , J. M. Morales , and Á. Porras This work was
supported by the Spanish Ministry of Science and Innovation
(AEI/10.13039/501100011033) through project PID2020-115460GB-I00. The work of
J. G. De la Varga was supported by the Spanish Ministry of Science and
Innovation training program for PhDs with fellowship number PRE2021-098958.
_Corresponding author: Salvador Pineda._ The authors are with Oasys Research
Group, University of Málaga Spain, 29071 Málaga, Spain (e-mail:
spineda@uma.es).

###### Abstract

The task of state estimation in active distribution systems faces a major
challenge due to the integration of different measurements with multiple
reporting rates. As a result, distribution systems are essentially
unobservable in real time, indicating the existence of multiple states that
result in identical values for the available measurements. Certain existing
approaches utilize historical data to infer the relationship between real-time
available measurements and the state. Other learning-based methods aim to
estimate the measurements acquired with a delay, generating pseudo-
measurements. Our paper presents a methodology that utilizes the outcome of an
unobservable state estimator to exploit information on the joint probability
distribution between real-time available measurements and delayed ones.
Through numerical simulations conducted on a realistic distribution grid with
insufficient real-time measurements, the proposed procedure showcases superior
performance compared to existing state forecasting approaches and those
relying on inferred pseudo-measurements.

###### Index Terms:

Real-time observability, active distribution networks, machine learning,
pseudo-measurements, state estimation.

##  I Introduction

The task of state estimation plays a critical role in power systems as it
enables accurate monitoring and assessment of the system’s operating
conditions. However, in distribution systems, state estimation faces a major
challenge due to the inherent lack of real-time observability, as certain
measurements can only be acquired with a delay [1]. Furthermore, an electric
power system may lose part of its input data due to topology changes, sensor
failures or communication errors [2], malicious attacks [3], and electrical
blackouts [4]. As a result, distribution systems are practically unobservable
in real time, indicating the existence of multiple states that result in
identical values for the available measurements. This effect is even
aggravated because distribution systems are historically measurement-scarce
and have to deal with multiple reporting rates and possible asynchronization
[5]. Meanwhile, the progress in distributed energy resources (DERs) and the
necessity for actively managing their operational tasks have resulted in the
emergence of active distribution networks (ADNs) [6]. Central to this paradigm
shift is the imperative for real-time monitoring, analysis, and control of
grid operations, for which the state estimation is a key component.

Under these circumstances, supervisory control and data acquisition (SCADA)
systems and feeder terminal units (FTUs) play an essential role in providing
real-time measurements to the distribution management system. Additionally,
micro-phasor measurement units ($\mu$-PMUs) [7] and advanced metering
infrastructure [8], which include smart meters (SMs) [9], also have great
potential to provide data for state estimation. Given the heterogeneous nature
and notably slow reporting times of this measurement set, most existing works
discard dynamic state estimation and favor a static state estimation framework
over a forecasting-aided one [5].

Hence, to obtain valid estimations in a real-time fashion, a static state
estimation process must incorporate additional information as available
measurements are usually insufficient to deliver observability. Many existing
approaches aim at generating the pseudo-measurements required to make the
system observable by leveraging known information [10], either employing
probabilistic and statistical methods or learning-based techniques. One
example of the former group is [11], where the first approach to incorporate
loads as pseudo-measurements in distribution system state estimation is
introduced. Another example of statistical methods to generate pseudo-
measurements is proposed in [12], where spatial and temporal correlations are
accounted for by means of a first-order vector autoregressive model. The
latter group of learning-based techniques has aroused great interest in
academia and industry, mainly through advances in Deep Learning and Neural
Networks (NN) models. In [13], an artificial NN is used to generate power
injection pseudo-measurements from a few real measurements.

However, previous models have some limitations, namely, the need for abundant
high-quality training data, the black-box nature, and the risk of physically
infeasible results [14]. These drawbacks could be partly averted by linking
physics-based prior knowledge to the models. A physics-conditioned
optimization problem is formulated in [15] to find the optimal choice for
injection pseudo-measurements given a few real-time measures, while
considering that the probability distribution for such unknown power
injections follows a previously trained Generative Adversarial Network model.
Another example can be found in [16], where an autoencoder framework is
physics-enhanced for pseudo-measurement generation. This hybrid learning
method employs Deep NNs to learn the state correlations while considering
physical flows in a power system.

Alternatively, in contrast to augmenting conventional methods with pseudo-
measurements, some references propose directly forecasting the state using the
available information. Authors in [17] simulate the behavior of power systems
in two steps: first predicting the state based on previous estimations,
followed by a filtering of the new set of measurements; both models,
implemented by two distinct NN. From the same family but from a different
category, a Deep NN is used in [18] to learn the mappings between real-time
measurements and the state, in combination with a Bayesian state estimation
algorithm for improved bad-data detection. State forecasting can also be
useful to estimate the state of observable systems. For instance, the authors
of [19] train a shallow NN to learn how to initialize the Gauss-Newton
algorithm by mapping available measurements to a point near the actual latent
state, using historical or simulation data.

Likewise, integrating the physics into the models can also be used to improve
state forecasting. In [20], a graph neural network is used to reconstruct
physical parameters such as admittances of power lines without the need for
complete observability. Once the model is tuned, the state can be directly
predicted through available measurements. Another NN is introduced in [21],
whose physics-informed approach is included by mirroring the structure of a
regularized and ‘locally linearized’ state estimation algorithm, showcasing
better performance and efficient training. Lastly, [22] proposes an original
method that assumes that the grid state vector is a smooth graph signal with
respect to the system admittance matrix. Subsequently, the authors develop a
regularized weighted least squares state estimator that takes into account the
graph structure of the network and does not require full observability.

Against this background, this paper introduces a novel data-driven static
state estimation methodology specifically designed to deal with a measurement
set with multiple reporting rates and insufficient real-time measurements. One
of the key advantages of our proposal is its unique utilization of an
unobservable state estimator’s outcome to extract valuable information
regarding the joint probability distribution between real-time available
measurements and delayed measurements. This innovative approach combines the
benefits of physics-informed methods while offering a distinct advantage over
existing approaches by providing a simple, transparent, and replicable
solution. To demonstrate the effectiveness of our proposal, we conduct
numerical simulations on a realistic case study with limited real-time
measurements. The results of these simulations showcase the improvement of our
proposed methodology when compared to existing state forecasting approaches
and those that rely on inferred pseudo-measurements.

The rest of the paper is organized as follows. In Section II we present the
state estimation formulation and its particularities when it comes to
unobservable systems. Afterward, our proposal and its advantages are
described. In Section III we provide the numerical results of a case study to
compare the performance of the proposed methodology with existing ones.
Finally, Section IV provides concluding remarks.

##  II Methodology

In this section, we begin by introducing the conventional static state
estimation problem and outlining its key components. Subsequently, we delve
into the challenges that arise when the observability of the entire network is
not achieved, mainly focusing on state forecasting methods and pseudo-
measurement generation techniques. Finally, we present our proposal that
effectively improves upon existing methods by leveraging the information on
the unobservable states provided by the system model equations and enhancing
the accuracy of the overall estimation process.

Using a generic notation, the relationship between the measurement vector
$z\in\mathbb{R}^{m}$ and the system state $x\in\mathbb{R}^{n}$ is
characterized as:

$z=h(x)+e$ (1)

where function $h(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is the
measurement function that depends on the physical laws governing the
relationship between the measured magnitudes and the state variables, and
$e\in\mathbb{R}^{m}$ denotes the measurement error vector, whose components
are typically assumed to follow normal distributions with zero mean and be
independent and uncorrelated with respect to the state variables [2].

If the number of measurements is greater than the number of unknown state
variables, and these measurements are sufficiently redundant and diverse, the
system satisfies the observability requirements and the state estimation
problem can be formulated as the following Weighted Least Square (WLS)
optimization problem:

$\hat{x}\in\arg\min_{x}\left(z-h\left(x\right)\right)^{T}W\left(z-h\left(x\right)\right)$
(2)

where $\hat{x}$ is the estimated state vector, $(\cdot)^{T}$ is the matrix
transposition operation, and the weight matrix
$W=diag\\{\sigma^{-2}_{1},\ldots,\sigma^{-2}_{m}\\}$) represents the user’s
confidence in the measured data, where $\sigma^{2}_{m}$ represents the error
variance of the $m$-$th$ measurement. For ease of notation, the solution of
problem (2) is denoted as $\hat{x}=\texttt{wls}(z)$. In the general case,
problem (2) results in a nonconvex optimization problem without constraints,
which conventionally has been solved using an iterative Newton-Raphson
algorithm [23]. Alternatively, problem (2) can also be directly solved by the
interior-point methods implemented in nonlinear optimization solvers [24].

As discussed in Section I, even when the numbers and types of measures
captured by the sensors are enough to make the system fully observable, some
of the measurements are not available in real time due to delays or
communication issues. Consequently, we denote the set of available real-time
measures as $z^{a}\in\mathbb{R}^{m_{a}}$, and the set of delayed measures as
$z^{d}\in\mathbb{R}^{m_{d}}$, and split the measurement function as follows:

$\displaystyle z^{a}=h^{a}(x)+e^{a}$ (3a) $\displaystyle z^{d}=h^{d}(x)+e^{d}$
(3b)

where $h^{a}(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}^{m_{a}}$,
$h^{d}(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}^{m_{d}}$,
$e^{a}\in\mathbb{R}^{m_{a}}$, $e^{d}\in\mathbb{R}^{m_{d}}$, and
$m=m_{a}+m_{d}$. In this work, we assume that the number and type of measures
available in real time are not enough to make the system observable [1, 5].
Furthermore, we consider we have access to a collection of historical data
comprising both available and delayed measurements. We denote this data set as
$\\{(z^{a}_{k},z^{d}_{k})\;\forall k\in\mathcal{K}\\}$. Although model (2)
cannot be used to infer the system state in real time because it is not
observable, we can always perform a retrospective state estimation using all
the historical measurements and data collected. That is, for each instance $k$
the retrospective system state denoted as $\hat{x}^{r}_{k}$ can be computed as
$\hat{x}^{r}_{k}=\texttt{wls}(z^{a}_{k},z^{d}_{k})$. Therefore, we assume we
have access to the enlarged historical data set
$\\{(z^{a}_{k},z^{d}_{k},\hat{x}^{r}_{k})\;\forall k\in\mathcal{K}\\}$. Under
this framework, the main goal is to find a $\omega$-parameterized function
$f(\cdot;\omega)$ that provides, as accurately as possible, the retrospective
system state $\hat{x}^{r}_{k}$ out of the measurements available in real time
$z^{a}_{k}$. Formally, we can estimate this function as follows:

$\omega^{*}\in\arg\min_{\omega}\sum_{k\in\mathcal{K}}(\hat{x}^{r}_{k}-f(z^{a}_{k};\omega))^{2}$
(4)

Once parameters $\omega$ are properly tuned, the system state of a new
instance $k^{\prime}$ can be directly determined as

$\hat{x}_{k^{\prime}}=f(z^{a}_{k^{\prime}};\omega^{*})$ (5)

Depending on the characteristics of function $f$, different strategies can be
derived. This paper examines existing strategies that fall within this general
framework while introducing novel approaches that surpass them in terms of
estimation accuracy. In the simplest case, we can use the historical data to
infer the joint probability distribution of the available measurements
$z^{a}_{k}$ and the retrospective states $\hat{x}^{r}_{k}$ using conventional
supervised regression techniques. These approaches are usually known as _state
forecasting_ and multiple examples can be found in the literature,
particularly since the tremendous boom of NN for prediction tasks [17, 19].

In most cases, state forecasting methods consider standard functions $f$ that
either simplify the underlying physical relationship between state variables
and measurements in order to obtain interpretable models, or try to
approximate function $h$ with elaborated machine learning approaches. For
instance, if $f$ belongs to the family of linear functions, model (4) reduces
to a linear regression problem that can be solved analytically to determine
the optimal parameters $\omega^{*}$. However, the assumption of a linear
relationship between state variables and measurements is highly restrictive
and can significantly compromise the accuracy of this approach. On the other
end of the spectrum, function $f$ can be represented by a neural network,
providing a more flexible and adaptable modeling approach at the expense of
losing interpretability, increasing sensitivity to parameter tuning, or
augmenting the risk of overfitting. Finally, some recent research works
propose physics-informed learning models to forecast the system state taking
into account the underlying electricity network. However, the promising
results in [18] are derived from the careful tuning of some of the design
parameters for the specific systems studied, while the method in [21] exploits
a NN with an embedded ‘locally linearized’ representation of the state-
estimation equations, with the consequent approximation error.

Alternatively, other methodologies aim at inferring the joint probability
distribution between available and delayed measurements using a function
$g(\cdot;\lambda)$, parameterized on $\lambda$. For the same historical data
set $\\{(z^{a}_{k},z^{d}_{k},\hat{x}^{r}_{k})\;\forall k\in\mathcal{K}\\}$,
the optimal value of $\lambda$ can be determined as follows:

$\lambda^{*}\in\arg\min_{\lambda}\sum_{k\in\mathcal{K}}(z^{d}_{k}-g(z^{a}_{k};\lambda))^{2}$
(6)

For an unseen instance $k^{\prime}$, the delayed measurements can be estimated
in real time as $\hat{z}^{d}_{k^{\prime}}=g(z^{a}_{k^{\prime}};\lambda^{*})$.
Since $\hat{z}^{d}_{k^{\prime}}$ are not actual measurements provided by
sensors, they are known as pseudo-measurements in the technical literature
[25]. By leveraging techniques for generating pseudo-measurements, the system
becomes observable and its state can be determined as
$\hat{x}_{k^{\prime}}=\texttt{wls}(z^{a}_{k^{\prime}},\hat{z}^{d}_{k^{\prime}})$.
Therefore, using a common mathematical framework, the relation between the
system state and the available measurements is characterized as

$\hat{x}_{k^{\prime}}=\texttt{wls}(z^{a}_{k^{\prime}},g(z^{a}_{k^{\prime}};\lambda^{*}))$
(7)

Again, selecting the proper function $g$ and solving problem (6) are the main
drawbacks of this strategy. Moreover, akin to the state forecasting approach
discussed earlier, most existing methodologies generate pseudo-measurements
either with a simple physical-agnostic approach or with complex physics-
informed strategies as discussed in Section I.

In this paper, we present a methodology aimed at improving the performance of
state forecasting approaches and pseudo-measurement generation techniques in
cases where the available real-time measurements $z^{a}$ are insufficient to
ensure system observability. Specifically, we seamlessly integrate the
measurement function $h$ into these techniques to enhance the estimation
accuracy of the system state variables or the pseudo-measurements,
respectively. To do so, for each instance $k$, we formulate the following
optimization problem

$\tilde{x}_{k}\in\arg\min_{x}\left(z^{a}_{k}-h^{a}\left(x\right)\right)^{T}W^{a}\left(z^{a}_{k}-h^{a}\left(x\right)\right)$
(8)

where $W^{a}$ denotes the weight matrix of available measurements only. Since
the number of available measurements is lower than the number of state
variables, problem (8) has multiple solutions and its optimal value is always
0. The optimization solver, therefore, delivers one solution satisfying
$z^{a}_{k}=h^{a}(\tilde{x}_{k})$, and thus, _carrying_ critical information on
the measurement function $h$. Once problem (8) is solved, we can generate
pseudo-measurements as $\tilde{z}^{d}_{k}=h^{d}(\tilde{x}_{k})$. This
procedure guarantees that there exists a state, namely, $\tilde{x}_{k}$ that
conforms to _all_ the measurements (both physical and pseudo-). Or in other
words, this set of measurements is consistent with $h$. Conversely, if the
pseudo-measurements $\tilde{z}^{d}_{k}$ are generated using the physics-
agnostics methods previously discussed, there is no guarantee of the existence
of state variables $x_{k}$ that simultaneously satisfy
$z^{a}_{k}=h^{a}(x_{k})$ and $g(z^{a}_{k};\lambda^{*})=h^{d}(x_{k})$.
Additionally, the proposed scheme seems straightforward and transparent in
contrast to the intricate physics-informed methods previously analyzed [15] or
[21]. For notation purposes, we denote the solution of problem (8) as
$\tilde{x}=\texttt{un}(z^{a})$.

The pseudo-measurements $\tilde{z}^{d}_{k}=h^{d}(\texttt{un}(z^{a}_{i}))$ we
propose can be used in different ways. For instance, we can use them as
additional features or explanatory variables in the state forecasting
procedure formulated in (4) and (5). By doing so, the optimal parameters
$\omega$ are obtained as

$\omega^{*}\in\arg\min_{\omega}\sum_{k\in\mathcal{K}}(\hat{x}^{r}_{k}-f(z^{a}_{k},h^{d}(\texttt{un}(z^{a}_{k}));\omega))^{2}$
(9)

and the system state of a new instance $k^{\prime}$ is computed as

$\hat{x}_{k^{\prime}}=f(z^{a}_{k^{\prime}},h^{d}(\texttt{un}(z^{a}_{k^{\prime}}));\omega^{*})$
(10)

Note that the only difference between approaches (4)-(5) and (9)-(10) is that
the proposed one also considers as explanatory variables or features of the
regression model the pseudo-measurements $\tilde{z}^{d}_{k}$. In that sense,
despite its apparent complexity, problem (9) is a least-squared method with
extra features obtained by solving problem (8) for each instance $k$ with a
non-convex optimization solver, and evaluating the known function $h^{d}$ on
the obtained system state. However, unlike approach (4)-(5), the estimation
model (9)-(10) accounts for the measurement function $h$ through problem
$\texttt{un}(\cdot)$ and through function $h^{d}$. Besides, this enhanced
state forecasting approach is simple, direct and easy to interpret unlike the
physics-informed black-box methodologies proposed in the literature.

Alternatively, the pseudo-measurements $\tilde{z}^{d}_{k}$ can also be used as
features to learn the delayed measurements $z^{d}_{k}$ as follows

$\lambda^{*}\in\arg\min_{\lambda}\sum_{k\in\mathcal{K}}(z^{d}_{k}-g(z^{a}_{k},h^{d}(\texttt{un}(z^{a}_{k}));\lambda))^{2}$
(11)

Similarly to (9), problem (11) is a conventional least-squared optimization
model that can be easily solved once the proposed pseudo-measurements are
computed. However, unlike problem (6), the one we propose incorporates
physical information about the underlying electricity network. For a new
instance $k^{\prime}$, the pseudo-measurements required to make the system
observable are obtained as

$\hat{z}^{d}_{k^{\prime}}=g(z^{a}_{k^{\prime}},h^{d}(\texttt{un}(z^{a}_{k^{\prime}}));\lambda^{*})$
(12)

While the results obtained from this enhanced pseudo-measurement generation
procedure can serve various purposes, their most intuitive application is as
measurements in the conventional state estimation problem (2). Consequently,
the state variables can be determined as

$\hat{x}_{k^{\prime}}=\texttt{wls}(z^{a}_{k^{\prime}},g(z^{a}_{k^{\prime}},h^{d}(\texttt{un}(z^{a}_{k^{\prime}}));\lambda^{*}))$
(13)

where the error variance $\sigma$ of the pseudo-measurements are obtained
through the residuals of the regression model (11).

In summary, this section describes four models, namely, (5), (7), (10), (13),
which basically involve mapping functions between the available measurements
$z^{a}$ and the state variables $x$. On the one hand, approaches (5) and (10)
represent state forecasting models that directly determine the system state
using the available measurements. However, while (5) relies solely on observed
data processed through physics-agnostic methods, the enhanced proposed model
(10) incorporates the system’s physical equations through the measurement
function $h$. On the other hand, approaches (7) and (13) aim at generating
pseudo-measurements to ensure system observability and determine the system
state using standard estimation techniques. Correspondingly, while model (7)
utilizes pseudo-measurements generated purely through physics-agnostic
techniques, the proposed methodology (13) does incorporate the physical
behavior of the system. Furthermore, in contrast to alternative physics-
informed methodologies prevalent in the technical literature, the models
presented in this paper exhibit a straightforward, direct approach that
facilitates effortless implementation and interpretation.

##  III Numerical Results

In the following, we provide results from a series of numerical experiments
designed to test the performance of our state estimation approaches against
alternative methods. We start by describing the setup of these experiments.

###  III-A Simulation setup

In the series of experiments, we utilize the well-known IEEE 33-bus
distribution system test case [26]. The data for this system is obtained from
the Matpower test cases library [27]. We have chosen the fully connected
configuration of this system for a partially meshed topology, which consists
of 33 nodes and 37 lines. Fig. 1 depicts the distribution network and the
location of the measuring devices for some of the analysis discussed in this
work.

Figure 1: IEEE-33 bus test case network and measurement set location.

To simulate a diverse set of operating conditions in ADNs with a high
penetration of renewables, prosumers, storage, and other DERs [28], the
voltage magnitude ($V_{i}$) is randomly chosen within a range of 0.95 p.u. to
1.05 p.u. Similarly, the voltage phase angle ($\theta_{i}$) is selected
between $-15\degree$ and $+15\degree$ [29]. In state estimation, the
measurement functions $h$ relate the set of measurable outputs to the system’s
state variables ($V_{i},\theta_{i}$) as follows

$\displaystyle h_{V_{i}}=V_{i}$ (14a) $\displaystyle
h_{\theta_{i}}=\theta_{i}$ (14b) $\displaystyle
h_{P_{ij}}=V_{i}V_{j}\left(G_{ij}\cos{\theta_{ij}}+B_{ij}\sin{\theta_{ij}}\right)-G_{ij}V_{i}^{2}$
(14c) $\displaystyle
h_{Q_{ij}}=V_{i}V_{j}\left(G_{ij}\sin{\theta_{ij}}-B_{ij}\cos{\theta_{ij}}\right)+V_{i}^{2}\left(B_{ij}-b_{ij}^{S}/2\right)$
(14d) $\displaystyle
h_{P_{i}}=V_{i}\sum_{j}V_{j}\left(G_{ij}\cos{\theta_{ij}}+B_{ij}\sin{\theta_{ij}}\right)$
(14e) $\displaystyle
h_{Q_{i}}=V_{i}\sum_{j}V_{j}\left(G_{ij}\sin{\theta_{ij}}-B_{ij}\cos{\theta_{ij}}\right)$
(14f)

where $G_{ij}$ and $B_{ij}$ are, respectively, the real and imaginary parts of
the bus admittance matrix, $b_{ij}^{S}$ is the shunt susceptance of the line
$ij$, and ${\theta_{ij}=\theta_{i}-\theta_{j}}$. The system measurements are
obtained by adding a Gaussian noise $\mathcal{N}(0,0.01^{2})$ to power
injections and flows and $\mathcal{N}(0,0.001^{2})$ to voltage magnitudes and
angles, as in [19].

In state estimation for unobservable systems, measurement functions (14)
become particularly critical. For instance, while each voltage phasor
measurement relates solely to an individual state variable, power flow
measurements implicate the four state variables of the two buses of the line,
and power injection measurements involve the state variables of all
interconnected buses. Therefore, having access to measurements with functions
that entail a high number of state variables offers richer information and
provides more clues about the system’s overall state, aiding in reconstructing
the unobservable variables.

The measurement set illustrated in Fig. 1 is detailed as follows.

  * •

Available measurements in real time ($m_{a}$):

    * –

One FTU provides measurements of voltage magnitude and outgoing active and
reactive power flows in the substation (red star).

    * –

Fifteen active and reactive power injections are available through SCADA
measurements (red squares).

    * –

Five PMUs provide voltage magnitude and angle measurements (red diamonds).

  * •

Delayed measurements ($m_{d}$):

    * –

Eighteen Smart Meters placed in the remaining nodes provide power injection
measurements (green circles). As previously stated, they are not available in
real time due to their slower reporting times [5].

In real time we only have access to the set $m_{a}=3+2\times(15+5)=43$ while
the measurements of $m_{d}=2\times(18)=36$ are delayed. Since $n=65$ (we fix
$\theta_{1}$ to 0), the system is unobservable in real time and its state
cannot be estimated through conventional techniques. However, the
retrospective state estimation is determined assuming that we have both sets
of available and delayed measurements, that is, $m=79$ measurements.

In Section II we describe two learning-based procedures to forecast the system
state using function $f$ and other two approaches to generate pseudo-
measurements through function $g$. Due to its simplicity and explainability,
we compare in this case study the performance of these four approaches
assuming that functions $f$ and $g$ belong to the family of linear functions.
To learn these functions, we randomly generate 10 000 instances. This set is
then divided into two subsets: a training dataset comprising 8 000 instances
and a test dataset consisting of 2 000 instances. By splitting the data in
this manner, we can effectively train and evaluate the performance of the four
approaches discussed.

The state estimation problems (2) and (8) are modeled in Pyomo 6.4.2 [30]
running in Python 3.10.13 and solved with Ipopt 3.12.8 [31]. In all cases, the
initial solution of the algorithm is the flat state, that is, all voltage
magnitudes and angles are set to 1 p.u. and $0\degree$, respectively. The
adjustment of the parameters of functions $f$ and $g$ are determined using the
linear regression function from the Scikit-learn package [32].

In order to demonstrate the performance of the proposed methodology, we
compare in this section the computational results of the following nine
approaches:

  * -

Benchmark (BN): This approach assumes that all measurements are available and
the system is then observable. Under the considered setup, however, this
situation is unrealistic and only used for benchmarking purposes.

  * -

Unobservable (UN): This approach consists of the solution of the optimization
problem (8) delivered by the non-convex optimization solver. Recall that, as
previously discussed, problem (8) has multiple solutions.

  * -

Flat (FL): This approach boils down to fixing the state variables to the flat
voltage profile.

  * -

State forecasting (SF): This strategy approximates the relationship between
the available measures and the state variables using a linear function and
therefore, implies the solution of (4) and (5).

  * -

Enhanced state forecasting (SF∗): This approach is similar to the previous one
but adds the additional enhanced features, i.e., $h^{d}(\texttt{un}(z^{a}))$,
described in Section II and therefore, it implies the solution of (9) and
(10).

  * -

1-NN: This is a naive approach that infers the system state as that
corresponding to the nearest neighbor measured as the Euclidean distance
between the available measurements.

  * -

20-NN: This is another naive forecasting approach where the system state is
computed as the average of the state variables of the 20 nearest neighbors
measured as the Euclidean distance between the available measurements.

  * -

Pseudo-measurements (PM): This approach first generates pseudo-measurements
through (6) and (7) and then solves the state estimation problem (2).

  * -

Enhanced pseudo-measurements (PM∗): This approach also generates pseudo-
measurements through (11) and (12), using the same additional features as SF∗.
Then it obtains the system state with (13).

###  III-B Enhancing effect of the proposed methodology

Before presenting the computational results of the series of experiments, we
illustrate in Fig. 2 the enhancing effect of the pseudo-measurements
$\tilde{z}^{d}$ obtained by solving the unobservable estimation problem (8).
In particular, we plot the measurement error of the active and reactive power
injection at node 19 for methods UN, PM, PM∗, 1-NN, and 20-NN in the test
dataset, with respect to the measurements $z^{d}$. For the sake of
illustration, the density plots have been depicted using a Gaussian kernel
density estimate with a factor of 5 for the standard deviation of the
smoothing kernel. As observed, although model (8) has multiple solutions, the
pseudo-measurements derived from the solution provided by the solver incur a
lower variance in measurement residuals than some naive learning techniques
such as 1-NN and 20-NN and in this node, even for the method PM. Additionally,
using the yielded pseudo-measurements as features improves the performance of
the regression model that aims at learning the delayed measures as a function
of the available ones. For both the active and reactive power, we can observe
that the residual probability distribution provided by PM∗ has a lower
variance than that obtained by PM.

The enhancing effect of the proposal is also shown in Table I, which contains
the mean and the standard deviation of the residuals of all delayed active
power injections. The results in this table indicate that the pseudo-
measurements obtained through the UN approach exhibit biased residuals.
However, using these pseudo-measurements as features in the regression model
of the PM∗ approach effectively rectifies this issue, yielding unbiased
estimations with the lowest variance for all delayed measurements
corresponding to active power injections. This behavior is consistently
replicated in all pseudo-measurements, but for brevity, only active power is
presented.

The reduction in variance obtained through method UN is obviously reliant on
the available measurement set. Solely relying on voltage phasor measurements
(magnitude and angle) provides limited information. These measurements only
constrain the corresponding state variables (bus voltages) to match their
measured values, leaving the remaining unobservable state variables
undetermined. In contrast, utilizing SCADA measurements rich in state variable
connections, like power flows and power injections, offers a significant
advantage since these measurements inherently involve the state variables of
multiple buses. Consequently, the state estimator leads to a more constrained
and informative solution, where the estimated state is not only consistent
with the provided measurements but also closer to the true system state. This
directly impacts the performance of the proposed approaches, SF* and PM*, as
evidenced by the numerical results presented in this section.

Figure 2: Measurement residual distributions of power injections at node 19. TABLE I: Active injection residuals statistics for the test set bus index | Mean, $\mu$ | Standard deviation, $\sigma$
---|---|---
UN | PM | PM∗ | UN | PM | PM∗
4 | 0.42 | 0 | 0 | 3.135 | 3.165 | 2.444
7 | -0.11 | 0 | 0 | 0.849 | 1.201 | 0.771
10 | -4.07 | 0 | 0 | 5.357 | 4.598 | 3.715
11 | 3.35 | 0 | 0 | 6.724 | 5.358 | 4.425
14 | -0.15 | 0 | 0 | 1.738 | 1.743 | 1.435
16 | -0.60 | 0 | 0 | 1.921 | 1.648 | 1.583
18 | 0 | 0 | 0 | 2.308 | 2.137 | 2.016
19 | -1.39 | 0 | 0 | 1.798 | 2.653 | 1.252
20 | 0.23 | 0 | 0 | 2.340 | 2.261 | 2.124
22 | -0.38 | 0 | 0 | 2.041 | 1.974 | 1.917
23 | 0.23 | 0 | 0 | 2.930 | 2.618 | 2.366
25 | -0.38 | 0 | 0 | 2.983 | 2.479 | 2.188
26 | -1.41 | 0 | 0 | 4.288 | 5.309 | 4.086
27 | -0.25 | 0 | 0 | 5.110 | 5.120 | 5.009
28 | -0.70 | 0 | 0 | 2.791 | 2.713 | 2.640
30 | -0.58 | 0 | 0 | 3.270 | 2.254 | 2.198
32 | -0.05 | 0 | 0 | 2.590 | 2.605 | 2.262

###  III-C Results and discussion

In order to compare the performance of the nine approaches presented in
Section III-A, we use the Root Mean Squared Error (RMSE) defined as

$\text{RMSE}=\sqrt{\frac{1}{K\cdot
R}\sum_{k=1}^{K}\sum_{r=1}^{R}(y_{k,r}-\hat{y}_{k,r})^{2}}$ (15)

where $K$ is the number of instances in the test dataset, $R$ is either the
number of buses or lines, depending on the magnitude of study; $y_{k,r}$
represents the true value of the magnitude, and $\hat{y}_{k,r}$ is the
estimated value by each approach. Table II collates the RMSE values for the
state variables (voltage magnitudes in p.u. and angles in degrees) and the
power variables (active and reactive power injections and flows in MW and
MVAr). By disregarding the unrealistic results of the benchmark approach, the
minimum error per magnitude is identified in bold.

TABLE II: RMSE values for the test set Method | $V$ | $\theta$ | $P$ | $Q$ | $P_{f}$ | $Q_{f}$
---|---|---|---|---|---|---
BN | 0.004 | 0.003 | 0.11 | 0.11 | 0.24 | 0.24
FL | 0.029 | 8.641 | 76.96 | 114.04 | 41.88 | 63.43
UN | 0.052 | 0.089 | 25.85 | 33.73 | 17.72 | 22.49
SF | 0.021 | 0.086 | 22.03 | 32.84 | 14.53 | 21.56
SF* | 0.018 | 0.085 | 20.06 | 31.77 | 13.40 | 20.96
PM | 0.085 | 0.097 | 21.79 | 32.76 | 15.87 | 22.24
PM* | 0.041 | 0.087 | 19.94 | 31.59 | 13.48 | 20.97
1NN | 0.038 | 0.170 | 63.55 | 99.33 | 38.01 | 59.36
20NN | 0.027 | 0.120 | 44.54 | 68.16 | 26.76 | 41.08

We begin our analysis of the results shown in Table II by considering the
outcomes obtained through the UN method first. Although this approach
estimates the state of an unobservable system, it significantly reduces
estimation errors in power variables and voltage angles, surpassing the
performance of the naive approaches 1-NN and 20-NN. Conversely, the voltage
magnitude errors are increased with respect to these methods. This conclusion
is aligned with the results provided in Fig. 2 and Table I and clearly
indicates that the state yielded by the non-unique solution problem (8)
contains valuable information to be exploited. In fact, the results in Table
II also show the improvement achieved by incorporating the pseudo-measurements
derived from the UN method as additional features in the state estimation
procedures since both SF∗ and PM∗ exhibit the lower estimation errors in all
variables than their respective standard versions, SF and PM. Furthermore, it
is worth noting that SF∗ and PM∗ exhibit the lowest estimation errors for
voltage and power injection variables, respectively. This makes sense as the
regression problems tackled by SF∗ and PM∗ are designed to respectively
minimize the squared error of the voltage and of a measurement set with high
proportion of power injection measurements. This reveals that incorporating a
solution of an unobservable state estimation routine as explanatory features
into state forecasting or pseudo-measurement generation procedures has a
positive effect on the state estimation accuracy.

Figure 3: Comparative of the absolute estimation error for the nine presented
approaches of the voltage magnitude (left) and the apparent power flow
(right).

Lastly, in order to provide a comprehensive view of the estimation errors,
Fig. 3 complements Table II with box plots corresponding to the distributions
of the estimation errors for the nine compared approaches, with ${K\cdot R}$
data points in each case. Specifically, the voltage magnitude and the apparent
power flow have been chosen as the key magnitudes. The voltage magnitude is
crucial to ensure a stable operation and prevent equipment damage, while the
precise estimation of apparent power flow is necessary to ensure that the
power system operates within its thermal and electrical limits, avoiding
overloads and potential failures. For ease of illustration, outliers have been
removed from the plots. The figure reveals that the enhancements introduced in
SF∗ and PM∗ not only reduce the mean error but also decrease its variance,
with respect to SF and PM, respectively. This observation highlights the
effectiveness of the proposed enhancing methodology in achieving more
consistent and accurate state estimations.

The second part of this analysis investigates how the proposed methodology
handles variations in active distribution network operating conditions. These
variations, often influenced by the penetration of renewables, storage, and
prosumers, can significantly impact performance. To evaluate this effect,
state variables are generated using different uniform probability distribution
ranges, as detailed in Table III for low, medium, and high variability
scenarios.

TABLE III: Variability scenarios Variability | Voltage Magnitude ($V$) | Voltage angle ($\theta$)
---|---|---
Low | $\mathcal{U}(0.975,1.025)$ | $\mathcal{U}(-7.5\degree,+7.5\degree)$
Medium | $\mathcal{U}(0.95,1.05)$ | $\mathcal{U}(-15.5\degree,+15.5\degree)$
High | $\mathcal{U}(0.925,1.075)$ | $\mathcal{U}(-22.5\degree,+22.5\degree)$

Beyond the variability of operating conditions, the performance of the
proposed approaches can also be influenced by the type and number of available
measurements, as discussed in Section III-B. To comprehensively evaluate the
robustness of our proposed methodology, we conduct additional experiments
where we vary the types and quantities of measurements used in the estimation
process. In particular, the number of PMU and SCADA measurement devices is
changed.

To quantify the performance of the proposed approaches across different
scenarios, we employ $\Delta$RMSE defined as the difference between the RMSE
of the standard method and the enhanced one, that is,

$\Delta\text{RMSE}=\text{RMSE}-\text{RMSE}^{*}$ (16)

Consequently, positive values of this metric indicate that the enhanced method
exhibits a lower error compared to the standard counterpart. Table IV
summarizes the performance comparison between PM and PM∗. Similar trends and
conclusions were observed for the SF method, which are omitted for brevity.
This table includes results for various variability scenarios (low, medium,
high) and different measurement setups (number of PMU and SCADA devices). The
table shows the difference in Root Mean Squared Error (RMSE) for both voltage
magnitudes ($V$) and apparent power flow ($S_{f}$). Interestingly, all the
values in Table IV are positive, signifying that the proposed method PM∗
consistently outperform PM across all scenarios.

TABLE IV: $\Delta$RMSE for varying variability and measurement data | Variability | # PMU | # SCADA | $\Delta\text{RMSE}\,(V)$ | $\Delta\text{RMSE}\,(Sf)$
---|---|---|---|---|---
PM vs. PM∗ | Low | 5 | 15 | 0.015 | 0.065
Medium | 5 | 15 | 0.044 | 0.729
High | 5 | 15 | 0.070 | 1.957
Medium | 15 | 5 | 0.002 | 0.070
Medium | 5 | 25 | 0.016 | 0.452

To assess the impact of variability on performance, we can examine the first
three rows. It is apparent that the enhanced methods consistently outperform
the baselines, with this advantage increasing as variability rises. This trend
suggests that the additional features introduced by the UN method become more
valuable as operating conditions become less predictable. In simpler terms,
when conditions are stable, the standard methods using linear relationships
perform well. However, under high variability, the additional features of the
proposed methods help them handle the increased complexity of the system,
leading to more accurate state estimation.

To understand how measurement type influences performance improvement, we
compare rows 2 and 4 of Table IV. These rows represent scenarios with the same
total number of measurements (43) but different ratios of voltage phasor (PMU)
and power injection (SCADA) devices. As evident, the scenario with more SCADA
devices (row 2) exhibits a higher improvement compared to the one with more
PMUs (row 4). This aligns with the proposed methods’ dependence on the UN
method. The UN method thrives on rich measurement sets that connect multiple
state variables across different buses, providing more informative data
compared to relying solely on isolated voltage measurements from PMUs.

Finally, we investigate the impact of the number of SCADA measurements,
considering their previously observed advantage over PMU measurements.
Interestingly, our analysis reveals an important trend when comparing rows 2
and 5 of Table IV. While both scenarios utilize SCADA devices, a medium amount
(15, as in row 2) provides greater improvement for both voltage magnitudes and
apparent power flow compared to a high amount (25, as in row 5). This can be
explained by the role SCADA measurements play in the pseudo-measurement
generation process within method PM. A larger set of SCADA injections
significantly improves the estimation accuracy of method PM itself, reducing
the need for the additional features introduced by the UN method.

##  IV Conclusions

Our paper presents a novel approach to address the challenge of state
estimation in distribution systems. This method leverages the outputs of an
unobservable state estimator to exploit information about the joint
probability distribution of real-time and delayed measurements. It enriches
learning-based methodologies by incorporating these outputs as additional
features, enabling either direct state variable forecasting or real-time
pseudo-measurement generation. Numerical simulations on a realistic
distribution network demonstrate the superiority of our proposed method across
various scenarios. In particular, it demonstrates a clear advantage in active
distribution networks with highly variable operating conditions. Additionally,
the performance improves when the unobservable state estimator utilizes
measurements encompassing a high number of state variables.

##  V Acknowledgments

The authors thank the Supercomputing and Bioinnovation Center (SCBI) of the
University of Málaga for their provision of computational resources (the
supercomputer Picasso) and technical support (www.scbi.uma.es/site).

## References

  * [1] A. Gómez-Expósito, C. Gómez-Quiles, and I. Džafić, “State Estimation in Two Time Scales for Smart Distribution Systems,” _IEEE Transactions on Smart Grid_ , vol. 6, no. 1, pp. 421–430, Jan. 2015.
  * [2] A. Abur and A. Gómez-Expósito, _Power system state estimation: theory and implementation_ , ser. Power engineering. New York, NY: Marcel Dekker, 2004.
  * [3] S. Gao, Z. He, X. Wei, Y. Liu, T. Huang, and J. Lei, “Bilevel Model for Protection-Branch Measurements-Based Topology Attack Against DC and AC State Estimations,” _IEEE Systems Journal_ , vol. 16, no. 4, pp. 5369–5379, 2022.
  * [4] Y. Wang, C. Chen, J. Wang, and R. Baldick, “Research on Resilience of Power Systems Under Natural Disasters–A Review,” _IEEE Transactions on Power Systems_ , vol. 31, no. 2, pp. 1604–1613, Mar. 2016.
  * [5] G. Cheng, Y. Lin, A. Abur, A. Gómez-Expósito, and W. Wu, “A Survey of Power System State Estimation using Multiple Data Sources: PMUS, SCADA, AMI, and Beyond,” _IEEE Transactions on Smart Grid_ , 2023\.
  * [6] T. A. Papadopoulos, E. O. Kontis, G. A. Barzegkar-Ntovom, and P. N. Papadopoulos, “A Three-Level Distributed Architecture for the Real-Time Monitoring of Modern Power Systems,” _IEEE Access_ , vol. 10, pp. 29 287–29 306, 2022.
  * [7] Y. Liu, L. Wu, and J. Li, “D-PMU based Applications for Emerging Active Distribution Systems: A Review,” _Electric Power Systems Research_ , vol. 179, p. 106063, 2020.
  * [8] R. R. Mohassel, A. Fung, F. Mohammadi, and K. Raahemifar, “A Survey on Advanced Metering Infrastructure,” _International Journal of Electrical Power & Energy Systems_, vol. 63, pp. 473–484, 2014.
  * [9] Y. Wang, Q. Chen, T. Hong, and C. Kang, “Review of Smart Meter Data Analytics: Applications, Methodologies, and Challenges,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 3, pp. 3125–3148, 2018.
  * [10] K. Dehghanpour, Z. Wang, J. Wang, Y. Yuan, and F. Bu, “A Survey on State Estimation Techniques and Challenges in Smart Distribution Systems,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 2, pp. 2312–2322, 2019.
  * [11] R. Singh, B. C. Pal, and R. A. Jabr, “Distribution System State Estimation through Gaussian Mixture Model of the Load as Pseudo-Measurement,” _IET Generation, Transmission & Distribution_, vol. 4, no. 1, pp. 50–59, Jan. 2010.
  * [12] J. Zhao, G. Zhang, Z. Y. Dong, and M. La Scala, “Robust Forecasting Aided Power System State Estimation Considering State Correlations,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 4, pp. 2658–2666, Jul. 2018.
  * [13] E. Manitsas, R. Singh, B. C. Pal, and G. Strbac, “Distribution System State Estimation Using an Artificial Neural Network Approach for Pseudo Measurement Modeling,” _IEEE Transactions on Power Systems_ , vol. 27, no. 4, pp. 1888–1896, Nov. 2012.
  * [14] B. Huang and J. Wang, “Applications of Physics-Informed Neural Networks in Power Systems–A Review,” _IEEE Transactions on Power Systems_ , vol. 38, no. 1, pp. 572–588, Jan. 2023.
  * [15] M. Kamal, W. Li, D. Deka, and H. Mohsenian-Rad, “Physics-Conditioned Generative Adversarial Networks for State Estimation in Active Power Distribution Systems with Low Observability,” in _2022 International Conference on Smart Grid Synchronized Measurements and Analytics (SGSMA)_ , May 2022, pp. 1–6.
  * [16] L. Wang, Q. Zhou, and S. Jin, “Physics-Guided Deep Learning for Power System State Estimation,” _Journal of Modern Power Systems and Clean Energy_ , vol. 8, no. 4, pp. 607–615, 2020.
  * [17] D. M. Vinod Kumar and S. C. Srivastava, “Power System State Forecasting Using Artificial Neural Networks,” _Electric Machines & Power Systems_, vol. 27, no. 6, pp. 653–664, May 1999.
  * [18] K. R. Mestav, J. Luengo-Rozas, and L. Tong, “Bayesian State Estimation for Unobservable Distribution Systems via Deep Learning,” _IEEE Transactions on Power Systems_ , vol. 34, no. 6, pp. 4910–4920, Nov. 2019.
  * [19] A. S. Zamzam, X. Fu, and N. D. Sidiropoulos, “Data-Driven Learning-Based Optimization for Distribution System State Estimation,” _IEEE Transactions on Power Systems_ , vol. 34, no. 6, pp. 4796–4805, Nov. 2019.
  * [20] L. Pagnier and M. Chertkov, “Physics-Informed Graphical Neural Network for Parameter & State Estimations in Power Systems,” Feb. 2021.
  * [21] L. Zhang, G. Wang, and G. B. Giannakis, “Real-Time Power System State Estimation and Forecasting via Deep Unrolled Neural Networks,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 15, pp. 4069–4077, Aug. 2019.
  * [22] L. Dabush, A. Kroizer, and T. Routtenberg, “State Estimation in Partially Observable Power Systems via Graph Signal Processing Tools,” _Sensors_ , vol. 23, no. 3, p. 1387, 2023.
  * [23] G. Wang, G. B. Giannakis, and J. Chen, “Robust and Scalable Power System State Estimation via Composite Optimization,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 6, pp. 6137–6147, Nov. 2019.
  * [24] E. Caro and A. Hernández, “State Estimation in Power Systems Based on a Mathematical Programming Approach,” in _Advances in Electric Power and Energy_. John Wiley & Sons, Ltd, 2020, pp. 23–57.
  * [25] A. Primadianto and C. N. Lu, “A Review on Distribution System State Estimation,” _IEEE Transactions on Power Systems_ , vol. 32, no. 5, pp. 3875–3883, 2016.
  * [26] P. Li, W. Wu, X. Wang, and B. Xu, “A Data-Driven Linear Optimal Power Flow Model for Distribution Networks,” _IEEE Transactions on Power Systems_ , vol. 38, no. 1, pp. 956–959, 2023.
  * [27] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: Steady-State operations, Planning, and Analysis Tools for Power Systems Research and Education,” _IEEE Transactions on power systems_ , vol. 26, no. 1, pp. 12–19, 2010.
  * [28] R. Hidalgo, C. Abbey, and G. Joós, “A Review of Active Distribution Networks Enabling Technologies,” in _IEEE PES General Meeting_ , 2010, pp. 1–9.
  * [29] J. A. D. Massignan, J. B. A. London, M. Bessani, C. D. Maciel, R. Z. Fannucchi, and V. Miranda, “Bayesian Inference Approach for Information Fusion in Distribution System State Estimation,” _IEEE Transactions on Smart Grid_ , vol. 13, no. 1, pp. 526–540, 2022.
  * [30] W. E. Hart, J. P. Watson, and D. L. Woodruff, “Pyomo: Modeling and Solving Mathematical Programs in Python,” _Mathematical Programming Computation_ , vol. 3, no. 3, pp. 219–260, 2011.
  * [31] A. Wächter and L. T. Biegler, “On the Implementation of an Interior-Point Filter Line-Search Algorithm for Large-Scale Nonlinear Programming,” _Mathematical programming_ , vol. 106, pp. 25–57, 2006.
  * [32] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine Learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.


# Asymptotic mass limit of large fully-heavy compact multiquarks

M. C. Gordillo<EMAIL_ADDRESS>Departamento de Sistemas Físicos, Químicos y
Naturales, Universidad Pablo de Olavide, Carretera de Utrera km 1, E-41013
Sevilla, Spain  Instituto Carlos I de Física Teórica y Computacional,
Universidad de Granada, E-18071 Granada, Spain.  J. M. Alcaraz-Pelegrina
<EMAIL_ADDRESS>Departamento de Física. Universidad de Córdoba, Campus de
Rabanales, Edif. C2 E-14071 Córdoba, Spain

###### Abstract

The properties of fully-heavy arrangements including a number of quarks
between 5 and 12 were calculated within the framework of a constituent quark
model by using a diffusion Monte Carlo technique. We considered only clusters
in which all the quarks had the same mass, and whose number of particles and
antiparticles were adequate to produce color singlets. All the multiquarks
were in their lowest possible values of $L^{2}$ and $S^{2}$ operators. This
means that we considered only color-spin wavefunctions that were antisymmetric
with respect to the interchange of any two quarks of the same type. We found
that in both all-$c$ and all-$b$ multiquarks, the mass per particle levels off
for arrangements with a number of quarks larger of equal than six. The
analysis of their structure implies that the fully-heavy multiquarks are
compact structures.

Protons and neutrons are the basic constituents of atomic nuclei. Quantum
chromodynamics (QCD) is the theory that describe them as a composite set of
quarks and gluons interacting through the strong force. However, QCD is not
limited to associations to the light quarks ($u,d$) that make up the nucleons,
but extends to other types of particles collectively called hadrons. Those
hadrons can include or being totally made of heavier quarks ($s,c,b$).

Unfortunately, as of today, it is impossible to solve analytically the QCD
equations and deduce the hadron spectrum. Among the phenomenological QCD-
inspired models designed to fill that gap the so-called quark model stands
out. It considers only the valence quarks and antiquarks within the hadrons
and was independently proposed by Murray Gell-Mann [1] and George Zweig [2].
Even though it was designed to account for the properties of mesons (one quark
and one antiquark) and baryons (three quarks), it also opens the door to
larger associations of quarks such as tetra and pentaquarks [3, 4]. We can
even have hexaquarks, such as the experimentally produced deuteron [5] and the
well-established $d^{*}(2380)$ resonance [6, 7, 8, 9, 10]. The quark model
does not impose, in principle, any limit to the upper size of those clusters
of quarks, and in this work we will use it to obtain the masses of all
possible fully-heavy multiquarks. Of all the possible compositions of those
clusters, we will stick to arrangements in which all the quarks have the same
mass. This means to consider (see below) sets up to 12 $c$ or $b$ quarks
and/or antiquarks. To do so, we have to solve the Schrödinger equation derived
from the Hamiltonian [11]:

$H=\sum_{i=1}^{N_{q}}\left(m_{i}+\frac{\vec{p}^{\,2}_{i}}{2m_{i}}\right)+\sum_{i<j}^{N_{q}}V(r_{ij})\,,$
(1)

where $N_{q}$ is the number of quarks, while $m_{i}$ and $\vec{p}_{i}$ are the
mass and momentum of the $i$ quark. This a non-relativistic approximation, and
it is expected to work best for the fully-heavy ensembles that will be
considering in this work. To produce experimentally those multiquarks is, in
principle, possible, as the discovery of the X(6900) (though to be a fully
$c$-tetraquark) attests [12]. $V(r_{ij})$, is a two-body potential that
depends only on the distance between quarks, $r_{ij}$, and can be written as
the sum of one-gluon exchange term given by [13, 14] :

$V_{\text{OGE}}(r_{ij})=\frac{1}{4}\alpha_{s}(\vec{\lambda}_{i}\cdot\vec{\lambda}_{j})[\frac{1}{r_{ij}}-\frac{2\pi}{3m_{i}m_{j}}\delta^{(3)}(r_{ij})(\vec{\sigma}_{i}\cdot\vec{\sigma}_{j})]\,,$
(2)

that includes both Coulomb and hyperfine terms, and the lineal confining
potential:

$V_{\text{CON}}(\vec{r}_{ij})=(b\,r_{ij}+\Delta)(\vec{\lambda}_{i}\cdot\vec{\lambda}_{j}).$
(3)

that approximates the contribution of multigluon exchanges. $\vec{\lambda}$
and $\vec{\sigma}$ are the Gell-Mann and Pauli matrices, respectively, and
account for the color and spin degrees of freedom. The Dirac delta function
was regularized in the standard way [15, 16, 17] in order to make possible the
calculations. The parameters needed to fully define the interaction were taken
from Refs. [15, 16], and were the same as the used in previous calculations
for smaller clusters [17, 18, 19, 20]. The masses of the hadrons computed with
this potential were found to be in good agreement with experimental data, when
available [17]. Since this non-relativistic approximation applies best to
heavy quarks, in order to describe light quarks (u, d or s) we would have to
include additional terms [21], something that will not be done in this work.

To solve the Schrödinger equation derived from the Hamiltonian in Eq. 1, we
resorted to a diffusion Monte Carlo (DMC) scheme [22, 23, 24, 25, 17]. This
will provide us with the desired masses of the ground states of the different
set of quarks. This method needs an initial approximation to the real many-
body wavefunction of the clusters, the trial function, that should include all
the information known a priori about the different systems. We chose the
expression [17]:

$\displaystyle\Psi({\bf r_{1},r_{2}},\ldots,{\bf
r_{n}},s_{1},s_{2},\ldots,s_{n},c_{1},c_{2},\ldots,c_{n})=$
$\displaystyle\Phi({\bf r_{1},r_{2}},\ldots,{\bf r_{n}})$
$\displaystyle\left[\chi_{s}(s_{1},s_{2},\ldots,s_{n})\bigotimes\chi_{c}(c_{1},c_{2},\ldots,c_{n})\right],$
(4)

where ${\bf r_{i}}$, $s_{i}$ and $c_{i}$ stand for the position, spin and
color of the particle $i$, that is inside a cluster of $n$ quarks. In this
work, we are going to consider only multiquark states that are eigenvectors of
the angular momentum operator, $L^{2}$ with eigenvalue $\ell$ = 0. This means
that $\Phi$ should depend on the distance between pairs of quarks and not on
their absolute positions. Following Ref. 17, we have used:

$\Phi({\bf r_{1},r_{2}},\ldots,{\bf
r_{n}})=\prod_{i=1}^{N_{q}}\exp(-a_{ij}r_{ij}),$ (5)

No other alternatives to the form of the radial part of the trial function
were considered in this work since, in principle, the DMC algorithm should be
able to correct its possible shortcomings and produce the exact masses of the
arrangements [24]. The $a_{ij}$ values were chosen in accordance to the
boundary conditions of the problem [17]. $\chi_{s}$ and $\chi_{c}$ are linear
combinations of the eigenvectos of the spin and color operators defined by:

$F^{2}=\left(\sum_{i=1}^{N_{q}}\frac{\lambda_{i}}{2}\right)^{2}$ (6)

and

$S^{2}=\left(\sum_{i=1}^{N_{q}}\frac{\sigma_{i}}{2}\right)^{2}.$ (7)

with eigenvalues $F^{2}$ = 0 (colorless functions) and $S$ = 0 or 1/2,
depending on whether the number of quarks in the multiquark is ever or odd,
respectively. Those are the lowest possible eigenvalues for the spin operator
and the only ones considered in this work. For instance, for the
$cccc\bar{c}\bar{c}\bar{c}\bar{c}$ ($c^{4}\bar{c}^{4}$) octaquark, we have 23
color and 14 spin functions meeting those criteria. This means 322
$\chi_{s}\bigotimes\chi_{c}$ possible combinations. That said, we have to
remember that since Eq. 5 is symmetric with respect to the exchange of any two
identical quarks, we have to produce spin-color combinations antisymmetric
with respect to those exchanges, as befits to a set of fermions as quarks are.
To do so, we apply the antisymmetry operator

$\mathcal{A}=\frac{1}{N}\sum_{{\alpha}=1}^{N}(-1)^{P}\mathcal{P_{\alpha}}$ (8)

to that color-spin set of functions. Here, $N$ is the number of possible
permutations of the set of quark indexes, $P$ is the order of the permutation,
and $\mathcal{P_{\alpha}}$ represents the matrices that define those
permutations. Once constructed the matrix derived from the operator in Eq. 8,
we have to check if we can find any eigenvector with eigenvalue equal to one.
If this is so, those combinations will be the input of the DMC calculation
[17]. For the octaquark, we have that of all the 322 color-spin functions,
only 2 are antisymmetric with respect to the interchange of all the pairs of
quarks and, separately, of all the pairs of antiquarks. The analysis of the
eigenvectors of the antisymmetry operator indicates that there are no
antisymmetric color-spin functions for structures in which any of the quark or
antiquark subsets contains more than 6 units. This means that the largest
possible fully heavy multiquark is the $c^{6}\bar{c}^{6}$ dodecaquark.
Moreover, neither the $c^{9}$ nonaquark nor the $c^{7}\bar{c}^{4}$
undecaquark, or their $b$-counterparts are viable structures. Independently,
the $c^{6}\bar{c}^{3}$ nonaquark is also impossible since no antisymmetric
color-spin combinations with respect to the interchange of any pair of $c$
quarks were found.

Table 1:  Masses of the mutiquarks considered in this work in MeV. The error bars are given in parenthesis.  | $c^{4}\bar{c}$ | $c^{3}\bar{c}^{3}$ | $c^{5}\bar{c}^{2}$
---|---|---|---
Mass | 8195(2) | 9614(2) | 11543(4)
| $c^{4}\bar{c}^{4}$ | $c^{5}\bar{c}^{5}$ | $c^{6}\bar{c}^{6}$
Mass | 13133(4) | 16539(4) | 19808(4)
| $b^{4}\bar{b}$ | $b^{3}\bar{b}^{3}$ | $b^{5}\bar{b}^{2}$
Mass | 24211(2) | 28822(2) | 33970(4)
| $b^{4}\bar{b}^{4}$ | $b^{5}\bar{b}^{5}$ | $b^{6}\bar{b}^{6}$
Mass | 38815(4) | 48599(4) | 58232(4)

The masses of the multiquarks obtained by the DMC method are given in Table 1.
As indicated above, all are colorless clusters with $S$=0 or 1/2 depending on
whether the total number of quarks is even or odd, respectively. We have to
stress that the color-spin functions used in the calculations are the
eigenvalues of the antisymmetry operator given in Eq. 8, with no quark
groupings different to those that put together identical particles. For
instance, in pentaquarks, we do not consider baryon+meson or
diquark+diquark+antiquark arrangements [26, 27], but a function that is
antisymmetric with respect to the exchange of any pair of the four quarks
considered to be undistinguishable. In any case, the results for that
particular multiquark are virtually identical to those of Ref. 28, in which
the same function is used. Those results validate our approach, that allows us
to dispense with Young-tableaux diagrams to calculate larger clusters.

Figure 1: Mass per particle in all the $c$-multiquarks considered in this
work. Where not shown, error bars are of the size of the symbols. The data for
the meson, baryon and tetraquark are taken from Ref. 17, while the upper
symbol for the hexaquark corresponds to the mass of the open charm hexaquark
given in Ref. 20. The dashed line represents the average mass for the open
charm hexaquark , and the hepta-, octa-, nona-, deca- and docedaquarks.

To better visualize the results in Table 1, we display the mass per particle
as a function of the number of particles in the cluster in Figs. 1 and 2. The
data not given in Table 1 are taken from Refs. 17 and 18. Something is
immediately apparent: from the open-charm hexaquark up, the mass per particle
of the clusters reaches a plateau both for $c$\- and $b$-multiquarks. This
basically means that to modify the number of quarks beyond six, we will have
to increase the mass of the system by a constant value of 1649 and 4854 MeV
per particle for $c$\- and $b$-multiquarks, respectively.

Figure 2: Same as in the previous figure, but for $b$-multiquarks. Error bars
are of the size of the symbols and not shown for simplicity. The dashed line
have the same meaning as in Fig. 1 but for the $b$-multiquarks. The source of
the data is the same as for their $c$-counterparts.  Figure 3: Radial
distribution functions for different particle pairs in $c$-multiquarks. Solid
lines, $c-c$ pairs (equivalent to $\bar{c}-\bar{c}$); symbols, $c-\bar{c}$
pairs. The hexaquak displayed correspond to the $c^{3}\bar{c}^{3}$ hidden
charm structure. All distributions are normalized to one.

The structure of the clusters can be deduced from the radial distribution
functions, depicted in Figs. 3 and 4. Those give us the probabilities of
having another particle at a particular distance of a given one. We show only
the more representative structures, the remaining ones being similar to those
displayed. First, we can see that all the clusters are compact structures,
i.e., the probability of finding another particle at distances beyond a
maximum of 2 fm goes rapidly to zero. In addition, in the majority of cases
there is very little difference between the probability of finding another
quark (solid lines) or an antiquark (symbols) for any particule at a given
distance. This is similar to what happens for smaller multiquarks [17, 18].
The only exception is the hidden-charm hexaquark, in which the $c-c$ and
$c-\bar{c}$ are noticiably different, and in which the first of them is
virtually identical to the corresponding to the $ccc$ baryon. The reason is
that in that system, the quarks and antiquarks group to produce a baryon and
an antibaryon glued together. The same happens with the $b^{3}\bar{b}^{3}$
system.

Figure 4: Comparison between the pair distribution functions for a
$c$-decaquark and a $b$-decaquark. Lines, $q-q$ pairs; symbols, $q-\bar{q}$
pairs.

In this work we have calculated the color-spin functions with an algorithm
that dispense with the need to use of Clebsch-Gordan coefficients. This is
necessary since the increase in the number of color-spin functions with the
number of quarks makes that approximation impossible. For instance, for an
heptaquark, we have 11 color and 14 spin functions that make a total of 154
combinations. This is to be compared with the 15 color-spin possibilities for
a pentaquark [28] or the 25 for an open-charm hexaquark [29, 30]. The use of
this technique in combination with a DMC algorithm, originally developed to
deal with many-body systems, allowed us obtain the masses of all possible
fully-heavy s-wave multiquarks. What we found is that, from a number of quarks
beyond six, the mass of those systems is linearly proportional to the number
of particles in the arrangements, i.e., in relative terms, there is no mass
penalty in producing progressively larger multiquarks, as it is in going from
a meson to a tetraquark. This means that, in mass terms, is equally probable
to have an open-charm hexaquark as to produce an heptaquark or octoquark.

###### Acknowledgements.

We acknowledge financial support from Ministerio de Ciencia e Innovación
MCIN/AEI/10.13039/501100011033 (Spain) under Grant No. PID2020-113565GB-C22
and from Junta de Andalucía group PAIDI-205. We also acknowledge the use of
the C3UPO computer facilities at the Universidad Pablo de Olavide.

## References

  * Gell-Mann [1964] M. Gell-Mann, A Schematic Model of Baryons and Mesons, Phys. Lett. 8, 214 (1964).
  * Zweig [1964] G. Zweig, Developments in the quark theory of hadrons, CERN-TH-412, NP-14146 (1964).
  * H.X.Chen  _et al._ [2016] H.X.Chen, W. Chen and X. Liu, and S.-L. Zhu, The hidden-charm pentaquark and tetraquark states, Phys. Rep. 639, 1 (2016).
  * H.X.Chen  _et al._ [2017] H.X.Chen, W. Chen, X. Liu and Y.R. Liu, and S.-L. Zhu, A review of the open charm and open bottom systems, Rep. Prog. Phys. 80, 076201 (2017).
  * Urey  _et al._ [1932] H. C. Urey, F. G. Brickwedde, and G. M. Murphy, A Hydrogen Isotope of Mass 2, Phys. Rev. 39, 164 (1932).
  * et al. [CELSIUS-WASA Collaboration(2009] M. B. et al. (CELSIUS-WASA Collaboration, Double-Pionic Fusion of Nuclear Systems and the “ABC” Effect:Approaching a Puzzle by Exclusive and Kinematically Complete Measurements, Phys. Rev. Lett. 102, 052301 (2009).
  * et al. [WASA-at COSY Collaboration(2011] P. A. et al. (WASA-at COSY Collaboration, Abashian-Booth-Crowe Effect in Basic Double-Pionic Fusion: A New Resonance?, Phys. Rev. Lett. 106, 242302 (2011).
  * et al. [WASA-at COSY Collaboration(2013a] P. A. et al. (WASA-at COSY Collaboration, Isospin decomposition of the basic double-pionic fusion in the region of the ABC effect, Phys. Lett. B 721, 229 (2013a).
  * et al. [WASA-at COSY Collaboration(2013b] P. A. et al. (WASA-at COSY Collaboration, Measurement of the pn $\rightarrow$ pp $\pi^{0}\pi^{-}$ reaction in isearch for the recently observed resonance structure in $d\pi^{0}\pi^{0}$ and $d\pi^{+}\pi^{-}$ systems, Phys. Rev. C 88, 055208 (2013b).
  * P. Adlarson et al. [WASA-at COSY Collaboration(2014] S. D. A. C. P. Adlarson et al. (WASA-at COSY Collaboration, Evidence for a New Resonance from Polarized Neutron-Proton Scattering, Phys. Rev. Lett. 112, 202301 (2014).
  * Valcarce  _et al._ [2005] A. Valcarce, H. Garcilazo, F Fernández, and P. González, Quark model study of few-baryon systems, Rep. Prog. Phys. 68, 965 (2005).
  * et al [LHCb collaboration] R. A. et al (LHCb collaboration), Observation of structure in the J/$\Psi$-pair mass spectrum, Sci. Bull. 65, 1983 (2020).
  * De Rújula  _et al._ [1975] A. De Rújula, H. Georgi, and S. L. Glashow, Hadron masses in a gauge theory, Phys. Rev. D 12, 147 (1975).
  * Bhaduri  _et al._ [1981] R. K. Bhaduri, L. E. Cohler, and Y. Nogami, A Unified Potential for Mesons and Baryons, Nuovo Cim. A 65, 376 (1981).
  * Semay and Silvestre-Brac [1994] C. Semay and B. Silvestre-Brac, Diquonia and potential models, Z. Phys. C 61, 271 (1994).
  * Silvestre-Brac [1996] B. Silvestre-Brac, Spectrum and static properties of heavy baryons, Few Body Syst. 20, 1 (1996).
  * Gordillo  _et al._ [2020] M. C. Gordillo, F De Soto, and J. Segovia, Diffusion Monte Carlo calculations of fully-heavy multiquark bound states, Phys. Rev. D 102, 111007 (2020).
  * Gordillo  _et al._ [2021] M. C. Gordillo, F. De Soto, and J. Segovia, Structure of the X(3872) as explained by a diffusion Monte Carlo calculation, Phys. Rev. D 104, 054036 (2021).
  * Gordillo  _et al._ [2022] M. C. Gordillo, F. De Soto, and J. Segovia, $x(3872)$’s excitation and its connection with production at hadron colliders, Phys. Rev. D 106, 094004 (2022).
  * Alcaraz-Pelegrina and Gordillo [2022] J. M. Alcaraz-Pelegrina and M. C. Gordillo, Diffusion monte carlo calculations of fully heavy compact hexaquarks, Phys. Rev. D 106, 114028 (2022).
  * Vijande  _et al._ [2005] J. Vijande, F. Fernández, and A. Valcarce, Constituent quark model study of the meson spectra, Journal of Physics G: Nuclear and Particle Physics 31, 481 (2005).
  * Kalos [1970] M. H. Kalos, Energy of a Boson Fluid with Lennard-Jones Potentials, Phys. Rev. A 2, 250 (1970).
  * Boronat and Casulleras [1994] J. Boronat and J. Casulleras, Monte Carlo analysis of an interatomic potential for He, Phys. Rev. B 49, 8920 (1994).
  * Hammond  _et al._ [1994] B. Hammond, W. Lester, and P. Reynolds,  _Monte Carlo Methods in ab Initio Quantum Chemistry_ (World Scientific, Singapore, 1994).
  * Sánchez-Baena  _et al._ [2018] J. Sánchez-Baena, J. Boronat, and F. Mazzanti, Diffusion Monte Carlo methods for spin-orbit-coupled ultracold Bose gases, Phys. Rev. A 98, 053632 (2018).
  * Yang  _et al._ [2022] G. Yang, J. Ping, and J. Segovia, Fully charm and bottom pentaquarks in a lattice-qcd inspired quark model, Phys. Rev. D 106, 014005 (2022).
  * Yan  _et al._ [2022] Y. Yan, Y. Wu, X. Hu, H. Huang, and J. Ping, Fully heavy pentaquarks in quark models, Phys. Rev. D 105, 014027 (2022).
  * An  _et al._ [2022] H.-T. An, S.-Q. Luo, Z.-W. Liu, and X. Liu, Fully heavy pentaquark states in constituent quark model, Phys. Rev. D 105, 074032 (2022).
  * Kima  _et al._ [2020] H. Kima, K.S. Kim, and M. Oka, Hexaquark picture for d∗(2380), Phys. Rev. D 102, 074023 (2020).
  * Lü  _et al._ [2022] Q. Lü, D.Y. Chen, and Y. Dong, Fulliy-heavy hexaquarks in a constituent quark model, (2022), arXiv:hep-ph/2208.03041 .