# Probabilistic Wind Power Forecasting with Missing Values via Adaptive
Quantile Regression

Honglin Wen Honglin Wen is with the Department of Electrical Engineering,
Shanghai Jiao Tong University, Shanghai, China, e-mail:
(linlin00@sjtu.edu.cn).

###### Abstract

Missing values challenge the probabilistic wind power forecasting at both
parameter estimation and operational forecasting stages. In this paper, we
illustrate that we are allowed to estimate forecasting functions for each
missing patterns conveniently, and propose an adaptive quantile regression
model whose parameters can adapt to missing patterns. For that, we
particularly design a feature extraction block within the quantile regression
model, where parameters are set as a function of missingness pattern and only
account for observed values. To avoid the quantile-crossing phenomena, we
design a multi-task model to ensure the monotonicity of quantiles, where
higher quantiles are derived by the addition between lower quantiles and non-
negative increments modeled by neural networks. The proposed approach is
distribution-free and applicable to both missing-at-random and missing-not-at-
random cases. Case studies demonstrate that the proposed approach achieves the
state-of-the-art in terms of the continuous ranked probability score.

###### Index Terms:

Probabilistic forecasting, machine learning, missing values, adaptive quantile
regression, quantile-crossing.

##  I Introduction

Probabilistic wind power forecasting is deemed as the workhorse to accommodate
wind power uncertainty in power system operations and electricity markets. It
leverages the information up to the current time and communicates the
probability distribution of wind power generation at a future time in the
forms of quantiles, prediction intervals, densities, etc [1]. It has been
adopted in several power system applications such as energy trading [2] as
well as reserve management [3], and has attracted increasing interest from the
industries [4].

Usually, probabilistic wind power forecasting models are developed in a data-
driven manner via parametric or non-parametric approaches [5]. The parametric
approach assumes that wind power generation follows some kind of distribution
(for instance Gaussian), and estimates the shape parameters via machine
learning methods, whereas the non-parametric approach is distribution-free.
Among the non-parametric approaches, quantile regression (QR) [6] is the most
popular one, since it is easy to use and achieves the state of the art in
several forecasting competitions. However, it requires estimating a separate
model for each quantile level, which often leads to embarrassing quantile-
crossing phenomena, i.e., higher quantiles are smaller than lower quantiles.
Recently, it has been proposed a continuous and distribution-free
probabilistic forecasting approach, which is allowed to predict the full
distribution at once by transforming the base distribution to the desired one
[7] and therefore avoids the quantile-crossing phenomena by nature.
Specifically, both the base distribution and transforms are learned via
machine learning.

Although forecasting approaches and products have developed a lot by
leveraging cutting-edge techniques [8, 9] as well as data sharing mechanisms
[10, 11, 12], the forecasting community has less investigated missing values
within datasets. In fact, missing values are inevitable in real-world data,
which may be caused by sensor failures and communication errors, for instance.
In modern statistical theory, missingness mechanisms can be classified into
missing completely at random (MCAR), missing at random (MAR), and missing not
at random (MNAR) cases, according to if missingness is dependent on observed
or missing values. For example, missingness caused by sensor failures usually
belongs to the MCAR mechanism (as it is irrelevant to observed or missing
values), whereas missingness caused by wind power curtailment belongs to the
MNAR mechanism (as it often occurs at high wind speeds, i.e., the missingness
depends on missing values). Intuitively, missing values challenge the
calculation of model prediction. For that, a natural idea is to impute missing
values before training and forecasting [13], which is referred to as the
“impute, then predict” strategy. It has also been proposed to perform the
imputation and forecasting tasks simultaneously based on deep learning
techniques [14], which is adopted in the DeepAR model [15]. However, it is
suggested that even optimal single imputed data lead to biased parameter
estimation and prediction, since the uncertainty about the missing features is
discarded [16]. Alternatively, one can leverage multiple imputations [17, 18]
on the training data, and then develop a family of models, but it brings
issues with tractability. Therefore, it remains an open issue to develop
probabilistic wind power forecasting approaches in the context of missing
values.

Compared to the inference with missing values that aims at parameter
estimation for probabilistic models or imputation [19], forecasting/prediction
with missing values focuses on the quality of forecasts. The seminar works can
date back to [20, 21], where authors investigate time series forecasting in
the context of missing values via auto-regressive moving average (ARMA) and
auto-regressive integrated moving average (ARIMA) models. They represent ARMA
and ARIMA models in the state-space form and address the calculation issues
aroused by missing values by skipping the state update. It is useful at both
model estimation and operational forecasting stages but is restricted to point
forecasts and linear models. A robust optimization approach has been proposed
in [22], which minimizes the worst-case loss when a proportion of features are
missing. Thus, it is applicable to both point and probabilistic forecasting by
setting the corresponding loss functions. However, it only focuses on the
model estimation stage and requires controlling the number of missing
features. In [23], it has been proposed a “universal imputation” strategy
based on the assumption that data are MAR, where missing values and targets
are treated equally, and the focus is on the joint distribution of features
and targets. After estimating the joint distribution, the probabilistic
forecasts can be derived by marginalization with respect to features. Though
it demonstrates better forecasting quality than the common used “impute, then
predict” strategy, it relies on the fully conditional specification technique,
the training time of which compromises its practical applications. A similar
idea can be also found in [24].

In fact, given a missingness pattern of features, one can obtain the Bayes-
optimal estimate for the forecasting function, e.g. mean function and quantile
function via the typically used forecasting paradigm [16]. It is shown that
multi-layer perceptrons (MLPs) can be Bayes consistent even in the MNAR cases
[16]. However, it is intractable to train a submodel for each missingness
pattern. For features of dimension $d$, it may require estimating $2^{d}$ sub-
models in the worst case [25]. Besides, the samples to train each submodel are
considerably reduced, as it treats each missingness pattern separately.
Therefore, it is appealing to develop models that are adaptive to several
missingness patterns [26]. In this work, we propose an adaptive quantile
regression approach by designing models with adaptive parameters and
leveraging multiplication with missingness indicator [27]. That is, we set the
parameters of the quantile regression model as a function of missing patterns.
The model mainly consists of two modules. One is responsible for feature
extraction by taking the original features with missing values as inputs,
whereas the other is designed for nonlinear mapping. By using masked weight
matrices in the feature extraction module, we are allowed to calculate model
prediction with observations. Specifically, the bias (a part of parameters) in
this module is set as the function of missingness patterns. Then, the latent
features are fed into the followed module, i.e., a group of MLPs, which yield
several quantile functions guided by corresponding pinball losses. To avoid
the quantile-crossing phenomena, we adopt a multi-task framework similar to
[28], where the higher quantiles are derived by the addition between lower
quantiles and non-negative increments modeled by neural networks. We validate
the proposed approach based on data from wind toolkit [29], where values are
removed according to designed missingness mechanisms (including MAR and MNAR
mechanisms). Case studies demonstrate that the proposed model achieves state-
of-the-art in terms of continuous ranked probability score (CRPS), especially
in MNAR cases.

