# Reducing Learning Difficulties: One-Step Two-Critic Deep Reinforcement
Learning for Inverter-based Volt-Var Control in Active Distribution Networks

Qiong Liu, Ye Guo, Lirong Deng, Haotian Liu, Dongyu Li, Hongbin Sun, and Wenqi
Huang This work was supported in part by the National Key R&D Program of China
(2020YFB0906000, 2020YFB0906005).Qiong Liu, Ye Guo are with the Tsinghua-
Berkeley Shenzhen Institute, Tsinghua University, Shenzhen, 518071, Guangdong,
China, e-mail<EMAIL_ADDRESS>Deng is with the Department of
Electrical Engineering, Shanghai University of Electric Power, Shanghai,
200000, ChinaHaotian Liu, Hongbin Sun are with the State Key Laboratory of
Power Systems, Department of Electrical Engineering, Tsinghua University,
Beijing 100084, ChinaDonyu Li is with the School of Cyber Science and
Technology, Beihang University, Beijing, 100191, China.Wenqi Huang is with the
Digital Grid Research Institute, China Southern Power Grid, Guangzhou,
China.Manuscript received XX, 2022; revised XX, 2022.

###### Abstract

A one-step two-critic deep reinforcement learning (OSTC-DRL) approach for
inverter-based volt-var control (IB-VVC) in active distribution networks is
proposed in this paper. First, the problem of IB-VVC is formulated as a one-
step Markov decision process, which reduces the difficulties of the DRL
learning task. Correspondingly, a one-step actor-critic DRL scheme is
designed, which has a simpler structure and avoids the problem of Q-value
over-estimation. Second, considering two objectives of VVC: minimizing power
loss and eliminating voltage violations, we utilize two critics to approximate
the rewards of two objectives separately, which reduces the difficulties of
the approximation tasks of each critic. OSTC-DRL under the simper structure
improves the approximation accuracy of critics, accelerates the convergence
process, and improves the control performance. The OSTC-DRL approach
cooperates well with any actor-critic DRL algorithms for the centralized IB-
VVC problem, and two centralized DRL algorithms were taken as examples. The
multi-agent OSTC-DRL approach is also developed and applied to the
decentralized IB-VVC problem. Extensive simulation experiments show that the
two proposed OSTC-DRL algorithms require fewer iteration times and return
better results than the recent DRL algorithms, and the multi-agent OSTC-DRL
algorithms work well for decentralized IB-VVC problems.

###### Index Terms:

Volt-Var control, deep reinforcement learning, actor-critic, active
distribution network.

##  I Introduction

On our efforts toward a carbon-neutral society, more and more distributed
energy resources will be integrated into active distribution networks (ADNs).
The increasing penetration of distributed generations (DGs) poses new
challenges to voltage regulation. To tackle this problem, Volt-Var control
(VVC) will be more and more important for active distribution networks. VVC
optimizes the output of the reactive power resources to eliminate voltage
violations and minimize power loss. Most DGs are inverter-based energy
resources (IB-ERs) that are possible to provide reactive power support
rapidly. There is an increasing interest to utilize these resources to achieve
the VVC [1, 2].

Recent methods for inverter-based VVC (IB-VVC) problems can be divided into
model-based and data-driven. Model-based methods solve the IB-VVC problems
based on a reliable model of the ADN. However, such a reliable model may be
difficult to acquire for distribution system operators [3]. As an alternative
solution, data-driven methods learn optimal actions from measurements
directly. In data-driven methods, deep reinforcement learning (DRL) methods
are intensively studied. Generally, actor-critic DRL algorithms are usually
implemented to deal with IB-VVC problems [4, 5, 6, 7, 8]. Actor-critic
algorithms train both the state-action value network (also named as critic)
and the policy network (also named as actor). They reduce the variance,
improve data efficiency, and accelerate the learning process [9]. Meanwhile,
DRL algorithms can achieve real-time decision-making by shifting the
computational expense from online optimization to offline training. However,
they may still suffer optimality problems especially in minimizing power loss.
The learning process may be slow, and the convergence tends to be unstable.

From the perspective of the general DRL field, the intrinsic reason for the
above problems is the unavoidable estimation error of critic networks
including overestimation or underestimation [10]. The overestimation error is
accumulated when using the temporal difference learning method to train critic
networks [11]. It may worsen the performance of the policy network or even
lead to divergent behavior. The underestimated bias will not be explicitly
accumulated through the policy update but may still degrade performance [12,
13]. Deep deterministic policy gradient (DDPG) uses a replay buffer and soft
target updates to alleviate the overestimation of critics [14]. Twin delayed
deep deterministic policy gradient (TD3) [11] and soft actor-critic (SAC) [15]
mitigates the overestimation bias and its accumulation by the technology of
clipped double-Q learning. Nevertheless, clipped double-Q learning leads to an
underestimation bias, which also degrades performance [12, 13]. Applying TD3
or SAC algorithms to VVC problems directly may encounter similar problems.

From the property of VVC problems, two objectives of VVC: minimizing power
loss and eliminating voltage violations, brings additional difficulties to the
learning process of DRL algorithms. Generally, the reward of DRL is designed
as the combination of power loss and voltage violation rate [16, 7]. The
weight ratio between power loss and voltage violations needs to be tuned
carefully, otherwise, the performance of DRL would degrade [17, 18]. To
alleviate the problem, a Lagrangian relaxation method is introduced to tune
the weight ratio between power loss and voltage violations in the reward
design [17]. To simply the learning task, the reward is the negative voltage
violation rate when voltage violations appear, and the negative of the power
loss when no voltage violation appears [18]. The critic network only needs to
learn the penalty of voltage violations when voltage violations appear. The
above reward design methods still mix the two objectives in one critic, and
degrade the estimation accuracy.

To address the estimation error of the critic network, papers [19, 20] derive
the critic based on the power flow model. It avoids the inaccuracy estimation
of critic value directly. However, the methods need a reliable model of ADNs,
and the error of the model would degrade the VVC performance.

Given the literature review above, we observe that two main problems may
degrade the performance of DRL for IB-VVC in ADNs: 1) The intrinsic problem is
the estimation errors of critic networks including the overestimation and
underestimation errors; 2) The two objectives of VVC bring additional
difficulties to the learning process of critic networks in DRL algorithms. To
address the two problems, we propose a one-step two-critic DRL (OSTC-DRL)
approach for IB-VVC in ADNs. For the first problem, we formulate the IB-VVC as
a one-step Markov decision process (MDP), and design a one-step actor-critic
DRL scheme to solve the one-step MDP. For the second problem, we use two
critic networks to approximate the rewards of power loss and voltage violation
rate separately. Based on the OSTC-DRL approach, we design two kinds of
centralized DRL algorithms which are OSTC with a deterministic policy (OSTC-
DP) derived from DDPG and OSTC soft actor-critic (OSTC-SAC) derived from SAC.
For decentralized VVC, we extend the OSTC-DRL approach to a multi-agent OSTC-
DRL approach and design two multi-agent DRL algorithms which are multi-agent
OSTC-DP and multi-agent OSTC-SAC. Our proposed approach is simple, stable, and
efficient. Compared with the existing DRL-based VVC algorithms, the main
contributions of this paper and the technical advancements are summarized as
follows:

  * 1.

We propose an OSTC-DRL approach for IB-VVC in ADNs that accelerates the
convergence rate and improves control performance by reducing the difficulty
of DRL learning tasks and decreasing the approximation difficulty of each
critic. It is a simple yet effective DRL approach for IB-VVC.

  * 2.

We analyze that general DRL approaches like DDPG or SAC for IB-VVC problems
complicate the problems and would lead to overestimation or underestimation
issues due to the incompatible between the approaches and the problems, and
the proposed OSTC-DRL approach addresses those issues by learning the policy
to maximize the recent reward rather than the infinity discounted accumulated
reward.

  * 3.

We analyze that the two objectives of IB-VVC would increase the learning
difficulties of one critic, and the proposed OSTC-DRL approach addresses the
issue by using two critics to learn the two objectives separately.

  * 4.

We extend the OSTC-DRL approach to the multi-agent OSTC-DRL approach for
decentralized IB-VVC. Multi-agent OSTC-DRL algorithms have a similar
performance to the centralized OSTC algorithms even when each actor executes
based on the local one bus information. The associated codes in this paper
will be shared on Github.

##  II Problem Formulation

Generally, the problem of IB-VVC is formulated as an MDP. MDP is designed for
formulating long horizon decisions or multi-period optimization problems.
However, the IB-VVC problem can be formulated as a single-period optimization
problem. Using MDP to formulate the problem is feasible but it increases the
complexity of optimization tasks.

For a single-period problem, it would be better to formulate the problem as a
one-step MDP, also known as the contextual bandit [21]. The one-step MDP is
defined by a tuple $(\mathcal{S},\mathcal{A},\mathcal{R})$. At each time step,
the agent observes a state $s\in\mathcal{S}$, and selects actions
$a\in\mathcal{A}$ based on its policy $\pi:\mathcal{S}\rightarrow\mathcal{A}$,
receiving a reward $r\in\mathcal{R}$. For the one-step MDP, the next state is
not included in the process.

For the IB-VVC problems, state space, action space, and reward function are
defined for one-step MDP as follows:

  * 1)

State: The state is $s=(P,Q,V,Q_{G})$. $Q_{G}$ is the vector of reactive power
generations of controllable devices, such as IB-ERs and SVCs. Since this paper
focuses on the IB-VVC, for simplicity, we assume the topology of the ADN does
not change in the control process, so the topology information is not added to
the state [16]. Adding $Q_{G}$ is to reflect the working condition of control
variables in ADNs. Otherwise, we cannot obtain load reactive power $Q_{D}$ and
inverter reactive power $Q_{G}$ with only the nodal reactive power $Q$.

  * 2)

Action: The action is implemented after the agent observe the state $s$. The
action is $a=Q_{G}$, where $Q_{G}$ are the reactive power outputs of all IB-
ERs and SVCs, and the range
$\left|Q_{G}\right|\leq\sqrt{S_{G}^{2}-\overline{P_{G}}^{2}}$ for IB-ERs and
$\underline{Q_{G}}\leq Q_{G}\leq\overline{Q_{G}}$ for SVC [22, 23].
$\overline{P_{G}}$ is the upper limit of active power generation of IB-ERs,
and $\overline{Q_{G}}$ and $\underline{Q_{G}}$ the upper and bottom limit of
reactive power generation of SVCs.

  * 3)

Reward: We assume the new state $s^{\prime}$ can be obtained after executing
action $a$ immediately. Reward is calculated based on $s^{\prime}$. VVC
problems have two objectives: minimizing active power loss, and eliminating
voltage violations. Hence, the reward consists of two terms: the negative of
active power loss $r_{p}$, and the negative of voltage violation rate $r_{v}$.
The reward of power loss $r_{p}$ is defined as:

$r_{p}=-\sum P^{\prime}.$ (1)

Similar to [18, 20], the reward of voltage violation $r_{v}$ is defined as

$r_{v}=-\sum\left[\max\left(V^{\prime}-\bar{V},0\right)+\max\left(\underline{V}-V^{\prime},0\right)\right].$
(2)

##  III One-Step Two-Critic DRL

This section introduces details of the OSTC-DRL approach. We first propose a
one-step DRL scheme to learn from the one-step MDP data, and then utilize the
two-critic to learn the reward of power loss and voltage violations
separately. Fig. 1 shows the framework of the one-step two-critic DRL
approach.

###  III-A One-Step Actor-Critic DRL

Generally, actor-critic DRL algorithms, such as DDPG, TD3, and SAC are used to
solve the IB-VVC problem. However, it complicates the problem and leads the
overestimation or underestimation issues because of the following reasons:

  * 1.

Actor-critic DRL algorithms usually solve long horizontal decision problems.
However, IB-VVC is a single-period optimization problem. Directly applying
those DRL algorithms in IB-VVC problems complicates the problem. It may need
more time and more data in the learning process.

  * 2.

The overestimation issue occurs because actor-critic DRL algorithms like DDPG
use temporal difference learning to update the critic network. The temporal
difference learning is

$y_{t}=r_{t}+\gamma Q_{\phi}(s_{t+1},a_{t+1}),\
a_{t+1}\sim\pi_{\theta}(s_{t+1}).$ (3)

The overestimation error occur in the term $Q_{\phi}(s_{t+1},a_{t+1})$ and the
overestimation error would accumulate in the temporal difference learning
[11].

  * 3.

To address the overestimation issue, clipped double-Q learning is proposed in
TD3 or SAC. Clipped double-Q learning is

$y_{t}=r_{t}+\gamma\min_{i=1,2}Q_{\phi_{i}}(s_{t+1},a_{t+1}),\
a_{t+1}\sim\pi_{\theta}(s_{t+1}).$ (4)

The term $\min_{i=1,2}Q_{\phi_{i}}(s_{t+1},a_{t+1})$ address the
overestimation issue but leads to the underestimation issue. The
underestimation issue would not accumulate in the learning process, but still
degrade the learning performance [12, 13].

To address the problems above, we design the one-step DRL approach where the
objective of one-step DRL is to learn a policy $\pi$ that maximizes the
expected reward,

$\pi^{*}=\arg\max_{\pi}E_{a\sim\pi}(r(s,a)),$ (5)

where $r=r_{p}+c_{v}r_{v}$, and $c_{v}$ is a constant. The learning task is
simpler compared with the ordinary DRL algorithm because the ordinary DRL
algorithms maximizes the expected infinite-horizon discounted accumulated
reward $E_{a\sim\pi}(\sum_{t=0}^{\infty}\gamma r_{t}(s,a))$, where $\gamma<1$
is a discount factor. Meanwhile, the term $Q_{\phi_{i}}(s_{t+1},a_{t+1})$ is
not in the process of learning critic networks, thus avoiding the
overestimation or underestimation issue successfully.

Actor-critic DRL algorithms use the state-action value function $Q(s,a)$
(known as the critic) to guide the update of the policy $\pi$ (known as the
actor). For one-step actor-critic DRL, the critic $Q(s,a)$ is the expected
reward for state $s$ and action $a$,

