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 % $Id: qps_ot.m 2385 2014-10-11 15:33:45Z ray $ 0097 % by Ray Zimmerman, PSERC Cornell 0098 % Copyright (c) 2010-2011 by Power System Engineering Research Center (PSERC) 0099 % 0100 % This file is part of MATPOWER. 0101 % See http://www.pserc.cornell.edu/matpower/ for more info. 0102 % 0103 % MATPOWER is free software: you can redistribute it and/or modify 0104 % it under the terms of the GNU General Public License as published 0105 % by the Free Software Foundation, either version 3 of the License, 0106 % or (at your option) any later version. 0107 % 0108 % MATPOWER is distributed in the hope that it will be useful, 0109 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0110 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0111 % GNU General Public License for more details. 0112 % 0113 % You should have received a copy of the GNU General Public License 0114 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0115 % 0116 % Additional permission under GNU GPL version 3 section 7 0117 % 0118 % If you modify MATPOWER, or any covered work, to interface with 0119 % other modules (such as MATLAB code and MEX-files) available in a 0120 % MATLAB(R) or comparable environment containing parts covered 0121 % under other licensing terms, the licensors of MATPOWER grant 0122 % you additional permission to convey the resulting work. 0123 0124 %% check for Optimization Toolbox 0125 % if ~have_fcn('quadprog') 0126 % error('qps_ot: requires the Optimization Toolbox'); 0127 % end 0128 0129 %%----- input argument handling ----- 0130 %% gather inputs 0131 if nargin == 1 && isstruct(H) %% problem struct 0132 p = H; 0133 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0134 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0135 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0136 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0137 if isfield(p, 'u'), u = p.u; else, u = []; end 0138 if isfield(p, 'l'), l = p.l; else, l = []; end 0139 if isfield(p, 'A'), A = p.A; else, A = []; end 0140 if isfield(p, 'c'), c = p.c; else, c = []; end 0141 if isfield(p, 'H'), H = p.H; else, H = []; end 0142 else %% individual args 0143 if nargin < 9 0144 opt = []; 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 0157 %% define nx, set default values for missing optional inputs 0158 if isempty(H) || ~any(any(H)) 0159 if isempty(A) && isempty(xmin) && isempty(xmax) 0160 error('qps_ot: LP problem must include constraints or variable bounds'); 0161 else 0162 if ~isempty(A) 0163 nx = size(A, 2); 0164 elseif ~isempty(xmin) 0165 nx = length(xmin); 0166 else % if ~isempty(xmax) 0167 nx = length(xmax); 0168 end 0169 end 0170 else 0171 nx = size(H, 1); 0172 end 0173 if isempty(c) 0174 c = zeros(nx, 1); 0175 end 0176 if ~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0177 (isempty(u) || all(u == Inf)) 0178 A = sparse(0,nx); %% no limits => no linear constraints 0179 end 0180 nA = size(A, 1); %% number of original linear constraints 0181 if isempty(u) %% By default, linear inequalities are ... 0182 u = Inf * ones(nA, 1); %% ... unbounded above and ... 0183 end 0184 if isempty(l) 0185 l = -Inf * ones(nA, 1); %% ... unbounded below. 0186 end 0187 if isempty(xmin) %% By default, optimization variables are ... 0188 xmin = -Inf * ones(nx, 1); %% ... unbounded below and ... 0189 end 0190 if isempty(xmax) 0191 xmax = Inf * ones(nx, 1); %% ... unbounded above. 0192 end 0193 if isempty(x0) 0194 x0 = zeros(nx, 1); 0195 end 0196 if isempty(H) || ~any(any(H)) 0197 isLP = 1; %% it's an LP 0198 else 0199 isLP = 0; %% nope, it's a QP 0200 end 0201 0202 %% default options 0203 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0204 verbose = opt.verbose; 0205 else 0206 verbose = 0; 0207 end 0208 0209 %% split up linear constraints 0210 ieq = find( abs(u-l) <= eps ); %% equality 0211 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0212 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0213 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0214 Ae = A(ieq, :); 0215 be = u(ieq); 0216 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0217 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0218 0219 %% grab some dimensions 0220 nlt = length(ilt); %% number of upper bounded linear inequalities 0221 ngt = length(igt); %% number of lower bounded linear inequalities 0222 nbx = length(ibx); %% number of doubly bounded linear inequalities 0223 0224 %% set up options 0225 if verbose > 1 0226 vrb = 'iter'; %% seems to be same as 'final' 0227 elseif verbose == 1 0228 vrb = 'final'; 0229 else 0230 vrb = 'off'; 0231 end 0232 if have_fcn('optimoptions') %% Optimization Tbx 6.3 + (R2013a +) 0233 %% could use optimset for everything, except some options are not 0234 %% recognized by optimset, only optimoptions, such as 0235 %% ot_opt.Algorithm = 'dual-simplex' 0236 if isLP 0237 ot_opt = optimoptions('linprog'); 0238 if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt) 0239 ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt); 0240 end 0241 else 0242 ot_opt = optimoptions('quadprog'); 0243 if have_fcn('quadprog_ls') 0244 ot_opt.Algorithm = 'interior-point-convex'; 0245 else 0246 ot_opt.LargeScale = 'off'; 0247 end 0248 if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt) 0249 ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt); 0250 end 0251 end 0252 ot_opt = optimoptions(ot_opt, 'Display', vrb); 0253 else %% need to use optimset() 0254 if isLP 0255 ot_opt = optimset('linprog'); 0256 if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt) 0257 ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt); 0258 end 0259 else 0260 ot_opt = optimset('quadprog'); 0261 if have_fcn('quadprog_ls') 0262 ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex'); 0263 else 0264 ot_opt = optimset(ot_opt, 'LargeScale', 'off'); 0265 end 0266 if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt) 0267 ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt); 0268 end 0269 end 0270 ot_opt = optimset(ot_opt, 'Display', vrb); 0271 end 0272 0273 %% call the solver 0274 if isLP 0275 [x, f, eflag, output, lam] = ... 0276 linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 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 mu_l(ieq(kl)) = -lam.eqlin(kl); 0300 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0301 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0302 0303 mu_u = zeros(nA, 1); 0304 mu_u(ieq(ku)) = lam.eqlin(ku); 0305 mu_u(ilt) = lam.ineqlin(1:nlt); 0306 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0307 0308 lambda = struct( ... 0309 'mu_l', mu_l, ... 0310 'mu_u', mu_u, ... 0311 'lower', lam.lower(1:nx), ... 0312 'upper', lam.upper(1:nx) ... 0313 ); 0314 end