fm_upfc.m

基于PSAT 软件的多目标最优潮流计算用于中小型电力系统的分析和管理

function  fm_upfc(flag)
% FM_UPFC define Unified Power Flow Controller - UPFC
%
% FM_UPFC(FLAG)
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%Author:    Hugo M. Ayres (revised by Federico Milano)
%Date:      18-Feb-2006
%Version:   1.0.0
%
%E-mail:    [e-mail]
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2006 Hugo M. Ayres & Federico Milano

global Upfc Syn Pod Bus DAE PV Line jay

bus1 = Upfc.bus1;
bus2 = Upfc.bus2;
V1 = DAE.V(bus1);
V2 = DAE.V(bus2);
a1 = DAE.a(bus1);
a2 = DAE.a(bus2);

kp = Upfc.Cp./(1-Upfc.Cp);
V = V1.*exp(jay.*a1)-V2.*exp(jay.*a2);
theta = angle(V)-pi/2;

vp = DAE.x(Upfc.vp);
vq = DAE.x(Upfc.vq);
iq = DAE.x(Upfc.iq);

Kr = Upfc.con(:,7);
Tr = Upfc.con(:,8);

vp_max = Upfc.con(:,9);
vp_min = Upfc.con(:,10);
vq_max = Upfc.con(:,11);
vq_min = Upfc.con(:,12);
iq_max = Upfc.con(:,13);
iq_min = Upfc.con(:,14);

Upfc.gamma = atan2(vq,vp)+theta-a1;
ss = sin(a1-a2+Upfc.gamma);
cc = cos(a1-a2+Upfc.gamma);

