# Multiple Model Reference Adaptive Control with Blending for Non-Square
Multivariable Systems

Alex Lovi, Barış Fidan, and Christopher Nielsen This work is supported by the
Natural Sciences and Engineering Research Council of Canada (NSERC). This work
has been submitted to the IEEE for possible publication. Copyright may be
transferred without notice, after which this version may no longer be
accessible.Department of Electrical and Computer Engineering, University of
Waterloo, Waterloo ON, N2L 3G1 Canada (e-mail: alovi@uwaterloo.ca).Department
of Electrical and Computer Engineering, University of Waterloo, Waterloo ON,
N2L 3G1 Canada (e-mail: cnielsen@uwaterloo.ca).Department of Mechanical and
Mechatronics Engineering, University of Waterloo, Waterloo ON, N2L 3G1 Canada
(e-mail: fidan@uwaterloo.ca).

###### Abstract

In this paper we develop a multiple model reference adaptive controller
(MMRAC) with blending. The systems under consideration are non-square, i.e.,
the number of inputs is not equal to the number of states; multi-input,
linear, time-invariant with uncertain parameters that lie inside of a known,
compact, and convex set. Moreover, the full state of the plant is available
for feedback. A multiple model online identification scheme for the plant’s
state and input matrices is developed that guarantees the estimated parameters
converge to the underlying plant model under the assumption of persistence of
excitation. Using an exact matching condition, the parameter estimates are
used in a control law such that the plant’s states asymptotically track the
reference signal generated by a state-space model reference. The control
architecture is proven to provide boundedness of all closed-loop signals and
to asymptotically drive the state tracking error to zero. Numerical
simulations illustrate the stability and efficacy of the proposed MMRAC
scheme.

###### Index Terms:

Adaptive control, Model reference adaptive control, Multiple model, Polytopic
uncertainty, Uncertain systems.

##  I Introduction

The use of multiple models to describe the dynamics of uncertain systems has
been studied [1, 2, 3, 4], and shown to improve the transient-time performance
[5], the steady-state tracking performance [4, 6], and robustness [7] when
compared to a single model approach. Multiple model control techniques can be
broadly divided into _switching_ control, which allows for the selection of a
best model in different dynamic situations [8, 7], and _blending_ control,
where the information from different models is combined to get a single
description of the system [9, 10, 3, 11]. In this paper we opt for a blending
technique since it allows for better closed-loop performance when compared to
any single model technique [12], and avoids possible undesirable transient-
time behavior that switching control may exhibit [13, 14].

Multiple model approaches with blending have various promising applications in
identification and control of time-varying (TV) systems, including adaptive
identification of MIMO, linear, periodically TV systems (with known period)
[15] and improvement of the transients and adaptation speed in adaptive
control of uncertain, TV, MIMO systems [9]. The stability and robustness
properties were studied in [3]. Mixing adaptive techniques have also been used
to achieve faster tracking for a class of nonlinear discrete-time systems
[16]. For the case of experimental results, the use of mixing control has
allowed to decrease overshoots, settling time, steady state error [17, 18],
and designing fault tolerant controllers [10]. Other applications include
multiple model estimation of power systems models [19], development of an
adaptive controller for partially-observed Boolean dynamical systems [20], and
distributed state estimation using a network of local sensors [21].

Model reference control has been rigorously and methodically studied for many
years and proposes promising applications with detailed design procedures.
When we consider systems with uncertainty, model reference adaptive
controllers (MRAC) are a versatile technique that achieves robust and
satisfactorily closed-loop performances [22, 23]. The combination of MRAC with
multiple model approaches to obtain a continuous input signal calculated using
all the identification errors from all the models is studied in [1, 24].
Similar work is present in [25], where adaptive identification models are
considered instead of fixed models. In [9], a similar identification scheme is
integrated with linear quadratic optimal controllers to design a multiple MRAC
(MMRAC) scheme for tracking reference signals generated by a linear reference
model. The asymptotic tracking of the MMRAC scheme for linear, time-invariant
(LTI), MIMO systems is proven and simulations are presented to evaluate and
validate the performance.

In this paper we consider non-square, multi-input, LTI systems with polytopic
parameter uncertainty, whose unknown plant models are in the interior of the
convex hull of a finite number of fixed models. Based on full state
measurements, we develop a parameter identification scheme that estimates a
weight vector which determines the convex combination that yields the state-
space representation of the uncertain system.

This paper extends our previous work [26] in several significant ways. The
systems under consideration do not have to be square, i.e., the number of
inputs can be less or equal to the number of states. Moreover, unlike in [26],
in this article the number of fixed models that define the polytopic
uncertainty of the plant can be arbitrary. This relaxation allows us to
develop a procedure to obtain a set of fixed corner systems to use in the
identification process. Furthermore, we provide sufficient conditions under
which the parameter estimates asymptotically converge to the plant’s uncertain
parameters. The simulations presented here point to faster convergence of the
parameter estimates to their true values, making the state error of the
proposed MMRAC scheme also converge faster compared to a single model MRAC
approach.

The problem formulation is presented in Section II, together with the
assumptions. In Section III, we present a selection process for the corner
systems. In Section IV, we present the identification scheme with its
stability analysis. We use these results in Section V to develop the MMRAC
scheme. A set of simulations are presented comparing the MMRAC scheme to the
single model case for uncertain systems tracking references in Section VI
using MATLAB and Simulink based simulations. We finalize with the conclusions
in Section VII.

#### Notation

If $S\subset\mathbb{R}^{n}$, then $\operatorname{conv}\left(S\right)$ denotes
the closed convex hull of $S$, the interior of $S$ is written
$\operatorname{int}(S)$. Let $\|\cdot\|$ denote the 2-norm of both vectors and
matrices. The kernel of a matrix $A$ is denoted by $\operatorname{Ker}A$.
Given two matrices $Q,R\in\mathbb{R}^{n\times n}$, we write $Q\succeq R$ if
$Q-R$ is positive semi-definite. For two signals $f_{1}$ and $f_{2}$ we write
$f_{1}\equiv f_{2}$, if there exist $\alpha$, $\beta>0$ such that for all
$t\geq 0$, ${\|f_{1}\left(t\right)-f_{2}\left(t\right)\|\leq\beta e^{-\alpha
t}}$. Further, if $e_{f}=f_{1}-f_{2}$ satisfies $\dot{e}_{f}(t)=-\alpha
e_{f}(t)$ then we write $f_{1}\stackrel{{\scriptstyle\alpha}}{{=}}f_{2}$. If
$x\in\mathbb{R}^{N}$, then $x_{i}$ denotes the $i$-th component, and
${\bar{x}=\left(x_{1},\cdots,x_{N-1}\right)\in\mathbb{R}^{N-1}}$. A signal
$\Phi:[0,\infty)\rightarrow\mathbb{R}^{k}$ is persistently exciting (PE) if
there exist $\alpha_{\Phi 1},\alpha_{\Phi 2},T_{\Phi}>0$ such that for every
$t\geq 0$

$\alpha_{\Phi
1}\mathbb{I}\preceq\int_{t}^{t+T_{\Phi}}\Phi(\tau)\Phi^{\top}(\tau)\mathrm{d}\tau\preceq\alpha_{\Phi
2}\mathbb{I}.$ (1)

##  II Problem Formulation

Consider the MIMO, LTI system

$\dot{x}_{p}\left(t\right)=A_{p}x_{p}\left(t\right)+B_{p}u\left(t\right),\>x_{p}\left(0\right)=x_{p0},\>t\geq
0,$ (2)

where $A_{p}\in\mathbb{R}^{n\times n}$ and $B_{p}\in\mathbb{R}^{n\times m}$
are unknown constant matrices, $B_{p}$ is full column rank,
$x_{p}\left(t\right)\in\mathbb{R}^{n}$, $u\left(t\right)\in\mathbb{R}^{m}$ are
the state of the system and the control input, respectively. We assume that
$x_{p}(t)$ is available for feedback and that the number states $n$ and inputs
$m$ are known.

###### Remark 1

The results of this article can be applied to systems with structured
uncertainties of the form
${\left(A_{p}\left(\eta\right),B_{p}\left(\eta\right)\right)}$, where
${A_{p}\left(\eta\right)\in\mathbb{R}^{n\times n}}$ and
${B_{p}\left(\eta\right)\in\mathbb{R}^{n\times m}}$ are unknown constant
matrices dependent on the uncertain parameter vector ${\eta\in
S\subset\mathbb{R}^{q}}$, belonging to a compact set.  $\triangle$

This article aims to design a controller such that the state $x_{p}$ of the
plant (2) asymptotically tracks the signal $x_{r}$ generated by the reference
model

$\dot{x}_{r}\left(t\right)=A_{r}x_{r}\left(t\right)+B_{r}r\left(t\right),\>x_{r}\left(0\right)=x_{r0},\>t\geq
0,$ (3)

where $A_{r}\in\mathbb{R}^{n\times n}$ and $B_{r}\in\mathbb{R}^{n\times m}$
are known constant matrices, $A_{r}$ is Hurwitz, and
$r:\left[0,\infty\right)\rightarrow\mathbb{R}^{m}$ is a known, bounded,
piecewise continuous reference signal. It is assumed that the plant (2) and
the reference model (3) satisfy the exact matching conditions [23], as stated
in the following assumption.

###### Assumption 1

There exist matrices $K^{*}\in\mathbb{R}^{m\times n}$ and
${L^{*}\in\mathbb{R}^{m\times m}}$ such that

$\displaystyle A_{p}+B_{p}K^{*}$ $\displaystyle=A_{r},$ (4) $\displaystyle
B_{p}L^{*}$ $\displaystyle=B_{r}.$ (5)

$\triangle$

Assumption 1 is a necessary and sufficient condition to guarantee the
existence of a solution to the tracking problem when the plant’s model (2) is
known [23]. In addition, (5) is equivalent to
$\operatorname{Im}B_{r}\subseteq\operatorname{Im}B_{p}$, which means that,
without loss of generality, we can assume that
$\operatorname{Im}B_{p}=\operatorname{Im}B_{r}$.

Our approach to solving the aforementioned control design problem involves
online parameter identification of the system matrices
$\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}\eqqcolon\Theta_{p}$ by defining a
compact, convex, uncertainty polytope such that for every extreme point of the
polytope, also referred to as corner, there is a fixed model

$\displaystyle\dot{x}_{i}(t)$ $\displaystyle=A_{i}x_{i}(t)+B_{i}u(t),$

with system matrices
$\begin{bmatrix}A_{i}&B_{i}\end{bmatrix}\eqqcolon\Theta_{i}$,
$i\in\left\\{1,\cdots,N\right\\}$. We define the set

$\mathcal{S}\coloneqq\left\\{\begin{bmatrix}A_{1}&B_{1}\end{bmatrix},\cdots,\begin{bmatrix}A_{N}&B_{N}\end{bmatrix}\right\\}$
(6)

consisting of every corner model.

###### Assumption 2

There exist a finite set $\mathcal{S}$ of known system matrices
$\begin{bmatrix}A_{i}&B_{i}\end{bmatrix}\in\mathbb{R}^{n\times(n+m)}$, such
that

  1. i)

${\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}\in\operatorname{int}\left(\operatorname{conv}\left(\mathcal{S}\right)\right)}$.

  2. ii)

Every convex combination of the $B_{i}$ is full column rank.

$\triangle$

The first item of Assumption 2 can be achieved from a system identification
process which we describe in Section III-B. We require Assumption 2 (ii) to
ensure that we do not have redundant inputs, or we do not drop rank of the
number of inputs, losing control authority of the system. Given a set
$\mathcal{S}$, there exist numerical methods to verify whether every convex
combinations of the $B_{i}$’s is full rank [27].

Every point of the polytope $\operatorname{conv}\left(\mathcal{S}\right)$ can
be expressed as a convex combination of the corner models
$\begin{bmatrix}A_{i}&B_{i}\end{bmatrix}$, which implies that the following
set is non-empty:

$\mathcal{W}\coloneqq\left\\{w\in\left[0,1\right]^{N}:\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}=\sum_{i=1}^{N}w_{i}\begin{bmatrix}A_{i}&B_{i}\end{bmatrix},\sum_{i=1}^{N}w_{i}=1\right\\}.$
(7)

Consequently, the problem of identifying the unknown system matrix
$\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$ is equivalent to the problem of
identifying a vector $w\in\mathcal{W}$. A preliminary study of the problem
with $N=n+1$, and $m=n$ was studied in [26]. In this paper unlike in [26], we
study the problem for any number of inputs $m>0$, and an arbitrary number $N$
of corner models.

##  III Corner Model Selection

Consider a set $\mathcal{S}$ that satisfies Assumption 2 and has $N$ elements.
In this section we develop a constructive procedure that, starting with
$\mathcal{S}$, produces a new set

$\displaystyle\mathcal{S}^{\prime}=\left\\{\begin{bmatrix}A_{1}^{\prime}&B_{1}^{\prime}\end{bmatrix},\cdots,\begin{bmatrix}A^{\prime}_{N^{\prime}}&B^{\prime}_{N^{\prime}}\end{bmatrix}\right\\},$
(8)

that also satisfies Assumption 2, and such that for every
$i\in\left\\{1,\cdots,N\right\\}$ there exist matrices
$K_{i}\in\mathbb{R}^{m\times n}$ and ${L_{i}\in\mathbb{R}^{m\times m}}$ such
that,

$\displaystyle A_{i}^{\prime}+B_{i}^{\prime}K_{i}$ $\displaystyle=A_{r},$ (9)
$\displaystyle B_{i}^{\prime}L_{i}$ $\displaystyle=B_{r}.$ (10)

###  III-A Satisfying the Matching Conditions for the Corner Models

