Nonlinear System Identification with Prior Knowledge of the Region of Attraction
Mohammad Khosravi Roy S. Smith
math.OC cs.SY eess.SP eess.SY


Hydrogen Non-Equilibrium Ionisation Effects in Coronal Mass Ejections
    P. Pagano<ref> A. Bemporad<ref> D.H. Mackay<ref>
          
    
=====================================================================

empty
empty


We consider the problem of nonlinear system identification when prior knowledge is available on the region of attraction (ROA) of an equilibrium point. We propose an identification method in the form of an optimization problem, minimizing the fitting error and guaranteeing the desired stability property. The problem is approached by joint identification the dynamics and a Lyapunov function verifying the stability property. In this setting, the hypothesis set is a reproducing kernel Hilbert space, and with respect to each point of the given subset of the ROA, the Lie derivative inequality of the Lyapunov function imposes a constraint. The problem is a non-convex infinite-dimensional optimization with infinite number of constraints. To obtain a tractable formulation, only  a suitably designed finite subset of the constraints are considered. The resulting problem admits a solution in form of a linear combination of the sections of the  kernel and its derivatives. An equivalent optimization problem with a quadratic cost function subject to linear and bilinear constraints is derived. A suitable change of variable gives a convex reformulation of the problem. To reduce the number of hyperparameters, the optimization problem is adapted to the case of diagonal kernels. The method is demonstrate  by means of an example.	




§ INTRODUCTION

The identification of nonlinear systems has received significant attention due to its potential in modeling various phenomena in science and engineering <cit.>.
Given the measurement data, the techniques of optimization, statistics, and system identification are to mathematically model the physical systems <cit.>. In many situations modeling involves more than fitting  nonlinear dynamics to the measurement data; one should include additional features as prior knowledge which are expected according to our understanding of the system.
For example, the system properties like stability, passivity and positivity are already considered for the identification of linear dynamics <cit.>.  

For nonlinear systems, stabilizability of the dynamics is considered as a part of the identification problem in <cit.>. Identification of a stable nonlinear dynamical system has been studied in <cit.> mainly motivated by imitation learning.
In <cit.>, hidden Markov models and Gaussian mixtures are used for modeling the dynamics.
A similar approach is presented in <cit.> with guaranteed global stability.
In <cit.>, a two-stage approach is presented where, first a parametric Lyapunov function as well as a model for the dynamics are learned, and then, the learned dynamics is stabilized using the Lyapunov function.
The approach presented in <cit.> models the dynamics as a  weakly nonlinear system which consists of a stable linear part for capturing the baseline behavior, and a nonlinear part to account for more complex phenomena, and  a phase variable for the coupling these two parts. It is shown in <cit.> that the derived model is stable and time-varying.
An identification method is introduced for learning a globally stable system in <cit.>. Similar to the current paper, the hypothesis space in <cit.> is a smooth vector-valued reproducing kernel Hilbert space  (SVRKHS) <cit.>.

Meanwhile, the stability condition in <cit.> is only imposed locally over the data points by forcing the eigenvalues of the corresponding Jacobian to be negative at sampling points.
 
In this paper, we propose a nonlinear system identification method designed to include the available knowledge on a subset of the region of attraction (ROA) of a stable equilibrium point. 
Assuming that this stability property can be verified by a quadratic Lyapunov function, the problem is then formulated as a joint nonparametric estimation over a hypothesis space for the unknown dynamics, characterized here by a SVRKHS, and also, the space of positive definite matrices in order to determine the Lyapunov function. 
The resulting formulation is a non-convex optimization problem over an infinite dimensional space with infinite number of bilinear constraints, arising from the Lie derivative of the Lyapunov function with respect the points of the given subset of ROA. 
In order to make the problem tractable, we first introduce a suitable finite subset of the given subset of ROA such that verifying the Lie derivative inequality on these points guarantees the desired stability property. Following this, we reformulate the problem into a  finite dimensional optimization problem with a quadratic cost function, and linear and bilinear constraints, modeling the stability of the system at the equilibrium point and in the given region. We prove that this problem admits a solution with a linear parametric representation in terms of the sections of the kernel as well as its derivatives. 
Using a non-obvious change of variables, we derive a convex reformulation of the problem. Following this, in order to mitigate the hyperparameter tuning issue, we present the case for diagonal kernels. The method is demonstrated numerically by means of an example. 



§ NOTATIONS AND PRELIMINARIES

