Empirical Differential Gramians for Nonlinear Model Reduction

TL;DR

This paper introduces an empirical method for nonlinear model reduction using differential Gramians, enabling balanced truncation along fixed trajectories without solving complex PDEs, demonstrated on an RL network.

Contribution

It presents a novel empirical approach to compute differential Gramians for nonlinear systems, simplifying the model reduction process along specific trajectories.

Findings

Balanced truncation is feasible without PDEs for fixed trajectories.

The method relies on impulse and initial state responses of the variational system.

Demonstrated successfully on a reinforcement learning network.

Abstract

In this paper, we present an empirical balanced truncation method for nonlinear systems with linear time-invariant input vector field components. First, we define differential reachability and observability Gramians. They are matrix valued functions of the state trajectory (i.e. the initial state and input trajectory) of the original nonlinear system, and it is difficult to find them as functions of the initial state and input. The main result of this paper is to show that for a fixed state trajectory, it is possible to compute the values of these Gramians by using impulse and initial state responses of the variational system. Therefore, balanced truncation is doable along the fixed state trajectory without solving nonlinear partial differential equations, differently from conventional nonlinear balancing methods. We further develop an approximation method, which only requires trajectories of the original nonlinear systems. Our methods are demonstrated by an RL network along a trajectory.

Equations

