# An Exact Poly-Time Membership-Queries Algorithm for Extraction a three-Layer
ReLU Network

Amit Daniely
The Hebrew University and Google Research Tel-Aviv
<EMAIL_ADDRESS>
&Elad Granot
The Hebrew University, Jerusalem
<EMAIL_ADDRESS>

###### Abstract

As machine learning increasingly becomes more prevalent in our everyday life,
many organizations offer neural-networks based services as a black-box. The
reasons for hiding a learning model may vary: e.g., preventing copying of its
behavior or keeping back an adversarial from reverse-engineering its mechanism
and revealing sensitive information about its training data.

However, even as a black-box, some information can still be discovered by
specific queries. In this work, we show a polynomial-time algorithm that uses
a polynomial number of queries to mimic precisely the behavior of a three-
layer neural network that uses ReLU activation.

##  1 Introduction

With the growth of the machine learning field in general and neural networks
model in specific, into new applicative disciplines, many commercial companies
offer machine learning services, allowing public use of pre-trained networks
as a black-box (e.g., IBM Watson [2] or Google AutoML [1]). Those networks
allow the user to query the model, and in many cases, return the exact output
of the network to allow the users to reason about the model output.

However, the reason to maintain the model as a black box is to keep the
intellectual-property of the companies and protect them from others to steal
their models. Moreover, sometimes the training phase uses sensitive data, and
exposing the model might expose this sensitive training data (as demonstrated
in [12]).

However, those models were shown to be sensitive to membership queries attacks
(e.g., [8]). Their work had showed that the architecture of the network could
be extracted under under some basic assumptions over the data. Fornasier et
al. [5] showed a provable algorithm that can extract two-hidden layers model
efficiently as long as the activation function is 3-times differentiable.
However, the widely used ReLU function [7] was left for further investigation.
In this work, we provide a provable algorithm that allows an adversary to
extract the exact parameters of a three-layer network in polynomial time and
polynomial amount of queries with respect to the (unknown) number of
parameters, given that the activation function is the ReLU function. Such
attempts were done before but returned only approximated function (e.g., [11,
9]) or used heuristics that might require an unknown amount of attempts to
succeed (e.g., [4]). The only exact proven algorithm, to our knowledge, was
done by Milli et al. [6] and managed to extract the parameters of a two-layer
network in a polynomial amount of queries, where the hidden layer’s width is
at most the input size and its weights are linearily independent. The idea
behind all those works is to use the fact that a ReLU-activated network is
piecewise-linear, and to inspect the critical points where the linearity
breaks.

This paper continues those line of works by showing:

  1. 1.

a polynomial time and a polynomial query complexities algorithm for a two-
layer neural network with finitely many hidden neurons; and

  2. 2.

a polynomial time and a polynomial query complexities algorithm for a three-
layer neural network with a first-layer right-invertible weight matrix.

The idea behind our algorithm is to obtain a network similar to the original
one, up to an affine transformation. We show that such a network exists and
explain how to find the affine transformations in order to make both models
identical. Furthermore, we keep our queries to the positive orthant, thus
allowing the algorithm to handle a second hidden layer, after peeling the
first layer up to the ReLU function. To our knowledge, this is a first time
where a polynomial-time algorithm is proposed to extract the weights of a
three-layer ReLU model. We hope that this framework would give insights for
further works about model-extraction, allowing to improve the guarantees for
more layers and architectures.

###  1.1 Related Work

Several works proposed to use membership queries to copy the behavior of a
learning model. Some of them (e.g., [11, 9]) suggested to train a new model,
using the black-box model as an annotator. In this line of works, random
samples where chosen and labeled using the black-box model, to create a new
training dataset. Then, an alternative learning model was used on this trainig
data. These works showed that given this technique one can obtain as accurate
prediction as the original model, without the need to collect nor annotate any
real-world data.

To obtain an exact copy of a neural-network, a prior such as the activation
functions was required. In Fornasier et al. [5], it was assumed three-times
differentiable activation functions. Under this assumtion, they approximated
the Hessian of the network on several points and showed those values can
reveal the weights and the activation functions themselves of a three-layer
network.

Unfortunately, this work does not hold for ReLU activation functions, which
are not differentiable. Other works tried to address the case of neural
networks with ReLU functions. All of them used the fact that such networks are
piecewise-linear, where the non-linear parts encode information on the weights
of the network. Rolnick and Kording [8] proved that under mild assumptions,
one can recover the architecture and the weights of any deep ReLU network.
However, their suggested algorithm may suffer from unknown query and time
complexity. In Milli et al. [6], a proven algorithm was able to extract two-
layer ReLU networks in efficient number of queries. They assumed that the rows
of the first weight-matrix are linearly independent, thus limited the possible
width of the network, and far from colinearity. In Carlini et al. [4], it was
suggested to use a heuristically algorithm, that was found in practice to
extract well ReLU networks. However, it was not guaranteed to return the
entire parameters. Moreover, the authors write that their method might require
exponential amount of time, even for the case of a two-layer neural network.
The reason is that in order find the sign of the weights, a brute-force
approach is required.

Our work extends the works above by showing an efficient provable algorithm
that guarantees to extract the weights of a three-layer ReLU neural-network,
using a polynomial amount of queries. Moreover, our algorithm also improves
the known extraction algorithms of two-layer models by allowing to extend the
width of the hidden layer to any number of neurons, while keeping polynomial
time and query complexities with respect to the number of neurons. Our
algorithm contributes to those aspects by, first, showing a simple method to
identify whether a given neuron is in the first layer or not, thus able to
proceed beyond the one-hidden layer case. Second, we consider a neural network
as a sum of a piecewise-linear function and an affine function. This new
definition allows us to use a new approach that finds an isomorphism to the
network without calculating the sign of the weights.

##  2 Results

###  2.1 Depth Two Networks

We consider a $2$-layer network111We assume that our input space is
$\mathbb{R}_{+}^{d}$ since we will need this assumption when we use our result
for depth two networks in order to recover depth three networks. Note that
this assumption is on one hand makes the reconstruction harder. Very similar
result is valid for the case that input space is $\mathbb{R}^{d}$.
$\mathcal{N}:\mathbb{R}_{+}^{d}\rightarrow\mathbb{R}$ as follow:

$\mathcal{N}(\textbf{x})=\sum_{j=1}^{d_{1}}u_{j}\sigma(\textbf{w}_{j}\textbf{x}+b_{j})$
(1)

where $\sigma:\mathbb{R}\rightarrow\mathbb{R}=\max(\cdot,0)$ is the ReLU
function, and for any $j\in[d_{1}]$, $\textbf{w}_{j}\in\mathbb{R}^{d}$,
$b_{j}\in\mathbb{R}$, and222Note that the assumption that the $u_{i}$’s are in
$\\{\pm 1\\}$ does not limit the expressivity of the class of networks that we
consider. Indeed, any neuron $\textbf{x}\mapsto
u\sigma(\textbf{w}\textbf{x}+b)$ calculates the same function as
$\textbf{x}\mapsto\frac{u}{|u|}\sigma(|u|\textbf{w}\textbf{x}+|u|b)$. Thus,
our assumption about the $u_{i}$’s is merely for the sake of standardization.
$u_{j}\in\\{\pm 1\\}$. We assume that the $\textbf{w}_{j}$’s, the $b_{j}$’s
and the $u_{j}$’s, along with the width $d_{1}$, are unknown to the user,
which has only black box access to $\mathcal{N}(\textbf{x})$, for any
$\textbf{x}\in\mathbb{R}_{+}^{d}$. We don’t make any further assumption on the
network weights, rather than a mild general position assumption, which is
outlined in section 3.2.

######  Theorem 1.

There is an algorithm that can reconstruct of a network as defined in(1) using
$\mathrm{poly}(d,d_{1})$ time and queries.

We note that by reconstruction we mean that the algorithm will find weights
$\textbf{w}^{\prime}_{0},\ldots,\textbf{w}^{\prime}_{d_{1}}\in\mathbb{R}^{d}$,
$b^{\prime}_{0},\ldots,b^{\prime}_{d_{1}}\in\mathbb{R}$, and
$u^{\prime}_{1},\ldots,u^{\prime}_{d_{1}}\in\\{\pm 1\\}$ such that

$\forall\textbf{x}\in\mathbb{R}_{+}^{d},\;\;\mathcal{N}(\textbf{x})=\textbf{w}^{\prime}_{0}\textbf{x}+b^{\prime}_{0}+\sum_{j=1}^{d_{1}}u^{\prime}_{j}\sigma(\textbf{w}^{\prime}_{j}\textbf{x}+b^{\prime}_{j}).$
(2)

###  2.2 Depth Three Networks

We consider a $3$-layer network
$\mathcal{N}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ as follow:

$\mathcal{N}(\textbf{x})=\textbf{u}^{T}\sigma(V\sigma(W\textbf{x}+\textbf{b})+\textbf{c})$
(3)

where $W\in\mathbb{R}^{d_{1}\times d}$, $\textbf{b}\in\mathbb{R}^{d_{1}}$,
$V\in\mathbb{R}^{d_{2}\times d_{1}}$, $\textbf{c}\in\mathbb{R}^{d_{2}}$,
$\textbf{u}\in\\{\pm 1\\}^{d_{2}}$ and
$\sigma:\mathbb{R}\rightarrow\mathbb{R}=\max(\cdot,0)$ is the ReLU function,
which is defined element-wise. We assume
$W,V,\textbf{u},\textbf{b},\textbf{c}$, along with $d_{1}$ and $d_{2}$, are
unknown to the user, which have only black box access to $\mathcal{N}(x)$ for
any $\textbf{x}\in\mathbb{R}^{d}$. We denote by $\textbf{w}_{j}$ the $j$th row
of $W$ and assume that it is of unit norm333We note that as the assumption
that $\textbf{u}\in\\{\pm 1\\}^{d_{2}}$, this assumption is merely for the
sake of standartization, and does not restrict the expressivity of the class
of networks that we consider. Indeed, any neuron of the form
$\textbf{x}\mapsto\sigma(\textbf{w}\textbf{x}+b)$ can be replaced by
$\textbf{x}\mapsto\|\textbf{w}\|\sigma\left(\frac{\textbf{w}}{\|\textbf{w}\|}\textbf{x}+\frac{b}{\|\textbf{w}\|}\right)$..

