optupfcconde1.m

matlab在电力系统中的应用

%===============================================================================================
% 程序OptUPFCConDel.m实现了一个含UPFC装置的单机无穷大系统的全状态反馈线性最优控制设计。
% 其中发电机用恒次暂太Eq1模型，忽略全网所有的电阻。
%===============================================================================================

clear;

global ku1 ku2  Cd Tj D Tr1 Tr2 Ta1 Ta2 PTs SB
global X1 X2 X3  X31 Xpg zz xl xT2 K x2
global Eq1s USs wN UCs PLs QLs PCs aEqs r1s a1s r2s a2s  UAs 



%===============================================================================================
%1. 含UPFC的单机—无穷大系统描述：
%===============================================================================================

%（1）发电机(G)参数：
SGN=300;PGN=240;QGN=180;UGN=10.5;wN=2*pi*50;
xd=1.32;xd1=0.3;%xd2=0.44;
Tj=6;D=5.0;
%(2)变压器（T1)参数：
ST1N=300;UhT1=242;UlT1=10.5;KT1N=10.5/242;usT1=14;
%(3)线路(L)参数：
l=300;mL=2;xL1=0.42;xL2=0.42;xL0=4*xL1;
%(4)变压器（T2)参数：
ST2N=280;UhT2=220;UlT2=121;KT2N=220/121;usT2=14;
%(5)无穷大系统（S)描述：
US=133.1;                       
%(6)UPFC参数：
x1=0.121;x2=24.2;n1=20;n2=10;Cd=1e0;UCN=25;
ku1=2/pi*n1;ku2=2/pi*n2;
Tr1=0.03;Tr2=0.03;Ta1=0.04;Ta2=0.04;

%==============================================================================================
%2. 标幺值换算：
%==============================================================================================

%（1）选择基准容量（SB），220kV电压等级基准电压（UB）换算交流系统参数为标幺值：
SB=240;UB=UGN/KT1N;
%发电机参数：
Ugn=UGN/KT1N/UB;
Pgn=PGN/SB;Qgn=QGN/SB;
Tj=Tj*SGN/SB;
xd=xd*UGN^2/SGN/KT1N^2*SB/UB^2;
xd1=xd1*UGN^2/SGN/KT1N^2*SB/UB^2;
%变压器T1:
xT1=usT1/100*UhT1^2/ST1N*SB/UB^2;
%线路：
xl=l/mL*xL1*SB/UB^2;
%变压器T2:
xT2=usT2/100*UhT2^2/ST2N*SB/UB^2;
%无穷大系统：
US=US*KT2N/UB;
%UPFC参数：
x1=x1/(UB^2/SB);  x2=x2/(UB^2/SB);
ku1=ku1/UB;       ku2=ku2/UB;         %U1=ku1*UC*cos(r1);U2=ku2*UC*cos(r2);

%=============================================================================================
%至此完成了含UPFC的单机无穷大系统的描述。
%=============================================================================================


%=============================================================================================
% 3解算含UPFC的潮流，确定稳态运行点：
%=============================================================================================

