Synthesizing Robust Domains of Attraction for State-Constrained Perturbed Polynomial Systems

TL;DR

This paper introduces a semi-definite programming approach to compute robust domains of attraction for perturbed polynomial systems with state constraints, ensuring trajectories remain within safe bounds despite uncertainties.

Contribution

It presents a novel SDP-based method that relaxes Zubov's equation to approximate the maximal robust domain of attraction for constrained systems.

Findings

The method guarantees the existence of solutions under certain conditions.

Solutions' sub-level sets inner-approximate the true domain of attraction.

Illustrative examples demonstrate the effectiveness of the approach.

Abstract

In this paper we propose a novel semi-definite programming based method to compute robust domains of attraction for state-constrained perturbed polynomial systems. A robust domain of attraction is a set of states such that every trajectory starting from it will approach an equilibrium while never violating a specified state constraint, regardless of the actual perturbation. The semi-definite program is constructed by relaxing a generalized Zubov's equation. The existence of solutions to the constructed semi-definite program is guaranteed and there exists a sequence of solutions such that their strict one sub-level sets inner-approximate the interior of the maximal robust domain of attraction in measure under appropriate assumptions. Some illustrative examples demonstrate the performance of our method.

Figures

Tables

Table 1. Table 1: Parameters that control the performance of the semi-definite program ( 22 ) on the examples presented in this section. δ , α 𝛿 𝛼 \delta,\alpha and R 𝑅 R are the constant values in ( 22 ), where 𝒳 ∞ = { 𝐱 ∣ q ​ ( 𝐱 ) < α } subscript 𝒳 conditional-set 𝐱 𝑞 𝐱 𝛼 \mathcal{X}_{\infty}=\{\bm{x}\mid q(\bm{x})<\alpha\} with q ​ ( 𝐱 ) = ‖ 𝐱 ‖ 2 𝑞 𝐱 superscript norm 𝐱 2 q(\bm{x})=\|\bm{x}\|^{2} and B ​ ( 𝟎 , R ) = { 𝐱 ∣ R − ‖ 𝐱 ‖ 2 ≥ 0 } 𝐵 0 𝑅 conditional-set 𝐱 𝑅 superscript norm 𝐱 2 0 B(\bm{0},R)=\{\bm{x}\mid R-\|\bm{x}\|^{2}\geq 0\} . k 𝑘 k , d s subscript 𝑑 𝑠 d_{s} and d s ′ subscript 𝑑 superscript 𝑠 ′ d_{s^{\prime}} denote the degree of the polynomials u k , { s 0 , s 1 , s 2 , i , i = 1 , … , m D , s 3 } u_{k},\{s_{0},s_{1},s_{2,i},i=1,\ldots,m_{D},s_{3}\} and { s 4 , j , s 5 , j , s 6 , j , s 7 , j , s 8 , j , s 9 , j , s 10 , l , j , j , l = 1 , … , n 𝒳 } formulae-sequence subscript 𝑠 4 𝑗 subscript 𝑠 5 𝑗 subscript 𝑠 6 𝑗 subscript 𝑠 7 𝑗 subscript 𝑠 8 𝑗 subscript 𝑠 9 𝑗 subscript 𝑠 10 𝑙 𝑗 𝑗 𝑙 1 … subscript 𝑛 𝒳 \{s_{4,j},s_{5,j},s_{6,j},s_{7,j},s_{8,j},s_{9,j},s_{10,l,j},j,l=1,\ldots,n_{\mathcal{X}}\} in ( 22 ), respectively; T 𝑇 T : computation times (seconds).

Ex.	$k$	$\delta$	$\alpha$	$R$	$d_s$	$d_{s^{\prime }}$	$T$
4.1	8	1	${10}^{-4}$	$1.01$	10	8	1.18
4.1	10	1	${10}^{-4}$	1.01	12	10	1.20
4.1	16	1	${10}^{-4}$	1.01	18	16	1.39
4.1	24	1	${10}^{-4}$	1.01	26	24	1.50
4.2	8	1	${10}^{-2}$	1.211	10	8	0.68
4.2	16	1	${10}^{-2}$	1.211	18	16	8.80
4.3	4	1	${10}^{-2}$	1.01	6	4	1.40
4.3	8	1	${10}^{-2}$	1.01	10	8	0.68
4.3	12	1	${10}^{-2}$	1.01	14	12	1.46
4.4	4	1	${10}^{-2}$	1.01	6	4	0.68
4.4	6	1	${10}^{-2}$	1.01	8	6	0.82
4.4	10	1	${10}^{-2}$	1.01	12	10	3.12
4.5	3	1	${10}^{-2}$	1.01	4	2	20.15
4.5	4	1	${10}^{-2}$	1.01	4	4	35.56
4.5	5	1	${10}^{-2}$	1.01	6	4	1257.20

Table 2. Table 2 : Relative error estimations of computed robust domains of attraction to the maximal robust domain of attraction as a function of the approximating polynomial degree for Example 4.1 .

$k$	$8$	$10$	$16$	$24$
error	13.2%	4.48%	3.41%	2.85%

Table 3. Table 3 : Relative volume error estimations of computed robust domains of attraction to the maximal robust domain of attraction as a function of the approximating polynomial degree for Example 4.2 .

$k$	$8$	$16$
error	7.99%	5.84%

Table 4. Table 4 : Relative volume error estimations of computed robust domains of attraction to the maximal robust domain of attraction as a function of the approximating polynomial degree for Example 4.3 .

$k$	$4$	$8$	$12$
error	25.02%	6.47%	1.75%

Table 5. Table 5 : Relative volume error estimations of computed robust domains of attraction to the maximal robust domain of attraction as a function of the approximating polynomial degree for Example 4.4 .

$k$	$4$	$6$	$10$
error	8.88%	6.94%	4.98%

Table 6. Table 6 : Relative volume error estimations of computed robust domains of attraction to the maximal robust domain of attraction on planes x 1 − x 2 subscript 𝑥 1 subscript 𝑥 2 x_{1}-x_{2} with x 3 = x 4 = x 5 = x 6 = x 7 = 0 subscript 𝑥 3 subscript 𝑥 4 subscript 𝑥 5 subscript 𝑥 6 subscript 𝑥 7 0 x_{3}=x_{4}=x_{5}=x_{6}=x_{7}=0 and x 1 − x 7 subscript 𝑥 1 subscript 𝑥 7 x_{1}-x_{7} with x 2 = x 3 = x 4 = x 5 = x 6 = 0 subscript 𝑥 2 subscript 𝑥 3 subscript 𝑥 4 subscript 𝑥 5 subscript 𝑥 6 0 x_{2}=x_{3}=x_{4}=x_{5}=x_{6}=0 as a function of the approximating polynomial degree for Example 4.5 .

$k$	$3$	$4$	$5$
error($x_1-x_2$)	32.17%	16.83%	10.57%
error($x_1-x_7$)	32.35%	20.40%	11.21%

Table 7. Table 7: An illustration for relative volume error estimations affected by δ 𝛿 \delta .

Ex. 4.1		$k$	$8$	$10$	$16$	$24$
	$\delta =3$	error	2.93%	2.11%	1.75%	1.73%
	$\delta =1$	error	$13.2\%$	4.48%	3.41%	2.85%

Table 8. Table 8: An illustration for relative volume error estimations affected by α 𝛼 \alpha .

Ex. 4.1		$k$	$8$	$10$	$16$	$24$
	$\alpha =0.01$	error	5.73%	4.06%	3.35%	2.76%
	$\alpha ={10}^{-4}$	error	13.2%	4.48%	3.41%	2.85%

Table 9. Table 9: An illustration for relative volume error estimations affected by R 𝑅 R .

Ex. 4.1		$k$	$8$	$10$	$16$	$24$
	$R=1.1$	error	14.55%	7.58%	3.75%	2.82%
	$R=1.01$	error	13.2%	4.48%	3.41%	2.85%

Equations

$$\dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t),\mathbf{d}(t)),$$

$$\mathcal{D}=\{\mathbf{d}(\cdot )\mid \mathbf{d}(\cdot ):[0,+\infty )\to Dismeasurable\}.$$

$$\Vert {\mathbf{\varphi }}_{{\mathbf{x}}_0}^{\mathbf{d}}(t)\Vert \le C{e}^{-\sigma t}\Vert {\mathbf{x}}_0\Vert$$

$$\mathcal{X}=\left\{\mathbf{x}\in {\mathbb{R}}^n|\underset{i=1}{\overset{n_{\mathcal{X}}}{\bigwedge }}[h_i(\mathbf{x})<1]\right\},$$

$${\mathcal{D}}_{ad}({\mathbf{x}}_0)=\{\mathbf{d}(\cdot )\in \mathcal{D}\mid {\mathbf{\varphi }}_{{\mathbf{x}}_0}^{\mathbf{d}}(t)\in \mathcal{X}\text{ for }t\in [0,\infty )\}.$$

$$\mathcal{R}:=\{{\mathbf{x}}_0\in {\mathbb{R}}^n\mid {\mathcal{D}}_{ad}({\mathbf{x}}_0)=\mathcal{D}\text{ and }\underset{t\to \infty }{\lim }{\mathbf{\varphi }}_{{\mathbf{x}}_0}^{\mathbf{d}}(t)=0\text{ for }\mathbf{d}(\cdot )\in \mathcal{D}\}.$$

$${\mathcal{D}}_{ad,\varrho }({\mathbf{x}}_0):=\{\mathbf{d}(\cdot )\in \mathcal{D}\mid \mathtt{dist}({\mathbf{\varphi }}_{{\mathbf{x}}_0}^{\mathbf{d}}(t),{\mathcal{X}}^c)>\varrho \text{ for }t\in [0,\infty )\}.$$

$${\mathcal{R}}_0:=\left\{{\mathbf{x}}_0\in {\mathbb{R}}^n|\begin{array}{rl}& \text{there exists }\varrho >0\text{ with }{\mathcal{D}}_{ad,\varrho }({\mathbf{x}}_0)=\mathcal{D}\text{ and}\\ & \text{ there exists a function }\beta (t)\text{ satisfying }\\ & \underset{t\to \infty }{\lim }\beta (t)=0\text{ with }\Vert {\mathbf{\varphi }}_{{\mathbf{x}}_0}^{\mathbf{d}}(t)\Vert \le \beta (t)\text{ for }\\ & t\in [0,\infty )\text{ and }\mathbf{d}(\cdot )\in \mathcal{D}.\end{array}\right\}.$$

$$B(0,R)=\{\mathbf{x}\in {\mathbb{R}}^n\mid h(\mathbf{x})\ge 0\},$$

$$\dot{\mathbf{x}}(t)=\mathbf{F}(\mathbf{x}(t),\mathbf{d}(t)),$$

