# Fast Adaptive Regression-based Model Predictive Control

Eslam Mostafa, Hussein A. Aly, Ahmed Elliethy Authors are with the Department
of Electrical and Computer Engineering, Military Technical College, Cairo
14627, Egypt (e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>a.s.elliethy@mtc.eg.edu).This paper has supplementary downloadable material
available provided by the authors. The material consists of a supplementary
document.Color versions of one or more of the figures in this paper are
available in the electronic version of this manuscript.

###### Abstract

Model predictive control (MPC) is an optimal control method that predicts the
future states of the system being controlled and estimates the optimal control
inputs that drive the predicted states to the required reference. The
computations of the MPC are performed at pre-determined sample instances over
a finite time horizon. The number of sample instances and the horizon length
determine the performance of the MPC and its computational cost. A long
horizon with a large sample count allows the MPC to better estimate the inputs
when the states have rapid changes over time, which results in better
performance but at the expense of high computational cost. However, this long
horizon is not always necessary, especially for slowly-varying states. In this
case, a short horizon with less sample count is preferable as the same MPC
performance can be obtained but at a fraction of the computational cost. In
this paper, we propose an adaptive regression-based MPC that predicts the best
minimum horizon length and the sample count from several features extracted
from the time-varying changes of the states. The proposed technique builds a
synthetic dataset using the system model and utilizes the dataset to train a
support vector regressor that performs the prediction. The proposed technique
is experimentally compared with several state-of-the-art techniques on both
linear and non-linear models. The proposed technique shows a superior
reduction in computational time with a reduction of about 35-65% compared with
the other techniques without introducing a noticeable loss in performance.

###### Index Terms:

Regression analysis, MPC, Control, Parametrization, Wavelet, SVR, Optimization

##  I Introduction

Model predictive control (MPC) [1] is an advanced control method that has been
widely used in many applications[2, 3, 4, 5, 6]. Given a discrete system with
its states defined at specific sampling instances, the MPC utilizes the
mathematical model of the system to predict its future states at each sample
instant over a finite future horizon time. The predicted states along with a
set of given system constraints are used to formulate an optimization problem
that is solved online in every control cycle to estimate the optimal control
inputs at each sample instant over the horizon. Only the estimated input at
the first sample instant is applied to the system, and the MPC repeats the
same process for all subsequent control cycles. For a long horizon with a
large number of sample instances, the predicted behavior of the system becomes
more intimate to the required reference [7], i.e., better control performance.
However, for such a long horizon, the MPC easily may not finalize the
computations involved in solving the online optimization problem within the
control cycle, and this results in a lagging control input to the system in
this case.

To avoid the aforementioned lagging control input problem, several studies in
the literature focus on speeding up the MPC by reducing its online
computations. For example, several techniques try to minimize the number of
sample instances within the horizon. In [8, 9, 10], the input and state
trajectories, that represent their variations over time, is parametrized with
some basis functions. The parameterization of the trajectories reduces the
degrees of freedom of the online optimization problem by calculating the
control inputs only at specific sampling instances in the horizon while
evaluating the rest using the parametrized version. In [11, 12, 13], non-
uniformly spaced sample instances are used such that smaller intervals between
sample instances are used with the near future of the horizon while larger
ones are used with the distant future. These techniques use a fixed horizon
but keep the number of sample instances relatively small. However, using a
fixed horizon is not optimal in all control scenarios. As shown in [14], a
short horizon is enough for controlling a vehicle on highways where speed
fluctuations are not fast, while a longer horizon is needed when driving in a
city due to higher fluctuations in speed and environmental variables.
Following this, the technique in [15] adaptively changes the horizon according
to the curvature value for the state trajectories. However, the technique uses
a heuristic rule to determine the horizon length, which results in a sub-
optimal horizon. In [16], a variable horizon MPC is achieved by defining
several fixed-horizon optimization problems with different horizon lengths.
The overall complexity is relatively reduced by utilizing the time-varying
move blocking technique [17] which fixes adjacent in-time decision variables
of the optimization problem or its derivatives to be constant over several
control cycles, and thus it results in sub-optimal control performance. In
[18], the horizon length is added as an extra degree of freedom to the MPC
formulation, and thus it incurs additional computations in the control cycle.
The technique in [19] deals with systems with non-linear dynamics and
incrementally adjusts the horizon length to its minimum possible value that
guarantees stabilization. This incremental adjusting for the horizon is
practically not suitable for fast applications.

Another line of research is the so-called explicit MPC [20, 21, 22], which
pre-computes the optimal control inputs offline as a function of the current
state and reference states. Thus, the online optimization reduces to a simple
search within the pre-computed values. The major drawback here is the
searching time for the optimal solution which quickly increases when the
number of states, horizon length, or the number of control inputs increases.
Thus the explicit MPC can only be applied to small problems. This drawback is
partially mitigated in [23] which uses a partial enumeration technique that
offline computes and stores in a table the optimal control inputs for only
frequently occurring constraint sets. The table is searched online for the
best control and updated to incorporate new constraint sets as the control
progresses in time.

Recently, machine learning techniques are employed for horizon prediction. In
[24, 25], a novel technique is proposed that uses reinforcement learning (RL)
[26] to predict the optimal horizon length using a policy function of the
current state. The policy of the RL is modeled as a neural network (NN) [27]
which is trained using the data collected during the operation of the system
being controlled. The training of the NN is performed online within each
control cycle, which adds more computations within the cycle. Another
technique in [28] trains a NN offline that is used to predict the optimal
horizon at run-time. However, these techniques perform the prediction solely
based on the instantaneous state of the system, without any consideration for
the future values of the states. Moreover, the techniques do not employ any
feature engineering for training the NN, and thus it easily over-fits the
training data. Additionally, the technique in [24] predicts the horizon length
only while the technique in [28] uses the move-blocking strategy [17] that
fixes the ratio between the number of sample instances and the time span of
the predicted horizon. Therefore, the predictions performed by these
techniques are not optimal in general.

In this paper, our goal is to accurately estimate both the best minimum
horizon length and the number of samples without introducing a noticeable loss
in the performance of the MPC. To this goal, we first propose a mathematical
formulation that relates the horizon length and the sample count with the
performance of the MPC, then we propose an efficient solution for it.
Specifically, we propose an adaptive regression-based MPC (ARMPC) that
predicts the best minimum horizon length and the sample count according to the
current and future variations exhibited in the reference states of the model.
To train the regression, we build a dataset by extracting several features
that capture the variations of the reference trajectories of the model over
future time with the associated best horizon and sample count in each
situation. In run time, we extract the same set of features which presented to
the regressor to predict the best minimum horizon length and the sample count.

Compared with previous techniques, our proposed ARMPC has several advantages.
First, it estimates both the best horizon length and the best number of
samples on the horizon, and this allows the proposed technique to provide more
reduction in computations. Second, our technique does not rely on raw values
of the reference states but employs feature engineering to extract several
distinctive features from the reference trajectories, and thus it avoids over-
fitting in the learning phase. Finally, these features are extracted not from
the instantaneous values of the states but over a future span of the horizon,
which allows a more accurate estimate for the horizon and the sample count.
These advantages are reflected in our experimental results where we compared
the proposed ARMPC with three different state-of-the-art techniques on both
linear and non-linear models. The proposed ARMPC shows a superior reduction in
computational time with a reduction of about 35-65% compared with the other
techniques without introducing a noticeable loss in performance. The source
code of the proposed technique will be available111https://github.com/ahmed-
elliethy/fast-regression-mpc online upon acceptance.

The remainder of this paper is organized as follows. Section II presents a
background for the MPC and briefly outlines its parametrization. Section III
presents motivation examples and formulates our problem statement. Section IV
introduces the proposed adaptive regression based MPC. Section V describes the
experimental setup and discusses the experimental results that evaluate our
proposed technique. Section VI summarizes our conclusion and presents our
future work.

##  II MPC Background

A continuous linear time-invariant system can be described in state space form
as [10]

$\dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t),$ (1)

where $t\in\mathbb{R}$ is a continuous time instant,
$\mathbf{x}(t)\in\mathbb{R}^{m\times 1}$ and
$\mathbf{u}(t)\in\mathbb{R}^{n\times 1}$ represent a vector of states and a
vector of system inputs at time $t$, respectively,
$\mathbf{A}\in\mathbb{R}^{m\times m}$ is the state matrix, and
$\mathbf{B}\in\mathbb{R}^{m\times n}$ is the input matrix. The system can be
discretized using any of the discretization methods (such as the Euler method)
by

$\dot{\mathbf{x}}_{k}=\frac{\mathbf{x}_{k+1}-\mathbf{x}_{k}}{t_{s}},$

where $k\in\mathbb{Z}$ is a discrete-time instant, $t_{s}$ is the sampling
time, and $\mathbf{x}_{k}=[x^{1}_{k},\dots,x^{m}_{k}]^{T}$ represents the
vector of states at $k$. With the discretization, (1) can be written as

$\mathbf{x}_{k+1}=\mathbf{A}_{d}\mathbf{x}_{k}+\mathbf{B}_{d}\mathbf{u}_{k},$
(2)

where $\mathbf{u}_{k}=[u^{1}_{k},\dots,u^{n}_{k}]^{T}$ represents the vector
of system inputs at $k$, $\mathbf{A}_{d}=\mathbf{I}+t_{s}\mathbf{A}$,
$\mathbf{B}_{d}=t_{s}\mathbf{B}$, and $\mathbf{I}$ is the identity matrix.

The goal of MPC is to estimate the optimal vector of system inputs for a fixed
horizon length of $\mathcal{T}$ time steps in control cycles. More clearly, in
the $c^{\text{th}}$ control cycle, the MPC estimates the optimal vector of
system inputs from the time instant $c$ to $c+\mathcal{T}-1$. Then only the
first control input (at the time instant $c$) is applied to the system and new
system states are predicted. In the next cycle, the MPC estimates the system
inputs using the newly predicted states, then only the first input is applied
again to the system. This loop will be repeated to find the optimum control
input in every control cycle.

To do so, in the $c^{\text{th}}$ control cycle, the MPC first expresses the
system states as a function of only the inputs and the initial state
$\mathbf{x}_{c}$ in the cycle as

$\displaystyle\left[\begin{array}[]{c}\mathbf{x}_{c+1}\\\ \mathbf{x}_{c+2}\\\
\vdots\\\ \mathbf{x}_{c+\mathcal{T}}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c}\mathbf{A}_{d}\mathbf{x}_{c}+\mathbf{B}_{d}\mathbf{u}_{c}\\\
\mathbf{A}_{d}\mathbf{x}_{c+1}+\mathbf{B}_{d}\mathbf{u}_{c+1}\\\ \vdots\\\
\mathbf{A}_{d}\mathbf{x}_{c+\mathcal{T}-1}+\mathbf{B}_{d}\mathbf{u}_{\mathcal{T}-1}\end{array}\right]$
(3) $\displaystyle=\mathbf{S}_{c}\mathbf{z}_{c}+\mathbf{n}_{c},$

where

$\displaystyle\mathbf{S}_{c}$
$\displaystyle=\left[\begin{array}[]{ccccc}\mathbf{B}_{d}&\mathbf{0}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{A}_{d}\mathbf{B}_{d}&\mathbf{B}_{d}&\mathbf{0}&\cdots&\mathbf{0}\\\
&\vdots&&\\\
\mathbf{A}_{d}^{\mathcal{T}-1}\mathbf{B}_{d}&\mathbf{A}_{d}^{\mathcal{T}-2}\mathbf{B}_{d}&\mathbf{A}_{d}^{\mathcal{T}-3}\mathbf{B}_{d}&\ldots&\mathbf{B}_{d}\end{array}\right],\mathbf{n}_{c}=\left[\begin{array}[]{l}\mathbf{A}_{d}\mathbf{x}_{c}\\\
\mathbf{A}_{d}^{2}\mathbf{x}_{c}\\\ \vdots\\\
\mathbf{A}_{d}^{\mathcal{T}}\mathbf{x}_{c}\end{array}\right],$

