annotate scripts/optimization/qp.m @ 8920:eb63fbe60fab

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