The set of natural numbers, the set of non-negative integers, the set of real numbers, the n-dimensional Euclidean space and the space of n by m real matrices are denoted by , ,  , ^n and ^n× m respectively.
The identity matrix and zero vector in the Euclidean space are denoted by  and  respectively.
The set of symmetric positive definite matrices in ^n× n is denoted by _++^n.
For any pair of symmetric matrices X,Y∈^n× n, we write X≽Y if X-Y∈_++^n.
Given ∈_++^n, ·_ is a norm on ^n defined as _:=^, for any ∈^n. When =, we drop subscript . The disk in ^n with center c and radius r>0 is denoted by (c,r) and defined as (c,r):={x∈^n |  x-c< r}.
For a vector α = (α_i)_i=1^n∈^n, we define |α|:=∑_i=1^nα_i≤ s and
for a function f of n variables we denote 
the partial derivative ∂^α_1/∂ x_1^α_1…∂^α_n/∂ x_n^α_n
f() by ∂_α f().
Similarly, if f is a function with k multivariable arguments, ∂^i_α f denotes the partial derivative of f with respect to the i^ argument.
The derivative operator is denoted by D, i.e., the derivative of f is shown by Df. 
The interior of set  is denoted by ^o.
Let  be a compact subset of ^n such that  is the closure of ^o.
Then, C^s(,^n) is defined as the Banach space of functions g:→^n where ∂_α g|_^o is well-defined and has a continuous extension to , for any α∈^n such that |α|≤ s. 
The norm on C^s(,^n) is defined as 

    g_C^s(,^n)=_α∈^n,|α|≤ ssup_∂_αf().

One can define C^2s(×,^n× n) similarly.
Let  be a set and  be a subset of . The indicator function of , denoted by _, is defined as
_(y) = 0, if y∈ and  _(y) = ∞, otherwise.


§ PROBLEM STATEMENT

Let  be an open domain in ^n and f:→^n be an unknown vector field defined on  which is C^2(,^n). Consider the corresponding dynamical system defined as

    ẋ = f(x),   x(0) = x_0,

where x_0∈ is the initial point. Denote the solution of (<ref>), at time instant t≥ 0, by x(t;x_0).
Let the origin be an asymptotically stable equilibrium of dynamical system (<ref>). 
Also, let the corresponding region of attraction (ROA) be denoted by , i.e., we have 

    := {x_0∈ | lim_t→∞x(t;x_0)=}.

Let Ω be a known inner approximation for the region of attraction of the origin. 
More precisely, Ω is a known compact set with non-empty interior such that ∈Ω⊆.

Consider a set of trajectories of system,   
like {x(·;x_0^i)  | 1≤ i≤}, where the corresponding initial points belong to , i.e.,  {x_0^1,…,x_0^}⊂.
For any i=1,…,, suppose that the i^ trajectory is sampled at time instants 
0≤ t^i_1<t^i_2<⋯<t^i_n_i where n_i∈.
Let x_k^i denote x(t_k^i;x_0^i), for 1≤ k≤ n_i.
Given these samples, one can estimate the time derivative of x(·;x^i_0) at the sampling time instants.
This estimation can be performed using various techniques, e.g. see <cit.> and the references therein, or by simply utilizing a nonlinear regression method and subsequently obtaining the derivatives numerically or analytically. 
Let these estimations be denoted by y_k^i, for 1≤ k≤ n_i.
One should note that y_k^i is an approximation of f(x_k^i).
Considering these samples of trajectories and their estimated derivatives, we get a set of data, denoted by , which contains
(x_k^i,y_k^i) pairs. For notation simplicity, we drop the superscripts and simply show set  as {(x_j,y_j) | 1≤ j≤ m}, where m=∑_i=1^n_i.


Given that origin is a stable equilibrium point of (<ref>) and 
the set Ω is provided as the prior knowledge about the region of attraction of the origin, the problem is to estimate the unknown vector field f, in a given class of functions ⊆ C^2(,^n), using
the set of data .

In the next section, we introduce a tractable formulation of this problem as a nonparametric estimation.
The formulation can be extended to the case of multiple equilibria and multiple regions of attraction.

	


§ MAIN RESULTS: IDENTIFICATION METHOD

We know that f satisfies the constraint that x= is an equilibrium point of (<ref>), 
i.e.,
f()=.
Moreover, we know that x= is stable and Ω is a subset of the corresponding region of attraction. 
Assume that these stability features of f can be verified by an unspecified quadratic Lyapunov function V(x) = 1/2x^Px where here P is a positive definite matrix.
More precisely, there exist an unknown ϵ>0 and an unknown positive definite matrix P≽ such that 

    DV(x) f(x) = x^Pf(x)≤ -ϵx^2,   ∀x∈Ω.


