1 | function varargout = dilate(F,w) |
---|
2 | % DILATE Derives a matrix dilation |
---|
3 | % |
---|
4 | % [G,H,M] = DILATE(X,w) where X is a symmetric variable derives the |
---|
5 | % decomposition and orthogonal complement used in a matrix dilation. |
---|
6 | % |
---|
7 | % X is decomposed as M(w)ŽG(x)M(w), and H(w) is an orthogonal complement to |
---|
8 | % M(w), with affine dependence in w. These matrices can be used to apply |
---|
9 | % the matrix dilation lemma to obtain a constraint affine in w |
---|
10 | % |
---|
11 | % M(w)ŽG(x)M(w) > 0 <==> existence of W s.t with G + W*H' + H*W' > 0 |
---|
12 | % |
---|
13 | % F = DILATE(F,w) where F is a SET object is used to construct the dilated |
---|
14 | % uncertain constraint, i.e an SDP constraint where polynomial dependence |
---|
15 | % w.r.t uncertain variables in all SDP constraints in F are converted |
---|
16 | % (conservatively) to affine dependence using the matrix dilation approach. |
---|
17 | % |
---|
18 | % See Yasuaki OISHI, A Region-Dividing Approach to Robust Semidefinite |
---|
19 | % Programming and Its Error Bound,DEPARTMENT OF MATHEMATICAL INFORMATICS |
---|
20 | % GRADUATE SCHOOL OF INFORMATION SCIENCE AND TECHNOLOGY THE UNIVERSITY OF |
---|
21 | % TOKYO BUNKYO-KU, TOKYO 113-8656, JAPAN , February 2006 |
---|
22 | % |
---|
23 | % See also ROBUSTIFY, SOLVEROBUST, UNCERTAIN |
---|
24 | |
---|
25 | % Author Johan Löfberg |
---|
26 | % $Id: dilate.m,v 1.2 2006/10/25 09:18:47 joloef Exp $ |
---|
27 | |
---|
28 | if isa(F,'sdpvar') |
---|
29 | [G,H,M] = matrix_dilate(F,w) |
---|
30 | varargout{1} = G; |
---|
31 | varargout{2} = H; |
---|
32 | varargout{3} = M; |
---|
33 | elseif isa(F,'lmi') |
---|
34 | Fnew = []; |
---|
35 | if nargin == 1 |
---|
36 | w = []; |
---|
37 | else |
---|
38 | w = w(:); |
---|
39 | end |
---|
40 | unc_declarations = is(F,'uncertain'); |
---|
41 | if any(unc_declarations) |
---|
42 | w = [w;recover(getvariables(sdpvar(F(find(unc_declarations)))))]; |
---|
43 | end |
---|
44 | for i = 1:length(F) |
---|
45 | if (is(F(i),'lmi') | ((length(sdpvar(F(i))) == 1) & is(F(i),'elementwise'))) & max(degree(sdpvar(F(i)),w))>0 |
---|
46 | [G,H,M] = matrix_dilate(sdpvar(F(i)),w); |
---|
47 | W = sdpvar(size(G,1),size(H,2)); |
---|
48 | Fnew = Fnew + set(G + W*H' + H*W' >= 0); |
---|
49 | else |
---|
50 | Fnew = Fnew + F(i); |
---|
51 | end |
---|
52 | end |
---|
53 | varargout{1} = Fnew; |
---|
54 | end |
---|
55 | |
---|
56 | |
---|
57 | function [G,H,M] = matrix_dilate(F,w) |
---|
58 | % Given an SDP F(x,w) with F polynomial in w, DILATE rewrites the problem |
---|
59 | % to F(x,w)=M(w)'G(x)M(w), and derives the orhogonal complements H(w) to |
---|
60 | % M(w), to be used in the dilated constraint G + W*H' + H*W' |
---|
61 | |
---|
62 | allvariables = getvariables(F); |
---|
63 | wvariables = getvariables(w); |
---|
64 | Fbasis = getbase(F); |
---|
65 | |
---|
66 | % Make sure it is ordered according to internal index |
---|
67 | w = recover(wvariables); |
---|
68 | |
---|
69 | % Degrees w.r.t the uncertain variables |
---|
70 | d = degree(recover(allvariables),w); |
---|
71 | |
---|
72 | % robustify code by removing unused variables |
---|
73 | w = w(find(d)); |
---|
74 | wvariables = wvariables(find(d)); |
---|
75 | |
---|
76 | % Degrees w.r.t the uncertain variables |
---|
77 | monomtable = yalmip('monomtable'); |
---|
78 | d = max(sum(monomtable(allvariables,wvariables),2)); |
---|
79 | n = size(F,1); |
---|
80 | |
---|
81 | % Sufficiently many monomials |
---|
82 | v = monolist(w,max(d)); |
---|
83 | % Some numerical format on these variables |
---|
84 | for i = 2:length(v) |
---|
85 | vvariables(i,1) = getvariables(v(i)); |
---|
86 | end |
---|
87 | monomtable = yalmip('monomtable'); |
---|
88 | wmonoms = [zeros(1,length(w));monomtable(vvariables(2:end),wvariables)]; |
---|
89 | |
---|
90 | % Find the linear certain term in the matrix |
---|
91 | linear_indicies = find(sum(monomtable(allvariables,wvariables),2) ==0); |
---|
92 | F0 = reshape(Fbasis(:,1) + Fbasis(:,1+linear_indicies)*recover(allvariables(linear_indicies)),n,n); |
---|
93 | |
---|
94 | % now find the matrix that multiplies with each monomial in M |
---|
95 | Fi = []; |
---|
96 | Fmonoms = monomtable(allvariables,wvariables); |
---|
97 | for i = 2:length(v) |
---|
98 | vmonoms = monomtable(vvariables(i),wvariables); |
---|
99 | index = findrows(Fmonoms,vmonoms); |
---|
100 | if isempty(index) |
---|
101 | Fi = [Fi zeros(n)]; |
---|
102 | else |
---|
103 | temp = monomtable(allvariables(index),:); |
---|
104 | temp(wvariables) = 0; |
---|
105 | base = reshape(Fbasis(:,1+index),n,n); |
---|
106 | if nnz(temp) == 0 |
---|
107 | Fi = [Fi base]; |
---|
108 | else |
---|
109 | Fi = [Fi base*recover(allvariables(find(temp)))]; |
---|
110 | end |
---|
111 | end |
---|
112 | end |
---|
113 | G = [F0 Fi/2;Fi'/2 zeros(n*(length(v)-1))]; |
---|
114 | |
---|
115 | % The outer factor |
---|
116 | M = kron(v,eye(n)); |
---|
117 | |
---|
118 | % Now create an orthogonal complement to M |
---|
119 | ii = []; |
---|
120 | jj = []; |
---|
121 | ss = []; |
---|
122 | for i = 2:length(v) |
---|
123 | monom = wmonoms(i,:); |
---|
124 | [dummy,index] = max(monom); |
---|
125 | monomnew = monom; |
---|
126 | monomnew(index) = monomnew(index) - 1; |
---|
127 | ii = [ii i]; |
---|
128 | jj = [jj i-1]; |
---|
129 | ss = [ss 1]; |
---|
130 | ii = [ii findrows(wmonoms,monomnew)]; |
---|
131 | jj = [jj i-1]; |
---|
132 | ss = [ss -recover(wvariables(index))]; |
---|
133 | end |
---|
134 | H = kron(sparse(ii,jj,ss),eye(n)); |
---|