$$\displaystyle\left\{\begin{array}[]{l}\dot{x}(t)=Ax(t)+Bu(t),\\ y(t)=Cx(t),\end{array}\right.$$

$$x(t)=e^{A(t-t_0)}x(t_0)+{\int }_{t_0}^{t}B{e}^{A(t-\tau )}u(\tau )d\tau .$$

$$G_c(t_0,t_f):={\int }_{t_0}^{t_f}{e}^{A(t-t_0)}B{B}^{\mathrm{T}}{e}^{A^{\mathrm{T}}(t-t_0)}dt,$$

$$G_o(t_0,t_f):={\int }_{t_0}^{t_f}{e}^{A^{\mathrm{T}}(t-t_0)}{C}^{\mathrm{T}}C{e}^{A(t-t_0)}dt.$$

$$T{G}_c(0,\infty )T^{\mathrm{T}}$$

$$=\mathrm{diag}\{{\sigma }_1,\dots ,{\sigma }_n\},{\sigma }_i\ge {\sigma }_{i+1}.$$

$$G_c(t_0,t_f)={\int }_{t_0}^{t_f}{x}_{\mathrm{Imp}}(t)x_{\mathrm{Imp}}^{\mathrm{T}}(t)dt.$$

$$\displaystyle y_{\rm Is}(t):=[\begin{array}[]{ccc}y_{{\rm Is},1}(t)&\cdots&y_{{\rm Is},n}(t)\end{array}]=Ce^{A(t-t_{0})}.$$

$$G_c(t_0,t_f)={\int }_{t_0}^{t_f}{y}_{\mathrm{Is}}^{\mathrm{T}}(t)y_{\mathrm{Is}}(t)dt.$$

$$\displaystyle\Sigma:\left\{\begin{array}[]{l}\dot{x}(t)=f(x(t))+Bu(t),\\ y(t)=h(x(t)),\end{array}\right.$$

$$\displaystyle d\Sigma:\left\{\begin{array}[]{l}\displaystyle\delta\dot{x}(t):=\frac{d\delta x(t)}{dt}=\frac{\partial f(\varphi_{t-t_{0}})}{\partial x}\delta x(t)+B\delta u(t),\\[5.69054pt] \displaystyle\delta y(t)=\frac{\partial h(\varphi_{t-t_{0}})}{\partial x}\delta x(t),\end{array}\right.$$

$$\frac{d}{dt}\frac{\partial {\phi }_{t-\tau }(x_{\tau },u)}{\partial x_{\tau }}$$

$$=\frac{\partial f({\phi }_{t-\tau }(x_{\tau },u))}{\partial x_{\tau }}$$

$$=\frac{\partial f({\phi }_{t-\tau }(x_{\tau },u))}{\partial {\phi }_{t-\tau }}\frac{\partial {\phi }_{t-\tau }(x_{\tau },u)}{\partial x_{\tau }}$$

$$\delta x(t)=\frac{\partial {\phi }_{t-t_0}(x_0,u)}{\partial x}\delta x_0+{\int }_{t_0}^t\frac{\partial {\phi }_{t-\tau }(x(\tau ),u)}{\partial x}B\delta u(\tau )d\tau .$$

$$\delta y(t)=$$

$$G_{\mathcal{R}}(t_0,t_f,x_0,u):={\int }_{t_0}^{t_f}\frac{\partial {\phi }_{t-t_0}}{\partial x}B{B}^{\mathrm{T}}\frac{{\partial }^{\mathrm{T}}{\phi }_{t-t_0}}{\partial x}dt,$$

$$G_{\mathcal{O}}(t_0,t_f,x_0,u)$$

$$:={\int }_{t_0}^{t_f}\frac{{\partial }^{\mathrm{T}}{\phi }_{t-t_0}}{\partial x}\frac{{\partial }^{\mathrm{T}}h({\phi }_{t-t_0})}{\partial {\phi }_{t-t_0}}\frac{\partial h({\phi }_{t-t_0})}{\partial {\phi }_{t-t_0}}\frac{\partial {\phi }_{t-t_0}}{\partial x}dt,$$

$$G_{\mathcal{R}}(t_0,t_f,x_0,u)={\int }_{t_0}^{t_f}$$

$$\Lambda =\mathrm{diag}\{{\sigma }_1,\dots ,{\sigma }_n\},{\sigma }_1\ge \cdots \ge {\sigma }_n>0$$

$$T_{\phi }{G}_{\mathcal{R}}(t_0,t_f,x_0,u)T_{\phi }^{\mathrm{T}}$$

$$\delta x_{\mathrm{Imp},i}(t)=\frac{\partial {\phi }_{t-t_0}(x_0,u)}{\partial x}{B}_i,$$

$$G_{\mathcal{R}}(t_0,t_f,x_0,u)={\int }_{t_0}^{t_f}\delta x_{\mathrm{Imp}}(t)\delta x_{\mathrm{Imp}}^{\mathrm{T}}(t)dt,$$

$$\displaystyle\hskip 14.22636pt\delta x_{\rm Imp}(t):=[\begin{array}[]{ccc}\delta x_{{\rm Imp},1}(t)&\cdots&\delta x_{{\rm Imp},m}(t)\end{array}].$$

$${\bar{G}}_{\mathcal{R}}(t_0,t_f,x_0,u)$$

$$={\int }_{t_0}^{t_f}\frac{\partial {\psi }_{t-t_0}}{\partial x}\frac{\partial f({\psi }_{t-t_0},u)}{\partial u}\frac{{\partial }^{\mathrm{T}}f({\psi }_{t-t_0},u)}{\partial u}\frac{{\partial }^{\mathrm{T}}{\psi }_{t-t_0}}{\partial x}dt.$$

$$\delta {\bar{x}}_{\mathrm{Imp}}(t)$$

$$={\int }_{t_0}^t\frac{\partial {\psi }_{t-\tau }(x(\tau ),u)}{\partial x}\frac{\partial f({\psi }_{\tau -t_0}(x_0,u),u)}{\partial u}{\delta }_D(\tau -t_0)d\tau$$

$$=\frac{\partial {\psi }_{t-t_0}}{\partial x}\frac{\partial f(x_0,u(t_0))}{\partial u}.$$

$$\delta y_{\mathrm{Is},i}(t)=\frac{\partial h({\phi }_{t-t_0}(x_0,u))}{\partial x}\frac{\partial {\phi }_{t-t_0}(x_0,u)}{\partial x}{e}_i^n,$$

$$G_{\mathcal{O}}(t_0,t_f,x_0,u)={\int }_{t_0}^{t_f}\delta y_{\mathrm{Is}}^{\mathrm{T}}(t)\delta y_{\mathrm{Is}}(t)dt,$$

$$\displaystyle\hskip 14.22636pt\delta y_{\rm Is}(t):=[\begin{array}[]{ccc}\delta y_{{\rm Is},1}(t)&\cdots&\delta y_{{\rm Is},n}(t)\end{array}].$$

$$d{\Sigma }_{(x_0,u)}(\delta x_0,\delta u)$$

$$:=\underset{s\to 0}{\lim }\frac{\Sigma (x_0+s\delta x_0,u+s\delta u)-\Sigma (x_0,u)}{s}$$

$$d{\Sigma }_{(x_0,u)}(\delta x_0,\delta u)$$

Full text

††thanks: A preliminary version of this paper is presented at the 20th IFAC World Congress, July 2017. Corresponding author Y. Kawano. Tel. +81 82 424 7582. Fax +81 82 422 7193.

,

1 Introduction

Along with the development of new technologies, control systems are becoming more complex and large-scale. To capture systems’ components which are essential for controller design and analysis, model order reduction techniques have been established, see e.g. [1]. In systems and control, typical methods are balanced truncation and moment matching [1, 43], and both of them have been extended to nonlinear systems [33, 8, 24, 3, 2, 17]. In contrast to successive theoretical developments, nonlinear model reduction methods still have computational challenges, since they require solutions to nonlinear partial differential equations (PDEs). There are few papers tackling this challenging problem such as [9, 31, 28, 32, 20]. As a data driven model order reduction method, proper orthogonal decomposition (POD) [16, 1] is often used in practice. However, POD is mainly proposed for non-control systems.

For linear time-invariant (LTI) systems, POD and balancing are connected based on the fact that the controllability and observability Gramians can be computed by using impulse and initial state responses, respectively. That is, balanced truncation of LTI systems can be performed by using empirical data. Applying linear empirical methods to nonlinear systems have attracted various research interests, see e.g., [27, 12, 6, 15, 11, 42]. Such methods are exploited in order to reduce the computational complexity of nonlinear controller design such as model predictive control [13, 5].

However, these empirical methods have been proposed only around a steady-state because the aforementioned nonlinear balancing method gives the same reduced order model as the linear balancing method at a steady-state. For analysis and control of nonlinear systems, a steady-state is not always important. For instance, in a trajectory tracking control problem, a reduced order model around the trajectory could be useful. Also, a limit cycle may be important, and analysis or stabilization of a limit cycle is interesting to research. In order to tackle such problems, it is worth developing empirical nonlinear model reduction methods, which are also applicable around a non steady-state. Recently, a connection between POD and nonlinear controllability functions is established by [19] in a stochastic setting. Empirical nonlinear observability Gramians have also been proposed [26, 30]. Nevertheless, neither of these two methods deals with both controllability and observability Gramians, and there is no direct connection between these two works.

In this paper, we propose an empirical balancing method for nonlinear systems with linear time-invariant input vector field components by utilizing its variational system. Since the variational system can be viewed as a linear time-varying system (LTV) along the trajectory of the nonlinear system, one can extend the concept of the controllability and observability Gramians of the LTV system [25, 40]. We call them the differential reachability and observability Gramians, respectively. They depend on the state trajectory of the nonlinear system. In general, it is not easy to obtain them as functions of the trajectory. Nevertheless, we show that their values at each fixed trajectory can be computed from the impulse and initial state responses of the variational system along this fixed trajectory. These obtained trajectory-wise Gramians are constant matrices, and thus one can compute balanced coordinates and a reduced order model in a similar manner as in the LTI case.

The proposed empirical balancing method requires the variational system model. For large-scale systems, computing it may be challenging. Therefore, we further develop approximation methods, which do not require the variational model. Our approach is based on the fact that the variational system is a state space representation of the Fréchet derivative of an operator defined by the nonlinear system, and we use its discretization approximation. For the observability Gramian, similar approximation methods are found in [26, 30]. However, there has been no corresponding controllability Gramian, which has been a bottleneck for developing the corresponding balancing method.

Similar nonlinear balanced realizations are found in flow balancing [38, 39, 37] and in differential balancing [24] but they are not empirical methods and require solutions to nonlinear PDEs. Moreover, [24] does not give the concept of a Gramian. A preliminary version of our work is found in [23]. In this paper, we further develop the discretization approximation methods. Moreover, we newly propose another differential balancing method for a class of nonlinear systems, which only requires the impulse responses of the variational system.

The remainder of this paper is organized as follows. In Section 2, we provide comprehensive back ground of linear balanced truncation in order to help understanding the whole picture of this paper. In Section 3, we define the differential reachability and observability Gramians and then a differentially balanced realization along a trajectory of the system. In Section 4, we show that the value of the differential reachability/observability Gramian can be computed by using the impulse/initial state responses of the variational system. Then, we develop approximation methods, which only require empirical data of the original nonlinear system. In Section 5, we study positive definiteness of the differential reachability and observability Gramians related with nonlinear local strong accessibility and local observability. Next, we propose another differential balancing method, which is further computationally oriented. In Section 6, an example demonstrates our method for an RL network. Finally in Section 6, we conclude the paper by summarizing our results.

2 Review of Linear Empirical Balancing

In this section, we summarize the results for balanced truncation of linear time-invariant (LTI) systems (for more details, see, e.g. [1, 42]) in order to help understanding the whole picture of this paper.

Consider the following SISO LTI system:

[TABLE]

where $x(t)\in{\mathbb{R}}^{n}$ and $u(t),y(t)\in{\mathbb{R}}$; $A\in{\mathbb{R}}^{n\times n}$, $B\in{\mathbb{R}}^{n}$, and $C^{\top}\in{\mathbb{R}}^{n}$. Its general solution is

[TABLE]

Based on the general solution, the controllability and observability Gramians are defined as

[TABLE]

They are positive definite for finite interval $[t_{0},t_{f}]$ for $t_{f}>t_{0}$ if and only if the system is controllable and observable.

Let assume that the system is exponentially stable. When $t_{0}=0$ and $t_{f}\to\infty$, it is known that the eigenvalues of the product $G_{o}(0,\infty)G_{c}(0,\infty)$ correspond to the Hankel singular values of the linear system. Furthermore, there is a change of coordinates $z=Tx$ such that

[TABLE]

In this coordinate, $z_{i}$ are sorted in descending order corresponding to the Hankel singular values $\sigma_{i}$. That is, $z_{i}$ is more important to capture the input-output behavior than $z_{i+1}$ if $\sigma_{i}>\sigma_{i+1}$. In balanced truncation, to approximate the input-output behavior by a reduced order model, the state variables corresponding to small Hankel singular values are truncated.

It is possible to compute the controllability/observability Gramian based on the impulse/initial state responses. From (1), the impulse response of the linear system is $x_{\rm Imp}(t)=Be^{A(t-t_{0})}$. From (2), one notices that

[TABLE]

Next, let $e^{n}_{i}\in{\mathbb{R}}^{n}$ denote the standard basis, i.e., whose $i$th element is $1$, and the other elements are zero, and let $y_{{\rm Is},i}(t)$ denote the corresponding output response. Then, we have

[TABLE]

Moreover, from (3), one notices that

[TABLE]

Therefore, balanced truncation can be achieved based on empirical data. In this paper, we consider to extend these results to nonlinear systems.

3 Differential Balancing along a Trajectory

We present an empirical balancing method for a nonlinear system with LTI input vector field components by using its variational system; the reason considering such a vector field is elaborated in Remark 4.1 below. The proposed empirical balancing method is based on two Gramians, which we call differential reachability and observability Gramians. They can be viewed as extensions of Gramians for linear time-varying (LTV) systems [25, 40] because the variational system can be viewed as an LTV system along a trajectory of the nonlinear system.

3.1 Preliminaries

Consider the following nonlinear system with LTI input vector field components (i.e., the input vector fields are constants):

[TABLE]

where $x(t)\in{\mathbb{R}}^{n}$, $u(t)\in{\mathbb{R}}^{m}$, and $y(t)\in{\mathbb{R}}^{p}$; $f:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ and $h:{\mathbb{R}}^{n}\to{\mathbb{R}}^{p}$ are of class $C^{2}$, and $B\in{\mathbb{R}}^{n\times m}$. Let $\varphi_{t-t_{0}}(x_{0},u)$ denote the state trajectory $x(t)$ of the system $\Sigma$ starting from $x(t_{0})=x_{0}\in{\mathbb{R}}^{n}$ for each choice of $u\in L_{2}^{m}[t_{0},\infty)$. Note that since $f$ is of class $C^{2}$, if $u$ is also of class $C^{2}$, then the solution $\varphi_{t-t_{0}}(x_{0},u)$ is a class $C^{2}$ function of $(t,x_{0})$ as long as it exists. Throughout the paper, we assume that $(\varphi_{t-t_{0}}(x_{0},u),u(t))$ are of class $C^{2}$ in a finite time interval $[t_{0},t_{f}]$.

In our method, we use the prolonged system [7] of the system $\Sigma$, which consists of the original system $\Sigma$ and its variational system $d\Sigma$ along $x(t)=\varphi_{t-t_{0}}(x_{0},u)$,

[TABLE]

where $\delta x(t)\in{\mathbb{R}}^{n}$, $\delta u(t)\in{\mathbb{R}}^{m}$ and $\delta y(t)\in{\mathbb{R}}^{p}$. In the time interval $[t_{0},t_{f}]$, the solution $\delta x(t)$ exists for any bounded input $\delta u(t)$ because the variational system $d\Sigma$ is an LTV system along $\varphi_{t-t_{0}}(x_{0},u)$.

Since the variational system is an LTV system, it is possible to extend the aforementioned linear empirical balancing method to a nonlinear system via the variational system. To this end, we compute the solution $\delta x(t)$ of $d\Sigma$. It follows from the chain rule that

[TABLE]

That is, $\partial\varphi_{t-\tau}(x_{\tau},u)/\partial x_{\tau}$ is the transition matrix of $\partial f(\varphi_{t-\tau})/\partial x$ as an LTV system. From the general solution of an LTV system, the solution $\delta x(t)$ to the variational system $d\Sigma$ starting from $\delta x(t_{0})=\delta x_{0}$ with input $\delta u(t)$ along the trajectory $\varphi_{t-t_{0}}(x_{0},u)$ is obtained as

[TABLE]

For the analysis, furthermore, we use a corresponding output when $\delta u\equiv 0$, namely

[TABLE]

3.2 Differential Balanced Realization

Inspired by results for LTI or LTV systems [25, 40], we define the differential reachability and observability Gramians from the variational systems as follows.

Definition 3.1**.**

For given $x_{0}\in{\mathbb{R}}^{n}$ and $u\in L_{2}^{m}[t_{0},t_{f}]$, the differential reachability Gramian is defined as

[TABLE]

where the arguments of $\varphi_{t-t_{0}}$ are $(x_{0},u)$.

Definition 3.2**.**

For given $x_{0}\in{\mathbb{R}}^{n}$ and $u\in L_{2}^{m}[t_{0},t_{f}]$, the differential observability Gramian is defined as

[TABLE]

where the arguments of $\varphi_{t-t_{0}}$ are $(x_{0},u)$.

Note that in the LTI case, they respectively reduce to the controllability Gramian (2) and observability Gramian (3). These differential Gramians exist in $[t_{0},t_{f}]$, $t_{f}>t_{0}$ from the assumption that the solution $\varphi_{t-t_{0}}(x_{0},u)$ exists and is of class $C^{2}$ in $[t_{0},t_{f}]$.

Remark 3.3**.**

Our differential Gramians can be viewed as extensions of Gramians for LTV systems [25, 40]. By substituting $t=t_{f}+t_{0}-\tau$ into (7), we have

[TABLE]

where $\varphi_{t_{f}-\tau}(x_{f},{\mathcal{F}}_{-}(u))$ is the backward trajectory of the system $\Sigma$ starting from $x(t_{f})=x_{f}$ with the input ${\mathcal{F}}_{-}(u)=u(t_{f}+t_{0}-\tau)\in L_{2}^{m}[t_{0},t_{f}]$. This is an extension of the reachability Gramian for an LTV system in [40] to nonlinear prolonged systems. Similarly, the differential observability Gramian is an extended concept of the observability Gramian for LTV systems. $\lhd$

In a similar manner as a standard procedure, one can define a balanced realization between the differential reachability and observability Gramians. Since these differential Gramians are defined as functions of $\varphi(x_{0},u)$, we define our balanced realization trajectory-wise as follows.

Definition 3.4**.**

Let the differential reachability Gramian $G_{{\mathcal{R}}}(t_{0},t_{f},x_{0},u)\in{\mathbb{R}}^{n\times n}$ and differential observability Gramian $G_{{\mathcal{O}}}(t_{0},t_{f},x_{0},u)\in{\mathbb{R}}^{n\times n}$ at fixed $\varphi_{t-t_{0}}(x_{0},u)$ be positive definite. A realization of the system $\Sigma$ is said to be a differentially balanced realization along $\varphi_{t-t_{0}}(x_{0},u)$ if there exists a constant diagonal matrix

[TABLE]

such that $G_{{\mathcal{R}}}(t_{0},t_{f},x,u)=G_{{\mathcal{O}}}(t_{0},t_{f},x,u)=\Lambda$.

It is possible to show that there always exists a differentially balanced realization along $\varphi_{t-t_{0}}(x_{0},u)$ if the differential Gramians are positive definite. Positive definiteness of them will be discussed in Section 5.1 related with local strong accessibility and local observability of the nonlinear system $\Sigma$.

Theorem 3.5**.**

Suppose that differential Gramians $G_{{\mathcal{R}}}(t_{0},t_{f},x_{0},u)$ and $G_{{\mathcal{O}}}(t_{0},t_{f},x_{0},u)$ at fixed $\varphi_{t-t_{0}}(x_{0},u)$ are positive definite. Then, there exists a non-singular matrix $T_{\varphi}\in{\mathbb{R}}^{n\times n}$ which achieves

[TABLE]

That is, a differentially balanced realization along $\varphi_{t-t_{0}}(x_{0},u)$ is obtained after a coordinate transformation $z=T_{\varphi}x$. $\lhd$

Since $G_{{\mathcal{R}}}(t_{0},t_{f},x_{0},u)$ and $G_{{\mathcal{O}}}(t_{0},t_{f},x_{0},u)$ are constant matrices, it is possible to prove Theorem 3.5 in a similar manner as for the LTI case [1]. As in the LTI case, one can compute a reduced order model by truncating the state variables $z_{k},z_{k+1},\dots,z_{n}$ corresponding to small $\sigma_{k},\sigma_{k+1},\dots,\sigma_{n}$. Clearly, a reduced order model changes for a different trajectory and time interval.

4 Empirical Methods

4.1 Empirical Differential Gramians

In the previous section, we defined a differentially balanced realization along a fixed trajectory $\varphi_{t-t_{0}}(x_{0},u)$. For computing the differential Gramians as functions of $\varphi_{t-t_{0}}(x_{0},u)$, or equivalently $(x_{0},u)$, one needs to solve nonlinear partial differential equations (nPDEs) as for similar nonlinear balancing methods [24, 38, 39, 37] in general. Hereafter, we focus on computing the values of the differential Gramians trajectory-wise.

First, we show that the differential reachability Gramian $G_{{\mathcal{R}}}(t_{0},t_{f},\allowbreak x_{0},u)$ along a fixed trajectory $\varphi_{t-t_{0}}(x_{0},u)$ can be computed by using an impulse response of the variational system $d\Sigma$. Let $\delta_{D}(\cdot)$ be Dirac’s delta function, and let $\delta x_{\rm Imp,i}(t)$ be the impulse response of the variational system $d\Sigma$ along the trajectory $\varphi_{t-t_{0}}(x_{0},u)$ with $\delta u(t)=e^{m}_{i}\delta_{D}(t-t_{0})$, where $e^{m}_{i}\in{\mathbb{R}}^{m}$ is the standard basis. Then, by substituting $\delta x_{0}=0$ and $u(t)=e^{m}_{i}\delta_{D}(t-t_{0})$ into (5), we have

[TABLE]

where $B_{i}$ is the $i$th column vector of $B$. From (7), we obtain

[TABLE]

Therefore, for each $x_{0}\in{\mathbb{R}}^{n}$ and $u\in L_{2}^{m}[t_{0},t_{f}]$, the value of the differential reachability Gramian $G_{{\mathcal{R}}}(t_{0},t_{f},x_{0},u)$ is obtained by using the impulse response of $d\Sigma$.

Remark 4.1**.**

The equality (9) does not hold if $B$ is not constant. Indeed, for the system $\dot{x}=f(x,u)$ and its trajectory $\psi_{t-t_{0}}(x_{0},u)$, the differential reachability Gramian is

[TABLE]

However, the impulse response of the corresponding variational system is

[TABLE]

The reachability Gramian and impulse response do not coincide with each other for non-constant $B$. $\lhd$

Next, we show that the differential observability Gramian $G_{{\mathcal{O}}}(t_{0},t_{f},x_{0},u)$ along a fixed trajectory $\varphi_{t-t_{0}}(x_{0},u)$ can be computed by using initial state responses. By substituting $\delta x_{0}=e^{n}_{i}$ and $\delta u=0$ into (6), one obtains the initial output response of $d\Sigma$ along $\varphi_{t-t_{0}}(x_{0},u)$ as

[TABLE]

From (8), we obtain

[TABLE]

Thus, for each $x_{0}\in{\mathbb{R}}^{n}$ and $u\in L_{2}^{m}[t_{0},t_{f}]$, the value of the differential observability Gramian $G_{{\mathcal{O}}}(t_{0},t_{f},x_{0},u)$ is obtained by using the initial state response of $d\Sigma$.

In summary, the value of the differential reachability/observability Gramian for given $x_{0}$ and $u$ is obtained by computing impulse/initial state responses of a variational system $d\Sigma$ along the trajectory $\varphi_{t-t_{0}}(x_{0},u)$. Therefore, trajectory-wise differential balanced truncation is doable based on empirical data.

4.2 Approximation of the Fréchet Derivative

The empirical approach in the previous subsection requires the variational system model in addition to the original system model. If the original nonlinear systems are large-scale, computing the variational system model may need an effort. Therefore, we present approximation methods not requiring the variational system model.

In order to be self-contained, we first introduce the Fréchet derivative of a nonlinear operator. Consider a nonlinear operator $\Sigma(x_{0},u):{\mathbb{R}}^{n}\times L_{2}[t_{0},t_{f}]\ni(x_{0},u)\mapsto(x_{f},y)\in{\mathbb{R}}^{n}\times L_{2}[t_{0},t_{f}]$ defined by the system $\Sigma$. A linear operator $d\Sigma_{(x_{0},u)}(\delta x_{0},\delta u)$ is said to be the Fréchet derivative if for each $x_{0}\in{\mathbb{R}}^{n}$ and $u\in L_{2}[t_{0},t_{f}]$, the following limit exists

[TABLE]

for all $\delta x_{0}\in{\mathbb{R}}^{n}$ and $\delta u\in L_{2}[t_{0},t_{f}]$. From its definition, the Fréchet derivative of nonlinear operator $\Sigma(x_{0},u)$ is given by the variational system $d\Sigma$. Therefore, by using approximation of the Fréchet derivative, one can approximately compute the impulse or initial state responses of the variational systems. A simple approximation is

[TABLE]

Since the nonlinear operator $\Sigma(x_{0},u)$ is given by the system $\Sigma$, a state space representation of the discretized approximation $d\Sigma_{(x_{0},u)}^{\rm app}(\delta x_{0},\delta u)$ is obtained as follows.

[TABLE]

Therefore, $\delta x(t)$ and $\delta y(t)$ can be approximately computed as $\delta x(t)\simeq(x^{2}(t)-x^{1}(t))/s$ and $\delta y(t)\simeq y_{v}(t)$, where $\delta x_{0}$ and $\delta u$ coincide with the differences of pairs of the initial states $(x^{2}(t_{0})-x^{1}(t_{0}))/s$ and inputs $(u^{2}-u^{1})/s$, respectively.

From the above discussion, an approximation of the impulse response (9) is obtained as

[TABLE]

Similar to the reachability Gramian, we need $m+1$ trajectories of the original nonlinear system in this computation by changing $x^{2}(t)$ depending on the choice of $\delta u$.

Next, an approximation of the initial state response (13) is

[TABLE]

In this computation, we need $n+1$ trajectories of the original nonlinear system.

In summary, the differential reachability and observability Gramians can be approximately computed by using $n+m+1$ trajectories of the original nonlinear system, where $x^{1}(t)$ is same for the approximations of both differential reachability and observability Gramians. Note that even if one does not have an exact model of a real-life system, one only needs the impulse and initial state responses. Therefore, it may be possible to compute an approximation of a differentially balanced realization along $\varphi_{t-t_{0}}(x_{0},u)$ by empirical data.

By applying our empirical methods, a change of coordinates $z=T_{\varphi}x$ for balanced realization is obtained, and $T_{\varphi}$ depends on $\varphi_{t-t_{0}}(x_{0},u)$. Still it is challenging to construct a reduced order model, which gives a good approximation for the whole trajectories because this essentially requires solving nPDEs. A potential solution to this problem is to apply deep learning techniques. After computing $T_{\varphi}$ for different choices of $\varphi_{t-t_{0}}(x_{0},u)$, a function fitting method gives a global nonlinear change of coordinates for model reduction. An obtained reduced order model gives a good approximation at least around $\varphi_{t-t_{0}}(x_{0},u)$ used for the computation of $T_{\varphi}$. We can take arbitrary many trajectories, thus resulting in an approximate global method for model reduction. Another potential solution is to employ a basic idea of proper orthogonal decomposition. First, we compute the summation of differential reachability/observability Gramian, e.g. $(1/r)\sum_{i=1}^{r}G_{{\mathcal{R}}}(t_{0},t_{f},x_{i},u_{i})$ for different choices of $\varphi_{t-t_{0}}(x_{i},u_{i})$, $i=1,\dots,r$. Then, we construct a linear change of coordinates which simultaneously diagonalize them and use this for truncation.

4.3 Literature Review

In literature, there are similar nonlinear balancing methods. We compare our methods with them.

First, another type of differential balancing method [24] employs the following differential controllability and observability functions $L_{{\mathcal{C}}}$ and $L_{{\mathcal{O}}}$.

[TABLE]

where $x(t_{0})=x_{0}\in{\mathbb{R}}^{n}$, $u\in L_{2}^{m}(-\infty,t_{0}]$, $\delta x(t_{0})=\delta x_{0}\in{\mathbb{R}}^{n}$ and $\delta x(-\infty)=0$.

[TABLE]

where $x(t_{0})=x_{0}\in{\mathbb{R}}^{n}$, $\delta x(t_{0})=\delta x_{0}\in{\mathbb{R}}^{n}$, $\delta x(\infty)=0$, $u(t)\equiv 0$, and $\delta u(t)\equiv 0$. Note that the differential controllability function gives the minimum energy to transfer the state of the prolonged system from $\delta x(-\infty)=0$ to $\delta x(t_{0})=\delta x_{0}$ for given $x(t_{0})=x_{0}$ and $u$. Therefore, it depends on $x_{0}$, $u$, and $\delta x_{0}$. A similar discussion holds for the differential observability function.

In fact, by using (6) and (8), the differential observability function and our differential observability Gramian are directly related as

[TABLE]

However, the differential reachability Gramian in (7) and the differential controllability function in (14) are different. This corresponds to the difference between reachability and controllability of LTV systems [40]. Reachability is the property to transfer the state from zero to an arbitrary terminal state, and controllability is the property to transfer the state from an arbitrary initial state to zero; they are the same properties for LTI systems. Based on the controllability Gramian of LTV systems, we define the differential controllability Gramian as

[TABLE]

where the arguments of $\varphi_{t_{0}-\tau}$ are $(x(\tau),u)$. If this differential controllability Gramian $G_{{\mathcal{C}}}(t_{0},x_{0},u)$ exists and is positive definite, the differential controllability function $L_{{\mathcal{C}}}(x_{0},u,\delta x_{0})$ can be described as

[TABLE]

The differential controllability Gramian is defined by using a backward trajectory of the nonlinear system $\Sigma$. In contrast, the differential reachability Gramian is based on a forward trajectory and is computationally oriented.

Relating with differential balancing, flow balancing is proposed by [38, 39, 37]. For flow balancing, the reachability and observability Gramians are defined on different time intervals, and the input is fixed for any initial state. Thus, the Gramians for flow balancing are defined as functions of the initial states. In contrast, our differential Gramians also depend on the input trajectory in addition to the initial state. Moreover, in order to achieve flow balancing, solutions to PDEs are required. Our methods may be applicable to develop empirical methods for flow balancing, which is included in future work.

The papers [27, 12, 6, 15, 11, 42] extend linear empirical balancing methods to nonlinear systems by focusing on a steady-state and attract a lot of research interests as computationally tractable nonlinear model reduction methods. Except [6], these methods can be viewed as our method with an approximation of the Fréchet derivative at a steady-state, and [6] gives an empirical method with differential controllability (not reachability) and observability Gramians. In other words, we provide interpretations of those methods in terms of the variational system and an approximation of the Fréchet derivative. For observability, similar Gramians as ours are found for non-control systems [26] and control systems [30]. However, those papers do not provide the explicit description of the Gramians by using the solution of the original system or an interpretation in terms of the Fréchet derivative and do not establish the corresponding controllability Gramian.

This is the first paper to develop empirical nonlinear balancing methods, which releases the requirement of $\varphi_{t-t_{0}}(x_{0},u)$ being a steady state. This relaxation is beneficial to enlarge the class of applications such as analysis and stabilization of a limit cycle and reducing computational complexity of trajectory tracking controller design for an arbitrary trajectory. Furthermore, as in [27, 12, 6, 15, 11, 42], one may use non-impulse or non-initial state responses for model reduction. These different choices of inputs or initial states enable us to deal with wider classes of model reduction problems such as in [14] although such methods may not be interpreted in terms of Gramians.

5 Further Analysis

In this section, we give some remarks for differential balancing proposed in this paper. First, we study positive definiteness of differential reachability and observability Gramians in terms of nonlinear local strong accessibility and local observability when $u\equiv 0$. Next, we show that for a specific class of systems, one can achieve another empirical differential balancing only by using the impulse responses of the variational system.

5.1 Positive Definiteness of Gramians along Autonomous System

The differentially balanced realization is defined for positive definite differential reachability and observability Gramians. In a specific case when $u\equiv 0$, the positive definiteness implies local accessibility and observability of the original nonlinear system $\Sigma$, and the converse is true for local observability; see e.g. [29] for the definitions of local strong accessibility and local observability.

Theorem 5.1**.**

Let $f(x)$ be of class $C^{\infty}$. Then, the system $\Sigma$ is locally strongly accessible if the differential reachability Gramian $G_{{\mathcal{R}}}(t_{1},t_{2},x_{0},u)$ is positive definite for any $x_{0}$ and $u\equiv 0$ for any subinterval $[t_{1},t_{2}]\subset[t_{0},t_{f}]$.

{pf}

For the sake of the simplicity of the discussion, we consider the single input case. Throughout the proof, we use the fact that the variational system $d\Sigma$ along $\varphi_{t-t_{0}}(x_{0},0)$ is an LTV system. The differential reachability Gramian $G_{{\mathcal{R}}}(t_{1},t_{2},x_{0},0)$ is nothing but the controllability Gramian [41] in the sense of LTV systems. For LTV systems, it has been shown in [34, 41] that the controllability Gramian is positive definite for any subinterval $[t_{1},t_{2}]$ if and only if the LTV system satisfies the Kalman-like controllability rank condition [34]; the discussion until here holds for the multiple-input system. In the single input case, the Kalman-like controllability rank condition [34] becomes as follows for the variational system: there exists $r>0$ such that

[TABLE]

where $ad_{f}^{0}B=f$ and $ad_{f}^{r}B:=[ad_{f}^{r-1}B\ B]:=(\partial B/\partial x)ad_{f}^{r-1}B-(\partial ad_{f}^{r-1}B/\partial x)B$, $i=1,2,\dots$; see [29]. In the multiple-input case, $[ad_{f}^{r-1}B_{i},B_{j}]$, $i,j=1,\dots,m$ are also needed to be taken into account. Condition (16) is a sufficient condition for local strong accessibility with Lie algebra; see e.g. [29]. $\lhd$

Remark 5.2**.**

Theorem 5.1 gives a sufficient condition for local strong accessibility in terms of the differential reachability Gramian. The gap between condition (16) and strong accessibility is that the other Lie brackets that appear in the local strong accessibility rank condition [29], e.g., $[B,[f,B]]$ are missing; for more details see e.g. [29]. To cover such Lie brackets, one needs to compute multiple differential reachability Gramians by changing inputs.

Now, we provide the sketch of the idea of using multiple differential reachability Gramians in the single input case. Consider two differential reachability Gramians $G_{{\mathcal{R}}}(t_{1},t_{2},x_{0},0)$ and $G_{{\mathcal{R}}}(t_{1},t_{2},x_{0},u_{1})$, where $u_{1}(t)=1$ for all $t\geq t_{0}$. From the results on the controllability analysis of LTV systems [34, 41], one can confirm that if there exists a non-zero constant vector $v\in{\mathbb{R}}^{n}$ such that

[TABLE]

for any subinterval $[t_{1},t_{2}]\subset[t_{0},t_{f}]$ $(t_{1}<t_{2})$ then

[TABLE]

To cover all Lie brackets that appear in the local strong accessibility rank condition, one needs to compute a large amount of differential reachability Gramians for different inputs. This could even be an infinite number. $\lhd$

For observability, we have a similar result.

Theorem 5.3**.**

Let $u\equiv 0$ and $\delta u\equiv 0$. Also let $f(x)$ and $h(x)$ be of class $C^{\infty}$. Suppose that the observability codistribution [29] of $\Sigma$ has a constant dimension. Then, the system $\Sigma$ is locally observable if and only if the differential observability Gramian $G_{{\mathcal{O}}}(t_{1},t_{2},x_{0},u)$ is positive definite for any $x_{0}$ and $u\equiv 0$ for any subintervals $[t_{1},t_{2}]\subset[t_{0},t_{f}]$.

{pf}

When the observability codistribution has a constant dimension, a system $\Sigma$ with $u\equiv 0$ is locally observable if and only if the nonlinear observability rank condition holds for all initial states [29]. One can confirm that this nonlinear observability rank condition is nothing but the Kalman-like observability rank condition  [34] for the variational system $d\Sigma$ along the trajectory $\varphi_{t-t_{0}}(x_{0},0)$ as an LTV systems. That is, the system $\Sigma$ is locally observable if and only if its variational system is differentially observable [34, 41] under the assumption for the rank of the observability codistribution. Furthermore, the LTV system is differentially observable if and only if its observability Gramian is positive definite for any subinterval $[t_{1},t_{2}]$ [41], which is nothing but $G_{{\mathcal{O}}}(t_{0},t_{f},x_{0},0)$. ∎

Remark 5.4**.**

In Theorem 5.3, a necessary and sufficient condition for the local observability of an autonomous system is provided, which is a sufficient condition for the local observability of a control system. For conservativeness, a similar discussion as Remark 5.2 holds. That is, to make the method less conservative, one may need to compute a large amount of differential observability Gramian with different inputs $u$. $\lhd$

The paper [30] gives a sufficient condition for local observability for non-zero $u$. As well known for LTV systems, the differential reachability and observability Gramians along $\varphi_{t-t_{0}}(x_{0},0)$ are positive definite if and only if the variational systems along $\varphi_{t-t_{0}}(x_{0},0)$ is completely controllable and observable, respectively. The above theorems connect complete controllability [34] and observability [34] of the variational system and nonlinear local strong accessibility and observability, respectively. In addition, the theorems provide an empirical method for checking nonlinear accessibility and observability because one can compute the differential reachability and observability Gramians along $\varphi_{t-t_{0}}(x_{0},u)$ by using impulse and initial state responses of the variational system, respectively.

5.2 Another Differential Balancing Method for Variationally Symmetric Systems

Balancing methods including ours require two Gramians in general. One is for controllability, and the other is for observability. However, for linear systems, there is a class of systems for which one Gramian is constructed from the other. Such systems are called symmetric [1, 35, 25]. Motivated by the results for symmetric systems, we develop another differential balancing method.

This symmetry concept is extended to nonlinear systems [18] and variational systems [22]. We further extend the latter symmetry concept.

Definition 5.5**.**

The system $\Sigma$ is said to be variationally symmetric if there exists a class $C^{1}$ and non-singular $S:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n\times n}$ such that

[TABLE]

hold.

Even though $B$ is constant, a variationally symmetric system can have a nonlinear output because $S$ is a function. If $S$ is constant, the output function should be linear for a system being variationally symmetric. Variational symmetry implies that after a change of coordinates $\delta z=S(x)\delta x$, the variational system becomes

[TABLE]

In the LTI case, the system (22) is called the dual system of the original system $\Sigma$, and the variational symmetry property is called symmetry. Many physical systems such as mechanical systems and RL circuits have this property; see e.g. [36].

For a symmetric LTI system, the controllability (observability) Gramian of the dual system corresponds to the observability (controllability) Gramian of the dual system [1, 35, 25]. Motivated by this fact, we consider to achieve model reduction based on the differential reachability Gramians of the original system and the system (22). For a variationally symmetric system, these two differential reachability Gramians are connected to each other. A similar relation holds for the differential observability Gramians. We leave to the reader this up.

Theorem 5.6**.**

For the variationally symmetric system $\Sigma$ with respect to $S$, the differential reachability Gramian of the system (22) satisfies

[TABLE]

for any $x_{0}\in{\mathbb{R}}^{n}$ and $u\in L_{2}^{m}[t_{0},t_{f}]$ if it exists, where the arguments of $\varphi_{t-t_{0}}$ are $(x_{0},u)$.

{pf}

Throughout this proof, we omit the arguments of $f$, $h$ and $S$, which are all $\varphi_{t-\tau}(x_{\tau},u)$. By using (4) and (18), compute

[TABLE]

Therefore, $S(\varphi_{t-\tau})(\partial\varphi_{t-\tau}/\partial x_{\tau})S^{-1}$ is the transition matrix of (22).

From (19), the differential reachability Gramian of (22) satisfies

[TABLE]

That completes the proof. ∎

In the linear case, the cross Gramian [25, 15] is useful for analysis of symmetric systems. However, the concept of a cross Gramian is missing in the differential case. Similar to the value of the differential reachability Gramian of the original nonlinear system, that of the system (22) can be computed by using the impulse response of the variational systems $d\Sigma$ for the variationally symmetric system. The applicability of the theory developed in this section is illustrated by an example in Section 6. In fact, the paper [42] gives an efficient empirical method for computing the linear controllability Gramian, but not for observablity Gramian. This method can be extended to the computation of the differential reachability Gramians. Therefore, differential balancing based on the only differential reachability Gramian can be proceed computationally efficiently.

6 Example

As a simple example of variationally symmetric system, we consider the following nonlinear RL circuit with nonlinear resisters.

[TABLE]

Note that due to the existence of $x_{i}^{2}$, we cannot use generalized differential balancing [21, 24] nor generalized incremental balancing [3].

This system is variationally symmetric with respect to the identity matrix because $\partial f(x)/\partial x$ is symmetric. Thus, $G_{{\mathcal{R}}}=G_{{\mathcal{R}}}^{*}$ holds, i.e., we only have to compute the differential reachability Gramian $G_{{\mathcal{R}}}$. We compute its value numerically based on the method in Section 4.2 with $s=0.01$ (we also try our numerical method with $s=1$, and the obtained reduced order model is similar as for $s=0.01$ in this example). For this computation, we need snapshots of the trajectories of the system. Since the system is a single input system, we need two trajectories; one is $x^{1}(t)$ around which a reduced oder model is constructed, and the other is $x^{2}(t)$ needed for the computation of the impulse response. For instance, we choose $x^{1}(t)$ as the trajectory starting from $x(0)=0$ with input $u=\sin(t)+\sin(3t)$. As a numerical computational method of snapshots, we use the forward Euler method with the step size $\Delta t=0.01$. Then, the differential reachability Gramian $G_{{\mathcal{R}}}$ is computed numerically in the time interval $[0,100]$. In the case when $G_{{\mathcal{R}}}=G_{{\mathcal{R}}}^{*}$, empirical differential balanced truncation is achieved by computing the eigenspace of $G_{{\mathcal{R}}}$ and truncating state variables corresponding to small eigenvalues. Figure 1 shows eigenvalues of $G_{{\mathcal{R}}}$. Figure 2 shows the output trajectories of the original system, 10-and 5-dimensional reduced order models. It can be observed that the trajectory of the 10th order model follows the trajectory of the original model really well. The 5th order model is still a good approximation except the time interval $[20,30]$ when the output trajectory changes significantly. Every process here is conduced by using Matlab 2019a on MacBook Pro (2.7 GHz Intel Core i7;16 GB 2133 MHz LPDDR3).

7 Conclusion and Future Work

7.1 Conclusion

In this paper, we have proposed a nonlinear empirical differential balancing method along a fixed state trajectory for nonlinear systems. Our method is based on the differential reachability and observability Gramians, which are functions of the state trajectory. The values of these Gramians at each trajectory is computable by using impulse and initial state responses of the variational system along the trajectory. We have also developed approximation methods for computing them, which only requires empirical data of the original nonlinear systems.

7.2 Possible Application

In [13, 5], empirical balancing at a stead-state is used to reduce the computational complexity of nonlinear model predictive control (MPC) [10, 4]. Our proposed empirical differential balancing method along a fixed state trajectory can be used to reduce the computational complexity around non steady states. In MPC, we repeatedly solve the following nonlinear optimal control problem.

[TABLE]

If the optimal control input $u$ in the time interval $[t_{0},t_{0}+\Delta t]$ is obtained, one can compute the state trajectory of the controlled system in this time interval. Along this trajectory, it is possible to achieve the proposed empirical differential balanced truncation. Then, we have a reduced order model. To compute the optimal control input in the next time interval $[t_{0}+\Delta t,t_{0}+2\Delta t]$, one can use the reduced order model. For this reduced order model, one can compute the corresponding cost function to (23). By solving the reduced order optimal control problem, one has the optimal control input for the reduced order model, which is an approximation of the optimal control input for the original system in the time interval $[t_{0}+\Delta t,t_{0}+2\Delta t]$. Thus, one can use this input for controlling the original system and have the state trajectory of the controlled original system in the time interval $[t_{0}+\Delta t,t_{0}+2\Delta t]$. Then, one can again exploit our empirical model reduction method for obtaining a reduced order model. By repeating this procedure, one can compute an approximation of the optimal control input in each time interval recursively.

Acknowledgment

This work of Y. Kawano was supported by JSPS KAKENHI Grant Number JP19K23517 and JST CREST Grant Number JPMJCR15K2, Japan.