QPS_GLPK Quadratic Program Solver based on QUADPROG/LINPROG.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, OPT)
A wrapper function providing a MATPOWER 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)
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
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, 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] = ...
qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_glpk(H, c, A, l, u)
x = qps_glpk(H, c, A, l, u, xmin, xmax)
x = qps_glpk(H, c, A, l, u, xmin, xmax, x0)
x = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_glpk(problem), where problem is a struct with fields:
H, c, A, l, u, xmin, xmax, x0, opt
all fields except 'c', 'A' and 'l' or 'u' are optional
x = qps_glpk(...)
[x, f] = qps_glpk(...)
[x, f, exitflag] = qps_glpk(...)
[x, f, exitflag, output] = qps_glpk(...)
[x, f, exitflag, output, lambda] = qps_glpk(...)
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] = qps_glpk(H, c, A, l, u, xmin, [], x0, opt);
See also GLPK.0001 function [x, f, eflag, output, lambda] = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_GLPK Quadratic Program Solver based on QUADPROG/LINPROG. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % A wrapper function providing a MATPOWER standardized interface for using 0006 % GLKP to solve the following LP (linear programming) problem: 0007 % 0008 % min 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 : dummy matrix (possibly sparse) of quadratic cost coefficients 0018 % for QP problems, which GLPK does not handle 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 (NOT USED) 0026 % OPT : optional options structure with the following fields, 0027 % all of which are also optional (default values shown in 0028 % parentheses) 0029 % verbose (0) - controls level of progress output displayed 0030 % 0 = no progress output 0031 % 1 = some progress output 0032 % 2 = verbose progress output 0033 % glpk_opt - options struct for GLPK, value in 0034 % verbose overrides these options 0035 % PROBLEM : The inputs can alternatively be supplied in a single 0036 % PROBLEM struct with fields corresponding to the input arguments 0037 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0038 % 0039 % Outputs: 0040 % X : solution vector 0041 % F : final objective function value 0042 % EXITFLAG : exit flag, 1 - optimal, <= 0 - infeasible, unbounded or other 0043 % OUTPUT : struct with errnum and status fields from GLPK output args 0044 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0045 % multipliers on the constraints, with fields: 0046 % mu_l - lower (left-hand) limit on linear constraints 0047 % mu_u - upper (right-hand) limit on linear constraints 0048 % lower - lower bound on optimization variables 0049 % upper - upper bound on optimization variables 0050 % 0051 % Note the calling syntax is almost identical to that of GLPK. The main 0052 % difference is that the linear constraints are specified with A, L, U 0053 % instead of A, B, Aeq, Beq. 0054 % 0055 % Calling syntax options: 0056 % [x, f, exitflag, output, lambda] = ... 0057 % qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt) 0058 % 0059 % x = qps_glpk(H, c, A, l, u) 0060 % x = qps_glpk(H, c, A, l, u, xmin, xmax) 0061 % x = qps_glpk(H, c, A, l, u, xmin, xmax, x0) 0062 % x = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt) 0063 % x = qps_glpk(problem), where problem is a struct with fields: 0064 % H, c, A, l, u, xmin, xmax, x0, opt 0065 % all fields except 'c', 'A' and 'l' or 'u' are optional 0066 % x = qps_glpk(...) 0067 % [x, f] = qps_glpk(...) 0068 % [x, f, exitflag] = qps_glpk(...) 0069 % [x, f, exitflag, output] = qps_glpk(...) 0070 % [x, f, exitflag, output, lambda] = qps_glpk(...) 0071 % 0072 % 0073 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0074 % H = [ 1003.1 4.3 6.3 5.9; 0075 % 4.3 2.2 2.1 3.9; 0076 % 6.3 2.1 3.5 4.8; 0077 % 5.9 3.9 4.8 10 ]; 0078 % c = zeros(4,1); 0079 % A = [ 1 1 1 1; 0080 % 0.17 0.11 0.10 0.18 ]; 0081 % l = [1; 0.10]; 0082 % u = [1; Inf]; 0083 % xmin = zeros(4,1); 0084 % x0 = [1; 0; 0; 1]; 0085 % opt = struct('verbose', 2); 0086 % [x, f, s, out, lambda] = qps_glpk(H, c, A, l, u, xmin, [], x0, opt); 0087 % 0088 % See also GLPK. 0089 0090 % MATPOWER 0091 % $Id: qps_glpk.m 2497 2014-12-17 19:56:20Z ray $ 0092 % by Ray Zimmerman, PSERC Cornell 0093 % Copyright (c) 2010-2014 by Power System Engineering Research Center (PSERC) 0094 % 0095 % This file is part of MATPOWER. 0096 % See http://www.pserc.cornell.edu/matpower/ for more info. 0097 % 0098 % MATPOWER is free software: you can redistribute it and/or modify 0099 % it under the terms of the GNU General Public License as published 0100 % by the Free Software Foundation, either version 3 of the License, 0101 % or (at your option) any later version. 0102 % 0103 % MATPOWER is distributed in the hope that it will be useful, 0104 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0105 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0106 % GNU General Public License for more details. 0107 % 0108 % You should have received a copy of the GNU General Public License 0109 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0110 % 0111 % Additional permission under GNU GPL version 3 section 7 0112 % 0113 % If you modify MATPOWER, or any covered work, to interface with 0114 % other modules (such as MATLAB code and MEX-files) available in a 0115 % MATLAB(R) or comparable environment containing parts covered 0116 % under other licensing terms, the licensors of MATPOWER grant 0117 % you additional permission to convey the resulting work. 0118 0119 %% check for Optimization Toolbox 0120 % if ~have_fcn('quadprog') 0121 % error('qps_glpk: requires the MEX interface to GLPK'); 0122 % end 0123 0124 %%----- input argument handling ----- 0125 %% gather inputs 0126 if nargin == 1 && isstruct(H) %% problem struct 0127 p = H; 0128 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0129 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0130 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0131 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0132 if isfield(p, 'u'), u = p.u; else, u = []; end 0133 if isfield(p, 'l'), l = p.l; else, l = []; end 0134 if isfield(p, 'A'), A = p.A; else, A = []; end 0135 if isfield(p, 'c'), c = p.c; else, c = []; end 0136 if isfield(p, 'H'), H = p.H; else, H = []; end 0137 else %% individual args 0138 if nargin < 9 0139 opt = []; 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 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('qps_glpk: 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 error('qps_glpk: GLPK handles only LP problems, not QP problems'); 0167 nx = size(H, 1); 0168 end 0169 if isempty(c) 0170 c = zeros(nx, 1); 0171 end 0172 if ~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 * ones(nA, 1); %% ... unbounded above and ... 0179 end 0180 if isempty(l) 0181 l = -Inf * ones(nA, 1); %% ... unbounded below. 0182 end 0183 if isempty(xmin) %% By default, optimization variables are ... 0184 xmin = -Inf * ones(nx, 1); %% ... unbounded below and ... 0185 end 0186 if isempty(xmax) 0187 xmax = Inf * ones(nx, 1); %% ... unbounded above. 0188 end 0189 if isempty(x0) 0190 x0 = zeros(nx, 1); 0191 end 0192 0193 %% default options 0194 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0195 verbose = opt.verbose; 0196 else 0197 verbose = 0; 0198 end 0199 0200 %% split up linear constraints 0201 ieq = find( abs(u-l) <= eps ); %% equality 0202 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0203 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0204 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0205 AA = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0206 bb = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx)]; 0207 0208 %% grab some dimensions 0209 nlt = length(ilt); %% number of upper bounded linear inequalities 0210 ngt = length(igt); %% number of lower bounded linear inequalities 0211 nbx = length(ibx); %% number of doubly bounded linear inequalities 0212 neq = length(ieq); %% number of equalities 0213 nie = nlt+ngt+2*nbx; %% number of inequalities 0214 0215 ctype = [repmat('S', neq, 1); repmat('U', nlt+ngt+2*nbx, 1)]; 0216 vtype = repmat('C', nx, 1); 0217 0218 %% set options struct for GLPK 0219 if ~isempty(opt) && isfield(opt, 'glpk_opt') && ~isempty(opt.glpk_opt) 0220 glpk_opt = glpk_options(opt.glpk_opt); 0221 else 0222 glpk_opt = glpk_options; 0223 end 0224 glpk_opt.msglev = verbose; 0225 0226 %% call the solver 0227 [x, f, errnum, extra] = ... 0228 glpk(c, AA, bb, xmin, xmax, ctype, vtype, 1, glpk_opt); 0229 0230 %% set exit flag 0231 if isfield(extra, 'status') 0232 output.errnum = errnum; 0233 output.status = extra.status; 0234 eflag = -errnum; 0235 if eflag == 0 && extra.status == 5 0236 eflag = 1; 0237 end 0238 else 0239 output.errnum = []; 0240 output.status = errnum; 0241 if have_fcn('octave') 0242 if errnum == 180 || errnum == 151 0243 eflag = 1; 0244 else 0245 eflag = 0; 0246 end 0247 else 0248 if errnum == 5 0249 eflag = 1; 0250 else 0251 eflag = 0; 0252 end 0253 end 0254 end 0255 0256 %% repackage lambdas 0257 if isempty(extra) || ~isfield(extra, 'lambda') || isempty(extra.lambda) 0258 lambda = struct( ... 0259 'mu_l', zeros(nA, 1), ... 0260 'mu_u', zeros(nA, 1), ... 0261 'lower', zeros(nx, 1), ... 0262 'upper', zeros(nx, 1) ... 0263 ); 0264 else 0265 lam.eqlin = extra.lambda(1:neq); 0266 lam.ineqlin = extra.lambda(neq+(1:nie)); 0267 lam.lower = extra.redcosts; 0268 lam.upper = -extra.redcosts; 0269 lam.lower(lam.lower < 0) = 0; 0270 lam.upper(lam.upper < 0) = 0; 0271 kl = find(lam.eqlin > 0); %% lower bound binding 0272 ku = find(lam.eqlin < 0); %% upper bound binding 0273 0274 mu_l = zeros(nA, 1); 0275 mu_l(ieq(kl)) = lam.eqlin(kl); 0276 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0277 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0278 0279 mu_u = zeros(nA, 1); 0280 mu_u(ieq(ku)) = -lam.eqlin(ku); 0281 mu_u(ilt) = -lam.ineqlin(1:nlt); 0282 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0283 0284 lambda = struct( ... 0285 'mu_l', mu_l, ... 0286 'mu_u', mu_u, ... 0287 'lower', lam.lower(1:nx), ... 0288 'upper', lam.upper(1:nx) ... 0289 ); 0290 end