In a nutshell, we mainly contribute a probabilistic forecasting approach with
missing values by designing quantile regression models with adaptive
parameters, which is applicable to both MAR and MNAR cases. The paper is
organized as follows. Section II describes the preliminaries of probabilistic
wind power forecasting and quantile regression. Section III formulates the
problem, whereas section IV presents the proposed approach. Section V presents
the setups of case studies and Section VI presents the results. Section VII
concludes this paper.

Notations: we denote random variables as uppercase letters (such as $Y$), and
their realizations as lowercase letters (such as $y$). We denote time as $t$
and use it as subscripts to represent random variables and realizations at
time $t$, for instance, $Y_{t}$ and $y_{t}$. Missing values are denoted as NA,
and the observations blurred with missing values are denoted as
$\tilde{y}_{t}\in\mathbb{R}\cup\texttt{NA}$. Missingness indicators are the
realizations of random variable $M_{t}$ and denoted as $m_{t}\in\\{0,1\\}$,
where $m_{t}=1$ implies $\tilde{y}_{t}=\texttt{NA}$ and $m_{t}=0$ implies
$\tilde{y}_{t}=y_{t}$.

##  II Preliminaries

###  II-A Probabilistic Wind Power Forecasting

Probabilistic Wind Power Forecasting aims at communicating the probabilistic
distribution of wind power generation at the future time with a lead time $k$,
i.e., $Y_{t+k}$. It often relies on a model $\mathcal{M}$ with parameters
$\boldsymbol{\theta}$, and leverages the information up to the current time
$t$, i.e., $\boldsymbol{x}_{t}$. The information $\boldsymbol{x}_{t}$ may
include weather and wind power generation at the previous time. In this work,
let us assume that $\boldsymbol{x}_{t}$ consists of lagged wind power
generation values of length $h$, i.e.,
$\boldsymbol{x}_{t}=[y_{t-h+1},y_{t-h+2},\cdots,y_{t}]^{\top}$. Denote the
cumulative distribution function (c.d.f) of $Y_{t+k}$ as $F_{t+k}(y)$. Then,
probabilistic wind power forecasting can be described as

$\hat{F}_{t+k}(y)=F_{t+k}(y|\boldsymbol{x}_{t};\mathcal{M},\boldsymbol{\theta}).$
(1)

The model $\mathcal{M}$ can be set as some kind of distributional model such
as logit-normal distribution [8]. Alternatively, it can be set as a group of
increasing quantiles. Denote the $\alpha$-th quantile of $Y_{t+k}$ as
$q_{t+k}^{\alpha}$, which is defined as

$q_{t+k}^{\alpha}=F^{-1}_{t+k}(\alpha).$ (2)

Then the forecast for distribution $F_{t+k}(y)$ can be also derived as

$\displaystyle\\{\hat{q}_{t+k}^{\alpha_{1}},\hat{q}_{t+k}^{\alpha_{2}},\cdots,\hat{q}_{t+k}^{\alpha_{m}}\\},\
\alpha_{1}<\alpha_{2}<\cdots<\alpha_{m},$

where $\hat{q}_{t+k}^{\alpha_{i}}$ is the estimated $\alpha_{i}$-th quantile.
Such quantiles can be derived via quantile regression [6]. Let
$g(\boldsymbol{x};\boldsymbol{\theta},\alpha)$ represent a quantile function,
such that

$q_{t+k}^{\alpha}=g(\boldsymbol{x}_{t};\boldsymbol{\theta},\alpha).$ (3)

In particular, it can be represented as the inner product of coefficient and
latent features learned via machine learning [30]. Denote the coefficient as
$\boldsymbol{w}$ and the features as
$\phi(\boldsymbol{x}_{t};\boldsymbol{\theta}_{\phi})$, i.e.,
$\boldsymbol{\theta}=\\{\boldsymbol{w},\boldsymbol{\theta}_{\phi}\\}$. We
describe the function as

$g(\boldsymbol{x}_{t};\boldsymbol{\theta},\alpha)=\boldsymbol{w}^{\top}\phi(\boldsymbol{x}_{t};\boldsymbol{\theta}_{\phi}),$
(4)

The parameter $\boldsymbol{\theta}$ can be estimated via machine learning
based on historical data $\\{(\boldsymbol{x}_{t},y_{t+k})|t=1,2,\cdots,n\\}$.
Concretely, they are estimated by minimizing the pinball loss $\mathcal{L}$,
i.e.,

$\displaystyle\mathcal{L}=\frac{1}{n}\sum_{t}\ell(y_{t+k},g(\boldsymbol{x}_{t};\boldsymbol{\theta},\alpha))$

where $\ell(y_{t+k},g(\boldsymbol{x}_{t};\boldsymbol{\theta},\alpha))$ is
defined as

$\displaystyle\ell(y_{t+k},g(\boldsymbol{x}_{t};\boldsymbol{\theta},\alpha))=$
(5)
$\displaystyle\max(\alpha(y_{t+k}-g(\boldsymbol{x}_{t};\boldsymbol{\theta},\alpha)),(\alpha-1)(y_{t+k}-g(\boldsymbol{x}_{t};\boldsymbol{\theta},\alpha))).$

Figure 1: The sketch of the proposed approach, where blank blocks in
$\tilde{\boldsymbol{x}}_{t}$ indicate missing values.

###  II-B Missingness Mechanism

In the modern statistical theory [19], missingness mechanisms can be
classified into three categories: MCAR, MAR, and MNAR. Intuitively, the
observations for both $\boldsymbol{x}_{t}$ and $y_{t+k}$ may contain missing
values. Therefore, we write the observations as $\tilde{\boldsymbol{x}}_{t}$
and $\tilde{y}_{t+k}$. We introduce a missingness indicator vector
$\boldsymbol{m}_{t}\in\\{0,1\\}^{d}$ for $\tilde{\boldsymbol{x}}_{t}$, which
is the realization of a random variable $\boldsymbol{M}_{t}$, s.t.,
$\tilde{x}_{t,i}=x_{t,i}$ when $m_{t,i}=0$ and $\tilde{x}_{t,i}=\texttt{NA}$
when $m_{t,i}=1$. The observed part of $\boldsymbol{x}_{t}$ is denoted as
$\boldsymbol{x}_{o,t}=o(\boldsymbol{x}_{t},\boldsymbol{m}_{t})$, whereas the
missing part of $\boldsymbol{x}_{t}$ is denoted as
$\boldsymbol{x}_{m,t}=o(\boldsymbol{x}_{t},1-\boldsymbol{m}_{t})$.

Taking the multivariable $\boldsymbol{x}_{t}$ as an example, the parametric
model for the joint distribution of the data sample and its mask is described
as

$f(\boldsymbol{x}_{t},\boldsymbol{m}_{t};\xi,\psi)=f(\boldsymbol{x}_{t};\xi)f(\boldsymbol{m}_{t}|\boldsymbol{x}_{t};\psi),$
(6)

