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; |
---|