Stability Analysis using Quadratic Constraints for Systems with Neural Network Controllers

Consider the feedback system consisting of a plant $G$ and state-feedback controller $\pi$ as shown in Figure 1. As a first step, we assume the plant $G$ is a linear, time-invariant (LTI) system defined by the following discrete-time model:

$\displaystyle x(k+1)$

$\displaystyle=A_{G}\ x(k)+B_{G}\ u(k),$

(2)

where $x(k)\in\mathbb{R}^{n_{G}}$ is the state, $u(k)\in\mathbb{R}^{n_{u}}$ is the input, $A_{G}\in\mathbb{R}^{n_{G}\times n_{G}}$ and $B_{G}\in\mathbb{R}^{n_{G}\times n_{u}}$. The controller $\pi:\mathbb{R}^{n_{G}}\rightarrow\mathbb{R}^{n_{u}}$ is an $\ell$-layer, feed-forward neural network (NN) defined as:

	$\displaystyle w^{0}(k)$	$\displaystyle=x(k),$		(3a)
	$\displaystyle w^{i}(k)$	$\displaystyle=\phi^{i}\left(\ W^{i}w^{i-1}(k)+b^{i}\ \right),\ i=1,\ldots,\ell,$		(3b)
	$\displaystyle u(k)$	$\displaystyle=W^{\ell+1}w^{\ell}(k)+b^{\ell+1},$		(3c)

where $w^{i}\in\mathbb{R}^{n_{i}}$ are the outputs (activations) from the $i^{th}$ layer and $n_{0}=n_{G}$. The operations for each layer are defined by a weight matrix $W^{i}\in\mathbb{R}^{n_{i}\times n_{i-1}}$, bias vector $b^{i}\in\mathbb{R}^{n_{i}}$, and activation function $\phi^{i}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}^{n_{i}}$. The activation function $\phi^{i}$ is applied element-wise, i.e.

$\displaystyle\phi^{i}(v):=\begin{bmatrix}\varphi(v_{1}),\cdots,\varphi(v_{n_{i}})\end{bmatrix}^{\top},$

(4)

where $\varphi:\mathbb{R}\rightarrow\mathbb{R}$ is the (scalar) activation function selected for the NN. Common choices for the scalar activation function include $\varphi(\nu):=\tanh(\nu)$, sigmoid $\varphi(\nu):=\frac{1}{1+e^{-\nu}}$, ReLU $\varphi(\nu):=\max(0,\nu)$, and leaky ReLU $\varphi(\nu):=\max(a\nu,\nu)$ with $a\in(0,1)$. We assume the activation $\varphi$ is identical in all layers; this can be relaxed with minor changes to the notation.

The state vector $x_{*}$ is an equilibrium point of the feedback system with input $u_{*}$ if the following conditions hold:

	$\displaystyle x_{*}$	$\displaystyle=A_{G}\ x_{*}+B_{G}\ u_{*},$		(5a)
	$\displaystyle u_{*}$	$\displaystyle=\pi(x_{*}).$		(5b)

Let $\chi(k;x_{0})$ denote the solution to the feedback system at time $k$ from initial condition $x(0)=x_{0}$. Our goal is to analyze asymptotic stability of the equilibrium point and to find the largest estimate of the region of attraction, defined below, using an ellipsoidal inner approximation.

It is useful to isolate the nonlinear activation functions from the linear operations of the NN as done in [8, 13]. Define $v^{i}$ as the input to the activation function $\phi^{i}$:

$\displaystyle v^{i}(k):=W^{i}w^{i-1}(k)+b^{i},\ i=1,\ldots,\ell.$

(7)

The nonlinear operation of the $i^{th}$ layer (3b) is thus expressed as $w^{i}(k)=\phi^{i}(v^{i}(k))$. Gather the inputs and outputs of all activation functions:

$\displaystyle v_{\phi}:=\begin{bmatrix}v^{1}\\ \vdots\\ v^{\ell}\end{bmatrix}\in\mathbb{R}^{n_{\phi}}\mbox{ and }w_{\phi}:=\begin{bmatrix}w^{1}\\ \vdots\\ w^{\ell}\end{bmatrix}\in\mathbb{R}^{n_{\phi}},$

(8)

