Mercurial > hg > octave-lyh
annotate scripts/optimization/glpk.m @ 8306:43795cf108d0
initial implementation of fsolve
remove old fsolve code
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Sun, 28 Sep 2008 21:09:35 +0200 |
parents | 854683691d7a |
children | eb63fbe60fab |
rev | line source |
---|---|
7017 | 1 ## Copyright (C) 2005, 2006, 2007 Nicolo' Giorgetti |
5233 | 2 ## |
3 ## This file is part of Octave. | |
4 ## | |
5 ## Octave is free software; you can redistribute it and/or modify it | |
6 ## under the terms of the GNU General Public License as published by | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
5233 | 9 ## |
10 ## Octave is distributed in the hope that it will be useful, but | |
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 ## General Public License for more details. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
5233 | 18 |
5244 | 19 ## -*- texinfo -*- |
20 ## @deftypefn {Function File} {[@var{xopt}, @var{fmin}, @var{status}, @var{extra}] =} glpk (@var{c}, @var{a}, @var{b}, @var{lb}, @var{ub}, @var{ctype}, @var{vartype}, @var{sense}, @var{param}) | |
21 ## Solve a linear program using the GNU GLPK library. Given three | |
22 ## arguments, @code{glpk} solves the following standard LP: | |
23 ## | |
6741 | 24 ## @iftex |
25 ## @tex | |
26 ## $$ | |
27 ## \min_x C^T x | |
28 ## $$ | |
29 ## @end tex | |
30 ## @end iftex | |
31 ## @ifnottex | |
5244 | 32 ## @example |
33 ## min C'*x | |
34 ## @end example | |
6741 | 35 ## @end ifnottex |
5244 | 36 ## |
5237 | 37 ## subject to |
5244 | 38 ## |
6741 | 39 ## @iftex |
40 ## @tex | |
41 ## $$ | |
42 ## Ax = b \qquad x \geq 0 | |
43 ## $$ | |
44 ## @end tex | |
45 ## @end iftex | |
46 ## @ifnottex | |
5244 | 47 ## @example |
48 ## @group | |
49 ## A*x = b | |
50 ## x >= 0 | |
51 ## @end group | |
52 ## @end example | |
6741 | 53 ## @end ifnottex |
5244 | 54 ## |
5237 | 55 ## but may also solve problems of the form |
5244 | 56 ## |
6741 | 57 ## @iftex |
58 ## @tex | |
59 ## $$ | |
60 ## [ \min_x | \max_x ] C^T x | |
61 ## $$ | |
62 ## @end tex | |
63 ## @end iftex | |
64 ## @ifnottex | |
5244 | 65 ## @example |
66 ## [ min | max ] C'*x | |
67 ## @end example | |
6741 | 68 ## @end ifnottex |
5244 | 69 ## |
5237 | 70 ## subject to |
5244 | 71 ## |
6741 | 72 ## @iftex |
73 ## @tex | |
74 ## $$ | |
75 ## Ax [ = | \leq | \geq ] b \qquad LB \leq x \leq UB | |
76 ## $$ | |
77 ## @end tex | |
78 ## @end iftex | |
79 ## @ifnottex | |
5244 | 80 ## @example |
81 ## @group | |
82 ## A*x [ "=" | "<=" | ">=" ] b | |
83 ## x >= LB | |
84 ## x <= UB | |
85 ## @end group | |
86 ## @end example | |
6741 | 87 ## @end ifnottex |
5244 | 88 ## |
89 ## Input arguments: | |
90 ## | |
91 ## @table @var | |
92 ## @item c | |
93 ## A column array containing the objective function coefficients. | |
94 ## | |
95 ## @item a | |
96 ## A matrix containing the constraints coefficients. | |
97 ## | |
98 ## @item b | |
99 ## A column array containing the right-hand side value for each constraint | |
100 ## in the constraint matrix. | |
101 ## | |
102 ## @item lb | |
103 ## An array containing the lower bound on each of the variables. If | |
104 ## @var{lb} is not supplied, the default lower bound for the variables is | |
105 ## zero. | |
106 ## | |
107 ## @item ub | |
108 ## An array containing the upper bound on each of the variables. If | |
109 ## @var{ub} is not supplied, the default upper bound is assumed to be | |
110 ## infinite. | |
111 ## | |
112 ## @item ctype | |
113 ## An array of characters containing the sense of each constraint in the | |
114 ## constraint matrix. Each element of the array may be one of the | |
115 ## following values | |
116 ## @table @code | |
117 ## @item "F" | |
7096 | 118 ## A free (unbounded) constraint (the constraint is ignored). |
5244 | 119 ## @item "U" |
7096 | 120 ## An inequality constraint with an upper bound (@code{A(i,:)*x <= b(i)}). |
5244 | 121 ## @item "S" |
7096 | 122 ## An equality constraint (@code{A(i,:)*x = b(i)}). |
5244 | 123 ## @item "L" |
7096 | 124 ## An inequality with a lower bound (@code{A(i,:)*x >= b(i)}). |
5244 | 125 ## @item "D" |
7096 | 126 ## An inequality constraint with both upper and lower bounds |
127 ## (@code{A(i,:)*x >= -b(i)} @emph{and} (@code{A(i,:)*x <= b(i)}). | |
5244 | 128 ## @end table |
129 ## | |
130 ## @item vartype | |
131 ## A column array containing the types of the variables. | |
132 ## @table @code | |
7096 | 133 ## @item "C" |
134 ## A continuous variable. | |
135 ## @item "I" | |
136 ## An integer variable. | |
5244 | 137 ## @end table |
138 ## | |
139 ## @item sense | |
140 ## If @var{sense} is 1, the problem is a minimization. If @var{sense} is | |
141 ## -1, the problem is a maximization. The default value is 1. | |
142 ## | |
143 ## @item param | |
144 ## A structure containing the following parameters used to define the | |
145 ## behavior of solver. Missing elements in the structure take on default | |
146 ## values, so you only need to set the elements that you wish to change | |
147 ## from the default. | |
148 ## | |
149 ## Integer parameters: | |
150 ## | |
151 ## @table @code | |
152 ## @item msglev (@code{LPX_K_MSGLEV}, default: 1) | |
153 ## Level of messages output by solver routines: | |
154 ## @table @asis | |
155 ## @item 0 | |
156 ## No output. | |
157 ## @item 1 | |
158 ## Error messages only. | |
159 ## @item 2 | |
160 ## Normal output . | |
161 ## @item 3 | |
162 ## Full output (includes informational messages). | |
163 ## @end table | |
164 ## | |
165 ## @item scale (@code{LPX_K_SCALE}, default: 1) | |
166 ## Scaling option: | |
167 ## @table @asis | |
168 ## @item 0 | |
7001 | 169 ## No scaling. |
5244 | 170 ## @item 1 |
7001 | 171 ## Equilibration scaling. |
5244 | 172 ## @item 2 |
173 ## Geometric mean scaling, then equilibration scaling. | |
174 ## @end table | |
175 ## | |
176 ## @item dual (@code{LPX_K_DUAL}, default: 0) | |
177 ## Dual simplex option: | |
178 ## @table @asis | |
179 ## @item 0 | |
180 ## Do not use the dual simplex. | |
181 ## @item 1 | |
182 ## If initial basic solution is dual feasible, use the dual simplex. | |
183 ## @end table | |
184 ## | |
185 ## @item price (@code{LPX_K_PRICE}, default: 1) | |
186 ## Pricing option (for both primal and dual simplex): | |
187 ## @table @asis | |
188 ## @item 0 | |
189 ## Textbook pricing. | |
190 ## @item 1 | |
191 ## Steepest edge pricing. | |
192 ## @end table | |
193 ## | |
194 ## @item round (@code{LPX_K_ROUND}, default: 0) | |
195 ## Solution rounding option: | |
196 ## @table @asis | |
197 ## @item 0 | |
198 ## Report all primal and dual values "as is". | |
199 ## @item 1 | |
200 ## Replace tiny primal and dual values by exact zero. | |
201 ## @end table | |
202 ## | |
203 ## @item itlim (@code{LPX_K_ITLIM}, default: -1) | |
204 ## Simplex iterations limit. If this value is positive, it is decreased by | |
205 ## one each time when one simplex iteration has been performed, and | |
206 ## reaching zero value signals the solver to stop the search. Negative | |
207 ## value means no iterations limit. | |
208 ## | |
209 ## @item itcnt (@code{LPX_K_OUTFRQ}, default: 200) | |
210 ## Output frequency, in iterations. This parameter specifies how | |
211 ## frequently the solver sends information about the solution to the | |
212 ## standard output. | |
213 ## | |
214 ## @item branch (@code{LPX_K_BRANCH}, default: 2) | |
215 ## Branching heuristic option (for MIP only): | |
216 ## @table @asis | |
217 ## @item 0 | |
218 ## Branch on the first variable. | |
219 ## @item 1 | |
220 ## Branch on the last variable. | |
221 ## @item 2 | |
222 ## Branch using a heuristic by Driebeck and Tomlin. | |
223 ## @end table | |
224 ## | |
225 ## @item btrack (@code{LPX_K_BTRACK}, default: 2) | |
226 ## Backtracking heuristic option (for MIP only): | |
227 ## @table @asis | |
228 ## @item 0 | |
229 ## Depth first search. | |
230 ## @item 1 | |
231 ## Breadth first search. | |
232 ## @item 2 | |
233 ## Backtrack using the best projection heuristic. | |
234 ## @end table | |
235 ## | |
236 ## @item presol (@code{LPX_K_PRESOL}, default: 1) | |
237 ## If this flag is set, the routine lpx_simplex solves the problem using | |
238 ## the built-in LP presolver. Otherwise the LP presolver is not used. | |
239 ## | |
240 ## @item lpsolver (default: 1) | |
241 ## Select which solver to use. If the problem is a MIP problem this flag | |
242 ## will be ignored. | |
243 ## @table @asis | |
244 ## @item 1 | |
245 ## Revised simplex method. | |
246 ## @item 2 | |
247 ## Interior point method. | |
248 ## @end table | |
249 ## @item save (default: 0) | |
7001 | 250 ## If this parameter is nonzero, save a copy of the problem in |
5244 | 251 ## CPLEX LP format to the file @file{"outpb.lp"}. There is currently no |
252 ## way to change the name of the output file. | |
253 ## @end table | |
254 ## | |
255 ## Real parameters: | |
256 ## | |
257 ## @table @code | |
258 ## @item relax (@code{LPX_K_RELAX}, default: 0.07) | |
259 ## Relaxation parameter used in the ratio test. If it is zero, the textbook | |
260 ## ratio test is used. If it is non-zero (should be positive), Harris' | |
261 ## two-pass ratio test is used. In the latter case on the first pass of the | |
262 ## ratio test basic variables (in the case of primal simplex) or reduced | |
263 ## costs of non-basic variables (in the case of dual simplex) are allowed | |
264 ## to slightly violate their bounds, but not more than | |
265 ## @code{relax*tolbnd} or @code{relax*toldj (thus, @code{relax} is a | |
6248 | 266 ## percentage of @code{tolbnd} or @code{toldj}}. |
5244 | 267 ## |
268 ## @item tolbnd (@code{LPX_K_TOLBND}, default: 10e-7) | |
5289 | 269 ## Relative tolerance used to check if the current basic solution is primal |
5244 | 270 ## feasible. It is not recommended that you change this parameter unless you |
271 ## have a detailed understanding of its purpose. | |
272 ## | |
273 ## @item toldj (@code{LPX_K_TOLDJ}, default: 10e-7) | |
274 ## Absolute tolerance used to check if the current basic solution is dual | |
275 ## feasible. It is not recommended that you change this parameter unless you | |
276 ## have a detailed understanding of its purpose. | |
277 ## | |
278 ## @item tolpiv (@code{LPX_K_TOLPIV}, default: 10e-9) | |
279 ## Relative tolerance used to choose eligible pivotal elements of the | |
280 ## simplex table. It is not recommended that you change this parameter unless you | |
281 ## have a detailed understanding of its purpose. | |
282 ## | |
5289 | 283 ## @item objll (@code{LPX_K_OBJLL}, default: -DBL_MAX) |
5244 | 284 ## Lower limit of the objective function. If on the phase II the objective |
285 ## function reaches this limit and continues decreasing, the solver stops | |
286 ## the search. This parameter is used in the dual simplex method only. | |
287 ## | |
288 ## @item objul (@code{LPX_K_OBJUL}, default: +DBL_MAX) | |
289 ## Upper limit of the objective function. If on the phase II the objective | |
290 ## function reaches this limit and continues increasing, the solver stops | |
291 ## the search. This parameter is used in the dual simplex only. | |
292 ## | |
293 ## @item tmlim (@code{LPX_K_TMLIM}, default: -1.0) | |
294 ## Searching time limit, in seconds. If this value is positive, it is | |
295 ## decreased each time when one simplex iteration has been performed by the | |
296 ## amount of time spent for the iteration, and reaching zero value signals | |
297 ## the solver to stop the search. Negative value means no time limit. | |
298 ## | |
299 ## @item outdly (@code{LPX_K_OUTDLY}, default: 0.0) | |
300 ## Output delay, in seconds. This parameter specifies how long the solver | |
301 ## should delay sending information about the solution to the standard | |
302 ## output. Non-positive value means no delay. | |
303 ## | |
304 ## @item tolint (@code{LPX_K_TOLINT}, default: 10e-5) | |
5289 | 305 ## Relative tolerance used to check if the current basic solution is integer |
5244 | 306 ## feasible. It is not recommended that you change this parameter unless |
307 ## you have a detailed understanding of its purpose. | |
308 ## | |
309 ## @item tolobj (@code{LPX_K_TOLOBJ}, default: 10e-7) | |
310 ## Relative tolerance used to check if the value of the objective function | |
311 ## is not better than in the best known integer feasible solution. It is | |
312 ## not recommended that you change this parameter unless you have a | |
313 ## detailed understanding of its purpose. | |
314 ## @end table | |
315 ## @end table | |
316 ## | |
317 ## Output values: | |
318 ## | |
319 ## @table @var | |
320 ## @item xopt | |
321 ## The optimizer (the value of the decision variables at the optimum). | |
322 ## @item fopt | |
323 ## The optimum value of the objective function. | |
324 ## @item status | |
325 ## Status of the optimization. | |
326 ## | |
327 ## Simplex Method: | |
328 ## @table @asis | |
329 ## @item 180 (@code{LPX_OPT}) | |
330 ## Solution is optimal. | |
331 ## @item 181 (@code{LPX_FEAS}) | |
332 ## Solution is feasible. | |
333 ## @item 182 (@code{LPX_INFEAS}) | |
334 ## Solution is infeasible. | |
335 ## @item 183 (@code{LPX_NOFEAS}) | |
336 ## Problem has no feasible solution. | |
337 ## @item 184 (@code{LPX_UNBND}) | |
338 ## Problem has no unbounded solution. | |
339 ## @item 185 (@code{LPX_UNDEF}) | |
340 ## Solution status is undefined. | |
341 ## @end table | |
342 ## Interior Point Method: | |
343 ## @table @asis | |
344 ## @item 150 (@code{LPX_T_UNDEF}) | |
345 ## The interior point method is undefined. | |
346 ## @item 151 (@code{LPX_T_OPT}) | |
347 ## The interior point method is optimal. | |
348 ## @end table | |
349 ## Mixed Integer Method: | |
350 ## @table @asis | |
351 ## @item 170 (@code{LPX_I_UNDEF}) | |
352 ## The status is undefined. | |
353 ## @item 171 (@code{LPX_I_OPT}) | |
354 ## The solution is integer optimal. | |
355 ## @item 172 (@code{LPX_I_FEAS}) | |
356 ## Solution integer feasible but its optimality has not been proven | |
357 ## @item 173 (@code{LPX_I_NOFEAS}) | |
358 ## No integer feasible solution. | |
359 ## @end table | |
360 ## @noindent | |
361 ## If an error occurs, @var{status} will contain one of the following | |
362 ## codes: | |
5232 | 363 ## |
5244 | 364 ## @table @asis |
365 ## @item 204 (@code{LPX_E_FAULT}) | |
366 ## Unable to start the search. | |
367 ## @item 205 (@code{LPX_E_OBJLL}) | |
368 ## Objective function lower limit reached. | |
369 ## @item 206 (@code{LPX_E_OBJUL}) | |
370 ## Objective function upper limit reached. | |
371 ## @item 207 (@code{LPX_E_ITLIM}) | |
372 ## Iterations limit exhausted. | |
373 ## @item 208 (@code{LPX_E_TMLIM}) | |
374 ## Time limit exhausted. | |
375 ## @item 209 (@code{LPX_E_NOFEAS}) | |
376 ## No feasible solution. | |
377 ## @item 210 (@code{LPX_E_INSTAB}) | |
378 ## Numerical instability. | |
379 ## @item 211 (@code{LPX_E_SING}) | |
380 ## Problems with basis matrix. | |
381 ## @item 212 (@code{LPX_E_NOCONV}) | |
382 ## No convergence (interior). | |
383 ## @item 213 (@code{LPX_E_NOPFS}) | |
384 ## No primal feasible solution (LP presolver). | |
385 ## @item 214 (@code{LPX_E_NODFS}) | |
386 ## No dual feasible solution (LP presolver). | |
387 ## @end table | |
388 ## @item extra | |
389 ## A data structure containing the following fields: | |
390 ## @table @code | |
391 ## @item lambda | |
392 ## Dual variables. | |
393 ## @item redcosts | |
394 ## Reduced Costs. | |
395 ## @item time | |
396 ## Time (in seconds) used for solving LP/MIP problem. | |
397 ## @item mem | |
6333 | 398 ## Memory (in bytes) used for solving LP/MIP problem (this is not |
399 ## available if the version of GLPK is 4.15 or later). | |
5244 | 400 ## @end table |
401 ## @end table | |
402 ## | |
403 ## Example: | |
404 ## | |
405 ## @example | |
406 ## @group | |
407 ## c = [10, 6, 4]'; | |
408 ## a = [ 1, 1, 1; | |
409 ## 10, 4, 5; | |
410 ## 2, 2, 6]; | |
411 ## b = [100, 600, 300]'; | |
412 ## lb = [0, 0, 0]'; | |
413 ## ub = []; | |
414 ## ctype = "UUU"; | |
415 ## vartype = "CCC"; | |
416 ## s = -1; | |
417 ## | |
418 ## param.msglev = 1; | |
419 ## param.itlim = 100; | |
420 ## | |
421 ## [xmin, fmin, status, extra] = ... | |
422 ## glpk (c, a, b, lb, ub, ctype, vartype, s, param); | |
423 ## @end group | |
424 ## @end example | |
425 ## @end deftypefn | |
5232 | 426 |
5233 | 427 ## Author: Nicolo' Giorgetti <giorgetti@dii.unisi.it> |
428 ## Adapted-by: jwe | |
5232 | 429 |
5237 | 430 function [xopt, fmin, status, extra] = glpk (c, a, b, lb, ub, ctype, vartype, sense, param) |
5232 | 431 |
5233 | 432 ## If there is no input output the version and syntax |
5237 | 433 if (nargin < 3 || nargin > 9) |
6046 | 434 print_usage (); |
5233 | 435 return; |
436 endif | |
5232 | 437 |
5233 | 438 if (all (size (c) > 1) || iscomplex (c) || ischar (c)) |
439 error ("C must be a real vector"); | |
440 return; | |
441 endif | |
442 nx = length (c); | |
443 ## Force column vector. | |
444 c = c(:); | |
5232 | 445 |
5237 | 446 ## 2) Matrix constraint |
5232 | 447 |
5233 | 448 if (isempty (a)) |
449 error ("A cannot be an empty matrix"); | |
450 return; | |
451 endif | |
452 [nc, nxa] = size(a); | |
453 if (! isreal (a) || nxa != nx) | |
454 error ("A must be a real valued %d by %d matrix", nc, nx); | |
455 return; | |
456 endif | |
457 | |
5237 | 458 ## 3) RHS |
5232 | 459 |
5233 | 460 if (isempty (b)) |
461 error ("B cannot be an empty vector"); | |
462 return; | |
463 endif | |
464 if (! isreal (b) || length (b) != nc) | |
465 error ("B must be a real valued %d by 1 vector", nc); | |
466 return; | |
467 endif | |
468 | |
5237 | 469 ## 4) Vector with the lower bound of each variable |
5232 | 470 |
5237 | 471 if (nargin > 3) |
5233 | 472 if (isempty (lb)) |
8036
854683691d7a
fix invalid memory read in glpk
Jaroslav Hajek <highegg@gmail.com>
parents:
7151
diff
changeset
|
473 lb = zeros (nx, 1); |
5233 | 474 elseif (! isreal (lb) || all (size (lb) > 1) || length (lb) != nx) |
475 error ("LB must be a real valued %d by 1 column vector", nx); | |
476 return; | |
477 endif | |
478 else | |
5237 | 479 lb = zeros (nx, 1); |
7151 | 480 endif |
5233 | 481 |
5237 | 482 ## 5) Vector with the upper bound of each variable |
5232 | 483 |
5237 | 484 if (nargin > 4) |
5233 | 485 if (isempty (ub)) |
486 ub = repmat (Inf, nx, 1); | |
487 elseif (! isreal (ub) || all (size (ub) > 1) || length (ub) != nx) | |
488 error ("UB must be a real valued %d by 1 column vector", nx); | |
489 return; | |
490 endif | |
491 else | |
492 ub = repmat (Inf, nx, 1); | |
7151 | 493 endif |
5232 | 494 |
5237 | 495 ## 6) Sense of each constraint |
5232 | 496 |
5237 | 497 if (nargin > 5) |
498 if (isempty (ctype)) | |
499 ctype = repmat ("S", nc, 1); | |
500 elseif (! ischar (ctype) || all (size (ctype) > 1) || length (ctype) != nc) | |
5244 | 501 error ("CTYPE must be a char valued vector of length %d", nc); |
5237 | 502 return; |
503 elseif (! all (ctype == "F" | ctype == "U" | ctype == "S" | |
504 | ctype == "L" | ctype == "D")) | |
505 error ("CTYPE must contain only F, U, S, L, or D"); | |
506 return; | |
507 endif | |
508 else | |
509 ctype = repmat ("S", nc, 1); | |
7151 | 510 endif |
5237 | 511 |
512 ## 7) Vector with the type of variables | |
513 | |
514 if (nargin > 6) | |
5289 | 515 if (isempty (vartype)) |
5233 | 516 vartype = repmat ("C", nx, 1); |
517 elseif (! ischar (vartype) || all (size (vartype) > 1) | |
518 || length (vartype) != nx) | |
5244 | 519 error ("VARTYPE must be a char valued vector of length %d", nx); |
5233 | 520 return; |
521 elseif (! all (vartype == "C" | vartype == "I")) | |
522 error ("VARTYPE must contain only C or I"); | |
523 return; | |
524 endif | |
525 else | |
526 ## As default we consider continuous vars | |
527 vartype = repmat ("C", nx, 1); | |
528 endif | |
5232 | 529 |
5289 | 530 ## 8) Sense of optimization |
531 | |
532 if (nargin > 7) | |
533 if (isempty (sense)) | |
534 sense = 1; | |
535 elseif (ischar (sense) || all (size (sense) > 1) || ! isreal (sense)) | |
536 error ("SENSE must be an integer value"); | |
537 elseif (sense >= 0) | |
538 sense = 1; | |
539 else | |
540 sense = -1; | |
541 endif | |
542 else | |
543 sense = 1; | |
544 endif | |
545 | |
546 ## 9) Parameters vector | |
5233 | 547 |
548 if (nargin > 8) | |
549 if (! isstruct (param)) | |
550 error ("PARAM must be a structure"); | |
551 return; | |
552 endif | |
553 else | |
554 param = struct (); | |
555 endif | |
5232 | 556 |
5233 | 557 [xopt, fmin, status, extra] = ... |
5237 | 558 __glpk__ (c, a, b, lb, ub, ctype, vartype, sense, param); |
5232 | 559 |
560 endfunction |