mp_lp.m

可进行电力系统多节点系统的优化潮流计算

function [xout, lambdaout, howout, success] = mp_lp(f,A,b,VLB,VUB,x0,N,verbosein,skip_bpmpd)
% MP_LP  Linear program solver.
%   Calls bp() from BPMPD package to solve linear program if available.
%   Otherwise uses the LP solver from the Optimization Toolbox. It looks for
%   linprog() first, then lp(). Calling syntax is equivalent to lp() from
%   version 1 of the Optimization Toolbox.

%   MATPOWER
%   $Id: mp_lp.m,v 1.10 2005/02/03 18:58:25 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   and Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2005 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

success = 0;
if nargin < 9
	skip_bpmpd = 0;
end
if skip_bpmpd | ~have_fcn('bpmpd')  %% no bpmpd, use LP solver from Opt Tbx
    if have_fcn('linprog')          %% use linprog 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
            lpopt = optimset('Display', 'iter');
        elseif verbosein == 1
            lpopt = optimset('Display', 'final');
        else
            lpopt = optimset('Display', 'off');
        end
        [xout,fval,exitflag,output,lambda] = linprog(f,Aieq,bieq,Aeq,beq,VLB,VUB,x0,lpopt);
        howout = exitflag;
        if exitflag == 1
            success = 1;
        end
        lambdaout = [   lambda.eqlin;
                        lambda.ineqlin;
                        lambda.lower;
                        lambda.upper    ];
    elseif have_fcn('lp')   %% use lp from Opt Tbx ver 1.x/2.x
        [xout, lambdaout, howout] = lp(f,A,b,VLB,VUB,x0,N,verbosein);
        if strcmp(howout, 'ok')
            success = 1;
        end
    else
        error('This function requires an LP solver. Please install the Optimization Toolbox or the BPMPD solver.');
    end
else            %% use bpmpd
    n = length(f);
    m = length(b);
    
    if nargin < 8
        verbosein = 0;
        if nargin < 7
            N = 0;
            if nargin < 6
                x0 = zeros(n,1);    % Until bpmpd features warm start, this is a dummy arg
                if nargin < 5
                    VUB = [];
                    if nargin < 4
                        VLB = [];
                    end
                end
            end
        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 = find(VLB > -1e15); % interpret limits <= -1e15 as unbounded
        lval = VLB(llist);
    else
        llist = [];
        lval = [];
    end
    
    if ~isempty(VUB)
        ulist = find(VUB < 1e15); % interpret limits >= 1e15 as unbounded
        uval = VUB(ulist);
    else
        ulist = [];
        uval = [];
    end
    
    if verbosein == -1
        prnlev = 0;
    else
        prnlev = 1;
    end
    if strcmp(computer, 'PCWIN')
        if prnlev
            fprintf('Windows version of BPMPD_MEX cannot print to screen.\n');
        end
        prnlev = 0;   % The DLL incarnation of bp was born mute and deaf,
    end               % probably because of acute shock after realizing its fate.
                      % Can't be allowed to try to speak or its universe crumbles.
    
    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([], 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
    if success
        % double-check feasibility
		bpmpd_bug_fatal = 0;
       	err_tol = 1e-4;
        nb = length(b);
        if ~isempty(VLB)
        	lb_violation = VLB - xout;
			if verbosein > 0
				fprintf('max variable lower bound violatation: %g\n', max(lb_violation));
			end
        else
        	lb_violation = zeros(size(xout));
        end
        if ~isempty(VUB)
            ub_violation = xout - VUB;
			if verbosein > 0
				fprintf('max variable upper bound violation: %g\n', max(ub_violation));
			end
        else
        	ub_violation = zeros(size(xout));
        end
        if N > 0
            eq_violation = abs( A(1:N,:) * xout - b(1:N) );
			if verbosein > 0
				fprintf('max equality constraint violation: %g\n', max(eq_violation));
			end
        else
        	eq_violation = zeros(N,1);
        end
        if N < nb
            ineq_violation = A(N+1:nb,:) * xout - b(N+1:nb);
			if verbosein > 0
				fprintf('max inequality constraint violation: %g\n', max(ineq_violation));
			end
        else
        	ineq_violation = zeros(nb-N,1);
        end

		if any( [ 	lb_violation;
					ub_violation;
					eq_violation;
					ineq_violation ] > err_tol)
			err_cnt = 0;
			if any( lb_violation > err_tol )
				err_cnt = err_cnt + 1;
				errs{err_cnt} = ...
					sprintf('variable lower bound violated by %g', ...
						max(lb_violation));
			end
			if any( ub_violation > err_tol )
				err_cnt = err_cnt + 1;
				errs{err_cnt} = ... 
					sprintf('variable upper bound violated by %g', ...
						max(ub_violation));
			end
			if any( eq_violation > err_tol )
				err_cnt = err_cnt + 1;
				errs{err_cnt} = ... 
					sprintf('equality constraint violated by %g', ...
						max(eq_violation));
			end
			if any( ineq_violation > err_tol )
				err_cnt = err_cnt + 1;
				errs{err_cnt} = ... 
					sprintf('inequality constraint violated by %g', ...
						max(ineq_violation));
			end
			fprintf('\nWARNING: This version of BPMPD_MEX has a bug which caused it to return\n');
			fprintf(  '         an incorrect (infeasible) solution for this particular problem.\n');
			for err_idx = 1:err_cnt
				fprintf('         %s\n', errs{err_idx});
			end
			if bpmpd_bug_fatal
				error(sprintf('%s\n%s', ...
					'To avoid this bug in BPMPD_MEX you will need', ...
					'to use a different QP solver for this problem.'));
			else
				fprintf('         Retrying with LP solver from Optimization Toolbox ...\n\n');
				skip_bpmpd = 1;		%% try again using another solver
				[xout, lambdaout, howout, success] = ...
					mp_lp(f,A,b,VLB,VUB,x0,N,verbosein,skip_bpmpd);
			end
		end
    end
end

return;