# Multiple Steady and Oscillatory Solutions in a Collapsible Channel Flow

Danyang<EMAIL_ADDRESS>Xiaoyu Luo  and Peter S. Stewart
School of Mathematics and Statistics, University of Glasgow, Glasgow, G12 8SQ
UK

(date; date)

###### Abstract

We study flow driven through a finite-length planar rigid channel by a fixed
upstream flux, where a segment of one wall is replaced by a pre-stressed
elastic beam subject to uniform external pressure. The steady and unsteady
systems are solved using a finite element method. Previous studies have shown
that the system can exhibit three steady states for some parameters (termed
the upper, intermediate and lower steady branches, respectively). Of these,
the intermediate branch is always unstable while the upper and lower steady
branches can (independently) become unstable to self-excited oscillations. We
show that for some parameter combinations the system is unstable to both upper
and lower branch oscillations simultaneously. However, we show that these two
instabilities eventually merge together for large enough Reynolds numbers,
exhibiting a nonlinear limit cycle which retains characteristics of both the
upper and lower branches of oscillations. Furthermore, we show that increasing
the beam pre-tension suppresses the region of multiple steady states but
preserves the onset of oscillations. Conversely, increasing the beam thickness
(a proxy for increasing bending stiffness) suppresses both multiple steady
states and the onset of oscillations.

###### keywords:

flow-vessel interaction; collapsible channel flows; stability.

##  1 Introduction

Fluid flow through flexible conduits in the human body can exhibit a wide
variety of interesting physiological phenomena [Heil and Hazel, 2011]; such
flows can be investigated experimentally using a Starling Resistor, where
fluid is driven through a segment of externally pressurised flexible tubing by
either a prescribed driving pressure or a prescribed upstream flux [e.g.
Bertram, 1982, 1986]. In particular, pressure driven flow can exhibit a
phenomenon known as ‘flow limitation’, where the flow rate along the tube does
not continue to increase as the driving increases: the large flow speeds
collapse the tube through the Bernoulli effect inhibiting the flow [Bertram
and Castles, 1999; Bertram and Tscherry, 2006]. Similarly, flux driven systems
can exhibit the associated phenomenon of ‘pressure drop limitation’, where the
pressure difference between the ends of the collapsible tube does not increase
with increases in flow rate [Bertram, 1986; Bertram _et al._ , 1990]. In
physiology, flow limitation can occur during tidal breathing, when the flow
rate of air being expelled from the lung becomes maximised during forced
expiration [Tantucci, 2013].

In some cases these highly collapsed steady configurations can co-exist with
other, more inflated, configurations of the vessel. Such multiplicity in
steady solutions is evident in Starling Resistor experiments, where the system
can exhibit hysteresis i.e. different steady configurations at the same
operating point dependent on the history and termed an ‘open-to-closed
transition’ [Bertram _et al._ , 1991; Bertram and Castles, 1999]. Such co-
existing steady states are also evident in theoretical models of flow in
collapsible tubes [Hazel and Heil, 2003; Heil and Boyle, 2010]. In particular,
lumped parameter models indicate that the collapsible tube can typically
exhibit three co-existing branches of steady solutions across a range of
parameters [Armitstead _et al._ , 1996]; these three branches have been termed
the upper branch (where the tube is mostly inflated), the lower branch (where
the tube is highly collapsed) with an intermediate branch between them. The
upper (lower) steady solutions is connected to the intermediate branch through
a saddle-node bifurcation at the upper (lower) branch limit point; the upper
and lower steady branches are stable to non-oscillatory perturbations and the
intermediate branch is always unstable [Armitstead _et al._ , 1996].

In physiology, flow limitation in the lung airways is often accompanied by
wheezing, associated with rapid fluttering of the airway wall [Grotberg and
Gavriely, 1989]. Similarly, Starling Resistor experiments investigating flow
limitation can sometimes exhibit spontaneous self-excited oscillations which
fall into distinct frequency bands [Bertram and Pedley, 1982; Bertram _et al._
, 1990, 1991]. Such self-excited oscillations are also evident in theoretical
models of flow in collapsible tubes, such as lumped parameter models [Bertram
and Pedley, 1982; Armitstead _et al._ , 1996], cross-sectionally averaged one-
dimensional models [Jensen, 1990] and full three-dimensional models [Heil and
Boyle, 2010; Whittaker _et al._ , 2010].

In this study we consider a theoretical model for a planar analog of the
Starling Resistor experiment, formed by removing a segment of one wall of a
rigid channel and inserting an elastic wall. This planar analog exhibits
multiple steady solutions in some parameter regimes [Luo and Pedley, 2000;
Heil, 2004; Stewart, 2017; Wang _et al._ , 2021; Herrada _et al._ , 2021],
with a three branch structure qualitatively similar to the collapsible tube
models. In addition, this collapsible channel system exhibits transition to
self-excited oscillations from both the lower (collapsed) branch of steady
solutions [Heil, 2004; Stewart, 2017] as well as the upper (inflated) branch
of steady solutions [Herrada _et al._ , 2021; Wang _et al._ , 2021]. Fully
developed upper branch oscillations exhibit an upstream propagating hump along
the compliant segment; this wave is suppressed by reflection at the upper
rigid segment and replaced by a new wave at the downstream end of the
compliant segment [Wang _et al._ , 2021]. By contrast, fully developed lower
branch oscillations exhibit a constriction at the downstream end of the
compliant segment which propagates up and downstream over a narrow range [Luo
and Pedley, 1996; Wang _et al._ , 2021].

In particular, we use a model where the wall takes the form of a thin pre-
stressed (massless) elastic beam with resistance to both bending and
stretching (Sec. 2) [Cai and Luo, 2003; Luo _et al._ , 2008; Wang _et al._ ,
2021]. We use this model to explore how the size of the region of multiple
steady states is influenced by the system parameters, particularly the beam
pre-tension and bending stiffness (Sec. 3). Finally, we also use this model to
explore the onset of self-excited oscillations from both the upper and lower
branches of steady solutions, showing that in some cases these two families of
oscillations merge together (Sec. 4).

##  2 The Model

Figure 1: The fluid-beam model.

We consider an incompressible Newtonian fluid (of density $\rho$ and viscosity
$\mu$) driven through a finite-length rigid channel with uniform width $D$. We
consider a parabolic inlet flow driven through the channel with flux $Q$
against a fixed outlet pressure $p_{0}$. One segment of the upper wall is
replaced by a plane strained elastic beam of thickness $h$ subject to a
uniform external pressure $p_{e}$. The corresponding lengths of the upstream
and downstream rigid and elastic segments are denoted $L_{u}$, $L_{d}$ and
$L$, respectively. We parameterize the domain using a Cartesian coordinate
system with origin at the intersection between the upstream rigid segment and
the compliant segment on the entirely rigid wall (see Fig. 1). Time is denoted
$t$.

We model the flexible wall as a massless elastic beam and denote the axial
pre-tension along the beam as $T$. The extensional and bending stiffnesses of
the beam are denoted as $EA$ and $EJ$, respectively. Here $E$ is the Young’s
modulus of the material while $A$ and $J$ are the cross-sectional area and the
second moment of inertia of the cross-sectional area of the beam,
respectively. The undeformed elastic beam is parameterized by the coordinate
$l$ (where $0\leqslant l\leqslant L$).

###  2.1 Governing Equations

To form non-dimensional variables we scale all lengths on $D$, velocities on
$Q/D$, time on $Q/D^{2}$ and pressures on the inertial scale according to

$\displaystyle p=\frac{\rho Q^{2}}{D^{2}}\tilde{p}+p_{0},$ (1)

where $p$ ($\tilde{p}$) denotes the dimensional (dimensionless) fluid
pressure.

The corresponding dimensionless parameters take the form

$\displaystyle\tilde{c}_{\lambda}$ $\displaystyle=\frac{(EA)D}{\rho Q^{2}},$
$\displaystyle\tilde{c}_{\kappa}$ $\displaystyle=\frac{EJ}{\rho Q^{2}D},$
$\displaystyle\tilde{T}$ $\displaystyle=\frac{TD}{\rho Q^{2}},$
$\displaystyle\tilde{R}e$ $\displaystyle=\frac{Q\rho}{\mu},$
$\displaystyle\tilde{p}_{e}$ $\displaystyle=\frac{(p_{e}-p_{0})D^{2}}{\rho
Q^{2}},$ (2)

where $\tilde{c}_{\lambda}$, $\tilde{c}_{\kappa}$ and $\tilde{T}$ are the
dimensionless extensional, bending stiffnesses and pre-tension of the elastic
beam, respectively, $\tilde{R}e$ is the Reynolds number and $\tilde{p}_{e}$ is
the external pressure. The corresponding dimensionless lengths of the
upstream, downstream and collapsible segments of the channel and the
dimensionless beam thickness are scaled as

$\displaystyle\tilde{L}_{u}$ $\displaystyle=\frac{L_{u}}{D},$
$\displaystyle\tilde{L}_{d}$ $\displaystyle=\frac{L_{d}}{D},$
$\displaystyle\tilde{L}$ $\displaystyle=\frac{L}{D},$ $\displaystyle\tilde{h}$
$\displaystyle=\frac{h}{D},$ (3)

respectively. We henceforth focus on the dimensionless quantities and drop the
tildes for simplicity.

####  2.1.1 Fluid Equations

The governing equations for the (two-dimensional) fluid are the incompressible
Navier-Stokes equations in the form,

$\nabla\cdot{\mathbf{u}}=0,\quad\frac{\partial{\mathbf{u}}}{\partial
t}+\left({\mathbf{u}}\cdot\nabla\right){\mathbf{u}}=\nabla\cdot{\boldsymbol{\sigma}},$
(4)

where ${\mathbf{u}}=(u,v)$ is the planar fluid velocity and the Newtonian
stress tensor ${\boldsymbol{\sigma}}$ takes the form

${\boldsymbol{\sigma}}=-p{\mathbf{I}}+Re^{-1}\left(\nabla{\mathbf{u}}+{\nabla{\mathbf{u}}}^{\rm{T}}\right),$
(5)

where $\mathbf{I}$ is the identity matrix and the superscript T represents the
matrix transpose.

####  2.1.2 Beam Equations

We denote the two components of beam deformation as
${\mathbf{x}_{b}}=(x_{b}(l,~{}t),y_{b}(l,~{}t))$ in terms of the beam
coordinate $l$; using a modified constitutive law for the massless beam, the
dimensionless governing equations for the beam can be written as (details see
Wang _et al._ [2021])

