1 | function [p,Bi,Ci,Pn,Pfinal] = PWA(h,Xdomain) |
---|
2 | % PWA Tries to create a PWA description |
---|
3 | % |
---|
4 | % [p,Bi,Ci,Pn,Pfinal] = PWA(h,X) |
---|
5 | % |
---|
6 | % Input |
---|
7 | % h : scalar SDPVAR object |
---|
8 | % X : SET object |
---|
9 | % |
---|
10 | % Output |
---|
11 | % |
---|
12 | % p : scalar SDPVAR object representing the PWA function |
---|
13 | % Bi,Ci,Pn,Pfinal : Data in MPT format |
---|
14 | % |
---|
15 | % The command tries to expand the nonlinear operators |
---|
16 | % (min,max,abs,...) used in the variable h, in order |
---|
17 | % to generate an epi-graph model. Given this epigraph model, |
---|
18 | % it is projected to the variables of interest and the |
---|
19 | % defining facets of the PWA function is extracted. The second |
---|
20 | % argument can be used to limit the domain of the PWA function. |
---|
21 | % If no second argument is supplied, the PWA function is created |
---|
22 | % over the domain -100 to 100. |
---|
23 | % |
---|
24 | % A new sdpvar object p is created, representing the same |
---|
25 | % function as h, but in a slightly different internal format. |
---|
26 | % Additionally, the PWA description in MPT format is created. |
---|
27 | % |
---|
28 | % The function is mainly inteded to be used for easy |
---|
29 | % plotting of convex PWA functions |
---|
30 | |
---|
31 | % Author Johan Löfberg |
---|
32 | % $Id: pwa.m,v 1.8 2006/02/07 10:58:39 joloef Exp $ |
---|
33 | |
---|
34 | |
---|
35 | t = sdpvar(1,1); |
---|
36 | [F,failure,cause] = expandmodel(set(h<t),[],sdpsettings('allowmilp',0)); |
---|
37 | if failure |
---|
38 | error(['Could not expand model (' cause ')']); |
---|
39 | return |
---|
40 | end |
---|
41 | |
---|
42 | % Figure out what the actual original variables are |
---|
43 | % note, by construction, they |
---|
44 | all_initial = getvariables(h); |
---|
45 | all_extended = yalmip('extvariables'); |
---|
46 | all_variables = getvariables(F); |
---|
47 | gen_here = getvariables(t); |
---|
48 | non_ext_in = setdiff(all_initial,all_extended); |
---|
49 | lifted = all_variables(all_variables>gen_here); |
---|
50 | vars = union(setdiff(setdiff(setdiff(all_variables,all_extended),gen_here),lifted),non_ext_in); |
---|
51 | nx = length(vars); |
---|
52 | |
---|
53 | X = recover(vars); |
---|
54 | |
---|
55 | if nargin == 1 |
---|
56 | Xdomain = set(-100 < X < 100); |
---|
57 | else |
---|
58 | Xdomain = set(-10000 < X < 10000)+Xdomain; |
---|
59 | end |
---|
60 | |
---|
61 | [Ai,Bi,Ci,Pn] = generate_pwa(F,t,X,Xdomain,nx); |
---|
62 | |
---|
63 | Pfinal = union(Pn); |
---|
64 | sol.Pn = Pn; |
---|
65 | sol.Bi = Bi; |
---|
66 | sol.Ci = Ci; |
---|
67 | sol.Ai = Ai; |
---|
68 | sol.Pfinal = Pfinal; |
---|
69 | p = pwf(sol,X,'convex'); |
---|
70 | % |
---|
71 | % binarys = recover(all_variables(find(ismember(all_variables,yalmip('binvariables'))))) |
---|
72 | % if length(binarys) > 0 |
---|
73 | % |
---|
74 | % Binary_Equalities = []; |
---|
75 | % Binary_Inequalities = []; |
---|
76 | % Mixed_Equalities = []; |
---|
77 | % top = 1; |
---|
78 | % for i = 1:length(F) |
---|
79 | % Fi = sdpvar(F(i)); |
---|
80 | % if is(F(i),'equality') |
---|
81 | % if all(ismember(getvariables(Fi),yalmip('binvariables'))) |
---|
82 | % Binary_Equalities = [Binary_Equalities;(top:top-1+prod(size(Fi)))']; |
---|
83 | % Mixed_Equalities = [Mixed_Equalities;(top:top-1+prod(size(Fi)))']; |
---|
84 | % end |
---|
85 | % else |
---|
86 | % if all(ismember(getvariables(Fi),yalmip('binvariables'))) |
---|
87 | % Binary_Inequalities = [Binary_Inequalities;(top:top-1+prod(size(Fi)))']; |
---|
88 | % end |
---|
89 | % end |
---|
90 | % top = top+prod(size(Fi))'; |
---|
91 | % end |
---|
92 | % P = sdpvar(F); |
---|
93 | % P_ineq = extsubsref(P,setdiff(1:length(P),[Binary_Equalities; Binary_Inequalities])) |
---|
94 | % P_binary_eq = extsubsref(P,Binary_Equalities);HK1 = getbase(P_binary_eq); |
---|
95 | % P_binary_ineq = extsubsref(P,Binary_Inequalities);HK2 = getbase(P_binary_ineq); |
---|
96 | % nbin = length(binarys); |
---|
97 | % enums = dec2decbin(0:2^nbin-1,nbin)' |
---|
98 | % if isempty(HK2) |
---|
99 | % HK2 = HK1*0; |
---|
100 | % end |
---|
101 | % for i = 1:size(enums,2) |
---|
102 | % if all(HK1*[1;enums(:,i)]==0) |
---|
103 | % if all(HK2*[1;enums(:,i)]>=0) |
---|
104 | % Pi = replace(P_ineq,binarys,enums(:,i)) |
---|
105 | % end |
---|
106 | % end |
---|
107 | % end |
---|
108 | % |
---|
109 | % else |
---|
110 | % [Ai,Bi,Ci,Pn] = generate_pwa(F,t,X,Xdomain,nx); |
---|
111 | % end |
---|
112 | % |
---|
113 | |
---|
114 | |
---|
115 | |
---|
116 | function [Ai,Bi,Ci,Pn] = generate_pwa(F,t,X,Xdomain,nx) |
---|
117 | |
---|
118 | % Project, but remember that we already expanded the constraints |
---|
119 | P = polytope(projection(F+set(t<10000)+Xdomain,[X;t],[],1));Xdomain = polytope(Xdomain); |
---|
120 | [H,K] = double(P); |
---|
121 | facets = find(H(:,end)<0); |
---|
122 | |
---|
123 | region = find(~H(:,end) & any(H(:,1:nx),2) ); |
---|
124 | Hr = H(region,1:nx); |
---|
125 | Kr = K(region,:); |
---|
126 | |
---|
127 | H = H(facets,:); |
---|
128 | K = K(facets); |
---|
129 | K = K./H(:,end); |
---|
130 | H = H./repmat(H(:,end),1,size(H,2)); |
---|
131 | |
---|
132 | nx = length(X); |
---|
133 | Pn = []; |
---|
134 | cib = [H(:,1:nx) K]; |
---|
135 | Ai = {}; |
---|
136 | Bi = cell(0); |
---|
137 | Ci = cell(0); |
---|
138 | if length(Kr > 0) |
---|
139 | Xdomain = intersect(Xdomain,polytope(Hr,Kr)); |
---|
140 | end |
---|
141 | [Hr,Kr] = double(Xdomain); |
---|
142 | |
---|
143 | for i = 1:length(K) |
---|
144 | j = setdiff(1:length(K),i); |
---|
145 | HiKi = repmat(cib(i,:),length(K)-1,1)-cib(j,:); |
---|
146 | Pi = polytope([HiKi(:,1:nx);Hr],[HiKi(:,end);Kr]); |
---|
147 | |
---|
148 | if isfulldim(Pi) |
---|
149 | Pn = [Pn Pi]; |
---|
150 | Bi{end+1} = -cib(i,1:end-1); |
---|
151 | Ci{end+1} = cib(i,end); |
---|
152 | Ai{end+1} = []; |
---|
153 | end |
---|
154 | end |
---|