and $\mathbf{z}_{c}\in\mathbb{R}^{n\mathcal{T}\times
1}=\left[\mathbf{u}_{c}^{T},\mathbf{u}_{c+1}^{T},\ldots,\mathbf{u}_{c+\mathcal{T}-1}^{T}\right]^{T}$.
Then, the MPC estimates the optimal vector of system inputs
$\mathbf{z}_{c}^{*}$ that minimizes a cost function $J$ in every control cycle
over $\mathcal{T}$ time steps as [10]

$\mathbf{z}_{c}^{*}=\arg\min_{\mathbf{z}_{c}}J(\mathbf{z}_{c},\mathcal{T}),$
(4)

where

$\displaystyle
J(\mathbf{z}_{c},\mathcal{T})=\sum_{k=c}^{c+\mathcal{T}-1}\left(\left[\mathbf{r}_{k}-\mathbf{x}_{k}\right]^{T}\mathbf{Q}\left[\mathbf{r}_{k}-\mathbf{x}_{k}\right]\right.\left.+\left[\mathbf{v}_{k}-\mathbf{u}_{k}\right]^{T}\mathbf{R}\left[\mathbf{v}_{k}-\mathbf{u}_{k}\right]\right),$
(5) s.t.
$\displaystyle\mathbf{x}_{k+1}=\mathbf{A}_{d}\mathbf{x}_{k}+\mathbf{B}_{d}\mathbf{u}_{k}\quad\forall
k=c,\ldots,c+\mathcal{T}-1,$
$\displaystyle\mathbf{x_{\text{min}}}\leq\mathbf{x}_{k}\leq\mathbf{x_{\text{max}}}\quad\forall
k=c,\ldots,c+\mathcal{T}-1,$
$\displaystyle\mathbf{u_{\text{min}}}\leq\mathbf{u}_{k}\leq\mathbf{u_{\text{max}}}\quad\forall
k=c,\ldots,c+\mathcal{T}-1,$
$\displaystyle\Delta\mathbf{u_{\text{min}}}\leq\Delta\mathbf{u}_{k}\leq\Delta\mathbf{u_{\text{max}}}\quad\forall
k=c,\ldots,c+\mathcal{T}-1.$

In the above cost function, $\mathbf{r}_{k}=[r^{1}_{k},\dots,r^{m}_{k}]^{T}$
and $\mathbf{v}_{k}=[v^{1}_{k},\dots,v^{n}_{k}]^{T}$ are the reference states
and reference inputs at $k$, respectively. $\mathbf{x_{\text{min}}}$, and
$\mathbf{x_{\text{max}}}$ are the minimum and maximum bounds of the states,
respectively. Similarly, $\mathbf{u}_{\text{min}}$ and
$\mathbf{u_{\text{max}}}$ are the minimum and maximum bounds of the inputs,
respectively. $\Delta\mathbf{u_{\text{min}}}$ and
$\Delta\mathbf{u_{\text{max}}}$ are the minimum and the maximum bounds of the
inputs rate of change. The matrices $Q$ and $R$ are the weight matrices for
states and inputs, respectively.

Within the sampling time $t_{s}$, the MPC controller should solve the
optimization problem in (4) to find the optimal control inputs and apply the
first control input to the system, at every control cycle. However, the
solution time of (4) may exceed $t_{s}$, especially for long horizon
$\mathcal{T}$, and thus the controller gives a lagging input to the system in
this case. In the following, we discuss the idea of parametrization that
reduces the computational time of solving (4).

### Parametrization

The optimization function in (4) is solved to find $\mathcal{T}$ vectors of
system inputs. By parametrizing the input trajectories with some basis
functions, the optimization (4) is solved but for a fewer number of control
inputs. Specifically, the system inputs are estimated only at specific sample
instances over the horizon in every control cycle, while the rest is evaluated
from the parameterized version. Without loss of generality, we build our
discussion here upon the parametrization technique in [10], which simply
linearizes the trajectories with line segments, as shown in Fig. 1.

Let $\mathcal{P}$ represents the number of sample instances over the horizon
$\mathcal{T}$ at which the control inputs are estimated from the optimization
(4) (such as the yellow points in Fig. 1). The distance between these sample
instances in the unit of time steps can be expressed as

$\Delta T=\frac{\mathcal{T}-1}{\mathcal{P}-1}.$

The input $\mathbf{u}_{k}$ can be derived by linear interpolation of the
inputs at the next sample $k_{n}$ and the previous one $k_{p}$ as

$\mathbf{u}_{k}=(1-w)\mathbf{u}_{k_{p}}+(w)\mathbf{u}_{k_{n}},$ (6)

where

$\displaystyle k_{p}$ $\displaystyle=\left\lfloor\frac{k}{\Delta
T}\right\rfloor,\;\;\;k_{n}=k_{p}+1,\;\;\textnormal{and }w=\frac{k}{\Delta
T}-k_{p}.$

Figure 1:  An example of $\mathcal{T}=5$ control inputs parameterized by
$\mathcal{P}=3$. The green circles represent the instances at which the
control inputs are computed without parametrization. After parametrization,
the computations are conducted only at the yellow circles, while the other red
points are linearly interpolated.

By modifying the dynamic equation (3) by $\mathbf{u}_{k}$ in (6), we get

$\mathbf{x}_{k+1}=\mathbf{A}_{d}\mathbf{x}_{k}+(1-w)\mathbf{B}_{d}\mathbf{u}_{k_{p}}+(w)\mathbf{B}_{d}\mathbf{u}_{k_{n}},$

thus, the cost function (5) is modified to be

$\displaystyle
J^{\prime}(\mathbf{z}_{c},\mathcal{T},\mathcal{P})=\sum_{k\in\Upsilon}\left(\left[\mathbf{r}_{k}-\mathbf{x}_{k}\right]^{T}\mathbf{Q}\left[\mathbf{r}_{k}-\mathbf{x}_{k}\right]\right.\left.+\left[\mathbf{v}_{k}-\mathbf{u}_{k}\right]^{T}\mathbf{R}\left[\mathbf{v}_{k}-\mathbf{u}_{k}\right]\right),$
(7)

s.t. the same constraints as in (5), where $\Upsilon=\\{c,c+\Delta
T,\dots,c+(\mathcal{P}-1)\Delta T\\}$. Now,

$\mathbf{z}_{c}^{*}=\arg\min_{\mathbf{z}_{c}}J^{\prime}(\mathbf{z}_{c},\mathcal{T},\mathcal{P}).$
(8)

This makes the optimization to be solved on $\mathcal{P}<\mathcal{T}$ steps.
Moreover, the dimensions of the matrices of the optimization problem are
reduced. Specifically, the dimension of $\mathbf{z}_{c}$ will be
($n\mathcal{P}\times 1$) instead of ($n\mathcal{T}\times 1$) and the matrix
$\mathbf{S}_{c}$ will have a dimension of ($m\mathcal{T}\times n\mathcal{P}$)
instead of ($m\mathcal{T}\times n\mathcal{T}$). Therefore, the parametrization
reduces the computational time required in each control cycle. However, this
reduction in computational time may be at the expense of the MPC performance,
especially for small $\mathcal{P}$. This is because the control inputs that
are evaluated using (6) may be different than the inputs estimated from the
original optimization in (4). This difference vanishes when $\mathcal{P}$ is
close to $\mathcal{T}$.

##  III Motivation and problem statement

As shown in (8), there are two factors that strongly affect both the
computational burden and the performance of the MPC, which are the horizon
length $\mathcal{T}$ and the sample count $\mathcal{P}$. To enhance the
performance of the MPC, both $\mathcal{T}$ and $\mathcal{P}$ are fixed to
specific large values for all control cycles. When the trajectories of states
have rapid changes over future time, the MPC takes the advantage of large
$\mathcal{T}$ and $\mathcal{P}$ for better preparing the appropriate inputs
$\mathbf{z}^{*}_{c}$. This allows the MPC to avoid the overshoot that may
occur at these rapid changes, which results in better performance but at the
expense of more computational cost. However, fixing $\mathcal{T}$ and
$\mathcal{P}$ are not optimal for all trajectories in every control cycle. For
example, if the trajectories do not show too many variations, then we can use
smaller values for $\mathcal{T}$ and $\mathcal{P}$ to solve (8) without
noticeable loss in control performance but at much less computational cost.
This means that it is better to set a variable value for $\mathcal{T}$ and
$\mathcal{P}$ for every control cycle according to the variations exhibited in
the trajectories. In the rest of the paper, we represent the variable horizon
length and sample count as $T_{c}$ and $P_{c}$, respectively, for the
$c^{\text{th}}$ control cycle.

In the following, we first present an experimental validation that illustrates
the effect of different fixed settings for $\mathcal{T}$ and $\mathcal{P}$ on
both the MPC’s computational time and performance followed by our problem
statement that mathematically formulates the problem of adaptive selection of
$T_{c}$ and $P_{c}$ in every control cycle.

###  III-A Motivating examples

We validate our claim on two different models. The first model is a linear
vehicle model that is built based on a simple bicycle model approximation of a
vehicle and we control the vehicle’s lateral position and orientation. The
second one is a more complex non-linear robot model that is linearized around
the operating point using Linear parameter varying (LPV)[29] which
encapsulates the nonlinearity of the robot model into a linear form. We
control the position and orientation of the robot. Both vehicle and robot
models are presented with more details in Sec. S.I and Sec. S.II, respectively
in the supplementary material. Figure 2 shows the reference state trajectories
for both the vehicle and the robot. As shown in the figure, two different sets
of reference state trajectories are examined, the first set has slow changes
(shown in black color) while the second one has rapid changes (shown in green
color).

Figure 2: The vehicle’s state reference trajectories are presented in (a) and
(b). The robot’s state reference trajectories are presented in (c), (d), and
(e). The black line represents a reference trajectory that is characterized by
slow changes, while the green line represents a reference trajectory that is
characterized by rapid changes.

We run our experiments using four different settings for $\mathcal{T}$ and
$\mathcal{P}$ of the MPC. To differentiate between these settings, we denote
MPC($a$,$b$) for a specific setting of $\mathcal{T}$ to $a$ and $\mathcal{P}$
to $b$. The four settings used in our experiments are MPC(40,40), MPC(5,5),
MPC(40,3), and MPC(40,25). Thus, the first two settings are concerned with
both long and short values for $\mathcal{T}$, respectively, without
parametrization, i.e., $\mathcal{P}=\mathcal{T}$. The other two settings are
concerned with a long $\mathcal{T}$, but with different settings for
$\mathcal{P}$. We evaluate the performance of the MPC with every setting with
both the average computational time and the average cost over all control
cycles. Assume that we have $H$ control cycles222We used the same $H$ for all
experiments., the average cost $E$ is computed as

$E=\frac{1}{H}\sum_{c=1}^{H}\frac{1}{\mathcal{T}}J^{\prime}(\mathbf{z}_{c}^{*},\mathcal{T},\mathcal{P}).$
(9)

