# Data-driven System Interconnections
and a Novel Data-enabled Internal Model Control

Yasaman Pedari, Jaeho Lee,
Yongsoon Eun, Hamid Ossareh  This work was supported by NASA under agreement
VT-80NSSC20M0213 and by National Science Foundation award CMMI-2238424. Y.
Pedari and H. Ossareh are with the Dept. of Elec. and Biomed. Eng., Univ. of
Vermont, Burlington, VT, USA (ypedari, h.ossareh@uvm.edu). J. Lee and Y. Eun
are with the Dept. of Elec. Eng. & Comp. Sci., DGIST, Daegu, South Korea
(yeun, jaeho.lee@dgist.ac.kr).

###### Abstract

Over the past two decades, there has been a growing interest in control
systems research to transition from model-based methods to data-driven
approaches. In this study, we aim to bridge a divide between conventional
model-based control and emerging data-driven paradigms grounded in Willems’
“fundamental lemma”. Specifically, we study how input/output data from two
separate systems can be manipulated to represent the behavior of
interconnected systems, either connected in series or through feedback. Using
these results, this paper introduces the Internal Behavior Control (IBC), a
new control strategy based on the well-known Internal Model Control (IMC) but
viewed under the lens of Behavioral System Theory. Similar to IMC, the IBC is
easy to tune and results in perfect tracking and disturbance rejection but,
unlike IMC, does not require a parametric model of the dynamics. We present
two approaches for IBC implementation: a component-by-component one and a
unified one. We compare the two approaches in terms of filter design,
computations, and memory requirements.

###### Index Terms:

Data-driven Control, Fundamental Lemma, Data-driven System Interconnections,
Internal Model Control

##  I Introduction

Recently, there has been a growing interest in data-driven control methods,
mainly due to a surge in available large datasets. These methods potentially
offer seamless end-to-end and flexible design, particularly for systems that
are difficult to model precisely. Following the trend towards data-driven
control, the work by Willems and coauthors on behavioral system theory (BST)
stands out [1]. They demonstrated that a sufficiently exciting trajectory of a
Linear Time-Invariant (LTI) system can be used to represent all its potential
trajectories. This concept termed as the “fundamental lemma”, highlights the
power of using data for learning and control.

In this study, we aim to bridge a divide between conventional model-based
control and emerging data-driven paradigms grounded in Willems’ fundamental
lemma. Specifically, we present how noise-free111Handeling noisy offline data
within the context of BST is an active area of research and is a topic for
future work. input/output data from multiple systems can be used to represent
the dynamics of an interconnected system in the BST framework. The existing
literature (e.g., [2, 3, 4]) has also studied the interconnection problem in
the context of BST, but not from an algorithmic standpoint, rather from a
purely theoretical standpoint as constraints on trajectories or latent
variables, or through parametric (i.e., state-space) representations of the
dynamics. Applications of [2, 3, 4] in controls have been reported in [5, 6,
7, 8]. Building on these foundations, the first contribution of this paper is
an algorithmic approach to convert individual system trajectories collected in
isolation into a data-driven representation of the interconnected system,
either in series or in feedback.

Using these data-driven interconnections, we propose a new control framework
called the Internal Behavior Control (IBC), which can be viewed as a data-
driven version of the classical Internal Model Control (IMC). The IMC uses a
model of the plant, as well as a causal approximate inverse, inside the
controller [9]. Here, we replace the model and the inverse with data-enabled
predictors. As we show, IBC inherits the desirable properties of the model-
based IMC: integral control, ease of tuning, and ability to naturally handle
Hammerstein-Weiner systems (e.g., actuator saturation). We offer two
approaches for implementing the proposed IBC. The first implementation is
component-by-component, where separate data-enabled predictors are leveraged
for the forward and inverse models. The second implementation combines the
data-driven forward and inverse predictors using the system interconnection
ideas discussed earlier. This effectively results in a single data-enabled
controller, which is more compact and more memory efficient than the first
implementation. The focus of this paper will be on single-input single-output
(SISO), stable, minimum-phase systems, though extension to the multi-input
multi-output case is not difficult.

Note that other model-based control strategies have also been recently viewed
under the lens of BST. Successful examples include the DeePC algorithm [10],
which stands parallel to the Model Predictive Control (MPC), and the design of
optimal Linear Quadratic Gaussian controllers directly from the subspace
matrices in [11]. These controllers may lack interpretability and may be hard
to tune. The proposed IBC does not suffer from these challenges. In summary,
this work mainly offers two contributions: $i$) methods and algorithms for
data-driven system interconnections; and $ii$) the IBC, which is a data-driven
control framework inspired by the classical IMC. To the best of our knowledge,
with the exception of the work [12] that studies a similar problem for control
of second-order Volterra systems, this is the first study of the IMC
controller in the context of BST.

This paper is organized as follows. Section II provides an overview of BST and
IMC. Data-driven interconnections and the IBC are presented in Sections III
and IV, respectively. The conclusions and future work are provided in Section
V.

The notation in this work is as follows. Given a system $G$, the integers
$n(G)$ and $L(G)$ denote the order and $L$-delay (i.e., the relative degree in
the SISO case), respectively. For a scalar signal $v(t)\in\mathbb{R}$, the
vector of values from time $t_{1}$ to $t_{2}$ is denoted by
$\mathbf{v}_{[t_{1},t_{2}]}=[v(t_{1}),\cdots,v(t_{2})]$. When there is no risk
of confusion, we omit the subscript and denote the vector by boldface letters,
i.e., $\mathbf{v}$. Given two vectors, $\mathbf{x}$ and $\mathbf{y}$, we
define col$(\mathbf{x},\mathbf{y})={[\mathbf{x}^{T}\;\;\mathbf{y}^{T}]}^{T}$.
For the vector $\mathbf{v}_{[t_{1},t_{2}]}$, the Hankel matrix of order $T$ is
defined as

$\mathscr{H}(\mathbf{v}_{[t_{1},t_{2}]})=\begin{bmatrix}v(t_{1})&v(t_{1}+1)&\cdots&v(t_{2}-T+1)\\\
v(t_{1}+1)&v(t_{1}+2)&\cdots&v(t_{2}-T)\\\ \vdots&\vdots&\ddots&\vdots\\\
v(t_{1}+T-1)&v(t_{1}+T)&\cdots&v(t_{2})\end{bmatrix}.$

##  II OVERVIEW OF BST AND IMC

###  II-A Overview of BST and Data-enabled predictions

This section provides a review of BST, for details see [13]. Consider a SISO
LTI system$\,G\,$with input $u\in\mathbb{R}$ and output $y\in\mathbb{R}$. For
any positive integer $T$, we define the “restricted behavior” of $G$, denoted
by$\,B|_{T}(G)$, as the set of all $T$-long
trajectories$\,\mathbf{w}=\text{col}(\mathbf{u},\mathbf{y})\in\mathbb{R}^{2T},\,$where$\,\mathbf{y}\in\mathbb{R}^{T}\,$is
the output trajectory corresponding to$\,\mathbf{u}\in\mathbb{R}^{T}$from some
initial condition.

Consider an offline-collected single $T_{d}$-long trajectory,
$\mathbf{w}^{d}=\text{col}(\mathbf{u}^{d},\mathbf{y}^{d})\in\mathbb{R}^{2T_{d}}$
of system $G$. Let $T_{p}$ be a positive integer to be later selected. We
construct the “data matrix”, $\mathcal{H}$, as

$\mathcal{H}=\begin{bmatrix}\mathcal{H}_{u}\\\
\mathcal{H}_{y}\end{bmatrix}\in\mathbb{R}^{2(T_{p}+1)\times(T_{d}-T_{p})}$ (1)

where

$\mathcal{H}_{u}=\mathscr{H}_{T_{p}+1}(\mathbf{u}^{d})=\begin{bmatrix}U_{p}\\\
U_{f}\end{bmatrix},\;\mathcal{H}_{y}=\mathscr{H}_{T_{p}+1}(\mathbf{y}^{d})=\begin{bmatrix}Y_{p}\\\
Y_{f}\end{bmatrix},$

with the partitioning $U_{p},Y_{p}\in\mathbb{R}^{T_{p}\times(T_{d}-T_{p})}$
and $U_{f},Y_{f}\in\mathbb{R}^{1\times(T_{d}-T_{p})}$. If $T_{p}\geq n(G)$ and
$\mathcal{H}$ satisfies the “low-rank condition”,

$\text{rank}(\mathcal{H})=T_{p}+1+n(G),$ (2)

then, any $(T_{p}+1)$-long trajectory of the system will belong to the column
space of $\mathcal{H}$222This result has become known as the “Generalized
Persistency of Excitation Condition” see Corollary 21 in [14]. This condition
formally involves the “lag” of the system, which in the SISO case is equal to
the order, $n$.. This condition allows us to use the data matrix for the
purpose of output prediction. Specifically, given the recent $T_{p}$ samples
of the input/output trajectory
($\mathbf{u}_{\text{ini}},\mathbf{y}_{\text{ini}}\in\mathbb{R}^{T_{p}}$), and
a “future” input ${u}_{\text{pred}}\in\mathbb{R}$, the unique future output,
$\mathbf{y}_{\text{pred}}\in\mathbb{R}$, can be predicted by solving:

$\mathcal{H}\mathbf{g}=\begin{bmatrix}\mathbf{u}_{\text{ini}}\\\
u_{\text{pred}}\\\ \mathbf{y}_{\text{ini}}\\\ y_{\text{pred}}\\\
\end{bmatrix}$ (3)

for the unknowns $\mathbf{g}$ and $y_{\text{pred}}$. The solution for vector
$\mathbf{g}$ is generally not unique. Often, the minimum-norm solution is
employed:

$\mathbf{g}^{*}=\left[\begin{array}[]{c}U_{p}\\\ U_{f}\\\
Y_{p}\end{array}\right]^{\dagger}\left[\begin{array}[]{c}\mathbf{u_{ini}}\\\
{u_{\text{pred}}}\\\ \mathbf{y_{ini}}\end{array}\right]$ (4)

and so

${y_{\text{pred}}}=Y_{f}\mathbf{g}^{*}$ (5)

###### Remark 1

Two assumptions in this review are not always necessary: First, the data
matrix does not require a single long trajectory; it can be formed from
independent trajectories; for details see [13]. Second, our predictors are
single-step predictors (i.e., $u_{\text{pred}},y_{\text{pred}}\in\mathbb{R}$).
However, longer prediction horizons may be used provided the data matrix
satisfies the low-rank condition.

###### Remark 2

Because of the use of a single step prediction (i.e.,
$y_{\text{pred}}\in\mathbb{R}$), the solution to (4) and (5) is equivalent to
the prediction obtained from a certainty-equivalent ARX model identified from
the offline data. For details, please see [15].

###  II-B Data-Enabled prediction of System Inverse

This section provides a review of data-enabled predictions for the system
inverse, as required for the IBC approach; for details see [16]. Given an
output trajectory, the goal is to compute the associated input using the
offline data. The ideas are similar to the previous section, except that we
predict the input instead of the output and $\mathcal{H}_{y}$ is now of a
different order. More specifically, the data matrix for the inverse system is
modified to

$\mathcal{H}_{\text{inv}}=\begin{bmatrix}\mathcal{H}_{u}\\\
\mathcal{H}_{y,\text{inv}}\end{bmatrix}\in\mathbb{R}^{(2T_{p}+2+L(G))\times(T_{d}-T_{p})}$
(6)

where $\mathcal{H}_{u}$ is as in (1) and $\mathcal{H}_{y,\text{inv}}$ is
defined as,

$\mathcal{H}_{y,\text{inv}}=\mathscr{H}_{T_{p}+L(G)+1}(\mathbf{y}^{d})=\begin{bmatrix}Y_{p}\\\
Y_{fL}\end{bmatrix}$

with the partitioning $Y_{p}\in\mathbb{R}^{T_{p}\times(T_{d}-T_{p})}$ and
$Y_{fL}\in\mathbb{R}^{(1+L(G))\times(T_{d}-T_{p})}$. Here, we have considered
different lengths for the offline data: $T_{d}$ for $u_{d}$ and $T_{d}+L(G)$
for $y_{d}$. The low-rank condition for the inverse case remains the same as
(2). If this condition is satisfied, $\mathcal{H}_{\text{inv}}$ can be used to
compute the input. To see this, at every timestep $t$, given the recent
$T_{p}+1+L(G)$ samples of the output and the same samples of the input without
the last $L(G)$ elements, we can uniquely compute the input $L(G)$ timesteps
in the past. Mathematically, given

$\mathbf{u}_{\text{ini}}^{\text{inv}}=\mathbf{u}_{[t-T_{p}-L(G),t-L(G)-1]},$
$\mathbf{y}_{\text{ini}}^{\text{inv}}=\mathbf{y}_{[t-T_{p}-L(G),t-L(G)-1]},\;\mathbf{y}_{\text{pred}}^{\text{inv}}=\mathbf{y}_{[t-L(G),t]},$

$\mathbf{u}_{\text{pred}}^{\text{inv}}=u(t-L(G))\in\mathbb{R}$ can be obtained
by solving:

$\begin{bmatrix}U_{p}\\\ U_{f}\\\ Y_{p}\\\
Y_{fL}\end{bmatrix}\mathbf{g}^{\text{inv}}=\begin{bmatrix}\mathbf{u}_{\text{ini}}^{\text{inv}}\\\
{u}_{\text{pred}}^{\text{inv}}\\\ \mathbf{y}_{\text{ini}}^{\text{inv}}\\\
{y}_{\text{pred}}^{\text{inv}}\end{bmatrix}$ (7)

for the unknowns
$\mathbf{g}^{\text{inv}}\in\mathbb{R}^{(T_{d}-T_{p})\times(2T_{p}+L(G)+2)}$
and ${u}_{\text{pred}}^{\text{inv}}$. As in the case of forward prediction,
the minimum-norm solution to $\mathbf{g}^{\text{inv}}$ is often employed,
which leads to:

${u}_{\text{pred}}^{\text{inv}}=U_{f}\begin{bmatrix}U_{p}\\\ Y_{p}\\\
Y_{fL}\end{bmatrix}^{\dagger}\begin{bmatrix}\mathbf{u}_{\text{ini}}^{\text{inv}}\\\
\mathbf{y}_{\text{ini}}^{\text{inv}}\\\
{y}_{\text{pred}}^{\text{inv}}\end{bmatrix}.$ (8)

###  II-C Overview of IMC

$\hat{G}^{-1}$$F$$G$$\hat{G}$$s_{2}$$r$$s_{1}$$u$$y$$\hat{y}$-+$e$- Figure 1:
Block diagram of the IMC structure

A block diagram of the IMC structure is shown in Fig. 1, where $G$ represents
the system, $\hat{G}$ is a model of the system, $\hat{G}^{-1}$ is the ideal
inverse of the model (to be revisited), and F is a filter to be defined. The
goal is to select $F$ to achieve a certain performance objective, e.g., to
minimize the $H_{2}$-norm of the tracking error for step setpoints. For open-
loop-stable minimum-phase systems, this $H_{2}$-norm is theoretically
minimized by setting $F=1$. In practice, $\hat{G}^{-1}$ is improper and cannot
be implemented. Moreover, even if it could be implemented, it would render the
feedback loop sensitive to uncertainties. This requires the augmentation of
$\hat{G}^{-1}$ with a low-pass filter whose order is at least the relative
degree of $\hat{G}$. If a description of $\hat{G}$’s uncertainties is
available, the time constant of the filter is selected to ensure robust
performance and stability [17], [18]. Otherwise, the time constant is treated
as a tuning parameter to directly trade-off performance and robustness,
rendering tuning an easy task. Note that input or output nonlinearities (i.e.,
Hammerstein-Wiener models [19]) can be naturally compensated for by including
them in the model path. In our study, we aim to substitute $\hat{G}$ and
$\hat{G}^{-1}$ with a data-enabled predictor while preserving IMC’s key
features, namely perfect steady-state tracking and disturbance rejection, ease
of tuning, and minimization of $H_{2}$ norm for step commands.

##  III DATA-DRIVEN SYSTEM INTERCONNECTIONS

In this section, we present how to create a data-driven representation of the
interconnection of two system – whether they are in series or in feedback –
using available input/output data from each individual system. As we show in
the next section, this allows us to view traditional model-based control
techniques, which frequently depend on system block interconnections, under
the lens of data-driven predictors.

As discussed in [8], in order to interconnect two systems, the output of one
system must match the input of the other. However, in general, the offline
input-output data may be collected in isolation, leading to mismatches. To
tackle this issue, we introduce the following function, the specific form of
which will be defined later:

$\mathbf{y}^{*}=\mathcal{Z}(\mathbf{u}^{*},\mathbf{u}^{d},\mathbf{y}^{d}):\mathbb{R}^{T}\times\mathbb{R}^{T_{d}}\times\mathbb{R}^{T_{d}}\xrightarrow{}\mathbb{R}^{T}$

In words, given an input sequence $\mathbf{u}^{*}$ , this function predicts
the output sequence, $\mathbf{y}^{*}$ , of the system starting from zero
initial conditions, while ensuring that
$\text{col}(\mathbf{u}^{*},\mathbf{y}^{*})$ is an admissible trajectory, i.e.,
compatible with$w^{d}=\text{col}(\mathbf{u}^{d},\mathbf{y}^{d})$. To this end,
we introduce a new input trajectory $\mathbf{u}^{\text{mod}}$, which is
identical to $\mathbf{u}^{*}$ except it is padded with $T_{p}$ zeros at the
beginning:

$\mathbf{u}^{\text{mod}}(t)=\begin{cases}0&\text{if }0\leq t<T_{p},\\\
\mathbf{u}^{*}(t-T_{p})&\text{if }T_{p}\leq t\leq T+T_{p}.\end{cases}$

We then define $\mathcal{Z}$ recursively, where $\mathbf{y}^{*}$ is computed
element by element. Specifically, using (4) and (5), ${y}^{*}(t)$is given by:

${y}^{*}(t)=y^{\text{mod}}(t+T_{P}),$

where $y^{\text{mod}}(t)$ is computed recursively as