$$\Vert \mathbf{F}({\mathbf{x}}_1,\mathbf{d})-\mathbf{F}({\mathbf{x}}_2,\mathbf{d})\Vert \le L_{\mathbf{f}}\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert$$

$${\mathcal{R}}_0=\left\{{\mathbf{x}}_0\in {\mathbb{R}}^n|\begin{array}{rl}& \text{there exists }\varrho >0\text{ with }{\mathcal{D}}_{ad,\varrho }({\mathbf{x}}_0)=\mathcal{D}\text{ and }\\ & \text{ there exists a function}\underset{t\to \infty }{\lim }\beta (t)=0\text{ with }\\ & \Vert {\mathbf{\psi }}_{{\mathbf{x}}_0}^{\mathbf{d}}(t)\Vert \le \beta (t)\text{ for }t\in [0,\infty )\text{ and }\mathbf{d}(\cdot )\in \mathcal{D}.\end{array}\right\},$$

$$\begin{split}l[1-h_{j}(\bm{x})]=\left\{\begin{array}[]{lll}1-h_{j}(\bm{x}),&\text{if }1-h_{j}(\bm{x})>0\\ 0,&\text{otherwise.}\\ \end{array}\right.\end{split}$$

$$g(\bm{x})=\left\{\begin{array}[]{lll}&\alpha-\alpha e^{-\frac{\|\bm{x}\|^{s}}{\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})}},\bm{x}\in\mathcal{X}_{\infty}\\ &q(\bm{x}),\bm{x}\in\mathbb{R}^{n}\setminus\mathcal{X}_{\infty}\end{array}\right.,$$

$${\mathcal{X}}_{\infty }=\{\mathbf{x}\in {\mathbb{R}}^n\mid q(\mathbf{x})<\alpha \}$$

$$B(0,\overline{r})\subset {\mathcal{X}}_{\infty }\text{ and }\partial B(0,\overline{r})\cap \partial {\mathcal{X}}_{\infty }=\varnothing$$

$$|g(\bm{x})-g(\bm{y})|\leq K\max\big{\{}\|\bm{x}\|^{s},\|\bm{y}\|^{s},\|\bm{x}\|^{s-2},\|\bm{y}\|^{s-2}\big{\}}\|\bm{x}-\bm{y}\|,\forall\bm{x},\bm{y}\in B(\bm{0},\overline{r}),$$

$$|g(\mathbf{x})-g(\mathbf{y})|\le K^{\prime }\Vert \mathbf{x}-\mathbf{y}\Vert .$$

$$|g(\mathbf{x})-g(\mathbf{y})|\le K^{\prime }\Vert \mathbf{x}-\mathbf{y}\Vert ,\forall \mathbf{y}\in B(\mathbf{x},{\sigma }^{\prime }).$$

$$|g(\mathbf{x})-g(\mathbf{y})|=|q(\mathbf{x})-q(\mathbf{y})|\le K_1\Vert \mathbf{x}-\mathbf{y}\Vert .$$

$$\begin{array}{rl}|g(\mathbf{x})-g(\mathbf{y})|& =|\alpha -\alpha +\alpha e^{-\frac{\Vert \mathbf{y}{\Vert }^s}{\text{dist}(\mathbf{y},{\mathcal{X}}_{\infty }^c)}}|\\ & =|\alpha e^{-\frac{\Vert \mathbf{y}{\Vert }^s}{\text{dist}(\mathbf{y},{\mathcal{X}}_{\infty }^c)}}|\le |\alpha e^{-\frac{\Vert \mathbf{y}{\Vert }^s}{\Vert \mathbf{x}-\mathbf{y}\Vert }}|\le \alpha \frac{\Vert \mathbf{x}-\mathbf{y}\Vert }{\Vert \mathbf{y}{\Vert }^s}.\end{array}$$

$$\begin{array}{rl}& {\int }_0^{\infty }g({\mathbf{\psi }}_{\mathbf{x}}^{\mathbf{d}}(t))dt\\ & ={\int }_0^{T}g({\mathbf{\psi }}_{\mathbf{x}}^{\mathbf{d}}(t))dt+{\int }_T^{\infty }g({\mathbf{\psi }}_{\mathbf{x}}^{\mathbf{d}}(t))dt\\ & ={\int }_0^{T}g({\mathbf{\psi }}_{\mathbf{x}}^{\mathbf{d}}(t))dt+{\int }_0^{\infty }g({\mathbf{\psi }}_{\mathbf{y}}^{\mathbf{d}}(t))dt-{\int }_0^{\infty }g(0)dt\\ & \le MT+{\int }_0^{\infty }{L}_g\Vert {\mathbf{\psi }}_{\mathbf{y}}^{\mathbf{d}}(t)\Vert dt\\ & \le MT+L_{g}C\Vert \mathbf{y}\Vert {\int }_0^{\infty }{e}^{-\sigma t}dt\\ & \le MT+\frac{L_{g}C\overline{r}}{\sigma },\end{array}$$

$$\begin{split}&|g(\bm{x})-g(\bm{y})|\\ &=\alpha|e^{-\frac{\|\bm{y}\|^{s}}{\texttt{dist}(\bm{y},\mathcal{X}_{\infty}^{c})}}-e^{-\frac{\|\bm{x}\|^{s}}{\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})}}|\\ &\leq\alpha\big{|}\frac{\|\bm{y}\|^{s}}{\texttt{dist}(\bm{y},\mathcal{X}_{\infty}^{c})}-\frac{\|\bm{x}\|^{s}}{\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})}\big{|}\\ &=\frac{\alpha}{\texttt{dist}(\bm{y},\mathcal{X}_{\infty}^{c})\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})}\big{|}\|\bm{y}\|^{s}\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})-\|\bm{x}\|^{s}\texttt{dist}(\bm{y},\mathcal{X}_{\infty}^{c})\big{|}\\ &\leq\frac{\alpha}{M_{2}^{2}}\big{|}\|\bm{y}\|^{s}\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})-\|\bm{x}\|^{s}\texttt{dist}(\bm{y},\mathcal{X}_{\infty}^{c})\big{|}\\ &\leq\frac{\alpha}{M_{2}^{2}}\Big{(}\big{|}\|\bm{y}\|^{s}\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})-\|\bm{x}\|^{s}\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})\big{|}\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\big{|}\|\bm{x}\|^{s}\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})-\|\bm{x}\|^{s}\texttt{dist}(\bm{y},\mathcal{X}_{\infty}^{c})\big{|}\Big{)}\\ &\leq\frac{\alpha}{M_{2}^{2}}\Big{(}M_{1}\big{|}\|\bm{x}\|^{s}-\|\bm{y}\|^{s}\big{|}+\|\bm{x}\|^{s}\big{|}\texttt{dist}(\bm{y},\mathcal{X}_{\infty}^{c})-\texttt{dist}(\bm{x},\mathcal{X}_{\infty}^{c})\big{|}\Big{)}\\ &\leq\frac{\alpha}{M_{2}^{2}}\Big{(}M_{1}M^{\prime}\max\{\|\bm{x}\|^{s-2},\|\bm{y}\|^{s-2}\}\big{|}\|\bm{x}\|-\|\bm{y}\|\big{|}+\|\bm{x}\|^{s}L_{\texttt{dist}}\big{|}\|\bm{y}\|-\|\bm{x}\|\big{|}\Big{)}\\ &\leq K\max\{\|\bm{x}\|^{s-2},\|\bm{y}\|^{s-2},\|\bm{x}\|^{s},\|\bm{y}\|^{s}\}\big{|}\|\bm{x}\|-\|\bm{y}\|\big{|},\end{split}$$

$$\begin{split}&\big{|}\|\bm{x}\|^{s}-\|\bm{y}\|^{s}\big{|}\\ &=s\|\bm{\xi}\|^{s-1}\cdot\big{|}\|\bm{x}\|-\|\bm{y}\|\big{|}\\ &=s\|\bm{\xi}\|^{s-2}\cdot\|\bm{\xi}\|\cdot\big{|}\|\bm{x}\|-\|\bm{y}\|\big{|}\\ &\leq M^{\prime}\max\{\|\bm{x}\|^{s-2},\|\bm{y}\|^{s-2}\}\big{|}\|\bm{x}\|-\|\bm{y}\|\big{|},\end{split}$$

$$\begin{split}V(\bm{x}):=\sup_{\bm{d}(\cdot)\in\mathcal{D}}\sup_{t\in[0,\infty)}\big{\{}\int_{0}^{t}g(\bm{\psi}_{\bm{x}}^{\bm{d}}(\tau))d\tau+h^{\prime}(\bm{\psi}_{\bm{x}}^{\bm{d}}(t))\big{\}}\end{split}$$

$$v(\bm{x}):=1-e^{-\delta V(\bm{x})}=\sup_{\bm{d}(\cdot)\in\mathcal{D}}\sup_{t\in[0,\infty)}\big{\{}1-e^{\delta\tilde{V}}\big{\}},$$

$$\begin{split}V(\bm{x})=\sup_{\bm{d}(\cdot)\in\mathcal{D}}\max\Big{\{}G(\bm{x},t,\bm{d}(\cdot))+&V(\bm{\psi}_{\bm{x}}^{\bm{d}}(t)),\sup_{\tau\in[0,t]}\{G(\bm{x},\tau,\bm{d}(\cdot))+h^{\prime}(\bm{\psi}_{\bm{x}}^{\bm{d}}(\tau))\}\Big{\}}.\end{split}$$

$$\begin{split}v(\bm{x})=\sup_{\bm{d}(\cdot)\in\mathcal{D}}\max\Big{\{}1+&(v(\bm{\psi}_{\bm{x}}^{\bm{d}}(t))-1)e^{-\delta G(\bm{x},t,\bm{d}(\cdot))},\\ &\sup_{\tau\in[0,t]}\{1-e^{-\delta G(\bm{x},\tau,\bm{d}(\cdot))-\delta h^{\prime}(\bm{\psi}_{\bm{x}}^{\bm{d}}(\tau))}\}\Big{\}}.\end{split}$$

$$\begin{split}&|V(\bm{x}_{0})-V(\bm{y})|\\ &\leq\sup_{\bm{d}(\cdot)\in\mathcal{D}}\sup_{t\in[0,\infty)}\Big{(}\int_{0}^{t}|g(\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(\tau))-g(\bm{\psi}_{\bm{y}}^{\bm{d}}(\tau))|d\tau+|h^{\prime}(\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(t))-h^{\prime}(\bm{\psi}_{\bm{y}}^{\bm{d}}(t))|\Big{)}\\ &\leq\sup_{\bm{d}(\cdot)\in\mathcal{D}}\int_{0}^{\infty}|g(\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(\tau))-g(\bm{\psi}_{\bm{y}}^{\bm{d}}(\tau))|d\tau+\sup_{\bm{d}(\cdot)\in\mathcal{D}}\sup_{t\in[0,\infty)}|h^{\prime}(\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(t))-h^{\prime}(\bm{\psi}_{\bm{y}}^{\bm{d}}(t))|.\end{split}$$