We will make the following assumptions on the network’s weights:

  1. 1.

$W$ is right invertible.

  2. 2.

The top layer has non-zero partial derivatives: For the second layer function
$F:\mathbb{R}^{d_{1}}\to\mathbb{R}$ given by
$F(\textbf{x})=\textbf{u}^{T}\sigma(V\textbf{x}+\textbf{c})$ we assume that
for any $\textbf{x}\in\mathbb{R}_{+}^{d_{1}}$ and $j\in[d_{1}]$, the
derivative of $F$ in the direction of $\mathbf{e}_{j}$ and $-\mathbf{e}_{j}$
is not zero. We note that if $d_{1}$ is large compared to $d$ ($d_{1}\geq 3d$
would be enough) this assumption is valid for most choices of $\textbf{u},V$
and c (see theorem 3). Hence, the non-zero partial derivatives assumption is
valid for a “generic" layer that increases the dimension.

  3. 3.

All the weights in the second layer are non-zero.

  4. 4.

The weights are in general position as detailed in section 3.2

######  Theorem 2.

There is an algorithm that can reconstruct of a network as defined in (3)
using $\mathrm{poly}(d,d_{1},d_{2})$ time and queries.

We note that by reconstruction we mean that the algorithm will find weights
$\textbf{v}^{\prime}_{0},\ldots,\textbf{v}^{\prime}_{d_{2}}\in\mathbb{R}^{d_{1}}$,
$c^{\prime}_{0},\ldots,c^{\prime}_{d_{2}}\in\mathbb{R}$, and
$u^{\prime}_{1},\ldots,u^{\prime}_{d_{2}}\in\\{\pm 1\\}$, as well as a matrix
$W^{\prime}\in\mathbb{R}^{d_{1}\times d}$ and a vector
$\textbf{b}^{\prime}\in\mathbb{R}^{d_{1}}$ such that

$\forall\textbf{x}\in\mathbb{R}^{d},\;\;\mathcal{N}(\textbf{x})=\textbf{v}^{\prime}_{0}\textbf{x}+c^{\prime}_{0}+\sum_{j=1}^{d_{2}}u^{\prime}_{j}\sigma\left(\textbf{v}^{\prime}_{j}\sigma\left(W^{\prime}\textbf{x}+\textbf{b}^{\prime}\right)+c^{\prime}_{j}\right).$

##  3 Proofs

###  3.1 Piecewise linear functions

We denote by $\mathbf{e}_{1},\ldots,\mathbf{e}_{d}$ the standard basis of
$\mathbb{R}^{d}$ and by $\mathbb{B}(\textbf{x},\delta)$ the open ball around
$\textbf{x}\in\mathbb{R}^{d}$ with radius $\delta>0$. Let
$f:\mathbb{R}^{d}\to\mathbb{R}$ be piecewise linear, with finitely many
pieces. A critical point is a point $\textbf{x}\in\mathbb{R}^{d}$ such that
for every $\delta>0$, $f$ is not affine in $\mathbb{B}(\textbf{x},\delta)$. A
critical hyperplane is an affine hyperplane $P$, whose intersection with the
set of critical points is a polytope of dimension $d-1$ exactly. Note that
there are finitely many critical hyperplanes for any piecewise linear
function, that any critical point belongs to at least one critical hyperplane,
and most444By most we mean all, except a set whose dimension is $d-2$.
critical points belongs to exactly one critical hyperplane. We will call such
points non-degenerate. Furthermore, we will say that a critical point x is
$\delta$-non-degenerate if exactly one critical hyperplane intersects
$\mathbb{B}(\textbf{x},\delta)$.

For the function $\mathcal{N}$ computed by a network such as (1) or (3), we
note that for any $j\in[d_{1}]$, the hyperplane
$P_{j}=\\{\textbf{x}:\textbf{w}_{j}\textbf{x}+b_{j}=0\\}$ is a critical
hyperplane. Also, if x is a critical point, then at least one of the neurons
is in a critical state (i.e. its input is $0$). In this case, we will say that
x is a critical point of that neuron.

We next describe a few simple algorithms that we will use frequently.

  * •

Reconstruction of an affine function:  We note that if
$f:\mathbb{B}(\textbf{x},2\delta)\to\mathbb{R}$ is affine, then it can be
reconstructed with $d+1$ queries in $\\\ B(\textbf{x},\delta)$ and $O(d)$
time. Indeed, for any $y\in\mathbb{B}(\textbf{x},\delta)$,

$f(\textbf{x}+\mathbf{y})=f(\textbf{x})+\sum_{i=1}^{d}y_{i}\frac{f(\textbf{x}+\delta\mathbf{e}_{i})-f(\textbf{x})}{\delta}$

Hence, $f$ can be fully reconstructed by querying
$f(\textbf{x}),f(\textbf{x}+\delta\mathbf{e}_{1}),\ldots,f(\textbf{x}+\delta\mathbf{e}_{d})$.

  * •

Reconstruction of critical points in one dimension: Assume that
$f:\mathbb{R}\to\mathbb{R}$ is piece-wise linear, with at most $k$ pieces and
has the following properties: (1) All its critical points are in
$\left(-\frac{1}{\delta},\frac{1}{\delta}\right)$, (2) each piece is of length
at least $\delta$, and (3) there are no two pieces that shares the same affine
function. Then, the critical points of $f$ can be reconstructed precisely
using $O\left(k\log(1/\delta)\right)$ queries and time. Indeed, the following
algorithm recovers the left-most critical point using
$O\left(\log(1/\delta)\right)$ queries and time.

    1. 1.

Set $x_{L}=-\frac{1}{\delta},x_{R}=\frac{1}{\delta}$

    2. 2.

For $j=1,\ldots,\lceil\log_{2}(2/\delta^{2})\rceil+1$:

      1. (a)

If the affine function in the left proximity of $\frac{x_{L}+x_{R}}{2}$ is the
same as the affine function in the proximity of $x_{L}$, set
$x_{L}=\frac{x_{L}+x_{R}}{2}$. Else, set $x_{R}=\frac{x_{L}+x_{R}}{2}$.

    3. 3.

Let $\Lambda_{L}$ and $\Lambda_{R}$ be the affine functions in the left
proximities of $x_{L}$ and $x_{R}$. Return the point $x$ for which
$\Lambda_{L}(x)=\Lambda_{R}(x)$.

To see that the algorithm works, note that if $x^{*}$ is the left-most
critical point then throughout the execution of the algorithm, we have that
$x_{L}\leq x^{*}<x_{R}$, and when the for loop terminates,
$|x_{L}-x_{R}|<\delta$. Hence, in the final step, we have that $x_{L}$ is in
the left-most piece, while $x_{R}$ is in the piece that is adjacent to the
left-most piece.

  * •

Reconstruction of a critical plane:  Let $f:\mathbb{R}^{d}\to\mathbb{R}$ be a
piece-wise linear function. Assume that x is a critical point. Assume
furthermore that $\mathbb{B}(\textbf{x},2\delta)$ intersects only one critical
space $P$. The space $P$ can be reconstructed using $O(d)$ time and queries as
follows: We can choose a random unit vector $\mathbf{e}$, and compute the
tangent spaces $\Lambda_{1},\Lambda_{2}$ of $f$ at
$\textbf{x}+\delta\mathbf{e}$ and $\textbf{x}-\delta\mathbf{e}$. Then, $P$ is
the null space of $\Lambda_{1}-\Lambda_{2}$.

###  3.2 General Position Assumption

In this section we describe our general position assumptions. These
assumptions are quantified by a parameter $\delta$, and the time and query
complexity of our algorithms will depend linearly on $\log(1/\delta)$. While
the description of the conditions is somewhat cumbersome, we argue that these
are very mild conditions. Indeed, if the weights are slightly perturbed, by a
noise of $2^{-d}$ magnitude, the condition will hold with $\delta$ that is
small enough to guarantee polynomial time and query complexity.

  1. 1.

Any critical point in the axes $\\{t\mathbf{e}_{i}:t\in\mathbb{R}\\}$ is
$\delta$-non-degenerate, and corresponds to a single neuron.

  2. 2.

For each $i\in[d]$, the piece-wise linear function
$t\mapsto\mathcal{N}(t\mathbf{e}_{i})$ has the property that there are no two
pieces that share the same affine function.

In the case of depth three network, we assume that the above condition apply
also for the sub-network that is defined by the top two layers.

###  3.3 Reconstruction of Depth Two Network

The strategy of the algorithm is given an access to a network of the form

$\mathcal{N}(\textbf{x})=\textbf{w}_{0}\textbf{x}+b_{0}+\sum_{j=1}^{d_{1}}u_{j}\sigma(\textbf{w}_{j}\textbf{x}+b_{j}),$

to extract one of the non-linear neurons, up to an affine map. Namely, to find
parameters $u,b$ and w such that for some $j$,
$u_{j}\sigma(\textbf{w}_{j}\textbf{x}+b_{j})-u\sigma(\textbf{w}\textbf{x}+b)$
is an affine function. This is enough for the reconstruction of the entire
network, as after finding such a neuron, we can subtract it from the current
network, and create a black-box access to a network with the same set of non-
linear neurons, except the $j$’th neuron. We can then continue and reconstruct
another neuron from the new network, and so on, until we are left with an
affine function, which can be easily reconstructed. We can assume that the
critical hyperplane of any non-linear neuron intersects at least one axis in
the positive orthant, as otherwise this neuron computes a linear function over
$R_{+}^{d}$, and we can therefore drop it and change the liner neuron
accordingly.