In Fig. 3, we plot both the average cost $E$ and the average computational
time for the four different settings. As shown in the figure, if the reference
state trajectories exhibit rapid changes, the average cost becomes very large
when using a short horizon or using a long horizon with a small number of
sample counts, but the computational time, in this case, is relatively small.
For the same trajectories, when using a long horizon or using a long horizon
with a sufficient number of sample counts, the average cost becomes smaller
but at the expense of more computational time. In the case of reference state
trajectories with slow changes, the average cost becomes small, i.e., good
performance, for all different settings used in the MPC even with short
horizons or a small number of samples. However, this good performance is
obtained in a much smaller computational time when using short horizons or a
small sample count.

From the results, we can conclude that despite increasing the horizon length
and sample count improving the control performance in general, it is not
required in all situations. In real scenarios, the state trajectories may
exhibit both rapid and slow variations. Therefore, the horizon length and the
sample count should not be fixed and should be adaptively selected according
to the variations encountered in the state trajectories in every control
cycle. With this adaptive selection, the computational cost of the MPC is
reduced, while it maintains its performance unaffected.

Figure 3: Plot of the average computational time and the average cost against
different settings of MPC for (a) the vehicle and (b) the robot. If the
reference state trajectories exhibit rapid changes, the average cost is large
when using a short horizon MPC(5,5) or using a long horizon with a small
sample count of MPC(40,3). For slow variations, the average cost becomes small
for all settings. However, the computational time is much smaller for short
horizons MPC(5,5) or for a small sample count MPC(40,3).

###  III-B Problem statement

A pictorial representation of our problem is shown in Fig. 4. According to the
situation of the trajectories in the $c^{\text{th}}$ control cycle, the values
of $T_{c}$ and $P_{c}$ should be selected as minimum as possible such that the
performance of the MPC is not deteriorated compared with its performance when
using large values for horizon length and sample count. Let $\mathcal{T}^{l}$
represents this large value, then our problem can be formulated as

$\displaystyle\\{T_{c}^{*},P_{c}^{*}\\}$ $\displaystyle=\min_{T,P}\\{T,P\\}$
(10) s.t.
$\displaystyle\ell\left(\frac{1}{T}J^{\prime}(\mathbf{z}_{c}^{*},T,P),\frac{1}{\mathcal{T}^{l}}J^{\prime}({\boldsymbol{\gamma}}^{*},\mathcal{T}^{l},\mathcal{T}^{l})\right)<\epsilon,$
$\displaystyle T\in\\{2,\dots,\mathcal{T}^{l}\\},\;P<T,$

where $\epsilon$ is a small positive number and ${\boldsymbol{\gamma}}^{*}$ is
the vector of optimal control inputs that is estimated when using
$\mathcal{T}^{l}$ for both the horizon length and the sample count. In our
formulation, we assess the loss in performance of MPC by $\ell$ which measures
the relative difference of the average of the cost function $J^{\prime}$ (7)
when using the values of $\\{T,P\\}$ instead of $\mathcal{T}^{l}$ for both
horizon length and the sample count, i.e.,

$\ell\left(a,b\right)=\frac{\lvert a-b\rvert}{b}.$ (11)

Figure 4: Problem statement architecture: In each $c^{th}$ control cycle, we
need to determine what is the most suitable $\\{T,P\\}$ to be used in MPC
optimization.

A possible solution for (10) is that, in every control cycle, we can start
with small values of $T$ and $P$ then iteratively increment both values until
$\\{T_{c}^{*},P_{c}^{*}\\}$ are estimated. However, this solution is not
practical in run time (i.e., when the controller is in action) as it involves
solving the MPC optimization (8) several times in every control cycle to
obtain $\\{T_{c}^{*},P_{c}^{*}\\}$. In the following section, we present our
proposed efficient approach to estimate $\\{T_{c}^{*},P_{c}^{*}\\}$.

##  IV Proposed adaptive regression-based MPC

In this section, we present our adaptive regression-based MPC (ARMPC) scheme
which adaptively estimate $\\{T_{c}^{*},P_{c}^{*}\\}$ in every control cycle
from (10) using a regression model. Our proposed approach is illustrated in
Fig. 5. Specifically, we build a dataset to train a regression model, then we
used the regression model to predict $\\{T_{c}^{*},P_{c}^{*}\\}$. The
predicted $\\{T_{c}^{*},P_{c}^{*}\\}$ are used to estimate the required
control inputs $\mathbf{z}^{*}_{c}$ from (8) as usual. In the following, we
discuss the feature extraction, the dataset creation, and the regression model
in more detail.

Figure 5: The overall architecture of the proposed regression-based MPC. At
each control cycle, reference state trajectories along with the current states
of the system being controlled are passed through a feature extractor, then
the extracted features are fed to the regression model to predict $T_{c}^{*}$
and $P_{c}^{*}$ that will be used by the MPC in this control cycle.

###  IV-A Feature Extraction

Let
${\boldsymbol{\alpha}}_{c}^{i}=[r^{i}_{c},\dots,r^{i}_{c+\mathcal{T}^{l}-1}]$
denotes the $i^{\text{th}}$ reference state trajectory in the $c^{\text{th}}$
control cycle. We first extract a feature vector
${{\boldsymbol{\beta}}_{c}^{i}}$ from each ${{\boldsymbol{\alpha}}_{c}^{i}}$.
We designed ${{\boldsymbol{\beta}}_{c}^{i}}$ to capture the variation in the
reference state trajectory as it is the key element that determines the best
minimum values of horizon length and sample count, as we indicated earlier.
Specifically, the vector ${{\boldsymbol{\beta}}_{c}^{i}}$ composed of

  * •

Curvature $\mathcal{C}({{\boldsymbol{\alpha}}_{c}^{i}})$ is a value that
quantifies the amount by which ${{\boldsymbol{\alpha}}_{c}^{i}}$ deviates from
being a straight line. Mathematically [30],

$\mathcal{C}({{\boldsymbol{\alpha}}_{c}^{i}})=[\mathcal{K}(r^{1}_{c}),\dots,\mathcal{K}(r^{m}_{c})],$

where

$\mathcal{K}(r^{i}_{c})=\frac{1}{\mathcal{T}^{l}}\sum_{k=c}^{c+\mathcal{T}^{l}-1}\frac{\lvert\nabla^{2}r_{k}^{i}\rvert}{\left(1+(\nabla
r_{k}^{i})^{2}\right)^{3/2}},$ (12)

and $\nabla r_{k}^{i}=r_{k+1}^{i}-r_{k}^{i}$. It is clear from (12), that the
curvature value is large for curved trajectories and this value decreases
whenever the trajectory shape gets closer to a straight line.

Figure 6: Wavelet decomposition.

  * •

Wavelet coefficients $\mathcal{W}({{\boldsymbol{\alpha}}_{c}^{i}})$ which are
extracted from ${{\boldsymbol{\alpha}}_{c}^{i}}$ using the wavelet
decomposition described in the pyramid architecture in Fig. 6. Specifically,
two sets of coefficients are computed from ${{\boldsymbol{\alpha}}_{c}^{i}}$:
approximation coefficients $\mathbf{a}^{c}_{i}$ and detail coefficients
$\mathbf{d}^{c}_{i}$ [31]. The approximation and the detail coefficients are
computed by convolving ${{\boldsymbol{\alpha}}_{c}^{i}}$ with a low pass
filter and a high pass filter, respectively, followed by dyadic decimation.
The same procedure is repeated $L$ times with the approximation coefficients
and all resultant coefficients are used to form
$\mathcal{W}({{\boldsymbol{\alpha}}_{c}^{i}})$. The wavelet decomposition is
used in our features because it localizes the trajectory in both time and
frequency [31]. Thus, a stretched wavelet helps capture the slowly varying
changes in the state trajectory while a compressed wavelet helps capture
abrupt changes in the state trajectory.

Thus,

${{\boldsymbol{\beta}}_{c}^{i}}=[\mathcal{C}({{\boldsymbol{\alpha}}_{c}^{i}}),\mathcal{W}({{\boldsymbol{\alpha}}_{c}^{i}})].$
(13)

Besides capturing the variations in the reference state trajectory, we augment
our features by the error vector $\mathbf{e}_{c}=[e^{1}_{c},\dots,e^{m}_{c}]$,
where $e^{i}_{c}=\lvert r^{i}_{c}-x^{i}_{c}\rvert$ represents the absolute
error between the values of the reference state and the current state at the
time instant $c$. The reason for taking the error into account is that the
error has a noticeable effect on the selection of horizon length and sample
count. When the error is large, it is preferred to use large values for both
horizon length and sample count to perfectly derive the states to their
reference.

After we compute ${{\boldsymbol{\beta}}_{c}^{i}}$ from each
${{\boldsymbol{\alpha}}_{c}^{i}}$, we concatenate the vectors
${{\boldsymbol{\beta}}_{c}^{i}}$ for all $i=1,\dots,m$ along with the error
vector $\mathbf{e}_{c}$ to form the complete feature vector
${{\boldsymbol{\beta}}_{c}}$, i.e.,

${{\boldsymbol{\beta}}_{c}}=[{\boldsymbol{\beta}}^{1}_{c},\dots,{\boldsymbol{\beta}}^{m}_{c},\mathbf{e}_{c}]^{T}.$
(14)

To generate the reference states trajectories for feature extraction, we
synthetically utilize the discrete state space equation (2) for a given model
by randomly varying the control inputs, while satisfying their constraints.

###  IV-B Dataset Creation

We build the dataset by associating each feature vector
${{\boldsymbol{\beta}}_{c}}$ of the $c^{\text{th}}$ control cycle with its
corresponding values of the best minimum horizon $T_{c}^{*}$ and sample count
$P_{c}^{*}$. Algorithm 1 illustrates how $\\{T_{c}^{*},P_{c}^{*}\\}$ are
obtained. Specifically, our dataset is built in two successive steps. First,
we estimate $T_{c}^{*}$, then we estimate $P_{c}^{*}$ for the specific
estimated $T_{c}^{*}$. To estimate $T_{c}^{*}$, we initially set both the
horizon length and the sample count to a small value $T$, then iteratively
increment $T$ and solve (8) with both horizon length and sample count set to
$T$. In each iteration, the loss in performance (11) of the MPC is computed.
If the loss falls below a threshold $\epsilon$, $T_{c}^{*}$ is estimated and
the iterations are stopped. Otherwise, $T_{c}^{*}$ is set to the maximum
allowable horizon length $\mathcal{T}^{l}$. Once we obtain $T_{c}^{*}$, we
estimate $P_{c}^{*}$ in the next step. We set the sample count to a small
value $P$, then iteratively increment $P$ and solve (8) with horizon length
equal to $T_{c}^{*}$ and sample count set to $P$. In each iteration, the loss
(11) is computed and if it falls below $\epsilon$, the value of $P_{c}^{*}$ is
returned.

Algorithm 1 Estimation of $T_{c}^{*}$ and $P_{c}^{*}$ for the $c^{\text{th}}$
control cycle for dataset creation.

Input: $\mathcal{T}^{l},\epsilon$
Output: $T_{c}^{*},P_{c}^{*}$

1:$T\leftarrow 2$

2:$a\leftarrow\frac{1}{\mathcal{T}^{l}}J^{\prime}({\boldsymbol{\gamma}}_{c}^{*},\mathcal{T}^{l},\mathcal{T}^{l})$

3:while $T<\mathcal{T}^{l}$ do

4: $b\leftarrow\frac{1}{T}J^{\prime}(\mathbf{z}_{c}^{*},T,T)$

5: if $\ell(b,a)<\epsilon$ then

6: $T_{c}^{*}\leftarrow T$

7: break

8: end if

9: $T\leftarrow T+1$

10:end while