where $n_{\phi}:=n_{1}+\cdots+n_{\ell}$, and define the combined nonlinearity $\phi:\mathbb{R}^{n_{\phi}}\rightarrow\mathbb{R}^{n_{\phi}}$ by stacking the activation functions:

$\displaystyle\phi(v_{\phi}):=\begin{bmatrix}\phi^{1}(v^{1})\\ \vdots\\ \phi^{\ell}(v^{\ell})\end{bmatrix}.$

(9)

Thus $w_{\phi}(k)=\phi(v_{\phi}(k))$, where the scalar activation function $\varphi$ is applied element-wise to each entry of $v_{\phi}$. Finally, the NN control policy $\pi$ defined in (3) can be rewritten as:

	$\displaystyle\begin{bmatrix}u(k)\\ v_{\phi}(k)\end{bmatrix}$	$\displaystyle=N\begin{bmatrix}x(k)\\ w_{\phi}(k)\\ 1\end{bmatrix}$		(10a)
	$\displaystyle w_{\phi}(k)$	$\displaystyle=\phi(v_{\phi}(k)).$		(10b)

The matrix $N$ depends on the weights and biases as follows, where the vertical and horizontal bars partition $N$ compatibly with the inputs $(x,w_{\phi},1)$ and outputs $(u,v_{\phi})$:

	$\displaystyle N$	$\displaystyle:=\left[\begin{array}[]{c|cccc|c}0&0&0&\cdots&W^{\ell+1}&b^{\ell+1}\\ \hline\cr W^{1}&0&\cdots&0&0&b^{1}\\ 0&W^{2}&\cdots&0&0&b^{2}\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\ 0&0&\cdots&W^{\ell}&0&b^{\ell}\end{array}\right]$		(11f)
		$\displaystyle:=\begin{bmatrix}N_{ux}&N_{uw}&N_{ub}\\ N_{vx}&N_{vw}&N_{vb}\end{bmatrix}.$		(11g)

This decomposition of the NN, depicted in Figure 2, isolates the activation functions in preparation for the stability analysis.

Suppose $(x_{*},u_{*})$ satisfies (5). Then $x_{*}$ can be propagated through the NN to obtain equilibrium values $v_{*}^{i}$, $w_{*}^{i}$ for the inputs/outputs of each activation function ($i=1,\ldots,\ell$), yielding $(v_{\phi},w_{\phi})=(v_{*},w_{*})$. Thus $(x_{*},u_{*},v_{*},w_{*})$ is an equilibrium point of (2) and (3) if:

	$\displaystyle x_{*}$	$\displaystyle=A_{G}\ x_{*}+B_{G}\ u_{*},$		(12a)
	$\displaystyle\begin{bmatrix}u_{*}\\ v_{*}\end{bmatrix}$	$\displaystyle=N\begin{bmatrix}x_{*}\\ w_{*}\\ 1\end{bmatrix},$		(12b)
	$\displaystyle w_{*}$	$\displaystyle=\phi(v_{*}).$		(12c)

The stability analysis relies on quadratic constraints (QCs) to bound the activation function. A typical constraint is the sector bound as defined next.

The interpretation of the sector $[\alpha,\beta]$ is that $\varphi$ lies between lines passing through the origin with slope $\alpha$ and $\beta$. Many activation functions are bounded in the sector $[0,1]$, e.g. $\tanh$ and ReLU. Figure 3 illustrates $\varphi(\nu)=\tanh(\nu)$ (blue solid) and the global sector defined by $[0,1]$ (red solid lines).

The global sector constraint is often too coarse for stability analysis, and a local sector constraint provides tighter bounds.

As an example, $\varphi(\nu):=\tanh(\nu)$ restricted to the interval $[-\bar{\nu},\bar{\nu}]$ satisfies the local sector bound $[\alpha,\beta]$ with $\alpha:=\tanh(\bar{\nu})/\bar{\nu}>0$ and $\beta:=1$. As shown in Figure 3 (green dashed lines), the local sector provides a tighter bound than the global sector. These bounds are valid for a symmetric interval around the origin with $\underline{\nu}=-\bar{\nu}$; non-symmetric intervals ($\underline{\nu}\neq-\bar{\nu}$) can be handled similarly.