where $\xi$ and $\psi$ represent the parameters of distribution for data and
mask. The data sample can be split into an observed part
$\boldsymbol{x}_{o,t}$ and a missing part $\boldsymbol{x}_{m,t}$, i.e.,
$\boldsymbol{x}_{t}=(\boldsymbol{x}_{o,t},\boldsymbol{x}_{m,t})$. Then the
missingness mechanism is referred to as MCAR if

$\displaystyle
f(\boldsymbol{m}_{t}|\boldsymbol{x}_{t};\psi)=f(\boldsymbol{m}_{t};\psi),$

which means the missingness is independent of the data sample. The missingness
mechanism is referred to as MAR if

$\displaystyle
f(\boldsymbol{m}_{t}|\boldsymbol{x}_{t};\psi)=f(\boldsymbol{m}_{t}|\boldsymbol{x}_{o,t};\psi).$

That is, the missingness is dependent on $\boldsymbol{x}_{o,t}$, yet
independent of $\boldsymbol{x}_{m,t}$. Similarly, the missingness mechanism is
referred to as MNAR if

$\displaystyle
f(\boldsymbol{m}_{t}|\boldsymbol{x}_{t};\psi)=f(\boldsymbol{m}_{t}|\boldsymbol{x}_{t}^{o},\boldsymbol{x}_{t}^{m};\psi),$

i.e., the missingness is dependent on both the observed and missing values. In
fact, both the real data and missing mechanisms are unavailable.
Alternatively, we are only allowed to get access to
$(\tilde{\boldsymbol{x}}_{t},\tilde{y}_{t+k})$. In this work, we aim at
developing a forecasting approach based on only observations without
assumptions on the missingness mechanisms.

##  III Methodology

In this section, we illustrate the approach to develop models directly with
observations $(\tilde{\boldsymbol{x}}_{t},\tilde{y}_{t+k})$. In what follows,
we first discuss the challenge within the probabilistic wind power forecasting
aroused by missing values, i.e., the discrete nature of NA. Then, we describe
the idea to deal with such an issue, that is, adapting the model parameters
with missingness patterns.

###  III-A Problem Formulation

As $(\tilde{\boldsymbol{x}}_{t},\tilde{y}_{t+k})$ has no contributions to the
parameter estimation if $\tilde{y}_{t+k}$ is missing, we only focus on samples
where $\tilde{y}_{t+k}$ is observed in the following context. Then, the
estimate of parameters is derived via

$\hat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}}\frac{1}{n}\sum_{t}\ell(\tilde{y}_{t+k},g(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta},\alpha)),$
(7)

though $\tilde{\boldsymbol{x}}_{t}$ contains missing values. Unfortunately,
most off-the-shelf machine learning methods do not work with the half-discrete
nature of $\mathbb{R}\cup\texttt{NA}$. Alternatively, we are allowed to
estimate a quantile function for each missingness pattern by selecting samples
with such missingness pattern, as suggested by [16]. For instance, for a
missingness pattern $\mathcal{P}_{i}\in\\{0,1\\}^{d}$, we denote the
corresponding quantile function as $g_{\mathcal{P}_{i}}$ with parameters
$\boldsymbol{\theta}_{i}$. It is expressed as

$g_{\mathcal{P}_{i}}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{i},\alpha)=g(\boldsymbol{x}_{o,t};\alpha,\boldsymbol{\theta}_{i},\boldsymbol{m}_{t}=\mathcal{P}_{i}).$
(8)

Therefore, the quantile function can be constructed as a combination of
functions

$g(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta},\alpha)=\left\\{\begin{aligned}
g_{\mathcal{P}_{1}}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{1},\alpha),\
\boldsymbol{m}_{t}=\mathcal{P}_{1},\\\
g_{\mathcal{P}_{2}}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{2},\alpha),\
\boldsymbol{m}_{t}=\mathcal{P}_{2},\\\ \cdots\\\
g_{\mathcal{P}_{2^{d}}}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{2^{d}},\alpha),\
\boldsymbol{m}_{t}=\mathcal{P}_{2^{d}}.\end{aligned}\right.$ (9)

Intuitively, it requires estimating $2^{d}$ sub-models, which scales
exponentially with the dimension and is therefore impractical.

###  III-B Method

To address the aforementioned tractability issue, we seek to develop a model
whose parameters can adapt to missingness patterns. It is described as

$\displaystyle
g(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}(\boldsymbol{m}_{t}),\alpha).$

For that, we expect the feature function $\phi(\cdot)$ can adapt to
missingness patterns, i.e.,

$g(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}(\boldsymbol{m}_{t}),\alpha)=\boldsymbol{w}^{\top}\phi(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{\phi}(\boldsymbol{m}_{t})).$
(10)

Specifically, we design the feature function $\phi(\cdot)$ as a composition of
$\phi_{1}(\cdot)$ and $\phi_{2}(\cdot)$, where $\phi_{1}(\cdot)$ is a linear
function with parameters $\boldsymbol{\theta}_{\phi_{1}}$ and
$\phi_{2}(\cdot)$ is a nonlinear function with parameters
$\boldsymbol{\theta}_{\phi_{2}}$. We have

$\phi(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{\phi}(\boldsymbol{m}_{t}))=\phi_{2}(\phi_{1}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{\phi_{1}}(\boldsymbol{m}_{t})),\boldsymbol{\theta}_{\phi_{2}}),$
(11)

where $\boldsymbol{\theta}_{\phi_{1}}$ is a function of $\boldsymbol{m}_{t}$.
Let $\boldsymbol{z}_{t}$ denote the output of
$\phi_{1}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{\phi_{1}}(\boldsymbol{m}_{t}))$,
i.e.,

$\displaystyle\boldsymbol{z}_{t}=\phi_{1}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{\phi_{1}}(\boldsymbol{m}_{t})).$

Obviously, NA will not contribute to $\boldsymbol{z}_{t}$. Thus, we define
$\phi_{1}(\cdot)$ as a function that outputs the linear combination of
observations in $\tilde{\boldsymbol{x}}_{t}$, which is described as

$\displaystyle z_{t,i}=$
$\displaystyle\sum_{j:m_{t,j}=0}\boldsymbol{W}_{\phi_{1}}(\boldsymbol{m}_{t})[i,j]\tilde{x}_{t,j}+\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})[i]$
(12) $\displaystyle=$
$\displaystyle\boldsymbol{W}_{\phi_{1}}(\boldsymbol{m}_{t})[i,:]{\rm
diag}(1-\boldsymbol{m}_{t})\tilde{\boldsymbol{x}}_{t}+\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})[i],$

where $\boldsymbol{W}_{\phi_{1}}(\boldsymbol{m}_{t})$ and
$\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})$ are the parameters of
$\phi_{1}(\cdot)$, i.e.,
$\boldsymbol{\theta}_{\phi_{1}}=\\{\boldsymbol{W}_{\phi_{1}}(\boldsymbol{m}_{t}),\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})\\}$,
and ${\rm diag}(\cdot)$ returns a square diagonal matrix with the elements of
the input vector. Then $\boldsymbol{z}_{t}$ is used as the input to
$\phi_{2}(\cdot)$, i.e.,