In the estimation problem, we need to minimize the fitting error, ∑_i=1^m y_i-f(x_i)_^2, subject to f∈, f() = and (<ref>).
There are two main issues: the correct choice of function class , and dealing with the (uncountable) infinite number of constraints introduced in (<ref>).
These issues are addressed in this section.



 §.§ From Infinite to Finite Number of Constraints

Since an optimization problem with an infinite-dimensional feasible set and infinite  number of imposed constraints can lead to intractability, particularly when the problem is non-convex as it is here, we need to introduce a finite number of suitable constraints implying (<ref>).
To this end, the notion of (α,β)-grid is introduced in the next definition.
One should note that, based on the discussed given below, the  (α,β)-grid is a suitably selected finite subset of Ω such that verifying stability condition 

on its elements implies the desired stability feature given in (<ref>).


Let {z_1,…,z_} be a finite subset of Ω\{} denoted by . We say   is an (α,β)-grid for Ω if 

    Ω⊆(⋃_i=1^(z_i,αz_i) )∪(,β).


The role of (α,β)-grid in the estimation problem is shown by the next theorem. 
Define L_1,f and L_2,f respectively as

    [                                       L_1,f                                          :=                               sup_x∈ΩDf(x),;                                       L_2,f                                          := sup_x∈Ωsup_h_1,h_2∈(,1)
    |D^2f(x)(h_1,h_2)|.;                                             ]
	
Since f is C^2(,^n), we have that L_1,f,L_2,f<∞.


	Let :={z_1,…,z_}⊂Ω\{} and P∈^n be such that 
	
    [   z_i^Pf(z_i)≤ -z_i^2,            ∀ i = 1,…,,; 1/2(PDf()+Df()^P) ≼ -,;                    P≽, ]

	Given ϵ∈(0,1), let α and β be real positive scalars where
	
    α≤(1+L_1,fP/ϵ+L_1,fP)^1/2   -1,
    	  β≤1-ϵ/L_2,fP,

	and  be an (α,β)-grid for Ω.
	Then, the following holds
	
    x^Pf(x)≤ -ϵx^2,   ∀x∈Ω.



Let x∈Ω. Since  is an (α,β)-grid, then by (<ref>), we know that x either belongs to  (,β), or it belongs to (z_i,αz_i), for some 1≤ i≤.

Case I: Assume that x∈(,β), i.e., x<β. Since, f is C^2(,_≥ 0^n) and f()=, from Taylor expansion at origin, we have that

    f(x) = Df()x + r(x),

where r:→^n is a C^2(,^n) function such that

    sup_x∈Ω, xr(x)/x^2≤ L_2,f.

Accordingly, one can easily see that

    x^Pf(x) =
    1/2x^(PDf()+Df()^P)x + x^P r(x).

Due to  (<ref>), (<ref>), and the Cauchy-Schwartz inequality, we have

    x^Pf(x) 
    ≤ 
    -x^2
    + 
    PL_2,fx^3.

According to (<ref>) and since x< β, we have

    x^Pf(x)≤ -ϵx^2.

Case II: Assume that  x∈(z_i,αz_i), i.e., x = z_i + a where a is a vector such that a<αz_i.
From the triangle inequality, we have that

    -ϵz_i^2 (1+α)^2
    ≤
    -ϵ(z_i + a)^2
    ≤ -ϵx^2.

Due to (<ref>), L_1,f is a Lipschitz constant for f.
Therefore, we have

    f(z_i + a)-f(z_i )≤ L_1,fa,

and

    f(z_i + a)≤ 
    L_1,f(z_i+a),

where the second inequality follows from f()= and the triangle inequality.
From (<ref>) and the Cauchy-Schwartz inequality,
we have 

    x^Pf(x) 
       =z_i^Pf(z_i) 
    +z_i^P(f(z_i + a)-f(z_i ))
    +a^Pf(z_i + a)
    
       ≤ -z_i^2 + Pz_if(z_i + a)-f(z_i )
    +
    Paf(z_i + a).

Note that Pa≤Pa and Pz_i≤Pz_i, for any i=1,…,. 
Therefore, from (<ref>) and (<ref>), 
we have 

    x^Pf(x) 
       ≤ 
    -z_i^2 + Pz_if(z_i + a)-f(z_i )
    +
    Paf(z_i + a)
       ≤ 
    -z_i^2 + Pz_i L_1,fa
    +
    Pa L_1,f(z_i+a).

