fm_rlc.m

电力系统的psat

function  fm_rlc(flag)
% FM_RLC define dynamic quasi-stationary phasors
%
% FM_RLC(FLAG)
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%
%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 RLC Bus DAE Settings

% Data Structure RLC.con:
%
%         col #1:      bus1 number
%         col #2:      bus2 number
%         col #3:      frequency base [Hz]
%         col #4:      power base [MVA]
%         col #5:      type   (1) RC
%                             (2) RL
%                             (3) RLC
%         col #6:      R  [p.u.] Resistance
%         col #7:      Xl [p.u.] Inductance
%         col #8:      Xc [p.u.] conductance

Wb = Settings.rad;

is1 = find(RLC.con(:,6) == 1);
is2 = find(RLC.con(:,6) == 2);
is3 = find(RLC.con(:,6) == 3);
bus1 = RLC.bus1;
bus2 = RLC.bus2;

b11 = bus1(is1);   b12 = bus2(is1);
b21 = bus1(is2);   b22 = bus2(is2);
b31 = bus1(is3);   b32 = bus2(is3);

vq = zeros(RLC.n,1);
vd = zeros(RLC.n,1);
iq = zeros(RLC.n,1);
id = zeros(RLC.n,1);

if ~isempty(is1)
    vd(is1) = DAE.x(RLC.vd(is1));
    vq(is1) = DAE.x(RLC.vq(is1));
elseif ~isempty(is2)
    id(is2) = DAE.x(RLC.id(is2));
    iq(is2) = DAE.x(RLC.iq(is2));
else
    vd(is3) = DAE.x(RLC.vd(is3));
    vq(is3) = DAE.x(RLC.vq(is3));
    id(is3) = DAE.x(RLC.id(is3));
    iq(is3) = DAE.x(RLC.iq(is3));
end

V1 = DAE.V(bus1);
V2 = DAE.V(bus2);
t1 = DAE.a(bus1);
t2 = DAE.a(bus2);

c1 = cos(t1);
c2 = cos(t2);
s1 = sin(t1);
s2 = sin(t2);
c12 = cos(t1-t2);
s12 = sin(t1-t2);

R  = RLC.con(:,7);
xl = RLC.con(:,8);
xc = RLC.con(:,9);

% Algebraic Equations:

switch flag
case 1

	for i = 1:length(is1)
        h = is1(i);
        DAE.gp(b11(i)) = DAE.gp(b11(i)) + (V1(h)^2-V1(h)*c1(h)*vd(h)-V1(h)*V2(h)*c12(h)-V1(h)*s1(h)*vq(h))/R(h);
        DAE.gq(b11(i)) = DAE.gq(b11(i)) + (-V1(h)*s1(h)*vd(h)-V1(h)*V2(h)*s12(h)+V1(h)*c1(h)*vq(h))/R(h);
        DAE.gp(b12(i)) = DAE.gp(b12(i)) + (-V1(h)*V2(h)*c12(h)+V2(h)*c2(h)*vd(h)+V2(h)^2+V2(h)*s2(h)*vq(h))/R(h);
        DAE.gq(b12(i)) = DAE.gq(b12(i)) + (V1(h)*V2(h)*s12(h)+V2(h)*s2(h)*vd(h)-V2(h)*c2(h)*vq(h))/R(h);
    end
    for i = 1:length(is2)
        h = is2(i);
        DAE.gp(b21(i)) = DAE.gp(b21(i)) + (V1(h)*V2(h)*s12(h)-V1(h)*c1(h)*R(h)*iq(h)+V1(h)*s1(h)*R(h)*id(h))/xl(h);
        DAE.gq(b21(i)) = DAE.gq(b21(i)) + (V1(h)^2-V1(h)*V2(h)*c12(h)-V1(h)*s1(h)*R(h)*iq(h)-V1(h)*c1(h)*R(h)*id(h))/xl(h);
        DAE.gp(b22(i)) = DAE.gp(b22(i)) + (-V1(h)*V2(h)*s12(h)+V2(h)*c2(h)*R(h)*iq(h)-V2(h)*s2(h)*R(h)*id(h))/xl(h);
        DAE.gq(b22(i)) = DAE.gq(b22(i)) + (-V1(h)*V2(h)*c12(h)+V2(h)^2+V2(h)*s2(h)*R(h)*iq(h)+V2(h)*c2(h)*R(h)*id(h))/xl(h);
    end
    for i = 1:length(is3)
        h = is3(i);

        RLC.dat(i,1) = (V1(h)*V2(h)*s12(h)-V1(h)*c1(h)*R(h)*id(h)+V1(h)*s1(h)*R(h)*iq(h))/(xl(h)-xc(h));
        RLC.dat(i,2) = (V1(h)^2-V1(h)*V2(h)*c12(h)-V1(h)*s1(h)*R(h)*id(h)-V1(h)*c1(h)*R(h)*iq(h))/(xl(h)-xc(h));
        RLC.dat(i,3) = (-V1(h)*V2(h)*s12(h)+V2(h)*c2(h)*R(h)*id(h)-V2(h)*s2(h)*R(h)*iq(h))/(xl(h)-xc(h));
        RLC.dat(i,4) = (-V1(h)*V2(h)*c12(h)+V2(h)^2+V2(h)*s2(h)*R(h)*id(h)+V2(h)*c2(h)*R(h)*iq(h))/(xl(h)-xc(h));

        %RLC.dat(i,1) = -V1(h)*s1(h)*id(h) + V1(h)*c1(h)*iq(h);
        %RLC.dat(i,2) =  V1(h)*c1(h)*id(h) + V1(h)*s1(h)*iq(h);
        %RLC.dat(i,3) =  V2(h)*s2(h)*id(h) - V2(h)*c2(h)*iq(h);
        %RLC.dat(i,4) = -V2(h)*c2(h)*id(h) - V2(h)*s2(h)*iq(h);

        DAE.gp(b31(i)) = DAE.gp(b31(i)) + RLC.dat(i,1);
        DAE.gq(b31(i)) = DAE.gq(b31(i)) + RLC.dat(i,2);
        DAE.gp(b32(i)) = DAE.gp(b32(i)) + RLC.dat(i,3);
        DAE.gq(b32(i)) = DAE.gq(b32(i)) + RLC.dat(i,4);
    end

