# Regularized Identification with Internal Positivity Side-Information

Mohammad Khosravi  Roy S. Smith

###### Abstract

In this paper, we present an impulse response identification scheme that
incorporates the internal positivity side-information of the system. The
realization theory of positive systems establishes specific criteria for the
existence of a positive realization for a given transfer function. These
transfer function criteria are translated to a set of suitable conditions on
the shape and structure of the impulse responses of positive systems.
Utilizing these conditions, the impulse response estimation problem is
formulated as a constrained optimization in a reproducing kernel Hilbert space
equipped with a stable kernel, and suitable constraints are imposed to encode
the internal positivity side-information. The optimization problem is
infinite-dimensional with an infinite number of constraints. An equivalent
finite-dimensional convex optimization in the form of a convex quadratic
program is derived. The resulting equivalent reformulation makes the proposed
approach suitable for numerical simulation and practical implementation. A
Monte Carlo numerical experiment evaluates the impact of incorporating the
internal positivity side-information in the proposed identification scheme.
The effectiveness of the proposed method is demonstrated using data from a
heating system experiment.

##  I Introduction

In various dynamical systems, the characteristic variables are constrained to
be non-negative or bounded, either by the nature of their definition or
according to the physics of the underlying system. For instance, charges in
RC-circuits, temperatures and thermal energies in buildings, mass flows in
compartmental systems, the population of certain species of animals or
bacteria, the concentration of pathogens, level of traffic and congestion in
networks and roads, prices of stocks and goods, the pressure of fluids, and
many other quantities of interest are always non-negative [1, 2, 3, 4]. In the
broad sense, a system described only by such non-negative variables is called
a _positive system_ [5]. Depending on our perspective, i.e., whether the
positivity feature is considered as an input-output property or a state-space
characteristic, we have two central notions of positivity in the system theory
literature [5]: internal positivity and external positivity, where the main
focus of the literature is on the former one. In externally positive systems,
the non-negativity of the input signal implies the same feature for the output
signal. Meanwhile, in internally positive systems, the state trajectory and
output signal are non-negative when the initial state and input signal are
non-negative.

Positive systems have received extensive attention in the past decades owing
to being omnipresent in various fields of science and their wide range of
applications [6, 7]. Luenberger pioneered the system theoretic approach to
positive systems with his seminal work [8] in the 1980s. Since then various
subjects of system theory aspects have been tackled, e.g., realization theory
[7], controllability and reachability [9], observability and observer design
[10], robust stability [11, 12], positive stabilization [13, 14], fault
detection and estimation [15], decentralized and distributed control [16, 17],
and, large scale positive systems and scalable control [18, 19].

Concerning the identification problem, there are two aspects when the
underlying system is known to be internally positive. First, respecting the
positivity property can be an essential issue in some applications, such as
implementing model predictive control. Accordingly, any accurate mathematical
modeling approach is expected to include this feature and construct an
internally positive system. The second aspect is informed system
identification [20, 21, 22, 23] and concerns utilizing the internal positivity
side-information and integrating features of this knowledge in the model to
improve the estimation accuracy. Indeed, disregarding the positivity
information can lead to models, which are not physically interpretable and
explainable, or behaviors in contradiction with our expectations [24, 6].
While positive systems have been extensively researched from various
viewpoints [19], their identification problem is not well studied, especially
with regard to these aspects. For instance, a set of conditions is introduced
in [25] for identifying compartmental models, which is a particular case of
internal positivity feature. In [26], assuming the output sequence of data is
a Poisson process, a maximum likelihood approach is presented for third-order
positive systems with distinct real poles. In [24], a particular situation is
considered, where state variable measurements are provided, in addition to
input-output data, and the stability and scalability issues are discussed.
Since internally positive systems are also externally positive, side-
information on internal positivity implies external positivity. Accordingly,
from the perspective of informed system identification, one can consider
external positivity as partial information to be integrated into the model. To
this end, one may employ the external positive system identification methods
[27, 28]. For example, since the externally positive systems are precisely
those with the non-negative impulse response, a kernel-based nonparametric
maximum a posteriori approach is introduced in [29, 28] for estimating a non-
negative finite impulse response (FIR). In this approach, the covariance of
the prior distribution is specified by the stable kernels, while the mean is
designed arbitrarily as an exponentially decaying FIR. Note that in the FIR
identification approaches with an external positivity constraint [27, 29, 28],
the complete information of internal positivity is not exploited. The kernel-
based methods [30], which resolve the issues of bias-variance trade-off,
robustness, and model order selection [31, 32, 33, 34], also provide a
suitable framework for the integration of various sorts of side-information
into the model [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]. These include the
stability of the system, the smoothness of the impulse response, time
constants and resonant frequencies [39]. Together with the realization theory
of positive systems [5], the kernel-based framework can provide a suitable
foundation for impulse response identification of positive systems.

This paper extends our previous work [21] and presents an identification
method that integrates the internal positivity side-information in the
estimated impulse response. From the realization theory of positive systems
[5], we know that the impulse response of an internally positive system has a
specific form, i.e., it has a dominant non-negative part, where the
corresponding transfer function has structured poles, and a residual part.
This specific form can be translated to a set of structural constraints on the
impulse response. Accordingly, the estimation problem is expressed in the form
of a constrained optimization in a stable reproducing kernel Hilbert space,
where suitable constraints are imposed to encode the internal positivity side-
information. Though this problem is initially formulated in an infinite-
dimensional space and with an infinite number of constraints, we derive an
equivalent finite-dimensional convex optimization in the form of a convex
quadratic program. We evaluate the impact of incorporating the internal
positivity side-information and assess the performance of the proposed
identification scheme through a Monte Carlo numerical experiment. The efficacy
of the proposed positive system identification technique is confirmed using
data from a thermal dynamics experiment.

##  II Notations

In this paper, the set of natural numbers, the set of integers, the set of
non-negative integers, the set of real numbers, the set of non-negative real
numbers, $n$-dimensional Euclidean space, and the set of $n$ by $m$ matrices
are respectively shown by ${\mathbb{N}}$, ${\mathbb{Z}}$, ${\mathbb{Z}}_{+}$,
${\mathbb{R}}$, ${\mathbb{R}}_{+}$, ${\mathbb{R}}^{n}$, and
${\mathbb{R}}^{n\times m}$. The positive orthant of ${\mathbb{R}}^{n}$ is
denoted by ${\mathbb{R}}_{+}^{n}$. The identity matrix and zero matrix are
denoted by $\mathbb{I}$ and ${\mathbbold{0}}_{n}$, respectively. Also, the
$n$-dimensional zero vector and the all-ones vector are denoted by $0_{n}$ and
$1_{n}$, respectively. When the dimension is clear from the context, we drop
the subscript. For $p\in[1,\infty)$, the $p$-norm of vector
${\mathrm{h}}=(h_{s})_{s=0}^{\infty}\in{\mathbb{R}}^{{\mathbb{Z}}_{+}}$ is
defined as $\|{\mathrm{h}}\|_{p}=(\sum_{s\geq 0}|h_{s}|^{p})^{\frac{1}{p}}$,
and the $\infty$-norm of ${\mathrm{h}}$ is defined as
$\|{\mathrm{h}}\|_{\infty}=\sup_{s\geq 0}|h_{s}|$. The space of vectors
${\mathrm{h}}$ with finite $p$-norm is denoted by $\ell^{p}$. With respect to
each ${\mathrm{u}}=(u_{s})_{s=0}^{\infty}\in\ell^{\infty}$ and
$t\in{\mathbb{Z}}$, the linear map
${\mathrm{L}}^{\\!{\mathrm{u}}}_{t}:\ell^{1}\to{\mathbb{R}}$ is defined as
${\mathrm{L}}^{\\!{\mathrm{u}}}_{t}({\mathrm{g}})=\sum_{s=0}^{\infty}g_{s}u_{t-s}$,
for any ${\mathrm{g}}=(g_{s})_{s=0}^{\infty}\in\ell^{1}$. Given a subset
${\mathcal{C}}\subset{\mathcal{X}}$, the function
$\delta_{{\mathcal{C}}}:{\mathcal{X}}\to\\{0,+\infty\\}$ is defined as
$\delta_{{\mathcal{C}}}(x)=0$, if $x\in{\mathcal{C}}$ and
$\delta_{{\mathcal{C}}}(x)=\infty$, otherwise. The set of polynomials in $x$
with maximum degree $n$ and real coefficients is denoted by
${\mathbb{R}}_{n}[x]$. For transfer function $G$, $r(G)$ denotes its spectral
radius.

##  III System Identification with Internal Positivity Side-Information

Let
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}:=({g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{t})_{t=0}^{\infty}$
be the impulse response of stable and causal system ${\mathcal{S}}$ and
$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z):=\sum_{t=0}^{\infty}{g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{t}z^{-t}$
be the corresponding transfer function. We call impulse response
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$, or equivalently the
system ${\mathcal{S}}$, _internally positive_ if there exists a realization
such that the state trajectory and the output remains non-negative given that
the initial state and the input are non-negative (see Definition 2). Suppose a
bounded signal ${\mathrm{u}}\in\ell^{\infty}$ is applied to the input of
system ${\mathcal{S}}$. Let $y_{t}$ denote the measured output at time instant
$t\in{\mathscr{T}}$, where ${\mathscr{T}}:=\\{t_{i}\ \\!|\\!\
i=0,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}\\!-\\!1\\}$, for a given
$n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}\in{\mathbb{N}}$. In other words, we
have

$y_{t}:={\mathrm{L}}^{\\!{\mathrm{u}}}_{t}({\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}})+w_{t},\qquad
t\in{\mathscr{T}},$ (1)

where $w_{t}$ denotes the uncertainty in the output measured at time instant
$t$, for $t\in{\mathscr{T}}$. Accordingly, we have a set of input-output
measurement data denoted be ${\mathscr{D}}$. Based on the introduced setting,
we introduce the following impulse response identification problem.

######  Problem.

Using data ${\mathscr{D}}$, estimate the impulse response of
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$, given the side-
information that ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is
internally positive.

In addressing this problem, the main concern is the appropriate integration of
the available internal positivity side-information into the impulse response
identification problem. To this end, we need to exploit suitable conditions
inducing the desired positivity feature. In the next section, these conditions
are discussed, and the estimation problem is formulated accordingly.

##  IV Mathematical Formulation for Internally Positive System Identification

In the realization theory of positive systems, sufficient conditions are
introduced under which the transfer function of a system admits a so-called
positive realization. We employ these conditions together with the notion of
stable reproducing kernel Hilbert spaces (RKHS) for bridging to the impulse
response identification of positive systems.

The following definition introduces the external positivity notion for impulse
response ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$, or equivalently
for system ${\mathcal{S}}$.

######  Definition 1 ([5]).

The impulse response ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is
said to be externally positive if for each $\tau\in{\mathbb{Z}}$ and for any
input signal ${\mathrm{u}}=(u_{t})_{t\in{\mathbb{Z}}}$, we know that the
(noiseless) output of the system is non-negative for all $t\leq\tau$, when we
have $u_{t}\geq 0$, for each $t\leq\tau$. The set of externally positive
impulse responses is denoted by $\overline{{\mathscr{P}}}$.

Since the condition introduced in Definition 1 is imposed for all
$\tau\in{\mathbb{Z}}$, one can see that externally positivity is equivalent to
the non-negativity of input signals implying the non-negativity of the output
signal of the system, i.e., if $u_{t}\geq 0$, for all $t\in{\mathbb{Z}}$, then
$y_{t}={\mathrm{L}}^{\\!{\mathrm{u}}}_{t}({\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}})\geq
0$, for each $t\in{\mathbb{Z}}$. The next theorem introduces a necessary and
sufficient condition for external positivity.

######  Theorem 1 ([5]).

The impulse response ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is
externally positive if and only if $g_{t}\geq 0$, for any
$t\in{\mathbb{Z}}_{+}$.

From Theorem 1, one can see that the external positive stable impulse
responses are exactly the non-negative ones, i.e., we have
$\overline{{\mathscr{P}}}=\ell^{1}\cap{\mathbb{R}}_{+}^{{\mathbb{Z}}_{+}}$.
Along with the external positivity notion, we have another notion of
positivity which is the main concern of the current paper and introduced in
the following definition.

######  Definition 2 ([5]).

The impulse response ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is
said to be internally positive, or simply positive, if there exists a
realization for system ${\mathcal{S}}$ as

$\left\\{\begin{array}[]{ccl}{\mathrm{x}}_{t+1}\\!\\!&\\!\\!\\!=\\!\\!\\!&\\!\\!{\mathrm{A}}{\mathrm{x}}_{t}+{\mathrm{b}}u_{t},\\\
y_{t}\\!\\!&\\!\\!\\!=\\!\\!\\!&\\!\\!{\mathrm{c}}{\mathrm{x}}_{t}+du_{t},\\\
\end{array}\right.\qquad\forall\ t\in{\mathbb{Z}},$ (2)

where ${\mathrm{A}}\in{\mathbb{R}}^{{n_{\mathrm{x}}}\times{n_{\mathrm{x}}}}$,
${\mathrm{b}}\in{\mathbb{R}}^{{n_{\mathrm{x}}}}$,
${\mathrm{c}}\in{\mathbb{R}}^{1\times{n_{\mathrm{x}}}}$, $d\in{\mathbb{R}}$
and ${n_{\mathrm{x}}}\in{\mathbb{N}}$, such that
${\mathrm{x}}_{0}\in{\mathbb{R}}_{+}^{{n_{\mathrm{x}}}}$ and
$u_{t}\in{\mathbb{R}}_{+}$, for all $t\geq 0$, implies that
${\mathrm{x}}_{t}\in{\mathbb{R}}_{+}^{{n_{\mathrm{x}}}}$ and
$y_{t}\in{\mathbb{R}}_{+}$, for each $t\geq 0$. The realization (2) with the
this property is called a _positive realization_ for
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$, or equivalently, for
$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$. Moreover, the set of internally
positive impulse responses is denoted by ${\mathscr{P}}$.

The internal positivity enforces a specific attribute on
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ according to the
_Kronecker’s theorem_ given below.

######  Theorem 2 ([46]).

With respect to impulse response ${\mathrm{g}}=(g_{t})_{t=0}^{\infty}$, define
Hankel operator $\mathrm{Hankel}({\mathrm{g}}):\ell^{\infty}\to\ell^{\infty}$
with entrywise representation in the standard basis of $\ell^{\infty}$ as
following

$\mathrm{Hankel}({\mathrm{g}}):=\begin{bmatrix}g_{0}&g_{1}&g_{2}&\ldots\\\
g_{1}&g_{2}&\reflectbox{$\ddots$}&\\\ g_{2}&\reflectbox{$\ddots$}&&\\\
\vdots&&&\\\
\end{bmatrix}=\begin{bmatrix}g_{i+j-2}\end{bmatrix}_{i,j=1}^{\infty}.$ (3)

Then, $G(z):=\sum_{t=0}^{\infty}g_{t}z^{-t}$ is a rational function if and
only if the rank of $\mathrm{Hankel}({\mathrm{g}})$ is finite, i.e., we have

$\mathrm{rank}\big{(}\mathrm{Hankel}({\mathrm{g}})\big{)}=\dim\Big{\\{}\mathrm{Hankel}({\mathrm{g}}){\mathrm{v}}\
\\!\big{|}\ \\!{\mathrm{v}}\in\ell^{\infty}\Big{\\}}<\infty.$

We call impulse response ${\mathrm{g}}$ _finite Hankel rank_ when the property
above is satisfied.

According to (2), we know that
$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)={\mathrm{c}}(z\mathbb{I}-{\mathrm{A}})^{-1}{\mathrm{b}}+d$
is a rational function. Therefore, due to Theorem 2, the internal positivity
of ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ implies that