11:$T_{c}^{*}\leftarrow T$

12:$P\leftarrow 2$

13:while $P<T_{c}^{*}$ do

14: $c\leftarrow\frac{1}{T_{c}^{*}}J^{\prime}(\mathbf{z}_{c}^{*},T_{c}^{*},P)$

15: if $\ell(c,a)<\epsilon$ then

16: $P_{c}^{*}\leftarrow P$

17: break

18: end if

19: $P\leftarrow P+1$

20:end while

21:$P_{c}^{*}\leftarrow P$

The Algorithm 1 has high computational cost as it involves solving the MPC
optimization (8) several times to obtain $\\{T_{c}^{*},P_{c}^{*}\\}$. However,
all these computations are performed offline, i.e., in the training time. But,
once the dataset is built, we predict the values of
$\\{T_{c}^{*},P_{c}^{*}\\}$ in run time using very simple calculations, thanks
to the regression model that we discuss next.

###  IV-C Regression Analysis

Our training dataset contains $N$ records, one for each control cycle. Every
record is composed of the feature vector ${{\boldsymbol{\beta}}_{c}}$ for the
$c^{\text{th}}$ control cycle with its corresponding values of
$\\{T_{c}^{*},P_{c}^{*}\\}$. Our goal here is to build two regression models,
the first is for predicting the horizon length and the second one is for
predicting the sample count, from the given feature vector. We used the non-
linear support vector regression (SVR) [32] technique for building our
regression models because of its robustness to outliers, excellent
generalization capability, and high prediction accuracy. The objective of the
SVR is to find the hyperplane
${\boldsymbol{\omega}}^{T}\phi({{\boldsymbol{\beta}}_{c}})+b$ that holds
maximum training observation within the margin $\tau$ (tolerance level), as
shown in Fig. 7, where $\phi({{\boldsymbol{\beta}}_{c}})$ is a transformation
that maps ${{\boldsymbol{\beta}}_{c}}$ to a higher-dimensional space.
Mathematically, the goal of the training of the SVR model is to find the best
parameters ${\boldsymbol{\omega}}^{*}$ for the hyperplane by solving

$\displaystyle\\{{\boldsymbol{\omega}}^{*},b^{*}\\}=$
$\displaystyle\arg\min_{{\boldsymbol{\omega}},b}\frac{1}{2}\|{\boldsymbol{\omega}}\|^{2}+C\sum_{c=1}^{N}(\zeta^{+}_{c}+\zeta^{-}_{c}),$
(15) s.t. $\displaystyle\forall
c:T_{c}^{*}-\left({\boldsymbol{\omega}}^{T}\phi({{\boldsymbol{\beta}}_{c}})+b\right)\leq\tau+\zeta^{+}_{c},$
$\displaystyle\forall
c:{\boldsymbol{\omega}}^{T}\phi({{\boldsymbol{\beta}}_{c}})+b-T_{c}^{*}\leq\tau+\zeta^{-}_{c},$
$\displaystyle\forall c:\zeta^{+}_{c},\zeta^{-}_{c}\geq 0,$

where $C>0$ is a regularization parameter that penalizes the number of
deviations larger than $\tau$. The $\zeta^{+}_{c}$ and $\zeta^{-}_{c}$ are
slack variables that allow the regression to have errors, as shown in Fig. 7.

Figure 7: The kernel function $\phi$ transforms the data into a higher
dimensional feature space to make it possible to find a linear hyperplane that
holds maximum number of training observation within the margin.

We used the same formulation (15) to train the second SVR model that predicts
the sample counts but with a slight modification for the constraints. Since we
estimate $P_{c}^{*}$ in our dataset creation for a specific $T_{c}^{*}$, as
shown in Sec. IV-B, we augment the feature vector ${{\boldsymbol{\beta}}_{c}}$
with $T_{c}^{*}$ and treat $P_{c}^{*}$ as the required output.

The training of the regression in (15) is solved using the dual problem form
as shown in [33] with a Gaussian kernel
$K(\mathbf{x},\mathbf{y})=\exp\left(-\|\mathbf{x}-\mathbf{y}\|^{2}\right)=\phi(\mathbf{x})^{T}\phi(\mathbf{y})$.
We used the cross-validation method [34] with grid search to estimate the
values of the hyper-parameters $C$, $\tau$, $\zeta^{+}_{c}$, and
$\zeta^{-}_{c}$. As we mentioned, we build two regression models, the first is
for predicting the best horizon length $T_{c}^{*}$, and the second one is for
predicting the best sample count $P_{c}^{*}$. The output from the training is
the best parameters $\\{{\boldsymbol{\omega}}^{*}_{t},b^{*}_{t}\\}$ and
$\\{{\boldsymbol{\omega}}^{*}_{p},b^{*}_{p}\\}$ of the hyperplanes of the SVRs
corresponding to both $T_{c}^{*}$ and $P_{c}^{*}$, respectively. The complete
control process of the proposed ARMPC is outlined in Algorithm 2.

Algorithm 2 Adaptive Regression-based Model Predictive Control

1:Create the dataset as in Algorithm 1.

2:Train two SVRs by solving Eq. (15) to estimate
$\\{{\boldsymbol{\omega}}^{*}_{t},b^{*}_{t}\\}$ and
$\\{{\boldsymbol{\omega}}^{*}_{p},b^{*}_{p}\\}$.

3:$c$ $\leftarrow$ 0\.

4:while $c<H$ do

5: Compute the feature vector ${{\boldsymbol{\beta}}_{c}^{i}}$ for each
${\boldsymbol{\alpha}}_{c}^{i}$ using Eq.(13).

6: Compute the error vector
$\mathbf{e}_{c}\leftarrow[e^{1}_{c},\dots,e^{m}_{c}]$, where $e^{i}_{c}=\lvert
r^{i}_{c}-x^{i}_{c}\rvert$.

7: Compute the feature vector ${{\boldsymbol{\beta}}_{c}}$ using Eq.(14).

8: $T_{c}^{*}$ $\leftarrow$
${\boldsymbol{\omega}}^{*T}_{t}\phi({{\boldsymbol{\beta}}_{c}})+b^{*}_{t}$.

9: $P_{c}^{*}$ $\leftarrow$
${\boldsymbol{\omega}}^{*T}_{p}\phi([{{\boldsymbol{\beta}}_{c}},T_{c}^{*}])+b^{*}_{p}$.

10: Solve
$\mathbf{z}_{c}^{*}\leftarrow\arg\min_{\mathbf{z}_{c}}J^{\prime}(\mathbf{z}_{c},T_{c}^{*},P_{c}^{*})$.

11: Get the first $n$ elements from $\mathbf{z}_{c}^{*}$ as the optimal input
$\mathbf{u}_{c}\leftarrow\mathbf{z}_{c}^{*}[1:n]$.

12: Apply the optimal input $\mathbf{u}_{c}$ and update states
$\mathbf{x}_{c+1}\leftarrow\mathbf{A}_{d}\mathbf{x}_{c}+\mathbf{B}_{d}\mathbf{u}_{c}$.

13: $c\leftarrow c+1$.

14:end while

##  V Experimental results

We evaluate our proposed ARMPC on linear vehicle and non-linear robot models
used in Sec. III-A. We are interested in controlling both position and
orientation for both models. For the vehicle, we control the lateral position
only, while for the robot, we control both longitudinal and lateral positions,
assuming the states of both models are measurable to the MPC.

Independent of our experiments, we build the dataset to train our regression
model as we discussed in Sec. IV-B. We set the parameters of Algorithm 1 as
$\epsilon=0.05$ and $T^{l}=40$. The generated training datasets contain
$N=2000000$ for both models. Also, we used the Daubechies wavelet with order 2
[35] to perform the wavelet decomposition with $L=3$ levels. The parameters of
our regression model are estimated with 5-fold cross-validation and their
values are $C=7.4129$, $\tau=0.62$, and $\zeta^{+}_{c}=\zeta^{-}_{c}=0.1$.
Figure S.3 in the supplementary material shows the reference spatial
trajectories for both the vehicle and the robot. Also, Fig. S.4 in the
supplementary material shows their reference state trajectories, which as in
the reality, they may contain both rapid and slow variations as shown in the
figure. In all conducted experiments, the output measurements from both models
are obtained with added noise to simulate the modeling errors and/or
disturbances that may occur in a real scenario. An extended Kalman filter [36]
is used in state estimation. All experiments are carried out for 40 seconds.
We used Matlab/Simulink software running under windows 10, with a PC (i7-8550U
CPU @1.80 GHz, 16 GB RAM).

In the following three subsections, we first present the performed experiments
that compare the proposed ARMPC with other state-of-the-art techniques. Then,
we present an ablation study that compares the proposed ARMPC with defeatured
versions obtained by introducing some modifications to the proposed technique.
Finlay, we discuss the behavior of the proposed ARMPC when encountering a non-
optimal control situation that may arise due to the in-feasibility of solving
the MPC optimization before the control cycle time is over or due to system
characteristics and constraints.

###  V-A Comparison with other techniques

We compare the proposed ARMPC with three different state-of-the-art MPC
techniques, which are the parametrized MPC (PMPC) [10], the adaptive dual MPC
(ADMPC) [15], and the adaptive neural network MPC (ANMPC) [28] techniques.
Additionally, we include the standard MPC (SMPC) with a long horizon as a
baseline for performance. We used the same notations in Sec. III-A to denote
both the PMBC and SMPC, which are MPC(40,$\mathcal{P}$) and MPC(40,40),
respectively. Our comparison is performed using two experiments. In the first
experiment, we tune the parameters of every technique to provide the best
performance so we can compare the amount of speed-up of each technique with
the baseline. In the second experiment, we tune the parameters to give the
same speed-up as our proposed ARMPC so we can compare the performance of each
technique. The parameters of all techniques for both experiments are listed in
Sec. S.III in the supplementary material. To ensure a fair comparison, the
parameters of the MPC’s optimization (7) and constraints settings in all
compared techniques are kept the same as shown in Table S.II and Table S.III,
respectively, in the supplementary material. We used the quadprog [37] convex
solver for solving the MPC’s optimization (7) for all techniques. For every
technique, we measure the time span from the start of every control cycle till
the estimation of the control inputs, then we report the average computational
time across all cycles.

Figure 8: The results of the first experiment in which we tune the parameters
of every technique to provide the best performance. The proposed ARMPC
provides a much smaller average computational time with comparable performance
to the MPC with long horizon MPC(40,40). Also, the proposed ARMPC shows a
35%-52% reduction in the average computational time compared with the other
techniques without loss in performance.

The results333A video that shows simulation for the vehicle controlled by the
proposed ARMPC is shown in Sec. S.IV in the supplementary material. of the
first experiment are shown in Fig. 8. Specifically, we plot both the average
cost $E$ in (9) and the average computational time for all techniques under
comparison. As shown in the figure, the proposed ARMPC provides a much smaller
average computational time with comparable performance to the MPC(40,40).
Specifically, the proposed ARMPC reduces the average computational time by 61%
and 64.5% compared with MPC(40,40), for the vehicle and the robot,
respectively. Also, the proposed ARMPC shows a 35%-52% reduction in the
average computational time compared with the other techniques without loss in
performance. This reduction in computational time is obtained because our
proposed ARMPC does not estimate the horizon length only and it does not rely
on raw values of the reference states but employs feature engineering to
prevent overfitting in the learning phase. Additionally, the proposed ARMPC
extracts these features over a future span of the horizon which leads to a
more accurate estimation. Conversely, both ADMPC and PMPC techniques fix the
sample count so it is required to use a large enough sample count to give good
performance in all situations which reflects in the overall average
computational time. Also, the ANMPC technique fixes the ratio between the
predicted sample count and the time span of the horizon, uses only the
instantaneous values of the states, and does not employ any feature
engineering for making predictions. These drawbacks make the technique provide
a poor estimation of the sample count which is reflected in the figure where
the ANMPC technique shows the poorest performance with high computational
time.