% Algebraic Jacobians:

case 2

    for i = 1:length(is1)
        k = is1(i);
        h = b11(i);
        j = b12(i);
        DAE.J12(h,h) = DAE.J12(h,h) + (2*V1(k)-c1(k)*vd(k)-V2(k)*c12(k)-s1(k)*vq(k))/R(k);
        DAE.J12(h,j) = DAE.J12(h,j) + -V1(k)*c12(k)/R(k);
        DAE.J11(h,h) = DAE.J11(h,h) + (V1(k)*s1(k)*vd(k)+V1(k)*V2(k)*s12(k)-V1(k)*c1(k)*vq(k))/R(k);
        DAE.J11(h,j) = DAE.J11(h,j) + -V1(k)*V2(k)*s12(k)/R(k);

        DAE.J22(h,h) = DAE.J22(h,h) + (-s1(k)*vd(k)-V2(k)*s12(k)+c1(k)*vq(k))/R(k);
        DAE.J22(h,j) = DAE.J22(h,j) +  -V1(k)*s12(k)/R(k);
        DAE.J21(h,h) = DAE.J21(h,h) + (-V1(k)*c1(k)*vd(k)-V1(k)*V2(k)*c12(k)-V1(k)*s1(k)*vq(k))/R(k);
        DAE.J21(h,j) = DAE.J21(h,j) + V1(k)*V2(k)*c12(k)/R(k);

        DAE.J12(j,h) = DAE.J12(j,h) + -V2(k)*c12(k)/R(k);
        DAE.J12(j,j) = DAE.J12(j,j) + (-V1(k)*c12(k)+c2(k)*vd(k)+2*V2(k)+s2(k)*vq(k))/R(k);
        DAE.J11(j,h) = DAE.J11(j,h) + V1(k)*V2(k)*s12(k)/R(k);
        DAE.J11(j,j) = DAE.J11(j,j) + (-V1(k)*V2(k)*s12(k)-V2(k)*s2(k)*vd(k)+V2(k)*c2(k)*vq(k))/R(k);

        DAE.J22(j,h) = DAE.J22(j,h) + V2(k)*s12(k)/R(k);
        DAE.J22(j,j) = DAE.J22(j,j) + (V1(k)*s12(k)+s2(k)*vd(k)-c2(k)*vq(k))/R(k);
        DAE.J21(j,h) = DAE.J21(j,h) + V1(k)*V2(k)*c12(k)/R(k);
        DAE.J21(j,j) = DAE.J21(j,j) + (-V1(k)*V2(k)*c12(k)+V2(k)*c2(k)*vd(k)+V2(k)*s2(k)*vq(k))/R(k);
    end
    for i = 1:length(is2)
        k = is2(i);
        h = b21(i);
        j = b22(i);
        DAE.J12(h,h) = DAE.J12(h,h) + (V2(k)*s12(k)-c1(k)*R(k)*iq(k)+s1(k)*R(k)*id(k))/xl(k);
        DAE.J11(h,h) = DAE.J11(h,h) + (V1(k)*V2(k)*c12(k)+V1(k)*s1(k)*R(k)*iq+V1(k)*c1(k)*R(k)*id(k))/xl(k);
        DAE.J22(h,h) = DAE.J22(h,h) + (2*V1(k)-V2(k)*c12(k)-s1(k)*R(k)*iq(k)-c1(k)*R(k)*id(k))/xl(k);
        DAE.J21(h,h) = DAE.J21(h,h) + (V1(k)*V2(k)*s12(k)-V1(k)*c1(k)*R(k)*iq(k)+V1(k)*s1(k)*R(k)*id(k))/xl(k);

        DAE.J12(h,j) = DAE.J12(h,j) + V1(k)*s12(k)/xl(k);
        DAE.J11(h,j) = DAE.J11(h,j) - V1(k)*V2(k)*c12(k)/xl(k);
        DAE.J22(h,j) = DAE.J22(h,j) - V1(k)*c12(k)/xl(k);
        DAE.J21(h,j) = DAE.J21(h,j) - V1(k)*V2(k)*s12(k)/xl(k);

        DAE.J12(j,j) = DAE.J12(j,j) + (-V1(k)*s12(k)+c2(k)*R(k)*iq(k)-s2(k)*R(k)*id(k))/xl(k);
        DAE.J11(j,j) = DAE.J11(j,j) + (V1(k)*V2(k)*c12(k)-V2(k)*s2(k)*R(k)*iq(k)-V2(k)*c2(k)*R(k)*id(k))/xl(k);
        DAE.J22(j,j) = DAE.J22(j,j) + (-V1(k)*c12(k)+2*V2(k)+s2(k)*R(k)*iq(k)+c2(k)*R(k)*id(k))/xl(k);
        DAE.J21(j,j) = DAE.J21(j,j) + (-V1(k)*V2(k)*s12(k)+V2(k)*c2(k)*R(k)*iq(k)-V2(k)*s2(k)*R(k)*id(k))/xl(k);

        DAE.J12(j,h) = DAE.J12(j,h) - V2(k)*s12(k)/xl(k);
        DAE.J11(j,h) = DAE.J11(j,h) - V1(k)*V2(k)*c12(k)/xl(k);
        DAE.J22(j,h) = DAE.J22(j,h) - V2(k)*c12(k)/xl(k);
        DAE.J21(j,h) = DAE.J21(j,h) + V1(k)*V2(k)*s12(k)/xl(k);
    end
    for i = 1:length(is3)
        k = is3(i);
        h = b31(i);
        j = b32(i);

        %DAE.J12(h,h) = DAE.J12(h,h) - s1(k)*id(k) + c1(k)*iq(k);
        %DAE.J11(h,h) = DAE.J11(h,h) - V1(k)*c1(k)*id(k) - V1(k)*s1(k)*iq(k);
        %DAE.J22(h,h) = DAE.J22(h,h) + c1(k)*id(k) + s1(k)*iq(k);
        %DAE.J21(h,h) = DAE.J21(h,h) - V1(k)*s1(k)*id(k) + V1(k)*c1(k)*iq(k);

        %DAE.J12(j,j) = DAE.J12(j,j) + s2(k)*id(k) - c2(k)*iq(k);
        %DAE.J11(j,j) = DAE.J11(j,j) + V2(k)*c2(k)*id(k) + V2(k)*s2(k)*iq(k);
        %DAE.J22(j,j) = DAE.J22(j,j) - c2(k)*id(k) - s2(k)*iq(k);
        %DAE.J21(j,j) = DAE.J21(j,j) + V2(k)*s2(k)*id(k) - V2(k)*c2(k)*iq(k);

        DAE.J12(h,h) = DAE.J12(h,h) + (V2(k)*s12(k)-c1(k)*R(k)*id(k)+s1(k)*R(k)*iq(k))/(xl(k)-xc(k));
        DAE.J11(h,h) = DAE.J11(h,h) + (V1(k)*V2(k)*c12(k)+V1(k)*s1(k)*R(k)*id+V1(k)*c1(k)*R(k)*iq(k))/(xl(k)-xc(k));
        DAE.J22(h,h) = DAE.J22(h,h) + (2*V1(k)-V2(k)*c12(k)-s1(k)*R(k)*id(k)-c1(k)*R(k)*iq(k))/(xl(k)-xc(k));
        DAE.J21(h,h) = DAE.J21(h,h) + (V1(k)*V2(k)*s12(k)-V1(k)*c1(k)*R(k)*id(k)+V1(k)*s1(k)*R(k)*iq(k))/(xl(k)-xc(k));

        DAE.J12(h,j) = DAE.J12(h,j) + V1(k)*s12(k)/(xl(k)-xc(k));
        DAE.J11(h,j) = DAE.J11(h,j) - V1(k)*V2(k)*c12(k)/(xl(k)-xc(k));
        DAE.J22(h,j) = DAE.J22(h,j) - V1(k)*c12(k)/(xl(k)-xc(k));
        DAE.J21(h,j) = DAE.J21(h,j) - V1(k)*V2(k)*s12(k)/(xl(k)-xc(k));

        DAE.J12(j,j) = DAE.J12(j,j) + (-V1(k)*s12(k)+c2(k)*R(k)*id(k)-s2(k)*R(k)*iq(k))/(xl(k)-xc(k));
        DAE.J11(j,j) = DAE.J11(j,j) + (V1(k)*V2(k)*c12(k)-V2(k)*s2(k)*R(k)*id(k)-V2(k)*c2(k)*R(k)*iq(k))/(xl(k)-xc(k));
        DAE.J22(j,j) = DAE.J22(j,j) + (-V1(k)*c12(k)+2*V2(k)+s2(k)*R(k)*id(k)+c2(k)*R(k)*iq(k))/(xl(k)-xc(k));
        DAE.J21(j,j) = DAE.J21(j,j) + (-V1(k)*V2(k)*s12(k)+V2(k)*c2(k)*R(k)*id(k)-V2(k)*s2(k)*R(k)*iq(k))/(xl(k)-xc(k));

        DAE.J12(j,h) = DAE.J12(j,h) - V2(k)*s12(k)/(xl(k)-xc(k));
        DAE.J11(j,h) = DAE.J11(j,h) - V1(k)*V2(k)*c12(k)/(xl(k)-xc(k));
        DAE.J22(j,h) = DAE.J22(j,h) - V2(k)*c12(k)/(xl(k)-xc(k));
        DAE.J21(j,h) = DAE.J21(j,h) + V1(k)*V2(k)*s12(k)/(xl(k)-xc(k));

    end

