--- /dev/null
+## Copyright (C) 2005-2012 Nicolo' Giorgetti
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @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})
+## Solve a linear program using the GNU @sc{glpk} library. Given three
+## arguments, @code{glpk} solves the following standard LP:
+## @tex
+## $$
+## \min_x C^T x
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## min C'*x
+## @end example
+##
+## @end ifnottex
+## subject to
+## @tex
+## $$
+## Ax = b \qquad x \geq 0
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+## A*x = b
+## x >= 0
+## @end group
+## @end example
+##
+## @end ifnottex
+## but may also solve problems of the form
+## @tex
+## $$
+## [ \min_x | \max_x ] C^T x
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## [ min | max ] C'*x
+## @end example
+##
+## @end ifnottex
+## subject to
+## @tex
+## $$
+## Ax [ = | \leq | \geq ] b \qquad LB \leq x \leq UB
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+## A*x [ "=" | "<=" | ">=" ] b
+## x >= LB
+## x <= UB
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## Input arguments:
+##
+## @table @var
+## @item c
+## A column array containing the objective function coefficients.
+##
+## @item A
+## A matrix containing the constraints coefficients.
+##
+## @item b
+## A column array containing the right-hand side value for each constraint
+## in the constraint matrix.
+##
+## @item lb
+## An array containing the lower bound on each of the variables. If
+## @var{lb} is not supplied, the default lower bound for the variables is
+## zero.
+##
+## @item ub
+## An array containing the upper bound on each of the variables. If
+## @var{ub} is not supplied, the default upper bound is assumed to be
+## infinite.
+##
+## @item ctype
+## An array of characters containing the sense of each constraint in the
+## constraint matrix. Each element of the array may be one of the
+## following values
+## @table @asis
+## @item "F"
+## A free (unbounded) constraint (the constraint is ignored).
+##
+## @item "U"
+## An inequality constraint with an upper bound (@code{A(i,:)*x <= b(i)}).
+##
+## @item "S"
+## An equality constraint (@code{A(i,:)*x = b(i)}).
+##
+## @item "L"
+## An inequality with a lower bound (@code{A(i,:)*x >= b(i)}).
+##
+## @item "D"
+## An inequality constraint with both upper and lower bounds
+## (@code{A(i,:)*x >= -b(i)} @emph{and} (@code{A(i,:)*x <= b(i)}).
+## @end table
+##
+## @item vartype
+## A column array containing the types of the variables.
+## @table @asis
+## @item "C"
+## A continuous variable.
+##
+## @item "I"
+## An integer variable.
+## @end table
+##
+## @item sense
+## If @var{sense} is 1, the problem is a minimization. If @var{sense} is
+## -1, the problem is a maximization. The default value is 1.
+##
+## @item param
+## A structure containing the following parameters used to define the
+## behavior of solver. Missing elements in the structure take on default
+## values, so you only need to set the elements that you wish to change
+## from the default.
+##
+## Integer parameters:
+##
+## @table @code
+## @item msglev (@w{@code{LPX_K_MSGLEV}}, default: 1)
+## Level of messages output by solver routines:
+## @table @asis
+## @item 0
+## No output.
+##
+## @item 1
+## Error messages only.
+##
+## @item 2
+## Normal output.
+##
+## @item 3
+## Full output (includes informational messages).
+## @end table
+##
+## @item scale (@w{@code{LPX_K_SCALE}}, default: 1)
+## Scaling option:
+## @table @asis
+## @item 0
+## No scaling.
+##
+## @item 1
+## Equilibration scaling.
+##
+## @item 2
+## Geometric mean scaling, then equilibration scaling.
+## @end table
+##
+## @item dual (@w{@code{LPX_K_DUAL}}, default: 0)
+## Dual simplex option:
+## @table @asis
+## @item 0
+## Do not use the dual simplex.
+##
+## @item 1
+## If initial basic solution is dual feasible, use the dual simplex.
+## @end table
+##
+## @item price (@w{@code{LPX_K_PRICE}}, default: 1)
+## Pricing option (for both primal and dual simplex):
+## @table @asis
+## @item 0
+## Textbook pricing.
+##
+## @item 1
+## Steepest edge pricing.
+## @end table
+##
+## @item round (@w{@code{LPX_K_ROUND}}, default: 0)
+## Solution rounding option:
+## @table @asis
+## @item 0
+## Report all primal and dual values "as is".
+##
+## @item 1
+## Replace tiny primal and dual values by exact zero.
+## @end table
+##
+## @item itlim (@w{@code{LPX_K_ITLIM}}, default: -1)
+## Simplex iterations limit. If this value is positive, it is decreased by
+## one each time when one simplex iteration has been performed, and
+## reaching zero value signals the solver to stop the search. Negative
+## value means no iterations limit.
+##
+## @item itcnt (@w{@code{LPX_K_OUTFRQ}}, default: 200)
+## Output frequency, in iterations. This parameter specifies how
+## frequently the solver sends information about the solution to the
+## standard output.
+##
+## @item branch (@w{@code{LPX_K_BRANCH}}, default: 2)
+## Branching heuristic option (for MIP only):
+## @table @asis
+## @item 0
+## Branch on the first variable.
+##
+## @item 1
+## Branch on the last variable.
+##
+## @item 2
+## Branch using a heuristic by Driebeck and Tomlin.
+## @end table
+##
+## @item btrack (@w{@code{LPX_K_BTRACK}}, default: 2)
+## Backtracking heuristic option (for MIP only):
+## @table @asis
+## @item 0
+## Depth first search.
+##
+## @item 1
+## Breadth first search.
+##
+## @item 2
+## Backtrack using the best projection heuristic.
+## @end table
+##
+## @item presol (@w{@code{LPX_K_PRESOL}}, default: 1)
+## If this flag is set, the routine lpx_simplex solves the problem using
+## the built-in LP presolver. Otherwise the LP presolver is not used.
+##
+## @item lpsolver (default: 1)
+## Select which solver to use. If the problem is a MIP problem this flag
+## will be ignored.
+## @table @asis
+## @item 1
+## Revised simplex method.
+##
+## @item 2
+## Interior point method.
+## @end table
+##
+## @item save (default: 0)
+## If this parameter is nonzero, save a copy of the problem in
+## CPLEX LP format to the file @file{"outpb.lp"}. There is currently no
+## way to change the name of the output file.
+## @end table
+##
+## Real parameters:
+##
+## @table @code
+## @item relax (@w{@code{LPX_K_RELAX}}, default: 0.07)
+## Relaxation parameter used in the ratio test. If it is zero, the textbook
+## ratio test is used. If it is non-zero (should be positive), Harris'
+## two-pass ratio test is used. In the latter case on the first pass of the
+## ratio test basic variables (in the case of primal simplex) or reduced
+## costs of non-basic variables (in the case of dual simplex) are allowed
+## to slightly violate their bounds, but not more than
+## @code{relax*tolbnd} or @code{relax*toldj (thus, @code{relax} is a
+## percentage of @code{tolbnd} or @code{toldj}}.
+##
+## @item tolbnd (@w{@code{LPX_K_TOLBND}}, default: 10e-7)
+## Relative tolerance used to check if the current basic solution is primal
+## feasible. It is not recommended that you change this parameter unless you
+## have a detailed understanding of its purpose.
+##
+## @item toldj (@w{@code{LPX_K_TOLDJ}}, default: 10e-7)
+## Absolute tolerance used to check if the current basic solution is dual
+## feasible. It is not recommended that you change this parameter unless you
+## have a detailed understanding of its purpose.
+##
+## @item tolpiv (@w{@code{LPX_K_TOLPIV}}, default: 10e-9)
+## Relative tolerance used to choose eligible pivotal elements of the
+## simplex table. It is not recommended that you change this parameter unless
+## you have a detailed understanding of its purpose.
+##
+## @item objll (@w{@code{LPX_K_OBJLL}}, default: -DBL_MAX)
+## Lower limit of the objective function. If on the phase II the objective
+## function reaches this limit and continues decreasing, the solver stops
+## the search. This parameter is used in the dual simplex method only.
+##
+## @item objul (@w{@code{LPX_K_OBJUL}}, default: +DBL_MAX)
+## Upper limit of the objective function. If on the phase II the objective
+## function reaches this limit and continues increasing, the solver stops
+## the search. This parameter is used in the dual simplex only.
+##
+## @item tmlim (@w{@code{LPX_K_TMLIM}}, default: -1.0)
+## Searching time limit, in seconds. If this value is positive, it is
+## decreased each time when one simplex iteration has been performed by the
+## amount of time spent for the iteration, and reaching zero value signals
+## the solver to stop the search. Negative value means no time limit.
+##
+## @item outdly (@w{@code{LPX_K_OUTDLY}}, default: 0.0)
+## Output delay, in seconds. This parameter specifies how long the solver
+## should delay sending information about the solution to the standard
+## output. Non-positive value means no delay.
+##
+## @item tolint (@w{@code{LPX_K_TOLINT}}, default: 10e-5)
+## Relative tolerance used to check if the current basic solution is integer
+## feasible. It is not recommended that you change this parameter unless
+## you have a detailed understanding of its purpose.
+##
+## @item tolobj (@w{@code{LPX_K_TOLOBJ}}, default: 10e-7)
+## Relative tolerance used to check if the value of the objective function
+## is not better than in the best known integer feasible solution. It is
+## not recommended that you change this parameter unless you have a
+## detailed understanding of its purpose.
+## @end table
+## @end table
+##
+## Output values:
+##
+## @table @var
+## @item xopt
+## The optimizer (the value of the decision variables at the optimum).
+##
+## @item fopt
+## The optimum value of the objective function.
+##
+## @item status
+## Status of the optimization.
+##
+## Simplex Method:
+## @table @asis
+## @item 180 (@w{@code{LPX_OPT}})
+## Solution is optimal.
+##
+## @item 181 (@w{@code{LPX_FEAS}})
+## Solution is feasible.
+##
+## @item 182 (@w{@code{LPX_INFEAS}})
+## Solution is infeasible.
+##
+## @item 183 (@w{@code{LPX_NOFEAS}})
+## Problem has no feasible solution.
+##
+## @item 184 (@w{@code{LPX_UNBND}})
+## Problem has no unbounded solution.
+##
+## @item 185 (@w{@code{LPX_UNDEF}})
+## Solution status is undefined.
+## @end table
+## Interior Point Method:
+## @table @asis
+## @item 150 (@w{@code{LPX_T_UNDEF}})
+## The interior point method is undefined.
+##
+## @item 151 (@w{@code{LPX_T_OPT}})
+## The interior point method is optimal.
+## @end table
+## Mixed Integer Method:
+## @table @asis
+## @item 170 (@w{@code{LPX_I_UNDEF}})
+## The status is undefined.
+##
+## @item 171 (@w{@code{LPX_I_OPT}})
+## The solution is integer optimal.
+##
+## @item 172 (@w{@code{LPX_I_FEAS}})
+## Solution integer feasible but its optimality has not been proven
+##
+## @item 173 (@w{@code{LPX_I_NOFEAS}})
+## No integer feasible solution.
+## @end table
+## @noindent
+## If an error occurs, @var{status} will contain one of the following
+## codes:
+##
+## @table @asis
+## @item 204 (@w{@code{LPX_E_FAULT}})
+## Unable to start the search.
+##
+## @item 205 (@w{@code{LPX_E_OBJLL}})
+## Objective function lower limit reached.
+##
+## @item 206 (@w{@code{LPX_E_OBJUL}})
+## Objective function upper limit reached.
+##
+## @item 207 (@w{@code{LPX_E_ITLIM}})
+## Iterations limit exhausted.
+##
+## @item 208 (@w{@code{LPX_E_TMLIM}})
+## Time limit exhausted.
+##
+## @item 209 (@w{@code{LPX_E_NOFEAS}})
+## No feasible solution.
+##
+## @item 210 (@w{@code{LPX_E_INSTAB}})
+## Numerical instability.
+##
+## @item 211 (@w{@code{LPX_E_SING}})
+## Problems with basis matrix.
+##
+## @item 212 (@w{@code{LPX_E_NOCONV}})
+## No convergence (interior).
+##
+## @item 213 (@w{@code{LPX_E_NOPFS}})
+## No primal feasible solution (LP presolver).
+##
+## @item 214 (@w{@code{LPX_E_NODFS}})
+## No dual feasible solution (LP presolver).
+## @end table
+##
+## @item extra
+## A data structure containing the following fields:
+## @table @code
+## @item lambda
+## Dual variables.
+##
+## @item redcosts
+## Reduced Costs.
+##
+## @item time
+## Time (in seconds) used for solving LP/MIP problem.
+##
+## @item mem
+## Memory (in bytes) used for solving LP/MIP problem (this is not
+## available if the version of @sc{glpk} is 4.15 or later).
+## @end table
+## @end table
+##
+## Example:
+##
+## @example
+## @group
+## c = [10, 6, 4]';
+## A = [ 1, 1, 1;
+## 10, 4, 5;
+## 2, 2, 6];
+## b = [100, 600, 300]';
+## lb = [0, 0, 0]';
+## ub = [];
+## ctype = "UUU";
+## vartype = "CCC";
+## s = -1;
+##
+## param.msglev = 1;
+## param.itlim = 100;
+##
+## [xmin, fmin, status, extra] = ...
+## glpk (c, A, b, lb, ub, ctype, vartype, s, param);
+## @end group
+## @end example
+## @end deftypefn
+
+## Author: Nicolo' Giorgetti <giorgetti@dii.unisi.it>
+## Adapted-by: jwe
+
+function [xopt, fmin, status, extra] = glpk (c, A, b, lb, ub, ctype, vartype, sense, param)
+
+ ## If there is no input output the version and syntax
+ if (nargin < 3 || nargin > 9)
+ print_usage ();
+ return;
+ endif
+
+ if (all (size (c) > 1) || iscomplex (c) || ischar (c))
+ error ("glpk:C must be a real vector");
+ return;
+ endif
+ nx = length (c);
+ ## Force column vector.
+ c = c(:);
+
+ ## 2) Matrix constraint
+
+ if (isempty (A))
+ error ("glpk: A cannot be an empty matrix");
+ return;
+ endif
+ [nc, nxa] = size(A);
+ if (! isreal (A) || nxa != nx)
+ error ("glpk: A must be a real valued %d by %d matrix", nc, nx);
+ return;
+ endif
+
+ ## 3) RHS
+
+ if (isempty (b))
+ error ("glpk: B cannot be an empty vector");
+ return;
+ endif
+ if (! isreal (b) || length (b) != nc)
+ error ("glpk: B must be a real valued %d by 1 vector", nc);
+ return;
+ endif
+
+ ## 4) Vector with the lower bound of each variable
+
+ if (nargin > 3)
+ if (isempty (lb))
+ lb = zeros (nx, 1);
+ elseif (! isreal (lb) || all (size (lb) > 1) || length (lb) != nx)
+ error ("glpk: LB must be a real valued %d by 1 column vector", nx);
+ return;
+ endif
+ else
+ lb = zeros (nx, 1);
+ endif
+
+ ## 5) Vector with the upper bound of each variable
+
+ if (nargin > 4)
+ if (isempty (ub))
+ ub = Inf (nx, 1);
+ elseif (! isreal (ub) || all (size (ub) > 1) || length (ub) != nx)
+ error ("glpk: UB must be a real valued %d by 1 column vector", nx);
+ return;
+ endif
+ else
+ ub = Inf (nx, 1);
+ endif
+
+ ## 6) Sense of each constraint
+
+ if (nargin > 5)
+ if (isempty (ctype))
+ ctype = repmat ("S", nc, 1);
+ elseif (! ischar (ctype) || all (size (ctype) > 1) || length (ctype) != nc)
+ error ("glpk: CTYPE must be a char valued vector of length %d", nc);
+ return;
+ elseif (! all (ctype == "F" | ctype == "U" | ctype == "S"
+ | ctype == "L" | ctype == "D"))
+ error ("glpk: CTYPE must contain only F, U, S, L, or D");
+ return;
+ endif
+ else
+ ctype = repmat ("S", nc, 1);
+ endif
+
+ ## 7) Vector with the type of variables
+
+ if (nargin > 6)
+ if (isempty (vartype))
+ vartype = repmat ("C", nx, 1);
+ elseif (! ischar (vartype) || all (size (vartype) > 1)
+ || length (vartype) != nx)
+ error ("glpk: VARTYPE must be a char valued vector of length %d", nx);
+ return;
+ elseif (! all (vartype == "C" | vartype == "I"))
+ error ("glpk: VARTYPE must contain only C or I");
+ return;
+ endif
+ else
+ ## As default we consider continuous vars
+ vartype = repmat ("C", nx, 1);
+ endif
+
+ ## 8) Sense of optimization
+
+ if (nargin > 7)
+ if (isempty (sense))
+ sense = 1;
+ elseif (ischar (sense) || all (size (sense) > 1) || ! isreal (sense))
+ error ("glpk: SENSE must be an integer value");
+ elseif (sense >= 0)
+ sense = 1;
+ else
+ sense = -1;
+ endif
+ else
+ sense = 1;
+ endif
+
+ ## 9) Parameters vector
+
+ if (nargin > 8)
+ if (! isstruct (param))
+ error ("glpk: PARAM must be a structure");
+ return;
+ endif
+ else
+ param = struct ();
+ endif
+
+ [xopt, fmin, status, extra] = ...
+ __glpk__ (c, A, b, lb, ub, ctype, vartype, sense, param);
+
+endfunction
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff