[37] | 1 | function pass = akerblad |
---|
| 2 | |
---|
| 3 | G=tf(1,[1 1 1]); |
---|
| 4 | Gd=c2d(G,0.1); |
---|
| 5 | [A,B]=ssdata(Gd); |
---|
| 6 | [n,m]=size(B); |
---|
| 7 | N=100; |
---|
| 8 | |
---|
| 9 | Ci=[1 0;0 1; 0 0]; |
---|
| 10 | Di=[0; 0; 1]; |
---|
| 11 | p=3; |
---|
| 12 | q=4; |
---|
| 13 | Ei=[1 3;-1 -3;0 0; 0 0]; |
---|
| 14 | Fi=[0;0;1; -1]; |
---|
| 15 | bi=[10;10;1;1]; |
---|
| 16 | P=[1 0;0 1;]; |
---|
| 17 | W=[1 0;0 1;]; |
---|
| 18 | E=[];F=[];C=[];D=[]; |
---|
| 19 | ci=[1;1;zeros(p,1)]; |
---|
| 20 | x0=[2;2]; |
---|
| 21 | b=[]; |
---|
| 22 | c=[]; |
---|
| 23 | for k=1:N, |
---|
| 24 | E(:,:,k)=Ei; |
---|
| 25 | F(:,:,k)=Fi; |
---|
| 26 | C(:,:,k)=Ci; |
---|
| 27 | D(:,:,k)=Di; |
---|
| 28 | b=[b;bi]; |
---|
| 29 | c=[c;ci]; |
---|
| 30 | end |
---|
| 31 | c=[c;1;1;zeros(n,1);1;zeros(n,1)]; |
---|
| 32 | d=sparse([],[],[],(N+1)*n,1); |
---|
| 33 | d(1:n)=x0; |
---|
| 34 | H=sparse([],[],[]); |
---|
| 35 | H(1:n,1+1:n+1)=eye(n); |
---|
| 36 | M=sparse([],[],[],N*q,N*(n+m+1)+1+n); |
---|
| 37 | L=sparse([],[],[]); |
---|
| 38 | |
---|
| 39 | for k=1:N, |
---|
| 40 | M(k*q+1-q:k*q,k*(n+m+1)-n-m:k*(n+m+1))=[zeros(q,1) E(:,:,k) F(:,:,k)]; |
---|
| 41 | L(k*(p+2)+1-p-2:k*(p+2),k*(n+m+1)-n-m:k*(n+m+1))=[-1 zeros(1,n+m);1 zeros(1,n+m);zeros(p,1) -2*C(:,:,k) -2*D(:,:,k)]; |
---|
| 42 | H(k*n+1:k*n+n,k*(n+m+1)+1-n-m:k*(1+n+m)+n+1)=[-A -B zeros(n,1) ... |
---|
| 43 | eye(n)]; |
---|
| 44 | end |
---|
| 45 | L(k*(p+2)+1:(k)*(p+2)+2+n+1+n,k*(n+m+1)+1:k*(n+m+1)+n+1)=[-1 ... |
---|
| 46 | zeros(1,n);1 zeros(1,n);zeros(n,1) -2*P;zeros(1,n+ ... |
---|
| 47 | 1);zeros(n,1) -W]; |
---|
| 48 | f=sparse([],[],[],N*(n+m+1)+n+1,1); |
---|
| 49 | f(1:n+m+1:N*(n+m+1)+n+1)=1; |
---|
| 50 | [yy,ss,rr,sol]=my_socp2(f,d,c,b,L,M,H,N); |
---|
| 51 | |
---|
| 52 | norm(yy) |
---|
| 53 | norm(ss) |
---|
| 54 | |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | function [y,s,r,sol]=my_socp2(f,b,c,d,A,B,H,N) |
---|
| 64 | %[t,x,s,r,gamma,lambda]=my_socp(f,c,d,A,B) |
---|
| 65 | %Solves the primal and dual LP+SOCP on the form |
---|
| 66 | %min f'y | max -d'gamma-c'lambda-b'rho |
---|
| 67 | %st s=c-Ay | st -A'lambda-B'gamma-H'rho=f |
---|
| 68 | % r=d-By | 0<gamma, ||lambda_1||<lambda_0 |
---|
| 69 | % Hy=b | |
---|
| 70 | % 0<r, ||s_1||<s_0 | |
---|
| 71 | [n,m]=size(A); |
---|
| 72 | [p,m]=size(B); |
---|
| 73 | [q,m]=size(H); |
---|
| 74 | |
---|
| 75 | tic |
---|
| 76 | yalmip('clear') |
---|
| 77 | |
---|
| 78 | y=sdpvar(m,1); |
---|
| 79 | s = c-A*y; |
---|
| 80 | r = d-B*y; |
---|
| 81 | F=lmi('0<r'); |
---|
| 82 | F=addlmi(F,'H*y=b'); |
---|
| 83 | for i=1:N |
---|
| 84 | F=addlmi(F,'||s(i*(2+3)-3:i*5)||<s(i*5-4)'); |
---|
| 85 | end |
---|
| 86 | F=addlmi(F,'||s((i+1)*(2+3)-3:(i+1)*5-1)||<s((i+1)*5-4)'); |
---|
| 87 | F=addlmi(F,'||s(n-1:n)||<s(n-2)'); |
---|
| 88 | toc |
---|
| 89 | sol=solvesdp(F,[],f'*y,sdpsettings('sedumi.beta',0.8)) |
---|
| 90 | |
---|
| 91 | y=double(y); s=double(s); r=double(r); |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | |
---|