Figure 9: Density distribution for the estimated sample count by the ARMPC and
ANMPC techniques. The proposed ARMPC provides a smaller number of occurrences
of a large sample count, while the ANMPC technique provides large number of
occurrences of a large sample count due to the wrong estimation.

In Fig. 9, we plot the histograms of the predicted sample count that is
obtained by the proposed ARMPC and the ANMPC techniques for the first
experiment. The histograms show that the proposed ARMPC provides a smaller
number of occurrences of a large sample count which reinforces the conclusion
drawn from Fig. 8 that the proposed ARMPC has a smaller average computational
time compared with the ANMPC technique. Note that, we omit the other
techniques from the histograms because these techniques do not estimate
variable sample count.

The results of the second experiment are shown in Fig. 10 where we plot both
$E$ and the average computational time for all techniques. As shown in the
figure, all techniques fail to provide comparable performance with the
proposed ARMPC for the same average computational time. For the PMPC
technique, we decrease its sample count till we reach the same average
computational time as our proposed ARMPC. However, this negatively affects the
performance, since there are some situations (states with large variations),
and the MPC needs enough samples to represent the horizon length. Also, for
the ADMPC technique, modifying the number of dense samples, the number of the
sparse sample or the heuristic threshold leads to speeding up the
computational time but at expense of the performance. The ANMPC technique, as
shown in the first experiment, provides a poor estimation of the sample count.
So to enforce reducing its computation, we limit the range of the sample
counts that are used in its training in this experiment and thus the technique
always predicts a small sample count to save computations. However, as shown
in the figure, the performance is negatively affected.

Figure 10: The results of the second experiment in which we tune the
parameters of each technique to give the same speed-up as our proposed ARMPC.
All techniques fail to provide comparable performance with the proposed ARMPC
for the same average computational time.

###  V-B Ablation study: comparison with defeatured versions

In this section, we study the effect of each feature proposed in Sec.IV-A on
the overall performance and computational time of the proposed ARMPC.
Specifically, we compare the proposed ARMPC with three de-featured versions.
Each de-featured version is obtained by dropping one of the proposed features
from the training of the SVR and consequently from the online prediction. We
denote (ARMPC-W), (ARMPC-C), and (ARMPC-E) for the de-featured versions of the
proposed ARMPC that is obtained by dropping the wavelet, the curvature, and
the error features, respectively. Additionally, we study the effect of
predicting only the horizon length (without the sample count). We denote
(ARMPC-P) for the de-featured version that trains only one SVR to predict the
horizon length only. For this version, we always set
$\mathcal{P}=\mathcal{T}$.

Figure 11: Plot of the average computational time and the average cost of the
proposed ARMPC in comparison with four defeatured versions. Drooping any of
the proposed features has a negative effect on the performance and increases
the computational time due to the wrong estimation of optimal
$\\{\mathcal{T},\mathcal{P}\\}$.

Figure 11 shows $E$ and the average computational time of the proposed ARMPC
in comparison with the four defeatured versions. As shown, predicting the
horizon length only without estimating the best minimum sample count as in
ARMPC-P increases the average computational time with no significant
enhancement in the performance. Thus, estimating both the sample count and the
horizon length has a very important impact on the computational time
reduction. Also, dropping any feature from the proposed features affects the
quality of the estimation which is negatively reflected in the performance.
This is because the dropped features make the rest indistinctive and in this
case, the dataset may have records with the same feature values but with two
different associated values of $\\{\mathcal{T},\mathcal{P}\\}$. For example,
suppose that there is a situation where the reference state trajectories have
slow variations and small curvature. If the instantaneous error is large in
this case, then Algorithm 1 will choose a large value for
$\\{\mathcal{T},\mathcal{P}\\}$ and vice-versa, even with such slow variations
and small curvature. Thus if the error feature is dropped, then several
records will appear in the dataset with similar values of the curvature and
wavelet features but with different associated
$\\{\mathcal{T},\mathcal{P}\\}$. This results in low fitting for the SVR and
will lead to non-optimal predictions. Similarly, dropping the wavelet features
decreases the performance as the remaining features can not capture whether
the state trajectories have rapid or slow variations. Therefore, a dataset
record may be constructed with a small value of ${\mathcal{T},\mathcal{P}}$ if
the instantaneous error is small regardless of the variations encountered in
the trajectories. So we can conclude that all parts of the proposed features
are very important for obtaining the best performance in the least amount of
computational time.

###  V-C Infeasibility and non-optimal control

In this subsection, we discuss the behavior of the proposed ARMPC when
encountering a non-optimal control situation that may arise due to the in-
feasibility of solving the MPC optimization before the control cycle time is
over or due to system characteristics and constraints.

When the in-feasibility happens because the optimization incurs high
computational costs due to a large horizon and can not find the optimal
solution before the control cycle time is over, the optimization returns a
sub-optimal solution. To compare our proposed ARMPC approach with the SMPC
approach, we conduct the following experiment on both the vehicle and the
robot models. We limit the maximum number of iterations $\mathcal{N}$ the MPC
solver allows to find the optimal solution and quantify the percentage of the
number of times both the ARMPC and the SMPC approaches converge to the optimal
solution from the total number of control cycles $H$. We denote this
percentage by $\sigma$. We repeat the same experiment by varying $\mathcal{N}$
and record the associated $\sigma$ for both ARMPC and SMPC approaches and plot
$\sigma$ against $\mathcal{N}$ in Fig. 12. As shown in the figure, when
$\mathcal{N}$ is large, both ARMPC and SMPC converge to the optimal solution
often. However, when $\mathcal{N}$ is smaller, the proposed ARMPC shows a
significantly higher percentage of converging to the optimal solution. This is
because the ARMPC approach estimates the optimal $\\{T_{c}^{*},P_{c}^{*}\\}$,
and therefore, it has smaller calculations involved in finding the optimal
solution. Thus, the proposed ARMPC approach shows a significant advantage over
SMPC in this case.

Figure 12: A plot of the percentage of converging to the optimal solution
$\sigma$ for different values of the MPC solver’s allowable maximum number of
iterations $\mathcal{N}$. For large values of $\mathcal{N}$, both proposed
ARMPC and SMPC approaches show high values of $\sigma$. However, For small
values of $\mathcal{N}$, ARMPC shows a significantly higher $\sigma$ than
SMPC. This is because the proposed ARMPC adaptively estimates the best minimum
horizon in every control cycle, and therefore, it has smaller calculations
involved in finding the optimal solution.

In the second case, when the in-feasibility happens due to system
characteristics and constraints (for example, when the constraints are too
tight), we compare the behavior of the ARMPC and the SMPC approaches through
the following experiment. We alter the two models’ characteristics by adding
more tightened constraints to the MPC optimization problem to force the
optimization problem to be infeasible. The added constraints are presented for
both the vehicle and the robot models in more detail in Table S.IV in the
supplementary material. In figure 13, we plot $\sigma$ for different values of
$\mathcal{N}$. The figure shows that both ARMPC and SMPC fail to get
acceptable $\sigma$ values. Specifically, neither has a value of $\sigma$
greater than 2% even with higher values of $\mathcal{N}$. This is simply
because the solution can not be found due to the MPC’s optimization problem
infeasibility.

Figure 13: A plot of the percentage of converging to the optimal solutions
$\sigma$ for different values of solver’s allowable maximum number of
iterations $\mathcal{N}$ in case of a tightly constrained optimization
problem. Both ARMPC and SMPC fail to show an acceptable percentage of optimal
solutions due to the infeasibility of the optimization problem.

##  VI Conclusion and future work

In this paper, we propose an adaptive regression-based MPC (ARMPC) technique
that predicts the best minimum horizon length and the sample count from
several features extracted from the reference state trajectories of the system
being controlled. The features are designed to capture the variation in the
trajectories by using the wavelet decomposition coefficients and the curvature
value, in addition to the instantaneous error between the reference state and
the current state. We conducted several experiments on both linear and non-
linear models to compare the proposed technique with three different state-of-
the-art techniques. The results show that the proposed technique provides a
superior reduction in computational time with a reduction of about 35-65%
compared with the other techniques without introducing a noticeable loss in
performance. Additionally, we showed experimentally that dropping any of the
proposed features makes our regression model not provide an accurate
estimation for the best minimum horizon length and the sample count which
affects both the performance and the computational time.

In the future, we plan to apply the proposed approach to non-linear MPC with
non-linear MPC solvers, such as the genetic algorithm. Another direction is to
apply machine learning techniques to the genetic algorithm to adaptively
select the best parameters according to some features extracted from the
reference and state trajectories.

