Mercurial > hg > octave-nkf

comparison scripts/optimization/qp.m @ 5289:b98bf1d70a0a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-21 14:43:14 by jwe]
author jwe
date Thu, 21 Apr 2005 14:43:20 +0000
parents
children 63cf9851f55c
comparison
equal deleted inserted replaced
5288:7434b60da83e 5289:b98bf1d70a0a
1 ## Copyright (C) 2000, 2001, 2004, 2005 Gabriele Pannocchia.
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
7 ## the Free Software Foundation; either version 2, or (at your option)
8 ## any later version.
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
16 ## along with Octave; see the file COPYING. If not, write to the Free
17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
18 ## 02111-1307, USA.
19
20 ## -*- texinfo -*-
21 ## @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{Ain}, @var{A_ub})
22 ## Solve the quadratic program
23 ## @ifinfo
24 ##
25 ## @example
26 ## min 0.5 x'*H*x + x'*q
27 ## x
28 ## @end example
29 ##
30 ## @end ifinfo
31 ## @iftex
32 ## @tex
33 ## @end tex
34 ## @end iftex
35 ## subject to
36 ## @ifinfo
37 ##
38 ## @example
39 ## A*x = b
40 ## lb <= x <= ub
41 ## A_lb <= Ain*x <= A_ub
42 ## @end example
43 ## @end ifinfo
44 ## @iftex
45 ## @tex
46 ## @end tex
47 ## @end iftex
48 ##
49 ## @noindent
50 ## using a null-space active-set method.
51 ##
52 ## Any bound (@var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb},
53 ## @var{A_ub}) may be set to the empty matrix (@code{[]}) if not
54 ## present. If the initial guess is feasible the algorithm is faster.
55 ##
56 ## The value @var{info} is a structure with the following fields:
57 ## @table @code
58 ## @item solveiter
59 ## The number of iterations required to find the solution.
60 ## @item info
61 ## An integer indicating the status of the solution, as follows:
62 ## @table @asis
63 ## @item 0
64 ## The problem is feasible and convex. Global solution found.
65 ## @item 1
66 ## The problem is not convex. Local solution found.
67 ## @item 2
68 ## The problem is not convex and unbounded.
69 ## @item 3
70 ## Maximum number of iterations reached.
71 ## @item 6
72 ## The problem is infeasible.
73 ## @end table
74 ## @end table
75 ## @end deftypefn
76
77 function [x, obj, INFO, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, A_in, A_ub)
78
79 if (nargin == 5 || nargin == 7 || nargin == 10)
80
81 ## Checking the quadratic penalty
82 n = issquare (H);
83 if (n == 0)
84 error ("qp: quadratic penalty matrix not square");
85 endif
86
87 n1 = issymmetric (H);
88 if (n1 == 0)
89 ## warning ("qp: quadratic penalty matrix not symmetric");
90 H = (H + H')/2;
91 endif
92
93 ## Checking the initial guess (if empty it is resized to the
94 ## right dimension and filled with 0)
95 if (isempty (x0))
96 x0 = zeros (n, 1);
97 elseif (length (x0) != n)
98 error ("qp: the initial guess has incorrect length");
99 endif
100
101 ## Linear penalty.
102 if (length (q) != n)
103 error ("qp: the linear term has incorrect length");
104 endif
105
106 ## Equality constraint matrices
107 if (isempty (A) || isempty(b))
108 n_eq = 0;
109 A = zeros (n_eq, n);
110 b = zeros (n_eq, 1);
111 else
112 [n_eq, n1] = size (A);
113 if (n1 != n)
114 error ("qp: equality constraint matrix has incorrect column dimension");
115 endif
116 if (length (b) != n_eq)
117 error ("qp: equality constraint matrix and vector have inconsistent dimension");
118 endif
119 endif
120
121 ## Bound constraints
122 Ain = zeros (0, n);
123 bin = zeros (0, 1);
124 n_in = 0;
125 if (nargin > 5)
126 if (! isempty (lb))
127 if (length(lb) != n)
128 error ("qp: lower bound has incorrect length");
129 else
130 Ain = [Ain; eye(n)];
131 bin = [bin; lb];
132 endif
133 endif
134
135 if (! isempty (ub))
136 if (length (ub) != n)
137 error ("qp: upper bound has incorrect length");
138 else
139 Ain = [Ain; -eye(n)];
140 bin = [bin; -ub];
141 endif
142 endif
143 endif
144
145 ## Inequality constraints
146 if (nargin > 7)
147 [dimA_in, n1] = size (A_in);
148 if (n1 != n)
149 error ("qp: inequality constraint matrix has incorrect column dimension");
150 else
151 if (! isempty (A_lb))
152 if (length (A_lb) != dimA_in)
153 error ("qp: inequality constraint matrix and lower bound vector inconsistent");
154 else
155 Ain = [Ain; A_in];
156 bin = [bin; A_lb];
157 endif
158 endif
159 if (! isempty (A_ub))
160 if (length (A_ub) != dimA_in)
161 error ("qp: inequality constraint matrix and upper bound vector inconsistent");
162 else
163 Ain = [Ain; -A_in];
164 bin = [bin; -A_ub];
165 endif
166 endif
167 endif
168 endif
169 n_in = length (bin);
170
171 ## Now we should have the following QP:
172 ##
173 ## min_x 0.5*x'*H*x + x'*q
174 ## s.t. A*x = b
175 ## Ain*x >= bin
176
177 ## Discard inequality constraints that have -Inf bounds since those
178 ## will never be active.
179 idx = isinf (bin) & bin < 0;
180 bin(idx) = [];
181 Ain(idx,:) = [];
182
183 ## Check if the initial guess is feasible.
184 rtol = sqrt (eps);
185
186 eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+norm (b)));
187 in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+norm (bin))));
188
189 info = 0;
190 if (eq_infeasible || in_infeasible)
191 ## The initial guess is not feasible.
192 ## First define xbar that is feasible with respect to the equality
193 ## constraints.
194 if (eq_infeasible)
195 if (rank (A) < n_eq)
196 error ("qp: equality constraint matrix must be full row rank")
197 endif
198 xbar = pinv (A) * b;
199 else
200 xbar = x0;
201 endif
202
203 ## Check if xbar is feasible with respect to the inequality
204 ## constraints also.
205 if (n_in > 0)
206 res = Ain * xbar - bin;
207 if (any (res < -rtol * (1 + norm (bin))))
208 ## xbar is not feasible with respect to the inequality
209 ## constraints. Compute a step in the null space of the
210 ## equality constraints, by solving a QP. If the slack is
211 ## small, we have a feasible initial guess. Otherwise, the
212 ## problem is infeasible.
213 if (n_eq > 0)
214 Z = null (A);
215 if (isempty (Z))
216 ## The problem is infeasible because A is square and full
217 ## rank, but xbar is not feasible.
218 info = 6;
219 endif
220 endif
221
222 if (info != 6)
223 ## Solve an LP with additional slack variables to find
224 ## a feasible starting point.
225 gamma = eye (n_in);
226 if (n_eq > 0)
227 Atmp = [Ain*Z, gamma];
228 btmp = -res;
229 else
230 Atmp = [Ain, gamma];
231 btmp = bin;
232 endif
233 ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)];
234 lb = [-Inf*ones(n-n_eq,1); zeros(n_in,1)];
235 ub = [];
236 ctype = repmat ("L", n_in, 1);
237 [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
238
239 if ((status == 180 || status == 181 || status == 151)
240 && all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp))))
241 ## We found a feasible starting point
242 if (n_eq > 0)
243 x0 = xbar + Z*P(1:n-n_eq)
244 else
245 x0 = P(1:n);
246 endif
247 else
248 ## The problem is infeasible
249 info = 6;
250 endif
251 endif
252 else
253 ## xbar is feasible. We use it a starting point.
254 x0 = xbar;
255 endif
256 else
257 ## xbar is feasible. We use it a starting point.
258 x0 = xbar;
259 endif
260 endif
261
262 if (info == 0)
263 ## The initial (or computed) guess is feasible.
264 ## We call the solver.
265 maxit = 100;
266 [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);
267 else
268 iter = 0;
269 x = x0;
270 lambda = [];
271 endif
272 obj = 0.5 * x' * H * x + q' * x;
273 INFO.solveiter = iter;
274 INFO.info = info;
275
276 else
277 usage ("[x, obj, info, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, Ain, A_ub)");
278 endif
279
280 endfunction