1 | function [pres,dres] = checklmi(F) |
---|
2 | % checklmi(F) Displays/calculates constraint residuals on set-object F |
---|
3 | % |
---|
4 | % [pres,dres] = checkset(F) |
---|
5 | % |
---|
6 | % pres : Primal constraint residuals |
---|
7 | % dres : Dual constraint residuals |
---|
8 | % |
---|
9 | % If no output argument is supplied, tabulated results are displayed |
---|
10 | % |
---|
11 | % Primal constraint residuals are calculated as: |
---|
12 | % |
---|
13 | % Semidefinite constraint F(x)>0 : min(eig(F)) |
---|
14 | % Element-wise constraint F(x)>0 : min(min(F)) |
---|
15 | % Equality constraint F==0 : -max(max(abs(F))) |
---|
16 | % Second order cone t>||x|| : t-||x|| |
---|
17 | % Integrality constraint on x : max(abs(x-round(x))) |
---|
18 | % Rank constraint rank(X) < r : r-rank(X) |
---|
19 | % Sum-of-square constraint : Minus value of largest (absolute value) coefficient |
---|
20 | % in the polynomial p-v'*v |
---|
21 | % |
---|
22 | % Dual constraints are evaluated similarily. |
---|
23 | % |
---|
24 | % See also SET, SOLVESDP, SOLVESOS, SOSD, DUAL |
---|
25 | |
---|
26 | % Author Johan Löfberg |
---|
27 | % $Id: checkset.m,v 1.16 2006/05/11 10:49:13 joloef Exp $ |
---|
28 | |
---|
29 | % Check if solution avaliable |
---|
30 | currsol = evalin('caller','yalmip(''getsolution'')'); |
---|
31 | if isempty(currsol) |
---|
32 | disp('No solution available.') |
---|
33 | return |
---|
34 | end |
---|
35 | |
---|
36 | nlmi = size(F.clauses,2); |
---|
37 | spaces = [' ']; |
---|
38 | if (nlmi == 0) |
---|
39 | if nargout == 0 |
---|
40 | disp('empty LMI') |
---|
41 | else |
---|
42 | pres = []; |
---|
43 | dres = []; |
---|
44 | end |
---|
45 | return |
---|
46 | end |
---|
47 | |
---|
48 | lmiinfo{1} = 'Matrix inequality'; |
---|
49 | lmiinfo{2} = 'Elementwise inequality'; |
---|
50 | lmiinfo{3} = 'Equality constraint'; |
---|
51 | lmiinfo{4} = 'Second order cone constraint'; |
---|
52 | lmiinfo{5} = 'Rotated Lorentz constraint'; |
---|
53 | lmiinfo{7} = 'Integer constraint'; |
---|
54 | lmiinfo{8} = 'Binary constraint'; |
---|
55 | lmiinfo{9} = 'KYP constraint'; |
---|
56 | lmiinfo{10} = 'Eigenvalue constraint'; |
---|
57 | lmiinfo{11} = 'SOS constraint'; |
---|
58 | |
---|
59 | header = {'ID','Constraint','Type','Primal residual','Dual residual','Tag'}; |
---|
60 | |
---|
61 | if nargout==0 |
---|
62 | disp(' '); |
---|
63 | end |
---|
64 | |
---|
65 | % Checkset is very slow for multiple SOS |
---|
66 | % The reason is that REPLACE has to be called |
---|
67 | % for every SOS. Instead, precalc on one vector |
---|
68 | p=[]; |
---|
69 | ParVar=[]; |
---|
70 | soscount=1; |
---|
71 | for j = 1:nlmi |
---|
72 | if F.clauses{j}.type==11 |
---|
73 | pi = F.clauses{j}.data; |
---|
74 | [v,ParVari] = sosd(pi); |
---|
75 | if isempty(v) |
---|
76 | p=[p;0]; |
---|
77 | else |
---|
78 | p=[p;pi]; |
---|
79 | ParVar=unique([ParVar(:);ParVari(:)]); |
---|
80 | end |
---|
81 | end |
---|
82 | end |
---|
83 | if ~isempty(ParVar) |
---|
84 | ParVar = recover(ParVar); |
---|
85 | p = replace(p,ParVar,double(ParVar)); |
---|
86 | end |
---|
87 | |
---|
88 | for j = 1:nlmi |
---|
89 | constraint_type = F.clauses{j}.type; |
---|
90 | if constraint_type~=11 |
---|
91 | F0 = double(F.clauses{j}.data); |
---|
92 | end |
---|
93 | if (constraint_type~=11) & any(isnan(F0(:))) |
---|
94 | res(j,1) = NaN; |
---|
95 | else |
---|
96 | switch F.clauses{j}.type |
---|
97 | case {1,9} |
---|
98 | res(j,1) = full(min(real(eig(F0)))); |
---|
99 | case 2 |
---|
100 | res(j,1) = full(min(min(F0))); |
---|
101 | case 3 |
---|
102 | res(j,1) = -full(max(max(abs(F0)))); |
---|
103 | case 4 |
---|
104 | res(j,1) = full(F0(1)-norm(F0(2:end))); |
---|
105 | case 5 |
---|
106 | res(j,1) = full(2*F0(1)*F0(2)-norm(F0(3:end))^2); |
---|
107 | case {7,8} |
---|
108 | res(j,1) = full(max(max(abs(F0-round(F0))))); |
---|
109 | case 11 |
---|
110 | if 0 |
---|
111 | p = F.clauses{j}.data; |
---|
112 | [v,ParVar] = sosd(p); |
---|
113 | if ~isempty(ParVar) |
---|
114 | ParVar = recover(ParVar); |
---|
115 | p = replace(p,ParVar,double(ParVar)); |
---|
116 | end |
---|
117 | if isempty(v) |
---|
118 | res(j,1)=nan; |
---|
119 | else |
---|
120 | h = p-v'*v; |
---|
121 | res(j,1) = full(max(max(abs(getbase(h))))); |
---|
122 | end |
---|
123 | else |
---|
124 | %p = F.clauses{j}.data; |
---|
125 | [v,ParVar] = sosd(F.clauses{j}.data); |
---|
126 | if isempty(v) |
---|
127 | res(j,1)=nan; |
---|
128 | else |
---|
129 | h = p(soscount)-v'*v; |
---|
130 | res(j,1) = full(max(max(abs(getbase(h))))); |
---|
131 | end |
---|
132 | soscount=soscount+1; |
---|
133 | end |
---|
134 | |
---|
135 | otherwise |
---|
136 | res(j,1) = nan; |
---|
137 | end |
---|
138 | end |
---|
139 | |
---|
140 | % Get the internal index |
---|
141 | LMIid = F.LMIid(j); |
---|
142 | dual = yalmip('dual',LMIid); |
---|
143 | if isempty(dual) | any(isnan(dual)) |
---|
144 | resdual(j,1) = NaN; |
---|
145 | else |
---|
146 | switch F.clauses{j}.type |
---|
147 | case {1,9} |
---|
148 | resdual(j,1) = min(eig(dual)); |
---|
149 | case 2 |
---|
150 | resdual(j,1) = min(min(dual)); |
---|
151 | case 3 |
---|
152 | resdual(j,1) = -max(max(abs(dual))); |
---|
153 | case 4 |
---|
154 | resdual(j,1) = dual(1)-norm(dual(2:end)); |
---|
155 | case 5 |
---|
156 | resdual(j,1) = 2*dual(1)*dual(2)-norm(dual(3:end))^2; |
---|
157 | case 7 |
---|
158 | resdual(j,1) = nan; |
---|
159 | otherwise |
---|
160 | gap = nan; |
---|
161 | end |
---|
162 | end |
---|
163 | % if isempty(dual) | any(isnan(dual)) | any(isnan(F0)) |
---|
164 | % gap = NaN; |
---|
165 | % else |
---|
166 | % switch F.clauses{j}.type |
---|
167 | % case {1,9} |
---|
168 | % gap = trace(F0*dual); |
---|
169 | % case {2,3} |
---|
170 | % gap = F0(:)'*dual(:); |
---|
171 | % case 4 |
---|
172 | % gap = nan; |
---|
173 | % case 5 |
---|
174 | % gap = nan; |
---|
175 | % case 7 |
---|
176 | % gap = nan; |
---|
177 | % otherwise |
---|
178 | % gap = nan; |
---|
179 | % end |
---|
180 | % end |
---|
181 | |
---|
182 | if nargout==0 |
---|
183 | data{j,1} = ['#' num2str(j)]; |
---|
184 | data{j,2} = F.clauses{j}.symbolic; |
---|
185 | data{j,3} = lmiinfo{F.clauses{j}.type}; |
---|
186 | data{j,4} = res(j,1); |
---|
187 | data{j,5} = resdual(j,1); |
---|
188 | % data{j,6} = gap; |
---|
189 | data{j,6} = F.clauses{j}.handle; |
---|
190 | |
---|
191 | |
---|
192 | if ~islinear(F.clauses{j}.data) |
---|
193 | if is(F.clauses{j}.data,'sigmonial') |
---|
194 | classification = 'sigmonial'; |
---|
195 | elseif is(F.clauses{j}.data,'bilinear') |
---|
196 | classification = 'bilinear'; |
---|
197 | elseif is(F.clauses{j}.data,'quadratic') |
---|
198 | classification = 'quadratic'; |
---|
199 | else |
---|
200 | classification = 'polynomial'; |
---|
201 | end |
---|
202 | data{j,3} = [data{j,3} ' (' classification ')']; |
---|
203 | end |
---|
204 | end |
---|
205 | end |
---|
206 | |
---|
207 | if nargout>0 |
---|
208 | pres = res; |
---|
209 | dres = resdual; |
---|
210 | else |
---|
211 | if length([data{:,6}])==0 |
---|
212 | header = {header{:,1:5}}; |
---|
213 | temp = {data{:,1:5}}; |
---|
214 | data = reshape(temp,length(temp)/5,5); |
---|
215 | end |
---|
216 | |
---|
217 | if all(isnan(resdual)) |
---|
218 | header = {header{:,[1 2 3 4]}}; |
---|
219 | temp = {data{:,[1 2 3 4]}}; |
---|
220 | data = reshape(temp,length(temp)/4,4); |
---|
221 | end |
---|
222 | table('',header,data) |
---|
223 | disp(' '); |
---|
224 | end |
---|