QPS_MIPS Quadratic Program Solver based on MIPS.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_MIPS(H, C, A, L, U, XMIN, XMAX, X0, OPT)
Uses the MATLAB Interior Point Solver (MIPS) 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
linsolver ('') - linear system solver for solving update steps
'' = default solver (currently same as '\')
'\' = built-in \ operator
'PARDISO' = PARDISO solver (if available)
feastol (1e-6) - termination tolerance for feasibility
condition
gradtol (1e-6) - termination tolerance for gradient
condition
comptol (1e-6) - termination tolerance for complementarity
condition
costtol (1e-6) - termination tolerance for cost condition
max_it (150) - maximum number of iterations
step_control (0) - set to 1 to enable step-size control
sc.red_it (20) - maximum number of step-size reductions if
step-control is on
cost_mult (1) - cost multiplier used to scale the objective
function for improved conditioning.
xi (0.99995) - constant used in alpha updates*
sigma (0.1) - centering parameter*
z0 (1) - used to initialize slack variables*
alpha_min (1e-8) - returns EXITFLAG = -1 if either alpha
parameter becomes smaller than this value*
rho_min (0.95) - lower bound on rho_t*
rho_max (1.05) - upper bound on rho_t*
mu_threshold (1e-5) - KT multipliers smaller than this value
for non-binding constraints are forced to zero
max_stepsize (1e10) - returns EXITFLAG = -1 if the 2-norm of
the reduced Newton step exceeds this value*
* see the corresponding Appendix in the manual for details
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_mips(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_mips(H, c, A, l, u)
x = qps_mips(H, c, A, l, u, xmin, xmax)
x = qps_mips(H, c, A, l, u, xmin, xmax, x0)
x = qps_mips(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_mips(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_mips(...)
[x, f] = qps_mips(...)
[x, f, exitflag] = qps_mips(...)
[x, f, exitflag, output] = qps_mips(...)
[x, f, exitflag, output, lambda] = qps_mips(...)
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_mips(H, c, A, l, u, xmin, [], x0, opt);
See also MIPS.0001 function [x, f, eflag, output, lambda] = qps_mips(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_MIPS Quadratic Program Solver based on MIPS. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_MIPS(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % Uses the MATLAB Interior Point Solver (MIPS) to solve the following 0006 % QP (quadratic programming) problem: 0007 % 0008 % min 1/2 X'*H*X + C'*X 0009 % X 0010 % 0011 % subject to 0012 % 0013 % L <= A*X <= U (linear constraints) 0014 % XMIN <= X <= XMAX (variable bounds) 0015 % 0016 % Inputs (all optional except H, C, A and L): 0017 % H : matrix (possibly sparse) of quadratic cost coefficients 0018 % C : vector of linear cost coefficients 0019 % A, L, U : define the optional linear constraints. Default 0020 % values for the elements of L and U are -Inf and Inf, 0021 % respectively. 0022 % XMIN, XMAX : optional lower and upper bounds on the 0023 % X variables, defaults are -Inf and Inf, respectively. 0024 % X0 : optional starting value of optimization vector X 0025 % OPT : optional options structure with the following fields, 0026 % all of which are also optional (default values shown in 0027 % parentheses) 0028 % verbose (0) - controls level of progress output displayed 0029 % 0 = no progress output 0030 % 1 = some progress output 0031 % 2 = verbose progress output 0032 % linsolver ('') - linear system solver for solving update steps 0033 % '' = default solver (currently same as '\') 0034 % '\' = built-in \ operator 0035 % 'PARDISO' = PARDISO solver (if available) 0036 % feastol (1e-6) - termination tolerance for feasibility 0037 % condition 0038 % gradtol (1e-6) - termination tolerance for gradient 0039 % condition 0040 % comptol (1e-6) - termination tolerance for complementarity 0041 % condition 0042 % costtol (1e-6) - termination tolerance for cost condition 0043 % max_it (150) - maximum number of iterations 0044 % step_control (0) - set to 1 to enable step-size control 0045 % sc.red_it (20) - maximum number of step-size reductions if 0046 % step-control is on 0047 % cost_mult (1) - cost multiplier used to scale the objective 0048 % function for improved conditioning. 0049 % xi (0.99995) - constant used in alpha updates* 0050 % sigma (0.1) - centering parameter* 0051 % z0 (1) - used to initialize slack variables* 0052 % alpha_min (1e-8) - returns EXITFLAG = -1 if either alpha 0053 % parameter becomes smaller than this value* 0054 % rho_min (0.95) - lower bound on rho_t* 0055 % rho_max (1.05) - upper bound on rho_t* 0056 % mu_threshold (1e-5) - KT multipliers smaller than this value 0057 % for non-binding constraints are forced to zero 0058 % max_stepsize (1e10) - returns EXITFLAG = -1 if the 2-norm of 0059 % the reduced Newton step exceeds this value* 0060 % * see the corresponding Appendix in the manual for details 0061 % PROBLEM : The inputs can alternatively be supplied in a single 0062 % PROBLEM struct with fields corresponding to the input arguments 0063 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0064 % 0065 % Outputs: 0066 % X : solution vector 0067 % F : final objective function value 0068 % EXITFLAG : exit flag 0069 % 1 = first order optimality conditions satisfied 0070 % 0 = maximum number of iterations reached 0071 % -1 = numerically failed 0072 % OUTPUT : output struct with the following fields: 0073 % iterations - number of iterations performed 0074 % hist - struct array with trajectories of the following: 0075 % feascond, gradcond, compcond, costcond, gamma, 0076 % stepsize, obj, alphap, alphad 0077 % message - exit message 0078 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0079 % multipliers on the constraints, with fields: 0080 % mu_l - lower (left-hand) limit on linear constraints 0081 % mu_u - upper (right-hand) limit on linear constraints 0082 % lower - lower bound on optimization variables 0083 % upper - upper bound on optimization variables 0084 % 0085 % Note the calling syntax is almost identical to that of QUADPROG 0086 % from MathWorks' Optimization Toolbox. The main difference is that 0087 % the linear constraints are specified with A, L, U instead of 0088 % A, B, Aeq, Beq. 0089 % 0090 % Calling syntax options: 0091 % [x, f, exitflag, output, lambda] = ... 0092 % qps_mips(H, c, A, l, u, xmin, xmax, x0, opt) 0093 % 0094 % x = qps_mips(H, c, A, l, u) 0095 % x = qps_mips(H, c, A, l, u, xmin, xmax) 0096 % x = qps_mips(H, c, A, l, u, xmin, xmax, x0) 0097 % x = qps_mips(H, c, A, l, u, xmin, xmax, x0, opt) 0098 % x = qps_mips(problem), where problem is a struct with fields: 0099 % H, c, A, l, u, xmin, xmax, x0, opt 0100 % all fields except 'c', 'A' and 'l' or 'u' are optional 0101 % x = qps_mips(...) 0102 % [x, f] = qps_mips(...) 0103 % [x, f, exitflag] = qps_mips(...) 0104 % [x, f, exitflag, output] = qps_mips(...) 0105 % [x, f, exitflag, output, lambda] = qps_mips(...) 0106 % 0107 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0108 % H = [ 1003.1 4.3 6.3 5.9; 0109 % 4.3 2.2 2.1 3.9; 0110 % 6.3 2.1 3.5 4.8; 0111 % 5.9 3.9 4.8 10 ]; 0112 % c = zeros(4,1); 0113 % A = [ 1 1 1 1; 0114 % 0.17 0.11 0.10 0.18 ]; 0115 % l = [1; 0.10]; 0116 % u = [1; Inf]; 0117 % xmin = zeros(4,1); 0118 % x0 = [1; 0; 0; 1]; 0119 % opt = struct('verbose', 2); 0120 % [x, f, s, out, lambda] = qps_mips(H, c, A, l, u, xmin, [], x0, opt); 0121 % 0122 % See also MIPS. 0123 0124 % MIPS 0125 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0126 % by Ray Zimmerman, PSERC Cornell 0127 % 0128 % $Id: qps_mips.m 2660 2015-03-20 02:10:18Z ray $ 0129 % 0130 % This file is part of MIPS. 0131 % Covered by the 3-clause BSD License (see LICENSE file for details). 0132 % See http://www.pserc.cornell.edu/matpower/ for more info. 0133 0134 %%----- input argument handling ----- 0135 %% gather inputs 0136 if nargin == 1 && isstruct(H) %% problem struct 0137 p = H; 0138 else %% individual args 0139 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u); 0140 if nargin > 5 0141 p.xmin = xmin; 0142 if nargin > 6 0143 p.xmax = xmax; 0144 if nargin > 7 0145 p.x0 = x0; 0146 if nargin > 8 0147 p.opt = opt; 0148 end 0149 end 0150 end 0151 end 0152 end 0153 0154 %% define nx, set default values for H and c 0155 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H)) 0156 if (~isfield(p, 'A') || isempty(p.A)) && ... 0157 (~isfield(p, 'xmin') || isempty(p.xmin)) && ... 0158 (~isfield(p, 'xmax') || isempty(p.xmax)) 0159 error('qps_mips: LP problem must include constraints or variable bounds'); 0160 else 0161 if isfield(p, 'A') && ~isempty(p.A) 0162 nx = size(p.A, 2); 0163 elseif isfield(p, 'xmin') && ~isempty(p.xmin) 0164 nx = length(p.xmin); 0165 else % if isfield(p, 'xmax') && ~isempty(p.xmax) 0166 nx = length(p.xmax); 0167 end 0168 end 0169 p.H = sparse(nx, nx); 0170 else 0171 nx = size(p.H, 1); 0172 end 0173 if ~isfield(p, 'c') || isempty(p.c) 0174 p.c = zeros(nx, 1); 0175 end 0176 if ~isfield(p, 'x0') || isempty(p.x0) 0177 p.x0 = zeros(nx, 1); 0178 end 0179 0180 %%----- run optimization ----- 0181 p.f_fcn = @(x)qp_f(x, p.H, p.c); 0182 [x, f, eflag, output, lambda] = mips(p); 0183 0184 %%----- objective function ----- 0185 function [f, df, d2f] = qp_f(x, H, c) 0186 f = 0.5 * x' * H * x + c' * x; 0187 if nargout > 1 0188 df = H * x + c; 0189 if nargout > 2 0190 d2f = H; 0191 end 0192 end