In the next proposition we show how to use the information of Assumption 1 and
a set $\mathcal{S}$ that satisfies Assumption 2 to obtain a new set
$\mathcal{S}^{\prime}$ such that for every element of $\mathcal{S}^{\prime}$
there exist matrices $K_{i}\in\mathbb{R}^{m\times n}$ and
${L_{i}\in\mathbb{R}^{m\times m}}$ that satisfy (9) and (10).

###### Proposition 1

Suppose that the plant (2) and reference model (3) are such that Assumption 1
is satisfied. If there exists a set $\mathcal{S}$ that satisfies Assumption 2,
then there exists a set $\mathcal{S}^{\prime}$ that also satisfies Assumption
2, and furthermore

  1. 1.

$\operatorname{conv}(\mathcal{S}^{\prime})\subseteq\operatorname{conv}(\mathcal{S})$.

  2. 2.

For each $i\in\left\\{1,\cdots,N\right\\}$, there exist matrices
$K_{i}\in\mathbb{R}^{m\times n}$ and ${L_{i}\in\mathbb{R}^{m\times m}}$ such
that (9) and (10) are satisfied.

###### Proof:

Let $(A_{p},B_{p})$ and $(A_{r},B_{r})$ be the matrix pairs that represent a
plant, and a reference model, respectively. Assume that they satisfy
Assumption 1. Moreover, assume that there exists a set $\mathcal{S}$ with $N$
elements that satisfies Assumption 2. The first, and trivial case, is if for
every $i=1,\cdots,N$, there exist matrices $K_{i}\in\mathbb{R}^{m\times n}$
and ${L_{i}\in\mathbb{R}^{m\times m}}$ such that every corner model
$[A^{\prime}_{i}~{}~{}B^{\prime}_{i}]:=[A_{i}~{}~{}B_{i}]$ satisfies (9) and
(10). In that case, we can define $\mathcal{S}^{\prime}=\mathcal{S}$, and we
have satisfied Assumptions 1 and 2.

For the case when there exists some $j\in\\{1,\cdots,N\\}$ such that
$[A^{\prime}_{j}~{}~{}B^{\prime}_{j}]:=[A_{j}~{}~{}B_{j}]\in\mathcal{S}$ does
not satisfy (9) or (10), consider the closed convex hull of $\mathcal{S}$,
i.e., $\operatorname{conv}\left(\mathcal{S}\right).$ The set
$\operatorname{conv}\left(\mathcal{S}\right)$ is a compact, convex polytope,
and it is non-empty since $\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$ is an
element of $\operatorname{conv}\left(\mathcal{S}\right)$ by Assumption 2 (i).
Next, take the set of matrices that satisfy Equations 4 and 5:

$\displaystyle\mathcal{T}\coloneqq$
$\displaystyle\big{\\{}\begin{bmatrix}A_{x}&B_{x}\end{bmatrix}\in\mathbb{R}^{n\times(n+m)}:$
$\displaystyle A_{x}+B_{x}K=A_{r},B_{x}L=B_{r},L\in\mathbb{R}^{m\times
m},K\in\mathbb{R}^{m\times n}\big{\\}}.$

The set $\mathcal{T}$ is a hyperplane, since it can be written as $n(n+m)$
linear equations (see Section 2.2.1 of [28]), therefore it is a convex
polytope. From Assumption 1, we see that if we replace $K$ and $L$ by $K^{*}$
and $L^{*}$, respectively, the element
$\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$ belongs to $\mathcal{T}$, which
means $\mathcal{T}$ is non-empty. The intersection set

$\mathcal{P}\coloneqq\operatorname{conv}\left(\mathcal{S}\right)\cap\mathcal{T}.$

is also non-empty, since $\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$ is an
element of both $\operatorname{conv}\left(\mathcal{S}\right)$ and
$\mathcal{T}$. Moreover, the intersection of a compact, convex polytope and a
convex polytope is another compact, convex polytope (see Section 2.3.1 of
[28]), which in turn implies that we can obtain $N^{\prime}$ corner models
$\begin{bmatrix}A_{i}^{\prime}&B_{i}^{\prime}\end{bmatrix}$ such that

$\mathcal{P}=\operatorname{conv}\Big{(}\underbrace{\left\\{\begin{bmatrix}A_{1}^{\prime}&B_{1}^{\prime}\end{bmatrix},\cdots,\begin{bmatrix}A^{\prime}_{N^{\prime}}&B^{\prime}_{N^{\prime}}\end{bmatrix}\right\\}}_{\eqqcolon\mathcal{S}^{\prime}}\Big{)}.$

Every convex combination of the matrices $B_{i}^{\prime}$ can be written as a
convex combination of the original $B_{i}$ matrices, and every convex
combination of the $B_{i}$’s is full column rank, which implies that every
convex combination of the $B_{i}^{\prime}$ is also full column rank. The
vertices, edges and faces of $\mathcal{P}$ are obtained by intersecting the
previous vertices, edges and faces from $\operatorname{conv}(\mathcal{S})$
with the set $\mathcal{T}$. Since $\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$ is
not on a vertex, edge or face of $\operatorname{conv}(\mathcal{S})$ it cannot
be on a vertex, edge or face of $\mathcal{P}$ and must be in the interior.
Hence, the set $\mathcal{S}^{\prime}$ satisfies Assumption 2, and also every
element of $\mathcal{S}^{\prime}$ satisfies (9) and (10).

∎

###### Remark 2

The results provided in Proposition 1 deal with systems that do not have the
same number of inputs as states; the special case of $n=m$ is solved
trivially, i.e. we have $\mathcal{S}=\mathcal{S}^{\prime}$. Since every convex
combination of $B_{i}$ is full column rank, that means that $B_{i}$ can be
inverted for all $i\in\left\\{1,\cdots,N\right\\}$, and the gains can be
obtained as

$\displaystyle K_{i}$ $\displaystyle=B_{i}^{-1}(A_{r}-A_{i})$ $\displaystyle
L_{i}$ $\displaystyle=B_{i}^{-1}B_{r}.$

$\triangle$

###  III-B Obtaining the sets $\mathcal{S}$ and $\mathcal{S}^{\prime}$

The modeling uncertainty in (2) is taken to be parametric uncertainty in the
matrices in this state-space system model. There is an implicit assumption
here that the state-space system model (2) is derived from physical laws
rather than from a state-space realization of an input-output system model. If
we consider the maximum range of values each entry of $A_{p}$, and $B_{p}$ can
take we get that we can bound $a_{ij}\in[a_{ij\min},a_{ij\max}]$ and
$b_{ij}\in[b_{ij\min},b_{ij\max}]$. When we consider the minimums and maximums
of every entry we write the following matrices

$\displaystyle A_{p\min}\leq A_{p}\leq A_{p\max},$ $\displaystyle
B_{p\min}\leq B_{p}\leq B_{p\max},$

where $A_{p\min}$ is a matrix where each entry $a_{ij}$ takes its minimum
value, and $A_{p\max}$ is a matrix where each entry $a_{ij}$ takes its maximum
value, $B_{p\min}$, and $B_{p\max}$ are matrices that take the minimums and
maximums of each entry $b_{ij}$, respectively, and the inequality is
considered entry-wise. Let $\mathcal{S}$ be the set of all possible system
matrices $\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$ each entry of which is
either the corresponding entry of $[A_{p\min}~{}~{}B_{p\min}]$ or the
corresponding entry of $[A_{p\max}~{}~{}B_{p\max}]$. Note that, by
construction, $\mathcal{S}$ satisfies Assumption 2 (i). The number of elements
in $\mathcal{S}$ is $N=2^{n(n+m)}$. If the entries can be parameterized by an
uncertainty vector $\eta\in\mathbb{R}^{q}$, then $\mathcal{S}$ can be defined
in terms of $2^{q}$ elements (see Remark 1).

The naïve corner model selection process described above guarantees that
$\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}\in\operatorname{int}\left(\operatorname{conv}\left(S\right)\right)$.
However, it is not evident that Assumption 2 (ii) is satisfied. The results
from [27] can be used to verify this condition. If Assumption 2 (ii) is
satisfied, then the last step is to use constructive procedure from the proof
of Proposition 1 to obtain $\mathcal{S}^{\prime}$.

###  III-C Example

In this example we consider a system with two states and one input. We assume
that the plant only has uncertainty in the $B_{p}$ matrix. The _unknown_ input
matrix to the system is $B_{p}=\begin{bmatrix}2\\\ 2\end{bmatrix}$. The
reference models input’s matrix is $B_{r}=\begin{bmatrix}10\\\
10\end{bmatrix}$. Note that by taking $L^{*}=5$ we can satisfy (5), and
Assumption 1 is satisfied. The polytopic uncertainty for $B_{p}$ is given by

$B_{p\min}=\begin{bmatrix}1\\\ 1\end{bmatrix},B_{p\max}=\begin{bmatrix}4\\\
5\end{bmatrix}.$

Using the selection procedure from Section III-B we can take all possible
combinations of the minimums and maximums of every entry of $B_{p\min}$ and
$B_{p\max}$ to get

$\displaystyle B_{1}=\begin{bmatrix}1\\\
1\end{bmatrix},B_{2}=\begin{bmatrix}1\\\
5\end{bmatrix},B_{3}=\begin{bmatrix}4\\\
5\end{bmatrix},B_{4}=\begin{bmatrix}4\\\ 1\end{bmatrix}.$

Since $A_{p}$ is completely known, we need to redefine the sets $\mathcal{S}$
and $\mathcal{T}$ as

$\displaystyle\mathcal{S}$
$\displaystyle=\left\\{B_{1},\cdots,B_{4}\right\\},$
$\displaystyle\mathcal{T}$ $\displaystyle=\big{\\{}B_{x}\in\mathbb{R}^{2\times
1}:B_{x}L=B_{r},L\in\mathbb{R}\big{\\}}.$

It is easily verifiable that
$B_{p}\in\operatorname{int}\left(\operatorname{conv}\left(\mathcal{S}\right)\right)$,
and that any convex combination of $B_{i}$ is full column rank. Nevertheless,
if we consider $B_{2}$, $B_{3}$, or $B_{4}$ we cannot satisfy (10). Using the
proof of Proposition 1 we solve for the set $\mathcal{P}$ graphically (closed
line segment going from $B_{1}$ to $B_{2}^{\prime}$), as shown in Figure 1 to
obtain

$\displaystyle\mathcal{P}=\operatorname{conv}(\mathcal{S})\cap\mathcal{T}=\operatorname{conv}\Bigg{(}\underbrace{\left\\{\begin{bmatrix}1\\\
1\end{bmatrix},\begin{bmatrix}4.5\\\
4.5\end{bmatrix}\right\\}}_{\mathcal{S}^{\prime}}\Bigg{)}.$

Note that
$\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}\in\operatorname{int}\left(\mathcal{P}\right)$.
In addition, every convex combination of $B_{1}$ and $B_{2}^{\prime}$ is full
column rank, and we can satisfy (10) with $L_{1}=10$, and
$L_{2^{\prime}}=20/9$.

Figure 1: Graphic illustration of the projection process to select corner
models.

##  IV Online Parameter Identification

Our multiple-model reference adaptive control (MMRAC) design will utilize a
blending-based multiple model parameter identification (MMPI) scheme that
generates estimates
$\begin{bmatrix}\hat{A}_{p}(t)&\hat{B}_{p}(t)\end{bmatrix}$ of
$\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$, as a weighted sum of the corner
models. In this section, we provide the design, stability and convergence
analysis of this MMPI scheme.

###  IV-A Multiple-Model Parameter Identifier Design

Filtering both sides of (2) by the linear filter $\dfrac{1}{s+\lambda}$, where
$\lambda>0$ is a design parameter, we obtain the parametric model

$\displaystyle z\left(t\right)$
$\displaystyle\stackrel{{\scriptstyle\lambda}}{{=}}\Theta_{p}\begin{bmatrix}\phi_{1}\left(t\right)\\\
\phi_{2}\left(t\right)\end{bmatrix}\eqqcolon\Theta_{p}\Phi(t),$ (11)

where $\phi_{1}\left(t\right),z\left(t\right)\in\mathbb{R}^{n}$,
$\phi_{2}\left(t\right)\in\mathbb{R}^{m}$ are generated by the filters

$\displaystyle\begin{bmatrix}\dot{\phi}_{1}(t)\\\
\dot{\phi}_{2}(t)\end{bmatrix}$
$\displaystyle=-\lambda\begin{bmatrix}\phi_{1}(t)\\\
\phi_{2}(t)\end{bmatrix}+\begin{bmatrix}x_{p}(t)\\\
u(t)\end{bmatrix},\>\begin{bmatrix}\phi_{1}(0)\\\
\phi_{2}(0)\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix},$ $\displaystyle
z(t)$ $\displaystyle=\dot{\phi}_{1}(t)=-\lambda\phi_{1}(t)+x_{p}(t).$ (12)

The relation (11) can be verified by taking time derivatives of both sides and
taking the difference, i.e., (12) implies for
$e_{z}(t)=z(t)-\Theta_{p}\Phi(t)$ that

$\displaystyle\dot{e}_{z}(t)$
$\displaystyle=\dot{z}(t)-\Theta_{p}\dot{\Phi}(t)$
$\displaystyle=-\lambda\dot{\phi}_{1}(t)+\dot{x}_{p}-\Theta_{p}\begin{bmatrix}\dot{\phi}_{1}(t)\\\
\dot{\phi}_{2}(t)\end{bmatrix}$
$\displaystyle=-\lambda\dot{\phi}_{1}(t)+\Theta_{p}\begin{bmatrix}x_{p}(t)\\\
u(t)\end{bmatrix}-\Theta_{p}\left(-\lambda\begin{bmatrix}\phi_{1}(t)\\\
\phi_{2}(t)\end{bmatrix}+\begin{bmatrix}x_{p}(t)\\\ u(t)\end{bmatrix}\right)$
$\displaystyle=-\lambda\left(-\lambda\phi_{1}(t)+x_{p}(t)\right)+\lambda\Theta_{p}\Phi(t)$
$\displaystyle=-\lambda\left(-\lambda\phi_{1}(t)+x_{p}(t)-\Theta_{p}\Phi(t)\right)$
$\displaystyle=-\lambda e_{z}(t),$ (13)