$\displaystyle
q_{t+k}^{\alpha}=\boldsymbol{w}^{\top}\phi_{2}(\boldsymbol{z}_{t};\boldsymbol{\theta}_{\phi_{2}}).$

##  IV Proposed Approach

With the main idea described in Section III, we now describe the proposed
approach in detail. The linear function $\phi_{1}(\cdot)$ is implemented
similar to the NeuMiss block proposed in [27], whereas the nonlinear function
$\phi_{2}(\cdot)$ is implemented as an MLP. In particular, we place the
quantile regression model in a multi-task framework and ensure the
monotonicity of quantiles. We sketch the proposed approach in Figure 1, which
consists of a feature extraction block and several MLPs correlated with
previous ones via addition successively.

###  IV-A Feature Extraction

In the feature extraction block, we replace NA in $\tilde{\boldsymbol{x}}_{t}$
with 0, and denote it as $\hat{\boldsymbol{x}}_{t}$, which operates as

$\displaystyle\hat{\boldsymbol{x}}_{t}=\tilde{\boldsymbol{x}}_{t}\odot(1-\boldsymbol{m}_{t}),$

where $\odot$ is an elementwise product operator. Then equation (12) can be
also rewritten compactly as

$\boldsymbol{z}_{t}=\boldsymbol{W}_{\phi_{1}}(\boldsymbol{m}_{t})\hat{\boldsymbol{x}}_{t}+\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})$
(13)

In particular, we use a special case of adaptive regression where
$\boldsymbol{W}_{\phi_{1}}(\boldsymbol{m}_{t})$ are static whereas
$\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})$ is a function of the
missingness patterns. Specifically, we set
$\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})$ as

$\boldsymbol{b}_{\phi_{1}}(\boldsymbol{m}_{t})=\boldsymbol{b}_{\phi_{1}}\odot\boldsymbol{m}_{t}.$
(14)

Then equation (13) can be rewritten as

$\boldsymbol{z}_{t}=\boldsymbol{W}_{\phi_{1}}\hat{\boldsymbol{x}}_{t}+\boldsymbol{b}_{\phi_{1}}\odot\boldsymbol{m}_{t}.$
(15)

In this work, we design the feature extraction block similar to the NeuMiss
model proposed in [27], which is illustrated in Figure 2. Consider that we
stack $l_{\phi_{1}}$ such blocks via skip connections. For the $l$-th block,
we denote its input as $\boldsymbol{h}^{l-1}_{\phi_{1}}$, and its output as
$\boldsymbol{h}^{l}_{\phi_{1}}$. Specially,
$\boldsymbol{h}^{0}_{\phi_{1}}=\hat{\boldsymbol{x}}_{t}$. For each block, we
have

$\boldsymbol{h}^{l}_{\phi_{1}}=\boldsymbol{W}_{\phi_{1}}^{l}\boldsymbol{h}^{l-1}_{\phi_{1}}\odot(1-\boldsymbol{m}_{t})+\boldsymbol{h}^{0}_{\phi_{1}},$
(16)

where $\boldsymbol{W}_{\phi_{1}}^{l}$ is the weight matrix in the $l$-th
block. The output of the last block is denoted as
$\boldsymbol{h}^{l_{\phi_{1}}}_{\phi_{1}}$. Then, the latent features
$\boldsymbol{z}_{t}$ is derived by adding
$\boldsymbol{h}^{l_{\phi_{1}}}_{\phi_{1}}$ with
$\boldsymbol{b}_{\phi_{1}}\odot\boldsymbol{m}_{t}$, i.e.,

$\boldsymbol{z}_{t}=\boldsymbol{h}^{l_{\phi_{1}}}_{\phi_{1}}+\boldsymbol{b}_{\phi_{1}}\odot\boldsymbol{m}_{t},$
(17)

which is further fed into nonlinear functions to yield quantiles.

Figure 2: The structure of feature extraction block.

###  IV-B Non-crossing Quantile Regression Neural Network

Specifically, we define a non-negative nonlinear function
$\phi_{2,\alpha_{i}}$ for each quantile level $\alpha_{i}$. Accordingly, we
respectively denote the coefficient and parameters as
$\boldsymbol{w}_{\alpha_{i}}$ and $\boldsymbol{\theta}_{\phi_{2},\alpha_{i}}$.
The quantile function $q_{t+k}^{\alpha_{1}}$ is derived by the composition of
$\phi_{1}(\cdot)$ and $\phi_{2,\alpha_{1}}(\cdot)$, i.e.,

$\begin{split}q_{t+k}^{\alpha_{1}}&=g(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{1},\alpha_{1})\\\
&=\boldsymbol{w}^{\top}_{\alpha_{1}}\phi_{2,\alpha_{\phi_{1}}}(\phi_{1}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{\phi_{1}});\boldsymbol{\theta}_{\phi_{2},\alpha_{1}}).\end{split}$
(18)

As for the quantile level $\alpha_{i}\ (i>1)$, we derive the quantile function
$q_{t+k}^{\alpha_{i}}$ by adding $q_{t+k}^{\alpha_{i-1}}$ with an increment,
which is set as the composition of $\phi_{1}(\cdot)$ and
$\phi_{2,\alpha_{i}}(\cdot)$. It is described as

$\begin{split}q_{t+k}^{\alpha_{i}}&=g(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{i},\alpha_{i})\\\
&=q_{t+k}^{\alpha_{i-1}}+\boldsymbol{w}^{\top}_{\alpha_{i}}\phi_{2,\alpha_{i}}(\phi_{1}(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{\phi_{1}});\boldsymbol{\theta}_{\phi_{2},\alpha_{i}}).\end{split}$
(19)

As $\phi_{2,\alpha_{i}}$ is non-negative, we have

$\displaystyle q_{t+k}^{\alpha_{i-1}}\leq q_{t+k}^{\alpha_{i}}.$

Each function $\phi_{2,\alpha_{i}}$ is implemented with an
$l_{\phi_{2}}$-layer MLP, where each layer takes
$\boldsymbol{h}^{l-1}_{\phi_{2},\alpha_{i}}$ as input, and output
$\boldsymbol{h}^{l}_{\phi_{2},\alpha_{i}}$. Specially, $\boldsymbol{z}_{t}$ is
denoted as $\boldsymbol{h}^{0}_{\phi_{2},\alpha_{i}}$. The $l$-th layer
operates as

$\boldsymbol{h}^{l}_{\phi_{2},\alpha_{i}}=\boldsymbol{W}_{\phi_{2},\alpha_{i}}^{l}\boldsymbol{h}^{l-1}_{\phi_{2},\alpha_{i}}+\boldsymbol{b}_{\phi_{2},\alpha_{i}}^{l},$
(20)