$y^{\text{mod}}(t)=\begin{cases}0&\text{if }0\leq t<T_{p}\\\
Y_{f}\begin{bmatrix}U_{P}\\\ U_{f}\\\
Y_{p}\end{bmatrix}^{\dagger}\begin{bmatrix}\mathbf{u}^{\text{mod}}_{[t-T_{p},t-1]}\\\
{u}^{\text{mod}}(t)\\\
\mathbf{y}^{\text{mod}}_{[t-T_{p},t-1]}\end{bmatrix}&\text{if }T_{p}\leq t\leq
T+T_{p}\end{cases}$

Here, the Hankel matrices $U_{p}$, $U_{f}$ , $Y_{p}$, and $Yf$ are constructed
using $\mathbf{u}^{d}$ and $\mathbf{y}^{d}$, see (1). Note that $\mathcal{Z}$
is only defined when the data matrix of $\mathbf{w}^{d}$ has rank
$T_{p}+1+n(G)$ for a selected $T_{p}\geq n(G)$. The above ideas can be used to
“remove” the effect of initial conditions from the offline output
$\mathbf{y}^{d}$, as shown in the theorem below. The proof is straightforward
and omitted.

$G_{1}$$G_{2}$$\mathbf{u}_{1}$$\mathbf{y}_{1}$$\mathbf{y}_{s}$$G_{s}$ Figure
2: Series interconnection given data from individual systems. In general,
since $\mathbf{y}_{1}\neq\mathbf{u}_{2}$, the interconnection is not
achievable. In the figure, the successful interconnection is shown with
$G_{2}$ taking $\mathbf{y}_{1}$ as input and generating $\mathbf{y}_{s}$ as
output.

###### Theorem 1

Consider a SISO LTI system $G$ and let $T_{p}\geq n(G)$. Suppose a
$T_{d}$-long trajectory
$\mathbf{w}^{d}=\text{col}(\mathbf{u}^{d},\mathbf{y}^{d})$ of $G$, whose data
matrix has rank $T_{p}+1+n(G)$, is given. The trajectory
$\mathbf{w}_{0}=\text{col}(\mathbf{u}^{d},\mathcal{Z}(\mathbf{u}^{d},\mathbf{u}^{d},\mathbf{y}^{d}))$
is a trajectory of $G$ starting from zero initial conditions. Furthermore, if
$\mathbf{w}^{d}$ is collected starting from zero initial conditions, then
$\mathbf{w}_{0}=\mathbf{w}^{d}$.

We now use $\mathcal{Z}$ to study the data-driven interconnection of two
systems in series (see Fig. 2), where their corresponding offline data may be
of different lengths and may be collected in isolation. The following Theorem
provides a mechanism to construct an admissible trajectory of the series
interconnection directly from the offline data of the individual systems. Such
trajectory can then be used for data-driven prediction of the series
interconnection using the previously described methods.

###### Theorem 2

Consider
$\mathbf{w}_{1}=\text{col}(\mathbf{u}_{1},\mathbf{y}_{1})\in\mathscr{B}|_{T_{1}}(G_{1})$
and
$\mathbf{w}_{2}=\text{col}(\mathbf{u}_{2},\mathbf{y}_{2})\in\mathscr{B}|_{T_{2}}(G_{2})$
as $T_{1}$ and $T_{2}$ long trajectories of SISO LTI systems $G_{1}$ and
$G_{2}$, respectively. If the data matrix of trajectory $\mathbf{w}_{2}$ has
rank $T_{p}+1+n(G_{2})$ for $T_{p}\geq n(G_{2})$, then the trajectory define
by

$\mathbf{w}_{s}=\begin{bmatrix}\mathbf{u}_{1}\\\
\mathbf{y}_{s}\end{bmatrix},\;\mathbf{y}_{s}=\mathcal{Z}(\mathbf{y}_{1},\mathbf{u}_{2},\mathbf{y}_{2})$

satisfies $\mathbf{w}_{s}\in\mathscr{B}|_{T_{1}}(G_{s})$, where $G_{s}$
denotes the LTI system that results from connecting systems $G_{1}$ and
$G_{2}$ in series.

###### Proof:

To interconnect $G_{1}$ and $G_{2}$ in series, $G_{1}$’s output must match
$G_{2}$’s input. However, the two trajectories, $\mathbf{w}_{1}$ and
$\mathbf{w}_{2}$ generally do not satisfy this condition as they may have been
collected in isolation. To force the output of $G_{1}$ to be equal to the
input of $G_{2}$, we use $\mathcal{Z}$ to generate $G_{2}$’s output,
$\mathbf{y}_{s}$, given $G_{1}$’s output,$\,\mathbf{y}_{1}$, as input. The
trajectory$\,\text{col}(\mathbf{u}_{1},\mathbf{y}_{s})$ would then represent a
valid trajectory of the two systems in series, where the second system starts
from zero initial conditions. An illustration of the proof is given in Fig. 2.
∎

$G_{1}$$G_{2}$$\mathbf{u}_{1}-\check{\mathbf{y}}_{2}$$\mathbf{u}_{1}$$\mathbf{y}_{1}$$\check{\mathbf{y}}_{2}$+$G_{f}$
Figure 3: Positive feedback interconnection given data from individual
systems. In general, $\mathbf{y}_{1}\neq\mathbf{u}_{2}$, so the
interconnection is not achievable. In the figure, the successful
interconnection is shown with G2 taking $\mathbf{y}_{1}$ as input and
generating $\check{\mathbf{y}}_{2}$ as output. The interconnected feedback
system takes the signal from the subtraction
$\mathbf{u}_{1}-\check{\mathbf{y}}_{2}$ as its input and produces
$\mathbf{y}_{1}$ as its output.

We now study a similar result for the feedback interconnection of two systems
(see Fig. 3).

###### Theorem 3

Consider
$\mathbf{w}_{1}=\text{col}(\mathbf{u}_{1},\mathbf{y}_{1})\in\mathscr{B}|_{T_{1}}(G_{1})$
and
$\mathbf{w}_{2}=\text{col}(\mathbf{u}_{2},\mathbf{y}_{2})\in\mathscr{B}|_{T_{2}}(G_{2})$
as $T_{1}$ and $T_{2}$ long trajectories of SISO LTI systems $G_{1}$ and
$G_{2}$, respectively. If the data matrix of trajectory $\mathbf{w}_{2}$ has
rank $T_{p}+1+n(G_{2})$ for $T_{p}\geq n(G_{2})$, then the trajectory define
by

$\mathbf{w}_{f}=\text{col}(\mathbf{u}_{1}-\check{\mathbf{y}}_{2},\mathbf{y}_{1}),\;\check{\mathbf{y}}_{2}=\mathcal{Z}(\mathbf{y}_{1},\mathbf{u}_{2},\mathbf{y}_{2})$

satisfies $\mathbf{w}_{f}\in\mathscr{B}|_{T_{1}}(G_{f})$, where $G_{f}$
denotes the LTI system that results from connecting systems $G_{1}$ and
$G_{2}$ through positive feedback (shown in Fig. 3).

###### Proof:

In a positive feedback configuration, the output of $G_{1}$ becomes the input
of $G_{2}$. Concurrently, $G_{1}$’s input is derived by subtracting $G_{2}$’s
output from the reference. However, the two trajectories, $\mathbf{w}_{1}$ and
$\mathbf{w}_{2}$ generally do not satisfy these conditions as they may have
been collected in isolation. To set the output of $G_{1}$ equal the $G_{2}$’s
input, we will feed $\mathbf{y}_{1}$ into $G_{2}$ to generate
$\check{\mathbf{y}}_{2}=\mathcal{Z}(\mathbf{y}_{1},\mathbf{u}_{2},\mathbf{y}_{2})$.
Next, by subtracting $\check{\mathbf{y}}_{2}$ from $\mathbf{u}_{1}$, we derive
the input to the feedback system, while the output is represented by
$\mathbf{y}_{1}$. The trajectory
$\text{col}(\mathbf{u}_{1}-\check{\mathbf{y}}_{2},\mathbf{y}_{1})$ would then
represent a valid trajectory of the two systems in positive feedback
interconnection. An illustration of the proof is given in Fig. 3. ∎

Note that both Theorems 2 and 3 involve regenerating the output of $G_{2}$,
which requires $\mathbf{w}_{2}$ to satisfy the low-rank condition.
Alternatively, the input and/or output of $G_{1}$ can be regenerated, in which
case w1 needs to satisfy the low-rank condition instead. Also, note that
Theorems 2 and 3 can also be easily extended to negative feedback, parallel
configurations, and more complicated structures. In the following section, we
will use these results to develop a data-driven parallel to IMC, called IBC.

##  IV INTERNAL BEHAVIOR CONTROL (IBC)

Referring back to Fig. 1, this section introduces the IBC, a data-driven
parallel to IMC. We present two IBC approaches: the first is referred to as
component-by-component (CBC), where we replace each block in Fig. 1 with a
data-driven predictor. The second is referred to as the unified approach,
where we use the results in Section III to combine the data-driven predictors
into a single one.

###  IV-A Component-by-Component (CBC) Approach