$Q^{\pi}(s,a)=\underset{a\sim\pi}{E}\left[r\mid s,a\right].$ (6)

The aim of the actor is selecting an action to maximize its critic value
$Q(s,a)$:

$\pi(a|s)=\arg\max_{a}Q\left(s,\pi(a|s)\right).$ (7)

In practical implementation, both the critic $Q\left(s,a\right)$ and the actor
$\pi(a|s)$ are approximated by neural networks $Q_{\phi}$ and $\pi_{\theta}$
with parameters of $\phi$ and $\theta$. $Q_{\phi}(s,a)$ is learned by
minimizing the MSE loss:

$L_{Q}(\phi)=\frac{1}{|B|}\sum_{\left(s,a,r\right)\in
B}\left(Q_{\phi}(s,a)-r\right)^{2},$ (8)

where $B$ is a mini-batches data $\left(s,a,r\right)$ sampling from the replay
buffer, and $|B|$ is the number of samples in the batch.

The actor $\pi_{\theta}$ is learned by maximizing the loss function:

$L_{\pi}({\theta})=\frac{1}{|B|}\sum_{s\in
B}Q_{\phi}\left(s,\mu_{\theta}(s)\right).$ (9)

Figure 1: The framework of one-step two-critic deep reinforcement learning.

###### Remark 1

The one-step DRL scheme is the simplest version of DDPG or TD3. Setting
$\gamma=0$ for DDPG or TD3 would derive the scheme directly. After setting
$\gamma=0$, the target network, clipped double-Q learning, delayed policy
updates, target policy smoothing, and the soft updating of the target network
is not necessary for DDPG or TD3. Compared with DDPG with 4 neural networks,
and TD3 with 6 neural networks, the scheme only needs 2 neural networks.
However, the scheme would have better performance because it is designed for
the IB-VVC problem specifically.

###  III-B Two-Critic Scheme for One-Step Actor-Critic DRL

(a) Two-bus system. (b) The functions $V(Q_{G})$, $P_{loss}(Q_{G})$.

Figure 2: The relationship between voltage $V$, power loss $P_{loss}$ and
reactive power generation $Q_{G}$ for a two-bus system.

VVC has two objectives: minimizing power loss, and eliminating voltage
violations. For such a problem with multiple objectives, existing DRL methods
integrate them into one objective and use one critic network to learn it.
Neural networks indeed have the universe approximation capability to
approximate any continuous function with arbitrarily small error with enough
number of hidden nodes theoretically [24]. However, using one critic network
to learn the integration of two objectives may degrade its approximation
accuracy and convergence speed in the learning process because of the
following properties of VVC problems:

  * 1)

The reward for active power loss $r_{p}$ and the reward for voltage violation
rate $r_{v}$ have different properties. As shown in Fig. 2, the relationship
between reactive power injection $Q_{G}$ and voltage $V$ at ADNs is close to
linear, whereas the relationship between reactive power injection $Q_{G}$ and
power loss $P_{loss}$ is strongly nonlinear. If we use one critic network to
approximate the reward containing two objectives, the critic network may have
to mix the two functions. In addition, the data in the reply buffer is dynamic
and does not satisfy the independent-identical-distribution condition, which
increases the difficulty further.

  * 2)

Eliminating voltage violations is more important than minimizing power loss,
so $c_{v}$ should be set large enough to penalize voltage violations. However,
the performance of the critic may suffer from the numerical stability problem
for a large $c_{v}$ [18].

  * 3)

Difficulties in learning the function of power loss and voltage violation may
be different. The difficult task costs more time or more data, and the
approximation accuracy may be lower compared with the easier task. One critic
may mix the two tasks, and the worse performance of the difficult task maybe
dominate in the learning process.

To address the problems above, we use two critic networks to approximate the
two objectives separately. The two-critic scheme decreases the approximation
difficulties of each critic, thus having a faster convergence rate and a small
approximation error. Correspondingly, the reward stored in MDP is designed as
$r=(r_{p},r_{v})$.

The critic for power loss $Q_{p}(s,a)$ and voltage violations $Q_{v}(s,a)$ are

$\begin{split}Q_{p}^{\pi}(s,a)&=\underset{a\sim\pi}{E}\left[r_{p}\mid
s,a\right]\\\ Q_{v}^{\pi}(s,a)&=\underset{a\sim\pi}{E}\left[r_{v}\mid
s,a\right].\end{split}$ (10)

In real application, the two critics $Q_{p}\left(s,a\right)$ and
$Q_{p}\left(s,a\right)$ are approximated by two neural networks $Q_{\phi_{p}}$
and $Q_{\phi_{v}}$ with parameters of $\phi_{v}$ and $\phi_{v}$.
$Q_{\phi_{p}}(s,a)$ and $Q_{\phi_{v}}(s,a)$ are learned by minimizing the MSE
losses,

$\begin{split}L_{Q_{p}}(\phi_{p})&=\frac{1}{|B|}\sum_{\left(s,a,r_{p}\right)\in
B}\left(Q_{\phi_{p}}(s,a)-r_{p}\right)^{2}\\\
L_{Q_{v}}(\phi_{v})&=\frac{1}{|B|}\sum_{\left(s,a,r_{v}\right)\in
B}\left(Q_{\phi_{v}}(s,a)-c_{v}r_{v}\right)^{2}.\end{split}$ (11)

The network actor is updated by maximizing the loss function,

$L_{\pi}({\theta})=\frac{1}{|B|}\sum_{s\in
B}\left(Q_{p_{\phi}}\left(s,\pi_{\theta}(s)\right)+Q_{v_{\phi}}\left(s,\pi_{\theta}(s)\right)\right).$
(12)

Algorithm 1 One-Step Two-Critic with Deterministic Policy (OSTC-DP) algorithm

1:Initial policy parameters $\theta$, Q-function parameters $\phi_{p}$,
$\phi_{v}$, empty replay buffer $\mathcal{D}$;

2:for each environment step do

3: Observe state $s$, and select action $a=$

4: $clip(\pi_{\theta}(s)+\epsilon,a_{Low},a_{High})$, where
$\epsilon\sim\mathcal{N}$;

5: Execute $a$ in the environment, and observe reward

6: $r_{p},r_{v}$;

7: Store $(s,a,r_{p},r_{v})$ in replay buffer;

8: if it’s time to update then

9: for $j$ in range (how many updates) do

10: Randomly sample a batch of transitions,

11: $B={(s,a,r_{p},r_{v})}$ from $\mathcal{D}$;

12: Update $Q_{p}$ and $Q_{v}$ by one step of

13: gradient descent using
$\nabla_{\phi_{p}}\frac{1}{|B|}\sum_{\left(s,a,r_{p}\right)\in
B}\left(Q_{\phi_{p}}(s,a)-r_{p}\right)^{2}$
$\nabla_{\phi_{v}}\frac{1}{|B|}\sum_{\left(s,a,r_{v}\right)\in
B}\left(Q_{\phi_{v}}(s,a)-r_{v}\right)^{2}$

14: Update policy by one step of gradient ascent

15: using  $\begin{split}\nabla_{\theta}\frac{1}{|B|}\sum_{s\in
B}\big{(}Q_{p_{\phi}}\left(s,\pi_{\theta}(s)\right)\\\
+c_{v}Q_{v_{\phi}}\left(s,\pi_{\theta}(s)\right)\big{)}\end{split}$

16: end for

17: end if

18:end for

OSTC-DRL is compatible well with any off-policy actor-critic algorithms. For
example, we design an OSTC-DP derived from DDPG to show the approach is
compatible well with deterministic policies. Algorithm 1 provides the detail
of the OSTC-DP. We also design OSTC-SAC derived from SAC to show the approach
is compatible well with stochastic policies. To obtain OSTC-SAC, we need to
make little modifications in steps 2 and 9 in Algorithm 1. In step 2, we need
to replace the deterministic policy
$\pi_{\theta}(s)=\tanh\left(\mu_{\theta}(s)\right)$ as the stochastic policy
$\pi_{\theta}(s,\xi)=\tanh\left(\mu_{\theta}(s)+\sigma_{\theta}(s)\odot\xi\right),\quad\xi\sim\mathcal{N}(0,I)$,
where $\mu_{\theta}$ , $\sigma_{\theta}$ are neural networks. In step 9, we
need to add the entropy regularization term,
$\nabla_{\theta}\frac{1}{|B|}\sum_{s\in
B}\big{(}Q_{p_{\phi}}\left(s,\pi_{\theta}(s)\right)+c_{v}Q_{v_{\phi}}\left(s,\pi_{\theta}(s)\right)-\alpha\log\pi_{\theta}\left(\pi_{\theta}(s)\mid
s\right)$. The temperature $\alpha$ can be adjusted by minimizing the loss
$L(\alpha)=\frac{1}{|B|}\sum_{s\in
B,a\sim\pi_{\theta}}[-\alpha\log\pi_{\theta}(a|s)-\alpha\mathcal{H}]$, where
$\mathcal{H}$ is the entropy target.

##  IV Extending OSTC-DRL to Multi-agent OSTC-DRL for Decentralized IB-VVC

The centralized OSTC-DRL approach requires massive real-time communication and
is fragile for single-point failure. For ADNs that each sub-area only can
acquire local information in real-time, and to enhance the robustness against
communication failures, we extend the OSTC-DRL to the multi-agent form for
decentralized IB-VVC. It is based on a centralized training decentralized
execution approach. In off-policy DRL algorithms, both actor networks and
critic networks are trained by the data sampling from the data buffer, so the
interaction data stored in the data buffer with a delayed time has little
influence in the training stage. In executing stage, each actor-network of the
sub-area just requires real-time local observation to make the decision. The
framework of multi-agent OSTC-DRL is shown in Fig. 3.

###  IV-A Formulating VVC as One-step Markov Game

Markov game is an extension of MDP for a multi-agent system. Similar to the
one-step MDP, we formulate the IB-VVC as a one-step Markov game where multiple
agents interact with the same environment. We use a tuple
$\left(\mathcal{S},\left[\mathcal{O}_{i}\right]_{n},\left[\mathcal{A}_{i}\right]_{n},\mathcal{R}\right)$
to describe a Markov game with $n$ agent, where $s\in\mathcal{S}$ is the full
state of environment, $o_{i}\in\mathcal{O}_{i}$ and $a_{i}\in\mathcal{A}_{i}$
are the local observations and actions for each agent $i$, and
$[r_{p},r_{v}]^{T}\in\mathcal{R}$ are the reward functions relating power loss
and voltage violation rate that are defined as
$\mathcal{S}\times\mathcal{A}_{1}\times\dots\mathcal{A}_{N}\times\mapsto\mathcal{R}_{p},\mathcal{R}_{v}$.
The goal of each agent is to maximize the expected reward
$J=E_{a_{i}\sim\pi_{i}(o_{i})}(r_{p}(s,a_{1},a_{2},\dots,a_{n})+c_{v}r_{v}(s,a_{1},a_{2},\dots,a_{n}))$.

The definition of state space and the reward is the same as the centralized
version in section II. The local observations of each agent is
$o_{i}=\left(P_{i},Q_{i},Q_{G,i},Q_{C,i}\right)$, where $o_{i}$ is the local
information of of $i^{th}$ subarea. The selection of local observation space
depends on the measurement conditions of each sub-area. It can be the
measurement of one sub-area and its neighbor areas, or just one bus of the IB-
ER or SVC installed. The local action of each agent is
${a}_{i,t}=\left(Q_{G,i},Q_{C,i}\right)$, where $Q_{G,i},Q_{C,i}$ are the
controllable reactive power resources in $i^{th}$ area, including IB-ERs and
SVCs.

Figure 3: The framework of multi-agent one-step two-critic deep reinforcement
learning.

###  IV-B Multi-Agent OSTC

As shown in Fig. 3, multi-agent OSTC-DRL contains two stages: centralized
training and decentralized execution. In the centralized training stage, the
off-policy multi-agent OSTC-DRL algorithm samples data from the data buffer to
train actors and critics. In the decentralized execution stage, the actor
generates actions and then applies the actions to the environment. The actions
generated by actors are based on real-time local measurements. The interaction
data will be stored in the data buffer, and the data storage can be at a slow
rate.

Since the data buffer can obtain the full state of ADNs and actions of each
actor at a slow rate, one critic $Q(s,a_{1},\dots,a_{n})$ is enough to
represent the state-action value of all actors. Similar to the OSTC-DRL, we
design two critics that are the critic of power loss
$Q_{p}(s,a_{1},\dots,a_{n})$ and the critic of voltage violation
$Q_{v}(s,a_{1},\dots,a_{n})$,

$\begin{split}Q_{p}^{\pi}(s,a_{1},\dots,a_{n})&=\underset{a_{i}\sim\pi_{i}}{E}\left[r_{p}\mid
s,a_{1},\dots,a_{n}\right]\\\
Q_{v}^{\pi}(s,a_{1},\dots,a_{n})&=\underset{a_{i}\sim\pi_{i}}{E}\left[r_{v}\mid
s,a_{1},\dots,a_{n}\right].\end{split}$ (13)

The critic networks $Q_{\phi_{p}}(s,a_{1},\dots,a_{n})$ and
$Q_{\phi_{v}}(s,a_{1},\dots,a_{n})$ are the same as the centralized OSTC-DRL,
which are learned by minimizing the MSE losses,

$\begin{split}L_{Q_{p}}(\phi_{p})&=\frac{1}{|B|}\sum_{\left(s,a_{1},\dots,a_{n},r_{p}\right)\in
B}\left(Q_{\phi_{p}}(s,a_{1},\dots,a_{n})-r_{p}\right)^{2}\\\
L_{Q_{v}}(\phi_{v})&=\frac{1}{|B|}\sum_{\left(s,a_{1},\dots,a_{n},r_{v}\right)\in
B}\left(Q_{\phi_{v}}(s,a_{1},\dots,a_{n})-r_{v}\right)^{2}.\end{split}$ (14)

The network actors are updated by the loss function,

$\begin{split}L_{\pi}(\theta_{1},\dots\theta_{n})=\frac{1}{|B|}\sum_{s\in
B}\big{(}Q_{p_{\phi}}\left(s,\pi_{\theta_{1}}(o_{1}),\dots,\pi_{\theta_{n}}(o_{n})\right)\\\
+c_{v}Q_{v_{\phi}}\left(s,\pi_{\theta_{1}}(o_{1}),\dots,\pi_{\theta_{n}}(o_{n})\right)\big{)}.\end{split}$
(15)

Algorithm 2 shows the details of multi-agent OSTC-DP. In step 6, we need to
collect the $(s,a_{1},\dots,a_{n},r_{p},r_{v})$ of the ADNs at the same time.
It is accessible because recently measurement devices can add a timestamp to
the measured data, and then the data can be updated to the center with a slow
time rate or a constant time interval decay. Of course, those data may be non-
synchronized, additional methods [25] can be used to preprocess the non-
synchronized data.

To obtain multi-agent OSTC-SAC, we need to make little modifications in steps
3 and 10 in Algorithm 2. In step 3, we need replace the policy as the
stochastic policy
$\pi_{\theta_{i}}(o_{i},\xi)=\tanh\left(\mu_{\theta_{i}}(o_{i})+\sigma_{\theta_{i}}(o_{i})\odot\xi_{i}\right),\quad\xi_{i}\sim\mathcal{N}(0,I)$,
where $\mu_{\theta_{i}}$ , $\sigma_{\theta_{i}}$ are neural networks. In step
10, we need add the entropy regularization term to the critic networks,
$L_{\pi}(\theta_{1},\dots\theta_{n})=\frac{1}{|B|}\sum_{s\in
B}\big{(}Q_{p_{\phi}}\left(s,\pi_{\theta_{1}}(o_{1}),\dots,\pi_{\theta_{n}}(o_{n})\right)+c_{v}Q_{v_{\phi}}\left(s,\pi_{\theta_{1}}(o_{1}),\dots,\pi_{\theta_{n}}(o_{n})\right)-\sum_{i}^{n}\alpha\log\pi_{\theta_{i}}\left(\pi_{\theta}(s)\mid
s\right)\big{)}$. The temperature $\alpha$ can be adjusted by minimizing the
loss $L(\alpha)=\frac{1}{|B|}\sum_{s\in
B,a\sim\pi_{\theta}}[-\alpha\log\pi_{\theta}(a|s)-\alpha\mathcal{H}]$, where
$\mathcal{H}$ is the entropy target. We name the algorithm derived from OSTC-
SAC as multi-agent OSTC-SAC.

Multi-agent OSTC-DP and multi-agent OSTC-SAC are flexible for the measurement
conditions of ADNs. The algorithms are designed for the measurement conditions
that each actor executes based on the local single bus information or sub-area
information.

Algorithm 2 Multi-agent One-Step Two-Critic with Deterministic Policies
(Multi-agent OSTC-DP) algorithm

1:Initial policy parameters $\theta_{1},\dots\theta_{n}$, Q-function
parameters $\phi_{1}$, $\phi_{2}$, empty replay buffer $\mathcal{D}$;

2:for each environment step do

3: for each agent $i$ do

4: Observe $o_{i}$, and select action $a_{i}=$

5: $clip(\pi_{\theta_{i}}(o_{i})+\epsilon_{i},a_{Low},a_{High})$, where
$\epsilon_{i}\sim\mathcal{N}$;

6: Execute $a_{i}$ in the environment;

7: end for

8: Collect reward $r_{p}$, $r_{v}$ of the ADN;

9: Store $(s,a_{1},\dots,a_{n},r_{p},r_{v})$ in replay buffer;

10: if it’s time to update then

11: for $j$ in range (how many updates) do

12: Randomly sample a batch of transitions, $B=$

13: $(s,a_{1},\dots,a_{n},r_{p},r_{v})$ from $\mathcal{D}$;

14: Update $Q_{p}$ and $Q_{v}$ by one step of gradient

15: descent using Equ. (14);

16: Update policies by one step of gradient ascent

17: using Equ. (15);

18: end for

19: end if

20:end for

##  V Simulation

Numerical simulation was conducted on 33-bus [26] and 69-bus [27] test
distribution networks to demonstrate the advantages of the proposed OSTC-DRL
approach. In the 33-bus test distribution network, 3 IB-ERs of 2 MVar reactive
power capacity and 1.5 MW active power were connected to buses 18, 22, and 25,
respectively, and 1 SVC of 2 MVar reactive power capacity was connected to bus
33. In the 69-bus test distribution network, 4 IB-ERs of 2 MVar reactive power
capacity and 1.5 MW active power were connected to buses 6, 24, 45, and 58
respectively, and 1 SVC of 2 MVar reactive power capacity was connected to bus
14. All load and generation levels were multiplied with the fluctuation ratio
of one day with 96 points extracted from paper [16] and a $20\%$ uniform
distribution noise to reflect the variance. The voltage limits for all buses
were set to be [0.95, 1.05]. The algorithms were implemented in Python. The
balanced power flow was solved by Pandapower [28] to simulate ADNs, and the
implementation of the DRL algorithms used PyTorch.

###  V-A Simulations for the Centralized OSTC-DRL Approach

We performed simulations to understand the contribution of each component:
one-step and two-critic, and show the superiority of the proposed centralized
OSTC-DRL approach. We designed 3 classes of simulation experiments.

  * Deterministic policy: 1) DDPG [14]; 2) One-step with deterministic policy (OS-DP) derived from the DDPG by applying the one-step DRL scheme; 3) Two-critic with deterministic policy (TC-DP) derived from the DDPG by applying the two-critic technology; 4) One-step two-critic DRL with deterministic policy (OSTC-DP) derived from the DDPG by applying the one-step DRL scheme and two-critic technology.

  * Stochastic policy: 5) SAC [29]; 6) One-step with stochastic policy (OS-SAC) derived from the SAC by applying the one-step DRL scheme; 7) Two-critic with stochastic policy (TC-SAC) derived from the SAC by applying the two-critic technology; 8) One-step two-critic DRL with stochastic policy (OSTC-SAC) derived from the SAC by applying the one-step DRL scheme and two-critic technology.

  * Model-based: 9) Model-based optimization method with accurate power flow model. Model-based optimization was solved by recalling PandaPower with the interior point solver.