$$|g({\mathbf{x}}^{\prime })-g({\mathbf{y}}^{\prime })|\le K\max \{\Vert {\mathbf{x}}^{\prime }{\Vert }^s,\Vert {\mathbf{y}}^{\prime }{\Vert }^s,\Vert {\mathbf{x}}^{\prime }{\Vert }^{s-2},\Vert {\mathbf{y}}^{\prime }{\Vert }^{s-2}\}\Vert {\mathbf{x}}^{\prime }-{\mathbf{y}}^{\prime }\Vert ,\forall {\mathbf{x}}^{\prime },{\mathbf{y}}^{\prime }\in B(0,\overline{r}).$$

Full text

1 Introduction

A robust domain of attraction of interest in this paper is a set of states from which the system will finally approach an equilibrium while never breaching a specified state constraint regardless of the actual perturbation. Computing it is a fundamental task in the analysis of dynamical systems such as power systems [1] and turbulence phenomena in fluid dynamics [5].

Existing approaches to approximating robust domains of attraction can be classified into non-Lyapunov and Lyapunov based categories. Non-Lyapunov based approaches include, but are not limited to, trajectory reversing methods [15], polynomial level-set methods [42] and reachable set computation based methods [23, 20, 13]. Contrasting with non-Lyapunov based methods, Lyapunov based methods are still dominant in estimating robust domains of attraction, e.g., [41, 22, 9, 34, 16]. Such methods are based on the search of a Lyapunov function $V(\bm{x})$ and a positive scalar $b$ such that the first-order Lie derivative of $V(\bm{x})$ is negative over the sub-level set $C=\{\bm{x}\mid V(\bm{x})\leq b\}$. Given such $V(\bm{x})$ and $b$, it can be shown that the connected component of $C$ containing the equilibrium is a robust domain of attraction. Generally, the search for Lyapunov functions is non-trivial for nonlinear systems due to the non-constructive nature of the Lyapunov theory, apart from some cases where the Jacobian matrix of the linearized system associated with the nonlinear system of interest is Hurwitz. However, with the advance of real algebraic geometry [31, 6] and polynomial optimization [33, 11, 19] in the last decades, especially the sum-of-squares (SOS) decomposition technique, finding a Lyapunov function which is decreasing over a given state constraint set can be reduced to a convex programming problem for polynomial systems. This results in a large amount of findings which adopt convex optimization based approaches to the search for polynomial Lyapunov functions, e.g., [30, 10, 4] and the references therein. However, if we return to the problem of estimating robust domains of attraction, it resorts to addressing a bilinear semi-definite program, e.g., [21, 37, 39, 16], which falls within the non-convex programming framework and is notoriously hard to solve [7]. Moreover, the existence of solutions to (bilinear) semi-definite programs is not explored theoretically in the existing literature.

Another way to compute Lyapunov functions and estimate robust domains of attraction is based on solving the Zubov’s equation [45], which is a Hamilton-Jacobi type partial differential equation. Zubov’s equation was originally inferred to describe the maximal domain of attraction for nonlinear systems free of perturbation inputs and state constraints. Recently, it was extended to perturbed nonlinear systems in [8] and further to state-constrained perturbed nonlinear systems in [18]. The appealing aspect of Zubov’s method in [18] is that it touches upon the problem of computing the maximal robust domain of attraction, whose interior is described exactly via the strict one sub-level set of the unique viscosity solution to a generalized Zubov’s equation. However, it is well-known that it is notoriously hard to solve the Zubov type equation generally [45, 18].

In this paper we propose a novel semi-definite programming based method to compute robust domains of attraction for state-constrained perturbed polynomial systems. We first customize the generalized Zubov’s equation in [18] to characterize the maximal robust domain of attraction based on Kirszbraun’s extension theorem for Lipschitz maps. Then we relax the customized Zubov’s equation into a system of inequalities and further encode these inequalities in the form of sum-of-squares constraints such that a robust domain of attraction can be generated via solving a semi-definite program. The semi-definite program falls within the convex programming framework and can be solved efficiently in polynomial time via interior-point methods, consequently providing a practical method for computing robust domains of attraction. Under appropriate assumptions, the existence of solutions to the constructed semi-definite program is guaranteed and there exists a sequence of solutions such that their strict one sub-level sets inner-approximate the interior of the maximal robust domain of attraction in measure. Some illustrative examples demonstrate the performance of our method.

The main contributions of this paper are summarized as follows.

A novel semi-definite programming based method is proposed to synthesize robust domains of attraction. The semi-definite program is constructed based on a customized Zubov’s equation of the form in [18]. Unlike [18] which reduces the problem of computing robust domains of attraction to a problem of solving a Zubov type equation, our approach reduces the robust domains of attraction generation problem to a semi-definite programming problem. 2. 2.

Under appropriate assumptions, the existence of solutions to the constructed semi-definite program is guaranteed and there exists a sequence of solutions such that their strict one sub-level sets inner-approximate the interior of the maximal robust domain of attraction in measure.

This paper is structured as follows. In Section 2 we introduce basic notations used throughout this paper and the robust domains of attraction generation problem of interest. In Section 3 we detail our semi-definite programming based method for computing robust domains of attraction for state-constrained perturbed polynomial systems. After evaluating this approach on five illustrative examples in Section 4, we conclude our paper in Section 5.

2 Preliminaries

In this section we first formulate the problem of generating robust domains of attraction in Subsection 2.1, and then introduce the concept of the maximal robust domain of uniform attraction in Subsection 2.2. The maximal robust domain of uniform attraction, which is the interior of the maximal robust domain of attraction, plays an important role for synthesizing robust domains of attraction.

The following basic notations will be used throughout the rest of this paper: $\mathbb{N}$ denotes the set of non-negative integers. $\mathbb{R}^{n}$ denotes the set of $n$-dimensional real vectors. $\mathbb{R}[\cdot]$ denotes the ring of polynomials in variables given by the argument. $\mathbb{R}_{k}[\cdot]$ denotes the set of real polynomials of degree at most $k$ in variables given by the argument, $k\in\mathbb{N}$. $C^{1}(\mathbb{R}^{n})$ denotes the set of continuously differentiable functions over $\mathbb{R}^{n}$. $\|\bm{x}\|$ denotes the 2-norm, i.e., $\|\bm{x}\|=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}$, where $\bm{x}=(x_{1},\ldots,x_{n})^{\top}$. Vectors are denoted by boldface letters. $\Delta^{\circ}$, $\Delta^{c}$ and $\overline{\Delta}$ denote the interior, the complement and the closure of a set $\Delta$, respectively. $B(\bm{x},r)$ denotes the closed ball around $\bm{x}$ with radius $r>0$ in $\mathbb{R}^{n}$.

2.1 Robust Domains of Attraction

A state-constrained perturbed dynamical system of interest in this paper is of the following form:

[TABLE]

where $\bm{x}(\cdot):[0,\infty)\rightarrow\mathcal{X}$, $\bm{d}(\cdot):[0,\infty)\rightarrow D$, $\mathcal{X}\subset\mathbb{R}^{n}$ is a bounded open set, $D=\{\bm{d}\in\mathbb{R}^{m}\mid\bigwedge_{i=1}^{m_{D}}[h_{i}^{D}(\bm{d})-1\leq 0]\}$ is a compact set in $\mathbb{R}^{m}$ with $h_{i}^{D}(\bm{d})\in\mathbb{R}[\bm{d}]$, and $\bm{f}(\bm{x},\bm{d})\in\mathbb{R}[\bm{x},\bm{d}]$, thus satisfying the local Lipschitz condition.

Denote the set of admissible perturbation inputs as

[TABLE]

As a consequence, for $\bm{x}_{0}\in\mathcal{X}$ and $\bm{d}(\cdot)\in\mathcal{D}$, there exists a unique absolutely continuous trajectory $\bm{\phi}_{\bm{x}_{0}}^{\bm{d}}(t)$ satisfying (1) a.e. with $\bm{\phi}_{\bm{x}_{0}}^{\bm{d}}(0)=\bm{x}_{0}$ for some time interval $[0,T]$ with $T>0$ [38].

Additionally, we have Assumption 2.1 for system (1) throughout this paper.{assumption}

1. (1)

$\bm{f}(\bm{0},\bm{d})=\bm{0},\forall\bm{d}\in D$, i.e., the fixed point $\bm{x}=\bm{0}$ is invariant under all perturbations. 2. (2)

there exist positive constants $C$, $\sigma$, $\overline{r}$ such that

[TABLE]

for $\bm{x}_{0}\in B(\bm{0},\overline{r})$ and $\bm{d}(\cdot)\in\mathcal{D}$, i.e., the equilibrium state $\bm{0}$ is uniformly locally exponentially stable for system (1). 3. (3)

The equilibrium $\bm{x}=\bm{0}$ resides in the interior of the state constraint set $\mathcal{X}$, i.e., there exists $\overline{r}>0$ such that $B(\bm{0},\overline{r})\subset\mathcal{X}$. We may without loss of generality assume that this $\overline{r}$ is the same as that in (2). In addition, we also assume that $\overline{r}$ is sufficiently small such that every trajectory starting within $B(\bm{0},\overline{r})$ will never leave the set $\mathcal{X}$. 4. (4)

The set $\mathcal{X}$ is of the following form

[TABLE]

where $h_{i}(\bm{x})\in\mathbb{R}[\bm{x}]$ with $h_{i}(\bm{x})\geq 0$ over $\mathbb{R}^{n}$ and $h_{i}(\bm{0})=0$. Note that this implies that $\partial\mathcal{X}\subseteq\cup_{i=1}^{n_{\mathcal{X}}}\{\bm{x}\in\mathcal{X}\mid h_{i}(\bm{x})=1\}$.

Remark 2.1**.**

*The assumption that the fixed point $\bm{x}=\bm{0}$ is invariant under all perturbations is relatively conservative. One of application scenarios of this assumption is that the steady state of the real system is fixed however, an uncertain (or, perturbed) model rather than a deterministic model was established to capture behaviors of the system due to incomplete information, e.g., [8, 39, 10]. *

Systems with locally exponentially stable equilibria are widely studied in the existing literature, e.g., [22]. Since it is not enough to know that the system will converge to an equilibrium eventually in many applications, there is a need to estimate how fast the system approaches $\bm{0}$. The concept of exponential stability can be used for this purpose [36].

The goal of this paper is to synthesize robust domains of attraction of the origin for system (1). The maximal robust domain of attraction, which is the set of states such that every possible trajectory starting from it will approach the origin while never leaving the state constraint set $\mathcal{X}$, is formally formulated in Definition 2.2.

Definition 2.2** ((Maximal) Robust Domain of Attraction).**