In the CBC approach, our objective is to replace the model-based blocks in the
IMC structure, specifically $\hat{G}^{-1}$ and $\hat{G}$, with data-driven
predictors. For this, we leverage the data-driven forward and inverse
predictors, presented in Section II-A and II-B, respectively.

####  IV-A1 CBC-IBC Forward Model

Given a $T_{d}$-long trajectory
$\mathbf{w}^{d}=\text{col}(\mathbf{u}^{d},\mathbf{y}^{d})$ of $\hat{G}$, whose
data matrix has rank $T_{p}+1+n(G)$ for $T_{p}\geq n(G)$, we can replace the
forward model of IMC, i.e., $\hat{G}$, with the data-enabled predictor given
in(4)–(5).

####  IV-A2 CBC-IBC Inverse Model

In the original IMC setup, an inverse model of the system, $\hat{G}^{-1}$ , is
also required to “cancel out” the system dynamics, effectively making the
controlled system (without plant/model mismatch) behave like an identity map
from $r$ to $y$. In the CBC-IBC setup, we leverage the inverse predictor
reviewed in Section II-B to replace $\hat{G}^{-1}$ with its data-driven
counterpart. The inverse predictor uses the same offline data as the forward
predictor. While the length of the offline input data remains unchanged, the
inverse predictor requires that the output sequence be $L(G)$ steps longer. In
contrast to the model-based IMC, the inverse predictor in IBC is already
causal and implementable because of the delay of $L(G)$ included in the
prediction, see Eq. (8), where $u^{pred}_{\text{inv}}=u(t-L(G))$

####  IV-A3 CBC-IBC Filter

In the traditional IMC framework, a low-pass filter with unity DC gain is
often employed. A choice is

$F(z)=\frac{1}{(\frac{\tau}{T_{s}}z+(1-\frac{\tau}{T_{s}}))^{L}}$ (9)

where $L$ and $T_{s}$ denote the $L$-delay of system $G$ and the sampling
period, respectively. This choice of order for the filter ensures that the
series interconnection of the filter with the system inverse becomes causal
and thus implementable. The scalar $\tau\in(0,1)$ is a free parameter, which
can be used to trade-off performance and robustness. Smaller $\tau$ values
yield faster but less stable responses, while larger values make the response
more robust but slower. In the CBC-IBC setup, the filter design differs from
that in IMC and the unified approach, discussed in the next subsection. This
is because a delay of $L(G)$ has already been incorporated into the inverse
system, as discussed above, and so the data-driven inverse is already
implementable. Algorithm 1 describes the controller operation in the CBC-IBC
at each timestep of the control loop. For details on the signal names, please
refer to Fig. 1.

Algorithm 1 CBC-IBC controller implementation

1:Setup: Choose $T_{p}\geq n(G)$. Make Hankel matrices $U_{p}$, $U_{f}$,
$Y_{p}$, $Y_{f}$, and $Y_{fl}$ through (1) and (7) ensuring that both the
forward and inverse data matrices have rank $T_{p}+1+n(G)$. Choose $\tau$ for
the filter $F$.

2:External Inputs: $r(t)$

3:Inputs from the previous loop: $\mathbf{u}_{[t-T_{p},t-1]}$,
$\check{\mathbf{y}}_{[t-T_{p},t-1]}$, $\mathbf{s}_{1,[t-T_{p}-L(G),t-1]}$,
$\mathbf{s}_{2,[t-T_{p}-L(G),t-L(G)-1]}$

4:Outputs: Control command $u(t)$

5:Compute $e(t)=y(t)-\hat{y}(t-1)$

6:Compute $s_{1}(t)=r(t)-e(t)$

7:Data-enabled prediction of System Inverse:
$s_{2}(t)=U_{f}\begin{bmatrix}Y_{p}\\\ U_{p}\\\
Y_{fL}\end{bmatrix}^{\dagger}\begin{bmatrix}\mathbf{s}_{1,[t-T_{p}-L(G),t-L(G)-1]}\\\
\mathbf{s}_{2,[t-T_{p}-L(G),t-L(G)-1]}\\\
\mathbf{s}_{1,[t-L(G),t]}\end{bmatrix}$

8:Compute, $u(t)$, by filtering $s_{2}(t)$ through $F$

9:Compute the date-enabled prediction (for the next loop):
$\hat{y}(t)=Y_{f}\begin{bmatrix}U_{p}\\\ Y_{p}\\\
U_{f}\end{bmatrix}^{\dagger}\begin{bmatrix}\mathbf{u}_{[t-T_{p},t-1]}\\\
\hat{\mathbf{y}}_{[t-T_{p},t-1]}\\\ u(t)\end{bmatrix}$

The properties of CBC-IBC are discussed in the following theorem.

###### Theorem 4

Suppose $G$ is an open-loop-stable minimum-phase SISO LTI system. The CBC-IBC
controller using filter $z^{L}F(z)$ performs identically to the classical IMC
controller using filter $F(z)$. Furthermore, the IBC controller results in
perfect steady-state tracking and disturbance rejection, and minimizes the
$H_{2}$ norm of the tracking error for step commands.

###### Proof:

We know that the model-based forward predictions in IMC are identical to their
data-driven counterparts in IBC. For the inverse prediction, the two differ
by$\,L(G)$ delays as discussed earlier. The statement of the theorem thus
follows. Finally, the tracking and disturbance rejection properties
and$\,H_{2}$ optimality are directly inherited from the model-based IMC. ∎

###  IV-B Unified Approach

It is known that the IMC controller can be represented as a single LTI
controller, $C$, as shown in Fig. 4. This controller is given by
$C(s)=\frac{\hat{G}^{-1}F}{1-F}$ , which has the same order as $G$, i.e.,
$n(C)=n(G)$. In this section, we apply the tools developed in section III to
create a unified IBC approach that uses a single data-driven predictor to
replace $C$, thus eliminating the necessity of separate data-enabled
predictors for the forward and inverse systems, as was required in the CBC
approach. The unified approach thus yields a more compact controller compared
to CBC-IBC. We now apply the results of Section III to create a data-driven
representation of $C$, which we achieve by finding an admissible trajectory
for $C$ from the input-output data of $G$. The results are summarized in the
following theorem.

$\hat{G}^{-1}$$F$$G$$\hat{G}$$r$$s_{3}$$u$$y$+-$C$ Figure 4: Simplified block
diagram of the IMC structure

###### Theorem 5

Given

  * •

$T_{d}$-long input-output data $\mathbf{u}^{d}$ and $\mathbf{y}^{d}$ from the
plant, $G$, where the data matrix has rank $T_{p}+1+n(G)$ for $T_{p}\geq
n(G)$,

  * •

and filter $F$ in (9),

then the trajectory defined by

$\mathbf{w}_{c}=\text{col}(\mathbf{y}^{d}-\bar{\mathbf{y}}^{d},\bar{\mathbf{u}}^{d}),\;\bar{\mathbf{y}}^{d}=f*\mathbf{y}^{d},\;\bar{\mathbf{u}}^{d}=f*\mathbf{u}^{d},$

satisfies $\mathbf{w}_{c}\in\mathscr{B}|_{T_{d}}(C)$, where $f$ is the impulse
response of the filter $F$, and ‘$*$’ denotes the convolution operator.

###### Proof:

The proof directly stems from Theorems 2 and 3 with the input-output sequences
for the inverse system set to$\,\mathbf{y}^{d}\,$and$\,\mathbf{u}^{d}$. ∎

Now, if the data matrix of $\mathbf{w}_{c}$ has rank $T_{p}+1+n(C)$ for
$T_{p}\geq n(C)$, then this data matrix can be used in the IBC unified
approach to make data-enabled predictions of $C$. Note that $n(C)=n(G)$.
Specifically, the data matrix of $\mathbf{w}_{c}$ is defined as,

$\mathcal{H}_{C}=\begin{bmatrix}\mathcal{H}_{C,u}\\\
\mathcal{H}_{C,y}\end{bmatrix}\in\mathbb{R}^{2(T_{p}+1)\times(T_{d}-T_{p})}$
(10)

where

$\mathcal{H}_{C,u}=\mathscr{H}_{T_{p}+1}(\mathbf{y}^{d}-\bar{\mathbf{y}}^{d})=\begin{bmatrix}E_{p}\\\
E_{f}\end{bmatrix},\;\mathcal{H}_{C,y}=\mathscr{H}_{T_{p}+1}(\bar{\mathbf{u}}^{d})=\begin{bmatrix}F_{p}\\\
F_{f}\end{bmatrix},$ (11)

with the partitioning $E_{p},F_{p}\in\mathbb{R}^{T_{p}\times(T_{d}-T_{p})}$
and $E_{f},F_{f}\in\mathbb{R}^{1\times(T_{d}-T_{p})}$. Data-enabled
predictions based on $\mathcal{H}_{C}$ will be used directly in the unified
IBC approach. Theorem 4 also applies to the unified approach, guaranteeing an
identical performance to the classical IMC controller using the same filter.
To conclude this section, Algorithm 2 presents the unified IBC controller at
each timestep of the control loop. For details of the signal names, please see
Fig. 4.

Algorithm 2 Unified-IBC controller implementation

1:External Inputs: $r(t)$

