QPS_CLP Quadratic Program Solver based on CLP - COIN-OR Linear Programming.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_CLP(H, C, A, L, U, XMIN, XMAX, X0, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_CLP(PROBLEM)
A wrapper function providing a standardized interface for using
CLP 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 (NOT USED)
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
clp_opt - options struct for CLP, 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 : exit flag, 1 - optimal, -1 - infeasible, -2 - unbounded
-3 - max iterations/time exceeded
OUTPUT : struct with fields
exitflag - raw CLP exit flag: 0 - optimal, 1 - infeasible,
2 - unbounded, 3 - max iterations/time exceeded
status - string with explanation of exitflag
(iter - depending on build of solver this may contain
the number of iterations)
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 CLP. 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_clp(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_clp(H, c, A, l, u)
x = qps_clp(H, c, A, l, u, xmin, xmax)
x = qps_clp(H, c, A, l, u, xmin, xmax, x0)
x = qps_clp(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_clp(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_clp(...)
[x, f] = qps_clp(...)
[x, f, exitflag] = qps_clp(...)
[x, f, exitflag, output] = qps_clp(...)
[x, f, exitflag, output, lambda] = qps_clp(...)
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_clp(H, c, A, l, u, xmin, [], x0, opt);
See also QPS_MASTER, CLP.0001 function [x, f, eflag, output, lambda] = qps_clp(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_CLP Quadratic Program Solver based on CLP - COIN-OR Linear Programming. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_CLP(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_CLP(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % CLP to solve the 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 (NOT USED) 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 % clp_opt - options struct for CLP, value in verbose 0034 % overrides these options 0035 % PROBLEM : The inputs can alternatively be supplied in a single 0036 % PROBLEM struct with fields corresponding to the input arguments 0037 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0038 % 0039 % Outputs: 0040 % X : solution vector 0041 % F : final objective function value 0042 % EXITFLAG : exit flag, 1 - optimal, -1 - infeasible, -2 - unbounded 0043 % -3 - max iterations/time exceeded 0044 % OUTPUT : struct with fields 0045 % exitflag - raw CLP exit flag: 0 - optimal, 1 - infeasible, 0046 % 2 - unbounded, 3 - max iterations/time exceeded 0047 % status - string with explanation of exitflag 0048 % (iter - depending on build of solver this may contain 0049 % the number of iterations) 0050 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0051 % multipliers on the constraints, with fields: 0052 % mu_l - lower (left-hand) limit on linear constraints 0053 % mu_u - upper (right-hand) limit on linear constraints 0054 % lower - lower bound on optimization variables 0055 % upper - upper bound on optimization variables 0056 % 0057 % Note the calling syntax is almost identical to that of CLP. The main 0058 % difference is that the linear constraints are specified with A, L, U 0059 % instead of A, B, Aeq, Beq. 0060 % 0061 % Calling syntax options: 0062 % [x, f, exitflag, output, lambda] = ... 0063 % qps_clp(H, c, A, l, u, xmin, xmax, x0, opt) 0064 % 0065 % x = qps_clp(H, c, A, l, u) 0066 % x = qps_clp(H, c, A, l, u, xmin, xmax) 0067 % x = qps_clp(H, c, A, l, u, xmin, xmax, x0) 0068 % x = qps_clp(H, c, A, l, u, xmin, xmax, x0, opt) 0069 % x = qps_clp(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_clp(...) 0073 % [x, f] = qps_clp(...) 0074 % [x, f, exitflag] = qps_clp(...) 0075 % [x, f, exitflag, output] = qps_clp(...) 0076 % [x, f, exitflag, output, lambda] = qps_clp(...) 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_clp(H, c, A, l, u, xmin, [], x0, opt); 0093 % 0094 % See also QPS_MASTER, CLP. 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_clp: requires the MEX interface to CLP'); 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_clp: 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 ~issparse(A) 0177 A = sparse(A); 0178 end 0179 if ~issparse(H) 0180 H = sparse(H); 0181 end 0182 0183 0184 %% default options 0185 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0186 verbose = opt.verbose; 0187 else 0188 verbose = 0; 0189 end 0190 0191 %% set up options struct for CLP 0192 if ~isempty(opt) && isfield(opt, 'clp_opt') && ~isempty(opt.clp_opt) 0193 clp_opt = clp_options(opt.clp_opt); 0194 else 0195 clp_opt = clp_options; 0196 end 0197 0198 if have_feature('opti_clp') %% use OPTI Toolbox verision's MEX interface 0199 clp_opt.display = verbose; 0200 0201 [x, f, exitflag, iter, lam] = clp(tril(H), c, A, l, u, xmin, xmax, clp_opt); 0202 0203 output.iter = iter; 0204 0205 %% repackage lambdas 0206 % if isempty(x) 0207 % x = NaN(nx, 1); 0208 % f = NaN; 0209 % end 0210 % if isempty(lam) 0211 % lambda = struct( ... 0212 % 'mu_l', zeros(nA, 1), ... 0213 % 'mu_u', zeros(nA, 1), ... 0214 % 'lower', zeros(nx, 1), ... 0215 % 'upper', zeros(nx, 1) ... 0216 % ); 0217 % else 0218 mu_l = lam.dual_row; 0219 mu_u = -lam.dual_row; 0220 lower = lam.dual_col; 0221 upper = -lam.dual_col; 0222 mu_l(mu_l < 0) = 0; 0223 mu_u(mu_u < 0) = 0; 0224 lower(lower < 0) = 0; 0225 upper(upper < 0) = 0; 0226 0227 lambda = struct( ... 0228 'mu_l', mu_l, ... 0229 'mu_u', mu_u, ... 0230 'lower', lower, ... 0231 'upper', upper ... 0232 ); 0233 % end 0234 else 0235 clp_opt.verbose = verbose; 0236 0237 %% split up linear constraints 0238 ieq = find( abs(u-l) <= eps ); %% equality 0239 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0240 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0241 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0242 Ae = A(ieq, :); 0243 be = u(ieq); 0244 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0245 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0246 0247 %% grab some dimensions 0248 nlt = length(ilt); %% number of upper bounded linear inequalities 0249 ngt = length(igt); %% number of lower bounded linear inequalities 0250 nbx = length(ibx); %% number of doubly bounded linear inequalities 0251 0252 %% call the solver 0253 [x, z, exitflag] = ... 0254 clp(H, c, Ai, bi, Ae, be, xmin, xmax, clp_opt); 0255 0256 %% repackage lambdas 0257 if isempty(x) 0258 x = NaN(nx, 1); 0259 f = NaN; 0260 else 0261 if isempty(H) || ~any(any(H)) 0262 f = c'*x; 0263 else 0264 f = 0.5 * x'*H*x + c'*x; 0265 end 0266 end 0267 if isempty(z) 0268 lambda = struct( ... 0269 'mu_l', zeros(nA, 1), ... 0270 'mu_u', zeros(nA, 1), ... 0271 'lower', zeros(nx, 1), ... 0272 'upper', zeros(nx, 1) ... 0273 ); 0274 else 0275 neq = length(be); 0276 nie = length(bi); 0277 lam.eqlin = z(1:neq); 0278 lam.ineqlin = z(neq+(1:nie)); 0279 %%----- MEXCLP DOES NOT RETURN MULTIPLIERS ON VARIABLE BOUNDS :-/ ----- 0280 lam.lower = NaN(nx, 1); 0281 lam.upper = NaN(nx, 1); 0282 kl = find(lam.eqlin > 0); %% lower bound binding 0283 ku = find(lam.eqlin < 0); %% upper bound binding 0284 0285 mu_l = zeros(nA, 1); 0286 mu_l(ieq(kl)) = lam.eqlin(kl); 0287 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0288 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0289 0290 mu_u = zeros(nA, 1); 0291 mu_u(ieq(ku)) = -lam.eqlin(ku); 0292 mu_u(ilt) = -lam.ineqlin(1:nlt); 0293 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0294 0295 lambda = struct( ... 0296 'mu_l', mu_l, ... 0297 'mu_u', mu_u, ... 0298 'lower', lam.lower, ... 0299 'upper', lam.upper ... 0300 ); 0301 end 0302 end 0303 0304 %% set eflag 0305 eflag = -exitflag; 0306 if eflag == 0 %% success 0307 eflag = 1; 0308 end 0309 0310 %% set status 0311 output.exitflag = exitflag; 0312 switch exitflag 0313 case 0 0314 output.status = 'optimal'; 0315 case 1 0316 output.status = 'primal infeasible'; 0317 case 2 0318 output.status = 'dual infeasible'; 0319 case 3 0320 output.status = 'max iterations or time exceeded'; 0321 otherwise 0322 output.status = 'unknown exit code'; 0323 end