Denote

[TABLE]

The maximal robust domain of attraction $\mathcal{R}$ is defined as

[TABLE]

*A robust domain of attraction $\Omega$ is a subset of the maximal robust domain of attraction $\mathcal{R}$, i.e., $\Omega\subseteq\mathcal{R}$. *

2.2 Robust Domains of Uniform Attraction

In order to relate robust domains of attraction to a Zubov type equation, a uniform version of the maximal robust domain of attraction is presented in [18]. In this subsection we introduce the maximal robust domain of uniform attraction.

To this end, we define the distance between a point $\bm{x}\in\mathbb{R}^{n}$ and a set $A\subset\mathbb{R}^{n}$ as $\mathtt{dist}(\bm{x},A):=\inf_{\bm{y}\in A}\|\bm{x}-\bm{y}\|$. Then, for $\varrho\geq 0$, we define the set of $\varrho-$admissible perturbation inputs as

[TABLE]

Note that $\mathcal{D}_{ad,0}(\bm{x}_{0})=\mathcal{D}_{ad}(\bm{x}_{0})$. The maximal robust domain of uniform attraction is then defined by

[TABLE]

The set $\mathcal{R}_{0}$ is a uniform version of $\mathcal{R}$ in the sense that for every initial state $\bm{x}_{0}\in\mathcal{R}_{0}$ the trajectories have a positive distance of at least $\varrho$ to $\mathcal{X}^{c}$ and converge towards $\bm{0}$ with a speed characterized by $\beta(t)$. Neither $\varrho$ or $\beta(t)$ depends on $\bm{d}(\cdot)\in\mathcal{D}$. The relation between $\mathcal{R}_{0}$ and $\mathcal{R}$ is uncovered in Lemma 2.3.

Lemma 2.3**.**

[18*]**

The following two statements hold:*

$\mathcal{R}_{0}$* is open.* 2. 2.

$\mathcal{R}_{0}=\mathcal{R}^{\circ}$.

$\mathcal{R}_{0}=\mathcal{R}^{\circ}$ implies that $\mathcal{R}_{0}$ and $\mathcal{R}$ coincide except for a set with void interior.

3 Computation of Robust Domains of Attraction

In this section we detail our method for synthesizing robust domains of attraction. Subsection 3.1 presents an auxiliary system, to which the global solution over $t\in[0,\infty)$ exists for every $\bm{x}\in\mathbb{R}^{n}$, and then Subsection 3.2 presents a customized Zubov’s equation for state-constrained perturbed polynomial systems as well as a semi-definite programming based method for synthesizing robust domains of attraction. Finally, we show in Subsection 3.3 that the constructed semi-definite program is able to generate a convergent sequence of robust domains of attraction to the maximal robust domain of uniform attraction in measure under appropriate assumptions.

3.1 System Reformulation

The system reformulation part in this subsection is similar to that in [43, 44]. Different from the present work, the problems of computing inner-approximations of backward reachable sets over finite time horizons and robust invariant sets over the infinite time horizon are respectively considered in [43, 44] based on relaxing Hamilton-Jacobi type partial differential equations. For self-containedness and ease of understanding, we also give it an appropriate description in this section. As $\bm{f}\in\mathbb{R}[\bm{x},\bm{d}]$ in system (1), $\bm{f}$ is only locally Lipschitz continuous over $\bm{x}$. Therefore, the existence of a global solution $\bm{\phi}_{\bm{x}_{0}}^{\bm{d}}(t)$ over $t\in[0,\infty)$ to system (1) is not guaranteed for any initial state $\bm{x}_{0}\in\mathbb{R}^{n}$ and any perturbation input $\bm{d}(\cdot)\in\mathcal{D}$ [12]. However, the existence of global solutions is a prerequisite for constructing the Zubov’s equation, to which the strict one sub-level set of the viscosity solution characterizes the maximal robust domain of uniform attraction. In this subsection we construct an auxiliary system, to which the global solution over $t\in[0,\infty)$ with any initial state $\bm{x}_{0}\in\mathbb{R}^{n}$ and any perturbation input $\bm{d}(\cdot)\in\mathcal{D}$ exists. Also, its solution coincides with the solution to system (1) over a compact set

[TABLE]

where $h(\bm{x})=R-\|\bm{x}\|^{2}$. The compact set $B(\bm{0},R)$ is chosen to satisfy $\mathcal{X}\subset B(\bm{0},R)$ and $\partial\mathcal{X}\cap\partial B(\bm{0},R)=\emptyset$. Such $R$ in (3) exists since $\mathcal{X}$ is a bounded set in $\mathbb{R}^{n}$. The set $B(\bm{0},R)$ in (3) plays three important roles in our semi-definite programming based approach, which will be shown in Subsection 3.2.

The condition $\mathcal{X}\subseteq B(\bm{0},R)$ guarantees that the maximal robust domain of uniform attraction $\mathcal{R}_{0}$ for system (1) can be exactly characterized by trajectories of the auxiliary system (4), as formulated in Proposition 3.2. 2. 2.

The condition $\partial\mathcal{X}\cap\partial B(\bm{0},R)=\emptyset$ assures that the strict one sub-level set of the approximating polynomial returned by solving the semi-definite program (22) in Subsection 3.2 is a robust domain of attraction. It is useful in justifying Theorem 3.12 in Subsection 3.2. 3. 3.

The condition that $h(\bm{x})$ is of the form $R-\|\bm{x}\|^{2}$, i.e., $h(\bm{x})=R-\|\bm{x}\|^{2}$, is used to guarantee the existence of solutions to the semi-definite program (22) in Subsection 3.2 under appropriate conditions. It is reflected in justifying Theorem 3.16 in Subsection 3.3.

The auxiliary system is of the following form:

[TABLE]

where $\bm{F}(\bm{x},\bm{d}):\mathbb{R}^{n}\times D\rightarrow\mathbb{R}^{n}$, which is globally Lipschitz continuous over $\bm{x}\in\mathbb{R}^{n}$ uniformly over $\bm{d}\in D$ with Lipschitz constant $L_{\bm{f}}$, i.e.,

[TABLE]

for $\bm{x}_{1},\bm{x}_{2}\in\mathbb{R}^{n}$ and $\bm{d}\in D$, where $L_{\bm{f}}$ is the Lipschitz constant of $\bm{f}$ over $B(\bm{0},R)$. Moreover, $\bm{F}(\bm{x},\bm{d})=\bm{f}(\bm{x},\bm{d})$ over $B(\bm{0},R)\times D$, implying that the trajectories governed by system (4) coincide with the ones generated by system (1) over the set $B(\bm{0},R)$.

The existence of system (4) is guaranteed by Kirszbraun’s theorem [14], which is recalled as Theorem 3.1.

Theorem 3.1** (Kirszbraun’s Theorem).**

*Let $A\subset\mathbb{R}^{k}$ be a set and $\bm{f}^{\prime}:A\rightarrow\mathbb{R}^{n}$ a function, where $k\geq 1$ is an integer. Suppose there exists $\gamma\geq 0$ such that $\|\bm{f}^{\prime}(\bm{z}_{1})-\bm{f}^{\prime}(\bm{z}_{2})\|\leq\gamma\|\bm{z}_{1}-\bm{z}_{2}\|$ for $\bm{z}_{1},\bm{z}_{2}\in A$. Then there is a function $\bm{F}^{\prime}:\mathbb{R}^{k}\rightarrow\mathbb{R}^{n}$ such that $\bm{F}^{\prime}(\bm{z})=\bm{f}^{\prime}(\bm{z})$ for $\bm{z}\in A$ and $\|\bm{F}^{\prime}(\bm{z}_{1})-\bm{F}^{\prime}(\bm{z}_{2})\|\leq\gamma\|\bm{z}_{1}-\bm{z}_{2}\|$ for $\bm{z}_{1},\bm{z}_{2}\in\mathbb{R}^{k}$. *

For instance, \bm{F}(\bm{x},\bm{d})=\inf_{\bm{y}\in B(\bm{0},R)}\big{(}\bm{f}(\bm{y},\bm{d})+\bm{z}L_{\bm{f}}\cdot\|\bm{x}-\bm{y}\|\big{)} satisfies (4), where $\bm{z}$ is an $n$-dimensional vector with each component equal to one.

Thus, for any pair $(\bm{d}(\cdot),\bm{x}_{0})\in\mathcal{D}\times\mathbb{R}^{n}$, there exists a unique absolutely continuous trajectory $\bm{x}(t)=\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(t)$ satisfying (4) a.e. with $\bm{x}(0)=\bm{x}_{0}$ for $t\in[0,\infty)$. This requirement is the basis of deriving the Zubov’s equation in [18]. Moreover, we have the following proposition stating that the sets $\mathcal{R}$ and $\mathcal{R}_{0}$ for system (1) coincide with the corresponding sets for system (4) as well.

Proposition 3.2**.**

$\mathcal{R}=\{\bm{x}_{0}\in\mathbb{R}^{n}\mid\mathcal{D}_{ad}(\bm{x}_{0})=\mathcal{D}\text{ and }\lim_{t\rightarrow\infty}\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(t)=0\text{ for }\bm{d}(\cdot)\\ \in\mathcal{D}\}$* and*

[TABLE]

*where $\mathcal{R}$ and $\mathcal{R}_{0}$ are respectively the maximal robust domain of attraction and the maximal robust domain of uniform attraction for system (1). *

Proof 3.3**.**

*Since $\mathcal{X}\subset B(\bm{0},R)$, $\bm{f}(\bm{x},\bm{d})=\bm{F}(\bm{x},\bm{d})$ over $\bm{x}\in\mathcal{X}$ and $\bm{d}\in D$, the trajectories for system (1) and (4) coincide in the state constraint set $\mathcal{X}$, it is obvious that the conclusion holds. *

3.2 Synthesizing Robust Domains of Attraction

In this subsection we first follow the procedure in [18] to characterize the maximal robust domain of uniform attraction $\mathcal{R}_{0}$ of system (4) as the strict one sub-level set of the viscosity solution to a generalized Zubov’s partial differential equation, which is specific to state-constrained perturbed polynomial systems. Based on this Zubov type equation, we construct a semi-definite program for generating robust domains of attraction.