2:Inputs from the previous loop:
$\mathbf{u}_{[t-T_{p},t-1]}$,$\mathbf{s}_{3,[t-T_{p},t-1]}$

3:Outputs: Control command $u(t)$

4:Setup: Set $T_{p}\geq n(G)$. Filter the input-output sequence
$\mathbf{u}^{d}$ and $\mathbf{y}^{d}$ through $F$ to obtain
$\bar{\mathbf{u}}^{d}$ and $\bar{\mathbf{y}}^{d}$, respectively. Then, using
the difference ${\mathbf{y}}^{d}-\bar{\mathbf{y}}^{d}$ as input and
$\bar{\mathbf{u}}^{d}$ as output, generate the Hankel matrices $E_{P}$,
$E_{f}$, $F_{P}$, and $F_{f}$ through Eq. (11), ensuring that the resulting
data matrix has rank $T_{p}+1+n(G)$.

5:Compute $s_{3}(t)=r(t)-y(t)$

6:Compute and return$\,u(t)=F_{f}\begin{bmatrix}E_{p}\\\ F_{p}\\\
E_{f}\end{bmatrix}^{\dagger}\begin{bmatrix}\mathbf{s}_{3,[t-T_{p},t-1]}\\\
\mathbf{u}_{[t-T_{p},t-1]}\\\ \mathbf{s}_{3,[t-T_{p},t]}\end{bmatrix}$

###### Remark 3

We recall from [17] that IMC and constraint-free MPC are essentially identical
in the sense that for each tuning of the MPC controller, there exists a unique
IMC filter, F, that results in the same control actions. A similar
relationship exists between their data-driven counterparts. Specifically, just
as IBC serves as a data-enabled parallel of IMC, DeePC [10] is a data-driven
parallel of MPC. Therefore, it is reasonable to infer that the same essential
equivalence exists between IBC and DeePC when no constraints are involved and
a similar performance is to be expected from both methods in the noise-free
setup.

###### Remark 4

The unified approach has less memory requirements than CBC, in terms of both
online and offline data. Regarding online data, CBC requires solving two data-
enabled predictions, one for the forward and one for the inverse model. The
forward prediction requires a minimum memory of the past $n(G)$ input-output
samples, while the inverse requires the past $n(G)+L(G)$ samples. The unified
approach, on the other hand, solves for only one data-enabled prediction in
every control loop but requires memory of the past $n(G)$ samples. Therefore,
the unified requires $n(G)+L(G)$ fewer samples and is thus more memory
efficient in terms of online data. Regarding offline data, CBC requires a
minimum of $T_{d}=3n+1+L$ to satisfy the low-rank condition, while unified
requires $T_{d}=3n+1$, so the unified approach requires $L$ fewer offline data
samples

##  V SIMULATION RESULTS

In this section, we illustrate the proposed IBC controller using an example.
Consider the second-order plant, $G(s)=\frac{10(s+1)}{(s+2)(s+4)}$. The
dynamics are discretized using the zero-order hold method with a sampling
period of $0.01$ s. The time constant, $\tau$, of the filter is set to $0.5$
s. In the CBC approach, the filter is advanced by a time shift of $L=1$ to
avoid unnecessary delays (consistent with Theorem 4). For the IBC controller
design, we collect offline data by feeding the system a random input sequence.
The collected output, $\mathbf{y}^{d}$, as well as the filtered input and
output signals used in the unified approach, are displayed in Fig. 5 Since
$n(G)=2$, $T_{p}$ is chosen as 2 in both approaches. The number of offline
data samples, $T_{d}$, is chosen as the minimum: $T_{d}=8$ for CBC, and
$T_{d}=7$ (see Remark 4). Both algorithms achieve perfect tracking as depicted
in Fig. 6. A step disturbance is applied to the plant input around time $13$
s. Both IBC approaches successfully recover from this disturbance. Note that
while both approaches yield identical performance, the unified approach is
more data efficient and solves a single data-enabled prediction in each
iteration, in contrast to the CBC approach, which solves two.

Figure 5: Noise-free raw and filtered input/output data collected offline.
Because the input is random, all Hankel matrices satisfy their respective low-
rank condition. Figure 6: Illustration of the effectiveness of IBC in both CBC
and unified approaches in reference tracking and disturbance rejection.

##  VI CONCLUSIONS AND FUTURE WORK

This paper first presented a data-driven formulation of series and feedback
interconnection of dynamical systems, where we described how data-driven
predictors can be designed for the interconnected system using data from each
individual system. Using these ideas, we introduced the Internal Behavior
Control (IBC), a new data-driven control strategy based on the well-known
Internal Model Control (IMC). Similar to IMC, the IBC is easy to tune and
results in perfect tracking and disturbance rejection but, unlike IMC, does
not require a parametric model. We outlined two approaches for IBC
implementation: CBC and Unified. We comapred the two approaches in terms of
filter design, computations, and memory requirements, where we showed that
unified is more memory efficient. There are numerous avenues for future
research. First, the focus of this paper was on open-loop stable, minimum-
phase SISO LTI systems. Future work will relax these assumptions to study
unstable, non-minimum phase, and/or MIMO systems in the context of IBC.
Second, we assumed here that the offline data is free of “noise” - measurement
noise, process disturbance, or nonlinearities, all of which make the data
inconsistent with an LTI structure. Investigation of noise in the offline data
and its impact on the IBC performance is another topic of future work. Third,
in traditional IMC, the filter is tuned using the knowledge of process
uncertainties. Future work will investigate uncertainty quantification of the
data-driven predictors and IBC filter tuning. Finally, we did not offer
specific guidelines on how to design the offline data to ensure that the data
matrices for the series and feedback interconnections satisfy their respective
rank conditions. Input design is another interesting topic for future work.

## References

  * [1] Jan C Willems and Jan W Polderman.  Introduction to mathematical systems theory: a behavioral approach, volume 26.  Springer Science & Business Media, 1997.
  * [2] Jan C Willems.  The behavioral approach to open and interconnected systems.  IEEE control systems magazine, 27(6):46–99, 2007.
  * [3] Jan C Willems.  Paradigms and puzzles in the theory of dynamical systems.  IEEE Transactions on automatic control, 36(3):259–294, 1991.
  * [4] Jan C Willems.  Models for dynamics.  Dynamics reported: a series in dynamical systems and their applications, pages 171–269, 1989.
  * [5] Thabiso Maupong and Paolo Rapisarda.  Data-driven control: the full interconnection case.  In 22nd International Symposium on Mathematical Theory of Networks and Systems, July 2016.
  * [6] Thabiso M Maupong and Paolo Rapisarda.  Data-driven control: A behavioral approach.  Systems & Control Letters, 101:37–43, 2017.
  * [7] Harry L Trentelman and Jan C Willems.  Control in a behavioral setting.  In Proceedings of 35th IEEE Conference on Decision and Control, volume 2, pages 1824–1829. IEEE, 1996.
  * [8] Jan C Willems.  On interconnections, control, and feedback.  IEEE Transactions on Automatic control, 42(3):326–339, 1997.
  * [9] Carlos E Garcia and Manfred Morari.  Internal model control. a unifying review and some new results.  Industrial & Engineering Chemistry Process Design and Development, 21(2):308–323, 1982.
  * [10] Jeremy Coulson, John Lygeros, and Florian Dörfler.  Data-enabled predictive control: In the shallows of the deepc.  In 2019 18th European Control Conference (ECC), pages 307–312. IEEE, 2019.
  * [11] Wouter Favoreel.  Subspace methods for identification and control of linear and bilinear systems.  PhD thesis, Katholiek Universiteit Leuven, 1999.
  * [12] Juan G Rueda-Escobedo and Johannes Schiffer.  Data-driven internal model control of second-order discrete volterra systems.  In 59th IEEE Conference on Decision and Control (CDC), pages 4572–4579. IEEE, 2020.
  * [13] Ivan Markovsky and Florian Dörfler.  Identifiability in the behavioral setting.  homepages. vub. ac. be/imarkovs/publications/identifiability. pdf, 2020.
  * [14] Ivan Markovsky and Florian Dörfler.  Identifiability in the behavioral setting.  IEEE Transactions on Automatic Control, 68(3):1667–1677, 2022.
  * [15] Florian Dörfler, Jeremy Coulson, and Ivan Markovsky.  Bridging direct and indirect data-driven control formulations via regularizations and relaxations.  IEEE Transactions on Automatic Control, 68(2):883–897, 2022.
  * [16] Yongsoon Eun, Jaeho Lee, and Hyungbo Shim.  Data-driven inverse of linear systems and application to disturbance observers.  In 2023 American Control Conference (ACC), pages 2806–2811. IEEE, 2023.
  * [17] Manfred Morari.  Robust process control.  Chemical engineering research & design, 65(6):462–479, 1987.
  * [18] Daniel E Rivera, Manfred Morari, and Sigurd Skogestad.  Internal model control: Pid controller design.  Industrial & engineering chemistry process design and development, 25(1):252–265, 1986.
  * [19] LA Aguirre, MCS Coelho, and MV Correa.  On the interpretation and practice of dynamical differences between hammerstein and wiener models.  IEE Proceedings-Control Theory and Applications, 152(4):349–356, 2005.


