[37] | 1 | function model = savemptmodel(model,Pfinal,Pn,Fi,Gi,details); |
---|
| 2 | |
---|
| 3 | if length(Fi)>0 |
---|
| 4 | if length(model) == 0 |
---|
| 5 | model{1} = fakemptmodel(Pfinal, Pn, Fi, Gi, details.Ai, details.Bi, details.Ci); |
---|
| 6 | [H,K] = double(Pfinal); |
---|
| 7 | model{1}.epicost = []; |
---|
| 8 | model{1}.convex = 1; |
---|
| 9 | else |
---|
| 10 | anyqp = nnz([details.Ai{:}])>0; |
---|
| 11 | for i = 1:length(model) |
---|
| 12 | anyqp = anyqp | nnz([model{i}.Ai{:}])>0; |
---|
| 13 | if anyqp |
---|
| 14 | break |
---|
| 15 | end |
---|
| 16 | end |
---|
| 17 | if anyqp |
---|
| 18 | model = savemptmodelqp(model,Pfinal,Pn,Fi,Gi,details); |
---|
| 19 | return |
---|
| 20 | else |
---|
| 21 | newmodel = fakemptmodel(Pfinal, Pn, Fi, Gi, details.Ai, details.Bi, details.Ci); |
---|
| 22 | newmodel.epicost = []; |
---|
| 23 | newmodel.convex = 1; |
---|
| 24 | |
---|
| 25 | replace = zeros(length(model),1); |
---|
| 26 | discard = 0; |
---|
| 27 | for i = 1:length(model) |
---|
| 28 | |
---|
| 29 | if Pfinal == model{i}.Pfinal |
---|
| 30 | |
---|
| 31 | if isempty(newmodel.epicost) |
---|
| 32 | B = reshape([details.Bi{:}]',size(details.Bi{1},2),[])'; |
---|
| 33 | c = reshape([details.Ci{:}]',[],1); |
---|
| 34 | [H,K] = double(Pfinal); |
---|
| 35 | newmodel.epicost = generate_epicost(H,K,B,c); |
---|
| 36 | end |
---|
| 37 | if isempty(model{i}.epicost) |
---|
| 38 | B = reshape([model{i}.Bi{:}]',size(model{i}.Bi{1},2),[])'; |
---|
| 39 | c = reshape([model{i}.Ci{:}]',[],1); |
---|
| 40 | [H,K] = double(model{i}.Pfinal); |
---|
| 41 | model{i}.epicost = generate_epicost(H,K,B,c); |
---|
| 42 | end |
---|
| 43 | |
---|
| 44 | if newmodel.epicost <= model{i}.epicost |
---|
| 45 | discard = 1; |
---|
| 46 | break |
---|
| 47 | end |
---|
| 48 | |
---|
| 49 | if newmodel.epicost >= model{i}.epicost |
---|
| 50 | replace(i,1) = 1; |
---|
| 51 | end |
---|
| 52 | |
---|
| 53 | elseif Pfinal >= model{i}.Pfinal |
---|
| 54 | |
---|
| 55 | if isempty(newmodel.epicost) |
---|
| 56 | B = reshape([details.Bi{:}]',size(details.Bi{1},2),[])'; |
---|
| 57 | c = reshape([details.Ci{:}]',[],1); |
---|
| 58 | [H,K] = double(Pfinal); |
---|
| 59 | newmodel.epicost = generate_epicost(H,K,B,c); |
---|
| 60 | end |
---|
| 61 | if isempty(model{i}.epicost) |
---|
| 62 | B = reshape([model{i}.Bi{:}]',size(model{i}.Bi{1},2),[])'; |
---|
| 63 | c = reshape([model{i}.Ci{:}]',[],1); |
---|
| 64 | [H,K] = double(model{i}.Pfinal); |
---|
| 65 | model{i}.epicost = generate_epicost(H,K,B,c); |
---|
| 66 | end |
---|
| 67 | if newmodel.epicost >= model{i}.epicost |
---|
| 68 | replace(i,1) = 1; |
---|
| 69 | end |
---|
| 70 | |
---|
| 71 | elseif Pfinal <= model{i}.Pfinal |
---|
| 72 | |
---|
| 73 | if isempty(newmodel.epicost) |
---|
| 74 | B = reshape([details.Bi{:}]',size(details.Bi{1},2),[])'; |
---|
| 75 | c = reshape([details.Ci{:}]',[],1); |
---|
| 76 | [H,K] = double(Pfinal); |
---|
| 77 | newmodel.epicost = generate_epicost(H,K,B,c); |
---|
| 78 | end |
---|
| 79 | if isempty(model{i}.epicost) |
---|
| 80 | B = reshape([model{i}.Bi{:}]',size(model{i}.Bi{1},2),[])'; |
---|
| 81 | c = reshape([model{i}.Ci{:}]',[],1); |
---|
| 82 | [H,K] = double(model{i}.Pfinal); |
---|
| 83 | model{i}.epicost = generate_epicost(H,K,B,c); |
---|
| 84 | end |
---|
| 85 | |
---|
| 86 | if newmodel.epicost <= model{i}.epicost |
---|
| 87 | discard = 1; |
---|
| 88 | break |
---|
| 89 | end |
---|
| 90 | end |
---|
| 91 | end |
---|
| 92 | if ~discard |
---|
| 93 | model = {model{find(~replace)},newmodel}; |
---|
| 94 | end |
---|
| 95 | end |
---|
| 96 | end |
---|
| 97 | end |
---|
| 98 | |
---|
| 99 | function model = savemptmodelqp(model,Pfinal,Pn,Fi,Gi,details); |
---|
| 100 | newmodel = fakemptmodel(Pfinal, Pn, Fi, Gi, details.Ai, details.Bi, details.Ci); |
---|
| 101 | newmodel.epicost = []; |
---|
| 102 | newmodel.convex = 1; |
---|
| 103 | |
---|
| 104 | replace = zeros(length(model),1); |
---|
| 105 | discard = 0; |
---|
| 106 | |
---|
| 107 | for i = 1:length(model) |
---|
| 108 | |
---|
| 109 | % Trivial pruning, why not... |
---|
| 110 | if quadraticLarger(details,model{i}) |
---|
| 111 | if Pfinal == model{i}.Pfinal |
---|
| 112 | discard = 1; |
---|
| 113 | break; |
---|
| 114 | elseif Pfinal <= model{i}.Pfinal |
---|
| 115 | discard = 1; |
---|
| 116 | break; |
---|
| 117 | end |
---|
| 118 | end |
---|
| 119 | |
---|
| 120 | if Pfinal==model{i}.Pfinal |
---|
| 121 | |
---|
| 122 | elseif Pfinal >= model{i}.Pfinal |
---|
| 123 | |
---|
| 124 | doreplace = 1; |
---|
| 125 | for k = 1:length(details.Ai) |
---|
| 126 | Qnew = [details.Ai{k} 0.5*details.Bi{k}' ;0.5*details.Bi{k} details.Ci{k}]; |
---|
| 127 | for j = 1:length(model{i}.Pfinal) |
---|
| 128 | Qold = [model{i}.Ai{j} 0.5*model{i}.Bi{j}';0.5*model{i}.Bi{j} model{i}.Ci{j}]; |
---|
| 129 | if ~all(real(eig(full(Qold-Qnew)))>=-1e-12) |
---|
| 130 | doreplace = 0; |
---|
| 131 | end |
---|
| 132 | end |
---|
| 133 | end |
---|
| 134 | if doreplace |
---|
| 135 | replace(i,1) = 1; |
---|
| 136 | end |
---|
| 137 | |
---|
| 138 | elseif Pfinal <= model{i}.Pfinal |
---|
| 139 | |
---|
| 140 | % if relaxationLarger(details,model{i},Pn,model{i}.Pn) |
---|
| 141 | % discard = 1; |
---|
| 142 | % break |
---|
| 143 | % end |
---|
| 144 | % |
---|
| 145 | if length(Pn)==1 |
---|
| 146 | discard = 1; |
---|
| 147 | Qnew = [details.Ai{1} 0.5*details.Bi{1}' ;0.5*details.Bi{1} details.Ci{1}]; |
---|
| 148 | for j = 1:length(model{i}.Pfinal) |
---|
| 149 | Qold = [model{i}.Ai{j} 0.5*model{i}.Bi{j}';0.5*model{i}.Bi{j} model{i}.Ci{j}]; |
---|
| 150 | if ~all(real(eig(full(Qnew-Qold)))>=-1e-12) |
---|
| 151 | discard = 0; |
---|
| 152 | end |
---|
| 153 | end |
---|
| 154 | if discard |
---|
| 155 | break |
---|
| 156 | end |
---|
| 157 | end |
---|
| 158 | % |
---|
| 159 | % if length(Pn)==1 |
---|
| 160 | % discard = 1; |
---|
| 161 | % Qnew = [details.Ai{1} 0.5*details.Bi{1}' ;0.5*details.Bi{1} details.Ci{1}]; |
---|
| 162 | % for j = 1:length(model{i}.Pn) |
---|
| 163 | % Qold = [model{i}.Ai{j} 0.5*model{i}.Bi{j}';0.5*model{i}.Bi{j} model{i}.Ci{j}]; |
---|
| 164 | % Pis = intersect(model{i}.Pn(j),newmodel.Pfinal); |
---|
| 165 | % if isfulldim(Pis) |
---|
| 166 | % [H,K] = double(Pis); |
---|
| 167 | % if ~quadraticLarger2(Qnew,Qold,H,K) |
---|
| 168 | % discard = 0; |
---|
| 169 | % end |
---|
| 170 | % end |
---|
| 171 | % end |
---|
| 172 | % if discard |
---|
| 173 | % break |
---|
| 174 | % end |
---|
| 175 | % end |
---|
| 176 | % |
---|
| 177 | end |
---|
| 178 | end |
---|
| 179 | |
---|
| 180 | if ~discard |
---|
| 181 | model = {model{find(~replace)},newmodel}; |
---|
| 182 | end |
---|
| 183 | |
---|
| 184 | function XbiggerY = relaxationLarger(X,Y,P1,P2) |
---|
| 185 | XbiggerY = 1; |
---|
| 186 | x = sdpvar(length(X.Ai{1}),1); |
---|
| 187 | for k = 1:length(X.Ai) |
---|
| 188 | Xq = [X.Ai{k} 0.5*X.Bi{k}' ;0.5*X.Bi{k} X.Ci{k}]; |
---|
| 189 | p1 = [x;1]'*Xq*[x;1]; |
---|
| 190 | for j = 1:length(Y.Ai) |
---|
| 191 | Yq = [Y.Ai{j} 0.5*Y.Bi{j}' ;0.5*Y.Bi{j} Y.Ci{j}]; |
---|
| 192 | if all(real(eig(full(Xq-Yq)))>0) |
---|
| 193 | else |
---|
| 194 | isc = intersect(P1(k),P2(j)); |
---|
| 195 | if isfulldim(isc) |
---|
| 196 | p2 = [x;1]'*Yq*[x;1]; |
---|
| 197 | [H,K] = double(isc); |
---|
| 198 | [xx,cc,vv]= solvemoment(set(H*x < K),p1-p2) |
---|
| 199 | end |
---|
| 200 | end |
---|
| 201 | % if ~all(real(eig(Xq-Yq))>=-1e-12) |
---|
| 202 | % XbiggerY = 0; |
---|
| 203 | % return |
---|
| 204 | % end |
---|
| 205 | end |
---|
| 206 | end |
---|
| 207 | |
---|
| 208 | |
---|
| 209 | function XbiggerY = quadraticLarger2(Q1,Q2,H,K) |
---|
| 210 | |
---|
| 211 | x = sdpvar(size(H,2),1); |
---|
| 212 | obj = [x;1]'*Q1*[x;1]-[x;1]'*Q2*[x;1]; |
---|
| 213 | |
---|
| 214 | %sol = solvemoment(set(H*x <= K),obj,[],2); |
---|
| 215 | sol = solvesdp(set(H*x <= K),obj,sdpsettings('solver','kktqp')); |
---|
| 216 | %if relaxdouble(obj) > -1e-5 |
---|
| 217 | if double(obj) > -1e-5 |
---|
| 218 | |
---|
| 219 | XbiggerY = 1; |
---|
| 220 | else |
---|
| 221 | XbiggerY = 0; |
---|
| 222 | end |
---|
| 223 | |
---|
| 224 | function XbiggerY = quadraticLarger(X,Y) |
---|
| 225 | XbiggerY = 1; |
---|
| 226 | for k = 1:length(X.Ai) |
---|
| 227 | Xq = [X.Ai{k} 0.5*X.Bi{k}' ;0.5*X.Bi{k} X.Ci{k}]; |
---|
| 228 | for j = 1:length(Y.Ai) |
---|
| 229 | Yq = [Y.Ai{j} 0.5*Y.Bi{j}' ;0.5*Y.Bi{j} Y.Ci{j}]; |
---|
| 230 | if ~all(real(eig(Xq-Yq))>=-1e-12) |
---|
| 231 | XbiggerY = 0; |
---|
| 232 | return |
---|
| 233 | end |
---|
| 234 | end |
---|
| 235 | end |
---|
| 236 | |
---|
| 237 | function epicost = generate_epicost(H,K,B,c) |
---|
| 238 | epicost = polytope([B -ones(size(B,1),1);H zeros(size(H,1),1)],[-c;K]); |
---|
| 239 | |
---|
| 240 | |
---|
| 241 | |
---|
| 242 | function C = fakemptmodel(Pfinal, Pn, Fi, Gi, Ai, Bi, Ci) |
---|
| 243 | |
---|
| 244 | dummystruct.dummyfield = []; |
---|
| 245 | |
---|
| 246 | C.sysStruct = dummystruct; |
---|
| 247 | C.probStruct = dummystruct; |
---|
| 248 | C.details.origSysStruct = dummystruct; |
---|
| 249 | C.details.origProbStruct = dummystruct; |
---|
| 250 | nr = length(Pn); |
---|
| 251 | |
---|
| 252 | C.Pfinal = Pfinal; |
---|
| 253 | C.Pn = Pn; |
---|
| 254 | C.Fi = Fi; |
---|
| 255 | C.Gi = Gi; |
---|
| 256 | if isempty(Ai), |
---|
| 257 | Ai = cell(1, nr); |
---|
| 258 | end |
---|
| 259 | C.Ai = Ai; |
---|
| 260 | C.Bi = Bi; |
---|
| 261 | C.Ci = Ci; |
---|
| 262 | C.dynamics = repmat(1, 1, nr); |
---|
| 263 | C.overlaps = 0; |
---|