gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/LPrelax.m
function [x2, duals_rlx, idx_workc, idx_bindc] = LPrelax(a, f, b, nequs, vlb, vub, idx_workc, mpopt); % MATPOWER Version 2.0 % by Deqiang (David) Gan, PSERC Cornell 12/12/97 % Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC) % See http://www.pserc.cornell.edu/ 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); atemp = full(atemp); [x2, duals] = 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); % if verbose, % print out something for the developer % fprintf('\n (%g, %g) \', ... % length(idx_workc), length(idx_bindc)); % end; 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;