clear;clc;clf;
%--------------------------------------------------------------------
%求解函數方程fun得到發電機功角delta、轉速W及交軸次暫態電勢Eq的初值
%--------------------------------------------------------------------
x0=[0;1.0;1.2;0;1.0];               %給定計算方程組的迭代初值
options=optimset('display','off');
disp('求解潮流方程以及有功、無功平衡方程后所得的各變量的初值:')
x=fsolve(@fun,x0,options)
disp('發電機功角delta的初值:')
delta0=x(1)
disp('發電機轉速W的初值:')
W0=x(2)
disp('發電機交軸次暫態電勢Eq的初值:')
Eq0=x(3)
%---------------------------------------------------------------------
%使用精細積分求解微分方程
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
Vref=1.1;Vo=1.0;XL=0.504;D=0.1;Ke=40;Tj=8;Pm=0.8;  
Tdo1=7.90;                %Tdo1---直軸暫態時間常數
Xd=0.982;Xd1=0.344;Xq=0.982;  %Xd---直軸電抗，Xd1---直軸暫態電抗，Xq---交軸電抗
X2=Xd1+XL+Xd1*XL/(Xd1+XL);
Xdsigma=Xd+XL+Xd*XL/(Xd+XL);
Xd1sigma=X2;
Xqsigma=Xq+XL+Xq*XL/(Xq+XL);
N=20;      
m=2^N;
tao=0.01;           
Delta_T=tao/m;   %積分時間區段
I=eye(3);        %單位陣
V0=[delta0;W0;Eq0];       %發電機功角、轉速及內電勢的積分初值
for n=1:5
    delta=V0(1,1);
    W=V0(2,1);
    Eq=V0(3,1);
    UGd0=Vo*Xq*sin(delta)/Xqsigma;
    UGq0=(1-(Xdsigma-Xd)/Xd1sigma)*Vo*cos(delta)+Eq*(Xdsigma-Xd)/Xd1sigma;
    UG0=sqrt(UGd0^2+UGq0^2);
    K1=Eq*Vo*cos(delta)/Xd1sigma+Vo^2*(Xd1sigma-Xqsigma)*cos(2*delta)/(Xd1sigma*Xqsigma);
    K2=Vo*sin(delta)/Xd1sigma;
    K3=Xd1sigma/Xdsigma;
    K4=(Xdsigma-Xd1sigma)*Vo*sin(delta)/Xd1sigma;
    K5=UGd0*Vo*Xq*cos(delta)/(UG0*Xqsigma)-UGq0*Vo*Xd1*sin(delta)/(UG0*Xd1sigma);
    K6=UGq0*((Xdsigma-Xd)/Xd1sigma)/UG0;
    Ddelta=(W-1)*100*pi;
    DW=(Pm-Eq*Vo*sin(delta)/X2-D*(W-1))/Tj;
    DEq=(-Ke*UG0-Xdsigma*Eq/Xd1sigma-(Xdsigma-Xd1sigma)*Vo*cos(delta)/Xd1sigma);
    A=[Ddelta;DW;DEq];
    H=[0,100*pi,0;-K1/Tj,-D/Tj,-K2/Tj;-(K4+Ke*K5)/Tdo1,0,-(1/K3+Ke*K6)/Tdo1];
    Ta=H*Delta_T+(H*Delta_T)^2*(I+(H*Delta_T)/3+(H*Delta_T)^2/12)/2;
    for k=0:N
        Ta=2*Ta+Ta*Ta;
    end
    T=I+Ta;
    A=T*A;
    V0=V0+A
end
Y0=V0
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
X3=Xd1+XL; 
Xd1sigma=X3;
for n=0:50
    delta=Y0(1,1);
    W=Y0(2,1);
    Eq=Y0(3,1);
    Ddelta=(W-1)*100*pi;
    DW=(Pm-Eq*Vo*sin(delta)/X3-D*(W-1))/Tj;
    DEq=(-Ke*UG0-Xdsigma*Eq/Xd1sigma-(Xdsigma-Xd1sigma)*Vo*cos(delta)/Xd1sigma);
    A=[Ddelta;DW;DEq];
    t=n*tao;
    subplot(2,2,1),plot(t,delta),hold on,grid on,title('發電機功角\delta')
    xlabel('t(s)'),ylabel('\delta'),text(0,delta0,'\fontsize{10}\bullet\fontname{宋體}功角初值')
    subplot(2,2,2),plot(t,W),hold on,grid on,title('發電機轉速W')
    xlabel('t(s)'),ylabel('W'),text(0,W0,'\fontsize{10}\bullet\fontname{宋體}轉速初值')
    subplot(2,2,3),plot(t,Eq),hold on,grid on,title('發電機內電勢Eq')
    xlabel('t(s)'),ylabel('Eq'),text(0,Eq0,'\fontsize{10}\bullet\fontname{宋體}內電勢初值')
    subplot(2,2,4),plot3(delta,W,Eq),hold on,grid off,title('三個狀態量的關系'),axis vis3d
    xlabel('delta'),ylabel('W'),zlabel('Eq'),text(delta0,W0,Eq0,'\fontsize{10}\bullet\fontname{宋體}初始值')
    view([-114,54])
    H=[0,100*pi,0;-K1/Tj,-D/Tj,-K2/Tj;-(K4+Ke*K5)/Tdo1,0,-(1/K3+Ke*K6)/Tdo1];
    Ta=H*Delta_T+(H*Delta_T)^2*(I+(H*Delta_T)/3+(H*Delta_T)^2/12)/2;
    for k=0:N
        Ta=2*Ta+Ta*Ta;
    end
    T=I+Ta;
    A=T*A;
    Y0=Y0+A
end