fm_dynidx.m

这是一个很适合研究和学习用的电力系统仿真软件

function fm_dynidx
% FM_DYNIDX define indices of state variables for components
%           which are not initialized during the power flow
%           analysis
%
% FM_DYNIDX
%
%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-2006 Federico Milano

global DAE Tg Syn Exc Tcsc Svc Upfc Sssc Statcom Jimma Fl Mixload

if Syn.n
  Syn.delta = zeros(Syn.n,1);
  Syn.omega = zeros(Syn.n,1);
  Syn.e1q = zeros(Syn.n,1);
  Syn.e1d = zeros(Syn.n,1);
  Syn.e2q = zeros(Syn.n,1);
  Syn.e2d = zeros(Syn.n,1);
  Syn.psiq = zeros(Syn.n,1);
  Syn.psid = zeros(Syn.n,1);
  for i = 1:Syn.n
    Syn.delta(i) = DAE.n + 1;
    Syn.omega(i) = DAE.n + 2;
    syn_ord = Syn.con(i,5);
    switch syn_ord
     case 2
      DAE.n = DAE.n+2;
     case 3
      Syn.e1q(i) =   DAE.n + 3;
      DAE.n = DAE.n+3;
     case 4
      Syn.e1q(i) =   DAE.n + 3;
      Syn.e1d(i) =   DAE.n + 4;
      DAE.n = DAE.n+4;
     case 5.1
      Syn.e1q(i) =   DAE.n + 3;
      Syn.e1d(i) =   DAE.n + 4;
      Syn.e2d(i) =   DAE.n + 5;
      DAE.n = DAE.n+5;
     case 5.2
      Syn.e1q(i) =   DAE.n + 3;
      Syn.e2q(i) =   DAE.n + 4;
      Syn.e2d(i) =   DAE.n + 5;
      DAE.n = DAE.n+5;
     case 5.3
      Syn.e1q(i) =   DAE.n + 3;
      Syn.psid(i) =   DAE.n + 4;
      Syn.psiq(i) =   DAE.n + 5;
      DAE.n = DAE.n+5;
     case 6
      Syn.e1q(i) =   DAE.n + 3;
      Syn.e1d(i) =   DAE.n + 4;
      Syn.e2q(i) =   DAE.n + 5;
      Syn.e2d(i) =   DAE.n + 6;
      DAE.n = DAE.n+6;
     case 8
      Syn.e1q(i)   = DAE.n + 3;
      Syn.e1d(i)   = DAE.n + 4;
      Syn.e2q(i)   = DAE.n + 5;
      Syn.e2d(i)   = DAE.n + 6;
      Syn.psiq(i)  = DAE.n + 7;
      Syn.psid(i)  = DAE.n + 8;
      DAE.n = DAE.n+8;
    end
  end
  Syn.vf = 2*ones(Syn.n, 1);
  Syn.pm = ones(Syn.n, 1);
end

if Exc.n
  Exc.vm  = zeros(Exc.n,1);
  Exc.vr1 = zeros(Exc.n,1);
  Exc.vr2 = zeros(Exc.n,1);
  Exc.vr3 = zeros(Exc.n,1);
  Exc.vf = zeros(Exc.n,1);
  for i = 1:Exc.n
    switch Exc.con(i,2)
     case 1
      Exc.vm(i)  = DAE.n + 1;
      Exc.vr1(i) = DAE.n + 2;
      Exc.vr2(i) = DAE.n + 3;
      Exc.vf(i) = DAE.n + 4;
      DAE.n = DAE.n + 4;
     case 2
      Exc.vm(i)  = DAE.n + 1;
      Exc.vr1(i) = DAE.n + 2;
      Exc.vr2(i) = DAE.n + 3;
      Exc.vf(i) = DAE.n + 4;
      DAE.n = DAE.n + 4;
     case 3
      Exc.vm(i)  = DAE.n + 1;
      Exc.vr3(i) = DAE.n + 2;
      Exc.vf(i) = DAE.n + 3;
      DAE.n = DAE.n + 3;
    end
  end
  Exc.vrif  = ones(Exc.n,1);
  Exc.vrif0 = ones(Exc.n,1);
end

Tg = dynidx(Tg);

global Oxl
if Oxl.n
  Oxl.v = DAE.n+[1:Oxl.n]';
  DAE.n = DAE.n + Oxl.n;
  Oxl.If = zeros(Oxl.n,1);
end

