# Parameter Refinement of a Ballbot and Predictive Control for Reference
Tracking with Linear Parameter-Varying Embedding

1st Dimitrios S. Karachalios This work was funded by the German Research
Foundation (DFG), project number 419290163. Institute for Electrical
Engineering in Medicine
University of Luebeck
Lübeck, Germany
email<EMAIL_ADDRESS>2nd Hossam S. Abbas  Institute
for Electrical Engineering in Medicine
University of Luebeck
Lübeck, Germany
email<EMAIL_ADDRESS>
###### Abstract

In this study, we implement a control method for stabilizing a ballbot that
simultaneously follows a reference. A ballbot is a robot balancing on a
spherical wheel where the single point of contact with the ground makes it
omnidirectional and highly maneuverable but with inherent instability. After
introducing the scheduling parameters, we start the analysis by embedding the
nonlinear dynamic model derived from first principles to a linear parameter-
varying (LPV) formulation. Continuously, and as an extension of a past study,
we refine the parameters of the nonlinear model that enhance significantly its
accuracy. The crucial advantages of the LPV formulation are that it consists
of a nonlinear predictor that can be used in model predictive control (MPC) by
retaining the convexity of the quadratic optimization problem with linear
constraints and further evades computational burdens that appear in other
nonlinear MPC methods with only a slight loss in performance. The LPVMPC
control method can be solved efficiently as a quadratic program (QP) that
provides timing that supports real-time implementation. Finally, to illustrate
the method, we test the control designs on a two-set point 1D non-smooth
reference with sudden changes, to a 2D nonstationary smooth reference known as
Lissajous curves, and to a single-set point 1D non-smooth reference where for
this case theoretical guarantees such as stability and recursive feasibility
are provided.

###### Index Terms:

Identification, Control, Ballbot, Model Predictive Control, Linear Parameter-
Varying, Quadratic Program

##  I Introduction

A ballbot Fig. 1 is an omnidirectional mobile robot that balances on a single
spherical wheel [1, 2, 3]. The single point of contact with the ground makes
this under-actuated system agile but challenging to control. The use cases of
a ballbot include healthcare, virtual assistants, etc., as service robots,
where maneuverability in busy environments while maintaining stability is
important. A ballbot system’s stability implies balancing and trajectory
tracking through a predetermined path. Linear control of a ballbot [1] is
proven insufficient, and the LPV model, being much closer to the nonlinear
model, has shown great potential in several applications, such as autonomous
vehicles [4].

The dynamic model structure is well-known for robotic systems based on
physical laws. The identification of the parameters can be obtained through
classical mechanics [5] and enhanced white modeling, but parameters that
describe friction can be determined only experimentally. This study will
provide a method that refines the parameters that define the nonlinear ballbot
model from the measurements derived in [1]. For controlling the ballbot,
often, PID controllers considered [6, 1] or with double-loop approaches [3].
One of the limitations of the aforementioned methods is the consideration of
input/state constraints.

Model predictive control (MPC) has been widely applied for reference tracking
with constraints [4]. It can generate optimal trajectories that steer the
robot in a constrained space. Since autonomously moving robots are safety-
critical systems, nonlinear MPC (NMPC) is gaining popularity due to its
ability to utilize high-fidelity nonlinear models, enabling more accurate and
precise control actions. Another advantage of MPC is that it can be leveraged
with planning algorithms, which is crucial for such autonomously moving
objects. The drawback of NMPC is the computational complexity that makes
online minimization of the underlying objective function over a nonconvex
manifold cumbersome. Therefore, a good alternative is to utilize the LPV
formulation that can embed the nonlinear model equivalently and introduce the
LPVMPC framework that can be solved efficiently as QP, with nearly the same
computational burden as in linear MPC, enabling real-time implementation.

Figure 1: The Ballbot constructed at the Institute for Electrical Engineering
in Medicine (IME) and side view on the $xz$-plane.

Contents: Sec. II contains preliminaries and definitions. In Sec. II-A, the
nonlinear ballbot model is presented. Continuously, in Sec II-B, the linear
parameter-varying (LPV) model is presented with technical aspects. In Sec.
II-C, we refer to other stabilization methods that will help tune our terminal
cost, and in Sec. II-D, the LPVMPC complete framework is introduced with all
the details along with algorithms for ease of implementation. Sec. III
contains results; in particular, in Sec.III-A, the nonlinear model is
identified. In Sec. III-B, three representative scenarios on reference
tracking illustrate the developed method together with results that can ensure
theoretical guarantees. Finally, a discussion followed by acknowledgment
concludes the study.

Notations: $I_{n}$ is the identity matrix of dimension $n$. The notation
$Q\succeq 0$ represents the positive semi-definiteness of a matrix $Q$. The
weighted norm $\|x\|_{Q}$ is defined as $\|x\|_{Q}^{2}=x^{\top}Qx$ and
similarly $\lVert x\rVert^{2}_{P}=x^{\top}Px,~{}\lVert
u\rVert^{2}_{R}=u^{\top}Ru$. The function $\texttt{diag}(\mathbf{x})$
constructs a diagonal square matrix from a vector $\mathbf{x}$. The set of
positive integers, including zero, is denoted by $\mathbb{Z}_{+}\cup\\{0\\}$;
$i|k$ denotes at time $k$ the $i^{\textrm{th}}$ prediction.

##  II Preliminaries and Definitions

###  II-A Nonlinear Model of a Ballbot

The ballbot model is thoroughly discussed regarding the Euler-Lagrange method
of formulating the dynamics in [1, 2]. A first attempt towards model discovery
was performed in [1], where a linearized model was used to balance the
ballbot. The ballbot is a $5$-DoF system, where three are the body’s
rotational motion, and the two are the translation in the $xy$-plane. We will
see that the body’s angular position $\theta$ is used mainly for balancing,
and the angular position of the ball $\phi$ is used mainly for tracking see
Fig. 1.

The symmetry of the ballbot control problem helps us focus independently on
just the $xz$-plane [1] to simulate and understand its stability and
controllability. The ballbot state space model in (1) gives us the two
observable and controllable states as $q=[\phi_{y},\theta_{y}]^{T}$, together
with $\tau_{y}$, being the control input torque of the virtual wheel actuator
of the ballbot in the $xz$-plane as in Fig. 1 (right).

The ballbot model consists of the following matrices $M,C,D,G$, that denote
the mass, Coriolis force (centrifugal forces), frictional torque, and
gravitational forces, respectively, together with the input vector
$\widetilde{B}$. The nonlinear mechanical model is

$\displaystyle
M(q)\ddot{q}+C(q,\dot{q})+D(\dot{q})+G(q)=\widetilde{B}\tau_{y},$ (1)

where the matrices are defined as:

$\displaystyle\begin{aligned} M&=\begin{bmatrix}b_{1}&-b_{2}+\ell
r_{b}\cos{\theta_{y}}\\\ -b_{2}+\ell
r_{b}\cos{\theta_{y}}&b_{3}\end{bmatrix},\\\ C&=\begin{bmatrix}-\ell
r_{b}\sin({\theta_{y})\dot{\theta}_{y}^{2}}\\\
0\end{bmatrix},~{}D=\begin{bmatrix}b_{4}\dot{\phi}_{y}\\\ 0\end{bmatrix},\\\
G&=\begin{bmatrix}0\\\ -\ell
g\sin({\theta_{y}})\end{bmatrix},~{}\widetilde{B}=\begin{bmatrix}\frac{r_{b}}{r_{w}}\\\
-\frac{r_{b}}{r_{w}}\end{bmatrix}.\end{aligned}$ (2)

The nonlinear state-space representation of the planar ballbot model is given
next in (3).

$\displaystyle\underbrace{\begin{bmatrix}\dot{q}\\\
\ddot{q}\end{bmatrix}}_{\dot{x}}=\underbrace{\begin{bmatrix}\dot{q}\\\
M^{-1}(-C-D-G-\epsilon)\end{bmatrix}+\begin{bmatrix}0\\\
M^{-1}\widetilde{B}\end{bmatrix}\tau_{y}}_{f(x,u)},$ (3)

where, $x=[\phi_{y},\theta_{y},\dot{\phi}_{y},\dot{\theta}_{y}]^{\top}$ is the
state vector of the system with $\phi_{y}$ the angular displacement,
$\theta_{y}$ the tilt angle, $\dot{\phi}_{y}$ the angular speed of the ballbot
and $\dot{\theta}_{y}$ the angular speed of the tilt angle. The input to each
plane, e.g., $\tau_{y}$, is the virtual wheel’s torque Fig. 1. The virtual
torques $(\tau_{x},\tau_{y},\tau_{z})$ can be converted to real torques of the
three motors $(\tau_{1},\tau_{2},\tau_{3})$ as in [1] by using (4), where for
the design in Fig. 1 it holds that the Zenith angle is
$\alpha=\pi/4=45^{\circ}$.

$\begin{bmatrix}\tau_{1}\\\ \tau_{2}\\\
\tau_{3}\end{bmatrix}=\frac{1}{3}\begin{bmatrix}\frac{2}{\cos{\alpha}}&0&\frac{1}{\sin{\alpha}}\\\
-\frac{1}{\cos{\alpha}}&\frac{\sqrt{3}}{\cos{\alpha}}&\frac{1}{\sin{\alpha}}\\\
-\frac{1}{\cos{\alpha}}&-\frac{\sqrt{3}}{\cos{\alpha}}&\frac{1}{\sin{\alpha}}\end{bmatrix}\begin{bmatrix}\tau_{x}\\\
\tau_{y}\\\ \tau_{z}\end{bmatrix}.$ (4)

Now that we have a well-defined nonlinear model, we can introduce an
equivalent representation with an LPV embedding.

###  II-B The Linear Parameter-Varying (LPV) Embedding

The linear parameter-varying (LPV) embedding [7, 8, 9, 10] can be seen as a
structured linearization over an adaptive operational point that depends on
the scheduling parameters. The nonlinear system in (3), can be embedded in the
continuous-time LPV representation by denoting with the subscript ”c” the
continuous operators that depend on the variable $\rho$ as $A_{c}(\rho)$,
$B_{c}(\rho)$ and define the continuous-time LPV model as

$\left\\{\begin{aligned}
\dot{x}(t)&=A_{c}(\rho(t))x(t)+B_{c}(\rho(t))u(t)\equiv f(x(t),u(t)),\\\
\rho(t)&=\sigma(x(t)),~{}x_{0}=x(0),~{}t\geq 0,\end{aligned}\right.$ (5)

where $\rho$ is the so-called scheduling parameter that belongs to the
scheduling set $\mathcal{P}\subset\mathbb{R}^{n_{\rho}}$ that defines the
range of the scheduling parameter, and it should be compact. Further, $\rho$
depends on the states through the mapping
$\sigma\colon\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{\rho}}$.

For the nonlinear system under consideration (on one plane of the ballbot),
the state dimension is $n_{x}=4$. We consider two scheduling parameters of the
LPV model as $\theta_{y}(t)$ and $\dot{\theta}_{y}(t)$, thus $n_{\rho}=2$ and
$\rho(t)=(\theta_{y}(t),\dot{\theta}_{y}(t))$. Therefore, in this case, the
functionality $\sigma$111Functionality:
$\sigma(x)=x^{\top}\left[\begin{array}[]{cccc}0&1&0&0\\\
0&0&0&1\end{array}\right]^{\top}=\left[\begin{array}[]{cc}\theta_{y}&\dot{\theta}_{y}\\\
\end{array}\right]=\rho$. can be seen as a linear function that chooses from
the state vector only some states. Finally, the input dimension is $n_{u}=1$
with $u(t)=\tau_{y}(t)$ representing the virtual torque of the virtual wheel
Fig. 1.

The considered LPV system’s matrices of the continuous-time state-space model
are given by (6), and the elements of the matrices $A_{c}$ and $B_{c}$ are
given in Table I.