For showing the generalized Zubov’s equation we first introduce a running cost $g(\bm{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}$ and a function $h^{\prime}(\bm{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}$ satisfying Assumption 3.2 as in [18]. {assumption}

The function $g(\bm{x})$ is a locally Lipschitz continuous function over $\bm{x}\in\mathbb{R}^{n}$ satisfying

1. (a)

$g(\bm{x})\geq 0$ with $g(\bm{0})=0$, $\forall\bm{x}\in\mathbb{R}^{n}$; 2. (b)

$\inf\{g(\bm{x})\mid\|\bm{x}\|\geq c\}>0$ for every $c>0$; 3. (c)

$\int_{0}^{\infty}g(\bm{\psi}_{\bm{x}}^{\bm{d}}(t))dt$ is finite if $t(\bm{x},\bm{d}(\cdot))$ is finite, where $t(\bm{x},\bm{d}(\cdot))=\inf\{t\geq 0\mid\bm{\psi}_{\bm{x}}^{\bm{d}}(t)\in B(\bm{0},\overline{r})\}$. 2. 2.

$h^{\prime}(\bm{x})=-\min_{j\in\{1,\ldots,n_{\mathcal{X}}\}}h^{\prime}_{j}(\bm{x})$ with $h^{\prime}_{j}(\bm{x})=\ln(l[1-h_{j}(\bm{x})])$ and

[TABLE]

with the convention $\ln 0=-\infty$.

The function $h^{\prime}(\bm{x})$ fulfills the requirement in (A4) in [18], i.e., $h^{\prime}(\bm{x})$ is locally Lipschitz continuous on $\mathcal{X}$, $h^{\prime}(\bm{x})=\infty$ iff $\bm{x}\notin\mathcal{X}$, and $\lim_{n\rightarrow\infty}h^{\prime}(\bm{x}_{n})=\infty$ when $\lim_{n\rightarrow\infty}\bm{x}_{n}=\bm{x}\notin\mathcal{X}$, $h^{\prime}(\bm{0})=0$. Throughout this paper, the function $g(\bm{x})$ in Assumption 3.2 is chosen as

[TABLE]

where $s$ is a sufficiently large positive scalar value, $\alpha>0$, and $q(\bm{x})\in\mathbb{R}[\bm{x}]$ with $q(\bm{x})>0$ for $\bm{x}\neq\bm{0}$ is a nonnegative polynomial such that

[TABLE]

is a nonempty robust domain of uniform attraction. The function $q(\bm{x})$ could be a (local) Lyapunov function and there are numerous methods for computing it, e.g., semi-definite programming based methods [30, 35]. In this paper we assume that $q(\bm{x})$ is given. Also, without loss of generality we assume that

[TABLE]

since $\overline{r}$ in Assumption 2.1 can be sufficiently small. The function $g(\bm{x})$ in (7) satisfies Assumption 3.2 as well as another important property shown in Lemma 3.4, which will lift the value functions in (11) and (12) shown later to be Lipschitz continuous.

Lemma 3.4**.**

The function $g(\bm{x})$ in (7) satisfies Assumption 3.2 and the following inequality,

[TABLE]

*where $B(\bm{0},\overline{r})$ is defined in Assumption 2.1, $K$ is some positive constant and $s$ is a sufficiently large positive scalar value. *

Proof 3.5**.**

We first prove that the function $g(\bm{x})$ in (7) is locally Lipschitz continuous over $\bm{x}\in\mathbb{R}^{n}$. It is obvious that for $\bm{x},\bm{y}\in\mathcal{X}_{\infty}^{\circ}$ or $\bm{x},\bm{y}\in\mathbb{R}^{n}\setminus\mathcal{X}_{\infty}$, there exists a constant $K^{\prime}$ such that

[TABLE]

In the following we just need to show that for $\bm{x}\in\partial\mathcal{X}_{\infty}$ (Since $\mathcal{X}_{\infty}$ is open, $\bm{x}\in\mathbb{R}^{n}\setminus\mathcal{X}_{\infty}$), there exist a neighborhood $B(\bm{x},\sigma^{\prime})$ of $\bm{x}$ and a constant $K^{\prime}>0$ such that

[TABLE]

Since $\bm{x}\in\partial\mathcal{X}_{\infty}$, $g(\bm{x})=q(\bm{x})=\alpha$. When $\bm{y}\in(\mathbb{R}^{n}\setminus\mathcal{X}_{\infty})\cap B(\bm{x},\sigma^{\prime})$, there exists a constant $K_{1}>0$,

[TABLE]

When $\bm{y}\in\mathcal{X}_{\infty}\cap B(\bm{x},\sigma^{\prime})$, we have

[TABLE]

The last inequality in (10) uses the fact that $e^{-z}\leq\frac{1}{z}$ for $z\geq 0$. Therefore, there exist a neighborhood $B(\bm{x},\sigma^{\prime})$ of $\bm{x}$ satisfying $\bm{0}\notin B(\bm{x},\sigma^{\prime})$ (since $\bm{x}\neq\bm{0}$) and a constant $K^{\prime}>0$, $|g(\bm{x})-g(\bm{y})|\leq K^{\prime}\|\bm{x}-\bm{y}\|\text{~{}holds},$ where $K^{\prime}\geq\max\{\max_{\bm{y}\in B(\bm{x},\sigma^{\prime})}\frac{\alpha}{\|\bm{y}\|^{s}},K_{1}\}$.

Below we show that the function $g(\bm{x})$ satisfies Assumption 3.2. It is trivial to prove that the function $g(\bm{x})$ satisfies conditions (a) and (b) in Assumption 3.2. Next we prove that the function $g(\bm{x})$ in (7) satisfies (c) in Assumption 3.2. Suppose that $T=t(\bm{x},\bm{d}(\cdot))<\infty$ and $\bm{y}=\bm{\psi}_{\bm{x}}^{\bm{d}}(T)$. Since $\bm{\psi}_{\bm{x}}^{\bm{d}}(\cdot):\mathbb{R}\rightarrow\mathbb{R}^{n}$ is continuous over $t$, $g(\bm{\psi}_{\bm{x}}^{\bm{d}}(t))$ is continuous over $t$ as well, implying that $g(\bm{\psi}_{\bm{x}}^{\bm{d}}(t))$ can attain the maximum $M$ over $t\in[0,T]$. Therefore, we have

[TABLE]

where $L_{g}$ is the Lipschitz constant of $g(\bm{x})$ over $\overline{\mathcal{X}}$, and $C,\overline{r},\sigma$ are defined in (2). Therefore, $\int_{0}^{\infty}g(\bm{\psi}_{\bm{x}}^{\bm{d}}(t))dt$ is finite if $t(\bm{x},\bm{d}(\cdot))<\infty$.

In the following we show that the function $g(\bm{x})$ satisfies (9) over $B(\bm{0},\overline{r})$.

Let $\bm{x},\bm{y}\in B(\bm{0},\overline{r})$, $M_{1}$ and $M_{2}$ are two positive constants such that $\texttt{dist}(\bm{z},\mathcal{X}_{\infty}^{c})\leq M_{1}$ and $M_{2}\leq\texttt{dist}(\bm{z},\mathcal{X}_{\infty}^{c})$ for $\bm{z}\in B(\bm{0},\overline{r})$ (Such $M_{1}$ and $M_{2}$ exist since $\texttt{dist}(\cdot,\mathcal{X}_{\infty}^{c}):B(\bm{0},\overline{r})\rightarrow(0,\infty)$ is Lipschitz continuous and $B(\bm{0},\overline{r})$ is a compact set with $B(\bm{0},\overline{r})\subset\mathcal{X}_{\infty}$ and $\partial B(0,\overline{r})\cap\partial\mathcal{X}_{\infty}=\emptyset$.), and $L_{\texttt{dist}}$ is the Lipschitz constant of the distance function $\texttt{dist}(\cdot,\mathcal{X}_{\infty}^{c})$ over $B(\bm{0},\overline{r})$, we have that

[TABLE]

where $M^{\prime}=s\max_{\bm{z}\in B(\bm{0},\overline{r})}\sum_{i=1}^{n}|z_{i}|$ and $K=\frac{\alpha\max\{M_{1}M^{\prime},L_{\texttt{dist}}\}}{M_{2}^{2}}$. The inequality \big{|}\|\bm{x}\|^{s}-\|\bm{y}\|^{s}\big{|}\leq M^{\prime}\max\{\|\bm{x}\|^{s-2},\|\bm{y}\|^{s-2}\}\big{|}\|\bm{x}\|-\|\bm{y}\|\big{|} is obtained in the following way:

[TABLE]

where $\|\bm{x}\|=\sqrt{\sum_{i=1}^{n}x_{i}^{2}},$ $\bm{\xi}=\lambda\bm{x}+(1-\lambda)\bm{y}$ and $\lambda$ is a constant falling within $(0,1)$.

*The proof is completed. *

Denote

[TABLE]

and

[TABLE]

where $\tilde{V}=-\int_{0}^{t}g(\bm{\psi}_{\bm{x}}^{\bm{d}}(\tau))d\tau-h^{\prime}(\bm{\psi}_{\bm{x}}^{\bm{d}}(t))$ and $\delta$ is some positive constant.

According to Theorem 3.1 in [18], we have the following conclusion that

$\mathcal{R}_{0}=\{\bm{x}\in\mathbb{R}^{n}\mid V(\bm{x})<\infty\}=\{\bm{x}\in\mathbb{R}^{n}\mid v(\bm{x})<1\}$. 2. 2.

$V(\bm{x})$ in (11) is continuous on $\mathcal{R}_{0}$. In addition, $\lim_{n\rightarrow\infty}V(\bm{x}_{n})=\infty$ if $\lim_{n\rightarrow\infty}\bm{x}_{n}=\bm{x}\notin\mathcal{R}_{0}$ or $\lim_{n\rightarrow\infty}\|\bm{x}_{n}\|=\infty.$

According to Proposition 4.2 in [18], $V(\bm{x})$ and $v(\bm{x})$ satisfy the following dynamic programming principle.

Lemma 3.6**.**

Assume that $G(\bm{x},t,\bm{d}(\cdot))=\int_{0}^{t}g(\bm{\psi}_{\bm{x}}^{\bm{d}}(\tau))d\tau$. Then the following assertions are satisfied:

for $\bm{x}\in\mathcal{R}_{0}$ and $t\geq 0$, we have:

[TABLE] 2. 2.

for $\bm{x}\in\mathbb{R}^{n}$ and $t\geq 0$, we have:

[TABLE]

We further exploit the Lipschitz continuity property of $V(\bm{x})$ and $v(\bm{x})$. The Lipschitz continuity property of $v(\bm{x})$ plays a key role in guaranteeing the existence of solutions to the constructed semi-definite program (22) theoretically, which will be introduced later.

Lemma 3.7**.**

Under Assumption 2.1 and Assumption 3.2, then

$V(\bm{x})$* in (11) is locally Lipschitz continuous over $\mathcal{R}_{0}$.* 2. 2.

$v(\bm{x})$* in (12) is locally Lipschitz continuous over $\mathbb{R}^{n}$.*

Proof 3.8**.**

1. Since $\sup_{\bm{d}(\cdot)\in\mathcal{D}}\|\bm{\psi}_{\bm{x}}^{\bm{d}}(t)-\bm{\psi}_{\bm{y}}^{\bm{d}}(t)\|\leq e^{L_{\bm{f}}t}\|\bm{x}-\bm{y}\|$ for $t\in[0,\infty)$ and $\bm{x},\bm{y}\in\mathbb{R}^{n}$, we have that for $\bm{x}_{0}\in\mathcal{R}_{0}$ and $t\in[0,\infty)$, there exist $\delta_{\bm{x}_{0},t}>0$ and $\rho>0$ such that $B(\bm{x}_{0},\delta_{\bm{x}_{0},t})\subseteq\mathcal{X}$ and $\mathtt{dist}(\bm{\psi}_{\bm{y}}^{\bm{d}}(\tau),\mathcal{X}^{c})\geq\frac{\rho}{2}$ for $\bm{y}\in B(\bm{x}_{0},\delta_{\bm{x}_{0},t})$, $\bm{d}(\cdot)\in\mathcal{D}$ and $\tau\in[0,t]$.

Therefore, for $\bm{y}\in B(\bm{x}_{0},\delta_{\bm{x}_{0}})$, where $B(\bm{x}_{0},\delta_{\bm{x}_{0}})\subseteq\mathcal{R}_{0}$, we obtain that

[TABLE]

According to Lemma 3.4,

[TABLE]

Thus, analogous to the proof of Proposition 4.2 in [8], we obtain that

[TABLE]

where $L_{S}$ is some positive constant.

As to $\sup_{\bm{d}(\cdot)\in\mathcal{D}}\sup_{t\in[0,\infty)}|h^{\prime}(\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(t))-h^{\prime}(\bm{\psi}_{\bm{y}}^{\bm{d}}(t))|$, following the proof of (ii) of Theorem 3.1 in [18] and the fact that $B(\bm{x}_{0},\delta_{\bm{x}_{0}})$ is compact (it indicates that there exists $T^{\prime}\in[0,\infty)$ such that $\bm{\psi}_{\bm{y}}^{\bm{d}}(t)\in B(\bm{0},\overline{r})$ for $t\geq T^{\prime}$, $\bm{d}(\cdot)\in\mathcal{D}$ and $\bm{y}\in B(\bm{x}_{0},\delta_{\bm{x}_{0}})$.), we have that there exists a non-negative constant $T$, which is independent of $\bm{d}(\cdot)\in\mathcal{D}$ and $\bm{y}\in B(\bm{x}_{0},\delta_{\bm{x}_{0}})$, such that

[TABLE]

where $T_{\bm{d}}\in[0,T]$. Taking $\delta_{0}=\min\{\delta_{\bm{x}_{0},T},\delta_{\bm{x}_{0}}\}$, we have that for $\bm{y}\in B(\bm{x}_{0},\delta_{0})$,

[TABLE]

where $L_{\bm{f}}$ is defined in (5) and $L_{h^{\prime}}$ is the Lipschitz constant of the function $h^{\prime}$ over the compact set $\overline{\Omega(B(\bm{x}_{0},\delta_{0}),[0,T])}$, which is the closure of the set of states visited by system (4) starting from $B(\bm{x}_{0},\delta_{0})$ within the time interval $[0,T]$ and thus is a subset of the set $\mathcal{X}$.

Combining (13) and (14), we have

[TABLE]

where $L=L_{h^{\prime}}e^{L_{\bm{f}}T}+L_{S}.$ Thus, $V(\bm{x})$ in (11) is Lipschitz continuous on the compact set $B(\bm{x}_{0},\delta_{0})$ and thus $V(\bm{x})$ is locally Lipschitz continuous on $\mathcal{R}_{0}$.

2. We will prove that $v(\bm{x})$ is locally Lipschitz continuous over $\mathbb{R}^{n}$ based on $0\leq v(\bm{x})\leq 1$ for $\bm{x}\in\mathbb{R}^{n}$ and the fact in Lemma 3.6 that for $\bm{x}\in\mathbb{R}^{n}$ and $t\geq 0$,

[TABLE]

For $\bm{x},\bm{y}\in\mathbb{R}^{n}$, we have that for any $\epsilon>0$, there exists a perturbation input $\bm{d}(\cdot)\in\mathcal{D}$ such that

[TABLE]

where $T>0$, $l>0$, $k>0$, $K>0$, and $L_{g}$ is the Lipschitz constant of the function $g(\bm{x})$ over the compact set $\Omega(B,[0,T])$, which is the closure of the set of states visited by system (4) starting from the set $B$ within the time interval $[0,T]$, and $B$ is a compact set covering $\bm{x}$ and $\bm{y}$.

*Thus, $v(\bm{x})$ is locally Lipschitz continuous over $\mathbb{R}^{n}$. *

Consequently, combining Theorem 4.4 in [18] and Lemma 3.7 we have:

Theorem 3.9**.**

The value function $V(\bm{x})$ in (11) is the unique locally Lipschitz continuous viscosity solution to the following equation

[TABLE]

Likewise, the value function $v(\bm{x})$ in (12) is the unique bounded and locally Lipschitz continuous viscosity solution to the generalized Zubov’s equation

[TABLE]

As a direct consequence of (16), we have that if a continuously differentiable function $u(\bm{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}$ satisfies (16), then $u(\bm{x})$ satisfies the constraints:

[TABLE]

Corollary 3.10**.**

*Assume a continuously differentiable function $u(\bm{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a solution to (17), then $v(\bm{x})\leq u(\bm{x})$ over $\bm{x}\in\mathbb{R}^{n}$ and consequently $\Omega=\{\bm{x}\mid u(\bm{x})<1\}\subset\mathcal{R}_{0}$ is a robust domain of uniform attraction. *

Proof 3.11**.**

It is obvious that (17) is equivalent to

[TABLE]

If $u(\bm{x})$ is a viscosity super-solution to (16), according to the comparison principle in Proposition 4.7 in [18], $v(\bm{x})\leq u(\bm{x})$ holds. Consequently, $\Omega=\{\bm{x}\mid u(\bm{x})<1\}\subset\mathcal{R}_{0}$ is a robust domain of uniform attraction. In the following we show that $u(\bm{x})$ is a viscosity super-solution to (16)

Let’s first recall the concept of viscosity super-solution to (16). A lower semicontinuous function $u_{l}(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a viscosity super-solution of (16) [18] if for all $\phi(\bm{x})\in C^{1}(\mathbb{R}^{n})$ such that $u_{l}(\bm{x})-\phi(\bm{x})$ has a local minimum at $\bm{x}_{0}$, we have

[TABLE]

Since $u(\bm{x})$ satisfies (18), we just show that $\inf_{\bm{d}\in D}\{-\frac{\partial\phi(\bm{x})}{\partial\bm{x}}\mid_{\bm{x}=\bm{x}_{0}}\bm{F}(\bm{x}_{0},\bm{d})-\delta g(\bm{x}_{0})(1-u(\bm{x}_{0}))\}\geq 0$, where $\phi(\bm{x})\in C^{1}(\mathbb{R}^{n})$ and $u(\bm{x})-\phi(\bm{x})$ has a local minimum at $\bm{x}_{0}$. Without loss of generality, we assume that $u(\bm{x}_{0})-\phi(\bm{x}_{0})=0$. There exists $\delta_{0}>0$ such that

[TABLE]

Suppose that $\inf_{\bm{d}\in D}\{-\frac{\partial\phi(\bm{x})}{\partial\bm{x}}\mid_{\bm{x}=\bm{x}_{0}}\bm{F}(\bm{x}_{0},\bm{d})-\delta g(\bm{x}_{0})(1-u(\bm{x}_{0}))\}\geq 0$ does not hold, i.e., $\inf_{\bm{d}\in D}\{-\frac{\partial\phi(\bm{x})}{\partial\bm{x}}\mid_{\bm{x}=\bm{x}_{0}}\bm{F}(\bm{x}_{0},\bm{d})-\delta g(\bm{x}_{0})(1-\phi(\bm{x}_{0}))\}<0$. Then there exists $\epsilon>0$ such that

[TABLE]

Further, there exists $\bm{d}_{1}\in D$ such that

[TABLE]

and consequently there exists $\delta^{\prime}>0$ with $\delta^{\prime}\leq\delta_{0}$ such that

[TABLE]

Since $\bm{\psi}_{\bm{x}_{0}}^{\bm{d}}(t)$ is absolutely continuous over $t$ for $\bm{d}(\cdot)\in\mathcal{D}$, there exists $\theta>0$ such that for $\tau\in[0,\theta]$,

[TABLE]

where $\bm{d}^{\prime}_{1}(\cdot)\in\mathcal{D}$ with $\bm{d}^{\prime}_{1}(\tau)=\bm{d}_{1}$ for $\tau\in[0,\theta]$. Therefore, we have that for $\tau\in[0,\theta]$,

[TABLE]

where $\bm{d}^{\prime}_{1}(\cdot)\in\mathcal{D}$ with $\bm{d}^{\prime}_{1}(\tau)=\bm{d}_{1}$ for $\tau\in[0,\theta]$.

By applying Gronwall’s inequality [17] to (20) together with (19) with the time interval [0, $\theta$], we have that

[TABLE]

where $G=\int_{0}^{\theta}g(\bm{\psi}_{\bm{x}_{0}}^{\bm{d}^{\prime}_{1}}(t))dt$. Therefore,

[TABLE]

which contradicts the fact that

[TABLE]

*Therefore, we conclude that $u(\bm{x})$ is a viscosity super-solution to (16). *

From Corollary 3.10 we observe that a robust domain of attraction can be found by solving (17) rather than (16). However, $u(\bm{x})$ is required to satisfy (17) over $\mathbb{R}^{n}$, which is a strong condition. This requirement renders the search for a continuously differentiable solution to (17) nonetheless nontrivial. Regarding this issue, we further relax this condition and restrict the search for a continuously differentiable function $u(\bm{x})$ in the compact set $B(\bm{0},R)\setminus\mathcal{X}_{\infty}$, where $B(\bm{0},R)$ is defined in (3). Also, since $\bm{F}(\bm{x},\bm{d})=\bm{f}(\bm{x},\bm{d})$ for $(\bm{x},\bm{d})\in B(\bm{0},R)\times D$ and $g(\bm{x})=q(\bm{x})$ for $\bm{x}\in B(\bm{0},R)\setminus\mathcal{X}_{\infty}$, we obtain the following system of constraints:

[TABLE]

Theorem 3.12**.**

*Let $u(\bm{x})$ be a continuously differentiable solution to (21) and $\delta$ be a positive value, then $\Omega=\{\bm{x}\in B(\bm{0},R)\mid u(\bm{x})<1\}$ is a robust domain of attraction. *

Proof 3.13**.**

Firstly, since $\bm{F}(\bm{x},\bm{d})=\bm{f}(\bm{x},\bm{d})$ for $(\bm{x},\bm{d})\in B(\bm{0},R)\times D$ and $g(\bm{x})=q(\bm{x})$ over $\bm{x}\in B(\bm{0},R)\setminus\mathcal{X}_{\infty}$, (21) is equivalent to

[TABLE]

Since $u(\bm{x})-1\geq 0$ for $\bm{x}\in B(\bm{0},R)\setminus\mathcal{X}$, $\Omega\subset\mathcal{X}$ holds. Next we prove that every possible trajectory initialized in the set $\Omega$ will approach the equilibrium state $\bm{0}$ eventually while never leaving the state constraint set $\mathcal{X}$.

Assume that there exist $\bm{y}\in\Omega$, a perturbation input $\bm{d}^{\prime}(\cdot)\in\mathcal{D}$ and $\tau>0$ such that

[TABLE]

and

[TABLE]

Obviously, $\bm{y}\notin\mathcal{X}_{\infty}$ and $\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(t)\notin\mathcal{X}_{\infty}$ for $t\in[0,\tau]$ since $\mathcal{X}_{\infty}$ is a robust domain of attraction. That is,

[TABLE]

Applying Gronwall’s inequality [17] to $-\frac{\partial u}{\partial\bm{x}}\bm{F}(\bm{x},\bm{d})-\delta g(\bm{x})(1-u(\bm{x}))\geq 0$ with the time interval [0, $\tau$], we have that

[TABLE]

where $G=\int_{0}^{\tau}g(\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(t))dt>0$. Therefore, $u(\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(\tau))<1$. However, since $\mathcal{X}\subseteq B(\bm{0},R)$ and $\partial\mathcal{X}\cap\partial B(\bm{0},R)=\emptyset$, $\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(\tau)\in B(\bm{0},R)\setminus\mathcal{X}$ holds and consequently $u(\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(\tau))\geq 1$. This is a contradiction. Thus, every possible trajectory initialized in the set $\Omega$ never leaves the set $\mathcal{X}$.

Lastly, we prove that every possible trajectory initialized in the set $\Omega$ approaches the equilibrium state $\bm{0}$ eventually. Since every possible trajectory initialized in the set $\mathcal{X}_{\infty}$ approaches the equilibrium state $\bm{0}$ eventually, it is enough to prove that every possible trajectory initialized in the set $\Omega\setminus\mathcal{X}_{\infty}$ will enter the set $\mathcal{X}_{\infty}$ in finite time. Assume that there exist $\bm{y}\in\Omega$ and a perturbation input $\bm{d}^{\prime}(\cdot)$ such that $\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(t)\notin\mathcal{X}_{\infty}$ for all $t\geq 0$. According to the second constraint and the third constraint in (21), we have $u(\bm{x})\geq 0$ for $\bm{x}\in B(\bm{0},R)\setminus\mathcal{X}_{\infty}$. Also, since $\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(t)\in B(\bm{0},R)$ for all $t\geq 0$, $u(\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(t))\geq 0$ holds for all $t\geq 0$. Moreover, applying Gronwall’s inequality [17] again to

[TABLE]

with the time interval [0, $\tau$] for $\tau>0$, we have that $u(\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(\tau))<1$. This implies that

[TABLE]

Also, since

[TABLE]

we obtain that

[TABLE]

Consequently, we conclude that

[TABLE]

contradicting the fact that $u(\bm{\psi}_{\bm{y}}^{\bm{d}^{\prime}}(t))\geq 0$ holds for $t\geq 0$. Therefore, every possible trajectory initialized in the set $\Omega$ will enter the set $\mathcal{X}_{\infty}$ in finite time and consequently will asymptotically approach the equilibrium state $\bm{0}$.

*Therefore, $\Omega$ is a robust domain of attraction. *

When $u(\bm{x})$ in (21) is a polynomial in $\mathbb{R}[\bm{x}]$, based on the sum-of-squares decomposition for multivariate polynomials, (21) is recast as the following semi-definite program:

According to Theorem 3.12, $\mathcal{R}_{u}=\{\bm{x}\in B(\bm{0},R)\mid u(\bm{x})<1\}$ is a robust domain of attraction, where $u(\bm{x})\in\mathbb{R}[\bm{x}]$ is the solution to (22).

Remark 3.14**.**

$\{\bm{x}\in B(\bm{0},R)\mid u(\bm{x})<1\}$* is still a robust domain of attraction if the origin $\bm{0}$ is asymptotically stable for (1) rather than uniformly locally exponentially stable, where $u(\bm{x})$ is the solution to (22). The proof of Theorem 3.12 does not require that the equilibrium state $\bm{0}$ is uniformly locally exponentially stable. *

3.3 Analysis of (22)

In this subsection we exploit some properties pertinent to (22) and show that there exist solutions to (22) under appropriate assumptions. Moreover, we show that there exists a sequence of solutions to (22) such that their strict one sub-level sets approximate the interior of the maximal robust domain of attraction in measure.

{assumption}

One of the polynomials defining the set $D$ is equal to $h_{i}^{D}:=\|\bm{d}\|^{2}-R_{D}$ for some constant $R_{D}\geq 0$.

As argued in [24], Assumption 3.3 is without loss of generality since the set $D$ is compact, the redundant constraint $\|\bm{d}\|^{2}-R_{D}-1\leq 0$ can always be incorporated into the description of $D$ for sufficiently large $R_{D}$.

We in the following show that given an arbitrary $\epsilon>0$, there exists a polynomial solution $p(\bm{x})$ to (22) such that $|p(\bm{x})-v(\bm{x})|<\epsilon$ holds for $\bm{x}\in B(\bm{0},R)$. Before this, we introduce a lemma from [26].

Lemma 3.15** (Lemma B.5 in [26]).**

Let $B(\bm{0},R)$ be a compact subset in $\mathbb{R}^{n}$ and $u(\bm{x}):B(\bm{0},R)\rightarrow\mathbb{R}$ be a locally Lipschitz function. If there exists a continuous function $\rho:B(\bm{0},R)\rightarrow\mathbb{R}$ such that for each $\bm{d}\in D$,

[TABLE]

where $\mathcal{L}u(\bm{x})=\nabla_{\bm{x}}u(\bm{x})\cdot\bm{f}(\bm{x},\bm{d})=\frac{\partial u}{\partial\bm{x}}\bm{f}(\bm{x},\bm{d})$ (recall that $\nabla_{\bm{x}}u(\bm{x})$ is defined a.e., since $u$ is locally Lipschitz.), then for any given $\epsilon>0$, there exists some smooth function $\psi(\bm{x})$ defined on $B(\bm{0},R)$ such that

[TABLE]

*over $\bm{x}\in B(\bm{0},R)$. *

Theorem 3.16**.**

Under Assumption 3.3, if $\delta$ is a positive value larger than or equal to one in (22), then for any $\epsilon>0$ there exists a polynomial solution $p(\bm{x})$ to (22) such that

[TABLE]

Proof 3.17**.**

When $\delta$ is a positive value larger than or equal to one, we have that $u(\bm{x})+(1-h_{j}(\bm{x}))-1\geq u(\bm{x})+(1-h_{j}(\bm{x}))^{\delta}-1\geq 0$ for $\bm{x}\in\overline{\mathcal{X}\setminus\mathcal{X}_{\infty}}$, i.e., $u(\bm{x})-h_{j}(\bm{x})\geq u(\bm{x})+(1-h_{j}(\bm{x}))^{\delta}-1\geq 0$ for $\bm{x}\in\overline{\mathcal{X}\setminus\mathcal{X}_{\infty}}$, $j=1,\ldots,n_{\mathcal{X}}$. Also, since $v(\bm{x})$ in (12) satisfies (16), we have $v(\bm{x})$ satisfies (21). Consequently, for any $\epsilon_{1}>0$, $v^{\prime}(\bm{x})=v(\bm{x})+\epsilon_{1}$ satisfies the following constraints:

[TABLE]

According to Lemma 3.7, $v(\bm{x})$ is locally Lipschitz continuous over $\bm{x}\in B(\bm{0},R)$. Therefore, $v^{\prime}(\bm{x})$ is locally Lipschitz continuous over $\bm{x}\in B(\bm{0},R)$ as well. According to Lemma 3.15, we have that for any $\epsilon_{2}<\frac{\epsilon_{1}\alpha\delta}{2}$ with $\epsilon_{2}>0$, there exists a continuous function $p^{\prime}(\bm{x})$ such that

[TABLE]

and

[TABLE]

Since $B(\bm{0},R)$ is compact, there exists a polynomial $p(\bm{x})$ of sufficiently high degree such that

[TABLE]

[TABLE]

Then we have

[TABLE]

and

[TABLE]

Thus, we have

[TABLE]

*The polynomial $p(\bm{x})$ is therefore strictly feasible in (22), which follows from the classical Putinar’s Positivstellensatz [33]. Since $\epsilon_{1}$ is arbitrary and $\epsilon_{2}<\epsilon_{1}\alpha\delta$, the conclusion in Theorem 3.16 holds. *

Given $(\epsilon_{k})_{k=1}^{\infty}$ with $\epsilon_{k}>0$ and $\lim_{k\rightarrow\infty}\epsilon_{k}=0$ with $k\in\mathbb{N}$, according to Theorem 3.16, there exists a sequence

[TABLE]

satisfying (22) such that $0\leq p_{k}(\bm{x})-v(\bm{x})<\epsilon_{k}.$ Denote

[TABLE]

we next show that $\mathcal{R}_{k,0}$ inner-approximates the interior of the maximal robust domain of attraction in measure with $k$ approaching infinity.

Theorem 3.18**.**

Let \big{(}p_{k}(\bm{x})\big{)}_{k=1}^{\infty} and $\mathcal{R}_{k,0}$ be the sequence in (23) and the set in (24) respectively. Then the set $\mathcal{R}_{k,0}$ converges to the interior of the maximal robust domain of attraction from inside in measure with $k$ tending towards infinity, i.e.,

[TABLE]

Proof 3.19**.**

*Following the proof of Theorem 3 in [25], we have that $\lim_{k\rightarrow\infty}\mu(\mathcal{R}_{0}\setminus\mathcal{R}_{k,0})=0$, where $\mathcal{R}_{0}=\{\bm{x}\in\mathbb{R}^{n}\mid v(\bm{x})<0\}$. According to Lemma 2.3, we have that $\lim_{k\rightarrow\infty}\mu(\mathcal{R}^{\circ}\setminus\mathcal{R}_{k,0})=0.$ *

4 Examples and Discussions

In this section we illustrate our approach with five examples. All computations were performed on an i7-P51s 2.6GHz CPU with 4GB RAM running Windows 10. For the numerical implementation, we formulate the sum-of-squares problem (22) using the Matlab package YALMIP [27] and employ Mosek of the academic version [29] as a semi-definite programming solver. The parameters that control the performance of our method are given in Table 1.

Example 4.1**.**

Consider a one-dimensional system adapted from [20], which is given by

[TABLE]

with $\mathcal{X}=\{x\in\mathbb{R}\mid x^{2}<1\}$ and $D=\{d\in\mathbb{R}\mid 100d^{2}-1\leq 0\}$.

The origin is a locally uniformly exponentially stable state. The maximal robust domain of attraction in this case is determined analytically as $\mathcal{R}=(-0.5,0.4)$. Theorem 3.12 indicates that the strict one sub-level set of the approximating polynomial $u$ computed by solving (22) is a robust domain of attraction. Plots of the computed robust domains of attraction for approximating polynomials of degree $k=8,10,16,24$ are presented in Fig. 1. The visualized results in Fig. 1 further confirm that the strict one sub-level set of the approximating polynomial $u(\bm{x})$ returned by solving (22) is indeed a robust domain of attraction. The relative volume errors, which are computed approximately by Monte Carlo integration, are also reported in Table 2. From Fig. 1 and Table 2, we observe fairly good tightness of the estimates since $k=10$.

Example 4.2**.**

The second example considers scaled version of the reversed-time Van der Pol oscillator free of perturbations [20] given by

[TABLE]

with $\mathcal{X}=\{\bm{x}\mid\frac{x^{2}+y^{2}}{1.21}<1\}$.

The origin is a locally uniformly exponentially stable state. For this example, there exists a limit cycle, which is the boundary of the maximal robust domain of attraction. The limit cycle is included in $\mathcal{X}$. Theorem 3.12 indicates that the strict one sub-level set of the approximating polynomial $u(\bm{x})$ computed by solving (22) is a robust domain of attraction. Plots of computed robust domains of attraction for approximating polynomials of degree $k=8,16$, are shown in Fig. 2. The visualized results in Fig. 2 further confirm that the strict one sub-level set of the approximating polynomial $u(\bm{x})$ returned by solving (22) is indeed a robust domain of attraction. In order to quantitatively assess the quality of computed robust domains of attraction, we use the simulation technique to synthesize an estimate $\tilde{\mathcal{R}}$ of the maximal robust domain of attraction by gridding the state space, and then compute the relative volume errors approximately according to the formula $(1-\frac{\text{number of grid states in }\mathcal{R}_{k}}{\text{number of grid states in }\tilde{\mathcal{R}}})\times 100\%$, where $\mathcal{R}_{k}$ is the robust domain of attraction formed by the approximating polynomial of degree $k$. The estimate $\tilde{\mathcal{R}}$ is shown in Fig. 2. The visualized results in Fig. 2 indicate that the estimate $\tilde{\mathcal{R}}$ approximates the maximal robust domain of attraction quite well. The relative volume errors are reported in Table 3. From Fig. 2 and Table 3 we observe fairly good tightness of estimates since $k=8$.

Example 4.3**.**

In this example we consider a chemical oscillator from [30]. The simplest, but chemically plausible trimolecular reaction is

[TABLE]

in which species $X$ is in dynamical equilibrium with species $A$ with a forward rate of reaction $k_{1}$ and a backward rate of reaction $k_{-1}$, and so on. Using the law of mass action, and non-dimensionalising the equations, we have

[TABLE]

where $x$, $y$ are the non-dimensional concentrations of $X$ and $Y$, and $a$, $b$ are non-negative constant parameters that depend on the concentrations of $A$ and $B$.

Like [30], we take $a=0.5$ and $b=0.5$. The system has a locally uniformly exponentially stable state $(1,0.5)$. Since $x$ and $y$ are the non-dimensional concentrations of $X$ and $Y$, they are naturally positive, i.e. $x>0$ and $y>0$. On the other hand, they should have upper bound on the concentrations. In this paper we impose the inequality constraint $(x-1)^{2}+4(y-0.5)^{2}<1$.

After translating the equilibrium $(1,0.5)$ to the origin $(0,0)$ and making $x_{1}=x,x_{2}=2y$, we obtain the equivalent system of interest in this example,

[TABLE]

with $\mathcal{X}=\{\bm{x}\mid x_{1}^{2}+x_{2}^{2}<1\}$.

The origin is a locally uniformly exponentially stable state for system (25). Theorem 3.12 indicates that the strict one sub-level set of the approximating polynomial $u(\bm{x})$ computed by solving (22) is a robust domain of attraction. Plots of computed robust domains of attraction for approximating polynomials of degree $k=4,8,12$, are shown in Fig. 3. The visualized results in Fig. 3 further confirm that the strict one sub-level set of the approximating polynomial $u(\bm{x})$ returned by solving (22) is indeed a robust domain of attraction. Like Example 4.2, we use the simulation technique to quantitatively assess the quality of computed robust domains of attraction. The relative volume errors are listed in Table 4. From Fig. 3 and Table 4 we observe fairly good tightness of estimates since $k=8$.

Example 4.4**.**

Consider a system from [18], whose dynamics are described by

[TABLE]

The origin for system (26) is locally uniformly exponentially stable. In this example $D=[4.9,5.1]$ and $\mathcal{X}=\{\bm{x}\mid x^{2}+y^{2}<1\}$. In order to fit (22), we reformulate (26) as the following equivalent system

[TABLE]

where $D=\{d\in\mathbb{R}\mid 100d^{2}-1\leq 0\}$ and $\mathcal{X}=\{\bm{x}\mid x^{2}+y^{2}<1\}$.

Theorem 3.12 indicates that the strict one sub-level set of the solution to (22) is a robust domain of attraction. Plots of computed robust domains of attraction $\mathcal{R}_{k}$, $k=4,6,10$, are shown in Fig. 4. We also give an estimate of the maximal robust domain of attraction by simulation methods and estimate the relative volume errors as in Example 4.2. From Table 5, which lists the relative volume errors, we observe fairly good tightness of the estimates since $k=4$.

Example 4.5**.**

Consider a seven-dimensional system, which is mainly employed to illustrate the scalability issue of our semi-definite programming based method in dealing with high dimensional system.

[TABLE]

where $D=\{d\in\mathbb{R}\mid d^{2}+0.75-1\leq 0\}$ and $\mathcal{X}=\{\bm{x}\mid\|\bm{x}\|^{2}<1\}$.

The equilibrium state $\bm{0}$ is locally uniformly exponentially stable. Theorem 3.12 indicates that the strict one sub-level set of the solution to (22) is a robust domain of attraction. Plots of computed robust domains of attraction for approximating polynomials of degree $k=3,4,5$, on planes $x_{1}-x_{2}$ with $x_{3}=x_{4}=x_{5}=x_{6}=x_{7}=0$ and $x_{1}-x_{7}$ with $x_{2}=x_{3}=x_{4}=x_{5}=x_{6}=0$ are shown in Fig. 5. In order to shed light on the accuracy of the computed domains of attraction, we use the simulation technique to synthesize coarse estimations of the maximal robust domain of attraction on planes $x_{1}-x_{2}$ with $x_{3}=x_{4}=x_{5}=x_{6}=x_{7}=0$ and $x_{1}-x_{7}$ with $x_{2}=x_{3}=x_{4}=x_{5}=x_{6}=0$ by taking initial states in the state spaces $\{\bm{x}\mid\|\bm{x}\|^{2}<1\wedge x_{3}=0\wedge x_{4}=0\wedge x_{5}=0\wedge x_{6}=0\wedge x_{7}=0\}$ and $\{\bm{x}\mid\|\bm{x}\|^{2}<1\wedge x_{2}=0\wedge x_{3}=0\wedge x_{4}=0\wedge x_{5}=0\wedge x_{6}=0\}$, respectively. They are the gray regions in Fig. 5. The relative volume errors on both planes $x_{1}-x_{2}$ with $x_{3}=x_{4}=x_{5}=x_{6}=x_{7}=0$ and $x_{1}-x_{7}$ with $x_{2}=x_{3}=x_{4}=x_{5}=x_{6}=0$, which are computed following the way in Example 4.2, are reported in Table 6.

Based on Examples 4.1$\sim$4.5, we conclude that approximating polynomials of higher degree would return less conservative robust domains of attraction. Although the size of the semi-definite program in (22) grows extremely fast with the number of state and perturbation variables and the degree of the polynomials in (22), it is worth emphasizing that we are dealing with nonlinear non-convex infinite-dimensional problems by solving a semi-definite programming problem, which is relatively simple to implement. Yet, despite the difficulty of the problems considered, the constructed semi-definite program (22) possess solutions whose strict one sub-level sets inner-approximate the interior of the maximal robust domain of attraction in measure under appropriate assumptions according to Theorem 3.18. In order to improve the scalability issue of our method and further apply it to higher dimensional systems, some techniques such as exploiting the algebraic structure [32] of the semi-definition programming (22) and using template polynomials such as (scaled-) diagonally-dominant-sums-of-squares polynomials [28, 2, 3] would facilitate such gains.

In the rest we further give a brief discussion on the other parameters that control the performance of the semi-definite program (22) in terms of relative volume error estimations of computed robust domains of attraction to the maximal robust domain of attraction. The parameters related to the degree of polynomials in Table 1 are already discussed either above or in the existing literature pertinent to the sum-of-squares optimization: polynomials of higher degree in (22) would return less conservative robust domains of attraction generally. Thus, we just discuss the parameters $\alpha$, $\delta$ and $R$. Finding an optimal combination of these parameters to obtain the least conservative estimates is not the focus of this discussion and would be discussed in detail in the future work. The parameters $\alpha$ and $R$ could reflect the size of the sets $\mathcal{X}_{\infty}$ and $B(\bm{0},R)$ in (22), respectively. The analysis on these three parameters is based on the other parameter values, i.e., the degrees of polynomials, listed in Table 1, and is summarized in Table 7 $\sim$ 9. “-” in Table 7 $\sim$ 9 means that (22) did not return a feasible solution and thus we did not obtain the relative volume error estimation. For the convenience of comparison we also add the relative volume error estimations listed in Table 2 $\sim$ 6 into Table 7 $\sim$ 9.

Table 7 $\sim$ 9 indicate that the parameters $\delta$, $\alpha$ and $R$ definitely influence the performance of (22) for some cases. Increasing $\delta$ and $\alpha$ and/or decreasing $R$ may result in less conservative estimates, but this does not hold always. This influence is weakening with the degree of approximating polynomials increasing.

5 Conclusion

In this paper a semi-definite programming based method was proposed for synthesizing robust domains of attraction for state-constrained perturbed polynomial systems. The semi-definite program, which falls within the convex programming framework and can be solved in polynomial time via interior-point methods, was constructed from a generalized Zubov’s equation. Under appropriate assumptions the existence of solutions to the constructed semi-definite program is guaranteed and there exists a sequence of solutions such that their strict one sub-level sets inner-approximate the interior of the maximal robust domain of attraction in measure. Finally, we evaluated the performance of the method on five examples.