TABLE I: Parameter setting for the reinforcement learning algorithm Algo. | Parameter | Value
---|---|---
Shared | Optimizer | Adam
Activation function | ReLU
Number of hidden layers | 2
Actor hidden layer neurons | {512, 512 }
Critic hidden layer neurons | {512, 512 }
Batch size | 128
Replay buffer size | $3\times 10^{4}$
Critic learning rate | $3\times 10^{-4}$
Actor learning rate | $1\times 10^{-4}$
Voltage violation penalty $c_{v}$ | 50
Initial random step | 960
Iterations per time step | 4
The deterministic policy | Exploration Policy | $\mathcal{N}(0,0.1)$
The stochastic policy | Entropy target | $-\dim(\mathcal{A})$
Temperature learning rate | $3\times 10^{-4}$
TABLE II: Quantified indices of the benchmarks in the final 50 episodes for
centralized DRL algorithms

Algorithm | Reward | $P_{\text{loss }}/\mathrm{MW}$ | VVR/p.u.2
---|---|---|---
33-bus | 69-bus | 33-bus | 69-bus | 33 -bus | 69-bus
Model-based | MBO1 | $-4.199$ | $-3.957$ | $4.199$ | $3.957$ | 0 | 0
Deterministic policy | DDPG | $-4.898$ | $-4.701$ | $4.898$ | $4.622$ | $4.882\text{e-}6$ | $1.580\text{e-}3$
OS-DP | $-4.867$ | $-4.867$ | $4.867$ | $4.778$ | $1.129\text{e-}5$ | $1.782\text{e-}3$
| TC-DP | $-4.289$ | $-4.217$ | $4.288$ | $4.140$ | $2.642\text{e-}5$ | $1.533\text{e-}3$
| OSTC-DP | -4.270 | -4.161 | $4.268$ | $4.104$ | $3.988\text{e-}5$ | $1.137\text{e-}3$
Stochastic policy | SAC | $-4.548$ | $-4.491$ | $4.548$ | $4.471$ | $0$ | $3.859\text{e-}4$
OS-SAC | $-4.368$ | $-4.383$ | $4.357$ | $4.339$ | $2.00\text{e-}4$ | $8.800\text{e-}4$
| TC-SAC | $-4.384$ | $-4.275$ | $4.373$ | $4.205$ | $2.184\text{e-}4$ | $1.396\text{e-}3$
| OSTC-SAC | -4.272 | -4.191 | $4.268$ | $4.121$ | $9.205\text{e-}5$ | $1.402\text{e-}3$

  * 1

MBO is the model-based optimization method using an accurate ADN model, which
can be seen as the optimal result. However, the accurate ADN model is not
available in real applications.

  * 2

VVR is the daily accumulation voltage violation rate.

Figure 4: Testing results of the training stage for the 33-bus test
distribution network. Figure 5: Testing results of the training stage for the
69-bus test distribution network.

(a) The contribution of one-step

(b) The contribution of two-critic

(c) The contribution of one-step two-critic

Figure 6: The reward error in the final 50 episodes (reward error = reward of
MBO - reward of the mentioned method)

The result of model-based optimization can be seen as optimal, which is a
baseline for the performance of DRL algorithms. We trained the DRL agent using
300 days of data. The parameter setting for four DRL algorithms is provided in
Table I. The discount factor $\gamma$ for DDPG, TC-DP, SAC, TC-SAC is 0.9. In
the training process, we tested the DRL algorithms in the same environment at
each step.

The testing results in the training process are shown in Figs. 4 and 5. The
reward, power loss, and voltage violation rate in those figures are the daily
accumulation values. To show the superiority of the proposed OSTC-DRL approach
over the traditional DRL approaches clearly, we omitted the results of OS-DP,
TC-DP, OS-SAC, and TC-SAC in Figs. 4 and 5. We made the following two
observations.

First, the proposed OSTC-DP and OSTC-SAC converged faster compared with DDPG
and SAC, see the learning trajectory from days 20-50 in Figs. 4 and 5.

Second, the proposed OSTC-DP and OSTC-SAC achieved higher rewards compared
with DDPG and SAC, see the learning trajectory from days 250-300 in Figs. 4
and 5. OSTC-DP and OSTC-SAC also had smaller power losses.

The voltage violation rates of all DRL algorithms had small fluctuations
around zeros. Increasing the voltage violation penalty can decrease the
voltage violation rate, but it just alleviates the fluctuations and cannot
avoid voltage violation completely. DRL algorithms learn by trial and error,
so in the training trajectory, DRL algorithms must trail both sides of the
voltage boundary many times to find the optimal solution. The alternative way
to address the voltage violation issues is by tightening the voltage limits.
For example, for the normal voltage operation interval $[0.95,1.05]$, we set
the objective voltage interval for DRL algorithms is $[0.955,1.045]$. Even
though there are slight voltage violations for the interval $[0.955,1.045]$ in
the learning process, there is no voltage violation for the interval
$[0.95,1.05]$.

To quantify the advantages of the proposed OSTC-DRL approach, Table II gives
the converged results of the 9 methods. We use the accuracy equation
$Acc=(R_{i}-R_{\text{MBO}})/R_{\text{MBO}}$ where $i$ represent OSTC-DP, DDPG,
OSTC-SAC, or SAC, and MBO represents model based optimization method using an
accurate power flow model. From the perspective of reward, the accuracy of
OSTC-DP, DDPG, OSTC-SAC, and SAC were $1.694\%$, $16.65\%$, $1.742\%$, and
$8.304\%$ in the 33-bus network, $5.133\%$, $18.79\%$, $5.891\%$, $13.47\%$ in
the 69-bus network. We can see that the accuracy of OSTC-DP was $9.831$ times
as DDPG in the 33-bus network, and $3.660$ times as DDPG in the 69-bus
network. The accuracy of OSTC-SAC was $4.766$ times as SAC in the 33-bus
network, and $2.286$ times as SAC in the 69-bus network. The results showed
that the OSTC-DRL approach can improve the performance of DRL algorithms
considerably for both deterministic policies and stochastic policies.

To show the contribution of the components of one-step and two-critic, Fig. 6
gives the ablation study results. The contribution of the component of one-
step is shown in Fig. 6(a). The one-step improved the performance of those
algorithms except that compares the results of DDPG and OS-DP for the 69-bus
network. The contribution of the component of the two-critic is shown in Fig.
6(b). Two-critic improved the performance of all of those algorithms. The
contribution of the combination of the one-step and the two-critic is shown in
Fig. 6(c). The algorithms with the OSTC-DRL approach achieved the best
performance among the algorithms.

Generally, SAC has better control performance than the DDPG because SAC has
two extra components: clipped double-Q learning, and entropy regularization
[29, 15]. However, our simulation results showed that the performances of TC-
DP and OSTC-DP were better than the corresponding of TC-SAC and OSTC-SAC,
respectively. The reasons may be as follows. “clipped double-Q learning”
mitigated the overestimation of Q value whereas bringing an underestimation
bias [13, 12]. Entropy regularization accelerated the learning process and
prevented the policy from prematurely converging to a bad local optimum.
However, it may bring additional regularization errors when the estimation
accuracy of the Q function is high enough. When the estimation accuracy of the
Q value was not high, clipped double-Q learning and entropy regularization
brought positive influence and led to the final results of SAC and OS-SAC over
the corresponding DDPG and OS-DP. However, when the estimation accuracy of the
Q value was high, clipped double-Q learning and entropy regularization brought
negative influence and led to the final results of TC-DP and OSTC-DP over the
corresponding TC-SAC and OSTC-SAC, respectively.

TABLE III: Quantified indices of the benchmarks in the final 50 episodes for
multi-agent DRL algorithms

