fm_hvdc.m

电力系统的psat

function  fm_hvdc(flag)
% FM_HVDC define HVDC transmission system
%
% FM_HVDC(FLAG)
%       FLAG = 0 initialization
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%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 Hvdc Bus DAE

% Data Structure: Hvdc.dat:
% col #1:  cosar
% col #2:  Vdr
% col #3:  Sr
% col #4:  Vdi
% col #5:  Si
% col #6:  cosgi
% col #7:  I0r
% col #8:  I0i
% col #9:  1/Rd
% col #10: Ld/Rd
% col #11: cosar_max
% col #12: cosar_min
% col #13: cosgi_max
% col #14: cosgi_min
% col #15: 1/KI
% col #16: Vn*In/Sn

Id = DAE.x(Hvdc.Id);
xr = DAE.x(Hvdc.xr);
xi = DAE.x(Hvdc.xi);

bus1 = Hvdc.bus1;
bus2 = Hvdc.bus2;
V1 = DAE.V(bus1);
V2 = DAE.V(bus2);

xcr = Hvdc.con(:,9);
xci = Hvdc.con(:,10);
ar = Hvdc.con(:,11);
ai = Hvdc.con(:,12);

Kp = Hvdc.con(:,14);

I0r_max = Hvdc.con(:,21);
I0r_min = Hvdc.con(:,22);
I0i_max = Hvdc.con(:,23);
I0i_min = Hvdc.con(:,24);

cosar = Hvdc.dat(:,1);
Vdr = Hvdc.dat(:,2);
Sr = Hvdc.dat(:,3);
Vdi = Hvdc.dat(:,4);
Si = Hvdc.dat(:,5);
cosgi = Hvdc.dat(:,6);
I0r = Hvdc.dat(:,7);
I0i = Hvdc.dat(:,8);
Rd = Hvdc.dat(:,9);
LdonRd = Hvdc.dat(:,10);
cosar_max = Hvdc.dat(:,11);
cosar_min = Hvdc.dat(:,12);
cosgi_max = Hvdc.dat(:,13);
cosgi_min = Hvdc.dat(:,14);
xr_max = Hvdc.dat(:,11);
xr_min = Hvdc.dat(:,12);
xi_max = Hvdc.dat(:,13);
xi_min = Hvdc.dat(:,14);
KIinv = Hvdc.dat(:,15);
VISn = Hvdc.dat(:,16);

