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-2015 by Power System Engineering Research Center (PSERC) 0097 % by Ray Zimmerman, PSERC Cornell 0098 % 0099 % $Id: qps_ot.m 2661 2015-03-20 17:02:46Z ray $ 0100 % 0101 % This file is part of MATPOWER. 0102 % Covered by the 3-clause BSD License (see LICENSE file for details). 0103 % See http://www.pserc.cornell.edu/matpower/ for more info. 0104 0105 %% check for Optimization Toolbox 0106 % if ~have_fcn('quadprog') 0107 % error('qps_ot: requires the Optimization Toolbox'); 0108 % end 0109 0110 %%----- input argument handling ----- 0111 %% gather inputs 0112 if nargin == 1 && isstruct(H) %% problem struct 0113 p = H; 0114 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0115 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0116 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0117 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0118 if isfield(p, 'u'), u = p.u; else, u = []; end 0119 if isfield(p, 'l'), l = p.l; else, l = []; end 0120 if isfield(p, 'A'), A = p.A; else, A = []; end 0121 if isfield(p, 'c'), c = p.c; else, c = []; end 0122 if isfield(p, 'H'), H = p.H; else, H = []; end 0123 else %% individual args 0124 if nargin < 9 0125 opt = []; 0126 if nargin < 8 0127 x0 = []; 0128 if nargin < 7 0129 xmax = []; 0130 if nargin < 6 0131 xmin = []; 0132 end 0133 end 0134 end 0135 end 0136 end 0137 0138 %% define nx, set default values for missing optional inputs 0139 if isempty(H) || ~any(any(H)) 0140 if isempty(A) && isempty(xmin) && isempty(xmax) 0141 error('qps_ot: LP problem must include constraints or variable bounds'); 0142 else 0143 if ~isempty(A) 0144 nx = size(A, 2); 0145 elseif ~isempty(xmin) 0146 nx = length(xmin); 0147 else % if ~isempty(xmax) 0148 nx = length(xmax); 0149 end 0150 end 0151 else 0152 nx = size(H, 1); 0153 end 0154 if isempty(c) 0155 c = zeros(nx, 1); 0156 end 0157 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0158 (isempty(u) || all(u == Inf))) 0159 A = sparse(0,nx); %% no limits => no linear constraints 0160 end 0161 nA = size(A, 1); %% number of original linear constraints 0162 if isempty(u) %% By default, linear inequalities are ... 0163 u = Inf(nA, 1); %% ... unbounded above and ... 0164 end 0165 if isempty(l) 0166 l = -Inf(nA, 1); %% ... unbounded below. 0167 end 0168 if isempty(xmin) %% By default, optimization variables are ... 0169 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0170 end 0171 if isempty(xmax) 0172 xmax = Inf(nx, 1); %% ... unbounded above. 0173 end 0174 if isempty(x0) 0175 x0 = zeros(nx, 1); 0176 end 0177 if isempty(H) || ~any(any(H)) 0178 isLP = 1; %% it's an LP 0179 else 0180 isLP = 0; %% nope, it's a QP 0181 end 0182 0183 %% default options 0184 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0185 verbose = opt.verbose; 0186 else 0187 verbose = 0; 0188 end 0189 0190 %% split up linear constraints 0191 ieq = find( abs(u-l) <= eps ); %% equality 0192 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0193 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0194 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0195 Ae = A(ieq, :); 0196 be = u(ieq); 0197 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0198 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0199 0200 %% grab some dimensions 0201 nlt = length(ilt); %% number of upper bounded linear inequalities 0202 ngt = length(igt); %% number of lower bounded linear inequalities 0203 nbx = length(ibx); %% number of doubly bounded linear inequalities 0204 0205 %% set up options 0206 if verbose > 1 0207 vrb = 'iter'; %% seems to be same as 'final' 0208 elseif verbose == 1 0209 vrb = 'final'; 0210 else 0211 vrb = 'off'; 0212 end 0213 if have_fcn('optimoptions') %% Optimization Tbx 6.3 + (R2013a +) 0214 %% could use optimset for everything, except some options are not 0215 %% recognized by optimset, only optimoptions, such as 0216 %% ot_opt.Algorithm = 'dual-simplex' 0217 if isLP 0218 ot_opt = optimoptions('linprog'); 0219 if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt) 0220 ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt); 0221 end 0222 else 0223 ot_opt = optimoptions('quadprog'); 0224 if have_fcn('quadprog_ls') 0225 ot_opt.Algorithm = 'interior-point-convex'; 0226 else 0227 ot_opt.LargeScale = 'off'; 0228 end 0229 if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt) 0230 ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt); 0231 end 0232 end 0233 ot_opt = optimoptions(ot_opt, 'Display', vrb); 0234 else %% need to use optimset() 0235 if isLP 0236 ot_opt = optimset('linprog'); 0237 if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt) 0238 ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt); 0239 end 0240 else 0241 ot_opt = optimset('quadprog'); 0242 if have_fcn('quadprog_ls') 0243 ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex'); 0244 else 0245 ot_opt = optimset(ot_opt, 'LargeScale', 'off'); 0246 end 0247 if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt) 0248 ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt); 0249 end 0250 end 0251 ot_opt = optimset(ot_opt, 'Display', vrb); 0252 end 0253 0254 %% call the solver 0255 if isLP 0256 [x, f, eflag, output, lam] = ... 0257 linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0258 else 0259 [x, f, eflag, output, lam] = ... 0260 quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0261 end 0262 0263 %% repackage lambdas 0264 if isempty(x) 0265 x = NaN(nx, 1); 0266 end 0267 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ... 0268 isempty(lam.lower) && isempty(lam.upper)) 0269 lambda = struct( ... 0270 'mu_l', NaN(nA, 1), ... 0271 'mu_u', NaN(nA, 1), ... 0272 'lower', NaN(nx, 1), ... 0273 'upper', NaN(nx, 1) ... 0274 ); 0275 else 0276 kl = find(lam.eqlin < 0); %% lower bound binding 0277 ku = find(lam.eqlin > 0); %% upper bound binding 0278 0279 mu_l = zeros(nA, 1); 0280 mu_l(ieq(kl)) = -lam.eqlin(kl); 0281 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0282 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0283 0284 mu_u = zeros(nA, 1); 0285 mu_u(ieq(ku)) = lam.eqlin(ku); 0286 mu_u(ilt) = lam.ineqlin(1:nlt); 0287 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0288 0289 lambda = struct( ... 0290 'mu_l', mu_l, ... 0291 'mu_u', mu_u, ... 0292 'lower', lam.lower(1:nx), ... 0293 'upper', lam.upper(1:nx) ... 0294 ); 0295 end