Denote the output of the $j$’th neuron before the activation by
$\mathcal{N}_{j}(\textbf{x})=\textbf{w}_{j}\textbf{x}+b_{j}$. Let
$\textbf{x}_{1},\textbf{x}_{2}\in\mathbb{R}_{+}^{d}$ be two points such that
exactly one neuron $j\in[d_{1}]$ changed its state (i.e. changed from active
to inactive, or vice versa) in the segment
$[\textbf{x}_{1},\textbf{x}_{2}]:=\\{\lambda\textbf{x}_{1}+(1-\lambda)\textbf{x}_{2}:\lambda\in[0,1]\\}$.
Moreover, assume that no neuron changes its state in neighborhoods of
$\textbf{x}_{1}$ and $\textbf{x}_{2}$, so that the change in the state happens
in the interior of $[x_{1},x_{2}]$. We note that finding such a pair of points
$\textbf{x}_{1},\textbf{x}_{2}$ can be done by considering a ray
$\ell(\rho)=\rho\mathbf{e}_{i}$, and seeking a critical point $\tilde{\rho}$
of the (one dimensional) function $N\circ\ell$. Under our general position
assumptions, for some $j\in[d]$, there is such a $\rho$ in $\mathbb{R}$, and
it can be found efficiently. Given such a $\rho$, and again under our general
position assumptions, we can take $\textbf{x}_{1}=\ell(\rho-\epsilon)$ and
$\textbf{x}_{2}=\ell(\rho+\epsilon)$, for small enough $\epsilon$.

We will explain next how given such two points, we can reconstruct the $j$’th
neuron, up to an affine function. First, the reconstruction of $u_{j}$ is
simple. Indeed, in the segment $[x_{1},x_{2}]$,
$\mathcal{N}^{\prime}(\textbf{x}):=\mathcal{N}(\textbf{x})-u_{j}\sigma(\textbf{w}_{j}\textbf{x}+b_{j})$
is affine, as no neuron, except the $j$’th neuron, changes its mode. Hence,
$\mathcal{N}(\textbf{x})=u_{j}\sigma(\textbf{w}_{j}\textbf{x}+b_{j})+\mathcal{N}^{\prime}(\textbf{x})$
is a sum of an affine function and the $j$’th neuron. In particular, it is
convex iff the $j$’th neuron is convex iff $u_{j}=1$. Hence, to reconstruct
$u_{j}$ we only need to check if the restriction of $N$ to $[x_{1},x_{2}]$ is
convex or concave.

We next explain how to reconstruct an affine map $\Lambda$ such that
$\sigma(\Lambda(\textbf{x}))-\sigma(\mathcal{N}_{j}(\textbf{x}))$ is affine.
Let $\Lambda_{1},\Lambda_{2}:\mathbb{R}^{d}\to\mathbb{R}$ be the affine maps
computed by the networks in the neighborhoods of $\textbf{x}_{1}$ and
$\textbf{x}_{2}$ respectively. We have that $\Lambda:=\Lambda_{1}-\Lambda_{2}$
is either $N_{j}$ or $-N_{j}$. Hence, we have that
$\sigma(\Lambda(\textbf{x}))$ is either $\sigma(\mathcal{N}_{j}(\textbf{x}))$
or
$\sigma(-\mathcal{N}_{j}(\textbf{x}))=\sigma(\mathcal{N}_{j}(\textbf{x}))-\mathcal{N}_{j}(\textbf{x})$.

We note that the query complexity of this algorithm is
$O\left(dd_{1}\log(1/\delta)\right)$. Indeed, to recover a neuron, we need, at
the worst case, to perform binary search over all the $d$ axes, which costs
$O\left(d\log(1/\delta)\right)$, and reconstructing two affine functions at
the cost of $O(d)$. The run-time of the algorithm is
$O\left(d^{2}d^{2}_{1}\log(1/\delta)\right)$, as each query costs $O(dd_{1})$
(since we need to evaluate the neurons that we subtracted).

###  3.4 Reconstruction of Depth Three Network

####  3.4.1 Extracting a Set that Contains the Critical Spaces of the First
Layer Neurons

The first step of the algorithm would be to find a list

$L=\left\\{\hat{P}_{1},\ldots,\hat{P}_{m}\right\\}$

of affine hyper-planes that contains the critical spaces of all neurons in the
first layer. This is very similar to our algorithm for reconstructing depth
two networks, and we briefly sketch it next. Let
$p_{1},\ldots,p_{m}\in\mathbb{R}^{d}$ be all the critical points of the
restriction of $\mathcal{N}$ to the line
$\\{t\mathbf{e}_{1}:t\in\mathbb{R}\\}$. Note that $m=O(d_{1}d_{2})$ (e.g.
[10]) and these points can be found efficiently. Let
$\hat{P}_{1},\ldots,\hat{P}_{m}$ be the critical hyperspaces corresponding to
these points. The space $\hat{P}_{i}$ can be efficiently found by $O(d)$
queries near $p_{i}$. Finally, since the second layer has non-zero partial
derivatives, for every neuron in the first layer, one of the $p_{i}$’s is in
its critical space, and hence $\hat{P}_{i}$ is its critical space.

####  3.4.2 Identifying First Layer Critical Spaces

The next step is to take the list

$L=\left\\{\hat{P}_{1},\ldots,\hat{P}_{m}\right\\}$

from the previous step, and remove all the hyperplans that are not critical
spaces of neurons in the first layer. After that, the list will contain
precisely the critical spaces of the neurons in the first layer.

In order to do so, it is enough to be able to efficiently decide, given the
list $L$, whether a given hyperplane $P$ is a critical space of a neuron in
the first layer. Suppose that this is not the case. Let
$\textbf{z}_{1},\textbf{z}_{2}\in P$ be two points that are separated by

$P_{1}=\left\\{\textbf{x}:\textbf{w}_{1}\textbf{x}+b_{1}=0\right\\}$

and slightly perturbed inside555Namely,
$\textbf{z}_{i}=\tilde{\textbf{z}}_{i}+\mathbf{e}$, where $\mathbf{e}$ is a
uniform vector in $B(\tilde{\textbf{z}}_{i},\delta)\cap P$. The point
$\tilde{\textbf{z}}_{1},\tilde{\textbf{z}}_{2}$ and $\delta$ are chosen such
that $B(\tilde{\textbf{z}}_{1},\delta)$ and $B(\tilde{\textbf{z}}_{1},\delta)$
are separated by $P_{1}$. $P$. We claim that at least one of
$\textbf{z}_{1},\textbf{z}_{2}$ is not a critical point of $\mathcal{N}$.
Indeed, since $P$ is not a critical space of a neuron in the first layer, and
the points are perturbed inside $P$, they are not critical points of any
neuron in the first layer w.p. 1. As for a neuron in the second layer, say

$h(\textbf{x})=\sum_{i=1}^{d_{1}}v_{i}\sigma(\textbf{w}_{i}\textbf{x}+b_{i})$

we have that the affine functions computed by the network in the proximity of
$\textbf{z}_{1}$ and $\textbf{z}_{2}$ are different. Indeed, these functions
are

$\Lambda_{j}(\textbf{x})=\left(\sum_{i=1}^{d_{1}}v_{i}\sigma^{\prime}(\textbf{w}_{i}\textbf{z}_{j}+b_{i})\textbf{w}_{i}\right)\textbf{x}+\left(\sum_{i=1}^{d_{1}}v_{i}\sigma^{\prime}(\textbf{w}_{i}\textbf{z}_{j}+b_{i})b_{i}\right),\;\;\;j=1,2$

Since $\textbf{z}_{1}$ and $\textbf{z}_{2}$ are on different sides of $P_{1}$,
the vectors

$\left(\sigma^{\prime}(\textbf{w}_{1}\textbf{z}_{j}+b_{1}),\ldots,\sigma^{\prime}(\textbf{w}_{d_{1}}\textbf{z}_{j}+b_{d_{1}})\right),\;\;\;j=1,2$

are different. Hence, since the $\textbf{w}_{1},\ldots,\textbf{w}_{d_{1}}$ are
independent (since $W$ is assumed to be right invertible) and
$v_{1},\ldots,v_{d_{1}}$ are non-zero we have that

$\left(\sum_{i=1}^{d_{1}}v_{i}\sigma^{\prime}(\textbf{w}_{i}\textbf{z}_{1}+b_{i})\textbf{w}_{i}\right)\neq\left(\sum_{i=1}^{d_{1}}v_{i}\sigma^{\prime}(\textbf{w}_{i}\textbf{z}_{2}+b_{i})\textbf{w}_{i}\right)$

and in particular $\Lambda_{1}\neq\Lambda_{2}$. It follows that the null-space
of at least one of these affine functions, say $\Lambda_{1}$, is not $P$.
Since $\textbf{z}_{1}$ is perturbed inside $P$, the probability that
$\textbf{z}_{1}$ is critical point of $h$ (i.e.
$h(\textbf{z}_{1})=\Lambda_{1}(\textbf{z}_{1})=0$) is $0$.

By the above discussion, together with the fact that if $P$ is a critical
space, then both $\textbf{z}_{1},\textbf{z}_{2}$ are critical point of
$\mathcal{N}$, we conclude that the following algorithm can verify whether $P$
is a critical space:

  1. 1.

For any $k\in[m]$, such that $\hat{P}_{k}\neq P$:

    1. (a)

Choose two perturbed points $\textbf{z}_{1},\textbf{z}_{2}$ in $P$ that are
separated by $\hat{P}_{k}$

    2. (b)

If one of these points is not critical, return "$P$ is not critical"

  2. 2.

Return "$P$ is critical"

####  3.4.3 Identifying Directions

Since the rows in $W$ assumed to have unit norm, the list of the critical
spaces of the neurons in the first layer, obtained in the previous step,
determines the weights up to sign. Namely, we can reconstruct a list

$L=\left\\{(\hat{\textbf{w}}_{1},\hat{b}_{1}),\ldots,(\hat{\textbf{w}}_{d_{1}}\hat{b}_{d_{1}})\right\\}$

that define precisely the neurons on the first layer, up so sign. It therefore
remains to recover the correct signs.

