QPS_OSQP Quadratic Program Solver based on OSQP.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_OSQP(H, C, A, L, U, XMIN, XMAX, X0, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_OSQP(PROBLEM)
A wrapper function providing a standardized interface for using
OSQP 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
3 = even more verbose progress output
osqp_opt - options struct for OSQP (see
https://osqp.org/docs/interfaces/solver_settings.html),
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 : OSQP exit flag
1 = converged
0 or negative values OSQP status value
(see OSQP documentation for details)
OUTPUT : OSQP results struct
(see OSQP 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_osqp(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_osqp(H, c, A, l, u)
x = qps_osqp(H, c, A, l, u, xmin, xmax)
x = qps_osqp(H, c, A, l, u, xmin, xmax, x0)
x = qps_osqp(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_osqp(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_osqp(...)
[x, f] = qps_osqp(...)
[x, f, exitflag] = qps_osqp(...)
[x, f, exitflag, output] = qps_osqp(...)
[x, f, exitflag, output, lambda] = qps_osqp(...)
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_osqp(H, c, A, l, u, xmin, [], x0, opt);
See also QPS_MASTER, OSQP.0001 function [x, f, eflag, output, lambda] = qps_osqp(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_OSQP Quadratic Program Solver based on OSQP. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_OSQP(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_OSQP(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % OSQP 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 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 % 3 = even more verbose progress output 0034 % osqp_opt - options struct for OSQP (see 0035 % https://osqp.org/docs/interfaces/solver_settings.html), 0036 % value in 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 : OSQP exit flag 0045 % 1 = converged 0046 % 0 or negative values OSQP status value 0047 % (see OSQP documentation for details) 0048 % OUTPUT : OSQP results struct 0049 % (see OSQP documentation for details) 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 QUADPROG 0058 % from MathWorks' Optimization Toolbox. The main difference is that 0059 % the linear constraints are specified with A, L, U instead of 0060 % A, B, Aeq, Beq. 0061 % 0062 % Calling syntax options: 0063 % [x, f, exitflag, output, lambda] = ... 0064 % qps_osqp(H, c, A, l, u, xmin, xmax, x0, opt) 0065 % 0066 % x = qps_osqp(H, c, A, l, u) 0067 % x = qps_osqp(H, c, A, l, u, xmin, xmax) 0068 % x = qps_osqp(H, c, A, l, u, xmin, xmax, x0) 0069 % x = qps_osqp(H, c, A, l, u, xmin, xmax, x0, opt) 0070 % x = qps_osqp(problem), where problem is a struct with fields: 0071 % H, c, A, l, u, xmin, xmax, x0, opt 0072 % all fields except 'c', 'A' and 'l' or 'u' are optional 0073 % x = qps_osqp(...) 0074 % [x, f] = qps_osqp(...) 0075 % [x, f, exitflag] = qps_osqp(...) 0076 % [x, f, exitflag, output] = qps_osqp(...) 0077 % [x, f, exitflag, output, lambda] = qps_osqp(...) 0078 % 0079 % 0080 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0081 % H = [ 1003.1 4.3 6.3 5.9; 0082 % 4.3 2.2 2.1 3.9; 0083 % 6.3 2.1 3.5 4.8; 0084 % 5.9 3.9 4.8 10 ]; 0085 % c = zeros(4,1); 0086 % A = [ 1 1 1 1; 0087 % 0.17 0.11 0.10 0.18 ]; 0088 % l = [1; 0.10]; 0089 % u = [1; Inf]; 0090 % xmin = zeros(4,1); 0091 % x0 = [1; 0; 0; 1]; 0092 % opt = struct('verbose', 2); 0093 % [x, f, s, out, lambda] = qps_osqp(H, c, A, l, u, xmin, [], x0, opt); 0094 % 0095 % See also QPS_MASTER, OSQP. 0096 0097 % MP-Opt-Model 0098 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0099 % by Ray Zimmerman, PSERC Cornell 0100 % 0101 % This file is part of MP-Opt-Model. 0102 % Covered by the 3-clause BSD License (see LICENSE file for details). 0103 % See https://github.com/MATPOWER/mp-opt-model for more info. 0104 0105 %% check for OSQP 0106 % if ~have_feature('osqp') 0107 % error('qps_osqp: requires OSQP (https://osqp.org)'); 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 lpqp = 'LP'; 0141 if isempty(A) && isempty(xmin) && isempty(xmax) 0142 error('qps_osqp: LP problem must include constraints or variable bounds'); 0143 else 0144 if ~isempty(A) 0145 nx = size(A, 2); 0146 elseif ~isempty(xmin) 0147 nx = length(xmin); 0148 else % if ~isempty(xmax) 0149 nx = length(xmax); 0150 end 0151 end 0152 else 0153 lpqp = 'QP'; 0154 nx = size(H, 1); 0155 end 0156 if isempty(c) 0157 c = zeros(nx, 1); 0158 end 0159 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0160 (isempty(u) || all(u == Inf))) 0161 A = sparse(0,nx); %% no limits => no linear constraints 0162 end 0163 nA = size(A, 1); %% number of original linear constraints 0164 if isempty(u) %% By default, linear inequalities are ... 0165 u = Inf(nA, 1); %% ... unbounded above and ... 0166 end 0167 if isempty(l) 0168 l = -Inf(nA, 1); %% ... unbounded below. 0169 end 0170 0171 %% default options 0172 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0173 verbose = opt.verbose; 0174 else 0175 verbose = 0; 0176 end 0177 0178 %% set up options struct for OSQP 0179 osqp_opt = struct(); 0180 if ~isempty(opt) && isfield(opt, 'osqp_opt') && ~isempty(opt.osqp_opt) 0181 osqp_opt = nested_struct_copy(osqp_opt, opt.osqp_opt); 0182 end 0183 if verbose > 1 0184 osqp_opt.verbose = 1; 0185 % if verbose > 2 0186 % osqp_opt. 0187 % else 0188 % osqp_opt. 0189 % end 0190 else 0191 osqp_opt.verbose = 0; 0192 end 0193 0194 %% add variable bounds to linear constraints 0195 if isempty(xmin) 0196 if isempty(xmax) 0197 k = []; 0198 else 0199 k = find(xmax < Inf); 0200 end 0201 else 0202 if isempty(xmax) 0203 k = find(xmin > -Inf); 0204 else 0205 k = find(xmin > -Inf | xmax < Inf); 0206 end 0207 end 0208 nv = length(k); 0209 if nv 0210 Av = sparse(1:nv, k, 1, nv, nx); 0211 if isempty(xmin) 0212 lv = -Inf(nv,1); 0213 else 0214 lv = xmin(k); 0215 end 0216 if isempty(xmax) 0217 uv = Inf(nv,1); 0218 else 0219 uv = xmax(k); 0220 end 0221 AA = [A; Av]; 0222 ll = [l; lv]; 0223 uu = [u; uv]; 0224 else 0225 AA = A; 0226 ll = l; 0227 uu = u; 0228 end 0229 0230 %% set up model 0231 o = osqp(); 0232 o.setup(H, c, AA, ll, uu, osqp_opt); 0233 if ~isempty(x0) 0234 o.warm_start('x', x0); 0235 end 0236 0237 %% solve model 0238 if verbose 0239 vn = osqpver; 0240 fprintf('OSQP Version %s -- %s solver\n', vn, lpqp); 0241 end 0242 res = o.solve(); 0243 if verbose 0244 fprintf('OSQP solution status: %s\n', res.info.status); 0245 end 0246 0247 %% extract results 0248 x = res.x; 0249 f = res.info.obj_val; 0250 eflag = res.info.status_val; 0251 if eflag > 1 0252 eflag = -10 * eflag; 0253 end 0254 output = res; 0255 0256 %% separate duals into binding lower & upper bounds 0257 mu_l = -res.y; 0258 mu_u = res.y; 0259 mu_l(mu_l < 0) = 0; 0260 mu_u(mu_u < 0) = 0; 0261 0262 %% variable bounds 0263 lam_lower = zeros(nx, 1); 0264 lam_lower(k) = mu_l(nA+1:end); 0265 lam_upper = zeros(nx, 1); 0266 lam_upper(k) = mu_u(nA+1:end); 0267 0268 %% package lambdas 0269 lambda = struct( ... 0270 'mu_l', mu_l(1:nA), ... 0271 'mu_u', mu_u(1:nA), ... 0272 'lower', lam_lower, ... 0273 'upper', lam_upper ... 0274 );