The local and global sector constraints above were defined to be centered at the point $(\nu,\varphi(\nu))=(0,0)$. The stability analysis will require offset sectors centered around an arbitrary point $(\nu_{*},\varphi(\nu_{*}))$ on the function. For example, $\varphi(\nu)=\tanh(\nu)$ satisfies the global sector bound (red solid) around the point $(\nu_{*},\tanh(\nu_{*}))$ with $[\alpha,\beta]=[0,1]$, as shown in Figure 4. It satisfies a tighter local sector bound (green dashed) when the input is restricted to $\nu\in[\underline{\nu},\bar{\nu}]$. An explicit expression for this local sector is $\beta=1$ and

$\displaystyle\alpha:=\min\left(\frac{\tanh(\bar{\nu})-\tanh(\nu_{*})}{\bar{\nu}-\nu_{*}},\frac{\tanh(\nu_{*})-\tanh(\underline{\nu})}{\nu_{*}-\underline{\nu}}\right).$

The local sector upper bound $\beta$ can be tightened further. This leads to the following definition of an offset local sector.

Offset local sector constraints can also be defined for the combined nonlinearity $\phi$, given by (9). Let $\underline{v},\bar{v},v_{*}\in\mathbb{R}^{n_{\phi}}$ be given with $\underline{v}\leq v_{*}\leq\bar{v}$. Assume that the activation input $v_{\phi}\in\mathbb{R}^{n_{\phi}}$ lies, element-wise, in the interval $[\underline{v},\bar{v}]$ and the $i^{th}$ input/output pair is $w_{\phi,i}=\varphi(v_{\phi,i})$. Further assume the scalar activation function satisfies the local sector $[\alpha_{i},\beta_{i}]$ around the point $v_{*,i}$ with the input restricted to $v_{\phi,i}\in[\underline{v}_{i},\bar{v}_{i}]$ for $i=1,\ldots,n_{\phi}$. The local sector bounds can be computed for $\varphi$ on the given interval either analytically (as above for $\tanh$) or numerically. These local sectors can be stacked into vectors $\alpha_{\phi},\beta_{\phi}\in\mathbb{R}^{n_{\phi}}$ that provide QCs satisfied by the combined nonlinearity $\phi$.

In order to apply the local sector and slope bounds in the stability analysis, we must first compute the bounds $\underline{v},\overline{v}\in\mathbb{R}^{n_{\phi}}$ on the activation input $v_{\phi}$. The process to compute the bounds is briefly discussed here with more details provided in [22]. Let $v_{*}^{1}$ be the equilibrium value at the first NN layer. Select $\underline{v}^{1}$, $\bar{v}^{1}\in\mathbb{R}^{n_{1}}$ with $\underline{v}^{1}\leq v_{*}^{1}\leq\bar{v}^{1}$. The assumed bounds on $v^{1}$ can be used to compute an interval $[\underline{w}^{1},\bar{w}^{1}]$ for the output $w^{1}=\phi^{1}(v^{1})$^†††For example, if $\varphi(\nu)=\tanh(\nu)$ then the input bound $\nu\in[-\bar{\nu},\bar{\nu}]$ implies the output bound $\varphi(\nu)\in[-\tanh(\bar{\nu}),\tanh(\bar{\nu})]$. which can then be used to compute bounds $[\underline{v}^{2},\bar{v}^{2}]$ on the input $v^{2}$ to the next activation function.^‡‡‡The next activation input is $v^{2}:=W^{2}w^{1}+b^{2}$. The largest value of the $i^{th}$ entry of this vector is obtained by solving the following optimization: $\displaystyle\bar{v}^{2}_{i}:=\max_{\underline{w}^{1}\leq w^{1}\leq\bar{w}^{1}}y^{\top}w^{1}+b_{i}^{2}$ (18) where $y^{\top}$ is the $i^{th}$ row of $W^{2}$. Define $c:=\frac{1}{2}(\bar{w}^{1}+\underline{w}^{1})$ and $r:=\frac{1}{2}(\bar{w}^{1}-\underline{w}^{1})$. The optimization can be rewritten as: $\displaystyle\bar{v}^{2}_{i}:=\left(y^{\top}c+b_{i}^{2}\right)+\max_{-r\leq\delta\leq r}y^{\top}\delta$ (19) This has the explicit solution $\bar{v}^{2}_{i}=y^{\top}c+b_{i}^{2}+\sum_{j=1}^{n_{1}}|y_{j}r_{j}|$. Similarly, the minimal value is $\underline{v}^{2}_{i}=y^{\top}c+b_{i}^{2}-\sum_{j=1}^{n_{1}}|y_{j}r_{j}|$. The intervals computed for $w^{1}$ and $v^{2}$ will contain their equilibrium value $w_{*}^{1}$ and $v_{*}^{2}$. This process can be propagated through all layers of the NN to obtain the bounds $\underline{v},\overline{v}\in\mathbb{R}^{n_{\phi}}$ for the activation function input $v_{\phi}$. The remainder of the paper will assume the local sector bounds have been computed as briefly summarized in the following property.