i.e., $e_{z}(t)$ is an exponentially decaying signal.

For each of the fixed models, define

$\displaystyle z_{i}\left(t\right)$
$\displaystyle\coloneqq\Theta_{i}\begin{bmatrix}\phi_{1}\left(t\right)\\\
\phi_{2}\left(t\right)\end{bmatrix}=\Theta_{i}\Phi(t),\>\forall
i\in\left\\{1,\cdots,N\right\\}.$ (14)

Then, the filtered state estimation error for each of the $N$ fixed models is
defined as

$\varepsilon_{i}\left(t\right)\coloneqq\dfrac{z\left(t\right)-z_{i}\left(t\right)}{m_{s}^{2}(t)},\forall
i\in\left\\{1,\cdots,N\right\\},$ (15)

where $m_{s}^{2}(t)\coloneqq 1+\alpha\|\Phi(t)\|^{2}$, $\alpha>0$, is a
normalization signal which guarantees that $\dfrac{\Phi(t)}{m_{s}(t)}$ is
bounded. By (7), (11)–(14), for any
$w=\left[w_{1},\dots,w_{N}\right]^{\top}\in\mathcal{W}$, we obtain

$\displaystyle z\left(t\right)$
$\displaystyle=\sum_{i=1}^{N}w_{i}z_{i}\left(t\right)+e_{z}(t).$ (16)

Using $\sum_{i=1}^{N}w_{i}=1$, (16) further implies

$\displaystyle
e_{z}(t)=\sum_{i=1}^{N}w_{i}z\left(t\right)-\sum_{i=1}^{N}w_{i}z_{i}\left(t\right)$
$\displaystyle=\sum_{i=1}^{N}w_{i}\left(z\left(t\right)-z_{i}\left(t\right)\right)$
$\displaystyle=\sum_{i=1}^{N}w_{i}\varepsilon_{i}\left(t\right)m_{s}^{2}(t).$
(17)

Adding $-\varepsilon_{N}\left(t\right)m_{s}^{2}(t)$ to both sides of (17), we
obtain

$\displaystyle e_{z}(t)-\varepsilon_{N}\left(t\right)m_{s}^{2}(t)$
$\displaystyle=\sum_{i=1}^{N-1}w_{i}\left(\varepsilon_{i}\left(t\right)-\varepsilon_{N}\left(t\right)\right)m_{s}^{2}(t),$

which implies, for any
$w=\left[w_{1},\dots,w_{N}\right]^{\top}\in\mathcal{W}$, that

$e_{z}(t)-\varepsilon_{N}\left(t\right)=\sum_{i=1}^{N-1}w_{i}\left(\varepsilon_{i}\left(t\right)-\varepsilon_{N}\left(t\right)\right).$
(18)

Defining the $n\times(N-1)$ time-varying matrix

$E\left(t\right)\coloneqq\begin{bmatrix}\varepsilon_{1}\left(t\right)-\varepsilon_{N}\left(t\right)&\cdots&\varepsilon_{N-1}\left(t\right)-\varepsilon_{N}\left(t\right)\end{bmatrix},$
(19)

we can rewrite (18) in matrix form as

$E\left(t\right)\bar{w}=e_{z}(t)-\varepsilon_{N}\left(t\right)$ (20)

for any
$\left[\bar{w}^{\top},w_{N}\right]^{\top}=\left[w_{1},\dots,w_{N-1},w_{N}\right]^{\top}\in\mathcal{W}$.
Equation 20 motivates using the following recursive adaptive law [22] to
generate estimate
$\hat{w}(t)=[\hat{\bar{w}}^{\top}\left(t\right),\hat{w}_{N}\left(t\right)]^{\top}\in\mathbb{R}^{N}$,
such that
$\lim_{t\rightarrow\infty}\sum_{i=1}^{N}\hat{w}_{i}(t)\begin{bmatrix}A_{i}&B_{i}\end{bmatrix}=\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$.

$\displaystyle\begin{split}\dot{\hat{\bar{w}}}\left(t\right)&=-\Gamma
E^{\top}\left(t\right)E\left(t\right)\hat{\bar{w}}\left(t\right)-\Gamma
E^{\top}\left(t\right)\varepsilon_{N}\left(t\right),\\\
\hat{w}_{N}\left(t\right)&=1-\sum_{i=1}^{N-1}\hat{w}_{i}\left(t\right),\end{split}$
(21)

where the tuning parameter $\Gamma\in\mathbb{R}^{(N-1)\times(N-1)}$ is a
symmetric positive definite matrix. Let $w^{*}$ be some element in
$\mathcal{W}$, and based on (21), the estimation error
$\tilde{\bar{w}}\left(t\right)\coloneqq\hat{\bar{w}}\left(t\right)-\bar{w}^{*}$
satisfies

$\displaystyle\dot{\tilde{\bar{w}}}\left(t\right)$
$\displaystyle=\dot{\hat{\bar{w}}}\left(t\right)-\dot{\bar{w}}^{*}=\dot{\hat{\bar{w}}}\left(t\right)$
$\displaystyle=-\Gamma
E^{\top}\left(t\right)E\left(t\right)\hat{\bar{w}}\left(t\right)-\Gamma
E^{\top}\left(t\right)\varepsilon_{N}\left(t\right).$ (22)

Pre-multiplying (20) by $\Gamma E^{\top}(t)$ and moving all the terms to the
left hand side, we obtain

$\Gamma E^{\top}(t)E\left(t\right)\bar{w}+\Gamma
E^{\top}(t)\varepsilon_{N}\left(t\right)-\Gamma E^{\top}(t)e_{z}(t)=0$ (23)

Adding (23) to (22), we further obtain

$\displaystyle\dot{\tilde{\bar{w}}}\left(t\right)$ $\displaystyle=-\Gamma
E^{\top}\left(t\right)E\left(t\right)\tilde{\bar{w}}\left(t\right)-\Gamma
E^{\top}(t)e_{z}(t).$ (24)

###  IV-B Stability and Convergence of the Identification Process

In this subsection, we establish the stability and convergence properties of
the estimation scheme (21) utilizing the following lemma.

######  Lemma 1 (See [22], Barbalat’s Lemma)

For a function $f:[0,\infty)\rightarrow\mathbb{R}$, if
$f,\>\dot{f}\in\mathcal{L}_{\infty}$ and ${f\in\mathcal{L}_{p}}$, for some
$p\in\left[1,\infty\right)$, then
$\lim\limits_{t\rightarrow\infty}f\left(t\right)=0$.

The main stability and convergence properties are established in the following
theorems and lemmas below.

###### Theorem 1

Consider the system (2) with definitions (15), and (19). Let Assumption 2
hold, $w^{*}=\left[\bar{w}^{*\top},w^{*}_{N}\right]^{\top}\in\mathcal{W}$ be
an arbitrary vector within the set defined in (7), and
$\tilde{\bar{w}}(t)\coloneqq\hat{\bar{w}}(t)-\bar{w}^{*}$. For any initial
condition $\hat{\bar{w}}(0)\in\mathbb{R}^{N-1}$, the estimation scheme (21)
guarantees that:

  1. (i)

$\hat{\bar{w}}$, $\tilde{\bar{w}}$, $\dot{\tilde{\bar{w}}}$ and $E$ are
bounded signals.

  2. (ii)

$E\tilde{\bar{w}}$ and $\dot{\tilde{\bar{w}}}$ are square integrable.

  3. (iii)

$\lim_{t\rightarrow\infty}E(t)\tilde{\bar{w}}(t)=0$.

  4. (iv)

$\hat{\bar{w}}(t)$ asymptotically converges to a constant vector
$\bar{\bar{w}}\in\mathbb{R}^{N-1}$.

###### Proof:

Consider the Lyapunov-like function

$V_{1}(\tilde{\bar{w}}(t),e_{z}(t))=\frac{1}{2}\tilde{\bar{w}}^{\top}(t)\Gamma^{-1}\tilde{\bar{w}}(t)+\frac{1}{2\lambda}e_{z}^{\top}(t)e_{z}(t).$
(25)

Taking the time derivative of (25) along (24), we have

$\displaystyle\dfrac{\mathrm{d}V_{1}(t)}{\mathrm{d}t}$
$\displaystyle=\tilde{\bar{w}}^{\top}(t)\Gamma^{-1}\dot{\tilde{\bar{w}}}(t)+\dfrac{1}{\lambda}e_{z}^{\top}(t)\dot{e}_{z}(t)$
$\displaystyle=-\tilde{\bar{w}}^{\top}(t)\Gamma^{-1}\Gamma
E^{\top}(t)E(t)\tilde{\bar{w}}(t)$
$\displaystyle-\tilde{\bar{w}}^{\top}(t)\Gamma^{-1}\Gamma
E^{\top}(t)e_{z}(t)+\dfrac{1}{\lambda}e_{z}^{\top}\dot{e}_{z}$
$\displaystyle=-\tilde{\bar{w}}^{\top}(t)E^{\top}(t)E(t)\tilde{\bar{w}}(t)-\tilde{\bar{w}}^{\top}(t)E^{\top}(t)e_{z}(t)+\dfrac{1}{\lambda}e_{z}^{\top}\dot{e}_{z}$
$\displaystyle=-\dfrac{1}{2}\|E\tilde{\bar{w}}\|^{2}+\dfrac{1}{2}\|e_{z}(t)\|^{2}-\dfrac{1}{2}\|E(t)\tilde{\bar{w}}(t)+e_{z}(t)\|^{2}.$
(26)

Since $e_{z}(t)$ is an exponentially decaying signal, Equation 26 implies that
$\dot{V}_{1}$ is bounded and there exists a time instant $t_{1}>0$ such that,
for all $t\geq t_{1}$, $\dot{V}_{1}(t)\leq 0$. Hence. $V_{1}$,
$\tilde{\bar{w}}$, and $\hat{\bar{w}}$ are bounded. Because of normalization
(15), $\varepsilon_{i}$ terms are bounded and hence $E$ is bounded, which
together with (24) also implies that $\dot{\tilde{\bar{w}}}$ is bounded,
finishing the proof of (i).

Since $V_{1}(t)$ is always positive, bounded, and decaying, the integral of
(26) for $t=0$ to $\infty$ is finite, and hence $E\tilde{\bar{w}}$ and
$E\tilde{\bar{w}}+e_{z}(t)$ are square integrable. Since $E$ is bounded, this,
together with (24), further implies that $\dot{\tilde{\bar{w}}}$ is square
integrable, completing the proof of (ii).

Items (i) and (ii) together with Barbalat’s Lemma imply (iii). Items (i) and
(ii) further imply that $\lim_{t\rightarrow\infty}\tilde{\bar{w}}(t)$ and,
hence, $\lim_{t\rightarrow\infty}\hat{\bar{w}}(t)$ exist and are finite,
proving (iv). ∎

Theorem 1 implies that $\|E(t)\tilde{\bar{w}}(t)\|$ asymptotically converges
to zero, but this does not mean that $\tilde{\bar{w}}(t)$ converges to the set
$\mathcal{W}$. We can now state, and prove, the main result of the
identification process.

###### Theorem 2

Consider the system (2) with definitions (15), and (19), and the estimation
scheme (21). Let Assumption 2 hold,
$w^{*}=\left[\bar{w}^{*\top},w^{*}_{N}\right]^{\top}\in\mathcal{W}$ be an
arbitrary vector within the set defined in (7),
$\tilde{\bar{w}}(t)\coloneqq\hat{\bar{w}}(t)-\bar{w}^{*}$, and $\Phi$ be
bounded and satisfy the PE condition (1). Then, for any initial condition
$\hat{\bar{w}}(0)\in\mathbb{R}^{N-1}$, the estimated system matrix
$\sum_{i=1}^{N}\hat{w}_{i}(t)\Theta_{i}$ asymptotically converges to
$\Theta_{p}$.

###### Proof:

Let
$Q(t)\coloneqq\sum_{i=1}^{N-1}(\tilde{w}_{i}(t)(\Theta_{N}-\Theta_{i}))\in\mathbb{R}^{n\times(n+m)}$.
It is easily verifiable that if $Q(t)=0$, then
$\sum_{i=1}^{N}\hat{w}_{i}(t)\Theta_{i}=\Theta_{p}$, and
$\hat{w}(t)\in\mathcal{W}$. Theorem 1 (iv) implies that $Q(t)$ asymptotically
converges to a constant matrix $\bar{Q}\in\mathbb{R}^{n\times(n+m)}$. Hence,
to establish that $\sum_{i=1}^{N}\hat{w}_{i}(t)\Theta_{i}$ asymptotically
converges to $\Theta_{p}$, we will show that $\bar{Q}=0$.

Expressing Theorem 1 (iii) in summation form, we get

$\displaystyle\lim_{t\rightarrow\infty}E(t)\tilde{\bar{w}}(t)$
$\displaystyle=\lim_{t\rightarrow\infty}Q(t)\dfrac{\Phi(t)}{1+\alpha\|\Phi(t)\|^{2}}$
$\displaystyle=\lim_{t\rightarrow\infty}\bar{Q}\dfrac{\Phi(t)}{1+\alpha\|\Phi(t)\|^{2}}=0.$
(27)

Since $\Phi$ is assumed to be bounded, (IV-B) implies that

$\displaystyle\lim_{t\rightarrow\infty}\bar{Q}\Phi(t)=0.$

