1 | function varargout = pwa_yalmip(varargin) |
---|
2 | %PWA_YALMIP Defines a piecewise function using data from MPT |
---|
3 | % |
---|
4 | %Only intended for internal use in YALMIP |
---|
5 | % |
---|
6 | % Three classes of PWA functions can be generated |
---|
7 | % |
---|
8 | % 1) Convexity aware epigraph version with a |
---|
9 | % convex domain and objective modelled as c'x <t. |
---|
10 | % Note, if used in non-concvex setting, a milp |
---|
11 | % model will be generated as a back-up. |
---|
12 | % |
---|
13 | % This is the function generated to describe value |
---|
14 | % function of a single mpLP problem. |
---|
15 | % |
---|
16 | % 2) Overlapping partitions with convex functions |
---|
17 | % in each partition. Requires one binary per region. |
---|
18 | % |
---|
19 | % This is the function used to describe the value |
---|
20 | % function generated when not removing overlaps in |
---|
21 | % binary mpLP. |
---|
22 | % |
---|
23 | % 3) Non-overlapping general pwa function. Requires |
---|
24 | % one binary per region. |
---|
25 | % |
---|
26 | % This is the function generated for value functions |
---|
27 | % when remove overlaps is done. It is also used to |
---|
28 | % describe the pwa optimizer. |
---|
29 | % |
---|
30 | % The input is a cell with structs. Each struct has |
---|
31 | % the fields Pn, Pfinal, Bi and Ci (or Fi and Gi) |
---|
32 | |
---|
33 | % Author Johan Löfberg |
---|
34 | % $Id: pwa_yalmip.m,v 1.3 2006/06/07 14:36:32 joloef Exp $ |
---|
35 | |
---|
36 | switch class(varargin{1}) |
---|
37 | |
---|
38 | case {'struct','cell'} % Should only be called internally |
---|
39 | |
---|
40 | if isa(varargin{2},'double') |
---|
41 | % Called from YALMIP to get double |
---|
42 | pwastruct = varargin{1}; |
---|
43 | x = varargin{2}; |
---|
44 | index = varargin{5}; |
---|
45 | |
---|
46 | val = inf; |
---|
47 | for i = 1:length(pwastruct) |
---|
48 | [ii,jj] = isinside(pwastruct{i}.Pn,x); |
---|
49 | if ii |
---|
50 | for k = 1:length(jj) |
---|
51 | val = min(val,pwastruct{i}.Bi{jj(k)}(index,:)*x+pwastruct{i}.Ci{jj(k)}(index)); |
---|
52 | end |
---|
53 | end |
---|
54 | end |
---|
55 | if isinf(val) |
---|
56 | val = nan; |
---|
57 | end |
---|
58 | varargout{1} = val; |
---|
59 | return |
---|
60 | end |
---|
61 | |
---|
62 | if nargin<3 |
---|
63 | pwaclass = 'general' |
---|
64 | end |
---|
65 | |
---|
66 | if isa(varargin{1},'struct') |
---|
67 | varargin{1} = {varargin{1}}; |
---|
68 | end |
---|
69 | |
---|
70 | % Put in standard format |
---|
71 | if ~isfield(varargin{1}{1},'Bi') |
---|
72 | if ~isfield(varargin{1}{1},'Fi') |
---|
73 | error('Wrong format on input to PWA (requires Bi or Fi)'); |
---|
74 | else |
---|
75 | for i = 1:length(varargin{1}) |
---|
76 | varargin{1}{i}.Bi = varargin{1}{i}.Fi |
---|
77 | varargin{1}{i}.Ci = varargin{1}{i}.Gi |
---|
78 | end |
---|
79 | end |
---|
80 | end |
---|
81 | |
---|
82 | if ~isfield(varargin{1}{1},'Pfinal') |
---|
83 | error('Wrong format on input to PWA (requires field Pn)'); |
---|
84 | end |
---|
85 | |
---|
86 | % Create binary variables already here to avoid generating |
---|
87 | % binary variables for the same variable several times |
---|
88 | switch varargin{3} |
---|
89 | case 'convex' |
---|
90 | % No binary needed |
---|
91 | varargin{end+1} = []; |
---|
92 | case 'convexoverlapping' |
---|
93 | % One binary per overlap |
---|
94 | varargin{end+1} = binvar(length(varargin{1}),1); |
---|
95 | case 'general' |
---|
96 | % one binary per region |
---|
97 | varargin{end+1} = binvar(length(varargin{1}{1}.Pn),1); |
---|
98 | otherwise |
---|
99 | end |
---|
100 | |
---|
101 | % Create one variable for each row |
---|
102 | % Inefficient but the only way currently in YALMIP |
---|
103 | varargout{1} = []; |
---|
104 | varargin{end+1} = 1; |
---|
105 | |
---|
106 | for i = 1:length(varargin{1}{1}.Ci{1}) |
---|
107 | varargin{end} = i; |
---|
108 | varargout{1} = [varargout{1};yalmip('addextendedvariable',mfilename,varargin{:})]; |
---|
109 | end |
---|
110 | |
---|
111 | case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph |
---|
112 | |
---|
113 | switch varargin{1} |
---|
114 | |
---|
115 | case 'graph' |
---|
116 | |
---|
117 | % Can only generate the first class of PWA functions |
---|
118 | t = varargin{2}; % The YALMIP variables modelling this pwa |
---|
119 | pwa_struct = varargin{3};% MPT structure |
---|
120 | x = varargin{4}; % Argument |
---|
121 | flag = varargin{5}; % Type of PWA function |
---|
122 | d = varargin{6}; % Binary for nonconvex cases |
---|
123 | index = varargin{7}; % Which row in Bix+Ci |
---|
124 | |
---|
125 | switch flag |
---|
126 | |
---|
127 | case 'convex' |
---|
128 | |
---|
129 | [H,K] = double(pwa_struct{1}.Pfinal); |
---|
130 | costs = reshape([pwa_struct{1}.Bi{:}]',length(x),[])'*x+reshape([pwa_struct{1}.Ci{:}]',[],1); |
---|
131 | F = set(H* x <= K) + set(costs< t); |
---|
132 | |
---|
133 | case 'convexoverlapping' |
---|
134 | |
---|
135 | % Local costs |
---|
136 | s = sdpvar(length(pwa_struct),1); |
---|
137 | |
---|
138 | % Some region, one defined as minimizer |
---|
139 | F = set(sum(d)==1); |
---|
140 | maxcost = -inf; |
---|
141 | mincost = inf; |
---|
142 | for i = 1:length(pwa_struct) |
---|
143 | % Reduce complexity of max function by |
---|
144 | % removing equal costs |
---|
145 | B = reshape([pwa_struct{i}.Bi{:}]',length(x),[])'; |
---|
146 | C = reshape([pwa_struct{i}.Ci{:}]',[],1); |
---|
147 | [ii,jj,kk] = unique(round(1e8*[B C])/1e8,'rows'); |
---|
148 | cost = reshape([pwa_struct{i}.Bi{jj}]',length(x),[])'*x+reshape([pwa_struct{i}.Ci{jj}]',[],1); |
---|
149 | [M,m] = derivebounds(cost); |
---|
150 | maxcost = max(maxcost,max(M)); |
---|
151 | mincost = min(mincost,min(m)); |
---|
152 | [H,K] = double(pwa_struct{i}.Pfinal); |
---|
153 | [M,m] = derivebounds(H*x - K); |
---|
154 | F = F + set(H*x-K <= (1+M)*(1-d(i))); |
---|
155 | end |
---|
156 | for i = 1:length(pwa_struct) |
---|
157 | |
---|
158 | B = reshape([pwa_struct{i}.Bi{:}]',length(x),[])'; |
---|
159 | C = reshape([pwa_struct{i}.Ci{:}]',[],1); |
---|
160 | [ii,jj,kk] = unique(round(1e8*[B C])/1e8,'rows'); |
---|
161 | |
---|
162 | cost = reshape([pwa_struct{i}.Bi{jj}]',length(x),[])'*x+reshape([pwa_struct{i}.Ci{jj}]',[],1); |
---|
163 | |
---|
164 | [M,m] = derivebounds(cost); |
---|
165 | F = F + set(cost < t + 2*(1+maxcost)*(1-d(i))); |
---|
166 | end |
---|
167 | [t_bounds] = yalmip('getbounds',getvariables(t)); |
---|
168 | bounds(t,max([t_bounds(1) mincost]),min([t_bounds(2) 3*maxcost])); |
---|
169 | |
---|
170 | |
---|
171 | case 'general' |
---|
172 | |
---|
173 | % In one region |
---|
174 | F = set(sum(d) == 1); |
---|
175 | |
---|
176 | % Extract the wanted row |
---|
177 | for i = 1:length(pwa_struct{1}.Pn) |
---|
178 | pwa_struct{1}.Bi{i} = pwa_struct{1}.Bi{i}(index,:); |
---|
179 | pwa_struct{1}.Ci{i} = pwa_struct{1}.Ci{i}(index); |
---|
180 | end |
---|
181 | |
---|
182 | % value in different regions |
---|
183 | costs = reshape([pwa_struct{1}.Bi{:}]',length(x),[])'*x+reshape([pwa_struct{1}.Ci{:}]',[],1); |
---|
184 | [M,m] = derivebounds(costs); |
---|
185 | |
---|
186 | % We know that t is always between min(m) and max(M). |
---|
187 | % However, user might know more |
---|
188 | [t_bounds] = yalmip('getbounds',getvariables(t)); |
---|
189 | bounds(t,max([t_bounds(1) min(m)]),min([t_bounds(2) max(M)])); |
---|
190 | |
---|
191 | % t equals some of the costs |
---|
192 | % Big-M : |costs-t| < max(costs)-min(t) < M-m |
---|
193 | F = F + set( -2*(M-m+1).*(1-d) <= costs-t <= 2*(1+M-m).*(1-d)); |
---|
194 | for i = 1:length(pwa_struct{1}.Pn) |
---|
195 | [H,K] = double(pwa_struct{1}.Pn(i)); |
---|
196 | [M,m] = derivebounds(H*x - K); |
---|
197 | F = F + set(H*x - K <= M*(1-d(i))); |
---|
198 | end |
---|
199 | |
---|
200 | otherwise |
---|
201 | |
---|
202 | varargout{1} = []; |
---|
203 | varargout{2} = struct('convexity','milp','monotoncity','none','definiteness','none'); |
---|
204 | varargout{3} = []; |
---|
205 | return |
---|
206 | end |
---|
207 | |
---|
208 | varargout{1} = F; |
---|
209 | varargout{2} = struct('convexity','convex','monotoncity','none','definiteness','none'); |
---|
210 | varargout{3} = x; |
---|
211 | |
---|
212 | case 'milp' |
---|
213 | |
---|
214 | % FIX : Should create case for overlapping convex PWAs, |
---|
215 | % used in a nonconvex fashion... |
---|
216 | |
---|
217 | % Can only generate the first class of PWA functions |
---|
218 | t = varargin{2}; % The YALMIP variables modelling this pwa |
---|
219 | pwa_struct = varargin{3};% MPT structure |
---|
220 | x = varargin{4}; % Argument |
---|
221 | flag = varargin{5}; % Type of PWA function |
---|
222 | d = varargin{6}; % Binary for nonconvex cases |
---|
223 | index = varargin{7}; % Which row in Bix+Ci |
---|
224 | |
---|
225 | switch flag |
---|
226 | |
---|
227 | case {'general','convex'} |
---|
228 | |
---|
229 | F = set(sum(d) == 1); |
---|
230 | |
---|
231 | % Extract the wanted row |
---|
232 | for i = 1:length(pwa_struct{1}.Pn) |
---|
233 | pwa_struct{1}.Bi{i} = pwa_struct{1}.Bi{i}(index,:); |
---|
234 | pwa_struct{1}.Ci{i} = pwa_struct{1}.Ci{i}(index); |
---|
235 | end |
---|
236 | |
---|
237 | % Cost in different regions |
---|
238 | costs = reshape([pwa_struct{1}.Bi{:}]',length(x),[])'*x+reshape([pwa_struct{1}.Ci{:}]',[],1); |
---|
239 | [M,m] = derivebounds(costs); |
---|
240 | |
---|
241 | % We know that t is always between min(m) and max(M). |
---|
242 | % However, user might know more |
---|
243 | [t_bounds] = yalmip('getbounds',getvariables(t)); |
---|
244 | bounds(t,max([t_bounds(1) min(m)]),min([t_bounds(2) max(M)])); |
---|
245 | |
---|
246 | % Might be called with a convex function, |
---|
247 | % but wants the complete MILP description |
---|
248 | % Example : Someone tries to maximize value |
---|
249 | % function |
---|
250 | if length(d)~=length(costs) |
---|
251 | d = binvar(length(costs),1); |
---|
252 | F = set(sum(d) == 1); |
---|
253 | end |
---|
254 | % t equals some of the costs |
---|
255 | % Big-M : |costs-t| < max(costs)-min(t) < M-m |
---|
256 | F = F + set( -2*(M-m+1).*(1-d) <= costs-t <= 2*(1+M-m).*(1-d)); |
---|
257 | for i = 1:length(pwa_struct{1}.Pn) |
---|
258 | [H,K] = double(pwa_struct{1}.Pn(i)); |
---|
259 | [M,m] = derivebounds(H*x - K); |
---|
260 | F = F + set(H*x - K <= (M+1)*(1-d(i))); |
---|
261 | end |
---|
262 | |
---|
263 | otherwise |
---|
264 | |
---|
265 | varargout{1} = []; |
---|
266 | varargout{2} = struct('convexity','convex','monotoncity','none','definiteness','none'); |
---|
267 | varargout{3} = []; |
---|
268 | return |
---|
269 | end |
---|
270 | |
---|
271 | varargout{1} = F; |
---|
272 | varargout{2} = struct('convexity','milp','monotoncity','milp','definiteness','milp'); |
---|
273 | varargout{3} = x; |
---|
274 | |
---|
275 | |
---|
276 | otherwise |
---|
277 | error('PWA_YALMIP called with CHAR argument?'); |
---|
278 | end |
---|
279 | otherwise |
---|
280 | error('Strange type on first argument in SDPVAR/NORM'); |
---|
281 | end |
---|