where $\boldsymbol{W}_{\phi_{2},\alpha_{i}}^{l}$ and
$\boldsymbol{b}_{\phi_{2},\alpha_{i}}^{l}$ respectively represent the weight
and bias in this layer. It is followed by a Relu function $\sigma(\cdot)$ ,
which operates as

$\sigma(\boldsymbol{h}^{l}_{\phi_{2},\alpha_{i}})=\max(\boldsymbol{h}^{l}_{\phi_{2},\alpha_{i}},\boldsymbol{0}),$
(21)

where $\max$ returns the maximum between
$\boldsymbol{h}^{l}_{\phi_{2},\alpha_{i}}$ and $\boldsymbol{0}$ elementwisely.

To estimate all the parameters, we minimize the total loss, which is defined
as

$\mathcal{L}=\frac{1}{n\times
m}\sum_{t}\sum_{i}\ell(\tilde{y}_{t+k},g(\tilde{\boldsymbol{x}}_{t};\boldsymbol{\theta}_{i},\alpha_{i})).$
(22)

##  V Case Study

In this section, we describe the data for case validation, which come from the
open-source wind toolkit [29]. Then we introduce experimental setups, where
missingness is simulated based on different mechanisms. After that, we
describe benchmark models and qualification metrics.

###  V-A Data Description

As we focus on very short-term wind power forecasting, here we use a wind
power measurement dataset collected from a wind farm located in South
Carolina. It is an hourly dataset where values are normalized by its capacity,
and spans from 2007 to 2013. In each case, we split the first $70\%$ of data
for training models, the following $10\%$ of data for validation, and the last
$20\%$ of data for genuine forecast verification.

###  V-B Experimental Setups

In this work, we simulate missingness based on both MAR and MNAR mechanisms.
For the MAR mechanism, we consider missingness that spread sporadically and in
blocks, though it is independent of the wind power generation value. As for
the MNAR mechanism, we consider a self-masking case, where values greater than
a threshold will be missing. The designed 3 cases are described as follows:

####  V-B1 Case 1

In this case, we randomly remove $20\%$ of data, which spread sporadically
over the dataset.

####  V-B2 Case 2

In this case, we remove data in blocks, which are randomly located over the
dataset. The length of each block is uniformly distributed between 5 and 30
steps, whereas the number of blocks is fixed as 300.

####  V-B3 Case 3

In this case, we remove data whose values are greater than 0.87.

In each case, we consider the lead time of 1, 2, and 3, and use wind power
generation values at previous time as features. As feature selection is not
the focus of this work, we choose it empirically as 6 lags. Sophisticated
feature selection approaches can be certainly used.

###  V-C Benchmark Models

Three types of models are considered as benchmarks: naive models, “impute,
then predict” strategy-based models, and “universal imputation” strategy-based
models. Besides, we set a quantile regression model trained on complete data
as a reference. We describe them as follows:

####  V-C1 Climatology

It is a naive model that estimates the empirical distribution of wind power
generation based on historical samples.

####  V-C2 IM-Gaussian

It is an “impute, then predict” strategy-based model. Missing values are
imputed via MissForest [31], based on which a parametric probabilistic
forecasting model with Gaussian distributional assumption is developed.

####  V-C3 IM-QR

It is an “impute, then predict” strategy-based model. Missing values are also
imputed via MissForest, based on which quantile regression models are
developed.

####  V-C4 DeepAR

It is a state-of-the-art “impute, then predict” strategy-based model where
imputation and forecasting are performed simultaneously.

####  V-C5 UI

It is a “universal imputation” strategy-based model proposed in [23] based on
a fully conditional specification [18]. For further descriptions, we refer
readers to [23].

####  V-C6 R-QR

It is a quantile regression model trained on the complete dataset, which
serves as a reference.

###  V-D Qualification Metrics

To assess the quality of forecasts, we verify the calibration and sharpness of
forecasts. Concretely, the calibration of predictive densities is assessed
with reliability diagrams. The sharpness of forecasts is assessed with the
average width of central prediction intervals, which reveals how predictive
densities concentrate the information. In addition, we assess the quality of
forecasts with a skill score, namely continuous ranked probability score
(CRPS). Given the lead time $k$, we denote the predictive c.d.f for wind power
generation at time $t+k$ as $\hat{F}_{t+k}$, and the real generation value as
$y_{t+k}$. The CRPS is calculated via

$\displaystyle{\rm
CRPS}(\hat{F}_{t+k},y_{t+k})=\int_{y}(\hat{F}_{t+k}(y)-\mathcal{I}(y-y_{t+k}))^{2}dy,$

where $\mathcal{I}(\cdot)$ is a step function. We report the average CRPS of
all test samples for each lead time. For further information on forecast
verification, we refer readers to [1].

Figure 3: 90% prediction intervals by the proposed model for 144 successive observations in case 1 when lead time is 1. TABLE I: The CRPS values of forecasts by the proposed and benchmark models with different lead times in case 1 (%). Lead time | Climatology | IM-Gaussian | IM-QR | DeepAR | UI | R-QR | Proposed
---|---|---|---|---|---|---|---
1 | 18.6 | 7.5 | 7.8 | 7.8 | 6.9 | 6.9 | 7.2
2 | 18.6 | 9.9 | 9.9 | 10.2 | 9.1 | 9.3 | 9.4
3 | 18.6 | 11.7 | 11.7 | 12.1 | 10.9 | 11.2 | 11.4

(a) Reliability diagram

(b) Prediction interval width

Figure 4: Reliability and sharpness of forecasts by the proposed and benchmark models with $k=1$ in case 1. TABLE II: The CRPS values of forecasts by the proposed and benchmark models with different lead times in case 2 (%). Lead time | Climatology | IM-Gaussian | IM-QR | DeepAR | UI | R-QR | Proposed
---|---|---|---|---|---|---|---
1 | 18.6 | 7.2 | 7.4 | 7.4 | 6.6 | 6.9 | 6.4
2 | 18.6 | 9.6 | 9.6 | 9.9 | 8.9 | 9.3 | 9.2
3 | 18.6 | 11.5 | 11.5 | 11.7 | 11.9 | 11.2 | 11.2

(a) Reliability diagram

(b) Prediction interval width

Figure 5: Reliability and sharpness of forecasts by the proposed and benchmark
models with $k=1$ in case 2.

##  VI Results

The results of the aforementioned three cases are presented in this section,
and followed by discussion.

###  VI-A Case 1

The CRPS values of forecasts by the proposed and benchmark models are
presented in Table I. Intuitively, The performance of climatology is the worst
among all models, as it only communicates forecasts via unconditional
empirical distribution. Unlike the common situations where QR often
outperforms Gaussian distributional models, The performance of IM-Gaussian and
IM-QR is quite close, which suggests that imputation has a complex impact on
the training of downstream forecasting models. To our surprise, the
performance of DeepAR is the worst among the three “impute, then predict”
strategy-based models. In fact, DeepAR imputes missing values via the
intermediate results of the recurrent neural network during training, which
may have a negative impact on forecasting. The performance of UI and the
proposed model is comparable to the reference model. Especially, the UI model
slightly outperforms the reference model, which may be due to its robustness
to overfitting. The 90% prediction intervals by the proposed model for 144
successive observations when lead time is 1 is shown in Figure 3.