% Differential Equations:

case 3

    if ~isempty(is1)
        DAE.f(RLC.vd(is1)) = (vq(is1)+xc(is1).*(V1(is1).*c1(is1)-vd(is1)-V2(is1).*c2(is1))./R(is1))*Wb;
        DAE.f(RLC.vq(is1)) = (-vd(is1)+xc(is1).*(V1(is1).*s1(is1)-vq(is1)-V2(is1).*s2(is1))./R(is1))*Wb;
    end
    if ~isempty(is2)
        DAE.f(RLC.id(is2)) = (iq(is2)+(V1(is2).*c1(is2)-V2(is2).*c2(is2)-R(is2).*id(is2))./xl(is2))*Wb;
        DAE.f(RLC.iq(is2)) = (-id(is2)+(V1(is2).*s1(is2)-V2(is2).*s2(is2)-R(is2).*iq(is2))./xl(is2))*Wb;
    end
    if ~isempty(is3)
        vdl = -V1(is3).*s1(is3)+V2(is3).*s2(is3)-R(is3).*id(is3)-vd(is3);
        vql =  V1(is3).*c1(is3)-V2(is3).*c2(is3)-R(is3).*iq(is3)-vq(is3);
        DAE.f(RLC.id(is3)) = ( iq(is3) + vdl./xl(is3))*Wb;
        DAE.f(RLC.iq(is3)) = (-id(is3) + vql./xl(is3))*Wb;
        DAE.f(RLC.vd(is3)) = ( vq(is3) + xc(is3).*id(is3))*Wb;
        DAE.f(RLC.vq(is3)) = (-vd(is3) + xc(is3).*iq(is3))*Wb;
    end