## References

  * [1] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control: Theory and practice-A survey,” _Automatica_ , vol. 25, no. 3, pp. 335–348, 1989\.
  * [2] D. Hrovat, S. Di Cairano, H. E. Tseng, and I. V. Kolmanovsky, “The development of model predictive control in automotive industry: A survey,” in _2012 IEEE International Conference on Control Applications_. IEEE, 2012, pp. 295–302.
  * [3] M. Bakošová and J. Oravec, “Robust MPC of an unstable chemical reactor using the nominal system optimization,” _Acta Chimica Slovaca_ , vol. 7, no. 2, pp. 87–93, 2014.
  * [4] C. Wang and Z. Song, “Convex model predictive control for rocket vertical landing,” in _2018 37th Chinese Control Conference (CCC)_ , 2018, pp. 9837–9842.
  * [5] D. Song, Q. Chang, S. Zheng, S. Yang, J. Yang, and Y. H. Joo, “Adaptive model predictive control for yaw system of variable-speed wind turbines,” _Journal of Modern Power Systems and Clean Energy_ , vol. 9, no. 1, pp. 219–224, 2020.
  * [6] E. Mostafa, A. Elliethy, and H. A. Aly, “Response-time enhancement of a hybrid controller of unmanned aerial vehicles using parametrization,” in _2022 13th International Conference on Electrical Engineering (ICEENG)_. IEEE, 2022, pp. 150–153.
  * [7] J. M. Maciejowski, _Predictive control: with constraints_. London: Pearson education, 2002.
  * [8] B. Khan and J. A. Rossiter, “Alternative parameterisation within predictive control: a systematic selection,” _International Journal of Control_ , vol. 86, no. 8, pp. 1397–1409, 2013.
  * [9] M. Muehlebach and R. D’Andrea, “A method for reducing the complexity of model predictive control in robotics applications,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 3, pp. 2516–2523, 2019.
  * [10] P. Hyatt, C. S. Williams, and M. D. Killpack, “Parameterized and GPU-parallelized real-time model predictive control for high degree of freedom robots,” _arXiv preprint arXiv:2001.04931_ , 2020.
  * [11] J. Funke, M. Brown, S. M. Erlien, and J. C. Gerdes, “Collision avoidance and stabilization for autonomous vehicles in emergency scenarios,” _IEEE Transactions on Control Systems Technology_ , vol. 25, no. 4, pp. 1204–1216, 2016\.
  * [12] P. Ogren and J. Robinson, “Receding horizon control of uavs using gradual dense-sparse discretizations,” in _AIAA Guidance, Navigation, and Control Conference_ , 2010, p. 8084.
  * [13] B.-A. Kim, Y. S. Son, S.-H. Lee, and C. C. Chung, “Model predictive control using dual prediction horizons for lateral control,” _IFAC Proceedings Volumes_ , vol. 46, no. 10, pp. 280–285, 2013.
  * [14] A. Rezaei and J. B. Burl, “Effects of time horizon on model predictive control for hybrid electric vehicles,” _IFAC-PapersOnLine_ , vol. 48, no. 15, pp. 252–256, 2015.
  * [15] M. Kim, D. Lee, J. Ahn, M. Kim, and J. Park, “Model predictive control method for autonomous vehicles using time-varying and non-uniformly spaced horizon,” _IEEE Access_ , vol. 9, pp. 86 475–86 487, 2021.
  * [16] R. C. Shekhar and J. M. Maciejowski, “Robust variable horizon MPC with move blocking,” _Systems & Control Letters_, vol. 61, no. 4, pp. 587–594, 2012\.
  * [17] R. Cagienard, P. Grieder, E. Kerrigan, and M. Morari, “Move blocking strategies in receding horizon control,” in _IEEE Conference on Decision and Control_ , vol. 2, 2004, pp. 2023–2028 Vol.2.
  * [18] P. O. Scokaert and D. Q. Mayne, “Min-max feedback model predictive control for constrained linear systems,” _IEEE Transactions on Automatic control_ , vol. 43, no. 8, pp. 1136–1142, 1998.
  * [19] A. J. Krener, “Adaptive horizon model predictive control,” _IFAC-PapersOnLine_ , vol. 51, no. 13, pp. 31–36, 2018.
  * [20] A. Bemporad, F. Borrelli, and M. Morari, “Piecewise linear optimal controllers for hybrid systems,” in _American Control Conference_ , vol. 2. IEEE, 2000, pp. 1190–1194.
  * [21] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “The explicit linear quadratic regulator for constrained systems,” _Automatica_ , vol. 38, no. 1, pp. 3–20, 2002.
  * [22] A. Alessio and A. Bemporad, “A survey on explicit model predictive control,” in _Nonlinear model predictive control_. New York: Springer, 2009, pp. 345–369.
  * [23] G. Pannocchia, J. B. Rawlings, and S. J. Wright, “Fast, large-scale model predictive control by partial enumeration,” _Automatica_ , vol. 43, no. 5, pp. 852–860, 2007.
  * [24] E. Bøhn, S. Gros, S. Moe, and T. A. Johansen, “Reinforcement learning of the prediction horizon in model predictive control,” _IFAC-PapersOnLine_ , vol. 54, no. 6, pp. 314–320, 2021.
  * [25] ——, “Optimization of the model predictive control meta-parameters through reinforcement learning,” _arXiv preprint arXiv:2111.04146_ , 2021.
  * [26] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. Massachusetts: MIT press, 2018.
  * [27] K. Suzuki, _Artificial neural networks: Architectures and applications_. Norderstedt, Germany: BoD–Books on Demand, 2013.
  * [28] M. S. M. Gardezi and A. Hasan, “Machine learning based adaptive prediction horizon in finite control set model predictive control,” _IEEE Access_ , vol. 6, pp. 32 392–32 400.
  * [29] M. Misin, “Model predictive controller of a UAV using the LPV approach,” Master’s thesis, Universitat Politècnica de Catalunya, 2020.
  * [30] L. Najman and P. Romon, _Modern approaches to discrete curvature_. New York: Springer, 2017, vol. 2184.
  * [31] S. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 11, no. 7, pp. 674–693, 1989.
  * [32] S. R. Gunn _et al._ , “Support vector machines for classification and regression,” _ISIS technical report_ , vol. 14, no. 1, pp. 5–16, 1998.
  * [33] V. Vapnik, _The nature of statistical learning theory_. Netherlands: Springer science & business media, 1999\.
  * [34] J. Chorowski, J. Wang, and J. M. Zurada, “Review and performance comparison of SVM-and ELM-based classifiers,” _Neurocomputing_ , vol. 128, pp. 507–516, 2014.
  * [35] I. Daubechies, _Ten lectures on wavelets_. Philadelphia: SIAM, 1992.
  * [36] Y. Kim and H. Bang, “Introduction to kalman filter and its applications,” in _Introduction and Implementations of the Kalman Filter_ , F. Govaers, Ed. Rijeka: IntechOpen, 2018, ch. 2.
  * [37] A. Geletu, “Solving optimization problems using the matlab optimization toolbox-a tutorial,” _TU-Ilmenau, Fakultät für Mathematik und Naturwissenschaften_ , 2007.


# Density of States and Magnetic Correlations at a Metal-Mott Insulator
Interface

M. Jiang1,2, G.G. Batrouni,3,4 and R.T. Scalettar1 1Physics Department,
University of California, Davis, California 95616, USA  2Department of
Mathematics, University of California, Davis, California 95616, USA  3INLN,
Université de Nice-Sophia Antipolis, CNRS; 1361 route des Lucioles, 06560
Valbonne, France  4Institut Universitaire de France

###### Abstract

The possibility of novel behavior at interfaces between strongly and weakly
correlated materials has come under increased study recently. In this paper,
we use determinant Quantum Monte Carlo to determine the inter-penetration of
metallic and Mott insulator physics across an interface in the two dimensional
Hubbard Hamiltonian. We quantify the behavior of the density of states at the
Fermi level and the short and long range antiferromagnetism as functions of
the distance from the interface and with different interaction strength,
temperature and hopping across the interface. Induced metallic behavior into
the insulator is evident over several lattice spacings, whereas
antiferromagnetic correlations remain small on the metallic side. At large
interface hopping, singlets form between the two boundary layers, shielding
the two systems from each other.

###### pacs:

71.10.Fd, 71.30.+h, 02.70.Uu

##  I Introduction

Over the last decade, advances in synthesis techniques have made possible the
construction of well-defined interfaces involving strongly correlated
materials, notably transition metal oxides Zubko11 ; Mannhart10 ; Munakata11 ;
Freeland10 . Indeed, the early discovery that the interface between two
insulating oxides, LaAlO3 and SrTiO3 is a high-mobility two dimensional
conductor Ohtomo02  or even a superconductor Thiel06  has emphasized that
novel physics can arise at such boundaries, beyond a simple interpolation
between materials properties on either side. Subsequently, numerous
heterostructure interfaces have been shown to exhibit unique phenomena that
are not present in their bulk constituents. In particular, Munakata et al.
presented a subtle proximity effect that arises between a normal metal and an
antiferromagnetic Mott insulator, which can be understood in the framework of
Ruderman-Kittel-Kasuya-Yoshida (RKKY) interactions Munakata11 .

Advances in characterization techniques are also expanding the horizon of the
field. It was shown recently that, by setting up a standing wave of the
incident X-rays and adjusting the position where the intensity is peaked
inside the sample, one can measure electronic excitations as a function of
location in the sample. This standing wave angle resolved photoemission
spectroscopy (SW-ARPES) gray10  method allows a depth selective probe of
buried layers and interfaces, and hence the construction of chemical and
electronic structure profiles within the sample.

Various numerical methods have been employed to model qualitative features of
these experiments. Early attempts include the study of LaAlO3 and SrTiO3
systems in which the electric fields arising from the La3+/Sr2+ charge
difference and the on-site interactions are treated within the Hartree-Fock
approximation Okamoto . A metallic interface between band and Mott insulators
in a quasi-one-dimensional lattice was explored with the Lanczos method
kancharla06 . Important insights into electronic properties at interfaces have
also been gained by dynamical mean field theory (DMFT) OkamotoDMFT . For
example, it has been suggested that the Kondo effect governing the interface
between metal and Mott insulator is inefficient so that Mott insulators are
impenetrable to the metal helmes08 . More specifically, the quasiparticle
weight decays as $1/x^{2}$ with distance $x$ from the metal, but the prefactor
of this decay was found to be very small. Zenia et al., however, emphasized
that even a small proximity effect can induce density of states to open up a
metallic channel inside an insulator sandwiched between to metals at
sufficiently low temperatures zenia09 , leading to perfect conductance. This
Fermi liquid is “fragile”: finite temperature, disorder, or frequency rapidly
return the behavior to that of a conventional N-I-N junction. Possible device
applications were suggested as a consequence of this sensitivity. Zubko11 ;
Mannhart10 ; Munakata11 ; Freeland10  In a related work Gutzwiller10 , a
Gutzwiller approximation approach was extended to inhomogeneous systems. The
decay length of the exponential fall-off of the penetration of metallic
character into the insulator region was shown to diverge as the metal-
insulator transition is approached.

A further interesting set of studies involves the dynamical response of
strongly correlated systems with an interface, for example the possibility of
the suppression of charge transport, current rectification, and the behavior
of holon-doublon pairs, which are important for devices silva10 . So far, such
dynamic phenomena have been explored primarily for one-dimensional systems,
where time-dependent Hartree Fock yonemitsu07  and density-matrix
renormalization group (DMRG) methods are especially effective white04 .

In this paper, we address the metal-insulator interface problem by using
determinant Quantum Monte Carlo (DQMC)blankenbecler81 ; white89 , a
numerically exact approach, to treat correlated electron models. Here we focus
on an inhomogeneous Hubbard model on a single two dimensional lattice divided
by a linear interface in two regions, one weakly correlated (with $U=0$), and
one at intermediate coupling (with $U$ non-zero). The chemical potential is
chosen to maintain particle-hole symmetry, so that all sites of the lattice
are on average half-filled. This choice avoids the fermion sign problem loh90
and enables the evaluation of magnetic correlations and the density of states
at low temperature.

Such DQMC simulations are at present limited to lattices of $400-1000$ sites
(depending on the interaction strength and temperature). By focussing on a
linear interface in a 2D lattice, we are able to explore systems with a fairly
large linear extent. Most of our results will be for $20\times 20$ lattices.
This enables us, for example, to evaluate the penetration depth across the
boundary, since the interface effects have sufficient room to heal before the
lattice edge. The behavior of the Hubbard model across a planar interface in a
3D material on systems of smaller linear extent has been explored by Euverte
et al euverte12 . The key finding is that up to an interface hopping $V$ which
is on the order of the bulk hybridization $t$, the effects of the interface
can extend well past the two layers immediately at the interface. That is, the
interface affects properties $3-4$ layers deep on both the strongly and weakly
correlated sides of the boundary. When $V\geq 2t$, however, there is a return
to the values characteristic of decoupled materials with $V=0$. This “revival”
of magnetic order on the insulating side is driven by the formation of
singlets between fermions on the two layers immediately adjacent to the
interface which acts to decouple the two materials.

This paper is organized as follows: In section II we explicitly write down the
Hamiltonian and provide a very brief overview of DQMC. Section III focuses on
the magnetic properties of the 2D lattice with a linear interface, both the
near-neighbor spin correlation and the structure factor. Here the interaction
strength $U$ and temperature $T$ are varied at fixed interface hopping $V$.
The key result is a determination of the penetration depth as a function of
$T$ and $U$. Section IV extends these results to the density of states.
Section V examines the variation with $V$, and Section VI summarizes our
results.

##  II Model and methodology

We consider the two-dimensional Hubbard model

$\begin{split}\hat{H}=&-t\sum\limits_{\langle
ij\rangle,L,\sigma}(c^{\dagger}_{iL\sigma}c_{jL\sigma}^{\vphantom{dagger}}+h.c.)-\mu\sum\limits_{i,L,\sigma}n_{iL\sigma}\\\
&-\sum\limits_{i,\langle
LL^{\prime}\rangle,\sigma}t_{LL^{\prime}}(c^{\dagger}_{iL\sigma}c_{iL^{\prime}\sigma}^{\vphantom{dagger}}+h.c.)\\\
&+\sum\limits_{i,L}U_{L}(n_{iL\uparrow}-\frac{1}{2})(n_{iL\downarrow}-\frac{1}{2})\end{split}$
(1)