[1]organization = Department of Industrial Engineering and Operations
Research, addressline=IIT Bombay, city=Mumbai, Maharashtra, postcode=400076,
country=India [2]organization = Department of Mechanical Engineering,
addressline=IIT Bombay, city=Mumbai, Maharashtra, postcode=400076,
country=India

# Online Universal Dirichlet Factor Portfolios

Purushottam Parthasarathy<EMAIL_ADDRESS>Avinash Bhardwaj
<EMAIL_ADDRESS>Manjesh K. Hanawal<EMAIL_ADDRESS>
###### Abstract

We revisit the online portfolio allocation problem and propose universal
portfolios that use factor weighing to produce portfolios that out-perform
uniform dirichlet allocation schemes. We show a few analytical results on the
lower bounds of portfolio growth when the returns are known to follow a factor
model. We also show analytically that factor weighted dirichlet sampled
portfolios dominate the wealth generated by uniformly sampled dirichlet
portfolios. We corroborate our analytical results with empirical studies on
equity markets that are known to be driven by factors.

###### keywords:

universal portfolio , portfolio optimization , online portfolio selection ,
dirichlet portfolio

##  1 Introduction

Portfolio selection is an important and well studied problem across the
disciplines of finance, economics, operations research and computer science.
While there are many flavors to approaching this problem, two major streams of
literature exist:

  * 1.

Single or Multi-period optimization (The best we can do in a single or multi
period from a mean-variance perspective)

  * 2.

Online portfolio selection (The best in terms of growth rate that we can do
over the long run)

Conceptually there is no strong argument to separate the two approaches, and
there is indeed some literature that combines the two[10]. In this work we
explore the state of the art for online portfolio selection and propose a
simple heuristic for asset allocation in factor driven markets that provides
an effective alternative to some of the complex techniques that have been
proposed in the literature.

###  1.1 Literature review

Cover[5] shows that one can come up with an adapted strategy for online
portfolio selection that asymptotically outperforms the best stock in the
portfolio. The portfolio that achieves this is defined as the performance
weighted average of all the possible portfolios in the n-simplex, assuming n
assets. In fact, Cover’s[5] central result is that the exponential growth rate
of a ‘non-anticipating’ fully adapted portfolio with no-lookahead bias (one
that is proposed) asymptotically approaches the growth rate of a constant re-
balanced portfolio with the benefit of hindsight. The update in each step
tilts towards the portfolios that have performed well so far and continuously
tracks a weighted average portfolio. Note that as [5] restricts the benchmark
portfolio to constant re-balanced portfolios, it becomes clear that although
[5]’s choice of a benchmark is reasonable within the context of what is
presented, without this restriction of constant re-balancing a portfolio with
the advantage of hindsight can achieve significantly higher returns[14]. In
subsequent work done as an extension to [5], Cover and Ordentlich[7] find a
tight error bound on the maximum deviation between the best constant re-
balanced portfolio in hindsight and the one they arrive at using side
information without the benefit of hindsight. Side information is modelled in
a discrete fashion here, selected from a set of k possible values. Based on
past association of a particular discrete value to a best portfolio, one can
put forward a best portfolio whenever k is revealed as side information again.
In related work [11] show that exponential multiplicative updates yield almost
the same performance as the universal portfolio and yet have less memory
requirements. In [11] the authors use a utility function that serves to trade-
off a performance measure and a distance measure. The performance is based on
$x_{t}$, the current relative price vector and the distance measure is chosen
to be the relative entropy between the new weight vector and the current one.
An approximation of the optimization objective function along with a fully
invested constraint leads to a recursive update algorithm that uses an
exponential update to the weight vector. In various empirical tests, the
algorithm achieves better out of sample performance than Cover’s[5] portfolio
and some other benchmarks.

As a sharp contrast to the universal portfolio type idea where best assets are
assumed to perform well in the future, some approaches use a contrarian
approach. The authors in [4] show that a forecast step that uses mean-
reversion can be used to get a performant portfolio. The authors calculate a
transfer(i $\to$ j) which depends upon a few factors that mimic mean reversion
and use lead-lag relationships. The authors show that this heuristic does well
empirically on a few datasets.

There is also a set of work that has been done that explores meta-algorithms
(MAs). The authors in [8] use daily S&P500 data for the past 21 years to show
that MAs outperform existing portfolio selection algorithms by several orders
of magnitude. As a good summary of available work, Li et al. [13] have written
a comprehensive review paper that describes many of the major techniques in
this area.

##  2 Constant re-balanced portfolios in factor driven markets

###  2.1 Constant re-balanced portfolios for a single factor market

In this section we derive a relationship between the growth rate of a constant
re-balanced portfolio in a market that is driven by a single factor. This
derivation serves to setup some of the preliminaries for results we show
later. Through the body of this work, there are m assets, n periods of
trading, and $b_{1}...b_{m}$ are portfolio allocations that drive returns of
the portfolio. The symbol t stands for transpose and not time. $X_{i}$ is the
$i^{th}$ relative price vector. A relative price is defined as the ratio of
the price in the next time period to the current one. The long term average
growth rate of a portfolio is given as:

$\frac{1}{n}\log(\prod_{i=1}^{n}b^{t}X_{i})$

Consider a single factor model as:

$\begin{bmatrix}X_{1}^{i}\\\ X_{2}^{i}\\\ .\\\ .\\\ X_{m}^{i}\\\
\end{bmatrix}=\begin{bmatrix}\beta_{1}^{t}\;R_{m}^{i}\\\
\beta_{2}^{t}\;R_{m}^{i}\\\ .\\\ .\\\ \beta_{m}^{t}\;R_{m}^{i}\\\
\end{bmatrix}=R_{m}^{i}\;\boldsymbol{\beta}$

where $R_{m}^{i}$ is a scalar that denotes single-factor (market) returns for
the $i^{th}$ period , and $\beta_{k}$ is a quantity that measures the factor
exposure of the $k^{th}$ asset to the factor return. Using the above
expression into the definition of the growth rate we have:

$G=\frac{1}{n}\log(\prod_{i=1}^{n}b^{t}R_{m}^{i}\boldsymbol{\beta})$

$G=\frac{1}{n}\log(\prod_{i=1}^{n}R_{m}^{i}\prod_{i=1}^{n}b^{t}\boldsymbol{\beta})$

$G=\frac{1}{n}\left[\log(\prod_{i=1}^{n}R_{m}^{i})+n\;\log
b^{t}\boldsymbol{\beta}\right]$

$G=\log
b^{t}\boldsymbol{\beta}+\frac{1}{n}\left[\log(\prod_{i=1}^{n}R_{m}^{i})\right]=\log
b^{t}\boldsymbol{\beta}+G_{m}$

which shows us that the excess average growth rate over the market growth rate
$G_{m}$ is given by the first term in the above expression. The maximum growth
rate is achieved at the vertex of the simplex which corresponds to the highest
$\beta$ asset, which is equivalent to holding investment at the asset with the
highest $\beta$.

###  2.2 Growth Rate of multi-factor portfolios

In this section we derive the relationship between the growth of the market
and the growth of a constant re-balanced portfolio when the portfolio
constituents can be assumed to follow a log-factor model. A stock market as
modelled in Cover[6] is represented as a vector of price-relatives that are
assumed to follow a distribution $F(x)$. Define $X_{i}$ as the one-period
price-relative for asset $i$, which is equal to $p^{end}_{i}/p^{start}_{i}$
where $p_{i}$’s are measured at the start and end of the time period (single
period model). A portfolio $b=(b_{1},b_{2},b_{3}...b_{m}),b_{i}\geq 0,\Sigma
b_{i}=1$ is the allocation of wealth to assets. The wealth that is generated
by engaging $b$ over the single period price-relatives is $S=b^{t}X$. As the
investment model involves reinvesting of the capital every period, we are
interested in the cumulative product of these single period gains and thereby
characterize the expected logarithm of this gain[6] as:

$W(b,F)=\int\log b^{t}x\;dF(x)=E\left[\log b^{t}X\right]$ (1)

The optimal $b^{*}$ that maximizes (1) is also the doubling rate of the
portfolio as shown in [6]. We now consider a 2-factor log returns factor model
of the market where known factors drive the growth of a market. Using this
setup, we derive the relationship between the growth of a constant re-balanced
portfolio as described above and the factor growth rates. For literature on
factor based models in general see [9] and [12] particularly for momentum
factor related effects. The factor model for asset $j$ we use is:

$R_{ji}=(r_{1i})^{\beta^{1}_{j}}(r_{2i})^{\beta^{2}_{j}}\quad j\in(1...m)$

where $r_{1i}$ and $r_{2i}$ are returns attributed to factor 1 and factor 2
and $\beta$ denotes the factor exposures. The average growth rate of this
portfolio is given by:

$GR_{avg}=\frac{1}{n}\,E\,\sum_{i=1}^{n}\log(\sum_{j=1}^{m}b_{j}R_{j})$

$GR_{avg}=\frac{1}{n}\,E\,\sum_{i=1}^{n}\log(\sum_{j=1}^{m}b_{j}r_{1i}^{\beta^{1}_{j}}r_{2i}^{\beta^{2}_{j}})$