$\mathrm{rank}\big{(}\mathrm{Hankel}({\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}})\big{)}<\infty.$
(4)

The realization (2) needs to have special structure introduced in the next
theorem.

######  Theorem 3 ([5]).

The impulse response ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is
internally positive if and only if there exist a realization as in (2) such
that the entries of ${\mathrm{A}},{\mathrm{b}},{\mathrm{c}}$ and $d$ are non-
negative.

The next corollary is easily concluded from Theorem 3 for the internally
positive impulse responses.

######  Theorem 4 ([5]).

If impulse response ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is
internally positive then we have

${g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{t}\geq 0,\qquad\forall
t\in{\mathbb{Z}}_{+}.$ (5)

Theorem 4 implies that any internally positive impulse response is externally
positive as well. Accordingly, due to (4), we have
${\mathscr{P}}\subset\underline{{\mathscr{P}}}$, where the set of impulse
responses $\underline{{\mathscr{P}}}$ is defined as

$\underline{{\mathscr{P}}}:=\overline{{\mathscr{P}}}\cap\Big{\\{}{\mathrm{g}}\\!\in\\!\ell^{1}\\!\
\big{|}\\!\ \mathrm{rank}(\mathrm{Hankel}({\mathrm{g}}))<\infty\Big{\\}}.$ (6)

Meanwhile, from the next example, which is a modified version of an example
given in [7], one can see that the inclusion in
${\mathscr{P}}\subset\underline{{\mathscr{P}}}$ is strict, i.e.,
${\mathscr{P}}\neq\underline{{\mathscr{P}}}$.

######  Example.

Let ${{\rho}}\in(0,1)$ and $\omega$ be an irrational real number. Define the
impulse response ${\mathrm{g}}=(g_{t})_{t=0}^{\infty}$ as

$g_{t}={{\rho}}^{t}(1+\cos(2\pi\omega t)),\qquad\forall t\in{\mathbb{Z}}_{+}.$
(7)

One can see that ${\mathrm{g}}$ is a non-negative impulse response with
following transfer function

$G(z)=\frac{1}{1-{{\rho}}z^{-1}}+\frac{1-{{\rho}}\cos w\ z^{-1}}{1-2\rho\cos
w\ z^{-1}+\rho^{2}z^{-2}}.$ (8)

Therefore, we have ${\mathrm{g}}\in\underline{{\mathscr{P}}}$. However, there
is no positive realization for the impulse response ${\mathrm{g}}$ [5].
Accordingly, due to Definition 2, we know that ${\mathrm{g}}$ is not
internally positive, i.e., ${\mathrm{g}}\notin{\mathscr{P}}$ and
${\mathscr{P}}\neq\underline{{\mathscr{P}}}$.

Based on the above discussion, conditions (4) and (5) are necessary but not
sufficient for ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ being an
internally positive impulse response. Using the following theorem, we derive
sufficient conditions for the internal positivity of ${\mathrm{g}}$ which are
used later to formulate the identification problem of internally positive
systems.

######  Theorem 5 ([5]).

Let ${\mathrm{g}}$ be a non-negative impulse response and $G$ be the
corresponding transfer function If $G$ is a strictly proper rational function
with a unique dominant pole $\rho\in(0,1)$, then there exists a positive
realization for $G$ .

With respect to each $\rho\in(0,1)$, define
${\mathscr{P}}_{\\!\rho}\subset\ell^{1}$ as the set of non-negative impulse
responses satisfying (4) such that we have

$\exists\ \\!{a}\\!\in\\!(0,\infty),\quad\lim_{t\to\infty}\rho^{-t}g_{t}={a},$
(9)

i.e., $\lim_{t\to\infty}\rho^{-t}g_{t}$ is well-defined and equal to a
positive real scalar ${a}$. Furthermore, we define ${\mathscr{P}}_{(0,1)}$ as
${\mathscr{P}}_{\\!(0,1)}=\cup_{\rho\in(0,1)}{\mathscr{P}}_{\\!\rho}$. Based
on Theorem 5, we have the following corollary for ${\mathscr{P}}_{\\!\rho}$
and ${\mathscr{P}}_{\\!(0,1)}$.

######  Corollary 6.

For any $\rho\in(0,1)$, each impulse response
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}\in{\mathscr{P}}_{\\!\rho}$
is internally positive, i.e., ${\mathscr{P}}_{\\!\rho}\subset{\mathscr{P}}$.
Moreover, we have ${\mathscr{P}}_{\\!(0,1)}\subset{\mathscr{P}}$.

###### Proof.

See Appendix A-A. ∎

Note that ${\mathscr{P}}_{\\!(0,1)}$ contains exactly the impulse responses
satisfying conditions (4), (5), and, (9). Hence, Corollary 6 says that any
impulse response in $\ell^{1}$ which satisfies these conditions is internally
positive. Accordingly, one can employ (4), (5) and (9) in the identification
problem to enforce internal positivity on the impulse response to be
estimated. The next theorem further highlights the importance of positive
systems ${\mathscr{P}}_{\\!(0,1)}$.

######  Theorem 7.

The set of impulse responses ${\mathscr{P}}_{\\!(0,1)}$ is dense in
${\mathscr{P}}$ with respect to $p$-norm topology, for any $p\in[1,\infty]$.

###### Proof.

See Appendix A-B. ∎

With respect to each
${\mathrm{g}}=(g_{t})_{t=0}^{\infty}\in{\mathscr{P}}_{(0,1)}$, one can define
impulse response ${\mathrm{h}}=(h_{t})_{t=0}^{\infty}$ such that
$h_{t}=g_{t}-{a}\rho^{t}$, for $t\in{\mathbb{Z}}_{+}$, where $\rho$ and ${a}$
are the positive scalars introduced in (9). Note that
$\lim_{t\to\infty}\rho^{-t}h_{t}=0$ and $\rho\in(0,1)$, which implies that
${\mathrm{h}}=(h_{t})_{t=0}^{\infty}$ is a stable impulse response dominated
by

${\mathrm{f}}_{\rho}=(f_{t})_{t=0}^{\infty}:=(\rho^{t})_{t=0}^{\infty}.$ (10)

Since, for each $t$, we have $g_{t}={a}\rho^{t}+h_{t}$, to identify the
internally positive impulse response ${\mathrm{g}}=(g_{t})_{t=0}^{\infty}$, we
need to estimate $\rho,{a}$ and the stable impulse response
${\mathrm{h}}=(h_{t})_{t=0}^{\infty}$, dominated by
${\mathrm{f}}_{\rho}=(\rho^{t})_{t=0}^{\infty}$, and meanwhile ensure that
${\mathrm{g}}=(g_{t})_{t=0}^{\infty}$ satisfies properties (4) and (5).
Accordingly, we need a suitable hypothesis space for
${\mathrm{h}}=(h_{t})_{t=0}^{\infty}$. To this end, we employ stable
reproducing kernel Hilbert spaces introduced below.

######  Definition 3 ([47, 48]).

The non-zero symmetric function
$\mathbbm{k}:{\mathbb{Z}}_{+}\times{\mathbb{Z}}_{+}\to{\mathbb{R}}$ is said to
be a _Mercer kernel_ if, for any $m\in{\mathbb{N}}$,
$t_{1},\ldots,t_{m}\in{\mathbb{T}}$ and $a_{1},\ldots,a_{m}\in{\mathbb{R}}$,
we have
$\sum_{i=1}^{m}\\!\sum_{j=1}^{m}\\!a_{i}\mathbbm{k}(t_{i},t_{j})a_{j}\geq 0$.
Moreover, the _section_ of kernel $\mathbbm{k}$ at $t\in{\mathbb{Z}}_{+}$ is
denoted by $\mathbbm{k}_{t}$ and defined as the function
$\mathbbm{k}(t,\cdot):{\mathbb{Z}}_{+}\to{\mathbb{R}}$. Furthermore, the
positive kernel $\mathbbm{k}$ is said to be _stable_ if, for any
${\mathrm{u}}=(u_{t})_{t\in{\mathbb{Z}}_{+}}\in\ell^{\infty}$, we have
$\sum_{t\in{\mathbb{Z}}_{+}}\\!\\!|\sum_{s\in{\mathbb{Z}}_{+}}\\!\\!u_{s}\mathbbm{k}(t,s)|<\infty$.

######  Theorem 8 ([47, 48]).

Given a Mercer kernel
$\mathbbm{k}:{\mathbb{Z}}_{+}\times{\mathbb{Z}}_{+}\to{\mathbb{R}}$, there
exists a _unique_ Hilbert space
${\mathcal{H}}_{\mathbbm{k}}\subseteq{\mathbb{R}}^{{\mathbb{Z}}_{+}}$ endowed
with inner product
${\langle{\cdot,\cdot}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}$ and norm
$\|\cdot\|_{{\mathcal{H}}_{\mathbbm{k}}}$, called a _RKHS with kernel_
$\mathbbm{k}$, such that, for each $t\in{\mathbb{Z}}_{+}$, we have

  * i)

$\mathbbm{k}_{t}\in{\mathcal{H}}_{\mathbbm{k}}$, and

  * ii)

${\langle{{\mathrm{g}},\mathbbm{k}_{t}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}=g_{t}$,
for all
${\mathrm{g}}=(g_{t})_{t\in{\mathbb{Z}}_{+}}\in{\mathcal{H}}_{\mathbbm{k}}$.

The second feature is called reproducing property. Moreover,
${\mathcal{H}}_{\mathbbm{k}}\subset\ell^{1}$ if and only if $\mathbbm{k}$ is a
stable kernel. In this case, ${\mathcal{H}}_{\mathbbm{k}}$ is said to be a
_stable RKHS_.

Given a stable kernel $\mathbbm{k}$, we take ${\mathcal{H}}_{\mathbbm{k}}$ as
the hypothesis space for the stable impulse response ${\mathrm{h}}$.
Considering the set of input-output data ${\mathscr{D}}$, we define the
empirical loss function
${\mathcal{E}}_{\rho}:{\mathbb{R}}\times{\mathcal{H}}_{\mathbbm{k}}\to{\mathbb{R}}_{+}$
as

${\mathcal{E}}_{\rho}({a},{\mathrm{h}}):=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\Big{(}y_{t_{i}}-{a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})\Big{)}^{2},$
(11)

where we assume the hyperparameter $\rho\in(0,1)$ is given. The estimation of
$\rho$ is will be discussed later. We formulate the identification problem
with internal positivity side-information as following regularized
optimization

$\begin{array}[]{cl}\operatornamewithlimits{min}\limits_{{a}\in{\mathbb{R}},{\mathrm{h}}\in{\mathcal{H}}_{\mathbbm{k}}}&{\mathcal{E}}_{\rho}({a},{\mathrm{h}})+\lambda\\!\
\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2},\\\
\text{s.t.}&h_{t}+{a}\rho^{t}\geq 0,\quad\forall t\geq 0,\\\
&\mathrm{rank}(\mathrm{Hankel}({\mathrm{h}}))<\infty,\\\
&{a}\geq{a}_{\min},\end{array}$ (12)

where ${a}_{\min}>0$ is a given lower-bound for ${a}$ to ensure that ${a}>0$,
and $\lambda>0$ is the regularization weight. Note that, similarly to the
standard problem formulation in the literature on the kernel-based impulse
response identification [33], the objective function in (12) is an empirical
loss function regularized with the RKHS norm of ${\mathrm{h}}$. This ensures
the stability of ${\mathrm{h}}$ and also allows incorporating other features
such as exponential decay and smoothness [39]. The next theorem says that the
solution of (12) leads to an internally positive estimation of impulse
response ${\mathrm{g}}$. Before proceeding to the theorem, we need to
introduce an assumption.

######  Assumption 1.

There exist $C\in{\mathbb{R}}_{+}$ and
$\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}\in(0,\rho)$ such that we have
$|\mathbbm{k}(t,t)|\leq C\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{2t}$, for any
$t\in{\mathbb{Z}}_{+}$.

######  Theorem 9.

Let Assumption 1 hold, ${a}$ and ${\mathrm{h}}=(h_{t})_{t=1}^{\infty}$ be a
solution pair for (12), and, the impulse response
${\mathrm{g}}=(g_{t})_{t=0}^{\infty}$ be defined as $g_{t}=h_{t}+{a}\rho^{t}$,
for any $t\in{\mathbb{Z}}_{+}$. Then, ${\mathrm{g}}$ is internally positive.

###### Proof.

Let $G$ be the transfer function which corresponds to ${\mathrm{g}}$. Due to
(12), the rank of Hankel operator $\mathrm{Hankel}({\mathrm{h}})$ is finite.
Subsequently, according to Theorem 2, the transfer function corresponding to
${\mathrm{h}}$, $H$, has finite order. On the other hand, we know that

$G(z)=\frac{{a}z^{-1}}{1-\rho z^{-1}}+H(z).$ (13)

Therefor, the order of $G$ is finite. Accordingly, due to Theorem 2, we know
that $\mathrm{rank}(\mathrm{Hankel}({\mathrm{g}}))<\infty$, i.e.,
${\mathrm{g}}$ satisfies (4). Furthermore, according to the first constraint
in (12), one can see that ${\mathrm{g}}$ is a non-negative impulse response
and (5) holds for ${\mathrm{g}}$. Moreover, from the reproducing property of
kernel, we know that $h_{t}={\langle{\mathbbm{k}_{t},{\mathrm{h}}}\rangle}$
and
$\|\mathbbm{k}_{t}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2}={\langle{\mathbbm{k}_{t},\mathbbm{k}_{t}}\rangle}=\mathbbm{k}(t,t)$,
for any $t\in{\mathbb{Z}}_{+}$. Consequently, due to Cauchy-Schwartz
inequality and Assumption 1, we have

$\begin{split}\\!\\!\\!\\!\\!|h_{t}|=|{\langle{\mathbbm{k}_{t},{\mathrm{h}}}\rangle}|&\leq\|\mathbbm{k}_{t}\|_{{\mathcal{H}}_{\mathbbm{k}}}\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}\\\
&\qquad\
=\mathbbm{k}(t,t)^{\frac{1}{2}}\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}\leq
C^{\frac{1}{2}}\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{t}.\end{split}$
(14)

for any $t\in{\mathbb{Z}}_{+}$. Following this, one can see that

$\begin{split}0&\leq\operatornamewithlimits{liminf}\limits_{t\to\infty}\rho^{-t}h_{t}\leq\operatornamewithlimits{limsup}\limits_{t\to\infty}\rho^{-t}h_{t}\leq\operatornamewithlimits{limsup}\limits_{t\to\infty}\rho^{-t}|h_{t}|\\\
&\leq\operatornamewithlimits{limsup}\limits_{t\to\infty}C^{\frac{1}{2}}\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{t}\rho^{-t}=0,\end{split}$
(15)

where the last equality is due to
$\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}\in(0,\rho)$. Hence,
$\lim_{t\to\infty}\rho^{-t}h_{t}$ is well-defined and we have
$\lim_{t\to\infty}\rho^{-t}h_{t}=0$. Subsequently, due to the definition of
${\mathrm{g}}$, it follows that $\lim_{t\to\infty}\rho^{-t}g_{t}={a}$ and
${\mathrm{g}}$ satisfies (9). Therefore, ${\mathrm{g}}$ belongs to
${\mathscr{P}}_{\\!\rho}$, and consequently, due to Corollary 6,
${\mathrm{g}}$ is internally positive. ∎