Here $c^{\dagger}_{iL\sigma}(c^{\phantom{\dagger}}_{iL\sigma})$ are fermionic
creation (destruction) operators of spin $\sigma$ at site $i$ in line $L$.
While $i,L$ can also be regarded as the $x$ and $y$ site labels $i_{x},i_{y}$,
the current notation emphasizes the broken rotational symmetry (the interface
is parallel to the $x$ axis). Our labeling convention is $L=\cdots-3,-2,-1$
for the $U=0$ lines and $L=1,2,3,\cdots$ for the $U\neq 0$ lines, so that the
interface connects $L=\pm 1$. $t$ and $t_{LL^{\prime}}$ are the intra and
inter-line nearest-neighbor hoppings. $t_{LL^{\prime}}=t$ except at the
interface where $t_{-1,1}=V$. The geometry is shown in Fig. 1.

Figure 1: (color online) Geometry of the two-dimensional lattice with a one-
dimensional interface (dashed red box). The metallic lines (negative L) have
$U=0$ while the insulating lines (positive L) have $U\neq 0$. The lattice has
periodic boundary conditions in the $x$ direction and open boundary conditions
in the $L$ direction. $V$ denotes the hopping/hybridization (heavy lines)
across the interface. We adopt the lattice size $N=20\times 20$ throughout the
paper, unless otherwise indicated.

The DQMC approach provides an exact solution to the finite temperature
properties of $H$ by decoupling the interaction term $U$ through the
introduction of an auxiliary Hubbard-Stratonovich (HS) field. The resulting
action is quadratic in the fermion operators, so that the trace over those
coordinates can be performed, leaving determinants (one for spin up and one
for spin down) of matrices whose dimensions are the spatial lattice size, and
whose entries depend on the HS field. The auxiliary field is sampled
stochastically, and the values of the up and down spin fermion Green’s
function for the configurations generated are used to construct the various
observables. The CPU time for the algorithm scales as the cube of the lattice
size and linearly footnote3  with inverse temperature $\beta$.

At low temperatures, the fermion determinants can become negative, precluding
their use as a probability for sampling the HS field. To avoid this “sign
problem”, we consider the half-filled case, $\mu=0$. Here particle-hole
symmetry (PHS) implies that $\langle n_{iL\sigma}\rangle=0.5$ for all
temperatures and interaction strengths, even if $U_{L}$ vary spatially with
$L$ and the hybridizations are not all equal. PHS can also be used to
demonstrate that the up and down spin determinants, although they can
individually become negative, always have the same sign. As a consequence,
their product is always positive, allowing the study of the low temperature
behavior at the metal-insulator interface.

Figure 2: (color online) Spin correlation $C_{L}(x)=\langle
S_{i,L}S_{i+x,L}\rangle$ as a function of separation $x$ for different lines
$L$. Diminishing antiferromagnetism of the insulator adjacent to the interface
is shown. The induced antiferromagnetic long-range orders in the metallic
region are evident in spite of their weaknesses. $U=0|4$ denotes the interface
setup.

##  III Magnetic properties

We gain quantitative information concerning the spin correlations in different
lines by measuring the spin correlation $C_{L}(x)=\langle
S^{-}_{i,L}S^{+}_{i+x,L}\rangle$, with
$S^{-}_{i,L}=c^{\phantom{\dagger}}_{i\uparrow}c^{\dagger}_{i\downarrow}$,
parallel to the interface. Figure 2 shows that antiferromagnetism in the
insulating layer $L=1$ immediately adjacent to the metal is diminished by
20-30% relative to larger positive $L$ which are deeper in the Mott insulator.
Data are shown for boundary hybridization $V=t$ and interaction $U=4t$ in the
insulator. Despite a somewhat smaller amplitude, the correlations in $L=1$
appear to remain long ranged. Induced antiferromagnetic long-range
orderfootnote1  in the metallic region is also evident, although its amplitude
is an order of magnitude smaller than that on the insulating side.

Figure 3: (color online) AF structure factor $S_{AF}(L)$ in each line $L$ for
different correlation strengths $U$ and temperatures $T$. The smooth
penetration of antiferromagnetic order into the metal and the diminishing AF
order in the insulator are clearly seen. Fitting $S_{AF}(L)$ to a hyperbolic
function openboundary , $a\tanh L/\lambda+b$ allows the extraction of a
penetration length $\lambda$, shown in the insets. $\lambda$ decreases with
larger correlation $U$ and lower temperature $T$. Further information

Compact information on the antiferromagnetic long-range order for each line
$L$ can be obtained by measuring the antiferromagnetic (AF) structure factor,

$S_{AF}(L)=\frac{1}{N}\sum\limits_{x}(-1)^{x}\cdot C_{L}(x)$ (2)

where $N$ is the linear lattice size. Figure 3 illustrates the $L$-dependent
AF structure factor for different correlation strengths $U$ and temperatures
$T$. The smooth penetration of antiferromagnetic order into the metal and the
diminishing AF order in the insulator are clearly seen. As might be expected,
the induced long range AF order in the metal adjacent to the interface is
stronger for correlated insulator with largerfootnote4  on-site repulsion $U$
(Fig. 3a) and/or lower temperature (Fig. 3b). The overall magnitude of
$S_{AF}$ in the metal is rather small, consistent with the small real-space
correlations in Fig. 2. The results of Figs. 2,3 are rather different from
previous work of Sherman et. al. Sherman  which found that the
antiferromagnetic order can penetrate into the metal to a depth of ten lattice
spacings, but are consistent with the shorter range effects described in
helmes08 ; zenia09 . Figure 3 indicates that by layer $L=-3$ the influence of
the contact with the insulator is minimal.

On the other hand, we find that contact with the metal has a substantially
greater effect on the insulator. Here the “Kondo proximity effect” diminishes
the long range AF order of the insulating lines adjacent to the interface.
This tendency to paramagnetism competes with the long range order induced by
the exchange energy $J\sim t^{2}/U$ deep in the insulator (Fig. 3a). To
quantify the proximity effects, we fit the curves of $S_{AF}(L,U,T)$ with the
hyperbolic form $S_{AF}(L,U,T)=a\tanh L/\lambda+b$. The penetration depth
$\lambda$ decreases with larger correlation $U$. Interestingly, $\lambda$
increases at higher temperature. It is possible that raising $T$ weakens the
magnetic order in the insulator, thereby enhancing the effects of the contact
with the metal.

Figure 4: (color online) The $L$ dependence of the density of states
$N(\omega)$ shows the evolution of the spectral properties across the metal-
insulator interface. Lines of sites within the Mott insulator ($L>0$ have
nonzero $U$) are characterized by the presence of Kondo resonance peaks split
by antiferromagnetic order at half-filling. However this AF gap between the
two peaks is dramatically weakened as $L\rightarrow 1$. On the other hand, the
densities of states associated with lines of sites on the metallic side of the
interface show little effect of contact with the insulator. Even for the line
$L=-1$ most immediately adjacent to the boundary there is only a small dimple
in $N(\omega)$ at $\omega=0$.

Figure 5: (color online) $L$-dependent density of states for different on-site
repulsions $U$ and temperatures $T$. The smooth penetration of metalicity into
the insulator and diminished metallic behaviors in the metal are clearly seen.
Similar to Fig. 3, we employ a hyperbolic fitting function and extract a
penetration depth $\lambda$ (insets), which we find decreases with larger
correlation $U$ and increases with higher temperature $T$.

##  IV Density of states

One criterion to distinguish metal from insulator is the single particle
density of states, which can be extracted from the local imaginary-time
dependent Green’s function
$G_{L}(\tau)=\sum\limits_{i\sigma}Tc_{iL\sigma}(\tau)c^{\dagger}_{iL\sigma}(0)$
from

$G_{L}(\tau)=\int_{-\infty}^{\infty}\
d\omega\frac{e^{-\omega\tau}}{e^{-\beta\omega}+1}\ N_{L}(\omega)$ (3)

and using the maximum entropy method maxent . The density of states is probed
experimentally with photoemission spectroscopy. In a translationally invariant
system one often also calculates the momentum dependent spectral function
$A(k,\omega)$, measured in angle-resolved photoemission spectroscopy (ARPES).

As mentioned in the introduction, a dramatic experimental achievement in the
past few years has been the development of standing wave ARPES gray10  which
has enabled the probing of strong correlation effects layer by layer in a
sample. Here we make some predictions for possible features to be seen in
SWARPES by obtaining the row dependent single particle density of states. As
with the row-by-row magnetic correlations, $N_{L}(\omega)$ shown in Fig. 4
characterizes the penetration of metalicity into the insulator. Kondo
proximity effects are evident in the evolution of $N_{L}(\omega)$. Figure 4
shows that the insulating lines for $L=6,5,4$, relatively far from the
interface, are characterized by the presence of a Kondo resonance peak split
by antiferromagnetic order at half-filling. However as the interface is
approached for $L=3,2,1$ the gap is increasingly filled in, evidence of the
coupling to the metallic half of the lattice. Meanwhile the metallic line
($L=-1$) immediately adjacent to the interface shows some influence of the
boundary, albeit a rather small one: the central peak at $\omega=0$ shows a
slight dip. This is completely gone for $L=-2$. Evidentally the metallic
behavior of $N(\omega)$ penetrates much further into the insulator than the
insulating physics does into the metal.

Further information on the spectral weight at Fermi level
$N_{L}(U,T)|_{\omega=0}$ as a function of on-site repulsion $U$ and
temperature $T$ is shown in Fig. 5. The penetration of metallic behavior into
the insulator and, conversely, the diminished spectral weight at the Fermi
surface in the metal are clearly seen. Similar to Fig. 3, we adopt a
hyperbolic function for fitting. The resulting penetration depths $\lambda$
are given in the inset, and decrease with larger correlation $U$ and lower
temperature $T$.

##  V Effect of variation of interface hopping

Figure 6: (color online) (a) Evolution of the antiferromagnetic (AF) structure
factors $S_{AF}(L)$ with interface hybridization $V$. Except the line of sites
$L=1$, adjacent to the metal, $S_{AF}(L>1)$ exhibits a nontrivial revival with
increasing the $V$. We attribute this affect to singlet formation at the
metal-insulator interface which leaves the lines $L>1$ decoupled from the
metal. (b) Further evidence of the strong magnetic coupling across the
interface. $U=0|4,\beta=10$.

We show in Fig. 6 the evolution of the line dependent AF structure factors
$S_{AF}(L)$ with interface hopping $V$. As expected, all lines $L>0$ with
nonzero on-site repulsion $U=4$ have large $S_{AF}$ when the metal and
insulator are decoupled ($V=0$). With increasing hybridization all
$S_{AF}(L>0)$ initially decrease. The most dramatic feature in Fig. 6(a) is
the different behaviors of $S_{AF}(L=1)$ and $S_{AF}(L>1)$. In the line $L=1$,
adjacent to the metal-insulator interface, $S_{AF}(L=1)$ monotonically
decreases with $V$. However, $S_{AF}(L>1)$ show nontrivial revivals with
increasing $V$ and ultimately recover to values characteristic of the
decoupled case $V=0$ when there is no contact with the metal whatsoever. This
behavior of $S_{AF}(L)$ is qualitatively similar to the multi-layer metal-
insulator interface in the three-dimensional Hubbard model euverte12 , and can
be attributed to “singlet” formation at the metal-insulator interface, which
both suppresses $S_{AF}(L=1)$ and also leaves the remaining lines $L>1$
decoupled from the metal.

