smanyx.m

机电模式初值的计算方法

%=======机电模式初值的计算方法=================
方法一：已检验
%利用简化二阶系统矩阵计算的特征值
one=eye(gn,gn);
x=2*pi;
Ge=ke/(one+te*x);
Ggov=0;
Gpss=0;
L1=one+x*tdp+k3+Ge*k6;
L2=-k4-Ge*k5;
LD=D+Ge*Gpss*k2/L1+Ggov;
Lk=k1+k2*L2/L1;
c=[0*one,one*100*pi;-inv(tj)*Lk,-inv(tj)*LD];%利用简化二阶系统矩阵计算的特征值
[v1p,de]=eig(c);
value_first=diag(de);
%---------------------
方法二：周二专已检验
%直接利用矩阵的特征值的平方根计算
one=eye(gn,gn);
x=2*pi;
Ge=ke/(one+te*x);
Ggov=0;
Gpss=0;
L1=one+x*tdp+k3+Ge*k6;
L2=-k4-Ge*k5;
LD=D+Ge*Gpss*k2/L1+Ggov;
Lk=k1+k2*L2/L1;
c=-(2*pi*50*Lk+j*2*pi*LD)/tj;
[votor_z,val]=eig(c);
valu=diag(val);
value_first=sqrt(valu);