######  Remark 1.

For $a$, ${\mathrm{h}}$ and ${\mathrm{g}}$ introduced in Theorem 9, we have

${\mathcal{E}}_{\rho}(a,{\mathrm{h}})=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\big{(}y_{t_{i}}-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{g}})\big{)}^{2},$
(16)

i.e., in the cost function of (12), the first term is the sum of squared
errors for the impulse response fitting when the dominant pole $\rho$ is
known.

##  V Towards a Tractable Solution

In this section, we investigate optimization problem (12) which was introduced
for impulse response identification with internal positivity. This
optimization problem is in an infinite-dimensional space with an infinite
number of constraints. In the following, we analyze this problem and provide a
tractable approach for deriving its solution.

Let ${\mathcal{V}}_{\mathbbm{k}}$ be the Hilbert space
${\mathbb{R}}\times{\mathcal{H}}_{\mathbbm{k}}$ which is endowed with inner
product
${\langle{\cdot,\cdot}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}:{\mathcal{V}}_{\mathbbm{k}}\times{\mathcal{V}}_{\mathbbm{k}}\to{\mathbb{R}}$
defined as following

${\langle{({a}_{1},{\mathrm{h}}_{1}),({a}_{2},{\mathrm{h}}_{2})}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}={a}_{1}{a}_{2}+{\langle{{\mathrm{h}}_{1},{\mathrm{h}}_{2}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}},$
(17)

for any ${a}_{1},{a}_{2}\in{\mathbb{R}}$ and
${\mathrm{h}}_{1},{\mathrm{h}}_{2}\in{\mathcal{H}}_{\mathbbm{k}}$. Also, let
${\mathscr{F}}\subseteq{\mathcal{H}}_{\mathbbm{k}}$ be the set of finite
Hankel rank impulse responses in ${\mathcal{H}}_{\mathbbm{k}}$, i.e.,

${\mathscr{F}}=\big{\\{}{\mathrm{h}}\in{\mathcal{H}}_{\mathbbm{k}}\\!\
\big{|}\\!\ \mathrm{rank}(\mathrm{Hankel}({\mathrm{h}}))<\infty\big{\\}}.$
(18)

We define function
${\mathcal{J}}_{\\!{\mathscr{F}}}:{\mathcal{V}}_{\mathbbm{k}}\to\mathbb{R}\cup\\{+\infty\\}$
as

${\mathcal{J}}_{\\!{\mathscr{F}}}({a},{\mathrm{h}})={\mathcal{E}}_{\rho}({a},{\mathrm{h}})+\sum_{s=0}^{\infty}\delta_{{{\mathscr{R}}}_{s}}({a},{\mathrm{h}})+\delta_{{\mathscr{F}}}({\mathrm{h}})+\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2},$
(19)

where ${{\mathscr{R}}}_{s}\subseteq{\mathcal{V}}_{\mathbbm{k}}$ is the
following set

${{\mathscr{R}}}_{s}:=\Big{\\{}\big{(}{a},(h_{s})_{s\in{\mathbb{Z}}_{+}}\\!\big{)}\in{\mathcal{V}}_{\mathbbm{k}}\\!\
\Big{|}\\!\ h_{s}+{a}\rho^{s}\geq 0,{a}\geq{a}_{\min}\Big{\\}},$ (20)

for $s\in{\mathbb{Z}}_{+}$. From the definition of
${\mathcal{J}}_{\\!{\mathscr{F}}}$, it follows easily that the optimization
problem (12) is equivalent to

$\inf_{({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}}\
\\!{\mathcal{J}}_{\\!{\mathscr{F}}}({a},{\mathrm{h}}).$ (21)

For $({a},{\mathrm{h}})=({a}_{\min},\mathbf{0})$, where $\mathbf{0}$ denotes
the zero vector in ${\mathcal{H}}_{\mathbbm{k}}$, one can easily see that

${\mathcal{J}}_{\\!{\mathscr{F}}}({a}_{\min},\mathbf{0})=\sum_{i=1}^{n}\big{(}y_{t_{i}}-{a}_{\min}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})\big{)}^{2}<\infty.$
(22)

Since, for any $({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}$, we have
${\mathcal{J}}_{\\!{\mathscr{F}}}({a},{\mathrm{h}})\geq 0$, it follows that
(21) is bounded. However, this arguments does not guarantee the existence of
solution for (21). In the following, we show that under mild conditions the
optimization problem (21) admits a solution when the kernel $\mathbbm{k}$
meets certain criteria. Let function
${\mathcal{J}}:{\mathcal{V}}_{\mathbbm{k}}\to\mathbb{R}\cup\\{+\infty\\}$ be
defined as following

${\mathcal{J}}({a},{\mathrm{h}})={\mathcal{E}}_{\rho}({a},{\mathrm{h}})+\sum_{s=0}^{\infty}\delta_{{{\mathscr{R}}}_{s}}({a},{\mathrm{h}})+\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2},$
(23)

and consider the optimization problem

$\inf_{({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}}\
\\!{\mathcal{J}}({a},{\mathrm{h}}).$ (24)

One can easily see that
${\mathcal{J}}_{\\!{\mathscr{F}}}={\mathcal{J}}+\delta_{{\mathscr{F}}}$, which
implies that

${\mathcal{J}}({a},{\mathrm{h}})\leq{\mathcal{J}}_{\\!{\mathscr{F}}}({a},{\mathrm{h}}),\qquad\forall\\!\
({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}.$ (25)

Consequently, if (24) has a solution $({a}^{*},{\mathrm{h}}^{*})$ such that
the operator $\mathrm{Hankel}({\mathrm{h}}^{*})$ is finite-rank, then
$({a}^{*},{\mathrm{h}}^{*})$ is a solution for (21) as well. In other words,
the identification problem with internal positivity side-information
introduced in (12) admits a solution. Hence, we need to study the solution
behavior of (24). To this end, we require several technical assumptions.

######  Assumption 2.

There exists $i\in\\{0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1\\}$
such that ${\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})\neq 0$.

Note that if Assumption 2 does not hold, then, for all
$i\in\\{0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1\\}$, we have
${\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})=0$, which means
that the dominant pole is not excited by the input signal ${\mathrm{u}}$.
Accordingly, this assumption essentially says that the input signal excites
the dominant pole.

######  Assumption 3.

There exists $\underline{t}\leq 0$ such that $u_{t}=0$, for $t<\underline{t}$.

Assumption 3 is a technical assumption and introduced mainly for the sake of
our mathematical arguments, i.e., to guarantee the continuity of operator
${\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}:{\mathcal{H}}_{\mathbbm{k}}\to{\mathbb{R}}$,
for $i\in\\{0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1\\}$. Indeed,
this assumption holds in realistic situations such as when the system is
initially at rest. Based on these assumptions, we can show the existence and
uniqueness for the solution of (24).

######  Theorem 10.

Under Assumptions 2 and 3, optimization problem (24) admits a unique solution,
i.e., there exists $({a}^{*},{\mathrm{h}}^{*})\in{\mathcal{V}}_{\mathbbm{k}}$
such that

${\mathcal{J}}({a}^{*},{\mathrm{h}}^{*})<{\mathcal{J}}({a},{\mathrm{h}}),$
(26)

for any
$({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}\backslash\big{\\{}({a}^{*},{\mathrm{h}}^{*})\big{\\}}$.

###### Proof.

Let set ${\mathscr{R}}\subset{\mathcal{V}}_{\mathbbm{k}}$ be defined as
${{\mathscr{R}}}=\bigcap_{s=0}^{\infty}{{\mathscr{R}}}_{s}$. Accordingly, one
can see that
$\sum_{s=0}^{\infty}\delta_{{{\mathscr{R}}}_{s}}=\delta_{{{\mathscr{R}}}}$,
and hence, we have

${\mathcal{J}}({a},{\mathrm{h}})={\mathcal{E}}_{\rho}({a},{\mathrm{h}})+\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2}+\delta_{{{\mathscr{R}}}}({a},{\mathrm{h}}).$
(27)

With respect to each $s\in{\mathbb{Z}}_{+}$, define set
${\mathcal{Q}}_{s}\subset{\mathcal{V}}_{\mathbbm{k}}$ as

${\mathcal{Q}}_{s}:=\Big{\\{}\big{(}{a},(h_{s})_{s\in{\mathbb{Z}}_{+}}\\!\big{)}\in{\mathcal{V}}_{\mathbbm{k}}\\!\
\Big{|}\\!\ h_{s}+{a}\rho^{s}\geq 0\Big{\\}}.$ (28)

For any $s\in{\mathbb{Z}}_{+}$, due to the reproducing property of the kernel,
we have

$\begin{split}h_{s}+{a}\rho^{s}&={\langle{{\mathrm{h}},\mathbbm{k}_{s}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}+{a}\rho^{s}={\langle{({a},{\mathrm{h}}),(\rho^{s},\mathbbm{k}_{s})}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}.\end{split}$
(29)

Therefore, we know that ${\mathcal{Q}}_{s}\subset{\mathcal{V}}_{\mathbbm{k}}$
is a half-space, and hence, it is a non-empty, closed and convex subset of
${\mathcal{V}}_{\mathbbm{k}}$, for all $s\in{\mathbb{Z}}_{+}$. Note that
$[{a}_{\min},\infty)\times{\mathcal{H}}_{\mathbbm{k}}$ is a non-empty, closed
and convex subsets of ${\mathcal{V}}_{\mathbbm{k}}$. One can see that

${{\mathscr{R}}}=\big{(}\operatornamewithlimits{\text{\scalebox{1.25}{$\cap$}}}\limits_{s=0}^{\infty}{\mathcal{Q}}_{s}\big{)}\\!\
\cap\\!\ \big{(}[{a}_{\min},\infty)\times{\mathcal{H}}_{\mathbbm{k}}\big{)},$
(30)

and also, we know that $({a}_{\min},\mathbf{0})$ belongs to
$[{a}_{\min},\infty)\times{\mathcal{H}}_{\mathbbm{k}}$ and
${\mathcal{Q}}_{s}$, for each $s\in{\mathbb{Z}}_{+}$. Therefore,
${{\mathscr{R}}}$ is a non-empty, closed and convex subset of
${\mathcal{V}}_{\mathbbm{k}}$. Consequently, it follows that
$\delta_{{{\mathscr{R}}}}:{\mathcal{V}}\to{\mathbb{R}}\cup\\{+\infty\\}$ is a
proper, convex and lower semi-continuous function, where we have
${\delta_{{{\mathscr{R}}}}({a}_{\min},\mathbf{0})=0}$. With respect to each
$i\in\\{0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1\\}$, define
$\varphi_{i}$ as

$\varphi_{i}:=\sum_{s=0}^{t_{i}-\underline{t}}u_{t_{i}-s}\mathbbm{k}_{s}=u_{t_{i}}\mathbbm{k}_{0}+u_{t_{i}-1}\mathbbm{k}_{1}+\ldots+u_{\underline{t}}\mathbbm{k}_{t_{i}-\underline{t}}.$
(31)

Since ${\mathcal{H}}_{\mathbbm{k}}$ is a linear space which contains the
sections of the kernel, we know that
$\varphi_{i}\in{\mathcal{H}}_{\mathbbm{k}}$, for each
$i\in\\{0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1\\}$. Moreover,
from Assumption 3, reproducing property of kernel and the linearity property
of inner product, it follows that

$\begin{split}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})&=\sum_{s=0}^{t_{i}-\underline{t}}{\langle{{\mathrm{h}},\mathbbm{k}_{s}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}u_{t_{i}-s}\\\
&=\langle{\mathrm{h}},\sum_{s=0}^{t_{i}-\underline{t}}\mathbbm{k}_{s}u_{t_{i}-s}\rangle_{{\mathcal{H}}_{\mathbbm{k}}}={\langle{{\mathrm{h}},\varphi_{i}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}.\end{split}$
(32)

Accordingly, we have

${a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}+{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})={a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}+{\langle{{\mathrm{h}},\psi_{i}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}={\langle{({a},{\mathrm{h}}),\psi_{i}}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}},$
(33)

where $\psi_{i}\in{\mathcal{V}}_{\mathbbm{k}}$ is the vector defined as
$\psi_{i}=({\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho}),\varphi_{i})$,
for $i=0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1$. Therefore, one
can see that

${\mathcal{J}}({a},{\mathrm{h}})\\!=\\!\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\\!\big{(}y_{t_{i}}\\!-\\!{\langle{({a},{\mathrm{h}}),\psi_{i}}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}\big{)}^{2}+\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2}\\!+\\!\delta_{{{\mathscr{R}}}}({a},{\mathrm{h}}),$

Accordingly, since $\delta_{{{\mathscr{R}}}}$ is proper, convex and lower
semi-continuous, and also, due to the fact that

${\mathcal{J}}({a}_{\min},\mathbf{0})=\sum_{i=1}^{n}\big{(}y_{t_{i}}-{a}_{\min}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})\big{)}^{2}<\infty,$
(34)

we know that ${\mathcal{J}}$ is a proper, convex and lower semi-continuous
function. This implies that (24) has a solution [49]. Define bilinear operator
${\mathrm{Q}}:{\mathcal{V}}_{\mathbbm{k}}\times{\mathcal{V}}_{\mathbbm{k}}\to{\mathbb{R}}$
as

$\begin{split}{\mathrm{Q}}\big{(}&({a}_{1},{\mathrm{h}}_{1}),({a}_{2},{\mathrm{h}}_{2})\big{)}\\\
=&\\!\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}{\langle{({a}_{1},{\mathrm{h}}_{1}),\psi_{i}}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}{\langle{\psi_{i},({a}_{2},{\mathrm{h}}_{2})}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}\\!+\\!\lambda{\langle{{\mathrm{h}}_{1},{\mathrm{h}}_{2}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}.\end{split}$
(35)

For any $({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}$, one can easily
that

${\mathrm{Q}}\big{(}({a},{\mathrm{h}}),({a},{\mathrm{h}})\big{)}=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}{\langle{({a},{\mathrm{h}}),\psi_{i}}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}^{2}\\!+\\!\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2}\geq
0.$ (36)

Moreover, since $\lambda$ is a positive real scalar, if
${\mathrm{Q}}\big{(}({a},{\mathrm{h}}),({a},{\mathrm{h}})\big{)}=0$, then, we
need to have ${\mathrm{h}}=0$ and
${\langle{({a},{\mathrm{h}}),\psi_{i}}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}=0$,
for all $i=0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1$, which
implies that
${a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})=0$.
Subsequently, due to Assumption 2, we have ${a}=0$, i.e.,
$({a},{\mathrm{h}})=(0,\mathbf{0})$. Based on this argument, we know that
${\mathrm{Q}}$ is a positive definite bilinear operator. Therefore, the
function $f:{\mathcal{V}}_{\mathbbm{k}}\to{\mathbb{R}}$, defined as
$f({\mathrm{v}})={\mathrm{Q}}({\mathrm{v}},{\mathrm{v}})$, for all
${\mathrm{v}}\in{\mathcal{V}}_{\mathbbm{k}}$, is strictly convex [49]. Note
that we have

${\mathcal{J}}({a},{\mathrm{h}})=f({a},{\mathrm{h}})-2{\mathrm{L}}({a},{\mathrm{h}})+\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}y_{t_{i}}^{2}+\delta_{{{\mathscr{R}}}}({a},{\mathrm{h}}),$
(37)

