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
max_it (0) - maximum number of iterations allowed
0 = use algorithm default
ot_opt - options struct for QUADPROG/LINPROG, values in
verbose and max_it override 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 % max_it (0) - maximum number of iterations allowed 0034 % 0 = use algorithm default 0035 % ot_opt - options struct for QUADPROG/LINPROG, values in 0036 % verbose and max_it override 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,v 1.13 2011/09/09 15:26:09 cvs Exp $ 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 0197 %% default options 0198 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0199 verbose = opt.verbose; 0200 else 0201 verbose = 0; 0202 end 0203 if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it) 0204 max_it = opt.max_it; 0205 else 0206 max_it = 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 ~isempty(opt) && isfield(opt, 'ot_opt') && ~isempty(opt.ot_opt) 0226 ot_opt = opt.ot_opt; 0227 else 0228 if isempty(H) || ~any(any(H)) 0229 ot_opt = optimset('linprog'); 0230 else 0231 ot_opt = optimset('quadprog'); 0232 if have_fcn('quadprog_ls') 0233 ot_opt = optimset(ot_opt, 'Algorithm', 'interior-point-convex'); 0234 else 0235 ot_opt = optimset(ot_opt, 'LargeScale', 'off'); 0236 end 0237 end 0238 end 0239 if max_it 0240 ot_opt = optimset(ot_opt, 'MaxIter', max_it); 0241 end 0242 if verbose > 1 0243 ot_opt = optimset(ot_opt, 'Display', 'iter'); %% seems to be same as 'final' 0244 elseif verbose == 1 0245 ot_opt = optimset(ot_opt, 'Display', 'final'); 0246 else 0247 ot_opt = optimset(ot_opt, 'Display', 'off'); 0248 end 0249 0250 %% call the solver 0251 if isempty(H) || ~any(any(H)) 0252 [x, f, eflag, output, lam] = ... 0253 linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0254 else 0255 [x, f, eflag, output, lam] = ... 0256 quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt); 0257 end 0258 0259 %% repackage lambdas 0260 kl = find(lam.eqlin < 0); %% lower bound binding 0261 ku = find(lam.eqlin > 0); %% upper bound binding 0262 0263 mu_l = zeros(nA, 1); 0264 mu_l(ieq(kl)) = -lam.eqlin(kl); 0265 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0266 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0267 0268 mu_u = zeros(nA, 1); 0269 mu_u(ieq(ku)) = lam.eqlin(ku); 0270 mu_u(ilt) = lam.ineqlin(1:nlt); 0271 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0272 0273 lambda = struct( ... 0274 'mu_l', mu_l, ... 0275 'mu_u', mu_u, ... 0276 'lower', lam.lower(1:nx), ... 0277 'upper', lam.upper(1:nx) ... 0278 );