In order to recover the correct sign of $(\hat{\textbf{w}}_{1},\hat{b}_{1})$,
we can simply do the following test: Choose a point x such that
$\hat{\textbf{w}}_{1}\textbf{x}+\hat{b}_{1}=0$, and query the network in the
points $\textbf{x}+\epsilon\textbf{z},\textbf{x}-\epsilon\textbf{z}$, for
small $\epsilon$, where z is a unit vector that has the property that is
orthogonal to $\hat{\textbf{w}}_{2},\ldots,\hat{\textbf{w}}_{d_{2}}$ but
$\hat{\textbf{w}}_{1}\textbf{z}>0$. Such a z exists since $W$ is right
invertible and hence $\textbf{w}_{1},\ldots,\textbf{w}_{d_{2}}$ are linearly
independent.

Now, when moving from x to either $\textbf{x}+\epsilon\textbf{z}$ or
$\textbf{x}-\epsilon\textbf{z}$, the value of all the neurons in the first
layer, possibly except the one that corresponds to
$(\hat{\textbf{w}}_{1},\hat{b}_{1})$, does not change. As for one that
corresponds to $(\hat{\textbf{w}}_{1},\hat{b}_{1})$, if its weights are
$(\hat{\textbf{w}}_{1},\hat{b}_{1})$, then its value changes when we move from
x to $\textbf{x}+\epsilon\textbf{z}$ but not when we move from x to
$\textbf{x}-\epsilon\textbf{z}$. On the other hand, if its weights are
$(-\hat{\textbf{w}}_{1},-\hat{b}_{1})$, then the value changes when we move
from x to $\textbf{x}-\epsilon\textbf{z}$ but not when we move from x to
$\textbf{x}+\epsilon\textbf{z}$.

Hence, in order to have the correct sign we can do the following: If
$\mathcal{N}(\textbf{x})\neq\mathcal{N}(\textbf{x}+\epsilon\textbf{z})$ then
keep $(\hat{\textbf{w}}_{1},\hat{b}_{1})$. Else, replace it with
$(-\hat{\textbf{w}}_{1},-\hat{b}_{1})$. The reason that this test work is the
above discussion, together with the assumption that the second layer has non-
zero partial derivatives.

####  3.4.4 Reconstruction of the top two layers

After having the weights of the first layer at hand, and since $W$ is assumed
to be right invertible, we can access the sub-network defined by the top two
layers directly. Namely, given $\textbf{x}\in\mathbb{R}_{+}^{d_{1}}$, we can
find $\textbf{z}\in\mathbb{R}^{d}$ such that

$\textbf{x}=\left(\sigma\left(\textbf{w}_{1}\textbf{z}+b_{1}\right),\ldots,\sigma\left(\textbf{w}_{d_{1}}\textbf{z}+b_{d_{1}}\right)\right)$

e.g., by taking $\textbf{z}=W^{-1}(\textbf{x}-\textbf{b})$ where $W^{-1}$ is
the right inverse of $W$. Now, $\mathcal{N}(\textbf{z})$ is exactly the value
of the top layer on the input x. Hence, the problem of reconstructing the top
two layers boils down to the problem of reconstructing a depth two network,
which we already solved.

We note that the query complexity of this algorithm is
$O\left(dd_{1}d_{2}\log(1/\delta)+d_{1}^{2}d_{2}^{2}\right)$. Indeed, as the
binary search of the first layer finds $m=O(d_{1}d_{2})$ points, on which the
identification of their layer is quadratic in the query complexity, we need,
at the worst case, $O\left(dd_{1}d_{2}\log(1/\delta)\right)$ queries for the
search and reconstruction, and additional $O\left(d_{1}^{2}d_{2}^{2}\right)$
queries for the identification. The run-time of the algorithm is
$O\left(d^{2}d_{1}d_{2}(d_{1}+d_{2})\log(1/\delta)\right)$, as each query
costs $O(d\max(d_{1},d_{2}))$ to evaluate. Finally, as we note that a slight
perturbation of the weights, by a noise of $2^{-d}$ magnitude, we get that
$\frac{1}{\delta}$ is polynomial in $d$ with probability $1$, then we can
guarantee polynomial time and query complexity in a general position case.

###  3.5 On the non-zero partial derivatives assumption

Consider a ReLU network

$\mathcal{N}(\textbf{x})=\sum_{j=1}^{d_{1}}u_{j}\sigma(\textbf{w}_{j}\textbf{x}+b_{j})$
(4)

and assume that for any $j$, $u_{j}\neq 0$ (otherwise the corresponding neuron
can be dropped). We have that

$\frac{\partial\mathcal{N}}{\partial
x_{i}}(\textbf{x})=\sum_{j=1}^{d_{1}}u_{j}\sigma^{\prime}(\textbf{w}_{j}\textbf{x}+b_{j})\textbf{w}_{j}(i)$

Now, if the weights are random, say that the $\textbf{w}_{j}$’s are
independent random variables such that $\textbf{w}_{j}(i)$ has a continuous
distribution, then w.p. 1 we have that for every non-zero vector,
$\textbf{z}\in\\{0,1\\}^{d_{1}}$ it holds that

$\sum_{j=1}^{d_{1}}z_{j}\textbf{w}_{j}(i)\neq 0$

and hence $\frac{\partial\mathcal{N}}{\partial x_{i}}(\textbf{x})\neq 0$,
unless the vector
$\Lambda(\textbf{x}):=\left(\textbf{w}_{1}\textbf{x}+b_{1},\ldots\textbf{w}_{d_{1}}\textbf{x}+b_{d_{1}}\right)$
is in the negative orthant $\mathbb{R}_{-}^{d_{1}}$. It follows the non-zero
partial derivatives assumption holds, provided if and only if the affine map
$\Lambda$ maps the positive orthant $\mathbb{R}_{+}^{d}$ to the complement of
the negative orthant $\mathbb{R}_{-}^{d_{1}}$. The following lemma shows that
if $d_{1}\gg d$, then this is often the case.

######  Lemma 1.

Assume that the pairs $(\textbf{w}_{j},b_{j})$ are independent and
symmetric666That is, for all $j\in[d_{1}]$, the distributions of
$(\textbf{w}_{j},b_{j})$ and $(-\textbf{w}_{j},-b_{j})$ are the same., then

$\Pr\left(\Lambda(\mathbb{R}^{d})\cap\mathbb{R}^{d_{1}}_{-}\neq\emptyset\right)\leq\frac{\left(\frac{ed_{1}}{d}\right)^{d+1}}{2^{d_{1}}}$

###### Proof.

We first note that the number of orthants that has non-negative intersection
with $\Lambda(\mathbb{R}^{d})$ is exactly the number of functions in the class

$H=\left\\{j\in[d_{1}]\mapsto\mathrm{sign}(\textbf{z}(j)):\textbf{z}\in\Lambda(\mathbb{R}^{d})\right\\}$

Since $\Lambda(\mathbb{R}^{d})$ is an affine space of dimension at most $d$,
$H$ has VC dimension at most $d+1$ (e.g. [3]). Hence, by the Sauer-Shelah
lemma (again, [3])

$|H|\leq\sum_{i=0}^{d+1}\binom{d_{1}}{i}\leq\left(\frac{ed_{1}}{d}\right)^{d+1}$

Finally, since the $(\textbf{w}_{j},b_{j})$’s and symmetric, the probability
that $\Lambda(\mathbb{R}^{d})$ intersects $\mathbb{R}_{-}^{d_{1}}$ is the same
as the probability that it intersects any other orthant. Since there are
$2^{d_{1}}$ orthants, and $\Lambda(\mathbb{R}^{d})$ intersects at most
$\left(\frac{ed_{1}}{d}\right)^{d+1}$ of them, it follows that
$\Pr\left(\Lambda(\mathbb{R}^{d})\cap\mathbb{R}^{d_{1}}_{-}\neq\emptyset\right)\leq\frac{\left(\frac{ed_{1}}{d}\right)^{d+1}}{2^{d_{1}}}$.
∎

All in all we have the following theorem

######  Theorem 3.

Assume that the pairs $(\textbf{w}_{j},b_{j})$ are independent, symmetric, and
has continuous marginals, w.p.
$1-\frac{\left(\frac{ed_{1}}{d}\right)^{d+1}}{2^{d_{1}}}$ we have that
$\frac{\partial\mathcal{N}}{\partial x_{i}}(\textbf{x})\neq 0$ for all
$\textbf{x}\in\mathbb{R}^{d}$ and $i\in[d]$.

##  4 Discussion and Future Work

This work joins a set of empirical and theoretical works showing that the
extraction of a ReLU network given membership-queries is possible. Here we
provide a constructive proof that three-layers model extraction can be done in
polynomial time. Our nonrestrictive assumptions make it feasible to construct
a fully connected, a convolutional, and many other architectures.

As empirical evidence show, we believe that it is possible to prove similar
results with even fewer assumptions for deeper models and more complex
architectures. We leave this challenge for future works.

