mp_qp.m

以上所传的内容为电力系统比较全面的分析程序,才用的是matlab编写m的文件

function [xout, lambdaout, howout, success] = mp_qp(H,f,A,b,VLB,VUB,x0,N,verbosein)
% MP_QP  Quadratic program solver.
%   Calls bp() from BPMPD package to solve quadratic program if available.
%   Otherwise uses the QP solver from the Optimization Toolbox. It looks for
%   quadprog() first, then qp(). Calling syntax is equivalent to qp() from
%   version 1 of the Optimization Toolbox.

%   MATPOWER
%   $Id: mp_qp.m,v 1.5 2004/09/01 17:42:51 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   and Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

success = 0;
if ~have_fcn('bpmpd')       %% bpmpd not available, use QP solver from Opt Tbx
    if have_fcn('quadprog') %% use quadprog from Opt Tbx ver 2.x+
        if N
            nb = length(b);
            Aeq = A(1:N,:);
            beq = b(1:N);
            Aieq = A(N+1:nb,:);
            bieq = b(N+1:nb);
        else
            Aeq = [];
            beq = [];
            Aieq = A;
            bieq = b;
        end
        if verbosein > 1
            qpopt = optimset('Display', 'iter');
        elseif verbosein == 1
            qpopt = optimset('Display', 'final');
        else
            qpopt = optimset('Display', 'off');
        end
        [xout,fval,exitflag,output,lambda] = quadprog(H,f,Aieq,bieq,Aeq,beq,VLB,VUB,x0,qpopt);
        howout = exitflag;
        if exitflag == 1
            success = 1;
        end
        lambdaout = [   lambda.eqlin;
                        lambda.ineqlin;
                        lambda.lower;
                        lambda.upper    ];
    elseif have_fcn('qp')   %% use qp from Opt Tbx ver 1.x/2.x
        [xout, lambdaout, howout] = qp(H,f,A,b,VLB,VUB,x0,N,verbosein);
        if strcmp(howout, 'ok')
            success = 1;
        end
    else
        error('This function requires a QP solver. Please install the Optimization Toolbox or the BPMPD solver.');
    end
else            %% use bpmpd
    n = length(f);
    m = length(b);
    
    if nargin < 9
        verbosein = 0;
        if nargin < 8
            N = 0;
            if nargin < 7
                x0 = zeros(n,1);    % Until bpmpd features warm start, this is a dummy arg
                if nargin < 6
                    VUB = [];
                    if nargin < 5
                        VLB = [];
                    end
                end
            end
        end
    end
    
    if ~isempty(H)
        if ~issparse(H)
            H = sparse(H);
        end
    end
    
    if ~issparse(A)
        A = sparse(A);
    end
    
    e = -ones(m,1);
    if N>0
        e(1:N,:) = zeros(N,1);
    end
    
    if ~isempty(VLB)
        llist = 1:n;
        lval = VLB;
    else
        llist = [];
        lval = [];
    end
    
    if ~isempty(VUB)
        ulist = 1:n;
        uval = VUB;
    else
        ulist = [];
        uval = [];
    end
    
    if verbosein == -1
        prnlev = 0;
    else
        prnlev = 1;
    end
    
    myopt = bpopt;
    %myopt(14)= 1e-1;   % TPIV1  first relative pivot tolerance (desired)
    %myopt(20)= 1e-10;  % TOPT1  stop if feasible and rel. dual gap less than this
    %myopt(22)= 1e-6;   % TFEAS1 relative primal feasibility tolerance
    myopt(23)= 1e-6;    % TFEAS2 relative dual feasibility tolerance
    %myopt(29)= 1e-10;  % TRESX  acceptable primal residual
    %myopt(30)= 1e-10;  % TRESY  acceptable dual residual
    %myopt(38)= 0;      % SMETHOD1 prescaling method
    
    [xout,y,s,w,howout] = bp(H, A, b, f, e, llist, lval, ...
                            ulist, uval, myopt, prnlev);
    
    ilt = find(w<=0);
    igt = find(w>0);
    mulow = zeros(n,1);
    muupp = zeros(n,1);
    muupp(ilt) = -w(ilt);
    mulow(igt) = w(igt);
    
    lambdaout = -y;
    if ~isempty(VLB)
        lambdaout = [lambdaout; mulow];
    end
    if ~isempty(VUB)
        lambdaout = [lambdaout; muupp];
    end
    
    % zero out lambdas smaller than a certain tolerance
    ii = find(abs(lambdaout)<1e-9);
    lambdaout(ii) = zeros(size(ii));
    
    % The next is necessary for proper operation of constr.m
    if strcmp(howout, 'infeasible primal')
        lambdaout = zeros(size(lambdaout));
    end
    if strcmp(howout, 'optimal solution')
        success = 1;
    end
end

return;