MIQPS_MATPOWER Mixed Integer Quadratic Program Solver for MATPOWER.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
MIQPS_MATPOWER(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
A common wrapper function for various QP solvers.
Solves 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),
'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer).
(MOSEK, GLPK, OT allow only 'C', 'B', or 'I')
OPT : optional options structure with the following fields,
all of which are also optional (default values shown in
parentheses)
alg ('DEFAULT') : determines which solver to use, can be either
a string (new-style) or a numerical alg code (old-style)
'DEFAULT' : (or 0) automatic, first available of CPLEX,
Gurobi, MOSEK, Opt Tbx (MILPs only), GLPK (MILPs only)
'CPLEX' : (or 500) CPLEX
'GLPK' : GLPK, (MILP problems only, i.e. empty H matrix)
'GUROBI' : (or 700) Gurobi
'MOSEK' : (or 600) MOSEK
'OT' : (or 300) Optimization Toolbox, INTLINPROG
(MILP problems only, i.e. empty H matrix)
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
cplex_opt - options struct for CPLEX
glpk_opt - options struct for GLPK
grb_opt - options struct for GBUROBI_MEX
intlinprog_opt - options struct for INTLINPROG
linprog_opt - options struct for LINPROG
mosek_opt - options struct for MOSEK
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 = converged
0 or negative values = algorithm specific failure codes
OUTPUT : output struct with the following fields:
alg - algorithm code of solver used
(others) - algorithm specific fields
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_matpower(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_matpower(H, c, A, l, u)
x = miqps_matpower(H, c, A, l, u, xmin, xmax)
x = miqps_matpower(H, c, A, l, u, xmin, xmax, x0)
x = miqps_matpower(H, c, A, l, u, xmin, xmax, x0, vtype)
x = miqps_matpower(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_matpower(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_matpower(...)
[x, f] = miqps_matpower(...)
[x, f, exitflag] = miqps_matpower(...)
[x, f, exitflag, output] = miqps_matpower(...)
[x, f, exitflag, output, lambda] = miqps_matpower(...)
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_matpower(H, c, A, l, u, xmin, [], x0, vtype, opt);0001 function [x, f, eflag, output, lambda] = miqps_matpower(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0002 %MIQPS_MATPOWER Mixed Integer Quadratic Program Solver for MATPOWER. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % MIQPS_MATPOWER(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT) 0005 % A common wrapper function for various QP solvers. 0006 % Solves the following QP (quadratic programming) problem: 0007 % 0008 % min 1/2 X'*H*X + C'*X 0009 % X 0010 % 0011 % subject to 0012 % 0013 % L <= A*X <= U (linear constraints) 0014 % XMIN <= X <= XMAX (variable bounds) 0015 % 0016 % Inputs (all optional except H, C, A and L): 0017 % H : matrix (possibly sparse) of quadratic cost coefficients 0018 % C : vector of linear cost coefficients 0019 % A, L, U : define the optional linear constraints. Default 0020 % values for the elements of L and U are -Inf and Inf, 0021 % respectively. 0022 % XMIN, XMAX : optional lower and upper bounds on the 0023 % X variables, defaults are -Inf and Inf, respectively. 0024 % X0 : optional starting value of optimization vector X 0025 % VTYPE : character string of length NX (number of elements in X), 0026 % or 1 (value applies to all variables in x), 0027 % allowed values are 'C' (continuous), 'B' (binary), 0028 % 'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer). 0029 % (MOSEK, GLPK, OT allow only 'C', 'B', or 'I') 0030 % OPT : optional options structure with the following fields, 0031 % all of which are also optional (default values shown in 0032 % parentheses) 0033 % alg ('DEFAULT') : determines which solver to use, can be either 0034 % a string (new-style) or a numerical alg code (old-style) 0035 % 'DEFAULT' : (or 0) automatic, first available of CPLEX, 0036 % Gurobi, MOSEK, Opt Tbx (MILPs only), GLPK (MILPs only) 0037 % 'CPLEX' : (or 500) CPLEX 0038 % 'GLPK' : GLPK, (MILP problems only, i.e. empty H matrix) 0039 % 'GUROBI' : (or 700) Gurobi 0040 % 'MOSEK' : (or 600) MOSEK 0041 % 'OT' : (or 300) Optimization Toolbox, INTLINPROG 0042 % (MILP problems only, i.e. empty H matrix) 0043 % verbose (0) - controls level of progress output displayed 0044 % 0 = no progress output 0045 % 1 = some progress output 0046 % 2 = verbose progress output 0047 % skip_prices (0) - flag that specifies whether or not to 0048 % skip the price computation stage, in which the problem 0049 % is re-solved for only the continuous variables, with all 0050 % others being constrained to their solved values 0051 % price_stage_warn_tol (1e-7) - tolerance on the objective fcn 0052 % value and primal variable relative match required to avoid 0053 % mis-match warning message 0054 % cplex_opt - options struct for CPLEX 0055 % glpk_opt - options struct for GLPK 0056 % grb_opt - options struct for GBUROBI_MEX 0057 % intlinprog_opt - options struct for INTLINPROG 0058 % linprog_opt - options struct for LINPROG 0059 % mosek_opt - options struct for MOSEK 0060 % PROBLEM : The inputs can alternatively be supplied in a single 0061 % PROBLEM struct with fields corresponding to the input arguments 0062 % described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt 0063 % 0064 % Outputs: 0065 % X : solution vector 0066 % F : final objective function value 0067 % EXITFLAG : exit flag 0068 % 1 = converged 0069 % 0 or negative values = algorithm specific failure codes 0070 % OUTPUT : output struct with the following fields: 0071 % alg - algorithm code of solver used 0072 % (others) - algorithm specific fields 0073 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0074 % multipliers on the constraints, with fields: 0075 % mu_l - lower (left-hand) limit on linear constraints 0076 % mu_u - upper (right-hand) limit on linear constraints 0077 % lower - lower bound on optimization variables 0078 % upper - upper bound on optimization variables 0079 % 0080 % Note the calling syntax is almost identical to that of QUADPROG 0081 % from MathWorks' Optimization Toolbox. The main difference is that 0082 % the linear constraints are specified with A, L, U instead of 0083 % A, B, Aeq, Beq. 0084 % 0085 % Calling syntax options: 0086 % [x, f, exitflag, output, lambda] = ... 0087 % miqps_matpower(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0088 % 0089 % x = miqps_matpower(H, c, A, l, u) 0090 % x = miqps_matpower(H, c, A, l, u, xmin, xmax) 0091 % x = miqps_matpower(H, c, A, l, u, xmin, xmax, x0) 0092 % x = miqps_matpower(H, c, A, l, u, xmin, xmax, x0, vtype) 0093 % x = miqps_matpower(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0094 % x = miqps_matpower(problem), where problem is a struct with fields: 0095 % H, c, A, l, u, xmin, xmax, x0, vtype, opt 0096 % all fields except 'c', 'A' and 'l' or 'u' are optional 0097 % x = miqps_matpower(...) 0098 % [x, f] = miqps_matpower(...) 0099 % [x, f, exitflag] = miqps_matpower(...) 0100 % [x, f, exitflag, output] = miqps_matpower(...) 0101 % [x, f, exitflag, output, lambda] = miqps_matpower(...) 0102 % 0103 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0104 % H = [ 1003.1 4.3 6.3 5.9; 0105 % 4.3 2.2 2.1 3.9; 0106 % 6.3 2.1 3.5 4.8; 0107 % 5.9 3.9 4.8 10 ]; 0108 % c = zeros(4,1); 0109 % A = [ 1 1 1 1; 0110 % 0.17 0.11 0.10 0.18 ]; 0111 % l = [1; 0.10]; 0112 % u = [1; Inf]; 0113 % xmin = zeros(4,1); 0114 % x0 = [1; 0; 0; 1]; 0115 % opt = struct('verbose', 2); 0116 % [x, f, s, out, lambda] = miqps_matpower(H, c, A, l, u, xmin, [], x0, vtype, opt); 0117 0118 % MATPOWER 0119 % Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC) 0120 % by Ray Zimmerman, PSERC Cornell 0121 % 0122 % This file is part of MATPOWER. 0123 % Covered by the 3-clause BSD License (see LICENSE file for details). 0124 % See https://matpower.org for more info. 0125 0126 %%----- input argument handling ----- 0127 %% gather inputs 0128 if nargin == 1 && isstruct(H) %% problem struct 0129 p = H; 0130 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0131 if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end 0132 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0133 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0134 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0135 if isfield(p, 'u'), u = p.u; else, u = []; end 0136 if isfield(p, 'l'), l = p.l; else, l = []; end 0137 if isfield(p, 'A'), A = p.A; else, A = []; end 0138 if isfield(p, 'c'), c = p.c; else, c = []; end 0139 if isfield(p, 'H'), H = p.H; else, H = []; end 0140 else %% individual args 0141 if nargin < 10 0142 opt = []; 0143 if nargin < 9 0144 vtype = []; 0145 if nargin < 8 0146 x0 = []; 0147 if nargin < 7 0148 xmax = []; 0149 if nargin < 6 0150 xmin = []; 0151 end 0152 end 0153 end 0154 end 0155 end 0156 end 0157 0158 %% default options 0159 if ~isempty(opt) && isfield(opt, 'alg') && ~isempty(opt.alg) 0160 alg = opt.alg; 0161 %% convert integer codes to string values 0162 if ~ischar(alg) 0163 switch alg 0164 case 0 0165 alg = 'DEFAULT'; 0166 case 300 0167 alg = 'OT'; 0168 case 500 0169 alg = 'CPLEX'; 0170 case 600 0171 alg = 'MOSEK'; 0172 case 700 0173 alg = 'GUROBI'; 0174 otherwise 0175 error('miqps_matpower: %d is not a valid algorithm code', alg); 0176 end 0177 end 0178 else 0179 alg = 'DEFAULT'; 0180 end 0181 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0182 verbose = opt.verbose; 0183 else 0184 verbose = 0; 0185 end 0186 if strcmp(alg, 'DEFAULT') 0187 if have_fcn('gurobi') %% use Gurobi by default, if available 0188 alg = 'GUROBI'; 0189 elseif have_fcn('cplex') %% if not, then CPLEX, if available 0190 alg = 'CPLEX'; 0191 elseif have_fcn('mosek') %% if not, then MOSEK, if available 0192 alg = 'MOSEK'; 0193 elseif isempty(H) || ~any(any(H)) %% if not, and linear objective 0194 if have_fcn('intlinprog') %% then Optimization Tbx, if available 0195 alg = 'OT'; 0196 elseif have_fcn('glpk') %% if not, and then GLPK, if available 0197 alg = 'GLPK'; 0198 end 0199 else 0200 error('miqps_matpower: no solvers available - requires CPLEX, Gurobi, MOSEK, INTLINPROG or GLPK'); 0201 end 0202 end 0203 0204 %%----- call the appropriate solver ----- 0205 switch alg 0206 case 'CPLEX' 0207 [x, f, eflag, output, lambda] = ... 0208 miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt); 0209 case 'GLPK' 0210 [x, f, eflag, output, lambda] = ... 0211 miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt); 0212 case 'GUROBI' 0213 [x, f, eflag, output, lambda] = ... 0214 miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt); 0215 case 'MOSEK' 0216 [x, f, eflag, output, lambda] = ... 0217 miqps_mosek(H, c, A, l, u, xmin, xmax, x0, vtype, opt); 0218 case 'OT' 0219 [x, f, eflag, output, lambda] = ... 0220 miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt); 0221 otherwise 0222 error('miqps_matpower: ''%s'' is not a valid algorithm code', alg); 0223 end 0224 if ~isfield(output, 'alg') || isempty(output.alg) 0225 output.alg = alg; 0226 end