昨天学习了一下SOSTOOLS工具箱编程的逻辑和语法,今天尝试理解给出的DEMO,为解决科学问题打基础
1.可行性
给定一个多项式
p
(
x
)
=
2
x
1
4
+
2
x
1
3
x
2
−
x
1
2
x
2
2
+
5
x
2
4
p(x)=2x_1^4+2x_1^3x_2-x_1^2x_2^2+5x_2^4
p(x)=2x14+2x13x2−x12x22+5x24,确定
p
(
x
)
≥
0
p(x)\geq0
p(x)≥0 是否可行。
可以用 sosineq 来解决这个问题,也可以用 findsos 来解决这个问题
clc;clear all;
syms x1 x2
prog1 = sosprogram([x1;x2]);
p = 2*x1^4+2*x1^3*x2-x1^2*x2^2+5*x2^4
prog1 = sosineq(prog1,p);
options.solver = 'mosek';
[prog1,info] = sossolve(prog1,options);
clc;clear all;
syms x1 x2
prog1 = sosprogram([x1;x2]);
p = 2*x1^4+2*x1^3*x2-x1^2*x2^2+5*x2^4
findsos(p,'rational')
2. 搜索非线性系统的Lyapunov函数
对于一个非线性系统
[
x
1
˙
x
2
˙
x
3
˙
]
=
[
−
x
1
3
−
x
1
x
3
2
−
x
2
−
x
1
2
x
2
−
x
3
−
3
x
3
x
3
2
+
1
+
3
x
1
2
x
3
]
\left[ \begin{matrix} \dot{x_1} \\ \dot{x_2} \\ \dot{x_3} \end{matrix} \right] = \left[ \begin{matrix} -x_1 ^3-x_1 x_3^2\\ -x_2-x_1 ^2x_2 \\ -x_3-\frac{3x_3}{x_3^2+1}+3x_1^2x_3 \end{matrix} \right]
x1˙x2˙x3˙
=
−x13−x1x32−x2−x12x2−x3−x32+13x3+3x12x3
如何得到Lyapunov函数呢?
如果这个非线性系统是稳定的,那么它的Lyapunov函数
V
(
x
)
V(x)
V(x)必须满足:
V
−
ϵ
(
x
1
2
+
x
2
2
+
x
3
2
)
≥
0
V-\epsilon(x_1^2+x_2^2+x_3^2)\geq0
V−ϵ(x12+x22+x32)≥0
−
∂
V
∂
x
1
x
1
˙
−
∂
V
∂
x
2
x
2
˙
−
∂
V
∂
x
3
x
3
˙
≥
0
-\frac{{\displaystyle \partial }V}{ {\displaystyle \partial }x_1}\dot{x_1}- \frac{{\displaystyle \partial }V}{ {\displaystyle \partial }x_2}\dot{x_2}- \frac{{\displaystyle \partial }V}{ {\displaystyle \partial }x_3}\dot{x_3}\geq0
−∂x1∂Vx1˙−∂x2∂Vx2˙−∂x3∂Vx3˙≥0
下面要构建SOSP问题,选择
ϵ
=
1
\epsilon=1
ϵ=1:
V
−
(
x
1
2
+
x
2
2
+
x
3
2
)
≥
0
V-(x_1^2+x_2^2+x_3^2)\geq0
V−(x12+x22+x32)≥0
−
∂
V
∂
x
1
(
x
3
2
+
1
)
x
1
˙
−
∂
V
∂
x
2
(
x
3
2
+
1
)
x
2
˙
−
∂
V
∂
x
3
(
x
3
2
+
1
)
x
3
˙
≥
0
-\frac{{\displaystyle \partial }V}{ {\displaystyle \partial }x_1}(x_3^2+1)\dot{x_1}- \frac{{\displaystyle \partial }V}{ {\displaystyle \partial }x_2}(x_3^2+1)\dot{x_2}- \frac{{\displaystyle \partial }V}{ {\displaystyle \partial }x_3}(x_3^2+1)\dot{x_3}\geq0
−∂x1∂V(x32+1)x1˙−∂x2∂V(x32+1)x2˙−∂x3∂V(x32+1)x3˙≥0
clc;clear all;
pvar x1 x2 x3
f = [(-x1^3-x1*x3^2)*(x3^2+1);
(-x2-x1^2*x2)*(x3^2+1);
(-x3^3-4*x3+3*x1^2*x3^3+3*x1^2*x3)];
prog = sosprogram([x1;x2;x3]);
[prog,V] = sospolyvar(prog,[x1^2;x2^2;x3^2],'wscoeff');
prog = sosineq(prog,V-(x1^2+x2^2+x3^2));
expr = -diff(V,x1)*f(1)-diff(V,x2)*f(2)-diff(V,x3)*f(3);
prog = sosineq(prog,expr);
options.solver = 'mosek'
prog = sossolve(prog,options)
SOLVE = sosgetsol(prog,V)