global Pss
if Pss.n
  Pss.v1 = zeros(Pss.n,1);
  Pss.v2 = zeros(Pss.n,1);
  Pss.v3 = zeros(Pss.n,1);
  Pss.va = zeros(Pss.n,1);
  Pss.vss = zeros(Pss.n,1);
  for i = 1:Pss.n
    switch Pss.con(i,2)
     case 1
      Pss.v1(i) = DAE.n + 1;
      Pss.vss(i) = DAE.n + 2;
      DAE.n = DAE.n + 2;
     case 2
      Pss.v1(i) = DAE.n + 1;
      Pss.v2(i) = DAE.n + 2;
      Pss.v3(i) = DAE.n + 3;
      Pss.vss(i) = DAE.n + 4;
      DAE.n = DAE.n + 4;
     case 3
      Pss.v1(i) = DAE.n + 1;
      Pss.v2(i) = DAE.n + 2;
      Pss.v3(i) = DAE.n + 3;
      Pss.vss(i) = DAE.n + 4;
      DAE.n = DAE.n + 4;
     case 4
      Pss.v1(i) = DAE.n + 1;
      Pss.v2(i) = DAE.n + 2;
      Pss.v3(i) = DAE.n + 3;
      Pss.va(i) = DAE.n + 4;
      Pss.vss(i) = DAE.n + 5;
      DAE.n = DAE.n + 5;
     case 5
      Pss.v1(i) = DAE.n + 1;
      Pss.v2(i) = DAE.n + 2;
      Pss.v3(i) = DAE.n + 3;
      Pss.va(i) = DAE.n + 4;
      Pss.vss(i) = DAE.n + 5;
      DAE.n = DAE.n + 5;
    end
  end
  Pss.Vs = zeros(Pss.n,1);
  Pss.s1 = zeros(Pss.n,1);
end

global Thload
if Thload.n
  Thload.T = DAE.n + [1:2:2*Thload.n]';
  Thload.G = DAE.n + [2:2:2*Thload.n]';
  DAE.n = DAE.n + 2*Thload.n;
end

Fl = dynidx(Fl);
Svc = dynidx(Svc);
Statcom = dynidx(Statcom);
Tcsc = dynidx(Tcsc);
Sssc = dynidx(Sssc);
Upfc = dynidx(Upfc);

global Mass
if Mass.n
  for i = 1:Mass.n
    Mass.delta_HP(i) = DAE.n + 1;
    Mass.omega_HP(i) = DAE.n + 2;
    Mass.delta_IP(i) = DAE.n + 3;
    Mass.omega_IP(i) = DAE.n + 4;
    Mass.delta_LP(i) = DAE.n + 5;
    Mass.omega_LP(i) = DAE.n + 6;
    Mass.delta_EX(i) = DAE.n + 7;
    Mass.omega_EX(i) = DAE.n + 8;
    DAE.n = DAE.n + 8;
  end
end

global SSR
if SSR.n
  for i = SSR.n
    SSR.Id(i)  = DAE.n + 1;
    SSR.Iq(i)  = DAE.n + 2;
    SSR.If(i)  = DAE.n + 3;
    SSR.Edc(i) = DAE.n + 4;
    SSR.Eqc(i) = DAE.n + 5;
    SSR.delta_HP(i) = DAE.n + 6;
    SSR.omega_HP(i) = DAE.n + 7;
    SSR.delta_IP(i) = DAE.n + 8;
    SSR.omega_IP(i) = DAE.n + 9;
    SSR.delta_LP(i) = DAE.n + 10;
    SSR.omega_LP(i) = DAE.n + 11;
    SSR.delta(i) = DAE.n + 12;
    SSR.omega(i) = DAE.n + 13;
    SSR.delta_EX(i) = DAE.n + 14;
    SSR.omega_EX(i) = DAE.n + 15;
    DAE.n = DAE.n + 15;
  end
  SSR.Tm = zeros(SSR.n,1);
  SSR.Efd = zeros(SSR.n,1);
end

global Sofc
if Sofc.n
  Sofc.Ik = zeros(Sofc.n,1);
  Sofc.Vk = zeros(Sofc.n,1);
  Sofc.pH2 = zeros(Sofc.n,1);
  Sofc.pH20 = zeros(Sofc.n,1);
  Sofc.pO2 = zeros(Sofc.n,1);
  Sofc.qH2 = zeros(Sofc.n,1);
  Sofc.m = zeros(Sofc.n,1);
  for i = 1:Sofc.n
    Sofc.Ik(i) = DAE.n + 1;
    Sofc.Vk(i) = DAE.n + 2;
    Sofc.pH2(i) = DAE.n + 3;
    Sofc.pH20(i) = DAE.n + 4;
    Sofc.pO2(i) = DAE.n + 5;
    Sofc.qH2(i) = DAE.n + 6;
    Sofc.m(i) = DAE.n + 7;
    DAE.n = DAE.n + 7;
  end