Since a<αz_i, one can conclude that

    x^Pf(x) 
    < -z_i^2 (1 - 2PL_1,fα 
    -PL_1,fα^2).

From (<ref>), one can see that

    1 - 2PL_1,fα 
    -PL_1,fα^2 ≥ϵ(1+α)^2.

Due to (<ref>), 
(<ref>),
and (<ref>), we have that

    x^Pf(x)
    <
    -ϵ(1+α)^2z_i^2
    ≤ -ϵx^2.

This concludes the proof.


	Theorem <ref> guarantees that in order to satisfy the infinite number of constraints given in (<ref>), it is enough to satisfy the finite number of constraints introduced in (<ref>), given a suitable (α,β)-grid for Ω. 


	One can verify that if (<ref>) holds for a positive-definite matrix P, then a scaled version of P satisfies the inequalities given in (<ref>). 
	This can be interpreted as the other direction of Theorem <ref>.	


 
	For any α,β>0, an (α,β)-grid exists for Ω.


Since 

    Ω⊆∪_z∈Ω,z(z,αz)∪(z,β)

and Ω is a compact set, this open cover  has a finite sub-cover.
Hence, there exist {z_1,…,z_}⊂Ω\{} such that (<ref>) holds.	



	For the given Ω, there are infinite choices for (α,β)-grids.
	Moreover, one can see that taking small values for α and β results in fine and large (α,β)-grid.  	






 §.§ Identifying the Dynamics in the Smooth Vector-valued Reproducing Kernel Hilbert Spaces

The function class taken for approximating the unknown vector field is a type of Hilbert spaces called smooth vector-valued reproducing kernel Hilbert spaces (SVRKHS) which are introduced below (see <cit.> for more details). Based on the suitable structure of SVRKHS, we will prove 
that the problem admits a solution with a specific finite linear parametric form. This allows reducing the optimization problem to the coefficients of this representation and subsequently, a tractable finite-dimensional optimization problem is obtained.

Let  be a compact subset of ^n with non-empty interior ^o such that  is the closure of ^o  
and  Ω⊂^o.

A Smooth Vector-valued Reproducing Kernel Hilbert Space (SVRKHS), denoted by , is a Hilbert space of functions g∈ C^s(,^n) such that 
for any ∈, 

we have sup_g∈,g_≤ 1g()<∞.


The function ∈ C^2s(×,^n× n) is an operator-valued positive-definite Mercer kernel <cit.> when for any m∈, ,,_1,…,_m∈ and _1,…,_m∈^n, we have (,)=(,)^ and
∑_1≤ i,j ≤ m_i^(_i,_j)_j ≥  0.

For any ∈, let _ denote the function defined by (,·):→^n.
This is called the section of kernel  at  or the feature map.

With respect to any Mercer kernel ∈ C^2s(×,^n× n), there exists a SVRKHS of functions g∈ C^s(,^n), denoted by  and endowed with inner product ··_ and norm ·_, such that for any (,)∈×^n and  for any α∈^n with |α|≤ s, we have

	
  i) ∂_α^1 _∈, and
	
  ii) g∂_α^1 __=^∂_αg(), for all g∈.
 
The second feature is called the reproducing property.

We suppose that the kernel  is suitably chosen such that function h:→^n, defined as h():=, belongs to . A simple example is (,) defined as (,) where  is a polynomial kernel. 
Also, let assume ∈^n_++ is a positive-definite 
finite-dimensional transformation such that
f∈, i.e.,  is a positive-definite change of coordinates on ^n which transforms the vector filed f to an element of .
More precisely, we know that g∈ where g:Ω→^n is defined as g() = ( f()), for any ∈Ω. For example  might be a scaling of the identity matrix. Note that this is mainly a technical assumption which is used later to simplify the mathematical arguments.

Let the fitting loss or the error function, denoted by _,, be the function _:→ defined as

    _(f) := ∑_i=1^y_i-f(x_i)_^2,
 
where ∈_++^n is an error weighting matrix. 
Additionally, we can consider a suitable kernel-based regularization due to f∈. More precisely, let the regularization function :→_≥ 0 be defined as ( f)  :=   f_^2.
The identification problem is now formulated as following 

    min_f∈ ∑_i=1^_,(f) + λ( f)