$\displaystyle\footnotesize A_{c}(\rho)=\begin{bmatrix}0&0&1&0\\\ 0&0&0&1\\\
0&A_{32}&A_{33}&A_{34}\\\
0&A_{42}&A_{43}&A_{44}\end{bmatrix},~{}B_{c}(\rho)=\begin{bmatrix}0\\\ 0\\\
B_{31}\\\ B_{41}\end{bmatrix},$ (6)

TABLE I: Elements of the operators $A_{c}(\rho)$ and $B_{c}(\rho)$, with their respective scheduling parameters given within the parentheses $A_{ij}(\rho)$ and $B_{ij}(\rho)$, and their scheduling signals are given by $\rho=(\theta_{y},\dot{\theta}_{y})$ as shown in the table. $d(\theta_{y})$ | $b_{1}b_{3}-(-b_{2}+\ell r_{b}\cos({\theta_{y}}))^{2}$
---|---
$A_{32}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $\frac{\ell g\operatorname{sinc}(\theta_{y})(b_{2}-\ell r_{b}\cos({\theta_{y}}))}{d(\theta_{y})}$
$A_{33}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $-\frac{b_{3}b_{4}}{d(\theta_{y})}$
$A_{34}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $\frac{b_{3}\ell r_{b}\sin({\theta_{y}})\dot{\theta}_{y}}{d(\theta_{y})}$
$A_{42}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $\frac{b_{1}\ell g\operatorname{sinc}(\theta_{y})}{d(\theta_{y})}$
$A_{43}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $-\frac{b_{4}(b_{2}-\ell r_{b}\cos{(\theta_{y})})}{d(\theta_{y})}$
$A_{44}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $\frac{\ell r_{b}\sin{(\theta_{y})}\dot{\theta}_{y}(b_{2}-\ell r_{b}\cos{(\theta_{y})})}{d(\theta_{y})}$
$B_{31}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $\frac{r_{b}(b_{3}-(b_{2}-\ell r_{b}\cos{(\theta_{y})})}{r_{w}d(\theta_{y})}$
$B_{41}\big{(}\theta_{y},\dot{\theta}_{y}\big{)}$ | $\frac{r_{b}((b_{2}-\ell r_{b}\cos{(\theta_{y})})-b_{1})}{r_{w}d(\theta_{y})}$

A discretization of the nonlinear continuous model is mandatory before setting
the framework that combines the LPV with Model Predictive Control (MPC).
Therefore, to discretize (5), many methods exist that differ in accuracy and
numerical properties (s.a., numerical stability). For instance, the Forward
Euler method is commonly used due to its ease of implementation, but it
usually endures low-level accuracy and, most importantly, inherits numerical
instability. Therefore, we evade such numerical pitfalls by utilizing the
Runge–Kutta $4^{\textrm{th}}$ order with a single-time step (RK4). Consider
the sampling time $t_{s}$ that spans the continuous time as
$t=t_{s}k,~{}k\in\mathbb{Z}_{+}\cup\\{0\\}$, and, it is denoted as
$x(t_{k})=x(kt_{s})=x_{k}$. We denote the angular velocity of the tilt angle
$\theta_{y}$ as $\dot{\theta}_{y}(t):=\omega_{y}(t)$ and the corresponding
discrete value is $\omega_{k}$. Similarly, we denote the angular velocity
$\dot{\phi}_{y}(t):=\varphi_{y}(t)$, thus, the complete discrete state vector
remains
$x_{k}=\left[\begin{array}[]{cccc}\phi_{k}&\theta_{k}&\varphi_{k}&\omega_{k}\end{array}\right]^{\top}$.

By denoting the continuous in-time LPV operator $f_{\text{LPV}}$ (i.e.,
$f_{\text{LPV}}$ is equivalent with the nonlinear operator $f$) of the ballbot
system, we can rewrite (5) as
$\dot{x}(t)=f_{\text{LPV}}(x(t),u(t)):=A_{c}(\sigma(x(t)))x(t)+B_{c}(\sigma(x(t)))u(t)$
with the input $u(t)=\tau_{y}(t)$. The Runge-Kutta discretization method of
order $4^{\textrm{th}}$ (RK4) with a single time step $t_{s}$ when zero-order
hold (ZOH) is considered remains:

$\text{RK4}~{}\left\\{\begin{aligned}
\kappa_{1}&=f_{\text{LPV}}\left(x_{k},u_{k}\right)\\\
\kappa_{2}&=f_{\text{LPV}}\left(x_{k}+(t_{s}/2)\kappa_{1},u_{k}\right)\\\
\kappa_{3}&=f_{\text{LPV}}\left(x_{k}+(t_{s}/2)\kappa_{2},u_{k}\right)\\\
\kappa_{4}&=f_{\text{LPV}}\left(x_{k}+t_{s}\kappa_{3},u_{k}\right)\\\
x_{k+1}&=x_{k}+(t_{s}/6)(\kappa_{1}+2\kappa_{2}+2\kappa_{3}+\kappa_{4}).\\\
\end{aligned}\right.$ (7)

Thus, we introduce the following remark to proceed with linear operations with
the LPV formulation, which will be essential for deriving a predictor that
retains the characteristics of the underlying optimization problem.

######  Remark 1 (The operator $f_{\text{LPV}}$ with a given scheduling
variable $\rho$ is linear)

When the scheduling parameter $\rho$ is given, the LPV operator remains linear
w.r.t. the state $x$ and input $u$. In such cases, the LPV continuous operator
will be denoted as:

$f_{\text{LPV}}^{(\rho)}(\rho,x(t),u(t)):=A_{c}(\rho)x(t)+B_{c}(\rho)u(t).$
(8)

Consequently, when the scheduling parameter can be considered known, we will
denote the LPV operator with a superscript $\rho$ as
$f_{\textrm{LPV}}^{(\rho)}$ that will indicate that $\rho$ is independent of
the state/input allowing linear actions on both. In particular, Remark 1
allows the LPV predictor to be inserted in a linear MPC framework and retain
all the characteristics of the optimization problem that stay invariant under
linear operations. For instance, the convexity of the optimization problem and
the linear constraints can be retained. At the same time, efficient algorithms
of QP can solve the underlying control problem efficiently.

###  II-C Stabilization and closed-loop parameters identification

To identify the parameters of the ballbot as a real physical system, it was
stabilized in [1] through a proportional integral derivative (PID) controller.
During stabilization in [1], a multi-harmonic excitation signal provided
enough measurements to identify a linear model that is defined with the
$p$-parameters provided in Table II and in the study [1]. The identified
linear model could stabilize the ballbot with an optimal feedback gain $K$
utilizing the solution of the linear quadratic regulator (LQR) approach.

TABLE II: Physical parameters of the ballbot model (2), with $p$ being the linearized model parameters from [1]. $p_{1}$ | $p_{2}$ | $p_{3}$ | $p_{4}$ | $p_{5}$
---|---|---|---|---
-342.6038 | -52.8301 | -1425.9 | -36.0734 | -9.1477
$p_{6}$ | $\ell[m]$ | $r_{b}[m]$ | $r_{w}[m]$ | $g[m/s^{2}]$
251.8476 | 0.2978 | 0.12 | 0.05 | 9.81

In this study, we will use the identified linear model from [1] for
discovering the $b$ parameters that define the refined values in (1) after
solving (offline) a nonlinear optimization problem as will be explained in
Sec. III-A. Consequently, the LPV embedding of the identified nonlinear model
will allow the LPV integration within the proposed MPC scheme as shown in Sec.
III-B that will handle the stabilization and reference tracking of the ballbot
simultaneously without cascading complex control designs based on the
linearized model.

###  II-D Model Predictive Control with LPV Embeddings

Model predictive control with linear parameter-varying embeddings implies that
the MPC optimally controls an underlying time-varying model based on the
scheduling parameters shown in the LPV formulation (5). The standard MPC form
begins with the cost function defined as $J_{k}(u)$ as seen in (9), where
$x_{i|k}$ denotes the discretized state vector at time $k$ and prediction $i$.
The $i$ varies in length $i=0,\ldots,N-1$ with $N$ the prediction horizon
length. The reference state vector is given as $x^{\mathrm{ref}}_{i|k}$.

$\displaystyle J_{k}(u)$ $\displaystyle=\sum_{i=0}^{N-1}(\lVert
x_{i|k}-x^{\mathrm{ref}}_{i|k}\rVert^{2}_{Q}+\lVert u_{i|k}\rVert^{2}_{R})+$
(9) $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\underbrace{(\lVert
x_{N|k}-x^{\mathrm{ref}}_{N|k}\rVert^{2}_{P}}_{\text{terminal cost}},$

where $Q\succeq 0,~{}R\succ 0,~{}P\succ 0$ are weighting symmetric matrices.
Theoretically, to provide closed-loop stabilization of the LPV model, for some
feedback gain $K(\rho)$ the $A(\rho)-B(\rho)K(\rho)$ should be Hurwitz for all
possible parameter variation of $\rho$. In this study, with the suitable
tuning of the LQR on the linearized model, we could provide a non-parametric
state feedback gain $K$ that satisfies the Hurwitz condition for the LPV over
a dense grid that covers the scheduling space $\mathcal{P}$. In addition, the
algebraic Riccati solution $P$ will be used to penalize the terminal cost of
the MPC.

In (9) and on the one hand, the explicit decision variables are the input
entries $u_{i|k},~{}i=0,\ldots,N-1$. Thus, the optimizer should provide a
solution in $u$ that minimizes the energy function (9). On the other hand, the
$x_{i|k},~{}i=0,\ldots,N$ consists of the set with the implicit decision
variables inserted through a predictor. Suppose the prediction $x_{i|k}$
depends nonlinearly on previous states. In that case, the energy function will
not remain quadratic or convex, inevitably leading to nonconvex optimization
formulation where the cost of solving such problems inherits high complexity.

To evade such a computational burden, by utilizing the LPV formulation
together with the property in Remark 1, when the scheduling variable is given,
the implicit decision variables $x_{i|k}$ are introduced linearly in the
energy function and retain the quadratic manifold where the minimization is
quite efficient through the quadratic program (QP) that will introduce next.

Importantly, good scheduling estimation will make the adaptive nature of the
LPV over a potentially nonlinear manifold quite accurate. Good predictions for
the scheduling variable lie at the heart of the LPV method, and methods that
introduce ways, e.g., scheduling tube control [10], can provide further
guarantees, such as recursive feasibility and stability.

Back to our aim, we want to minimize the quadratic cost function $J$ in (9)
under the linear equality constraints that satisfy the LPV (given scheduling)
predictor model together with the following linear inequalities given in
(LABEL:eq:constraints).

$\displaystyle\begin{aligned}
&\mathcal{X}=\\{x_{k}\in\mathbb{R}^{n_{{\mathrm{x}}}}~{}|~{}G^{x}x_{k}\leq
h^{x}\\},\\\ &\mathcal{U}=\\{u_{k}\in\mathbb{R}^{n_{u}}~{}|~{}G^{u}u_{k}\leq
h^{u}\\},\end{aligned}$ (10)

where, the state/input inequality constraints emerge as

$\displaystyle G^{x}$
$\displaystyle=\left[\begin{array}[]{c}+I_{n_{\mathrm{x}}}\\\
-I_{n_{\mathrm{x}}}\end{array}\right]\in\mathbb{R}^{(2n_{\mathrm{x}})\times
n_{\mathrm{x}}},~{}h^{x}=\left[\begin{array}[]{c}+x_{\text{max}}\\\
-x_{\text{min}}\end{array}\right]\in\mathbb{R}^{2n_{\mathrm{x}}},$
$\displaystyle G^{u}$ $\displaystyle=\left[\begin{array}[]{c}+I_{n_{u}}\\\
-I_{n_{u}}\end{array}\right]\in\mathbb{R}^{(2n_{u})\times
n_{u}},~{}h^{u}=\left[\begin{array}[]{c}+u_{\text{max}}\\\
-u_{\text{min}}\end{array}\right]\in\mathbb{R}^{2n_{u}}.$

The superiority of the MPC control strategy comes exactly at the point where
constraints can be introduced. MPC casts the control problem into optimization
with the major feature, among other methods, that can handle constraints. MPC
can be extended easily to include trajectory planning of the ballbot when it
moves in a structured environment with obstacles. Tangent planes can play a
significant role in keeping the optimization problem as QP, like the work in
[4].

The QP optimization problem is formulated together with the LPV predictor at
the ”present” time $k$ as the function
$\mathrm{QP}(\rho_{i|k},x_{k},x^{\mathrm{ref}}_{i+k})$ with the input
arguments; 1) the given (from previous time) scheduling signal $\rho_{i|k}$;
2) the initial conditions $x_{k}$; 3) the reference $x_{i+k}^{\textrm{ref}}$.
The Quadratic Program’s (QP) functionality is defined next in Algorithm 1.