end

global Cac
if Cac.n
  Cac.q1 = zeros(Cac.n,1);
  for i = 1:Cac.n
    Cac.q1(i) = DAE.n + 1;
    DAE.n = DAE.n + 1;
  end
end

global Cluster
if Cluster.n
  Cluster.Vs = zeros(Cluster.n,1);
  for i = 1:Cluster.n
    Cluster.Vs(i) = DAE.n + 1;
    DAE.n = DAE.n + 1;
  end
end

global Exload
if Exload.n
  Exload.xp = zeros(Exload.n,1);
  Exload.xq = zeros(Exload.n,1);
  for i = 1:Exload.n
    Exload.xp(i) = DAE.n + 1;
    Exload.xq(i) = DAE.n + 2;
    DAE.n = DAE.n + 2;
  end
end

global Wind
if Wind.n
  Wind.vw = zeros(Wind.n,1);
  for i = 1:Wind.n
    Wind.vw(i) = DAE.n + 1;
    DAE.n = DAE.n + 1;
  end
end

global Cswt
if Cswt.n
  Cswt.omega_wr = zeros(Cswt.n,1);
  Cswt.omega_m = zeros(Cswt.n,1);
  Cswt.gamma = zeros(Cswt.n,1);
  Cswt.e1r = zeros(Cswt.n,1);
  Cswt.e1m = zeros(Cswt.n,1);
  for i = 1:Cswt.n
    Cswt.omega_wr(i) = DAE.n + 1;
    Cswt.omega_m(i) = DAE.n + 2;
    Cswt.gamma(i) = DAE.n + 3;
    Cswt.e1r(i) = DAE.n + 4;
    Cswt.e1m(i) = DAE.n + 5;
    DAE.n = DAE.n + 5;
  end
end

global Dfig
if Dfig.n
  Dfig.theta_p = zeros(Dfig.n,1);
  Dfig.omega_m = zeros(Dfig.n,1);
  Dfig.idr = zeros(Dfig.n,1);
  Dfig.iqr = zeros(Dfig.n,1);
  for i = 1:Dfig.n
    Dfig.omega_m(i) = DAE.n + 1;
    Dfig.theta_p(i) = DAE.n + 2;
    Dfig.idr(i) = DAE.n + 3;
    Dfig.iqr(i) = DAE.n + 4;
    DAE.n = DAE.n + 4;
  end
end

global Ddsg
if Ddsg.n
  Ddsg.theta_p = zeros(Ddsg.n,1);
  Ddsg.omega_m = zeros(Ddsg.n,1);
  Ddsg.ids = zeros(Ddsg.n,1);
  Ddsg.iqs = zeros(Ddsg.n,1);
  Ddsg.idc = zeros(Ddsg.n,1);
  for i = 1:Ddsg.n
    Ddsg.omega_m(i) = DAE.n + 1;
    Ddsg.theta_p(i) = DAE.n + 2;
    Ddsg.ids(i) = DAE.n + 3;
    Ddsg.iqs(i) = DAE.n + 4;
    Ddsg.idc(i) = DAE.n + 5;
    DAE.n = DAE.n + 5;
  end
end

global Busfreq
if Busfreq.n
  Busfreq.x = zeros(Busfreq.n,1);
  Busfreq.w = zeros(Busfreq.n,1);
  for i = 1:Busfreq.n
    Busfreq.x(i) = DAE.n + 1;
    Busfreq.w(i) = DAE.n + 2;
    DAE.n = DAE.n + 2;
  end
end

global Pmu
if Pmu.n
  Pmu.vm = zeros(Pmu.n,1);
  Pmu.thetam = zeros(Pmu.n,1);
  for i = 1:Pmu.n
    Pmu.vm(i) = DAE.n + 1;
    Pmu.thetam(i) = DAE.n + 2;
    DAE.n = DAE.n + 2;
  end
end

global Pod
if Pod.n
  Pod.v1 = zeros(Pod.n,1);
  Pod.v2 = zeros(Pod.n,1);
  Pod.v3 = zeros(Pod.n,1);
  Pod.Vs = zeros(Pod.n,1);
  for i = 1:Pod.n
    Pod.v1(i) = DAE.n + 1;
    Pod.v2(i) = DAE.n + 2;
    Pod.v3(i) = DAE.n + 3;
    Pod.Vs(i) = DAE.n + 4;
    DAE.n = DAE.n + 4;
  end
end

Jimma = dynidx(Jimma);
Mixload = dynidx(Mixload);