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