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