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
max_it (0) - maximum number of iterations allowed
0 = use algorithm default
ipopt_opt - options struct for IPOPT, values in
verbose and max_it override 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, lam] = qps_ipopt(H, c, A, l, u, xmin, [], x0, opt);
See also IPOPT, IPOPT_OPTIONS.
https://projects.coin-or.org/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 % max_it (0) - maximum number of iterations allowed 0032 % 0 = use algorithm default 0033 % ipopt_opt - options struct for IPOPT, values in 0034 % verbose and max_it override 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 = first order optimality conditions satisfied 0044 % 0 = maximum number of iterations reached 0045 % -1 = numerically failed 0046 % OUTPUT : output struct with the following fields: 0047 % iterations - number of iterations performed 0048 % hist - struct array with trajectories of the following: 0049 % feascond, gradcond, compcond, costcond, gamma, 0050 % stepsize, obj, alphap, alphad 0051 % message - exit message 0052 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0053 % multipliers on the constraints, with fields: 0054 % mu_l - lower (left-hand) limit on linear constraints 0055 % mu_u - upper (right-hand) limit on linear constraints 0056 % lower - lower bound on optimization variables 0057 % upper - upper bound on optimization variables 0058 % 0059 % Note the calling syntax is almost identical to that of QUADPROG 0060 % from MathWorks' Optimization Toolbox. The main difference is that 0061 % the linear constraints are specified with A, L, U instead of 0062 % A, B, Aeq, Beq. 0063 % 0064 % Calling syntax options: 0065 % [x, f, exitflag, output, lambda] = ... 0066 % qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 0067 % 0068 % x = qps_ipopt(H, c, A, l, u) 0069 % x = qps_ipopt(H, c, A, l, u, xmin, xmax) 0070 % x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0) 0071 % x = qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt) 0072 % x = qps_ipopt(problem), where problem is a struct with fields: 0073 % H, c, A, l, u, xmin, xmax, x0, opt 0074 % all fields except 'c', 'A' and 'l' or 'u' are optional 0075 % x = qps_ipopt(...) 0076 % [x, f] = qps_ipopt(...) 0077 % [x, f, exitflag] = qps_ipopt(...) 0078 % [x, f, exitflag, output] = qps_ipopt(...) 0079 % [x, f, exitflag, output, lambda] = qps_ipopt(...) 0080 % 0081 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0082 % H = [ 1003.1 4.3 6.3 5.9; 0083 % 4.3 2.2 2.1 3.9; 0084 % 6.3 2.1 3.5 4.8; 0085 % 5.9 3.9 4.8 10 ]; 0086 % c = zeros(4,1); 0087 % A = [ 1 1 1 1; 0088 % 0.17 0.11 0.10 0.18 ]; 0089 % l = [1; 0.10]; 0090 % u = [1; Inf]; 0091 % xmin = zeros(4,1); 0092 % x0 = [1; 0; 0; 1]; 0093 % opt = struct('verbose', 2); 0094 % [x, f, s, out, lam] = qps_ipopt(H, c, A, l, u, xmin, [], x0, opt); 0095 % 0096 % See also IPOPT, IPOPT_OPTIONS. 0097 % https://projects.coin-or.org/Ipopt/. 0098 0099 % MATPOWER 0100 % $Id: qps_ipopt.m,v 1.5 2010/12/15 18:40:42 cvs Exp $ 0101 % by Ray Zimmerman, PSERC Cornell 0102 % Copyright (c) 2010 by Power System Engineering Research Center (PSERC) 0103 % 0104 % This file is part of MATPOWER. 0105 % See http://www.pserc.cornell.edu/matpower/ for more info. 0106 % 0107 % MATPOWER is free software: you can redistribute it and/or modify 0108 % it under the terms of the GNU General Public License as published 0109 % by the Free Software Foundation, either version 3 of the License, 0110 % or (at your option) any later version. 0111 % 0112 % MATPOWER is distributed in the hope that it will be useful, 0113 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0114 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0115 % GNU General Public License for more details. 0116 % 0117 % You should have received a copy of the GNU General Public License 0118 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0119 % 0120 % Additional permission under GNU GPL version 3 section 7 0121 % 0122 % If you modify MATPOWER, or any covered work, to interface with 0123 % other modules (such as MATLAB code and MEX-files) available in a 0124 % MATLAB(R) or comparable environment containing parts covered 0125 % under other licensing terms, the licensors of MATPOWER grant 0126 % you additional permission to convey the resulting work. 0127 0128 %% check for IPOPT 0129 % if ~have_fcn('ipopt') 0130 % error('qps_ipopt: requires IPOPT (https://projects.coin-or.org/Ipopt/)'); 0131 % end 0132 0133 %%----- input argument handling ----- 0134 %% gather inputs 0135 if nargin == 1 && isstruct(H) %% problem struct 0136 p = H; 0137 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0138 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0139 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0140 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0141 if isfield(p, 'u'), u = p.u; else, u = []; end 0142 if isfield(p, 'l'), l = p.l; else, l = []; end 0143 if isfield(p, 'A'), A = p.A; else, A = []; end 0144 if isfield(p, 'c'), c = p.c; else, c = []; end 0145 if isfield(p, 'H'), H = p.H; else, H = []; end 0146 else %% individual args 0147 if nargin < 9 0148 opt = []; 0149 if nargin < 8 0150 x0 = []; 0151 if nargin < 7 0152 xmax = []; 0153 if nargin < 6 0154 xmin = []; 0155 end 0156 end 0157 end 0158 end 0159 end 0160 0161 %% define nx, set default values for missing optional inputs 0162 if isempty(H) || ~any(any(H)) 0163 if isempty(A) && isempty(xmin) && isempty(xmax) 0164 error('qps_ipopt: LP problem must include constraints or variable bounds'); 0165 else 0166 if ~isempty(A) 0167 nx = size(A, 2); 0168 elseif ~isempty(xmin) 0169 nx = length(xmin); 0170 else % if ~isempty(xmax) 0171 nx = length(xmax); 0172 end 0173 end 0174 H = sparse(nx,nx); 0175 else 0176 nx = size(H, 1); 0177 end 0178 if isempty(c) 0179 c = zeros(nx, 1); 0180 end 0181 if ~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0182 (isempty(u) || all(u == Inf)) 0183 A = sparse(0,nx); %% no limits => no linear constraints 0184 end 0185 nA = size(A, 1); %% number of original linear constraints 0186 if nA 0187 if isempty(u) %% By default, linear inequalities are ... 0188 u = Inf * ones(nA, 1); %% ... unbounded above and ... 0189 end 0190 if isempty(l) 0191 l = -Inf * ones(nA, 1); %% ... unbounded below. 0192 end 0193 end 0194 if isempty(x0) 0195 x0 = zeros(nx, 1); 0196 end 0197 0198 %% default options 0199 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0200 verbose = opt.verbose; 0201 else 0202 verbose = 0; 0203 end 0204 if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it) 0205 max_it = opt.max_it; 0206 else 0207 max_it = 0; 0208 end 0209 0210 %% make sure args are sparse/full as expected by IPOPT 0211 if ~isempty(H) 0212 if ~issparse(H) 0213 H = sparse(H); 0214 end 0215 end 0216 if ~issparse(A) 0217 A = sparse(A); 0218 end 0219 0220 0221 %%----- run optimization ----- 0222 %% set options struct for IPOPT 0223 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt) 0224 options.ipopt = ipopt_options(opt.ipopt_opt); 0225 else 0226 options.ipopt = ipopt_options; 0227 end 0228 options.ipopt.jac_c_constant = 'yes'; 0229 options.ipopt.jac_d_constant = 'yes'; 0230 options.ipopt.hessian_constant = 'yes'; 0231 options.ipopt.least_square_init_primal = 'yes'; 0232 options.ipopt.least_square_init_duals = 'yes'; 0233 % options.ipopt.mehrotra_algorithm = 'yes'; %% default 'no' 0234 if verbose 0235 options.ipopt.print_level = min(12, verbose*2+1); 0236 else 0237 options.ipopt.print_level = 0; 0238 end 0239 if max_it 0240 options.ipopt.max_iter = max_it; 0241 end 0242 0243 %% define variable and constraint bounds, if given 0244 if nA 0245 options.cu = u; 0246 options.cl = l; 0247 end 0248 if ~isempty(xmin) 0249 options.lb = xmin; 0250 end 0251 if ~isempty(xmax) 0252 options.ub = xmax; 0253 end 0254 0255 %% assign function handles 0256 funcs.objective = @(x) 0.5 * x' * H * x + c' * x; 0257 funcs.gradient = @(x) H * x + c; 0258 funcs.constraints = @(x) A * x; 0259 funcs.jacobian = @(x) A; 0260 funcs.jacobianstructure = @() A; 0261 funcs.hessian = @(x, sigma, lambda) tril(H); 0262 funcs.hessianstructure = @() tril(H); 0263 0264 %% run the optimization 0265 [x, info] = ipopt(x0,funcs,options); 0266 0267 if info.status == 0 || info.status == 1 0268 eflag = 1; 0269 else 0270 eflag = 0; 0271 end 0272 if isfield(info, 'iter') 0273 output.iterations = info.iter; 0274 end 0275 output.info = info.status; 0276 f = funcs.objective(x); 0277 0278 %% repackage lambdas 0279 kl = find(info.lambda < 0); %% lower bound binding 0280 ku = find(info.lambda > 0); %% upper bound binding 0281 mu_l = zeros(nA, 1); 0282 mu_l(kl) = -info.lambda(kl); 0283 mu_u = zeros(nA, 1); 0284 mu_u(ku) = info.lambda(ku); 0285 0286 lambda = struct( ... 0287 'mu_l', mu_l, ... 0288 'mu_u', mu_u, ... 0289 'lower', info.zl, ... 0290 'upper', info.zu );