$\displaystyle
c_{\kappa}\kappa\frac{\partial\left(\lambda\kappa\right)}{\partial
l}+c_{\lambda}\frac{\partial\lambda}{\partial l}+\lambda{\sigma_{1}}=0,$ (6)
$\displaystyle-c_{\kappa}\frac{\partial}{\partial
l}\left(\frac{1}{\lambda}\frac{\partial(\lambda\kappa)}{\partial
l}\right)+c_{\lambda}\lambda\kappa(\lambda-1)+\lambda\kappa
T+\lambda{\sigma_{2}}-\lambda p_{e}=0,$ (7) $\displaystyle\frac{\partial
x_{b}}{\partial l}=\lambda\cos\theta,\quad\frac{\partial y_{b}}{\partial
l}=\lambda\sin\theta,\quad\frac{\partial\theta}{\partial l}=\lambda\kappa,$
(8)

where $\theta$ is the angle between the rigid wall and the tangent vector of
the deformed beam (see Fig. 1), $\kappa$ is the curvature of the beam and
$\lambda$ is the principal stretch of the beam, which can each be computed in
terms of the beam deformation as

$\displaystyle\kappa$ $\displaystyle=\frac{1}{\lambda^{3}}\left(\frac{\partial
x_{b}}{\partial l}\frac{\partial^{2}y_{b}}{\partial l^{2}}-\frac{\partial
y_{b}}{\partial l}\frac{\partial^{2}x_{b}}{\partial l^{2}}\right),$
$\displaystyle\lambda$ $\displaystyle=\sqrt{\left(\frac{\partial
x_{b}}{\partial l}\right)^{2}+\left(\frac{\partial y_{b}}{\partial
l}\right)^{2}},$ $\displaystyle(0\leq l\leq L).$ (9)

In addition, ${\sigma_{1}}$ and ${\sigma_{2}}$ are the tangent and normal
components of the fluid traction on the beam,

$\sigma_{1}=\left(-{\boldsymbol{\sigma}}\hat{{\mathbf{n}}}\right)\cdot\hat{{\mathbf{t}}},\quad\sigma_{2}=\left(-{\boldsymbol{\sigma}}\hat{{\mathbf{n}}}\right)\cdot\hat{{\mathbf{n}}},$
(10)

where $\hat{{\mathbf{t}}}$ and $\hat{{\mathbf{n}}}$ are the tangent and normal
unit vectors of the deformed beam (see Fig. 1).

####  2.1.3 Boundary Conditions

We prescribe a parabolic inlet flow with unit flux in the form

$\displaystyle u$ $\displaystyle=6y(1-y),$ $\displaystyle v$
$\displaystyle=0,$ $\displaystyle(x$ $\displaystyle=-L_{u},0\leqslant
y\leqslant 1).$ (11)

We assume the no-slip condition along the rigid walls as well as continuity of
velocity between the elastic beam and the fluid in the form

$\displaystyle{\mathbf{u}}$ $\displaystyle={\mathbf{0}},~{}~{}\quad$
$\displaystyle(y=0;~{}~{}y=1,-L_{u}\leq x\leq 0,~{}L\leq x\leq L+L_{d}),$ (12)
$\displaystyle{\mathbf{u}}$
$\displaystyle={\mathbf{u}_{b}}=\frac{\partial{\mathbf{x}_{b}}}{\partial
t},\quad$ $\displaystyle({\mathbf{x}}\in\partial\Omega_{b}).$ (13)

The two ends of the elastic beam are fixed to the rigid wall in the form

$\displaystyle x_{b}(0,t)=0,\quad y_{b}(0,t)=1,\quad x_{b}(L,t)=L,\quad
y_{b}(L,t)=1,\quad\theta(0,t)=\theta(L,t)=0.$ (14)

###  2.2 The Finite Element Method

A finite-element method is used to solve the coupled fluid-beam system. We
divide the fluid domain into three sections, denoted as A, B and C for the
upstream, compliant and downstream compartments (Fig. 1), respectively [Luo
and Pedley, 1996; Luo _et al._ , 2008]. In section B, we use an adaptive mesh
that consists of rotational spines that connect fixed nodes in the bottom wall
with nodes in the deformable beam [Cai and Luo, 2003]. Then nodes are seeded
along these spines covering region B. Each spine can rotate around its fixed
node on the rigid wall. Hence, nodes in section B can move along the
rotational spine as the beam is deformed. Further details of the numerical
method are provide elsewhere [Luo _et al._ , 2008; Hao _et al._ , 2016]. A
mesh of 36657 elements is used for the numerical solutions in this study with
time-step $\Delta t=0.01$. Convergence tests of grid- and time-step-
independence were carried out between three meshes and two choices of time-
step (for details see Wang _et al._ [2021]).

###  2.3 Parameter Choices

In our previous studies we modeled the elasticity of the beam by varying the
extensional stiffness $c_{\lambda}$ [Luo _et al._ , 2008; Wang _et al._ ,
2021], which is proportional to the bending stiffness of the beam since
$c_{\kappa}=\left(h^{2}/12\right)c_{\lambda}$, and fixed the other beam
parameters including the pre-tension ($T=0$), beam thickness ($h=0.01$) and
external pressure ($p_{e}=1.95$). Conversely, in this study we fix
$c_{\lambda}=1600$ and examine the role of these other parameters on both the
steady (Sec. 3) and unsteady (Sec. 4) behaviour of the system. In particular
we note that increasing the beam thickness allows us to alter the bending
stiffness of the structure independently of the extensional stiffness.
Following Luo _et al._ [2008], we fix the geometry of the channel according to
$L_{u}=L=5$, $L_{d}=30$.

##  3 Multiple Steady Solutions

We discuss the predictions of the steady system considering the role of
increasing pre-tension (Sec. 3.1) and beam thickness (Sec. 3.2).

###  3.1 Role of Pre-tension

Figure 2: Steady solutions for fixed beam thickness $h=0.01$: (a) maximal
($y_{max}$, solid line) and minimal ($y_{min}$, dash-dot line) beam position
against $Re$ for $T=0$ (color black), $T=1$ (color red) and $T=5$ (color blue)
for $p_{e}=1.1$; (b), (c) The steady beam shapes for $p_{e}=1.1$ with $Re=340$
(solid line) and $Re=360$ (dashed line) for $T=0$ and $T=5$, respectively.
(d)-(f) The upper and lower limit points in parameter space spanned by the
external pressure $p_{e}$ and Reynolds number $Re$ for $T=0$, $T=0.1$ and
$T=1$, respectively. Note the vertical solid and dashed lines in panel (a) are
at $Re=340$ and $360$, respectively. The region with multiple steady solutions
is shaded in blue in panels (d)-(f).

In order to investigate the role of pre-tension $T$ in setting the steady beam
shape, Fig. 2 summarises the steady solutions of the maximal and minimal beam
positions (Fig. 2a), several typical beam shapes (Fig. 2b) and the upper and
lower limit points (Fig. 2c) with different values of pre-tension with fixed
wall thickness $h=0.01$. In particular, the maximal ($y_{max}$, solid line)
and minimal ($y_{min}$, dash-dot line) beam deflection on the $y-$direction as
a function of Reynolds number $Re$ is shown in Fig. 2(a) for $T=0$ (black),
$T=1$ (red) and $T=5$ (blue) at fixed external pressure $p_{e}=1.1$. For low
Reynolds numbers the beam is fully inflated (i.e. $y_{min}=1$); this is the
upper branch of steady solutions. As the Reynolds number increases the elastic
beam becomes increasingly collapsed though the Bernoulli effect. Notably, for
$T=0$ the system enters a region with three steady states for $354.56\lesssim
Re\lesssim 362.78$ similar to our previous model [Wang _et al._ , 2021]. As
the Reynolds number further increases, the system again exhibits a unique
steady state for $Re>362.78$; this branch is the lower branch of steady
solutions. The upper and lower branches of steady solutions are connected by
an intermediate branch, which they intersect at the upper and lower limit
points, respectively.

A similar region with multiple steady states is observed for $T=1$
($355.79\lesssim Re\lesssim 358.75$), although the beam pre-tension has
significantly narrowed the range of Reynolds numbers for which it is evident
for this choice of external pressure. The region with multiple steady states
vanishes entirely for large pre-tension ($T=5$, blue lines in Fig. 2a).

To illustrate these steady configurations in detail, Fig. 2(b,c) illustrates
the possible steady wall shapes at both $Re=340$ and $Re=360$ for $T=0$ and
$T=5$, respectively. At $Re=340$, located on the upper branch for $T=0$ (see
the solid vertical line in Fig. 2a), the steady wall shape is fully inflated
and is unique (solid line in Fig. 2b). Conversely, at $Re=360$, located in the
region with multiple steady states for $T=0$ (see the dashed vertical line in
Fig. 2a), there are three possible steady wall shapes (Fig. 2b). Whereas for
$T=5$ (Fig. 2c), the system has a unique steady state for both Reynolds
numbers and the wall is significantly less deflected due to the larger pre-
tension.

Fig. 2(d-f) summarises the region with multiple steady states for $T=0$ (Fig.
2d), $T=0.1$ (Fig. 2e) and $T=1$ (Fig. 2f) in parameter space spanned by the
external pressure and Reynolds number. In each case we see a triangular region
with multiple steady states with external pressure below a critical value and
Reynolds number above a critical value. As $T$ increases, this critical point
is gradually pushed to larger Reynolds number and lower external pressures.
However, the region with multiple steady states (Fig. 2f) remains of
approximately the same width relative to the critical point. In summary, large
pre-tension suppresses multiple steady states of the system by shifting the
critical point for multiple solutions across the parameter space.

###  3.2 Role of Beam Thickness

Figure 3: Steady solutions for $T=0$, showing (a) the maximal ($y_{max}$,
solid line) and minimal ($y_{min}$, dash-dot line) beam position against $Re$
for $h=0.01$ (color black), $h=0.086$ (color red) and $h=0.122$ (color blue);
the upper and lower limit points in parameter space spanned by the external
pressure $p_{e}$ and Reynolds number $Re$ for $h=0.01$ (b), $h=0.027$ (c) and
$h=0.086$ (d), respectively. In panels (b)-(d), the region with multiple
steady solutions is shaded in blue.

