Mercurial > hg > octave-nkf

annotate scripts/optimization/qp.m @ 10060:8f51a90eb8d1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement default opts query and register opts for qp
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 05 Jan 2010 07:28:42 +0100
parents 665ad34efeed
children 17ce2a700a97
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 -*-
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
21 ## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H})
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
22 ## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q})
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
23 ## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b})
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
24 ## @deftypefnx {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})
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
25 ## @deftypefnx {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})
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
26 ## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@dots{}, @var{options})
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
27 ## Solve the quadratic program
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
28 ## @tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
29 ## $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
30 ## \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
31 ## $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
32 ## @end tex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
33 ## @ifnottex
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
34 ##
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
35 ## @example
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
36 ## @group
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
37 ## min 0.5 x'*H*x + x'*q
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
38 ## x
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
39 ## @end group
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
40 ## @end example
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
41 ##
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
42 ## @end ifnottex
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
43 ## subject to
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
44 ## @tex
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
45 ## $$
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
46 ## 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
47 ## $$
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
48 ## @end tex
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
49 ## @ifnottex
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
50 ##
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
51 ## @example
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
52 ## @group
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
53 ## A*x = b
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
54 ## lb <= x <= ub
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
55 ## 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
56 ## @end group
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
57 ## @end example
6741
00116015904d [project @ 2007-06-18 16:07:14 by jwe]
jwe
parents: 6526
diff changeset
58 ## @end ifnottex
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
59 ##
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
60 ## @noindent
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
61 ## using a null-space active-set method.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
62 ##
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
63 ## 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
64 ## @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
65 ## 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
66 ##
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
67 ## @table @var
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
68 ## @item options
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
69 ## An optional structure containing the following
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
70 ## parameter(s) used to define the behavior of the solver. Missing elements
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
71 ## in the structure take on default values, so you only need to set the
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
72 ## elements that you wish to change from the default.
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
73 ##
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
74 ## @table @code
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
75 ## @item MaxIter (default: 200)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
76 ## Maximum number of iterations.
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
77 ## @end table
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
78 ## @end table
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
79 ##
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
80 ## 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
81 ## @table @code
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
82 ## @item solveiter
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
83 ## The number of iterations required to find the solution.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
84 ## @item info
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
85 ## An integer indicating the status of the solution, as follows:
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
86 ## @table @asis
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
87 ## @item 0
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
88 ## The problem is feasible and convex. Global solution found.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
89 ## @item 1
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
90 ## The problem is not convex. Local solution found.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
91 ## @item 2
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
92 ## The problem is not convex and unbounded.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
93 ## @item 3
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
94 ## Maximum number of iterations reached.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
95 ## @item 6
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
96 ## The problem is infeasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
97 ## @end table
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
98 ## @end table
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
99 ## @end deftypefn
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
100
10060
8f51a90eb8d1 implement default opts query and register opts for qp
Jaroslav Hajek <highegg@gmail.com>
parents: 10059
diff changeset
101 ## PKG_ADD: __all_opts__ ("qp");
8f51a90eb8d1 implement default opts query and register opts for qp
Jaroslav Hajek <highegg@gmail.com>
parents: 10059
diff changeset
102
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
103 function [x, obj, INFO, lambda] = qp (x0, H, varargin)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
104
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
105 nargs = nargin;
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
106
10060
8f51a90eb8d1 implement default opts query and register opts for qp
Jaroslav Hajek <highegg@gmail.com>
parents: 10059
diff changeset
107 if (nargin == 1 && ischar (x0) && strcmp (x0, 'defaults'))
8f51a90eb8d1 implement default opts query and register opts for qp
Jaroslav Hajek <highegg@gmail.com>
parents: 10059
diff changeset
108 x = optimset ("MaxIter", 200);
8f51a90eb8d1 implement default opts query and register opts for qp
Jaroslav Hajek <highegg@gmail.com>
parents: 10059
diff changeset
109 return;
8f51a90eb8d1 implement default opts query and register opts for qp
Jaroslav Hajek <highegg@gmail.com>
parents: 10059
diff changeset
110 endif
8f51a90eb8d1 implement default opts query and register opts for qp
Jaroslav Hajek <highegg@gmail.com>
parents: 10059
diff changeset
111
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
112 if (nargs > 2 && isstruct (varargin{end}))
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
113 options = varargin{end};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
114 nargs--;
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
115 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
116 options = struct ();
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
117 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
118
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
119 if (nargs >= 3)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
120 q = varargin{1};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
121 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
122 q = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
123 endif
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
124
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
125 if (nargs >= 5)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
126 A = varargin{2};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
127 b = varargin{3};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
128 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
129 A = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
130 b = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
131 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
132
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
133 if (nargs >= 7)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
134 lb = varargin{4};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
135 ub = varargin{5};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
136 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
137 lb = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
138 ub = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
139 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
140
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
141 if (nargs == 10)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
142 A_lb = varargin{6};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
143 A_in = varargin{7};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
144 A_ub = varargin{8};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
145 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
146 A_lb = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
147 A_in = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
148 A_ub = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
149 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
150
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
151 if (nargs == 2 || nargs == 3 || nargs == 5 || nargs == 7 || nargs == 10)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
152
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
153 maxit = optimget (options, "MaxIter", 200)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
154
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
155 ## Checking the quadratic penalty
9872
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
156 if (! issquare (H))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
157 error ("qp: quadratic penalty matrix not square");
9872
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
158 elseif (! ishermitian (H))
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
159 ## warning ("qp: quadratic penalty matrix not hermitian");
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
160 H = (H + H')/2;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
161 endif
9872
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
162 n = rows (H);
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
163
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
164 ## 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
165 ## right dimension and filled with 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
166 if (isempty (x0))
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
167 x0 = zeros (n, 1);
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
168 elseif (numel (x0) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
169 error ("qp: the initial guess has incorrect length");
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
170 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
171
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
172 ## Linear penalty.
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
173 if (isempty (q))
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
174 q = zeros (n, 1);
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
175 elseif (numel (q) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
176 error ("qp: the linear term has incorrect length");
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
177 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
178
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
179 ## Equality constraint matrices
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
180 if (isempty (A) || isempty (b))
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
181 A = zeros (0, n);
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
182 b = zeros (0, 1);
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
183 n_eq = 0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
184 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
185 [n_eq, n1] = size (A);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
186 if (n1 != n)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
187 error ("qp: equality constraint matrix has incorrect column dimension");
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
188 endif
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
189 if (numel (b) != n_eq)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
190 error ("qp: equality constraint matrix and vector have inconsistent dimension");
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
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
193
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
194 ## Bound constraints
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
195 Ain = zeros (0, n);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
196 bin = zeros (0, 1);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
197 n_in = 0;
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
198 if (nargs > 5)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
199 if (! isempty (lb))
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
200 if (numel (lb) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
201 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
202 elseif (isempty (ub))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
203 Ain = [Ain; eye(n)];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
204 bin = [bin; lb];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
205 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
206 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
207
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
208 if (! isempty (ub))
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
209 if (numel (ub) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
210 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
211 elseif (isempty (lb))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
212 Ain = [Ain; -eye(n)];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
213 bin = [bin; -ub];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
214 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
215 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
216
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
217 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
218 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
219 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
220 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
221 ## 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
222 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
223 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
224 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
225 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
226 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
227 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
228 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
229 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
230 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
231 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
232 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
233 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
234 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
235 endif
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
236 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
237
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
238 ## Inequality constraints
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
239 if (nargs > 7)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
240 [dimA_in, n1] = size (A_in);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
241 if (n1 != n)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
242 error ("qp: inequality constraint matrix has incorrect column dimension");
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
243 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
244 if (! isempty (A_lb))
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
245 if (numel (A_lb) != dimA_in)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
246 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
247 elseif (isempty (A_ub))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
248 Ain = [Ain; A_in];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
249 bin = [bin; A_lb];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
250 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
251 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
252 if (! isempty (A_ub))
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
253 if (numel (A_ub) != dimA_in)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
254 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
255 elseif (isempty (A_lb))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
256 Ain = [Ain; -A_in];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
257 bin = [bin; -A_ub];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
258 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
259 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
260
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
261 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
262 rtol = sqrt (eps);
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8280
diff changeset
263 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
264 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
265 ## 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
266 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
267 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
268 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
269 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
270 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
271 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
272 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
273 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
274 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
275 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
276 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
277 endif
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
278 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
279 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
280
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
281 ## Now we should have the following QP:
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
282 ##
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
283 ## min_x 0.5*x'*H*x + x'*q
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
284 ## s.t. A*x = b
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
285 ## Ain*x >= bin
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
286
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
287 ## Discard inequality constraints that have -Inf bounds since those
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
288 ## will never be active.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
289 idx = isinf (bin) & bin < 0;
6523
589c8dbba916 [project @ 2007-04-12 19:06:44 by jwe]
jwe
parents: 6417
diff changeset
290
6526
2a715c6409a5 [project @ 2007-04-13 20:31:50 by jwe]
jwe
parents: 6523
diff changeset
291 bin(idx) = [];
2a715c6409a5 [project @ 2007-04-13 20:31:50 by jwe]
jwe
parents: 6523
diff changeset
292 Ain(idx,:) = [];
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
293
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
294 n_in = numel (bin);
5310
2fbcdc356fc7 [project @ 2005-04-27 15:38:53 by jwe]
jwe
parents: 5307
diff changeset
295
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
296 ## Check if the initial guess is feasible.
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
297 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
298 || isa (b, "single"))
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
299 rtol = sqrt (eps ("single"));
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
300 else
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
301 rtol = sqrt (eps);
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
302 endif
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
303
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
304 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
305 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
306
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
307 info = 0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
308 if (eq_infeasible || in_infeasible)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
309 ## The initial guess is not feasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
310 ## 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
311 ## constraints.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
312 if (eq_infeasible)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
313 if (rank (A) < n_eq)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
314 error ("qp: equality constraint matrix must be full row rank")
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
315 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
316 xbar = pinv (A) * b;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
317 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
318 xbar = x0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
319 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
320
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
321 ## Check if xbar is feasible with respect to the inequality
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
322 ## constraints also.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
323 if (n_in > 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
324 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
325 if (any (res < -rtol * (1 + abs (bin))))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
326 ## xbar is not feasible with respect to the inequality
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
327 ## constraints. Compute a step in the null space of the
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
328 ## equality constraints, by solving a QP. If the slack is
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
329 ## small, we have a feasible initial guess. Otherwise, the
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
330 ## problem is infeasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
331 if (n_eq > 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
332 Z = null (A);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
333 if (isempty (Z))
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
334 ## The problem is infeasible because A is square and full
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
335 ## rank, but xbar is not feasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
336 info = 6;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
337 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
338 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
339
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
340 if (info != 6)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
341 ## Solve an LP with additional slack variables to find
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
342 ## a feasible starting point.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
343 gamma = eye (n_in);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
344 if (n_eq > 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
345 Atmp = [Ain*Z, gamma];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
346 btmp = -res;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
347 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
348 Atmp = [Ain, gamma];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
349 btmp = bin;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
350 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
351 ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
352 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
353 ub = [];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
354 ctype = repmat ("L", n_in, 1);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
355 [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
356 if ((status == 180 || status == 181 || status == 151)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
357 && 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
358 ## We found a feasible starting point
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
359 if (n_eq > 0)
5306
63cf9851f55c [project @ 2005-04-26 14:41:40 by jwe]
jwe
parents: 5289
diff changeset
360 x0 = xbar + Z*P(1:n-n_eq);
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
361 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
362 x0 = P(1:n);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
363 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
364 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
365 ## The problem is infeasible
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
366 info = 6;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
367 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
368 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
369 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
370 ## xbar is feasible. We use it a starting point.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
371 x0 = xbar;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
372 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
373 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
374 ## xbar is feasible. We use it a starting point.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
375 x0 = xbar;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
376 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
377 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
378
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
379 if (info == 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
380 ## The initial (or computed) guess is feasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
381 ## We call the solver.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
382 [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
383 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
384 iter = 0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
385 x = x0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
386 lambda = [];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
387 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
388 obj = 0.5 * x' * H * x + q' * x;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
389 INFO.solveiter = iter;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
390 INFO.info = info;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
391
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
392 else
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5310
diff changeset
393 print_usage ();
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
394 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
395
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
396 endfunction