where ${\mathrm{L}}:{\mathcal{V}}_{\mathbbm{k}}\to{\mathbb{R}}$ is the bounded
linear operator defined as

${\mathrm{L}}({a},{\mathrm{h}})=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}y_{t_{i}}{\langle{({a},{\mathrm{h}}),\psi_{i}}\rangle}_{{\mathcal{V}}_{\mathbbm{k}}}.$
(38)

Therefore, since $f$ is strictly convex, ${\mathrm{L}}$ is linear, and,
$\delta_{{{\mathscr{R}}}}$ is convex, we know that
${\mathcal{J}}:{\mathcal{V}}_{\mathbbm{k}}\to{\mathbb{R}}$ is a strictly
convex function, and consequently, the the solution of optimization problem
(24) is unique [49]. This concludes the proof. ∎

Due to Theorem 10, we know that the convex program (24) has a unique solution
$({a}^{*},{\mathrm{h}}^{*})$. Meanwhile, one should note that (24) is an
infinite dimensional optimization problem with infinite number of constraints.
Thus, obtaining the solution $({a}^{*},{\mathrm{h}}^{*})$ is not
straightforward. On the other hand, since ${\mathrm{h}}^{*}$ belongs to the
set of stable responses ${\mathcal{H}}_{\mathbbm{k}}$ dominated by
${\mathrm{f}}_{\rho}$, one may intuitively expect that
${\mathrm{h}}^{*}\in{{\mathscr{R}}}_{m}$, for large enough
$m\in{\mathbb{Z}}_{+}$. In other words, the solution to the optimization
problem (24) is determined by a finite number of constraints and the remaining
constraints are unnecessary. In order to formalize this idea, let function
${\mathcal{J}}_{m}:{\mathcal{V}}_{\mathbbm{k}}\to\mathbb{R}\cup\\{+\infty\\}$
be defined as

${\mathcal{J}}_{m}({a},{\mathrm{h}})={\mathcal{E}}_{\rho}({a},{\mathrm{h}})+\sum_{s=0}^{m}\delta_{{{\mathscr{R}}}_{s}}({a},{\mathrm{h}})+\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2},$
(39)

and, consider the following program

$\inf_{({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}}\
\\!{\mathcal{J}}_{m}({a},{\mathrm{h}}).$ (40)

Note that (40) is equivalent to

$\\!\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{cl}\operatornamewithlimits{min}\limits_{{a}\in{\mathbb{R}},{\mathrm{h}}\in{\mathcal{H}}_{\mathbbm{k}}}\\!\\!\\!&\\!\\!\\!\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\\!\big{(}y_{t_{i}}-{a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})\big{)}^{2}\\!\\!+\\!\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2},\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\\
\text{s.t.}\\!\\!\\!&\\!\\!\\!h_{t}+{a}\rho^{t}\geq 0,\quad\forall\\!\
t\in\\{0,1,\ldots,m\\},\\\ \\!\\!\\!&\\!\\!\\!{a}\geq{a}_{\min}.\end{array}$
(41)

The next theorem guarantees the existence and uniqueness of solution for
optimization problem (40), or equivalently, for program (41).

######  Theorem 11.

Under the assumptions of Theorem 10, for each $m\in{\mathbb{Z}}_{+}$, problem
(40) admits unique solution
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$.

###### Proof.

Define set
${\mathscr{P}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\subset{\mathcal{V}}_{\mathbbm{k}}$
as
${\mathscr{P}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}=\bigcap_{s=0}^{m}{{\mathscr{R}}}_{s}$.
By replacing ${\mathscr{P}}$ with
${\mathscr{P}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}$ in the proof of
Theorem 10 and then repeating the same steps, the claim follows. ∎

Once the existence and uniqueness for the solution of (40) are established by
Theorem 11, a reasonable concern is the asymptotic behavior of
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$,
especially with respect to $({a}^{*},{\mathrm{h}}^{*})$. The next theorem
reveals this link saying that the solution
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$
coincides with $({a}^{*},{\mathrm{h}}^{*})$ when $m$ is large enough.

######  Theorem 12.

Let Assumptions 2 and 3 hold. Then, the following statements hold:
_i)_ Under Assumption 1, there is a non-negative integer $m$ such that
${a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}={a}^{*}$ and
${\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}={\mathrm{h}}^{*}$.
_ii)_ For a non-negative integer $m$, one has
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=({a}^{*},{\mathrm{h}}^{*})$
if and only if
${\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})<\infty$.
_iii)_ If
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=({a}^{*},{\mathrm{h}}^{*})$,
for a non-negative integer $m$, then
${a}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}}={a}^{*}$ and
${\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}}={\mathrm{h}}^{*}$,
for every $\overline{m}\geq m$.

###### Proof.

Part i) For any ${a}\geq{a}_{\min}$ and any $m\in{\mathbb{Z}}_{+}$, we have

${\mathcal{J}}_{m}({a},\mathbf{0})={\mathcal{J}}({a},\mathbf{0})=\sum_{i=1}^{n}\big{(}y_{t_{i}}-{a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})\big{)}^{2}\in[0,\infty).$
(42)

Let ${a}_{0}$ be
${\mathrm{argmin}}_{{a}\geq{a}_{\min}}{\mathcal{J}}({a},\mathbf{0})$, and,
define $C_{0}$ as $C_{0}:={\mathcal{J}}({a}_{0},\mathbf{0})$, which is equal
to ${\mathcal{J}}_{m}({a}_{0},\mathbf{0})$, for any $m\in{\mathbb{Z}}_{+}$. If
$C_{0}=0$, then one can easily see that $({a}_{0},\mathbf{0})$ is a solution
for optimization problems (24) and (40), where, because of their uniqueness,
$({a}_{0},\mathbf{0})=({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=({a}^{*},{\mathrm{h}}^{*})$.
Now, we consider the case $C_{0}>0$. Define $m_{0}$ as following

$m_{0}=\min\bigg{\\{}m\in{\mathbb{Z}}_{+}\\!\ \bigg{|}\\!\
m\geq\frac{1}{2}\frac{\ln(C_{0}C)-\ln({a}_{\min}^{2}\lambda)}{\ln(\rho)-\ln(\rho_{{\scalebox{0.5}{$\mathrm{d}$}}})}\bigg{\\}}.$
(43)

Based on the definition of $m_{0}$, one can easily see that

${a}_{\min}\rho^{s}-\lambda^{-\frac{1}{2}}C_{0}^{\frac{1}{2}}C^{\frac{1}{2}}\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{s}\geq
0,$ (44)

for any $s\geq m_{0}$. Let $m$ be an arbitrary integer such that $m\geq
m_{0}$, and, consider the convex program (40) with unique solution
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$.
We know that
$h_{s}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}+{a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\rho^{s}\geq
0$, for each $s=0,1,\ldots,m$. On the other hand, for $s>m$, from the
reproducing property and the Cauchy-Schwartz inequality, it follows that

$\begin{split}h_{s}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}+{a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\rho^{s}&\geq{\langle{{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},\mathbbm{k}_{s}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}+{a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\rho^{s}\\\
&\geq-\|{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\|_{{\mathcal{H}}_{\mathbbm{k}}}\mathbbm{k}(s,s)^{\frac{1}{2}}+{a}_{\min}\rho^{s}.\end{split}$
(45)

Note that due to
${\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})\leq{\mathcal{J}}_{m}({a}_{0},\mathbf{0})$,
one has

$\lambda\|{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2}\leq{\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})\leq
C_{0},$ (46)

which implies that
$\|{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\|_{{\mathcal{H}}_{\mathbbm{k}}}\leq\lambda^{-\frac{1}{2}}C_{0}^{\frac{1}{2}}$.
Hence, according to (44), (45) and Assumption 1, we have

$h_{s}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}+{a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\rho^{s}\geq{a}_{\min}\rho^{s}-\lambda^{-\frac{1}{2}}C_{0}^{\frac{1}{2}}C^{\frac{1}{2}}\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{s}\geq
0,$ (47)

which implies that
$h_{s}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}+{a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\rho^{s}\geq
0$, for all $s\in{\mathbb{Z}}_{+}$. Therefore, we have $\sum_{s\geq
1}\delta_{{{\mathscr{R}}}_{s}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=0$,
and subsequently, it follows that
${\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})={\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$.
On the other hand, according to the definition of $({a}^{*},{\mathrm{h}}^{*})$
and
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$,
we know that
${\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})\leq{\mathcal{J}}_{m}({a}^{*},{\mathrm{h}}^{*})$
and
${\mathcal{J}}({a}^{*},{\mathrm{h}}^{*})\leq{\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$.
Since
${\mathcal{J}}_{m}({a},{\mathrm{h}})\leq{\mathcal{J}}({a},{\mathrm{h}})$, for
all $({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}$, one can see that

${\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})\leq{\mathcal{J}}_{m}({a}^{*},{\mathrm{h}}^{*})\leq{\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}).$
(48)

Hence, we have
${\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})={\mathcal{J}}_{m}({a}^{*},{\mathrm{h}}^{*})$,
and, subsequently, due to Theorem 11, one has
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=({a}^{*},{\mathrm{h}}^{*})$.
This concludes the proof of part i).
Part ii) Consider the case that
${\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})<\infty$.
This implies that $\sum_{s\geq
1}\delta_{{{\mathscr{R}}}_{s}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=0$,
and consequently,
${\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})={\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$.
One can easily see that,

${\mathcal{J}}_{m}({a},{\mathrm{h}})\leq{\mathcal{J}}_{\overline{m}}({a},{\mathrm{h}})\leq{\mathcal{J}}({a},{\mathrm{h}}),$
(49)

for any $m\leq\overline{m}$ and each
$({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}$. Accordingly, due to the
definition of $({a}^{*},{\mathrm{h}}^{*})$ and
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$,
we have

$\begin{split}{\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})&\leq{\mathcal{J}}_{m}({a}^{*},{\mathrm{h}}^{*})\leq{\mathcal{J}}({a}^{*},{\mathrm{h}}^{*})\\\
&\ \
\qquad\leq{\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})={\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}),\end{split}$
(50)

which implies that $({a}^{*},{\mathrm{h}}^{*})$ is a solution for (40). Since,
this solution is unique according to Theorem 11, we need to have
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=({a}^{*},{\mathrm{h}}^{*})$.
The converse is straightforward, and concludes the proof of part ii).
Part iii) From the previous part, we know that
${\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})<\infty$.
Consequently, we have $\sum_{s\geq
1}\delta_{{{\mathscr{R}}}_{s}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=0$
which implies that

${\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})={\mathcal{J}}_{\overline{m}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})={\mathcal{J}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}).$
(51)

Accordingly, due to the definition of
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$
and
$({a}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}})$,
we have

$\begin{split}{\mathcal{J}}_{\overline{m}}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})&={\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})\\\
&\leq{\mathcal{J}}_{m}({a}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}})\leq{\mathcal{J}}_{\overline{m}}({a}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}}).\end{split}$
(52)

Therefore,
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})$
is the unique solution of optimization problem
$\inf_{({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}}\
\\!{\mathcal{J}}_{\overline{m}}({a},{\mathrm{h}})$, and consequently, we have

$({a}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}\overline{m}\mathrm{)}$}}}})=({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=({a}^{*},{\mathrm{h}}^{*}),$
(53)

where the second equality holds by assumption. This concludes the proof of
part iii) and the proof of Theorem 12. ∎

The following observation is a direct result of Theorem 12.

######  Corollary 13.

Under the assumptions of Theorem 12, there exists a non-negative integer
$m^{*}$ such that
$({a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=({a}^{*},{\mathrm{h}}^{*})$
if and only if $m\geq m^{*}$. Indeed, for any $m<m^{*}$, there exists
$s\in{\mathbb{Z}}_{+}$ such that
$h_{s}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}+{a}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}\rho^{s}<0$.
This implies that

$m^{*}=\min\big{\\{}m\in{\mathbb{Z}}_{+}\\!\ \big{|}\\!\
{\mathcal{J}}({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})<\infty\big{\\}}.$
(54)

Moreover, for $m_{0}$ introduced in (43), we have $m^{*}\leq m_{0}$.

######  Corollary 14.

Let $m\in{\mathbb{Z}}_{+}$ be such that the impulse response
${\mathrm{g}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}:=(h_{s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}+{a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\rho^{s})_{s=0}^{\infty}$
is not non-negative, i.e., $m<m^{*}$. Then, based on the proof of Theorem 12,
one can see that there exists $\underline{s}\in{\mathbb{Z}}_{+}$ such that
$g_{\underline{s}}<0$ and $s<m_{0}$, where $m_{0}$ is introduced in (43).

Due to Theorem 12 and Corollary 13, it suffices to consider only a finite
number of constraints in optimization problem (24) as in (40). The feasible
set of this optimization problem is infinite-dimension which makes the problem
intractable at the current format. In the remainder of this section, we derive
a practical heuristic for obtaining the solution of (12). To this end, we need
to introduce the _representer theorem_ and an additional definition.

######  Theorem 15 (Representer Theorem, [50, 51]).

Let $e:{\mathbb{R}}^{n}\to{\mathbb{R}}\cup\\{+\infty\\}$ be a given function,
${\mathcal{H}}$ be a Hilbert space endowed with inner product
${\langle{\cdot,\cdot}\rangle}_{{\mathcal{H}}}$, and
${\mathrm{w}}_{1},\ldots,{\mathrm{w}}_{n}$ be given vectors in
${\mathcal{H}}$. Consider the following optimization problem

$\min_{{\mathrm{w}}\in{\mathcal{H}}}\
e({\langle{{\mathrm{w}},{\mathrm{w}}_{1}}\rangle}_{{\mathcal{H}}},\ldots,{\langle{{\mathrm{w}},{\mathrm{w}}_{n}}\rangle}_{{\mathcal{H}}})+\kappa(\|{\mathrm{w}}\|_{{\mathcal{H}}}),$
(55)

where $\kappa:{\mathbb{R}}_{+}\to{\mathbb{R}}$ is an increasing function. If
(55) has a solution, there also exists a solution in the linear span of
${\mathrm{w}}_{1},\ldots,{\mathrm{w}}_{n}$.

Before proceeding to the next theorem, we define the matrices
${\mathrm{O}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}\times
n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}$,
${\mathrm{L}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}\times(m+1)}$
and ${\mathrm{K}}\in{\mathbb{R}}^{(m+1)\times(m+1)}$ as following

$\\!\\!\\!\\!\\!\\!\begin{array}[]{ll}{\mathrm{O}}(i,j)={\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i-1}}({\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{j-1}}(\mathbbm{k})),&1\\!\leq\\!i,j\\!\leq\\!n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}},\\\
{\mathrm{L}}(i,j)={\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i-1}}(\mathbbm{k}_{j-1}),&1\\!\leq\\!i\\!\leq\\!n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}},1\\!\leq\\!j\\!\leq\\!m\\!+\\!1,\\\
{\mathrm{K}}(i,j)=\mathbbm{k}(i-1,j-1),&1\\!\leq\\!i,j\\!\leq\\!m\\!+\\!1.\\\
\end{array}\\!\\!\\!\\!\\!\\!$ (56)

Also, the vectors
${\mathrm{y}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}$,
${\mathrm{b}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}$ and
${\mathrm{c}}\in{\mathbb{R}}^{m+1}$ are defined respectively as

${\mathrm{y}}:=\begin{bmatrix}y_{t_{i}}\end{bmatrix}_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\,,\quad{\mathrm{b}}:=\begin{bmatrix}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})\end{bmatrix}_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\,,$
(57)

and