This section uses a Lyapunov function and the offset local sector to compute an inner approximation for the ROA of the feedback system of $G$ and $\pi$. To simplify notation, the interval bound on $v^{1}$ is assumed to be symmetrical about $v^{1}_{*}$, i.e. $\underline{v}^{1}=2v^{1}_{*}-\bar{v}^{1}$ so that $\bar{v}^{1}-v_{*}^{1}=v_{*}^{1}-\underline{v}^{1}$. This can be relaxed to handle non-symmetrical intervals with minor notational changes.

Consider the uncertain feedback system in Figure 5, consisting of an uncertain plant $F_{u}(G,\Delta)$ and a NN controller $\pi$ as defined by (3). The uncertain plant $F_{u}(G,\Delta)$ is an interconnection of a nominal plant $G$ and a perturbation $\Delta$. The nominal plant $G$ is defined by the following equations:

	$\displaystyle x(k+1)$	$\displaystyle=A_{G}\ x(k)+B_{G1}\ q(k)+B_{G2}\ u(k),$		(24a)
	$\displaystyle p(k)$	$\displaystyle=C_{G}\ x(k)+D_{G1}\ q(k)+D_{G2}\ u(k),$		(24b)

where $x(k)\in\mathbb{R}^{n_{G}}$ is the state, $u(k)\in\mathbb{R}^{n_{u}}$ is the control input, $p(k)\in\mathbb{R}^{n_{p}}$ and $q(k)\in\mathbb{R}^{n_{q}}$ are the input and output of $\Delta$, $A_{G}\in\mathbb{R}^{n_{G}\times n_{G}}$, $B_{G1}\in\mathbb{R}^{n_{G}\times n_{q}}$, $B_{G2}\in\mathbb{R}^{n_{G}\times n_{u}}$, $C_{G}\in\mathbb{R}^{n_{p}\times n_{G}}$, $D_{G1}\in\mathbb{R}^{n_{p}\times n_{q}}$, and $D_{G2}\in\mathbb{R}^{n_{p}\times n_{u}}$. The perturbation is a bounded, causal operator $\Delta:\ell^{n_{p}}_{2e}\rightarrow\ell^{n_{q}}_{2e}$. The nominal plant $G$ and perturbation $\Delta$ form the interconnection $F_{u}(G,\Delta)$ through the constraint

$\displaystyle q(\cdot)=\Delta(p(\cdot)).$

(25)

Denote the set of perturbations to be considered as $\mathcal{S}$.

Let $\chi(k;x_{0},\Delta)$ denote the solution to the feedback system of $F_{u}(G,\Delta)$ and $\pi$ with $\Delta\in\mathcal{S}$ at time $k$ from the initial condition $x(0)=x_{0}$^¶¶¶An input/output model is used for the perturbation $\Delta$ so that its internal state and initial condition is not explicitly considered.. Define the robust ROA associated with $x_{*}$ as follows.

The objective is to prove the uncertain feedback system is asymptotically stable and, if so, to find the largest estimate of the robust ROA using an ellipsoidal inner approximation.

The perturbation can represent various types of uncertainty [14], [15], including saturation, time delay, unmodeled dynamics, and slope-restricted nonlinearities. The input-output relationship of $\Delta$ is characterized with an integral quadratic constraint (IQC), which consists of a ‘virtual’ filter $\Psi_{\Delta}$ applied to the input $p$ and output $q$ of $\Delta$ and a constraint on the output $r$ of $\Psi_{\Delta}$. The filter $\Psi_{\Delta}$ is an LTI system of the form:

	$\displaystyle\psi(k+1)$	$\displaystyle=A_{\Psi}\ \psi(k)+B_{\Psi 1}\ p(k)+B_{\Psi 2}\ q(k)$		(27a)
	$\displaystyle r(k)$	$\displaystyle=C_{\Psi}\ \psi(k)+D_{\Psi 1}\ p(k)+D_{\Psi 2}\ q(k)$		(27b)
	$\displaystyle\psi(0)$	$\displaystyle=0$		(27c)

