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