${\mathrm{c}}:=\begin{bmatrix}\rho^{j}\end{bmatrix}_{j=0}^{m}\,.$ (58)

######  Theorem 16.

Let Assumption 3 hold. Then, for any non-negative integer $m$, there exists
${\mathrm{x}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}=[x_{0}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},\ldots,x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}]^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m+1}$ such
that the unique solution of (41),
${\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}=(h_{t}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})_{t=0}^{\infty}$,
admits the following parametric representation

$h^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}_{t}=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}(\mathbbm{k}_{t})+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\mathbbm{k}_{s}(t).$
(59)

Moreover,
$({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{x}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
is the solution of the following convex quadratic program

$\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{cl}\\!\\!\\!\operatornamewithlimits{min}\limits_{{a}\in{\mathbb{R}},\\!\
{\mathrm{x}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m+1}}\\!\\!\\!\\!&\\!\\!\big{\|}{\mathrm{y}}-{\mathrm{b}}{a}-\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\end{bmatrix}{\mathrm{x}}\big{\|}^{2}\\!+\\!\lambda{\mathrm{x}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\\!\\!\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\\\
{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$ }}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\
\end{bmatrix}{\mathrm{x}},\\!\\!\\\
\\!\\!\\!\mathrm{s.t.}\\!\\!\\!\\!&\\!\\!\begin{bmatrix}{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\
\end{bmatrix}{\mathrm{x}}+{\mathrm{c}}{a}\geq 0,\\\
\\!\\!\\!&\\!\\!{a}\geq{a}_{\min}.\end{array}\\!\\!\\!\\!\\!\\!\\!\\!$ (60)

###### Proof.

For $s=0,\ldots,m$, let ${\mathcal{A}}_{s}$ be the following set

${\mathcal{A}}_{s}=\big{\\{}({a},x)\in{\mathbb{R}}^{2}\\!\ \big{|}\\!\
x+{a}\rho^{s}\geq 0,\ {a}\geq{a}_{\min}\big{\\}},$ (61)

and, define the function
$e:{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m+1}\to{\mathbb{R}}\cup\\{+\infty\\}$
such that for any
${\mathrm{x}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m+1}$
we have

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!e&(x_{0},\ldots,x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m})=\\\
\\!&\min_{{a}\in{\mathbb{R}}}\Bigg{[}\\!\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}(y_{t_{i}}-{a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})-x_{i})\\!+\\!\sum_{s=0}^{m}\delta_{{\mathcal{A}}_{s}}\\!({a},x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s})\Bigg{]}.\\!\\!\\!\\!\\!\end{split}$
(62)

Also, for $i=0,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1$ and
$s=0,\ldots,m$, let $\varphi_{i}$ and
$\varphi_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}$ be defined respectively
as in (31) and
$\varphi_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}=\mathbbm{k}_{s}$. Due to
the reproducing property and (32), we know that
${\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})={\langle{{\mathrm{h}},\varphi_{i}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}$
and
$h_{s}={\langle{{\mathrm{h}},\mathbbm{k}_{s}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}$,
for any ${\mathrm{h}}=(h_{s})_{s=0}^{\infty}$. Accordingly, due to (62), one
can easily see that (40) is equivalent to the following optimization problem

$\min_{{\mathrm{h}}\in{\mathcal{H}}_{\mathbbm{k}}}e({\langle{{\mathrm{h}},\varphi_{0}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}},\ldots,{\langle{{\mathrm{h}},\varphi_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}})+\lambda\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2},$
(63)

with unique solution
${\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$. Therefore, due to
Theorem 15, we know that
${\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$ belongs to the
span of
$\varphi_{0},\varphi_{1},\ldots,\varphi_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m}$,
i.e., there exists
${\mathrm{x}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}=[x_{0}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},\ldots,x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}]^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m+1}$ such
that

${\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\varphi_{i}=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\varphi_{i}+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\mathbbm{k}_{s}.$
(64)

Due to reproducing property, we know that
$h_{t}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}={\langle{{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},\mathbbm{k}_{t}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}$,
for any $t\in{\mathbb{Z}}_{+}$. Accordingly, from (32), (64), the linearity
property of inner product, and the reproducing property, we have

$\begin{split}h^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}_{t}&=\big{\langle}\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\varphi_{i}+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\mathbbm{k}_{s},\\!\,\mathbbm{k}_{t}\big{\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}\\\
&=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\langle{\varphi_{i},\mathbbm{k}_{t}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\langle{\mathbbm{k}_{s},\mathbbm{k}_{t}}\rangle}_{{\mathcal{H}}_{\mathbbm{k}}}\\\
&=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}(\mathbbm{k}_{t})+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\mathbbm{k}_{s}(t).\end{split}$
(65)

Moreover, for $j=0,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1$, we have

$\begin{split}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{j}}({\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})&={\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{j}}(\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\varphi_{i}+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\mathbbm{k}_{s})\\\
&=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{j}}(\varphi_{i})+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{j}}(\mathbbm{k}_{s})\\\
&=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{j}}({\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}(\mathbbm{k}))+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{j}}(\mathbbm{k}_{s}).\end{split}$
(66)

Considering optimization problem (41), which is equivalent to (40), we replace
${\mathrm{h}}$ with the parametric form given in (64). Then, due to (56),
(65), (66), and the definition of vectors ${\mathrm{b}}$, ${\mathrm{c}}$ and
${\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$, the optimization
problem (60) follows. This concludes the proof. ∎

######  Remark 2.

Let the system be initially at rest and the sampling times be
${\mathscr{T}}=\\{0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1\\}$.
With respect to each $n_{1},n_{2}\in{\mathbb{Z}}_{+}$, we define matrix
${\mathrm{K}}_{n_{1},n_{2}}\in{\mathbb{R}}^{n_{1}\times n_{2}}$ such that
${\mathrm{K}}_{n_{1},n_{2}}(i,j)={\mathbbm{k}(i-1,j-1)}$, for
$i=1,\ldots,n_{1}$ and $j=1,\ldots,n_{2}$. Then, one can easily see that

${\mathrm{O}}={\mathrm{T}}_{{\mathrm{u}}}{\mathrm{K}}_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}},n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}{\mathrm{T}}_{{\mathrm{u}}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\,,\quad{\mathrm{L}}={\mathrm{T}}_{{\mathrm{u}}}{\mathrm{K}}_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}},m+1}\,,$
(67)

and

${\mathrm{K}}={\mathrm{K}}_{m+1,m+1}\,,$ (68)

where ${\mathrm{T}}_{{\mathrm{u}}}\in{\mathbb{R}}^{n\times n}$ is the Toeplitz
matrix defined as
${\mathrm{T}}_{{\mathrm{u}}}=[u_{i-j}]_{i,j=1}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}$.

Theorem 16 offers a practical way to solve problem (39). Due to Theorem 12 and
Corollary 13, we know that this solution coincides with the solution of (23)
provided that $m\geq m^{*}$. Nevertheless, compared to (23), the main
optimization problem (19) has an additional constraint on the rank of
resulting Hankel operator being finite. In the remainder of this section, we
fill this gap by employing the notion of finite Hankel rank kernels.

######  Definition 4.

We call the Mercer kernel $\mathbbm{k}:{\mathbb{Z}}_{+}\times{\mathbb{Z}}_{+}$
a _finite Hankel rank_ if, for each $s\in{\mathbb{Z}}_{+}$, the section of
kernel at $s$ is a finite Hankel rank impulse response, i.e.,

$\mathrm{rank}(\mathrm{Hankel}(\mathbbm{k}_{s}))<\infty,\quad\forall\\!\
s\in{\mathbb{Z}}_{+}.$ (69)

The commonly employed stable kernels in the literature [33, 52] are
_Tuned/Correlated_ (TC), _Diagonal/Correlated_ (TC), and _Stable Spline_ (TC),
which are respectively denoted by
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}$,
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{DC}$}}}}$ and
$\mathbbm{k}_{\SS}$, and defined as follows:

  * •

Tuned/Correlated (TC) kernel:

$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}(s,t)=\beta^{\max(s,t)},\quad\forall
s,t\in{\mathbb{Z}}_{+},$ (70)

where $\beta\in[0,1)$ is the hyperparameter of
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}$,

  * •

Diagonal/Correlated (DC) kernel:

$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{DC}$}}}}(s,t)=\beta^{\frac{s+t}{2}}\gamma^{|s-t|},\quad\forall
s,t\in{\mathbb{Z}}_{+},$ (71)

where $\beta\in[0,1)$ and $\gamma\in[-1,1]$ are the hyperparameters of
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{DC}$}}}}$, and,

  * •

Stable Spline (SS) kernel:

$\mathbbm{k}_{\SS}(s,t)=\frac{\beta^{s+t+\max(s,t)}}{2}-\frac{\beta^{3\max(s,t)}}{6},\quad\forall
s,t\in{\mathbb{Z}}_{+},$ (72)

where $\rho\in[0,1)$ is the hyperparameter of $\mathbbm{k}_{\SS}$.

As shown by the next theorem, these common kernels are finite Hankel rank.

######  Theorem 17.

Every finite support kernel,
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{DC}$}}}}$,
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}$ and
$\mathbbm{k}_{\SS}$ is finite Hankel rank.

###### Proof.

See Appendix A-C. ∎

Based on the notion of finite Hankel rank kernel and our previous discussion,
we can show when the solution of our estimation problem (12) can be obtained
by solving (60).

######  Theorem 18.

Under the assumptions of Theorem 12, if kernel $\mathbbm{k}$ is finite Hankel
rank, then the unique solution of (40) satisfies

$\mathrm{rank}\big{(}\mathrm{Hankel}({\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})\big{)}<\infty.$
(73)

Moreover,
$({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
is a solution of (12) provided that $m\geq m^{*}$.

###### Proof.

Due to (31) and (64), we know that

$\begin{split}{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}x_{i}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\sum_{s=0}^{t_{i}-\underline{t}}u_{t_{i}-s}\mathbbm{k}_{s}+\sum_{s=0}^{m}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\mathbbm{k}_{s}.\end{split}$
(74)

Rearranging the terms in (74), one can see that there exists real scalars
$\overline{x}_{0},\ldots,\overline{x}_{\overline{t}}$, where
$\overline{t}=\max\\{m,t_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}-\underline{t}\\}$,
such that we have

${\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}=\sum_{s=0}^{\overline{t}}\overline{x}_{s}\mathbbm{k}_{s}=\overline{x}_{0}\mathbbm{k}_{0}+\ldots+\overline{x}_{\overline{t}}\mathbbm{k}_{\overline{t}}.$
(75)

Therefore, we have

$\mathrm{Hankel}({\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}})=\sum_{s=0}^{\overline{t}}\overline{x}_{s}\mathrm{Hankel}(\mathbbm{k}_{s}),$
(76)

and subsequently, it follows that

$\mathrm{rank}(\mathrm{Hankel}({\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}}))\leq\sum_{s=0}^{\overline{t}}\mathrm{rank}(\mathrm{Hankel}(\mathbbm{k}_{s}))<\infty,$
(77)

i.e., (73) holds. For $m\geq m^{*}$, we know that
$({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
is the solution of (24). Accordingly, due to (73), we have

${\mathcal{J}}_{\\!{\mathscr{F}}}({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})={\mathcal{J}}({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})<\infty.$
(78)

On the other hand, for any $({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}$,
one can see that

${\mathcal{J}}({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})\leq{\mathcal{J}}({a},{\mathrm{h}})\leq{\mathcal{J}}_{\\!{\mathscr{F}}}({a},{\mathrm{h}}).$
(79)

Hence, due to (78), we have
${\mathcal{J}}_{\\!{\mathscr{F}}}({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})\leq{\mathcal{J}}_{\\!{\mathscr{F}}}({a},{\mathrm{h}})$,
for any $({a},{\mathrm{h}})\in{\mathcal{V}}_{\mathbbm{k}}$, i.e.,
$({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
is a solution of (12). This concludes the proof. ∎

Based on the above discussion, in order to solve the optimization problem
(12), it suffices to find the solution of the quadratic program (60) where
$m\geq m^{*}$ and $\mathbbm{k}$ is a given finite Hankel rank kernel such as
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}$,
$\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{DC}$}}}}$, or
$\mathbbm{k}_{\SS}$. Corollary 13 provides a bound for $m^{*}$. In some
special cases, we can provide a more practical bound.

######  Theorem 19.

Let the assumptions of Theorem 12 hold. If $\mathbbm{k}$ is either the TC
kernel (70) or the DC kernel (71), then we have $m^{*}\leq
t_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}-\underline{t}+1$.

###### Proof.

Let $m=t_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}-\underline{t}+1$, and
consider
$({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$,
the unique solution of (40). For $t=0,1,\ldots,m$, we know that
$h_{t}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}+{a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\rho^{t}\geq
0$. One the other hand, due to the definition of the TC kernel and (74), we
have
$h_{t}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}=\beta^{t-m}h_{m}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$,
for any $t>m$. Note that, due to Assumption 1, we have $\beta^{s}\leq
C\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{s}$, for all $s\in{\mathbb{Z}}_{+}$.
This implies that $\beta\leq\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{2}$.
Therefore, since $\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{2}<\rho$, we have

$\begin{split}h_{t}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}+{a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\rho^{t}&=\beta^{t-m}h_{m}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}+{a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\rho^{m}\rho^{t-m}\\\
&\geq\beta^{t-m}h_{m}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}+\beta^{t-m}{a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\rho^{m}\geq
0,\end{split}$ (80)

where the last equality is due to
$h_{m}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}+{a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\rho^{m}\geq
0$. Therefore,
$({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
is feasible for (24), and subsequently, we have
${\mathcal{J}}({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})<\infty$.
Hence, according to Theorem 12, we know that
$({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})=({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
which implies that $m^{*}\leq
m=t_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}-\underline{t}+1$. Based on a
similar argument, one can show the same result for the case of DC kernel. This
concludes the proof. ∎

######  Remark 3.

For the special cases $\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{TC}$}}}}$
and $\mathbbm{k}_{\text{{\scalebox{0.75}{$\mathrm{DC}$}}}}$, one can relax
Assumption 1 to
$|\mathbbm{k}(t,t)|\leq\rho_{{\scalebox{0.5}{$\mathrm{d}$}}}^{t}$, for all
$t\in{\mathbb{Z}}_{+}$.

Based on the above discussion, we present a practical implementation scheme
for identifying a linear internally positive dynamical system. The algorithm
and further details are provided in the next section.

##  VI Numerical Implementation Algorithm

Based on the discussion in Sections IV and V, to identify system
${\mathcal{S}}$ with internal positivity side-information, we need to solve
convex quadratic program (60). This optimization problem can be solved using
standard off-the-shelf solvers such as Matlab’s quadprog or cvx supported by
MOSEK [53]. Note that (60) depends on non-negative integer parameter $m$ which
is supposed be larger or equal to the parameter $m^{*}$ introduced in
Corollary 13. One possible approach is to set the value of $m$ to $m_{0}$,
which is introduced in (43) and guaranteed to have property. Also, one may
take initially $m$ equal to
$n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-\underline{t}$ and iteratively
increase it until $m$ exceeds $m^{*}$ and a non-negative impulse response is
obtained. This is of special interest when a suitable quadratic program solver
with a warm-starting feature is available [54]. In this iterative approach, at
each iteration $m$, one should check whether the estimated impulse response
${\mathrm{g}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}=(g^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}_{s})_{s=0}^{\infty}$
is non-negative, or equivalently, $m\geq m^{*}$. According to Corollary 14,
for verifying this stopping condition, it is enough to see whether
$g_{s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\geq 0$ holds, for
$s<m_{0}$. Due to (43), we know that $m_{0}$ depends logarithmically on the
parameters of problem, and thus, the size of $m_{0}$ is not prohibitive in the
practical examples. According to Theorem 19, when the TC kernel is employed
and the initial value of $m$ is set to
$n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-\underline{t}$, the introduced
iterative scheme takes only a single iteration. The outline of this approach
is summarized in Algorithm 1.

Algorithm 1 System Identification with Internal Positivity Side-Information

1:Input: Set of data ${\mathscr{D}}$, finite Hankel rank kernel stable
$\mathbbm{k}$, dominant pole $\rho$, regularization weight $\lambda$, and,
$\Delta_{m}\in{\mathbb{N}}$.

2:$m\ \reflectbox{$\to$}\
t_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}-\underline{t}+1$

