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