To investigate the role of increased bending stiffness on the steady
configuration of the beam, Fig. 3 demonstrates the steady solutions with
various beam thickness, plotting the maximal ($y_{max}$) and minimal
($y_{min}$) beam deflection against the Reynolds number $Re$ for three beam
thicknesses ($h=0.01,~{}0.086$, and $0.122$) with fixed pre-tension $T=0$ and
external pressure $p_{e}=1.1$ in Fig. 3(a). Similar to Fig. 2, the system
shows multiple steady states for $h=0.01$ ($354.56\lesssim Re\lesssim
362.78$); this region of multiple solutions narrows with increased beam
thickness $h=0.086$ ($354.87\lesssim Re\lesssim 355.8$). Additionally, the
steady system is unique for sufficiently large beam thickness $h=0.122$. Fig.
3(b-d) summarises the corresponding upper and lower limit points in parameter
space spanned by the external pressure $p_{e}$ and Reynolds number $Re$ for
the same three beam thicknesses, showing that the region with multiple steady
states narrows as $h$ increases, but the critical point for the existence of
multiple solutions does not vary much across the parameters tested. In
summary, this figure demonstrates that increasing the beam thickness (and
hence the bending stiffness) eventually suppresses multiple steady states;
however, this is achieved by narrowing the tongue while holding the critical
point fixed.

##  4 Self-excited Oscillations

In order to test the stability of the system for a given parameter combination
we apply a small increment to the steady solution (here we use the steady
solution corresponding to a $1\%$ increase in the extensional stiffness
$c_{\lambda}$). The system is deemed stable if the unsteady solution converges
to its corresponding steady solution, and unstable if the perturbation grows
[Drazin, 2002]. The state between stable and unstable is termed as neutrally
stable. In this model the perturbation grows in an oscillatory manner and
saturates into a finite amplitude limit cycle. In this study, we report the
dynamics of this oscillatory limit cycle and ignore the initial transient. In
particular, we present phase portraits of the oscillation as a function of the
wall pressure at the upstream and downstream ends of the compliant segment
(e.g. Fig. 5c-f, i-l and Fig. 6b-e) and compute the fluid pressure on the wall
at the channel midpoint time-averaged over a period of oscillation (e.g. Fig.
5a,g, Fig. 6a, Fig. 8 and Fig. 9). In this section, we focus on the unsteady
solutions for fixed extensional stiffness ($c_{\lambda}=1600$), testing the
stability of both the upper and lower steady branches. In particular, we
investigate oscillatory solutions with no pre-tension and low beam thickness
(Sec. 4.1), large beam thickness (Sec. 4.2) and large pre-tension (Sec. 4.3).

###  4.1 Multiple Oscillatory Solutions with Low Bending Stiffness

Figure 4: An overview of the parameter space spanned by the Reynolds number
and external pressure at $T=0$ and $h=0.01$. The computed neutrally stable
points are denoted as filled black circles, and the upper and lower limit
points are denoted as filled black squares connected with dashed and solid
lines. The region where the system is stable is shaded in grey and the region
with multiple steady solutions is shaded in blue. (b), (c) show the zoom-in of
the region near $p_{e}=2.25$ and $p_{e}=1.48$, respectively. Note the unsteady
solutions of the operating points for $p_{e}=2.25$: U1 ($Re=168.85$), U2
($Re=171$), U3 ($Re=173$), U4 ($Re=174.5$); and for $p_{e}=1.48$: L1 (U1a)
($Re=266.8$), L2 ($Re=271$), L3 ($Re=274$) are shown in Fig. 5.

Figure 5: Unsteady solutions of the upper and lower steady branches for
$p_{e}=2.25$ and $p_{e}=1.48$. (a, g) Bifurcation diagram of midpoint wall
pressure as a function of the Reynolds number for $p_{e}=2.25$ and
$p_{e}=1.48$, respectively; (b, h) Oscillation period as a function of the
Reynolds number for $p_{e}=2.25$ and $p_{e}=1.48$, respectively; (c-f, i-l)
Phase portraits in the space spanned by the wall pressure measured at the
upstream and downstream ends of the compliant segment for $p_{e}=2.25$ and
$p_{e}=1.48$, respectively.

Figure 6: Unsteady bifurcation diagram of the system for $T=0$, $h=0.01$ and
$p_{e}=1.1$, showing (a) The midpoint wall pressure $p_{mid}$ as a function of
the Reynolds number $Re$. The steady midpoint pressure $p_{mid}$ (solid
(stable) and dashed (unstable) lines), the time-averaged midpoint wall
pressure $p_{mid}^{avg}$ (triangles), the upper and lower branch limit points
(squares); (b-e) Phase portraits in space spanned by the wall pressure
measured at the upstream and downstream ends of the compliant segment at
operating points i-iv labelled in (a). The corresponding values of the upper,
intermediate and lower steady branches are denoted by open circles.

Figure 7: Self-excited oscillations arising at operating point iii ($Re=358$,
$p_{e}=1.1$, $T=0$ and $h=0.01$ labelled in Fig. 6a), showing (a) time history
of the wall pressure measured at upstream ($p_{w}(0,t)$) and downstream
($p_{w}(L,t)$) ends of the compliant segment; (b-f) streamlines and pressure
contours of the flow at five equally-spaced time instances over a period of
oscillation, labelled in (a).

We first analyse unsteady simulations of perturbations to the upper and lower
steady branches with no pre-tension ($T=0$) and a thin beam ($h=0.01$). Here
we extend our previous analysis [Wang _et al._ , 2021] to examine the
parameter space spanned by the external pressure and Reynolds numbers, as
shown in Fig. 4 (Note that in Wang _et al._ [2021] we held the external
pressure $p_{e}=1.95$ throughout). The steady behaviour of the system is
similar to the cases presented in Fig. 2(d-f) and 3(b-d), where three steady
states can co-exist in a narrow tongue between the upper and lower limit
points (blue shaded region marked by filled symbols). Testing the stability of
these steady solutions we find that both the upper and lower steady branches
can independently become unstable to oscillations (each via a supercritical
Hopf bifurcation) in the neighbourhood of the region of parameter space with
multiple steady states. The computed neutrally stable points from the upper
and lower steady branches are marked as filled circles, and the neutral
stability curves are denoted as dot-dashed lines connecting these neutral
points. In this parameter space the upper branch is unstable within a narrow
tongue to the left of the trace of the upper branch limit points. This upper
branch neutral curve intersects the trace of the upper branch limit points at
a co-dimension 2 (Takens-Bogdanov) point [Glendinning, 1994; Strogatz, 2018]
at $p_{e}\approx 2.25$, $Re\approx 174.5$ (we term this the upper Takens-
Bogdanov point). For external pressures less than this upper Takens-Bogdanov
point the upper branch oscillation is eventually restabilised as the Reynolds
number increases via an interaction between the oscillatory limit cycle and
the upper branch limit point [see details in Wang _et al._ , 2021]. However,
this unstable tongue extends to slightly larger external pressures than those
which admit see multiple steady states. For external pressures greater than
the upper Takens-Bogdanov point the upper branch oscillation instead
restabilises via a second Hopf bifurcation as the Reynolds number increases.
This interaction is explored further in Fig. 5, where we examine the time-
averaged midpoint pressure of fully developed upper branch oscillations (Fig.
5a), their corresponding period (Fig. 5b) and illustrate several limit cycles
of oscillation (Fig. 5c-f). As in Wang _et al._ [2021], this upper branch
oscillation is associated with an overall decrease in the time-averaged
midpoint pressure compared to the steady state (Fig. 5a), and the period of
oscillation increases as the upper limit point is approached (Fig. 5b). The
oscillation develops a rather complicated limit cycle (Fig. 5c-f), similar to
the upper branch oscillations reported in Wang _et al._ [2021].

