QPS_IPOPT Quadratic Program Solver based on IPOPT.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_IPOPT(H, C, A, L, U, XMIN, XMAX, X0, OPT)
Uses 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 = first order optimality conditions satisfied
0 = maximum number of iterations reached
-1 = numerically failed
OUTPUT : output struct with the following fields:
iterations - number of iterations performed
hist - struct array with trajectories of the following:
feascond, gradcond, compcond, costcond, gamma,
stepsize, obj, alphap, alphad
message - exit message
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 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 % Uses IPOPT to solve the following QP (quadratic programming) problem: 0006 % 0007 % min 1/2 X'*H*X + C'*X 0008 % X 0009 % 0010 % subject to 0011 % 0012 % L <= A*X <= U (linear constraints) 0013 % XMIN <= X <= XMAX (variable bounds) 0014 % 0015 % Inputs (all optional except H, C, A and L): 0016 % H : matrix (possibly sparse) of quadratic cost coefficients 0017 % C : vector of linear cost coefficients 0018 % A, L, U : define the optional linear constraints. Default 0019 % values for the elements of L and U are -Inf and Inf, 0020 % respectively. 0021 % XMIN, XMAX : optional lower and upper bounds on the 0022 % X variables, defaults are -Inf and Inf, respectively. 0023 % X0 : optional starting value of optimization vector X 0024 % OPT : optional options structure with the following fields, 0025 % all of which are also optional (default values shown in 0026 % parentheses) 0027 % verbose (0) - controls level of progress output displayed 0028 % 0 = no progress output 0029 % 1 = some progress output 0030 % 2 = verbose progress output 0031 % ipopt_opt - options struct for IPOPT, value in verbose 0032 % overrides these options 0033 % PROBLEM : The inputs can alternatively be supplied in a single 0034 % PROBLEM struct with fields corresponding to the input arguments 0035 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0036 % 0037 % Outputs: 0038 % X : solution vector 0039 % F : final objective function value 0040 % EXITFLAG : exit flag 0041 % 1 = first order optimality conditions satisfied 0042 % 0 = maximum number of iterations reached 0043 % -1 = numerically failed 0044 % OUTPUT : output struct with the following fields: 0045 % iterations - number of iterations performed 0046 % hist - struct array with trajectories of the following: 0047 % feascond, gradcond, compcond, costcond, gamma, 0048 % stepsize, obj, alphap, alphad 0049 % message - exit message 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_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 0065 % 0066 % x = qps_ipopt(H, c, A, l, u) 0067 % x = qps_ipopt(H, c, A, l, u, xmin, xmax) 0068 % x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0) 0069 % x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 0070 % x = qps_ipopt(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_ipopt(...) 0074 % [x, f] = qps_ipopt(...) 0075 % [x, f, exitflag] = qps_ipopt(...) 0076 % [x, f, exitflag, output] = qps_ipopt(...) 0077 % [x, f, exitflag, output, lambda] = qps_ipopt(...) 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_ipopt(H, c, A, l, u, xmin, [], x0, opt); 0093 % 0094 % See also IPOPT, IPOPT_OPTIONS. 0095 % https://github.com/coin-or/Ipopt. 0096 0097 % MATPOWER 0098 % Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC) 0099 % by Ray Zimmerman, PSERC Cornell 0100 % 0101 % This file is part of MATPOWER. 0102 % Covered by the 3-clause BSD License (see LICENSE file for details). 0103 % See https://matpower.org for more info. 0104 0105 %% check for IPOPT 0106 % if ~have_fcn('ipopt') 0107 % error('qps_ipopt: requires IPOPT (https://github.com/coin-or/Ipopt)'); 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 if isempty(A) && isempty(xmin) && isempty(xmax) 0141 error('qps_ipopt: LP problem must include constraints or variable bounds'); 0142 else 0143 if ~isempty(A) 0144 nx = size(A, 2); 0145 elseif ~isempty(xmin) 0146 nx = length(xmin); 0147 else % if ~isempty(xmax) 0148 nx = length(xmax); 0149 end 0150 end 0151 H = sparse(nx,nx); 0152 else 0153 nx = size(H, 1); 0154 end 0155 if isempty(c) 0156 c = zeros(nx, 1); 0157 end 0158 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0159 (isempty(u) || all(u == Inf))) 0160 A = sparse(0,nx); %% no limits => no linear constraints 0161 end 0162 nA = size(A, 1); %% number of original linear constraints 0163 if nA 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 end 0171 if isempty(x0) 0172 x0 = zeros(nx, 1); 0173 end 0174 0175 %% default options 0176 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0177 verbose = opt.verbose; 0178 else 0179 verbose = 0; 0180 end 0181 0182 %% make sure args are sparse/full as expected by IPOPT 0183 if ~isempty(H) 0184 if ~issparse(H) 0185 H = sparse(H); 0186 end 0187 end 0188 if ~issparse(A) 0189 A = sparse(A); 0190 end 0191 0192 0193 %%----- run optimization ----- 0194 %% set options struct for IPOPT 0195 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt) 0196 options.ipopt = ipopt_options(opt.ipopt_opt); 0197 else 0198 options.ipopt = ipopt_options; 0199 end 0200 options.ipopt.jac_c_constant = 'yes'; 0201 options.ipopt.jac_d_constant = 'yes'; 0202 options.ipopt.hessian_constant = 'yes'; 0203 options.ipopt.least_square_init_primal = 'yes'; 0204 options.ipopt.least_square_init_duals = 'yes'; 0205 % options.ipopt.mehrotra_algorithm = 'yes'; %% default 'no' 0206 if verbose 0207 options.ipopt.print_level = min(12, verbose*2+1); 0208 else 0209 options.ipopt.print_level = 0; 0210 end 0211 0212 %% define variable and constraint bounds, if given 0213 if nA 0214 options.cu = u; 0215 options.cl = l; 0216 end 0217 if ~isempty(xmin) 0218 options.lb = xmin; 0219 end 0220 if ~isempty(xmax) 0221 options.ub = xmax; 0222 end 0223 0224 %% assign function handles 0225 funcs.objective = @(x) 0.5 * x' * H * x + c' * x; 0226 funcs.gradient = @(x) H * x + c; 0227 funcs.constraints = @(x) A * x; 0228 funcs.jacobian = @(x) A; 0229 funcs.jacobianstructure = @() A; 0230 funcs.hessian = @(x, sigma, lambda) tril(H); 0231 funcs.hessianstructure = @() tril(H); 0232 0233 %% run the optimization 0234 [x, info] = ipopt(x0,funcs,options); 0235 0236 if info.status == 0 || info.status == 1 0237 eflag = 1; 0238 else 0239 eflag = 0; 0240 end 0241 if isfield(info, 'iter') 0242 output.iterations = info.iter; 0243 end 0244 output.info = info.status; 0245 f = funcs.objective(x); 0246 0247 %% check for empty results (in case optimization failed) 0248 if isempty(info.lambda) 0249 lam = NaN(nA, 1); 0250 else 0251 lam = info.lambda; 0252 end 0253 if isempty(info.zl) && ~isempty(xmin) 0254 zl = NaN(nx, 1); 0255 else 0256 zl = info.zl; 0257 end 0258 if isempty(info.zu) && ~isempty(xmax) 0259 zu = NaN(nx, 1); 0260 else 0261 zu = info.zu; 0262 end 0263 0264 %% repackage lambdas 0265 kl = find(lam < 0); %% lower bound binding 0266 ku = find(lam > 0); %% upper bound binding 0267 mu_l = zeros(nA, 1); 0268 mu_l(kl) = -lam(kl); 0269 mu_u = zeros(nA, 1); 0270 mu_u(ku) = lam(ku); 0271 0272 lambda = struct( ... 0273 'mu_l', mu_l, ... 0274 'mu_u', mu_u, ... 0275 'lower', zl, ... 0276 'upper', zu );