Mercurial > hg > octave-lyh

comparison 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
comparison
equal deleted inserted replaced
10058:1be07aac495d 10059:665ad34efeed
16 ## You should have received a copy of the GNU General Public License 16 ## You should have received a copy of the GNU General Public License
17 ## along with Octave; see the file COPYING. If not, see 17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>. 18 ## <http://www.gnu.org/licenses/>.
19 19
20 ## -*- texinfo -*- 20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb}, @var{A_in}, @var{A_ub}) 21 ## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H})
22 ## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q})
23 ## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b})
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})
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})
26 ## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@dots{}, @var{options})
22 ## Solve the quadratic program 27 ## Solve the quadratic program
23 ## @tex 28 ## @tex
24 ## $$ 29 ## $$
25 ## \min_x {1 \over 2} x^T H x + x^T q 30 ## \min_x {1 \over 2} x^T H x + x^T q
26 ## $$ 31 ## $$
56 ## using a null-space active-set method. 61 ## using a null-space active-set method.
57 ## 62 ##
58 ## Any bound (@var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb}, 63 ## Any bound (@var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb},
59 ## @var{A_ub}) may be set to the empty matrix (@code{[]}) if not 64 ## @var{A_ub}) may be set to the empty matrix (@code{[]}) if not
60 ## present. If the initial guess is feasible the algorithm is faster. 65 ## present. If the initial guess is feasible the algorithm is faster.
66 ##
67 ## @table @var
68 ## @item options
69 ## An optional structure containing the following
70 ## parameter(s) used to define the behavior of the solver. Missing elements
71 ## in the structure take on default values, so you only need to set the
72 ## elements that you wish to change from the default.
73 ##
74 ## @table @code
75 ## @item MaxIter (default: 200)
76 ## Maximum number of iterations.
77 ## @end table
78 ## @end table
61 ## 79 ##
62 ## The value @var{info} is a structure with the following fields: 80 ## The value @var{info} is a structure with the following fields:
63 ## @table @code 81 ## @table @code
64 ## @item solveiter 82 ## @item solveiter
65 ## The number of iterations required to find the solution. 83 ## The number of iterations required to find the solution.
78 ## The problem is infeasible. 96 ## The problem is infeasible.
79 ## @end table 97 ## @end table
80 ## @end table 98 ## @end table
81 ## @end deftypefn 99 ## @end deftypefn
82 100
83 function [x, obj, INFO, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, A_in, A_ub) 101 function [x, obj, INFO, lambda] = qp (x0, H, varargin)
84 102
85 if (nargin == 5 || nargin == 7 || nargin == 10) 103 nargs = nargin;
104
105 if (nargs > 2 && isstruct (varargin{end}))
106 options = varargin{end};
107 nargs--;
108 else
109 options = struct ();
110 endif
111
112 if (nargs >= 3)
113 q = varargin{1};
114 else
115 q = [];
116 endif
117
118 if (nargs >= 5)
119 A = varargin{2};
120 b = varargin{3};
121 else
122 A = [];
123 b = [];
124 endif
125
126 if (nargs >= 7)
127 lb = varargin{4};
128 ub = varargin{5};
129 else
130 lb = [];
131 ub = [];
132 endif
133
134 if (nargs == 10)
135 A_lb = varargin{6};
136 A_in = varargin{7};
137 A_ub = varargin{8};
138 else
139 A_lb = [];
140 A_in = [];
141 A_ub = [];
142 endif
143
144 if (nargs == 2 || nargs == 3 || nargs == 5 || nargs == 7 || nargs == 10)
145
146 maxit = optimget (options, "MaxIter", 200)
86 147
87 ## Checking the quadratic penalty 148 ## Checking the quadratic penalty
88 if (! issquare (H)) 149 if (! issquare (H))
89 error ("qp: quadratic penalty matrix not square"); 150 error ("qp: quadratic penalty matrix not square");
90 elseif (! ishermitian (H)) 151 elseif (! ishermitian (H))
95 156
96 ## Checking the initial guess (if empty it is resized to the 157 ## Checking the initial guess (if empty it is resized to the
97 ## right dimension and filled with 0) 158 ## right dimension and filled with 0)
98 if (isempty (x0)) 159 if (isempty (x0))
99 x0 = zeros (n, 1); 160 x0 = zeros (n, 1);
100 elseif (length (x0) != n) 161 elseif (numel (x0) != n)
101 error ("qp: the initial guess has incorrect length"); 162 error ("qp: the initial guess has incorrect length");
102 endif 163 endif
103 164
104 ## Linear penalty. 165 ## Linear penalty.
105 if (length (q) != n) 166 if (isempty (q))
167 q = zeros (n, 1);
168 elseif (numel (q) != n)
106 error ("qp: the linear term has incorrect length"); 169 error ("qp: the linear term has incorrect length");
107 endif 170 endif
108 171
109 ## Equality constraint matrices 172 ## Equality constraint matrices
110 if (isempty (A) || isempty(b)) 173 if (isempty (A) || isempty (b))
174 A = zeros (0, n);
175 b = zeros (0, 1);
111 n_eq = 0; 176 n_eq = 0;
112 A = zeros (n_eq, n);
113 b = zeros (n_eq, 1);
114 else 177 else
115 [n_eq, n1] = size (A); 178 [n_eq, n1] = size (A);
116 if (n1 != n) 179 if (n1 != n)
117 error ("qp: equality constraint matrix has incorrect column dimension"); 180 error ("qp: equality constraint matrix has incorrect column dimension");
118 endif 181 endif
119 if (length (b) != n_eq) 182 if (numel (b) != n_eq)
120 error ("qp: equality constraint matrix and vector have inconsistent dimension"); 183 error ("qp: equality constraint matrix and vector have inconsistent dimension");
121 endif 184 endif
122 endif 185 endif
123 186
124 ## Bound constraints 187 ## Bound constraints
125 Ain = zeros (0, n); 188 Ain = zeros (0, n);
126 bin = zeros (0, 1); 189 bin = zeros (0, 1);
127 n_in = 0; 190 n_in = 0;
128 if (nargin > 5) 191 if (nargs > 5)
129 if (! isempty (lb)) 192 if (! isempty (lb))
130 if (length(lb) != n) 193 if (numel (lb) != n)
131 error ("qp: lower bound has incorrect length"); 194 error ("qp: lower bound has incorrect length");
132 elseif (isempty (ub)) 195 elseif (isempty (ub))
133 Ain = [Ain; eye(n)]; 196 Ain = [Ain; eye(n)];
134 bin = [bin; lb]; 197 bin = [bin; lb];
135 endif 198 endif
136 endif 199 endif
137 200
138 if (! isempty (ub)) 201 if (! isempty (ub))
139 if (length (ub) != n) 202 if (numel (ub) != n)
140 error ("qp: upper bound has incorrect length"); 203 error ("qp: upper bound has incorrect length");
141 elseif (isempty (lb)) 204 elseif (isempty (lb))
142 Ain = [Ain; -eye(n)]; 205 Ain = [Ain; -eye(n)];
143 bin = [bin; -ub]; 206 bin = [bin; -ub];
144 endif 207 endif
164 endfor 227 endfor
165 endif 228 endif
166 endif 229 endif
167 230
168 ## Inequality constraints 231 ## Inequality constraints
169 if (nargin > 7) 232 if (nargs > 7)
170 [dimA_in, n1] = size (A_in); 233 [dimA_in, n1] = size (A_in);
171 if (n1 != n) 234 if (n1 != n)
172 error ("qp: inequality constraint matrix has incorrect column dimension"); 235 error ("qp: inequality constraint matrix has incorrect column dimension");
173 else 236 else
174 if (! isempty (A_lb)) 237 if (! isempty (A_lb))
175 if (length (A_lb) != dimA_in) 238 if (numel (A_lb) != dimA_in)
176 error ("qp: inequality constraint matrix and lower bound vector inconsistent"); 239 error ("qp: inequality constraint matrix and lower bound vector inconsistent");
177 elseif (isempty (A_ub)) 240 elseif (isempty (A_ub))
178 Ain = [Ain; A_in]; 241 Ain = [Ain; A_in];
179 bin = [bin; A_lb]; 242 bin = [bin; A_lb];
180 endif 243 endif
181 endif 244 endif
182 if (! isempty (A_ub)) 245 if (! isempty (A_ub))
183 if (length (A_ub) != dimA_in) 246 if (numel (A_ub) != dimA_in)
184 error ("qp: inequality constraint matrix and upper bound vector inconsistent"); 247 error ("qp: inequality constraint matrix and upper bound vector inconsistent");
185 elseif (isempty (A_lb)) 248 elseif (isempty (A_lb))
186 Ain = [Ain; -A_in]; 249 Ain = [Ain; -A_in];
187 bin = [bin; -A_ub]; 250 bin = [bin; -A_ub];
188 endif 251 endif
219 idx = isinf (bin) & bin < 0; 282 idx = isinf (bin) & bin < 0;
220 283
221 bin(idx) = []; 284 bin(idx) = [];
222 Ain(idx,:) = []; 285 Ain(idx,:) = [];
223 286
224 n_in = length (bin); 287 n_in = numel (bin);
225 288
226 ## Check if the initial guess is feasible. 289 ## Check if the initial guess is feasible.
227 if (isa (x0, "single") || isa (H, "single") || isa (q, "single") || isa (A, "single") 290 if (isa (x0, "single") || isa (H, "single") || isa (q, "single") || isa (A, "single")
228 || isa (b, "single")) 291 || isa (b, "single"))
229 rtol = sqrt (eps ("single")); 292 rtol = sqrt (eps ("single"));
305 x0 = xbar; 368 x0 = xbar;
306 endif 369 endif
307 endif 370 endif
308 371
309 if (info == 0) 372 if (info == 0)
310 ## FIXME -- make maxit a user option.
311 ## The initial (or computed) guess is feasible. 373 ## The initial (or computed) guess is feasible.
312 ## We call the solver. 374 ## We call the solver.
313 maxit = 200;
314 [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit); 375 [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);
315 else 376 else
316 iter = 0; 377 iter = 0;
317 x = x0; 378 x = x0;
318 lambda = []; 379 lambda = [];