1 | function test_sos_peyrl |
---|
2 | |
---|
3 | sdpvar x1 x2 x3 x4; |
---|
4 | x = [x1;x2;x3;x4]; |
---|
5 | n = length(x); |
---|
6 | |
---|
7 | % System matrices: |
---|
8 | Ao = [0 1 0 0; |
---|
9 | 0 0 1 0; |
---|
10 | 0 0 0 1; |
---|
11 | -1 0 -2 0]; |
---|
12 | |
---|
13 | bo = [0 0 0 1]'; |
---|
14 | |
---|
15 | fo = [-1 -1 -2 -1]; |
---|
16 | |
---|
17 | % Transformation-Matrix: |
---|
18 | Pinv = [1 -1 -1 0; |
---|
19 | -1 1 -1 0; |
---|
20 | 0 -1 1 0; |
---|
21 | 0 2 0 -1]; |
---|
22 | |
---|
23 | A = Ao;%inv(Pinv)*Ao*Pinv; |
---|
24 | b = bo;%inv(Pinv)*bo; |
---|
25 | |
---|
26 | f = fo;%fo*Pinv; |
---|
27 | |
---|
28 | % System equations: |
---|
29 | f1 = A*x+b*f*x; |
---|
30 | f2 = A*x+b; |
---|
31 | f3 = A*x-b; |
---|
32 | |
---|
33 | % Cell boundaries: |
---|
34 | g11 = f*x+1; |
---|
35 | g12 = -f*x+1; |
---|
36 | g21 = f*x-1; |
---|
37 | g31 = -f*x-1; |
---|
38 | h120 = f*x-1; |
---|
39 | h130 = f*x+1; |
---|
40 | |
---|
41 | |
---|
42 | % The Lyapunov functions Vi(x): |
---|
43 | disp('Constructing Lyapunov functions') |
---|
44 | |
---|
45 | N = 4; |
---|
46 | |
---|
47 | % The Lyapunov functions Vi(x): |
---|
48 | lm = monolist(x,N,2); |
---|
49 | p1 = sdpvar(size(lm,1),1); define(p1); |
---|
50 | V1 = p1'*lm; |
---|
51 | |
---|
52 | lm = monolist(x,N); |
---|
53 | p2 = sdpvar(size(lm,1),1); define(p2); |
---|
54 | V2 = p2'*lm; |
---|
55 | p3 = sdpvar(size(lm,1),1); define(p3); |
---|
56 | V3 = p3'*lm; |
---|
57 | |
---|
58 | |
---|
59 | sdpvar eps1 eps2 eps3 eps4 eps5 eps6 |
---|
60 | pd1 = eps1*sum(x.^N); |
---|
61 | pd2 = eps2*sum(x.^N); |
---|
62 | pd3 = eps3*sum(x.^N); |
---|
63 | pd4 = eps4*sum(x.^N); |
---|
64 | pd5 = eps5*sum(x.^N); |
---|
65 | pd6 = eps6*sum(x.^N); |
---|
66 | |
---|
67 | % Create the aijk, bijk SOS-polys and cij: |
---|
68 | lm = monolist(x,N/2-1); |
---|
69 | A11 = sdpvar(size(lm,1)); (A11); a11 = lm'*A11*lm; |
---|
70 | A12 = sdpvar(size(lm,1)); (A12); a12 = lm'*A12*lm; |
---|
71 | A21 = sdpvar(size(lm,1)); (A21); a21 = lm'*A21*lm; |
---|
72 | A31 = sdpvar(size(lm,1)); (A31); a31 = lm'*A31*lm; |
---|
73 | |
---|
74 | B11 = sdpvar(size(lm,1)); (B11); b11 = lm'*B11*lm; |
---|
75 | B12 = sdpvar(size(lm,1)); (B12); b12 = lm'*B12*lm; |
---|
76 | B21 = sdpvar(size(lm,1)); (B21); b21 = lm'*B21*lm; |
---|
77 | B31 = sdpvar(size(lm,1)); (B31); b31 = lm'*B31*lm; |
---|
78 | |
---|
79 | lm = monolist(x,N-1); |
---|
80 | C12 = sdpvar(size(lm,1)); (C12); c12 = lm'*C12*lm; |
---|
81 | C13 = sdpvar(size(lm,1)); (C13); c13 = lm'*C13*lm; |
---|
82 | |
---|
83 | % Constraints: |
---|
84 | F = set(sos(V1 - a11*g11 - a12*g12)); |
---|
85 | F = F + set(sos(V2 - a21*g21 )); |
---|
86 | F = F + set(sos(V3 - a31*g31 )); |
---|
87 | |
---|
88 | F = F + set(sos(-jacobian(V1,x)*f1 - b11*g11 - b12*g12 )); |
---|
89 | F = F + set(sos(-jacobian(V2,x)*f2 - b21*g21) ); |
---|
90 | F = F + set(sos(-jacobian(V3,x)*f3 - b31*g31) ); |
---|
91 | |
---|
92 | F = F + set(coefficients(V1+c12*h120-V2,x)==0); |
---|
93 | F = F + set(coefficients(V1+c13*h130-V3,x)==0); |
---|
94 | |
---|
95 | |
---|
96 | F = F + set(eps1>0.1) + set(eps2>0.1) + set(eps3>0.1); |
---|
97 | F = F + set(eps4>0.1) + set(eps5>0.1) + set(eps6>0.1); |
---|
98 | |
---|
99 | F = F + set(A11>0) + set(A12>0) + set(A21>0) + set(A31>0) + set(B11>0) + ... |
---|
100 | set(B12>0) + set(B21>0) + set(B31>0); |
---|
101 | |
---|
102 | % Call solver: |
---|
103 | parametric=recover(setdiff(depends(F),depends(x))); |
---|
104 | [sol,v,Q,residual,model]=solvesos(F,[],sdpsettings('sos.post',1,'sedumi.free',0),parametric); |
---|
105 | |
---|
106 | mbg_asserttrue(sol.problem==0 | sol.problem == 4); |
---|
107 | mbg_asserttrue(norm(residual) < 1e-5); |
---|