where λ>0 is the regularization weight and  is the set of smooth vector fields such that x=0 is a stable equilibrium point and  attractive in the region Ω, and also, for any f∈, we have that f∈. 
Note that (<ref>)
is a non-convex optimization problem with an infinite-dimensional feasible set and infinite  number of constraints.
In the followings, we show that this problem has a tractable reformulation.

Due to Theorem <ref>, for imposing the stability feature given in (<ref>), it is sufficient to take a suitable (α,β)-grid, like = {z_1,…,z_}, and  solve optimization problem over the grid,

    [                     _f∈ C^2(,^2), P∈_++^n               _i=1^y_i-f(x_i)_^2 + λ f^2_;                                      s.t. z_i^Pf(z_i)≤ -z_i^2, 
            ∀ i = 1,…,,;                    1/2(PDf()+Df()^P) ≼ -,;                                   f() = ,;                                       f∈,;                                       P≽. ]

The existence of such a grid is guaranteed by Theorem <ref>.
Rewriting  optimization problem (<ref>) in terms of g, one has

    [                                     _g∈, P≽               _i=1^y_i-^-1g(x_i)_^2 + λg^2_;                                        s.t. z_i^P^-1g(z_i)≤ -z_i^2,         ∀ i = 1,…,,;                1/2(P^-1Dg()+Dg()^^-1P) ≼ -,;                                     g() = . ]

The problem (<ref>) in a non-convex infinite-dimensional optimization and therefore, it is not tractable. However, in order to address this issue, we derive a finite dimensional problem equivalent to (<ref>).

With respect to a given ∈_++^n, we define _ as

    _ := {g∈ |    z_i^P^-1g(z_i)
    ≤ -z_i^2, 
     ∀ i = 1,…,,
    
       1/2(P^-1D
    g()+Dg()^^-P) ≼  -,
    
       g() = }.


For any ≽, the set _ is a non-empty, closed and convex subset of .	


Since  and  are positive definite matrices, all of the eigenvalues of matrix M defined as M:=1/2(^-1+^-) are strictly larger than zero.	
Let the function g_γ:→^n be defined as g_γ()=-2γ^-1 where γ is a positive real scalar smaller than smallest eigenvalue of M. Since -1/2γ g_γ∈, we know that g_γ∈.
We have

    _i^^-1g_γ(_i)=
    -2γ^-1_i^M_i≤ -z_i^2,

and 

    1/2(^-1Dg_γ()+Dg_γ()^^-1) = -2γ^-1M≼ -.

Moreover, we know that g_γ()=.	Therefore, g_γ∈_ and thus, _ is non-empty.
The convexity of _ is due to the linear dependency of the left-hand sides of the constraints with respect to g. 
Now, let (g_k)_k=1^∞∈_ be a sequence converging to g∈.

Let  be an arbitrary vector in ^n.
Then, for any ∈ and any α∈^n, due to the reproducing property and the Cauchy-Schwarz inequality, we have

    |^(∂_αg()-∂_αg_k())| 
       =
    |g-g_k∂_α^1 __|
    
       ≤g-g_k_∂_α^1 __.

This shows that lim_k→∞∂_αg_k()=∂_αg(). 
Therefore, we have 

    lim_k→∞g_k()=g(),
    
        lim_k→∞Dg_k()=Dg(),

and

    lim_k→∞∂_αg_k(_i) =∂_αg(_i),     ∀ i=1,…,.

Since for any k∈, g_k satisfies the constraints and the left-hand sides of the constraints depend linearly on g_k and g,  it follows that, the constraints are also satisfied by g, i.e., g∈_. Hence, _ is a closed subset of .

For ease of notation, 
define p:=+, m := 1+p + n, and also, set
_0 =, and _i = _i-, for 1+≤ i ≤ p. 

For any ≽ and λ>0, the optimization problem 

    _g∈_ ∑_i=1^y_i-^-1g(x_i)_^2 + λg^2_,

has a unique solution, denoted by g^*_.
Moreover, there exist vectors {_i}_i=0^p and {_j}_j=1^n such that, g^*_,
the solution of (<ref>), is in the  following form

    g^*_ 
    = 
    ∑_i=0^p__i_i
    +
    ∑_j=1^n∂_j^1__j.



Define :→∪{+∞} as

    (g) 
    := ∑_i=1^y_i-^-1g(x_i)_^2  + λg^2_ +  __(g),