Algorithm 1 The quadratic program (QP) based-LPVRK4

0: The time step $k\in\mathbb{Z}_{+}$, the weighted costs $Q,~{}R,~{}P$, the
triplet $(\rho_{i|k},x_{k},x^{\mathrm{ref}}_{k+i})$, for $i=0,\ldots,N-1$, the
horizon length $N\in\mathbb{Z}_{+}$ and the sampling time
$t_{s}\in\mathbb{R}_{+}$.

0: Optimal numerically design $u_{i|k},~{}i=0,\ldots,N-1$.Minimization of the
quadratic cost
function:$\underset{u_{i|k}}{min}\sum_{i=0}^{N-1}\Big{(}\lVert\hat{x}_{i|k}-x^{\mathrm{ref}}_{k+i}\rVert^{2}_{Q}+\lVert
u_{i|k}\rVert^{2}_{R}\Big{)}+\underbrace{\lVert\hat{x}_{N}-x^{\mathrm{ref}}_{k+N}\rVert_{P}^{2}}_{\text{terminal
cost}}$Subject to:

1: $\hat{x}_{0|k}=x_{0|k}=x_{k}$,

2: for i=0,…,N-1 do

3:
$\kappa_{1}=f_{\textrm{LPV}}^{(\rho_{i|k})}\left(\rho_{i|k},\hat{x}_{i|k},u_{i|k}\right)$,

4:
$\kappa_{2}=f_{\textrm{LPV}}^{(\rho_{i|k})}\left(\rho_{i|k},\hat{x}_{i|k}+(t_{s}/2)\kappa_{1},u_{i|k}\right)$,

5:
$\kappa_{3}=f_{\textrm{LPV}}^{(\rho_{i|k})}\left(\rho_{i|k},\hat{x}_{i|k}+(t_{s}/2)\kappa_{2},u_{i|k}\right)$,

6:
$\kappa_{4}=f_{\textrm{LPV}}^{(\rho_{i|k})}\left(\rho_{i|k},\hat{x}_{i|k}+t_{s}\kappa_{3},u_{i|k}\right)$,

7: Update:
$\hat{x}_{i+1|k}=\hat{x}_{i|k}+(t_{s}/6)(\kappa_{1}+2\kappa_{2}+2\kappa_{3}+\kappa_{4})$,

8: end for

9: Satisfy: $\hat{x}_{i|k}\in\mathcal{X},\forall i=1,\dots,N$,

10: Satisfy: $u_{i|k}\in\mathcal{U},\forall i=0,1,\dots,N-1$.

Next, we discuss Algorithm (1) step by step when it is called at the
simulation time instance $k\in\mathbb{Z}_{+}$:

  * •

At time $k$, the input to the QP consists of the matrices $Q,~{}R,~{}P$, the
initial condition $x_{k}$, the reference $x_{k+i}^{\textrm{ref}}$, and the
scheduling variable that is either available from previous time $k-1$ or
initialized from initial conditions $x_{k}$. In particular, when the
simulation starts at time $k=0$, and there is no available scheduling
prediction from the past, the scheduling variable is initialized as
$\rho_{i|0}=\sigma(x_{0}),~{}i=0,\ldots,N-1$.

  * •

The optimization problem consists of the quadratic cost function (9) that is
updated with the implicit decision variables $\hat{x}_{i|k}$ as a function of
the explicit decision variables $u_{i|k}$ with the LPV predictor. This can be
done with linear operations since the scheduling variable
$\rho_{i|k},~{}i=0,\ldots,N-1$ is considered known, and the linear property
explained in Remark 1. The RK4 in lines (2-7) is the discretization scheme
that predicts the evolution.

  * •

The QP can be solved together with the constraints in lines (9-10),
efficiently from [11] and provide the optimized input design at each time $k$
numerically as $u_{i|k},~{}i=0,\ldots,N-1$.

Algorithm 2 The Ballbot QP-based LPVMPC algorithm

0: Initial conditions $x_{0}$, with reference
$x^{\textrm{ref}}_{k},k\in\mathbb{Z}_{+}\cup\\{0\\}$

0: The control input $u_{k},~{}k=0,1,\ldots$, that drives the nonlinear system
to the reference by satisfying constraints.Algorithm steps

1: Initialize the scheduling parameter for $k=0$ as next:
$\rho_{i|0}=(\theta_{0},\omega_{0}),~{}i=0,\ldots,N$

2: while $k=0,1,\ldots$ do

3: Solve the QP in Algorithm 1
$u_{i|k}\leftarrow\texttt{QP}(\rho_{i|k},x_{k},x_{i+k}^{\text{ref}}),~{}i=0,\ldots,N-1\vspace{-3mm}$

4: Update the scheduling parameters from the designed input
$u_{i|k},~{}i=0,\ldots,N-1$ with (7) as
$x_{i+1|k}=\left[\begin{array}[]{c}\phi_{i+1|k}\\\ {\theta_{i+1|k}}\\\
\varphi_{i+1|k}\\\
{\omega_{i+1|k}}\end{array}\right]\underset{\text{RK4}}{\overset{\eqref{eq:dlpv}}{=}}f_{\textrm{LPV}}\left(\left[\begin{array}[]{c}\phi_{i|k}\\\
{\theta_{i|k}}\\\ \varphi_{i|k}\\\
{\omega_{i|k}}\end{array}\right],u_{i|k}\right)$

5: Update $\rho_{i|k}=\left(\theta_{i|k},\omega_{i|k}\right),~{}i=0,\ldots,N$

6: Apply $u_{k}=u_{0|k}$ to the continuous system (5)

7: $k\leftarrow k+1$

8: Measure $x_{k}$

9: Update $\rho_{i|k}\leftarrow\rho_{i+1|k-1},~{}i=0,\ldots,N-1$

10: end while

As we have introduced in detail the QP minimization problem that is tailored
with the LPV formulation in Algorithm 1, we are ready to proceed by stating
the method in Algorithm 2 that is connected with our application on
stabilizing and tracking a given reference for the ballbot. Next, every step
in the resulting Algorithm 2 is explained.

  * •

The input to Algorithm 2 consists of the initial conditions $x_{0}$ and the
full state reference $x^{\textrm{ref}}$. For the ballbot system, the state
reference contains the angular displacement of the ball $\phi_{k}$. In
particular, $\theta_{k}$ is the tilt angle, and its reference target is zero,
thus enforcing stabilization in that operating (unstable equilibrium) point.
Finally, for the rest of the states in $x^{\text{ref}}$ that concern the
angular velocities, again, we enforce them to zero as the aim is to reach the
reference $\phi^{\textrm{ref}}$ and then the ballbot to remain there with
$\theta=\dot{\theta}=\dot{\phi}=0$ (steady state) and $\phi$ at the target
point.

  * •

In Step 1, initialization at any instant $k$ of the scheduling parameters can
be done either from the previous time as
$\rho_{i|k}=\sigma(x_{i+1|k-1}),~{}i=0,\ldots,N-1$ or from the initial
conditions as explained in Algorithm 1.

  * •

In Steps 2-3, we can call the QP Algorithm 1 to provide the optimally designed
input $u_{i|k},~{}i=0,\ldots,N-1$. Note that by having $N$ inputs, we can
compute $N+1$ state estimation through the LPV predictor.

  * •

In Step 4, given $x_{k}$ and $u_{i|k}$, we can compute the true evolution
response from the nonlinear system after using the exact/embedded nonlinear
discrete predictor (7). The scheduling parameters to be used in Step 5 can be
revealed in Step 4. The hat notation is omitted as these computations result
in the system’s true response, assuming that no other disturbance or
measurement noise will affect the system.

  * •

In Step 6, the one-step implementation of the input to the nonlinear system
(5) is done by simulating with ode45 (Runge Kutta 45 with adaptive time step)
where between two consequent time instances that differ only one sampling time
$t_{s}$, zero-order hold (ZOH) (constant input of the sampling time) is
considered. This high-fidelity simulation ensures good agreement with the true
response of the physical plant.

  * •

In Step 7, we slide the prediction horizon by one (sliding the prediction
window with length $N$ by one sampling time $t_{s}$)

  * •

In Step 8, we measure the new initial conditions $x_{k}$.

  * •

Finally, in Step 9, we update the new scheduling variable at time $k$ from the
previous simulation time $k-1$, and the while loop in Steps 2-10 runs to
provide the next control inputs $u_{k}$ for arbitrary $k$.

##  III Results and Discussion

###  III-A Refinement of the ballbot physical parameters

A linearized continuous model at $x_{e}=[0,0,0,0]^{T}$ can be derived through
the LPV evaluated at $x_{e}$ as in (11).

$\displaystyle A_{c}(0,0)$ $\displaystyle=\begin{bmatrix}0&0&1&0\\\ 0&0&0&1\\\
0&\underbrace{A_{32}(0,0)}_{p_{1}}&\underbrace{A_{33}(0,0)}_{p_{2}}&\underbrace{A_{34}(0,0)}_{0}\\\
0&\underbrace{A_{42}(0,0)}_{p_{4}}&\underbrace{A_{43}(0,0)}_{p_{5}}&\underbrace{A_{44}(0,0)}_{0}\end{bmatrix},~{}$
(11) $\displaystyle B_{c}(0,0)$ $\displaystyle=\begin{bmatrix}0\\\ 0\\\
\underbrace{B_{3}(0,0)}_{p_{3}}\\\
\underbrace{B_{4}(0,0)}_{p_{6}}\end{bmatrix},~{}d(0)=b_{1}b_{3}-(-b_{2}+\ell
r_{b})^{2}.$

By introducing the functions
$h_{i}(b):\mathbb{R}^{4}\rightarrow\mathbb{R},~{}i=1,~{}\ldots,6$ with suport
the parameterized vector $b=(b_{1},b_{2},b_{3},b_{4})$, we can define:

$\displaystyle A_{32}(0,0)$
$\displaystyle:=h_{1}(b)=\frac{gl\left(b_{2}-l\mathrm{r_{b}}\right)}{d(0)},$
(12) $\displaystyle A_{33}(0,0)$
$\displaystyle:=h_{2}(b)=-\frac{b_{3}\,b_{4}}{d(0)},$ $\displaystyle
A_{42}(0,0)$ $\displaystyle:=h_{4}(b)=\frac{b_{1}\,g\,l}{d(0)},$
$\displaystyle A_{43}(0,0)$
$\displaystyle:=h_{5}(b)=-\frac{b_{4}\,\left(b_{2}-l\,\mathrm{r_{b}}\right)}{d(0)},$
$\displaystyle A_{34}(0,0)$ $\displaystyle:=0,~{}A_{44}(0,0):=0,$
$\displaystyle B_{3}(0,0)$
$\displaystyle:=h_{3}(b)=\frac{b_{3}\,\mathrm{r_{b}}-\mathrm{r_{b}}\,\left(b_{2}-l\,\mathrm{r_{b}}\right)}{\mathrm{r_{w}}\,d(0)},$
$\displaystyle B_{4}(0,0)$
$\displaystyle:=h_{6}(b)=\frac{\mathrm{r_{b}}\,\left(b_{2}-l\,\mathrm{r_{b}}\right)-b_{1}\mathrm{r_{b}}}{\mathrm{r_{w}}\,d(0)}.$

