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.
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 % 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 % https://projects.coin-or.org/Ipopt/. 0096 0097 % MATPOWER 0098 % $Id: qps_ipopt.m 2340 2014-06-27 19:15:10Z ray $ 0099 % by Ray Zimmerman, PSERC Cornell 0100 % Copyright (c) 2010 by Power System Engineering Research Center (PSERC) 0101 % 0102 % This file is part of MATPOWER. 0103 % See http://www.pserc.cornell.edu/matpower/ for more info. 0104 % 0105 % MATPOWER is free software: you can redistribute it and/or modify 0106 % it under the terms of the GNU General Public License as published 0107 % by the Free Software Foundation, either version 3 of the License, 0108 % or (at your option) any later version. 0109 % 0110 % MATPOWER is distributed in the hope that it will be useful, 0111 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0112 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0113 % GNU General Public License for more details. 0114 % 0115 % You should have received a copy of the GNU General Public License 0116 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0117 % 0118 % Additional permission under GNU GPL version 3 section 7 0119 % 0120 % If you modify MATPOWER, or any covered work, to interface with 0121 % other modules (such as MATLAB code and MEX-files) available in a 0122 % MATLAB(R) or comparable environment containing parts covered 0123 % under other licensing terms, the licensors of MATPOWER grant 0124 % you additional permission to convey the resulting work. 0125 0126 %% check for IPOPT 0127 % if ~have_fcn('ipopt') 0128 % error('qps_ipopt: requires IPOPT (https://projects.coin-or.org/Ipopt/)'); 0129 % end 0130 0131 %%----- input argument handling ----- 0132 %% gather inputs 0133 if nargin == 1 && isstruct(H) %% problem struct 0134 p = H; 0135 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0136 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0137 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0138 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0139 if isfield(p, 'u'), u = p.u; else, u = []; end 0140 if isfield(p, 'l'), l = p.l; else, l = []; end 0141 if isfield(p, 'A'), A = p.A; else, A = []; end 0142 if isfield(p, 'c'), c = p.c; else, c = []; end 0143 if isfield(p, 'H'), H = p.H; else, H = []; end 0144 else %% individual args 0145 if nargin < 9 0146 opt = []; 0147 if nargin < 8 0148 x0 = []; 0149 if nargin < 7 0150 xmax = []; 0151 if nargin < 6 0152 xmin = []; 0153 end 0154 end 0155 end 0156 end 0157 end 0158 0159 %% define nx, set default values for missing optional inputs 0160 if isempty(H) || ~any(any(H)) 0161 if isempty(A) && isempty(xmin) && isempty(xmax) 0162 error('qps_ipopt: LP problem must include constraints or variable bounds'); 0163 else 0164 if ~isempty(A) 0165 nx = size(A, 2); 0166 elseif ~isempty(xmin) 0167 nx = length(xmin); 0168 else % if ~isempty(xmax) 0169 nx = length(xmax); 0170 end 0171 end 0172 H = sparse(nx,nx); 0173 else 0174 nx = size(H, 1); 0175 end 0176 if isempty(c) 0177 c = zeros(nx, 1); 0178 end 0179 if ~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0180 (isempty(u) || all(u == Inf)) 0181 A = sparse(0,nx); %% no limits => no linear constraints 0182 end 0183 nA = size(A, 1); %% number of original linear constraints 0184 if nA 0185 if isempty(u) %% By default, linear inequalities are ... 0186 u = Inf * ones(nA, 1); %% ... unbounded above and ... 0187 end 0188 if isempty(l) 0189 l = -Inf * ones(nA, 1); %% ... unbounded below. 0190 end 0191 end 0192 if isempty(x0) 0193 x0 = zeros(nx, 1); 0194 end 0195 0196 %% default options 0197 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0198 verbose = opt.verbose; 0199 else 0200 verbose = 0; 0201 end 0202 0203 %% make sure args are sparse/full as expected by IPOPT 0204 if ~isempty(H) 0205 if ~issparse(H) 0206 H = sparse(H); 0207 end 0208 end 0209 if ~issparse(A) 0210 A = sparse(A); 0211 end 0212 0213 0214 %%----- run optimization ----- 0215 %% set options struct for IPOPT 0216 if ~isempty(opt) && isfield(opt, 'ipopt_opt') && ~isempty(opt.ipopt_opt) 0217 options.ipopt = ipopt_options(opt.ipopt_opt); 0218 else 0219 options.ipopt = ipopt_options; 0220 end 0221 options.ipopt.jac_c_constant = 'yes'; 0222 options.ipopt.jac_d_constant = 'yes'; 0223 options.ipopt.hessian_constant = 'yes'; 0224 options.ipopt.least_square_init_primal = 'yes'; 0225 options.ipopt.least_square_init_duals = 'yes'; 0226 % options.ipopt.mehrotra_algorithm = 'yes'; %% default 'no' 0227 if verbose 0228 options.ipopt.print_level = min(12, verbose*2+1); 0229 else 0230 options.ipopt.print_level = 0; 0231 end 0232 0233 %% define variable and constraint bounds, if given 0234 if nA 0235 options.cu = u; 0236 options.cl = l; 0237 end 0238 if ~isempty(xmin) 0239 options.lb = xmin; 0240 end 0241 if ~isempty(xmax) 0242 options.ub = xmax; 0243 end 0244 0245 %% assign function handles 0246 funcs.objective = @(x) 0.5 * x' * H * x + c' * x; 0247 funcs.gradient = @(x) H * x + c; 0248 funcs.constraints = @(x) A * x; 0249 funcs.jacobian = @(x) A; 0250 funcs.jacobianstructure = @() A; 0251 funcs.hessian = @(x, sigma, lambda) tril(H); 0252 funcs.hessianstructure = @() tril(H); 0253 0254 %% run the optimization 0255 [x, info] = ipopt(x0,funcs,options); 0256 0257 if info.status == 0 || info.status == 1 0258 eflag = 1; 0259 else 0260 eflag = 0; 0261 end 0262 if isfield(info, 'iter') 0263 output.iterations = info.iter; 0264 end 0265 output.info = info.status; 0266 f = funcs.objective(x); 0267 0268 %% check for empty results (in case optimization failed) 0269 if isempty(info.lambda) 0270 lam = NaN(nA, 1); 0271 else 0272 lam = info.lambda; 0273 end 0274 if isempty(info.zl) && ~isempty(xmin) 0275 zl = NaN(nx, 1); 0276 else 0277 zl = info.zl; 0278 end 0279 if isempty(info.zu) && ~isempty(xmax) 0280 zu = NaN(nx, 1); 0281 else 0282 zu = info.zu; 0283 end 0284 0285 %% repackage lambdas 0286 kl = find(lam < 0); %% lower bound binding 0287 ku = find(lam > 0); %% upper bound binding 0288 mu_l = zeros(nA, 1); 0289 mu_l(kl) = -lam(kl); 0290 mu_u = zeros(nA, 1); 0291 mu_u(ku) = lam(ku); 0292 0293 lambda = struct( ... 0294 'mu_l', mu_l, ... 0295 'mu_u', mu_u, ... 0296 'lower', zl, ... 0297 'upper', zu );