wolffd@0: function [x, fval, status] = qplcprog(H, f, Ai, bi, Ae, be, lb, ub, varargin)
wolffd@0: %QPLCPROG Positive-semidefinite quadratic programming based on linear complementary programming.
wolffd@0: %
wolffd@0: %	x = qplcprog(H, f);
wolffd@0: %	[x, fval, status] = qplcprog(H, f, Ai, bi, Ae, be, lb, ub, varargin);
wolffd@0: %
wolffd@0: % Inputs:
wolffd@0: %	H - Hessian for objective function (n x n matrix)
wolffd@0: %
wolffd@0: % Optional Inputs:
wolffd@0: %	f - Vector for objective function, default is 0 (n x 1 vector)
wolffd@0: %	Ai - Linear inequality constraint matrix, default is [] (ni x n matrix)
wolffd@0: %	bi - Linear inequality constraint vector, default is [] (ni x 1 vector)
wolffd@0: %	Ae - Linear equality constraint matrix, default is [] (ne x n matrix)
wolffd@0: %	be - Linear equality constraint vector, default is [] (ne x 1 vector)
wolffd@0: %	lb - Lower-bound constraint, default is 0 (n x 1 vector)
wolffd@0: %	ub - Upper-bound constraint, default is [] (n x 1 vector)
wolffd@0: %
wolffd@0: %	maxiter - LCP maximum iterations, default is 100000
wolffd@0: %	tiebreak - Method to break ties for pivot selection (default is 'first'). Options are:
wolffd@0: %			'first'		Select pivot with lowest index.
wolffd@0: %			'last'		Select pivot with highest index.
wolffd@0: %			'random'	Select a pivot at random.
wolffd@0: %	tolpiv - LCP pivot selection tolerance, default is 1.0e-9
wolffd@0: %	tolcon - QP constraint tolerance, default is 1.0e-6
wolffd@0: %
wolffd@0: % Outputs:
wolffd@0: %	x - Solution vector (n x 1 vector)
wolffd@0: %	fval - Value of objective function at solution (scalar)
wolffd@0: %	status - Status of optimization with values (integer)
wolffd@0: %		 6		Ray termination
wolffd@0: %		 1		Convergence to a solution (normal termination)
wolffd@0: %		 0		Termination prior to convergence
wolffd@0: %		-2		Infeasible problem
wolffd@0: %		-3		Bad pivot to an infeasible solution
wolffd@0: %
wolffd@0: % Notes:
wolffd@0: %	Depending upon the inputs, various types of positive-semidefinite quadratic programs can be
wolffd@0: %	solved. The general problem is
wolffd@0: %
wolffd@0: %		min  0.5 * x' * H * x + f' * x
wolffd@0: %		 x
wolffd@0: %
wolffd@0: %	subject to any or all of the following constraints
wolffd@0: %
wolffd@0: %		Ai * x <= bi
wolffd@0: %		Ae * x = be
wolffd@0: %		lb <= x <= ub
wolffd@0: %
wolffd@0: %	with the current requirement that lb must be specified and finite.
wolffd@0: %
wolffd@0: %	It is assumed that H is positive-semidefinite. Note that if H = 0, qplcp solves a linear
wolffd@0: %	program.
wolffd@0: 
wolffd@0: % Copyright 2008 The MathWorks, Inc.
wolffd@0: % $Revision: 1.1.6.3 $   $ Date: $
wolffd@0: 
wolffd@0: if nargin < 1
wolffd@0:    error('Finance:qplcprog:MissingInputArgument', ...
wolffd@0:       'Missing required Hessian matrix for QP.');
wolffd@0: end
wolffd@0: H = full(H);
wolffd@0: n = size(H,1);
wolffd@0: 
wolffd@0: if nargin < 2 || isempty(f)
wolffd@0:    f = zeros(n,1);
wolffd@0: else
wolffd@0:    f = full(f);
wolffd@0:    f = f(:);
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin == 3
wolffd@0:    error('Finance:qplcprog:IncompleteConstraintSpecification', ...
wolffd@0:       'Incomplete specification for linear inequality constraints.');
wolffd@0: end
wolffd@0: if nargin < 3
wolffd@0:    Ai = [];
wolffd@0:    bi = [];
wolffd@0: else
wolffd@0:    Ai = full(Ai);
wolffd@0:    bi = full(bi);
wolffd@0:    bi = bi(:);
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin == 5
wolffd@0:    error('Finance:qplcprog:IncompleteConstraintSpecification', ...
wolffd@0:       'Incomplete specification for linear equality constraints.');
wolffd@0: end
wolffd@0: if nargin < 5
wolffd@0:    Ae = [];
wolffd@0:    be = [];
wolffd@0: else
wolffd@0:    Ae = full(Ae);
wolffd@0:    be = full(be);
wolffd@0:    be = be(:);
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin < 7 || isempty(lb)
wolffd@0:    error('Finance:qplcprog:MissingBoundConstraint', ...
wolffd@0:       'Lower-bound constraint required.');
wolffd@0: else
wolffd@0:    lb = full(lb);
wolffd@0:    lb = lb(:);
wolffd@0:    if any(~isfinite(lb))
wolffd@0:       error('Finance:qplcprog:MissingBoundConstraint', ...
wolffd@0:          'Finite lower-bound constraint required.');
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin < 8
wolffd@0:    ub = [];
wolffd@0: else
wolffd@0:    ub = full(ub);
wolffd@0:    ub = ub(:);
wolffd@0:    if all(ub == Inf)
wolffd@0:       ub = [];
wolffd@0:    end
wolffd@0: end
wolffd@0: if any(~isfinite(ub))
wolffd@0:    error('Finance:qplcprog:InvalidBoundConstraint', ...
wolffd@0:       'If specified, finite upper-bound constraint required.');
wolffd@0: end
wolffd@0: 
wolffd@0: % process name-value pairs (if any)
wolffd@0: 
wolffd@0: if nargin > 8
wolffd@0:    if mod(nargin-8,2) ~= 0
wolffd@0:       error('Finance:qplcprog:InvalidNameValuePair', ...
wolffd@0:          'Invalid name-value pairs with either a missing name or value.');
wolffd@0:    end
wolffd@0: 
wolffd@0:    names = { 'maxiter', 'tiebreak', 'tolcon', 'tolpiv' };			% names
wolffd@0:    values = { 10000, 'first', 1.0e-4, 1.0e-9 };					% default values
wolffd@0:    try
wolffd@0:       [maxiter, tiebreak, tolcon, tolpiv] = parsepvpairs(names, values, varargin{:});
wolffd@0: 
wolffd@0:    catch E
wolffd@0:       E.throw
wolffd@0:    end
wolffd@0: else
wolffd@0:    maxiter = 100000;
wolffd@0:    tiebreak = 'first';
wolffd@0:    tolpiv = 1.0e-9;
wolffd@0:    tolcon = 1.0e-6;
wolffd@0: end
wolffd@0: 
wolffd@0: % check arguments
wolffd@0: 
wolffd@0: if maxiter <= 0
wolffd@0:    error('Finance:qplcprog:NonPositiveInteger', ...
wolffd@0:       'Maximum number of iterations (maxiter) must be a positive integer.');
wolffd@0: end
wolffd@0: 
wolffd@0: if tolcon < 2*eps
wolffd@0:    error('Finance:qplcprog:InvalidTolerance', ...
wolffd@0:       'Unrealistically small constraint tolerance (tolcon) specified.');
wolffd@0: end
wolffd@0: 
wolffd@0: if tolpiv < 2*eps
wolffd@0:    error('Finance:qplcprog:InvalidTolerance', ...
wolffd@0:       'Unrealistically small pivot tolerance (tolpiv) specified.');
wolffd@0: end
wolffd@0: 
wolffd@0: % set translation if necessary
wolffd@0: 
wolffd@0: if norm(lb) ~= 0
wolffd@0:    fx = f + H*lb;
wolffd@0:    if ~isempty(Ai) && ~isempty(bi)
wolffd@0:       bxi = bi - Ai*lb;
wolffd@0:    end
wolffd@0:    if ~isempty(Ae) && ~isempty(be)
wolffd@0:       bxe = be - Ae*lb;
wolffd@0:    end
wolffd@0:    if ~isempty(ub)
wolffd@0:       uxb = ub - lb;
wolffd@0:    end
wolffd@0:    translate = true;
wolffd@0: else
wolffd@0:    fx = f;
wolffd@0:    if ~isempty(Ai) && ~isempty(bi)
wolffd@0:       bxi = bi;
wolffd@0:    end
wolffd@0:    if ~isempty(Ae) && ~isempty(be)
wolffd@0:       bxe = be;
wolffd@0:    end
wolffd@0:    if ~isempty(ub)
wolffd@0:       uxb = ub;
wolffd@0:    end
wolffd@0:    translate = false;
wolffd@0: end
wolffd@0: 
wolffd@0: % set up constraints
wolffd@0: 
wolffd@0: if ~isempty(Ai) && ~isempty(bi)
wolffd@0:    if ~isempty(Ae) && ~isempty(be)
wolffd@0:       if ~isempty(ub)
wolffd@0:          A = [ Ai; Ae; -Ae; eye(n,n) ];
wolffd@0:          b = [ bxi; bxe; -bxe; uxb ];
wolffd@0:       else
wolffd@0:          A = [ Ai; Ae; -Ae ];
wolffd@0:          b = [ bxi; bxe; -bxe ];
wolffd@0:       end
wolffd@0:    else
wolffd@0:       if ~isempty(ub)
wolffd@0:          A = [ Ai; eye(n,n) ];
wolffd@0:          b = [ bxi; uxb ];
wolffd@0:       else
wolffd@0:          A = Ai;
wolffd@0:          b = bxi;
wolffd@0:       end
wolffd@0:    end
wolffd@0: else
wolffd@0:    if ~isempty(Ae) && ~isempty(be)
wolffd@0:       if ~isempty(ub)
wolffd@0:          A = [ Ae; -Ae; eye(n,n) ];
wolffd@0:          b = [ bxe; -bxe; uxb ];
wolffd@0:       else
wolffd@0:          A = [ Ae; -Ae ];
wolffd@0:          b = [ bxe; -bxe ];
wolffd@0:       end
wolffd@0:    else
wolffd@0:       if ~isempty(ub)
wolffd@0:          A = eye(n,n);
wolffd@0:          b = uxb;
wolffd@0:       else
wolffd@0:          error('Finance:qplcprog:InvalidConstraints', ...
wolffd@0:             'Insufficient or invalid constraints.');
wolffd@0:       end
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: % set up lcp problem
wolffd@0: 
wolffd@0: ns = size(A,1);
wolffd@0: 
wolffd@0: M = [ zeros(ns,ns) -A; A' H ];
wolffd@0: q = [ b; fx ];
wolffd@0: 
wolffd@0: % solve lcp problem
wolffd@0: 
wolffd@0: if isempty(varargin)
wolffd@0:    [z, w, status] = lcprog(M, q); %#ok
wolffd@0: else
wolffd@0:    [z, w, status] = lcprog(M, q, 'maxiter', maxiter, 'tiebreak', tiebreak, 'tolpiv', tolpiv); %#ok
wolffd@0: end
wolffd@0: 
wolffd@0: % solve qp problem
wolffd@0: 
wolffd@0: if status > 0
wolffd@0:    x = z(end-n+1:end);
wolffd@0:    if translate
wolffd@0:       x = x + lb;
wolffd@0:    end
wolffd@0:    fval = 0.5*x'*H*x + f'*x;
wolffd@0: 
wolffd@0:    if qplcprogrealitycheck(Ai, bi, Ae, be, lb, ub, x)
wolffd@0:       status = -2;
wolffd@0:       x = NaN(n,1);
wolffd@0:       fval = NaN;
wolffd@0:    end
wolffd@0: else
wolffd@0:    x = NaN(n,1);
wolffd@0:    fval = NaN;
wolffd@0: end
wolffd@0: 
wolffd@0: function notok = qplcprogrealitycheck(Ai, bi, Ae, be, lb, ub, x, tolcon)
wolffd@0: %QPLCPROGREALITYCHECK - Make sure candidate solution x is feasible.
wolffd@0: 
wolffd@0: if nargin < 7
wolffd@0:    error('Finance:qplcprog:MissingInputArgument', ...
wolffd@0:       'Missing required input arguments Ai, bi, Ae, be, lb, ub, x.');
wolffd@0: end
wolffd@0: if nargin < 8 || isempty(tolcon)
wolffd@0:    tolcon = 1.0e-6;
wolffd@0: end
wolffd@0: 
wolffd@0: notok = false;
wolffd@0: 
wolffd@0: if ~isempty(Ai)
wolffd@0:    u = max(Ai*x - bi,0);
wolffd@0:    if u > tolcon
wolffd@0:       warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0:          'Candidate solution violates inequality constraints.');
wolffd@0:       notok = true;
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: if ~isempty(Ae)
wolffd@0:    if norm(be - Ae*x,inf) > tolcon
wolffd@0:       warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0:          'Candidate solution violates equality constraints.');
wolffd@0:       notok = true;
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: if ~isempty(lb)
wolffd@0:    u = x - lb;
wolffd@0:    if min(u) < -tolcon
wolffd@0:       warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0:          'Candidate solution violates lower-bound constraints.');
wolffd@0:       notok = true;
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: if ~isempty(ub)
wolffd@0:    u = ub - x;
wolffd@0:    if min(u) < -tolcon
wolffd@0:       warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0:          'Candidate solution violates upper-bound constraints.');
wolffd@0:       notok = true;
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % [EOF]