To infer the nonlinear $b$-parameters, given the $p$-parameters values, we
must solve the nonlinear system of equations described in the set of equations
(12) for the unknown vector $b=(b_{1},b_{2},b_{3},b_{4})$ and with the
functionality $h_{i}(b):\mathbb{R}^{4}\rightarrow\mathbb{R},~{}i=1,\ldots,6$.
The problem remains to find $b^{*}$ that satisfies $p_{i}-h_{i}(b^{*})\approx
0$ for all $i=1,\ldots,6$. Thus, we define the operator
$F(b):\mathbb{R}^{4}\rightarrow\mathbb{R}^{6}$ along with the Jacobian next:

$\displaystyle F(b):=\left[\begin{array}[]{cc}p_{1}-h_{1}(b)\\\\[2.84526pt]
p_{2}-h_{2}(b)\\\\[2.84526pt] \vdots\\\\[2.84526pt] p_{6}-h_{6}(b)\\\
\end{array}\right],~{}\nabla F(b)=\left[\begin{array}[]{ccc}\frac{\partial
F_{1}}{\partial b_{1}}&\cdots&\frac{\partial F_{1}}{\partial
b_{4}}\\\\[2.84526pt] \frac{\partial F_{2}}{\partial
b_{1}}&\cdots&\frac{\partial F_{2}}{\partial b_{4}}\\\\[2.84526pt]
\vdots\\\\[2.84526pt] \frac{\partial F_{6}}{\partial
b_{1}}&\cdots&\frac{\partial F_{6}}{\partial b_{4}}\\\ \end{array}\right].$

The Newton scheme for updating the parameters $b$ with initialization $b_{0}$
and for $k\in\mathbb{Z}_{+}\cup\\{0\\}$ is

$b_{k+1}=b_{k}-\left[\nabla F(b_{k})\right]^{-1}F(b_{k}),~{}k=0,\ldots$ (13)

Upon convergence of (13), the optimal solution vector $b^{*}$ can be obtained
to satisfy $F(b^{*})\approx 0$ (with the smallest possible residual). In such
problems, where the nonlinear parameterized manifold is nonconvex, the
iterative scheme in (13) remains divergent when random initialization is
considered, or a solution can be algebraically correct without explaining the
physical law. To tackle such a problem, initialization of the $b$-parameters
should be done by respecting the engineering regularity of the problem. Thus,
some parameters $b_{i}$ can be initialized with physical meaning, e.g.,
$b_{1},~{}b_{2},~{}b_{3}$ being positive to explain the high-fidelity
knowledge of the physical plant [1] where $b_{4}$ as concerns friction can
remain arbitrary but also with a reasonable small value.

In Fig. 2, the iterative scheme (13) converged after five steps and gave the
optimal solution with a residual error that stagnated to the value $0.4330$.
This error cannot be improved further as the analysis relies on the identified
parameters $p$ from the study [1], where the dynamics are corrupted with noise
and/or slight nonlinear behavior from the physical plant.

Continuously, in Fig. 3, the simulation indicates the expected performance
after substituting the obtained parameters $b$ in Table III and applying a
multi-harmonic input to the stabilized plant. In particular, in Fig. 3, the
nonlinear model compared to the linearized model stays close to the dynamical
evolution for all states when the dynamics are close to the linearization
operation point $x_{e}=0$ as expected. The discrepancy increases for larger
deviations in $\phi$ from the origin where the linearized cannot be trusted.
Finally, in Fig. 3, the comparison between the LPV and the nonlinear model
certifies that the LPV is an equivalent embedding of the nonlinear system
where differences are buried in machine precision error. Finally, the
discrete-time implementation with the RK4 in (7) with a ZOH also certifies a
good accuracy compared to the ode45.

Figure 2: Convergence of the Newton scheme with residual error $\|F(b^{*})\|=0.4330$. The optimal parameter vector $b^{*}$ is shown in Tab. III. TABLE III: The refined $b$ parameters of the nonlinear model. $b^{*}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$
---|---|---|---|---
Optimal | 0.002483 | 0.059325 | 0.143093 | -0.07436
Figure 3: Models simulations; the linearized model given in dashed blue line;
non-linear model simulation given in black; LPV in dashed red lines; all done
using ode45 on MATLAB. Finally, with green is the RK4 method utilizing (7).
The input $\tau_{y}$ is a multiharmonic signal.

###  III-B LPVMPC for Reference Tracking

Reference tracking for a ballbot refers to the stabilized motion by the LPVMPC
method that can maneuver over any given state reference $x^{\mathrm{ref}}$. We
will consider three study scenarios. The 1st scenario will be a 1D reference
tracking of a nonsmooth discontinuous trajectory that will ”shock” the
optimization problem due to the sudden changes. The 2nd scenario will be a 2D
smooth reference as Lissajous curves, assuming the ballbot can be driven
independently in the $xz$ and $yz$ planes. Finally, a 3rd scenario will
illustrate the reference tracking of a single set point with a scheme of the
LPVMPC that can guarantee the stability of the closed-loop system and the
recursive feasibility of its optimization problem. These properties are
explained below. Tuning of the prediction horizon depends on the physical
system’s operational bandwidth and can vary significantly across applications.
The ballbot is a robotic system that operates around $1$(Hz); thus, a horizon
of about $1$ (s) is reasonable.

####  III-B1 1st Scenario: Nonsmooth 1D reference

In Fig. 4, we consider the following 2-set point nonsmooth reference. The
ballbot should start from the origin and at $t=1$ s should roll for one circle
$2\pi$ rad, then stay there for $2$ s and for the remaining $1$ s to return to
the origin; thus, it will roll a total distance of approximately $(2\cdot
2\pi\cdot r_{b}\approx 1.5~{}\text{m})$ within $4$ s.

At the beginning and near the origin, positive virtual torque applies on the
ballbot to decrease the tilt angle $\theta$ that instantly will move the
ballbot to the negative $\phi$, where immediately switching direction (normal
behavior of such non-minimal phase system, similar to a bike where to turn,
first a small maneuver to the opposite direction is needed). Continuously the
ballbot increases speed in the $\phi$ positive direction while approaching the
target $2\pi$ with satisfying accuracy and faster than the linear one after
$\approx 1.3$ s, then stays there and at around $2^{+}$ s where the prediction
horizon of time length $N\cdot t_{s}=20\cdot 0.05=1$ s receives the
information of the sudden change in reference from $2\pi$ to $0$ rad, the
ballbot slightly overshoots in $\phi$ so to change the tilt angle again and
starts to accelerate in the opposite direction. After less than a second, the
ballbot has reached its origin without overshooting, which outlines the good
performance seen mainly in nonlinear MPC frameworks.

In Fig. 4, we compare the MPC performances between the LPV and the linearized
model. As expected, the LPVMPC is generally faster at reaching the reference
with almost similar computational complexity as the linear MPC and the
potential to handle better strong maneuvering as the adaptive scheduling
variable $\rho_{i|k}$ provides linearizations over a trajectory of the
nonlinear ballbot model instead of a single fixed point for the linearized
model. The quadratic costs for the 1st scenario are defined as
$Q=\texttt{diag}([200,1,0.1,0.1])$, $R=1000$, and the terminal cost $P$ is the
Riccati solution from the LQR on the linearized model from [1] with the same
tuning of $Q$ and $R$ that will penalize the ”tail” of the horizon. The
sampling time is $t_{s}=0.05$ (s), with horizon length $N=20$, and the MPC
solution is obtained under the input/state hard constraints introduced in
Table IV and highlighted with the gray background in Fig. 4.

TABLE IV: MPC Constraints Variable | Lower bound | Upper bound | Units
---|---|---|---
$\phi_{k}$ | $-\infty$ | $+\infty$ | [rad]
$\theta_{k}$ | $-\pi/3$ | $+\pi/3$ | [rad]
$\dot{\phi}_{k}$ | $-10\pi$ | $+10\pi$ | $[\frac{\rm rad}{\rm s}]$
$\dot{\theta}_{k}$ | $-2\pi$ | $+2\pi$ | $[\frac{\rm rad}{\rm s}]$
$u$ | $-1.5$ | $+1.5$ | [Nm]

In the 1D reference tracking, as shown in Fig. 4, traveling in $y$ direction,
the ballbot system also balances itself simultaneously, as shown with $\theta$
being between ranges of $0\,rad$ to $0.5\,rad$ while traveling a distance of
$\phi$, which is referenced for $2\pi\,rad$, and angular velocities of both
are shown in $\dot{\theta}$, and $\dot{\phi}$ respectively. They show the
agile and robust movements of the ballbot. Also, in correspondence to the
reference and torque, $\tau_{y}$ is applied to the virtual wheel actuator. The
virtual torque is shown as discrete steps, i.e., ZOH, and the energy function
$J$ to be minimized to get the optimal control input $\tau_{y}$.

Figure 4: Ballbot reference tracking with the angular displacement $\phi$
traveling from the origin to $2\pi$ rad and back. A comparison between the
linear and LPV MPC frameworks is illustrated. $\theta$ measuring balancing of
ballbot, $\dot{\phi}$, $\dot{\theta}$ measuring the angular velocities,
$\tau_{y}$ is the virtual torque, $J$ is the energy cost function.

####  III-B2 2nd Scenario: Lissajous 2D smooth curves

Solving the LPVMPC reference tracking problem in parallel for the $yz$-plane
and $xz$-plane, we can drive the ballbot over the $xy$ plane with reference as
Lissajous curves in Fig. 5. The ballbot starts from
$(x,y)=(\phi_{x}r_{b},\phi_{y}r_{b})=(0,0)$, and following the path to global
coordinates as: $X^{\textrm{ref}}=\phi_{x}r_{b}=0.12\cdot 2\pi\sin(0.3t)$ and
$Y^{\textrm{ref}}=\phi_{y}r_{b}=0.12\cdot 2\pi\sin(0.4t)$, under the
constraints in Table IV. In that case, the reference is smoother without
sudden changes. Therefore, we can increase the weight for matching the $\phi$
without inheriting infeasibility problems; thus, the quadratic cost for the 2D
reference tracking in each direction can be set as
$Q=\texttt{diag}([1000,1,0.1,0.1])$. All the other quadratic weights and
parameter specifications are considered the same as in the 1st scenario. In
Fig. 5, the LPVMPC control strategy accurately drives the ballbot to the
reference within an area of around $2~{}m^{2}$ and for $70$ s.

Figure 5: Stabilization and Lissajous curve trajectory tracking of the ballbot
in 2D using the LPVMPC algorithm solved in both $xz$ and $yz$ planes.

####  III-B3 3rd Scenario: Single set point reference and LPVMPC with
guarantees

In our last scenario, regulated as in the 1st (e.g., quadratic costs, etc.),
we want to illustrate the crucial advantage of considering LPV models for MPC
tasks, as it can provide theoretical guarantees. On the one hand, when a
linearized model is considered instead of the underlying nonlinear, achieving
theoretical guarantees within MPC could be a straightforward task. Still, it
will remain quite conservative and thus not practically useful. On the other
hand, in the case of NMPC, guarantees can also be provided at the expense of
computational complexity. Consequently, an excellent alternative for providing
theoretical guarantees is the LPV formulation within MPC, which can be
considered an excellent trade-off between the complexity and conservativeness
of the previous two general approaches.

To prove theoretical guarantees, in Algorithm 1, we have to introduce
additionally the terminal equality constrain as

$\mathcal{X}_{t}=\\{\hat{x}_{N|k}\in\mathbb{R}^{n_{\textrm{x}}}~{}|~{}\hat{x}_{N|k}=x^{\text{ref}}_{N+k}\\}.$
(14)

Thus, the state constraints $\mathcal{X}$ in Algorithm 1 are enhanced with
$\hat{x}_{N|k}\in\mathcal{X}_{t}$. Note that in the presence of the terminal
constraint $\mathcal{X}_{t}$, the terminal cost remains zero; thus, in that
case, no use of the weight $P$ is needed. In Fig. 6, the dynamics are
illustrated when considering the terminal constraint (cyan line) in (14),
where we lose some performance compared with no guarantees (red line). Thus,
with active (14) and making use of the Theorem 1 in [12] (page: 6), all the
conditions are satisfied in our 3rd scenario, which leads to the following
conclusions; i) the LPVMPC optimization problem is recursively feasible; ii)
the closed-loop system satisfies the constraints $u_{i|k}\in\mathcal{U}$,
$\hat{x}_{i|k}\in\mathcal{X}$, and $\hat{x}_{N|k}\in\mathcal{X}_{t}$; iii)
$(u^{\text{ref}},x^{\text{ref}})=(0,\pi)$ is an exponentially stable (forced)
equilibrium of the closed-loop system. The proof of these properties follows
the same reasoning as the proof of Theorem 1 in [12].

Figure 6: Single-set point reference i.e., $\pi$, for the LPVMPC framework with and without theoretical guarantees. TABLE V: Comparison of computation times for a single QP using YALMIP with LPVMPC. The simulations are performed on a Dell Latitude 5590 laptop with an Intel(R) Core(TM) i7-8650U CPU and 16 GB of RAM. The scenarios are implemented in MATLAB, utilizing the YALMIP toolbox [11], with an optimality tolerance of $10^{-8}$. LPVMPC with $t_{s}=0.05$ s | Mean | Standard deviation
---|---|---
1st Scenario | 0.0086 s | 0.0104 s
2nd Scenario | 0.0087 s | 0.0132 s
3rd Scenario | 0.0058 s | 0.0025 s

##  IV Conclusion

In this study, we were concerned with the reference tracking of a ballbot by
utilizing the LPV embedding within the MPC framework. An important step was
the refinement of the parameters of the ballbot system from measurements of a
past study. Thus, having the nonlinear model allowed the nonlinear embedding
of the LPV form. An improved discretization method based on the well-known
Runge-Kutta $4^{\textrm{th}}$ order introduced in the LPVMPC framework
potentially improves the accuracy issue that can arise in the forward Euler
method. The results indicated that applying the LPVMPC-based control method
for simultaneously balancing and reference tracking can achieve real-time
implementation as the average timing of solving the QP problem is usually
$(10)$-times less than the sampling time under consideration.

Utilizing such a control design that preserves the nonlinear behavior of the
model is advantageous and practical, unlike methods that work with linearized
models or try to solve very complex optimization problems tailored directly to
the nonlinear, which might be impossible online. The method is illustrated
with practical reference tracking scenarios where theoretical guarantees can
be provided in the single-set point case. Proving stability and recursive
feasibility in the case of multiple reference points is left to future
investigation. Finally, as the theoretical establishment of such a robotic
system has matured, implementing this method in the physical ballbot in a real
environment is left for our immediate future endeavors.

## Acknowledgement

We thank Mr. Ing.-Ievgen Zhavzharov, who constructed the ballbot in Fig. 1.

## References

  * [1] M. Studt, I. Zhavzharov, and H. S. Abbas, “Parameter identification and LQR/MPC balancing control of a ballbot,” in _2022 European Control Conference (ECC)_ , 2022, pp. 1315–1321.
  * [2] P. Fankhauser and C. Gwerder, “Modeling and control of a ballbot,” 2010. [Online]. Available: https://api.semanticscholar.org/CorpusID:58555644
  * [3] D. B. Pham, H. Kim, J. Kim, and S.-G. Lee, “Balancing and transferring control of a ball segway using a double-loop approach [applications of control],” _IEEE Control Systems Magazine_ , vol. 38, no. 2, pp. 15–37, 2018.
  * [4] M. Nezami, D. S. Karachalios, G. Schildbach, and H. S. Abbas, “On the design of nonlinear MPC and LPVMPC for obstacle avoidance in autonomous driving*,” in _2023 9th International Conference on Control, Decision and Information Technologies (CoDIT)_ , 2023, pp. 1–6.
  * [5] B. Siciliano, L. Sciavicco, L. Villani, and G. Oriolo, _Robotics: Modelling, Planning and Control_ , 1st ed. Springer Publishing Company, Incorporated, 2008.
  * [6] S. Puychaison and B. S., “Mouse type ballbot identification and control using a convex-concave optimization,” in _In Proceedings of the International Automatic Control Conference (CACS 2019)_ , National Taiwan Ocean University, Keelung,Taiwan, Nov. 2019, pp. 1–6.
  * [7] P. S. G. Cisneros, S. Voss, and H. Werner, “Efficient nonlinear model predictive control via quasi-LPV representation,” in _2016 IEEE 55th Conference on Decision and Control (CDC)_ , 2016, pp. 3216–3221.
  * [8] M. M. Morato, J. E. Normey-Rico, and O. Sename, “Model predictive control design for linear parameter varying systems: A survey,” _Annual Reviews in Control_ , vol. 49, pp. 64–80, 2020.
  * [9] H. S. Abbas, R. Tóth, N. Meskin, J. Mohammadpour, and J. Hanema, “A robust MPC for input-output LPV models,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 12, pp. 4183–4188, 2016.
  * [10] H. S. Abbas, “Linear parameter-varying model predictive control for nonlinear systems using general polytopic tubes,” _Automatica_ , vol. 160, p. 111432, 2024. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S000510982300599X
  * [11] J. Löfberg, “YALMIP : A toolbox for modeling and optimization in matlab,” in _In Proceedings of the CACSD Conference_ , Taipei, Taiwan, 2004.
  * [12] C. Verhoek, J. Berberich, S. Haesaert, R. Tóth, and H. S. Abbas, “A linear parameter-varying approach to data predictive control,” 2023. [Online]. Available: https://arxiv.org/abs/2311.07140


# Realization of the all-optical phase modulator, filter, splitter, and self-
consistent logic gates based on assembled magneto-optical heterostructures

Jie Xu1,2, Yun You3, Fengwen Kang4, Sanshui Xiao5, Lujun Hong6, Yun Shen6,
Yamei Luo1,2, and Kosmas L. Tsakmakidis7 1School of Medical Information and
Engineering, Southwest Medical University, Luzhou 646000, China  2Medical
Engineering & Medical Informatics Integration and Transformational Medicine of
Luzhou Key Laboratory, Luzhou 646000, China  3School of Science, East China
Jiaotong University, Nanchang 330013, China  4College of Materials Science and
Engineering, Sichuan University, Chengdu 610065, China  5DTU Fotonik,
Department of Photonics Engineering, Technical University of Denmark, DK-2800
Kgs. Lyngby, Denmark  6Institute of Space Science and Technology, Nanchang
University, Nanchang 330031, China  7Section of Condensed Matter Physics,
Department of Physics, National and Kapodistrian University of Athens,
Panepistimioupolis, GR-157 84 Athens, Greece  Corresponding authors: J. Xu
(xujie011451@163.com), Y. Luo (luoluoeryan@126.com), and K. L. Tsakmakidis
<EMAIL_ADDRESS>
###### Abstract

All-optical computing has recently emerged as a vibrant research field in
response to the energy crisis and the growing demand for information
processing. However, the efficiency of subwavelength-scale all-optical devices
remains relatively low due to challenges such as back-scattering reflections
and strict surface roughness. Furthermore, achieving multifunctionality
through the reassembly of all-optical structures has thus far been rarely
accomplished. One promising approach to address these issues is the
utilization of one-way edge modes. In this study, we propose four types of
deep-subwavelength ($\sim 10^{-2}\lambda_{0}$, where $\lambda_{0}$ is the
wavelength in vacuum) all-optical functional devices: a phase modulator, a
filter, a splitter, and logic gates. These devices are based on robust one-way
modes but do not require an external magnetic field, which can allow for
flexible assembly. In particular, we investigate a phase modulation range
spanning from $-\pi$ to $\pi$, a perfect filter that divides the input port’s
one-way region into two output one-way regions with equal bandwidth, a multi-
frequency splitter with an equal splitting ratio (e.g., 50/50), and self-
consistent logic gates. We validate these theoretical findings through
comprehensive full-wave numerical simulations. Our findings may find
applications in minimal optical calculations and integrated optical circuits.

##  I Introduction

Optical communication is known for its low loss [1], parallel calculation
capability [2], and high speed. As a result, all-optical devices are
attracting increasing attention. Regular all-optical devices include, but are
not limited to, optical filters, sensors, splitters, switches,
frequency/phase/amplitude modulators, couplers, and amplifiers. Most of these
devices rely on the electro-optic (EO) effect, thermo-optic (TO) effect, or
nonlinear optical (NLO) effect. Among these optical effects, the linear EO
effect has been extensively studied in fields such as optical sensing and
optical frequency combs [3] due to its ultrafast response [4]. The TO effect
is commonly employed in the design of optical switches [5], while the NLO
effect is primarily used for light control with light, enabling applications
such as laser generation and optical switching [6]. However, all-optical
devices based on EO, TO, or NLO suffer from high power consumption due to
reflection induced by reciprocity/symmetry and fabrication errors.

In recent years, there has been a significant focus on one-way edge modes due
to their potential for breaking the diffraction limit (similar to SPPs) [7, 8,
9], robust processing/computing capabilities [10], and low loss. One common
approach to achieve one-way edge modes is by utilizing photonic crystals
(PhCs). The key lies in engineering a special surface that is bounded by
materials with zero and nonzero Chern numbers. According to the bulk-edge
correspondence [11, 12], this ensures the existence of one-way edge mode(s) on
this surface. The robustness of one-way modes in PhCs to imperfections or
bends has been proposed theoretically and verified experimentally by several
groups [13, 14, 15]. However, PhCs still pose challenges in terms of their
relatively complex manufacturing process and the large device size compared to
the wavelength of the guiding wave ($\lambda_{0}$). This poses a challenge for
achieving subwavelength or even deep/ultra-subwavelength optical devices.

Recently, our group has proposed several subwavelength magneto-optical (MO)
one-way structures and discovered interesting phenomena/devices, such as
unidirectional/bidirectional slow light [16, 17], a perfect optical buffer
with zero phase shift [18], and all-optical logic gates (LGs) [19]. The most
significant advantage of the MO one-way structure is its extremely simple
manufacturing process, requiring only two joint layers and not having strict
requirements for surface roughness. Ultra-subwavelength all-optical devices
can be achieved based on the MO structure, regardless of the nonlocal effect,
which holds true in common scenarios, especially when the wavenumber is
relatively small (e.g., $k<100k_{0}$) [20, 21].

Here, we propose the design of assembled MO structures and several interesting
deep-subwavelength ($\sim 10^{-2}\lambda_{0}$) all-optical functional devices,
including phase modulators, filters, splitters, and self-consistent logic
gates (LGs). We study the propagation characteristics of surface
magnetoplasmons (SMPs) [22] in theory, and demonstrate the low-loss properties
of one-way SMPs through full-wave numerical simulations. Furthermore, we
verify the functionality of these all-optical devices. The key to realizing
the all-optical phase modulator, splitter, and self-consistent LGs lies in the
utilization of one-way index-near-zero (INZ) modes. Unlike other INZ modes
found in epsilon-near-zero (ENZ) [23, 24, 25], mu-near-zero (MNZ) [26], or
epsilon-and-mu-near-zero (EMNZ) [27] metamaterials, the INZ modes in our work
exhibit a single propagation direction and can be easily adjusted within a
continuous band.