3:while stopping/exiting condition is not met do

4: Calculate vectors ${\mathrm{y}}$, ${\mathrm{b}}$ and ${\mathrm{c}}$ based
on (57) and (58).

5: Obtain matrices ${\mathrm{O}}$, ${\mathrm{L}}$ and ${\mathrm{K}}$ as in
(56), or by Remark 2

6: Solve quadratic program (60) for
${a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$ and
${\mathrm{x}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$.

7: Obtain ${\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$ based on
(59), or equivalently, (74).

8: ${\mathrm{g}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\
\reflectbox{$\to$}\
({a}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}\rho^{s}+h_{s}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})_{s=0}^{\infty}$

9: if ${\mathrm{g}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$ is non-
negative, then exit the loop,

10: else$\quad m\ \reflectbox{$\to$}\ m+\Delta_{m}$

11: end

12:end

13:Input: Internally positive impulse response ${\mathrm{g}}^{*}$, and also
${\mathrm{x}}^{*}$, ${a}^{*}$, ${\mathrm{h}}^{*}$.

In order to initialize Algorithm 1, in addition to the set of data
${\mathscr{D}}$, we need a suitable kernel $\mathbbm{k}$ and also an
estimation of the dominant pole $\rho$ and the regularization weight
$\lambda$. In general, deciding on the type of kernel depends on the shape and
smoothness of the impulse response to be identified. Once the type of kernel
$\mathbbm{k}$ is set, we need to estimate the vector of hyperparameters
$\theta_{\mathbbm{k}}$ which characterizes $\mathbbm{k}$. Accordingly, vector
$\theta$ defined as $\theta:=[\rho,\lambda,\theta_{\mathbbm{k}}]\in\Theta$ is
the overall vector of hyperparameters to be determined where $\Theta$ denotes
the space of feasible hyperparameters. For estimating $\theta$, we employ a
_cross-validation_ scheme endowed with a _Bayesian optimization_ heuristic
[55]. In this regard, we split the index set of data
$\\{0,1,\ldots,n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1\\}$ into two
mutually disjoint subsets ${\mathcal{I}}_{\text{T}}$ and
${\mathcal{I}}_{\text{V}}$ respectively for _training_ and _validation_. We
define the model evaluation metric $v:\Theta\to{\mathbb{R}}$ based on the
prediction error of the validation data as

$v(\theta)=\frac{1}{|{\mathcal{I}}_{\text{V}}|}\sum_{i\in{\mathcal{I}}_{\text{V}}}\big{(}y_{t_{i}}-{a}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{f}}_{\rho})-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})\big{)}^{2},$
(81)

where $({a},{\mathrm{h}})$ is the result of proposed identification scheme
using the training data corresponding to indices ${\mathcal{I}}_{\text{T}}$,
and also to the hyperparameter $\theta$. We then estimate the hyperparameters
$\theta$ as $\hat{\theta}:={\mathrm{argmin}}_{\theta\in\Theta}\\!\ v(\theta)$.
Since the dependency of $v(\theta)$ on $\theta$ has a black-box oracle form,
we need to employ Bayesian optimization algorithms such as _GP-LCB_ or similar
alternatives [55]. These heuristics are readily available in Matlab’s bayesopt
function.

##  VII Further Internal Positivity Side-Information and Extensions

In Section IV, we employed positive system realization theory to formulate the
identification problem with internal positivity side-information. The
resulting optimization problem, (12), is formulated using the fact that the
transfer function of system ${\mathcal{S}}$ is in the following form

$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)=F^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)+H^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)=\frac{a}{1-\rho
z^{-1}}+H^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z).$ (82)

The transfer function $G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ has a dominant
part $F^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z):=a(1-\rho z^{-1})^{-1}$ with
$\rho\in(0,1)$, and a suppressed part
$H^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$. Given that the impulse response of
the system satisfies specific properties, this formulation can be further
extended to the following cases

$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)=F^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)+H^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)=\frac{N^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z^{-1})}{(1-\rho
z^{-1})^{n}}+H^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z),$ (83)

and

$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)=F^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)+H^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)=\frac{N^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z^{-1})}{1-\rho^{n}z^{-n}}+H^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z),$
(84)

where $N^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is a polynomial with degree
less than $n$, i.e.,
$N^{{\scalebox{0.5}{$({\mathcal{S}})$}}}\in{\mathbb{R}}_{n-1}[z^{-1}]$. Note
that according to [7, Theorem 9], the transfer function of a positive system
with non-zero spectral radius is in form of (83) or (84), which are
generalized forms of (82). These cases correspond to the situations where, in
addition to internal positivity, we have further information on the dominant
part of the impulse response of the system. In this section, we discuss these
extensions.

###  VII-A Non-simple Unique Dominant Pole

We first discuss the extension which corresponds to (83), i.e., the transfer
function has a unique dominant pole with multiplicity $n$ where $n$ can be
larger than one. The mathematical proofs are omitted since they are similar to
the ones given in Section V.

With respect to each $n\in{\mathbb{N}}$ and $\rho\in(0,1)$, let
${\mathscr{F}}_{\\!{{\alpha}},n}$ be the following set of impulse responses

$\begin{split}{\mathscr{F}}_{\\!{{\alpha}},n}&=\bigg{\\{}{\mathrm{f}}=(f_{t})_{t=0}^{\infty}\in\ell^{1}\\!\
\Big{|}\\!\
\operatornamewithlimits{lim}\limits_{t\to\infty}t^{-n+1}\rho^{-t}f_{t}>0,\\\
&\qquad\qquad(1-\rho z^{-1})^{n}\\!\
\sum_{t=0}^{\infty}f_{t}z^{-t}\in{\mathbb{R}}_{n-1}[z^{-1}]\bigg{\\}}.\end{split}$
(85)

Also, define the impulse response sets ${\mathscr{P}}_{\\!{{\alpha}},n}$ and
${\mathscr{P}}_{\\!(0,1),n}$ respectively as

$\begin{split}{\mathscr{P}}_{\\!{{\alpha}},n}&=\Big{\\{}{\mathrm{g}}={\mathrm{f}}+{\mathrm{h}}\in\underline{{\mathscr{P}}}\\!\
\Big{|}\\!\ {\mathrm{f}}\in{\mathscr{F}}_{\\!{{\alpha}},n},\\\ &\qquad\quad\
{\mathrm{h}}=(h_{t})_{t=0}^{\infty}\in\ell_{1},\lim_{t\to\infty}\rho^{-t}h_{t}=0\Big{\\}},\end{split}$
(86)

and
${\mathscr{P}}_{\\!(0,1),n}=\cup_{\rho\in(0,1)}{\mathscr{P}}_{\\!{{\alpha}},n}$.
According to [7, Theorem 11] and based on an argument similar to the proof of
Corollary 6, one can show that, for any $\rho\in(0,1)$, the impulse responses
in ${\mathscr{P}}_{\\!{{\alpha}},n}$ are internally positive. Indeed,
${\mathscr{P}}_{\\!{{\alpha}},n}$ is exactly the set of impulse responses of
positive systems with dominant pole structure as in (84). Thus, to identify
impulse response ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ with
internal positivity side-information in the sense that
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}\in{\mathscr{P}}_{\\!{{\alpha}},n}$
we need to estimate ${\mathrm{f}}$ and ${\mathrm{h}}$ with properties given in
(86). One can see that each
${\mathrm{f}}=(f_{t})_{t=0}^{\infty}\in{\mathscr{F}}_{\\!{{\alpha}},n}$ is
uniquely characterized in terms of real positive number $a$ and vector
${\mathrm{a}}=[a_{j}]_{j=0}^{n-2}\in{\mathbb{R}}^{n-2}$ as

$f_{t}=at^{n-1}\rho^{t}+\operatornamewithlimits{\sum}\limits_{j=0}^{n-2}a_{j}t^{j}\rho^{t},\qquad\forall
t\geq 0.$ (87)

Subsequently, we reintroduce the empirical loss function
${\mathcal{E}}_{\rho,n}:{\mathbb{R}}\times{\mathbb{R}}^{n-2}\times{\mathcal{H}}_{\mathbbm{k}}\to{\mathbb{R}}_{+}$
as following

$\begin{split}&{\mathcal{E}}_{\rho,n}(a,{\mathrm{a}},{\mathrm{h}}):=\\\
&\quad\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\Big{[}y_{t_{i}}-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}\Big{(}a{\mathrm{f}}_{\rho,n-1}+\sum_{j=0}^{n-2}a_{j}{\mathrm{f}}_{\rho,j}\Big{)}-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})\Big{]}^{2}\\!\\!,\end{split}$
(88)

where, for $j=0,\ldots,n-1$, the impulse response
${\mathrm{f}}_{{{\alpha}},j}$ is defined as
${\mathrm{f}}_{{{\alpha}},j}=(t^{j}\rho^{t})_{t=0}^{\infty}$. According to
(88), the identification problem (12) is updated to the following optimization
problem

$\begin{array}[]{cl}\operatornamewithlimits{min}\limits_{\begin{subarray}{c}a\in{\mathbb{R}},{\mathrm{h}}\in{\mathcal{H}}_{\mathbbm{k}}\\\
{\mathrm{a}}\in{\mathbb{R}}^{n-1}\end{subarray}}&{\mathcal{E}}_{\rho,n}(a,{\mathrm{a}},{\mathrm{h}})+\lambda\\!\
\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2}+\varepsilon\|{\mathrm{a}}\|^{2},\\\
\text{s.t.}&h_{t}+\rho^{t}\Big{[}at^{n-1}\\!+\\!\operatornamewithlimits{\sum}\limits_{j=0}^{n-2}a_{j}t^{j}\Big{]}\geq
0,\ \forall t\geq 0,\\\
&\mathrm{rank}(\mathrm{Hankel}({\mathrm{h}}))<\infty,\\\ &a\geq
a_{\min},\end{array}$ (89)

where $a_{\min}>0$ is a given lower-bound for $a$ and $\varepsilon>0$ is a
regularization weight. Based on the same line of argument as in Section V, the
next proposition presents a finite dimensional convex quadratic program
equivalent to (89). Before proceeding to the proposition, we define matrices
${\mathrm{B}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}\times
n}$ and ${\mathrm{C}}\in{\mathbb{R}}^{(m+1)\times n}$ as follows

$\\!\\!\\!\\!\\!\\!\begin{array}[]{ll}{\mathrm{B}}(i,j)={\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i-1}}({\mathrm{f}}_{\rho,j-1}),&1\\!\leq\\!i\\!\leq\\!n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}},1\\!\leq\\!j\\!\leq\\!n,\\\
{\mathrm{C}}(i,j)=(i-1)^{j-1}\rho^{i-1},&1\\!\leq\\!i\\!\leq\\!m\\!+\\!1,1\\!\leq\\!j\\!\leq\\!n.\end{array}\\!\\!\\!\\!\\!\\!$
(90)

######  Theorem 20.

Let Assumptions 1, 2 and 3 hold and $\mathbbm{k}$ be a finite Hankel rank
kernel. Also, with respect to each $m\in{\mathbb{Z}}_{+}$, let
$(a^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{a}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{x}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
be the solution of following convex quadratic program

$\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{cl}\\!\\!\\!\operatornamewithlimits{min}\limits_{\begin{subarray}{c}{\mathrm{x}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m+1}\\\
a\in{\mathbb{R}},{\mathrm{a}}\in{\mathbb{R}}^{n-1}\end{subarray}}\\!\\!\\!\\!&\\!\Big{\|}{\mathrm{y}}\\!-\\!{\mathrm{B}}\\!\begin{bmatrix}{\mathrm{a}}\\\
a\end{bmatrix}\\!\\!-\\!\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\end{bmatrix}\\!{\mathrm{x}}\Big{\|}^{2}\\!+\\!\lambda{\mathrm{x}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\\!\\!\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\\\
{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$ }}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\
\end{bmatrix}{\mathrm{x}}\\!+\\!\varepsilon\|{\mathrm{a}}\|^{2}\\!\\!,\\\
\\!\\!\\!\mathrm{s.t.}\\!\\!\\!\\!&\\!\begin{bmatrix}{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\
\end{bmatrix}{\mathrm{x}}+{\mathrm{C}}\\!\begin{bmatrix}{\mathrm{a}}\\\
a\end{bmatrix}\geq 0,\\\ \\!\\!\\!&\\!\\!a\geq
a_{\min},\end{array}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$ (91)

and, the impulse response
${\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$ is defined
according to (59). Then, there exists $m^{*}$ such that
$(a^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{a}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}},{\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
is a solution of (89), for any $m\geq m^{*}$. Moreover, for any
$m_{1},m_{2}\geq m^{*}$, we have
$(a^{{{\scalebox{0.63}{$\mathrm{(}m_{1}\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}m_{1}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m_{1}\mathrm{)}$}}}})=(a^{{{\scalebox{0.63}{$\mathrm{(}m_{2}\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}m_{2}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m_{2}\mathrm{)}$}}}})$.

###  VII-B Multiple Simple Dominant Poles

In this section, we introduce the extension corresponding to the case in (84),
i.e., the dominant part of the transfer function has specially structured
multiple simple dominant poles.

For any $n\in{\mathbb{N}}$ and any $\rho\in(0,1)$, defined the impulse
response sets
${\mathscr{F}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}$,
${\mathscr{P}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}$
and ${\mathscr{P}}_{\\!(0,1)}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}$
respectively as

$\begin{split}{\mathscr{F}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}&=\bigg{\\{}{\mathrm{f}}=(f_{t})_{t=0}^{\infty}\in\underline{{\mathscr{P}}}\\!\
\Big{|}\\!\
\operatornamewithlimits{liminf}\limits_{t\to\infty}\rho^{-t}f_{t}>0,\\\
&\qquad\qquad(1-\rho^{n}z^{-n})\sum_{t=0}^{\infty}f_{t}z^{-t}\in{\mathbb{R}}_{n-1}[z^{-1}]\bigg{\\}},\\\
{\mathscr{P}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}&=\bigg{\\{}{\mathrm{g}}={\mathrm{f}}+{\mathrm{h}}\in\underline{{\mathscr{P}}}\\!\
\Big{|}\\!\
{\mathrm{f}}\in{\mathscr{F}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}},\\\
&\qquad\qquad{\mathrm{h}}=(h_{t})_{t=0}^{\infty}\in\ell_{1},\operatornamewithlimits{liminf}\limits_{t\to\infty}\rho^{-t}|h_{t}|=0\bigg{\\}},\end{split}$
(92)

