fm_flows.m

一个较好的MATLAB潮流程序

function [Fij, Jij, Hij, Fji, Jji, Hji] = fm_flows(type, mu1, mu2)
% FM_FLOWS compute current, active power and apparent power flow
%          in transmission lines and the associated Jacobian and
%          Hessian matrices
%
%[FIJ, JIJ, HIJ, FJI, JJI, HJI] = FM_FLOWS(TYPE,MU1,MU2)
%
% TYPE -> type of flows computations: (1) current flows
%                                     (2) active power flows
%                                     (3) apparent power flows
%
% MU1 -> column vector (Line.n,1) of Lagrangian multipliers for (i->j)
% MU2 -> column vector (Line.n,1) of Lagrangian multipliers for (j->i)
%
% FIJ -> column vector (Line.n,1) of flows from bus "i" to bus "j"
% JIJ -> Jacobian matrix (Line.n,Bus.n) from bus "i" to bus "j"
% HIJ -> Hessian matrix (Bus.n,Bus.n) from bus "i" to bus "j"
%
% FJI -> column vector (Line.n,1) of flows from bus "j" to bus "i"
% JJI -> Jacobian matrix (Line.n,Bus.n) from bus "j" to bus "i"
% HJI -> Hessian matrix (Bus.n,Bus.n) from bus "j" to bus "i"
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global Line Bus DAE jay

% ==========================================================================
% line flow solution
% ==========================================================================

tps = Line.con(:,11).*exp(jay*Line.con(:,12)*pi/180);
VV = DAE.V.*exp(jay*DAE.a);
r = Line.con(:,8);
rx = Line.con(:,9);
chrg = Line.con(:,10)/2;
z = r + jay*rx;
y = 1./z;
g = real(y);
b = imag(y);
nl = [1:Line.n];

switch type
case 1

    % ==========================================================================
    % AC Current Flows, Jacobian and Hessian
    % ==========================================================================

    Fij = ((VV(Line.from) - tps.*VV(Line.to)).*y + ...
        VV(Line.from).*(jay*chrg))./(tps.*conj(tps));
    Fij = abs(Fij).^2;
    Fji = (VV(Line.to) - VV(Line.from)./tps).*y + VV(Line.to).*(jay*chrg);
    Fji = abs(Fji).^2;

    ati = imag(tps);
    atr = real(tps);
    cij = cos(DAE.a(Line.from)-DAE.a(Line.to));
    sij = sin(DAE.a(Line.from)-DAE.a(Line.to));

    g2 = g.^2;
    b2 = b.^2;
    at1 = ati.^2+atr.^2;
    at2 = at1.*at1;
    Vij = DAE.V(Line.from).*DAE.V(Line.to);
    Vi = DAE.V(Line.from);
    Vj = DAE.V(Line.to);

    a1 = ati.*b.*chrg;                  a2 = ati.*g2;
    a3 = atr.*g.*chrg;                  a4 = atr.*g2;
    a5 = ati.*g.*chrg;                  a6 = atr.*b2;
    a7 = atr.*b.*chrg;                  a8 = ati.*b2;
    c1 = a4+a5+a6+a7;                   c2 = a1+a2-a3+a8;
    c3 = a1+a2+a3+a8;                   c4 = a4-a5+a6+a7;
    d1 = 2*((chrg+b).^2+g2)./at2;       d2 = 2*at1.*(b2+g2)./at2;
    d4 = 2*((chrg+b).^2+g2)./at1;       d3 = 2*(g2+b2)./at1;
    e1 = 2*(c2.*sij+c1.*cij)./at2;      e2 = 2*(c2.*cij-c1.*sij)./at2;
    e3 = 2*(c3.*sij+c4.*cij)./at1;      e4 = 2*(c3.*cij-c4.*sij)./at1;

    Jij = [sparse(nl,Line.from,-Vij.*e2,Line.n,Bus.n) + ...
           sparse(nl,Line.to,   Vij.*e2,Line.n,Bus.n), ...
           sparse(nl,Line.from, Vi.*d1 - Vj.*e1,Line.n,Bus.n) + ...
           sparse(nl,Line.to,   Vj.*d2 - Vi.*e1,Line.n,Bus.n)];
    Jji = [sparse(nl,Line.from,-Vij.*e4,Line.n,Bus.n) + ...
           sparse(nl,Line.to,   Vij.*e4,Line.n,Bus.n), ...
           sparse(nl,Line.from, Vi.*d3 - Vj.*e3,Line.n,Bus.n) + ...
           sparse(nl,Line.to,   Vj.*at1.*d4 - Vi.*e3,Line.n,Bus.n)];

    a1 =  mu1.*Vij.*e1;
    a2 = -mu1.*Vj.*e2;
    a3 = -mu1.*Vi.*e2;
    a4 =  mu1.*d1;
    a5 = -mu1.*e1;
    a6 =  mu1.*d2;

    Hij = [sparse(Line.from, Line.from, a1, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,  -a1, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from,-a1, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,   a1, Bus.n, Bus.n), ...
           sparse(Line.from, Line.from, a2, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,   a3, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from,-a2, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,  -a3, Bus.n, Bus.n); ...
           sparse(Line.from, Line.from, a2, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,  -a2, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from, a3, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,  -a3, Bus.n, Bus.n), ...
           sparse(Line.from, Line.from, a4, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,   a5, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from, a5, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,   a6, Bus.n, Bus.n)];

    a1 =  mu2.*Vij.*e3;
    a2 = -mu2.*Vj.*e4;
    a3 = -mu2.*Vi.*e4;
    a4 =  mu2.*d3;
    a5 = -mu2.*e3;
    a6 =  mu2.*d4.*at1;

    Hji = [sparse(Line.from, Line.from, a1, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,  -a1, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from,-a1, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,   a1, Bus.n, Bus.n), ...
           sparse(Line.from, Line.from, a2, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,   a3, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from,-a2, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,  -a3, Bus.n, Bus.n); ...
           sparse(Line.from, Line.from, a2, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,  -a2, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from, a3, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,  -a3, Bus.n, Bus.n), ...
           sparse(Line.from, Line.from, a4, Bus.n, Bus.n) + ...
           sparse(Line.from, Line.to,   a5, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.from, a5, Bus.n, Bus.n) + ...
           sparse(Line.to,   Line.to,   a6, Bus.n, Bus.n)];

case 2

    % ==========================================================================
    % AC active power flows in each line from bus i to bus j
    % ==========================================================================

    Fij = VV(Line.from).*conj(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
        VV(Line.from).*(jay*chrg))./(tps.*conj(tps)));
    Fij = real(Fij).^2;
    Fji = VV(Line.to).*conj((VV(Line.to) - VV(Line.from)./tps).*y ...
                            + VV(Line.to).*(jay*chrg));
    Fji = real(Fji).^2;

    ati = imag(tps);
    atr = real(tps);
    cji = cos(DAE.a(Line.to)-DAE.a(Line.from));
    sji = sin(DAE.a(Line.to)-DAE.a(Line.from));

    g2 = g.^2;                         b2 = b.^2;
    at2 = ati.^2+atr.^2;               at4 = at2.^2;

    a1 = atr.*b.*cji;        a2 = atr.*g.*sji;
    a3 = ati.*b.*sji;        a4 = ati.*g.*cji;
    a5 = atr.*b.*sji;        a6 = atr.*g.*cji;