We present the reliability diagrams and prediction interval width of forecasts
when $k=1$ in Figure 4. As shown, the reliability diagrams of “impute, then
predict” strategy-based models deviate from the ideal case in relatively
further distances. The prediction interval widths of the proposed and
benchmark models are comparable, all of which are smaller than that of the
reference model. The reliability of DeepAR is the worst among all models,
whereas its prediction interval widths are also the smallest. The reliability
diagrams of UI and the proposed model are close to the ideal case, which
reveals that they lead to little biases in the forecasts. Surprisingly, the
reliability of the reference model is worse than the proposed and UI model,
while the reference model yields larger prediction interval widths. It may be
caused by the overfitting of the reference model.

TABLE III: The CRPS values of forecasts by the proposed and benchmark models with different lead times in case 3 (percentage). Lead time | Climatology | IM-Gaussian | IM-QR | DeepAR | UI | R-QR | Proposed
---|---|---|---|---|---|---|---
1 | 15.4 | 6.6 | 6.8 | 7.1 | 6.4 | 6.0 | 6.0
2 | 15.4 | 8.6 | 8.7 | 9.3 | 8.3 | 7.9 | 8.0
3 | 15.4 | 10.1 | 9.9 | 10.8 | 9.5 | 9.4 | 9.4

(a) Reliability diagram

(b) Prediction interval width

Figure 6: Reliability and sharpness of forecasts by the proposed and benchmark
models with $k=1$ in case 3.

###  VI-B Case 2

The CRPS values of forecasts by the proposed and benchmark models are
presented in Table II. In this case, the differences in CRPS values among all
models are smaller than those in case 1. Different from the case 1,
missingness occurs in blocks in this case, leading to more samples with
complete observations. Therefore, missing values have less impact on the
quality of forecasts. Among the “impute, then predict” strategy-based models,
the performance of DeepAR is still the worst, though the difference between
DeepAR and IM-Gaussian/IM-QR is smaller compared to that in case 1. By
contrast, the performance of the proposed and UI models is still better than
that of “impute, then predict” strategy-based models, which is comparable to
the reference model. It suggests that the proposed and UI models are
applicable to both sporadic and block missingness.

The reliability diagrams and prediction interval width of forecasts when $k=1$
are presented in Figure 5. In this case, the reliability of DeepAR is still
the worst among all models. It is seen that the reliability diagram of the
proposed model fluctuates around the ideal case, though close to the ideal
case. It is caused by the monotonicity constraint on the proposed model. Such
constraints ensure that higher quantiles are no smaller than lower quantiles,
but influence the parameter estimation on the other hand. Further analysis is
contained in the following subsection. In general, the performance of the
proposed model in reliability and sharpness is quite good.

###  VI-C Case 3

We present the CRPS values of forecasts by the proposed and benchmark models
in Table III. Compared to cases 1 and 2, all models achieve better performance
in CRPS, as high wind power generation values are excluded, which leads to
larger forecast errors. Still, the quality of forecasts by the DeepAR model is
the worst among all models. The performance of UI and the proposed models is
better than that of “impute, then predict” strategy-based models. Especially
the proposed model achieves the best performance among all models, which
reveals that the proposed model is applicable to MNAR cases.

The reliability diagrams and prediction interval widths are shown in Figure 6.
Unlike in cases 1 and 2, the reliability of UI model deviates from the ideal
case to a large extent, which may be due to the fact that fully conditional
specification relies on MAR assumption. In the MNAR case, such an assumption
leads to more biases in the parameter estimation. By contrast, the proposed
model is more applicable to MNAR cases.

###  VI-D Discussion on Non-crossing Constraints

(a) Reliability diagram

(b) Prediction interval width

Figure 7: Reliability and sharpness of forecasts by the proposed and regular quantile regression model with $k=1$ based on the setting in case 2. TABLE IV: The training time of the proposed and benchmark models in case 1 (minutes). Models | Climatology | IM-Gaussian | IM-QR | DeepAR | UI | R-QR | Proposed
---|---|---|---|---|---|---|---
Training time | - | 32 | 10 | 64 | 41 | 10 | 6

As described in Section IV.B, we place hard constraints on quantiles to ensure
monotonicity by setting higher quantiles as the addition of lower quantiles
and non-negative increments. Certainly, it avoids the embarrassing quantile
crossing phenomena. However, such constraints also influence parameter
estimation. For illustration, we present the reliability diagrams and
prediction interval widths of regular quantile regression and the proposed
model based on the setting in case 2 in Figure 7. As shown, the regular
quantile regression model leads to a quantile-crossing phenomenon, while the
prediction interval widths are wider than those of the regular quantile
regression model. In addition, as the proposed model constructs quantiles via
addition operation, it can be only calculated sequentially, while the regular
quantile regression model is allowed to yield quantiles in parallel.

###  VI-E Training Time

We present the training time of all models in Table IV, where time spent on
imputation for the “impute, then predict” strategy-based models is not
included. It is seen that the proposed model is time-efficient in training. In
particular, the training time of UI model (which is based on fully conditional
specification here) will increase linearly with the dimension, whereas that of
the proposed model scales sublinearly thanks to the deep learning techniques.

##  VII Conclusion

In this work, we propose an adaptive quantile regression approach for
probabilistic wind power forecasting with missing values within the typically
used conditional distribution modeling framework. It is based on deep neural
network models, and contains a linear mapping whose bias is adaptive to
missingness patterns by design and a nonlinear mapping that is responsible for
yielding quantiles. It is applicable to both missing at-random and missing
not-at-random cases. In particular, higher quantiles are derived by the
addition between lower quantiles and non-negative increments, which avoids the
embarrassing quantile-crossing phenomena. Case studies demonstrate that the
proposed approach achieves state-of-the-art in terms of CRPS in both missing
at-random and missing not-at-random cases, and is time-efficient in training.
However, we only consider adaptive bias in this work, more efforts can be
taken into developing adaptive weights in the future. Besides, as case studies
suggest the hard constraints on quantiles influence the parameter estimation
while they ensure the monotonicity, advanced non-crossing quantile regression
methods can be investigated.

