%===============================================================================================
% 程序OptUPFCConDel.m實現(xiàn)了一個含UPFC裝置的單機(jī)無窮大系統(tǒng)的全狀態(tài)反饋線性最優(yōu)控制設(shè)計。
% 其中發(fā)電機(jī)用恒次暫太Eq1模型，忽略全網(wǎng)所有的電阻。
%===============================================================================================

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的單機(jī)—無窮大系統(tǒng)描述：
%===============================================================================================

%（1）發(fā)電機(jī)(G)參數(shù)：
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)參數(shù)：
ST1N=300;UhT1=242;UlT1=10.5;KT1N=10.5/242;usT1=14;
%(3)線路(L)參數(shù)：
l=300;mL=2;xL1=0.42;xL2=0.42;xL0=4*xL1;
%(4)變壓器（T2)參數(shù)：
ST2N=280;UhT2=220;UlT2=121;KT2N=220/121;usT2=14;
%(5)無窮大系統(tǒng)（S)描述：
US=133.1;                       
%(6)UPFC參數(shù)：
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. 標(biāo)幺值換算：
%==============================================================================================

%（1）選擇基準(zhǔn)容量（SB），220kV電壓等級基準(zhǔn)電壓（UB）換算交流系統(tǒng)參數(shù)為標(biāo)幺值：
SB=240;UB=UGN/KT1N;
%發(fā)電機(jī)參數(shù)：
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;
%無窮大系統(tǒng)：
US=US*KT2N/UB;
%UPFC參數(shù)：
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的單機(jī)無窮大系統(tǒng)的描述。
%=============================================================================================

%=============================================================================================
% 3解算含UPFC的潮流，確定穩(wěn)態(tài)運行點：
%=============================================================================================

%(1)等效系統(tǒng)參數(shù)：
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)給定正常運行時所需的已知各變量值：
%發(fā)電機(jī)(額定運行參數(shù)）：
%UGsf=Ugn;
PGs=Pgn*0.7;
PTs=PGs
%無窮大系統(tǒng)：
USs=US;
%UPFC：
UCs=UCN;
UAs=0.97262201770602;aAs=0.63386860225573;Eqs=1.68975237101889;QLs=1*0.24013979496714;
%(3)潮流計算：
%(a)設(shè)定初始值：
xx0=[0.78384993430858
  -0.16035512884840
   0.57608195543976
   0.11952502934240
   1.42433660831351];

%(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('UPFCflow3',xx0,options,[],X1,X2,X3,x2,UAs,aAs,Eqs,PGs,QLs,USs)
f = UPFCflow3(x,X1,X2,X3,x2,UAs,aAs,Eqs,PGs,QLs,USs)
%(c)變量向量化，完成全網(wǎng)潮流計算：
%交流側(cè)：
V1s=x(1)+x(3)*j;U1s=abs(V1s);a1s=angle(V1s);
V2s=x(2)+x(4)*j;U2s=abs(V2s);a2s=angle(V2s);
%%VGs=UGs*cos(x(5))+UGs*sin(x(5))*j;UGs=abs(VGs);aGs=angle(VGs);
VAs=UAs*(cos(aAs)+j*sin(aAs));UAxs=real(VAs);UAys=imag(VAs);
I1s=(VAs-V1s)/X1/j;I2s=(VAs-V2s-USs)/X2/j;I3s=I1s+I2s;
I2xs=real(I2s);I2ys=imag(I2s);
VBs=USs+I2s*j*(xl+xT2);
UBxs=real(VBs);UBys=imag(VBs);
aEqs=x(5);VEqs=Eqs*(cos(aEqs)+j*sin(aEqs));
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);QLsf=imag(SLs);PCs=real(S1s)+real(S2s);

%直流側(cè):死區(qū)角應(yīng)滿足條件0<=r1,r2<=90*pi/180;
r1s=acos(U1s/ku1/UCs);r2s=acos(U2s/ku2/UCs);%th1s=a1s-aAs;th2s=a2s-aAs;

%============================================================================================
%至此完成了含UPFC的單機(jī)無窮大系統(tǒng)穩(wěn)態(tài)運行點的確定。
%============================================================================================


%============================================================================================
%4. 調(diào)用用戶定義函數(shù)UPFCModel.m在穩(wěn)態(tài)運行點建立線性化的基本狀態(tài)空間方程：
%============================================================================================
zz=-X3*(X3-X31)/(X1*X2*(X3-X31+X3*X31*Xpg));
%%(1)基本代數(shù)方程線性化偏差化系數(shù)：

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) 微分方程有側(cè)未知量偏差方程系數(shù)（變量依次為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);
%*dUC;

hPp1=(cPp3*cAx1*cEg2*cAD3+cPp3*cAx3*cAD3+cPp4*cAy1*cEg2*cAD3+cPp1*cAD3+cPp4*cAy3*cAD3);
%*dr1;
hPp2=(cPp3*cAx1*cEg3+cPp2+cPp3*cAx4+cPp4*cAy4+cPp4*cAy1*cEg3);
%*da1