% Jacobians of state and service variables:

case 4

    if isempty(is1), is1 = 0; end
    if isempty(is2), is2 = 0; end

    for i = 1:RLC.n

        if ~isempty(find(is1 == i))

            DAE.Fx(RLC.vd(i),RLC.vd(i)) = -xc(i)/R(i)*Wb;
            DAE.Fx(RLC.vd(i),RLC.vq(i)) = Wb;
            DAE.Fx(RLC.vq(i),RLC.vd(i)) = -Wb;
            DAE.Fx(RLC.vq(i),RLC.vq(i)) = -xc(i)/R(i)*Wb;

            DAE.Fy(RLC.vd(i),bus1(i) + Bus.n) = xc(i)*c1(i)/R(i)*Wb;
            DAE.Fy(RLC.vd(i),bus2(i) + Bus.n) = -xc(i)*c2(i)/R(i)*Wb;
            DAE.Fy(RLC.vd(i),bus1(i)) = -xc(i)*V1(i)*s1(i)/R(i)*Wb;
            DAE.Fy(RLC.vd(i),bus2(i)) = xc(i)*V2(i)*s2(i)/R(i)*Wb;
            DAE.Fy(RLC.vq(i),bus1(i) + Bus.n) = xc(i)*s1(i)/R(i)*Wb;
            DAE.Fy(RLC.vq(i),bus2(i) + Bus.n) = -xc(i)*s2(i)/R(i)*Wb;
            DAE.Fy(RLC.vq(i),bus1(i)) = xc(i)*V1(i)*c1(i)/R(i)*Wb;
            DAE.Fy(RLC.vq(i),bus2(i)) = -xc(i)*V2(i)*c2(i)/R(i)*Wb;

            DAE.Gx(bus1(i),RLC.vd(i)) = -V1(i)*c1(i)/R(i);
            DAE.Gx(bus1(i),RLC.vq(i)) = -V1(i)*s1(i)/R(i);
            DAE.Gx(bus1(i) + Bus.n,RLC.vq(i)) = V1(i)*c1(i)/R(i);
            DAE.Gx(bus1(i) + Bus.n,RLC.vd(i)) = -V1(i)*s1(i)/R(i);
            DAE.Gx(bus2(i),RLC.vd(i)) = V2(i)*c2(i)/R(i);
            DAE.Gx(bus2(i),RLC.vq(i)) = V2(i)*s2(i)/R(i);
            DAE.Gx(bus2(i) + Bus.n,RLC.vq(i)) = -V2(i)*c2(i)/R(i);
            DAE.Gx(bus2(i) + Bus.n,RLC.vd(i)) = V2(i)*s2(i)/R(i);

        elseif ~isempty(find(is2 == i))

            DAE.Fx(RLC.id(i),RLC.id(i)) = -R(i)/xl(i)*Wb;
            DAE.Fx(RLC.id(i),RLC.iq(i)) = Wb;
            DAE.Fx(RLC.iq(i),RLC.id(i)) = -Wb;
            DAE.Fx(RLC.iq(i),RLC.iq(i)) = -R(i)/xl(i)*Wb;

            DAE.Fy(RLC.id(i),bus1(i)) = -V1(i)*s1(i)/xl(i)*Wb;
            DAE.Fy(RLC.id(i),bus2(i)) = V2(i)*s2(i)/xl(i)*Wb;
            DAE.Fy(RLC.id(i),bus1(i) + Bus.n) = c1(i)/xl(i)*Wb;
            DAE.Fy(RLC.id(i),bus2(i) + Bus.n) = -c2(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus1(i)) = V1(i)*c1(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus2(i)) = -V2(i)*c2(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus1(i) + Bus.n) = s1(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus2(i) + Bus.n) = -s2(i)/xl(i)*Wb;

            DAE.Gx(bus1(i),RLC.id(i)) = V1(i)*s1(i)*R(i)/xl(i);
            DAE.Gx(bus1(i),RLC.iq(i)) = -V1(i)*c1(i)*R(i)/xl(i);
            DAE.Gx(bus1(i) + Bus.n,RLC.iq(i)) = -V1(i)*s1(i)*R(i)/xl(i);
            DAE.Gx(bus1(i) + Bus.n,RLC.id(i)) = -V1(i)*c1(i)*R(i)/xl(i);
            DAE.Gx(bus2(i),RLC.id(i)) = -V2(i)*s2(i)*R(i)/xl(i);
            DAE.Gx(bus2(i),RLC.iq(i)) = V2(i)*c2(i)*R(i)/xl(i);
            DAE.Gx(bus2(i) + Bus.n,RLC.iq(i)) = V2(i)*s2(i)*R(i)/xl(i);
            DAE.Gx(bus2(i) + Bus.n,RLC.id(i)) = V2(i)*c2(i)*R(i)/xl(i);

        else % case 3

            DAE.Fx(RLC.id(i),RLC.id(i)) = -R(i)/xl(i)*Wb;
            DAE.Fx(RLC.id(i),RLC.iq(i)) =  Wb;
            DAE.Fx(RLC.id(i),RLC.vd(i)) = -Wb/xl(i);
            DAE.Fx(RLC.iq(i),RLC.id(i)) = -Wb;
            DAE.Fx(RLC.iq(i),RLC.iq(i)) = -R(i)/xl(i)*Wb;
            DAE.Fx(RLC.iq(i),RLC.vq(i)) = -Wb/xl(i);
            DAE.Fx(RLC.vd(i),RLC.vq(i)) =  Wb;
            DAE.Fx(RLC.vq(i),RLC.vd(i)) = -Wb;
            DAE.Fx(RLC.vd(i),RLC.id(i)) = xc(i)*Wb;
            DAE.Fx(RLC.vq(i),RLC.iq(i)) = xc(i)*Wb;

            DAE.Fy(RLC.id(i),bus1(i)) = -V1(i)*c1(i)/xl(i)*Wb;
            DAE.Fy(RLC.id(i),bus2(i)) =  V2(i)*c2(i)/xl(i)*Wb;
            DAE.Fy(RLC.id(i),bus1(i) + Bus.n) = -s1(i)/xl(i)*Wb;
            DAE.Fy(RLC.id(i),bus2(i) + Bus.n) =  s2(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus1(i)) = -V1(i)*s1(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus2(i)) =  V2(i)*s2(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus1(i) + Bus.n) =  c1(i)/xl(i)*Wb;
            DAE.Fy(RLC.iq(i),bus2(i) + Bus.n) = -c2(i)/xl(i)*Wb;

            %DAE.Gx(bus1(i),RLC.iq(i)) = V1(i)*c1(i);
            %DAE.Gx(bus1(i),RLC.id(i)) = -V1(i)*s1(i);
            %DAE.Gx(bus1(i) + Bus.n,RLC.id(i)) = V1(i)*c1(i);
            %DAE.Gx(bus1(i) + Bus.n,RLC.iq(i)) = V1(i)*s1(i);
            %DAE.Gx(bus2(i),RLC.iq(i)) = -V2(i)*c2(i);
            %DAE.Gx(bus2(i),RLC.id(i)) = V2(i)*s2(i);
            %DAE.Gx(bus2(i) + Bus.n,RLC.id(i)) = -V2(i)*c2(i);
            %DAE.Gx(bus2(i) + Bus.n,RLC.iq(i)) = -V2(i)*s2(i);

            DAE.Gx(bus1(i),RLC.iq(i)) = V1(i)*s1(i)*R(i)/(xl(i)-xc(i));
            DAE.Gx(bus1(i),RLC.id(i)) = -V1(i)*c1(i)*R(i)/(xl(i)-xc(i));
            DAE.Gx(bus1(i) + Bus.n,RLC.id(i)) = -V1(i)*s1(i)*R(i)/(xl(i)-xc(i));
            DAE.Gx(bus1(i) + Bus.n,RLC.iq(i)) = -V1(i)*c1(i)*R(i)/(xl(i)-xc(i));
            DAE.Gx(bus2(i),RLC.iq(i)) = -V2(i)*s2(i)*R(i)/(xl(i)-xc(i));
            DAE.Gx(bus2(i),RLC.id(i)) = V2(i)*c2(i)*R(i)/(xl(i)-xc(i));
            DAE.Gx(bus2(i) + Bus.n,RLC.id(i)) = V2(i)*s2(i)*R(i)/(xl(i)-xc(i));
            DAE.Gx(bus2(i) + Bus.n,RLC.iq(i)) = V2(i)*c2(i)*R(i)/(xl(i)-xc(i));

        end
    end

end