MIQPS_GLPK Mixed Integer Linear Program Solver based on GLPK - GNU Linear Programming Kit.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
MIQPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_GLPK(PROBLEM)
A wrapper function providing a standardized interface for using
GLKP to solve the following LP (linear programming) problem:
min C'*X
X
subject to
L <= A*X <= U (linear constraints)
XMIN <= X <= XMAX (variable bounds)
Inputs (all optional except H, C, A and L):
H : dummy matrix (possibly sparse) of quadratic cost coefficients
for QP problems, which GLPK does not handle
C : vector of linear cost coefficients
A, L, U : define the optional linear constraints. Default
values for the elements of L and U are -Inf and Inf,
respectively.
XMIN, XMAX : optional lower and upper bounds on the
X variables, defaults are -Inf and Inf, respectively.
X0 : optional starting value of optimization vector X (NOT USED)
VTYPE : character string of length NX (number of elements in X),
or 1 (value applies to all variables in x),
allowed values are 'C' (continuous), 'B' (binary) or
'I' (integer).
OPT : optional options structure with the following fields,
all of which are also optional (default values shown in
parentheses)
verbose (0) - controls level of progress output displayed
0 = no progress output
1 = some progress output
2 = verbose progress output
skip_prices (0) - flag that specifies whether or not to
skip the price computation stage, in which the problem
is re-solved for only the continuous variables, with all
others being constrained to their solved values
price_stage_warn_tol (1e-7) - tolerance on the objective fcn
value and primal variable relative match required to avoid
mis-match warning message
glpk_opt - options struct for GLPK, value in verbose
overrides these options
PROBLEM : The inputs can alternatively be supplied in a single
PROBLEM struct with fields corresponding to the input arguments
described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt
Outputs:
X : solution vector
F : final objective function value
EXITFLAG : exit flag, 1 - optimal, <= 0 - infeasible, unbounded or other
OUTPUT : struct with errnum and status fields from GLPK output args
LAMBDA : struct containing the Langrange and Kuhn-Tucker
multipliers on the constraints, with fields:
mu_l - lower (left-hand) limit on linear constraints
mu_u - upper (right-hand) limit on linear constraints
lower - lower bound on optimization variables
upper - upper bound on optimization variables
Note the calling syntax is almost identical to that of GLPK. The main
difference is that the linear constraints are specified with A, L, U
instead of A, B, Aeq, Beq.
Calling syntax options:
[x, f, exitflag, output, lambda] = ...
miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_glpk(H, c, A, l, u)
x = miqps_glpk(H, c, A, l, u, xmin, xmax)
x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0)
x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype)
x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_glpk(problem), where problem is a struct with fields:
H, c, A, l, u, xmin, xmax, x0, vtype, opt
all fields except 'c', 'A' and 'l' or 'u' are optional
x = miqps_glpk(...)
[x, f] = miqps_glpk(...)
[x, f, exitflag] = miqps_glpk(...)
[x, f, exitflag, output] = miqps_glpk(...)
[x, f, exitflag, output, lambda] = miqps_glpk(...)
Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
H = [ 1003.1 4.3 6.3 5.9;
4.3 2.2 2.1 3.9;
6.3 2.1 3.5 4.8;
5.9 3.9 4.8 10 ];
c = zeros(4,1);
A = [ 1 1 1 1;
0.17 0.11 0.10 0.18 ];
l = [1; 0.10];
u = [1; Inf];
xmin = zeros(4,1);
x0 = [1; 0; 0; 1];
opt = struct('verbose', 2);
[x, f, s, out, lambda] = miqps_glpk(H, c, A, l, u, xmin, [], x0, vtype, opt);
See also MIQPS_MASTER, GLPK.0001 function [x, f, eflag, output, lambda] = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0002 %MIQPS_GLPK Mixed Integer Linear Program Solver based on GLPK - GNU Linear Programming Kit. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % MIQPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_GLPK(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % GLKP to solve the following LP (linear programming) problem: 0008 % 0009 % min C'*X 0010 % X 0011 % 0012 % subject to 0013 % 0014 % L <= A*X <= U (linear constraints) 0015 % XMIN <= X <= XMAX (variable bounds) 0016 % 0017 % Inputs (all optional except H, C, A and L): 0018 % H : dummy matrix (possibly sparse) of quadratic cost coefficients 0019 % for QP problems, which GLPK does not handle 0020 % C : vector of linear cost coefficients 0021 % A, L, U : define the optional linear constraints. Default 0022 % values for the elements of L and U are -Inf and Inf, 0023 % respectively. 0024 % XMIN, XMAX : optional lower and upper bounds on the 0025 % X variables, defaults are -Inf and Inf, respectively. 0026 % X0 : optional starting value of optimization vector X (NOT USED) 0027 % VTYPE : character string of length NX (number of elements in X), 0028 % or 1 (value applies to all variables in x), 0029 % allowed values are 'C' (continuous), 'B' (binary) or 0030 % 'I' (integer). 0031 % OPT : optional options structure with the following fields, 0032 % all of which are also optional (default values shown in 0033 % parentheses) 0034 % verbose (0) - controls level of progress output displayed 0035 % 0 = no progress output 0036 % 1 = some progress output 0037 % 2 = verbose progress output 0038 % skip_prices (0) - flag that specifies whether or not to 0039 % skip the price computation stage, in which the problem 0040 % is re-solved for only the continuous variables, with all 0041 % others being constrained to their solved values 0042 % price_stage_warn_tol (1e-7) - tolerance on the objective fcn 0043 % value and primal variable relative match required to avoid 0044 % mis-match warning message 0045 % glpk_opt - options struct for GLPK, value in verbose 0046 % overrides these options 0047 % PROBLEM : The inputs can alternatively be supplied in a single 0048 % PROBLEM struct with fields corresponding to the input arguments 0049 % described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt 0050 % 0051 % Outputs: 0052 % X : solution vector 0053 % F : final objective function value 0054 % EXITFLAG : exit flag, 1 - optimal, <= 0 - infeasible, unbounded or other 0055 % OUTPUT : struct with errnum and status fields from GLPK output args 0056 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0057 % multipliers on the constraints, with fields: 0058 % mu_l - lower (left-hand) limit on linear constraints 0059 % mu_u - upper (right-hand) limit on linear constraints 0060 % lower - lower bound on optimization variables 0061 % upper - upper bound on optimization variables 0062 % 0063 % Note the calling syntax is almost identical to that of GLPK. The main 0064 % difference is that the linear constraints are specified with A, L, U 0065 % instead of A, B, Aeq, Beq. 0066 % 0067 % Calling syntax options: 0068 % [x, f, exitflag, output, lambda] = ... 0069 % miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0070 % 0071 % x = miqps_glpk(H, c, A, l, u) 0072 % x = miqps_glpk(H, c, A, l, u, xmin, xmax) 0073 % x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0) 0074 % x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype) 0075 % x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0076 % x = miqps_glpk(problem), where problem is a struct with fields: 0077 % H, c, A, l, u, xmin, xmax, x0, vtype, opt 0078 % all fields except 'c', 'A' and 'l' or 'u' are optional 0079 % x = miqps_glpk(...) 0080 % [x, f] = miqps_glpk(...) 0081 % [x, f, exitflag] = miqps_glpk(...) 0082 % [x, f, exitflag, output] = miqps_glpk(...) 0083 % [x, f, exitflag, output, lambda] = miqps_glpk(...) 0084 % 0085 % 0086 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0087 % H = [ 1003.1 4.3 6.3 5.9; 0088 % 4.3 2.2 2.1 3.9; 0089 % 6.3 2.1 3.5 4.8; 0090 % 5.9 3.9 4.8 10 ]; 0091 % c = zeros(4,1); 0092 % A = [ 1 1 1 1; 0093 % 0.17 0.11 0.10 0.18 ]; 0094 % l = [1; 0.10]; 0095 % u = [1; Inf]; 0096 % xmin = zeros(4,1); 0097 % x0 = [1; 0; 0; 1]; 0098 % opt = struct('verbose', 2); 0099 % [x, f, s, out, lambda] = miqps_glpk(H, c, A, l, u, xmin, [], x0, vtype, opt); 0100 % 0101 % See also MIQPS_MASTER, GLPK. 0102 0103 % MP-Opt-Model 0104 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0105 % by Ray Zimmerman, PSERC Cornell 0106 % 0107 % This file is part of MP-Opt-Model. 0108 % Covered by the 3-clause BSD License (see LICENSE file for details). 0109 % See https://github.com/MATPOWER/mp-opt-model for more info. 0110 0111 %% check for Optimization Toolbox 0112 % if ~have_feature('quadprog') 0113 % error('miqps_glpk: requires the MEX interface to GLPK'); 0114 % end 0115 0116 %%----- input argument handling ----- 0117 %% gather inputs 0118 if nargin == 1 && isstruct(H) %% problem struct 0119 p = H; 0120 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0121 if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end 0122 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0123 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0124 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0125 if isfield(p, 'u'), u = p.u; else, u = []; end 0126 if isfield(p, 'l'), l = p.l; else, l = []; end 0127 if isfield(p, 'A'), A = p.A; else, A = []; end 0128 if isfield(p, 'c'), c = p.c; else, c = []; end 0129 if isfield(p, 'H'), H = p.H; else, H = []; end 0130 else %% individual args 0131 if nargin < 10 0132 opt = []; 0133 if nargin < 9 0134 vtype = []; 0135 if nargin < 8 0136 x0 = []; 0137 if nargin < 7 0138 xmax = []; 0139 if nargin < 6 0140 xmin = []; 0141 end 0142 end 0143 end 0144 end 0145 end 0146 end 0147 0148 %% define nx, set default values for missing optional inputs 0149 if isempty(H) || ~any(any(H)) 0150 if isempty(A) && isempty(xmin) && isempty(xmax) 0151 error('miqps_glpk: LP problem must include constraints or variable bounds'); 0152 else 0153 if ~isempty(A) 0154 nx = size(A, 2); 0155 elseif ~isempty(xmin) 0156 nx = length(xmin); 0157 else % if ~isempty(xmax) 0158 nx = length(xmax); 0159 end 0160 end 0161 else 0162 error('miqps_glpk: GLPK handles only LP problems, not QP problems'); 0163 nx = size(H, 1); 0164 end 0165 if isempty(c) 0166 c = zeros(nx, 1); 0167 end 0168 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0169 (isempty(u) || all(u == Inf))) 0170 A = sparse(0,nx); %% no limits => no linear constraints 0171 end 0172 nA = size(A, 1); %% number of original linear constraints 0173 if isempty(u) %% By default, linear inequalities are ... 0174 u = Inf(nA, 1); %% ... unbounded above and ... 0175 end 0176 if isempty(l) 0177 l = -Inf(nA, 1); %% ... unbounded below. 0178 end 0179 if isempty(xmin) %% By default, optimization variables are ... 0180 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0181 end 0182 if isempty(xmax) 0183 xmax = Inf(nx, 1); %% ... unbounded above. 0184 end 0185 if isempty(x0) 0186 x0 = zeros(nx, 1); 0187 end 0188 0189 %% default options 0190 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0191 verbose = opt.verbose; 0192 else 0193 verbose = 0; 0194 end 0195 0196 %% split up linear constraints 0197 ieq = find( abs(u-l) <= eps ); %% equality 0198 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0199 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0200 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0201 AA = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0202 bb = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx)]; 0203 0204 %% grab some dimensions 0205 nlt = length(ilt); %% number of upper bounded linear inequalities 0206 ngt = length(igt); %% number of lower bounded linear inequalities 0207 nbx = length(ibx); %% number of doubly bounded linear inequalities 0208 neq = length(ieq); %% number of equalities 0209 nie = nlt+ngt+2*nbx; %% number of inequalities 0210 0211 ctype = [repmat('S', neq, 1); repmat('U', nlt+ngt+2*nbx, 1)]; 0212 0213 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I')) 0214 mi = 0; 0215 vtype = repmat('C', nx, 1); 0216 else 0217 mi = 1; 0218 %% expand vtype to nx elements if necessary 0219 if length(vtype) == 1 && nx > 1 0220 vtype = char(vtype * ones(nx, 1)); 0221 elseif size(vtype, 2) > 1 %% make sure it's a col vector 0222 vtype = vtype'; 0223 end 0224 end 0225 %% convert 'B' variables to 'I' and clip bounds to [0, 1] 0226 k = find(vtype == 'B'); 0227 if ~isempty(k) 0228 kk = find(xmax(k) > 1); 0229 xmax(k(kk)) = 1; 0230 kk = find(xmin(k) < 0); 0231 xmin(k(kk)) = 0; 0232 vtype(k) = 'I'; 0233 end 0234 0235 %% set options struct for GLPK 0236 if ~isempty(opt) && isfield(opt, 'glpk_opt') && ~isempty(opt.glpk_opt) 0237 glpk_opt = glpk_options(opt.glpk_opt); 0238 else 0239 glpk_opt = glpk_options; 0240 end 0241 glpk_opt.msglev = verbose; 0242 0243 %% call the solver 0244 [x, f, errnum, extra] = ... 0245 glpk(c, AA, bb, xmin, xmax, ctype, vtype, 1, glpk_opt); 0246 0247 %% set exit flag 0248 if isfield(extra, 'status') %% status found in extra.status 0249 output.errnum = errnum; 0250 output.status = extra.status; 0251 eflag = -errnum; 0252 if eflag == 0 && extra.status == 5 0253 eflag = 1; 0254 end 0255 else %% status found in errnum 0256 output.errnum = []; 0257 output.status = errnum; 0258 if have_feature('octave') 0259 if errnum == 180 || errnum == 151 || errnum == 171 0260 eflag = 1; 0261 else 0262 eflag = 0; 0263 end 0264 else 0265 if errnum == 5 0266 eflag = 1; 0267 else 0268 eflag = 0; 0269 end 0270 end 0271 end 0272 0273 %% repackage lambdas 0274 if isempty(extra) || ~isfield(extra, 'lambda') || isempty(extra.lambda) 0275 lambda = struct( ... 0276 'mu_l', zeros(nA, 1), ... 0277 'mu_u', zeros(nA, 1), ... 0278 'lower', zeros(nx, 1), ... 0279 'upper', zeros(nx, 1) ... 0280 ); 0281 else 0282 lam.eqlin = extra.lambda(1:neq); 0283 lam.ineqlin = extra.lambda(neq+(1:nie)); 0284 lam.lower = extra.redcosts; 0285 lam.upper = -extra.redcosts; 0286 lam.lower(lam.lower < 0) = 0; 0287 lam.upper(lam.upper < 0) = 0; 0288 kl = find(lam.eqlin > 0); %% lower bound binding 0289 ku = find(lam.eqlin < 0); %% upper bound binding 0290 0291 mu_l = zeros(nA, 1); 0292 mu_l(ieq(kl)) = lam.eqlin(kl); 0293 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0294 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0295 0296 mu_u = zeros(nA, 1); 0297 mu_u(ieq(ku)) = -lam.eqlin(ku); 0298 mu_u(ilt) = -lam.ineqlin(1:nlt); 0299 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0300 0301 lambda = struct( ... 0302 'mu_l', mu_l, ... 0303 'mu_u', mu_u, ... 0304 'lower', lam.lower(1:nx), ... 0305 'upper', lam.upper(1:nx) ... 0306 ); 0307 end 0308 0309 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices) 0310 if verbose 0311 fprintf('--- Integer stage complete, starting price computation stage ---\n'); 0312 end 0313 if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol) 0314 tol = opt.price_stage_warn_tol; 0315 else 0316 tol = 1e-7; 0317 end 0318 k = find(vtype == 'I' | vtype == 'B'); 0319 x(k) = round(x(k)); 0320 xmin(k) = x(k); 0321 xmax(k) = x(k); 0322 x0 = x; 0323 opt.glpk_opt.lpsolver = 1; %% simplex 0324 opt.glpk_opt.dual = 0; %% primal simplex 0325 if have_feature('octave') && have_feature('octave', 'vnum') >= 3.007 0326 opt.glpk_opt.dual = 1; %% primal simplex 0327 end 0328 0329 [x_, f_, eflag_, output_, lambda] = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt); 0330 if eflag ~= eflag_ 0331 error('miqps_glpk: EXITFLAG from price computation stage = %d', eflag_); 0332 end 0333 if abs(f - f_)/max(abs(f), 1) > tol 0334 warning('miqps_glpk: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1)); 0335 end 0336 xn = x; 0337 xn(abs(xn)<1) = 1; 0338 [mx, k] = max(abs(x - x_) ./ xn); 0339 if mx > tol 0340 warning('miqps_glpk: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k)); 0341 end 0342 output.price_stage = output_; 0343 end