## References

  * [1] Google Cloud AutoML.  URL https://cloud.google.com/automl.
  * [2] IBM Watson.  URL https://www.ibm.com/watson.
  * Anthony and Bartlet [1999] M. Anthony and P. Bartlet.  _Neural Network Learning: Theoretical Foundations_.  Cambridge University Press, 1999.
  * Carlini et al. [2020] Nicholas Carlini, Matthew Jagielski, and Ilya Mironov.  Cryptanalytic Extraction of Neural Network Models.  _arXiv preprint_ , (arXiv: 2003.04884), 2020.
  * Fornasier et al. [2019] Massimo Fornasier, Timo Klock, and Michael Rauchensteiner.  Robust and Resource Efficient Identification of Two Hidden Layer Neural Networks.  _arXiv preprint_ , (arXiv: 1907.00485), 2019.
  * Milli et al. [2019] Smitha Milli, Ludwig Schmidt, Anca D. Dragan, and Moritz Hardt.  Model Reconstruction from Model Explanations.  In _Proceedings of the Conference on Fairness, Accountability, and Transparency - FAT* ’19_ , pages 1–9, New York, New York, USA, jul 2019. ACM Press.  ISBN 9781450361255.  doi: 10.1145/3287560.3287562.
  * Nair and Hinton [2010] Vinod Nair and Geoffrey E Hinton.  Rectified linear units improve Restricted Boltzmann machines.  In _ICML 2010 - Proceedings, 27th International Conference on Machine Learning_ , pages 807–814, 2010.  ISBN 9781605589077.
  * Rolnick and Kording [2019] David Rolnick and Konrad P. Kording.  Reverse-Engineering Deep ReLU Networks.  _arXiv preprint_ , (arXiv: 1910.00744), oct 2019.
  * Shi et al. [2017] Yi Shi, Yalin Sagduyu, and Alexander Grushin.  How to steal a machine learning classifier with deep learning.  _2017 IEEE International Symposium on Technologies for Homeland Security, HST 2017_ , (April 2017), 2017.  doi: 10.1109/THS.2017.7943475.
  * Telgarsky [2016] Matus Telgarsky.  benefits of depth in neural networks.  In _Conference on Learning Theory_ , pages 1517–1539, 2016.
  * Tramer et al. [2016] Florian Tramer, Fan Zhang, Ari Juels, Michael K. Reiter, and Thomas Ristenpart.  Stealing Machine Learning Models via Prediction APIs.  _USENIX Security Symposium_ , 2016.  ISSN 00414131.
  * Zhang et al. [2019] Yuheng Zhang, Ruoxi Jia, Hengzhi Pei, Wenxiao Wang, Bo Li, and Dawn Song.  The Secret Revealer: Generative Model-Inversion Attacks Against Deep Neural Networks.  _arXiv preprint_ , (arXiv: 1911.07135), nov 2019.


# Dual Particle Output Feedback Control based on Lyapunov drifts for nonlinear
systems

Emilien Flayac, Karim Dahia, Bruno Hérissé, and Frédéric Jean

###### Abstract

This paper presents a dual receding horizon output feedback controller for a
general non linear stochastic system with imperfect information. The novelty
of this controller is that stabilization is treated, inside the optimization
problem, as a negative drift constraint on the control that is taken from the
theory of stability of Markov chains. The dual effect is then created by
maximizing information over the stabilizing controls which makes the global
algorithm easier to tune than our previous algorithm. We use a particle filter
for state estimation to handle nonlinearities and multimodality. The
performance of this method is demonstrated on the challenging problem of
terrain aided navigation.

##  I INTRODUCTION

Stochastic model predictive control (SMPC) is a widespread technique to deal
with control problems where the state is subject to stochastic disturbances,
hard constraints on the input and potentially soft constraints on the state.
Its principle is to apply a receding horizon strategy based on the resolution
of a finite horizon stochastic optimal control problem. Although, when only
partial and noisy information on the state is available through some
observations, classic SMPC combined with a state estimator may lead to
overcautious controls or even destabilizing ones. It is due to the fact that,
in general, the control influences the observations and then state estimation
in addition to guiding the system in a standard way. This property is known as
the dual effect property of the control. Stochastic optimal control problems
with imperfect information are much harder than their counterpart in the
perfect information case. The main reason is that to be optimal with partial
information, one needs to anticipate the information that will be available.
Consequently, the optimal law possesses the property of dual effect. As it is
computationally intractable, suboptimal outputfeedback control laws are
computed instead, with the idea to keep the dual effect property. Therefore,
they are called dual controllers.

Usually, dual controllers are either implicit or explicit. The design of
implicit dual controls, which we do not address in this paper, is based on the
idea to approximate the Bellman equation of the problem. See [1] for a review
on implicit dual control and [2] for an example of use of a particle filter in
implicit dual control. On the contrary, in explicit dual control, one modifies
the original problem to incorporate an explicit excitation of the system to
maintain the dual effect. Explicit dual effect can be included in the
optimization problem as a constraint or in the cost. Including it as a
constraint generally leads to controllers with persistent excitation,
constraints on information or approximation with scenario trees. Including it
in the cost leads to integrated experiment design, where a measure of
information is added to the original cost. In [3], a persistent excitation
controller is presented but it is supposed that the system is linearizable
which is not a suitable assumptions in our framework. Controllers with
constraints on the information may lead to infeasibilty problems that are also
hard to anticipate in nonlinear cases. In [4], this infeasibility issue is
addressed in the deterministic framework. Scenario tree methods are
computationally demanding and usually rely on Kalman filters inside the
optimization problem, like in [5], [6] which are not adapted to multimodal
cases.

In integrated experiment design, a new term based on the Fisher information
matrix (FIM) is added to the original cost as classically done in optimal
design [7]. However, if the original cost to minimize conditions the stability
of the system, then the resulting trade off may destabilize it.

That is why, in this paper, we propose an output feedback dual SMPC based on
integrated experiment design for information probing and on a Lyapunov
constraint on the control to ensure stability. Output feedback is obtained by
setting the initial condition of the optimization problem as particles from a
particle filter that is able to deal with state estimation for arbitrary
systems.

The paper is organized as follows. Section II recalls some important
probabilistic notations. Section III recalls the bases of SMPC with imperfect
information and describes our approach. Section IV recalls the principle of
particle filtering and section V describes the particular optimization problem
that is solved in a receding horizon manner. In [8] and [9], stochastic output
feedback stability is proved but a separation principle was assumed. Actually,
proving stochastic output feedback stability in the general case, when
observability depends on the control, and when the only efficient estimators
are particle filters, has never been proved to the best of our knowledge.
However, we demonstrate the efficiency of our method on a challenging
application in section VI.

##  II NOTATION

Let $(\Omega,\mathcal{F},P)$ be a probability space. In the following, random
variables refer to $\mathcal{F}$-measurable functions defined on $\Omega$. For
$i\in\mathbb{N}$, $\mathcal{B}(\mathbb{R}^{i})$ denotes the set of Borel sets
of $\mathbb{R}^{i}$. For a random variable X and a probability distribution,
$X\sim p$ means that $p$ is the probability law of $X$. $P(\cdot|\cdot)$ and
$E(\cdot|\cdot)$ denotes conditional probability and expectation. $P_{x}$ and
$E_{x}$ denotes especially probability and expectation conditionally to
$X_{0}=x$, for $x\in\mathbb{R}^{i}$. $P_{p}$ and $E_{p}$ denotes probability
and expectation conditionally to $X_{0}\sim p$.

##  III PROBLEM SETUP

###  III-A Description of the system

We consider a discrete-time controlled stochastic dynamical system
${(X_{k})}_{k\in\mathbb{N}}$ valued in $\mathbb{R}^{n_{x}}$ described by the
following equation, $\forall k\in\mathbb{N}$:

$\displaystyle X_{k+1}=f(X_{k},U_{k},{\xi}_{k}),$ $\displaystyle X_{0}\sim
p_{0},$ (1)

where:

  * •

$p_{0}$ is the initial probability law.

  * •

${(U_{k})}_{k\in\mathbb{N}}$ is the control process valued in
$\mathbb{R}^{n_{u}}$.

  * •

${(\xi_{k})}_{k\in\mathbb{N}}$ are i.i.d. random variables valued in
$\mathbb{R}^{n_{\xi}}$ distributed according to $p_{\xi}$. For each
$k\in\mathbb{N}$, $\xi_{k}$ represents an external disturbance on the
dynamics.

  * •

f :
$\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\xi}}\longrightarrow\mathbb{R}^{n_{x}}$
is measurable.

We also assume that the state of the system is only available through some
observations represented by a stochastic process ${(Y_{k})}_{k\in\mathbb{N}}$
valued in $\mathbb{R}^{n_{y}}$ which verifies, $\forall k\in\mathbb{N}$:

$\displaystyle Y_{k}=h(X_{k},\eta_{k}),$ (2)

where:

  * •

${(\eta_{k})}_{k\in\mathbb{N}}$ are i.i.d. random variables valued in
$\mathbb{R}^{n_{\eta}}$ distributed according to $p_{\eta}$. For each
$k\in\mathbb{N}$, $\eta_{k}$ represents an external disturbance on the
observations.

  * •

$h:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{\eta}}\longrightarrow\mathbb{R}^{n_{y}}$
is measurable.

For $k\in\mathbb{N}$, we define the information vector $I_{k}$ as follow:

$\displaystyle I_{k}$
$\displaystyle=(Y_{0},U_{0},\dots,Y_{k-1},U_{k-1},Y_{k}),$ (3)

From $I_{k}$, one can derive two important quantities in stochastic control
with imperfect information that are the conditional distribution of $X_{k}$
given $I_{k}$ denoted by $\mu_{k}$ and the conditional distribution of
$X_{k+i}$ given $(I_{k},U_{k},\dots,U_{k+i-1})$ for any $i\in\mathbb{N}^{*}$,
denoted by $\mu_{k+i|k}$. Moreover, we denote by $K$ the Markov kernel defined
by equation (1) and we assume that the conditional distribution defined by
equation (2) has a density with respect to the Lebesgue measure such that
there exists a likelihood function denoted by $\rho$. Therefore, for
$k\in\mathbb{N}$, $A\in\mathcal{B}(\mathbb{R}^{n_{x}})$ and
$B\in\mathcal{B}(\mathbb{R}^{n_{y}})$ :

$\displaystyle P(X_{k+1}\in A|X_{k}=x_{k},U_{k}=u_{k})$
$\displaystyle=K(A,x_{k},u_{k}),$ $\displaystyle P(Y_{k}\in B|X_{k}=x_{k})$
$\displaystyle=\int_{B}\rho(y_{k},x_{k})\mathrm{d}y_{k}.$

$\mu_{0}$ is supposed to be known. Thus, $\forall(k,i)\in\mathbb{N}^{2}$,
$\mu_{k}$ and $\mu_{k+i|k}$ verify the following nonlinear filtering equations
that can be summed up to:

$\displaystyle\mu_{k+1}$ $\displaystyle=F\left(\mu_{k},Y_{k+1},U_{k}\right),$
(4) $\displaystyle\mu_{k+i+1|k}$
$\displaystyle=G\left(\mu_{k+i|k},U_{k+i}\right).$ (5)

