QPS_IPOPT Quadratic Program Solver based on IPOPT.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_IPOPT(H, C, A, L, U, XMIN, XMAX, X0, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_IPOPT(PROBLEM)
A wrapper function providing a standardized interface for using
IPOPT 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
ipopt_opt - options struct for IPOPT, 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 = converged
0 = failed to converge
OUTPUT : output struct with the following fields:
status - see IPOPT documentation for INFO.status
https://coin-or.github.io/Ipopt/IpReturnCodes__inc_8h_source.html
iterations - number of iterations performed (INFO.iter)
cpu - see IPOPT documentation for INFO.cpu
eval - see IPOPT documentation for INFO.eval
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_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_ipopt(H, c, A, l, u)
x = qps_ipopt(H, c, A, l, u, xmin, xmax)
x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0)
x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_ipopt(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_ipopt(...)
[x, f] = qps_ipopt(...)
[x, f, exitflag] = qps_ipopt(...)
[x, f, exitflag, output] = qps_ipopt(...)
[x, f, exitflag, output, lambda] = qps_ipopt(...)
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_ipopt(H, c, A, l, u, xmin, [], x0, opt);
See also QPS_MASTER, IPOPT, IPOPT_OPTIONS.
https://github.com/coin-or/Ipopt.0001 function [x, f, eflag, output, lambda] = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_IPOPT Quadratic Program Solver based on IPOPT. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_IPOPT(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_IPOPT(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % IPOPT 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 % ipopt_opt - options struct for IPOPT, 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 0043 % 1 = converged 0044 % 0 = failed to converge 0045 % OUTPUT : output struct with the following fields: 0046 % status - see IPOPT documentation for INFO.status 0047 % https://coin-or.github.io/Ipopt/IpReturnCodes__inc_8h_source.html 0048 % iterations - number of iterations performed (INFO.iter) 0049 % cpu - see IPOPT documentation for INFO.cpu 0050 % eval - see IPOPT documentation for INFO.eval 0051 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0052 % multipliers on the constraints, with fields: 0053 % mu_l - lower (left-hand) limit on linear constraints 0054 % mu_u - upper (right-hand) limit on linear constraints 0055 % lower - lower bound on optimization variables 0056 % upper - upper bound on optimization variables 0057 % 0058 % Note the calling syntax is almost identical to that of QUADPROG 0059 % from MathWorks' Optimization Toolbox. The main difference is that 0060 % the linear constraints are specified with A, L, U instead of 0061 % A, B, Aeq, Beq. 0062 % 0063 % Calling syntax options: 0064 % [x, f, exitflag, output, lambda] = ... 0065 % qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 0066 % 0067 % x = qps_ipopt(H, c, A, l, u) 0068 % x = qps_ipopt(H, c, A, l, u, xmin, xmax) 0069 % x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0) 0070 % x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 0071 % x = qps_ipopt(problem), where problem is a struct with fields: 0072 % H, c, A, l, u, xmin, xmax, x0, opt 0073 % all fields except 'c', 'A' and 'l' or 'u' are optional 0074 % x = qps_ipopt(...) 0075 % [x, f] = qps_ipopt(...) 0076 % [x, f, exitflag] = qps_ipopt(...) 0077 % [x, f, exitflag, output] = qps_ipopt(...) 0078 % [x, f, exitflag, output, lambda] = qps_ipopt(...) 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_ipopt(H, c, A, l, u, xmin, [], x0, opt); 0094 % 0095 % See also QPS_MASTER, IPOPT, IPOPT_OPTIONS. 0096 % https://github.com/coin-or/Ipopt. 0097 0098 % MP-Opt-Model 0099 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0100 % by Ray Zimmerman, PSERC Cornell 0101 % 0102 % This file is part of MP-Opt-Model. 0103 % Covered by the 3-clause BSD License (see LICENSE file for details). 0104 % See https://github.com/MATPOWER/mp-opt-model for more info. 0105 0106 %% check for IPOPT 0107 % if ~have_feature('ipopt') 0108 % error('qps_ipopt: requires IPOPT (https://github.com/coin-or/Ipopt)'); 0109 % end 0110 0111 %%----- input argument handling ----- 0112 %% gather inputs 0113 if nargin == 1 && isstruct(H) %% problem struct 0114 p = H; 0115 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0116 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0117 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0118 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0119 if isfield(p, 'u'), u = p.u; else, u = []; end 0120 if isfield(p, 'l'), l = p.l; else, l = []; end 0121 if isfield(p, 'A'), A = p.A; else, A = []; end 0122 if isfield(p, 'c'), c = p.c; else, c = []; end 0123 if isfield(p, 'H'), H = p.H; else, H = []; end 0124 else %% individual args 0125 if nargin < 9 0126 opt = []; 0127 if nargin < 8 0128 x0 = []; 0129 if nargin < 7 0130 xmax = []; 0131 if nargin < 6 0132 xmin = []; 0133 end 0134 end 0135 end 0136 end 0137 end 0138 0139 %% define nx, set default values for missing optional inputs 0140 if isempty(H) || ~any(any(H)) 0141 if isempty(A) && isempty(xmin) && isempty(xmax) 0142 error('qps_ipopt: 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 H = sparse(nx,nx); 0153 else 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 nA 0165 if isempty(u) %% By default, linear inequalities are ... 0166 u = Inf(nA, 1); %% ... unbounded above and ... 0167 end 0168 if isempty(l) 0169 l = -Inf(nA, 1); %% ... unbounded below. 0170 end 0171 end 0172 if isempty(x0) 0173 x0 = zeros(nx, 1); 0174 end 0175 0176 %% default options 0177 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0178 verbose = opt.verbose; 0179 else 0180 verbose = 0; 0181 end 0182 0183 %% make sure args are sparse/full as expected by IPOPT 0184 if ~isempty(H) 0185 if ~issparse(H) 0186 H = sparse(H); 0187 end 0188 end 0189 if ~issparse(A) 0190 A = sparse(A); 0191 end 0192 0193 0194 %%----- run optimization ----- 0195 %% set options struct for IPOPT 0196 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt) 0197 options.ipopt = ipopt_options(opt.ipopt_opt); 0198 else 0199 options.ipopt = ipopt_options; 0200 end 0201 options.ipopt.jac_c_constant = 'yes'; 0202 options.ipopt.jac_d_constant = 'yes'; 0203 options.ipopt.hessian_constant = 'yes'; 0204 options.ipopt.least_square_init_primal = 'yes'; 0205 options.ipopt.least_square_init_duals = 'yes'; 0206 % options.ipopt.mehrotra_algorithm = 'yes'; %% default 'no' 0207 if verbose 0208 options.ipopt.print_level = min(12, verbose*2+1); 0209 else 0210 options.ipopt.print_level = 0; 0211 end 0212 0213 %% define variable and constraint bounds, if given 0214 if nA 0215 options.cu = u; 0216 options.cl = l; 0217 end 0218 if ~isempty(xmin) 0219 options.lb = xmin; 0220 end 0221 if ~isempty(xmax) 0222 options.ub = xmax; 0223 end 0224 0225 %% assign function handles 0226 funcs.objective = @(x) 0.5 * x' * H * x + c' * x; 0227 funcs.gradient = @(x) H * x + c; 0228 funcs.constraints = @(x) A * x; 0229 funcs.jacobian = @(x) A; 0230 funcs.jacobianstructure = @() A; 0231 funcs.hessian = @(x, sigma, lambda) tril(H); 0232 funcs.hessianstructure = @() tril(H); 0233 0234 %% run the optimization 0235 [x, info] = ipopt(x0,funcs,options); 0236 0237 if info.status == 0 || info.status == 1 0238 eflag = 1; 0239 else 0240 eflag = 0; 0241 end 0242 output = struct('status', info.status); 0243 if isfield(info, 'iter') 0244 output.iterations = info.iter; 0245 else 0246 output.iterations = []; 0247 end 0248 if isfield(info, 'cpu') 0249 output.cpu = info.cpu; 0250 end 0251 if isfield(info, 'eval') 0252 output.eval = info.eval; 0253 end 0254 f = funcs.objective(x); 0255 0256 %% check for empty results (in case optimization failed) 0257 if isempty(info.lambda) 0258 lam = NaN(nA, 1); 0259 else 0260 lam = info.lambda; 0261 end 0262 if isempty(info.zl) && ~isempty(xmin) 0263 zl = NaN(nx, 1); 0264 else 0265 zl = info.zl; 0266 end 0267 if isempty(info.zu) && ~isempty(xmax) 0268 zu = NaN(nx, 1); 0269 else 0270 zu = info.zu; 0271 end 0272 0273 %% repackage lambdas 0274 kl = find(lam < 0); %% lower bound binding 0275 ku = find(lam > 0); %% upper bound binding 0276 mu_l = zeros(nA, 1); 0277 mu_l(kl) = -lam(kl); 0278 mu_u = zeros(nA, 1); 0279 mu_u(ku) = lam(ku); 0280 0281 lambda = struct( ... 0282 'mu_l', mu_l, ... 0283 'mu_u', mu_u, ... 0284 'lower', zl, ... 0285 'upper', zu );