switch flag
 case 0 % Initialization

  pi = 3.14159265358979;
  
  Hvdc.dat(:,1) = cosar;
  Hvdc.dat(:,2) = Vdr;
  Hvdc.dat(:,3) = Sr;
  Hvdc.dat(:,4) = Vdi;
  Hvdc.dat(:,5) = Si;
  Hvdc.dat(:,6) = cosgi;
  Hvdc.dat(:,7) = I0r;
  Hvdc.dat(:,8) = I0i;
  Hvdc.dat(:,9) = 1./Hvdc.con(:,15);
  Hvdc.dat(:,10) = Hvdc.con(:,16)./Hvdc.con(:,15);
  Hvdc.dat(:,11) = cos(pi*Hvdc.con(:,18)/180);
  Hvdc.dat(:,12) = cos(pi*Hvdc.con(:,17)/180);
  Hvdc.dat(:,13) = cos(pi*Hvdc.con(:,20)/180);
  Hvdc.dat(:,14) = cos(pi*Hvdc.con(:,19)/180);
  Hvdc.dat(:,15) = 1./Hvdc.con(:,13);
  Hvdc.dat(:,16) = Hvdc.con(:,7).*Hvdc.con(:,8)./Hvdc.con(:,3);
  
  Vdr = ones(Hvdc.n,1);
  Sr = 1.3505*Hvdc.dat(:,16).*V1;
  Vdi = ones(Hvdc.n,1);
  Si = 1.3505*Hvdc.dat(:,16).*V2;
  
  I0r = I0r_max;
  I0i = I0i_min;
  cosar = (Hvdc.dat(:,11)-Hvdc.dat(:,12))/2;
  cosgi = Hvdc.dat(:,14);
  
 case 1 % algebraic equations
  
  DAE.gp = DAE.gp + sparse(Hvdc.bus1,1,VISn.*Vdr.*Id,Bus.n,1);
  DAE.gp = DAE.gp - sparse(Hvdc.bus2,1,VISn.*Vdi.*Id,Bus.n,1);
  DAE.gq = DAE.gq + sparse(Hvdc.bus1,1,VISn.*Id.*sqrt(1.8238*V1.*V1-Vdr.*Vdr),Bus.n,1);
  DAE.gq = DAE.gq + sparse(Hvdc.bus2,1,VISn.*Id.*sqrt(1.8238*V2.*V2-Vdi.*Vdi),Bus.n,1);
  
 case 2 % algebraic Jacobians
  
  a = 3040988094387517/2251799813685248;
  b = 95493/100000;
  c = 9247608590206622002595365425289;
  d = 5070602400912917605986812821504;
  e = 10141204801825835211973625643008;
  g = 18495217180413244005190730850578;
  h = 1/4503599627370496;
  
  Qi = (c*(VISn.*V2.*Id).^2-d*(VISn.*Vdi.*Id).^2).^0.5;
  Qr = (c*(VISn.*V1.*Id).^2-d*(VISn.*Vdr.*Id).^2).^0.5;
  
  j12_1 = VISn.*(a*xr + a*Kp.*(I0r-Id) + a*V1.*Kp./ar).*Id;
  j22_1 = h./Qr.*(g*(VISn.^2).*V1.*(Id.^2)-e*(VISn.^2)*Vdr.*(Id.^2).*(cosar+a*V1.*Kp./ar));
  j12_2 = -VISn.*(cosgi-a*Kp.*I0i).*Id;
  j22_2 = h./Qi.*(g*(VISn.^2).*V2.*(Id.^2)-e*(VISn.^2).*Vdi.*(Id.^2).*(cosgi-a*Kp.*I0i));
  
  i = find(I0r >= I0r_max | I0r <= I0r_min);
  if i,
    j12_1(i) = VISn(i).*cosar(i).*Id(i);
    j22_1(i) = h./Qr(i).*(g*(VISn(i).^2).*V1(i).*(Id(i).^2) - e*(VISn(i).^2).*Vdr(i).*(Id(i).^2).*cosar(i));
  end
  i = find(I0i >= I0i_max | I0i <= I0i_min);
  if i,
    j12_2(i) = -VISn(i).*cosgi(i).*Id(i);
    j22_2(i) = h./Qi(i).*(g*(VISn(i).^2).*V2(i).*(Id(i).^2) - e*(VISn(i).^2).*Vdi(i).*(Id(i).^2).*cosgi(i));
  end
  i = find(cosar >= cosar_max | cosar <= cosar_min);
  if i,
    j12_1(i) = a*VISn(i).*cosar(i).*Id(i);
    j22_1(i) = h./Qr(i).*(g*(VISn(i).^2).*V1(i).*(Id(i).^2) - 2*a*(VISn(i).^2).*Vdr(i).*(Id(i).^2).*cosar(i));
  end
  i = find(cosgi >= cosgi_max | cosgi <= cosgi_min);
  if i,
    j12_2(i) = -a*VISn(i).*cosgi(i).*Id(i);
    j22_2(i) = h./Qi(i).*(g*(VISn(i).^2).*V2(i).*(Id(i).^2) - 2*a*(VISn(i).^2).*Vdi(i).*(Id(i).^2).*cosgi(i));
  end
  
  DAE.J12 = DAE.J12 + sparse(Hvdc.bus1,Hvdc.bus1,j12_1,Bus.n,Bus.n);
  DAE.J22 = DAE.J12 + sparse(Hvdc.bus1,Hvdc.bus1,j22_1,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(Hvdc.bus2,Hvdc.bus2,j12_2,Bus.n,Bus.n);
  DAE.J22 = DAE.J12 + sparse(Hvdc.bus2,Hvdc.bus2,j22_2,Bus.n,Bus.n);
  
 case 3  % differential equations
  
  I0r = V1./ar;
  I0r = min(I0r,I0r_max);
  I0r = max(I0r,I0r_min);
  
  I0i = V2./ai;
  I0i = min(I0i,I0i_max);
  I0i = max(I0i,I0i_min);
  
  cosar = xr+Kp.*(I0r-Id);
  cosar = min(cosar,cosar_max);
  cosar = max(cosar,cosar_min);
  
  cosgi = xi+Kp.*(Id-I0i);
  cosgi = min(cosgi,cosgi_max);
  cosgi = max(cosgi,cosgi_min);
  
  a = 3040988094387517/2251799813685248;
  b = 95493/100000;
  Vdr = a.*V1.*cosar-b.*xcr.*Id;
  Sr = a.*VISn.*V1.*Id;
  Vdi = a.*V2.*cosgi-b.*xci.*Id;
  Si = a.*VISn.*V2.*Id;
  
  DAE.f(Hvdc.Id) = (Rd.*(Vdr-Vdi)-Id)./LdonRd;
  DAE.f(Hvdc.xr) = (I0r-Id)./KIinv;
  DAE.f(Hvdc.xi) = (Id-I0i)./KIinv;
  
  rectifier = 1;
  for i = 1:Hvdc.n
    if DAE.x(Hvdc.xr(i)) >= xr_max(i)
      DAE.x(Hvdc.xr(i)) = xr_max(i);
      if DAE.f(Hvdc.xr(i)) > 0; DAE.f(Hvdc.xr(i)) = 0; end
    elseif DAE.x(Hvdc.xr(i)) <= xr_min(i)
      DAE.x(Hvdc.xr(i)) = xr_min(i);
      if DAE.f(Hvdc.xr(i)) < 0; DAE.f(Hvdc.xr(i)) = 0; end
    else
      DAE.x(Hvdc.xi(i)) = xi_min(i);
      DAE.f(Hvdc.xi(i)) = 0;
      cosgi(i) = cosgi_min(i);
      Vdi(i) = a*V2(i)*cosgi(i)-b*xci(i)*Id(i);
      Si(i) = a*VISn(i)*V2(i)*Id(i);
      rectifier = 0;
    end
    if rectifier
      if DAE.x(Hvdc.xi(i)) >= xi_max(i)
        DAE.x(Hvdc.xi(i)) = xi_max(i);
        if DAE.f(Hvdc.xi(i)) > 0; DAE.f(Hvdc.xi(i)) = 0; end
      elseif DAE.x(Hvdc.xi(i)) <= xi_min(i)
        DAE.x(Hvdc.xi(i)) = xi_min(i);
        if DAE.f(Hvdc.xi(i)) < 0; DAE.f(Hvdc.xi(i)) = 0; end
      else
        DAE.x(Hvdc.xr(i)) = xr_min(i);
        DAE.f(Hvdc.xi(i)) = 0;
        cosar(i) = cosar_min(i);
        Vdr(i) = a*V1(i)*cosar(i)-b*xcr(i)*Id(i);
        Sr(i) = a*VISn(i)*V1(i)*Id(i);
      end
    end
  end
  
  Hvdc.dat(:,1) = cosar;
  Hvdc.dat(:,2) = Vdr;
  Hvdc.dat(:,3) = Sr;
  Hvdc.dat(:,4) = Vdi;
  Hvdc.dat(:,5) = Si;
  Hvdc.dat(:,6) = cosgi;
  Hvdc.dat(:,7) = I0r;
  Hvdc.dat(:,8) = I0i;
  
 case 4  % Jacobians of state variables
  
  a = 3040988094387517/2251799813685248;
  b = 95493/100000;
  c = 9247608590206622002595365425289;
  d = 5070602400912917605986812821504;
  e = 10141204801825835211973625643008;
  g = 18495217180413244005190730850578;
  h = 1/4503599627370496;
  m = -3040988094387517;
  
  for i = 1:Hvdc.n
    
    Qi = (c*VISn(i)^2*V2(i)^2*Id(i)^2-d*VISn(i)^2*Vdi(i)^2*Id(i)^2)^(1/2);
    Qr = (c*VISn(i)^2*V1(i)^2*Id(i)^2-d*VISn(i)^2*Vdr(i)^2*Id(i)^2)^(1/2);
    
    DAE.Fx(Hvdc.Id(i),Hvdc.Id(i)) = (Rd(i)*(-a*V1(i)*Kp(i)-b*xcr(i)-a*Kp(i)*V2(i)+b*xci(i))-1)/LdonRd(i);
    DAE.Fx(Hvdc.Id(i),Hvdc.xr(i)) = a*Rd(i)*V1(i)/LdonRd(i);
    DAE.Fx(Hvdc.Id(i),Hvdc.xi(i)) = -a*Rd(i)*V2(i)/LdonRd(i);
    DAE.Fx(Hvdc.xr(i),Hvdc.Id(i)) = -1/KIinv(i);
    DAE.Fx(Hvdc.xi(i),Hvdc.Id(i)) = 1/KIinv(i);
    
    DAE.Fy(Hvdc.xr(i),bus1(i) + Bus.n) = 1/ar(i)/KIinv(i);
    DAE.Fy(Hvdc.Id(i),bus1(i) + Bus.n) = Rd(i)*(cosar(i)+a*V1(i)*Kp(i)/ar(i))/LdonRd(i);
    DAE.Fy(Hvdc.xi(i),bus2(i) + Bus.n) = -1/ai(i)/KIinv(i);
    DAE.Fy(Hvdc.Id(i),bus2(i) + Bus.n) = Rd(i)*(-cosgi(i)+a*Kp(i)*I0i(i))/LdonRd(i);
    
    DAE.Gx(bus1(i),Hvdc.Id(i)) = VISn(i)*(-a*V1(i)*Kp(i)-b*xcr(i))*Id(i)+VISn(i)*Vdr(i);
    DAE.Gx(bus1(i),Hvdc.xr(i)) = a*VISn(i)*V1(i)*Id(i);
    
    DAE.Gx(bus1(i) + Bus.n,Hvdc.Id(i)) = h/Qr*(g*VISn(i)^2*V1(i)^2*Id(i)-e*VISn(i)^2*Vdr(i)*Id(i)^2* ...
                                               (-a*V1(i)*Kp(i)-b*xcr(i))-e*VISn(i)^2*Vdr(i)^2*Id(i));
    DAE.Gx(bus1(i) + Bus.n,Hvdc.xr(i)) = m/Qr*VISn(i)^2*Vdr(i)*Id(i)^2*V1(i);
    
    DAE.Gx(bus2(i),Hvdc.Id(i)) = -VISn(i)*(a*Kp(i)*V2(i)-b*xci(i))*Id(i)-VISn(i)*Vdi(i);
    DAE.Gx(bus2(i),Hvdc.xi(i)) = -a*VISn(i)*V2(i)*Id(i);
    
    DAE.Gx(bus2(i) + Bus.n,Hvdc.Id(i)) = h/Qi*(g*VISn(i)^2*V2(i)^2*Id(i)-e*VISn(i)^2*Vdi(i)*Id(i)^2* ...
                                               (a*Kp(i)*V2(i)-b*xci(i))-e*VISn(i)^2*Vdi(i)^2*Id(i));
    DAE.Gx(bus2(i) + Bus.n,Hvdc.xi(i)) = m/Qi*VISn(i)^2*Vdi(i)*Id(i)^2*V2(i);
    
    if I0r(i) >= I0r_max(i) | I0r(i) <= I0r_min(i)
      DAE.Fy(Hvdc.xr(i),bus1(i) + Bus.n) = 0;
      DAE.Fy(Hvdc.Id(i),bus1(i) + Bus.n) = a*cosar(i)*Rd(i)/LdonRd(i);
    end
    
    if I0i(i) >= I0i_max(i) | I0i(i) <= I0i_min(i)
      DAE.Fy(Hvdc.xi(i),bus2(i) + Bus.n) = 0;
      DAE.Fy(Hvdc.Id(i),bus2(i) + Bus.n) = -a*cosgi(i)*Rd(i)/LdonRd(i);
    end
    
    if cosar(i) >= cosar_max(i) | cosar(i) <= cosar_min(i)
      DAE.Fx(Hvdc.Id(i),Hvdc.Id(i)) = (Rd(i)*(-b*xcr(i)-a*Kp(i)*V2(i)+b*xci(i))-1)/LdonRd(i);
      %DAE.Fx(Hvdc.Id(i),Hvdc.xr(i)) = 0;
      DAE.Gx(bus1(i),Hvdc.Id(i)) = VISn(i)*(-b*xcr(i))*Id(i)+VISn(i)*Vdr(i);
      %DAE.Gx(bus1(i),Hvdc.xr(i)) = 0;
      %DAE.Gx(bus1(i) + Bus.n,Hvdc.xr(i)) = 0;
      DAE.Gx(bus1(i) + Bus.n,Hvdc.Id(i)) = h/Qr*(g*VISn(i)^2*V1(i)^2*Id(i) - e*VISn(i)^2*Vdr(i)*Id(i)^2* ...
                                                 (-b*xcr(i)) - e*VISn(i)^2*Vdr(i)^2*Id(i));
      DAE.Fy(Hvdc.Id(i),bus1(i) + Bus.n) = Rd(i)*a*cosar(i)/LdonRd(i);
    end
    
    if cosgi(i) >= cosgi_max(i) | cosgi(i) <= cosgi_min(i)
      DAE.Fx(Hvdc.Id(i),Hvdc.Id(i)) = (Rd(i)*(-a*V1(i)*Kp(i)-b*xcr(i)+b*xci(i))-1)/LdonRd(i);
      %DAE.Fx(Hvdc.Id(i),Hvdc.xi(i)) = 0;
      DAE.Gx(bus2(i),Hvdc.Id(i)) = -VISn(i)*(-b*xci(i))*Id(i)-VISn(i)*Vdi(i);
      %DAE.Gx(bus2(i),Hvdc.xi(i)) = 0;
      %DAE.Gx(bus2(i) + Bus.n,Hvdc.xi(i)) = 0;
      DAE.Gx(bus2(i) + Bus.n,Hvdc.Id(i)) = h/Qi*(g*VISn(i)^2*V2(i)^2*Id(i)-e*VISn(i)^2*Vdi(i)*Id(i)^2* ...
                                                 (-b*xci(i))-e*VISn(i)^2*Vdi(i)^2*Id(i));
      DAE.Fy(Hvdc.Id(i),bus2(i) + Bus.n) = -Rd(i)*a*cosgi(i)/LdonRd(i);
    end
    
    if (DAE.x(Hvdc.xr(i)) >= xr_max(i) | DAE.x(Hvdc.xr(i)) <= xr_min(i)) & (DAE.f(Hvdc.xr(i)) == 0)
      DAE.Fx(Hvdc.Id(i),Hvdc.xr(i)) = 0;
      DAE.Fx(Hvdc.xr(i),Hvdc.Id(i)) = 0;
      DAE.Fx(Hvdc.xr(i),Hvdc.xr(i)) = -1;
      DAE.Fy(Hvdc.xr(i),bus1(i) + Bus.n) = 0;
      DAE.Gx(bus1(i),Hvdc.xr(i)) = 0;
      DAE.Gx(bus1(i) + Bus.n,Hvdc.xr(i)) = 0;
    else
      DAE.Fx(Hvdc.xr(i),Hvdc.xr(i)) = 0;
    end
    if (DAE.x(Hvdc.xi(i)) >= xi_max(i) | DAE.x(Hvdc.xi(i)) <= xi_min(i)) & (DAE.f(Hvdc.xi(i)) == 0)
      DAE.Fx(Hvdc.Id(i),Hvdc.xi(i)) = 0;
      DAE.Fx(Hvdc.xi(i),Hvdc.Id(i)) = 0;
      DAE.Fx(Hvdc.xi(i),Hvdc.xi(i)) = -1;
      DAE.Fy(Hvdc.xi(i),bus2(i) + Bus.n) = 0;
      DAE.Gx(bus2(i),Hvdc.xi(i)) = 0;
      DAE.Gx(bus2(i) + Bus.n,Hvdc.xi(i)) = 0;
    else
      DAE.Fx(Hvdc.xi(i),Hvdc.xi(i)) = 0;
    end
  end
  
 case 5 % non-windup limiters
  
  Ac_idx = 2*Bus.n+DAE.n;
  
  for i = 1:Hvdc.n
    if (DAE.x(Hvdc.xr(i)) >= xr_max(i) | DAE.x(Hvdc.xr(i)) <= xr_min(i)) & (DAE.f(Hvdc.xr(i)) == 0)
      k = Hvdc.xr(i);
      DAE.tn(k) = 0;
      DAE.Ac(:,k) = zeros(Ac_idx,1);
      DAE.Ac(k,:) = zeros(1,Ac_idx);
      DAE.Ac(k,k) = DAE.Fx(k,k);
    end
    if (DAE.x(Hvdc.xi(i)) >= xi_max(i) | DAE.x(Hvdc.xi(i)) <= xi_min(i)) & (DAE.f(Hvdc.xi(i)) == 0)
      k = Hvdc.xi(i);
      DAE.tn(k) = 0;
      DAE.Ac(:,k) = zeros(Ac_idx,1);
      DAE.Ac(k,:) = zeros(1,Ac_idx);
      DAE.Ac(k,k) = DAE.Fx(k,k);
    end
  end
  
end