lprelax.asv

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

function [x2, duals_rlx, idx_workc, idx_bindc] = LPrelax(a, f, b, nequs, vlb, vub, idx_workc, mpopt);

%   MATPOWER
%   $Id: LPrelax.m,v 1.7 2004/09/01 16:17:05 ray Exp $
%   by Deqiang (David) Gan, PSERC Cornell & Zhejiang University
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

%% options
alg     = mpopt(11);
verbose = mpopt(31);

if opf_slvr(alg) == 1
    idx_workc = find(b < 0.001);
end

converged = 0;
while converged == 0

    atemp = a(idx_workc, :);
    btemp = b(idx_workc);

    [x2, duals] = mp_lp(f, atemp, btemp, vlb, vub, [], nequs, -1);

    diffs = b - a * x2;                 % diffs should be normalized by what means? under development
    idx_bindc = find(diffs < 1.0e-8);
    
    idx_violated = find(diffs < -1.0e-8);

    if isempty(idx_violated) == 1
        converged = 1;
    else
        flag = zeros(length(b), 1);         % set up flag from scratch
        flag(idx_workc) = ones(length(idx_workc), 1);   % enforce historical working constraints

        idx_add = find(diffs < 0.001);
        flag(idx_add) = ones(length(idx_add), 1);   % enforce violating constraints

        flag(1:nequs) = ones(nequs, 1);         % enforce original equality constraints
        idx_workc_new = find(flag);

        if length(idx_workc) == length(idx_workc_new)   % safeguard step
            if isempty(find(idx_workc - idx_workc_new)) ==1
                converged = 1;
            end
        end
        idx_workc = idx_workc_new;

    end

end

duals_rlx = zeros(length(b), 1);
duals_rlx(idx_workc) = duals(1:length(btemp));


return;