Home > matpower6.0 > qps_mips.m

qps_mips

PURPOSE ^

QPS_MIPS Quadratic Program Solver based on MIPS.

SYNOPSIS ^

function [x, f, eflag, output, lambda] = qps_mips(H, c, A, l, u, xmin, xmax, x0, opt)

DESCRIPTION ^

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 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 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-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 http://www.pserc.cornell.edu/matpower/ 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

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005