for any g∈. 
According to Theorem <ref>, _ is a non-empty, closed, and convex set, and therefore  __ is a proper lower-semicontinuous convex function <cit.>. 
Let g_γ be the element of  introduced in the proof of Theorem <ref>.
Since g_γ∈_, we have __(g_γ)=0.
Also, we know that ∑_i=1^y_i+2γ^-1^-1x_i_^2<∞.
Therefore, ∑_i=1^y_i-^-1g(x_i)_^2 is a proper and continuous convex function with respect to g. Since λ>0 and g_γ_<∞, we have that  is a proper and lower-semicontinuous strongly convex function.
Therefore, min_g∈(g) has a unique (finite) solution <cit.>, which means that (<ref>) admits a unique solution with finite cost.  
Define set ⊆ as

    := {∑_i=0^p__i_i
    +
    ∑_j=1^n∂_j^1__j 
    |_0,…,_p,_1,…,_n∈^n}.

This is a finite-dimensional subspace of  and consequently, it is a closed subspace. Hence, one can decompose g_^* as g_^* = g_^∥+g_^⊥ where  g_^∥∈ and g_^⊥∈^⊥.
Therefore,  for any ∈^n and for any i=0,…,p, from the reproducing property, we have ^ g_^⊥(_i) = g_^⊥__i_=0 and subsequently, we have ^ g_^∥(_i) = ^ g_^*(_i). Accordingly, one can conclude that g_^∥(_i) = g_^*(_i), for any i=0,…,p.
Similarly, due to the reproducing property, we have

    [                              ^Dg_^⊥()                                     = [
    ^∂^1_1 g_^⊥(), 
    …, 
    ^∂^1_n g_^⊥() ];                                                                           =       [
    g^⊥_∂^1_1__, 
    …, g^⊥_∂^1_n__];                                     =                                    ^, ]

which shows that Dg_^⊥() is zero and subsequently, Dg_^∥() =Dg_^*(). As g_^*∈_, it follows that g_^∥∈_. Also, we have 

    ∑_i=1^y_i-^-1g^∥(x_i)_^2 = ∑_i=1^y_i-^-1g(x_i)_^2,

and 

    g_^2 = g^∥_^2 +g^⊥_^2 ≥g^∥_^2.
 

We need to have g^⊥_^2=0, otherwise g^∥ is a feasible solution with objective value strictly smaller than minimum of the objective function. This means that g = g^∥∈ and g^*_ has the form given in (<ref>).
 
For simplicity of notation, we define _i = _i-p-1, for m-n≤ i ≤ m-1, and vector ∈^nm  as
:=[_0^…_m-1^]^.

In the next theorem, we introduce a finite dimensional version of (<ref>).
First, we need to introduce required notations.
Define matrix ∈^mn× mn and _i∈^n× mn respectively as
:=[_i_1,i_2]_i_1=0,i_2=0^m-1,m-1 and 
_i:=[_i,0,…,_i,m-1], for any i=0,…,m, 
where _i_1,i_2 is given, for any i_1,i_2=0,…,m, as

    _i_1,i_2  :=  
    (_i_2,_i_1),
          
    0≤ i_1,i_2≤ p,
    ∂^1_i_2-p(,_i_1),
           
    0≤ i_1≤ p<i_2< m,
    ∂^2_i_1-p(_i_2,),
           
    0≤ i_2≤ p<i_1< m,
    ∂^2_i_1-p∂^1_i_2-p(,),
           
    p+1≤ i_1,i_2< m.

Define the function :→^n× nm as

    () :=
    [__0(), 
    …, 
    __p(),
    ∂_1^1_(), 
    …, 
    ∂_n^1_() 
    ],     ∀x∈.

Based on this definition, one can easily see that _i= (_i), for any i=0,…,m.
Due to the reproducing property, we have the following proposition.

	Let ∈^nm and g:→^n be defined such that g()=(), for any x∈. 
	Then, we have that g∈ and 
	g_^2 =  ^.
	

Consider the following optimization problem

    [                           _≽, ∈^nm
    ∈^n× n                    _i=1^y_i-^-1_i_^2 + λ^;                                      s.t. _i^^-1_i≤ -_i^2, 
            ∀ i = p-+1,…,p,;       = 
    [ _m-n+1      …   _m-1     _m ],;                    1/2(P^-1D+D^^-1P) ≼ -,;                                    _0 = . ]

Then, for each solution of (<ref>), one can find a solution for (<ref>) with same cost value.


One can restate optimization problem (<ref>) in the following form

    _P≽ (_g∈_ _i=1^y_i-^-1g(x_i)_^2 + λg^2_).

From Theorem <ref>, we know that the solution of the inner problem is in the form of g()=(), for a vector ∈^mn. Utilizing the reproducing property of the kernel and according to (<ref>), substituting g()=() in (<ref>), results in optimization problem (<ref>). Solving (<ref>), we obtain ,,, and subsequently, g. 