The expressions of $F$ and $G$ are classic and can be found in [1]. A control
policy at time $k$ is then a map, denoted by $\pi_{k}$, that maps a
conditional distribution $\mu_{k}$ to a control $U_{k}$. We denote a sequence
of control policies by:

$\displaystyle\pi_{i:j}:=(\pi_{i},\dots,\pi_{j}),\;\text{for}\;i\leq j.$

###  III-B Dual stochastic model predictive control

Stochastic model predictive control is a widely used method for designing
controllers in the presence of possibly unbounded disturbances in nonlinear
dynamics as the one described in equation (1). It consists in solving a
finite-horizon discrete stochastic optimal problem, to only apply the first
control policy of the optimal sequence and to solve the problem again starting
from the new state of the system. See [10] for a review on general stochastic
MPC.

A MPC scheme is defined by the following features:

  * •

A time horizon, $T\in\mathbb{N}^{*}$.

  * •

A family of set of constraints on the control,
$\mathbb{U}_{i}\subset\mathbb{R}^{n_{u}}$, $\forall i=0,..,T-1$.

  * •

$\forall i=0,..,T-1$, an instantaneous cost $g_{i}$:
$\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\xi}}\longrightarrow\mathbb{R}$
and a final cost $g_{T}$: $\mathbb{R}^{n_{x}}\longrightarrow\mathbb{R}$.

The choice of the objective functions $g_{i}$ and $g_{T}$ and of the control
constraints $\mathbb{U}_{i}$ is a matter of modelling. Indeed, one often has
an economic cost to minimize that comes from practical considerations, such as
a price or a fuel consumption. Moreover, the control must be designed to
attain some target in the state space: it is the guiding problem. There are
two classical ways to address this issue in MPC:

  * •

Adding a new term to the cost to enforce stability in some sens. In this case,
the general cost is decomposed in the following way, for $i=0,..,T$:

$\displaystyle g_{i}=g_{i}^{\textrm{stab}}+g_{i}^{\textrm{eco}},$ (6)

Thus, $g_{i}$ and $g_{T}$ realize a compromise between convergence and
economic costs. For instance, in the LQG case,
$g_{i}(x,u,\xi)=x^{T}M_{x}x+u^{T}M_{u}u$ where $M_{x}$ and $M_{u}$ are
positive definite matrices. The first term drives the state of the system to
zero and the second term penalizes high controls. The compromise is dealt with
by tuning the matrices $M_{x}$ and $M_{u}$.

  * •

Adding a drift constraint on the first control, $U_{0}$, that enforces the
decreasing of some Lyapunov-like function, during the first time step only,
such that, for $i=0,..,T$:

$\displaystyle g_{i}=g_{i}^{\textrm{eco}},$ (7) $\displaystyle\text{a negative
drift condition on }U_{0}.\;$

Actually, since only $U_{0}$ is applied on the system, the Lyapunov function
decreases along the whole trajectory and then stability is obtained. It is
also known as Lyapunov Economic MPC, see [11] for a review in the
deterministic setting. In the stochastic setting, it has been applied with
output feedback for continuous-time nonlinear systems in [9] and for a
discrete-time linear system with bounded controls in [8].

In the presence of partial information represented by equation (2), the
theoretical receding horizon problem to solve is more complex than in the
perfect information case, as the the new state of the system is $\mu_{k}$
[12], and it evolves according to equation (4). Thus, the problem to solve at
each time k denoted by $(P^{k}_{C})$ can be written as follows:

$\begin{array}[]{rrclcc}\displaystyle\underset{\pi_{0:T-1}}{\text{min}}&\lx@intercol
E_{\mu_{k}}\left[\sum_{i=0}^{T-1}E[g_{i}(X_{i},U_{i},\xi_{i})|I_{i}]+E[g_{T}(X_{T})|I_{T}]\right]\hfil\lx@intercol\\\
\textrm{s.t.}&\tilde{\mu}_{i+1}&=&F\left(\tilde{\mu}_{i},Y_{i+1},U_{i}\right),\\\
&U_{i}&=&\pi_{i}(\tilde{\mu}_{i}),\\\ &U_{i}&\in&\mathbb{U}_{i},\;\forall
i=0,..,T-1,\\\ &\tilde{\mu}_{0}&=&\mu_{k}.\end{array}$

In $(P^{k}_{C})$, the conditional expectation w.r.t. the available information
$I_{i}$ emphasizes the fact that the observations, are propagated from $k$ to
$k+T$ and depend on the control. Therefore, as recalled in [1] and [13], one
can see from the Dynamic Programming Principle that the optimal control
policies of $(P^{k}_{C})$ generally have the dual effect property [14]. As
optimal policies are intractable in the general nonlinear case due to the
curse of dimensionality, suboptimal policies that are constructed instead.
They are called dual controllers as they preserve the dual effect property by:

  * •

controlling the system in the standard way defined by equation (6) or (7).

  * •

explicitely or implicitely probing information to improve the quality of the
future observations.

In this article, we focus on explicit dual controllers that take information
into account inside the cost. This technique is called integrated experiment
design. In fact, the new cost denoted by $g_{i}^{\textrm{ex}}$ realizes a
trade off between the original costs $g_{i}$ and $g_{T}$ and a measure of
information denoted by $g_{i}^{\textrm{info}}$ such that, for $i=0,..,T$:

$\displaystyle g_{i}^{\textrm{ex}}=g_{i}+g_{i}^{\textrm{info}},$ (8)

An important issue appears when, additionally, the guiding objective is
incorporated in the cost as described in equation (6). Consequently, (8)
represents the sum of three terms that usually contain weights that must be
tuned. By construction, the weights affect the convergence of the system to
the guiding goal which makes their tuning difficult. Indeed, if too much
importance is given to $g_{i}^{\textrm{info}}$ then probing will be efficient
but the system may get stuck far from the target. On the contrary, if too much
importance is given to $g_{i}^{\textrm{stab}}$ then the probing effect will
not be sufficient and output feedback performance may be poor. It is difficult
to know which case will occur a priori depending of the value of the weights
because the optimal costs are not known at first sight.

That is why, in this paper, we present an output-feedback explicit dual
controller based on the minimization of the Fisher information matrix subject
to a negative mean drift condition for nonlinear discrete-time systems coupled
with a particle filter for state estimation. In other words, we have chosen to
model our problem with equations (7) and (8).

Furthermore, in explicit dual control, one does not have to propagate the
information inside the dynamics with equation (4) anymore because it is dealt
with in a different way. It is sufficient to propagate $\mu_{i|0}$ with
equation (5) which is much simpler. Besides, one wants to solve in practice a
finite dimensional optimization problem so one looks for control values and
not policies anymore. Then, the general optimization problem associated with
our method, denoted by $(P^{k}_{\textrm{ex}})$, can be written as follows:

$\begin{array}[]{rrclcc}\displaystyle\underset{u_{0:T-1}}{\text{min}}&\lx@intercol
E_{\mu_{k}}\left[E[\sum_{i=0}^{T-1}g_{i}^{\textrm{ex}}(X_{i},u_{i},\xi_{i})+g_{T}^{\textrm{ex}}(X_{T})|I_{0}]\right]\hfil\lx@intercol\\\
\textrm{s.t.}&\tilde{\mu}_{i+1|0}&=&G\left(\tilde{\mu}_{i|0},u_{i}\right),\\\
&u_{i}&\in&\mathbb{U}_{i},\;\forall i=0,..,T-1,\\\
&\tilde{\mu}_{0}&=&\mu_{k}.\\\ &\text{negative}&\text{drift}&\text{condition
on}\;u_{0}.\\\ \end{array}$

##  IV STATE ESTIMATION

From the definition of $(P^{k}_{\textrm{ex}})$, it is clear that dual SMPC
requires the computation of $\mu_{k}$. However, in the general nonlinear case,
$\mu_{k}$ cannot be computed explicitly so approximation methods are needed.
Kalman filters are widespread and easy-to-compute approximations of the
posterior distribution but they may fail in the presence of high
nonlinearities and multimodality. That is why we use particle filters that are
known to handle these difficulties at the price of a higher computational
cost. A particle filter approximates the posterior distribution $\mu_{k}$ by a
set of N particles, ${\left(x^{l}_{k}\right)}_{l=1,..,N}$ valued in
$\mathbb{R}^{n_{x}}$, associated with nonnegative and normalized weights
${\left({\omega}^{l}_{k}\right)}_{l=1,..,N}$. The approximate distribution,
denoted by, $\mu_{k}^{N}$, is then defined as follows:

$\displaystyle\mu_{k}^{N}=\textstyle\sum_{l=1}^{N}\omega_{k}^{l}\delta_{x^{l}_{k}},$
(9)

Algorithm 1 described the steps of a particle filtering algorithm with
adaptive resampling. Furthermore, we define the conditional expectation of the
state, denoted by $\hat{x}_{k}$ w.r.t. the information $I_{k}$ and its
particle approximation, the empirical mean of the filter $\hat{x}_{k}^{N}$
i.e:

$\displaystyle\hat{x}_{k}=E[X_{k}|I_{k}],$
$\displaystyle\hat{x}_{k}^{N}=\frac{1}{N}\textstyle\sum_{l=1}^{N}\omega^{l}_{k}x^{l}_{k}.$

##  V CONTROL POLICY

The goal of this section is to present our new explicit dual receding horizon
control scheme. The main idea is to prioritize the guiding goal by adding a
stabilizing constraint on the first control of the finite horizon optimization
problem. This constraint consists in forcing the decreasing of a Lyapunov-like
function during the first time step starting from the current estimator the
state. It is equivalent to stabilizing an estimator of the state. The success
of such a technique depends highly on the estimation error. Indeed, if the
estimation error is high or even diverges, driving the estimator to the target
does not imply driving the true state to the target. It is a well known
problem in the deterministic setting [15]. In [9] and [8], output feedback
stability is proved but it was assumed that the estimation error converges
uniformly w.r.t. the control. This assumption is not realistic anymore in our
setting because the observations do depend on the control and the estimation
error may diverge if the dual effect is not taken into account. Accordingly,
our method consists in solving a version of $(P^{k}_{\textrm{ex}})$ where
$g_{i}^{\textrm{ex}}$ is based on the Fisher information matrix and the drift
condition comes from the theory of Markov chains stability.