## References

  * [1] T. Gneiting and M. Katzfuss, “Probabilistic forecasting,” _Annual Review of Statistics and Its Application_ , vol. 1, pp. 125–151, 2014.
  * [2] P. Pinson, C. Chevallier, and G. N. Kariniotakis, “Trading wind generation from short-term probabilistic forecasts of wind power,” _IEEE Transactions on Power Systems_ , vol. 22, no. 3, pp. 1148–1156, 2007.
  * [3] M. A. Matos and R. J. Bessa, “Setting the operating reserve using probabilistic wind power forecasts,” _IEEE transactions on power systems_ , vol. 26, no. 2, pp. 594–603, 2010.
  * [4] S. E. Haupt, M. G. Casado, M. Davidson, J. Dobschinski, P. Du, M. Lange, T. Miller, C. Mohrlen, A. Motley, R. Pestana _et al._ , “The use of probabilistic forecasts: Applying them in theory and practice,” _IEEE Power and Energy Magazine_ , vol. 17, no. 6, pp. 46–57, 2019.
  * [5] C. Sweeney, R. J. Bessa, J. Browell, and P. Pinson, “The future of forecasting for renewable energy,” _Wiley Interdisciplinary Reviews: Energy and Environment_ , vol. 9, no. 2, p. e365, 2020.
  * [6] R. Koenker and K. F. Hallock, “Quantile regression,” _Journal of economic perspectives_ , vol. 15, no. 4, pp. 143–156, 2001.
  * [7] H. Wen, P. Pinson, J. Ma, J. Gu, and Z. Jin, “Continuous and distribution-free probabilistic wind power forecasting: A conditional normalizing flow approach,” _IEEE Transactions on Sustainable Energy_ , vol. 13, no. 4, pp. 2250–2263, 2022.
  * [8] H. Wen, J. Ma, J. Gu, L. Yuan, and Z. Jin, “Sparse variational gaussian process based day-ahead probabilistic wind power forecasting,” _IEEE Transactions on Sustainable Energy_ , vol. 13, no. 2, pp. 957–970, 2022.
  * [9] H. Zhang, Y. Liu, J. Yan, S. Han, L. Li, and Q. Long, “Improved deep mixture density network for regional wind power probabilistic forecasting,” _IEEE Transactions on Power Systems_ , vol. 35, no. 4, pp. 2549–2560, 2020\.
  * [10] L. Cavalcante, R. J. Bessa, M. Reis, and J. Browell, “Lasso vector autoregression structures for very short-term wind power forecasting,” _Wind Energy_ , vol. 20, no. 4, pp. 657–675, 2017.
  * [11] C. Goncalves, R. J. Bessa, and P. Pinson, “Privacy-preserving distributed learning for renewable energy forecasting,” _IEEE Transactions on Sustainable Energy_ , vol. 12, no. 3, pp. 1777–1787, 2021.
  * [12] C. Goncalves, P. Pinson, and R. J. Bessa, “Towards data markets in renewable energy forecasting,” _IEEE Transactions on Sustainable Energy_ , vol. 12, no. 1, pp. 533–542, 2020.
  * [13] R. Tawn, J. Browell, and I. Dinwoodie, “Missing data in wind farm time series: Properties and effect on forecasts,” _Electric Power Systems Research_ , vol. 189, p. 106640, 2020.
  * [14] W. Cao, D. Wang, J. Li, H. Zhou, L. Li, and Y. Li, “Brits: Bidirectional recurrent imputation for time series,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [15] D. Salinas, V. Flunkert, J. Gasthaus, and T. Januschowski, “Deepar: Probabilistic forecasting with autoregressive recurrent networks,” _International Journal of Forecasting_ , vol. 36, no. 3, pp. 1181–1191, 2020\.
  * [16] J. Josse, N. Prost, E. Scornet, and G. Varoquaux, “On the consistency of supervised learning with missing values,” _arXiv preprint arXiv:1902.06931_ , 2019.
  * [17] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the em algorithm,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 39, no. 1, pp. 1–22, 1977\.
  * [18] S. Van Buuren, J. P. Brand, C. G. Groothuis-Oudshoorn, and D. B. Rubin, “Fully conditional specification in multivariate imputation,” _Journal of statistical computation and simulation_ , vol. 76, no. 12, pp. 1049–1064, 2006\.
  * [19] R. J. Little and D. B. Rubin, _Statistical analysis with missing data_. John Wiley & Sons, 2019, vol. 793\.
  * [20] R. H. Jones, “Maximum likelihood fitting of arma models to time series with missing observations,” _Technometrics_ , vol. 22, no. 3, pp. 389–395, 1980\.
  * [21] R. Kohn and C. F. Ansley, “Estimation, prediction, and interpolation for arima models with missing data,” _Journal of the American statistical Association_ , vol. 81, no. 395, pp. 751–761, 1986.
  * [22] A. Stratigakos, P. Andrianesis, A. Michiorri, and G. Kariniotakis, “Towards resilient energy forecasting: A robust optimization approach,” 2022.
  * [23] H. Wen, P. Pinson, J. Gu, and Z. Jin, “Wind energy forecasting with missing values within a fully conditional specification framework,” _International Journal of Forecasting_ , 2023.
  * [24] N. B. Ipsen, P.-A. Mattei, and J. Frellsen, “How to deal with missing data in supervised deep learning?” in _ICLR 2022-10th International Conference on Learning Representations_ , 2022.
  * [25] P. R. Rosenbaum and D. B. Rubin, “Reducing bias in observational studies using subclassification on the propensity score,” _Journal of the American statistical Association_ , vol. 79, no. 387, pp. 516–524, 1984.
  * [26] D. Bertsimas, A. Delarue, and J. Pauphilet, “Prediction with missing data,” _stat_ , vol. 1050, p. 7, 2021.
  * [27] M. Le Morvan, J. Josse, T. Moreau, E. Scornet, and G. Varoquaux, “Neumiss networks: differentiable programming for supervised learning with missing values.” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 5980–5990, 2020.
  * [28] Y. Park, D. Maddix, F.-X. Aubet, K. Kan, J. Gasthaus, and Y. Wang, “Learning quantile functions without quantile crossing for distribution-free time series forecasting,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2022, pp. 8127–8150.
  * [29] C. Draxl, A. Clifton, B.-M. Hodge, and J. McCaa, “The wind integration national dataset (wind) toolkit,” _Applied Energy_ , vol. 151, pp. 355–366, 2015.
  * [30] C. M. Bishop and N. M. Nasrabadi, _Pattern recognition and machine learning_. Springer, 2006, vol. 4, no. 4.
  * [31] D. J. Stekhoven and P. Bühlmann, “Missforest—non-parametric missing value imputation for mixed-type data,” _Bioinformatics_ , vol. 28, no. 1, pp. 112–118, 2012.


§.§ Problem Formulation

We consider finding the pose of a MAMR, with frame {$\mathcal{A}$} within our workspace  {$\mathcal{W}$} (see Fig. <ref>). Its position is denoted as $ \textbf{p} \in \mathbb{R}^3$ in {$\mathcal{W}$} and attitude as rotation matrix $ R \in SO(3)$ of the MAMR frame {$\mathcal{A}$} relative to {$\mathcal{W}$}. Additionally, the MAMR's linear velocity is denoted by $\textbf{V} \in \mathbb{R}^3$ expressed in {$\mathcal{W}$}. The MAMR's angular velocity expressed in {$\mathcal{W}$} relative to {$\mathcal{A}$} is represented by $\boldsymbol{\Omega} \in \mathbb{R}^3$.