switch flag
 
 case 0
  
  % reset transmission line reactance and admittance matrix
  Line.con(Upfc.line,9) = Line.con(Upfc.line,9) + Upfc.xcs;
  fm_y;

  % if vp = 0, gamma is as follows:
  Upfc.gamma = pi/2+theta-a1;

  DAE.x(Upfc.vp) = 0;
  DAE.x(Upfc.vq) = kp.*V1.*sin(a1-a2)./sin(a1-a2+Upfc.gamma);
  DAE.x(Upfc.iq) = Bus.Qg(Upfc.bus1)./V1;
  
  idx = find(DAE.x(Upfc.vp) > vp_max);
  if idx, upfcwarn(idx,' Vp is over its max limit.'), end
  idx = find(DAE.x(Upfc.vp) < vp_min);
  if idx, upfcwarn(idx,' Vp is under its min limit.'), end

  idx = find(DAE.x(Upfc.vq) > vq_max);
  if idx, upfcwarn(idx,' Vq is over its max limit.'), end
  idx = find(DAE.x(Upfc.vq) < vq_min);
  if idx, upfcwarn(idx,' Vq is under its min limit.'), end

  idx = find(DAE.x(Upfc.iq) > iq_max);
  if idx, upfcwarn(idx,' Ish is over its max limit.'), end
  idx = find(DAE.x(Upfc.iq) < iq_min);
  if idx, upfcwarn(idx,' Ish is under its min limit.'), end
  
  DAE.x(Upfc.vp) = max(DAE.x(Upfc.vp),vp_min);
  DAE.x(Upfc.vp) = min(DAE.x(Upfc.vp),vp_max);
  
  DAE.x(Upfc.vq) = max(DAE.x(Upfc.vq),vq_min);
  DAE.x(Upfc.vq) = min(DAE.x(Upfc.vq),vq_max);
  
  DAE.x(Upfc.iq) = max(DAE.x(Upfc.iq),iq_min);
  DAE.x(Upfc.iq) = min(DAE.x(Upfc.iq),iq_max);

  % reference voltages
  Upfc.vp0 = DAE.x(Upfc.vp);
  Upfc.vq0 = DAE.x(Upfc.vq);
  Upfc.Vref = DAE.x(Upfc.iq)./Kr + V1;   

  % eliminate PV components used for initializing UPFC's
  global Syn
  for i = 1:Upfc.n
    idxg = find(Syn.bus == Upfc.bus1(i));
    if ~isempty(idxg)
      upfcwarn(i,[' UPFC cannot be connected at the same bus as ' ...
                 'synchronous machines.'])
      continue
    end
    idx = findbus(PV,Upfc.bus1(i));
    PV = remove(PV,idx);
    if isempty(idx)
      upfcwarn(i,' no PV generator found at the bus.')
    end
  end
    
  fm_disp('Initialization of UPFC completed.')  
   
 case 1 % algebraic equations
  
  c1 = Upfc.y.*sqrt(vp.*vp+vq.*vq);
  
  P1 =  c1.*V2.*ss;
  Q1 =  c1.*V1.*cos(Upfc.gamma)-iq.*V1;
  Q2 = -c1.*V2.*cc;
  
  DAE.gp = DAE.gp + sparse(bus1,1,P1,Bus.n,1);
  DAE.gp = DAE.gp - sparse(bus2,1,P1,Bus.n,1);
  DAE.gq = DAE.gq + sparse(bus1,1,Q1,Bus.n,1);
  DAE.gq = DAE.gq + sparse(bus2,1,Q2,Bus.n,1); 
  
 case 2 % Jacobians of active & reactive powers

  c1 = Upfc.y.*sqrt(vp.*vp+vq.*vq);
  U = max(V1.^2+V2.^2-2.*V1.*V2.*cos(a1-a2),1e-6);
  
  L0 = V1.*V2.*cos(a1-a2);
  L1 = (U-V2.^2+L0)./U;
  L2 = (U-V1.^2+L0)./U;
  L5 = sin(a1-a2)./U;
  L3 = -V2.*L5;
  L4 =  V1.*L5;  
  
  P1a1 =  c1.*L1.*V2.*cc;
  P1a2 =  c1.*(L2-1).*V2.*cc;
  P2a1 = -c1.*L1.*V2.*cc;  
  P2a2 = -c1.*(L2-1).*V2.*cc;  
  P1v1 =  c1.*L3.*V2.*cc;
  P1v2 =  c1.*(ss+V2.*L4.*cc);
  P2v1 = -c1.*L3.*V2.*cc;
  P2v2 = -c1.*(ss+V2.*L4.*cc);  
  Q1a1 = -c1.*(L1-1).*V1.*sin(Upfc.gamma);
  Q1a2 = -c1.*L2.*V1.*sin(Upfc.gamma);
  Q2a1 =  c1.*L1.*V2.*ss;
  Q2a2 =  c1.*(L2-1).*V2.*ss;
  Q1v1 =  c1.*(cos(Upfc.gamma)-V1.*L3.*sin(Upfc.gamma))-iq; 
  Q1v2 = -c1.*L4.*V1.*sin(Upfc.gamma);
  Q2v1 =  c1.*L3.*V2.*ss;
  Q2v2 = -c1.*(cc-V2.*L4.*ss);

  DAE.J11 = DAE.J11 + sparse(bus1,bus1,P1a1,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(bus1,bus2,P1a2,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(bus2,bus1,P2a1,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(bus2,bus2,P2a2,Bus.n,Bus.n);  
  DAE.J12 = DAE.J12 + sparse(bus1,bus1,P1v1,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(bus1,bus2,P1v2,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(bus2,bus1,P2v1,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(bus2,bus2,P2v2,Bus.n,Bus.n);  
  DAE.J21 = DAE.J21 + sparse(bus1,bus1,Q1a1,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(bus1,bus2,Q1a2,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(bus2,bus1,Q2a1,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(bus2,bus2,Q2a2,Bus.n,Bus.n);  
  DAE.J22 = DAE.J22 + sparse(bus1,bus1,Q1v1,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(bus1,bus2,Q1v2,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(bus2,bus1,Q2v1,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(bus2,bus2,Q2v2,Bus.n,Bus.n);
  
 case 3 % differential equations

  vp0 = zeros(Upfc.n,1);
  vq0 = zeros(Upfc.n,1);
  
  ty1 = find(Upfc.con(:,2) == 1);
  ty2 = find(Upfc.con(:,2) == 2);
  
  tyvp = Upfc.con(:,15);
  tyvq = Upfc.con(:,16);
  tyiq = Upfc.con(:,17);
     
  if ty1  
    vp0(ty1) = Upfc.vp0(ty1);
    vq0(ty1) = Upfc.vq0(ty1);
  end
  if ty2
    vq0(ty2) = kp(ty2).*V1(ty2).*sin(a1(ty2)-a2(ty2))./ss(ty2);
    vp0(ty2) = Upfc.vp0(ty2);
  end
  
  Vref = Upfc.Vref;
 
  if Upfc.pod
    Vs = DAE.x(Pod.Vs(Pod.upfc));
    vp0(Upfc.pod) = vp0(Upfc.pod) + tyvp(Upfc.pod).*Vs;
    vq0(Upfc.pod) = vq0(Upfc.pod) + tyvq(Upfc.pod).*Vs;
    Vref(Upfc.pod) = Vref(Upfc.pod) + tyiq(Upfc.pod).*Vs;
  end
  
  u = ~(vp >= vp_max & DAE.f(Upfc.vp) > 0) & ...
      ~(vp <= vp_min & DAE.f(Upfc.vp) < 0);
  DAE.f(Upfc.vp) = u.*(vp0-vp)./Tr; 
    
  u = ~(vq >= vq_max & DAE.f(Upfc.vq) > 0) & ...
      ~(vq <= vq_min & DAE.f(Upfc.vq) < 0);
  DAE.f(Upfc.vq) = u.*(vq0-vq)./Tr;
  
  u = ~(iq >= iq_max & DAE.f(Upfc.iq) > 0) & ...
      ~(iq <= iq_min & DAE.f(Upfc.iq) < 0);
  DAE.f(Upfc.iq) = u.*(Kr.*(Vref-V1)-iq)./Tr;
  
  DAE.x(Upfc.vp) = min(vp_max,DAE.x(Upfc.vp));
  DAE.x(Upfc.vp) = max(vp_min,DAE.x(Upfc.vp));
  DAE.x(Upfc.vq) = min(vq_max,DAE.x(Upfc.vq));
  DAE.x(Upfc.vq) = max(vq_min,DAE.x(Upfc.vq));      
  DAE.x(Upfc.iq) = min(iq_max,DAE.x(Upfc.iq));
  DAE.x(Upfc.iq) = max(iq_min,DAE.x(Upfc.iq));

 case 4 % Jacobians of state variables

  c1 = Upfc.y.*V1;
  vs = sqrt(vp.*vp+vq.*vq);
  c2 = Upfc.y.*vs;
  U = max(V1.^2+V2.^2-2.*V1.*V2.*cos(a1-a2),1e-6);
  
  P1r =  c1.*V2.*ss;
  P2r = -c1.*V2.*ss;
  Q1r =  c1.*V1.*cos(Upfc.gamma);
  Q2r = -c1.*V2.*cc;
  
  P1a =  c2.*V2.*cc;
  P2a = -c2.*V2.*cc;
  Q1a = -c2.*V1.*sin(Upfc.gamma);  
  Q2a =  c2.*V2.*ss;  
  
  P1vp = P1r.*vp./(V1.*vs) - P1a.*vq./vs.^2;
  P2vp = P2r.*vp./(V1.*vs) - P2a.*vq./vs.^2;
  Q1vp = Q1r.*vp./(V1.*vs) - Q1a.*vq./vs.^2;
  Q2vp = Q2r.*vp./(V1.*vs) - Q2a.*vq./vs.^2;  
  
  P1vq = P1r.*vq./(V1.*vs) + P1a.*vp./vs.^2;
  P2vq = P2r.*vq./(V1.*vs) + P2a.*vp./vs.^2;
  Q1vq = Q1r.*vq./(V1.*vs) + Q1a.*vp./vs.^2;
  Q2vq = Q2r.*vq./(V1.*vs) + Q2a.*vp./vs.^2;

  ty1 = Upfc.con(:,2) == 1;
  ty2 = Upfc.con(:,2) == 2;

  tyvp = Upfc.con(:,15);
  tyvq = Upfc.con(:,16);
  tyiq = Upfc.con(:,17);
  
  DAE.Fx = DAE.Fx + sparse(Upfc.vp,Upfc.vp,-1./Tr,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Upfc.vq,Upfc.vq,-1./Tr,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Upfc.iq,Upfc.iq,-1./Tr,DAE.n,DAE.n);
  
  ap = (vp <= vp_max & vp >= vp_min);
  aq = (vq <= vq_max & vq >= vq_min);
  ac = (iq <= iq_max & iq >= iq_min);

  if Upfc.pod
    z = Pod.u(Pod.upfc)./Tr(Upfc.pod);
    DAE.Fx(Upfc.vp(Upfc.pod),Pod.Vs(Pod.upfc)) = z.*tyvp(Upfc.pod).*ap(Upfc.pod);
    DAE.Fx(Upfc.vq(Upfc.pod),Pod.Vs(Pod.upfc)) = z.*tyvq(Upfc.pod).*aq(Upfc.pod);
    DAE.Fx(Upfc.iq(Upfc.pod),Pod.Vs(Pod.upfc)) = z.*Kr(Upfc.pod).*tyiq(Upfc.pod).*ac(Upfc.pod);
  end

  b = find(ap);
  if b
    DAE.Gx = DAE.Gx + sparse(bus1(b),Upfc.vp(b),P1vp(b),2*Bus.n,DAE.n);  
    DAE.Gx = DAE.Gx + sparse(bus2(b),Upfc.vp(b),P2vp(b),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus1(b)+Bus.n,Upfc.vp(b),Q1vp(b),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus2(b)+Bus.n,Upfc.vp(b),Q2vp(b),2*Bus.n,DAE.n);
  end
  
  b = find(aq);
  if b
    a = find(aq.*ty2);
    if a
      Factor(a) = Upfc.Cp(a)./(1-Upfc.Cp(a));
      F1(a) = Factor(a).*V1(a).*V2(a).*sin(a1(a)-a2(a))./U(a).^.5;
      F2(a)=-F1(a);
      F3(a) = Factor(a).*(V1(a)-V2(a).*cos(a1(a)-a2(a)))./U(a).^.5;
      F4(a) = Factor(a).*(V2(a)-V1(a).*cos(a1(a)-a2(a)))./U(a).^.5;     
      DAE.Fy = DAE.Fy + sparse(Upfc.vq(a),bus1(a), F1(a)./Tr(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Upfc.vq(a),bus2(a), F2(a)./Tr(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Upfc.vq(a),bus1(a)+Bus.n, F3(a)./Tr(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Upfc.vq(a),bus2(a)+Bus.n, F4(a)./Tr(a),DAE.n,2*Bus.n);
    end
    DAE.Gx = DAE.Gx + sparse(bus1(b),Upfc.vq(b),P1vq(b),2*Bus.n,DAE.n);  
    DAE.Gx = DAE.Gx + sparse(bus2(b),Upfc.vq(b),P2vq(b),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus1(b)+Bus.n,Upfc.vq(b),Q1vq(b),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus2(b)+Bus.n,Upfc.vq(b),Q2vq(b),2*Bus.n,DAE.n); 
  end
  
  b = find(ac);
  if b
    DAE.Gx = DAE.Gx - sparse(Bus.n+bus1(b),Upfc.iq(b),V1(b),2*Bus.n,DAE.n);
    DAE.Fy = DAE.Fy - sparse(Upfc.iq(b),Bus.n+bus1(b),Kr(b)./Tr(b),DAE.n,2*Bus.n);
  end  
  
 case 5 % non-windup limiters

  a = find((vp >= vp_max | vp <= vp_min) & DAE.f(Upfc.vp) == 0);
  if ~isempty(a)
    DAE.tn(Upfc.vp(a)) = 0;
    DAE.Ac(Upfc.vp(a),:) = 0;
    DAE.Ac(:,Upfc.vp(a)) = 0;
    DAE.Ac(Upfc.vp(a),Upfc.vp(a)) = speye(length(a));
  end
  
  a = find((vq >= vq_max | vq <= vq_min) & DAE.f(Upfc.vq) == 0);
  if ~isempty(a)
    DAE.tn(Upfc.vq(a)) = 0;
    DAE.Ac(Upfc.vq(a),:) = 0;
    DAE.Ac(:,Upfc.vq(a)) = 0;
    DAE.Ac(Upfc.vq(a),Upfc.vq(a)) = speye(length(a));
  end  

  a = find((iq >= iq_max | iq <= iq_min) & DAE.f(Upfc.iq) == 0);
  if ~isempty(a)
    DAE.tn(Upfc.iq(a)) = 0;
    DAE.Ac(Upfc.iq(a),:) = 0;
    DAE.Ac(:,Upfc.iq(a)) = 0;
    DAE.Ac(Upfc.iq(a),Upfc.iq(a)) = speye(length(a));
  end
end

% function for creating warning messages
     function upfcwarn(idx, msg)
global Upfc Varname
fm_disp(strcat('Warning: UPFC #',int2str(idx), ...
	       ' at bus <',Varname.bus(Upfc.bus1(idx)), ...
	       '>: ',msg))