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 MATPOWER 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 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_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 MATPOWER 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 https://v8doc.sas.com/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-2016, Power Systems Engineering Research Center (PSERC) 0126 % by Ray Zimmerman, PSERC Cornell 0127 % 0128 % This file is part of MIPS. 0129 % Covered by the 3-clause BSD License (see LICENSE file for details). 0130 % See https://github.com/MATPOWER/mips for more info. 0131 0132 %%----- input argument handling ----- 0133 %% gather inputs 0134 if nargin == 1 && isstruct(H) %% problem struct 0135 p = H; 0136 else %% individual args 0137 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u); 0138 if nargin > 5 0139 p.xmin = xmin; 0140 if nargin > 6 0141 p.xmax = xmax; 0142 if nargin > 7 0143 p.x0 = x0; 0144 if nargin > 8 0145 p.opt = opt; 0146 end 0147 end 0148 end 0149 end 0150 end 0151 0152 %% define nx, set default values for H and c 0153 if ~isfield(p, 'H') || isempty(p.H) || ~any(any(p.H)) 0154 if (~isfield(p, 'A') || isempty(p.A)) && ... 0155 (~isfield(p, 'xmin') || isempty(p.xmin)) && ... 0156 (~isfield(p, 'xmax') || isempty(p.xmax)) 0157 error('qps_mips: LP problem must include constraints or variable bounds'); 0158 else 0159 if isfield(p, 'A') && ~isempty(p.A) 0160 nx = size(p.A, 2); 0161 elseif isfield(p, 'xmin') && ~isempty(p.xmin) 0162 nx = length(p.xmin); 0163 else % if isfield(p, 'xmax') && ~isempty(p.xmax) 0164 nx = length(p.xmax); 0165 end 0166 end 0167 p.H = sparse(nx, nx); 0168 else 0169 nx = size(p.H, 1); 0170 end 0171 if ~isfield(p, 'c') || isempty(p.c) 0172 p.c = zeros(nx, 1); 0173 end 0174 if ~isfield(p, 'x0') || isempty(p.x0) 0175 p.x0 = zeros(nx, 1); 0176 end 0177 0178 %%----- run optimization ----- 0179 p.f_fcn = @(x)qp_f(x, p.H, p.c); 0180 [x, f, eflag, output, lambda] = mips(p); 0181 0182 %%----- objective function ----- 0183 function [f, df, d2f] = qp_f(x, H, c) 0184 f = 0.5 * x' * H * x + c' * x; 0185 if nargout > 1 0186 df = H * x + c; 0187 if nargout > 2 0188 d2f = H; 0189 end 0190 end