octave_packages/optim-1.2.0/linprog.m

   1 ## Copyright (C) 2009 Luca Favatella <slackydeb@gmail.com>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn{Function File} {@var{x} =} linprog (@var{f}, @var{A}, @var{b})
  18 ## @deftypefnx{Function File} {@var{x} =} linprog (@var{f}, @var{A}, @var{b}, @var{Aeq}, @var{beq})
  19 ## @deftypefnx{Function File} {@var{x} =} linprog (@var{f}, @var{A}, @var{b}, @var{Aeq}, @var{beq}, @var{lb}, @var{ub})
  20 ## @deftypefnx{Function File} {[@var{x}, @var{fval}] =} linprog (@dots{})
  21 ## Solve a linear problem.
  22 ##
  23 ## Finds
  24 ##
  25 ## @example
  26 ## min (f' * x)
  27 ## @end example
  28 ##
  29 ## (both f and x are column vectors) subject to
  30 ##
  31 ## @example
  32 ## @group
  33 ## A   * x <= b
  34 ## Aeq * x  = beq
  35 ## lb <= x <= ub
  36 ## @end group
  37 ## @end example
  38 ##
  39 ## If not specified, @var{Aeq} and @var{beq} default to empty matrices.
  40 ##
  41 ## If not specified, the lower bound @var{lb} defaults to minus infinite
  42 ## and the upper bound @var{ub} defaults to infinite.
  43 ##
  44 ## @seealso{glpk}
  45 ## @end deftypefn
  46
  47 function [x fval] = linprog (f, A, b,
  48                              Aeq = [], beq = [],
  49                              lb = [], ub = [])
  50
  51   if (((nargin != 3) && (nargin != 5) && (nargin != 7)) ||
  52       (nargout > 2))
  53     print_usage ();
  54   endif
  55
  56   nr_f = rows(f);
  57
  58   # Sanitize A and b
  59   if (isempty (A) && isempty (b))
  60     A = zeros (0, nr_f);
  61     b = zeros (rows (A), 1);
  62   endif
  63
  64   nr_A = rows (A);
  65
  66   if (columns (f) != 1)
  67     error ("f must be a column vector");
  68   elseif (columns (A) != nr_f)
  69     error ("columns (A) != rows (f)");
  70   elseif (size (b) != [nr_A 1])
  71     error ("size (b) != [(rows (A)) 1]");
  72   else
  73
  74     ## Sanitize Aeq
  75     if (isempty (Aeq))
  76       Aeq = zeros (0, nr_f);
  77     endif
  78     if (columns (Aeq) != nr_f)
  79       error ("columns (Aeq) != rows (f)");
  80     endif
  81
  82     ## Sanitize beq
  83     if (isempty (beq))
  84       beq = zeros (0, 1);
  85     endif
  86     nr_Aeq = rows (Aeq);
  87     if (size (beq) != [nr_Aeq 1])
  88       error ("size (beq) != [(rows (Aeq)) 1]");
  89     endif
  90
  91     ## Sanitize lb
  92     if (isempty (lb))
  93       lb = - Inf (nr_f, 1);
  94     endif
  95     if (size (lb) != [nr_f 1])
  96       error ("size (lb) != [(rows (f)) 1]");
  97     endif
  98
  99     ## Sanitize ub
 100     if (isempty (ub))
 101       ub = Inf (nr_f, 1);
 102     endif
 103     if (size (ub) != [nr_f 1])
 104       error ("size (ub) != [(rows (f)) 1]");
 105     endif
 106
 107
 108     ## Call glpk
 109     ctype = [(repmat ("U", nr_A, 1));
 110              (repmat ("S", nr_Aeq, 1))];
 111     [x(1:nr_f, 1) fval(1, 1)] = glpk (f, [A; Aeq], [b; beq], lb, ub, ctype);
 112
 113   endif
 114
 115 endfunction
 116
 117 %!test
 118 %! f = [1; -1];
 119 %! A = [];
 120 %! b = [];
 121 %! Aeq = [1, 0];
 122 %! beq = [2];
 123 %! lb = [0; Inf];
 124 %! ub = [-Inf; 0];
 125 %! x_exp = [2; 0];
 126 %! assert (linprog (f, A, b, Aeq, beq, lb, ub), x_exp);
 127
 128 %!shared f, A, b, lb, ub, x_exp, fval_exp
 129 %! f  = [21 25 31 34  23 19 32  36 27 25 19]';
 130 %!
 131 %! A1 = [ 1  0  0  0   1  0  0   1  0  0  0;
 132 %!        0  1  0  0   0  1  0   0  1  0  0;
 133 %!        0  0  1  0   0  0  0   0  0  1  0;
 134 %!        0  0  0  1   0  0  1   0  0  0  1];
 135 %! A2 = [ 1  1  1  1   0  0  0   0  0  0  0;
 136 %!        0  0  0  0   1  1  1   0  0  0  0;
 137 %!        0  0  0  0   0  0  0   1  1  1  1];
 138 %! A  = [-A1; A2];
 139 %!
 140 %! b1 = [40; 50; 50; 70];
 141 %! b2 = [100; 60; 50];
 142 %! b  = [-b1; b2];
 143 %!
 144 %! lb = zeros (rows (f), 1);
 145 %! ub =   Inf (rows (f), 1);
 146 %!
 147 %! x_exp    = [40 0 50 10 0 50 10 0 0 0 50]';
 148 %! fval_exp = f' * x_exp;
 149 %!
 150 %!test
 151 %! Aeq = [];
 152 %! beq = [];
 153 %! [x_obs fval_obs] = linprog (f, A, b, Aeq, beq, lb, ub);
 154 %! assert ([x_obs; fval_obs], [x_exp; fval_exp]);
 155 %!
 156 %!test
 157 %! Aeq = zeros (1, rows (f));
 158 %! beq = 0;
 159 %! assert (linprog (f, A, b, Aeq, beq, lb, ub), x_exp);