##  II Tunable INZ modes and all-optical phase modulator

Figure 1: Assembled all-optical communication systems and functional devices
utilizing YIG with remanence, along with PEC and PMC walls.

Yttrium iron garnet (YIG) stands out as one of the most captivating MO
materials, finding widespread applications in optical communication. These
applications encompass optical isolators [28, 29], band-pass filters [30],
high-efficiency optomagnonic micron-sized resonators [31], robust LGs [19],
and more. Benefiting from the unique magnetization characteristics of YIG,
residual magnetism persists even after removing the bias magnetic field. This
property holds great potential for the realization of reliable all-optical
communication, as depicted in Fig. 1. In this paper, we demonstrate that by
utilizing YIG with remanence and PMC walls, a wide range of all-optical
functional devices, including phase modulators, filters, splitters, and LGs,
can be achieved. Importantly, these devices can be finely controlled and
integrated into a system, as illustrated in Fig. 1.

To investigate the impact of PMC walls on SMPs, we first construct four MO
heterostructures surrounded by PEC and/or PMC walls [32]. These structures are
named as follows: the PEC-air-YIG-PEC (EDYE) structure, the PEC-air-YIG-PMC
(EDYM) structure, the PMC-air-YIG-PEC (MDYE) structure, and the PMC-air-YIG-
PMC (MDYM) structure, as illustrated in the inset pictures in Fig. 2. It
should be noted that the YIG material with remanence are used in this study,
which means that all the proposed functional devices in this work do not
require precise control of an external magnetic field. Our theoretical
analysis reveals the following forms for the dispersion relationship of these
structures:

$\alpha_{r}X_{1}+\frac{\mu_{r2}}{\mu_{r1}}k-\alpha_{d}\mu_{vr}X_{2}=0$ (1)

where $X_{1}$ and $X_{2}$ are parameters that directly depend on the lower and
upper boundary conditions of the structures, respectively. For the four
structures, their values are as follows:

$\mathrm{Lower\ boundary}\begin{aligned}
&\begin{cases}&X_{1}=\frac{1}{tanh(\alpha_{r}d_{2})}\qquad(\mathrm{PEC})\\\
&X_{1}=\frac{\alpha_{r}\cdot
tanh(\alpha_{r}d_{2})-\frac{\mu_{r2}}{\mu_{r1}}k}{\alpha_{r}-\frac{\mu_{r2}}{\mu_{r1}}k\cdot
tanh(\alpha_{r}d_{2})}\qquad(\mathrm{PMC})\end{cases}\\\ \end{aligned}$ (2)
$\mathrm{Upper\ boundary}\begin{aligned}
&\begin{cases}&X_{2}=-\frac{1}{tanh(\alpha_{d}d_{1})}\qquad(\mathrm{PEC})\\\
&X_{2}=-tanh(\alpha_{d}d_{1})\qquad(\mathrm{PMC})\end{cases}\\\ \end{aligned}$
(3)

where $\mu_{r1}=1+i\frac{\nu\omega_{r}}{\omega(1+\nu^{2})}$ and
$\mu_{r2}=-\frac{\omega_{r}}{\omega(1+\nu^{2})}$ are the nondiagonal elements
of the tensor relative permeability of the YIG layers [33, 7]. $\alpha_{r}$
and $\alpha_{d}$ represent the attenuation factors of the surface wave in the
YIG layer and air, respectively [34]. $d_{1}$ and $d_{2}$ are the thicknesses
of the air and YIG layers. $\omega_{r}$ is the characteristic circular
frequency, which is denoted as $\omega_{m}$ when the YIG reaches saturation
magnetization. In this paper, we assume $\omega_{m}=2\pi\times 5\times 10^{9}$
rad/s and $\omega_{r}=2\pi\times 3.587\times 10^{9}$ rad/s. It can be observed
that as $k\to\pm\infty$, $X_{1}\to 1$ and $X_{2}\to-1$. Additionally, Eq. (1)
tends to

$\displaystyle\begin{cases}&1+\frac{\mu_{r2}}{\mu_{r1}}+\mu_{vr}=0\iff\omega_{sp}^{+}=\omega_{r}\qquad(k\to+\infty)\\\
&-1+\frac{\mu_{r2}}{\mu_{r1}}-\mu_{vr}=0\iff\omega_{sp}^{-}=0.5\omega_{r}\qquad(k\to-\infty)\end{cases}$
(4)

According to Equation (4), we can expect, at most, one asymptotic frequency
for $k>0$, denoted as $\omega_{sp}^{+}$, and one asymptotic frequency for
$k<0$, denoted as $\omega_{sp}^{-}$. Figure 2 illustrates the four dispersion
curves of SMPs in these structures as $d_{1}=d_{2}=0.03\lambda_{m}$.
Interestingly, when the lower PEC wall in the original EDYE structure is
replaced with a PMC wall (see Fig. 2(b)), the one-way region abruptly changes
from the higher region, $\omega_{sp}^{-}<\omega<\omega_{sp}^{+}$ (see Fig.
2(a)), to the lower region, $0<\omega<\omega_{sp}^{-}$, and simultaneously the
propagation direction changes oppositely. Similarly, when the upper PEC wall
in the EDYE structure is replaced with a PMC wall (Fig. 2(c)), the asymptotic
frequencies and the one-way region remain the same, but the one-way SMPs
consistently have positive values of k, i.e., $k>0$. Furthermore, when both
the upper and lower PEC walls are replaced by PMC walls (Fig. 2(d)), the one-
way region is completely disrupted due to the emergence of modes based on
total internal reflection (TIR).

In the inset of Fig. 2, we present simulation results with a setting of
$\nu=10^{-3}\omega_{m}$, which shows good agreement with our theoretical
analysis considering lossless materials ($\nu=0$). The corresponding
transmission coefficients in these simulations are approximately $62.20\%$
(effective propagation length $L_{eff}\simeq 1.63\lambda_{0}$), $87.73\%$
($L_{eff}\simeq 2.17\lambda_{0}$), and $86.82\%$ ($L_{eff}\simeq
5.06\lambda_{0}$). Thus, the proposed deep-subwavelength structures exhibit
truly low loss.

Figure 2: Dispersion diagrams of ’XDYX’ structures where the first ’X’ can be
either (a, b) PEC or (c, d) PMC, while the second ’X’ can also be either (a,
c) PEC or (b, d) PMC. The red line and black dashed lines represent the
dispersion curves of the SMPs. The cyan dashed area is the bulk zone of the
bulk YIG. The other parameters are $d_{1}=d_{2}=0.03\lambda_{m}$,
$\omega_{r}=0.7174\omega_{m}$, and $\varepsilon_{d}=1$.

Furthermore, one of the most significant differences between Fig. 2(a) and
Figs. 1(b,c) is the presence of an INZ mode, marked by a green circle in Fig.
2(a). In our previous work, we proposed a method to control the thickness
parameters in order to adjust the INZ modes in magnetized-YIG systems.
Similarly, in the current study, we can employ the same approach to achieve
tunable INZ modes. The advantage of integrability in this system is undeniable
compared to previous systems.

In Fig. 3, we investigate the impact of the thicknesses of YIG ($d_{2}$) and
air ($d_{1}$) on the INZ modes in two cases. In the first case, we assume
$d_{1}=d_{2}$ ($<0.1\lambda_{m}$, subwavelength). As shown in Figs. 2(a) and
2(b), when we increase $d_{1}$, the INZ frequency gradually decreases and
remains close to $\omega=0.5\omega_{m}$. The red line in Fig. 3(a) represents
a fitting curve with the following form:

$\bar{\omega}=5.4848\cdot\bar{d_{1}}^{2}-0.09412\cdot\bar{d_{1}}+0.5085$ (5)

where $\bar{\omega}=\omega/\omega_{m}$ and $\bar{d_{1}}=d_{1}/\lambda_{m}$
($\lambda_{m}=2\pi c/\omega_{m}$). Equation (5) allows us to estimate the INZ
frequency ($f_{\mathrm{INZ}}$) in the $d_{1}=d_{2}$ case. For example, when
$d_{1}=0.045\lambda_{m}$, according to Eq. (5), the INZ frequency should be
approximately $0.4932\omega_{m}$. The inset of Fig. 3(a) further confirms that
when $d_{1}=d_{2}=0.045\lambda_{m}$, a point source with $f=0.4932f_{m}$ can
excite SMPs with a sufficiently large effective wavenumber ($k_{eff}$) and
exhibit no phase shift during propagation. In the second case, we assume
$d_{2}=2d_{1}$. Similar to the first case, we plot $f_{\mathrm{INZ}}$ as a
function of $d_{1}$ in Fig. 3(c). The corresponding fitting equation is as
follows:

$\bar{\omega}=102.45\cdot\bar{d_{1}}^{3}-25.783\cdot\bar{d_{1}}^{2}+0.3716\cdot\bar{d_{1}}+0.5828$
(6)

In comparison to the red line in Fig. 3(a), the red line in Fig. 3(c) exhibits
a sharper slope, resulting in greater changes in $f_{\mathrm{INZ}}$ as $d_{1}$
increases. It is important to emphasize that, in theory, it is possible to
engineer the fitting curve by designing a suitable ratio between $d_{1}$ and
$d_{2}$, such as establishing a nearly linear relationship between
$f_{\mathrm{INZ}}$ and $d_{1}$ within the range of
$0.05\lambda_{m}<d_{1}<0.1\lambda_{m}$. This linear relationship could have
potential applications in phase modulators, optical filters, and frequency
combs.

Figure 3: The operating frequencies of INZ modes as a function of $d_{1}$ in
the cases of (a) $d_{1}=d_{2}$ and (c) $d_{2}=2d_{1}$. The dispersion curves
of one-way SMPs as the thicknesses changing (b) in case (a), and (d) in case
(c). Figure 4:  All-optical phase modulator relied on INZ modes. (a) The phase
angle of $E_{z}$ of the output INZ modes as a function of $L_{x}$. (b)-(e)
Electric filed distributions of the FEM simulations, and $L_{x}$ are equal to
(b) 40 mm, (c) 50 mm, (d) 60 mm, and (e) 70 mm, respectively. (f) The
schematic picture of the INZ modes-based phase modulator. Parameters are
$d_{1}=d_{2}=0.03\lambda_{m}$, $\nu=10^{-3}\omega_{m}$, and $f=0.5007f_{m}$.

It is well known that an INZ wave can travel with no phase shift, which
provides a potential approach to building a phase modulator. By designing a
joint structure that contains two one-way waveguides with the same one-way
region, and if the second waveguide supports the one-way INZ mode, it is
theoretically possible to control the output phase by manipulating the joint
position or the lengths of the waveguides. As analyzed in Fig. 2, the EDYE
structure possesses the same one-way region as the MDYE structure, and an INZ
mode is present in the EDYE structure. Therefore, as shown in Fig. 4(f), we
propose a novel phase modulator called the ’MDYE-EDYE’ structure. We further
investigate how the joint position $L_{x}$ affects the output phase of the
electric field ($E_{z}$). In Fig. 4(a), a clear linear relationship between
the phase (arg($E_{z}$)) and $L_{x}$ is observed, suggesting that such a
simple structure can achieve precise phase modulation.

