MIQPS_OT Mixed Integer Linear Program Solver based on INTLINPROG.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
MIQPS_OT(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_OT(PROBLEM)
A wrapper function providing a standardized interface for using
QUADPROG or LINPROG from the Optimization Toolbox to solve the
following QP (quadratic programming) problem:
min 1/2 X'*H*X + 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 : matrix (possibly sparse) of quadratic cost coefficients
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
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
intlinprog_opt - options struct for INTLINPROG, value in verbose
overrides these options
linprog_opt - options struct for LINPROG, 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 : QUADPROG/LINPROG exit flag
(see QUADPROG and LINPROG documentation for details)
OUTPUT : QUADPROG/LINPROG output struct
(see QUADPROG and LINPROG documentation for details)
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 QUADPROG
from MathWorks' Optimization Toolbox. 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_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_ot(H, c, A, l, u)
x = miqps_ot(H, c, A, l, u, xmin, xmax)
x = miqps_ot(H, c, A, l, u, xmin, xmax, x0)
x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype)
x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_ot(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_ot(...)
[x, f] = miqps_ot(...)
[x, f, exitflag] = miqps_ot(...)
[x, f, exitflag, output] = miqps_ot(...)
[x, f, exitflag, output, lambda] = miqps_ot(...)
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_ot(H, c, A, l, u, xmin, [], x0, vtype, opt);
See also MIQPS_MASTER, INTLINPROG, QUADPROG, LINPROG.0001 function [x, f, eflag, output, lambda] = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0002 %MIQPS_OT Mixed Integer Linear Program Solver based on INTLINPROG. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % MIQPS_OT(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_OT(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % QUADPROG or LINPROG from the Optimization Toolbox to solve the 0008 % following QP (quadratic programming) problem: 0009 % 0010 % min 1/2 X'*H*X + C'*X 0011 % X 0012 % 0013 % subject to 0014 % 0015 % L <= A*X <= U (linear constraints) 0016 % XMIN <= X <= XMAX (variable bounds) 0017 % 0018 % Inputs (all optional except H, C, A and L): 0019 % H : matrix (possibly sparse) of quadratic cost coefficients 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 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 % intlinprog_opt - options struct for INTLINPROG, value in verbose 0046 % overrides these options 0047 % linprog_opt - options struct for LINPROG, value in verbose 0048 % overrides these options 0049 % PROBLEM : The inputs can alternatively be supplied in a single 0050 % PROBLEM struct with fields corresponding to the input arguments 0051 % described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt 0052 % 0053 % Outputs: 0054 % X : solution vector 0055 % F : final objective function value 0056 % EXITFLAG : QUADPROG/LINPROG exit flag 0057 % (see QUADPROG and LINPROG documentation for details) 0058 % OUTPUT : QUADPROG/LINPROG output struct 0059 % (see QUADPROG and LINPROG documentation for details) 0060 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0061 % multipliers on the constraints, with fields: 0062 % mu_l - lower (left-hand) limit on linear constraints 0063 % mu_u - upper (right-hand) limit on linear constraints 0064 % lower - lower bound on optimization variables 0065 % upper - upper bound on optimization variables 0066 % 0067 % Note the calling syntax is almost identical to that of QUADPROG 0068 % from MathWorks' Optimization Toolbox. The main difference is that 0069 % the linear constraints are specified with A, L, U instead of 0070 % A, B, Aeq, Beq. 0071 % 0072 % Calling syntax options: 0073 % [x, f, exitflag, output, lambda] = ... 0074 % miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0075 % 0076 % x = miqps_ot(H, c, A, l, u) 0077 % x = miqps_ot(H, c, A, l, u, xmin, xmax) 0078 % x = miqps_ot(H, c, A, l, u, xmin, xmax, x0) 0079 % x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype) 0080 % x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0081 % x = miqps_ot(problem), where problem is a struct with fields: 0082 % H, c, A, l, u, xmin, xmax, x0, vtype, opt 0083 % all fields except 'c', 'A' and 'l' or 'u' are optional 0084 % x = miqps_ot(...) 0085 % [x, f] = miqps_ot(...) 0086 % [x, f, exitflag] = miqps_ot(...) 0087 % [x, f, exitflag, output] = miqps_ot(...) 0088 % [x, f, exitflag, output, lambda] = miqps_ot(...) 0089 % 0090 % 0091 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0092 % H = [ 1003.1 4.3 6.3 5.9; 0093 % 4.3 2.2 2.1 3.9; 0094 % 6.3 2.1 3.5 4.8; 0095 % 5.9 3.9 4.8 10 ]; 0096 % c = zeros(4,1); 0097 % A = [ 1 1 1 1; 0098 % 0.17 0.11 0.10 0.18 ]; 0099 % l = [1; 0.10]; 0100 % u = [1; Inf]; 0101 % xmin = zeros(4,1); 0102 % x0 = [1; 0; 0; 1]; 0103 % opt = struct('verbose', 2); 0104 % [x, f, s, out, lambda] = miqps_ot(H, c, A, l, u, xmin, [], x0, vtype, opt); 0105 % 0106 % See also MIQPS_MASTER, INTLINPROG, QUADPROG, LINPROG. 0107 0108 % MP-Opt-Model 0109 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0110 % by Ray Zimmerman, PSERC Cornell 0111 % 0112 % This file is part of MP-Opt-Model. 0113 % Covered by the 3-clause BSD License (see LICENSE file for details). 0114 % See https://github.com/MATPOWER/mp-opt-model for more info. 0115 0116 %% check for Optimization Toolbox 0117 % if ~have_feature('quadprog') 0118 % error('miqps_ot: requires the Optimization Toolbox'); 0119 % end 0120 0121 %%----- input argument handling ----- 0122 %% gather inputs 0123 if nargin == 1 && isstruct(H) %% problem struct 0124 p = H; 0125 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0126 if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end 0127 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0128 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0129 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0130 if isfield(p, 'u'), u = p.u; else, u = []; end 0131 if isfield(p, 'l'), l = p.l; else, l = []; end 0132 if isfield(p, 'A'), A = p.A; else, A = []; end 0133 if isfield(p, 'c'), c = p.c; else, c = []; end 0134 if isfield(p, 'H'), H = p.H; else, H = []; end 0135 else %% individual args 0136 if nargin < 10 0137 opt = []; 0138 if nargin < 9 0139 vtype = []; 0140 if nargin < 8 0141 x0 = []; 0142 if nargin < 7 0143 xmax = []; 0144 if nargin < 6 0145 xmin = []; 0146 end 0147 end 0148 end 0149 end 0150 end 0151 end 0152 0153 %% define nx, set default values for missing optional inputs 0154 if isempty(H) || ~any(any(H)) 0155 if isempty(A) && isempty(xmin) && isempty(xmax) 0156 error('miqps_ot: LP problem must include constraints or variable bounds'); 0157 else 0158 if ~isempty(A) 0159 nx = size(A, 2); 0160 elseif ~isempty(xmin) 0161 nx = length(xmin); 0162 else % if ~isempty(xmax) 0163 nx = length(xmax); 0164 end 0165 end 0166 else 0167 nx = size(H, 1); 0168 end 0169 if isempty(c) 0170 c = zeros(nx, 1); 0171 end 0172 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0173 (isempty(u) || all(u == Inf))) 0174 A = sparse(0,nx); %% no limits => no linear constraints 0175 end 0176 nA = size(A, 1); %% number of original linear constraints 0177 if isempty(u) %% By default, linear inequalities are ... 0178 u = Inf(nA, 1); %% ... unbounded above and ... 0179 end 0180 if isempty(l) 0181 l = -Inf(nA, 1); %% ... unbounded below. 0182 end 0183 if isempty(xmin) %% By default, optimization variables are ... 0184 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0185 end 0186 if isempty(xmax) 0187 xmax = Inf(nx, 1); %% ... unbounded above. 0188 end 0189 if isempty(x0) 0190 x0 = zeros(nx, 1); 0191 end 0192 if isempty(H) || ~any(any(H)) 0193 isLP = 1; %% it's an LP 0194 else 0195 isLP = 0; %% nope, it's a QP 0196 end 0197 0198 %% default options 0199 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0200 verbose = opt.verbose; 0201 else 0202 verbose = 0; 0203 end 0204 0205 %% split up linear constraints 0206 ieq = find( abs(u-l) <= eps ); %% equality 0207 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0208 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0209 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0210 Ae = A(ieq, :); 0211 be = u(ieq); 0212 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0213 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0214 0215 %% grab some dimensions 0216 nlt = length(ilt); %% number of upper bounded linear inequalities 0217 ngt = length(igt); %% number of lower bounded linear inequalities 0218 nbx = length(ibx); %% number of doubly bounded linear inequalities 0219 0220 %% mixed integer? 0221 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I')) 0222 mi = 0; 0223 vtype = repmat('C', 1, nx); 0224 else 0225 mi = 1; 0226 %% expand vtype to nx elements if necessary 0227 if length(vtype) == 1 && nx > 1 0228 vtype = char(vtype * ones(1, nx)); 0229 end 0230 end 0231 %% check bounds for 'B' variables to make sure they are between 0 and 1 0232 k = find(vtype == 'B'); 0233 if ~isempty(k) 0234 kk = find(xmax(k) > 1); 0235 xmax(k(kk)) = 1; 0236 kk = find(xmin(k) < 0); 0237 xmin(k(kk)) = 0; 0238 end 0239 intcon = find(vtype == 'B' | vtype == 'I'); 0240 0241 %% set up options 0242 if verbose > 1 0243 vrb = 'iter'; %% seems to be same as 'final' 0244 elseif verbose == 1 0245 vrb = 'final'; 0246 else 0247 vrb = 'off'; 0248 end 0249 0250 if isLP 0251 if mi 0252 ot_opt = optimoptions('intlinprog'); 0253 if ~isempty(opt) && isfield(opt, 'intlinprog_opt') && ~isempty(opt.intlinprog_opt) 0254 ot_opt = nested_struct_copy(ot_opt, opt.intlinprog_opt); 0255 end 0256 else 0257 ot_opt = optimoptions('linprog'); 0258 if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt) 0259 ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt); 0260 end 0261 end 0262 else 0263 if mi 0264 error('miqps_ot: MIQP problems not supported.'); 0265 end 0266 ot_opt = optimoptions('quadprog'); 0267 if have_feature('quadprog_ls') 0268 ot_opt.Algorithm = 'interior-point-convex'; 0269 else 0270 ot_opt.LargeScale = 'off'; 0271 end 0272 if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt) 0273 ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt); 0274 end 0275 end 0276 ot_opt = optimoptions(ot_opt, 'Display', vrb); 0277 0278 %% call the solver 0279 if isLP 0280 if mi 0281 [x, f, eflag, output] = ... 0282 intlinprog(c, intcon, Ai, bi, Ae, be, xmin, xmax, ot_opt); 0283 lam = []; 0284 else 0285 [x, f, eflag, output, lam] = ... 0286 linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0287 end 0288 else 0289 [x, f, eflag, output, lam] = ... 0290 quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0291 end 0292 0293 %% repackage lambdas 0294 if isempty(x) 0295 x = NaN(nx, 1); 0296 end 0297 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ... 0298 isempty(lam.lower) && isempty(lam.upper)) 0299 lambda = struct( ... 0300 'mu_l', NaN(nA, 1), ... 0301 'mu_u', NaN(nA, 1), ... 0302 'lower', NaN(nx, 1), ... 0303 'upper', NaN(nx, 1) ... 0304 ); 0305 else 0306 kl = find(lam.eqlin < 0); %% lower bound binding 0307 ku = find(lam.eqlin > 0); %% upper bound binding 0308 0309 mu_l = zeros(nA, 1); 0310 mu_l(ieq(kl)) = -lam.eqlin(kl); 0311 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0312 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0313 0314 mu_u = zeros(nA, 1); 0315 mu_u(ieq(ku)) = lam.eqlin(ku); 0316 mu_u(ilt) = lam.ineqlin(1:nlt); 0317 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0318 0319 lambda = struct( ... 0320 'mu_l', mu_l, ... 0321 'mu_u', mu_u, ... 0322 'lower', lam.lower(1:nx), ... 0323 'upper', lam.upper(1:nx) ... 0324 ); 0325 end 0326 0327 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices) 0328 if verbose 0329 fprintf('--- Integer stage complete, starting price computation stage ---\n'); 0330 end 0331 if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol) 0332 tol = opt.price_stage_warn_tol; 0333 else 0334 tol = 1e-7; 0335 end 0336 k = intcon; 0337 x(k) = round(x(k)); 0338 xmin(k) = x(k); 0339 xmax(k) = x(k); 0340 x0 = x; 0341 % opt.linprog_opt.Algorithm = 'dual-simplex'; %% dual-simplex 0342 0343 [x_, f_, eflag_, output_, lambda] = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt); 0344 % output 0345 % output.message 0346 if eflag ~= eflag_ 0347 error('miqps_ot: EXITFLAG from price computation stage = %d', eflag_); 0348 end 0349 if abs(f - f_)/max(abs(f), 1) > tol 0350 warning('miqps_ot: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1)); 0351 end 0352 xn = x; 0353 xn(abs(xn)<1) = 1; 0354 [mx, k] = max(abs(x - x_) ./ xn); 0355 if mx > tol 0356 warning('miqps_ot: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k)); 0357 end 0358 output.price_stage = output_; 0359 % output_ 0360 % output_.message 0361 end