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:    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 Upfc Bus jay DAE Settings

Vdc = max(DAE.x(Upfc.Vdc),1e-6);
Vmac = DAE.x(Upfc.Vmac);
m = DAE.x(Upfc.m);
Vmdc = DAE.x(Upfc.Vmdc);
alpha = DAE.x(Upfc.alpha);
x1 = DAE.x(Upfc.x1);
Imd = DAE.x(Upfc.Imd);
x2 = DAE.x(Upfc.x2);
Imq = DAE.x(Upfc.Imq);

bus1 = Upfc.bus1;
bus2 = Upfc.bus2;

V1 = DAE.V(bus1);
theta1 = DAE.a(bus1);
V2 = DAE.V(bus2);
theta2 = DAE.a(bus2);

omega = Settings.rad;

Rdc = Upfc.con(:,6);
Cdc = Upfc.con(:,7);
Rsh = Upfc.con(:,8);
Xsh = Upfc.con(:,9);
Kmac = Upfc.con(:,10);
Tmac = Upfc.con(:,11);
Kmdc = Upfc.con(:,12);
Tmdc = Upfc.con(:,13);
Vrefac = Upfc.con(:,14);
Vrefdc = Upfc.con(:,15);

Kpdc = Upfc.con(:,16);
Kidc = Upfc.con(:,17);
Kpac = Upfc.con(:,18);
Kiac = Upfc.con(:,19);

Rt = Upfc.con(:,20);
Xt = Upfc.con(:,21);
Rse = Upfc.con(:,22);
Xse = Upfc.con(:,23);
Kp = Upfc.con(:,24);
Ki = Upfc.con(:,25);
Pref = Upfc.con(:,26);
Qref = Upfc.con(:,27);

a1r = Upfc.dat(:,1);
a2r = Upfc.dat(:,2);
a1i = Upfc.dat(:,3);
a2i = Upfc.dat(:,4);

alpha_max = Upfc.dat(:,5);
alpha_min = Upfc.dat(:,6);
m_max = Upfc.dat(:,7);
m_min = Upfc.dat(:,8);
I_max = Upfc.con(:,28);
I_min = Upfc.con(:,29);

Gsh = Upfc.dat(:,9);
Bsh = Upfc.dat(:,10);
K = Upfc.dat(:,11);

