[37] | 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); |
---|