Furthermore, the lower branch of steady solutions is unstable within a tongue
to the right of the region with multiple steady states. This tongue is
stabilised for large external pressure ($p_{e}\gtrsim 2.12$), and tracks close
to the curve of lower branch limit points, merging at a second (lower) Takens-
Bogdanov point at $Re\approx 266.8$, $p_{e}\approx 1.48$. The neighbourhood of
this lower Takens-Bogdanov point is explored in Fig. 5(g-l), plotting the
time-averaged midpoint pressure (Fig. 5g), period of oscillation (Fig. 5h) and
several limit cycles (Fig. 5i-l) for fixed external pressure $p_{e}=1.48$.
During oscillation the midpoint pressure is again less than the corresponding
steady configuration. Just beyond the critical Reynolds number for lower
branch instability, the oscillation grows and saturates into an elliptical
limit cycle of larger period than the upper branch oscillations (open circles
in Fig. 5h) but of much smaller amplitude (open circles in Fig. 5g). A phase
portrait for point L1 is shown in Fig. 5(i). For Reynolds numbers sufficiently
close to critical ($Re\gtrsim 266.8$) this small-amplitude limit cycle is
maintained over the lifetime of our numerical simulations. However, for larger
Reynolds numbers ($Re\gtrsim 270.3$) the system visits this limit cycle only
transiently and eventually saturates into a limit cycle of larger amplitude
(open triangles in Fig. 5g) and shorter period (open triangles in Fig. 5h).
These nonlinear oscillations are a direct continuation of the oscillations
bifurcating from the upper branch steady state, displaying analogous phase
portraits (Fig. 5j,k). Furthermore, across the range of Reynolds numbers which
exhibit the small-amplitude limit cycle growing from the lower steady branch
($266.8\lesssim Re\lesssim 270$), it turns out that the system also exhibits a
saturated limit cycle analogous to the upper branch oscillations for these
points as well (amplitude and period shown as open triangles in Fig. 5(g,h)
and a typical phase portrait is shown in Fig. 5(l). Hence, this system appears
to exhibit hysteresis in the oscillatory behaviour across a small range of
parameters with two possible branches of fully developed oscillations (open
circles and open triangles in Fig. 5h); these two branches merge together at
$Re\approx 270.3$ and then continue to large Reynolds numbers (open triangles
in Fig. 5h). There are two possibilities to explain this observation: either
the system exhibits two co-existing limit cycles across a range of parameters,
or these reported lower branch (small-amplitude) limit cycles are in fact long
transients which eventually grow and saturate to a limit cycle analogous to
the upper branch oscillations (note that we see no evidence of this transition
for simulations close to the critical Reynolds number over the long time
intervals considered in our simulations). Either way, for large enough
Reynolds numbers the system exhibits a merged family of oscillations (open
triangles in Fig. 5h) which grow from the lower steady branch but are a direct
continuation of nonlinear oscillations bifurcating from the upper steady
branch.

It emerges that for external pressures less than the lower Takens-Bogdanov
point, this merged family of oscillations becomes more prominent with a
typical example for $p_{e}=1.1$ shown in Fig. 6, showing the bifurcation
diagram (Fig. 6a) and four typical limit cycles in the neighbourhood of the
region with multiple steady states (Fig. 6b-e). As the Reynolds number
increases the system becomes unstable along the upper branch of steady
solutions, exhibiting a decrease in the time-averaged midpoint pressure with
an elaborate limit cycle around the upper branch steady state (Fig. 6b). As
the Reynolds number increases into the region with multiple steady states the
oscillation encompass all three steady solutions (Fig. 6c, d), and there is
now no interaction with the intermediate steady state (a prominent feature of
purely upper branch oscillations, as discussed in Wang _et al._ [2021]). As
the Reynolds number continues to increase into the region with only a lower
steady state, the oscillation exhibits a qualitatively similar limit cycle
(Fig. 6e).

To further explore the nature of these merged oscillations, Fig. 7 examines
the oscillation in more detail at a point in parameter space which exhibits
multiple steady states, showing time-trace of the wall pressure at the two
ends of the compliant segment (Fig. 7a) and several snapshots of the flow-
field and pressure contours in the channel (Fig. 7b-e). The fully developed
(nonlinear) oscillation retains many of the characteristics of the upper and
lower branch oscillations from which it has merged (these individual
oscillations are discussed in Wang _et al._ [2021]). Similar to the upper
branch oscillations, the limit cycle involves an upstream propagating hump in
the wall profile which is suppressed by interaction with the upstream rigid
segment and replaced by a new upstream propagating hump originating at the
downstream end of the compliant segment. However, like the lower branch
oscillations, the wall profile exhibits significant indentation into the
channel (and so much lower fluid pressures). This greater wall indentation
leads to the formation of a low pressure region near the end of the compliant
segment (Fig. 7d), which creates a (weak) vorticity wave propagating along the
downstream rigid segment (Fig. 7e, f), similar to oscillations about a
collapsed steady state [Luo and Pedley, 1996]. Interestingly, the wall
pressures at the two ends of the compliant segment are much closer to being in
phase than the upper or lower branch oscillations alone. In summary, the
figure illustrates the flow profile of our new family of self-excited
oscillations, showing that it retains characteristics of both the upper and
lower oscillations that have been reported elsewhere [Wang _et al._ , 2021].

###  4.2 Oscillations with Large Beam Thickness

Figure 8: Bifurcation diagrams of midpoint wall pressure $p_{mid}$ as a
function of the Reynolds number $Re$ at $T=0$, $p_{e}=1.1$ for (a) $h=0.086$
and (b) $h=0.094$ and $h=0.122$. The steady midpoint pressure $p_{mid}$ is
denoted as solid (stable) and dashed (unstable) lines. The time-averaged
midpoint wall pressure $p_{mid}^{avg}$ is denoted as filled triangles, the
upper and lower branch limit points are denoted as squares.

To evaluate the role of the beam thickness, $h$, in the unsteady response of
the system, Fig. 8 presents bifurcation diagrams showing the onset of self-
excited oscillations, plotting the time-averaged midpoint pressure alongside
the corresponding steady midpoint pressure as a function of the Reynolds
number for three beam thicknesses ($h=0.086$, 0.094 and 0.012) with zero pre-
tension and fixed external pressure ($p_{e}=1.1$). For beam thickness
$h=0.086$ (Fig. 8a), the system exhibits multiple steady states in the range
$354.87\lesssim Re\lesssim 355.8$. In this case the oscillations from the
upper and lower steady branches merge into one family of oscillations across
the region with multiple steady states (similar to Fig. 6). When the beam
thickness is increased to $h=0.094$, the multiple steady states vanish but the
system still exhibits transition to self-excited oscillations at $Re\approx
344.07$ (denoted as the dashed line), which grow in amplitude as the Reynolds
number increases (Fig. 8b). However, as the beam thickness is further
increased to $h=0.122$ the system no longer becomes unstable to oscillations
and the steady system has a unique solution for all external pressures (Fig.
8b). In summary, increasing the thickness of the beam (and hence the bending
stiffness) suppresses the onset of self-excited oscillations as well as the
region with multiple steady states.

###  4.3 Oscillations with Large Pre-tension

Figure 9: Bifurcation diagrams of midpoint wall pressure $p_{mid}$ as a
function of the Reynolds number $Re$ at $h=0.01$, $p_{e}=1.1$ for (a) $T=1$
and (b) $T=5$. The steady midpoint pressure $p_{mid}$ is denoted as solid
(stable) and dashed (unstable) lines. The time-averaged midpoint wall pressure
$p_{mid}^{avg}$ is denoted as filled triangles, the upper and lower branch
limit points are denoted as squares.

In order to evaluate the role of pre-tension $T$ in the onset of self-excited
oscillations, Fig. 9 plots a bifurcation diagram showing the time-averaged
midpoint wall pressure compared to the steady midpoint pressure as a function
of the Reynolds number for large pre-tension $T=1$ (Fig. 9a) and $T=5$ (Fig.
9b) for fixed external pressure ($p_{e}=1.1$) and beam thickness ($h=0.01$).
For $T=1$, the steady system exhibits multiple steady solutions (Fig. 2) and
the oscillations initiated close to the upper and lower steady branches merge
into one family of oscillations across the region with multiple steady
solutions (Fig. 9a). However, for much larger pre-tension ($T=5$) the system
exhibits a unique steady state for all Reynolds numbers but this still becomes
unstable to self-excited oscillations at $Re\approx 338.3$. As before, the
time-averaged midpoint pressure is lower than the corresponding steady value
and decreases with the Reynolds number (Fig. 9b). Therefore, the onset of
oscillations is preserved with increasing pre-tension despite the loss of
multiple steady states.

##  5 Discussion

In this study we revisit a theoretical model for flow in a planar collapsible
channel where the flexible wall is modelled as a pre-stressed elastic beam
using a modified nonlinear constitutive law [Wang _et al._ , 2021],
investigating the influence of the external pressure, beam pre-tension and
thickness (a proxy for bending stiffness) on the steady and unsteady behaviour
of the system. The model was solved numerically using the finite element
method.

Similar to the Starling Resistor experiments [Bertram _et al._ , 1990; Bertram
and Castles, 1999] and previous models of flow in collapsible channels and
tubes [Armitstead _et al._ , 1996; Heil, 2004; Stewart, 2017], our model
predicts that the steady system can exhibit multiple co-existing states,
consisting of upper, intermediate and lower steady branches. The model
predicts that this region of multiple steady states is suppressed by
increasing either the wall pre-stress (Fig. 2) or the wall bending stiffness
(Fig. 3). In the former the critical point for multiple steady states is
postponed to larger Reynolds numbers and lower external pressures as the pre-
stress increases (Fig. 2d-f), whereas in the latter the critical point does
not move significantly but the region with multiple steady states narrows as
the bending stiffness increases (Fig. 3b-d).

Previous studies have indicated that self-excited oscillations can
(independently) grow from either the upper and the lower branches of steady
solutions [Wang _et al._ , 2021] in the neighbourhood of the region with
multiple steady states (Fig. 4, 5). Our model predicts that these two families
of oscillations eventually merge into a single family of oscillatory behaviour
for sufficiently low external pressure (Fig. 6). This new merged family of
oscillations retains characteristics of both the upper and lower branch
oscillations (Fig. 7). This new family of oscillations is preserved as the
beam pre-stress increases, despite the loss of multiple steady states (Fig.
9). Conversely, this merged family of oscillations is suppressed by increasing
the beam thickness (i.e. increasing the bending stiffness of the beam, Fig.
8).

## Acknowledgements

We gratefully acknowledge funding from the Chinese Scholarship Council (DYW)
and UK Engineering and Physical Sciences Research Council grants EP/S020950,
EP/S030875 and EP/N014642 (XYL and PSS).

## References

  * Armitstead _et al._ [1996] Armitstead, J. P., Bertram, C. D. and Jensen, O. E. [1996] “A study of the bifurcation behaviour of a model of flow through a collapsible tube,” _Bull. Math. Biol._ 58(4), 611–641.
  * Bertram [1982] Bertram, C. D. [1982] “Two models of instability in a thick-walled collapsible tube conveying a flow,” _J. Biomech._ 15(3), 223–224.
  * Bertram [1986] Bertram, C. D. [1986] “Unstable equilibrium behaviour in collapsible tubes,” _J. Biomech._ 19(1), 61–69.
  * Bertram and Castles [1999] Bertram, C. D. and Castles, R. J. [1999] “Flow limitation in uniform thick-walled collapsible tubes,” _J. Fluids Struct._ 13(3), 399–418.
  * Bertram and Pedley [1982] Bertram, C. D. and Pedley, T. J. [1982] “A mathematical model of unsteady collapsible tube behaviour,” _J. Biomech._ 15(1), 39–50.
  * Bertram _et al._ [1990] Bertram, C. D., Raymond, C. J. and Pedley, T. J. [1990] “Mapping of instabilities for flow through collapsed tubes of differing length,” _J. Fluids Struct._ 4(2), 125–153.
  * Bertram _et al._ [1991] Bertram, C. D., Raymond, C. J. and Pedley, T. J. [1991] “Application of nonlinear dynamics concepts to the analysis of self-excited oscillations of a collapsible tube conveying a fluid,” _J. Fluids Struct._ 5(4), 391–426.
  * Bertram and Tscherry [2006] Bertram, C. D. and Tscherry, J. [2006] “The onset of flow-rate limitation and flow-induced oscillations in collapsible tubes,” _J. Fluids Struct._ 22(8), 1029–1045.
  * Cai and Luo [2003] Cai, Z. X. and Luo, X. Y. [2003] “A fluid–beam model for flow in a collapsible channel,” _J. Fluids Struct._ 17(1), 125–146.
  * Drazin [2002] Drazin, P. G. [2002] _Introduction to hydrodynamic stability_ (Cambridge University, Cambridge).
  * Glendinning [1994] Glendinning, P. [1994] _Stability, instability and chaos: an introduction to the theory of nonlinear differential equations_ (Cambridge University, Cambridge).
  * Grotberg and Gavriely [1989] Grotberg, J. B. and Gavriely, N. [1989] “Flutter in collapsible tubes: a theoretical model of wheezes,” _J. Appl. Physiol._ 66(5), 2262–2273.
  * Hao _et al._ [2016] Hao, Y. J., Cai, Z. X., Roper, S. and Luo, X. Y. [2016] “An Arnoldi-frontal approach for the stability analysis of flows in a collapsible channel,” _Int. J. Appl. Mech._ 8(06), 1650073.
  * Hazel and Heil [2003] Hazel, A. L. and Heil, M. [2003] “Steady finite-reynolds-number flows in three-dimensional collapsible tubes,” _J. Fluid Mech._ 486, 79–103.
  * Heil [2004] Heil, M. [2004] “An efficient solver for the fully coupled solution of large-displacement fluid–structure interaction problems,” _Comput. Methods Appl. Mech. Eng._ 193(1-2), 1–23.
  * Heil and Boyle [2010] Heil, M. and Boyle, J. [2010] “Self-excited oscillations in three-dimensional collapsible tubes: simulating their onset and large-amplitude oscillations,” _J. Fluid Mech._ 652, 405–426.
  * Heil and Hazel [2011] Heil, M. and Hazel, A. L. [2011] “Fluid-structure interaction in internal physiological flows,” _Annu. Rev. Fluid Mech._ 43, 141–162.
  * Herrada _et al._ [2021] Herrada, M. A., Blanco-Trejo, S., Eggers, J. and Stewart, P. S. [2021] “Global stability analysis of flexible channel flow with a hyperelastic wall,” Revised to 2021.
  * Jensen [1990] Jensen, O. E. [1990] “Instabilities of flow in a collapsed tube,” _J. Fluid Mech._ 220, 623–659.
  * Luo _et al._ [2008] Luo, X. Y., Cai, Z. X., Li, W. G. and Pedley, T. J. [2008] “The cascade structure of linear instability in collapsible channel flows,” _J. Fluid Mech._ 600, 45–76.
  * Luo and Pedley [1996] Luo, X. Y. and Pedley, T. J. [1996] “A numerical simulation of unsteady flow in a two-dimensional collapsible channel,” _J. Fluid Mech._ 314, 191–225.
  * Luo and Pedley [2000] Luo, X. Y. and Pedley, T. J. [2000] “Multiple solutions and flow limitation in collapsible channel flows,” _J. Fluid Mech._ 420, 301–324.
  * Stewart [2017] Stewart, P. S. [2017] “Instabilities in flexible channel flow with large external pressure,” _J. Fluid Mech._ 825, 922–960.
  * Strogatz [2018] Strogatz, S. H. [2018] _Nonlinear dynamics and chaos_ (CRC press, Florida).
  * Tantucci [2013] Tantucci, C. [2013] “Expiratory flow limitation definition, mechanisms, methods, and significance,” _Pulm. Med._ 2013, 1–6.
  * Wang _et al._ [2021] Wang, D., Luo, X. and Stewart, P. [2021] “Energetics of collapsible channel flow with a nonlinear fluid-beam model,” _arXiv e-prints_ , arXiv:2106.16234, URL https://arxiv.org/abs/2106.16234.
  * Whittaker _et al._ [2010] Whittaker, R. J., Heil, M., Jensen, O. E. and Waters, S. L. [2010] “Predicting the onset of high-frequency self-excited oscillations in elastic-walled tubes,” _Proc. R. Soc. A_ 466(2124), 3635–3657.


# MultiCast: Zero-Shot Multivariate Time Series Forecasting Using LLMs

Georgios Chatzigeorgakidis  “Athena” Research Center
Greece
<EMAIL_ADDRESS>Konstantinos Lentzos  “Athena” Research Center
Greece
<EMAIL_ADDRESS>Dimitrios Skoutas  “Athena” Research Center
Greece
<EMAIL_ADDRESS>
###### Abstract

Predicting future values in multivariate time series is vital across various
domains. This work explores the use of large language models (LLMs) for this
task. However, LLMs typically handle one-dimensional data. We introduce
MultiCast, a zero-shot LLM-based approach for multivariate time series
forecasting. It allows LLMs to receive multivariate time series as input,
through three novel token multiplexing solutions that effectively reduce
dimensionality while preserving key repetitive patterns. Additionally, a
quantization scheme helps LLMs to better learn these patterns, while
significantly reducing token use for practical applications. We showcase the
performance of our approach in terms of RMSE and execution time against state-
of-the-art approaches on three real-world datasets.

###### Index Terms:

large language models, multivariate time series, forecasting

##  I Introduction

A time series is a sequence of data points, typically recorded at successive
equally spaced intervals of time. These data points can represent various
measurements, observations, or readings taken over time, such as temperature
readings, stock prices, sales figures, or sensor readings. Time series
analysis involves studying the patterns, trends, and relationships present in
the data to understand its behavior over time [1]. Time series forecasting
predicts future values of a time series based on its past observations.

Traditional time series forecasting methods have demonstrated considerable
efficacy over the years and continue to maintain relevance and widespread
adoption in contemporary practice [2]. In general, these methods can be
categorized into linear [3], [4] and non-linear models [5], [6].

Arguably, the most popular traditional time series method is AutoRegressive
Integrated Moving Average (ARIMA) [7]. ARIMA consists of three independent
components; (i) the AutoRegressive (AR) component assumes that the current
value of a time series is a linear combination of its past values, with the
addition of a white-noise term; (ii) the Moving Average (MA) component assumes
that the current value of a time series variable is a linear combination of
past white-noise terms, with no dependence on past values of the variable
itself; (iii) the integrated (I) component incorporates differencing to make
the time series stationary, allowing for the modeling of nonstationary time
series data.

Machine learning and, in particular, deep learning has emerged as a
transformative approach in the field of time series forecasting, offering new
advances [8, 9, 10, 11]. Moreover, pre-training has been used in deep
learning, to significantly accelerate the training process and increase
performance [12]. In domains such as computer vision and Natural Language
Processing (NLP), pre-training facilitates scaling of performance with the
availability of data. However, in the context of time series modeling, access
to sizable pretraining datasets is often limited.

Large Language Models (LLMs) have emerged as a popular tool for Natural
Language Processing (NLP) tasks, and have received considerable attention in
recent years. LLMs are pretrained models, trained on vast amounts of text
data. Their ability to learn rich representations of language has drawn the
attention of the scientific community over the past few years. Specifically,
LLMs are quite capable of capturing syntactic, semantic, and contextual
information [13]. Another interesting aspect of LLMs are emergent abilities
[14], which are capabilities that are not explicitly programmed or designed,
but rather spontaneously emerge from the complex internal processes of the
models. In the past few years, scientists have focused on leveraging the LLMs’
potential to solve problems from other fields than NLP. In particular, in time
series forecasting, by taking advantage of pre-learned representations of
language, LLMs can potentially capture temporal relationships and time series
dynamics [15]. However, most works have focused on univariate time series
forecasting, requiring either fine-tuning [16], or a few-shot prompting
approach [17] (i.e., providing a few examples via prompting to guide the
model’s behavior for a specific task).

In this work, we examine the utility of LLMs for multivariate time series
forecasting via zero-shot prompting (i.e., no additional examples are
provided). To the best of our knowledge, ours is the first work that addresses
this problem.

Our contributions are summarized as follows:

  * •

We introduce three dimensional multiplexing techniques to combine all
dimensions into a single string, passed to an LLM as input.

  * •

We employ SAX quantization on the time series to facilitate inference by the
model and to significantly reduce the computational cost and token usage.

  * •

We present an experimental evaluation against existing traditional, machine
learning, and LLM-based methods for time series forecasting.

##  II Related Work

LLMs have been applied into many different domains and contexts such as
healthcare [18], [19], [20], financial modeling [21], [22], [23], and
education and research [24], [25], as well as in time series data [26] for
many different tasks and application domains [27].

The authors of TIME-LLM [28], introduce a reprogramming framework aimed at
adapting LLMs for time series forecasting without altering their pre-trained
structure. TIME-LLM reprograms input time series into text prototype
representations that suit LLMs’ capabilities. By introducing Prompt-as-Prefix
(PaP), which enriches the input context with natural language instructions,
the reprogrammed input is then processed by the frozen LLM. The output is
projected to generate time series forecasts updating only lightweight input
transformation and output projection parameters, while the backbone language
model remains frozen. Scenarios for both short- and long-term are addressed,
as well as few- and one-shot learning.

LLMTIME [15] is the first approach to apply zero-shot forecasting on time
series using LLMs. The authors argue that the output of LLMs when predicting
digit-by-digit follows a multimodal distribution, which fits well in the case
of time series. To apply forecasting, the time series values are tokenized and
rescaled to a predefined number of digits to use fewer tokens. Then, to apply
forecasting, the time series with their tokenized values separated by commas
are passed to the model. Notably, the model’s output is limited to producing
only digits and commas (i.e., $[0-9,]$). At each time step, a predefined
number of samples is drawn and the final forecast is built using the median of
all samples after descaling the outputted values.

Despite the potential of LLMs for time series forecasting, there are several
limitations that need to be addressed.

  * •

No multivariate support: Most current approaches using LLMs for time series
forecasting focus on univariate time series data. This limitation restricts
the applicability of LLMs to certain types of time series data.

  * •

Fine-tuning requirement: Fine-tuning can be time-consuming and computationally
expensive, particularly for large models. It also requires a substantial
amount of training data, which may not always be available.

  * •

Number of tokens required: LLMs are extremely large models, capable of
efficiently running on computers equipped with GPUs of high capacity in RAM.
Thus, their broad availability depends on services that host such models,
which usually charge queries by token. Consequently, very large queries (e.g.,
a large time series in our context) would be rather expensive to run.

##  III MultiCast

In the following, we describe our approach to zero-shot multivariate time
series forecasting using LLMs. First, we go through the three separate token
multiplexing approaches that we propose. Then, we describe our approach to
reducing complexity using the SAX representation.

###  III-A Dimensional Multiplexing

The dimensional multiplexing process takes place after each dimension has been
rescaled to avoid decimals. Then, each digit is treated separately. An example
of this process is illustrated in the top row of Figure 1. Of course,
depending on the LLM used, its tokenizer must be adapted accordingly, as
discussed in [15]. After multiplexing, the tokens are replaced with their
corresponding corpus id before being passed onto the model for inference. When
the model produces the output, this process is reversed to obtain the final
result. We introduce three separate dimensional multiplexing techniques,
namely (i) digit-interleaving, (ii) value-interleaving, and (iii) value-
concatenation.

Figure 1: The three token multiplexing techniques.

####  III-A1 Digit-Interleaving

After each dimension has been rescaled, the Digit-Interleaving (DI)
multiplexing technique places the digits of each dimension per timestamp
interchangeably. This is exemplified in Figure 1a. Consider a 2-dimensional
time series. Specifically, $d_{1}=[1.7,2.6,...]$ and $d_{2}=[2.3,3.1,...]$ are
the two dimensions (i.e., we only show the first two timestamps for brevity).
After rescaling, the dimensions become $d_{1}=[17,26,...]$ and
$d_{2}=[23,31,...]$, respectively. Then, as described previously, each digit
is considered a separate token. Before being assigned the corresponding corpus
id, tokens are interchangeably placed per dimension for each timestamp,
reducing the dimensions to 1. The resulting series in the example would be
$d=[1273,2361,...]$. Then, each digit (token) and comma are assigned with the
corresponding id. This technique attempts to take advantage of the fact that,
in many multivariate time series, all values are correlated and similarly
scaled. Such an example are z-normalized series, which have zero mean with
values differing a few standard deviations from it. In such a case, the left-
wise digits per dimension will be all placed first; since the model is
producing the output token-by-token, this can help it infer the correct
scaling of the series. More formally, DI multiplexing can be formulated as
follows.

$\begin{split}I_{d}=\\{t_{111}...t_{d11}\quad t_{11b}...t_{d1b}\\}\quad
t_{c}\quad\\{t_{1n1}...t_{dn1}\quad t_{1nb}...t_{dnb}\\}\end{split}$ (1)

where $d$ is the number of dimensions, $b$ the predefined number of digits per
timestamp, and $n$ the time series length.

####  III-A2 Value-Interleaving

Figure 1b shows the Value-Interleaving (VI) dimensional multiplexing
technique. This time, instead of interchangeably placing the digits per
timestamp and dimension, we place the whole values of each dimension per
timestamp one after the other. Thus, in the example, the 1-dimensional result
will be $d=[1723,2631,...]$. Intuitively, this technique is more suitable in
cases where the dimensions of the series are on a different scale. We expect
the model to be able to distinguish between the different dimensions
–especially when they differ in scale–, and manage to internally demultiplex
the input before inference. The VI multiplexing can be formulated as follows.

$\begin{split}I_{t}s=\\{t_{111}...t_{11b}\quad t_{d11}...t_{d1b}\\}\quad
t_{c}\quad\\{t_{1n1}...t_{1nb}\quad t_{dn1}...t_{dnb}\\}\end{split}$ (2)

where $d$ is the number of dimensions, $b$ the predefined number of digits per
timestamp, and $n$ the time series length.

####  III-A3 Value-Concatenation

Finally, Figure 1b shows the Value-Concatenation (VC) dimensional multiplexing
technique which is an extension of the value-interleaving technique; for each
timestamp, we now place the values of each dimension separated by commas, thus
considering them as different values (e.g., in the figure, the 1-dimensional
result will be $d=[17,23,26,31,...]$. We expect this to further faciliate the
internal demultiplexing by the model before detecting any patterns. The VC
multiplexing can be formulated as follows.

$\begin{split}I_{t}d=\\{t_{111}...t_{11b}\\}\hskip 0.80005ptt_{c}\hskip
0.80005pt\\{t_{d11}...t_{d1b}\\}\hskip 0.80005ptt_{c}\hskip
0.80005pt\\{t_{1n1}...t_{1nb}\\}\hskip 0.80005ptt_{c}\hskip
0.80005pt\\{t_{dn1}...t_{dnb}\\}\end{split}$ (3)

where $d$ is the number of dimensions, $b$ the predefined number of digits per
timestamp and $n$ the length of the time series. Of course, in all cases, upon
receiving the multiplexed output from the model, the tokens must be properly
decoded, demultiplexed, and brought back to their initial scale for each
dimension, depending on the selected technique. A significant advantage of
this multiplexing technique against forecasting each dimension separately is
the fact that multivariate time series tend to have high interdimensional
correlations (e.g., temperature and humidity in weather data). We expect that
providing them altogether in the model can lead to the detection of such
interdimensional patterns, yielding better results.

###  III-B Quantization Using SAX

The Symbolic Aggregate approXimation (SAX) is a multi-resolution
representation of a time series introduced in [29]. It can be derived from its
Piecewise Aggregate Approximation (PAA) [30, 31] by quantizing the PAA
segments on the $v$-axis. A time series is first transformed into a PAA
representation of $w$ segments with real-valued coefficients. To obtain a
$SAX$ word for a time series, these coefficients are discretized along the
value axis using breakpoints assuming a $\mathcal{N}(0,1)$ Gaussian
distribution that enables the generation of equiprobable symbols for a given
cardinality. Although bitwise representations were used for these symbols in
the original paper, other encoding types are also possible. Two such popular
alternatives are using alphabetical characters or digits for each symbol.

Forecasting time series is an inherently difficult task due to the nature of
the data. This is also the case for zero-shot foreasting using LLMs, since
(i.e., as also described in [15]) they have to infer a sequence of tokens for
each timestamp, thus simulating a multi-modal distribution. This becomes even
harder when applying the above-mentioned dimensional multiplexing techniques.
Also, for large time series, such a process becomes significantly more
computationally intensive; plus, it requires many tokens, which, depending on
the application, can be rather expensive to infer according to currently LLM
pricing policies. To alleviate these issues, we quantize the time series
across all dimensions in both axes using the SAX representation, before
applying tokenization. We support two different quantization types, either
using an alphabetical or a digital SAX alphabet. Now, each value per timestamp
is consisted of only one token instead of multiple. For example, the time
series in Figure 1 could become $d_{1}^{sax}=[a,b,...]$ and
$d_{2}^{sax}=[b,c,...]$ after alphabetical quantization. We expect that it
will be easier for the model to detect patterns when dealing only with one
token per timestamp.

##  IV Experimental Evaluation

This section presents the results of our experiments. We first explain how we
set up our tests and assess the suggested methods.

###  IV-A Experimental Setup

####  IV-A1 System

We used Python and the Hugging Face API111https://huggingface.co/. The
experiments were run on a server with an AMD Ryzen Threadripper 3960X 24-core
CPU and 256GB memory. The experiments were run on CPU.

####  IV-A2 Datasets

We employ three real-world multivariate time series datasets.

Gas Rate: This is a 2-dimensional dataset containing carbon dioxide (CO2)
emissions. The first dimension contains the input CO2 measurements (ft3/min)
in a gas furnace. The second dimension contains the output CO2 percentage. The
dataset is obtained from the darts library222https://unit8co.github.io/darts.
Of course, the two dimensions are correlated, which makes this dataset ideal
for multivariate forecasting.

Electricity: This multivariate time series is part of the Electricity
Transformer Dataset (ETDataset)333https://github.com/zhouhaoyi/ETDataset. It
contains hourly measurements of various metrics, which were resampled on a
3-day basis, for a total of 242 timestamps. From this dataset, we extracted 3
dimensions of electricity measurements, specifically the High UseFul Load
(HUFL), High UseLess Load (HULL), and Oil Temperature (OT). Again, the
dimensions are correlated; specifically, OT is used as a target variable in
regression problems.

Weather: The weather dataset was generated by the Max Planck
Institute444https://www.bgc-jena.mpg.de/wetter/ and contains 21 weather-
related metrics obtained from a weather station located in Germany. From the
21 variables, we extracted the air temperatures (Tlog) measured in Celsius
degrees, the water vapor concentration (H2OC) measured in mmol/mol, the
saturation water vapor pressure (VPmax), measured in mbar, and the potential
temperature (Tpot) measured in Kelvin degrees. Again, being weather-related,
all dimensions are correlated.

TABLE I: Datasets. Dataset | Dimensions | Length
---|---|---
Gas Rate | 2 | 296
Electricity | 3 | 242
Weather | 4 | 217
TABLE II: Parameters. Parameter | Range
---|---
Dimensions |  2, 3, 4
Number of samples |  5, 10, 20
SAX segment length |  3, 6, 9
SAX alphabet size |  5, 10, 20

####  IV-A3 Competitors

We evaluate the following methods:

  * •

MultiCast (DI): MultiCast using the digit-interleaving dimensional
multiplexing method.

  * •

MultiCast (VI): MultiCast using the value-interleaving dimensional
multiplexing method on the same value.

  * •

MultiCast (VC): MultiCast using the value-concatenation dimensional
multiplexing method on consecutive values.

  * •

LLMTIME: The state-of-the-art in LLM-based zero-shot time series forecasting
(i.e., applied in each dimension separately).

  * •

ARIMA: Autoregressive Integrated Moving Average (ARIMA) is one of the most
widely used univariate time series forecasting methods.

  * •

LSTM [32]: Termed as Long-Short-Term Memory (LSTM), LSTMs are Recurrent Neural
Networks (RNNs) designed to handle the vanishing gradient problem. This
ability allows LSTMs to learn and remember information over time, making them
ideal for time series forecasting. LSTMs have been used successfully for
multivariate time series forecasting [33, 34].

####  IV-A4 Parameters

The parameters utilized in our experimental assessment are listed in Table II.
For each parameter, we performed tuning tests to establish their ranges and
default values, which are highlighted in bold within the table. More
specifically, the dimensions parameter corresponds to the dimensionality of
each dataset; the number of samples only applies to the LLM-related models and
is the number of inference values taken for each timestamp; the SAX segment
length is the level of quantization on the x-axis, which determines the level
of compression of a time series; the SAX alphabet size is the level of
quantization on the y-axis, as performed by the SAX method. Regarding the LSTM
parametrization, we performed a grid search, which yielded a 1-hidden-level
network of 128 units and a dropout rate of 0.2. It was trained for 30 epochs
using the Adam [35] optimizer with the Mean Squared Error (MAE) as loss
function.

####  IV-A5 Metrics

In accordance with standard practices in time series forecasting, the Root
Mean Squared Error (RMSE) metric was employed to evaluate our methods. RMSE is
formulated as $\sqrt{\Sigma_{i=1}^{n}{({y_{i}-\hat{y}_{i}})^{2}}/{n}}$, where
$y_{i}$ is the actual value, $\hat{y}_{i}$ is the predicted value at timestamp
$i$ and $n$ is the number of timestamps on which forecasting was applied.

###  IV-B LLM Model Selection

MultiCast can be used with any LLM to apply multivariate time series
forecasting. In the following, we evaluate its accuracy using LLaMA2 (i.e.,
the 7B parameter variant) and Phi-2 [36] as back-end models. LLaMA2 is one of
the most popular LLMs that achieves good performance with fewer parameters.
Phi-2 is a math-oriented LLM (i.e., 2.7B parameters), tailored to solving math
problems. Table III lists the forecast RMSE for the Gas Rate data set in both
dimensions for both LLMs. In both cases, the VI variant of MultiCast was used.
LLaMA2 achieves better performance (i.e., approx. twice as good) in all cases.
This can be attributed to the fewer parameters of the Phi-2 model; while it is
math-oriented and quite capable of solving complex problems described by
textual prompting, it seems to not properly detect the patterns in the series,
leading to larger errors.

TABLE III: LLM model comparison. Model | Dimension
---|---
GasRate | CO2
MultiCast (LLaMA2 / 7B) | 1.154 | 2.71
MultiCast (Phi-2 / 2.7B) | 2.106 | 4.676

Figures 2a and b depict two indicative examples of forecasting the first
dimension of the Gas Rate dataset using the LLaMA2 and Phi-2 models,
respectively. Clearly, the LLaMA2 model performs better, being able to
properly follow the upward trend of the time series and even infer two local
maxima of the original time series. Phi-2, on the other hand, fails to
accurately forecast the time series in this dimension; while it seems to
successfully detect the upward trend, its entire output is shifted 1 to 2
units on the y-axis. Since LLaMA2 seems to perform significantly better in all
cases, for the rest of the experiential evaluation, we will be using LLaMA2 as
the back-end model for MultiCast.

(a) LLaMA2 GasRate prediction

(b) Phi-2 GasRate prediction

Figure 2: Comparison of the two models.

###  IV-C Forecasting Accuracy

Next, we compare the prediction precision in terms of RMSE of all MultiCast
variants against the rest of the competitor approaches. Table IV lists the
results for the Gas Rate dataset. To better comprehend the insights behind the
results and acquire knowledge regarding the differences in forecasting ability
between the LLM-based models and the rest of the competition, for each
dimension, we denote the first best overall performance using bold font and
the second best using italic font. Interestingly, for the GasRate dimension,
the best overall approach was LLMTIME ($0.703$), followed by MultiCast (DI)
with $0.781$. The LLM-based approaches all seem to cope well with detecting
the underlying patterns for this dimension, thus producing good results. The
case is different for the second dimension (CO2), where the conventional
methods seem to yield a better overall performance, with ARIMA being the best
(2.63). MultiCast (VI) was the second-best overall and the best LLM-based
performer (2.71).

TABLE IV: Forecasting RMSE for the Gas Rate dataset. Model | Dimension
---|---
GasRate | CO2
MultiCast (DI) | 0.781 | 4.639
MultiCast (VI) | 1.154 | 2.71
MultiCast (VC) | 0.965 | 3.626
LLMTIME | 0.703 | 2.75
ARIMA | 0.92 | 2.63
LSTM | 1.122 | 3.89
TABLE V: Forecasting RMSE for the Electricity dataset. Model | Dimension
---|---
HUFL | HULL | OT
MultiCast (DI) | 5.914 | 1.444 | 9.198
MultiCast (VI) | 8.63 | 1.882 | 13.752
MultiCast (VC) | 2.424 | 1.913 | 10.230
LLMTIME | 4.299 | 1.432 | 7.543
ARIMA | 7.063 | 1.572 | 4.181
LSTM | 4.892 | 1.43 | 8.740
TABLE VI: Forecasting RMSE for the Weather dataset. Model | Dimension
---|---
Tlog | H2OC | VPmax | Tpot
MultiCast (DI) | 3.711 | 2.43 | 3.025 | 6.888
MultiCast (VI) | 3.26 | 2.122 | 2.387 | 11.352
MultiCast (VC) | 4.983 | 3.819 | 5.776 | 5.993
LLMTIME | 3.14 | 1.746 | 4.044 | 6.981
ARIMA | 3.324 | 2.686 | 4.331 | 6.067
LSTM | 3.524 | 1.796 | 2.708 | 5.559

Figure 3 depicts two indicative forecast outputs of the best MultiCast
approach (DI) for the first dimension of the Gas Rate data set against the
corresponding ARIMA result. Both seem to yield a good result here; MultiCast
seems to properly detect a continuously upward trend in the time series;
however, the result seems to have larger variance than that of the original
time series. On the other hand, the ARIMA approach does not clearly follow the
upward trend; however, its variance seems to be on par with the one of the
original time series.

(a) Multicast for the gas rate dataset (GasRate).

(b) ARIMA for the gas rate dataset (GasRate).

Figure 3: MultiCast (DI) versus ARIMA for the GasRate dimension.

Table V lists the results for the Electricity dataset. For the HUFL dimension,
MultiCast (VC) seems to yield significantly better RMSE than the rest of the
approaches. However, the rest of the MultiCast variants do not cope as well.
For the HULL dimension, all approaches seem to produce good results, with
LLMTIME achieving the best RMSE. Finally, for the OT dimension, ARIMA performs
significantly better than the competition. The MultiCast approaches do not
perform well. LLMTIME is the best among LLM-based models. This suggests a
possible drop in the performance of MultiCast as the dimensionality of the
time series increases since there is the extra step of demultiplexing the
input that the LLMs must infer. However, the error in the best LLM-based model
(9.198) is very close to that of the LSTM model (8.740).

Figure 4 illustrates an indicative example of the MultiCast (VC) forecast
output (Figure 4a) against the LSTM (Figure 4b) for the HUFL dimension of the
electricity data set. Clearly, MultiCast manages to correctly infer both the
trend and variance of the time series. On the other hand, the LSTM seems to
perform rather poorly, falsely yielding a non-existent linear upward trend.

(a) Multicast for the electricity dataset (HUFL).

(b) LSTM for the electricity dataset (HUFL).

Figure 4: MultiCast (VC) versus LSTM for the HUFL dimension.

The RMSE results for the Weather dataset are listed in Table VI. LLMTIME
achieves the best performance in the Tlog dimension, though, all approaches
except MultiCast (VC) are close. This is also the case for the H2OC dimension.
For the VPmax dimension, the best overall approach was MultiCast (VI), with
MultiCast (VC) again performing worse than the rest. However, this is reversed
in the Tpot dimension, where the MultiCast variant (VC) yields the best
performance among all LLM-based approaches. LSTM is the better performer in
this dimension. Notice that the degradation in forecasting accuracy for more
dimensions is not present in this case; the MultiCast variants are all either
close to, or outperform the rest in all dimensions. Another key takeaway here
is that the optimal multiplexing method differs from dimension to dimension
and from dataset to dataset. A comprehensive analysis on which dataset
characteristics cause this behavior is an interesting future work.

As in the rest of the cases, an indicative example of MultiCast against a
conventional method is illustrated in Figure 5. Clearly, the DI variant of
MultiCast (Figure 5a) yields better results than ARIMA (Figure 5b) here, able
to accurately estimate the upward trend and fluctuation at the end of the time
series.

(a) Multicast for the weather dataset (Tlog).

(b) ARIMA for the weather dataset (Tlog).

Figure 5: MultiCast (VI) versus ARIMA for the Tlog dimension.

Overall, we notice a trade-off when using MultiCast for multivariate time
series forecasting, as opposed to LLMTIME. Forecasting each dimension
separately using LLMTIME will completely ignore the interdimensional
correlations, which is not desirable in such scenarios. On the other hand,
MultiCast poses an additional challenge to the LLM models, which now have to
also infer the demultiplexing of the dimensions. Both cases hinder the
accuracy of the obtained result. Having in mind the interesting aspect of
emergent abilities, we argue that using very large LLMs (e.g., GPT-4, Gemini)
will further improve MultiCast’s performance.

###  IV-D Increasing Number of Samples

Table VII lists the accuracy in terms of RMSE of all LLM-based models for an
increasing number of samples. As a reminder, all LLM-based models draw several
samples of the values of each timestamp, and the final estimated value is
derived by computing the median among all samples. The LLMTIME approach seems
to produce better results for 5 and 10 samples. Interestingly, the error
worsens for 20 samples. This could be because the inherent variance of the
produced series tends to be averaged out as we draw more samples. However,
this is not the case for the MultiCast method; all three approaches seem to
produce better results for more samples, and the MultiCast DI variant achieves
the best performance for 20 samples. A drawback of drawing many samples is the
performance deterioration in execution time (i.e., each execution time is
listed below each corresponding RMSE in Table VII). Notice that, in all cases,
the execution time doubles when the number of samples is doubled, which is
expected since the model must infer twice as many tokens. Interestingly, the
LLMTIME requires slightly less total time (i.e., sum of time needed per
dimension) than its MultiCast counterparts, since the latter also need to
infer the multiplexing/demultiplexing of the tokens.

TABLE VII: Performance for an increasing number of samples. Method | Number of samples
---|---
5 | 10 | 20
MultiCast (DI) |  0.781 |  |
1036 sec |  0.762 |  |
2050 sec |  0.592 |  |
4159 sec |  |  |
MultiCast (VI) |  0.965 |  |
1041 sec |  1.302 |  |
2068 sec |  0.877 |  |
4131 sec |  |  |
MultiCast (VC) |  1.154 |  |
1168 sec |  0.704 |  |
2468 sec |  0.63 |  |
4981 sec |  |  |
LLMTIME |  0.703 |  |
1023 sec |  0.606 |  |
1939 sec |  0.842 |  |
3684 sec |  |  |

###  IV-E Performance Boost Using SAX

(a) 3 SAX segments.

(b) 6 SAX segments.

(c) 9 SAX segments.

Figure 6: Forecasting for various SAX segments (CO2%).

(a) 5 SAX symbols.

(b) 10 SAX symbols.

(c) 20 SAX symbols.

Figure 7: Forecasting for different SAX alphabet sizes (CO2%). Figure 8:
Forecasting using digits instead of letters as symbols.

Next, we will show the results obtained when quantization is applied using the
SAX method, as described in Section III-A. Specifically, we evaluate the
effects of increasing the length of the SAX segment and the size of the
alphabet on the performance of zero-shot time series forecasting using
MultiCast.

####  IV-E1 Increasing SAX Segment length

Table VIII lists the results for an increasing number of SAX segments for the
CO2$\%$ dimension of the Gas Rate dataset, in terms of RMSE and execution
time. We also list the results for using a different kind of SAX quantization;
either using alphabetical characters or digits to encode SAX words.
Compressing the time series significantly facilitates the inference process
since now the model has to generate only one symbol per timestamp, instead of
three or more. This is reflected in the execution times shown in Table VIII;
inference after applying SAX compression is more than an order of magnitude
faster, from 52 seconds in the best case (i.e., using 9 SAX segments) to 1168
seconds, when no quantization is applied. The large difference in performance
can have a big impact on forecasting tasks that are run on CPU, which may
often be the case in scenarios where access to a GPU with large enough memory
to fit an LLM is not possible.

TABLE VIII: Increasing SAX segment length. Method | SAX Segment Length
---|---
3 | 6 | 9
MultiCast SAX (alphabetical) |  1.089 |  |
148 sec |  0.983 |  |
77 sec |  0.888 |  |
54 sec |  |  |
MultiCast SAX (digital) |  0.992 |  |
156 sec |  0.99 |  |
71 sec |  0.912 |  |
52 sec |  |  |
MultiCast |  0.781
1168 sec |  |  |

As expected, quantizing the time series leads to a loss of information. Again,
this is reflected in the RMSE scores for the SAX approaches in Table VIII,
which are worse than when no quantization is applied. However, this may not
always be the case; having to infer only one symbol per timestamp is easier
for the LLM. Patterns, if they exist, will be easier to detect and guess. The
higher RMSE scores in these cases can be attributed to the quantization that
SAX applies on both axes. However, the final result, when plotted, could
properly follow the initial time series. This effect is illustrated in Figures
6b and 6c, for the CO2$\%$ dimension of the Gas Rate dataset. On the other
hand, in this case, MultiCast using 3 SAX segments managed to detect the
initial upward trend (Figure 6a), but the result worsened afterwards.

Figure 8 shows an indicative example of the prediction result when applying
SAX quantization using digits to encode symbols, for the CO2$\%$ dimension of
the Gas Rate dataset. It is easily noticeable that the resulting line (red in
the figure) closely follows the initial time series in this dimension. This
could suggest that it may be easier for the LLMs to detect patterns in time
series represented by numbers rather than alphabetical characters.

####  IV-E2 Increasing SAX Alphabet Size

In the following, we evaluate the performance of MultiCast, in terms of RMSE
and execution time, when increasing the size of the SAX alphabet. Table IX
lists the results. We should note that for digits we can only go up to an
alphabet of size 10. Again, the non-quantized MultiCast yields better
performance in terms of RMSE, but is significantly slower. Increasing the size
of the alphabet does not seem to affect the execution time. Also, in terms of
RMSE, larger alphabet sizes produced higher errors, possibly due to the
increase in complexity that the use of more symbols induces.

TABLE IX: Increasing SAX alphabet size. Method | SAX Alphabet Size
---|---
5 | 10 | 20
MultiCast SAX (alphabetical) |  0.983 |  |
77 sec |  1.198 |  |
81 sec |  1.273 |  |
83 sec |  |  |
MultiCast SAX (digital) |  0.99 |  |
71 sec |  1.21 |  |
75 sec | N/A |  |
MultiCast |  0.781
1168 sec |  |  |

Finally, Figure 7 shows an indicative forecasting example for 5, 10 and 20 SAX
symbols for the CO2$\%$ dimension of the Gas Rate dataset. The drop in RMSE
scores is also reflected here; only in the case of using five symbols does the
forecasted time series follow the trend of the original.

##  V Conclusions

In this paper, we presented MultiCast, an approach that leverages LLMs for
zero-shot multivariate time series forecasting. To make this model work with
multiple dimensions, we proposed three token multiplexing solutions that
reduce the dimensionality of the time series to one. This allows the time
series to be compatible with the input of an LLM, while still preserving its
ability to detect repetitive patterns. Additionally, we presented a
quantization solution that aims to facilitate the learning of existing
patterns in the series by LLMs. This solution also significantly reduces the
execution time. In our comprehensive experimental analysis using three real-
world datasets, we found that the use of LLMs for multivariate zero-shot time
series forecasting shows promise and offers a significant advantage compared
to other similar methods available in the literature: No expert knowledge or
time and resource-consuming parameter search processes are required. In the
future, we plan to expand our research on employing LLMs for zero-shot
solutions on other similar time series-related tasks, such as imputation,
anomaly detection, and change point detection, as well as evaluate MultiCast’s
inference performance using more LLMs as back-end models and further improving
it in more dimensions.

## Acknowledgment

This work was partially funded by the EU Horizon Europe projects STELAR
(101070122) and DT4GS (101056799).

## References

  * [1] R. H. Shumway, D. S. Stoffer, and D. S. Stoffer, _Time series analysis and its applications_. Springer, 2000, vol. 3.
  * [2] A. Zeng, M. Chen, L. Zhang, and Q. Xu, “Are transformers effective for time series forecasting?” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 37, no. 9, pp. 11 121–11 128, Jun. 2023.
  * [3] Z. Liu, Z. Zhu, J. Gao, and C. Xu, “Forecast methods for time series data: A survey,” _IEEE Access_ , vol. 9, pp. 91 896–91 912, 2021.
  * [4] J. G. De Gooijer and R. J. Hyndman, “25 years of time series forecasting,” _International Journal of Forecasting_ , vol. 22, no. 3, pp. 443–473, 2006\.
  * [5] J. Kuvulmaz, S. Usanmaz, and S. N. Engin, “Time-series forecasting by means of linear and nonlinear models,” in _MICAI 2005: Advances in Artificial Intelligence_ , A. Gelbukh, Á. de Albornoz, and H. Terashima-Marín, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005, pp. 504–513.
  * [6] C. Cheng, A. Sa-Ngasoongsong, O. Beyca, T. Le, H. Yang, Z. Kong, and S. Bukkapatnam, “Time series forecasting for nonlinear and non-stationary processes: a review and comparative study,” _IIE Transactions_ , vol. 47, no. 10, pp. 1053–1071, 2015.
  * [7] G. E. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung, _Time series analysis: forecasting and control_. John Wiley & Sons, 2015.
  * [8] K. e. a. Benidis, “Deep learning for time series forecasting: Tutorial and literature survey,” _ACM Comput. Surv._ , vol. 55, no. 6, dec 2022. [Online]. Available: https://doi.org/10.1145/3533382
  * [9] J. F. Torres, D. Hadjout, A. Sebaa, F. Martínez-Álvarez, and A. Troncoso, “Deep learning for time series forecasting: A survey,” _Big Data_ , vol. 9, no. 1, pp. 3–21, 2021.
  * [10] A. Mahmoud and A. Mohammed, _A Survey on Deep Learning for Time-Series Forecasting_. Springer International Publishing, 2021, pp. 365–392.
  * [11] S. Du, T. Li, Y. Yang, and S.-J. Horng, “Multivariate time series forecasting via attention-based encoder–decoder framework,” _Neurocomputing_ , vol. 388, pp. 269–279, 2020.
  * [12] J. Jiang, Y. Shu, J. Wang, and M. Long, “Transferability in deep learning: A survey,” 2022.
  * [13] W. X. Zhao, K. Zhou, J. Li, T. Tang, X. Wang, Y. Hou, Y. Min, B. Zhang, J. Zhang, Z. Dong _et al._ , “A survey of large language models,” _arXiv preprint arXiv:2303.18223_ , 2023.
  * [14] J. Wei, Y. Tay, R. Bommasani, C. Raffel, B. Zoph, S. Borgeaud, D. Yogatama, M. Bosma, D. Zhou, D. Metzler _et al._ , “Emergent abilities of large language models,” _arXiv preprint arXiv:2206.07682_ , 2022.
  * [15] N. Gruver, M. Finzi, S. Qiu, and A. G. Wilson, “Large language models are zero-shot time series forecasters,” in _NeurIPS_ , 2023.
  * [16] C. Chang, W.-C. Peng, and T.-F. Chen, “Llm4ts: Two-stage fine-tuning for time-series forecasting with pre-trained llms,” _arXiv preprint arXiv:2308.08469_ , 2023.
  * [17] H. Xue and F. D. Salim, “Promptcast: A new prompt-based learning paradigm for time series forecasting,” _IEEE Transactions on Knowledge and Data Engineering_ , 2023.
  * [18] P. Schmiedmayer, A. Rao, P. Zagar, V. Ravi, A. Zahedivash, A. Fereydooni, and O. Aalami, “Llm on fhir–demystifying health records,” _arXiv preprint arXiv:2402.01711_ , 2024.
  * [19] A. J. Nashwan, A. A. AbuJaber, and A. AbuJaber, “Harnessing the power of large language models (llms) for electronic health records (ehrs) optimization,” _Cureus_ , vol. 15, no. 7, 2023.
  * [20] M. Wornow, A. Lozano, D. Dash, J. Jindal, K. W. Mahaffey, and N. H. Shah, “Zero-shot clinical trial patient matching with llms,” _arXiv preprint arXiv:2402.05125_ , 2024.
  * [21] Y. Li, S. Wang, H. Ding, and H. Chen, “Large language models in finance: A survey,” in _Proceedings of the Fourth ACM International Conference on AI in Finance_ , 2023, pp. 374–382.
  * [22] H. Yang, X.-Y. Liu, and C. D. Wang, “Fingpt: Open-source financial large language models,” _arXiv preprint arXiv:2306.06031_ , 2023.
  * [23] H. Zhao, Z. Liu, Z. Wu, Y. Li, T. Yang, P. Shu, S. Xu, H. Dai, L. Zhao, G. Mai _et al._ , “Revolutionizing finance with llms: An overview of applications and insights,” _arXiv preprint arXiv:2401.11641_ , 2024.
  * [24] M. Hosseini, C. A. Gao, D. M. Liebovitz, A. M. Carvalho, F. S. Ahmad, Y. Luo, N. MacDonald, K. L. Holmes, and A. Kho, “An exploratory survey about using chatgpt in education, healthcare, and research,” _medRxiv_ , pp. 2023–03, 2023.
  * [25] S. Moore, R. Tong, A. Singh, Z. Liu, X. Hu, Y. Lu, J. Liang, C. Cao, H. Khosravi, P. Denny _et al._ , “Empowering education with llms-the next-gen interface and content generation,” in _International Conference on Artificial Intelligence in Education_. Springer, 2023, pp. 32–37.
  * [26] Y. Jiang, Z. Pan, X. Zhang, S. Garg, A. Schneider, Y. Nevmyvaka, and D. Song, “Empowering time series analysis with large language models: A survey,” 2024\.
  * [27] X. Zhang, R. R. Chowdhury, R. K. Gupta, and J. Shang, “Large language models for time series: A survey,” 2024.
  * [28] M. Jin, S. Wang, L. Ma, Z. Chu, J. Y. Zhang, X. Shi, P.-Y. Chen, Y. Liang, Y.-F. Li, S. Pan, and Q. Wen, “Time-LLM: Time series forecasting by reprogramming large language models,” in _The Twelfth International Conference on Learning Representations_ , 2024. [Online]. Available: https://openreview.net/forum?id=Unb5CVPtae
  * [29] J. Shieh and E. J. Keogh, “ _i_ SAX: indexing and mining terabyte sized time series,” in _SIGKDD_ , 2008, pp. 623–631.
  * [30] E. J. Keogh, K. Chakrabarti, M. J. Pazzani, and S. Mehrotra, “Dimensionality reduction for fast similarity search in large time series databases,” _Knowl. Inf. Syst._ , vol. 3, no. 3, pp. 263–286, 2001.
  * [31] B. Yi and C. Faloutsos, “Fast time sequence indexing for arbitrary Lp norms,” in _VLDB_ , 2000, pp. 385–394.
  * [32] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, pp. 1735–80, 12 1997.
  * [33] S. Alhirmizy and B. Qader, “Multivariate time series forecasting with lstm for madrid, spain pollution,” in _2019 international conference on computing and information science and technology and their applications (ICCISTA)_. IEEE, 2019, pp. 1–5.
  * [34] J. Ju and F.-A. Liu, “Multivariate time series data prediction based on att-lstm network,” _Applied sciences_ , vol. 11, no. 20, p. 9373, 2021.
  * [35] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [36] M. Javaheripi, S. Bubeck, M. Abdin, J. Aneja, S. Bubeck, C. C. T. Mendes, W. Chen, A. Del Giorno, R. Eldan, S. Gopi _et al._ , “Phi-2: The surprising power of small language models,” 2023.