SOSP估计ROA理论

学习完SOSTOOLS工具箱后，开始学习平方和规划（Sum-of-Square Programming，SOSP）理论估计非线性系统的吸引域（Region of Attraction, ROA）。因为，非线性系统的全局渐近稳定性很难得到验证，一般都是得到一个局部围绕平衡点的吸引域。

SOS的定义

对于有 $n$ 个实变量的多元多项式 $p(x)$，其中 $x$ 是分量为 $x_1,...,x_n$ 的向量。如果存在多项式 $f_i(x)$ 使得： $$p(x)=\sum _{i=1}^M{f}_i^2(x)\text{\hspace{1em}(1)}$$那么多元多项式 $p(x)$ 就被称为 SOS。

SOS多项式的集合表示为 $\sum [x]$，如果 $p(x)$ 是SOS，那么它是半正定的（平方的和都是非负的）。

对于 $p(x)$ 来说，它的阶数就是最高阶单项式的阶数，比如
$$p(x)=x_1^6+x_2^3{x}_3^4$$阶数为 $3+4=7$。

如果 $p(x)$ 是SOS，可以写成：
$$p(x)=Z^T(x)Q_{0}Z(x)=Z^T(x)(Q_0+\sum _{i=1}^N{\lambda }_i{M}_i)Z(x)\text{\hspace{1em}(2)}$$
其中 $Q_0>0$ 是Gram矩阵。
因此，如果存在 ${\lambda }_i$ 使得
$$Q=Q_0+\sum _{i=1}^N{\lambda }_i{M}_i\ge 0\text{\hspace{1em}(3)}$$
则 $p(x)$ 是SOS，这是一个LMI的可行性问题。

综上，我们可以总结一下：并不是所有的非负的多项式都是SOS，但是所有的SOS都是非负的多项式，也就是SOS是非负多项式的充分不必要条件。但是检查SOS分解是否存在在计算上可以通过多项式的复杂性进行分解，通过此来寻找非线性系统的Lyapunov函数。

广义S-procedure

S-procedure的目的就是将一些不是凸约束的问题转化成LMI约束。
考虑 $p_0(x)$，$x\in D$
$$p(x)\ge 0:{\mathbb{R}}^n⟹\mathbb{R}$$为了限制$x$的域，需要$p_i(x)\ge 0,...,p_m(x)\ge 0$
则有：
∩ i = 1 m { x ∈ R n ∣ p i ( x ) ≥ 0 } ⊆ { x ∈ R n ∣ p 0 ( x ) ≥ 0 } (4) \cap_{i=1}^m\{x\in\mathbb{R}^n|p_i{(x)}\geq0\}\subseteq\{x\in\mathbb{R}^n|p_0(x)\geq0\}\tag{4} ∩i=1m​{x∈Rn∣pi​(x)≥0}⊆{x∈Rn∣p0​(x)≥0}(4)
存在半正定函数 $q_1(x),...,q_m(x)$使得：
$$p_0(x)-\sum _{i=1}^m{q}_i(x)p_i(x)\ge 0\text{\hspace{1em}(5)}$$
我的理解大概就是：先判断$p_0(x)-{\sum }_{i=1}^m{q}_i(x)p_i(x)\ge 0$成立，同时，$p_i(x)\ge 0,...,p_m(x)\ge 0$，那么可以得到 $p(x)\ge 0$

吸引域估计

考虑以下形式的自治非线性系统：$\dot{x}=f(x(t))$，假设 $f(0)=0$。如果 $\xi (x_0,t)$是在$t$时刻的解，初始条件是$\xi (x_0,0)=x_0$。吸引域是 { x 0 ∈ R n : l i m t → ∞ ξ ( x 0 , t ) = 0 } \{x_0\in\mathbb{R}^n: lim_{t\rightarrow\infty}\xi(x_0,t)=0\} {x0​∈Rn:limt→∞​ξ(x0​,t)=0}。对于$t>0$，$x_0\in \mathcal{M}$意味着$\xi (x_0,t)\in \mathcal{M}$，那么，集合 $\mathcal{M}$ 在 $f(x(t))$下是不变的。
Theorem 1（Lyapunov稳定性的充分条件） 考虑以下形式的自治非线性系统：$\dot{x}=f(x(t))$，$x$是状态向量，$f$ 是局部 Lipschitz 函数，描述系统动力学。该系统局部稳定的条件是：存在连续可微函数$V$，使得下面两个条件成立：
(1) 对于所有$x\ne 0$，有$V(x)>0$，同时$V(0)=0$
(2) $\dot{V}(x)\le 0$，即 $\nabla xV(x)f(x)\le 0$
此外，如果 $\nabla xV(x)f(x)<0$，平衡点是局部渐近稳定的。

具有多项式向量场的Lyapunov函数可以通过上述两个条件放宽到SOS约束来构造，假设 ${\phi }_i$ 是正定多项式定义为：
$${\phi }_i(x)=\sum _{j=1}^n{ϵ}_{ij}{x}_j^2\text{\hspace{1em}(6 )}$$ 其中，${ϵ}_{ij}$是比较小的正实数。可以使用SOSP在原点建立平衡的局部渐近稳定性

Program 1：给定具有多项式向量场的ODE系统，如果能找到$V(x)$，使得$V(0)=0$，并且 $q(x)\in \sum [x]$，
$$V(x)-{\phi }_1(x)\in \sum [x]\text{\hspace{1em}(7)}$$ $$-{\nabla }_{x}V(x)f(x)-{\phi }_2(x)+q(x)g(x)\in \sum [x]\text{\hspace{1em}(8)}$$则可以检查原点平衡的局部渐近稳定性。
其中 $(7)$ 是由 $V(x)>0$，$V(0)=0$ 得到，$(8)$ 是由 $\nabla xV(x)f(x)\le 0$ 和 S-procedure得到。