Algorithm | Reward | $P_{\text{loss }}/\mathrm{MW}$ | VVR/p.u.
---|---|---|---
33-bus | 69-bus | 33-bus | 69-bus | 33 -bus | 69-bus
Deterministic policy | OSTC-DP | $-4.270$ | $-4.161$ | $4.268$ | $4.104$ | $3.988\text{e-}5$ | $1.137\text{e-}3$
MA-OSTC-DP-sub1 | $-4.264$ | $-4.142$ | $4.262$ | $4.115$ | $2.787\text{e-}5$ | $5.364\text{e-}4$
MA-OSTC-DP-local2 | $-4.301$ | $-4.397$ | $4.288$ | $4.269$ | $2.624\text{e-}4$ | $2.568\text{e-}3$
Stochastic policy | OSTC-SAC | $-4.272$ | $-4.191$ | $4.268$ | $4.121$ | $9.205\text{e-}5$ | $1.402\text{e-}3$
MA-OSTC-SAC-sub | $-4.343$ | $-4.234$ | $4.317$ | $4.128$ | $5.284\text{e-}4$ | $2.117\text{e-}3$
MA-OSTC-SAC-local | $-4.370$ | $-4.422$ | $4.331$ | $4.282$ | $7.828\text{e-}4$ | $2.797\text{e-}3$

  * 1

“MA” means multi-agent. ”sub” represents the sub-area bus information. “MA-
OSTC-DP-sub” means that the MA-OSTC-DP-sub algorithm works on that each actor
executes based on the sub-area bus information.

  * 2

“local” represents the local one bus information. “MA-OSTC-DP-sub” means that
the MA-OSTC-DP-sub algorithm works on that each actor executes based on the
local one bus information.

###  V-B Simulation for the Multi-Agent OSTC-DP Algorithm

The proposed multi-agent algorithms were designed for measurement conditions
of ADNs in that all the measurements are uploaded to the center at a slow rate
and each agent only can obtain its sub-area measurements in real-time. We
tested the multi-agent algorithms on the two measurement conditions: 1) each
agent can obtain the local one bus information, and 2) each agent can obtain
the sub-area bus information. This setting was to show the multi-agent OSTC-
DRL is flexible to different measurement situations.

For the first measurement conditions, the 33-bus test distribution network was
divided into 4 subareas and the 69-bus test distribution network was divided
into 5 subareas. Each subarea contained one bus in which the controllable
devices have been installed. For the second measurement condition, we divided
the distribution networks into 4 sub-areas for both case 33 and 69 systems. In
the 33-bus test distribution network, the sub-areas were $[7,8,\dots,18]$,
$[19,20,21,22]$, $[23,24,25]$ and $[26,27,\dots,33]$. In the 69-bus test
distribution network, the sub-areas were $[2,3,\dots,10]$, $[11,12,\dots,26]$,
$[36,37,\dots,45]$ and $[53,54,\dots,64]$. The partitioning is flexible. Some
buses can belong to two partitionings concurrently, or not belong to any
partitioning.

Correspondingly, we extended OSTC-DP to multi-agent OSTC-DP-local and multi-
agent OSTC-DP-sub for two measurement conditions. Also, we extended OSTC-SAC
to multi-agent OSTC-SAC-local and multi-agent OSTC-SAC-sub. The parameters of
DRL algorithms were the same as the corresponding subsection V-A except for
the number of actors. We tested the performance of DRL algorithms at each step
in the training process. Table III shows the qualified indices of the 6
algorithms in the final 50 episodes.

After enough time to learn, for deterministic policies, multi-agent OSTC-DP-
sub achieved similar performance as OSTC-DP, while the performance of multi-
agent OSTC-DP-local was slightly worse than the performance of OSTC-DP. For
stochastic policies, the rank of the performance of three algorithms from high
to low was OSTC-SAC, multi-agent OSTC-SAC-sub, and multi-agent OSTC-SAC-local.
Noting that the performances of all multi-agent OSTC DRL algorithms are better
than the performances of DDPG and SAC. Those results showed that the multi-
agent OSTC-DP and OSTC-SAC algorithms are robust for the information obtained
by each actor in the execution stage, even when each actor just can obtain its
local one bus information. It is reasonable because the voltage information of
one bus is influenced by other buses, so it can reflect the global information
partially. Meanwhile, the partial information degrades the control performance
of DRL slightly.

##  VI Conclusion

In this paper, we have proposed an OSTC-DRL approach for IB-VVC in ADNs. Based
on the OSTC-DRL approach, we designed two DRL algorithms that are OSTC-DP and
OSTC-SAC. We also extended the approach to the multi-agent OSTC-DRL approach
for decentralized IB-VVC problems. We designed multi-agent OSTC-DP and multi-
agent OSTC-SAC algorithms. Simulation results showed the contributions of one-
step and two-critic separately, and the OSTC-DRL approach has improved the VVC
performance considerably compared with the state-of-the-art DRL algorithms for
IB-VVC in ADNs. After extending to the multi-agent DRL algorithms, they
achieved nearly equal or slight degradation performance as the OSTC-DRL
algorithms.

The proposed OSTC-DRL approach focuses on single-period optimization problems
and the action space is continuous. However, for ADNs embedded with capacity
banks, on-load tap changers, and storage devices, the actions contain both
continuous and discrete, and the optimization task should consider the long
horizontal process. Therefore, we would extend our algorithm to mixed-integer
multi-period optimization problems in future works.

## References

  * [1] H. Liu, C. Zhang, Q. Chai, K. Meng, Q. Guo, and Z. Y. Dong, “Robust regional coordination of inverter-based volt/var control via multi-agent deep reinforcement learning,” _IEEE Trans. on Smart Grid_ , vol. 12, no. 6, pp. 5420–5433, Nov. 2021.
  * [2] Y. Tang, K. Dvijotham, and S. Low, “Real-time optimal power flow,” _IEEE Trans. on Smart Grid_ , vol. 8, no. 6, pp. 2963–2973, Nov. 2017.
  * [3] R. Albert, I. Albert, and G. L. Nakarado, “Structural vulnerability of the north american power grid,” _Phys. Rev. E_ , vol. 69, p. 025103, Feb 2004\.
  * [4] H. Liu and W. Wu, “Bi-level off-policy reinforcement learning for volt/var control involving continuous and discrete devices,” _arXiv preprint arXiv:2104.05902_ , 2021.
  * [5] D. Cao, W. Hu, J. Zhao, Q. Huang, Z. Chen, and F. Blaabjerg, “A multi-agent deep reinforcement learning based voltage regulation using coordinated PV inverters,” _IEEE Trans. on Power Syst._ , vol. 35, no. 5, pp. 4120–4123, Sept. 2020.
  * [6] P. Kou, D. Liang, C. Wang, Z. Wu, and L. Gao, “Safe deep reinforcement learning-based constrained optimal control scheme for active distribution networks,” _Appl. Energy_ , vol. 264, p. 114772, Apr. 2020.
  * [7] Y. Gao, W. Wang, and N. Yu, “Consensus multi-agent reinforcement learning for volt-var control in power distribution networks,” _IEEE Trans. on Smart Grid_ , vol. 12, no. 4, pp. 3594–3604, July 2021.
  * [8] H. Liu and W. Wu, “Online multi-agent reinforcement learning for decentralized inverter-based volt-var control,” _IEEE Trans. on Smart Grid_ , vol. 12, no. 4, pp. 2980–2990, July 2021.
  * [9] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [10] S. Thrun and A. Schwartz, “Issues in using function approximation for reinforcement learning,” in _Proceedings of the 1993 Connectionist Models Summer School_ , vol. 6. Hillsdale, NJ, USA: Lawrence Erlbaum, 1993.
  * [11] S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _Int. Conf. Mach. Learn. (ICML)_. Stockholm, Sweden: PMLR, 2018, pp. 1587–1596.
  * [12] L. Pan, Q. Cai, and L. Huang, “Softmax deep double deterministic policy gradients,” _arXiv preprint arXiv:2010.09177_ , 2020.
  * [13] K. Ciosek, Q. Vuong, R. Loftin, and K. Hofmann, “Better exploration with optimistic actor-critic,” _arXiv preprint arXiv:1910.12807_ , 2019.
  * [14] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning.” in _Int. Conf. Learn. Representations (ICLR)_ , Vancouver, BC, Canada, 2016.
  * [15] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _Int. Conf. Mach. Learn. (ICML)_. Stockholm, Sweden: PMLR, 2018, pp. 1861–1870.
  * [16] H. Liu and W. Wu, “Two-stage deep reinforcement learning for inverter-based volt-var control in active distribution networks,” _IEEE Trans. on Smart Grid_ , vol. 12, no. 3, pp. 2037–2047, May 2021.
  * [17] W. Wang, N. Yu, Y. Gao, and J. Shi, “Safe off-policy deep reinforcement learning algorithm for volt-var control in power distribution systems,” _IEEE Trans. on Smart Grid_ , vol. 11, no. 4, pp. 3008–3018, July 2020.
  * [18] Y. Zhang, X. Wang, J. Wang, and Y. Zhang, “Deep reinforcement learning based volt-var optimization in smart distribution systems,” _IEEE Trans. on Smart Grid_ , vol. 12, no. 1, pp. 361–371, Jan. 2020.
  * [19] Z. Yan and Y. Xu, “Real-time optimal power flow: A lagrangian based deep reinforcement learning approach,” _IEEE Trans. on Power Syst._ , vol. 35, no. 4, pp. 3270–3273, July 2020.
  * [20] X. Sun and J. Qiu, “Two-stage volt/var control in active distribution networks with multi-agent deep reinforcement learning method,” _IEEE Trans. on Smart Grid_ , vol. 12, no. 4, pp. 2903–2912, July 2021.
  * [21] M. Majzoubi, C. Zhang, R. Chari, A. Krishnamurthy, J. Langford, and A. Slivkins, “Efficient contextual bandits with continuous actions,” in _Adv. Neural Inf. Process. Syst. (NeurIPS)_ , H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, Eds., vol. 33, Vancouver, Canada, 2020, pp. 349–360.
  * [22] Q. Yang, G. Wang, A. Sadeghi, G. B. Giannakis, and J. Sun, “Two-timescale voltage control in distribution grids using deep reinforcement learning,” _IEEE Trans. on Smart Grid_ , vol. 11, no. 3, pp. 2313–2323, May 2020.
  * [23] V. Kekatos, G. Wang, A. J. Conejo, and G. B. Giannakis, “Stochastic reactive power management in microgrids with renewables,” _IEEE Trans. on Power Syst._ , vol. 30, no. 6, pp. 3386–3395, Nov. 2015.
  * [24] T. Chen and H. Chen, “Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems,” _IEEE Trans. Neural Netw._ , vol. 6, no. 4, pp. 911–917, Jul. 1995\.
  * [25] A. Alimardani, F. Therrien, D. Atanackovic, J. Jatskevich, and E. Vaahedi, “Distribution system state estimation based on nonsynchronized smart meters,” _IEEE Trans. on Smart Grid_ , vol. 6, no. 6, pp. 2919–2928, Nov. 2015.
  * [26] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Power Energy Mag._ , vol. 9, no. 4, pp. 101–102, Apr. 1989.
  * [27] D. Das, “Optimal placement of capacitors in radial distribution system using a fuzzy-GA method,” _Int. J. Electr. Power Energy Syst._ , vol. 30, no. 6-7, pp. 361–367, July 2008.
  * [28] L. Thurner, A. Scheidler, F. Schäfer, J. Menke, J. Dollichon, F. Meier, S. Meinecke, and M. Braun, “pandapower — an open-source python tool for convenient modeling, analysis, and optimization of electric power systems,” _IEEE Trans. on Power Syst._ , vol. 33, no. 6, pp. 6510–6521, Nov. 2018.
  * [29] T. Haarnoja, A. Zhou, K. Hartikainen, G. Tucker, S. Ha, J. Tan, V. Kumar, H. Zhu, A. Gupta, P. Abbeel _et al._ , “Soft actor-critic algorithms and applications,” _arXiv preprint arXiv:1812.05905_ , 2018.


# Gravity aided navigation using Viterbi map matching algorithm

Wenchao Li
School of Science
RMIT University
Melbourne, VIC, Australia
<EMAIL_ADDRESS>
& Christopher Gilliam
School of Engineering
RMIT University
Melbourne, VIC, Australia

& Xuezhi Wang
School of Science
RMIT University
Melbourne, VIC, Australia

& Allison Kealy
School of Science
RMIT University
Melbourne, VIC, Australia

& Andrew D. Greentree
Australian Research Council Centre of
Excellence for Nanoscale BioPhotonics
School of Science
RMIT University
Melbourne, VIC, Australia

& Bill Moran
Department of Electrical and
Electronic Engineering
University of Melbourne
Melbourne, VIC, Australia

###### Abstract

In GNSS-denied environments, aiding a vehicle’s inertial navigation system
(INS) is crucial to reducing the accumulated navigation drift caused by sensor
errors (e.g. bias and noise). One potential solution is to use measurements of
gravity as an aiding source. The measurements are matched to a geo-referenced
map of Earth’s gravity in order to estimate the vehicle’s position. In this
paper, we propose a novel formulation of the map matching problem using a
hidden Markov model (HMM). Specifically, we treat the spatial cells of the map
as the hidden states of the HMM and present a Viterbi style algorithm to
estimate the most likely sequence of states, i.e. most likely sequence of
vehicle positions, that results in the sequence of observed gravity
measurements. Using a realistic gravity map, we demonstrate the accuracy of
our Viterbi map matching algorithm in a navigation scenario and illustrate its
robustness compared to existing methods.

_Keywords_ Gravity aided navigation $\cdot$ Gravity map matching $\cdot$
Hidden Markov Model $\cdot$ Maximum a posterior (MAP) estimation

##  1 Introduction