We describe our state in the special euclidean group $SE(3)$, i.e. the group of homogenous transformations with entries in $\mathbb{R}^3$ associated with the Lie algebra, $\mathfrak{se}(3)$ of dimension 6. The main goal is to estimate the homogenous transformation matrix from the MAMR reference frame {$\mathcal{A}$} to the global frame {$\mathcal{W}$} (see Fig. <ref>).

\begin{equation*}
T = ^\mathcal{W}T_\mathcal{A} : \{A\} \rightarrow \{W\}
\end{equation*}

Therefore, the dynamics model can be represented as

\begin{equation} \label{eq:SE3-system}
\dot{T} = T \begin{bmatrix}
(\boldsymbol{\Omega} + \textbf{b} + \boldsymbol{\delta})_\times & \textbf{V} \\
0 & 1
\end{bmatrix}
\end{equation}

with $(\boldsymbol{\Omega} + \textbf{b} + \boldsymbol{\delta})$ the measured angular velocity from the gyroscope including its bias $\textbf{b}$ and noise $\boldsymbol{\delta}$. Additionally, $(\cdot)_\times$ denotes the skew-symmetric matrix associated with the cross product by itself.

Representation of the world reference frame {$\mathcal{W}$} and MAMR reference frame {$\mathcal{A}$}, together with gravity vector $G$ in green, and magnetic field measurements $B_i$ in orange for m EPMs. In purple is the state to estimate.

§.§ Measurement Model

We consider our MAMR to be under $m$ EPMs actuation, and to be fitted with an accelerometer and a 3D HE sensor, providing two types of measurements: acceleration, and magnetic field vector.

Considering that gravitational acceleration ($\textbf{g}$) dominates over linear accelerations as per common approach in literature <cit.>, the accelerometer measurement can be represented as (see Fig. <ref> in green)

\begin{equation}
\textbf{G} = R^T\textbf{g}
\end{equation}

where $R^T$ denotes the transpose of the MAMR's rotation matrix.

The magnetic field vector generated by an EPM$_j$ (with $j = 1, ..., m$) can be modeled as a dipole

\begin{equation} \label{eq:dipole_model}
\textbf{B}_{j} :=  \textbf{B}(\mu_{j},\textbf{r}_{j}) = \frac{\mu_0 |\mu_{j}|}{4 \pi |\textbf{r}_{j}|^3} \left(3 \hat{\textbf{r}}_{j} \hat{\textbf{r}}_{j}^T - I \right)\hat \mu_{j}
\end{equation}

with $\textbf{r}_{j}$ the distance between {$\mathcal{A}$} and EPM$_j$, and $\mu_{j}$ the EPM's magnetic moment in {$\mathcal{A}$}. This assumption is valid for far-enough distances from the EPMs and is commonly employed in other magnetic localization works <cit.>. Assuming that there are no metal objects in the workspace, the measured magnetic field $\textbf{B}$ equals the sum of the magnetic fields generated by each EPM.

\begin{equation} \label{eq:mag_model}
\textbf{B} = \sum_{j = 1}^{m} \textbf{B}_{j}
\end{equation}

Given the absence of the Earth's magnetic field measurement, a minimum of two magnetic field measurements for different configurations of the $m$ EPMs are needed for observability (see Fig. <ref> in orange, and Section <ref>). This is a valid assumption for systems where the magnetic field changes much quicker than the MAMR's pose, such as static or quasi-static systems. This being so, assuming null mean Gaussian measurement noises <cit.>, the measurement model can be expressed as follows. In addition to $n$ measurements of the magnetic field, their norm $\|\textbf{B}_i\|$ was also included. Unlike the full magnetic field measurement vector, which contains information on both position and orientation, the magnetic field norm is dependent only on the MAMR's position. When multiple measurements are present, the addition of the magnetic field norm increased convergence speed.

\begin{equation} \label{eq:meas_model}
\textbf{h} = \begin{bmatrix}
\|\textbf{B}_1\| \\
\vdots \\
\| \textbf{B}_i \| \\
\textbf{B}_1 \\
\vdots \\
\textbf{B}_i \\
\textbf{G}
\end{bmatrix}, \qquad i = 2, \ldots, n
\end{equation}

§.§ Extended Kalman Filter

Extended Kalman Filters (EKF) in $SO(3)$ and $SE(3)$ have been widely used and proved effective <cit.>. For the sake of summary, only the EKF equations are explicitly described here. Further detail on the formulation of EKF can be found in <cit.> and <cit.>.

The discrete dynamics of the estimated state can be described as

\begin{align} \label{eq:dynamics}
\hat{T}_{k+1} = \hat{T}_k \text{exp}(\text{K}_k\tilde{\textbf{y}}_kt) \\
\tilde{\textbf{y}}_k = \textbf{y}_k - \textcolor{black}{\textbf{h}}(\hat{T}_k)
\end{align}

with time-step $k = 0, t, 2t, ...$, $K_k$ the gain defined by the standard EKF prediction and update steps below, $\text{exp}(\cdot)$ the exponential map of $SE(3)$, $\textbf{h}$ the measurement model defined in eq.(<ref>), and $\textbf{y}_k$ the sensors' outputs in the measurement model format, i.e. the norm of the magnetic field, followed by the magnetic field and gravity.

§.§.§ Prediction
This step sees the propagation of the state covariance matrix $P_k \in \mathbb{R}^{6 \times 6}$ as

\begin{equation*}
P_k = F_k\overline{P}_{k-1}F_k^T + G_kQ_nG_k^T
\end{equation*}

with $P_k = \text{diag}(P_{k_{p}}, P_{k_{R}})$, where $P_{k_{p}}$ and $P_{k_{R}}$ denote the state covariance matrix of the position and orientation respectively. Additionally, input noise is considered as a null-mean Gaussian distribution with constant covariance $Q_n \in \mathbb{R}^{6\times6}$. Lastly, $F_k = \text{exp}(A_k t)$ and $G_k = T_k \frac{\partial}{\partial A_k}\text{exp}(A_k t)$, with $A_k$ defined by the Lie algebra as matrix $A_k = [\Omega_\times \quad V; 0 \quad 0]$.

§.§.§ Update
The second step sees the computation of the gain $K_k$ used in the update of the state as shown in eq. (<ref>) through

\begin{align*}
S_k & = H_kP_kH_k^T + R_n \\
K_k & = P_kH_k^TS_k^{-1} \\
\overline{P}_k & = P_k - K_kS_kK_k^T \\
\end{align*}

where $H_k = \frac{\partial \textcolor{black}{\textbf{h}_k}}{\partial T_k} $. Additionally, measurement noise is considered as a null-mean Gaussian distribution with constant covariance matrix $R_n \in \mathbb{R}^{m\times m}$ - $\textcolor{black}{\textbf{h}} \in \mathbb{R}^m$.

§.§ Error metrics:
The observer's performance was assessed through two different error metrics: one for the estimation of the MAMR's position and one for the MAMR's attitude.

e_p = p - p̂

e_R = tr(I - R̂^TR)