Mercurial > hg > octave-nkf
annotate scripts/optimization/qp.m @ 8507:cadc73247d65
style fixes
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 13 Jan 2009 14:08:36 -0500 |
parents | 5ee11a81688e |
children | eb63fbe60fab |
rev | line source |
---|---|
7017 | 1 ## Copyright (C) 2000, 2001, 2004, 2005, 2006, 2007 Gabriele Pannocchia. |
5289 | 2 ## |
3 ## This file is part of Octave. | |
4 ## | |
5 ## Octave is free software; you can redistribute it and/or modify it | |
6 ## under the terms of the GNU General Public License as published by | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
5289 | 9 ## |
10 ## Octave is distributed in the hope that it will be useful, but | |
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 ## General Public License for more details. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
5289 | 18 |
19 ## -*- texinfo -*- | |
6741 | 20 ## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb}, @var{A_in}, @var{A_ub}) |
5289 | 21 ## Solve the quadratic program |
6741 | 22 ## @iftex |
23 ## @tex | |
24 ## $$ | |
25 ## \min_x {1 \over 2} x^T H x + x^T q | |
26 ## $$ | |
27 ## @end tex | |
28 ## @end iftex | |
29 ## @ifnottex | |
5289 | 30 ## |
31 ## @example | |
32 ## min 0.5 x'*H*x + x'*q | |
33 ## x | |
34 ## @end example | |
35 ## | |
6741 | 36 ## @end ifnottex |
37 ## subject to | |
5289 | 38 ## @iftex |
39 ## @tex | |
6741 | 40 ## $$ |
41 ## Ax = b \qquad lb \leq x \leq ub \qquad A_{lb} \leq A_{in} \leq A_{ub} | |
42 ## $$ | |
5289 | 43 ## @end tex |
44 ## @end iftex | |
6741 | 45 ## @ifnottex |
5289 | 46 ## |
47 ## @example | |
48 ## A*x = b | |
49 ## lb <= x <= ub | |
6741 | 50 ## A_lb <= A_in*x <= A_ub |
5289 | 51 ## @end example |
6741 | 52 ## @end ifnottex |
5289 | 53 ## |
54 ## @noindent | |
55 ## using a null-space active-set method. | |
56 ## | |
57 ## Any bound (@var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb}, | |
58 ## @var{A_ub}) may be set to the empty matrix (@code{[]}) if not | |
59 ## present. If the initial guess is feasible the algorithm is faster. | |
60 ## | |
61 ## The value @var{info} is a structure with the following fields: | |
62 ## @table @code | |
63 ## @item solveiter | |
64 ## The number of iterations required to find the solution. | |
65 ## @item info | |
66 ## An integer indicating the status of the solution, as follows: | |
67 ## @table @asis | |
68 ## @item 0 | |
69 ## The problem is feasible and convex. Global solution found. | |
70 ## @item 1 | |
71 ## The problem is not convex. Local solution found. | |
72 ## @item 2 | |
73 ## The problem is not convex and unbounded. | |
74 ## @item 3 | |
75 ## Maximum number of iterations reached. | |
76 ## @item 6 | |
77 ## The problem is infeasible. | |
78 ## @end table | |
79 ## @end table | |
80 ## @end deftypefn | |
81 | |
82 function [x, obj, INFO, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, A_in, A_ub) | |
83 | |
84 if (nargin == 5 || nargin == 7 || nargin == 10) | |
85 | |
86 ## Checking the quadratic penalty | |
87 n = issquare (H); | |
88 if (n == 0) | |
89 error ("qp: quadratic penalty matrix not square"); | |
90 endif | |
91 | |
92 n1 = issymmetric (H); | |
93 if (n1 == 0) | |
94 ## warning ("qp: quadratic penalty matrix not symmetric"); | |
95 H = (H + H')/2; | |
96 endif | |
97 | |
98 ## Checking the initial guess (if empty it is resized to the | |
99 ## right dimension and filled with 0) | |
100 if (isempty (x0)) | |
101 x0 = zeros (n, 1); | |
102 elseif (length (x0) != n) | |
103 error ("qp: the initial guess has incorrect length"); | |
104 endif | |
105 | |
106 ## Linear penalty. | |
107 if (length (q) != n) | |
108 error ("qp: the linear term has incorrect length"); | |
109 endif | |
110 | |
111 ## Equality constraint matrices | |
112 if (isempty (A) || isempty(b)) | |
113 n_eq = 0; | |
114 A = zeros (n_eq, n); | |
115 b = zeros (n_eq, 1); | |
116 else | |
117 [n_eq, n1] = size (A); | |
118 if (n1 != n) | |
119 error ("qp: equality constraint matrix has incorrect column dimension"); | |
120 endif | |
121 if (length (b) != n_eq) | |
122 error ("qp: equality constraint matrix and vector have inconsistent dimension"); | |
123 endif | |
124 endif | |
125 | |
126 ## Bound constraints | |
127 Ain = zeros (0, n); | |
128 bin = zeros (0, 1); | |
129 n_in = 0; | |
130 if (nargin > 5) | |
131 if (! isempty (lb)) | |
132 if (length(lb) != n) | |
133 error ("qp: lower bound has incorrect length"); | |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
134 elseif (isempty (ub)) |
5289 | 135 Ain = [Ain; eye(n)]; |
136 bin = [bin; lb]; | |
137 endif | |
138 endif | |
139 | |
140 if (! isempty (ub)) | |
141 if (length (ub) != n) | |
142 error ("qp: upper bound has incorrect length"); | |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
143 elseif (isempty (lb)) |
5289 | 144 Ain = [Ain; -eye(n)]; |
145 bin = [bin; -ub]; | |
146 endif | |
147 endif | |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
148 |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
149 if (! isempty (lb) && ! isempty (ub)) |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
150 rtol = sqrt (eps); |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
151 for i = 1:n |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
152 if (abs(lb (i) - ub(i)) < rtol*(1 + max (abs (lb(i) + ub(i))))) |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
153 ## These are actually an equality constraint |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
154 tmprow = zeros(1,n); |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
155 tmprow(i) = 1; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
156 A = [A;tmprow]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
157 b = [b; 0.5*(lb(i) + ub(i))]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
158 n_eq = n_eq + 1; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
159 else |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
160 tmprow = zeros(1,n); |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
161 tmprow(i) = 1; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
162 Ain = [Ain; tmprow; -tmprow]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
163 bin = [bin; lb(i); -ub(i)]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
164 n_in = n_in + 2; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
165 endif |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
166 endfor |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
167 endif |
5289 | 168 endif |
169 | |
170 ## Inequality constraints | |
171 if (nargin > 7) | |
172 [dimA_in, n1] = size (A_in); | |
173 if (n1 != n) | |
174 error ("qp: inequality constraint matrix has incorrect column dimension"); | |
175 else | |
176 if (! isempty (A_lb)) | |
177 if (length (A_lb) != dimA_in) | |
178 error ("qp: inequality constraint matrix and lower bound vector inconsistent"); | |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
179 elseif (isempty (A_ub)) |
5289 | 180 Ain = [Ain; A_in]; |
181 bin = [bin; A_lb]; | |
182 endif | |
183 endif | |
184 if (! isempty (A_ub)) | |
185 if (length (A_ub) != dimA_in) | |
186 error ("qp: inequality constraint matrix and upper bound vector inconsistent"); | |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
187 elseif (isempty (A_lb)) |
5289 | 188 Ain = [Ain; -A_in]; |
189 bin = [bin; -A_ub]; | |
190 endif | |
191 endif | |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
192 |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
193 if (! isempty (A_lb) && ! isempty (A_ub)) |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
194 rtol = sqrt (eps); |
8507 | 195 for i = 1:dimA_in |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
196 if (abs (A_lb(i) - A_ub(i)) < rtol*(1 + max (abs (A_lb(i) + A_ub(i))))) |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
197 ## These are actually an equality constraint |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
198 tmprow = A_in(i,:); |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
199 A = [A;tmprow]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
200 b = [b; 0.5*(A_lb(i) + A_ub(i))]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
201 n_eq = n_eq + 1; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
202 else |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
203 tmprow = A_in(i,:); |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
204 Ain = [Ain; tmprow; -tmprow]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
205 bin = [bin; A_lb(i); -A_ub(i)]; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
206 n_in = n_in + 2; |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
207 endif |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
208 endfor |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
209 endif |
5289 | 210 endif |
211 endif | |
212 | |
213 ## Now we should have the following QP: | |
214 ## | |
215 ## min_x 0.5*x'*H*x + x'*q | |
216 ## s.t. A*x = b | |
217 ## Ain*x >= bin | |
218 | |
219 ## Discard inequality constraints that have -Inf bounds since those | |
220 ## will never be active. | |
221 idx = isinf (bin) & bin < 0; | |
6523 | 222 |
6526 | 223 bin(idx) = []; |
224 Ain(idx,:) = []; | |
5289 | 225 |
5310 | 226 n_in = length (bin); |
227 | |
5289 | 228 ## Check if the initial guess is feasible. |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
229 if (isa (x0, "single") || isa (H, "single") || isa (q, "single") || isa (A, "single") |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
230 || isa (b, "single")) |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
231 rtol = sqrt (eps ("single")); |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
232 else |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
233 rtol = sqrt (eps); |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
234 endif |
5289 | 235 |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
236 eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+abs (b))); |
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
237 in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+abs (bin)))); |
5289 | 238 |
239 info = 0; | |
240 if (eq_infeasible || in_infeasible) | |
241 ## The initial guess is not feasible. | |
242 ## First define xbar that is feasible with respect to the equality | |
243 ## constraints. | |
244 if (eq_infeasible) | |
245 if (rank (A) < n_eq) | |
246 error ("qp: equality constraint matrix must be full row rank") | |
247 endif | |
248 xbar = pinv (A) * b; | |
249 else | |
250 xbar = x0; | |
251 endif | |
252 | |
253 ## Check if xbar is feasible with respect to the inequality | |
254 ## constraints also. | |
255 if (n_in > 0) | |
256 res = Ain * xbar - bin; | |
8280
5ee11a81688e
qp.m: convert b <= x <= b and b <= A*x <= b to equality constraints
Gabriele Pannocchia <g.pannocchia@ing.unipi.it>
parents:
7795
diff
changeset
|
257 if (any (res < -rtol * (1 + abs (bin)))) |
5289 | 258 ## xbar is not feasible with respect to the inequality |
259 ## constraints. Compute a step in the null space of the | |
260 ## equality constraints, by solving a QP. If the slack is | |
261 ## small, we have a feasible initial guess. Otherwise, the | |
262 ## problem is infeasible. | |
263 if (n_eq > 0) | |
264 Z = null (A); | |
265 if (isempty (Z)) | |
266 ## The problem is infeasible because A is square and full | |
267 ## rank, but xbar is not feasible. | |
268 info = 6; | |
269 endif | |
270 endif | |
271 | |
272 if (info != 6) | |
273 ## Solve an LP with additional slack variables to find | |
274 ## a feasible starting point. | |
275 gamma = eye (n_in); | |
276 if (n_eq > 0) | |
277 Atmp = [Ain*Z, gamma]; | |
278 btmp = -res; | |
279 else | |
280 Atmp = [Ain, gamma]; | |
281 btmp = bin; | |
282 endif | |
283 ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)]; | |
284 lb = [-Inf*ones(n-n_eq,1); zeros(n_in,1)]; | |
285 ub = []; | |
286 ctype = repmat ("L", n_in, 1); | |
287 [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype); | |
288 if ((status == 180 || status == 181 || status == 151) | |
289 && all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp)))) | |
290 ## We found a feasible starting point | |
291 if (n_eq > 0) | |
5306 | 292 x0 = xbar + Z*P(1:n-n_eq); |
5289 | 293 else |
294 x0 = P(1:n); | |
295 endif | |
296 else | |
297 ## The problem is infeasible | |
298 info = 6; | |
299 endif | |
300 endif | |
301 else | |
302 ## xbar is feasible. We use it a starting point. | |
303 x0 = xbar; | |
304 endif | |
305 else | |
306 ## xbar is feasible. We use it a starting point. | |
307 x0 = xbar; | |
308 endif | |
309 endif | |
310 | |
311 if (info == 0) | |
6848 | 312 ## FIXME -- make maxit a user option. |
5289 | 313 ## The initial (or computed) guess is feasible. |
314 ## We call the solver. | |
6848 | 315 maxit = 200; |
5289 | 316 [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit); |
317 else | |
318 iter = 0; | |
319 x = x0; | |
320 lambda = []; | |
321 endif | |
322 obj = 0.5 * x' * H * x + q' * x; | |
323 INFO.solveiter = iter; | |
324 INFO.info = info; | |
325 | |
326 else | |
6046 | 327 print_usage (); |
5289 | 328 endif |
329 | |
330 endfunction |