and
${\mathscr{P}}_{\\!(0,1)}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}=\cup_{\rho\in(0,1)}{\mathscr{P}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}$.
One can easily see that for $n=1$ these sets coincide with
${\mathscr{P}}_{\\!{{\alpha}}}$ and ${\mathscr{P}}_{\\!(0,1)}$. Due to [7,
Theorem 12] and by following the same lines of argument as in the proof of
Corollary 6, we can show that
${\mathscr{P}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}$
contains internally positive impulse responses, for any $\rho\in(0,1)$.
Accordingly, the identification with internally positivity side-information in
the sense that impulse response belongs to
${\mathscr{P}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}$
translates to the estimation of ${\mathrm{f}}$ and ${\mathrm{h}}$ with
properties given in (92). Note that with respect to each
${\mathrm{f}}=(f_{t})_{t=0}^{\infty}\in{\mathscr{F}}_{\\!{{\alpha}}}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}$,
there exist real scalars
$a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{0},\ldots,a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{n-1}$
and
$a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{0},\ldots,a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{n-1}$
such that, for any $t\in{\mathbb{Z}}_{+}$, we have

$\begin{split}f_{t}&=\Big{(}\sum_{k=0}^{n-1}(a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{k}+\mathrm{j}a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{k})\rho^{t}\omega^{kt}\Big{)},\\\
0&=\mathrm{imag}\Big{(}\sum_{k=0}^{n-1}(a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{k}+\mathrm{j}a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{k})\rho^{t}\omega^{kt}\Big{)},\end{split}$
(93)

where $\omega=\mathrm{e}^{\frac{2\pi}{n}}$, i.e., ${\mathrm{f}}$ is uniquely
characterized in terms of vectors
${\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}=[a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{k}]_{k=0}^{n-1}$
and
${\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}=[a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{k}]_{k=0}^{n-1}$.
Hence, we can reintroduce the empirical loss function
${\mathcal{E}}_{\rho}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\times{\mathcal{H}}_{\mathbbm{k}}\to{\mathbb{R}}_{+}$
as follows

$\begin{split}{\mathcal{E}}_{\rho}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}&({\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}},{\mathrm{h}}):=\\\
&\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\Big{[}y_{t_{i}}-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}\Big{(}\sum_{k=0}^{n-1}(a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}-a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}})\Big{)}-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{h}})\Big{]}^{2}\\!\\!,\end{split}$

where
${\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}$
and
${\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}$
are defined as

${\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}=\big{(}\rho^{t}\\!\
(\omega^{kt})\big{)}_{t=0}^{\infty},\quad{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}=\big{(}\rho^{t}\\!\
\mathrm{imag}(\omega^{kt})\big{)}_{t=0}^{\infty},$ (94)

for $k=0,\ldots,n-1$. Therefore, the identification problem (12) is modified
to

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{cl}\operatornamewithlimits{min}\limits_{\begin{subarray}{c}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\in{\mathbb{R}}^{n}\\\
{\mathrm{h}}\in{\mathcal{H}}_{\mathbbm{k}}\end{subarray}}&{\mathcal{E}}_{\rho}^{{{\scalebox{0.63}{$\mathrm{(}n\mathrm{)}$}}}}({\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}},{\mathrm{h}})\\!+\\!\lambda\\!\
\|{\mathrm{h}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2}\\!+\\!\varepsilon\|{\mathrm{E}}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}\|^{2}\\!+\\!\varepsilon\|{\mathrm{E}}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\|^{2}\\!\\!,\\\
\text{s.t.}&h_{t}+\operatornamewithlimits{\sum}\limits_{k=0}^{n-1}(a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}-a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}})\geq
0,\quad\forall t\geq 0,\\\
&\operatornamewithlimits{\sum}\limits_{k=0}^{n-1}(a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}+a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}})=0,\quad\forall
t\geq 0,\\\
&\operatornamewithlimits{liminf}\limits_{t\to\infty}\rho^{-t}\operatornamewithlimits{\sum}\limits_{k=0}^{n-1}(a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}-a^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{k}{\mathrm{f}}_{{{\alpha}},k}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}})\geq
a_{\min},\\\ &\mathrm{rank}(\mathrm{Hankel}({\mathrm{h}}))<\infty,\\\
\end{array}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(95)

where ${a}_{\min}>0$ is a small positive real scalar, $\varepsilon>0$ is a
regularization weight and ${\mathrm{E}}\in{\mathbb{R}}^{n\times n}$ is defined
as ${\mathrm{E}}=\mathrm{diag}(0,1,1,\ldots,1)$. By an argument similar to
Section V, we can derive an equivalent finite dimensional convex quadratic
program for (95). To this end, define matrices
${\mathrm{V}}_{m}\in{\mathbb{R}}^{m\times n}$,
${\mathrm{B}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}\times
n}$, and,
${\mathrm{B}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}\times
n}$, respectively as

$\begin{array}[]{ll}{\mathrm{V}}_{m}(i,j)=\omega^{(i-1)(j-1)},&1\\!\leq\\!i\\!\leq\\!m,1\\!\leq\\!j\\!\leq\\!n,\\\
{\mathrm{B}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}(i,j)={\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i-1}}\\!({\mathrm{f}}_{\rho,j-1}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}),&1\\!\leq\\!i\\!\leq\\!n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}},1\\!\leq\\!j\\!\leq\\!n,\\\
{\mathrm{B}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}(i,j)={\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i-1}}\\!({\mathrm{f}}_{\rho,j-1}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}),&1\\!\leq\\!i\\!\leq\\!n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}},1\\!\leq\\!j\\!\leq\\!n,\\\
\end{array}$ (96)

for $m\in{\mathbb{Z}}_{+}\cup\\{\infty\\}$. Also, let
${\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{m}$,
${\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{m}$ and
${\mathrm{D}}_{m}$ be defined respectively as $({\mathrm{V}}_{m})$,
$\mathrm{imag}({\mathrm{V}}_{m})$, and
$\mathrm{diag}(1,\rho,\ldots,\rho^{m-1})$. One can see that the first
constraint in (95) is equivalent to

${\mathrm{h}}+{\mathrm{D}}_{\infty}{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{\infty}\\!\
{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}-{\mathrm{D}}_{\infty}{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{\infty}\\!\
{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\geq 0.$
(97)

The second constraint in (95) is

${\mathrm{D}}_{\infty}{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{\infty}\\!\
{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}+{\mathrm{D}}_{\infty}{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{\infty}\\!\
{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}=0,$ (98)

which implies that
${\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{\infty}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}+{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{\infty}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}=0$
due to $\rho>0$. As $\omega^{n}=1$, we know that ${\mathrm{V}}_{m}$ is an
$n$-periodic Vandermonde matrix. Therefore, (98) is equivalent to
${\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}+{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}=0$.
Similarly, one can show that the third constraint in (95) is equivalent to
${\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}+{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\geq
a_{\min}1_{n}$. Based on the discussion above and similar to those in Section
V, the next proposition presents the finite-dimensional convex quadratic
program equivalent to (95).

######  Theorem 21.

Let the assumptions of Proposition 20 hold. For any $m\in{\mathbb{Z}}_{+}$,
let
$({\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r},m\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i},m\mathrm{)}$}}}},{\mathrm{x}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
be the solution of following convex quadratic program

$\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{cl}\\!\\!\\!\operatornamewithlimits{min}\limits_{\begin{subarray}{c}{\mathrm{x}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+m+1}\\\
{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}\\!,\\!\
{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\in{\mathbb{R}}^{n}\end{subarray}}\\!\\!\\!\\!&\\!\Big{\|}{\mathrm{y}}\\!-\\!{\mathrm{B}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}\\!-\\!{\mathrm{B}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\\!-\\!\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\end{bmatrix}\\!{\mathrm{x}}\Big{\|}^{2}\\!\\!\\!+\\!\lambda{\mathrm{x}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\\!\\!\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\\\
{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$ }}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\
\end{bmatrix}\\!{\mathrm{x}}\\\ &\qquad+\\!\
\varepsilon\|{\mathrm{E}}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}\|^{2}+\varepsilon\|{\mathrm{E}}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\|^{2}\\!\\!,\\\
\\!\\!\\!\mathrm{s.t.}\\!\\!\\!\\!&\\!\begin{bmatrix}{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\
\end{bmatrix}\\!{\mathrm{x}}+{\mathrm{D}}_{m+1}\\!\begin{bmatrix}{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{m+1}\\!\\!&\\!\\!-{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{m+1}\\!\end{bmatrix}\\!\\!\begin{bmatrix}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}\\\
{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\end{bmatrix}\\!\\!\geq\\!0,\\\
\\!\\!\\!&\\!\\!{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}+{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}=0,\\\
\\!\\!\\!&\\!\\!{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r}\mathrm{)}$}}}}+{\mathrm{V}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}_{n}{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i}\mathrm{)}$}}}}\geq
a_{\min}1_{n},\end{array}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$ (99)

and, define the impulse response
${\mathrm{h}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}}$ by (59). Then,
there exists $m^{*}$ such that
$({\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{r},m\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}\mathrm{i},m\mathrm{)}$}}}},{\mathrm{x}}^{{\scalebox{0.63}{$\mathrm{(}m\mathrm{)}$}}})$
is a solution of (89), for each $m\geq m^{*}$. Moreover, for any
$m_{1},m_{2}\geq m^{*}$, we have
$({\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}r,m_{1}\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}i,m_{1}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m_{1}\mathrm{)}$}}}})=({\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}r,m_{2}\mathrm{)}$}}}},{\mathrm{a}}^{{{\scalebox{0.63}{$\mathrm{(}i,m_{2}\mathrm{)}$}}}},{\mathrm{h}}^{{{\scalebox{0.63}{$\mathrm{(}m_{2}\mathrm{)}$}}}})$.

###  VII-C Zero Spectral Radius

According to [7, Theorem 8], when the transfer function
$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ has zero spectral radius,
$r(G^{{\scalebox{0.5}{$({\mathcal{S}})$}}})=0$, the impulse response
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is internally positive if
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is non-negative. Based on
the proof of [7, Theorem 8], we know that
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ belongs to $c_{00}$, the
space of finitely non-zero impulse responses, i.e., there exists
$n_{{\mathrm{g}}}\in{\mathbb{Z}}_{+}$ such that
$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)=\sum_{t=0}^{n_{{\mathrm{g}}}-1}g_{t}z^{-t}$
and $g_{t}=0$, for all $t\geq n_{{\mathrm{g}}}$. Comparing to the other cases
of internal positivity side-information, the case of zero spectral radius
provides the weakest information. Indeed, this knowledge only says the
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is a non-negative and
finitely non-zero sequence and provides no further information about the
behavior of the impulse response of system ${\mathcal{S}}$.

Based on the discussion above, the identification problem with the internal
positivity side-information and the extra knowledge
$r(G^{{\scalebox{0.5}{$({\mathcal{S}})$}}})=0$ can be easily formulated as
following

$\begin{array}[]{cl}\operatornamewithlimits{min}\limits_{{\mathrm{g}}\in{\mathcal{H}}_{\mathbbm{k}}}&\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}-1}\big{(}y_{t_{i}}-{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}({\mathrm{g}})\big{)}^{2}+\lambda\\!\
\|{\mathrm{g}}\|_{{\mathcal{H}}_{\mathbbm{k}}}^{2},\\\ \text{s.t.}&g_{t}\geq
0,\quad\forall t=0,1,\ldots,n_{{\mathrm{g}}}-1,\\\ \end{array}$ (100)

where $\mathbbm{k}$ is a kernel zero on
${\mathbb{Z}}_{+}^{2}\backslash\\{0,\ldots,n_{{\mathrm{g}}}-1\\}^{2}$. One can
show that for ${\mathrm{g}}^{*}=(g_{t}^{*})_{t=0}^{\infty}$, the solution of
(100), and for $t=0,1,\ldots,n_{{\mathrm{g}}}$, we have

$g_{t}=\sum_{i=0}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}}x_{i}{\mathrm{L}}^{\\!{\mathrm{u}}}_{t_{i}}(\mathbbm{k}_{t})+\sum_{s=0}^{n_{{\mathrm{g}}}-1}x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+s}\mathbbm{k}_{s}(t),$
(101)

where
${\mathrm{x}}=[x_{0},\ldots,x_{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+n_{{\mathrm{g}}}-1}]^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!$ is the solution of following _convex quadratic program_

$\begin{array}[]{cl}\operatornamewithlimits{min}\limits_{{\mathrm{x}}\in{\mathbb{R}}^{n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}+n_{{\mathrm{g}}}}}&\big{\|}{\mathrm{y}}-\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\end{bmatrix}{\mathrm{x}}\big{\|}^{2}+\lambda{\mathrm{x}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\begin{bmatrix}{\mathrm{O}}\\!\\!&\\!\\!{\mathrm{L}}\\\
{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$ }}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\
\end{bmatrix}{\mathrm{x}},\\\
\mathrm{s.t.}&\begin{bmatrix}{\mathrm{L}}^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\\!\\!&\\!\\!{\mathrm{K}}\\\ \end{bmatrix}{\mathrm{x}}\geq
0.\end{array}$ (102)

Also, one can solve (100) directly for
$[g_{0},\ldots,g_{n_{{\mathrm{g}}}-1}]^{{\scalebox{0.63}{$\mathsf{T}$
}}}\\!\in{\mathbb{R}}^{n_{{\mathrm{g}}}}$ at the cost of a matrix inversion
which can be both imprecise and computationally demanding, especially when
$n_{{\mathrm{g}}}$ is large.

##  VIII Numerical Experiments

In this section, we provide numerical and experimental examples to verify the
efficacy and performance of the proposed method for impulse response
identification with internal positivity side-information. The first example
concerns the impact of incorporating internal positivity side-information on
the estimation quality and providing a comparative analysis for the proposed
identification scheme through a Monte Carlo analysis. The second example
concerns the efficacy of the proposed identification scheme on a set of data
collected from an experimental heating system.

###  VIII-A Monte Carlo Experiment

