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