By the AM-GM inequality and by the property that log is a monotone function:

$GR_{avg}\geq\frac{1}{n}\,E\,\sum_{i=1}^{n}\log(m\prod_{j=1}^{m}(b_{j}r_{1i}^{\beta^{1}_{j}}r_{2i}^{\beta^{2}_{j}})^{\frac{1}{m}})$

$GR_{avg}\geq\log
m+\frac{1}{mn}E\left[\sum_{i=1}^{n}\sum_{j=1}^{m}\log(b_{j}r_{1i}^{\beta^{1}_{j}}r_{2i}^{\beta^{2}_{j}})\right]$

$GR_{avg}\geq\log m+\frac{1}{mn}E\left[\sum_{i=1}^{n}\sum_{j=1}^{m}\log
b_{j}+\beta^{1}_{j}\log r_{1i}+\beta^{2}_{j}\log r_{2i}\right]$

$\displaystyle GR_{avg}$ $\displaystyle\geq\log
m+\frac{1}{m}\sum_{j=1}^{m}\log
b_{j}+\frac{1}{m}\sum_{j=1}^{m}\beta^{1}_{j}E\left[\frac{1}{n}\sum_{i=1}^{n}\log(r_{1i})\right]$
$\displaystyle+\frac{1}{m}\sum_{j=1}^{m}\beta^{2}_{j}E\left[\frac{1}{n}\sum_{i=1}^{n}\log(r_{2i})\right]$

$GR_{avg}\geq\log m+\frac{1}{m}\sum_{j=1}^{m}\log
b_{j}+\frac{1}{m}\sum_{j=1}^{m}\beta^{1}_{j}G_{1}+\frac{1}{m}\sum_{j=1}^{m}\beta^{2}_{j}G_{2}$

$GR_{avg}\geq\beta_{1}^{AM}G_{1}+\beta_{2}^{AM}G_{2}-\log\frac{1}{m\,b^{GM}}$

where AM stands for arithmetic mean and GM stands for geometric mean and
$G_{1}$ and $G_{2}$ are buy and hold returns from factor 1 and factor 2
respectively. The main insight from the result being that the growth rate of a
portfolio that rebalances regularly can do better than the underlying factor
growth rates, depending on the average $\beta$ exposure of the assets
underlying the portfolio and how diversified the portfolio is. The last term
denotes a portfolio-drag term that vanishes for a $1/m$ portfolio.

###  2.3 Universality of Dirichlet factor portfolios

Next we present a modification of Cover’s[5] universal portfolio selection
process to incorporate economic information. We focus on the alpha parameter
of the dirichlet sampling on the unit simplex. Recall that alpha is a
concentration parameter of the dirichlet distribution, in the sense that a
diffuse alpha leads to portfolios with no prior knowledge. On the other
extreme, an alpha where one component is very large and the others very small
would tend to bias the portfolio towards that particular asset. We propose
using dirichlet($\alpha=\beta$) distributions, where $\beta$ are factor
loadings of assets, to create universal portfolios that perform better than
the benchmark universal portfolio that [5] suggests. To show that the
portfolio is universal, recall that Cover and Ordentlich[7] in result (24) in
their paper show that the total wealth generated by the dirichlet
$(1,1,1...1)$ distribution is given as:

$W_{\mathbf{1}}=E_{F(x),b}\left[\int_{b\in
D_{\mathbf{1}}}S_{n}(b,X^{n})\;d\mu(b)\right]$

where

$S_{n}(b,X^{n})=\prod_{i=1}^{n}b^{t}X_{i}$

The technique that we propose involves the same process, and samples
repeatedly from a dirichlet distribution ($\alpha=\beta$) that is driven by
factor exposures that each asset has to the underlying factors:

$W_{\mathbf{\beta}}=E_{F(x),b}\left[\int_{b\in
D_{\mathbf{\beta}}}S_{n}(b,X^{n})\;d\mu(b)\right]$

Keeping in mind that Cover and Ordentlich[7] show that the dirichlet
$(1,1...,1)$ process leads to a portfolio that is universal, we are interested
in analyzing the properties of the ratio of these wealth processes, a ratio
where the denominator is the naive dirichlet wealth with no information and
the numerator is the new informed dirichlet wealth with knowledge of the
factor $\beta$s:

$R_{n}=\frac{E_{F(x),b}\left[\int_{b\in
D_{\mathbf{\beta}}}S_{n}(b,X^{n})\;d\mu(b)\right]}{E_{F(x),b}\left[\int_{b\in
D_{\mathbf{1}}}S_{n}(b,X^{n})\;d\mu(b)\right]}$

$R_{n}=\frac{\int_{b\in
D_{\mathbf{\beta}}}E_{F(x),b}\left[S_{n}(b,X^{n})\right]\;d\mu(b)}{\int_{b\in
D_{\mathbf{1}}}E_{F(x),b}\left[S_{n}(b,X^{n})\right]\;d\mu(b)}$

$R_{n}=\frac{\int_{b\in
D_{\mathbf{\beta}}}E_{F(x),b}\left[\prod_{i=1}^{n}b^{t}X_{i}\right]\;d\mu(b)}{\int_{b\in
D_{\mathbf{1}}}E_{F(x),b}\left[\prod_{i=1}^{n}b^{t}X_{i}\right]\;d\mu(b)}$

with the goal of explicitly showing that this ratio converges to a number
larger then one, which would show that the wealth gains from the
Dirchlet($\beta_{1},\beta_{2},\beta_{3}...\beta_{m}$) in expectation dominates
the gains from the dirichlet(1,1,1…1) portfolio. Recall that the one-factor
model led to a relationship between the asset return and market as follows:

$\begin{bmatrix}X_{1}\\\ X_{2}\\\ .\\\ .\\\ X_{m}\\\
\end{bmatrix}=\begin{bmatrix}\beta_{1}^{t}\;r_{m}\\\ \beta_{2}^{t}\;r_{m}\\\
.\\\ .\\\ \beta_{m}^{t}\;r_{m}\\\ \end{bmatrix}$

This leads to the ratio being evaluated as:

$R_{n}=\frac{\int_{b\in
D_{\mathbf{\beta}}}E_{F(x),b}\left[\prod_{i=1}^{n}R_{m}^{i}\prod_{i=1}^{n}b_{i}^{t}\beta\right]\;d\mu(b)}{\int_{b\in
D_{\mathbf{1}}}E_{F(x),b}\left[\prod_{i=1}^{n}R_{m}^{i}\prod_{i=1}^{n}b_{i}^{t}\beta\right]\;d\mu(b)}$

$R_{n}=\frac{\int_{b\in
D_{\mathbf{\beta}}}E_{F(x),b}\left[\prod_{i=1}^{n}R_{m}^{i}\right]\prod_{i=1}^{n}E_{F(x),b}\left[b_{i}^{t}\right]\beta\;d\mu(b)}{\int_{b\in
D_{\mathbf{1}}}E_{F(x),b}\left[\prod_{i=1}^{n}R_{m}^{i}\right]\prod_{i=1}^{n}E_{F(x),b}\left[b_{i}^{t}\right]\beta\;d\mu(b)}$

$R_{n}=\frac{\int_{\Delta}E_{F(x),b^{\beta}}\left[\prod_{i=1}^{n}R_{m}^{i}\right]\prod_{i=1}^{n}E_{F(x),b^{\beta}}\left[b_{i}^{\beta}\right]^{t}\beta\;d\mu(b)}{\int_{\Delta}E_{F(x),b^{1}}\left[\prod_{i=1}^{n}R_{m}^{i}\right]\prod_{i=1}^{n}E_{F(x),b^{1}}\left[b_{i}^{1}\right]^{t}\beta\;d\mu(b)}$

$R_{n}=\frac{\int_{\Delta}\left[m^{n}(\sum_{j=1}^{m}\beta_{j}^{2})^{n}\right]\;d\mu(b)}{\int_{\Delta}\left[(\sum_{j=1}^{m}\beta_{j})^{2n}\right]\;d\mu(b)}$

which is a ratio of two integrals over the unit (non-negative) simplex where
the numerator always dominates the denominator:

$F_{n}=\frac{m^{n}(\sum_{j=1}^{m}\beta_{j}^{2})^{n}}{(\sum_{j=1}^{m}\beta_{j})^{2n}}\geq
1$ (2)

By Cauchy-Schwartz $F_{n}$ is greater than 1 (see Lemma in Appendix).
Consequently, $R_{n}$ is lower bounded by 1 and we have shown that the new
factor tilted portfolio proposed is also universal.

##  3 Empirical Results: Dirichlet factor portfolios

We highlight the practical consequences of our results by showing four use
cases across two separate equity market settings, the Nasdaq 100 and the FTSE
100. For both asset universes, we show that it is possible to construct factor
dirichlet portfolios for factors that are already known in the literature to
have an impact on the cross-section of returns for that asset class.

