fm_sae1.m

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

function fm_sae1(flag)
% FM_SAE1 define a subtransmission area equivalent
%         with three loads and three LTCs. This model
%         uses one equivalent state variables presenting
%         a "mean" value of the LTCs' tap ratios.
%
% FM_SAE1(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
%
%see also FM_SAE2 and FM_SAE3
%
%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 Bus DAE SAE1

for i = 1:SAE1.n
   if isempty(SAE1.m)
      m = 1;
   else
      m = DAE.x(SAE1.m(i));
   end
   ha = Bus.int(round(SAE1.con(i,1)));
   hb = Bus.int(round(SAE1.con(i,2)));
   xt1 = SAE1.con(i,6);
   xt2 = SAE1.con(i,7);
   xt3 = SAE1.con(i,8);
   xa0 = SAE1.con(i,9);
   xb0 = SAE1.con(i,10);
   if xt1 == 0; xt1 = 0.2; end
   if xt2 == 0; xt2 = 0.2; end
   if xt3 == 0; xt3 = 0.2; end
   if xa0 == 0; xa0 = 0.2; end
   if xb0 == 0; xb0 = 0.2; end
   a1 = SAE1.con(i,11);
   b1 = SAE1.con(i,12);
   a2 = SAE1.con(i,13);
   b2 = SAE1.con(i,14);
   a3 = SAE1.con(i,15);
   b3 = SAE1.con(i,16);
   h1 = SAE1.con(i,17);
   k1 = SAE1.con(i,18);
   vrif1 = SAE1.con(i,19);
   h2 = SAE1.con(i,20);
   k2 = SAE1.con(i,21);
   vrif2 = SAE1.con(i,22);
   h3 = SAE1.con(i,23);
   k3 = SAE1.con(i,24);
   vrif3 = SAE1.con(i,25);
   mmax1 = SAE1.con(i,26);
   mmin1 = SAE1.con(i,27);
   mmax2 = SAE1.con(i,28);
   mmin2 = SAE1.con(i,29);
   mmax3 = SAE1.con(i,30);
   mmin3 = SAE1.con(i,31);

   m0 = (0.9519 + 0.9474 + 0.9492)/3;
   mmax = (mmax1+mmax2+mmax3)/3;
   mmin = (mmin1+mmin2+mmin3)/3;

   %  coefficienti nel caso mi = m + Bi
   % B1 = 0.9519 - m0;
   % B2 = 0.9474 - m0;
   % B3 = 0.9492 - m0;

   %  coefficienti nel caso m1 = m per qualunque mi
   B1 = 0;
   B2 = 0;
   B3 = 0;

   va = DAE.V(ha);
   vb = DAE.V(hb);
   delta = DAE.a(ha);
   theta = DAE.a(hb);
   if isempty(SAE1.m)
      hm = 1;
   else
      hm = SAE1.m(i);
   end
   t1 = m+B1;
   t2 = 1/t1;
   t4 = 1/xt1;
   t5 = t1*t1;
   t6 = 1/t5;
   t8 = t4+a1*t6;
   t9 = 1/t8;
   t10 = t4*t9;
   t12 = m+B2;
   t13 = 1/t12;
   t15 = 1/xt2;
   t16 = t12*t12;
   t17 = 1/t16;
   t19 = t15+a2*t17;
   t20 = 1/t19;
   t21 = t15*t20;
   t23 = m+B3;
   t24 = 1/t23;
   t26 = 1/xt3;
   t27 = t23*t23;
   t28 = 1/t27;
   t30 = t26+a3*t28;
   t31 = 1/t30;
   t32 = t26*t31;
   t34 = b1*t2*t10+b2*t13*t21+b3*t24*t32;
   t35 = 1/xa0;
   t36 = 1/xb0;
   t37 = xt1*xt1;
   t38 = 1/t37;
   t40 = xt2*xt2;
   t41 = 1/t40;
   t43 = xt3*xt3;
   t44 = 1/t43;
   t46 = t4+t26+t15+t35+t36-t38*t9-t41*t20-t44*t31;
   t47 = 1/t46;
   t48 = t34*t47;
   t49 = va*t35;
   t51 = va*vb;
   t52 = -delta+theta;
   t53 = sin(t52);
   t54 = t51*t53;
   t55 = t47*t35;
   t56 = t55*t36;
   t57 = t54*t56;
   t59 = vb*t36;
   t62 = xa0*xa0;
   t63 = 1/t62;
   t65 = t35-t63*t47;
   t66 = va*va;
   t68 = cos(t52);
   t70 = t56*t51*t68;
   t72 = xb0*xb0;
   t73 = 1/t72;
   t75 = t36-t73*t47;
   t76 = vb*vb;
   t81 = vb*t53*t56;
   t84 = va*t53*t56;
   t87 = t55*t59*t68;
   t91 = t55*t36*va*t68;
   t97 = t49+t59;
   t98 = t97*t47;
   t100 = t2*t4*t9;
   t106 = t13*t15*t20;
   t112 = t24*t26*t31;
   t117 = t46*t46;
   t118 = 1/t117;
   t119 = t97*t118;
   t121 = t8*t8;
   t122 = 1/t121;
   t128 = t19*t19;
   t129 = 1/t128;
   t135 = t30*t30;
   t136 = 1/t135;
   t142 = -2.0*t38*t122*a1/t5/t1-2.0*t41*t129*a2/t16/t12-2.0*t44*t136*a3/t27/t23;
   t148 = t5*t5;
   t149 = 1/t148;
   t152 = t4*t122*a1;
   t162 = t16*t16;
   t163 = 1/t162;
   t166 = t15*t129*a2;
   t176 = t27*t27;
   t177 = 1/t176;
   t180 = t26*t136*a3;
   t198 = (-b1*t6*t10+2.0*b1*t149*t152-b2*t17*t21+2.0*b2*t163*t166-b3*t28* ...
      t32+2.0*b3*t177*t180)*t47;
   t200 = t34*t118;
   t203 = t118*t35;
   t206 = t54*t203*t36*t142;
   t214 = t203*t36*t51*t68*t142;

   switch flag

   case 0

      if length(SAE1.con(1,:)) > 28
          xeq = zeros(5,1);
          A =  [1 0 1 0 0; 1 0 0 1 0; 1 0 0 0 1; 0 1 1 0 0; 0 1 0 1 0; ...
                  0 1 0 0 1; 0 0 1 1 0; 0 0 0 1 1; 0 0 1 0 1; 1 1 0 0 0];
          x1 = SAE1.con(i,35) + SAE1.con(i,32);
          x2 = SAE1.con(i,35) + SAE1.con(i,33) + SAE1.con(i,37);
          x3 = SAE1.con(i,35) + SAE1.con(i,34) + SAE1.con(i,37) + SAE1.con(i,38);
          x4 = SAE1.con(i,36) + SAE1.con(i,32) + SAE1.con(i,38) + SAE1.con(i,37);
          x5 = SAE1.con(i,36) + SAE1.con(i,33) + SAE1.con(i,38);
          x6 = SAE1.con(i,36) + SAE1.con(i,34);
          x7 = SAE1.con(i,32) + SAE1.con(i,37) + SAE1.con(i,33);
          x8 = SAE1.con(i,34) + SAE1.con(i,38) + SAE1.con(i,33);
          x9 = SAE1.con(i,32) + SAE1.con(i,37) + SAE1.con(i,38) + SAE1.con(i,34);
          x10= SAE1.con(i,35) + SAE1.con(i,37) + SAE1.con(i,38) + SAE1.con(i,36);
          xreal = [x1; x2; x3; x4; x5; x6; x7; x8; x9; x10];
          xeq = A\xreal;
          SAE1.con(i,9) = xeq(1);
          SAE1.con(i,10) = xeq(2);
          SAE1.con(i,6) = xeq(3);
          SAE1.con(i,7) = xeq(4);
          SAE1.con(i,8) = xeq(5);
      end

   case 1
      Pa = t48*t49-t57;
      Pb = t48*t59+t57;
      Qa = t65*t66-t70;
      Qb = t75*t76-t70;

      DAE.gp(ha) = Pa + DAE.gp(ha);
      DAE.gq(ha) = Qa + DAE.gq(ha);
      DAE.gp(hb) = Pb + DAE.gp(hb);
      DAE.gq(hb) = Qb + DAE.gq(hb);

      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

   case 2

      DAE.J12(ha,ha) = DAE.J12(ha,ha) + t48*t35-t81;
      DAE.J12(ha,hb) = DAE.J12(ha,hb) - t84;
      DAE.J11(ha,ha) = DAE.J11(ha,ha) + t70;
      DAE.J11(ha,hb) = DAE.J11(ha,hb) - t70;
      DAE.J22(ha,ha) = DAE.J22(ha,ha) + 2.0*t65*va-t87;
      DAE.J22(ha,hb) = DAE.J22(ha,hb) - t91;
      DAE.J21(ha,ha) = DAE.J21(ha,ha) - t57;
      DAE.J21(ha,hb) = DAE.J21(ha,hb) + t57;
      DAE.J12(hb,ha) = DAE.J12(hb,ha) + t81;
      DAE.J12(hb,hb) = DAE.J12(hb,hb) + t48*t36+t84;
      DAE.J11(hb,ha) = DAE.J11(hb,ha) - t70;
      DAE.J11(hb,hb) = DAE.J11(hb,hb) + t70;
      DAE.J22(hb,ha) = DAE.J22(hb,ha) - t87;
      DAE.J22(hb,hb) = DAE.J22(hb,hb) + 2.0*t75*vb-t91;
      DAE.J21(hb,ha) = DAE.J21(hb,ha) -t57;
      DAE.J21(hb,hb) = DAE.J21(hb,hb) +  t57;

      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

   case 3

      %   DAE.f :
      DAE.f(hm) = -h1*t1/3+k1*(t98*t100-vrif1)/3-h2*t12/3+k2*(t98*t106-vrif2)/3-h3* ...
         t23/3+k3*(t98*t112-vrif3)/3;

      if (m >= mmax)
         m = mmax;
         DAE.x(hm) = m;
         if(DAE.f(hm) > 0); DAE.f(hm) = 0;end
      end
      if (m <= mmin)
         m = mmin;
         DAE.x(hm) = m;
         if(DAE.f(hm) < 0); DAE.f(hm) = 0;end
      end

      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

   case 4

      DAE.Fx(hm,hm) = -h1/3+k1*(-t119*t2*t10*t142-t98*t6*t4*t9+2.0*t98*t149*t152)/3-h2 ...
         /3+k2*(-t119*t13*t21*t142-t98*t17*t15*t20+2.0*t98*t163*t166)/3-h3/3+k3*(-t119* ...
         t24*t32*t142-t98*t28*t26*t31+2.0*t98*t177*t180)/3;
      DAE.Gx(ha,hm) = t198*t49-t200*t49*t142+t206;
      DAE.Gx(ha+Bus.n,hm) = t63*t118*t142*t66+t214;
      DAE.Gx(hb,hm) = t198*t59-t200*t59*t142-t206;
      DAE.Gx(hb+Bus.n,hm) = t73*t118*t142*t76+t214;
      DAE.Fy(hm,Bus.n+ha) = k1*t35*t47*t100/3+k2*t35*t47*t106/3+k3*t35*t47*t112/3;
      DAE.Fy(hm,Bus.n+hb) = k1*t36*t47*t100/3+k2*t36*t47*t106/3+k3*t36*t47*t112/3;

      if ((m >= mmax | m <= mmin) & DAE.f(hm) == 0)
         DAE.Fx(hm,hm) = -1;
         DAE.Fy(hm,:) = zeros(1,2*Bus.n);
      end

      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

   case 5

      if ((m >= mmax | m <= mmin) & DAE.f(hm) == 0)
         DAE.tn(hm) = 0;
         DAE.Ac(hm,:) = zeros(1,DAE.n+2*Bus.n);
         DAE.Ac(hm,hm) = DAE.Fx(hm,hm);
      end

   end
end