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