Since the advent of Global Navigation Satellite Systems (GNSS) in the 1960s,
GNSS-based positioning has become ubiquitous in the areas of navigation,
guidance and control. However, in a growing number of situations GNSS is
either unavailable, e.g. in underwater environments (Wang et al., 2017a), or
deliberately jammed or spoofed (Hemann et al., 2016). In such GNSS-denied (or
contested) situations, a platform’s navigation performance is governed by on-
board sensors – chief among them being the inertial navigation system (INS),
which integrates inertial measurements from accelerometers and gyroscopes to
obtain position. Unfortunately, even with high precision INSs, very small
errors in the inertial sensors result in the build up of large navigation
errors over time (Titterton & Weston, 2004; Groves, 2013). Thus, an on-board
aiding source that provides a position fix is crucial to removing these
accumulated navigation errors and retaining confidence in the navigation
accuracy. A promising candidate, and the focus of this paper, is the use of
gravitational information to obtain a position fix (Han et al., 2016; Liu, Li
& Jing, 2019; Wang et al., 2022).

Gravity aided navigation originates from the early 1990’s (Affleck &
Jircitano, 1990; Jircitano & Dosch, 1991). The central concept is that an on-
board sensor measures elements of the gravitational vector (or gravity
gradient tensor) whilst the platform is in motion and these measurements are
then matched to a geo-referenced map of the Earth’s gravitational field to
determine a position. The advantage of such an approach is twofold: First,
gravity is an inherent physical property of the Earth, thus it is immune to
environmental interference, jamming and spoofing. Second, the measurement of
gravity is a passive operation that does not require external information,
which is advantageous in sensitive situations (Muniraj & Farhood, 2017). The
challenge in this approach however is how to match the measurements to a
position in the map. The design of this matching procedure, known as the map
matching algorithm, is key to the success of gravity aided navigation.

The map matching algorithm needs to take into account: 1) the sensor noise
when measuring gravity; 2) the uncertainty on the spatial location when a
measurement is obtained; 3) possible ambiguity in the map resulting in a
measurement being matched to multiple locations. In the literature, map
matching algorithms tackle these problems in one of two ways. The first
category of algorithms consider the problem as a single point matching
operation. Building on the framework introduced in the Sandia terrain-aided
navigation system (SITAN) (Hostetler & Andreas, 1983; Bergman, 1997), these
techniques use the geo-referenced map as a look-up function and compute a
predicted gravitational measurement using the INS estimated position. The
predicted measurement is then used in the prediction step of a statistical
filter. Due to the non-linear relationship between the estimation states and
gravitational measurements, the extended Kalman filters (EKF) is often used to
perform the estimation (Wei et al., 2017; Lee et al., 2015). This category of
algorithms take into account the sensor noise and spatial uncertainty but do
not consider the map ambiguity directly. The second category of map matching
algorithms directly considers the structure of the map. They approach the
problem as a sequential matching operation. The classical example is the
iterative closest contour point (ICCP) algorithm proposed in (Kamgar-Parsi &
Kamgar-Parsi, 1999) and utilised in (Wang et al., 2016; Han et al., 2017; Liu,
Li, Lin & Jing, 2019). The ICCP is based on the observation that a single,
scalar, gravity measurement creates an iso-contour of similar values in the
gravity map space. Thus, given a sequence of gravity measurements, the
matching problem is posed in terms of fitting a trajectory to a set of iso-
contours based on initial position estimates and sensor measurements. To make
the problem well posed, a regularization term is introduced to control the
shape of the trajectory. However, linking the kinematic constraints of the
platform’s motion to the regularization term is not straightforward.

In this paper, we propose a novel sequential map matching algorithm based on
formulating the matching operation as a hidden Markov model (HMM).
Specifically, we consider the spatial cells of the map to be the set of all
possible hidden states in the model and the platform’s trajectory to be a
particular sequence of these states. The observations from each state are the
gravitational measurements and the output and state transition probabilities
relate to the sensor noise and platform movement models, respectively. Thus,
the map matching problem is equivalent to determining this sequence of states
given a sequence of gravitational measurements. The advantage to this
formulation is twofold: First, our proposed approach works directly in the map
space and incorporates the modelling of the sensor noise and platform motion
with the underlying structure of the gravity map. Second, this type of HMM
problem can be solved efficiently using a dynamic programming method known as
the Viterbi algorithm (Viterbi, 1967; Ren & Karimi, 2009). Given the HMM
formulation, the Viterbi algorithm estimates the most likely sequence of
states that results in the sequence of observations and this sequence of
states is optimal in the maximum a posterior (MAP) sense. Accordingly, we
present a Viterbi map matching algorithm that outputs the most likely
trajectory of the platform given a sequence of gravitational measurements.
Furthermore, to tackle limited map resolution, we implement a two-layer,
coarse-to-fine, estimation scheme to achieve sub-cell accuracy. We validate
our algorithm hsing the ultra-high resolution, non-parametric, gravity maps
presented in (Hirt et al., 2013) and demonstrate that it outperforms existing
algorithms in terms of robustness and navigational accuracy. In particular, we
demonstrate that our algorithm is able to tackle varying map spatial
resolution as well as varying sensor noise levels.

##  2 Problem Formulation

Suppose that the actual position of a vehicle at time $k$ is
$\mathbf{s}_{k}=[x_{k},y_{k}]^{T}$, the INS position estimate is
$\mathbf{s}_{k}^{INS}=[x_{k}^{INS},y_{k}^{INS}]^{T}$ and the position
corrected via map matching algorithm is
$\mathbf{s}_{k}^{C\text{-}INS}=[x_{k}^{C\text{-}INS},y_{k}^{C\text{-}INS}]^{T}$.
The gravity measurement at time $k$ and position $\mathbf{s}_{k}$ is modeled
as

$\displaystyle z_{k}=g(\mathbf{s}_{k})+\varepsilon_{k}$ (1)

where $g(\mathbf{s}_{k})$ describes the true Earth’s gravity at position
$\mathbf{s}_{k}$ and $\varepsilon_{k}\sim\mathcal{N}(0,\sigma_{z}^{2})$ is
noise with $0$ mean and standard deviation $\sigma_{z}$. Without loss of
generality, in the simulation, the true Earth gravity is obtained from the
gravity map without noise for the purpose of demonstrating the algorithm.
Additionally, we assume that the INS provides velocity measurements
$\mathbf{v}_{k}$ in the navigation frame and it is modeled as:

$\displaystyle\mathbf{v}_{k}=\mathbf{v}^{0}_{k}+\bm{\phi}_{k}$ (2)

where $\mathbf{v}_{k}^{0}$ is the ground truth, $\bm{\phi}_{k}$ is the noise
term comprising with a bias $\mathbf{b}_{k}$ and an independent zero mean
Gaussian noise with co-variance $\sigma_{v}^{2}\mathbf{I}_{2}$ and
$\mathbf{I}_{2}$ is the $2\times 2$ identity matrix.

We assume that the INS position estimates are corrected after every segment of
$T$ gravity measurements. Without loss of generality, we assume that the INS
position estimates and measurements are recorded from $k+1$. Then after $T$
sampling times, i.e., at time $k+T$, we have following available segments of
measurements and INS position estimates:

$\displaystyle\mathbf{Z}_{k}=\\{z_{k+1},\cdots,z_{k+T}\\},\;\mathbf{V}_{k}=\\{\mathbf{v}_{k+1},\cdots,\mathbf{v}_{k+T}\\},\;\mathbf{S}^{INS}_{k}=\\{\mathbf{s}_{k+1}^{INS},\cdots,\mathbf{s}_{k+T}^{INS}\\}$

The corrected position sequence using
$(\mathbf{Z}_{k},\mathbf{S}^{INS}_{k},\mathbf{V}_{k})$ and the proposed map
matching algorithm is denoted by
$\mathbf{S}_{k}^{C\text{-}INS}=\\{\mathbf{s}_{k+1}^{C\text{-}INS},\cdots,\mathbf{s}_{k+T}^{C\text{-}INS}\\}$.

In this work, since the map matching algorithm is implemented every $T$ time
steps, therefore we actually are interested in how the algorithm works within
time slot $\\{k+1,k+2,\cdots,k+T\\}$. For simplicity, where there is no
confusion, time stamp $k$ in this time slot is dropped from now on and we
write $t=1,\cdots,T$.

##  3 Viterbi Map Matching Algorithm

###  3.1 Constructing a grid on the gravity map

The gravity map can be viewed as a grid of pixels with each pixel representing
an area equivalent to the resolution of the map. These pixels will henceforth
be referred to as _cells_. The resolution of the gravity map, i.e. the
dimensions of the cell, is $R_{x}\times R_{y}$. Since each cell on the gravity
map can be associated with a location, we simply let the central position of
each cell be its location. Denote the collection of all cells’ locations on
the interested area of the map by $\bm{\mathcal{C}}$. We define $g_{m}(\cdot)$
as the “look-up function” for finding the value from the gravity map using the
given location, e.g. $g_{m}(\mathbf{s}_{t}^{INS})$ is the gravity value
associated to the cell that $\mathbf{s}_{t}^{INS}$ lies in.

In fact, at a specific sampling time $t$, we are only interested in a limited
area around the INS position estimate $\mathbf{s}_{t}^{INS}$ since the true
position should be within that region. Therefore, given an estimate
$\mathbf{s}_{t}^{INS}$, we can find a cell, or equivalently a pixel, on the
geo-referenced map of gravity such that $\mathbf{s}_{t}^{INS}$ lies in. It
should be noted that, since the cell is predefined and constructed according
to the pixel of the map, and the estimate $\mathbf{s}_{t}^{INS}$ is arbitrary,
therefore $\mathbf{s}_{t}^{INS}$ can be any point within the underlying cell
and is unnecessary to be the center position of this cell. Then a grid of
cells around the cell containing $\mathbf{s}_{t}^{INS}$ can be constructed and
denoted by $\mathbf{B}_{t}$. Let the size of the grid be $nR_{x}\times
nR_{y}$, where $n$ is a relatively small odd integer, and the positions of the
cells within this grid be denoted by $\mathbf{c}^{j}_{t}$ for
$j=1,\cdots,n^{2}$. We also use $\mathbf{c}^{j}_{t}$ as the labels of the
cells. The cells within the grid are arranged as shown in Fig. 1, where the
central cell of the grid, labelled as $\mathbf{c}^{\bar{n}}_{t}$ with
$\bar{n}=\frac{n^{2}+1}{2}$, is the one that $\mathbf{s}_{t}^{INS}$ lies in.
Its position is

$\displaystyle\mathbf{c}^{\bar{n}}_{t}\triangleq[x^{\bar{n}}_{t},y^{\bar{n}}_{t}]^{T}=\arg\min_{[x,y]^{T}\in\bm{\mathcal{C}}}\left\|[x,y]^{T}-\mathbf{s}_{t}^{INS}\right\|_{2}$
(3)

Essentially, (3) means that $\mathbf{c}^{\bar{n}}_{t}$ is the cell that the
INS position estimate $\mathbf{s}_{t}^{INS}$ lies in at time $t$. As a result,
each cell is associated to a gravity value according to the gravity map. Based
on $\mathbf{c}^{\bar{n}}_{t}$ and $n$, the positions of other cells in the
grid $\mathbf{B}_{k}$ are

$\displaystyle\mathbf{c}^{j}_{t}$
$\displaystyle=\mathbf{c}^{\bar{n}}_{t}-[m_{1}R_{x},m_{2}R_{y}]^{T},\;\;\text{where}\;\;m_{1},m_{2}=-n,\cdots,n\;\;\text{and}\;\;j=1,\cdots,n^{2}$
(4)

Figure 1: Illustrative example of a $3\times 3$ grid of cells, at arbitrary
time $t$, from a gravity map $g_{m}$. Each cell has a spatial location
$\mathbf{c}_{t}^{j}=[x_{t}^{j},x_{t}^{j}]^{T}$ and an associated gravity value
$g_{m}(\mathbf{c}_{t}^{j})$ at this location. The central cell of the grid,
$\mathbf{c}^{5}_{t}$, is the one that the INS position estimate
$\mathbf{s}_{t}^{INS}$ lies in and can be determined via (3), and the
positions of other cells can then be calculated by (4).

###  3.2 Viterbi based map-matching algorithm

At each time $t=1,2,\cdots,T$, we construct the grid $\mathbf{B}_{t}$ and its
associated cells $\mathbf{c}_{t}^{j_{t}}$, where the subscript $t$ is again
the time index and superscript $j_{t}$ is the cell label
$j_{t}=1,\cdots,n^{2}$. Then the problem becomes how to select a cell in each
$\mathbf{B}_{t}$ such that this segment of cells is the optimal estimate of
the actual trajectory. For a segment of cells,
$\\{\mathbf{c}^{j_{1}}_{1},\cdots,\mathbf{c}^{j_{T}}_{T}\\}$, with
corresponding referenced gravity values,
$\\{g_{m}(\mathbf{c}^{j_{1}}_{1}),\cdots,g_{m}(\mathbf{c}^{j_{T}}_{T})\\}$,
gravity measurements, $\mathbf{Z}_{k}$, INS information, $\mathbf{V}_{k}$ and
$\mathbf{S}^{INS}_{k}$, the joint posterior probability is

$\displaystyle
p\left(\\{\mathbf{c}^{j_{1}}_{1},\cdots,\mathbf{c}^{j_{T}}_{T}\\}\Big{|}\mathbf{Z}_{k},\mathbf{S}^{INS}_{k},\mathbf{V}_{k}\right)=p\left(\\{g_{m}(\mathbf{c}^{j_{1}}_{1}),\cdots,g_{m}(\mathbf{c}^{j_{T}}_{T})\\}\Big{|}\mathbf{Z}_{k},\mathbf{S}^{INS}_{k},\mathbf{V}_{k}\right)$
(5)

The most probable trajectory from $t=1$ to $t=T$ is given by the MAP estimator

$\displaystyle\\{j_{1}^{*},\cdots,j_{T}^{*}\\}=\arg\max_{{\bm{\mathcal{I}}^{\bigtimes}}}p\left(\\{\mathbf{c}^{j_{1}}_{1},\cdots,\mathbf{c}_{T}^{j_{T}}\\}\Big{|}\mathbf{Z}_{k},\mathbf{S}^{INS}_{k},\mathbf{V}_{k}\right)$
(6)

