octave_packages/m/optimization/glpk.m

   1 ## Copyright (C) 2005-2012 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 3 of the License, or (at
   8 ## your option) 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, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  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 @sc{glpk} library.  Given three
  22 ## arguments, @code{glpk} solves the following standard LP:
  23 ## @tex
  24 ## $$
  25 ##   \min_x C^T x
  26 ## $$
  27 ## @end tex
  28 ## @ifnottex
  29 ##
  30 ## @example
  31 ## min C'*x
  32 ## @end example
  33 ##
  34 ## @end ifnottex
  35 ## subject to
  36 ## @tex
  37 ## $$
  38 ##   Ax = b \qquad x \geq 0
  39 ## $$
  40 ## @end tex
  41 ## @ifnottex
  42 ##
  43 ## @example
  44 ## @group
  45 ## A*x  = b
  46 ##   x >= 0
  47 ## @end group
  48 ## @end example
  49 ##
  50 ## @end ifnottex
  51 ## but may also solve problems of the form
  52 ## @tex
  53 ## $$
  54 ##   [ \min_x | \max_x ] C^T x
  55 ## $$
  56 ## @end tex
  57 ## @ifnottex
  58 ##
  59 ## @example
  60 ## [ min | max ] C'*x
  61 ## @end example
  62 ##
  63 ## @end ifnottex
  64 ## subject to
  65 ## @tex
  66 ## $$
  67 ##  Ax [ = | \leq | \geq ] b \qquad LB \leq x \leq UB
  68 ## $$
  69 ## @end tex
  70 ## @ifnottex
  71 ##
  72 ## @example
  73 ## @group
  74 ## A*x [ "=" | "<=" | ">=" ] b
  75 ##   x >= LB
  76 ##   x <= UB
  77 ## @end group
  78 ## @end example
  79 ##
  80 ## @end ifnottex
  81 ##
  82 ## Input arguments:
  83 ##
  84 ## @table @var
  85 ## @item c
  86 ## A column array containing the objective function coefficients.
  87 ##
  88 ## @item A
  89 ## A matrix containing the constraints coefficients.
  90 ##
  91 ## @item b
  92 ## A column array containing the right-hand side value for each constraint
  93 ## in the constraint matrix.
  94 ##
  95 ## @item lb
  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
  98 ## zero.
  99 ##
 100 ## @item ub
 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
 103 ## infinite.
 104 ##
 105 ## @item ctype
 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
 108 ## following values
 109 ## @table @asis
 110 ## @item "F"
 111 ## A free (unbounded) constraint (the constraint is ignored).
 112 ##
 113 ## @item "U"
 114 ## An inequality constraint with an upper bound (@code{A(i,:)*x <= b(i)}).
 115 ##
 116 ## @item "S"
 117 ## An equality constraint (@code{A(i,:)*x = b(i)}).
 118 ##
 119 ## @item "L"
 120 ## An inequality with a lower bound (@code{A(i,:)*x >= b(i)}).
 121 ##
 122 ## @item "D"
 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)}).
 125 ## @end table
 126 ##
 127 ## @item vartype
 128 ## A column array containing the types of the variables.
 129 ## @table @asis
 130 ## @item "C"
 131 ## A continuous variable.
 132 ##
 133 ## @item "I"
 134 ## An integer variable.
 135 ## @end table
 136 ##
 137 ## @item sense
 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.
 140 ##
 141 ## @item param
 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
 145 ## from the default.
 146 ##
 147 ## Integer parameters:
 148 ##
 149 ## @table @code
 150 ## @item msglev (@w{@code{LPX_K_MSGLEV}}, default: 1)
 151 ## Level of messages output by solver routines:
 152 ## @table @asis
 153 ## @item 0
 154 ## No output.
 155 ##
 156 ## @item 1
 157 ## Error messages only.
 158 ##
 159 ## @item 2
 160 ## Normal output.
 161 ##
 162 ## @item 3
 163 ## Full output (includes informational messages).
 164 ## @end table
 165 ##
 166 ## @item scale (@w{@code{LPX_K_SCALE}}, default: 1)
 167 ## Scaling option:
 168 ## @table @asis
 169 ## @item 0
 170 ## No scaling.
 171 ##
 172 ## @item 1
 173 ## Equilibration scaling.
 174 ##
 175 ## @item 2
 176 ## Geometric mean scaling, then equilibration scaling.
 177 ## @end table
 178 ##
 179 ## @item dual    (@w{@code{LPX_K_DUAL}}, default: 0)
 180 ## Dual simplex option:
 181 ## @table @asis
 182 ## @item 0
 183 ## Do not use the dual simplex.
 184 ##
 185 ## @item 1
 186 ## If initial basic solution is dual feasible, use the dual simplex.
 187 ## @end table
 188 ##
 189 ## @item price   (@w{@code{LPX_K_PRICE}}, default: 1)
 190 ## Pricing option (for both primal and dual simplex):
 191 ## @table @asis
 192 ## @item 0
 193 ## Textbook pricing.
 194 ##
 195 ## @item 1
 196 ## Steepest edge pricing.
 197 ## @end table
 198 ##
 199 ## @item round   (@w{@code{LPX_K_ROUND}}, default: 0)
 200 ## Solution rounding option:
 201 ## @table @asis
 202 ## @item 0
 203 ## Report all primal and dual values "as is".
 204 ##
 205 ## @item 1
 206 ## Replace tiny primal and dual values by exact zero.
 207 ## @end table
 208 ##
 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.
 214 ##
 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
 218 ## standard output.
 219 ##
 220 ## @item branch (@w{@code{LPX_K_BRANCH}}, default: 2)
 221 ## Branching heuristic option (for MIP only):
 222 ## @table @asis
 223 ## @item 0
 224 ## Branch on the first variable.
 225 ##
 226 ## @item 1
 227 ## Branch on the last variable.
 228 ##
 229 ## @item 2
 230 ## Branch using a heuristic by Driebeck and Tomlin.
 231 ## @end table
 232 ##
 233 ## @item btrack (@w{@code{LPX_K_BTRACK}}, default: 2)
 234 ## Backtracking heuristic option (for MIP only):
 235 ## @table @asis
 236 ## @item 0
 237 ## Depth first search.
 238 ##
 239 ## @item 1
 240 ## Breadth first search.
 241 ##
 242 ## @item 2
 243 ## Backtrack using the best projection heuristic.
 244 ## @end table
 245 ##
 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.
 249 ##
 250 ## @item lpsolver (default: 1)
 251 ## Select which solver to use.  If the problem is a MIP problem this flag
 252 ## will be ignored.
 253 ## @table @asis
 254 ## @item 1
 255 ## Revised simplex method.
 256 ##
 257 ## @item 2
 258 ## Interior point method.
 259 ## @end table
 260 ##
 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.
 265 ## @end table
 266 ##
 267 ## Real parameters:
 268 ##
 269 ## @table @code
 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}}.
 279 ##
 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.
 284 ##
 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.
 289 ##
 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.
 294 ##
 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.
 299 ##
 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.
 304 ##
 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.
 310 ##
 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.
 315 ##
 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.
 320 ##
 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.
 326 ## @end table
 327 ## @end table
 328 ##
 329 ## Output values:
 330 ##
 331 ## @table @var
 332 ## @item xopt
 333 ## The optimizer (the value of the decision variables at the optimum).
 334 ##
 335 ## @item fopt
 336 ## The optimum value of the objective function.
 337 ##
 338 ## @item status
 339 ## Status of the optimization.
 340 ##
 341 ## Simplex Method:
 342 ## @table @asis
 343 ## @item 180 (@w{@code{LPX_OPT}})
 344 ## Solution is optimal.
 345 ##
 346 ## @item 181 (@w{@code{LPX_FEAS}})
 347 ## Solution is feasible.
 348 ##
 349 ## @item 182 (@w{@code{LPX_INFEAS}})
 350 ## Solution is infeasible.
 351 ##
 352 ## @item 183 (@w{@code{LPX_NOFEAS}})
 353 ## Problem has no feasible solution.
 354 ##
 355 ## @item 184 (@w{@code{LPX_UNBND}})
 356 ## Problem has no unbounded solution.
 357 ##
 358 ## @item 185 (@w{@code{LPX_UNDEF}})
 359 ## Solution status is undefined.
 360 ## @end table
 361 ## Interior Point Method:
 362 ## @table @asis
 363 ## @item 150 (@w{@code{LPX_T_UNDEF}})
 364 ## The interior point method is undefined.
 365 ##
 366 ## @item 151 (@w{@code{LPX_T_OPT}})
 367 ## The interior point method is optimal.
 368 ## @end table
 369 ## Mixed Integer Method:
 370 ## @table @asis
 371 ## @item 170 (@w{@code{LPX_I_UNDEF}})
 372 ## The status is undefined.
 373 ##
 374 ## @item 171 (@w{@code{LPX_I_OPT}})
 375 ## The solution is integer optimal.
 376 ##
 377 ## @item 172 (@w{@code{LPX_I_FEAS}})
 378 ## Solution integer feasible but its optimality has not been proven
 379 ##
 380 ## @item 173 (@w{@code{LPX_I_NOFEAS}})
 381 ## No integer feasible solution.
 382 ## @end table
 383 ## @noindent
 384 ## If an error occurs, @var{status} will contain one of the following
 385 ## codes:
 386 ##
 387 ## @table @asis
 388 ## @item 204 (@w{@code{LPX_E_FAULT}})
 389 ## Unable to start the search.
 390 ##
 391 ## @item 205 (@w{@code{LPX_E_OBJLL}})
 392 ## Objective function lower limit reached.
 393 ##
 394 ## @item 206 (@w{@code{LPX_E_OBJUL}})
 395 ## Objective function upper limit reached.
 396 ##
 397 ## @item 207 (@w{@code{LPX_E_ITLIM}})
 398 ## Iterations limit exhausted.
 399 ##
 400 ## @item 208 (@w{@code{LPX_E_TMLIM}})
 401 ## Time limit exhausted.
 402 ##
 403 ## @item 209 (@w{@code{LPX_E_NOFEAS}})
 404 ## No feasible solution.
 405 ##
 406 ## @item 210 (@w{@code{LPX_E_INSTAB}})
 407 ## Numerical instability.
 408 ##
 409 ## @item 211 (@w{@code{LPX_E_SING}})
 410 ## Problems with basis matrix.
 411 ##
 412 ## @item 212 (@w{@code{LPX_E_NOCONV}})
 413 ## No convergence (interior).
 414 ##
 415 ## @item 213 (@w{@code{LPX_E_NOPFS}})
 416 ## No primal feasible solution (LP presolver).
 417 ##
 418 ## @item 214 (@w{@code{LPX_E_NODFS}})
 419 ## No dual feasible solution (LP presolver).
 420 ## @end table
 421 ##
 422 ## @item extra
 423 ## A data structure containing the following fields:
 424 ## @table @code
 425 ## @item lambda
 426 ## Dual variables.
 427 ##
 428 ## @item redcosts
 429 ## Reduced Costs.
 430 ##
 431 ## @item time
 432 ## Time (in seconds) used for solving LP/MIP problem.
 433 ##
 434 ## @item mem
 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).
 437 ## @end table
 438 ## @end table
 439 ##
 440 ## Example:
 441 ##
 442 ## @example
 443 ## @group
 444 ## c = [10, 6, 4]';
 445 ## A = [ 1, 1, 1;
 446 ##      10, 4, 5;
 447 ##       2, 2, 6];
 448 ## b = [100, 600, 300]';
 449 ## lb = [0, 0, 0]';
 450 ## ub = [];
 451 ## ctype = "UUU";
 452 ## vartype = "CCC";
 453 ## s = -1;
 454 ##
 455 ## param.msglev = 1;
 456 ## param.itlim = 100;
 457 ##
 458 ## [xmin, fmin, status, extra] = ...
 459 ##    glpk (c, A, b, lb, ub, ctype, vartype, s, param);
 460 ## @end group
 461 ## @end example
 462 ## @end deftypefn
 463
 464 ## Author: Nicolo' Giorgetti <giorgetti@dii.unisi.it>
 465 ## Adapted-by: jwe
 466
 467 function [xopt, fmin, status, extra] = glpk (c, A, b, lb, ub, ctype, vartype, sense, param)
 468
 469   ## If there is no input output the version and syntax
 470   if (nargin < 3 || nargin > 9)
 471     print_usage ();
 472     return;
 473   endif
 474
 475   if (all (size (c) > 1) || iscomplex (c) || ischar (c))
 476     error ("glpk:C must be a real vector");
 477     return;
 478   endif
 479   nx = length (c);
 480   ## Force column vector.
 481   c = c(:);
 482
 483   ## 2) Matrix constraint
 484
 485   if (isempty (A))
 486     error ("glpk: A cannot be an empty matrix");
 487     return;
 488   endif
 489   [nc, nxa] = size(A);
 490   if (! isreal (A) || nxa != nx)
 491     error ("glpk: A must be a real valued %d by %d matrix", nc, nx);
 492     return;
 493   endif
 494
 495   ## 3) RHS
 496
 497   if (isempty (b))
 498     error ("glpk: B cannot be an empty vector");
 499     return;
 500   endif
 501   if (! isreal (b) || length (b) != nc)
 502     error ("glpk: B must be a real valued %d by 1 vector", nc);
 503     return;
 504   endif
 505
 506   ## 4) Vector with the lower bound of each variable
 507
 508   if (nargin > 3)
 509     if (isempty (lb))
 510       lb = zeros (nx, 1);
 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);
 513       return;
 514     endif
 515   else
 516     lb = zeros (nx, 1);
 517   endif
 518
 519   ## 5) Vector with the upper bound of each variable
 520
 521   if (nargin > 4)
 522     if (isempty (ub))
 523       ub = Inf (nx, 1);
 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);
 526       return;
 527     endif
 528   else
 529     ub = Inf (nx, 1);
 530   endif
 531
 532   ## 6) Sense of each constraint
 533
 534   if (nargin > 5)
 535     if (isempty (ctype))
 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);
 539       return;
 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");
 543       return;
 544     endif
 545   else
 546     ctype = repmat ("S", nc, 1);
 547   endif
 548
 549   ## 7) Vector with the type of variables
 550
 551   if (nargin > 6)
 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);
 557       return;
 558     elseif (! all (vartype == "C" | vartype == "I"))
 559       error ("glpk: VARTYPE must contain only C or I");
 560       return;
 561     endif
 562   else
 563     ## As default we consider continuous vars
 564     vartype = repmat ("C", nx, 1);
 565   endif
 566
 567   ## 8) Sense of optimization
 568
 569   if (nargin > 7)
 570     if (isempty (sense))
 571       sense = 1;
 572     elseif (ischar (sense) || all (size (sense) > 1) || ! isreal (sense))
 573       error ("glpk: SENSE must be an integer value");
 574     elseif (sense >= 0)
 575       sense = 1;
 576     else
 577       sense = -1;
 578     endif
 579   else
 580     sense = 1;
 581   endif
 582
 583   ## 9) Parameters vector
 584
 585   if (nargin > 8)
 586     if (! isstruct (param))
 587       error ("glpk: PARAM must be a structure");
 588       return;
 589     endif
 590   else
 591     param = struct ();
 592   endif
 593
 594   [xopt, fmin, status, extra] = ...
 595     __glpk__ (c, A, b, lb, ub, ctype, vartype, sense, param);
 596
 597 endfunction