Due to Theorem <ref>, in order to solve (<ref>), it is sufficient to find the solutions of (<ref>), which is a finite-dimensional optimization. 
One can see that  (<ref>) is a finite dimensional optimization problem with a convex cost and, linear and bilinear constraints, and subsequently,  it is not a convex optimization.  

Note that the kernel  is characterized by a number of constants called hyperparameters. Here, this is implicitly assumed and for the sake of more transparent discussion, we have dropped this dependency in the notations. 
The hyperparameters are required to be estimated based on the data. 
This is commonly done using a cross-validation routine. 
The hyperparameter estimation  is essentially a computationally demanding procedure, especially when the kernel has a large number of hyperparameters, which can be the case when operator-valued kernels are used. 
The non-convexity of problem increases the computational complexity of the hyperparameters estimation to the point of potential intractability.
These issues will be addressed in the following section.





 §.§ Diagonal Kernels

In order to alleviate the issue of having a large number of hyperparameters, we take the kernel  as a diagonal kernel in the form of (,)=(,), where ∈ C^2s(×,) is a scalar valued Mercer kernel. 
Define :→^m as 

    ():=[    __0()        …    __p() ∂_1^1_()        … ∂_n^1_() ]^,     ∀x∈.

By defining the matrix ∈^n× m as
:=[_0…_m-1], one can see that for the unique  solution of (<ref>), we have g^*_ =. 
Similar to  and {_i}_i=0^m, we define matrices ∈^m× m and _i∈^m× 1 respectively as
:=[_i_1,i_2]_i_1=0,i_2=0^m,m and 
_i:=[_i,0,…,_i,m]^, for any i=1,…,m, 
where, for each i_1,i_2=0,…,m, _i_1,i_2  is defined similarly to (<ref>) but based on . Also, let ∈^m× n be the Jacobian or derivative of  at =, i.e.,

    :=
    [ ∂_1()     … ∂_n() ].

Accordingly, we have Dg()=.
Due to the reproducing property, we have the following proposition.

Let ∈^m× n and g:→^n be defined such that g()=(), for any x∈. 
Then, we have that g∈ and 
g_^2 =  (^).
	

Based on the discussion above, analogous to 
(<ref>), we can introduce the following  optimization problem

    [                                   _≽
     ∈^n× m                    _i=1^y_i-^-1_i_^2 + λ (^);                                         s.t. _i^^-1_i≤ -_i^2,             ∀ i = p-+1,…,p,;                      1/2(P^-1 + ^^^-1P) ≼ -,;                                       _0 = . ]

One should note that, the complexity of this optimization problem, in terms of the number of variables and constraints, is significantly lower than (<ref>), especially when the dimension of the state space, n, is large.





 §.§ Towards a Convex Formulation

In the estimation problem (<ref>), and subsequently (<ref>), the matrices W and  are introduced as arbitrary positive definite matrices. One can see that the mathematical arguments (up to Theorem   <ref>) only require the fact that  and  do not depend on g. This provides the opportunity of choosing them such that a change of variables lead to a convex formulation. In fact, we set :=^2 and :=. Note that for any ∈^n, one has _^2=. Accordingly, we have

    ∑_i=1^y_i-^-1g(x_i)_^2 =
    ∑_i=1^y_i-g(x_i)^2.
 
Therefore, optimization problem (<ref>) can be modified to give

    [                                  _g∈, P≽               _i=1^Py_i-g(x_i)^2 + λg^2_;                                     s.t. z_i^ g(z_i)≤ -z_i^2,         ∀ i = 1,…,,;                     1/2(Dg()+Dg()^) ≼ -,;                                  g() = . ]

One can see that (<ref>) is a convex optimization problem which by  Theorem <ref> has a unique solution of the form

    g^* 
    = 
    _i=0^p__i_i
    +
    _j=1^n∂_j^1__j.

By substituting the solution (<ref>) into (<ref>), we obtain a finite problem (analogous to (<ref>)) as

    [                                     _≽, ∈^nm
    ∈^n× n                                  _i=1^y_i-_i^2 + λ^;                                                s.t.           _i^_i≤ -_i^2,         ∀ i = p -  + 1,…,p,; = 
    [ _m - n + 1          …       _m-1         _m ],;                                      1/2(D+D^) ≼ -,;                                              _0 = . ]