Consider system ${\mathcal{S}}$ described with the impulse response
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}=({g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{t})_{t=0}^{\infty}$
defined as

${g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{t}={{\rho}}^{t}(1+\beta^{t}\cos(2\pi\omega
t)),\qquad\forall\\!\ t\in{\mathbb{Z}}_{+},$ (103)

where ${{\rho}}$ and $\beta$ are real scalars in $(0,1)$, and $\omega$ is an
irrational real number in $(0,1)$. One can easily see that
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is a non-negative impulse
response with transfer function

$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)\\!=\\!\frac{1}{1\\!-\\!{{\rho}}z^{-1}}+\frac{1-{{\rho}}\beta\cos
w\ z^{-1}}{1-2{{\rho}}\beta\cos w\ z^{-1}+{{\rho}}^{2}\beta^{2}z^{-2}}.$ (104)

Therefore, according to Corollary 6, we know that
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is internally positive.

####  VIII-A1 Simulation Configuration

In this numerical experiment, we set ${{\rho}}=0.98$, $\beta=0.92$, and
$\omega=\frac{1}{10}\pi^{2}$. Using Matlab’s idinput function, we generate a
set of $120$ random binary input signals, each with length of
$n_{{\scalebox{0.5}{$\\!{\mathscr{D}}$}}}=200$. The system is initially at
rest. The input signals are applied to the system and the corresponding
noiseless output is obtained. We consider three signal-to-noise ratio (SNR)
levels of $10$ dB, $20$ dB, and $30$ dB. With respect to each of these SNR
levels and each output signal, we generate a zero-mean white Gaussian signal
as the additive measurement uncertainty. The resulting noisy output is
measured at time instants $t_{i}=i$, for $i=0,1,\ldots,199$. Accordingly, with
respect to each of the mentioned SNR levels, we have $120$ sets of input-
output data as follows4

$\begin{split}{\mathscr{D}}_{i}^{{{\scalebox{0.63}{$\mathrm{(}10\text{dB}\mathrm{)}$}}}}&\\!=\\!\big{\\{}(u_{s}^{{{\scalebox{0.63}{$\mathrm{(}i\mathrm{)}$}}}},y_{s}^{{{\scalebox{0.63}{$\mathrm{(}i,10\text{dB}\mathrm{)}$}}}})\big{|}s\\!=\\!0,\ldots,199\big{\\}},\
\ i\\!=\\!1,\ldots,120,\\\
{\mathscr{D}}_{i}^{{{\scalebox{0.63}{$\mathrm{(}20\text{dB}\mathrm{)}$}}}}&\\!=\\!\big{\\{}(u_{s}^{{{\scalebox{0.63}{$\mathrm{(}i\mathrm{)}$}}}},y_{s}^{{{\scalebox{0.63}{$\mathrm{(}i,20\text{dB}\mathrm{)}$}}}})\big{|}s\\!=\\!0,\ldots,199\big{\\}},\
\ i\\!=\\!1,\ldots,120,\\\
{\mathscr{D}}_{i}^{{{\scalebox{0.63}{$\mathrm{(}30\text{dB}\mathrm{)}$}}}}&\\!=\\!\big{\\{}(u_{s}^{{{\scalebox{0.63}{$\mathrm{(}i\mathrm{)}$}}}},y_{s}^{{{\scalebox{0.63}{$\mathrm{(}i,30\text{dB}\mathrm{)}$}}}})\big{|}s\\!=\\!0,\ldots,199\big{\\}},\
\ i\\!=\\!1,\ldots,120,\end{split}$

where the superscript refers to the SNR level in the corresponding set of
data.

####  VIII-A2 Comparison Methods

For estimating the impulse response of system, we utilize the input-output
data sets and the following identification methods:

  * A.

The first method is based on the subspace approach implemented by Matlab’s
n4sid and using the true order of system ${\mathcal{S}}$. Once we obtain an
initial estimation
$\tilde{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}1\mathrm{)}$}}}}$, the
result is projected on the positive orthant by
$\hat{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}1\mathrm{)}$}}}}=\max(\tilde{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}1\mathrm{)}$}}}},0)$,
where the max-operation is performed coordinate-wise.

  * B.

In the second method, we employ the least squares approach, and then, similar
method A, the projection on the positive orthant is applied to estimate a non-
negative finite impulse response. More precisely, the estimation
$\hat{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}2\mathrm{)}$}}}}$ is
obtained by
$\hat{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}2\mathrm{)}$}}}}=\max(\tilde{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}2\mathrm{)}$}}}},0)$,
where

$\tilde{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}2\mathrm{)}$}}}}:={\mathrm{argmin}}_{{\mathrm{g}}\in{\mathbb{R}}^{n_{{\mathrm{g}}}}}\|{\mathrm{T}}_{{\mathrm{u}}}{\mathrm{g}}-{\mathrm{y}}\|^{2},$
(105)

and, vector ${\mathrm{y}}$ and Toeplitz matrix ${\mathrm{T}}_{{\mathrm{u}}}$
are respectively defined in (57) and Remark 2.

  * C.

This method is based on a constrained least squares approach, where the
external positivity feature is enforced by setting the feasible set to the
positive orthant. In other words, the impulse response is estimated as

$\hat{{\mathrm{g}}}^{{{\scalebox{0.63}{$\mathrm{(}3\mathrm{)}$}}}}:={\mathrm{argmin}}_{{\mathrm{g}}\in{\mathbb{R}}_{+}^{n_{{\mathrm{g}}}}}\|{\mathrm{T}}_{{\mathrm{u}}}{\mathrm{g}}-{\mathrm{y}}\|^{2}.$
(106)

  * D.

& E. The fourth and fifth methods are respectively similar to the second and
third approaches, but with an additional kernel-based regularization term
included in the corresponding optimization problems. For method D, one can
employ Matlab’s impulseest function and then taking the positive part of the
resulting FIR. Also, method E essentially corresponds to (100), or
equivalently (102).

  * F.

The sixth method is a Bayesian FIR estimation scheme for externally positive
systems [28, 29]. This scheme is based on maximum a posteriori estimation,
where the employed prior is a maximum entropy distribution with support on
positive orthant and kernel-based covariance.

  * G.

The last method is the scheme proposed in this paper and summarized in
Algorithm 1 (see Section VI).

One should note that in all of the mentioned methods, the resulting impulse
responses are non-negative. In order to have a fair comparison, in the kernel-
based methods D, E, F, and G, the same kernel type (71) is employed.

####  VIII-A3 Evaluation Metrics and Results

For evaluating the performances of these methods, we compare the resulting
bias-variance trade-offs, as shown in Table I.

Method A B C D E F G 10 dB Bias$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 2.14 2.86
1.60 0.80 0.77 0.62 0.27 Var$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 6.54 34.7 11.8
1.77 1.74 1.74 0.51 MSE$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 11.1 42.9 14.8 2.40
2.34 2.06 0.58 20 dB Bias$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 2.01 0.78 0.48
0.47 0.38 0.37 0.089 Var$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 1.26 3.90 1.69 0.99
0.85 0.84 0.114 MSE$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 5.32 4.50 1.93 1.22 0.99
0.97 0.122 30 dB Bias$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 2.01 0.27 0.23 0.22
0.24 0.25 0.031 Var$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 0.17 0.58 0.32 0.33 0.55
0.54 0.015 MSE$(\hat{{\mathrm{g}}})\\!\\!\\!\\!$ 4.20 0.65 0.37 0.38 0.61 0.60
0.016

Table I: The bias, variance and MSE resulting from the identification methods
listed in Section VIII-A. The last column corresponds to the proposed approach
which integrates internal positivity.

Moreover, for further quantitative comparison of the estimated impulse
responses, we use _coefficient of determination_ , or _R-squared_ , which is
denoted by $\mathrm{fit}$ and defined as

$\mathrm{fit}(\hat{{\mathrm{g}}})=100\times\left(1-\frac{\|\hat{{\mathrm{g}}}-{\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}\|_{2}}{\|{\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}\|_{2}}\right),$
(107)

where $\hat{{\mathrm{g}}}$ is the estimated impulse response. Figure 1
compares the resulting quality of fit for the different SNR levels.

Figure 1: Box-plots of the R-squared metric for the estimation results of the
different methods and SNR levels. All fits for method A are below $85$%.

####  VIII-A4 Discussion

We have several observations from the Monte Carlo numerical experiment. As
shown in Figure 1, all methods are outperformed by the proposed identification
scheme, which maximally incorporates the internal positivity side-information.
Indeed, the side-information helps excluding spurious model candidates and
subsequently increases the accuracy of the estimation. The bias-variance
results presented in Table I confirm this fact. While each of these methods
estimates a non-negative impulse response and partially integrates positivity
side-information, we can see that the level of integrated side-information is
less than that of the proposed method. The proposed approach incorporates this
information in the model maximally. Comparing methods B and C, one can see the
former one is a two-step procedure where the estimation is performed in the
first step and non-negativity of the impulse response is obtained in the
second step, while the latter approach is a single-step procedure that
considers impulse response non-negativity during the estimation. On the other
hand, according to the fitting results shown in Figure 1, C performs better
than method B. For methods D and E, we have similar arguments. This
observation highlights the importance of jointly considering the positivity
with the impulse response estimation, as done by the proposed method. Method A
knows the actual order of the system. However, according to the results
presented in Figure 1 and Table I, one can see that positivity is a more
advantageous and stronger side-information for impulse response estimation,
especially when it is incorporated with its maximum strength like in the
proposed scheme. Finally, one can see that the kernel-based methods D, E, F,
and G have better estimation performance comparing to the methods A, B, and C,
which is expected [33].

Figure 2: The experimental system (left) and the corresponding block-diagram
schematic (right).

###  VIII-B Heating System Experiment

In this example, we verify the efficacy of the proposed identification scheme
on a set of data collected from an experimental nonlinear heating system [56].

Figure 2 shows the experiment configuration and the corresponding control and
measurement schematic. In this experiment, a metal plate is heated up by a
$300$-watt Halogen lamp mounted almost $5$cm above the center of the plate. On
the other side of the plate, a thermocouple is placed, measuring the
temperature. The thermocouple is connected to a computer via an A/D board for
sampling and recording the temperature measurements. The lamp is supplied by a
thyristor-based power amplifier driven by a D/A board and controlled by a
computer. The sampling time for control and data acquisition is
$T_{\text{s}}=2$s. Accordingly, we have a nonlinear discrete-time system from
the input of D/A board to the output of A/D board, which is intuitively close
to a linear positive system. The nonlinearity of the system is mainly due to
the power amplifier [56]. Moreover, the system is subject to delay and
disturbances from the ambient. To make the identification problem more
challenging, we disregard the nonlinearity and external disturbances issues.

The system is actuated by a piece-wise constant input and the output of system
is measured for $801$ samples. The collected input-output data is shown in
Figure 3, which is also available in the DAISY database [57].

Figure 3: The figure demonstrates the measurement data corresponding to the
experiments.

We employ this data for identifying the system using the identification
methods discussed in Section VIII-A, and then, compare the results. Since the
lamp has failed near to the end of experiment and the tail of data has less
fidelity, we discard the last $101$ samples ($200$ seconds). We split the data
into a _training set_ , to be used for identification, and, a _test set_ ,
which is the base for comparing the identified models. The samples training
set contains the first $500$ measurement samples, and the next $200$ data
points belong to the test set. The quality of estimated model is evaluated
based on the _R-squared_ metric, which assess the prediction precision on the
test data, and, defined as follows

$\mathrm{fit}(\hat{{\mathrm{g}}})=100\times\left(1-\Bigg{[}\frac{\sum_{501\leq
i\leq 700}(y_{t_{i}}-\hat{y}_{t_{i}})^{2}}{\sum_{501\leq i\leq
700}(y_{t_{i}}-\overline{y})^{2}}\Bigg{]}^{\frac{1}{2}}\right),$ (108)

where $\hat{y}_{s}$ denotes the predicted output for time instant $s$ and
$\overline{y}$ is the average of output measurements in the test set. In [56],
a Hammerstein model is derived where the static nonlinear block is a
sinusoidal map derived by curve-fitting, and the linear block is an ARX model
with estimated coefficients. We denote this method by N+ARX. In the kernel-
based methods D, E, F, and G, we have employed the TC kernel (70). Also, for
the methods which are estimating an FIR, we have set $n_{{\mathrm{g}}}=200$.
We evaluate the R-squared metric on the test data for this model and the ones
estimated by the above methods. Table II reports fitting results where one can
see that the proposed method provides more accurate fit.

Method N+ARX A B C D E F G Fit [%] 48.5 83.4 81.2 80.4 81.7 89.7 85.8 92.2

Table II: R-squared metric evaluations of test data for different
identification methods.

This is also confirmed by Figure 4 which compares the test data with the
output signals predicted by methods N+ARX, E, F, and G. It seems that for
obtaining models with more accurate predictions, one should identify a
nonlinear dynamics.

Figure 4: The figure compares the test measurement data and the predicted
values.

##  IX Conclusion

In this paper, we have considered the problem of impulse response
identification when side-information is available on the internal positivity
of the system. We have employed the realization theory of positive systems to
introduce the identification scheme in which the positivity side-information
is integrated into the identified model. The resulting formulation is in the
form of a constrained optimization over a reproducing kernel Hilbert space
endowed with a stable kernel where the constraints are suitably designed to
incorporate the positivity side-information in the solution. We have borrowed
techniques and tools from optimization theory in normed spaces to derive an
equivalent finite-dimensional convex quadratic program. This gives a
computationally tractable identification scheme that incorporates the internal
positivity side-information and has the well-known advantageous features of
kernel-based methods. We have performed a Monte Carlo numerical experiment to
compare the performance of proposed approach with FIR identification methods
considering only the external positivity feature. This has empirically studied
the impact of integrating positivity side-information in terms of estimation
bias, variance, and mean squared error. The results show that the proposed
identification approach, which integrates internal positivity, outperforms the
schemes considering only external positivity. We have observed that
incorporating internal positivity side-information reduces the estimation bias
and variance. This observation is expected since FIR external positivity
implies the weakest form of information about an internally positive system
and fails to exploit the complete information of internal positivity. We have
further evaluated the effectiveness of the proposed identification scheme
using data from a heating system experiment.

##  Appendix A Appendix

###  A-A Proof of Corollary 6

Let
${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}\in{\mathscr{P}}_{\\!{{\rho}}}$
and $G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ be the corresponding transfer
function. We know that ${\mathrm{g}}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ is
a non-negative impulse response which satisfies (4). Accordingly, due to
Theorem 2, there exist ${n_{\mathrm{x}}}\in{\mathbb{N}}$,
${\mathrm{A}}\in{\mathbb{R}}^{{n_{\mathrm{x}}}\times{n_{\mathrm{x}}}}$,
${\mathrm{b}}\in{\mathbb{R}}^{{n_{\mathrm{x}}}}$,
${\mathrm{c}}\in{\mathbb{R}}^{1\times{n_{\mathrm{x}}}}$ and $d\in{\mathbb{R}}$
such that
$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)={\mathrm{c}}(z\mathbb{I}-{\mathrm{A}})^{-1}{\mathrm{b}}+d$.
Note that since $d={g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{0}$, we know
that $d\geq 0$. Consider impulse response
${\mathrm{g}}=(g_{t})_{t=0}^{\infty}$ where $g_{0}=0$ and
$g_{t}={g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{t}$, for $t\geq 1$, and let
$G$ be the transfer function corresponding to ${\mathrm{g}}$. One can easily
see that ${\mathrm{g}}$ is non-negative, and also, we have
$G(z)=G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)-d={\mathrm{c}}(z\mathbb{I}-{\mathrm{A}})^{-1}{\mathrm{b}}$
which is a strictly proper rational transfer function. Since there exists
${a}>0$ such that
$\lim_{t\to\infty}{{\rho}}^{-t}{g}^{{\scalebox{0.5}{$({\mathcal{S}})$}}}_{t}={a}$,
we know that $\lim_{t\to\infty}{{\rho}}^{-t}(g_{t}-{a}{{\rho}}^{t})=0$.
Therefore, the spectral radius of the rational transfer function
$G(z)-{a}(1-{{\rho}}z^{-1})^{-1}$ is less than ${{\rho}}$, and consequently,
${{\rho}}$ is the unique dominant pole of $G(z)$. Hence, according to Theorem
5, $G(z)$ admits a positive realization, i.e., there exist
${m_{\mathrm{x}}}\in{\mathbb{N}}$,
${\mathrm{A}}_{+}\in{\mathbb{R}}_{+}^{{m_{\mathrm{x}}}\times{m_{\mathrm{x}}}}$,
${\mathrm{b}}_{+}\in{\mathbb{R}}_{+}^{{m_{\mathrm{x}}}}$ and
${\mathrm{c}}_{+}\in{\mathbb{R}}_{+}^{1\times{m_{\mathrm{x}}}}$ such that
$G(z)={\mathrm{c}}_{+}(z\mathbb{I}-{\mathrm{A}}_{+})^{-1}{\mathrm{b}}_{+}$.
Therefore, we have
$G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}(z)={\mathrm{c}}_{+}(z\mathbb{I}-{\mathrm{A}}_{+})^{-1}{\mathrm{b}}_{+}+d$
which says that $G^{{\scalebox{0.5}{$({\mathcal{S}})$}}}$ has a positive