|
5233
|
1 ## Copyright (C) 2005 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 2, or (at your option) |
|
|
8 ## 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, write to the Free |
|
|
17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA |
|
|
18 ## 02111-1307, USA. |
|
|
19 |
|
5232
|
20 ## GLPK - An Octave Interface for the GNU GLPK library |
|
|
21 ## |
|
|
22 ## This routine calls the glpk library to solve an LP/MIP problem. A typical |
|
|
23 ## LP problem has following structure: |
|
|
24 ## |
|
|
25 ## [min|max] C'x |
|
|
26 ## s.t. |
|
|
27 ## Ax ["="|"<="|">="] b |
|
|
28 ## {x <= UB} |
|
|
29 ## {x >= LB} |
|
|
30 ## |
|
|
31 ## The calling syntax is: |
|
5233
|
32 ## [XOPT,FOPT,STATUS,EXTRA]=glpk(SENSE,C,... |
|
|
33 ## A,B,CTYPE,LB,UB,... |
|
|
34 ## VARTYPE,PARAM,LPSOLVER,SAVE) |
|
5232
|
35 ## |
|
|
36 ## For a quick reference to the syntax just type glpk at command prompt. |
|
|
37 ## |
|
|
38 ## The minimum number of input arguments is 4 (SENSE,C,A,B). In this case we |
|
|
39 ## assume all the constraints are '<=' and all the variables are continuous. |
|
|
40 ## |
|
|
41 ## --- INPUTS --- |
|
|
42 ## |
|
|
43 ## SENSE: indicates whether the problem is a minimization |
|
|
44 ## or maximization problem. |
|
|
45 ## SENSE = 1 minimize |
|
|
46 ## SENSE = -1 maximize. |
|
|
47 ## |
|
|
48 ## C: A column array containing the objective function |
|
|
49 ## coefficients. |
|
|
50 ## |
|
|
51 ## A: A matrix containing the constraints coefficients. |
|
|
52 ## |
|
|
53 ## B: A column array containing the right-hand side value for |
|
|
54 ## each constraint in the constraint matrix. |
|
|
55 ## |
|
|
56 ## CTYPE: A column array containing the sense of each constraint |
|
|
57 ## in the constraint matrix. |
|
|
58 ## CTYPE(i) = 'F' Free (unbounded) variable |
|
|
59 ## CTYPE(i) = 'U' "<=" Variable with upper bound |
|
|
60 ## CTYPE(i) = 'S' "=" Fixed Variable |
|
|
61 ## CTYPE(i) = 'L' ">=" Variable with lower bound |
|
|
62 ## CTYPE(i) = 'D' Double-bounded variable |
|
|
63 ## (This is case sensitive). |
|
|
64 ## |
|
|
65 ## LB: An array of at least length numcols containing the lower |
|
|
66 ## bound on each of the variables. |
|
|
67 ## |
|
|
68 ## UB: An array of at least length numcols containing the upper |
|
|
69 ## bound on each of the variables. |
|
|
70 ## |
|
|
71 ## VARTYPE: A column array containing the types of the variables. |
|
|
72 ## VARTYPE(i) = 'C' continuous variable |
|
|
73 ## VARTYPE(i) = 'I' Integer variable |
|
|
74 ## (This is case sensitive). |
|
|
75 ## |
|
|
76 ## PARAM: A structure containing some parameters used to define |
|
|
77 ## the behavior of solver. For more details type |
|
|
78 ## HELP GLPKPARAMS. |
|
|
79 ## |
|
|
80 ## LPSOLVER: Selects which solver using to solve LP problems. |
|
|
81 ## LPSOLVER=1 Revised Simplex Method |
|
|
82 ## LPSOLVER=2 Interior Point Method |
|
|
83 ## If the problem is a MIP problem this flag will be ignored. |
|
|
84 ## |
|
|
85 ## SAVE: Saves a copy of the problem if SAVE<>0. |
|
|
86 ## The file name can not be specified and defaults to "outpb.lp". |
|
|
87 ## The output file is CPLEX LP format. |
|
|
88 ## |
|
|
89 ## --- OUTPUTS --- |
|
|
90 ## |
|
|
91 ## XOPT: The optimizer. |
|
|
92 ## |
|
|
93 ## FOPT: The optimum. |
|
|
94 ## |
|
|
95 ## STATUS: Status of the optimization. |
|
5233
|
96 ## |
|
5232
|
97 ## - Simplex Method - |
|
5233
|
98 ## Value Code |
|
|
99 ## 180 LPX_OPT solution is optimal |
|
5232
|
100 ## 181 LPX_FEAS solution is feasible |
|
|
101 ## 182 LPX_INFEAS solution is infeasible |
|
|
102 ## 183 LPX_NOFEAS problem has no feasible solution |
|
|
103 ## 184 LPX_UNBND problem has no unbounded solution |
|
|
104 ## 185 LPX_UNDEF solution status is undefined |
|
|
105 ## |
|
|
106 ## - Interior Point Method - |
|
|
107 ## Value Code |
|
|
108 ## 150 LPX_T_UNDEF the interior point method is undefined |
|
|
109 ## 151 LPX_T_OPT the interior point method is optimal |
|
|
110 ## * Note that additional status codes may appear in |
|
|
111 ## the future versions of this routine * |
|
5233
|
112 ## |
|
5232
|
113 ## - Mixed Integer Method - |
|
|
114 ## Value Code |
|
|
115 ## 170 LPX_I_UNDEF the status is undefined |
|
|
116 ## 171 LPX_I_OPT the solution is integer optimal |
|
|
117 ## 172 LPX_I_FEAS solution integer feasible but |
|
|
118 ## its optimality has not been proven |
|
|
119 ## 173 LPX_I_NOFEAS no integer feasible solution |
|
5233
|
120 ## |
|
5232
|
121 ## EXTRA: A data structure containing the following fields: |
|
|
122 ## LAMBDA Dual variables |
|
|
123 ## REDCOSTS Reduced Costs |
|
|
124 ## TIME Time (in seconds) used for solving LP/MIP problem in seconds. |
|
|
125 ## MEM Memory (in bytes) used for solving LP/MIP problem. |
|
|
126 ## |
|
|
127 ## |
|
5233
|
128 ## In case of error the glpk returns one of the |
|
5232
|
129 ## following codes (these codes are in STATUS). For more informations on |
|
|
130 ## the causes of these codes refer to the GLPK reference manual. |
|
|
131 ## |
|
|
132 ## Value Code |
|
|
133 ## 204 LPX_E_FAULT unable to start the search |
|
|
134 ## 205 LPX_E_OBJLL objective function lower limit reached |
|
|
135 ## 206 LPX_E_OBJUL objective function upper limit reached |
|
|
136 ## 207 LPX_E_ITLIM iterations limit exhausted |
|
|
137 ## 208 LPX_E_TMLIM time limit exhausted |
|
|
138 ## 209 LPX_E_NOFEAS no feasible solution |
|
|
139 ## 210 LPX_E_INSTAB numerical instability |
|
|
140 ## 211 LPX_E_SING problems with basis matrix |
|
|
141 ## 212 LPX_E_NOCONV no convergence (interior) |
|
|
142 ## 213 LPX_E_NOPFS no primal feas. sol. (LP presolver) |
|
|
143 ## 214 LPX_E_NODFS no dual feas. sol. (LP presolver) |
|
|
144 |
|
5233
|
145 ## Author: Nicolo' Giorgetti <giorgetti@dii.unisi.it> |
|
|
146 ## Adapted-by: jwe |
|
5232
|
147 |
|
5233
|
148 function [xopt, fmin, status, extra] = glpk (sense, c, a, b, ctype, lb, ub, vartype, param, lpsolver, savepb) |
|
5232
|
149 |
|
5233
|
150 ## If there is no input output the version and syntax |
|
|
151 if (nargin < 4 || nargin > 11) |
|
|
152 usage ("[xopt, fopt, status, extra] = glpk (sense, c, a, b, ctype, lb, ub, vartype, param, lpsolver, savepb"); |
|
|
153 return; |
|
|
154 endif |
|
5232
|
155 |
|
5233
|
156 if (all (size (c) > 1) || iscomplex (c) || ischar (c)) |
|
|
157 error ("C must be a real vector"); |
|
|
158 return; |
|
|
159 endif |
|
|
160 nx = length (c); |
|
|
161 ## Force column vector. |
|
|
162 c = c(:); |
|
5232
|
163 |
|
5233
|
164 ## 3) Matrix constraint |
|
5232
|
165 |
|
5233
|
166 if (isempty (a)) |
|
|
167 error ("A cannot be an empty matrix"); |
|
|
168 return; |
|
|
169 endif |
|
|
170 [nc, nxa] = size(a); |
|
|
171 if (! isreal (a) || nxa != nx) |
|
|
172 error ("A must be a real valued %d by %d matrix", nc, nx); |
|
|
173 return; |
|
|
174 endif |
|
|
175 |
|
|
176 ## 4) RHS |
|
5232
|
177 |
|
5233
|
178 if (isempty (b)) |
|
|
179 error ("B cannot be an empty vector"); |
|
|
180 return; |
|
|
181 endif |
|
|
182 if (! isreal (b) || length (b) != nc) |
|
|
183 error ("B must be a real valued %d by 1 vector", nc); |
|
|
184 return; |
|
|
185 endif |
|
|
186 |
|
|
187 ## 5) Sense of each constraint |
|
5232
|
188 |
|
5233
|
189 if (nargin > 4) |
|
|
190 if (isempty (ctype)) |
|
|
191 ctype = repmat ("U", nc, 1); |
|
|
192 elseif (! ischar (ctype) || all (size (ctype) > 1) || length (ctype) != nc) |
|
|
193 error ("CTYPE must be a char valued %d by 1 column vector", nc); |
|
|
194 return; |
|
|
195 elseif (! all (ctype == "F" | ctype == "U" | ctype == "S" |
|
|
196 | ctype == "L" | ctype == "D")) |
|
|
197 error ("CTYPE must contain only F, U, S, L, or D"); |
|
|
198 return; |
|
|
199 endif |
|
|
200 else |
|
|
201 ctype = repmat ("U", nc, 1); |
|
|
202 end |
|
|
203 |
|
|
204 ## 6) Vector with the lower bound of each variable |
|
5232
|
205 |
|
5233
|
206 if (nargin > 5) |
|
|
207 if (isempty (lb)) |
|
|
208 lb = repmat (-Inf, nx, 1); |
|
|
209 elseif (! isreal (lb) || all (size (lb) > 1) || length (lb) != nx) |
|
|
210 error ("LB must be a real valued %d by 1 column vector", nx); |
|
|
211 return; |
|
|
212 endif |
|
|
213 else |
|
|
214 lb = repmat (-Inf, nx, 1); |
|
|
215 end |
|
|
216 |
|
|
217 ## 7) Vector with the upper bound of each variable |
|
5232
|
218 |
|
5233
|
219 if (nargin > 6) |
|
|
220 if (isempty (ub)) |
|
|
221 ub = repmat (Inf, nx, 1); |
|
|
222 elseif (! isreal (ub) || all (size (ub) > 1) || length (ub) != nx) |
|
|
223 error ("UB must be a real valued %d by 1 column vector", nx); |
|
|
224 return; |
|
|
225 endif |
|
|
226 else |
|
|
227 ub = repmat (Inf, nx, 1); |
|
|
228 end |
|
5232
|
229 |
|
5233
|
230 ## 8) Vector with the type of variables |
|
5232
|
231 |
|
5233
|
232 if (nargin > 7) |
|
|
233 if isempty (vartype) |
|
|
234 vartype = repmat ("C", nx, 1); |
|
|
235 elseif (! ischar (vartype) || all (size (vartype) > 1) |
|
|
236 || length (vartype) != nx) |
|
|
237 error ("VARTYPE must be a char valued %d by 1 column vector", nx); |
|
|
238 return; |
|
|
239 elseif (! all (vartype == "C" | vartype == "I")) |
|
|
240 error ("VARTYPE must contain only C or I"); |
|
|
241 return; |
|
|
242 endif |
|
|
243 else |
|
|
244 ## As default we consider continuous vars |
|
|
245 vartype = repmat ("C", nx, 1); |
|
|
246 endif |
|
5232
|
247 |
|
5233
|
248 ## 9) Parameters vector |
|
|
249 |
|
|
250 if (nargin > 8) |
|
|
251 if (! isstruct (param)) |
|
|
252 error ("PARAM must be a structure"); |
|
|
253 return; |
|
|
254 endif |
|
|
255 else |
|
|
256 param = struct (); |
|
|
257 endif |
|
5232
|
258 |
|
5233
|
259 ## 10) Select solver method: simplex or interior point |
|
|
260 |
|
|
261 if (nargin > 9) |
|
|
262 if (! (isreal (lpsolver) && isscalar (lpsolver))) |
|
|
263 error ("LPSOLVER must be a real scalar value"); |
|
|
264 return; |
|
|
265 endif |
|
|
266 else |
|
|
267 lpsolver = 1; |
|
|
268 endif |
|
|
269 |
|
|
270 ## 11) Save the problem |
|
5232
|
271 |
|
5233
|
272 if (nargin > 10) |
|
|
273 if (! (isreal (savepb) && isscalar (lpsolver))) |
|
|
274 error ("LPSOLVER must be a real scalar value"); |
|
|
275 return; |
|
|
276 endif |
|
|
277 else |
|
|
278 savepb = 0; |
|
|
279 end |
|
5232
|
280 |
|
5233
|
281 [xopt, fmin, status, extra] = ... |
|
|
282 __glpk__ (sense, c, a, b, ctype, lb, ub, vartype, param, lpsolver, savepb); |
|
5232
|
283 |
|
|
284 endfunction |
