QPS_OT Quadratic Program Solver based on QUADPROG/LINPROG.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_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
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
linprog_opt - options struct for LINPROG, value in
verbose overrides these options
quadprog_opt - options struct for QUADPROG, 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 : 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] = ...
qps_ot(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_ot(H, c, A, l, u)
x = qps_ot(H, c, A, l, u, xmin, xmax)
x = qps_ot(H, c, A, l, u, xmin, xmax, x0)
x = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_ot(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_ot(...)
[x, f] = qps_ot(...)
[x, f, exitflag] = qps_ot(...)
[x, f, exitflag, output] = qps_ot(...)
[x, f, exitflag, output, lambda] = qps_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] = qps_ot(H, c, A, l, u, xmin, [], x0, opt);
See also QUADPROG, LINPROG.0001 function [x, f, eflag, output, lambda] = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_OT Quadratic Program Solver based on QUADPROG/LINPROG. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_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 % 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 % linprog_opt - options struct for LINPROG, value in 0034 % verbose overrides these options 0035 % quadprog_opt - options struct for QUADPROG, value in 0036 % verbose overrides these options 0037 % PROBLEM : The inputs can alternatively be supplied in a single 0038 % PROBLEM struct with fields corresponding to the input arguments 0039 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0040 % 0041 % Outputs: 0042 % X : solution vector 0043 % F : final objective function value 0044 % EXITFLAG : QUADPROG/LINPROG exit flag 0045 % (see QUADPROG and LINPROG documentation for details) 0046 % OUTPUT : QUADPROG/LINPROG output struct 0047 % (see QUADPROG and LINPROG documentation for details) 0048 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0049 % multipliers on the constraints, with fields: 0050 % mu_l - lower (left-hand) limit on linear constraints 0051 % mu_u - upper (right-hand) limit on linear constraints 0052 % lower - lower bound on optimization variables 0053 % upper - upper bound on optimization variables 0054 % 0055 % Note the calling syntax is almost identical to that of QUADPROG 0056 % from MathWorks' Optimization Toolbox. The main difference is that 0057 % the linear constraints are specified with A, L, U instead of 0058 % A, B, Aeq, Beq. 0059 % 0060 % Calling syntax options: 0061 % [x, f, exitflag, output, lambda] = ... 0062 % qps_ot(H, c, A, l, u, xmin, xmax, x0, opt) 0063 % 0064 % x = qps_ot(H, c, A, l, u) 0065 % x = qps_ot(H, c, A, l, u, xmin, xmax) 0066 % x = qps_ot(H, c, A, l, u, xmin, xmax, x0) 0067 % x = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt) 0068 % x = qps_ot(problem), where problem is a struct with fields: 0069 % H, c, A, l, u, xmin, xmax, x0, opt 0070 % all fields except 'c', 'A' and 'l' or 'u' are optional 0071 % x = qps_ot(...) 0072 % [x, f] = qps_ot(...) 0073 % [x, f, exitflag] = qps_ot(...) 0074 % [x, f, exitflag, output] = qps_ot(...) 0075 % [x, f, exitflag, output, lambda] = qps_ot(...) 0076 % 0077 % 0078 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0079 % H = [ 1003.1 4.3 6.3 5.9; 0080 % 4.3 2.2 2.1 3.9; 0081 % 6.3 2.1 3.5 4.8; 0082 % 5.9 3.9 4.8 10 ]; 0083 % c = zeros(4,1); 0084 % A = [ 1 1 1 1; 0085 % 0.17 0.11 0.10 0.18 ]; 0086 % l = [1; 0.10]; 0087 % u = [1; Inf]; 0088 % xmin = zeros(4,1); 0089 % x0 = [1; 0; 0; 1]; 0090 % opt = struct('verbose', 2); 0091 % [x, f, s, out, lambda] = qps_ot(H, c, A, l, u, xmin, [], x0, opt); 0092 % 0093 % See also QUADPROG, LINPROG. 0094 0095 % MATPOWER 0096 % Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC) 0097 % by Ray Zimmerman, PSERC Cornell 0098 % 0099 % This file is part of MATPOWER. 0100 % Covered by the 3-clause BSD License (see LICENSE file for details). 0101 % See http://www.pserc.cornell.edu/matpower/ for more info. 0102 0103 %% check for Optimization Toolbox 0104 % if ~have_fcn('quadprog') 0105 % error('qps_ot: requires the Optimization Toolbox'); 0106 % end 0107 0108 %%----- input argument handling ----- 0109 %% gather inputs 0110 if nargin == 1 && isstruct(H) %% problem struct 0111 p = H; 0112 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0113 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0114 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0115 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0116 if isfield(p, 'u'), u = p.u; else, u = []; end 0117 if isfield(p, 'l'), l = p.l; else, l = []; end 0118 if isfield(p, 'A'), A = p.A; else, A = []; end 0119 if isfield(p, 'c'), c = p.c; else, c = []; end 0120 if isfield(p, 'H'), H = p.H; else, H = []; end 0121 else %% individual args 0122 if nargin < 9 0123 opt = []; 0124 if nargin < 8 0125 x0 = []; 0126 if nargin < 7 0127 xmax = []; 0128 if nargin < 6 0129 xmin = []; 0130 end 0131 end 0132 end 0133 end 0134 end 0135 0136 %% define nx, set default values for missing optional inputs 0137 if isempty(H) || ~any(any(H)) 0138 if isempty(A) && isempty(xmin) && isempty(xmax) 0139 error('qps_ot: LP problem must include constraints or variable bounds'); 0140 else 0141 if ~isempty(A) 0142 nx = size(A, 2); 0143 elseif ~isempty(xmin) 0144 nx = length(xmin); 0145 else % if ~isempty(xmax) 0146 nx = length(xmax); 0147 end 0148 end 0149 else 0150 nx = size(H, 1); 0151 end 0152 if isempty(c) 0153 c = zeros(nx, 1); 0154 end 0155 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0156 (isempty(u) || all(u == Inf))) 0157 A = sparse(0,nx); %% no limits => no linear constraints 0158 end 0159 nA = size(A, 1); %% number of original linear constraints 0160 if isempty(u) %% By default, linear inequalities are ... 0161 u = Inf(nA, 1); %% ... unbounded above and ... 0162 end 0163 if isempty(l) 0164 l = -Inf(nA, 1); %% ... unbounded below. 0165 end 0166 if isempty(xmin) %% By default, optimization variables are ... 0167 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0168 end 0169 if isempty(xmax) 0170 xmax = Inf(nx, 1); %% ... unbounded above. 0171 end 0172 if isempty(x0) 0173 x0 = zeros(nx, 1); 0174 end 0175 if isempty(H) || ~any(any(H)) 0176 isLP = 1; %% it's an LP 0177 else 0178 isLP = 0; %% nope, it's a QP 0179 end 0180 0181 %% default options 0182 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0183 verbose = opt.verbose; 0184 else 0185 verbose = 0; 0186 end 0187 %% Matlab or Octave 0188 matlab = have_fcn('matlab'); 0189 otver = have_fcn('quadprog', 'vnum'); 0190 0191 %% split up linear constraints 0192 ieq = find( abs(u-l) <= eps ); %% equality 0193 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0194 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0195 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0196 Ae = A(ieq, :); 0197 be = u(ieq); 0198 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0199 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0200 0201 %% grab some dimensions 0202 nlt = length(ilt); %% number of upper bounded linear inequalities 0203 ngt = length(igt); %% number of lower bounded linear inequalities 0204 nbx = length(ibx); %% number of doubly bounded linear inequalities 0205 0206 %% set up options 0207 if verbose > 1 0208 vrb = 'iter'; %% seems to be same as 'final' 0209 elseif verbose == 1 0210 vrb = 'final'; 0211 else 0212 vrb = 'off'; 0213 end 0214 if have_fcn('optimoptions') %% Optimization Tbx 6.3 + (R2013a +) 0215 %% could use optimset for everything, except some options are not 0216 %% recognized by optimset, only optimoptions, such as 0217 %% ot_opt.Algorithm = 'dual-simplex' 0218 if isLP 0219 ot_opt = optimoptions('linprog'); 0220 if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt) 0221 ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt); 0222 end 0223 else 0224 ot_opt = optimoptions('quadprog'); 0225 if have_fcn('quadprog_ls') 0226 ot_opt.Algorithm = 'interior-point-convex'; 0227 else 0228 ot_opt.LargeScale = 'off'; 0229 end 0230 if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt) 0231 ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt); 0232 end 0233 end 0234 ot_opt = optimoptions(ot_opt, 'Display', vrb); 0235 else %% need to use optimset() 0236 if isLP 0237 if matlab 0238 ot_opt = optimset('linprog'); 0239 else 0240 ot_opt = optimset(); 0241 end 0242 if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt) 0243 ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt); 0244 end 0245 else 0246 if matlab 0247 ot_opt = optimset('quadprog'); 0248 if have_fcn('quadprog_ls') 0249 ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex'); 0250 else 0251 ot_opt = optimset(ot_opt, 'LargeScale', 'off'); 0252 end 0253 else 0254 ot_opt = optimset(); 0255 end 0256 if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt) 0257 ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt); 0258 end 0259 end 0260 ot_opt = optimset(ot_opt, 'Display', vrb); 0261 end 0262 0263 %% call the solver 0264 if isLP 0265 if matlab 0266 [x, f, eflag, output, lam] = ... 0267 linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0268 else 0269 % don't use linprog under Octave (using GLPK directly is recommended) 0270 % [x, f] = linprog(c, Ai, bi, Ae, be, xmin, xmax); 0271 % eflag = []; 0272 % output = []; 0273 % lam = []; 0274 [x, f, eflag, output, lam] = ... 0275 quadprog(sparse(nx,nx), c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0276 end 0277 else 0278 [x, f, eflag, output, lam] = ... 0279 quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0280 end 0281 0282 %% repackage lambdas 0283 if isempty(x) 0284 x = NaN(nx, 1); 0285 end 0286 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ... 0287 isempty(lam.lower) && isempty(lam.upper)) 0288 lambda = struct( ... 0289 'mu_l', NaN(nA, 1), ... 0290 'mu_u', NaN(nA, 1), ... 0291 'lower', NaN(nx, 1), ... 0292 'upper', NaN(nx, 1) ... 0293 ); 0294 else 0295 kl = find(lam.eqlin < 0); %% lower bound binding 0296 ku = find(lam.eqlin > 0); %% upper bound binding 0297 0298 mu_l = zeros(nA, 1); 0299 % %% workaround for Octave optim 1.5.0 and earlier, which 0300 % %% has opposite sign convention for equality multipliers 0301 % if ~matlab && otver <= 1.005 0302 % mu_l(ieq(ku)) = lam.eqlin(ku); 0303 % else 0304 mu_l(ieq(kl)) = -lam.eqlin(kl); 0305 % end 0306 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0307 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0308 0309 mu_u = zeros(nA, 1); 0310 % %% workaround for Octave optim 1.5.0 and earlier, which 0311 % %% has opposite sign convention for equality multipliers 0312 % if ~matlab && otver <= 1.005 0313 % mu_u(ieq(kl)) = -lam.eqlin(kl); 0314 % else 0315 mu_u(ieq(ku)) = lam.eqlin(ku); 0316 % end 0317 mu_u(ilt) = lam.ineqlin(1:nlt); 0318 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0319 0320 %% workaround for Octave optim 1.5.0 and earlier, which 0321 %% has opposite sign convention for equality multipliers 0322 % if ~matlab && otver <= 1.005 0323 %% there are also issues with variable bounds that are 0324 %% converted to equalities, and maybe other issues 0325 % end 0326 0327 lambda = struct( ... 0328 'mu_l', mu_l, ... 0329 'mu_u', mu_u, ... 0330 'lower', lam.lower(1:nx), ... 0331 'upper', lam.upper(1:nx) ... 0332 ); 0333 end