% 編寫建立自激過程微分方程的M—函數(shù)
% 將該M—函數(shù)定義為sh_ge_se_ex_ode(shunt_generator_self_excited_ode)
function dydt=tongbufadianji(t,y)
% 首先指定全局變量
%global   a1 a2 a3 a4 Rf  
% 下面輸入電機基本數(shù)據(jù)：
r=2.9069E-03,Rfd=5.9013E-04,Rkd=1.1900E-02,Rkq=2.0081E-02;Ufd=24;w=314,Ll=3.0892E-04,Lmd=3.2164E-03,Lmq=9.7153E-04,Llfd=3.0712E-04,Llkd=4.9076E-04,Llkq=1.0365E-03,
Laa0=1/3*(Lmd+Lmq)+Ll;Laa2=1/3*(Lmd-Lmq);Mab0=1/2*Laa0;;Mab2=Laa2;Mafd0=2/3*Lmd,Makd0=2/3*Lmd,Makq0=Lmq,Lfd=Llfd+Lmd,Lkd=Llkd+Lmd,Lkq=Llkq+Lmq,Mfkd=Lmd;
L=[ -(Laa0+Laa2*cos(2*w*t)),            -(-Mab0+Mab2*cos(2*(w*t+2*pi/3))),  -(-Mab0+Mab2*cos(2*(w*t-2*pi/3))),  Mafd0*cos(w*t),	        Makd0*cos(w*t),          -Makq0*sin(w*t); 
    -(-Mab0+Mab2*cos(2*(w*t+2*pi/3))), -(Laa0+Laa2*cos(2*(w*t-2*pi/3))),   -(-Mab0+Mab2*cos(2*(w*t))),         Mafd0*cos(w*t-2*pi/3),	Makd0*cos(w*t-2*pi/3),  -Makq0*sin(w*t-2*pi/3);
    -(-Mab0+Mab2*cos(2*(w*t-2*pi/3))), -(-Mab0+Mab2*cos(2*(w*t))),         -(Laa0+Laa2*cos(2*(w*t+2*pi/3))),   Mafd0*cos(w*t+2*pi/3),	Makd0*cos(w*t+2*pi/3),  -Makq0*sin(w*t+2*pi/3);
    -Mafd0*cos(w*t),                   -Mafd0*cos(w*t-2*pi/3),             -Mafd0*cos(w*t+2*pi/3),             Lfd,                     Mfkd,                     0;
    -Makd0*cos(w*t),                   -Makd0*cos(w*t-2*pi/3),             -Makd0*cos(w*t+2*pi/3),             Mfkd,                    Lkd,                      0;
     Makq0*sin(w*t),                    Makq0*sin(w*t-2*pi/3),              Makq0*sin(w*t+2*pi/3),             0,                       0,                        Lkq]                        
    
    
G=[ 2*Laa2*sin(2*w*t),               2*Mab2*sin(2*(w*t+2*pi/3)),     2*Mab2*sin(2*(w*t-2*pi/3)),       -Mafd0*sin(w*t),	        -Makd0*sin(w*t),          -Makq0*cos(w*t); 
    2*Mab2*sin(2*(w*t+2*pi/3)),     2*Laa2*sin(2*(w*t-2*pi/3)),     2*Mab2*sin(2*(w*t)),              -Mafd0*sin(w*t-2*pi/3),	-Makd0*sin(w*t-2*pi/3),  -Makq0*cos(w*t-2*pi/3);
    2*Mab2*sin(2*(w*t-2*pi/3)),     2*Mab2*sin(2*w*t),            2*Laa2*sin(2*(w*t+2*pi/3)),       -Mafd0*sin(w*t+2*pi/3),	-Makd0*sin(w*t+2*pi/3),  -Makq0*cos(w*t+2*pi/3);
    Mafd0*sin(w*t),                  Mafd0*sin(w*t-2*pi/3),           Mafd0*sin(w*t+2*pi/3),             0,                         0,                        0;
    Makd0*sin(w*t),                  Makd0*sin(w*t-2*pi/3),           Makd0*sin(w*t+2*pi/3),             0,                         0,                        0;
    Makq0*cos(w*t),                  Makq0*cos(w*t-2*pi/3),           Makq0*cos(w*t+2*pi/3),             0,                         0,                        0]
    
% 下面進行曲線擬合：

R=[-r,   0,   0,  0,    0,   0;
    0,  -r,   0,  0,    0,   0;
    0,   0,  -r,  0,    0,   0;
    0,   0,   0,  Rfd,  0,   0;
    0,   0,   0,  0,    Rkd, 0;
    0,   0,   0,  0,    0,   Rkq]
Uabc=[100,0,0,Ufd,0,0]';
y(1)=0;
dydt=L\(Uabc-w*G*y-R*y);





% 下面列寫標準形式微分方程并將感應(yīng)電勢用勵磁電流的擬合多項式函數(shù)表示：
% 這里的a(1)、a(2)、a(3)、a(4)是指上一步的多項式擬合系數(shù)a(i)
%dydt=[1/Lf*(a(1)*iff^3+a(2)*iff^2+a(3)*iff^1+a(4))-Rf/Lf*iff];
% 給全局變量賦值
%a1=a(1);a2=a(2);a3=a(3);a4=a(4);