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