One such factor is size, where there is abundant literature to show that small
companies outperform their larger counterparts[9]. As a simple proxy for size,
we use the traded price at any point in time as an indication of size. While
this can be critiqued, it serves as a simple calculation keeping our
information set requirement within the pricing domain. Note that to get to a
set of weights that be used for the dirichlet sampling, a few cross-sectional
normalizations are also needed. As a starting example, we use the
transformations shown in Algorithm 1 to obtain an alpha value for our
dirichlet simulation at every step.

Algorithm 1 Calculating portfolio weights for SizeFactor based dirichlet
portfolio

while $t\neq T$ do

$\alpha_{jt}\leftarrow 1/p_{jt}$

$\alpha_{jt}\leftarrow\frac{\alpha_{jt}}{min(\alpha_{jt})}$

Simulate N-Managers dirichlet portfolios from Dirichlet($\alpha=\alpha_{t}$)

Calculate Universal Portfolios from factor biased managers and update the
wealth process

$w_{jt}\leftarrow\frac{\int(Sbdb)}{\int(Sdb)}$

$t\leftarrow t+1$

end while

The result of this algorithm is a portfolio that shows favorable
characteristics. We show the results for a few such factor portfolios in
Figure 1 for Nasdaq 100 constituents from August 2007 to June 2023. Results
are repeated for FTSE 100 constituents in Figure 2.

###  3.1 Nasdaq 100 Constituents

In this example we used a selection of 71 assets (listed cash equity) from the
US stock market that are a subset of the Nasdaq 100 index constituents.
Tickers that were part of the index as of June 6th, 2023 and had availability
of daily price time-series data for at least 4,000 (a sufficiently large
number) trading days were chosen as the asset universe. The tickers were
retrieved from the Nasdaq web page [3] on the same date, 101 tickers were
filtered using the above criteria to yield 71 tickers with at least the
minimum threshold of price history. Daily price data was downloaded from
EODHistoricalData’s (a commercial service) end of day API [1] . We down-
sampled the daily price data to a weekly time-series and used the weekly
median price as the entry and exit price for the trading simulation.

Table 1: Factor based Dirichlet Alpha | FactorType | Dirichlet Alpha
---|---|---
1 | Size | $\left[\frac{\frac{1}{p}}{min(\frac{1}{p})})\right]$
2 | Mom | $truncate\left[exp(winsorize(Z(r)))\right]$
3 | Sharpe | $truncate\left[winsorize(Z(S(r)))\right]$
4 | Compound | $\sqrt{MomSharpe}$

Table 1 shows the computational functional form behind each of the factors
shown in the charts. As we have access only to price data for this study, we
limit the experiment to only price-based factors. As an example, the momentum
factor is calculated as a function of the cumulative return of the asset (r).
In this case, we cross-sectionally z-score the cumulative return (r) to date,
winsorize the z-score between -6 and 6, exponentiate it and then truncate the
result to be always greater than 1. These back of the envelope functions are
meant to be easy to calculate and at the same time intend to mimic the factor
target. They also need to adhere to the allowed values for a dirichlet
distribution $\alpha$ parameter input.

$S$ in the sharpe factor stands for the rolling sharpe ratio, which is
computed using a 10 week span exponentially weighted moving average of the
returns divided by the 10 week span rolling standard deviation of the returns.
Rest of the function symbols are same as that in the Momentum factor. In each
case the most recent available metric is used. We also illustrate a compound
factor that is the geometric mean of the momentum and sharpe factor, that
opens up the possibility of constructing functional forms that depend on a
variety of factors.

The results are shown in Figure 1. Note that this back-test has an inherent
look-ahead bias in that it uses assets that survive the period between August
2007 and June 2023. Hence the asset universe comprises of higher quality
assets than the Nasdaq at that point in time. Since we are comparing only the
relative performance of the class of naive-dirichlet portfolios to the factor
biased dirichlet portfolios, this chart serves its purpose as a relative
comparison and is not indicative of absolute performance. The back-test shows
out-performance both for the momentum and size dirichlet portfolios in terms
of the growth rate of the portfolio.

Figure 2 shows the same effects for the constituents of the FTSE, a Europe
based index.

Figure 1: Backtest Results: Nasdaq 100 constituents with various dirichlet
factors

###  3.2 FTSE 100 Constituents

Now we use a selection of 87 cash equity assets from the stocks that are
listed on the London Stock Exchange (LSE) and are part of the FTSE 100 as of
June 22nd 2023. The criteria for selection were the same as in the previous
example with Nasdaq 100 stocks. Tickers that were part of the FTSE index as of
June 22nd, 2023 and had availability of daily price time-series data for 4,000
trading days as of July 6th 2023 were chosen as the asset universe. In the
FTSE dataset, the size effects are quite large. To control for this effect and
not to let it influence our results too much, we remove assets that had a
minimum price of less than $1 through the time period considered. Only one
such asset was removed from the dataset in this case. The tickers were
retrieved from the LSE web page [2]. Trading conditions and factors are the
same as what is described in the section on Nasdaq constituents. Similar
outsize performance is observed for factor dirichlet portfolios, with the size
based dirichlet portfolio dominating other allocations.

Figure 2: Backtest Results: FTSE 100 constituents with various dirichlet
factors

##  4 Conclusion

In this work, we’ve proposed a generative process to create robust, high
sharpe, diversified, universal portfolios that potentially outperform the
naive uninformed dirichlet universal portfolio allocation. An analytical proof
has been shown that the portfolios are indeed universal and they dominate
portfolios that arise out of a uniform dirichlet sampling. Particular examples
of this technique have been demonstrated on real market data from the U.S. and
European equities market. We believe that this technique can be used by small
and large investment managers, robo-advisors and algorithmic trading outfits
as a target portfolio for incremental alpha generation. Further work could
involve online learning of the parameter for dirichlet inputs based on recent
performance of factors in the market, and being able to blend them as shown in
the composite portfolio. Another line of work could include blending price and
non-price factors to create mixture portfolio processes.

## References

  * [1] Eod historical data.  https://eodhistoricaldata.com.  Accessed: 2023-07-24.
  * [2] Ftse100.  https://www.londonstockexchange .com/indices/ftse- 100/constituents/table.  Accessed: 2023-07-24.
  * [3] Nasdaq100.  https://www.nasdaq.com/market-activity/quotes/nasdaq- ndx-index.  Accessed: 2023-07-24.
  * Borodin et al. [2003] Allan Borodin, Ran El-Yaniv, and Vincent Gogan.  Can we learn to beat the best stock.  _Advances in Neural Information Processing Systems_ , 16, 2003.
  * Cover [1991] Thomas M Cover.  Universal portfolios.  _Mathematical finance_ , 1(1):1–29, 1991.
  * Cover [1999] Thomas M Cover.  _Elements of information theory_.  John Wiley & Sons, 1999.
  * Cover and Ordentlich [1996] Thomas M Cover and Erik Ordentlich.  Universal portfolios with side information.  _IEEE Transactions on Information Theory_ , 42(2):348–363, 1996.
  * Das and Banerjee [2011] Puja Das and Arindam Banerjee.  Meta optimization and its application to portfolio selection.  In _Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 1163–1171, 2011.
  * Fama and French [1993] Eugene F Fama and Kenneth R French.  Common risk factors in the returns on stocks and bonds.  _Journal of financial economics_ , 33(1):3–56, 1993.
  * Gyorfi et al. [2007] Laszlo Gyorfi, Andras Urban, and Istvan Vajda.  Kernel-based semi-log-optimal empirical portfolio selection strategies.  _International Journal of Theoretical and Applied Finance_ , 10(03):505–516, 2007.
  * Helmbold et al. [1998] David P Helmbold, Robert E Schapire, Yoram Singer, and Manfred K Warmuth.  On-line portfolio selection using multiplicative updates.  _Mathematical Finance_ , 8(4):325–347, 1998.
  * Jegadeesh and Titman [1993] Narasimhan Jegadeesh and Sheridan Titman.  Returns to buying winners and selling losers: Implications for stock market efficiency.  _The Journal of finance_ , 48(1):65–91, 1993\.
  * Li and Hoi [2014] Bin Li and Steven CH Hoi.  Online portfolio selection: A survey.  _ACM Computing Surveys (CSUR)_ , 46(3):1–36, 2014\.
  * Singer [1997] Yoram Singer.  Switching portfolios.  _International Journal of Neural Systems_ , 8(04):445–455, 1997.

## Appendix

Lemma to show (2) using Cauchy-Schwartz Inequality:

$\begin{Vmatrix}u.v\end{Vmatrix}^{2}\leq\begin{Vmatrix}u.u\end{Vmatrix}\begin{Vmatrix}v.v\end{Vmatrix}$

Setting $u$ as $\mathbf{1}$ and $v$ as $\beta$, both in $\mathbf{R^{m}}$

$\begin{Vmatrix}\mathbf{1}.\beta\end{Vmatrix}^{2}\leq\begin{Vmatrix}\mathbf{1}.\mathbf{1}\end{Vmatrix}\begin{Vmatrix}\beta.\beta\end{Vmatrix}$

$\begin{Vmatrix}\mathbf{1}.\beta\end{Vmatrix}^{2}\leq
m\sum_{i=1}^{m}\beta^{2}$

$(\sum_{i=1}^{m}\beta_{j})^{2}\leq m\sum_{i=1}^{m}\beta_{j}^{2}$

∎