and hence

$\displaystyle\lim_{t\rightarrow\infty}\bar{Q}\Phi(t)\Phi^{\top}(t)\bar{Q}^{\top}=0.$
(28)

Since $\Phi$ satisfies the PE condition (1), Equation (IV-B) implies that
$\bar{Q}=0$, completing the proof.

∎

The recursive adaptive law (21) guarantees that
$\sum_{i=1}^{N}\hat{w}_{i}\left(t\right)=1$, but not that
$\hat{w}\left(t\right)\in[0,1]^{N}$. Since the set $\left[0,1\right]^{N}$ is
convex, the projection of $\hat{w}\left(t\right)$ into $\left[0,1\right]^{N}$
is well-defined.

Define the compact set

$\Pi\coloneqq\left\\{\hat{\bar{w}}\in\left[0,1\right]^{N-1}:\sum_{i=1}^{N-1}\hat{\bar{w}}_{i}\leq
1\right\\},$ (29)

and let
$\operatorname{Pr}_{\Pi,\hat{\bar{w}}}:\mathbb{R}^{N-1}\rightarrow\Pi\subset\mathbb{R}^{N-1}$
denote the parameter projection operator [22]. Choose
$\hat{w}(0)\in\operatorname{int}(\Pi)$, then the recursive adaptive algorithm
(21) with parameter projection is as follows:

$\displaystyle\dot{\hat{\bar{w}}}\left(t\right)$
$\displaystyle=\operatorname{Pr}_{\Pi,\hat{\bar{w}}}\left(-\Gamma\left(E^{\top}\left(t\right)E\left(t\right)\hat{\bar{w}}\left(t\right)+E^{\top}\left(t\right)\varepsilon_{N}\left(t\right)\right)\right),$
$\displaystyle\hat{w}_{N}\left(t\right)$
$\displaystyle=1-\sum_{i=1}^{N-1}\hat{w}_{i}\left(t\right).$ (30)

The parameter projection operator enforces that $\Pi$ is a positively
invariant subset for the dynamics (30).

###### Corollary 1

The gradient adaptive law with parameter projection (30) has all the
properties established in Theorem 1 for the gradient adaptive law (21).
Furthermore, we get the following properties:

  1. (i)

If $\hat{w}\left(0\right)\in\Pi$, then $\hat{w}\left(t\right)\in\Pi$, $\forall
t\geq 0$.

  2. (ii)

If $\hat{w}\left(0\right)\in\operatorname{int}\Pi$, then $\hat{w}(t)$
asymptotically converges to the set $\mathcal{W}$.

###### Proof:

The proof follows applying Theorem 3.10.1 of [22] to (21) combined with
Theorems 1 and 2. ∎

##  V Multiple Model Reference Adaptive Control

In this section we combine the system identification scheme of Section IV with
a MMRAC controller to achieve asymptotic tracking of the reference system (3).

###  V-A Multiple Model Reference Adaptive Control Design

We use the recursive adaptive algorithm with parameter projection (30) to
design a MMRAC scheme to achieve the adaptive state tracking control task
stated in Section II. We utilize Theorem 1 (iii) to construct the proposed
MMRAC scheme.

If $w\in\mathcal{W}$, then multiplying both sides of (10) by $w_{i}$ and
summing over $i$ yields

$\displaystyle\sum_{i=1}^{N}w_{i}B_{i}L_{i}$
$\displaystyle=\sum_{i=1}^{N}w_{i}B_{r}=B_{r},$

which implies, together with (5) from Assumption 1, that

$B_{p}L^{*}=\sum_{i=1}^{N}w_{i}B_{i}L_{i}.$ (31)

Applying the same steps on (9) we get

$\displaystyle\begin{split}A_{r}=\sum_{i=1}^{N}w_{i}A_{r}&=\sum_{i=1}^{N}w_{i}\big{(}A_{i}+B_{i}K_{i}\big{)}\\\
&=A_{p}+\sum_{i=1}^{N}w_{i}B_{i}K_{i}.\end{split}$ (32)

Comparing (32) to (4) we get that

$B_{p}K^{*}=\sum_{i=1}^{N}w_{i}B_{i}K_{i}.$ (33)

Equations (31) and (33) motive us to generate estimates of the gains $K^{*}$
and $L^{*}$ using the estimates $\hat{w}\left(t\right)$, keeping in mind the
rank supposition in Assumption 2, as

$\displaystyle\hat{K}\left(t\right)$
$\displaystyle=\hat{B}_{p}^{\dagger}(t)\sum_{i=1}^{N}\hat{w}_{i}\left(t\right)B_{i}K_{i},$
(34) $\displaystyle\hat{L}\left(t\right)$
$\displaystyle=\hat{B}_{p}^{\dagger}(t)\sum_{i=1}^{N}\hat{w}_{i}\left(t\right)B_{i}L_{i},$
(35)

where $\hat{B}_{p}^{\dagger}(t)$ is the Moore-Penrose inverse of

$\hat{B}_{p}(t)=\sum_{i=1}^{N}\hat{w}_{i}\left(t\right)B_{i}.$ (36)

Equations 33, 31, 34 and 35, together with Corollary 1 further imply the
following:

###### Corollary 2

Consider the system (2) with definitions (15), (19), and the reference model
(3). If $\Phi(t)$ satisfies the PE condition with constants $\alpha_{\Phi 1}$,
$\alpha_{\Phi 2}$, and $T_{\Phi}$, then the estimates $\hat{K}\left(t\right)$
and $\hat{L}\left(t\right)$ defined in (34), and (35) are bounded, and
asymptotically converge to $K^{*}$ and $L^{*}$, respectively.

###### Proof:

Since $\Pi$ is compact, we get that $\hat{B}_{p}(t)$ belongs to a compact set,
which means the Moore-Penrose inverse of $\hat{B}_{p}(t)$ exists, and it is
bounded. Furthermore, from Assumption 2, we have that $\hat{B}_{p}(t)$ is full
column rank, which means that we can write

$\hat{B}_{p}^{\dagger}(t)=\left(\hat{B}_{p}^{\top}(t)\hat{B}_{p}(t)\right)^{-1}\hat{B}_{p}^{\top}(t).$
(37)

Combining (37) with Equations 34 and 35 we get that $\hat{K}\left(t\right)$,
and $\hat{L}\left(t\right)$ are bounded. The proof to show asymptotic
convergence is the same as the proof for Theorem 1 (iii). ∎

The control law we will consider to achieve asymptotic tracking of the
reference model is

$u\left(t\right)\coloneqq\hat{K}\left(t\right)x_{p}\left(t\right)+\hat{L}\left(t\right)r\left(t\right).$
(38)

###  V-B Stability Analysis of the Entire System

The main result of the paper is now presented.

###### Theorem 3

Consider the plant 2 and the reference model 3. If $\Phi(t)$ satisfies the PE
condition (1), and Assumptions 2 and 1 hold, then the MMRAC scheme 12, 15, 14,
30, 34, 35, 36, 37 and 38 guarantees that for any

  1. (i)

initial conditions of the plant (2),

  2. (ii)

initial conditions of the reference model (3), and

  3. (iii)

piecewise continuous and bounded reference signal
${r:\left[0,\infty\right)\rightarrow\mathbb{R}^{m}}$ in (3),

all closed-loop signals are bounded and $x_{p}\left(t\right)$ asymptotically
converges to $x_{r}\left(t\right)$.

###### Proof:

Let $x_{p}(0)=x_{p0}$ and $x_{r}(0)=x_{r0}$ be arbitrary initial plant and
reference model states, and $r(t)$ be any known, bounded, piecewise continuous
reference signal. Substituting (38) into (2) we get

$\displaystyle\dot{x}_{p}\left(t\right)$
$\displaystyle=A_{p}x_{p}\left(t\right)+B_{p}u\left(t\right)$
$\displaystyle=A_{p}x_{p}\left(t\right)+B_{p}\hat{K}\left(t\right)x_{p}\left(t\right)+B_{p}\hat{L}\left(t\right)r\left(t\right).$

Adding and subtracting $B_{p}K^{*}x_{p}\left(t\right)$ and
$B_{p}L^{*}r\left(t\right)$, defining
$\tilde{K}\left(t\right)\coloneqq\hat{K}\left(t\right)-K^{*}$,
$\tilde{L}\left(t\right)\coloneqq\hat{L}\left(t\right)-L^{*}$, and using
Assumption 1 we get

$\displaystyle\dot{x}_{p}\left(t\right)=A_{r}x_{p}\left(t\right)+B_{r}r\left(t\right)+B_{p}\tilde{K}\left(t\right)x_{p}\left(t\right)+B_{p}\tilde{L}\left(t\right)r\left(t\right).$

For the tracking error $e\left(t\right)\coloneqq
x_{p}\left(t\right)-x_{r}\left(t\right)$, this implies

$\displaystyle\dot{e}\left(t\right)$
$\displaystyle=A_{r}e\left(t\right)+B_{p}\tilde{K}\left(t\right)x_{p}\left(t\right)+B_{p}\tilde{L}\left(t\right)r\left(t\right)$
$\displaystyle=\left(A_{r}+B_{p}\tilde{K}(t)\right)e(t)+B_{p}\tilde{K}(t)x_{r}(t)+B_{p}\tilde{L}(t)r(t).$
(39)

Let $Q\in\mathbb{R}^{n\times n}$ be a fixed, symmetric, and positive definite
matrix. Then we can define $P\in\mathbb{R}^{n\times n}$ to be the unique
positive definite and symmetric solution of ${PA_{r}+A_{r}^{\top}P+Q=0}$.
Consider the positive definite function

$\displaystyle V\left(e\left(t\right)\right)$
$\displaystyle=e^{\top}\left(t\right)Pe\left(t\right).$ (40)

Taking the derivative of (40) along solutions of (39) we get

$\displaystyle\dot{V}\left(e\left(t\right)\right)$
$\displaystyle=2e^{\top}\left(t\right)P\dot{e}\left(t\right)$
$\displaystyle=2e^{\top}\left(t\right)PA_{r}e\left(t\right)+2e^{\top}\left(t\right)PB_{p}\tilde{K}\left(t\right)x_{p}\left(t\right)$
$\displaystyle\qquad+2e^{\top}\left(t\right)PB_{p}\tilde{L}\left(t\right)r\left(t\right).$
(41)

The first term satisfies

$\displaystyle 2e^{\top}\left(t\right)PA_{r}e\left(t\right)$
$\displaystyle=e^{\top}\left(t\right)PA_{r}e\left(t\right)+e^{\top}\left(t\right)PA_{r}e\left(t\right)$
$\displaystyle=e^{\top}\left(t\right)\left(PA_{r}+A_{r}^{\top}P\right)e\left(t\right)$
$\displaystyle=-e^{\top}\left(t\right)Qe\left(t\right).$

Hence, (41) can be rewritten as as

$\displaystyle\dot{V}\left(e\left(t\right)\right)=$
$\displaystyle-e^{\top}\left(t\right)Qe\left(t\right)+2e^{\top}\left(t\right)PB_{p}\tilde{K}\left(t\right)x_{p}\left(t\right)$
$\displaystyle\qquad+2e^{\top}\left(t\right)PB_{p}\tilde{L}\left(t\right)r\left(t\right).$

Substituting $x_{p}(t)=e(t)+x_{r}(t)$, we get

$\displaystyle\dot{V}\left(e\left(t\right)\right)=$
$\displaystyle-e^{\top}\left(t\right)\left(Q-PB_{p}\tilde{K}\left(t\right)\right)e\left(t\right)$
$\displaystyle+2e^{\top}\left(t\right)PB_{p}\tilde{K}\left(t\right)x_{r}\left(t\right)$
$\displaystyle+2e^{\top}\left(t\right)PB_{p}\tilde{L}\left(t\right)r\left(t\right).$

Defining

$c_{1}(t)\coloneqq\|PB_{p}\tilde{K}(t)\|,$

and

$c_{2}(t)\coloneqq
2\|PB_{p}\tilde{K}(t)\|\|x_{r}(t)\|+2\|PB_{p}\tilde{L}(t)\|\|r(t)\|,$

we have

$\dot{V}(e(t))\leq-\left(\lambda_{\textrm{min}}(Q)-c_{1}(t)\right)\|e(t)\|^{2}+c_{2}(t)\|e(t)\|,$
(42)

where $\lambda_{\textrm{min}}(Q)>0$ is the minimum eigenvalue of $Q$. Note
that (39) does not have a finite escape time; from Corollary 2 we have that
$\hat{K}(t)$ and $\hat{L}(t)$ are continuous and bounded, which means that
$\tilde{K}(t)$ and $\tilde{L}(t)$ are also continuous and bounded, hence (39)
may only go to infinity as time goes to infinity. Moreover, if $\Phi(t)$ is
PE, then $\tilde{K}(t)$ and $\tilde{L}(t)$ converge to zero asymptotically,
which lets us conclude that $c_{1}(t)$ is bounded, and there exists $t_{1}\geq
0$ such that $\lambda_{\textrm{min}}(Q)>c_{1}(t)$ for all $t\geq t_{1}$. We
can define

$\displaystyle\bar{c}_{1}=\sup\\{c_{1}(t),t\geq t_{1}\\},$ (43)
$\displaystyle\bar{c}_{2}=\sup\\{c_{2}(t),t\geq t_{1}\\}.$ (44)

This implies that for all $t\geq t_{1}$ we get that if

$\displaystyle\|e(t)\|>\dfrac{\bar{c}_{2}}{\lambda_{\min}Q-\bar{c}_{1}},$

then $\dot{V}<0$, and we can conclude that $e$ is bounded. This further
implies that $x_{p}$ is bounded, which finally implies that $u$ is bounded,
showing that all signals in the closed-loop system are bounded. Combining this
with Corollary 2 we get