The “metallic” line structure factors $S_{AF}(L)$ retain their small values
for the entire range of $V$. As with the real space spin correlation
functions, there is a greater effect of the metal on the insulator
(suppressing magnetic order) than of the insulator on the metal (inducing
magnetic order).

Figure 6(b) provides further evidence of the strong magnetic coupling across
the interface. The dominant feature is the rapid increase of antiferromagnetic
correlation across the metal-insulator interface $C_{\bot}(L=-1,L^{\prime}=1)$
as increasing the hopping $V$. Note also that magnetic correlation
$C_{\bot}(L=1,L^{\prime}=2)$ become smaller with larger hopping $V$ due to the
“singlet” formation at the interface $L=-1,1$.

##  VI Conclusion

We have used the numerically exact finite-temperature determinant Quantum
Monte Carlo (DQMC) method to study a metal-insulator interface in a two-
dimensional square lattice. By investigating the long-range antiferromagnetic
order and density of states, we demonstrated that the metallic behavior
penetrates into the insulator for several lattice spacings, with a penetration
depth $\lambda$ which decreases with increasing on-site repulsion, $U$, on the
insulator side of the interface. $\lambda$ is also (somewhat more weakly)
temperature dependent, decreasing as $T$ is lowered towards the critical
temperature $T_{c}=0$ for long range magnetic order on the insulating side.
The penetration length $\lambda$ thus shows an opposite temperature dependence
to that of the spin correlation length, $\xi$, which grows as $T$ is lowered.

The insulator-induced antiferromagnetic long-range order in the metal is
predominantly limited to the line immediately adjacent to the interface. That
is, magnetic characteristics of strong correlation appear to penetrate less
deeply into the metallic side of an interface than does weak correlation
physics penetrate into the insulator. This is consistent with previous results
for a planar interface in a 3D latticeeuverte12 . We note, however, that in
the 3D geometry, the effect of the contact with the insulator on the in-plane
conductivity $\sigma$ in the $U=0$ half of the lattice can extend beyond the
contact layer. This is consistent with our results for the density of states
at the Fermi level $N(\omega=0)$, which are shifted from their $U=0$ values
for layer $L=-2$ in addition to layer $L=-1$.

In the past several years, studies of the effect of spatially varying
densities and interaction strengths have been motivated by experiments on
trapped atomic gasesjaksch98 ; batrouni02 ; stoferle06 ; chin06 . In such
systems the spatial variation, eg a quadratic confining potential $V_{\rm
trap}(r)=V_{T}(r/l)^{2}$, has an explicit length scale $l$. This complicates
the determination of intrinsic length scales associated with the response of
the interacting fermions themselves. Our choice here of a sharp (scale free)
interface between metal and insulator ($U_{L}=0,L=-1,-2,-3...$ and
$U_{L}=U>0,L=1,2,3...$) allows us to attribute the lengths characterizing the
relaxation of properties on either side of the interface solely to the
fermionic correlations. This choice is, of course, also more appropriate to
the solid state context of a sharp interface between two materials.

One further motivation to study the metal-insulator interface is its close
relation to the question of “orbitally selective Mott transitions”
(OSMT)liebsch04 ; arita05 ; koga04 ; ferrero05 ; demedici05 ; knecht05 . Here
the central question is whether orbitals with different degree of electronic
correlation, coupled together by interorbital hybridization, necessarily
undergo the Mott transition simultaneously. The layer index ‘L’ in the
Hamiltonian considered in this paper bears a formal similarity to the orbital
degree of freedom in the OSMT, although of course the details of the coupling
via hybridization are rather different in the two cases. It is evident from
our data that a layer which shows the hallmarks of an antiferromagnetic Mott
insulator can coexist with layers which have the characteristic behavior of a
paramagnetic metal. That such coexistence is possible is similar to the
conclusion ultimately reached in numerical studies of the OSMT.

In this paper, the density of states and magnetic correlations near an
interface between $U=0$ and $U\neq 0$ regions have been found to be more or
less smooth interpolations between the bulk metal and Mott insulator. On the
other hand, Quantum Monte Carlo simulations have observed novel phasesnovel
such as spin liquids, to arise in models where the energy scales are poised at
the boundary between a semi-metal and an antiferromagnet. An interesting
extension of the present work would be, therefore, to bring such more general
regions into contact and study the properties at the interface.

##  VII Acknowledgements

We acknowledge support from the National Scence Foundation under grant NSF-
PIF-1005502. This work was also supported under ARO Award W911NF0710576 with
funds from the DARPA OLE Program and by the CNRS-UC Davis EPOCAL LIA joint
research grant. We are grateful to Bo Deans for useful input.

## References

  * (1) P. Zubko, S. Gariglio, M. Gabay, P. Ghosez, and J.M. Triscone, Annu. Rev. Condens. Matter Phys 2, 141 (2011), and references therein.
  * (2) J. Mannhart and D.G. Schlom, Science 327, 1607 (2010).
  * (3) K. Munakata, T.H. Geballe, and M.R. Beasley, Phys. Rev. B 84, 161405(R) (2011).
  * (4) J. W. Freeland, J. Chakhalian, A.V. Boris, J.M. Tonnerre, J.J. Kavich1, P. Yordanov, S. Grenier, P. Zschack1, E. Karapetrova, P. Popovich, H. N. Lee, and B. Keimer, Phys. Rev. B 81, 094414 (2010).
  * (5) A. Ohtomo, D. A. Muller, J.L. Grazul, and H.Y. Hwang, Nature 419, 378 (2002).
  * (6) S. Thiel, G. Hammerl, A. Schmehl, C.W. Schneider, and J. Mannhart, Science 313, 1942 (2006).
  * (7) A.X. Gray, C. Papp, B. Balke, S.-H. Yang, M. Huijben, E. Rotenberg, A. Bostwick, S. Ueda, Y. Yamashita, K. Kobayashi, E.M. Gullikson, J.B. Kortright, F.M.F. DeGroot, G. Rijnders, D.H.A. Blank, R. Ramesh, C.S. Fadley, Phys. Rev. B82, 205116 (2010).
  * (8) S. Okamoto and A. J. Millis, Nature (London) 428, 630 (2004); Phys. Rev. B70, 075101 (2004).
  * (9) S.S. Kancharla and E. Dagotto, Phys. Rev. B74, 195427 (2006).
  * (10) S. Okamoto and A. J. Millis, Phys. Rev. B70, 241104(R) (2004).
  * (11) R. W. Helmes, T. A. Costi, and A. Rosch, Phys. Rev. Lett. 101, 066802 (2008).
  * (12) H. Zenia, J. K. Freericks, H. R. Krishnamurthy, and T. Pruschke, Phys. Rev. Lett. 103, 116402 (2009).
  * (13) G. Borghi, M. Fabrizio, and E. Tosatti, Phys. Rev. B 81, 115134 (2010).
  * (14) Luis G.G.V. Dias da Silva, K.A. Al-Hassanieh, A.E. Feiguin, F.A. Reboredo2, and E. Dagotto Phys. Rev. B81, 125113 (2010).
  * (15) S. Yonemitsu, N. Maeshima, and T. Hasegawa, Phys. Rev. B76, 235118 (2007).
  * (16) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
  * (17) R. Blankenbecler, D.J. Scalapino, and R.L. Sugar, Phys. Rev. D24, 2278 (1981).
  * (18) S.R. White, D.J. Scalapino, R.L. Sugar, E.Y. Loh, Jr., J.E. Gubernatis, and R.T. Scalettar, Phys. Rev. B40, 506 (1989).
  * (19) E. Y. Loh, Jr., J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B41, 9301 (1990)
  * (20) A. Euverte, F. Hebert, S. Chiesa, R.T. Scalettar, and G.G. Batrouni, Phys. Rev. Lett. 108, 246401 (2012).
  * (21) The linear scaling in $\beta$ is restricted to cases when the matrix condition number does not become overly large: At large interaction strengths $U$ the necessity of recomputing the green’s function from scratch due to numerical instabilities can introduce $\beta^{2}$ terms. These can also arise in computation of the imaginary time dependent Greens function.
  * (22) ‘Long range order’ might be a bit of a misnomer. Although the correlations remain non-zero at long distances, it is likely this reflects an imposed order from the external contact with the insulator, as opposed to order mediated through the metallic layer itself.
  * (23) The increase of $S_{AF}$ at the two open boundaries is left out of our fitting. This phenomenon where “surface” magnetic correlations are larger than the “bulk”, despite the smaller number of neighbors, occurs because of the reduced quantum fluctuations (hopping) at the “surface”. See, for example, F. Zhang, S. Thevuthasan, R.T. Scalettar, R.R.P. Singh, and C.S. Fadley, Phys. Rev. B51, 12468 (1995), and references therein.
  * (24) Actually, magnetic correlations in the Hubbard model peak at $U\sim 8t$, ultimately declining when $U$ becomes very large owing to the decreasing exchange energy $J\sim t^{2}/U$ at large $U$.
  * (25) A. Sherman and N. Voropajeva, Int. J. Mod. Phys. B 24, 979 (2010).
  * (26) J.E. Gubernatis, M. Jarrell, R.N. Silver, and D.S. Sivia, Phys. Rev. B 44, 6011 (1991); A. W. Sandvik, Phys. Rev. B 57, 10287 (1998).
  * (27) D. Jaksch, C. Bruder, J.I. Cirac, C.W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * (28) G.G. Batrouni, V. Rousseau, R.T. Scalettar, M. Rigol, A. Muramatsu, P.J.H. Denteneer, and M. Troyer, Phys. Rev. Lett. 89, 117203 (2002).
  * (29) T. Stöferle, H. Moritz, K. Günter, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 96, 030401 (2006).
  * (30) J.K. Chin, D.E. Miller, Y. Liu, C. Stan, W. Setiawan, C. Sanner, K. Xu, and W. Ketterle, Nature (London) 443, 961 (2006).
  * (31) A. Liebsch, Phys. Rev. B70, 165103 (2004).
  * (32) R. Arita and K. Held, Phys. Rev. B72, 201102(R) (2005).
  * (33) A. Koga, N. Kawakami, T.M. Rice, and M. Sigrist, Phys. Rev. Lett. 92, 216402 (2004).
  * (34) M. Ferrero, F. Becca, M. Fabrizio, and M. Capone, Phys. Rev. B72, 205126 (2005).
  * (35) L. de’Medici, A. Georges, and S. Biermann, Phys. Rev. B72, 205124 (2005).
  * (36) C. Knecht, N. Blumer, and P.G.J. van Dongen, Phys. Rev. B72, 081103(R) (2005).
  * (37) An early example of these phenemena is the bond ordered wave phase between the spin density and charge density wave phases of the extended Hubbard model. See P. Sengupta, A.W. Sandvik, and D.K. Campbell, Phys. Rev. B65 155113 (2002). Recent spin liquid examples in Hubbard models are described in M. Aidelsburger, M. Itala, S. Nasciembéne, S. Trotzky, Y.A. Chen, and I. Bloch, Phys. Rev. Lett. 107, 255301 (2011); H. Morita, A. Watanabe, and M. Imada, J. Phys. Soc. Jpn. 71, 2109 (2002); Z. Meng, T. Lang, S. Wessel, F. Assaad, and A. Muramatsu, Nature (London) 464, 847 (2010); and Chia-Chen Chang and R.T. Scalettar, Phys. Rev. Lett. 109, 026404 (2012).