[37] | 1 | function varargout = pwq_yalmip(varargin) |
---|
| 2 | %PWQ_YALMIP Defines a piecewise quadratic function using data from MPT |
---|
| 3 | % |
---|
| 4 | %Only intended for internal use in YALMIP |
---|
| 5 | % |
---|
| 6 | % Currently only a container for PWQ functions. Can not be |
---|
| 7 | % used in actual optmization problem. |
---|
| 8 | |
---|
| 9 | % Author Johan Löfberg |
---|
| 10 | % $Id: pwq_yalmip.m,v 1.2 2006/05/14 12:39:57 joloef Exp $ |
---|
| 11 | |
---|
| 12 | switch class(varargin{1}) |
---|
| 13 | |
---|
| 14 | case {'struct','cell'} % Should only be called internally |
---|
| 15 | |
---|
| 16 | if isa(varargin{2},'double') |
---|
| 17 | % Called from YALMIP to get double |
---|
| 18 | pwastruct = varargin{1}; |
---|
| 19 | x = varargin{2}; |
---|
| 20 | index = varargin{5}; |
---|
| 21 | |
---|
| 22 | val = inf; |
---|
| 23 | for i = 1:length(pwastruct) |
---|
| 24 | [ii,jj] = isinside(pwastruct{i}.Pn,x); |
---|
| 25 | if ii |
---|
| 26 | for k = 1:length(jj) |
---|
| 27 | Q = pwastruct{i}.Ai{jj(k)}; |
---|
| 28 | val = min(val,x'*Q*x + pwastruct{i}.Bi{jj(k)}(index,:)*x+pwastruct{i}.Ci{jj(k)}(index)); |
---|
| 29 | end |
---|
| 30 | end |
---|
| 31 | end |
---|
| 32 | if isinf(val) |
---|
| 33 | val = nan; |
---|
| 34 | end |
---|
| 35 | varargout{1} = val; |
---|
| 36 | return |
---|
| 37 | end |
---|
| 38 | |
---|
| 39 | if nargin<3 |
---|
| 40 | pwaclass = 'general' |
---|
| 41 | end |
---|
| 42 | |
---|
| 43 | if isa(varargin{1},'struct') |
---|
| 44 | varargin{1} = {varargin{1}}; |
---|
| 45 | end |
---|
| 46 | |
---|
| 47 | % Put in standard format |
---|
| 48 | if ~isfield(varargin{1}{1},'Bi') |
---|
| 49 | if ~isfield(varargin{1}{1},'Fi') |
---|
| 50 | error('Wrong format on input to PWQ (requires Bi or Fi)'); |
---|
| 51 | else |
---|
| 52 | for i = 1:length(varargin{1}) |
---|
| 53 | varargin{1}{1}.Ai = cell(1, varargin{1}{i}.Fi); |
---|
| 54 | varargin{1}{i}.Bi = varargin{1}{i}.Fi |
---|
| 55 | varargin{1}{i}.Ci = varargin{1}{i}.Gi |
---|
| 56 | end |
---|
| 57 | end |
---|
| 58 | end |
---|
| 59 | |
---|
| 60 | if ~isfield(varargin{1}{1},'Pfinal') |
---|
| 61 | error('Wrong format on input to PWQ (requires field Pn)'); |
---|
| 62 | end |
---|
| 63 | |
---|
| 64 | % This will be a container for binary variables in the furture |
---|
| 65 | varargin{end+1} = binvar(length(varargin{1}{1}.Pn),1); |
---|
| 66 | |
---|
| 67 | % Create one variable for each row |
---|
| 68 | % Inefficient but the only way currently in YALMIP |
---|
| 69 | varargout{1} = []; |
---|
| 70 | varargin{end+1} = 1; |
---|
| 71 | |
---|
| 72 | for i = 1:length(varargin{1}{1}.Ci{1}) |
---|
| 73 | varargin{end} = i; |
---|
| 74 | varargout{1} = [varargout{1};yalmip('addextendedvariable',mfilename,varargin{:})]; |
---|
| 75 | end |
---|
| 76 | |
---|
| 77 | case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | switch varargin{1} |
---|
| 81 | |
---|
| 82 | case 'graph' |
---|
| 83 | |
---|
| 84 | varargout{1} = []; |
---|
| 85 | varargout{2} = struct('convexity','milp','monotoncity','none','definiteness','none'); |
---|
| 86 | varargout{3} = []; |
---|
| 87 | |
---|
| 88 | case 'milp' |
---|
| 89 | |
---|
| 90 | % FIX : Should create case for overlapping convex PWAs, |
---|
| 91 | % used in a nonconvex fashion... |
---|
| 92 | |
---|
| 93 | % Can only generate the first class of PWA functions |
---|
| 94 | t = varargin{2}; % The YALMIP variables modelling this pwa |
---|
| 95 | pwq_struct = varargin{3};% MPT structure |
---|
| 96 | x = varargin{4}; % Argument |
---|
| 97 | flag = varargin{5}; % Type of PWA function |
---|
| 98 | d = varargin{6}; % Binary for nonconvex cases |
---|
| 99 | index = varargin{7}; % Which row in Bix+Ci |
---|
| 100 | |
---|
| 101 | switch flag |
---|
| 102 | |
---|
| 103 | case {'general','convex'} |
---|
| 104 | |
---|
| 105 | if length(d)==1 |
---|
| 106 | % Don't introduce any binary variables when |
---|
| 107 | % there is only one quadratic cost |
---|
| 108 | F = set([]); |
---|
| 109 | cost = x'*pwq_struct{1}.Ai{1}*x+pwq_struct{1}.Bi{1}*x+pwq_struct{1}.Ci{1}; |
---|
| 110 | else |
---|
| 111 | n = length(x); |
---|
| 112 | m = length(d); |
---|
| 113 | z = sdpvar(n,m); |
---|
| 114 | F = set(sum(d) == 1); |
---|
| 115 | cost = 0; |
---|
| 116 | [Mx,mx] = derivebounds(x); |
---|
| 117 | for i = 1:m |
---|
| 118 | bounds(z(:,i),mx,Mx); |
---|
| 119 | F = F + set(-(Mx-mx)*(1-d(i)) <= z(:,i)-x <= (Mx-mx)*(1-d(i))); |
---|
| 120 | F = F + set(5*mx*d(i) <= z(:,i) <= 5*Mx*d(i)); |
---|
| 121 | [H,K] = double(pwq_struct{1}.Pn(i)); |
---|
| 122 | [M,m] = derivebounds(H*x-K); |
---|
| 123 | F = F + set(H*x-K <= M*(1-d(i))); |
---|
| 124 | cost = cost + z(:,i)'*pwq_struct{1}.Ai{i}*z(:,i)+pwq_struct{1}.Bi{i}*z(:,i)+d(i)*pwq_struct{1}.Ci{i}; |
---|
| 125 | end |
---|
| 126 | end |
---|
| 127 | |
---|
| 128 | varargout{1} = F; |
---|
| 129 | varargout{2} = struct('convexity','milp','monotoncity','milp','definiteness','milp'); |
---|
| 130 | varargout{2}.replacer = cost; |
---|
| 131 | varargout{3} = x; |
---|
| 132 | |
---|
| 133 | case {'convexoverlapping'} |
---|
| 134 | |
---|
| 135 | n = length(x); |
---|
| 136 | cost = 0; |
---|
| 137 | r = binvar(length(pwq_struct),1); |
---|
| 138 | d = {}; |
---|
| 139 | % % Some overlapping function is active |
---|
| 140 | % % (it will automatically be the smallest one) |
---|
| 141 | F = set(sum(r) == 1); |
---|
| 142 | for j = 1:length(pwq_struct) |
---|
| 143 | |
---|
| 144 | if length(pwq_struct{j}.Pn) == 1 |
---|
| 145 | d{j} = r(j); |
---|
| 146 | else |
---|
| 147 | d{j} = binvar(length(pwq_struct{j}.Pn),1); |
---|
| 148 | F = F + set(sum(d{j}) == r(j)); |
---|
| 149 | end |
---|
| 150 | m = length(d{j}); |
---|
| 151 | z = sdpvar(n,m,'full'); |
---|
| 152 | |
---|
| 153 | [Mx,mx] = derivebounds(x); |
---|
| 154 | for i = 1:m |
---|
| 155 | bounds(z(:,i),mx,Mx); |
---|
| 156 | F = F + set(-(Mx-mx)*(1-d{j}(i)) <= z(:,i)-x <= (Mx-mx)*(1-d{j}(i))); |
---|
| 157 | F = F + set(mx*d{j}(i) <= z(:,i) <= Mx*d{j}(i)); |
---|
| 158 | [H,K] = double(pwq_struct{j}.Pn(i)); |
---|
| 159 | [M,m] = derivebounds(H*x-K); |
---|
| 160 | F = F + set(H*x-K <= M*(1-d{j}(i))); |
---|
| 161 | cost = cost + z(:,i)'*pwq_struct{j}.Ai{i}*z(:,i)+pwq_struct{j}.Bi{i}*z(:,i)+d{j}(i)*pwq_struct{j}.Ci{i}; |
---|
| 162 | end |
---|
| 163 | end |
---|
| 164 | |
---|
| 165 | varargout{1} = F; |
---|
| 166 | varargout{2} = struct('convexity','milp','monotoncity','milp','definiteness','milp'); |
---|
| 167 | varargout{2}.replacer = cost; |
---|
| 168 | varargout{3} = x; |
---|
| 169 | |
---|
| 170 | |
---|
| 171 | otherwise |
---|
| 172 | |
---|
| 173 | varargout{1} = []; |
---|
| 174 | varargout{2} = struct('convexity','convex','monotoncity','none','definiteness','none'); |
---|
| 175 | varargout{3} = []; |
---|
| 176 | return |
---|
| 177 | end |
---|
| 178 | |
---|
| 179 | otherwise |
---|
| 180 | error('PWA_YALMIP called with CHAR argument?'); |
---|
| 181 | end |
---|
| 182 | |
---|
| 183 | %disp('*************************************************************************') |
---|
| 184 | %disp('PWQ functions can not be used yet.') |
---|
| 185 | %disp('Are you perhaps trying to use a quadratic value function? Not possible...'); |
---|
| 186 | %disp('*************************************************************************') |
---|
| 187 | %disp(' '); |
---|
| 188 | %error('PWQ_YALMIP still not implemented'); |
---|
| 189 | |
---|
| 190 | otherwise |
---|
| 191 | error('Strange type on first argument in PWQ_YALMIP'); |
---|
| 192 | end |
---|