$\displaystyle\lim_{t\rightarrow\infty}\left(\dot{V}(e(t))+e^{\top}\left(t\right)Qe\left(t\right)\right)=0.$

This implies that $e(t)$ converges to $0$, asymptotically, i.e., $x_{p}(t)$
asymptotically converges to $x_{r}(t)$. ∎

##  VI Simulations

Figure 2: State of the system $x_{p}(t)$, and state of the reference model
$x_{r}(t)$. Figure 3: Control efforts for MMRAC and a single model MRAC.
Figure 4: Euclidean norm of the tracking error for MMRAC and a single model
MRAC. Figure 5: Semilog plot for the norm 2 of the tracking errors, and the
linear regressions for MMRAC and a single model. Figure 6: Estimation of each
entry of $A_{p}$. The dashed line is the true value. Figure 7: Estimation of
each entry of $B_{p}$. The dashed line is the true value.

In this section, we illustrate the behavior and performance of the proposed
MMRAC scheme through a set of simulation tests performed on an uncertain model
(2), a reference model (3), and a set $\mathcal{S}$ (6) that defines the
polytopic uncertainty. We compare the results with simulations using a single
model for MRAC.

Consider the uncertain system with $n=3$, $m=2$ given by the matrices

$A_{p}=\begin{bmatrix}-4.725&-6.275&-2.175\\\ -0.925&-3.85&0.35\\\
-3.65&-8.125&-2.825\end{bmatrix},\>B_{p}=\begin{bmatrix}-0.575&-2.2\\\
-0.45&0.575\\\ -1.025&-1.625\end{bmatrix}.$

The reference model (3) is

$A_{r}=\begin{bmatrix}-1&0&0\\\ 0&-1&0\\\
1&1&-1\end{bmatrix},\>B_{r}=\begin{bmatrix}1&0\\\ 0&1\\\ 1&1\end{bmatrix}.$

It can be verified that Assumption 1 is satisfied by defining the matrices

$\displaystyle K^{*}=\begin{bmatrix}-3.16&-7.48&-0.36\\\
-0.87&-0.89&-0.89\end{bmatrix},L^{*}=\begin{bmatrix}-0.44&-1.66\\\
-0.34&0.44\end{bmatrix}.$

Considering the following $5$ matrix pairs

$\displaystyle A_{1}=\begin{bmatrix}-0.75&0.25&0.25\\\ -3.75&-4.75&-3.75\\\
-2.5&-2.5&-4.5\end{bmatrix},B_{1}=\begin{bmatrix}0.25&-0.5\\\ 1.25&2.5\\\
1.5&2\end{bmatrix},$ $\displaystyle A_{2}=\begin{bmatrix}-16&-30&-5\\\
-3.5&-8.5&-1.5\\\ -17.5&-36.5&-7.5\end{bmatrix},B_{2}=\begin{bmatrix}-5&-10\\\
-1&-2.5\\\ -6&-12.5\end{bmatrix},$ $\displaystyle
A_{3}=\begin{bmatrix}2&0&-1\\\ -0.5&-1.5&-0.5\\\
-0.5&0.5&-2.5\end{bmatrix},B_{3}=\begin{bmatrix}-1&1\\\ -0.5&0\\\
-1.5&1\end{bmatrix},$ $\displaystyle A_{4}=\begin{bmatrix}-1.5&-0.75&-0.75\\\
-0.25&-0.875&1.125\\\
0.25&0.375&-0.625\end{bmatrix},B_{4}=\begin{bmatrix}0.25&0.25\\\
0.125&-0.375\\\ 0.375&-0.125\end{bmatrix},$ $\displaystyle
A_{5}=\begin{bmatrix}-4&-1&-5\\\ 5&-2&8\\\
3&-1&2\end{bmatrix},B_{5}=\begin{bmatrix}2&-1\\\ -3&2\\\ -1&1\end{bmatrix},$

we can define the set

$\displaystyle\mathcal{S}=\left\\{\begin{bmatrix}A_{1}&B_{1}\end{bmatrix},\cdots,\begin{bmatrix}A_{5}&B_{5}\end{bmatrix}\right\\}.$

The set $\mathcal{S}$ satisfies Assumption 2. The design parameters for the
estimation scheme and the controller are $\Gamma=2I_{5\times 5}$,
$\lambda=0.5$, and $\alpha=0.01$, with the initial condition
$\hat{w}(0)=\begin{bmatrix}0.2&0.15&0.15&0.1&0.4\end{bmatrix}^{\top}$. The
input to the reference model is taken to be
$r_{1}\left(t\right)=r_{2}\left(t\right)=\sin{\left(t\right)}+0.5\sin{2t}$.
With these definitions the full controller is

$\displaystyle\begin{bmatrix}\dot{\phi}_{1}(t)\\\
\dot{\phi}_{2}(t)\end{bmatrix}$
$\displaystyle=-\lambda\begin{bmatrix}\phi_{1}(t)\\\
\phi_{2}(t)\end{bmatrix}+\begin{bmatrix}x_{p}(t)\\\ u(t)\end{bmatrix},$
$\displaystyle\dot{\hat{\bar{w}}}\left(t\right)$
$\displaystyle=\operatorname{Pr}_{\Pi,\hat{\bar{w}}}\left(-\Gamma\left(E^{\top}\left(t\right)E\left(t\right)\hat{\bar{w}}\left(t\right)+E^{\top}\left(t\right)\varepsilon_{5}\left(t\right)\right)\right),$
$\displaystyle\hat{w}_{5}\left(t\right)$
$\displaystyle=1-\sum_{i=1}^{4}\hat{w}_{i}\left(t\right),$
$\displaystyle\hat{K}\left(t\right)$
$\displaystyle=\hat{B}_{p}^{\dagger}(t)\sum_{i=1}^{5}\hat{w}_{i}\left(t\right)B_{i}K_{i},$
$\displaystyle\hat{L}\left(t\right)$
$\displaystyle=\hat{B}_{p}^{\dagger}(t)\sum_{i=1}^{5}\hat{w}_{i}\left(t\right)B_{i}L_{i},$
$\displaystyle u\left(t\right)$
$\displaystyle=\hat{K}\left(t\right)x_{p}\left(t\right)+\hat{L}\left(t\right)r\left(t\right).$

The simulations of the MMRAC are compared to a direct adaptive control
technique using a single model MRAC (see Chapter 9 of [23]). The simulations
consider the same initial condition for the system, the same reference, and
the initial gains are calculated as

$\displaystyle\hat{K}(0)=\hat{B}_{p}^{\dagger}(0)\sum_{i=0}^{5}w_{i}(0)K_{i},$
$\displaystyle\hat{L}(0)=\hat{B}_{p}^{\dagger}(0)\sum_{i=0}^{5}w_{i}(0)L_{i}.$

In Figure 2 we see the time evolution of each of the states for the reference
model, the states generated by the MMRAC, and the states generated by MRAC.
Convergence is achieved for both controllers, with a clear advantage of the
multiple-model technique while using approximately the same control effort, as
depicted in Figure 3. The convergence speed is a key factor to consider, and
there is a clear advantage of using multiple models when we analyze the
dynamics of the error. In Figure 4 we compare the time evolution of the norm
of the error for the case of MMRAC, and MRAC. To better visualize the
advantage in convergence speed, we present the norm of the error in log-scale
in Figure 5, as we perform a linear regression on each of the signals. We see
that the convergence speed of the MMRAC scheme is approximately two order of
magnitude faster than the single model approach (the slope of the linear
regressions are -0.0333 and -0.0103 for MMRAC and a single model,
respectively), matching with the conjecture in [29]. Finally, as mentioned in
Section IV, we achieve asymptotic convergence of the estimated system’s
matrices. In Figures 6 and 7 we see the asymptotic convergence of each entry
of the matrix pairs
$\begin{bmatrix}\hat{A}_{p}(t)&\hat{B}_{p}(t)\end{bmatrix}$ to
$\begin{bmatrix}A_{p}&B_{p}\end{bmatrix}$, respectively.

##  VII Conclusions

In this article we developed a multiple model reference adaptive control
(MMRAC) scheme for multi-input, linear, time-invariant systems with uncertain
parameters that lie inside a known compact and convex polytope. The controller
performs online parameter identification of the system matrices as a convex
combination of an arbitrary number of fixed model, one at each extreme point
of the convex polytope of uncertainty. The identification is proven to be
asymptotically stable, and sufficient conditions for perfect identification of
the uncertain system matrices are provided. The tracking controller of the
MMRAC guarantees that all close-loop signals are bounded, and that the
difference between the plant states and the ones generated by a linear
reference model asymptotically converge to zero. To verify the effectiveness
of the proposed MMRAC scheme we compare it to a single model MRAC scheme
through MATLAB and Simulink simulations. For the simulation examples, we
observe that the convergence speed is approximately two orders of magnitude
faster than using a single model, with similar control effort. Future lines of
research include generalizing the proposed scheme to linear, time-varying
systems and nonlinear systems.

## References

  * [1] K. S. Narendra and Z. Han, “A new approach to adaptive control using multiple models,” _International Journal of Adaptive Control and Signal Processing_ , vol. 26, no. 8, pp. 778–799, Aug. 2012.
  * [2] A. S. Morse, “A simple example of an adaptive control system,” _Communications in Information and Systems_ , vol. 11, no. 2, pp. 105–118, 2011.
  * [3] M. Kuipers and P. Ioannou, “Multiple model adaptive control with mixing,” _IEEE Transactions on Automatic Control_ , vol. 55, no. 8, pp. 1822–1836, Aug. 2010.
  * [4] J. P. Hespanha, D. Liberzon, and A. S. Morse, “Overcoming the limitations of adaptive control by means of logic-based switching,” _Systems & Control Letters_, vol. 49, no. 1, pp. 49–65, May 2003.
  * [5] G. M. Prasad, V. Kedia, and A. S. Rao, “Multi-model predictive control (mmpc) for non-linear systems with time delay: An experimental investigation,” in _IEEE International Conference on Measurement, Instrumentation, Control and Automation_ , 2020, pp. 1–5.
  * [6] L. Vu and D. Liberzon, “Supervisory control of uncertain linear time-varying systems,” _IEEE Transactions on Automatic Control_ , vol. 56, no. 1, pp. 27–42, 2011.
  * [7] J. P. Hespanha and A. S. Morse, “Switching between stabilizing controllers,” _Automatica_ , vol. 38, no. 11, pp. 1905–1917, Nov. 2002.
  * [8] J. P. Hespanha, D. Liberzon, A. S. Morse, B. D. O. Anderson, T. S. Brinsmead, and F. De Bruyne, “Multiple model adaptive control. Part 2: switching,” _International Journal of Robust and Nonlinear Control_ , vol. 11, no. 5, pp. 479–496, Apr. 2001.
  * [9] H. Zengin, N. Zengin, B. Fidan, and A. Khajepour, “Blending based multiple-model adaptive control of multivariable systems with application to lateral vehicle motion control,” _European Journal of Control_ , vol. 58, pp. 1–10, Mar. 2021.
  * [10] K. Büyükkabasakal, B. Fidan, and A. Savran, “Mixing adaptive fault tolerant control of quadrotor UAV,” _Asian Journal of Control_ , vol. 19, no. 4, pp. 1441–1454, Jul. 2017.
  * [11] J. L. Mancilla-Aguilar and R. A. García, “An algorithm for the robust exponential stabilization of a class of switched systems,” _International Journal of Robust and Nonlinear Control_ , vol. 25, no. 13, pp. 2062–2082, Sep. 2015.
  * [12] G. Evensen, _Data Assimilation: The Ensemble Kalman Filter_ , 2nd ed. Berlin, Heidelberg: Springer-Verlag, 2009.
  * [13] A. Dehghani, B. D. O. Anderson, and A. Lanzon, “Unfalsified adaptive control: A new controller implementation and some remarks,” in _IEEE European Control Conference_ , 2007, pp. 709–716.
  * [14] S. Baldi, G. Battistelli, E. Mosca, and P. Tesi, “Multi-model unfalsified adaptive switching supervisory control,” _Automatica_ , vol. 46, no. 2, pp. 249–259, 2010.
  * [15] K. S. Narendra and K. Esfandiari, “Adaptive identification and control of linear periodic systems using second-level adaptation,” _International Journal of Adaptive Control and Signal Processing_ , vol. 33, no. 6, pp. 956–971, Jun. 2019.
  * [16] Y. Zhang, X. Wang, and Z. Wang, “Adaptive multiple model control for a class of nonlinear discrete time systems: second-level adaption design approach,” _International Journal of Control_ , vol. 96, no. 2, pp. 497–507, 2023.
  * [17] V. K. Pandey, I. Kar, and C. Mahanta, “Controller design for a class of nonlinear mimo coupled system using multiple models and second level adaptation,” _ISA Transactions_ , vol. 69, pp. 256–272, 2017.
  * [18] L. Dutta and D. Kumar D., “Adaptive model predictive control design using multiple model second level adaptation for parameter estimation of two‐degree freedom of helicopter model,” _International Journal of Robust and Nonlinear Control_ , vol. 31, no. 8, pp. 3248–3278, 2021.
  * [19] K. Moffat and C. Tomlin, “The Multiple Model Adaptive Power System State Estimator,” in _IEEE Conference on Decision and Control_ , Dec. 2021, pp. 3525–3530.
  * [20] M. Imani and U. Braga-Neto, “Multiple Model Adaptive controller for Partially-Observed Boolean Dynamical Systems,” in _American Control Conference_ , May 2017, pp. 1103–1108.
  * [21] S. Wang, W. Ren, and J. Chen, “Fully distributed state estimation with multiple model approach,” in _IEEE Conference on Decision and Control_ , Dec. 2016, pp. 2920–2925.
  * [22] P. Ioannou and B. Fidan, _Adaptive Control Tutorial_. SIAM, 2006.
  * [23] G. Tao, _Adaptive Control Design and Analysis_. Wiley, 2003.
  * [24] Z. Han and K. S. Narendra, “New concepts in adaptive control using multiple models,” _IEEE Transactions on Automatic Control_ , vol. 57, no. 1, pp. 78–89, Jan. 2012.
  * [25] N. Ahmadian, A. Khosravi, and P. Sarhadi, “A New Approach to Adaptive Control of Multi-Input Multi-Output Systems Using Multiple Models,” _ASME Journal of Dynamic Systems, Measurement, and Control_ , vol. 137, no. 9, p. 091009, Sep. 2015.
  * [26] A. Lovi, B. Fidan, and C. Nielsen, “Multiple model reference adaptive tracking control of multivariable systems with blending,” in _IEEE Conference on Decision and Control_ , 2022, pp. 1362–1367.
  * [27] B. Kolodziejczak and T. Szulc, “Convex combinations of matrices — full rank characterization,” _Linear Algebra and its Applications_ , vol. 287, no. 1-3, pp. 215–222, 1999.
  * [28] S. P. Boyd and L. Vandenberghe, _Convex Optimization_ , 2nd ed. Cambridge, UK: Cambridge University Press, 2006.
  * [29] K. S. Narendra, W. Yu, and C. Wei, “Stability, robustness, and performance issues in second level adaptation,” in _American Control Conference_ , 2014, pp. 2377–2382.