%(1)等效系统参数：
X1=x1;X2=x2+xl+xT2;    X3=xd+xT1;X31=xd1+xT1;
Xpg=1/X1+1/X2+1/X3;    Xpt=1/X1+1/X2+1/xT1;
%(2)给定正常运行时所需的已知各变量值：
%发电机(额定运行参数）：
UGs=Ugn;
PGs=Pgn;
PTs=PGs;
%无穷大系统：
USs=US;
%UPFC：
UCs=UCN;
Q1s=-.3;U2s=0.2;a2s=(-36.7+180)*pi/180;
f2s=U2s*cos(a2s);g2s=U2s*sin(a2s);
%(3)潮流计算：
%(a)设定初始值：
f10=1;g10=0.1;fA0=1;gA0=0.1;aG0=50*pi/180;
xx0=[f10;g10;fA0;gA0;aG0];
xx0=[ 0.83158657851831
   0.51024630890743
   0.83141637433853
   0.51023210010221
   0.66549236517005];

%(b)计算UPFC变量：
options=foptions;
options(1)=1;options(2)=1e-10;options(3)=1e-12;
options(7)=1;options(5)=1;options(14)=1e5;                                                                                                                                                                                                                                                                                                                                                                                                                                                                 ;
x = fsolve('UPFCflow1',xx0,options,[],X1,X2,Xpt,xT1,USs,UGs,PGs,Q1s,f2s,g2s);
f = UPFCflow1(x,X1,X2,Xpt,xT1,USs,UGs,PGs,Q1s,f2s,g2s);
%(c)变量向量化，完成全网潮流计算：
%交流侧：
V1s=x(1)+x(2)*j;U1s=abs(V1s);a1s=angle(V1s);
V2s=f2s+g2s*j;U2s=abs(V2s);a2s=angle(V2s);
VGs=UGs*cos(x(5))+UGs*sin(x(5))*j;UGs=abs(VGs);aGs=angle(VGs);
VAs=x(3)+x(4)*j;UAs=abs(VAs);aAs=angle(VAs);UAxs=x(3);UAys=x(4);
I1s=(VAs-V1s)/X1/j;I2s=(VAs-V2s-USs)/X2/j;Is=I1s+I2s;
I2xs=real(I2s);I2ys=imag(I2s);aI2s=angle(I2s);
VBs=USs+I2s*j*(xl+xT2);
UBxs=real(VBs);UBys=imag(VBs);
VEqs=VAs+Is*X3*j;Eqs=abs(VEqs);aEqs=angle(VEqs);
Eq1s=(Eqs+(X3-X31)/X31*UAs*cos(aEqs-aAs))*X31/X3;VEq1s=Eq1s*(cos(aEqs)+j*sin(aEqs));
VA1s=VAs-V2s-I2s*x2*j;UA1s=abs(VA1s);aA1s=angle(VA1s);
SGs=VGs*conj(Is);S1s=V1s*conj(I1s);S2s=V2s*conj(I2s);SLs=VBs*conj(I2s);
PLs=real(SLs);QLs=imag(SLs);PCs=real(S1s)+real(S2s);

%直流侧:死区角应满足条件0<=r1,r2<=90*pi/180;
r1s=acos(U1s/ku1/UCs);r2s=acos(U2s/ku2/UCs);%th1s=a1s-aAs;th2s=a2s-aAs;

save zyUPFCdata0_1 wN Tj D Tr1 Tr2 Ta1 Ta2 SB Cd UCs


%============================================================================================
%至此完成了含UPFC的单机无穷大系统稳态运行点的确定。
%============================================================================================
%Eq1ss=(-Eqs+zz*(X1*U2s*cos(a2s-aEqs)+X2*U1s*cos(a1s-aEqs)+X1*USs*cos(aEqs)))/(zz*X1*X2*X3*Xpg/(X3-X31));


%============================================================================================
%4. 调用用户定义函数UPFCModel.m在稳态运行点建立线性化的基本状态空间方程：
%============================================================================================
zz=-X3*(X3-X31)/(X1*X2*(X3-X31+X3*X31*Xpg));
%%(1)基本代数方程线性化偏差化系数：

cAD1=ku1*cos(r1s);                                      cAD3=-ku1*UCs*sin(r1s);
cAD2=ku2*cos(r2s);                                      cAD4=-ku2*UCs*sin(r2s);

cEg1=zz*(X1*U2s*sin(a2s-aEqs)+X2*U1s*sin(a1s-aEqs)-X1*USs*sin(aEqs));
cEg2=zz*X2*cos(a1s-aEqs);                               cEg3=-zz*X2*U1s*sin(a1s-aEqs);
cEg4=zz*X1*cos(a2s-aEqs);                               cEg5=-zz*X1*U2s*sin(a2s-aEqs);

cAx1=cos(aEqs)/X3/Xpg;      cAx3=cos(a1s)/X1/Xpg;       cAx5=cos(a2s)/X2/Xpg;
cAx2=-Eqs*sin(aEqs)/X3/Xpg; cAx4=-U1s*sin(a1s)/X1/Xpg;  cAx6=-U2s*sin(a2s)/X2/Xpg;

cAy1=sin(aEqs)/X3/Xpg;      cAy3=sin(a1s)/X1/Xpg;       cAy5=sin(a2s)/X2/Xpg;
cAy2=Eqs*cos(aEqs)/X3/Xpg;  cAy4=U1s*cos(a1s)/X1/Xpg;   cAy6=U2s*cos(a2s)/X2/Xpg;

cPG1=(UAxs*sin(aEqs)-UAys*cos(aEqs))/X3;                cPG3=Eqs*sin(aEqs)/X3;
cPG2=Eqs*(UAxs*cos(aEqs)+UAys*sin(aEqs))/X3;            cPG4=-Eqs*cos(aEqs)/X3;

cPp1=(UAys*cos(a1s)-UAxs*sin(a1s))/X1;                  cPp3=-U1s*sin(a1s)/X1;
cPp2=-U1s*(sin(a1s)*UAys+cos(a1s)*UAxs)/X1;             cPp4=U1s*cos(a1s)/X1;

cPs1=(UAys*cos(a2s)+(USs-UAxs)*sin(a2s))/X2;            cPs3=-U2s*sin(a2s)/X2;
cPs2=-U2s*(sin(a2s)*UAys+(UAxs-USs)*cos(a2s))/X2;       cPs4=U2s*cos(a2s)/X2;
%%%%%%%%
cA1=UAxs/sqrt(UAxs^2+UAys^2);                           cA2=UAys/sqrt(UAxs^2+UAys^2);

cIx1=1/X2;            cIx2=-sin(a2s)/X2;                cIx3=-U2s*cos(a2s)/X2;

cIy1=-1/X2;            cIy2=cos(a2s)/X2;                 cIy3=-U2s*sin(a2s)/X2;

cBx1=-(xl+xT2);                                         cBy1=xl+xT2;

cPl1=I2xs;      cPl2=I2ys;        cPl3=UBxs;           cPl4=UBys;

cQl1=-I2ys;      cQl2=I2xs;         cQl3=UBys;           cQl4=-UBxs;

save zyUPFCdata1_1 cAD1 cAD2 cAD3 cAD4 cEg1 cEg2 cEg3 cEg4 cEg5 cAx1 cAx2 cAx3 cAx4... 
cAx5 cAx6 cAy1 cAy2 cAy3 cAy4 cAy5 cAy6 cPG1 cPG2 cPG3 cPG4 cPp1 cPp2 cPp3 cPp4... 
cPs1 cPs2 cPs3 cPs4 cA1 cA2 cIx1 cIx2 cIx3 cIy1 cIy2 cIy3 cBx1 cBy1 cPl1 cPl2 cPl3 cPl4 ...
cQl1 cQl2 cQl3 cQl4;

%%(2) 微分方程有侧未知量偏差方程系数（变量依次为dr1,da1,dr2,da2,daEq,dUC）:

hPG1=(cPG3*cAx1*cEg2*cAD3+cPG4*cAy1*cEg2*cAD3+cPG3*cAx3*cAD3+cPG4*cAy3*cAD3+cPG1*cEg2*cAD3);
%*dr1
hPG2=(cPG4*cAy1*cEg3+cPG1*cEg3+cPG3*cAx1*cEg3+cPG4*cAy4+cPG3*cAx4);
%*da1;
hPG3=(cPG1*cEg4*cAD4+cPG4*cAy1*cEg4*cAD4+cPG3*cAx5*cAD4+cPG4*cAy5*cAD4+cPG3*cAx1*cEg4*cAD4);
%*dr2;
hPG4=(cPG1*cEg5+cPG4*cAy1*cEg5+cPG3*cAx6+cPG4*cAy6+cPG3*cAx1*cEg5);
%*da2;
hPG5=(cPG4*cAy1*cEg1+cPG3*cAx1*cEg1+cPG1*cEg1+cPG2+cPG3*cAx2+cPG4*cAy2);
%*daEq;
hPG6=(cPG4*cAy1*cEg2*cAD1+cPG1*cEg2*cAD1+cPG3*cAx1*cEg4*cAD2+cPG3*cAx1*cEg2*cAD1+cPG4*cAy3*cAD1...
   +cPG4*cAy1*cEg4*cAD2+cPG4*cAy5*cAD2+cPG3*cAx5*cAD2+cPG1*cEg4*cAD2+cPG3*cAx3*cAD1);