==========================k1-k6===============================================
%从计算结果读入初值(初值应考虑到是不是机电模式--虚部是不是在1.26--15.7之间）
value_first=[eig1;eig2];
[num_all,z]=size(value_first);
num_no=0;
num_yes=0;
while (num_yes+num_no<num_all)
   num_yes=num_yes+1;
   %初值考虑到是不是机电模式--虚部是不是在1.26--15.7之间
   if abs(imag(value_first(num_yes)))<1.26|abs(imag(value_first(num_yes)))>15.7
      value_first(num_yes)=[];num_no=num_no+1;num_yes=num_yes-1;
   end
end
%是机电模式的特征值初值
%倪以信SMA法
no=zeros(gn,gn);
one=eye(gn,gn);
%收敛限制
ersp=1e-6;
%机电模式的个数
[n,n1]=size(value_first);
mm=zeros(n,1);%计算收敛前计算的次数
for k=1:n
value=value_first(k);
l=1;
while l==1
%各控制器的传递函数
Ge=ke/(one+te*value);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%计算简化二阶的系统阵
L1s=(one+value*tdp)+k3+Ge*(k6+Gpss*k2);
L2s=-k4-Ge*(k5-Gpss*k1);
Ks=k1+k2*inv(L1s)*L2s;
Ds=D+Ggov+k2*inv(L1s)*Ge*Gpss;
Ar=[no,one*100*pi;-inv(tj)*Ks,-inv(tj)*Ds];
[vot,val]=eig(Ar);
value0=diag(val);
[min_erro,h_s]=min(abs(value-value0));
if min_erro<ersp
      l=0;
end
votor=vot(:,h_s);
value=value0(h_s);
mm(k,1)=mm(k,1)+1;
end %while
%取出真实的计算结果
value_last(k,1)=value;
votor_last(:,k)=votor;
end %for



%================等效二阶系统框图k1-k15模型====================================
%从计算结果读入初值(初值应考虑到是不是机电模式--虚部是不是在1.26--15.7之间）
value_first=[eig1;eig2;eig3;];
[num_all,z]=size(value_first);
num_no=0;
num_yes=0;
while (num_yes+num_no<num_all)
   num_yes=num_yes+1;
   %初值考虑到是不是机电模式--虚部是不是在1.26--15.7之间
   if abs(imag(value_first(num_yes)))<1.26|abs(imag(value_first(num_yes)))>15.7
      value_first(num_yes)=[];num_no=num_no+1;num_yes=num_yes-1;
   end
end
%是机电模式的特征值初值
%倪以信SMA法
no=zeros(gn,gn);
one=eye(gn,gn);
%收敛限制
ersp=1e-6;
%机电模式的个数
[n,n1]=size(value_first);
mm=zeros(n,1);%计算收敛前计算的次数
for k=1:n
value=value_first(k);
l=1;
while l==1&mm(k)<50
%各控制器的传递函数
Ge=ke/(one+te*value);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%计算简化二阶的系统阵
Gd=one/(one+value*tqpp);
Gq=one/(one+value*tdpp);
%Gtd=inv(Gd*(one+value*tqp));
Gtq=one/(Gq*(one+value*tdp));
Gdpp=Gd/(one+k11*Gd);
%Gdpp=Gd/(one+(k11+k17*Gtd)*Gd);
Gqpp=Gq/(one+(k9+(k3+k6*Ge)*Gtq)*Gq);
Gqd=k15;
%Gqd=k15+k16*Gtd;
Gdq=k14+(k13+k7*Ge)*Gtq;
Gdo=k10;
%Gdo=k10+k18*Gtd;
Gqo=k8+(k4+k5*Ge)*Gtq;
Gw=Gpss*Ge*Gtq;
L_s=one/(one-Gdpp*Gqd*Gqpp*Gdq);
L_d=D+Ggov+L_s*k2*Gqpp*Gw-L_s*k12*Gdpp*Gqd*Gqpp*Gw;
L_o=k1+k2*L_s*(Gqpp*Gdq*Gdpp*Gdo-Gqpp*Gqo)+k12*L_s*(Gdpp*Gqd*Gqpp*Gqo-Gdpp*Gdo);
Ar=[no,one*100*pi;-inv(tj)*L_o,-inv(tj)*L_d];
[vot,val]=eig(Ar);
value0=diag(val);
[min_erro,h_s]=min(abs(value-value0));
if min_erro<ersp
      l=0;
end
votor=vot(:,h_s);
value=value0(h_s);
mm(k,1)=mm(k,1)+1;
end %while
%取出真实的计算结果
value_last(k,1)=value;
votor_last(:,k)=votor;
end %for

%================等效二阶系统法SMA框图k1-k18模型====================================
%从计算结果读入初值(初值应考虑到是不是机电模式--虚部是不是在1.26--15.7之间）
value_first=[eig1;eig2;eig3;];
[num_all,z]=size(value_first);
num_no=0;
num_yes=0;
while (num_yes+num_no<num_all)
   num_yes=num_yes+1;
   %初值考虑到是不是机电模式--虚部是不是在1.26--15.7之间
   if abs(imag(value_first(num_yes)))<1.26|abs(imag(value_first(num_yes)))>15.7
      value_first(num_yes)=[];num_no=num_no+1;num_yes=num_yes-1;
   end
end
%是机电模式的特征值初值
%倪以信SMA法
no=zeros(gn,gn);
one=eye(gn,gn);
%收敛限制
ersp=1e-6;
%机电模式的个数
[n,n1]=size(value_first);
mm=zeros(n,1);%计算收敛前计算的次数
for k=1:n
value=value_first(k);
l=1;
while l==1&mm(k)<50
%各控制器的传递函数
Ge=ke/(one+te*value);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%计算简化二阶的系统阵
Gd=one/(one+value*tqpp);
Gq=one/(one+value*tdpp);
%
Gtd=inv(Gd*(one+value*tqp));
Gtq=one/(Gq*(one+value*tdp));
%Gdpp=Gd/(one+k11*Gd);
Gdpp=Gd/(one+(k11+k17*Gtd)*Gd);
Gqpp=Gq/(one+(k9+(k3+k6*Ge)*Gtq)*Gq);
%Gqd=k15;
Gqd=k15+k16*Gtd;
Gdq=k14+(k13+k7*Ge)*Gtq;
%Gdo=k10;
Gdo=k10+k18*Gtd;
Gqo=k8+(k4+k5*Ge)*Gtq;
Gw=Gpss*Ge*Gtq;
L_s=one/(one-Gdpp*Gqd*Gqpp*Gdq);
L_d=D+Ggov+L_s*k2*Gqpp*Gw-L_s*k12*Gdpp*Gqd*Gqpp*Gw;
L_o=k1+k2*L_s*(Gqpp*Gdq*Gdpp*Gdo-Gqpp*Gqo)+k12*L_s*(Gdpp*Gqd*Gqpp*Gqo-Gdpp*Gdo);
Ar=[no,one*100*pi;-inv(tj)*L_o,-inv(tj)*L_d];
[vot,val]=eig(Ar);
value0=diag(val);
[min_erro,h_s]=min(abs(value-value0));
if min_erro<ersp
      l=0;
end
votor=vot(:,h_s);
value=value0(h_s);
mm(k,1)=mm(k,1)+1;
end %while
%取出真实的计算结果
value_last(k,1)=value;
votor_last(:,k)=votor;
end %for

                        以下内容没有用到
=====================k1-k6二阶迭代算法========================================
one=eye(gn,gn);
x=2*pi;
ges=ke/(one+te*x);
L1=x*tdp+k3+ges*k6;%GE(S);
L2=-k4-ges*k5;%GE(S);
LD=D;%GE*GPSS*Ls;
Lk=k1+k2*inv(L1)*L2;
c=[0*one,one*100*pi;-inv(tj)*Lk,-inv(tj)*LD];%要不要改成k1，d;
[v1p,de]=eig(c);
%---------------------
%从计算结果读入初值(初值应考虑到是不是机电模式--虚部是不是在1.26--15.7之间）
value_first=diag(de);
votor_first=v1p;
%倪以信SMA法
no=zeros(gn,gn);
one2=eye(2*gn,2*gn);
%收敛限制
ersp=1e-6;
[n,n1]=size(value_first);
mm=zeros(n,1);
for k=1:n-2
value=value_first(k);
votor=votor_first(:,k);
l=1;
while l==1%&mm<100
%各控制器的传递函数
Ge=ke/(one+te*value);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%各控制系统的传函导数
Gep=-ke*te/(one+te*value)^2;
Ggovp=0;
Gpssp=0;
kp=0;
aw=0;
%计算简化二阶的系统阵
L1s=value*tdp+k3+Ge*(k6+Gpss*k2);
L2s=-k4-Ge*(k5-Gpss*k1);
Ks=k1+k2*inv(L1s)*L2s;
Ds=D+Ggov+k2*inv(L1s)*Ge*Gpss;
Ar=[no,one*100*pi;-inv(tj)*Ks,-inv(tj)*Ds];
%计算系统阵的导数阵
L1sp=tdp+Gep*(k6-Gpss*kp*k2)-Ge*Gpssp*kp*k2;
L2sp=-Gep*(k5-Gpss*kp*k1)+Ge*Gpssp*kp*k1;
Dsp=Ggovp+k2*inv(L1s)*(-L1sp*inv(L1s)*Ge*Gpss+Gep*Gpss+Ge*Gpssp)*aw;
Ksp=k2*inv(L1s)*L2sp-k2*L1sp*(inv(L1s))^2*L2s;
Arp=[no,no;-inv(tj)*Ksp,-inv(tj)*Dsp];
%计算A*
Ar_k=inv(one2-Arp)*(Ar-Arp*value);
%计算第K+1次特征值
value_k=votor'*Ar_k*votor/(votor'*votor);
%计算特征向量
%各控制器的传递函数
Ge=ke/(one+te*value_k);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%计算k+1简化二阶的系统阵
L1s=value*tdp+k3+Ge*(k6+Gpss*k2);
L2s=-k4-Ge*(k5-Gpss*k1);
Ks=k1+k2*inv(L1s)*L2s;
Ds=D+Ggov+k2*inv(L1s)*Ge*Gpss;
Ar=[no,one*100*pi;-inv(tj)*Ks,-inv(tj)*Ds];
%88888
svy=value_k*one2-Ar;
sy=inv(svy)*votor;
votor_k=sy/max(sy);
%=======
if abs(value_k-value)<ersp&max(abs(votor_k-votor))<ersp
   l=0;
end
votor=votor_k;
value=value_k;
mm(k,1)=mm(k,1)+1;
end %while
%取出真实的计算结果
value_last(k,1)=value;
votor_last(:,k)=votor;
end %for

============================k1-k15==================================================
%从计算结果读入初值(初值应考虑到是不是机电模式--虚部是不是在1.26--15.7之间）
value_first=[eig1;eig2];
[num_all,z]=size(value_first);
num_no=0;
num_yes=0;
while (num_yes+num_no<num_all)
   num_yes=num_yes+1;
   %初值考虑到是不是机电模式--虚部是不是在1.26--15.7之间
   if abs(imag(value_first(num_yes)))<1.26|abs(imag(value_first(num_yes)))>15.7
      value_first(num_yes)=[];num_no=num_no+1;num_yes=num_yes-1;
   end
end
%是机电模式的特征值初值
%倪以信SMA法
no=zeros(gn,gn);
one=eye(gn,gn);
%收敛限制
ersp=1e-6;
%机电模式的个数
[n,n1]=size(value_first);
mm=zeros(n,1);%计算收敛前计算的次数
for k=1:n
value=value_first(k);
l=1;
while l==1&mm(k)<50
%各控制器的传递函数
Ge=ke/(one+te*value);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%计算简化二阶的系统阵
L_s=(one+value*tdpp)/(one+value*tdp);
L_w=Gpss*Ge*L_s;
L_9=one/(one+value*tdpp+k9);
L_36=-k3-k6*Ge;
L_9s=L_9/(one-L_9*L_36*L_s);
L_m=-L_9s*k2;
L_713=-k13-k7*Ge;
L_14=-k14+L_713*L_s;
L_11=one/(one+k11+value*tqpp);
L_15=-L_9s*k15*L_11;
L_e=L_15/(one-L_15*L_14);
L_45=-k4-k5*Ge;
L_8=L_45*L_s-k8;
L_10=-L_11*k10;
L_d=L_w*(L_m-L_e*k12)-D-Ggov;
L_o=L_8*(L_m-L_e*k12)-L_10*k12-k1;
Ar=[no,one*100*pi;inv(tj)*L_o,inv(tj)*L_d];
[vot,val]=eig(Ar);
value0=diag(val);
[min_erro,h_s]=min(abs(value-value0));
if min_erro<ersp
      l=0;
end
votor=vot(:,h_s);
value=value0(h_s);
mm(k,1)=mm(k,1)+1;
end %while
%取出真实的计算结果
value_last(k,1)=value;
votor_last(:,k)=votor;
end %for


============================k1-k18没用==================================================
%从计算结果读入初值(初值应考虑到是不是机电模式--虚部是不是在1.26--15.7之间）
value_first=[eig1;eig2];
[num_all,z]=size(value_first);
num_no=0;
num_yes=0;
while (num_yes+num_no<num_all)
   num_yes=num_yes+1;
   %初值考虑到是不是机电模式--虚部是不是在1.26--15.7之间
   if abs(imag(value_first(num_yes)))<1.26|abs(imag(value_first(num_yes)))>15.7
      value_first(num_yes)=[];num_no=num_no+1;num_yes=num_yes-1;
   end
end
%是机电模式的特征值初值
%倪以信SMA法
no=zeros(gn,gn);
one=eye(gn,gn);
%收敛限制
ersp=1e-6;
%机电模式的个数
[n,n1]=size(value_first);
mm=zeros(n,1);%计算收敛前计算的次数
for k=1:n
value=value_first(k);
l=1;
while l==1&mm(k)<50
%各控制器的传递函数
Ge=ke/(one+te*value);%励磁系统
Ggov=0;%原动机系统
Gpss=0;%PSS系统
%计算简化二阶的系统阵
L_q=(one+value*tqpp)/(one+value*tqp);
k15=k15+k16*L_q;
k11=k11+k17*L_q;
k10=k10+k18*L_q;
L_11=one/(one+k11+value*tqpp);
L_9=one/(one+value*tdpp+k9);
L_s=(one+value*tdpp)/(one+value*tdp);
L_713=-k13-k7*Ge;
L_36=-k3-k6*Ge;
L_45=-k4-k5*Ge;
L_w=Gpss*Ge*L_s;
L_9s=L_9/(one-L_9*L_36*L_s);
L_8=L_45*L_s-k8;
L_14=-k14+L_713*L_s;
L_m=-L_9s*k2;
L_15=-L_9s*k15*L_11;
L_e=L_15/(one-L_15*L_14);
L_10=-L_11*k10;
L_d=L_w*(L_m-L_e*k12)-D-Ggov;
L_o=L_8*(L_m-L_e*k12)-L_10*k12-k1;
Ar=[no,one*100*pi; inv(tj)*L_o,inv(tj)*L_d];
[vot,val]=eig(Ar);
value0=diag(val);
[min_erro,h_s]=min(abs(value-value0));
if min_erro<ersp
      l=0;
end
votor=vot(:,h_s);
value=value0(h_s);
mm(k,1)=mm(k,1)+1;
end %while
%取出真实的计算结果
value_last(k,1)=value;
votor_last(:,k)=votor;
end %for





G_t=(one+value*tdpp)/(one+value*tdp);
G_ed=k14+k13*G_t+k7*Ge*G_t;
G_eq=k3*G_t+k6*Ge*G_t+k9;
G_8=k8+k5*Ge*G_t+k4*G_t;
G_w=Gpss*Ge*G_t;
H_eq=one/(one+G_eq+value*tdpp);
H_ed=one/(one+value*tqpp+k11-k15*H_eq);
J_ed=-H_eq*k2+k12;
M_ed=J_ed*H_ed;
L_d=G_w*H_eq*(k2-k15*M_ed)+D+Ggov;
L_o=k1-k10*M_ed-G_8*H_eq*(k2-k15*M_ed);
Ar=[no,one*100*pi;-inv(tj)*L_o,-inv(tj)*L_d];