#  Towards Multiphase Clocking in Single-Flux Quantum Systems

Rassul Bairamkulov and Giovanni De Micheli Lausanne, Switzerland
<EMAIL_ADDRESS><EMAIL_ADDRESS>Integrated Systems
Laboratory, EPFL

###### Abstract

Rapid single-flux quantum (RSFQ) is one of the most advanced superconductive
electronics technologies. SFQ systems operate at tens of gigahertz with up to
three orders of magnitude smaller power as compared to CMOS. In conventional
SFQ systems, most gates require clock signal. Each gate should have the fanins
with equal logic depth, necessitating insertion of path-balancing (PB) DFFs,
incurring prohibitive area penalty.

Multiphase clocking is the effective method for reducing the path-balancing
overhead at the cost of reduced throughput. However, existing tools are not
directly applicable for technology mapping of multiphase systems. To overcome
this limitation, in this work, we propose a technology mapping tool for
multiphase systems. Our contribution is threefold. First, we formulate a phase
assignment as a Constraint Programming with Satisfiability (CP-SAT) problem,
to determine the phase of each element within the network. Second, we
formulate the path balancing problem as a CP-SAT to optimize the number of
DFFs within an asynchronous datapath. Finally, we integrate these methods into
a technology mapping flow to convert a logic network into a multiphase SFQ
circuit. In our case studies, by using seven phases, the size of the circuit
(expressed as the number of Josephson junctions) is reduced, on average, by
59.94 % as compared to the dual (fast-slow) clocking method, while
outperforming the state-of-the-art single-phase SFQ mapping tools.

††This work was supported by the SNF grant “Supercool: Design methods and
tools for superconducting electronics” under Grant 200021-1920981. The authors
thank Professor Peter Beerel from the University of Southern California for
the insightful discussions.

a) b) c)

Figure 1:  a) An example of a CMOS circuit. b) Equivalent SFQ circuit with a
splitter and two path balancing DFFs. c) A four-phase SFQ realization with no
path balancing DFFs. The numbers indicate the phase corresponding to each
clocked gate.

##  I Introduction

Superconductive electronics is one of the most promising beyond-CMOS
technologies, with Rapid Single-Flux Quantum (RSFQ) [1] as the most prominent
superconductive technology. RSFQ systems consistently achieve operating
frequencies on the order of tens of gigahertz [2, 3], with particular
structures accelerating to hundreds of gigahertz [4, 5]. Furthermore, the
operating power of the RSFQ systems is two to three orders of magnitude
smaller than CMOS, even considering the refrigeration power [6]. Due to the
low power, low noise, and high speed RSFQ systems are utilized in high-
resolution sensors [7] and high-performance wireless communication systems
[8]. In addition, these advantages make the RSFQ a promising candidate for
large-scale stationary computing systems (e.g., data centers) [6], as well as
low-power computing in space [9].

Elevating complexity of RSFQ systems is however a challenging task. Due to the
major differences between RSFQ and CMOS technologies, standard design
methodologies and tools are not applicable to RSFQ. In CMOS, transistors
encode information by voltage levels. In contrast, in RSFQ, Josephson
Junctions (JJ) encode information using single flux quantum (SFQ) pulses [10].
111 For an in-depth discussion of the SFQ technology we refer an interested
reader to [9, 10]. The absence or presence of a SFQ pulse encode,
respectively, a logic $\mathtt{0}$ or $\mathtt{1}$. The flow of pulses within
a logic network requires synchronization. While a few structures, such as
splitter and $\mathtt{OR}$, can operate asynchronously, most logic gates in
the original RSFQ cell library, such as $\mathtt{AND}$, $\mathtt{NOT}$ and
$\mathtt{XOR}$, require clock signal [1]. This feature requires the data to be
pipelined at the gate level.

Two features further complicate scaling the RSFQ systems. Due to the quantized
nature of SFQ pulses, most RSFQ logic gates are limited to fanout of one. A
special gate called splitter is necessary to duplicate a signal [9, 10], as
illustrated in Fig. 1b. In addition, path balancing is required to ensure the
correct order of data propagation within the network, as indicated by four
path balancing D-flip-flops (DFF) in Fig. 1b. The number of path balancing
DFFs and splitters can be prohibitively large, degrading the area and yield.
The complexity of the clock tree synthesis is also increased due to the large
number of clocked elements. Despite the advances in technology mapping [11,
12], the overhead of path balancing and clock distribution remains large,
motivating the investigation of alternative RSFQ architectures.

Alternative clocking methods have been explored in the literature to alleviate
the issue of path balancing. In dual clock methodology [13], a logic circuit
is partitioned into separate clocking domains. The gates within a single
domain are clocked at high frequency $f_{hi}$, while cross-domain
communication is performed via the nondestructive readout (NDRO) registers
operating at the low frequency $f_{lo}=f_{hi}/N$. In [13], the number of path
balancing DFFs is reduced by 35% with $N=5$ [13]. The throughput of the system
is however reduced by a factor of 5, establishing a tradeoff between the area
and the throughtput of the system. The primary limitation of the dual clock
method is the use of the large registers to repeat the data signals and
additional $\mathtt{AND}$ gates to suppress the spurious pulses within the
system [14]. In addition, similar to full path balancing, distributing the
high frequency clock signal makes the timing closure more challenging to
achieve.

Multiphase clocking has recently been proposed for designing area-efficient
SFQ systems [15]. Using several phase-shifted clock signals allows data to
propagate to the subsequent stages without DFFs. For example, by using four-
phase clocking, the path balancing DFFs can be completely removed, as shown in
Fig. 1c. Using multiple phases was demonstrated effective in reducing the SFQ
circuit area, thanks to the reduction in DFF count (i.e., smaller logic
network) and fewer clocked elements yielding smaller clock network size. With
only two phases, the number of DFFs and area are reduced by, respectively, 55%
and 41% [15], despite the overhead of generating and distributing separate
clock signals. Additional phases further reduce the system area. Furthermore,
in an $n$-phase system, the clock signals have $n$ times lower frequency,
simplifying the clock distribution process. Similar to dual clock method, the
area reduction is achieved by sacrificing the throughput. An $n$-phase system
has $n$ times smaller throughput as compared to a single phase system.

Although the multiphase clocking has been shown to reduce the path balancing
overhead, existing technology mapping tools offer limited support for the
multiphase systems. For example, path balancing DFFs cannot be efficiently
distributed within the complex datapaths utilizing asynchronous SFQ gates. To
bridge this gap, in this paper, we propose a novel methodology to apply the
multiphase clocking to SFQ circuits. Our contribution is threefold,

  1. 1.

we extend the DFF insertion methodology to assign the phase to each gate
within the circuit, including unclocked elements.

  2. 2.

we formulate PB DFF insertion as a Constraint Programming with SAT (CP-SAT)
problem to satisfy the timing constraints with minimum number of PB DFFs. We
improve the scalability by identifying independent paths allowing the
asynchronous paths to be processed separately.

  3. 3.

we integrate phase assignment and DFF insertion into the technology mapping
flow to realize an arbitrary logic network with multiple clock phases.

With only two phases, our flow outperforms the state-of-the-art single-phase
methods. Compared with the dual clocking method, our technology mapping flow
achieves up to 82% smaller size, while offering equal throughput.

The rest of the paper is organized as follows. The principles of SFQ
technology and multiphase clocking are reviewed in Section II. The CP-SAT-
based method for phase assignment is presented in Section III. PB DFF
insertion are described in Section IV. Experimental results and comparison
with the prior works are provided in Section V, followed by conclusions in
Section VI.

##  II Background

A logic network $\mathcal{N}=(\mathcal{V},\mathcal{E})$ is a directed acyclic
graph (DAG) where $\mathcal{V}$ is a set of nodes and
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is a set of edges. The set
of nodes $\mathcal{V}=\mathcal{I}\cup\mathcal{O}\cup\mathcal{G}$ consists of
three disjoint subsets representing, respectively, the primary inputs (PI),
primary outputs (PO) and gates within the network. A set of fanins $FI(g)$
(fanouts $FO(g)$) of gate $g\in\mathcal{V}$ denotes the nodes connected to $g$
via and incoming (outgoing) edge.

Technology mapping refers to the process of transforming an arbitrary logic
network $\mathcal{N}$ into a technology-specific representation
$\mathcal{N}^{\prime}$ using a particular set of primitives. The technology
mapping aims to optimize the target characteristics of the mapped network
$\mathcal{N}^{\prime}$, such as area, delay, or yield The nature of a
technology greatly influences the technology mapping process. In the next
subsection, the properties of the SFQ technology are reviewed.

###  II-A Single-Flux Quantum technology

Rapid Single-flux Quantum (RSFQ) is a cryogenic superconductive computing
logic family based on Josephson junctions (JJ). RSFQ gates consist of
superconductive loops storing quantized magnetic flux. The information is
transferred between the superconductive loops in form of SFQ pulses with the
area of $\Phi_{0}=\hbar/2e\approx 2.07$ mV$\cdot$ps [9]. Based on
synchronization mechanisms, the SFQ logic gates can be divided into three
major categories [16].

  * •

Asynchronous input, Asynchronous output ($\mathtt{AA}$) components process the
input information immediately upon arrival. The most common gates in this
category are splitter and merger. A splitter produces two pulses for each
incoming pulse. A merger directs signals from two input branches into one
output branch, effectively performing an $\mathtt{OR}$ function.

  * •

Asynchronous input, Synchronous output ($\mathtt{AS}$) elements process the
input information immediately upon arrival and release the output
synchronously after the arrival of the clock signal. The components of this
type include D-flip-flop (DFF), inverter ($\mathtt{NOT}$) and exclusive-or
($\mathtt{XOR}$).

  * •

Synchronous input, Asynchronous output ($\mathtt{SA}$) elements require the
input signals to arrive simultaneously to operate correctly. The result of the
computation is released immediately after processing. Assuming inputs arrive
simultaneously, a CB can be tuned to produce at most a single output pulse,
producing an $\mathtt{OR}$ element [17]. Furthermore, by adjusting the JJ size
and bias current, the $\mathtt{OR}$ structure can be transformed into
$\mathtt{AND}$ element, requiring both signals to arrive simultaneously to
produce an output pulse.

In conventional RSFQ [1], the $\mathtt{OR}$ and $\mathtt{AND}$ gates
incorporate the DFFs at the inputs to ensure the simultaneous release of the
data pulses. However, the simultaneity can be achieved by other $\mathtt{AS}$
gates, such as $\mathtt{NOT}$ and $\mathtt{XOR}$ [16], allowing richer
functionality to be realized in fewer clock cycles.

To correctly consider the type of each element, we divide the set of logic
gates into three disjoint subsets
$\mathcal{G}=\mathcal{G}_{\mathtt{AA}}\cup\mathcal{G}_{\mathtt{AS}}\cup\mathcal{G}_{\mathtt{SA}}$,
where each subset represents the elements of the corresponding category. Note
that only gates in the subset $\mathcal{G}_{\mathtt{AS}}$ require clock
signal.

###  II-B Multiphase clocking

Multiphase clocking (MPC) is a clocking technique based on several clock
signals with equal clock period. The seminal work by Li et al. [15]
demonstrated MPC as an effective remedy for path balancing in SFQ systems. A
$n$-phase system utilizes $n$ periodic signals $\\{t_{0},\cdots,t_{n-1}\\}$
operating at the same frequency and different phases. Each clocked element
$g\in\mathcal{G}_{\mathtt{AS}}$ within the network is synchronized by only one
clock signal at phase $\varphi(g)$. The epoch $S(g)$ of a gate $g$ is defined
as the number of clock cycles separating the gate $g$ from the PIs, as
illustrated in Fig. 2. The clock signals are ordered by phase
$\varphi\in\\{0,\cdots,n-1\\}$, i.e., during any epoch, the clock signal
$t_{i}$ arrives before clock signal $t_{j}$ if $i<j$. For consistency in I/O
timing, the PIs are placed at the same epoch $0$ and the POs are placed at the
same epoch. For example, the POs are Fig. 2, the fanouts $\mathtt{x}$ and
$\mathtt{y}$ are both placed at epoch 2. Observe that a path balancing DFF is
inserted before $\mathtt{x}$ to equalize the epochs of both POs. For
convenience, we define a _stage_ $\sigma(g)$ of a gate $g$ as

