Mercurial > octave

annotate scripts/optimization/qp.m @ 10059:665ad34efeed

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
qp.m: handle optimset options
author Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
date Tue, 05 Jan 2010 01:18:50 -0500
parents 72d6e0de76c7
children 8f51a90eb8d1
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
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
101 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
102
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
103 nargs = nargin;
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 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
106 options = varargin{end};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
107 nargs--;
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
108 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
109 options = struct ();
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
110 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
111
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
112 if (nargs >= 3)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
113 q = varargin{1};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
114 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
115 q = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
116 endif
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
117
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
118 if (nargs >= 5)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
119 A = varargin{2};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
120 b = varargin{3};
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 A = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
123 b = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
124 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
125
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
126 if (nargs >= 7)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
127 lb = varargin{4};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
128 ub = varargin{5};
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
129 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
130 lb = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
131 ub = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
132 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
133
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
134 if (nargs == 10)
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
135 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
136 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
137 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
138 else
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
139 A_lb = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
140 A_in = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
141 A_ub = [];
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
142 endif
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
143
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
144 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
145
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
146 maxit = optimget (options, "MaxIter", 200)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
147
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
148 ## Checking the quadratic penalty
9872
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
149 if (! issquare (H))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
150 error ("qp: quadratic penalty matrix not square");
9872
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
151 elseif (! ishermitian (H))
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
152 ## warning ("qp: quadratic penalty matrix not hermitian");
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
153 H = (H + H')/2;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
154 endif
9872
72d6e0de76c7 fix qp, condest and krylov
Jaroslav Hajek <highegg@gmail.com>
parents: 9211
diff changeset
155 n = rows (H);
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
156
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
157 ## 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
158 ## right dimension and filled with 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
159 if (isempty (x0))
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
160 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
161 elseif (numel (x0) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
162 error ("qp: the initial guess has incorrect length");
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
163 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
164
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
165 ## Linear penalty.
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
166 if (isempty (q))
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
167 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
168 elseif (numel (q) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
169 error ("qp: the linear term 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 ## 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
173 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
174 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
175 b = zeros (0, 1);
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
176 n_eq = 0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
177 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
178 [n_eq, n1] = size (A);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
179 if (n1 != n)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
180 error ("qp: equality constraint matrix has incorrect column dimension");
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
181 endif
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
182 if (numel (b) != n_eq)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
183 error ("qp: equality constraint matrix and vector have inconsistent dimension");
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 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
186
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
187 ## Bound constraints
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
188 Ain = zeros (0, n);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
189 bin = zeros (0, 1);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
190 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
191 if (nargs > 5)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
192 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
193 if (numel (lb) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
194 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
195 elseif (isempty (ub))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
196 Ain = [Ain; eye(n)];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
197 bin = [bin; lb];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
198 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
199 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
200
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
201 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
202 if (numel (ub) != n)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
203 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
204 elseif (isempty (lb))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
205 Ain = [Ain; -eye(n)];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
206 bin = [bin; -ub];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
207 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
208 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
209
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 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
211 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
212 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
213 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
214 ## 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
215 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
216 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
217 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
218 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
219 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
220 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
221 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
222 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
223 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
224 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
225 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
226 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
227 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
228 endif
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
229 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
230
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
231 ## Inequality constraints
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
232 if (nargs > 7)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
233 [dimA_in, n1] = size (A_in);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
234 if (n1 != n)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
235 error ("qp: inequality constraint matrix has incorrect column dimension");
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
236 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
237 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
238 if (numel (A_lb) != dimA_in)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
239 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
240 elseif (isempty (A_ub))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
241 Ain = [Ain; A_in];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
242 bin = [bin; A_lb];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
243 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
244 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
245 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
246 if (numel (A_ub) != dimA_in)
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
247 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
248 elseif (isempty (A_lb))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
249 Ain = [Ain; -A_in];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
250 bin = [bin; -A_ub];
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 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
253
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
254 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
255 rtol = sqrt (eps);
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8280
diff changeset
256 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
257 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
258 ## 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
259 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
260 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
261 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
262 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
263 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
264 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
265 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
266 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
267 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
268 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
269 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
270 endif
5289
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 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
273
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
274 ## Now we should have the following QP:
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
275 ##
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
276 ## min_x 0.5*x'*H*x + x'*q
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
277 ## s.t. A*x = b
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
278 ## Ain*x >= bin
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
279
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
280 ## Discard inequality constraints that have -Inf bounds since those
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
281 ## will never be active.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
282 idx = isinf (bin) & bin < 0;
6523
589c8dbba916 [project @ 2007-04-12 19:06:44 by jwe]
jwe
parents: 6417
diff changeset
283
6526
2a715c6409a5 [project @ 2007-04-13 20:31:50 by jwe]
jwe
parents: 6523
diff changeset
284 bin(idx) = [];
2a715c6409a5 [project @ 2007-04-13 20:31:50 by jwe]
jwe
parents: 6523
diff changeset
285 Ain(idx,:) = [];
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
286
10059
665ad34efeed qp.m: handle optimset options
Joshua Redstone <redstone@gmail.com>, John W. Eaton <jwe@octave.org>
parents: 9872
diff changeset
287 n_in = numel (bin);
5310
2fbcdc356fc7 [project @ 2005-04-27 15:38:53 by jwe]
jwe
parents: 5307
diff changeset
288
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
289 ## Check if the initial guess is feasible.
7795
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
290 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
291 || isa (b, "single"))
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
292 rtol = sqrt (eps ("single"));
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
293 else
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
294 rtol = sqrt (eps);
df9519e9990c Handle single precision eps values
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
295 endif
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
296
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
297 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
298 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
299
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
300 info = 0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
301 if (eq_infeasible || in_infeasible)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
302 ## The initial guess is not feasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
303 ## 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
304 ## constraints.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
305 if (eq_infeasible)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
306 if (rank (A) < n_eq)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
307 error ("qp: equality constraint matrix must be full row rank")
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
308 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
309 xbar = pinv (A) * b;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
310 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
311 xbar = x0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
312 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
313
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
314 ## Check if xbar is feasible with respect to the inequality
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
315 ## constraints also.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
316 if (n_in > 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
317 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
318 if (any (res < -rtol * (1 + abs (bin))))
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
319 ## xbar is not feasible with respect to the inequality
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
320 ## constraints. Compute a step in the null space of the
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
321 ## equality constraints, by solving a QP. If the slack is
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
322 ## small, we have a feasible initial guess. Otherwise, the
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
323 ## problem is infeasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
324 if (n_eq > 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
325 Z = null (A);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
326 if (isempty (Z))
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
327 ## The problem is infeasible because A is square and full
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
328 ## rank, but xbar is not feasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
329 info = 6;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
330 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
331 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
332
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
333 if (info != 6)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
334 ## Solve an LP with additional slack variables to find
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
335 ## a feasible starting point.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
336 gamma = eye (n_in);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
337 if (n_eq > 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
338 Atmp = [Ain*Z, gamma];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
339 btmp = -res;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
340 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
341 Atmp = [Ain, gamma];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
342 btmp = bin;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
343 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
344 ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
345 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
346 ub = [];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
347 ctype = repmat ("L", n_in, 1);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
348 [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
349 if ((status == 180 || status == 181 || status == 151)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
350 && 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
351 ## We found a feasible starting point
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
352 if (n_eq > 0)
5306
63cf9851f55c [project @ 2005-04-26 14:41:40 by jwe]
jwe
parents: 5289
diff changeset
353 x0 = xbar + Z*P(1:n-n_eq);
5289
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
354 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
355 x0 = P(1:n);
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
356 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
357 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
358 ## The problem is infeasible
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
359 info = 6;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
360 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
361 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
362 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
363 ## xbar is feasible. We use it a starting point.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
364 x0 = xbar;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
365 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
366 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
367 ## xbar is feasible. We use it a starting point.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
368 x0 = xbar;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
369 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
370 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
371
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
372 if (info == 0)
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
373 ## The initial (or computed) guess is feasible.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
374 ## We call the solver.
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
375 [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
376 else
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
377 iter = 0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
378 x = x0;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
379 lambda = [];
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
380 endif
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
381 obj = 0.5 * x' * H * x + q' * x;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
382 INFO.solveiter = iter;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
383 INFO.info = info;
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
384
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
385 else
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5310
diff changeset
386 print_usage ();
5289
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
b98bf1d70a0a [project @ 2005-04-21 14:43:14 by jwe]
jwe
parents:
diff changeset
389 endfunction