where $\bm{\mathcal{I}}=\\{1,\cdots,n^{2}\\}$,
$\bm{\mathcal{I}}^{\bigtimes}=\varprod_{i=1}^{T}\bm{\mathcal{I}}$, notation
$\varprod$ is Cartesian product and
$\\{\mathbf{c}^{j_{1}^{*}}_{1},\cdots,\mathbf{c}^{j_{T}^{*}}_{T}\\}$ is the
optimal path (locations) of the vehicle estimated from the given measurements.
It should be noted that $\mathbf{c}^{j_{t}}_{t}$ is actually the central
position of the each cells of the map; that is, the estimated path for time
interval ${1,\cdots,T}$ passes the center of each cell. Finally, the estimated
position sequence is then used to reset the corresponding segment of this
period estimated from INS.

The exhaustive searching process is required to find an optimal solution to
(6), where the computation load of the order of $n^{T}$ can explode quickly if
$n$ and/or $T$ getting large. In this work, we intend to find an efficient
algorithm to solve this problem.

Note that, for arbitrary sequence, $\\{j_{1},\cdots,j_{T}\\}$

$\displaystyle
p\left(\mathbf{c}_{t}^{j_{t}}\Big{|}\mathbf{c}_{t-1}^{j_{t-1}},\cdots,\mathbf{c}_{1}^{j_{1}}\right)=p\left(\mathbf{c}_{t}^{j_{t}}\Big{|}\mathbf{c}_{t-1}^{j_{t-1}}\right),\;\;\;\forall
t\in\\{2,\cdots,T\\}\;\;\text{and}\;\;j_{t}\in\bm{\mathcal{I}},$ (7)

which implies that the selection of $\mathbf{c}_{t}^{j_{t}}$ is only dependent
on that of $\mathbf{c}_{t-1}^{j_{t-1}}$. This is a typical Markov process
(Stroock, 2013). The sequence
$\\{\mathbf{c}_{1}^{j_{1}},\cdots,\mathbf{c}_{T}^{j_{T}}\\}$ with the
measurements can be modelled by an HMM. The Viterbi algorithm (Viterbi, 1967)
is available to compute (6) and provides a most likely sequence of positions.

At each time $t$, the likelihood function of $z_{t}$ given
$\mathbf{c}_{t}^{j_{t}}$ and $g_{m}(\mathbf{c}_{t}^{j_{t}})$ is denoted by

$\displaystyle
p_{z}\left(z_{t}\Big{|}\mathbf{c}_{t}^{j_{t}}\right)=p_{z}\left(z_{t}\Big{|}g_{m}(\mathbf{c}_{t}^{j_{t}})\right)$
(8)

where $p_{z}(\cdot)$ is the probability density function of the gravity
measurement. For each successive pair $\\{t-1,t\\}$ with $t=2,\cdots,T$, the
measurement likelihood function of $\mathbf{c}^{j_{t}}_{t}$ given
$\mathbf{v}_{t-1}$, the velocity measurement at time ${t-1}$, and
$\mathbf{c}_{t}^{j_{t}}$, is

$\displaystyle
p_{v}\left(\mathbf{c}^{j_{t}}_{t}\Big{|}\mathbf{v}_{t-1},\mathbf{c}_{t-1}^{j_{t-1}}\right)$
(9)

where $p_{v}(\cdot)$ is the probability density function of the velocity
measurement.

Then, given $\mathbf{c}^{j_{t-1}}_{t-1}$ and measurement $z_{t}$, the
posterior density of $\mathbf{c}^{j_{t}}_{t}$ is