From Fig. 4(a), it is evident that to achieve a phase shift range of
$[-\pi,\pi]$, the position shift $\Delta L_{x}$ should be approximately 20 mm,
whereas the theoretical value in the lossless condition is around 19.3125 mm
($k\simeq 3.1068k_{m}$). To clearly demonstrate the linear relationship
between the phase angle of $E_{z}$ and $L_{x}$, we perform four full-wave
simulations in a lossy condition with $\nu=10^{-3}\omega_{m}$, where $L_{x}$
is set to 40 mm (Fig. 4(b)), 50 mm (Fig. 4(c)), 60 mm (Fig. 4(d)), and 70 mm
(Fig. 4(e)). Consequently, the output $E_{z}$ in Figs. 3(b) and 3(d) are both
positive with nearly the same phase, as well as in Figs. 3(c) and 3(e).
Therefore, it is reasonable to believe that even in a real lossy condition,
our proposed ’MDYE-EDYE’ structure can function as a precise all-optical phase
modulator, offering advantages such as robustness against backscattering
reflection and a straightforward manufacturing process. Compared to previous
all-optical phase modulators based on Mach-Zehnder Interferometers [35, 36],
metasurfaces [37, 38], and/or metamaterials [39], our proposed structure
exhibits comparable performance and simplicity.

##  III Perfect filter, all-optical splitter, and self-consistent all-optical
logic gates based on INZ modes

Furthermore, through careful assembly, the remanence-based MO system can
function as a perfect filter, effectively separating electromagnetic (EM)
waves within a specific frequency band into two separate channels with
identical bandwidth. Fig. 5(b) illustrates the Y-shaped filter, which consists
of three arms. The horizontal arm serves as the input port, while the other
two arms, namely the ’EDYM’ and ’EDYE’ structures (previously studied in the
previous subsection, Figs. 1(a) and 1(b)), function as the output ports. It is
important to note that the magnetization direction in the ’EDYM’ structure has
been reversed to support the forward propagation of EM waves. The input port
is referred to as the ’EYYE’ structure, with the interconnected YIG layers
having opposite magnetization directions. As shown in Fig. 5(a), the ’EYYE’
structure exhibits a significantly wide one-way region, $0<\omega<\omega_{r}$,
which corresponds to the combined one-way regions of the ’EDYE’ structure
($0.5\omega_{r}<\omega<\omega_{r}$) and the ’EDYM’ structure
($0<\omega<0.5\omega_{r}$).

Figure 5: (a) Dispersion diagram of the YIG-YIG structure. (b) The schematic
picture and (c) simulations of the perfect optical filter with frequencies are
$f=0.5f_{m}$ (marked by blue circle in (a)) and $f=0.2f_{m}$ (marked by black
circle in (a)), respectively.

Benefiting from the unique relationship between the one-way regions of the
input port (’A’) and the output ports (’B’ and ’C’), along with the robustness
of one-way EM waves, it can be theoretically guaranteed that the excited one-
way EM wave from the input port will propagate into either of the output
ports. To verify the performance of the filter, we conducted two finite
element method (FEM) simulations in this structure with frequencies set at
$f=0.5f_{m}$ (the second image in Fig. 5(c)) and $f=0.2f_{m}$ (the third image
in Fig. 5(c)). As a result, the excited wave from ’A’ propagated to the
junction point and decisively injected into either ’B’ or ’C’, resulting in an
extremely high contrast ratio between the output ports. Therefore, the
simulation results align with our theoretical analysis. The proposed perfect
filter has the potential to serve crucial roles such as an optical switch in
an integrated optical circuit. By directly changing the working frequency, the
state of the communication system can be switched between ’work’ and ’stop’.

Figure 6: (a) Dispersion diagrams of three arms of a novel tunable splitter
shown in (c). (b) Relationship between transmittances of two output arms and
the operating frequency within the one-way region (colored area in (a)). (c)
The schematic diagram and simulations for splitting ratio $\eta=0/100$ (second
picture) and $\eta=50/50$ (last picture). Other parameters are
$d_{1}=d_{2}=0.01\lambda_{m}$ and $\nu=10^{-3}\omega_{m}$.

Besides, it is worth noting that the ’EDYE’, ’MDYE’, and ’EYYE’ structures
have overlapping one-way regions, but the wavenumbers of the one-way SMPs
differ. As a result, one-way SMPs with different frequencies propagating in
these structures will have different effective refractive indexes ($n_{eff}$)
and may induce different splitting ratios ($\eta$). Fig. 6(c) illustrates the
designed Y-shaped splitter, consisting of the ’MDYE’ (’A’), ’EYYE’ (’B’), and
’EDYE’ (’C’) structures. Figure 6(a) presents the dispersion curves of SMPs in
the vicinity of the overlapping one-way region
($0.5\omega_{r}\leq\omega<\omega_{r}$). It is important to note that we set
$d_{1}=d_{2}=0.01\lambda_{m}$ to engineer an SMP with an infinite wavenumber
at $\omega=0.5\omega_{r}$ (the lower limit of the shaded region in Fig. 6(a))
in the ’EDYE’ structure, which implies $n_{eff}\to\infty$ and consequently
zero transmittance in the ’A-C’ channel should be observed under this
condition. We performed FEM simulations in the designed splitter, with the
frequency changing within the one-way region, and calculated the
transmittance. As depicted in Fig. 6(b), the red and black lines represent the
transmittance of the ’A-B’ and ’A-C’ channels, respectively. As expected, when
$\omega=0.5\omega_{r}=0.3587\omega_{m}$, nearly 100% of the energy traveled
through the ’A-B’ channel (as shown in the second picture in Fig. 6(c)). More
importantly, as the frequency gradually increases to $0.4447\omega_{m}$, the
splitting ratio can approach 50/50 (as shown in the last picture in Fig.
6(c)), which is a highly desirable value in the field of optical
communication. Another noteworthy case occurs when $\omega\simeq
0.45\omega_{m}$ (indicated by the right arrow in Fig. 6(b)). In this scenario,
the EM wave cannot pass through the ’A-C’ channel either. This behavior can be
explained by the theory of total internal reflection: 1) the effective indices
of ’A’, ’B’, and ’C’ are approximately 19.4082 ($k\simeq 10.6745k_{m}$),
29.4398 ($k\simeq 16.1919k_{m}$), and 6.8467 ($k\simeq 3.7657k_{m}$),
respectively; 2) the incident angle is $\pi/3$, which is clearly larger than
the ’critical angle’ (equal to $arcsin(6.8467/19.4082)$) of the ’A-C’ channel.
When the frequency is further increased to approximately $0.6456\omega_{m}$, a
50/50 splitting ratio will be achieved (indicated by the black left arrow in
Fig. 6(b)). Therefore, our proposed remanence-based all-optical splitter has
the capability to maintain a constant splitting ratio for different
frequencies. For example, it can achieve $\eta=50/50$ at
$\omega=0.4447\omega_{m}$ and/or $\omega=0.6456\omega_{m}$, which we believe
is extremely important for parallel optical computing.

Figure 7:  (a) Schematic picture of all-optical LGs based on YIG with
remanence. (b) Logic operations and FEM simulations for three different
inputs, i.e. (’1’, ’0’), (’0’, ’1’), and (’1’, ’1’). (c) The truth table for
OR, AND, and NAND gates according to simulation results shown in (b).

In our previous work[19], we proposed a possible method to achieve all-optical
LGs with an extremely high contrast ratio. However, there was an unresolved
issue in the previously proposed LGs, which involved inconsistency between the
logic of the input and output ports. In this study, we propose another
Y-shaped system based on INZ modes, which can function as all-optical LGs such
as OR, AND, and NAND gates. Following a consistent positive logic convention,
the presence of EM energy is treated as logic ’1’ and is utilized in both the
input and output ports, while the phase angle of $E_{z}$ serves as the
reference quantity in the output port. Figure 7(a) illustrates the schematic
diagram of the designed system. One can observe that this structure is very
similar to the one shown in Fig. 6(c), except for the exchange of YIG layers
in arm ’B’. As discussed earlier, all three arms share the same one-way
region, $0.5\omega_{r}<\omega<\omega_{r}$. In contrast to arm ’B’ shown in
Fig. 6, the arm in Fig. 7 can sustain backward-propagating SMPs (as shown in
the inset of Fig. 7(a)). Consequently, two one-way channels are successfully
engineered based on the splitter model shown in Fig. 6. In contrast to the
working frequency used in our previous work, we choose the INZ frequency in
the output port as the working frequency in order to preserve the phase
information in the output section.

The proposed system can be theoretically treated as a combination of basic LGs
such as OR, AND, and NAND gates. For the OR gate, we utilize positive logic,
where any presence of EM energy is considered logic ’1’. Based on the robust
one-way channels, any input EM signal (logic ’1’) from either arm ’A’ or arm
’B’ can unidirectionally propagate through arm ’C’ (also logic ’1’). Hence,
the system naturally functions as an OR gate.

To implement an AND gate, we introduce a reference quantity, the phase angle
($\theta$) of $E_{z}$, to determine the state of the output port.
Specifically, we assume that the output EM signal with $\pi/2<\theta<\pi/2$
represents logic ’1’. As illustrated in Fig. 7(b), three input combinations,
[’1’, ’0’], [’0’, ’1’], and [’1’, ’1’], yield three corresponding outputs,
logic ’0’, logic ’0’, and logic ’1’. It should be noted that we set
$\theta\simeq 1.4$ rad for input-1 to engineer an appropriate output phase
angle. Therefore, in practical applications, the input $\theta$ needs to be
tailored for different inputs. In Fig. 7, the inputs $[\theta_{1},\theta_{2}]$
are [-1.9 rad, ], [, 1.4 rad], and [1.4 rad, 1.4 rad], respectively.

Figure 7(c) presents the truth table for the OR, AND, and NAND gates. It is
evident that the outputs are inverted for the AND and NAND gates. Therefore,
based on the analysis of the AND gate, the system can easily achieve NAND
operation by assuming that the output EM signal with $\pi/2<\theta<\pi/2$
represents logic ’0’. Furthermore, other basic LGs can be straightforwardly
extended. With improved integration and consistent logic, the proposed
Y-shaped LGs enhance the potential applications of one-way mode-based LGs in
the all-optical calculation of integrated optical circuits.

Figure 8:  Realization of 3D all-optical phase modulator, perfect filter, and
self-consistent LGs.

The above-mentioned all-optical phase modulator, filter, and LGs, which are
based on YIG with remanence, can be straightforwardly extended from 2D models
to 3D structures[34]. As depicted in the inset of Fig. 8(a), two PEC walls are
introduced in the z direction to confine the guided wave within the device.
The simulation result shown in Fig. 8(a) exhibits similarities to the
corresponding 2D case (Fig. 4(e)). Figures 7(b)-(e) illustrate the simulation
outcomes for the 3D perfect filter, 3D splitter, and 3D LGs, respectively, and
all the results are consistent with their 2D counterparts. Consequently, our
proposed all-optical structures hold great potential for practical device
implementation. Due to the characteristic of one-way propagation, the
manufacturing process should be significantly simpler compared to other
competing technologies. More importantly, as depicted in Fig. 1, our proposed
all-optical phase modulator, filter, splitter, and LGs can be easily assembled
on a chip. Therefore, it holds significant promise for enabling flexible all-
optical computation and communication.

##  IV Conclusion

In conclusion, we have proposed a series of all-optical devices based on YIG
with remanence. The all-optical phase modulator utilizes INZ modes and allows
for controllable phase modulation in a linear relationship with the length of
the boundary. The INZ frequency can be adjusted by manipulating the thickness
parameters, providing tunability. Moreover, we have designed a perfect filter,
which separates the input port’s one-way region into two parts and creates two
one-way channels with equal bandwidth. Additionally, we have introduced an
all-optical splitter that divides the input wave into two output ports with a
customizable and precise splitting ratio, such as 50/50. Furthermore, we have
proposed all-optical logic gates based on one-way SMPs, where consistent
positive logic and the phase angle of the output field are used to determine
the logic states. Basic logic gates, including OR, AND, and NAND, have been
achieved based on this principle. The feasibility and performance of the
proposed devices have been validated in both 2D and 3D simulations. Our
proposed all-optical devices, leveraging MO heterostructures, hold potential
for optical calculations, including parallel computing in integrated optical
circuits.

