fm_sssc.m

一个较好的MATLAB潮流程序

  t36 = mse.*Vdc;
  t37 = t36.*t3;
  t38 = t2.*t37;
  t40 = t36.*t7;
  t41 = t6.*t40;
  t43 = a1r.*t32-a1i.*t34+t38./4-t41./4;
  t44 = a0r.*t1;
  t46 = a0i.*t1;
  t48 = -t44.*t4+t46.*t8;
  t55 = a0i.*t19+a0r.*t28;
  t59 = t2.*t40;
  t61 = t6.*t37;
  t63 = a1i.*t32+a1r.*t34+t59./4+t61./4;
  t66 = -t46.*t4-t44.*t8;
  t69 = 1./Cdc;
  t71 = 1./Vdc;
  t77 = Vdc.^2;
  t78 = 1./t77;
  t84 = t30.*t48./4+t55.*t66./4;
  t86 = t69.*t71;
  t88 = t30.^2;
  t89 = t55.^2;
  t101 = -t44.*t16+t46.*t25;
  t109 = -t46.*t16-t44.*t25;
  t116 = t30.*t101./4+t55.*t109./4;
  t124 = t44.*t40+t46.*t37;
  t130 = t46.*t40-t44.*t37;
  t137 = t30.*t124./4+t55.*t130./4;
  t141 = a1r.*V1;
  t143 = a1i.*V1;
  t147 = a0r.*V1;
  t149 = a0i.*V1;
  t151 = -t147.*t21-t149.*t11;
  t159 = -t149.*t21+t147.*t11;
  t166 = t30.*t151+t55.*t159;
  t176 = a0r.*t11-a0i.*t21;
  t184 = a0i.*t11+a0r.*t21;
  t191 = t30.*t176+t55.*t184;
  t195 = a1r.*V2;
  t197 = a1i.*V2;
  t201 = a0r.*V2;
  t203 = a0i.*V2;
  t205 = t201.*t23+t203.*t13;
  t213 = t203.*t23-t201.*t13;
  t220 = t30.*t205+t55.*t213;
  t230 = -a0r.*t13+a0i.*t23;
  t238 = -a0i.*t13-a0r.*t23;
  t245 = t30.*t230+t55.*t238;
  t251 = t12.*t48./4+t22.*t66./4;
  t252 = Kmac.*t251;
  t253 = 1./Tmac;
  t257 = t12.*t101./4+t22.*t109./4;
  t259 = Kmac.*t257.*t253;
  t262 = t12.*t124./4+t22.*t130./4;
  t263 = Kmac.*t262;
  t270 = Kmac.*(-t22.*t30+t12.*t151+t12.*t55+t22.*t159);
  t277 = Kmac.*(t11.*t30+t12.*t176+t21.*t55+t22.*t184);
  t282 = Kmac.*(t12.*t205+t22.*t213);
  t287 = Kmac.*(t12.*t230+t22.*t238);
  t289 = K.*T1;
  t290 = t289.*t253;
  t291 = 1./T2;
  t295 = (-K+t290).*t291;
  t309 = 2.*Kmac.*t84;
  t312 = 2.*Kmac.*t116.*t253;
  t313 = 2.*Kmac.*t137;
  t315 = 2.*Kmac.*t166;
  t317 = 2.*Kmac.*t191;
  t319 = 2.*Kmac.*t220;
  t321 = 2.*Kmac.*t245;
  t338 = 1./Tmdc;

  DAE.Fx = DAE.Fx + sparse(Sssc.Vdc,Sssc.Vdc,((t2.*t4-t6.*t8).*t30./4+t43.*t48./4+(t2.*t8+t6.*t4).*t55./4+t63.*t66./4).*t69.*t71-(t43.*t30+t63.*t55).*t69.*t78-1./Rdc.*t69-2.*Rac.*t84.*t86+Rac.*(t88+t89).*t69.*t78,DAE.n,DAE.n);
  %DAE.Fx = DAE.Fx + sparse(Sssc.Vdc,Sssc.m,((t2.*t16-t6.*t25).*t30./4+t43.*t101./4+(t2.*t25+t6.*t16).*t55./4+t63.*t109./4).*t69.*t71-2.*Rac.*t116.*t86,DAE.n,DAE.n);
  %DAE.Fx = DAE.Fx + sparse(Sssc.Vdc,Sssc.beta,((-t59-t61).*t30./4+t43.*t124./4+(t38-t41).*t55./4+t63.*t130./4).*t69.*t71-2.*Rac.*t137.*t86,DAE.n,DAE.n);
  DAE.Fy = DAE.Fy + sparse(Sssc.Vdc,bus1,((-t141.*t21-t143.*t11).*t30+t43.*t151+(-t143.*t21+t141.*t11).*t55+t63.*t159).*t69.*t71-2.*Rac.*t166.*t86,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Sssc.Vdc,bus1+Bus.n,((a1r.*t11-a1i.*t21).*t30+t43.*t176+(a1i.*t11+a1r.*t21).*t55+t63.*t184).*t69.*t71-2.*Rac.*t191.*t86,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Sssc.Vdc,bus2,((t195.*t23+t197.*t13).*t30+t43.*t205+(t197.*t23-t195.*t13).*t55+t63.*t213).*t69.*t71-2.*Rac.*t220.*t86,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Sssc.Vdc,bus2+Bus.n,((-a1r.*t13+a1i.*t23).*t30+t43.*t230+(-a1i.*t13-a1r.*t23).*t55+t63.*t238).*t69.*t71-2.*Rac.*t245.*t86,DAE.n,2*Bus.n);


  c1 = find(c_type == 1);
  if ~isempty(c1)
    DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.Vdc,t252.*t253,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.Pmac,-t253,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.m,t259,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.beta,t263.*t253,DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus1,t270.*t253,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus1+Bus.n,t277.*t253,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus2,t282.*t253,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus2+Bus.n,t287.*t253,DAE.n,2*Bus.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Vdc,-t290.*t252.*t291,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Pmac,t295,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.m,(-Kd-t289.*t259).*t291,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.beta,-t290.*t263.*t291,DAE.n,DAE.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1,-t290.*t270.*t291,DAE.n,2*Bus.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1+Bus.n,-t290.*t277.*t291,DAE.n,2*Bus.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2,-t290.*t282.*t291,DAE.n,2*Bus.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2+Bus.n,-t290.*t287.*t291, ...
    %                         DAE.n,2*Bus.n);

    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.m,(-Kd-t289.*t259).*t291, ...
    %                         DAE.n,DAE.n);
  end


  c2 = find(c_type == 2);
  if ~isempty(c2)
    DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.Vdc,t309.*t253,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.Pmac,-t253,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.m,t312,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.beta,t313.*t253,DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus1,t315.*t253,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus1+Bus.n,t317.*t253,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus2,t319.*t253,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.Pmac,bus2+Bus.n,t321.*t253,DAE.n,2*Bus.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Vdc,-t290.*t309.*t291,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Pmac,t295,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.m,(-Kd-t289.*t312).*t291,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.beta,-t290.*t313.*t291,DAE.n,DAE.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1,-t290.*t315.*t291,DAE.n,2*Bus.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1+Bus.n,-t290.*t317.*t291,DAE.n,2*Bus.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2,-t290.*t319.*t291,DAE.n,2*Bus.n);
    %DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2+Bus.n,-t290.*t321.*t291,DAE.n,2*Bus.n);

    %DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.m,(-Kd-t289.*t312).*t291,DAE.n,DAE.n);

  end

  DAE.Fx = DAE.Fx + sparse(Sssc.Vmdc,Sssc.Vdc,Kmdc.*t338,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sssc.Vmdc,Sssc.Vmdc,-t338,DAE.n,DAE.n);
  %DAE.Fx = DAE.Fx + sparse(Sssc.beta,Sssc.Vdc,-Kp.*Kmdc.*t338,DAE.n,DAE.n);
  %DAE.Fx = DAE.Fx + sparse(Sssc.beta,Sssc.Vmdc,Kp.*t338-Ki,DAE.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus1,Sssc.Vdc,t251,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus1,Sssc.m,t257,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus1,Sssc.beta,t262,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus1+Bus.n,Sssc.Vdc,t22.*t48./4-t12.*t66./4,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus1+Bus.n,Sssc.m,t22.*t101./4-t12.*t109./4,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus1+Bus.n,Sssc.beta,t22.*t124./4-t12.*t130./4,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus2,Sssc.Vdc,-t14.*t48./4-t24.*t66./4,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus2,Sssc.m,-t14.*t101./4-t24.*t109./4,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus2,Sssc.beta,-t14.*t124./4-t24.*t130./4,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus2+Bus.n,Sssc.Vdc,-t24.*t48./4+t14.*t66./4,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus2+Bus.n,Sssc.m,-t24.*t101./4+t14.*t109./4,2*Bus.n,DAE.n);
  %DAE.Gx = DAE.Gx + sparse(bus2+Bus.n,Sssc.beta,-t24.*t124./4+t14.*t130./4,2*Bus.n,DAE.n);

  % m
  k = find(mse <= m_max & mse >= m_min & (DAE.f(Sssc.m) ~= 0 | DAE.t == -1));
  if ~isempty(k)

    DAE.Fx = DAE.Fx + sparse(Sssc.Vdc,Sssc.m,((t2.*t16-t6.*t25).*t30./4+t43.*t101./4+(t2.*t25+t6.*t16).*t55./4+t63.*t109./4).*t69.*t71-2.*Rac.*t116.*t86,DAE.n,DAE.n);

   DAE.Gx = DAE.Gx + sparse(bus1,Sssc.m,t257,2*Bus.n,DAE.n);
   DAE.Gx = DAE.Gx + sparse(bus1+Bus.n,Sssc.m,t22.*t101./4-t12.*t109./4,2*Bus.n,DAE.n);
   DAE.Gx = DAE.Gx + sparse(bus2,Sssc.m,-t14.*t101./4-t24.*t109./4,2*Bus.n,DAE.n);
   DAE.Gx = DAE.Gx + sparse(bus2+Bus.n,Sssc.m,-t24.*t101./4+t14.*t109./4,2*Bus.n,DAE.n);


   c1 = find(c_type(k) == 1);
    if ~isempty(c1)
      k = k(c1);

      DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.m,t259,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Vdc,-t290.*t252.*t291,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Pmac,t295,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.m,(-Kd-t289.*t259).*t291,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.beta,-t290.*t263.*t291,DAE.n,DAE.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1,-t290.*t270.*t291,DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1+Bus.n,-t290.*t277.*t291,DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2,-t290.*t282.*t291,DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2+Bus.n,-t290.*t287.*t291, ...
                               DAE.n,2*Bus.n);

    end

    c2 = find(c_type(k) == 2);
    if ~isempty(c2)
      k = k(c2);

      DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.m,t312,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Vdc,-t290.*t309.*t291,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.Pmac,t295,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.m,(-Kd-t289.*t312).*t291,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Sssc.m,Sssc.beta,-t290.*t313.*t291,DAE.n,DAE.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1,-t290.*t315.*t291,DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus1+Bus.n,-t290.*t317.*t291,DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2,-t290.*t319.*t291,DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.m,bus2+Bus.n,-t290.*t321.*t291,DAE.n,2*Bus.n);

    end

  end

  % beta
  k = find(beta <= b_max & beta >= b_min & (DAE.f(Sssc.beta) ~= 0 | DAE.t == -1));
  if ~isempty(k)

    DAE.Fx = DAE.Fx + sparse(Sssc.Vdc,Sssc.beta,((-t59-t61).*t30./4+t43.*t124./4+(t38-t41).*t55./4+t63.*t130./4).*t69.*t71-2.*Rac.*t137.*t86,DAE.n,DAE.n);

    DAE.Fx = DAE.Fx + sparse(Sssc.beta,Sssc.Vdc,-Kp.*Kmdc.*t338,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Sssc.beta,Sssc.Vmdc,Kp.*t338-Ki,DAE.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus1,Sssc.beta,t262,2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus1+Bus.n,Sssc.beta,t22.*t124./4-t12.*t130./4,2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus2,Sssc.beta,-t14.*t124./4-t24.*t130./4,2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(bus2+Bus.n,Sssc.beta,-t24.*t124./4+t14.*t130./4,2*Bus.n,DAE.n);

    c1 = find(c_type(k) == 1);
    if ~isempty(c1)
      k = k(c1);

      DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.beta,t263.*t253,DAE.n,DAE.n);
      %idx = find(diag(DAE.Fx(Sssc.m,Sssc.beta)));
      %if ~isempty(idx)
      %  for i = 1:length(idx)
      %    DAE.Fx(Sssc.m(idx(i)),Sssc.beta(idx(i)));
      %  end
      %end

    end

    c2 = find(c_type(k) == 2);
    if ~isempty(c2)
      k = k(c2);

      DAE.Fx = DAE.Fx + sparse(Sssc.Pmac,Sssc.beta,t313.*t253,DAE.n,DAE.n);
      %idx = find(diag(DAE.Fx(Sssc.m,Sssc.beta)));
      %if ~isempty(idx)
      %  for i = 1:length(idx)
      %    DAE.Fx(Sssc.m(idx(i)),Sssc.beta(idx(i)));
      %  end
      %end

    end
  end
  k = find(beta >= b_max & beta <= b_min & ...
           (DAE.f(Sssc.beta) == 0 | DAE.t == -1));
  if ~isempty(k)
    DAE.Fx = DAE.Fx + sparse(Sssc.beta(k),Sssc.beta(k),-1,DAE.n,DAE.n);
    idx = find(diag(DAE.Fx(Sssc.m(k),Sssc.beta(k))));
    if ~isempty(idx)
      for i = 1:length(idx)
        DAE.Fx(Sssc.m(idx(k(i))),Sssc.beta(idx(k(i)))) = 0;
      end
    end
  end


 case 5 % Hard limits

  global Settings

  idx = find((mse >= m_max | mse <= m_min) & DAE.f(Sssc.m) == 0);
  if ~isempty(idx)
    k = Sssc.m(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = -eye(length(idx));
    else
      DAE.Ac(k,k) = -speye(length(idx));
    end
  end

  idx = find((beta >= b_max | beta <= b_min) & DAE.f(Sssc.beta) == 0);
  if ~isempty(idx)
    k = Sssc.beta(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = -eye(length(idx));
    else
      DAE.Ac(k,k) = -speye(length(idx));
    end
  end

end