1 | function F = robustify_lp_conic(F_xw,Zmodel,x,w) |
---|
2 | |
---|
3 | % Creates robustified version of the uncertain set of linear inequalities |
---|
4 | % s.t A(w)*x <= b(w) for all F(w) >= 0 where F(w) is a conic set, here |
---|
5 | % given in SeDuMi (and YALMIP internal) numerical format. |
---|
6 | % |
---|
7 | % Based on Robust Optimization - Methodology and Applications. A. Ben-Tal |
---|
8 | % and A. Nemerovskii. Mathematical Programming (Series B), 92:453-480, 2002 |
---|
9 | % |
---|
10 | % Note, there are some sign errors in the paper. |
---|
11 | |
---|
12 | if length(F_xw) == 0 |
---|
13 | F = []; |
---|
14 | return |
---|
15 | end |
---|
16 | |
---|
17 | X = sdpvar(F_xw); |
---|
18 | b = []; |
---|
19 | A = []; |
---|
20 | % Some pre-calc |
---|
21 | xw = [x;w]; |
---|
22 | xind = find(ismembc(getvariables(xw),getvariables(x))); |
---|
23 | wind = find(ismembc(getvariables(xw),getvariables(w))); |
---|
24 | [Qs,cs,fs,dummy,nonquadratic] = vecquaddecomp(X,xw); |
---|
25 | for i = 1:length(X) |
---|
26 | %[Q,c,f,dummy,nonquadratic] = quaddecomp(X(i),xw); |
---|
27 | Q = Qs{i}; |
---|
28 | c = cs{i}; |
---|
29 | f = fs{i}; |
---|
30 | if nonquadratic |
---|
31 | error('Constraints can be at most quadratic, with the linear term uncertain'); |
---|
32 | end |
---|
33 | Q_ww = Q(wind,wind); |
---|
34 | Q_xw = Q(xind,wind); |
---|
35 | Q_xx = Q(xind,xind); |
---|
36 | c_x = c(xind); |
---|
37 | c_w = c(wind); |
---|
38 | |
---|
39 | b = [b;f + c_w'*w]; |
---|
40 | A = [A;-(c_x + 2*Q_xw*w)']; |
---|
41 | end |
---|
42 | |
---|
43 | % Try to find variables that only have simple bound constraints. These |
---|
44 | % variables can explicitly be optimized and thus speed up the construction, |
---|
45 | % and allow a model with fewer variables. |
---|
46 | [Zmodel2,lower,upper] = find_simple_variable_bounds(Zmodel); |
---|
47 | |
---|
48 | % Partition the uncertain variables |
---|
49 | simple_w = find( ~isinf(lower) & ~isinf(upper)); |
---|
50 | general_w = find( isinf(lower) | isinf(upper)); |
---|
51 | simple_w = recover(simple_w); |
---|
52 | general_w = recover(general_w); |
---|
53 | |
---|
54 | |
---|
55 | % Linear uncertain constraint is (Bbetai*x + cdi) >= 0 for all w, or |
---|
56 | % (bi' + (Bi*w)')*x + (ci'*w + di). |
---|
57 | cd = b; |
---|
58 | Bbeta = -A; |
---|
59 | |
---|
60 | F = set([]); |
---|
61 | top = 1; |
---|
62 | |
---|
63 | % To speed up the construction, compute the ci vectors for all constraints |
---|
64 | % in one call ci_basis = [c1 c2 ...] |
---|
65 | ci_basis = basis(cd',w); |
---|
66 | |
---|
67 | for i = 1:length(b) |
---|
68 | cdi = cd(i); |
---|
69 | Bbetai = Bbeta(i,:); |
---|
70 | |
---|
71 | if (nnz(ci_basis(:,i))==0) & isa(Bbetai,'double') |
---|
72 | % This constraint row is constant |
---|
73 | F = F + set(Bbetai*x + cdi >= 0); |
---|
74 | else |
---|
75 | |
---|
76 | if isempty(general_w) |
---|
77 | |
---|
78 | ci = ci_basis(:,i); |
---|
79 | |
---|
80 | di = basis(cdi,0); |
---|
81 | if isa(Bbetai,'double') |
---|
82 | Bi = zeros(1,length(w)); |
---|
83 | else |
---|
84 | Bi = basis(Bbetai,w)'; |
---|
85 | end |
---|
86 | bi = basis(Bbeta(i,:),0)'; |
---|
87 | % Scale to -1,1 uncertainty |
---|
88 | T = diag((upper-lower))/2; |
---|
89 | e = (upper+lower)/2; |
---|
90 | if nnz(Bi) == 0 |
---|
91 | F = F + set(bi'*x + (di-e'*T*ci) - norm(T*ci,1) > 0); |
---|
92 | else |
---|
93 | non_zeroBirow = find(sum(abs(Bi'),2)); |
---|
94 | zeroBirow = find(sum(abs(Bi'),2) == 0); |
---|
95 | t = sdpvar(length(non_zeroBirow),1); |
---|
96 | F = F + set((bi'-e'*T*Bi')*x + (di-e'*T*ci) - norm(T(zeroBirow,:)*ci,1) - sum(t) >= 0) + set(-t < T(non_zeroBirow,:)*(ci+Bi'*x) < t); |
---|
97 | end |
---|
98 | |
---|
99 | else |
---|
100 | |
---|
101 | lhs1 = 0; |
---|
102 | lhs2 = 0; |
---|
103 | top = 1; |
---|
104 | |
---|
105 | if Zmodel.K.f > 0 |
---|
106 | zeta = sdpvar(Zmodel.K.f,1); |
---|
107 | F = F + set(zeta >= 0); |
---|
108 | lhs1 = lhs1 + Zmodel.F_struc(top:top + Zmodel.K.f-1,2:end)'*zeta; |
---|
109 | lhs2 = lhs2 - Zmodel.F_struc(top:top + Zmodel.K.f-1,1)'*zeta; |
---|
110 | top = top + Zmodel.K.f; |
---|
111 | end |
---|
112 | |
---|
113 | if Zmodel.K.l > 0 |
---|
114 | zeta = sdpvar(Zmodel.K.l,1); |
---|
115 | F = F + set(zeta >= 0); |
---|
116 | lhs1 = lhs1 + Zmodel.F_struc(top:top + Zmodel.K.l-1,2:end)'*zeta; |
---|
117 | lhs2 = lhs2 - Zmodel.F_struc(top:top + Zmodel.K.l-1,1)'*zeta; |
---|
118 | top = top + Zmodel.K.l; |
---|
119 | end |
---|
120 | |
---|
121 | if Zmodel.K.q(1) > 0 |
---|
122 | for j = 1:length(Zmodel.K.q) |
---|
123 | zeta = sdpvar(Zmodel.K.q(j),1); |
---|
124 | F = F + set(cone(zeta)); |
---|
125 | lhs1 = lhs1 + Zmodel.F_struc(top:top + Zmodel.K.q(j)-1,2:end)'*zeta(:); |
---|
126 | lhs2 = lhs2 - Zmodel.F_struc(top:top + Zmodel.K.q(j)-1,1)'*zeta(:); |
---|
127 | top = top + Zmodel.K.q(j); |
---|
128 | end |
---|
129 | end |
---|
130 | |
---|
131 | if Zmodel.K.s(1) > 0 |
---|
132 | for j = 1:length(Zmodel.K.s) |
---|
133 | zeta = sdpvar(Zmodel.K.s(j)); |
---|
134 | F = F + set(zeta >= 0); |
---|
135 | lhs1 = lhs1 + Zmodel.F_struc(top:top + Zmodel.K.s(j)^2-1,2:end)'*zeta(:); |
---|
136 | lhs2 = lhs2 - Zmodel.F_struc(top:top + Zmodel.K.s(j)^2-1,1)'*zeta(:); |
---|
137 | top = top + Zmodel.K.s(j)^2; |
---|
138 | end |
---|
139 | end |
---|
140 | |
---|
141 | % if isempty(simple_w) |
---|
142 | ci = basis(cd(i),w); |
---|
143 | di = basis(cd(i),0); |
---|
144 | Bi = basis(Bbeta(i,:),w)'; |
---|
145 | bi = basis(Bbeta(i,:),0)'; |
---|
146 | F = F + set(lhs1 == Bi'*x + ci); |
---|
147 | F = F + set(lhs2 >= - (bi'*x + di)); |
---|
148 | end |
---|
149 | end |
---|
150 | end |
---|
151 | |
---|
152 | |
---|
153 | function b = basis(p,w) |
---|
154 | |
---|
155 | if isequal(w,0) |
---|
156 | b = getbasematrix(p,0); |
---|
157 | else |
---|
158 | n = length(w); |
---|
159 | if isequal(getbase(w),[zeros(n,1) eye(n)]) |
---|
160 | b = []; |
---|
161 | lmi_variables = getvariables(w); |
---|
162 | for i = 1:length(w) |
---|
163 | b = [b ; getbasematrix(p,lmi_variables(i))]; |
---|
164 | end |
---|
165 | else |
---|
166 | b = []; |
---|
167 | for i = 1:length(w) |
---|
168 | b = [b ; getbasematrix(p,getvariables(w(i)))]; |
---|
169 | end |
---|
170 | end |
---|
171 | end |
---|
172 | b = full(b); |
---|
173 | |
---|
174 | |
---|
175 | function [p,lower,upper] = find_simple_variable_bounds(p); |
---|
176 | |
---|
177 | if any(p.K.q > 0) | any(p.K.s > 0) |
---|
178 | lower = -inf(length(p.c),1); |
---|
179 | upper = inf(length(p.c),1); |
---|
180 | else |
---|
181 | [lower,upper,used_rows] = findulb(p.F_struc,p.K); |
---|
182 | used_rows = used_rows(~any(full(p.F_struc(used_rows,1+find(p.variabletype~=0))),2)); |
---|
183 | if ~isempty(used_rows) |
---|
184 | p_temp = p; |
---|
185 | p_temp.F_struc(p.K.f+used_rows,:)=[]; |
---|
186 | p_temp.K.l = p.K.l - length(used_rows); |
---|
187 | |
---|
188 | if size(p.F_struc,1) > 0 |
---|
189 | % These variables are still used in some other constraints |
---|
190 | still_used = find(sum(abs(p_temp.F_struc(:,2:end)),1) > 0); |
---|
191 | if ~isempty(still_used) |
---|
192 | % we have to keep these variables |
---|
193 | lower(still_used) = -inf; |
---|
194 | upper(still_used) = inf; |
---|
195 | % They are used here |
---|
196 | keep_rows = find(sum(abs(p.F_struc(:,1+still_used)),2) > 0); |
---|
197 | if any(keep_rows>(p.K.l + p.K.f)) |
---|
198 | % error('Tell johan to fix the SDP case in find_simple_variable_bounds!') |
---|
199 | end |
---|
200 | used_rows = used_rows + p.K.f; |
---|
201 | used_rows = setdiff(used_rows,keep_rows); |
---|
202 | p.F_struc(used_rows,:)=[]; |
---|
203 | p.K.l = p.K.l - nnz(used_rows>p.K.f); |
---|
204 | p.K.f = p.K.f - nnz(used_rows<=p.K.f); |
---|
205 | else |
---|
206 | p = p_temp; |
---|
207 | end |
---|
208 | else |
---|
209 | p = p_temp; |
---|
210 | end |
---|
211 | end |
---|
212 | end |
---|
213 | |
---|
214 | |
---|
215 | |
---|
216 | |
---|
217 | |
---|
218 | |
---|
219 | |
---|
220 | |
---|
221 | |
---|
222 | |
---|
223 | |
---|