In the case of diagonal kernels, the modified version of (<ref>) is the following

    [                                _∈^n ×  m , ≽                       ∑_i=1^y_i-_i^2 + λ (^);                                         s.t. _i^_i≤ -_i^2 , 
            ∀ i = p -  + 1,…, p,;                                1/2(+^^) ≼ -,;                                       _0 = . ]


Based on the above discussion, we have the following proposition.

The optimization problems (<ref>) and (<ref>) are feasible convex programs with strongly convex objective functions. Therefore, each of (<ref>) and (<ref>) admits a unique solution.


Due to convexity of optimization problems (<ref>) and (<ref>), one can use off-the-shelf optimization solvers, like CVX <cit.>, in order to obtain the solutions.


One should note that once  and g are obtained, the vector field f is calculated as ^-1g.


	


§ NUMERICAL EXPERIMENTS

In this section, we discuss a numerical example.
To this end, consider the dynamical system
defined as 

    [                               ẋ_1                                 =     -5x_2-4x_1  +x_1x_2^2-6x_1^3,;                               ẋ_2                                 = -20x_1-4x_2 + 4 x_1^2x_2 + x_2^3.;                                   ]

For system (<ref>), (x_1,x_2)= is a stable equilibrium point, which is attracting in region Ω defined as Ω:=(,1.5).
Assume that this information is provided as prior knowledge. 
Moreover, we consider two trajectories of the system, starting from  _0=(1,-1) and  _0=(-1,-1)), and take samples from each of them at different 19 locations with an additive measurement noise of (0,σ^2) where σ^2 = 0.001. 
For identifying (<ref>), we take two approaches:  1) we utilize the prior knowledge on the ROA and solve  (<ref>), 
2) only the stability of the equilibrium point is considered and we solve a modified version of (<ref>) where = and the grid constraints are removed.
Let the corresponding solutions be denoted by f̂ and f̌, respectively. 
The main difference of these two approaches is the inclusion of the prior knowledge on the region of attraction in the estimation method. Accordingly, the comparison of f̂ and f̌ can reflect the impact and the potential leverage of using the prior knowledge of the ROA on the estimation. 
For , we take a uniform polar grid of size =300 inside Ω.
Given these settings, we obtain estimations f̂ and f̌. The results are shown in Figure <ref>. The calculated coefficient of determination in the unit square, also known as R squared, for f̂ and f̌ is 93.6% and  80.4%, respectively.




Let _f, _f̂ and _f̌ denote the trajectories generated from f, f̂ and f̌, respectively.
We consider initial points _0=(-1,1),(-1.45,0),(1,1.5). The corresponding trajectories are shown in Figure <ref>. For point _0=(-1,1)∈Ω, all of the trajectories goes to the equilibrium point (0,0). One can see that trajectory _f̂ stays close to the trajectory of true system, _f, while _f̌ deviates from _f significantly.
Starting from _0=(-1.45,0)∈Ω, trajectories  _f and _f̌ stay close to each other and converge to (0,0), meanwhile _f̌ diverges. This confirms that the prior knowledge is satisfied by the estimated vector field f̂.
Finally, if _0=(-1.45,0)∉ROA, trajectory _f as well as trajectory _f̂ diverge. However, _f̌ converges to (0,0) which is not expected. 


	



§ CONCLUSION

We have discussed nonlinear system identification  when in addition to the measurement data, prior knowledge is available on a subset of the region of attraction (ROA) of an equilibrium point. The proposed identification method is an optimization problem minimizing the fitting error and guaranteeing the desired stability property. The resulting problem is a joint identification of the dynamics as well as a Lyapunov function for the stability property. 
Due to the functional hypothesis space for the dynamics and the Lie derivative inequalities for the stability, a bilinear infinite-dimensional optimization problem with infinite number of constraints is obtained. 
To get a tractable formulation, we consider a sufficient suitable finite subset of the constraints. The resulting problem admits a solution in form of a linear combination of the sections of the kernel and its derivatives. Subsequently, by a change of variable, we obtain a convex reformulation of the problem. Moreover, for reducing the number of hyperparameters, the optimization problem is adapted to the case of diagonal kernels. 
We have verified the approach and illustrated the results on an example. 
The fitting error for the estimation is low, and the estimated vector field confirms the expected behavior in the given subset of the ROA. 
In order to assess the impact of the prior knowledge, we have compared the method with an approach where the prior knowledge is not exploited. The comparison confirms the significance of the impact of the prior knowledge on the precision of the estimation and on the global behavior of the estimated system.








	

IEEEtran