$\displaystyle
p\left(\mathbf{c}^{j_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{t-1},\mathbf{c}^{j_{t-1}}_{t-1}\right)\propto$
$\displaystyle
p\left(z_{t}\Big{|}\mathbf{c}_{t}^{j_{t}},\mathbf{v}_{t-1},\mathbf{c}^{j_{t-1}}_{t-1}\right)p\left(\mathbf{c}^{j_{t}}_{t}\Big{|}\mathbf{v}_{t-1},\mathbf{c}^{j_{t-1}}_{t-1}\right)$
$\displaystyle=$ $\displaystyle
p_{z}\left(z_{t}\Big{|}g_{m}(\mathbf{c}^{j_{t}}_{t})\right)p_{v}\left(\mathbf{c}^{j_{t}}_{t}\Big{|}\mathbf{v}_{t-1},\mathbf{c}_{t-1}^{j_{t-1}}\right)$
(10)

We shall be considering this, in particular, in the case when the transition
is a single time step and $j_{1}=j_{t-1}$ and $j_{2}=j_{t}$, for a pair
$\\{j_{1},j_{2}\\}\in\bm{\mathcal{I}}\varprod\bm{\mathcal{I}}$, we use the
notation

$\displaystyle\left\\{j_{1}\rightarrow
j_{2},p\left(\mathbf{c}^{j_{2}}_{2}\Big{|}z_{2},\mathbf{v}_{1},\mathbf{c}^{j_{1}}_{1}\right)\right\\}$
(11)

for the potential path from $\mathbf{c}^{j_{1}}_{1}$ to
$\mathbf{c}^{j_{2}}_{2}$ and its associated conditional probability.

Similarly, for any pair
$\\{j_{t-1},j_{t}\\}\in\bm{\mathcal{I}}\varprod\bm{\mathcal{I}}$, we can have
such a path with its posterior probability as

$\displaystyle\left\\{j_{t-1}\rightarrow
j_{t},p\left(\mathbf{c}^{j_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{t-1},\mathbf{c}^{j_{t-1}}_{t-1}\right)\right\\}$
(12)

Then, starting from $t=2$, for a fixed $j_{1}$, a local optimal path with
greatest posterior probability is given by

$\displaystyle\left\\{j_{1}\rightarrow
j^{*}_{2}\Big{|}\;j^{*}_{2}=\arg\max_{j_{2}\in\bm{\mathcal{I}}}p\left(\mathbf{c}^{j_{2}}_{2}\Big{|}z_{2},\mathbf{v}_{1},\mathbf{c}^{j_{1}}_{1}\right)\right\\}$
(13)

Similarly, at $t=3$, for pair $\\{j^{*}_{2},j_{3}\\}$, where $j^{*}_{2}$ is
obtained from (13) and therefore fixed and $j_{3}\in\bm{\mathcal{I}}$, the
following local optimal path is

$\displaystyle\left\\{j^{*}_{2}\rightarrow
j^{*}_{3}\Big{|}\;j^{*}_{3}=\arg\max_{j_{3}\in\bm{\mathcal{I}}}p\left(\mathbf{c}^{j_{3}}_{3}\Big{|}z_{3},\mathbf{v}_{2},\mathbf{c}^{j^{*}_{2}}_{2}\right)\right\\}$
(14)

Iteratively, for any pair $\\{j^{*}_{t-1},j_{t}\\}$,
$j_{t}\in\bm{\mathcal{I}}$, the following optimization is done:

$\displaystyle\left\\{j^{*}_{t-1}\rightarrow
j^{*}_{t}\Big{|}\;j^{*}_{t}=\arg\max_{j_{t}\in\bm{\mathcal{I}}}p\left(\mathbf{c}^{j_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{t-1},\mathbf{c}^{j^{*}_{t-1}}_{t-1}\right)\right\\}$
(15)

As a result, for any fixed $j_{1}\in\bm{\mathcal{I}}$, there is a candidate
path that the vehicle potentially travels:

$\displaystyle\left\\{j_{1}\rightarrow j^{*}_{2}\rightarrow\cdots\rightarrow
j^{*}_{t}\rightarrow\cdots\rightarrow j^{*}_{T-1}\rightarrow
j^{*}_{T}\right\\}\triangleq\bm{\mathcal{J}}^{*}_{j_{1}}$ (16)

where

$\displaystyle
j^{*}_{t}=\begin{cases}\displaystyle\arg\max_{j_{t}\in\bm{\mathcal{I}}}p\left(\mathbf{c}^{j_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{1},\mathbf{c}^{j_{1}}_{1}\right)&t=2\\\
\displaystyle\arg\max_{j_{t}\in\bm{\mathcal{I}}}p\left(\mathbf{c}^{j_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{t-1},\mathbf{c}^{j^{*}_{t-1}}_{t-1}\right)&t=3,\cdots,T\end{cases}$
(17)

Since $\bm{\mathcal{J}}^{*}_{j_{1}}$ depends on the selection of $j_{1}$, we
are able to collect the potential paths corresponding to all possible values
of $j_{1}$, i.e.
$\bm{\mathcal{J}}^{*-}=\\{\bm{\mathcal{J}}^{*}_{1},\cdots,\bm{\mathcal{J}}^{*}_{n^{2}}\\}$.
Then,

$\displaystyle j^{*}_{1}$
$\displaystyle=\arg\max_{\bm{\mathcal{J}}^{*-}}p_{z}\left(z_{1}\Big{|}g_{m}(\mathbf{c}_{1}^{j_{1}})\right)p\left(\mathbf{c}^{j^{*}_{2}}_{2}\Big{|}z_{2},\mathbf{v}_{1},\mathbf{c}^{j_{1}}_{1}\right)\prod_{t=3}^{T}p\left(\mathbf{c}^{j^{*}_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{t-1},\mathbf{c}^{j^{*}_{t-1}}_{t-1}\right)$
(18)

and the complete estimated path is

$\displaystyle\left\\{j^{*}_{1}\rightarrow
j^{*}_{2}\rightarrow\cdots\rightarrow j^{*}_{t}\rightarrow\cdots\rightarrow
j^{*}_{T-1}\rightarrow j^{*}_{T}\right\\}\triangleq\bm{\mathcal{J}}^{*}$ (19)

Eventually, the estimated path $\bm{\mathcal{J}}^{*}$ is used to reset the INS
estimated position directly. An example of how the algorithm works is shown in
Fig. 2.

Figure 2: An example of the Viterbi map matching algorithm. The time sequence
is $t=1,2,3,4$ and the dimension of each grid is $5\times 5$. The blue lines
are the INS estimated positions while the red ones are the actual positions.
The black dashed lines indicate two candidate trajectories satisfying (16),
i.e. $\bm{\mathcal{J}}^{*}_{1}=\\{1\rightarrow 1\rightarrow 1\rightarrow 1\\}$
and $\bm{\mathcal{J}}^{*}_{25}=\\{25\rightarrow 25\rightarrow 25\rightarrow
25\\}$. All trajectories shown in this figure are for demonstration only.

In practice, given a measurement $z_{t}$ at a specific time $t$, we may find
more than one cell within a grid $\mathbf{B}_{t}$ that has the same likelihood
value using (8). Therefore multiple optimal trajectories may be identified
because of that they can have the same maximum posterior value. If this
happens we choose the trajectory that is closest to the INS estimated position
for simplicity.

###  3.3 An enhanced two-layer algorithm

In Section 3.2, a grid-based Viterbi algorithm is proposed to find a path used
to correct $\mathbf{S}^{INS}_{k}$. However, it should be noted that the
proposed algorithm only returns a path that passes through the center of each
cell. As a result, the deviation between this path and the actual trajectory
is highly dependent on the resolution of the gravity map since the actual
trajectory can be any point of each cell. In this section, an enhanced method
is presented to improve the accuracy of the estimated path.

The key idea of this enhanced algorithm is to generate refined sub-cells for
each cell and use these refined sub-cells to match the actual trajectory. This
strategy takes effect when constructing the grid $\mathbf{B}_{t}$. An example
of a refined cell is given in Fig. 3.

Figure 3: Illustration of the refined grid. The cell $\mathbf{c}_{t}^{5}$ is
an example to show how the refined sub-cells are constructed, where $n=3$ and
$o=3$.

Suppose that each cell is divided into $o\times o$ sub-cells, the centers of
which, for $\mathbf{c}_{t}^{j_{t}}$, are denoted by
$\mathbf{c}_{t}^{j_{t},l_{t}}$ where $l_{t}\in\bm{\mathcal{O}}$ and
$\bm{\mathcal{O}}=\\{1,2,\cdots,o^{2}\\}$. We also use
$\mathbf{c}_{t}^{j_{t},l_{t}}$ as label of $l_{t}$-th sub-cell of $j_{t}$-th
cell at time $t$. Since the resolution of gravity map is fixed, therefore the
sub-cells $\mathbf{c}_{t}^{j_{t},l_{t}}$ refine cells $\mathbf{c}_{t}^{j_{t}}$
to provide better resolution on positions and do not improve the actual
resolution of the gravity map. In the other word, the cell
$\mathbf{c}_{t}^{j_{t}}$ within $\mathbf{B}_{t}$ is divided into multiple sub-
cells $\mathbf{c}_{t}^{j_{t},l_{t}}$. Then we have the following relation, for
any given $j_{t}$:

$\displaystyle
g_{m}(\mathbf{c}_{t}^{j_{t}})=g_{m}(\mathbf{c}_{t}^{j_{t},l_{t}}),\;\;\forall
l_{t}\in\bm{\mathcal{O}}$ (20)

and

$\displaystyle
p_{z}\left(z_{t}\Big{|}g_{m}(\mathbf{c}_{t}^{j_{t},l_{t}})\right)=p_{z}\left(z_{t}\Big{|}g_{m}(\mathbf{c}_{t}^{j_{t}})\right),\;\;\forall
l_{t}\in\bm{\mathcal{O}}$ (21)

In similar vein to the (15), (16) and (17), the potential path for a given
$\\{j_{1},l_{1}\\}$ is

$\displaystyle\left\\{\\{j_{1},l_{1}\\}\rightarrow\\{j^{*}_{2},l^{*}_{2}\\}\rightarrow\cdots\rightarrow\\{j^{*}_{T-1},l^{*}_{T-1}\\}\rightarrow\\{j^{*}_{T},l^{*}_{T}\\}\right\\}\triangleq\bm{\mathcal{L}}_{j_{1},l_{1}}^{*}$
(22)

where $l^{*}_{t}$ in $\\{j^{*}_{t},l^{*}_{t}\\}$ is the optimal sub-cell and
$j^{*}_{t}$ is the associated cell. $\\{j^{*}_{t},l^{*}_{t}\\}$ is obtained
using $\\{j^{*}_{t-1},l^{*}_{t-1}\\}$ and it can be calculated from (17) via
replacing $j_{t}$ by $\\{j_{t},l_{t}\\}$.

Let
$\bm{\mathcal{L}}^{*-}=\\{\bm{\mathcal{L}}^{*}_{1,1},\bm{\mathcal{L}}^{*}_{1,2}\cdots,\bm{\mathcal{L}}^{*}_{n^{2},o^{2}}\\}$.
Then $\\{j^{*}_{1},l^{*}_{1}\\}$ in the estimated path

$\displaystyle\bm{\mathcal{L}}^{*}=\left\\{\\{j^{*}_{1},l^{*}_{1}\\}\rightarrow\\{j^{*}_{2},l^{*}_{2}\\}\rightarrow\cdots\rightarrow\\{j^{*}_{T-1},l^{*}_{T-1}\\}\rightarrow\\{j^{*}_{T},l^{*}_{T}\\}\right\\}$
(23)

is obtained by

$\displaystyle\\{j^{*}_{1},l^{*}_{1}\\}=\arg\max_{\bm{\mathcal{L}}^{*-}}p_{z}\left(z_{1}\Big{|}g_{m}(\mathbf{c}_{1}^{j_{1},l_{1})}\right)p\left(\mathbf{c}^{j^{*}_{2},l^{*}_{2}}_{2}\Big{|}z_{2},\mathbf{v}_{1},\mathbf{c}^{j_{1},l_{1}}_{1}\right)\prod_{t=3}^{T}p\left(\mathbf{c}^{j^{*}_{t},l^{*}_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{t-1},\mathbf{c}^{j^{*}_{t-1},l^{*}_{t-1}}_{t-1}\right)$
(24)

###  3.4 Reducing the Computation

In the algorithms presented above, the computational overhead increases
dramatically as more cells are considered in the algorithm. To mitigate this
issue, we constrain the number of $j_{t}$ in each time $t$ via a threshold.
Relative likelihood values to select the available $j_{t}$ are calculated as
follows:

$\displaystyle\bm{\mathcal{A}}_{t}=\left\\{j_{t}\left|\frac{p_{z}\left(z_{t}\Big{|}g_{m}(\mathbf{c}_{t}^{j_{t}})\right)}{\max_{j_{t}\in\bm{\mathcal{I}}}p_{z}\left(z_{t}\Big{|}g_{m}(\mathbf{c}_{t}^{j_{t}})\right)}\geq\alpha,\forall
j_{t}\in\bm{\mathcal{I}}\right.\right\\}$ (25)

where $\alpha\in[0,1]$. (25) implies that only those cells with relatively
high likelihood values are chosen. In other words, cells that are very
unlikely to form part of the path are discarded.

An example of higher likelihood cells from $\mathbf{c}_{t}^{j_{t}}$,
$j_{t}\in\bm{\mathcal{I}}$, is given in Fig. 4, where we assume that $T=3$ and
the size of the grid at each time stamp is $n=11$. The values of
$p_{z}\left(z_{t}\Big{|}g_{m}(\mathbf{c}_{1}^{j_{t}})\right)$, $t=1,2,3$, are
represented in different colors. At time $1$, there are three higher
likelihood cells. Similarly, at time 2 and 3, there are $3$ and $2$ grids
available respectively. As a result, $\bm{\mathcal{A}}_{t}$, $t=1,2,3$, is
obtained accordingly, e.g. $\bm{\mathcal{A}}_{1}=\\{8,18,29,40,61,83,94\\}$.
The computational complexity is reduced from $O((T-1)n^{2})$ to
$O(\sum_{t=1}^{T-1}\\#\bm{\mathcal{A}}_{t}\times\\#\bm{\mathcal{A}}_{t+1})$,
where $\\#\bm{\mathcal{A}}_{t}$ is the number of elements in
$\bm{\mathcal{A}}_{t}$.

Figure 4: Illustration of likelihood values of cells within the grids at
successive time stamp $t=1,2,3$. The windows size is $n=11\times 11$. For
arbitrary cell $j_{t}$ at time $t$, the likelihood values are calculated from
$p_{z}(z_{t}|g_{m}(\mathbf{c}_{t}^{j_{t}}))$ and normalized via
${p_{z}(z_{t}|g_{m}(\mathbf{c}_{t}^{j_{t}}))}\big{/}{\max_{j^{\prime}_{t}\in\bm{\mathcal{I}}}\left\\{p_{z}(z_{t}|g_{m}(\mathbf{c}_{t}^{j^{\prime}_{t}}))\right\\}}$.
There is no sub-cell for each cell to simplify the illustration. In each grid,
there are only up to $7$ cells having relative high likelihood values while
other cells’ values are close to $0$. If (25) is applied and let $\alpha=0.1$,
the maximum number of possible trajectories is $7^{3}$; otherwise, in the
algorithm without applying (25), the maximum number of possible trajectories
is $121^{3}$.

In practice, to simplify the problem, $\alpha$ can be set to be a fixed value,
such as $0.1$ or $0.2$.

Similarly, we use indexes in
${\bm{\mathcal{A}}_{1},\cdots,\bm{\mathcal{A}}_{T}}$ to find multiple
candidate paths using the same method as (22) by replacing the associated
indexes sequence. i.e.

$\displaystyle\left\\{\\{j_{1},l_{1}\\}\rightarrow\\{j^{*}_{2},l^{*}_{2}\\}\rightarrow\cdots\rightarrow\\{j^{*}_{T-1},l^{*}_{T-1}\\}\rightarrow\\{j^{*}_{T},l^{*}_{T}\\}\right\\}\triangleq\bm{\mathcal{R}}_{j_{1},l_{1}}^{*}$
(26)

where $j_{1}\in\bm{\mathcal{A}}_{1}$ and $j^{*}_{t}\in\bm{\mathcal{A}}_{t}$
for $t=2,\cdots,T$.

Let
$\bm{\mathcal{R}}^{*-}=\left\\{\bm{\mathcal{R}}_{j_{1},l_{1}}^{*}\big{|}\;\forall
j_{1}\in\bm{\mathcal{A}}_{1}\right\\}$ contain these candidate paths. Then
(24) can be rewritten as, where $j_{1}\in\bm{\mathcal{A}}_{1}$,

$\displaystyle\\{j^{*}_{1},l^{*}_{1}\\}=\arg\max_{\bm{\mathcal{R}}^{*-}}p_{z}\left(z_{1}\Big{|}g_{m}(\mathbf{c}_{1}^{j_{1},l_{1}})\right)p\left(\mathbf{c}^{j^{*}_{2},l^{*}_{2}}_{2}\Big{|}z_{2},\mathbf{v}_{1},\mathbf{c}^{j_{1},l_{1}}_{1}\right)\prod_{t=3}^{T}p\left(\mathbf{c}^{j^{*}_{t},l^{*}_{t}}_{t}\Big{|}z_{t},\mathbf{v}_{t-1},\mathbf{c}^{j^{*}_{t-1},l^{*}_{t-1}}_{t-1}\right)$
(27)

The estimated path is, then, as follows:

$\displaystyle\bm{\mathcal{R}}^{*}=\left\\{\\{j^{*}_{1},l^{*}_{1}\\}\rightarrow\\{j^{*}_{2},l^{*}_{2}\\}\rightarrow\cdots\rightarrow\\{j^{*}_{T-1},l^{*}_{T-1}\\}\rightarrow\\{j^{*}_{T},l^{*}_{T}\\}\right\\},\;\;\text{where}\;\;j_{t}\in\bm{\mathcal{A}}_{t}$
(28)

In summary, the algorithm proposed in Sections 3.2 with enhanced strategies
presented in 3.3 and 3.4 is given in Algorithm 1.

Result: The estimated trajectory

1 Initialization $k\leftarrow 1$

2while  _ $k<\text{Running time}$_ do

3 Obtain position via INS

4 Push measurements and INS estimated position into vector $\mathbf{Z}_{k}$,
$\mathbf{S}^{INS}_{k}$ and $\mathbf{V}_{k}$ respectively.

5 if  _$k\mod{T}==0$_ then

6 Construct grids $\mathbf{B}_{t}$ using $n$, $o$, (3) and sub-cell method
proposed in Section 3.3 for each $t=1,\cdots,T$

7 Find all cells $\bm{\mathcal{A}}_{t}$ satisfying (25) for each
$t=1,\cdots,T$

8 for  _$j_{1}\in\bm{\mathcal{A}}_{1}$ _ do

9 for  _$l_{1}\leftarrow 1$ to $o$ _ do

10 for  _$t\leftarrow 1$ to $T$ _ do

11 Find ordered sequence $\bm{\mathcal{R}}_{j_{1}}^{*}$ according to (26)

12 end for

13 Push $\bm{\mathcal{R}}_{j_{1}}^{*}$ into $\bm{\mathcal{R}}^{*-}$

14 end for

15

16 end for

17 Find the estimated path $\bm{\mathcal{R}}^{*}$ using
$\bm{\mathcal{R}}^{*-}$ via (28)

18 Replace the INS estimated trajectory using $\bm{\mathcal{R}}^{*}$

19 Clear $\mathbf{Z}_{k}$, $\mathbf{S}^{INS}_{k}$ and $\mathbf{V}_{k}$

20 end if

21 $k\leftarrow k+1$

22 end while

Algorithm 1 The Viterbi-based map matching algorithm

##  4 Simulation

In the first part of the simulation, we assume that a vehicle (aircraft)
travels across Australia, approximately from Melbourne to Sydney, at a
constant velocity $7.54^{\circ}$/h, roughly $837$ Km/h ($1^{\circ}\approx
111\text{Km}$ (Scales et al., 2017)). The total flight distance is about $814$
Km. The gravity map of the region of interest has been obtained from the
GGMplus software (Hirt et al., 2013) with a variety of resolutions, where the
component taken from the map is gravitational acceleration. The map is
represented in grids that have specified resolutions with the grid equally
spaced in terms of geodetic (GRS80) latitude and longitude. The map with lower
resolution is created by down-sampling the higher one. An illustrative example
is shown in Fig.5a

(a)

(b)

Figure 5: Illustration of the scenarios used in the simulation. The colored
map is the reference gravitational acceleration map of Victoria and New south
Wales (5a); and Western Australia (5b). The INS estimated trajectory without
aiding deviates gradually from the true trajectory. and its resolution is
$\frac{1}{200}^{\circ}$.

To simplify the simulation, the noise in the velocity measurement, i.e. (2),
from the IMU is modelled as
$\bm{\phi}_{k}=\mathbf{b}_{k}+\mathcal{N}(0,\sigma_{v}^{2}\mathbf{I}_{2})$,
where $\mathbf{b}_{k}$ is bias vector. The measured gravitational acceleration
is assumed to be corrupted by additive Gaussian noise
$\mathcal{N}(0,\sigma_{z}^{2})$, where $\sigma_{z}$ is the standard deviation.
As mentioned earlier, the true Earth gravity is obtained from gravity map
without noise for the purpose of demonstrating the algorithm.

A Monte Carlo (MC) simulation of 500 runs is carried out for each scenario and
simulation performance at time $k$ is measured by the root-mean-square error,
defined as

$\displaystyle\text{Error}_{k}=\frac{1}{nMC}\sum_{m=1}^{nMT}\text{Dist}\left(\mathbf{s}_{k,m},\mathbf{s}_{k,m}^{C\text{-}INS}\right)$
(29)

where $nMC$ is the number of MC, $\mathbf{s}_{k,m}$ and
$\mathbf{s}_{k,m}^{C\text{-}INS}$ are the position and updated position
respectively of the target at time $k$ and $m$-th simulation and
$\text{Dist}\left(\cdot,\cdot\right)$ is the Haversine distance (Km) between
two points given their longitude and latitude (Gade, 2010).

The algorithms are evaluated by different parameters. In the simulation, the
bias of the velocity is set to be $1^{\circ}/hr$ and $\sigma_{v}={9\times
10^{-6}}^{\circ}$/s, roughly $1$ m/s. In the simulation, since we assume that
the inertial sensors used in the vehicle INS exhibits white noise and bias,
such that the INS computed position deviates from ground truth at the end of
journey as shown in Fig. 5. The standard deviation of the gravitational
acceleration measurement is $\sigma_{z}=1,2$ mGal. The INS reports positions
and velocities every $12$ s.

Firstly, the performances of the Viterbi based map matching algorithm (VBMP)
proposed in Section 3.2 and the computation reduced algorithm (RVBMP) proposed
in Section 3.4 are compared in terms of the time efficiency. In the
simulation, standard deviation $\sigma_{z}=1$ mGal, the map resolution is
$1/100^{\circ}$, the correction rate, i.e. $T$, is set to be $6$, window size
is $n=13$ and the sub-cell algorithm proposed in 3.3 is not applied. It also
can be noticed that the value of $\alpha$ can potentially have impact on the
performance of RVBMP; that is, the larger the $\alpha$ is, the more cells will
be included in $\bm{\mathcal{A}}_{t}$ and more time will be required to find a
solution. In considering this, the performance of RVBMP using different values
of $\alpha$ are compared with VBMP. When $\alpha=0$, set
$\bm{\mathcal{A}}_{t}$ will contain all cells and the algorithm becomes VBMP.
The simulation results are given in Table 1, where TCR stands for time
consumption ratio between the time consumption using a specific $\alpha$ and a
referenced time, i.e. $\frac{\text{Time}_{\alpha}}{\text{Time}_{\text{ref}}}$,
and $\text{Time}_{\alpha=0.4}$ is chosen as the reference in Table 1, Mean and
Std. Dev. are sample mean and sample standard deviation of $\text{Error}_{k}$,
i.e.

$\displaystyle\text{Mean}=\frac{1}{L}\sum_{k=1}^{L}\text{Error}_{k},\;\;\;\text{Std.
Dev.}=\sqrt{\frac{1}{L-1}\sum_{k=1}^{L}\left(\text{Error}_{k}-\text{Mean}\right)}$

where $L$ is the total number of sampling points over the trajectory.

Table 1: Performance comparison between RVBMP and RVBMP using different values
of $\alpha$. The scenario is shown in Fig. 5a. VBMP and RVBMP
---
$\alpha$ | Mean (Km) | Std. Dev. (Km) | TCR
$0$ (VBMP) | $1.2802$ | $0.5011$ | $35.3817$
$0.05$ | $1.2890$ | $0.5101$ | $2.6162$
$0.1$ | $1.2786$ | $0.5219$ | $2.0119$
$0.2$ | $1.5114$ | $0.5470$ | $1.5049$
$0.3$ | $1.7097$ | $0.5655$ | $1.1986$
$0.4$ | $2.3764$ | $0.8590$ | $1$

It can be seen that the RVBMP performs better than VBMP in terms of time
efficiency and they have the similar tracking performance when $\alpha$ is
less than $0.1$.

Secondly, the performance of RVBMP and the enhanced two-layer algorithm
(RVBMP-2) in Section 3.3 are compared via standard deviation $\sigma_{z}=1$
mGal and two map resolutions. The correction rate, i.e. $T$, is set to be $6$,
window size is $n=13$ and the size of a sub-cell is $o=7$ for RVBMP-2. The
results are given in Fig. 6 and Table 2. The definitions of Mean and Std. Dev
are same as those in Table 1 and they are calculated over those estimated
results that are not divergent. The time consumption ratio (TCR) is given to
compare the time efficiency of these two algorithms. The TCR in this table is
defined as: $\frac{\text{Time}_{RVBMP-2}}{\text{Time}_{RVBMP}}$.

From Fig. 6 and Table 2, we can see that the RVBMP-2 outperforms RVBMP in both
resolutions, which indicates that the sub-cell can improve the RVBMP. It
should be noticed that the sub-cell strategy works better when the map
resolution is relative low. This is because of that the high map resolution
implies that high position resolution as well so that the sub-cell strategy
has less influence on the performance.

(a)

(b)

Figure 6: Simulation results. The scenario is shown in Fig. 5a. LABEL:sub@sim00:1 The map resolution is $1/100^{\circ}$; LABEL:sub@sim00:2 The map resolution is $1/200^{\circ}$. Table 2: Performance comparison between RVBMP and RVBMP-2 using two different map resolutions. The scenario is shown in Fig. 5a. Map resolution = $1/100^{\circ}$ |
---|---
Algorithm | Mean (Km) | Std. Dev. (Km) | TCR
RVBMP | $1.4358$ | $0.5715$ | $1$
RVBMP-2 | $0.9190$ | $0.4004$ | $24.9$
Map resolution = $1/200^{\circ}$ |
Algorithm | Mean (Km) | Std. Dev. (Km) | TCR
RVBMP | $0.5846$ | $0.2639$ | $1$
RVBMP-2 | $0.5266$ | $0.2538$ | $28.5$

Finally, to further demonstrate the efficiency of RVBMP-2, the performance of
RVBMP-2 and ICCP (Kamgar-Parsi & Kamgar-Parsi, 1999; Han et al., 2017) are
compared using two different noise levels, i.e. $\sigma_{z}$. The map
resolution is $1/200^{\circ}$, window size is $n=13$ and the size of sub-cell
is $o=5$. The simulation results are shown in Fig. 7-8 and summarized in Table
3, where the success rate represents the percentage of times the algorithm
followed the true trajectory without divergence, and the definitions of Mean
and Std. Dev. are as in Table 2.

(a)

(b)

Figure 7: Simulation results of ICCP and proposed algorithm when the
correction rate is $4$. The scenario is shown in Fig. 5a. LABEL:sub@sim1:1
$\sigma_{z}=1\;\text{mGal}$; LABEL:sub@sim1:2 $\sigma_{z}=2\;\text{mGal}$.

(a)

(b)

Figure 8: Simulation results of ICCP and proposed algorithm when the
correction rate is $6$. The scenario is shown in Fig. 5a. LABEL:sub@sim2:1
$\sigma_{z}=1\;\text{mGal}$; LABEL:sub@sim2:2 $\sigma_{z}=2\;\text{mGal}$.
Table 3: Performance comparison between RVBMP-2 and ICCP using two different
standard deviations of the noise, $\sigma_{z}$ and correction rate. The
scenario is shown in Fig. 5a. Correction rate = $4$
---
$\sigma_{z}$ (mGal) | Algorithm | Success rate | Mean (Km) | Std. Dev. (Km)
1 | ICCP | $92\%$ | $1.3424$ | $0.6642$
RVBMP-2 | $100\%$ | $0.8300$ | $0.3557$
2 | ICCP | $54\%$ | $2.0253$ | $1.1920$
RVBMP-2 | $100\%$ | $1.4238$ | $0.4146$
Correction rate = $6$
$\sigma_{z}$ (mGal) | Algorithm | Success rate | Mean (Km) | Std. Dev. (Km)
1 | ICCP | $75\%$ | $0.8675$ | $0.4297$
RVBMP-2 | $100\%$ | $0.5765$ | $0.2722$
2 | ICCP | $46\%$ | $1.1533$ | $0.5148$
RVBMP-2 | $100\%$ | $0.8478$ | $0.3466$

From Table 3 and Fig. 7-8, we see that RVBMP-2 algorithm outperforms ICCP both
in success rate and tracking performance.

In the second part of the simulation, a reference gravity map that covers a
major area of Western Australia (WA), Australia, see Fig. 5b, is used to
further demonstrate the performance of the proposed algorithm. The major
challenge of this scenario is that the gravitational acceleration is
relatively flat; that is, the variation in gravity across distinct
pixels/cells is, on the whole, smaller than that in the regions in the earlier
simulations. This flatness of the gravity map suggests that aiding with it
will perform less well than in the previous results.

In the simulation, the start point of the vehicle is near Perth, and the end
point is located at ($[-24.427111,124.978610]$). The distance is about $1206$
Km. The remaining parameters are similar to those in Table 3.

Table 4: Performance comparison between RVBMP-2 and ICCP using two different
standard deviations of the noise, $\sigma_{z}$ and the correction rate. The
scenario is shown in Fig. 5b. Correction rate = $4$
---
$\sigma_{z}$ (mGal) | Algorithm | Success rate | Mean (Km) | Std. Dev. (Km)
1 | ICCP | $23\%$ | $1.5880$ | $0.8657$
RVBMP-2 | $100\%$ | $1.4861$ | $0.6183$
2 | ICCP | $2\%$ | $2.4581$ | $1.0113$
RVBMP-2 | $100\%$ | $2.3573$ | $0.8325$
Correction rate = $6$
$\sigma_{z}$ (mGal) | Algorithm | Success rate | Mean (Km) | Std. Dev. (Km)
1 | ICCP | $22\%$ | $1.0967$ | $0.4457$
RVBMP-2 | $100\%$ | $1.0679$ | $0.4336$
2 | ICCP | $2\%$ | $1.6697$ | $0.7551$
RVBMP-2 | $99\%$ | $1.6003$ | $0.5975$

In the Table 4, we can see that, in terms of accuracy, the performances of
RVBMP-2 and ICCP are similar because of the map is relatively flat and
provides less information to the algorithm. However, the former has better a
success rate than the latter.

##  5 Summary

In this paper, we have presented a Viterbi-inspired map matching algorithm
based on a hidden Markov model (HMM) formalisation of the map matching
problem. Our algorithm determines the optimal trajectory of a platform given a
sequence of gravitational measurements, a geo-reference gravity map and INS
velocity information. This trajectory is optimal in the maximum a posterior
sense and is used to correct the platform’s INS. The simulation results
demonstrate that the algorithm is efficient in correcting the INS estimated
trajectory. In addition, we proposed an enhanced 2-layer addition to the
algorithm to improve the performance when the map resolution is low. We showed
by simulations that our proposed algorithm outperforms existing map matching
algorithms in terms of navigational accuracy and is also robust to varying
levels of sensor noise. In terms of future work, we have validated our
algorithm using gravitational measurements however it can be potentially
generalized to other aiding resources, such as terrain elevation and magnetic
field, that have similar maps. Additionally, we currently use a
straightforward method when combining the output of the algorithm with the INS
position – the INS position is simply reset with the corrected position.
Accordingly, another area of future work is to develop a more sophisticated
fusion process.

##  6 Conflict of interest

None

## References

  * (1)
  * Affleck & Jircitano (1990) Affleck, C. A. & Jircitano, A. (1990), Passive gravity gradiometer navigation system, in ‘IEEE Symposium on Position Location and Navigation’, pp. 60–66.
  * Bergman (1997) Bergman, N. (1997), ‘A bayesian approach to terrain-aided navigation’, IFAC Proceedings Volumes 30(11), 1457–1462.
  * Gade (2010) Gade, K. (2010), ‘A non-singular horizontal position representation’, The journal of navigation 63(3), 395–417.
  * Groves (2013) Groves, P. D. (2013), Principles of GNSS, inertial, and multisensor integrated navigation systems, Artech House.
  * Han et al. (2016) Han, Y., Wang, B., Deng, Z. & Fu, M. (2016), ‘An improved TERCOM-based algorithm for gravity-aided navigation’, IEEE Sensors Journal 16(8), 2537–2544.
  * Han et al. (2017) Han, Y., Wang, B., Deng, Z. & Fu, M. (2017), ‘A combined matching algorithm for underwater gravity-aided navigation’, IEEE/ASME Transactions on Mechatronics 23(1), 233–241.
  * Hemann et al. (2016) Hemann, G., Singh, S. & Kaess, M. (2016), Long-range GPS-denied aerial inertial navigation with LIDAR localization, in ‘2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)’, IEEE, pp. 1659–1666.
  * Hirt et al. (2013) Hirt, C., Claessens, S., Fecher, T., Kuhn, M., Pail, R. & Rexer, M. (2013), ‘New ultrahigh-resolution picture of Earth’s gravity field’, Geophysical research letters 40(16), 4279–4283.
  * Hostetler & Andreas (1983) Hostetler, L. & Andreas, R. (1983), ‘Nonlinear Kalman filtering techniques for terrain-aided navigation’, IEEE Transactions on Automatic Control 28(3), 315–323.
  * Jircitano & Dosch (1991) Jircitano, A. & Dosch, D. (1991), Gravity aided inertial navigation system (GAINS), in ‘Proc. 47th Annual Meeting of The Institute of Navigation’, pp. 221–229.
  * Kamgar-Parsi & Kamgar-Parsi (1999) Kamgar-Parsi, B. & Kamgar-Parsi, B. (1999), Vehicle localization on gravity maps, in ‘Unmanned Ground Vehicle Technology’, Vol. 3693, International Society for Optics and Photonics, pp. 182–191.
  * Lee et al. (2015) Lee, J., Kwon, J. & Yu, M. (2015), ‘Performance evaluation and requirements assessment for gravity gradient referenced navigation’, Sensors 15(7), 16833–16847.
  * Liu, Li & Jing (2019) Liu, F., Li, F. & Jing, X. (2019), ‘Navigability analysis of local gravity map with projection pursuit-based selection method by using gravitation field algorithm’, IEEE Access 7, 75873–75889.
  * Liu, Li, Lin & Jing (2019) Liu, F., Li, F., Lin, N. & Jing, X. (2019), ‘Gravity aided positioning based on real-time ICCP with optimized matching sequence length’, IEEE Access 7, 97440–97456.
  * Muniraj & Farhood (2017) Muniraj, D. & Farhood, M. (2017), A framework for detection of sensor attacks on small unmanned aircraft systems, in ‘2017 International Conference on Unmanned Aircraft Systems (ICUAS)’, IEEE, pp. 1189–1198.
  * Ren & Karimi (2009) Ren, M. & Karimi, H. A. (2009), ‘A hidden Markov model-based map-matching algorithm for wheelchair navigation’, The Journal of Navigation 62(3), 383–395.
  * Scales et al. (2017) Scales, K. L., Hazen, E. L., Jacox, M. G., Edwards, C. A., Boustany, A. M., Oliver, M. J. & Bograd, S. J. (2017), ‘Scale of inference: on the sensitivity of habitat models for wide-ranging marine predators to the resolution of environmental data’, Ecography 40(1), 210–220.
  * Stroock (2013) Stroock, D. W. (2013), An introduction to Markov processes, Vol. 230, Springer Science & Business Media.
  * Titterton & Weston (2004) Titterton, D. & Weston, J. (2004), Strapdown Inertial Navigation Technology, Institution of Engineering and Technology.
  * Viterbi (1967) Viterbi, A. (1967), ‘Error bounds for convolutional codes and an asymptotically optimum decoding algorithm’, IEEE Transactions on Information Theory 13(2), 260–269.
  * Wang et al. (2016) Wang, B., Yu, L., Deng, Z. & Fu, M. (2016), ‘A particle filter-based matching algorithm with gravity sample vector for underwater gravity aided navigation’, IEEE/ASME Transactions on Mechatronics 21(3), 1399–1408.
  * Wang et al. (2017a) Wang, H., Wu, L., Chai, H., Bao, L. & Wang, Y. (2017a), ‘Location accuracy of INS/gravity-integrated navigation system on the basis of ocean experiment and simulation’, Sensors 17(12), 2961.
  * Wang et al. (2022) Wang, X., Gilliam, C., Kealy, A., Close, J. & Moran, B. (2022), ‘Probabilistic map matching for robust inertial navigation aiding’, arXiv .
  * Wei et al. (2017) Wei, E., Dong, C., Liu, J., Yang, Y., Tang, S., Gong, G. & Deng, Z. (2017), ‘A robust solution of integrated SITAN with TERCOM algorithm: weight-reducing iteration technique for underwater vehicles’ gravity-aided inertial navigation system’, NAVIGATION, Journal of the Institute of Navigation 64(1), 111–122.