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