对于具有有理形式的系统
$$\dot{x}=\frac{n(x)}{d(x)},d(x)>0$$
将约束$(8)$乘以$d(x)$得到了可以使用SOS分解的形式，那么就可以使用 Program 1 来检查原点的稳定性，可以得到：
Program 2
$$V(x)-{\phi }_1(x)\in \sum [x]\text{\hspace{1em}(9)}$$ $$-{\nabla }_{x}V(x)n(x)-{\phi }_2(x)d(x)+q(x)g(x)d(x)\in \sum [x]\text{\hspace{1em}(10)}$$那么在建立平衡点的稳定性之后，就可以尝试去估计ROA了。
Lemma 1 假设 $\gamma >0$ , 并存在连续可微函数 $V:{\mathbb{R}}^n\to \mathbb{R}$ 使得：
(1) $V(0)=0$, 且对所有非零 $x\in {\mathbb{R}}^n$ 有 $V(x)>0$
(2) M V , γ = { x ∈ R n : V ( x ) ≤ γ } \mathcal{M}_{V,\gamma}=\{x\in\mathbb{R}^n:V(x)\leq\gamma\} MV,γ​={x∈Rn:V(x)≤γ}是有界的
(3) M V , γ \ { 0 } ⊂ { x ∈ R n : ∇ x V ( x ) f ( x ) ≤ 0 } \mathcal{M}_{V,\gamma} \backslash \{0\}\subset\{x\in\mathbb{R}^n: \nabla xV(x)f(x)\leq0\} MV,γ​\{0}⊂{x∈Rn:∇xV(x)f(x)≤0}
那么，对于所有的$x_0\in {\mathcal{M}}_{V,\gamma }$，非线性系统的解 $\xi (x_0,t)$ 存在于 $[0,\infty ]$ 上，且对于 $t\ge 0$，$li{m}_{t\to \infty }\xi (x_0,t)=0$。
如果 Lemma 1 满足，则 ${\mathcal{M}}_{V,\gamma }$ 是平衡点处的 ROA 不变子集。对于在 ${\mathcal{M}}_{V,\gamma }$ 内的任何初始条件 $x_0$，非线性系统的开环解都收敛于原点。
$(3)$中两个集合均由不等式描述，通过 S-procedure 来建立包含关系。如果 $\phi (x)$是正定的，$q(x)$ 是半正定的，同时 $$q(x)(V(x)-\gamma )\ge (\phi (x)+{\nabla }_{x}Vf(x))$$那么 $(3)$ 成立。

所以，要寻找满足 Lemma 1 非负条件的 Lyapunov函数的问题可以使用 SOS 松弛方法，将其转化为 SDP 问题。通过找到使得稳定性条件成立的 Lyapunov 函数的最大水平集 $\gamma$ 来估计吸引域的大小。

Program 3 有理系统 $\dot{x}=\frac{n(x)}{d(x)},d(x)>0$ 平衡点周围的吸引域可以通过找到使得 $\beta$ 最大化的 $V$, $q_1\in \sum [x]$, $q_2\in \sum [x]$ 来估计，这样：
$$max\beta$$ $$V(x)-{\phi }_1(x)\in \sum [x]\text{\hspace{1em}(11)}$$$$q_1(x)(s(x)-\beta )-(V(x)-\gamma )\in \sum [x]\text{\hspace{1em}(12)}$$$$-{\nabla }_{x}V(x)n(x)-{\phi }_2(x)d(x)+q_2(x)d(x)(V(x)-\gamma )\in \sum [x]\text{\hspace{1em}(13)}$$

条件 $(11)$ 确保 Lemma 1.1 满足，保证 $V(x)$的正定性
条件 $(12)$ 确保 M s , β = { s ( x ) ≤ β } ⊂ M V , γ = { V ( x ) ≤ γ } \mathcal{M}_{s,\beta}=\{s(x)\leq\beta\}\subset\mathcal{M}_{V,\gamma}=\{V(x)\leq\gamma\} Ms,β​={s(x)≤β}⊂MV,γ​={V(x)≤γ}，其中， $s(x)$ 是给定的形状函数
条件 $(13)$ 确保 Lemma 1.3 满足，强迫 $\dot{V}(x)$ 为负值。
因此要使得 $\beta$ 最大，来使得 ROA 的面积最大

算法： Program 3 可以通过重复以下步骤来迭代求解，直到 $\beta$ 不再改进。
(1) $\gamma$ step：初始化 $V$，通过找到 $q_2(x)$ 使得 (13) 成立来最大化 $\gamma$
(2) $\beta$ step：保持 $V$ $\gamma$ $q_2(x)$ 固定，并通过找到 $q_1(x)$使得 (12) 成立来最大化 $\beta$
(3) $V$ step：保持 $\gamma$ $\beta$ $q_1(x)$ $q_2(x)$ 固定，并找到 $V$ 使得 (11)-(13) 成立

$\gamma$：通过 $\gamma$ 的二分法来求解约束(13)，直到 $\gamma$的下估计值和上估计值之间的差异达到某个指定阈值。
$\beta$：通过 $\beta$ 的二分法来求解约束(12)，直到 $\beta$的下估计值和上估计值之间的差异达到某个指定阈值。
$V$ $q_1(x)$ $q_2(x)$ 不需要包含常数项。

论文中建议的多项式阶数如下：
$$deg(V)\ge deg({\phi }_1)$$ $$deg(q_1)+deg(s)\ge deg(V)$$ $$deg(q_2)+deg(d)\ge deg(n)-1$$