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