where $\psi(k)\in\mathbb{R}^{n_{\psi}}$ is the state, $r(k)\in\mathbb{R}^{n_{r}}$ is the output, and $A_{\Psi}$ is a Schur matrix. The state matrices have compatible dimensions. The dynamics of $\Psi_{\Delta}$ can be compactly denoted by $\left[\begin{array}[]{c|cc}A_{\Psi}&B_{\Psi 1}&B_{\Psi 2}\\ \hline\cr C_{\Psi}&D_{\Psi 1}&D_{\Psi 2}\end{array}\right]$. By $(p_{*},q_{*})=0$ from Assumption 1, the equilibrium state $\psi_{*}\in\mathbb{R}^{n_{\psi}}$ of (27) is also zero.

The Lyapunov analysis in the next subsection makes use of time-domain IQCs as defined next:

The notation $\Delta\in$ IQC$(\Psi_{\Delta},M_{\Delta})$ indicates that $\Delta$ satisfies the IQC defined by $\Psi_{\Delta}$ and $M_{\Delta}$. Therefore, the precise relation (25), for analysis, is replaced by the constraint (28) on $r$. The QC proposed in Lemma 1 is a special instance of a time-domain IQCs. Specifically, Lemma 1 defines a QC that holds at each time step $k$ and hence the inequality also holds summing over any finite horizons. This is referred to as the offset local sector IQC.

The time-domain IQCs, as defined here, hold on any finite horizon $N\geq 0$. These are typically called “hard IQCs” [14]. IQCs can also be defined in the frequency domain and equivalently expressed as time-domain constraints over an infinite horizon ($N=\infty$). These are called soft IQCs. Although this paper focuses on the use of hard IQCs, it is possible to also incorporate soft IQCs [21, 20, 24, 25].

Let $\zeta:=\left[\begin{smallmatrix}x\\ \psi\end{smallmatrix}\right]\in\mathbb{R}^{n_{\zeta}}$ define the extended state vector, $n_{\zeta}=n_{G}+n_{\psi}$, whose dynamics are

	$\displaystyle\zeta(k+1)$	$\displaystyle=\mathcal{A}\ \zeta(k)+\mathcal{B}\begin{bmatrix}q(k)\\ u(k)\end{bmatrix}$		(29a)
	$\displaystyle r(k)$	$\displaystyle=\mathcal{C}\ \zeta(k)+\mathcal{D}\begin{bmatrix}q(k)\\ u(k)\end{bmatrix}$		(29b)
	$\displaystyle u(k)$	$\displaystyle=\pi(x(k))$		(29c)

where the state-space matrices are

	$\displaystyle\mathcal{A}=\begin{bmatrix}A_{G}&0\\ B_{\Psi 1}C_{G}&A_{\Psi}\end{bmatrix},\ \mathcal{B}=\begin{bmatrix}B_{G1}&B_{G2}\\ B_{\Psi 1}D_{G1}+B_{\Psi 2}&B_{\Psi 1}D_{G2}\end{bmatrix},$
	$\displaystyle\mathcal{C}=\begin{bmatrix}D_{\Psi 1}C_{G}&C_{\Psi}\end{bmatrix},\ \mathcal{D}=\begin{bmatrix}D_{\Psi 1}D_{G1}+D_{\Psi 2}&D_{\Psi 1}D_{G2}\end{bmatrix}.$

Let $\zeta_{*}:=\left[\begin{smallmatrix}x_{*}\\ \psi_{*}\end{smallmatrix}\right]=0$ define the equilibrium point of the extended system (29). Since IQCs implicitly constrain the input $p$ of the extended system (29), the response of the extended system subject to IQCs “covers” the behaviors of the original uncertain feedback system. The following theorem provides a method for inner-approximating the robust ROA by performing analysis on the extended system subject to IQCs.