$\sigma(g)=nS(g)+\varphi(g).$ (1)

Figure 2:  An example two-phase network. Each PI is placed at stage 0 or 1, at
epoch 0. The POs are placed at stages 4 and 5, corresponding to epoch 2.

Multiphase clocking relaxes the path balancing requirements of a clocked
system. With $n$ phases, the stage difference $\Delta\sigma(i,j)$ between any
pair of sequentially adjacent elements $(i,j)$ should be no greater than $n$.
Compared with single-phase systems where a $q$-stage difference is balanced by
$q$ DFFs, the same difference within an $n$-phase system can be balanced with
only $\left\lfloor\frac{q}{n}\right\rfloor$ DFFs, greatly reducing the path
balancing overhead in imbalanced paths.

To realize the potential of multiphase clocking in SFQ systems, we propose the
two-step DFF insertion methodology based on CP-SAT. During the phase
assignment step, the SFQ gates are assigned a phase to reduce the number of
DFFs within the system, as described in Section III. During the DFF insertion
step, optimal DFF insertion is determined for each datapath, as described in
Section IV.

##  III Phase assignment

The phase assignment is the procedure of determining the stage (i.e., epoch
and phase) of each gate within the network. An ILP-based methodology for phase
assignment to reduce the number of path balancing DFFs is proposed in [15],

$\min\limits_{\sigma(g)\ \forall
g\in\mathcal{G}}\sum\limits_{(i,j)\in\mathcal{E}}\left\lfloor\frac{\sigma(j)-\sigma(i)}{n}\right\rfloor,$
(2) Subject to:  $\sigma(i)<\sigma(j)\hskip
10.00002pt\forall(i,j)\in\mathcal{E},$ (3)
$\left\lfloor\frac{\sigma(i)}{n}\right\rfloor=0\hskip 10.00002pt\forall
i\in\mathcal{I},$ (4)
$\left\lfloor\frac{\sigma(i)}{n}\right\rfloor=\left\lfloor\frac{\sigma(j)}{n}\right\rfloor\hskip
10.00002pt\forall i,j\in\mathcal{O}.$ (5)

This formulation is highly effective when considering the systems consisting
of only clocked gates, i.e. $\mathcal{G}=\mathcal{G}_{\mathtt{AS}}$,
$\mathcal{G}_{\mathtt{AA}}=\mathcal{G}_{\mathtt{SA}}=\varnothing$ This
formulation however offers limited support for the unclocked gates, such as
mergers and $\mathtt{SA}$ $\mathtt{AND}$ gate. Unclocked elements can greatly
improve the quality of SFQ technology mapping by enriching the functionality
available within a single clock cycle. In [18], for example, the area of the
circuit was reduced by up to 54% by utilizing asynchronous elements. However,
the unclocked gates impose timing constraints different from the clocked
gates. To extend the support of this ILP model to unclocked gates, additional
constraints are introduced in this section.

While the notions of phase, epoch, and stage are primarily relevant for the
clocked gates, it is however important to also assign a stage to the
$\mathtt{AA}$ and $\mathtt{SA}$ elements, since a stage of an unclocked
element guides the DFF insertion. Consider two networks containing an
asynchronous datapath illustrated in Fig. 3. Both networks realize the same
function, however, the possible locations of path balancing DFFs are
different. Assigning a stage to an unclocked gate is important to correctly
place the path balancing DFFs.

Figure 3:  The effect of phase assigned to splitters and mergers. Datapaths
(a) and (b) execute the same function. The phases assigned to the splitters
and mergers (denoted as $\mathtt{S}$ and $\mathtt{M}$, respectively) are
however different. The potential DFF locations (depicted as yellow rectangles)
are therefore different in the two datapaths.

###  III-A $\mathtt{SA}$ element constraints

Recall that the inputs of the $\mathtt{SA}$ gates require simultaneous arrival
of data pulses to the gate. To ensure the simultaneous arrival, the
$\mathtt{AS}$ gate should be placed immediately before the $\mathtt{SA}$ gate,
at the same clock stage, as shown in Fig. 4a. If the stage of $\mathtt{AS}$
gate is earlier than the $\mathtt{SA}$ gate, a path balancing DFF should be
inserted. In addition, the $\mathtt{SA}$ gate cannot be placed at the same
stage as the fanin, if the fanin is not an $\mathtt{AS}$ gate. For example, if
the fanin of the $\mathtt{SA}$ gate has fanout greater than one, a splitter is
inserted after the $\mathtt{AS}$ gate. Thus, the fanin type of the
$\mathtt{SA}$ gate is a splitter ($\mathtt{AA}$ gate). The propagation delay
incurred by the $\mathtt{AA}$ gates, such as splitters, may desynchronize the
inputs arriving to the $\mathtt{SA}$ gate [18], producing a data hazard, as
illustrated in Fig. 4b. To avoid this condition, the minimum stage of gate $g$
is increased by one, as shown in Fig. 4c.

a) b) c)

Figure 4:  Timing constraints of an $\mathtt{SA}$ gate. a) The gate should be
preceded by an $\mathtt{AS}$ gate at the same stage. b) Example of an invalid
assignment. The $\mathtt{AS}$ gate is preceded by an $\mathtt{AA}$ gate at the
same stage. c) The issue is resolved by shifting the $\mathtt{SA}$ gate to the
next stage and inserting two path balancing DFFs.

###  III-B Asynchronous paths

The sharing of DFFs by a complex unclocked path has been highlighed as early
as 1991 in [19], where sharing of DFFs is maximized to minimize the circuit
area. While estimating the minimum number of DFFs in special cases, such as
splitter tree, is relatively simple, more complex asynchronous paths, however,
require more sophisticated processing and, hence, prohibitive runtime. During
the phase assignment stage we therefore approximate the number of DFFs as

$C\left(\sigma\left(g\right)\right)=\left\lfloor\frac{\max\limits_{a\in
FO\left(g\right)}\left(\sigma\left(a\right)\right)-\sigma\left(g\right)}{n}\right\rfloor.$
(6)

The exact number of DFFs is determined during the DFF insertion stage
described in Section IV. The combined optimization problem is formulated as

$\min\limits_{\sigma\left(g\right)\hskip 7.00002pt\forall
g\in\mathcal{G}}\sum\limits_{\left(i,j\right)\in\mathcal{E}}\left\lfloor\frac{\sigma\left(j\right)-\sigma\left(i\right)+\left(j\in\mathcal{G}_{\mathtt{SA}}\right)}{n}\right\rfloor,$
(7) Subject to:  $\sigma\left(i\right)<\sigma\left(j\right)\hskip
10.00002pt\forall\left(i,j\right)\in\mathcal{E},j\in\mathcal{G}_{\mathtt{AS}},$
(8) $\sigma\left(i\right)\leq\sigma\left(j\right)\hskip
10.00002pt\forall\left(i,j\right)\in\mathcal{E},j\notin\mathcal{G}_{\mathtt{AS}},$
(9)
$\sigma\left(g\right)\geq\sigma\left(a\right)+\left(a\notin\mathcal{G}_{\mathtt{AS}}\right)\hskip
10.00002pt\forall g\in\mathcal{G}_{\mathtt{SA}},a\in FI\left(g\right),$ (10)
$\left\lfloor\frac{\sigma\left(i\right)}{n}\right\rfloor=0\hskip
10.00002pt\forall i\in\mathcal{I},$ (11)
$\left\lfloor\frac{\sigma\left(i\right)}{n}\right\rfloor=\left\lfloor\frac{\sigma\left(j\right)}{n}\right\rfloor\hskip
10.00002pt\forall i,j\in\mathcal{O}.$ (12)

Each gate $g$ with $\\#FO(g)>1$ requires a splitter tree to copy the signal to
each of the fanouts. After completing the phase assignment, a splitter tree
that maximizes sharing of path balancing DFFs can be produced. First, the
fanouts $FO(g)$ of a gate $g$ are sorted by phase in the ascending order,
producing a sequence $A=[a_{0},a_{1},\cdots,a_{k}]$,
$\sigma(a_{i})\leq\sigma(a_{j})\iff i<j$. Next, the splitters
$[s_{0},...,s_{k-1}]$ are inserted at phases $a_{0},\cdots,a_{k-1}$. Each
splitter $s_{i}$ has a splitter $s_{i-1}$ as a fanin, except $s_{0}$ whose
fanin is the gate $g$. Each splitter $s_{i}$ has the $a_{i}$ and $s_{i+1}$ as
a fanout, except $s_{k-1}$ whose fanouts are $a_{k-1}$ and $a_{k}$. Using this
method, the splitters are placed as late as possible within the network,
maximizing sharing of the datapath.

##  IV DFF placement

The complex asynchronous datapaths greatly complicate the DFF insertion
process. Consider for example the datapath shown in Fig. 3. The optimal
placement of the DFFs cannot be easily inferred from the datapath topology. In
this section, we propose a CP-SAT based placement methodology for determining
the optimal location of the DFFs. Our methodology consists of two stages.
First, we identify independent datapaths within an unbalanced logic network,
as described in Subsection IV-A. Next, for each independent datapaths, we
formulate a CP-SAT problem where the timing constraints are satisfied using
the minimum number of DFFs, as described in Subsection IV-B.

###  IV-A Independent datapaths

Recall that only a single clocked gate can be placed at the same stage along a
datapath. The $\mathtt{AS}$ and $\mathtt{SA}$ elements therefore prevent a DFF
from being placed at the same stage. Furthermore, any timing constraint
imposed by a DFF placed before the $\mathtt{AS}$/$\mathtt{SA}$ element, will
be dominated by the timing constraint imposed by the
$\mathtt{AS}$/$\mathtt{SA}$ element itself at any location after the
$\mathtt{AS}$/$\mathtt{SA}$ element. This feature allows us to represent a
network-wide DFF insertion problem as multiple smaller problems. This
partitioning capability brings two major advantages to the DFF insertion
process. First, while the network-wide DFF insertion problem can be solved in
$O(e^{N})$ time, the complexity of a partitioned problem is
$O(e^{k}N_{path})$, where $N$ is the network size, $N_{path}\propto N$ is the
number of independent datapaths and $k$ is the size of the largest partition
(typically, $k\ll N$). Second, partitioning facilitates multithreaded
processing, allowing faster processing on a multicore system.

The proposed DFF insertion method is applied to a single independent datapath
constrained by the synchronous gates at the input and the output, as
illustrated in Fig. 5a. The datapath $P$ can be described as a portion of the
logic network $P=(I,A,O)$, where set $A\subseteq\mathcal{G}_{\mathtt{AA}}$
contains the $\mathtt{AA}$ gates within the datapath and sets
$I,O\subseteq\mathcal{G}_{\mathtt{AS}}\cup\mathcal{G}_{\mathtt{SA}}$ describe,
respectively, the $\mathtt{AS}$ and $\mathtt{SA}$ gates at the input and
output of the datapath. In Fig. 5a, the three independent paths are

$P_{1}=\left(I=\\{\mathtt{A},\mathtt{B},\mathtt{D}\\},A=\\{\mathtt{1},\mathtt{2},\mathtt{3},\mathtt{4},\mathtt{5}\\},O=\\{\mathtt{W},\mathtt{X},\mathtt{Y},\mathtt{Z}\\}\right),$
$P_{2}=\left(I=\\{\mathtt{E}\\},A=\varnothing,O=\\{\mathtt{W}\\}\right),$
$P_{3}=\left(I=\\{\mathtt{C},\mathtt{F}\\},A=\\{\mathtt{6}\\},O=\\{\mathtt{Z}\\}\right).$

###  IV-B DFF insertion

After identifying each independent path within the network, we determine the
potential DFF sites for subsequent DFF insertion. To uniquely identify each
potential DFF site, we define $d=(g_{d}^{i},g_{d}^{o},\sigma(d))$, where
$g_{d}^{i}\in I\cup A$ and $g_{d}^{o}\in A\cup O$ are the fanin and fanout
elements of the DFFs and $\sigma(d)$ is the stage of the DFF site. We define a
chain $Q=(d_{min},\cdots,d_{max})$ as a sequence of sequentially adjacent DFF
locations situated between $d_{min}$ and $d_{max}$. The examples of DFF sites
and a chain are shown in Fig 5b. For convenience, we define the length of
chain $\Delta\sigma(Q)$ as the stage difference between $d_{min}$ and
$d_{max}$.

For each DFF site, we introduce a binary variable $\delta(d)$ equal to
$\mathtt{1}$ if the DFF is placed at $d$ and $\mathtt{0}$ otherwise. The
problem of minimizing the number of path balancing DFFs can therefore be
formulated as a CP-SAT problem minimizing

$\sum\limits_{d\in P}\delta(d).$ (13)

Several constraints describe the valid placement of the DFFs within an
independent path. First, for each stage, only a single DFF can be placed along
a chain. In Fig. 5b, for example, the stage $\sigma=5$ has five potential DFF
locations, $a,b,c,d$, and $e$. The DFFs here should however satisfy the
constraint

$\sum\limits_{d\in Q,\sigma_{d}=i}\delta(d)\leq 1.$ (14)

The combinations of DFF sites satisfying this constraint are $\\{a,c,e\\}$,
$\\{a,d\\}$, $\\{b,e\\}$, since only these sets do not place multiple DFFs
along any chain. Other conflicting DFF sites can be found at stages 3 and 6.

Second, the distance between any pair of clocked elements along a chain should
not exceed $n$ in an $n$-phase system. This requirement can be formulated as a
condition

$\bigvee\limits_{d\in Q,\Delta\sigma(Q)=n}\delta(d)=\mathtt{1}$ (15)

stating that at least one DFF should be placed along any chain of length $n$.

