1 ## Copyright (C) 2005-2012 Nicolo' Giorgetti
3 ## This file is part of Octave.
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 3 of the License, or (at
8 ## your option) any later version.
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.
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
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 @sc{glpk} library. Given three
22 ## arguments, @code{glpk} solves the following standard LP:
38 ## Ax = b \qquad x \geq 0
51 ## but may also solve problems of the form
54 ## [ \min_x | \max_x ] C^T x
67 ## Ax [ = | \leq | \geq ] b \qquad LB \leq x \leq UB
74 ## A*x [ "=" | "<=" | ">=" ] b
86 ## A column array containing the objective function coefficients.
89 ## A matrix containing the constraints coefficients.
92 ## A column array containing the right-hand side value for each constraint
93 ## in the constraint matrix.
96 ## An array containing the lower bound on each of the variables. If
97 ## @var{lb} is not supplied, the default lower bound for the variables is
101 ## An array containing the upper bound on each of the variables. If
102 ## @var{ub} is not supplied, the default upper bound is assumed to be
106 ## An array of characters containing the sense of each constraint in the
107 ## constraint matrix. Each element of the array may be one of the
111 ## A free (unbounded) constraint (the constraint is ignored).
114 ## An inequality constraint with an upper bound (@code{A(i,:)*x <= b(i)}).
117 ## An equality constraint (@code{A(i,:)*x = b(i)}).
120 ## An inequality with a lower bound (@code{A(i,:)*x >= b(i)}).
123 ## An inequality constraint with both upper and lower bounds
124 ## (@code{A(i,:)*x >= -b(i)} @emph{and} (@code{A(i,:)*x <= b(i)}).
128 ## A column array containing the types of the variables.
131 ## A continuous variable.
134 ## An integer variable.
138 ## If @var{sense} is 1, the problem is a minimization. If @var{sense} is
139 ## -1, the problem is a maximization. The default value is 1.
142 ## A structure containing the following parameters used to define the
143 ## behavior of solver. Missing elements in the structure take on default
144 ## values, so you only need to set the elements that you wish to change
147 ## Integer parameters:
150 ## @item msglev (@w{@code{LPX_K_MSGLEV}}, default: 1)
151 ## Level of messages output by solver routines:
157 ## Error messages only.
163 ## Full output (includes informational messages).
166 ## @item scale (@w{@code{LPX_K_SCALE}}, default: 1)
173 ## Equilibration scaling.
176 ## Geometric mean scaling, then equilibration scaling.
179 ## @item dual (@w{@code{LPX_K_DUAL}}, default: 0)
180 ## Dual simplex option:
183 ## Do not use the dual simplex.
186 ## If initial basic solution is dual feasible, use the dual simplex.
189 ## @item price (@w{@code{LPX_K_PRICE}}, default: 1)
190 ## Pricing option (for both primal and dual simplex):
196 ## Steepest edge pricing.
199 ## @item round (@w{@code{LPX_K_ROUND}}, default: 0)
200 ## Solution rounding option:
203 ## Report all primal and dual values "as is".
206 ## Replace tiny primal and dual values by exact zero.
209 ## @item itlim (@w{@code{LPX_K_ITLIM}}, default: -1)
210 ## Simplex iterations limit. If this value is positive, it is decreased by
211 ## one each time when one simplex iteration has been performed, and
212 ## reaching zero value signals the solver to stop the search. Negative
213 ## value means no iterations limit.
215 ## @item itcnt (@w{@code{LPX_K_OUTFRQ}}, default: 200)
216 ## Output frequency, in iterations. This parameter specifies how
217 ## frequently the solver sends information about the solution to the
220 ## @item branch (@w{@code{LPX_K_BRANCH}}, default: 2)
221 ## Branching heuristic option (for MIP only):
224 ## Branch on the first variable.
227 ## Branch on the last variable.
230 ## Branch using a heuristic by Driebeck and Tomlin.
233 ## @item btrack (@w{@code{LPX_K_BTRACK}}, default: 2)
234 ## Backtracking heuristic option (for MIP only):
237 ## Depth first search.
240 ## Breadth first search.
243 ## Backtrack using the best projection heuristic.
246 ## @item presol (@w{@code{LPX_K_PRESOL}}, default: 1)
247 ## If this flag is set, the routine lpx_simplex solves the problem using
248 ## the built-in LP presolver. Otherwise the LP presolver is not used.
250 ## @item lpsolver (default: 1)
251 ## Select which solver to use. If the problem is a MIP problem this flag
255 ## Revised simplex method.
258 ## Interior point method.
261 ## @item save (default: 0)
262 ## If this parameter is nonzero, save a copy of the problem in
263 ## CPLEX LP format to the file @file{"outpb.lp"}. There is currently no
264 ## way to change the name of the output file.
270 ## @item relax (@w{@code{LPX_K_RELAX}}, default: 0.07)
271 ## Relaxation parameter used in the ratio test. If it is zero, the textbook
272 ## ratio test is used. If it is non-zero (should be positive), Harris'
273 ## two-pass ratio test is used. In the latter case on the first pass of the
274 ## ratio test basic variables (in the case of primal simplex) or reduced
275 ## costs of non-basic variables (in the case of dual simplex) are allowed
276 ## to slightly violate their bounds, but not more than
277 ## @code{relax*tolbnd} or @code{relax*toldj (thus, @code{relax} is a
278 ## percentage of @code{tolbnd} or @code{toldj}}.
280 ## @item tolbnd (@w{@code{LPX_K_TOLBND}}, default: 10e-7)
281 ## Relative tolerance used to check if the current basic solution is primal
282 ## feasible. It is not recommended that you change this parameter unless you
283 ## have a detailed understanding of its purpose.
285 ## @item toldj (@w{@code{LPX_K_TOLDJ}}, default: 10e-7)
286 ## Absolute tolerance used to check if the current basic solution is dual
287 ## feasible. It is not recommended that you change this parameter unless you
288 ## have a detailed understanding of its purpose.
290 ## @item tolpiv (@w{@code{LPX_K_TOLPIV}}, default: 10e-9)
291 ## Relative tolerance used to choose eligible pivotal elements of the
292 ## simplex table. It is not recommended that you change this parameter unless
293 ## you have a detailed understanding of its purpose.
295 ## @item objll (@w{@code{LPX_K_OBJLL}}, default: -DBL_MAX)
296 ## Lower limit of the objective function. If on the phase II the objective
297 ## function reaches this limit and continues decreasing, the solver stops
298 ## the search. This parameter is used in the dual simplex method only.
300 ## @item objul (@w{@code{LPX_K_OBJUL}}, default: +DBL_MAX)
301 ## Upper limit of the objective function. If on the phase II the objective
302 ## function reaches this limit and continues increasing, the solver stops
303 ## the search. This parameter is used in the dual simplex only.
305 ## @item tmlim (@w{@code{LPX_K_TMLIM}}, default: -1.0)
306 ## Searching time limit, in seconds. If this value is positive, it is
307 ## decreased each time when one simplex iteration has been performed by the
308 ## amount of time spent for the iteration, and reaching zero value signals
309 ## the solver to stop the search. Negative value means no time limit.
311 ## @item outdly (@w{@code{LPX_K_OUTDLY}}, default: 0.0)
312 ## Output delay, in seconds. This parameter specifies how long the solver
313 ## should delay sending information about the solution to the standard
314 ## output. Non-positive value means no delay.
316 ## @item tolint (@w{@code{LPX_K_TOLINT}}, default: 10e-5)
317 ## Relative tolerance used to check if the current basic solution is integer
318 ## feasible. It is not recommended that you change this parameter unless
319 ## you have a detailed understanding of its purpose.
321 ## @item tolobj (@w{@code{LPX_K_TOLOBJ}}, default: 10e-7)
322 ## Relative tolerance used to check if the value of the objective function
323 ## is not better than in the best known integer feasible solution. It is
324 ## not recommended that you change this parameter unless you have a
325 ## detailed understanding of its purpose.
333 ## The optimizer (the value of the decision variables at the optimum).
336 ## The optimum value of the objective function.
339 ## Status of the optimization.
343 ## @item 180 (@w{@code{LPX_OPT}})
344 ## Solution is optimal.
346 ## @item 181 (@w{@code{LPX_FEAS}})
347 ## Solution is feasible.
349 ## @item 182 (@w{@code{LPX_INFEAS}})
350 ## Solution is infeasible.
352 ## @item 183 (@w{@code{LPX_NOFEAS}})
353 ## Problem has no feasible solution.
355 ## @item 184 (@w{@code{LPX_UNBND}})
356 ## Problem has no unbounded solution.
358 ## @item 185 (@w{@code{LPX_UNDEF}})
359 ## Solution status is undefined.
361 ## Interior Point Method:
363 ## @item 150 (@w{@code{LPX_T_UNDEF}})
364 ## The interior point method is undefined.
366 ## @item 151 (@w{@code{LPX_T_OPT}})
367 ## The interior point method is optimal.
369 ## Mixed Integer Method:
371 ## @item 170 (@w{@code{LPX_I_UNDEF}})
372 ## The status is undefined.
374 ## @item 171 (@w{@code{LPX_I_OPT}})
375 ## The solution is integer optimal.
377 ## @item 172 (@w{@code{LPX_I_FEAS}})
378 ## Solution integer feasible but its optimality has not been proven
380 ## @item 173 (@w{@code{LPX_I_NOFEAS}})
381 ## No integer feasible solution.
384 ## If an error occurs, @var{status} will contain one of the following
388 ## @item 204 (@w{@code{LPX_E_FAULT}})
389 ## Unable to start the search.
391 ## @item 205 (@w{@code{LPX_E_OBJLL}})
392 ## Objective function lower limit reached.
394 ## @item 206 (@w{@code{LPX_E_OBJUL}})
395 ## Objective function upper limit reached.
397 ## @item 207 (@w{@code{LPX_E_ITLIM}})
398 ## Iterations limit exhausted.
400 ## @item 208 (@w{@code{LPX_E_TMLIM}})
401 ## Time limit exhausted.
403 ## @item 209 (@w{@code{LPX_E_NOFEAS}})
404 ## No feasible solution.
406 ## @item 210 (@w{@code{LPX_E_INSTAB}})
407 ## Numerical instability.
409 ## @item 211 (@w{@code{LPX_E_SING}})
410 ## Problems with basis matrix.
412 ## @item 212 (@w{@code{LPX_E_NOCONV}})
413 ## No convergence (interior).
415 ## @item 213 (@w{@code{LPX_E_NOPFS}})
416 ## No primal feasible solution (LP presolver).
418 ## @item 214 (@w{@code{LPX_E_NODFS}})
419 ## No dual feasible solution (LP presolver).
423 ## A data structure containing the following fields:
432 ## Time (in seconds) used for solving LP/MIP problem.
435 ## Memory (in bytes) used for solving LP/MIP problem (this is not
436 ## available if the version of @sc{glpk} is 4.15 or later).
448 ## b = [100, 600, 300]';
456 ## param.itlim = 100;
458 ## [xmin, fmin, status, extra] = ...
459 ## glpk (c, A, b, lb, ub, ctype, vartype, s, param);
464 ## Author: Nicolo' Giorgetti <giorgetti@dii.unisi.it>
467 function [xopt, fmin, status, extra] = glpk (c, A, b, lb, ub, ctype, vartype, sense, param)
469 ## If there is no input output the version and syntax
470 if (nargin < 3 || nargin > 9)
475 if (all (size (c) > 1) || iscomplex (c) || ischar (c))
476 error ("glpk:C must be a real vector");
480 ## Force column vector.
483 ## 2) Matrix constraint
486 error ("glpk: A cannot be an empty matrix");
490 if (! isreal (A) || nxa != nx)
491 error ("glpk: A must be a real valued %d by %d matrix", nc, nx);
498 error ("glpk: B cannot be an empty vector");
501 if (! isreal (b) || length (b) != nc)
502 error ("glpk: B must be a real valued %d by 1 vector", nc);
506 ## 4) Vector with the lower bound of each variable
511 elseif (! isreal (lb) || all (size (lb) > 1) || length (lb) != nx)
512 error ("glpk: LB must be a real valued %d by 1 column vector", nx);
519 ## 5) Vector with the upper bound of each variable
524 elseif (! isreal (ub) || all (size (ub) > 1) || length (ub) != nx)
525 error ("glpk: UB must be a real valued %d by 1 column vector", nx);
532 ## 6) Sense of each constraint
536 ctype = repmat ("S", nc, 1);
537 elseif (! ischar (ctype) || all (size (ctype) > 1) || length (ctype) != nc)
538 error ("glpk: CTYPE must be a char valued vector of length %d", nc);
540 elseif (! all (ctype == "F" | ctype == "U" | ctype == "S"
541 | ctype == "L" | ctype == "D"))
542 error ("glpk: CTYPE must contain only F, U, S, L, or D");
546 ctype = repmat ("S", nc, 1);
549 ## 7) Vector with the type of variables
552 if (isempty (vartype))
553 vartype = repmat ("C", nx, 1);
554 elseif (! ischar (vartype) || all (size (vartype) > 1)
555 || length (vartype) != nx)
556 error ("glpk: VARTYPE must be a char valued vector of length %d", nx);
558 elseif (! all (vartype == "C" | vartype == "I"))
559 error ("glpk: VARTYPE must contain only C or I");
563 ## As default we consider continuous vars
564 vartype = repmat ("C", nx, 1);
567 ## 8) Sense of optimization
572 elseif (ischar (sense) || all (size (sense) > 1) || ! isreal (sense))
573 error ("glpk: SENSE must be an integer value");
583 ## 9) Parameters vector
586 if (! isstruct (param))
587 error ("glpk: PARAM must be a structure");
594 [xopt, fmin, status, extra] = ...
595 __glpk__ (c, A, b, lb, ub, ctype, vartype, sense, param);