switch flag
 case 0   % Initialization

  k0 = 1/4*6^(0.5);

  Zsh2 = Rsh.*Rsh+Xsh.*Xsh;
  Gsh = Rsh./Zsh2;
  Bsh = -Xsh./Zsh2;
  Upfc.dat(:,9) = Gsh;
  Upfc.dat(:,10) = Bsh;
  Upfc.dat(:,11) = omega.*(Rt+Rse)./(Xt+Xse); % K

  Z = Rt+Rse+jay*(Xt+Xse);
  a1 = (Rse+jay*Xse)/Z;
  a2 = (Rt+jay*Xt)/Z;

  Upfc.dat(:,1) = real(a1);
  Upfc.dat(:,2) = real(a2);
  Upfc.dat(:,3) = imag(a1);
  Upfc.dat(:,4) = imag(a2);

  % m_max
  Upfc.dat(:,7) = 1.5*(Vrefac+0.05*I_max)./Vrefdc/k0;
  m_max = 100; %Upfc.dat(:,11);

  % m_min
  Upfc.dat(:,8) = 0.5*(Vrefac-0.05*I_min)./Vrefdc/k0;
  m_min = -100; %Upfc.dat(:,12);

  % alpha_max
  CC = Vrefac.*Vrefac.*Gsh-Vrefdc.*Vrefdc./Rdc-Rsh.*I_max.*I_max;
  BB = -k0.*m_max.*Vrefdc.*Vrefac.*Gsh;
  AA = -k0.*m_max.*Vrefdc.*Vrefac.*Bsh;
  d = AA.*AA+BB.*BB;
  b = BB.*CC./d;
  c = (CC.*CC-AA.*AA)./d;
  Upfc.dat(:,5) = 5*acos(b+sqrt(b.*b-c));

  % alpha_min
  CC = Vrefac.*Vrefac.*Gsh-Vrefdc.*Vrefdc./Rdc-Rsh.*I_min.*I_min;
  BB = -k0.*m_min.*Vrefdc.*Vrefac.*Gsh;
  AA = -k0.*m_min.*Vrefdc.*Vrefac.*Bsh;
  d = AA.*AA+BB.*BB;
  b = BB.*CC./d;
  c = (CC.*CC-AA.*AA)./d;
  Upfc.dat(:,6) = -5*acos(b+sqrt(b.*b-c));

 case 1 % algebraic equations

  t1 = V1.^2;
  t3 = sqrt(6);
  t5 = t3.*m.*Vdc;
  t6 = V1.*Gsh;
  t7 = cos(alpha);
  t11 = V1.*Bsh;
  t12 = sin(alpha);
  t16 = cos(theta1);
  t17 = V1.*t16;
  t18 = Rse+Rt;
  t19 = cos(theta2);
  t20 = V2.*t19;
  t21 = Xse.^2;
  t22 = t21.*x1;
  t24 = sqrt(2);
  t25 = Kp.*t24;
  t26 = t25.*Pref;
  t29 = Kp.^2;
  t30 = t21.*t29;
  t31 = t30.*t24;
  t33 = Pref.*Imd.*V2;
  t36 = theta2-theta1;
  t37 = sin(t36);
  t39 = V1.*t37.*V2;
  t40 = omega.*Xt;
  t41 = Kp.*Qref;
  t45 = Xt.*t29;
  t46 = Qref.^2;
  t50 = Pref.^2;
  t54 = cos(t36);
  t55 = V1.*t54;
  t56 = V2.^2;
  t57 = t56.*V2;
  t58 = omega.^2;
  t62 = Imq.^2;
  t63 = t62.*t56;
  t65 = Xt.^2;
  t66 = t65.*t29;
  t68 = x2.^2;
  t70 = t56.*Xse;
  t73 = Imd.^2;
  t74 = t56.*t73;
  t77 = x1.^2;
  t81 = t65.*t58;
  t83 = t21.*t58;
  t87 = t24.*V1;
  t89 = t87.*t37.*t56;
  t90 = Kp.*Imq;
  t94 = t87.*t37;
  t95 = t56.*omega;
  t100 = Xse.*omega;
  t107 = 2.*t22.*V2.*t26-2.*t31.*t33-4.*t39.*t40.*t41+4.*t45.*t46.*Xse+4.*t45.*t50.*Xse-4.*t55.*t57.*t58+t30.*t63+t66.*t63+2.*Xt.*t68.*t70+t30.*t74+t66.*t74+2.*Xt.*t77.*t70+t81.*t74+t83.*t74+t83.*t63+t81.*t63+2.*t89.*t40.*t90-2.*t94.*t95.*Xse.*x2-4.*t39.*t100.*t41+2.*t89.*t100.*t90;
  t108 = t65.*x2;
  t110 = t25.*Qref;
  t113 = Xt.*x2;
  t115 = Xse.*Kp;
  t116 = t24.*Qref;
  t120 = t66.*t24;
  t122 = Qref.*Imq.*V2;
  t125 = t45.*t24;
  t127 = Imq.*V2;
  t134 = t21.*x2;
  t140 = t65.*x1;
  t141 = t95.*Imq;
  t144 = Xt.*t58;
  t145 = t63.*Xse;
  t150 = t95.*Imd;
  t153 = t74.*Xse;
  t158 = t24.*t57;
  t166 = Xt.*x1;
  t167 = t166.*t56;
  t171 = Xt.*Kp;
  t173 = Imd.*V2;
  t177 = t24.*Pref;
  t182 = t87.*t54;
  t183 = t56.*t58;
  t190 = 2.*t108.*V2.*t110+4.*t113.*V2.*t115.*t116-2.*t120.*t122-4.*t125.*Qref.*Xse.*t127+2.*t1.*t56.*t58+2.*t134.*V2.*t110-2.*t31.*t122-2.*t140.*t141+2.*t144.*t145-2.*t22.*t141+2.*t108.*t150+2.*t144.*t153+2.*t134.*t150-2.*t158.*t58.*Xse.*Imq-2.*t158.*t144.*Imq-4.*t167.*t100.*Imq+4.*t171.*t116.*t100.*t173-4.*t171.*t177.*t100.*t127+2.*t182.*t183.*Xt.*Imq-2.*t120.*t33;
  t197 = t65.*Kp.*t24;
  t199 = Pref.*omega.*t127;
  t206 = t113.*t56;
  t225 = t21.*Kp.*t24;
  t235 = Qref.*omega.*t173;
  t240 = t56.*Kp;
  t241 = t240.*Imd;
  t245 = t87.*t54.*t56;
  t246 = Kp.*Imd;
  t250 = t115.*Imd;
  t255 = 2.*t182.*t183.*Xse.*Imq-2.*t197.*t199-2.*t94.*t183.*Xse.*Imd+4.*t206.*t100.*Imd-2.*t94.*t183.*Xt.*Imd+t65.*t68.*t56+2.*t30.*t46+2.*t30.*t50+2.*t66.*t50+t21.*t77.*t56-2.*t225.*t199+t65.*t77.*t56+2.*t66.*t46+t21.*t68.*t56+2.*t225.*t235+2.*t197.*t235-2.*t140.*t241+2.*t245.*t40.*t246-4.*t167.*t250+2.*t45.*t153;
  t263 = t240.*Imq;
  t275 = t55.*V2;
  t276 = Kp.*Pref;
  t291 = t56.^2;
  t297 = t158.*omega;
  t322 = 2.*t158.*t40.*x1+2.*t291.*t58+4.*t95.*t171.*Pref-2.*t297.*t171.*Imd+2.*t158.*t100.*x1+4.*t166.*V2.*t115.*t177-2.*t297.*t250+2.*t140.*V2.*t26-2.*t108.*t263+4.*t95.*t115.*Pref-4.*t125.*Pref.*Xse.*t173;
  t330 = sqrt((t107+t190+t255-2.*t22.*t241-4.*t206.*t115.*Imq+2.*t45.*t145-2.*t134.*t263-2.*t94.*t95.*t113+2.*t245.*t100.*t246-2.*t182.*t95.*t166-4.*t275.*t40.*t276-2.*t182.*t95.*Xse.*x1-4.*t275.*t100.*t276+t322)./t58./t56);
  t331 = t330.*t24;
  t332 = Xt+Xse;
  t333 = 1./V2;
  t340 = 1./omega;
  t354 = atan((t94-t332.*(x2+Kp.*(t116.*t333-Imq)+omega.*Imd).*t340)./(t182-t24.*V2-t332.*(x1+Kp.*(t177.*t333-Imd)-omega.*Imq).*t340));
  t355 = theta2-t354;
  t356 = cos(t355);
  t359 = t17-t20-t331.*t356./2;
  t361 = t18.^2;
  t362 = t332.^2;
  t364 = 1./(t361+t362);
  t366 = sin(theta1);
  t367 = V1.*t366;
  t368 = sin(theta2);
  t369 = V2.*t368;
  t370 = sin(t355);
  t373 = t367-t369-t331.*t370./2;
  t376 = t18.*t359.*t364+t332.*t373.*t364;
  t382 = -t332.*t359.*t364+t18.*t373.*t364;
  DAE.gp = DAE.gp + sparse(bus1,1, t1.*Gsh-t5.*t6.*t7./4-t5.*t11.*t12./4+t17.*t376+t367.*t382,Bus.n,1);
  DAE.gq = DAE.gq + sparse(bus1,1, -t1.*Bsh+t5.*t11.*t7./4-t5.*t6.*t12./4+t367.*t376-t17.*t382,Bus.n,1);
  DAE.gp = DAE.gp + sparse(bus2,1, -t20.*t376-t369.*t382,Bus.n,1);
  DAE.gq = DAE.gq + sparse(bus2,1, -t369.*t376+t20.*t382,Bus.n,1);

 case 2 % Algebraic Jacobians

  t1 = sin(theta1);
  t2 = V1.*t1;
  t3 = Rse+Rt;
  t4 = cos(theta1);
  t5 = V1.*t4;
  t6 = cos(theta2);
  t7 = V2.*t6;
  t8 = Xse.^2;
  t9 = omega.^2;
  t10 = t8.*t9;
  t11 = Imd.^2;
  t12 = V2.^2;
  t13 = t11.*t12;
  t15 = Imq.^2;
  t16 = t15.*t12;
  t18 = Xt.^2;
  t19 = t18.*t9;
  t21 = Kp.^2;
  t22 = Xt.*t21;
  t23 = Qref.^2;
  t27 = Pref.^2;
  t31 = theta2-theta1;
  t32 = cos(t31);
  t33 = V1.*t32;
  t34 = t12.*V2;
  t35 = t34.*t9;
  t38 = t8.*t21;
  t40 = t18.*t21;
  t42 = x2.^2;
  t43 = Xt.*t42;
  t44 = t12.*Xse;
  t49 = x1.^2;
  t50 = Xt.*t49;
  t54 = t8.*x2;
  t55 = t12.*Kp;
  t56 = t55.*Imq;
  t59 = sqrt(2);
  t60 = t59.*V1;
  t61 = sin(t31);
  t62 = t60.*t61;
  t63 = t12.*omega;
  t64 = Xt.*x2;
  t65 = t63.*t64;
  t68 = t60.*t32;
  t69 = Xt.*x1;
  t70 = t63.*t69;
  t73 = t12.*t9;
  t74 = Xse.*Imd;
  t75 = t73.*t74;
  t78 = t33.*V2;
  t79 = omega.*Xt;
  t80 = Kp.*Pref;
  t81 = t79.*t80;
  t85 = t60.*t32.*t12;
  t86 = Kp.*Imd;
  t87 = t79.*t86;
  t90 = t64.*t12;
  t91 = Xse.*omega;
  t92 = t91.*Imd;
  t95 = t10.*t13+t10.*t16+t19.*t16+4.*t22.*t23.*Xse+4.*t22.*t27.*Xse-4.*t33.*t35+t38.*t16+t40.*t16+2.*t43.*t44+t38.*t13+t40.*t13+2.*t50.*t44+t19.*t13-2.*t54.*t56-2.*t62.*t65-2.*t68.*t70-2.*t62.*t75-4.*t78.*t81+2.*t85.*t87+4.*t90.*t92;
  t96 = Xt.*Imd;
  t97 = t73.*t96;
  t100 = Xse.*x1;
  t101 = t63.*t100;
  t104 = t91.*t80;
  t107 = t91.*t86;
  t110 = t59.*t34;
  t111 = t79.*x1;
  t114 = Xt.*Kp;
  t115 = t114.*Pref;
  t118 = t110.*omega;
  t119 = t114.*Imd;
  t122 = t91.*x1;
  t126 = t8.*Kp.*t59;
  t127 = Pref.*omega;
  t128 = Imq.*V2;
  t129 = t127.*t128;
  t132 = Xse.*Kp;
  t133 = t132.*Pref;
  t136 = t132.*Imd;
  t139 = t18.*x1;
  t141 = Kp.*t59;
  t142 = t141.*Pref;
  t145 = t12.^2;
  t148 = t18.*t42;
  t152 = t8.*t42;
  t156 = V1.^2;
  t162 = t8.*t49;
  t164 = -2.*t62.*t97-2.*t68.*t101-4.*t78.*t104+2.*t85.*t107+2.*t110.*t111+4.*t63.*t115-2.*t118.*t119+2.*t110.*t122-2.*t126.*t129+4.*t63.*t133-2.*t118.*t136+2.*t139.*V2.*t142+2.*t145.*t9+t148.*t12+2.*t38.*t23+t152.*t12+2.*t38.*t27+2.*t156.*t12.*t9+2.*t40.*t27+t162.*t12;
  t166 = t18.*t49;
  t170 = t69.*V2;
  t171 = t59.*Pref;
  t175 = t40.*t59;
  t177 = Pref.*Imd.*V2;
  t180 = t22.*t59;
  t181 = Pref.*Xse;
  t182 = Imd.*V2;
  t186 = t8.*x1;
  t190 = t38.*t59;
  t194 = t60.*t61.*t12;
  t195 = Kp.*Imq;
  t196 = t79.*t195;
  t199 = Xse.*x2;
  t200 = t63.*t199;
  t203 = t64.*V2;
  t204 = t59.*Qref;
  t209 = Qref.*Imq.*V2;
  t212 = Qref.*Xse;
  t217 = t141.*Qref;
  t222 = t63.*Imq;
  t225 = t18.*x2;
  t228 = t55.*Imd;
  t231 = t69.*t12;
  t236 = Xt.*t9;
  t237 = t16.*Xse;
  t240 = t166.*t12+2.*t40.*t23+4.*t170.*t132.*t171-2.*t175.*t177-4.*t180.*t181.*t182+2.*t186.*V2.*t142-2.*t190.*t177+2.*t194.*t196-2.*t62.*t200+4.*t203.*t132.*t204-2.*t175.*t209-4.*t180.*t212.*t128+2.*t54.*V2.*t217-2.*t190.*t209-2.*t139.*t222-2.*t225.*t56-2.*t139.*t228-4.*t231.*t136-2.*t186.*t222+2.*t236.*t237;
  t241 = Qref.*omega;
  t242 = t241.*t182;
  t245 = t63.*Imd;
  t248 = t13.*Xse;
  t253 = t9.*Xse;
  t254 = t253.*Imq;
  t257 = t236.*Imq;
  t264 = t91.*Imq;
  t268 = t18.*Kp.*t59;
  t271 = V1.*t61;
  t272 = t271.*V2;
  t273 = Kp.*Qref;
  t274 = t79.*t273;
  t282 = t91.*t195;
  t288 = t91.*t273;
  t291 = t132.*Imq;
  t296 = Xse.*Imq;
  t297 = t73.*t296;
  t300 = Xt.*Imq;
  t301 = t73.*t300;
  t310 = 2.*t22.*t248-2.*t186.*t228+2.*t194.*t282+2.*t225.*V2.*t217-4.*t272.*t288-4.*t90.*t291+2.*t22.*t237+2.*t68.*t297+2.*t68.*t301-4.*t114.*t171.*t91.*t128+2.*t268.*t242;
  t314 = 1./t9;
  t315 = (t95+t164+t240+2.*t126.*t242+2.*t225.*t245+2.*t236.*t248+2.*t54.*t245-2.*t110.*t254-2.*t110.*t257+4.*t114.*t204.*t91.*t182-4.*t231.*t264-2.*t268.*t129-4.*t272.*t274+t310).*t314;