###  V-A Foster-Lyapunov drift in case of perfect information

Perfect information is met when full knowledge of the state $X_{k}$ is
available. In the time homogeneous case, the control is computed from a state
feedback control policy i.e. a measurable maps, denoted by $\alpha$, that maps
a state $X_{k}$ to a control $U_{k}$. Thus, for any $\alpha$, one can define
the corresponding closed loop system as follows, $\forall k\in\mathbb{N}$:

$\displaystyle X_{k+1}=f(X_{k},\alpha(X_{k}),{\xi}_{k}),$ $\displaystyle
X_{0}\sim p_{0}.$ (10)

A state feedback control policy is said to be admissible if $\forall
x\in\mathbb{R}^{n},\alpha(x)\in\mathbb{U}_{0}$. Therefore, equation (10)
defines also a time-homogeneous Markov chain whose stability can be studied
via the classical theory of negative drifts conditions discussed in [16] and
recalled in [17]. In proposition 1, we focus on geometric drifts conditions
that are closely related to Lyapunov conditions for exponential stability for
continuous-time processes.

Algorithm 1 Particle filter with adaptive resampling

1:Create a sample of $N$ particles $x^{l}_{0}$ according to the law $\mu_{0}$
and initialize the weights $\omega_{0}^{l}$ with $\frac{1}{N}$

2:for $i=0,1,2\dots$ do

3: Prediction:

4: Given a control $u_{i}$ and the particles
${\left(x^{l}_{i}\right)}_{l=1,..,N}$, compute the predicted particles by
drawing samples from $K$ i.e.  $\displaystyle x^{l}_{i+1|i}\sim
K(\mathrm{d}x_{i+1|i},x_{i}^{l},u_{i}),\;\text{for}\;l=1,..,N.$

5: Correction:

6: Get the new observation $Y_{i+1}$

7: Compute the updated weights
${\left(\tilde{\omega}^{l}_{i}\right)}_{l=1,..,N}$ thanks to the likelihood
function $\rho$:
$\displaystyle\tilde{\omega}^{l}_{i}={\omega}^{l}_{i}\rho(Y_{i+1},x^{l}_{i+1|i})$

8: if Resampling then

9: Draw the a posteriori particles ${\left(x^{l}_{i+1}\right)}_{l=1,..,N}$
from the set ${\left(x^{l}_{i+1|i}\right)}_{l=1,..,N}$ and
${\left(\tilde{\omega}^{l}_{i}\right)}_{l=1,..,N}$ using a resampling
technique and set $\omega^{l}_{i+1}=\frac{1}{N}$

10: else

11: Set $x^{l}_{i+1}=x^{l}_{i+1|i}$ and
${\omega}^{l}_{i+1}=\tilde{\omega}^{l}_{i}/\sum_{l=1}^{N}\tilde{\omega}^{l}_{i}$

12: end if

13:end for

###### Proposition 1