Finally, recall that the $\mathtt{SA}$ gates should be preceded by an
$\mathtt{AS}$ gate to function correctly. Thus, if an $\mathtt{SA}$ element is
not directly preceded by an $\mathtt{AS}$ gate, a DFF should be placed before
the $\mathtt{SA}$ element,

$\delta(d)=\mathtt{1}\hskip 10.00002pt\forall
g_{d}^{o}\in\mathcal{G}_{\mathtt{SA}},\sigma(g_{d})=\sigma_{d}.$ (16)

The DFF sites before the $\mathtt{SA}$ gates are shaded with diagonal hatch
pattern in Fig. 5b.

Equations (13)-(16) constitute a CP-SAT problem where the conditions (14)-(16)
are satisfied with the minimum number of DFFs (13).

Figure 5: Example of an asynchronous datapath within a four-phase network
($n=4$). The black triangles, and black curved shapes represent the
$\mathtt{AS}$ and $\mathtt{SA}$ elements, respectively. Red circles denote the
$\mathtt{AA}$ elements. a) Three independent paths are drawn with red solid,
black solid, and black dotted arrows. b) DFF sites shown as rectangles within
the red independent path. The grey DFF sites represent the chain $\mathtt{Q}$
of length $n$. Diagonally shaded DFF sites precede the $\mathtt{SA}$ gates.
DFFs should therefore be placed at these sites.

##  V Experimental results

We integrate the DFF placement methodology into the technology mapping flow
for SFQ compound gate circuits. The SFQ circuits are synthesized with
$\mathtt{mockturtle}$ using a database of precomputed compound gate structures
[18]. The circuits are next decomposed into the functional primitives while
removing all DFFs within the circuit. We next apply our CP-SAT-based phase
assignment and DFF insertion procedures using Google OR-Tools [20].

We apply our mapping flow to synthesize a subset of EPFL [21] and ISCAS [22]
benchmark circuits. For each benchmark circuit, the total runtime does not
exceed 60.64 seconds. Phase assignment is a major component of the total
runtime. Optimal phase assignment has been successfully determined in all
circuits. In most cases, the phase assignment is completed in under 5 seconds,
with the longest runtime (42 seconds) required by the $\mathtt{voter}$
circuit. Naturally, larger circuits, such as $\mathtt{s13207}$,
$\mathtt{voter}$ and $\mathtt{priority}$ require longer runtime for phase
assignment. Applying our tool to larger circuits would require additional
processing time. Note that although finding the optimal phase assignment may
require prohibitive runtime, a feasible solution is often available much
earlier. Therefore, by relaxing the optimality requirement, the our framework
can support larger circuits. 222Mapping of $\mathtt{hyp}$ (the largest circuit
in EPFL suite) has been completed within one hour with only 15 and 20 minutes
allocated to phase assignment and DFF insertion, respectively (not presented
due to space constraints).

TABLE I: Comparison of multiphase clocking with dual clocking method [13] for different throughputs | 1/7 throughput | 1/12 throughput
---|---|---
| DCM | Multiphase | Change | Runtime | DCM | Multiphase | Change | Runtime
| #DFF | #JJ | #DFF | #JJ | #DFF | #JJ | (s) | #DFF | #JJ | #DFF | #JJ | #DFF | #JJ | (s)
$\mathtt{int2float}$ | 117 | 7’770 | 217 | 5’136 | +85% | –34% | 4.42 | 39 | 5’140 | 46 | 3’939 | +18% | –23% | 4.13
$\mathtt{priority}$ | 8’562 | 257’252 | 3’285 | 45’094 | –62% | –82% | 54.01 | 4’225 | 158’568 | 1’775 | 34’524 | –58% | –78% | 30.33
$\mathtt{voter}$ | 7’204 | 447’044 | 2’180 | 162’804 | –70% | –64% | 60.64 | 3’732 | 355’144 | 1’568 | 158’520 | –58% | –55% | 48.94
$\mathtt{c432}$ | 224 | 10’734 | 342 | 5’116 | +53% | –52% | 14.92 | 118 | 7’124 | 240 | 4’402 | +103% | –38% | 10.14
$\mathtt{c880}$ | 362 | 14’658 | 254 | 6’190 | –30% | –58% | 7.56 | 187 | 9’483 | 119 | 5’245 | –36% | –45% | 7.25
$\mathtt{c1908}$ | 282 | 13’169 | 125 | 3’529 | –56% | –73% | 10.22 | 144 | 8’739 | 69 | 3’137 | –52% | –64% | 8.78
$\mathtt{c3540}$ | 776 | 43’437 | 589 | 17’016 | –24% | –61% | 11.83 | 282 | 26’897 | 440 | 15’973 | +56% | –41% | 9.29
$\mathtt{c1355}$ | 193 | 8’739 | 46 | 4’515 | –76% | –48% | 4.18 | 119 | 6’149 | 44 | 4’501 | –63% | –27% | 3.89
$\mathtt{s13207}$ | 1’795 | 106’346 | 1’837 | 42’382 | +2% | –60% | 27.87 | 571 | 60’766 | 1’082 | 37’097 | +89% | –39% | 25.90
$\mathtt{s5378}$ | 645 | 50’766 | 808 | 21’761 | +25% | –57% | 7.87 | 255 | 34’053 | 368 | 18’681 | +44% | –45% | 7.07
$\mathtt{s382}$ | 56 | 4’448 | 89 | 2’411 | +59% | –46% | 4.17 | 9 | 2’750 | 48 | 2’124 | +433% | –23% | 3.93
Geomean | 557.04 | 29’868.20 | 413.48 | 11’965.54 | –25.8% | –59.9% | 12.02 | 227.85 | 19’565.38 | 229.57 | 10’467.10 | +0.8% | –47.5% | 10.03

Figure 6:  Normalized average number of DFFs (left) and JJs (right) with
different number of phases. The lines in the background represent normalized
number of DFFs and JJs for individual benchmarks.

We first compare our results with dual clocking method [14]. Two experiments
are reported where the frequency of the slow clock is 7 and 12 times smaller
than the frequency of the fast clock. For fair comparison, we use 7 and 12
phases to equalize the throughput. The results are shown in Table I. Note that
unlike [14], the DFFs within the $\mathtt{AND}$ and $\mathtt{OR}$ gate are
explicitly counted in our work due to the use of the unclocked ($\mathtt{SA}$)
$\mathtt{AND}$ and $\mathtt{OR}$ elements. Furthermore, additional elements,
such as pulse repeaters and NDRO flip-flops are needed in DCM. We therefore
focus on the number of JJs, commonly used in SFQ technology as an area metric.
The number of JJs is reduced, on average, by 59.9% when using 7-phase
clocking. The largest reduction is observed in the $\mathtt{priority}$
benchmark circuit containing many imbalanced paths that require long chains of
path balancing DFFs to be inserted. With multiphase clocking, these chains can
be balanced by approximately 7 times fewer DFFs, contributing to the area
savings.

For 12-phase clocking, the area is reduced by 47.5%. Observe that the smallest
reduction in area (and the increase in DFF count) are observed in circuits
with relatively small depth. These circuits therefore have fewer imbalanced
paths that could benefit from multiphase clocking. Despite relatively large
number of DFFs, we still manage to reduce the total JJ count due to the use of
unclocked elements ($\mathtt{AA}$ and $\mathtt{SA}$) supported by our tool.

The Fig. 6 (left) describes the relationship between the number of DFFs and
the number of phases. The numbers of DFFs are normalized with respect to a
single phase system mapped with our flow. Consistent with [15] the number of
DFFs reduces with multiple phases. The area savings due to an additional phase
gradually diminishes as more phases are added. We compare our results against
PBMap [11] and gate compounding [18], the state-of-the-art single-phase SFQ
mapping algorithms. With our mapping flow, two phases are sufficient, on
average, to outperform gate compounding. With three phases our flow yields
fewer DFFs than PBMap. The Fig. 6 (right) describes the relationship between
the circuit area (expressed as the number of JJs) and the number of phases.
Our flow outperforms PBMap with a single phase, primarily due to the use of
asynchronous gates (e.g., merger instead of an $\mathtt{OR}$ gate). A two-
phase network synthesized with our flow requires fewer JJs than a single-phase
compound-gate network.

##  VI Conclusions

RSFQ technology presents a remarkable opportunity to achieve unprecedented
performance and energy efficiency of computing systems. Realizing its full
potential however requires overcoming major technological issues, including
path balancing. Multiphase clocking can substantially reduce the path
balancing overhead in SFQ systems. Existing path balancing tools however offer
limited support of multiphase clocking in systems with asynchronous SFQ gates.
In this work, we presented a technology mapping and path balancing methodology
for multiphase SFQ systems based on constraint programming with satisfiability
(CP-SAT). An SFQ logic network is initially synthesized with
$\mathtt{mockturtle}$ logic synthesis library. Next, the circuit is decomposed
into primitives and the DFFs are removed from the network. Each primitive gate
is assigned a phase using the CP-SAT while satisfying the special timing
constraints imposed by the unclocked SFQ elements. Finally, using the novel
CP-SAT formulation, the minimum number of path balancing DFFs satisfying the
timing constraints is determined. In the experimental results, we showed an
average of 59% reduction in the number of JJs when compared to dual clocking
method. Furthermore, with only two phases, our methodology yields smaller
networks than the state-of-the-art single-phase SFQ mapping techniques.

## References

  * [1] K. Likharev, O. Mukhanov, and V. Semenov,  “Resistive Single Flux Quantum Logic for the Josephson-Junction Digital Technology,”  Proc. SQUID, Vol. 85, June 1985.
  * [2] T. Kawaguchi, M. Tanaka, K. Takagi, and N. Takagi,  “Demonstration of an 8-Bit SFQ Carry Look-Ahead Adder Using Clockless Logic Cells,”  Proc. ISEC, July 2015.
  * [3] Z. J. Deng, N. Yoshikawa, S. R. Whiteley, and T. Van Duzer,  “Data-Driven Self-Timed RSFQ High-Speed Test System,”  IEEE TASC, Vol. 7, No. 4, 1997.
  * [4] W. Chen _et al._ ,  “Rapid Single Flux Quantum T-Flip Flop Operating up to 770 GHz,”  IEEE TASC, Vol. 9, No. 2, 1999.
  * [5] Q. P. Herr, A. D. Smith, and M. S. Wire,  “High Speed Data Link between Digital Superconductor Chips,”  Appl. Phys. Lett., Vol. 80, No. 17, 2002.
  * [6] D. S. Holmes, A. L. Ripple, and M. A. Manheimer,  “Energy-Efficient Superconducting computing—Power Budgets and Requirements,”  IEEE TASC, Vol. 23, No. 3, 2013.
  * [7] K. Sternickel and A. I. Braginski,  “Biomagnetism Using SQUIDs: Status and Perspectives,”  Supercond. Sci. Technol., Vol. 19, No. 3, 2006.
  * [8] W. H. Tuttlebee,  Software Defined Radio: Enabling Technologies,  John Wiley and Sons, 2002.
  * [9] G. Krylov and E. G. Friedman,  Single Flux Quantum Integrated Circuit Design,  Springer, 2022.
  * [10] P. Bunyk, K. Likharev, and D. Zinoviev,  “RSFQ Technology: Physics and Devices,”  Int. J. High Speed Electron. Syst., Vol. 11, No. 01, 2001.
  * [11] G. Pasandi and M. Pedram,  “PBMap: A Path Balancing Technology Mapping Algorithm for Single Flux Quantum Logic Circuits,”  IEEE TASC, Vol. 29, No. 4, June 2019.
  * [12] N. Kito, K. Takagi, and N. Takagi,  “Logic-Depth-Aware Technology Mapping Method for RSFQ Logic Circuits With Special RSFQ Gates,”  IEEE TASC, Vol. 32, No. 4, 2021.
  * [13] G. Pasandi and M. Pedram,  “An Efficient Pipelined Architecture for Superconducting Single Flux Quantum Logic Circuits Utilizing Dual Clocks,”  IEEE TASC, Vol. 30, No. 2, 2020.
  * [14] G. Pasandi and M. Pedram,  “Depth-Bounded Graph Partitioning Algorithm and Dual Clocking Method for Realization of Superconducting SFQ Circuits,”  ACM JETC, Vol. 17, No. 1, October 2020.
  * [15] X. Li, M. Pan, T. Liu, and P. A. Beerel,  “Multi-Phase Clocking for Multi-Threaded Gate-Level-Pipelined Superconductive Logic,”  Proc. ISVLSI, pp. 62–67, 2022.
  * [16] R. Bairamkulov and G. De Micheli,  “Compound Logic Gates for Pipeline Depth Minimization in SFQ Integrated Systems,”  Proc. GLSVLSI, 2023.
  * [17] O. Mukhanov, V. Semenov, and K. Likharev,  “Ultimate Performance of the RSFQ Logic Circuits,”  IEEE Trans. Magn., Vol. 23, No. 2, 1987.
  * [18] R. Bairamkulov, A. Tempia Calvino, and G. De Micheli,  “Synthesis of SFQ Circuits with Compound Gates,”  Proc. VLSI-SoC, 2023.
  * [19] C. E. Leiserson and J. B. Saxe,  “Retiming Synchronous Circuitry,”  Algorithmica, Vol. 6, No. 1–6, June 1991.
  * [20] L. Perron,  “Operations Research and Constraint Programming at Google,”  Proc. CP, 2011.
  * [21] M. Soeken _et al._ ,  “The EPFL Logic Synthesis Libraries,”  arXiv Preprint arXiv:1805.05121v3, 2018.
  * [22] M. C. Hansen, H. Yalcin, and J. P. Hayes,  “Unveiling the ISCAS-85 Benchmarks: A Case Study in Reverse Engineering,”  IEEE Des. Test. Comput., Vol. 16, No. 3, 1999.