1######################################################################## 2## 3## Copyright (C) 2005-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {[@var{xopt}, @var{fmin}, @var{errnum}, @var{extra}] =} glpk (@var{c}, @var{A}, @var{b}, @var{lb}, @var{ub}, @var{ctype}, @var{vartype}, @var{sense}, @var{param}) 28## Solve a linear program using the GNU @sc{glpk} library. 29## 30## Given three arguments, @code{glpk} solves the following standard LP: 31## @tex 32## $$ 33## \min_x C^T x 34## $$ 35## @end tex 36## @ifnottex 37## 38## @example 39## min C'*x 40## @end example 41## 42## @end ifnottex 43## subject to 44## @tex 45## $$ 46## Ax = b \qquad x \geq 0 47## $$ 48## @end tex 49## @ifnottex 50## 51## @example 52## @group 53## A*x = b 54## x >= 0 55## @end group 56## @end example 57## 58## @end ifnottex 59## but may also solve problems of the form 60## @tex 61## $$ 62## [ \min_x | \max_x ] C^T x 63## $$ 64## @end tex 65## @ifnottex 66## 67## @example 68## [ min | max ] C'*x 69## @end example 70## 71## @end ifnottex 72## subject to 73## @tex 74## $$ 75## Ax [ = | \leq | \geq ] b \qquad LB \leq x \leq UB 76## $$ 77## @end tex 78## @ifnottex 79## 80## @example 81## @group 82## A*x [ "=" | "<=" | ">=" ] b 83## x >= LB 84## x <= UB 85## @end group 86## @end example 87## 88## @end ifnottex 89## 90## Input arguments: 91## 92## @table @var 93## @item c 94## A column array containing the objective function coefficients. 95## 96## @item A 97## A matrix containing the constraints coefficients. 98## 99## @item b 100## A column array containing the right-hand side value for each constraint in 101## the constraint matrix. 102## 103## @item lb 104## An array containing the lower bound on each of the variables. If @var{lb} 105## is not supplied, the default lower bound for the variables is zero. 106## 107## @item ub 108## An array containing the upper bound on each of the variables. If @var{ub} 109## is not supplied, the default upper bound is assumed to be infinite. 110## 111## @item ctype 112## An array of characters containing the sense of each constraint in the 113## constraint matrix. Each element of the array may be one of the following 114## values 115## 116## @table @asis 117## @item @qcode{"F"} 118## A free (unbounded) constraint (the constraint is ignored). 119## 120## @item @qcode{"U"} 121## An inequality constraint with an upper bound (@code{A(i,:)*x <= b(i)}). 122## 123## @item @qcode{"S"} 124## An equality constraint (@code{A(i,:)*x = b(i)}). 125## 126## @item @qcode{"L"} 127## An inequality with a lower bound (@code{A(i,:)*x >= b(i)}). 128## 129## @item @qcode{"D"} 130## An inequality constraint with both upper and lower bounds 131## (@code{A(i,:)*x >= -b(i)}) @emph{and} (@code{A(i,:)*x <= b(i)}). 132## @end table 133## 134## @item vartype 135## A column array containing the types of the variables. 136## 137## @table @asis 138## @item @qcode{"C"} 139## A continuous variable. 140## 141## @item @qcode{"I"} 142## An integer variable. 143## @end table 144## 145## @item sense 146## If @var{sense} is 1, the problem is a minimization. If @var{sense} is -1, 147## the problem is a maximization. The default value is 1. 148## 149## @item param 150## A structure containing the following parameters used to define the 151## behavior of solver. Missing elements in the structure take on default 152## values, so you only need to set the elements that you wish to change from 153## the default. 154## 155## Integer parameters: 156## 157## @table @code 158## @item msglev (default: 1) 159## Level of messages output by solver routines: 160## 161## @table @asis 162## @item 0 (@w{@code{GLP_MSG_OFF}}) 163## No output. 164## 165## @item 1 (@w{@code{GLP_MSG_ERR}}) 166## Error and warning messages only. 167## 168## @item 2 (@w{@code{GLP_MSG_ON}}) 169## Normal output. 170## 171## @item 3 (@w{@code{GLP_MSG_ALL}}) 172## Full output (includes informational messages). 173## @end table 174## 175## @item scale (default: 16) 176## Scaling option. The values can be combined with the bitwise OR operator and 177## may be the following: 178## 179## @table @asis 180## @item 1 (@w{@code{GLP_SF_GM}}) 181## Geometric mean scaling. 182## 183## @item 16 (@w{@code{GLP_SF_EQ}}) 184## Equilibration scaling. 185## 186## @item 32 (@w{@code{GLP_SF_2N}}) 187## Round scale factors to power of two. 188## 189## @item 64 (@w{@code{GLP_SF_SKIP}}) 190## Skip if problem is well scaled. 191## @end table 192## 193## Alternatively, a value of 128 (@w{@env{GLP_SF_AUTO}}) may be also 194## specified, in which case the routine chooses the scaling options 195## automatically. 196## 197## @item dual (default: 1) 198## Simplex method option: 199## 200## @table @asis 201## @item 1 (@w{@code{GLP_PRIMAL}}) 202## Use two-phase primal simplex. 203## 204## @item 2 (@w{@code{GLP_DUALP}}) 205## Use two-phase dual simplex, and if it fails, switch to the primal simplex. 206## 207## @item 3 (@w{@code{GLP_DUAL}}) 208## Use two-phase dual simplex. 209## @end table 210## 211## @item price (default: 34) 212## Pricing option (for both primal and dual simplex): 213## 214## @table @asis 215## @item 17 (@w{@code{GLP_PT_STD}}) 216## Textbook pricing. 217## 218## @item 34 (@w{@code{GLP_PT_PSE}}) 219## Steepest edge pricing. 220## @end table 221## 222## @item itlim (default: intmax) 223## Simplex iterations limit. It is decreased by one each time when one simplex 224## iteration has been performed, and reaching zero value signals the solver to 225## stop the search. 226## 227## @item outfrq (default: 200) 228## Output frequency, in iterations. This parameter specifies how frequently 229## the solver sends information about the solution to the standard output. 230## 231## @item branch (default: 4) 232## Branching technique option (for MIP only): 233## 234## @table @asis 235## @item 1 (@w{@code{GLP_BR_FFV}}) 236## First fractional variable. 237## 238## @item 2 (@w{@code{GLP_BR_LFV}}) 239## Last fractional variable. 240## 241## @item 3 (@w{@code{GLP_BR_MFV}}) 242## Most fractional variable. 243## 244## @item 4 (@w{@code{GLP_BR_DTH}}) 245## Heuristic by @nospell{Driebeck and Tomlin}. 246## 247## @item 5 (@w{@code{GLP_BR_PCH}}) 248## Hybrid @nospell{pseudocost} heuristic. 249## @end table 250## 251## @item btrack (default: 4) 252## Backtracking technique option (for MIP only): 253## 254## @table @asis 255## @item 1 (@w{@code{GLP_BT_DFS}}) 256## Depth first search. 257## 258## @item 2 (@w{@code{GLP_BT_BFS}}) 259## Breadth first search. 260## 261## @item 3 (@w{@code{GLP_BT_BLB}}) 262## Best local bound. 263## 264## @item 4 (@w{@code{GLP_BT_BPH}}) 265## Best projection heuristic. 266## @end table 267## 268## @item presol (default: 1) 269## If this flag is set, the simplex solver uses the built-in LP presolver. 270## Otherwise the LP presolver is not used. 271## 272## @item lpsolver (default: 1) 273## Select which solver to use. If the problem is a MIP problem this flag 274## will be ignored. 275## 276## @table @asis 277## @item 1 278## Revised simplex method. 279## 280## @item 2 281## Interior point method. 282## @end table 283## 284## @item rtest (default: 34) 285## Ratio test technique: 286## 287## @table @asis 288## @item 17 (@w{@code{GLP_RT_STD}}) 289## Standard ("textbook"). 290## 291## @item 34 (@w{@code{GLP_RT_HAR}}) 292## Harris' two-pass ratio test. 293## @end table 294## 295## @item tmlim (default: intmax) 296## Searching time limit, in milliseconds. 297## 298## @item outdly (default: 0) 299## Output delay, in seconds. This parameter specifies how long the solver 300## should delay sending information about the solution to the standard output. 301## 302## @item save (default: 0) 303## If this parameter is nonzero, save a copy of the problem in @nospell{CPLEX} 304## LP format to the file @file{"outpb.lp"}. There is currently no way to 305## change the name of the output file. 306## @end table 307## 308## Real parameters: 309## 310## @table @code 311## @item tolbnd (default: 1e-7) 312## Relative tolerance used to check if the current basic solution is primal 313## feasible. It is not recommended that you change this parameter unless you 314## have a detailed understanding of its purpose. 315## 316## @item toldj (default: 1e-7) 317## Absolute tolerance used to check if the current basic solution is dual 318## feasible. It is not recommended that you change this parameter unless you 319## have a detailed understanding of its purpose. 320## 321## @item tolpiv (default: 1e-10) 322## Relative tolerance used to choose eligible pivotal elements of the simplex 323## table. It is not recommended that you change this parameter unless you have 324## a detailed understanding of its purpose. 325## 326## @item objll (default: -DBL_MAX) 327## Lower limit of the objective function. If the objective function reaches 328## this limit and continues decreasing, the solver stops the search. This 329## parameter is used in the dual simplex method only. 330## 331## @item objul (default: +DBL_MAX) 332## Upper limit of the objective function. If the objective function reaches 333## this limit and continues increasing, the solver stops the search. This 334## parameter is used in the dual simplex only. 335## 336## @item tolint (default: 1e-5) 337## Relative tolerance used to check if the current basic solution is integer 338## feasible. It is not recommended that you change this parameter unless you 339## have a detailed understanding of its purpose. 340## 341## @item tolobj (default: 1e-7) 342## Relative tolerance used to check if the value of the objective function is 343## not better than in the best known integer feasible solution. It is not 344## recommended that you change this parameter unless you have a detailed 345## understanding of its purpose. 346## @end table 347## @end table 348## 349## Output values: 350## 351## @table @var 352## @item xopt 353## The optimizer (the value of the decision variables at the optimum). 354## 355## @item fopt 356## The optimum value of the objective function. 357## 358## @item errnum 359## Error code. 360## 361## @table @asis 362## @item 0 363## No error. 364## 365## @item 1 (@w{@code{GLP_EBADB}}) 366## Invalid basis. 367## 368## @item 2 (@w{@code{GLP_ESING}}) 369## Singular matrix. 370## 371## @item 3 (@w{@code{GLP_ECOND}}) 372## Ill-conditioned matrix. 373## 374## @item 4 (@w{@code{GLP_EBOUND}}) 375## Invalid bounds. 376## 377## @item 5 (@w{@code{GLP_EFAIL}}) 378## Solver failed. 379## 380## @item 6 (@w{@code{GLP_EOBJLL}}) 381## Objective function lower limit reached. 382## 383## @item 7 (@w{@code{GLP_EOBJUL}}) 384## Objective function upper limit reached. 385## 386## @item 8 (@w{@code{GLP_EITLIM}}) 387## Iterations limit exhausted. 388## 389## @item 9 (@w{@code{GLP_ETMLIM}}) 390## Time limit exhausted. 391## 392## @item 10 (@w{@code{GLP_ENOPFS}}) 393## No primal feasible solution. 394## 395## @item 11 (@w{@code{GLP_ENODFS}}) 396## No dual feasible solution. 397## 398## @item 12 (@w{@code{GLP_EROOT}}) 399## Root LP optimum not provided. 400## 401## @item 13 (@w{@code{GLP_ESTOP}}) 402## Search terminated by application. 403## 404## @item 14 (@w{@code{GLP_EMIPGAP}}) 405## Relative MIP gap tolerance reached. 406## 407## @item 15 (@w{@code{GLP_ENOFEAS}}) 408## No primal/dual feasible solution. 409## 410## @item 16 (@w{@code{GLP_ENOCVG}}) 411## No convergence. 412## 413## @item 17 (@w{@code{GLP_EINSTAB}}) 414## Numerical instability. 415## 416## @item 18 (@w{@code{GLP_EDATA}}) 417## Invalid data. 418## 419## @item 19 (@w{@code{GLP_ERANGE}}) 420## Result out of range. 421## @end table 422## 423## @item extra 424## A data structure containing the following fields: 425## 426## @table @code 427## @item lambda 428## Dual variables. 429## 430## @item redcosts 431## Reduced Costs. 432## 433## @item time 434## Time (in seconds) used for solving LP/MIP problem. 435## 436## @item status 437## Status of the optimization. 438## 439## @table @asis 440## @item 1 (@w{@code{GLP_UNDEF}}) 441## Solution status is undefined. 442## 443## @item 2 (@w{@code{GLP_FEAS}}) 444## Solution is feasible. 445## 446## @item 3 (@w{@code{GLP_INFEAS}}) 447## Solution is infeasible. 448## 449## @item 4 (@w{@code{GLP_NOFEAS}}) 450## Problem has no feasible solution. 451## 452## @item 5 (@w{@code{GLP_OPT}}) 453## Solution is optimal. 454## 455## @item 6 (@w{@code{GLP_UNBND}}) 456## Problem has no unbounded solution. 457## @end table 458## @end table 459## @end table 460## 461## Example: 462## 463## @example 464## @group 465## c = [10, 6, 4]'; 466## A = [ 1, 1, 1; 467## 10, 4, 5; 468## 2, 2, 6]; 469## b = [100, 600, 300]'; 470## lb = [0, 0, 0]'; 471## ub = []; 472## ctype = "UUU"; 473## vartype = "CCC"; 474## s = -1; 475## 476## param.msglev = 1; 477## param.itlim = 100; 478## 479## [xmin, fmin, status, extra] = ... 480## glpk (c, A, b, lb, ub, ctype, vartype, s, param); 481## @end group 482## @end example 483## @end deftypefn 484 485function [xopt, fmin, errnum, extra] = glpk (c, A, b, lb, ub, ctype, vartype, sense, param) 486 487 ## If there is no input output the version and syntax 488 if (nargin < 3 || nargin > 9) 489 print_usage (); 490 endif 491 492 if (! isvector (c) || iscomplex (c) || ischar (c) || any (! isfinite (c))) 493 error ("glpk: C must be a real vector with finite values"); 494 endif 495 nx = length (c); 496 ## Force column vector. 497 c = c(:); 498 499 ## 2) Matrix constraint 500 501 if (isempty (A)) 502 error ("glpk: A cannot be an empty matrix"); 503 endif 504 if (! isreal (A)) 505 error ("glpk: A must be real valued, not %s", typeinfo (A)); 506 endif 507 if (any (! isfinite (A(:)))) 508 error ("glpk: The values in A must be finite"); 509 endif 510 511 [nc, nxa] = size (A); 512 if (nxa != nx) 513 error ("glpk: A must be %d-by-%d, not %d-by-%d", 514 nc, nx, rows (A), columns (A)); 515 endif 516 517 ## 3) RHS 518 519 if (isempty (b)) 520 error ("glpk: B cannot be an empty vector"); 521 endif 522 if (! isreal (b) || length (b) != nc) 523 error ("glpk: B must be a real-valued %d-by-1 vector", nc); 524 endif 525 if (any (! isfinite (b(:)))) 526 error ("glpk: The values in B must be finite"); 527 endif 528 529 ## 4) Vector with the lower bound of each variable 530 531 if (nargin > 3) 532 if (isempty (lb)) 533 lb = zeros (nx, 1); 534 elseif (! isreal (lb) || all (size (lb) > 1) || length (lb) != nx 535 || any (isnan (lb))) 536 error ("glpk: LB must be a real-valued %d-by-1 column vector", nx); 537 endif 538 else 539 lb = zeros (nx, 1); 540 endif 541 542 ## 5) Vector with the upper bound of each variable 543 544 if (nargin > 4) 545 if (isempty (ub)) 546 ub = Inf (nx, 1); 547 elseif (! isreal (ub) || all (size (ub) > 1) || length (ub) != nx 548 || any (isnan (ub))) 549 error ("glpk: UB must be a real-valued %d-by-1 column vector", nx); 550 endif 551 else 552 ub = Inf (nx, 1); 553 endif 554 555 ## 6) Sense of each constraint 556 557 if (nargin > 5) 558 if (isempty (ctype)) 559 ctype = repmat ("S", nc, 1); 560 elseif (! ischar (ctype) || all (size (ctype) > 1) || length (ctype) != nc) 561 error ("glpk: CTYPE must be a char vector of length %d", nc); 562 elseif (! all (ctype == "F" | ctype == "U" | ctype == "S" 563 | ctype == "L" | ctype == "D")) 564 error ("glpk: CTYPE must contain only F, U, S, L, or D"); 565 endif 566 else 567 ctype = repmat ("S", nc, 1); 568 endif 569 570 ## 7) Vector with the type of variables 571 572 if (nargin > 6) 573 if (isempty (vartype)) 574 vartype = repmat ("C", nx, 1); 575 elseif (! ischar (vartype) || all (size (vartype) > 1) 576 || length (vartype) != nx) 577 error ("glpk: VARTYPE must be a char vector of length %d", nx); 578 elseif (! all (vartype == "C" | vartype == "I")) 579 error ("glpk: VARTYPE must contain only C or I"); 580 endif 581 else 582 ## As default we consider continuous vars 583 vartype = repmat ("C", nx, 1); 584 endif 585 586 ## 8) Sense of optimization 587 588 if (nargin > 7) 589 if (isempty (sense)) 590 sense = 1; 591 elseif (ischar (sense) || all (size (sense) > 1) || ! isreal (sense) 592 || any (! isfinite (sense))) 593 error ("glpk: SENSE must be an integer value"); 594 elseif (sense >= 0) 595 sense = 1; 596 else 597 sense = -1; 598 endif 599 else 600 sense = 1; 601 endif 602 603 ## 9) Parameters vector 604 605 if (nargin > 8) 606 if (! isstruct (param)) 607 error ("glpk: PARAM must be a structure"); 608 endif 609 else 610 param = struct (); 611 endif 612 613 [xopt, fmin, errnum, extra] = ... 614 __glpk__ (c, A, b, lb, ub, ctype, vartype, sense, param); 615 616endfunction 617 618 619%!testif HAVE_GLPK 620%! sense = -1; 621%! c = [10, 6, 4]'; 622%! A = [1, 1, 1; 10, 4, 5; 2, 2, 6]; 623%! b = [100, 600, 300]'; 624%! ctype = ['U', 'U', 'U']'; 625%! lb = [0, 0, 0]'; 626%! ub = []; 627%! vartype = ['C', 'C', 'C']'; 628%! param.msglev = 0; 629%! param.lpsolver = 1; 630%! [xmin, fmin, errnum, extra] = glpk (c, A, b, lb, ub, ctype, vartype, ... 631%! sense, param); 632%! assert (fmin, c' * xmin); 633%! for i = 1:3 634%! assert (A(i,:) * xmin <= b(i)); 635%! endfor 636 637%!testif HAVE_GLPK 638%! sense = 1; 639%! c = [-1, -1]'; 640%! A = [-2, 5; 2, -2]; 641%! b = [5, 1]'; 642%! ctype = ['U', 'U']'; 643%! lb = [0, 0]'; 644%! ub = []; 645%! vartype = ['I', 'I']'; 646%! param.msglev = 0; 647%! [xmin, fmin, errnum, extra] = glpk (c, A, b, lb, ub, ctype, vartype, ... 648%! sense, param); 649%! assert (fmin, c' * xmin); 650%! for i = 1:2 651%! assert (A(i,:) * xmin <= b(i)); 652%! endfor 653 654 655%!testif HAVE_GLPK 656%! sense = 1; 657%! c = [0, 0, 0, -1, -1]'; 658%! A = [-2, 0, 0, 1, 0; 0, 1, 0, 0, 2; 0, 0, 1, 3, 2]; 659%! b = [4, 12, 18]'; 660%! ctype = ['S', 'S', 'S']'; 661%! lb = [0, 0, 0, 0, 0]'; 662%! ub = []; 663%! vartype = ['C', 'C', 'C', 'C', 'C']'; 664%! param.msglev = 0; 665%! [xmin, fmin, errnum, extra] = glpk (c, A, b, lb, ub, ctype, vartype, ... 666%! sense, param); 667%! assert (fmin, c' * xmin); 668%! assert (A * xmin, b); 669 670%!error<C .* finite values> glpk(NaN, 2, 3) 671%!error<A must be finite> glpk(1, NaN, 3) 672%!error<B must be finite> glpk(1, 2, NaN) 673%!error<LB must be a real-valued> glpk(1, 2, 3, NaN) 674%!error<UB must be a real-valued> glpk(1, 2, 3, 4, NaN) 675%!error<SENSE must be .* integer> glpk(1, 2, 3, 4, 5, "F", "C", NaN) 676