Suppose that there exist $b>0$ and $\lambda_{min}\in[0,1[$, a measurable
function $V:\mathbb{R}^{n_{x}}\longrightarrow[0,+\infty[$, a compact set
$C\subset\mathbb{R}^{n_{x}}$ and an admissible state feedback control policy
$\alpha$ such that $E[V(f(x,\alpha(x),{\xi}_{0})]\leq\lambda_{min}V(x)$,
$\forall x\notin C$ and $\textrm{sup}_{x\in C}E_{x}[V(X_{1})]=b$

Then, $\forall\lambda\in[\lambda_{min},1[$,

$\displaystyle E_{x}[V(X_{k})]\leq\lambda^{k}V(x)+b(1-\lambda)^{-1},\forall
k\in\mathbb{N},\;\forall x\in\mathbb{R}^{n},$

where $X_{k}$ is computed with equation (10).

Proposition 1 is a slightly different reformulation of proposition 1 in [17],
but its proof follows from the one in [17]. In particular, we consider that
the parameter $\lambda$ can be chosen arbitrarily in $[\lambda_{min},1[$. In
practice, $\lambda$ is a parameter to tune that determines the convergence
speed of the system. Moreover, as explained in [17], if proposition 1 is
verified for a norm-like function V then for $r>0$, $P_{x}(\|X_{k}\|>r)$
decreases as the inverse of V so the distribution of state concentrates itself
around $0$. That is why, in the rest of the paper, we suppose that our guiding
goal is to drive the system (1) to 0. To do that, we also suppose that the
assumptions in proposition 1 are fulfilled and notably that the system (1) can
be stabilized with perfect information with some admissible state feedback
control policy $\alpha$.

###  V-B Receding horizon policy

In our dual MPC scheme, the information is quantified by the Fisher
information Matrix (FIM), denoted by $J_{i}$. Recursive computation of the FIM
can be found in [18]. Thus, explicit dual effect is created by minimizing some
functions of the FIM, $g_{i}^{\textrm{fish}}$ and $g_{T}^{\textrm{fish}}$ such
that, in $(P^{k}_{\textrm{ex}})$ :

$\displaystyle g_{i}^{\textrm{ex}}(X_{i},U_{i},\xi_{i})$
$\displaystyle=g_{i}^{\textrm{eco}}(X_{i},U_{i},\xi_{i})+g_{i}^{\textrm{fish}}(J_{i}),$
(11)

In $(P^{k}_{\textrm{ex}})$, the drift condition is taken from proposition 1
and, applied to $\hat{x}_{k}$ only when $\hat{x}_{k}\notin C$ such that, for
some $\lambda\in[\lambda_{min},1[$:

$\displaystyle E_{\hat{x}_{k}}[V(f(\hat{x}_{k},u_{0},{\xi}_{0})]\leq\lambda
V(\hat{x}_{k}),\;\textrm{when}\;\hat{x}_{k}\notin C.$

It is important to notice that the admissibility of the drift constraint is
guaranteed by the existence of the stabilizing admissible state feedback
policy $\alpha$.

A good approximation technique to approach $({P}_{\textrm{ex}}^{k})$ in
practice is the scenario approach [19]. It appears that, in output feedback
MPC, there is a synergy between the scenario approach and particle filtering.
In fact, the initial condition for each independent scenario is chosen as a
particle from the current set of particles. It improves global performance
compared to a similar technique involving a Kalman filer in which, the
initials conditions are always drawn according to a unimodal law. This method
was already used in [20] and [13]. The approximation, denoted by
$\left({P}_{\textrm{ex}}^{k,N}\right)$, can be defined,
$\forall\lambda\in[\lambda_{min},1[$, by:

$\begin{array}[]{rrclcc}\displaystyle\underset{u_{0}\cdots
u_{T-1}}{\text{min}}&\lx@intercol\sum_{l=1}^{N_{s}}\omega_{k}^{l}\left(\sum_{i=0}^{T-1}{g}_{i}^{\textrm{ex}}\left(X^{l}_{i},u_{i},\xi_{i}^{l}\right)+{g}_{T}^{\textrm{ex}}\left(X^{l}_{T}\right)\right)\hfil\lx@intercol\\\
\\\ \textrm{s.t.}&X^{l}_{i+1}&=&f(X^{l}_{i},u_{i},{\xi}^{l}_{i}),\\\
&u_{i}&\in&\mathbb{U}_{i},\\\ &X^{l}_{0}&=&{x}^{l}_{k},\;\forall
l=1,..,N_{s},\;\forall i=0,..,T-1,\\\ \end{array}$
$\frac{1}{N_{{dr}}}\textstyle\sum_{\ell=1}^{N_{{dr}}}V(f(\hat{x}_{k}^{N},u_{0},\tilde{\xi}_{0}^{\ell})\leq\lambda
V(\hat{x}_{k}^{N}),\;\textrm{when}\;\hat{x}_{k}^{N}\notin C,\\\ $

where:

  * •

$N_{s}<N$ is the number of scenarios considered. It is supposed to be less
than the number of particles for computational reasons, so $N_{s}$ particles
must be extracted from the original set.

  * •

$N_{dr}\in\mathbb{N}^{*}$ is the size of the sample used to approximate the
expectation in the drift constraint.

  * •

${({\xi}^{l}_{i})}_{i=0,..,T-1,l=1,..,N_{s}}$ and
${(\tilde{\xi}_{0}^{\ell})}_{\ell=1,..,N_{dr}}$ are i.i.d. random variables
sampled from $p_{\xi}$.

Additionally, when $\hat{x}_{k}^{N}\in C$ the drift condition is not
necessarily feasible so we decided to apply the stabilizing policy $\alpha$ to
$\hat{x}_{k}^{N}$, such that $U_{k}=\alpha(\hat{x}_{k}^{N})$, as it is done
for classical certainty equivalence controllers. This means that if
$\hat{x}_{k}^{N}$ enters $C$ there is no probing anymore. It is not a problem
as one can see $C$ as the target set.

###  V-C Global algorithm

The complete Output feedback algorithm is summarized in Algorithm 2. This
method has two main advantages compared to the one we presented in [13].

First, in [13], both the stability properties and the information were
incorporated as two terms of the cost. Therefore, the weights between the
terms were difficult to tune and, especially in receding horizon control where
they deeply influence the convergence of the system. In particular, it was
compulsory to decrease the weights on $g_{i}^{\textrm{fish}}$ and
$g_{T}^{\textrm{fish}}$ with time otherwise the system converged to a point
that was far from the target. Moreover, tuning the decreasing of this weight
was also complicated a priori. In our new method, the stability properties are
much less influenced by the cost because of the drift condition which is a
constraint. The most important parameter to tune is $\lambda$, and, in
principle, it influences only the convergence speed of the system and not its
qualitative properties of stability.

Secondly, from a numerical point of view, as stated in [9], the stability
properties are contained in one constraint and are easier to achieve in
practice. Indeed, the convergence of the system depends little on the quality
of the solution of the optimization problem and much more on its admissibility
which is easy to obtain with classical solvers.

Algorithm 2 Fisher/Lyapunov Output Feedback Control

1:Create a sample of $N$ particles ${x}^{l}_{0}$ according to the law
$\mu_{0}$ and initialize the weights $\omega_{0}^{l}$ with $\frac{1}{N}$.

2:for $k=0,1,2,\dots$ do

3: if $\hat{x}_{k}^{N}\notin C$ then

4: Solve $\left({P}_{\textrm{ex}}^{k,N}\right)$ starting from the particles
${x}^{l}_{k}$ and the weights $\omega_{k}^{l}$.

5: Get an optimal sequence $(u_{0}^{*},\dots,u_{T-1}^{*})$.

6: Set $U_{k}=u_{0}^{*}$.

7: else

8: Set $U_{k}=\alpha(\hat{x}_{k}^{N})$

9: end if

10: Compute the a posteriori particles $x^{l}_{k+1}$ and weights
$\omega_{k+1}^{l}$ given $U_{k}$ according to Algorithm 1

11:end for

##  VI APPLICATION AND NUMERICAL RESULTS

###  VI-A Description of the application

Algorithm 2 has been applied to the guidance and localization of a drone by
terrain-based navigation (TBN). The guiding goal is to drive a drone in a 3D
space from an uncertain initial condition $X_{0}$ to a compact set centered
around $0$. If the original target is not 0 then a translation can be made to
center the problem around 0. It is assumed that the drone can be described by
3 positions and 3 speeds,
$X_{k}={(x_{k},y_{k},z_{k},v^{x}_{k},v^{y}_{k},v^{z}_{k})}^{T}$ and the
control of 3 accelerations $U_{k}={(u^{x}_{k},u^{y}_{k},u^{z}_{k})}^{T}$.
$(x_{k},y_{k})$ represents a horizontal position and $z_{k}$ an altitude. Its
dynamics is linear with bounded controls such that, for $k\in\mathbb{N}$ :

$\displaystyle X_{k+1}=AX_{k}+BU_{k}+\xi_{k},$ $\displaystyle\|U_{k}\|\leq
U_{max},$ (12)

where $U_{max}>0$ and $A\in\mathbb{R}^{n_{x}\times n_{x}}$ and
$B\in\mathbb{R}^{n_{x}\times n_{u}}$ correspond to a discrete-time second
order system with damping on the speed. Because of the constraint on the
control, classical linear controllers are not feasible and the closed loop
system must be nonlinear. The guiding problem is addressed with a drift
constraint taken from proposition 8 in [17] such that $V(X)=e^{\|X\|}$ and C
is a ball centered around 0 whose radius depends only on the disturbances of
the dynamics. To guarantee the existence of an admissible stabilizing state
feedback policy, $\alpha$, the maximum control $U_{max}$ must be sufficiently
high. See [17] for the precise definition of $\alpha$, $U_{max}$ and $C$.

The speed is supposed to be measured. Concerning the position, the paradigm of
TBN is that only a measure of the difference between the altitude of the
drone, $z_{k}$, and the altitude of the corresponding vertical point on the
ground is available. The ground is represented by a map, $h_{m}$, that maps a
horizontal position $(x,y)$ to the corresponding height of the terrain. In
practice, $h_{m}$ is determined by a smooth interpolation of data points so it
is highly nonlinear. Therefore, the observation equation is:

$\displaystyle
h(X_{k},\eta_{k})={(z_{k}-h_{m}(x_{k},y_{k}),v^{x}_{k},v^{y}_{k},v^{z}_{k})}^{T}+\eta_{k}.$

Because of the potential multimodality of $h_{m}$, state estimation must be
dealt with by a particle filter. Moreover, it appears very naturally that dual
control is required in this application. Indeed, the quality of the
observations depends on the control and more precisely on the area that is
flied over by the drone. Let us assume that the drone flies over a flat area
with constant altitude then one measurement of height matches a whole
horizontal area and the estimation error on $(x_{k},y_{k})$ is of the order of
magnitude of the size of the area. On the contrary, if the drone flies over a
rough terrain, then one measurement of height corresponds to a smaller area on
the ground and the estimation error is reduced. Finally, an intuitive good
dual control strategy is to go toward the target avoiding to fly over
uninformative areas.

###  VI-B Numerical results

Figure 1 represents the horizontal projection of a trajectory obtained by
algorithm 2 with the terrain map in its background. The scenario has been
chosen such that, if the system goes in straight line to the target then it
flies over a flat area. We can see that in figure 1 that, as expected, the
system makes a detour to avoid the flat area, so that the particles tighten
around the true position, and finally reaches the target. In this run, the
algorithm has been stopped when the target enters the compact C.

Figure 1: Plot of one trajectory obtained by Fisher/Lyapunov particle control
and of the particles from the particle filter

Figure 2 shows RMSE in horizontal position after 30 Monte Carlo simulations of
three different policies: a constrained LQ MPC with no information probing,
the policy described in this article and the policy described in [13] where
stability was forced by the minimization of the distance to the target.

Figure 2: Plot of the RMSE in horizontal position for 3 policies: without the
FIM (blue), with the FIM and the drift (red), with the FIM and the distance to
the target (green)

The simulations were run in MATLAB and the optimization problems were solved
using the modelling language AMPL and the solver Ipopt.

##  VII CONCLUSION

In this paper, we have presented a new explicit dual output feedback
stochastic MPC for nonlinear systems. Its principle is to choose at each time
step the controls that maximizes the information over the control that forces
the mean decreasing of some Lyapunov-like function for discrete-time nonlinear
stochastic systems. In particular, it does not involve penalization of the
guiding goal in the cost which is a classical feature of MPC. Output feedback
is obtained by coupling the resolution of an optimization problem with a
particle filter. The method is applied to terrain aided navigation and appears
to be easier to tune that one of our previous method.

## References

  * [1] Ali Mesbah.  Stochastic model predictive control with active uncertainty learning: A Survey on dual control.  Annual Reviews in Control, November 2017.
  * [2] David S. Bayard and Alan Schumitzky.  Implicit dual control based on particle filtering and forward dynamic programming.  International Journal of Adaptive Control and Signal Processing, 2008.
  * [3] Giancarlo Marafioti, Robert R. Bitmead, and Morten Hovd.  Persistently exciting model predictive control: Persistently exciting model predictive control.  International Journal of Adaptive Control and Signal Processing, 28(6):536–552, June 2014.
  * [4] D. Telen, B. Houska, M. Vallerio, F. Logist, and J. Van Impe.  A study of integrated experiment design for NMPC applied to the Droop model.  Chemical Engineering Science, 160:370–383, March 2017.
  * [5] Kristian G. Hanssen and Bjarne Foss.  Scenario Based Implicit Dual Model.  In 5th IFAC Conference on Nonlinear Model Predictive Control NMPC 2015, volume 48, pages 416–241, Seville, Spain, September 2015\. Elsevier.
  * [6] Sankaranarayanan Subramanian, Sergio Lucia, and Sebastian Engell.  Economic Multi-stage Output Feedback NMPC using the Unscented Kalman Filter.  IFAC-PapersOnLine, 48(8):38–43, 2015.
  * [7] Valerii V. Fedorov and Peter Hackl.  Model-oriented design of experiments, volume 125.  Springer Science & Business Media, 2012.
  * [8] Peter Hokayem, Eugenio Cinquemani, Debasish Chatterjee, Federico Ramponi, and John Lygeros.  Stochastic receding horizon control with output feedback and bounded controls.  Automatica, 48(1):77–88, 2012.
  * [9] Tyler Homer and Prashant Mhaskar.  Output-Feedback Lyapunov-based Predictive Control of Stochastic Nonlinear Systems.  IEEE Transactions on Automatic Control, pages 1–1, 2017.
  * [10] Ali Mesbah.  Stochastic Model Predictive Control: An Overview and Perspectives for Future Research.  IEEE Control Systems, 36(6):30–44, December 2016.
  * [11] Matthew Ellis, Helen Durand, and Panagiotis D. Christofides.  A tutorial review of economic model predictive control methods.  Journal of Process Control, 24(8):1156–1178, August 2014.
  * [12] Dimitri P. Bertsekas.  Dynamic programming and optimal control. Vol. 1: […].  Number 3 in Athena scientific optimization and computation series. Athena Scientific, Belmont, Mass, 3. ed edition, 2005.  OCLC: 835987011.
  * [13] Emilien Flayac, Karim Dahia, Bruno Hérissé, and Frédéric Jean.  Nonlinear Fisher Particle Output Feedback Control and its application to Terrain Aided Navigation.  In 56th IEEE Conference on Decision and Control, Melbourne, December 2017.
  * [14] Yaakov Bar-Shalom and Edison Tse.  Dual Effect, Certainty Equivalence, and Separation in Stochastic Control.  IEEE Transactions on Automatic Control, 19(5):494–500, 1974.
  * [15] Vincent Andrieu and Laurent Praly.  A unifying point of view on output feedback designs for global asymptotic stabilization.  Automatica, 45(8):1789–1798, August 2009.
  * [16] Sean P. Meyn and Richard L. Tweedie.  Markov Chains and Stochastic stability.  Cambridge University Press, Cambridge, 2nd edition, 2009.
  * [17] Debasish Chatterjee and John Lygeros.  On stability and performance of stochastic predictive control techniques.  IEEE Transactions on Automatic Control, 60(2):509–514, 2015.
  * [18] Petr Tichavsky, Carlos H. Muravchik, and Arye Nehorai.  Posterior Cramér-Rao bounds for discrete-time nonlinear filtering.  IEEE Transactions on signal processing, 46(5):1386–1396, 1998.
  * [19] Giuseppe C. Calafiore and Marco C. Campi.  The Scenario Approach to Robust Control Design.  IEEE Transactions on Automatic Control, 51(5):742–753, May 2006\.
  * [20] Martin A. Sehr and Robert R. Bitmead.  Particle Model Predictive Control: Tractable Stochastic Nonlinear Output-Feedback MPC.  arXiv preprint arXiv:1612.00505, 2016.