## Acknowledgement

This work was supported by the National Natural Science Foundation of Sichuan
Province (No. 2023NSFSC1309), and the open fund of Luzhou Key Laboratory of
Intelligent Control and Application of Electronic Devices (No. ZK202210),
Sichuan Science and Technology Program (No. 2022YFS0616), the Science and
Technology Strategic Cooperation Programs of Luzhou Municipal People’s
Government and Southwest Medical University (No. 2019LZXNYDJ18). J.X. and Y.L.
thanks for the support of the Innovation Laboratory of Advanced Medical
Material & Physical Diagnosis and Treatment Technology. K.L.T. was supported
by the General Secretariat for Research and Technology (GSRT) and the Hellenic
Foundation for Research and Innovation (HFRI) under Grant No. 4509.

## References

  * [1] Limin Tong, Rafael R. Gattass, Jonathan B. Ashcom, Sailing He, Jingyi Lou, Mengyan Shen, Iva Maxwell, and Eric Mazur.  Subwavelength-diameter silica wires for low-loss optical wave guiding.  Nature, 426(6968):816–819, 2003.
  * [2] J. Feldmann, N. Youngblood, M. Karpov, H. Gehring, X. Li, M. Stappers, M. Le Gallo, X. Fu, A. Lukashchuk, A. S. Raja, J. Liu, C. D. Wright, A. Sebastian, T. J. Kippenberg, W. H. P. Pernice, and H. Bhaskaran.  Parallel convolutional processing using an integrated photonic tensor core.  Nature, 589(7840):52–58, 2021.
  * [3] Rongjin Zhuang, Kai Ni, Guanhao Wu, Ting Hao, Longzhao Lu, Yang Li, and Qian Zhou.  Electro-optic frequency combs: Theory, characteristics, and applications.  Laser & Photonics Reviews, page 2200353, 2023.
  * [4] Viola Valentina Vogler-Neuling, Artemios Karvounis, Andrea Morandi, Helena Weigand, Eric Dénervaud, and Rachel Grange.  Photonic Assemblies of Randomly Oriented Nanocrystals for Engineered Nonlinear and Electro-Optic Effects.  ACS Photonics, 9(7):2193–2203, 2022.
  * [5] Xiaojun Chen, Jiao Lin, and Ke Wang.  A Review of Silicon-Based Integrated Optical Switches.  Laser & Photonics Reviews, 17(4):2200571, 2023.
  * [6] Mengxin Ren, Wei Cai, and Jingjun Xu.  Tailorable Dynamics in Nonlinear Optical Metasurfaces.  Advanced Materials, 32(3):1806317, 2020.
  * [7] Jie Xu, Xiaohua Deng, Hang Zhang, Chiaho Wu, Martijn Wubs, Sanshui Xiao, and Linfang Shen.  Ultra-subwavelength focusing and giant magnetic-field enhancement in a low-loss one-way waveguide based on remanence.  Journal of Optics, 22(2):025003, 2020.
  * [8] K. L. Tsakmakidis, L. Shen, S. A. Schulz, X. Zheng, J. Upham, X. Deng, H. Altug, A. F. Vakakis, and R. W. Boyd.  Breaking lorentz reciprocity to overcome the time-bandwidth limit in physics and engineering.  Science, 356(6344):1260–1264, 2017.
  * [9] Shiqing Li, Kosmas L. Tsakmakidis, Qian Shen, Hang Zhang, Jinhua Yan, Shulin Sun, and Linfang Shen.  Broadband unidirectional guided-wave-driven metasurfaces for arbitrary wavefront control, 2023.
  * [10] Farzad Zangeneh-Nejad, Dimitrios L Sounas, Andrea Alù, and Romain Fleury.  Analogue computing with metamaterials.  Nature Reviews Materials, 6(3):207–225, 2021.
  * [11] Guo-Jing Tang, Xin-Tao He, Fu-Long Shi, Jian-Wei Liu, Xiao-Dong Chen, and Jian-Wen Dong.  Topological photonic crystals: physics, designs, and applications.  Laser & Photonics Reviews, 16(4):2100300, 2022.
  * [12] Guancong Ma, Meng Xiao, and Che Ting Chan.  Topological phases in acoustic and mechanical systems.  Nature Reviews Physics, 1(4):281–294, 2019.
  * [13] Yin Poo, Rui-xin Wu, Zhifang Lin, Yan Yang, and C. T. Chan.  Experimental Realization of Self-Guiding Unidirectional Electromagnetic Edge States.  Physical Review Letters, 106(9):093903, 2011.
  * [14] Peiheng Zhou, Gui-Geng Liu, Xin Ren, Yihao Yang, Haoran Xue, Lei Bi, Longjiang Deng, Yidong Chong, and Baile Zhang.  Photonic amorphous topological insulator.  Light: Science & Applications, 9(1):133, 2020.
  * [15] Mudi Wang, Ruo-Yang Zhang, Lei Zhang, Dongyang Wang, Qinghua Guo, Zhao-Qing Zhang, and Che Ting Chan.  Topological one-way large-area waveguide states in magnetic photonic crystals.  Physical Review Letters, 126(6):067401, 2021.
  * [16] Jie Xu, Panpan He, Delong Feng, Kangle Yong, Lujun Hong, Yun Shen, and Yun Zhou.  Slow wave and truly rainbow trapping in a one-way terahertz waveguide.  Optics Express, 29(7):11328–11341, 2021.
  * [17] Jie Xu, Qian Shen, Kai Yuan, Xiaohua Deng, Yun Shen, Hang Zhang, Chiaho Wu, Sanshui Xiao, and Linfang Shen.  Trapping and releasing bidirectional rainbow at terahertz frequencies.  Optics Communications, 473:125999, 2020.
  * [18] Yun Zhou, Panpan He, Sanshui Xiao, Fengwen Kang, Lujun Hong, Yun Shen, Yamei Luo, and Jie Xu.  Realization of tunable index-near-zero modes in nonreciprocal magneto-optical heterostructures.  Optics Express, 30(15):27259–27272, 2022.
  * [19] Jie Xu, Fengwen Kang, Yamei Luo, Sanshui Xiao, and Kosmas L Tsakmakidis.  All-optical digital logic based on unidirectional modes.  Advanced Optical Materials, 11(1):2201836, 2023.
  * [20] S Ali Hassani Gangaraj and Francesco Monticone.  Do truly unidirectional surface plasmon-polaritons exist?  Optica, 6(9):1158–1165, 2019.
  * [21] Kosmas L. Tsakmakidis, Konstantinos Baskourelos, and Tomasz Stefański.  Topological, nonreciprocal, and multiresonant slow light beyond the time-bandwidth limit.  Applied Physics Letters, 119(19):190501, 2021.
  * [22] JJ Brion, RF Wallis, A Hartstein, and E Burstein.  Theory of surface magnetoplasmons in semiconductors.  Physical Review Letters, 28(22):1455, 1972.
  * [23] Yiyu Zhou, M Zahirul Alam, Mohammad Karimi, Jeremy Upham, Orad Reshef, Cong Liu, Alan E Willner, and Robert W Boyd.  Broadband frequency translation through time refraction in an epsilon-near-zero material.  Nature communications, 11(1):1–7, 2020.
  * [24] Xinxiang Niu, Xiaoyong Hu, Saisai Chu, and Qihuang Gong.  Epsilon-near-zero photonics: a new platform for integrated devices.  Advanced Optical Materials, 6(10):1701292, 2018.
  * [25] Mário Silveirinha and Nader Engheta.  Tunneling of Electromagnetic Energy through Subwavelength Channels and Bends using $\varepsilon$ -Near-Zero Materials.  Physical Review Letters, 97(15):157403, 2006.
  * [26] João S. Marcos, Mário G. Silveirinha, and Nader Engheta.  $\mu$ -near-zero supercoupling.  Physical Review B, 91(19):195112, 2015.
  * [27] Ahmed M. Mahmoud and Nader Engheta.  Wave–matter interactions in epsilon-and-mu-near-zero structures.  Nature Communications, 5(1):5638, 2014.
  * [28] Takian Fakhrul, Stana Tazlaru, Lukáš Beran, Yan Zhang, Martin Veis, and Caroline A. Ross.  Magneto-Optical Bi:YIG Films with High Figure of Merit for Nonreciprocal Photonics.  Advanced Optical Materials, 7(13):1900056, 2019.
  * [29] Jiachang Liang, Yan Li, Tingge Dai, Yuejun Zhang, Xiaowei Zhang, Hongjun Liu, and Pengjun Wang.  On-chip Ce:YIG/Si Mach–Zehnder optical isolator with low power consumption.  Optics Express, 31(5):8375–8383, 2023.
  * [30] Maksym Popov, Yuzan Xiong, Igor Zavislyak, Hryhorii Chumak, Oleksandr Fedorchuk, Sujoy Saha, Rao Bidthanapally, Hongwei Qu, Michael R Page, and Gopalan Srinivasan.  Y-type hexagonal ferrite-based band-pass filter with dual magnetic and electric field tunability.  Scientific Reports, 13(1):1179, 2023.
  * [31] E. Almpanis, G. P. Zouros, P. A. Pantazopoulos, K. L. Tsakmakidis, N. Papanikolaou, and N. Stefanou.  Spherical optomagnonic microresonators: Triple-resonant photon transitions between Zeeman-split Mie modes.  Physical Review B, 101(5):054412, 2020.
  * [32] Constantine A Balanis.  Advanced Engineering Electromagnetics.  John Wiley & Sons, 2012.
  * [33] Qian Shen, Xiaodong Zheng, Hang Zhang, Yun You, and Linfang Shen.  Large-area unidirectional surface magnetoplasmons using uniaxial -near-zero material.  Optics Letters, 46(23):5978, 2021.
  * [34] Lujun Hong, Yazhou Wang, Yun Shen, Xiaohua Deng, Kai Yuan, Sanshui Xiao, and Jie Xu.  Broadband energy squeezing and tunneling based on unidirectional modes.  Optical Materials Express, 11(9):2975–2984, 2021.
  * [35] C Sturm, D Tanese, H S Nguyen, H Flayac, E Galopin, A Lemaıtre, I Sagnes, D Solnyshkov, A Amo, G Malpuech, and J Bloch.  All-optical phase modulation in a cavity-polariton Mach–Zehnder interferometer.  NATURE COMMUNICATIONS, 2014.
  * [36] Chong Zhang, Paul A. Morton, Jacob B. Khurgin, Jon D. Peters, and John E. Bowers.  Ultralinear heterogeneously integrated ring-assisted Mach–Zehnder interferometer modulator on silicon.  Optica, 3(12):1483, 2016.
  * [37] Ziqi Miao, Qiong Wu, Xin Li, Qiong He, Kun Ding, Zhenghua An, Yuanbo Zhang, and Lei Zhou.  Widely Tunable Terahertz Phase Modulation with Gate-Controlled Graphene Metasurfaces.  Physical Review X, 5(4):041027, 2015.
  * [38] Weihao Yang, Jun Qin, Jiawei Long, Wei Yan, Yucong Yang, Chaoyang Li, En Li, Juejun Hu, Longjiang Deng, Qingyang Du, and Lei Bi.  A self-biased non-reciprocal magnetic metasurface for bidirectional phase modulation.  Nature Electronics, 6(3):225–234, 2023.
  * [39] Seojoo Lee, Soojeong Baek, Teun-Teun Kim, Hyukjoon Cho, Sangha Lee, Ji-Hun Kang, and Bumki Min.  Metamaterials for Enhanced Optical Responses and their Application to Active Control of Terahertz Waves.  Advanced Materials, 32(35):2000250, 2020.