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 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_ipopt(H, c, A, l, u, xmin, [], x0, opt);
See also IPOPT, IPOPT_OPTIONS.
http://www.coin-or.org/projects/Ipopt.xml.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 http://www.jmu.edu/docs/sasdoc/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 % http://www.coin-or.org/projects/Ipopt.xml. 0096 0097 % MATPOWER 0098 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0099 % by Ray Zimmerman, PSERC Cornell 0100 % 0101 % $Id: qps_ipopt.m 2661 2015-03-20 17:02:46Z ray $ 0102 % 0103 % This file is part of MATPOWER. 0104 % Covered by the 3-clause BSD License (see LICENSE file for details). 0105 % See http://www.pserc.cornell.edu/matpower/ for more info. 0106 0107 %% check for IPOPT 0108 % if ~have_fcn('ipopt') 0109 % error('qps_ipopt: requires IPOPT (http://www.coin-or.org/projects/Ipopt.xml)'); 0110 % end 0111 0112 %%----- input argument handling ----- 0113 %% gather inputs 0114 if nargin == 1 && isstruct(H) %% problem struct 0115 p = H; 0116 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0117 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0118 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0119 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0120 if isfield(p, 'u'), u = p.u; else, u = []; end 0121 if isfield(p, 'l'), l = p.l; else, l = []; end 0122 if isfield(p, 'A'), A = p.A; else, A = []; end 0123 if isfield(p, 'c'), c = p.c; else, c = []; end 0124 if isfield(p, 'H'), H = p.H; else, H = []; end 0125 else %% individual args 0126 if nargin < 9 0127 opt = []; 0128 if nargin < 8 0129 x0 = []; 0130 if nargin < 7 0131 xmax = []; 0132 if nargin < 6 0133 xmin = []; 0134 end 0135 end 0136 end 0137 end 0138 end 0139 0140 %% define nx, set default values for missing optional inputs 0141 if isempty(H) || ~any(any(H)) 0142 if isempty(A) && isempty(xmin) && isempty(xmax) 0143 error('qps_ipopt: LP problem must include constraints or variable bounds'); 0144 else 0145 if ~isempty(A) 0146 nx = size(A, 2); 0147 elseif ~isempty(xmin) 0148 nx = length(xmin); 0149 else % if ~isempty(xmax) 0150 nx = length(xmax); 0151 end 0152 end 0153 H = sparse(nx,nx); 0154 else 0155 nx = size(H, 1); 0156 end 0157 if isempty(c) 0158 c = zeros(nx, 1); 0159 end 0160 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0161 (isempty(u) || all(u == Inf))) 0162 A = sparse(0,nx); %% no limits => no linear constraints 0163 end 0164 nA = size(A, 1); %% number of original linear constraints 0165 if nA 0166 if isempty(u) %% By default, linear inequalities are ... 0167 u = Inf(nA, 1); %% ... unbounded above and ... 0168 end 0169 if isempty(l) 0170 l = -Inf(nA, 1); %% ... unbounded below. 0171 end 0172 end 0173 if isempty(x0) 0174 x0 = zeros(nx, 1); 0175 end 0176 0177 %% default options 0178 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0179 verbose = opt.verbose; 0180 else 0181 verbose = 0; 0182 end 0183 0184 %% make sure args are sparse/full as expected by IPOPT 0185 if ~isempty(H) 0186 if ~issparse(H) 0187 H = sparse(H); 0188 end 0189 end 0190 if ~issparse(A) 0191 A = sparse(A); 0192 end 0193 0194 0195 %%----- run optimization ----- 0196 %% set options struct for IPOPT 0197 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt) 0198 options.ipopt = ipopt_options(opt.ipopt_opt); 0199 else 0200 options.ipopt = ipopt_options; 0201 end 0202 options.ipopt.jac_c_constant = 'yes'; 0203 options.ipopt.jac_d_constant = 'yes'; 0204 options.ipopt.hessian_constant = 'yes'; 0205 options.ipopt.least_square_init_primal = 'yes'; 0206 options.ipopt.least_square_init_duals = 'yes'; 0207 % options.ipopt.mehrotra_algorithm = 'yes'; %% default 'no' 0208 if verbose 0209 options.ipopt.print_level = min(12, verbose*2+1); 0210 else 0211 options.ipopt.print_level = 0; 0212 end 0213 0214 %% define variable and constraint bounds, if given 0215 if nA 0216 options.cu = u; 0217 options.cl = l; 0218 end 0219 if ~isempty(xmin) 0220 options.lb = xmin; 0221 end 0222 if ~isempty(xmax) 0223 options.ub = xmax; 0224 end 0225 0226 %% assign function handles 0227 funcs.objective = @(x) 0.5 * x' * H * x + c' * x; 0228 funcs.gradient = @(x) H * x + c; 0229 funcs.constraints = @(x) A * x; 0230 funcs.jacobian = @(x) A; 0231 funcs.jacobianstructure = @() A; 0232 funcs.hessian = @(x, sigma, lambda) tril(H); 0233 funcs.hessianstructure = @() tril(H); 0234 0235 %% run the optimization 0236 [x, info] = ipopt(x0,funcs,options); 0237 0238 if info.status == 0 || info.status == 1 0239 eflag = 1; 0240 else 0241 eflag = 0; 0242 end 0243 if isfield(info, 'iter') 0244 output.iterations = info.iter; 0245 end 0246 output.info = info.status; 0247 f = funcs.objective(x); 0248 0249 %% check for empty results (in case optimization failed) 0250 if isempty(info.lambda) 0251 lam = NaN(nA, 1); 0252 else 0253 lam = info.lambda; 0254 end 0255 if isempty(info.zl) && ~isempty(xmin) 0256 zl = NaN(nx, 1); 0257 else 0258 zl = info.zl; 0259 end 0260 if isempty(info.zu) && ~isempty(xmax) 0261 zu = NaN(nx, 1); 0262 else 0263 zu = info.zu; 0264 end 0265 0266 %% repackage lambdas 0267 kl = find(lam < 0); %% lower bound binding 0268 ku = find(lam > 0); %% upper bound binding 0269 mu_l = zeros(nA, 1); 0270 mu_l(kl) = -lam(kl); 0271 mu_u = zeros(nA, 1); 0272 mu_u(ku) = lam(ku); 0273 0274 lambda = struct( ... 0275 'mu_l', mu_l, ... 0276 'mu_u', mu_u, ... 0277 'lower', zl, ... 0278 'upper', zu );