a.m

计算发电机励磁系统的机断电压功率,能够判断系统的稳定性和能观能控性能

clc
  clear
%潮流计算
p=0.8;Q=0.47;u=1;
xd=1.226;xpd=0.190;
xq=0.830;xe=0.500;
Td0=8.940;D=1;
Tj=5;ke=1;Te=0.05;
WN=314;km=1;Tm=0.5;
s=p+j*Q;I=conj(s/u);q=angle(I);
ut=u+I*(j*xe);qt=angle(ut);
Eqd=ut+j*xq*I;o=angle(Eqd);
id=abs(I)*sin(o-q);iq=abs(I)*cos(o-q);
ud=abs(u)*sin(o);uq=abs(u)*cos(o);
utd=abs(ut)*sin(o-qt);utq=abs(ut)*cos(o-qt);
epq=utq+xd'*id;
xpdm=xpd+xe;xqm=xq+xe;
%求线性化系数k1--k6
disp('线性化系数k1--k6')
k1=(epq*uq)/xpdm+((uq*uq-ud*ud)*(xd'-xq))/(xpdm*xqm)
k2=ud/xpdm
k3=(xd-xpd)/xpdm
k4=(ud*(xd-xpd))/xpdm
k5=(xq*utd*uq)/(abs(ut)*xqm)-(xpd*utq*ud)/(abs(ut)*xpdm)
k6=(utq*xe)/(abs(ut)*xpdm)
%求原系统矩阵A,B,C,D
disp('原系统各矩阵:')
A=[-(D/Tj),-(k1/Tj),1/Tj,0,-(k2/Tj);WN,0,0,0,0;0,0,-(1/Tm),0,0;0,0,0,-(1/Te),0;0,-(k4/Td0),0,1/Td0,-(1/(Td0*k3))]
B=[0,0;0,0;0,km/Tm;ke/Te,0;0,0]
C=[0,k5,0,0,k6;0,k1,0,0,k2]
D=0;
%求解原系统传递函数
sys=ss(A,B,C,D);
disp('原系统传递函数为：')
H=tf(sys)
%基于阶跃型输入无静差,构造增广系统（求系统矩阵）
     %确定内摸矩阵（Ac,Bc,Cc)
     disp('系统伺服系统传递函数矩阵为：')
      m=zeros(2);n=zeros(5,2);
      Ac=m
      Bc=[1,0;0,1]
      Cc=eye(2)
   %求增广系统Am,Bm，Mm(参考输入矩阵）
 disp('增广系统传递函数矩阵为：') 
 Am=[A,B;-Bc*C,Ac]
 Bm=[n;Bc]
 Cm=[C m]
%判断增广系统的稳定性
[b,c]=eig(Am);
disp('增广系统各极点')
r=real(diag(c))
re=(r>=0);
if re==0
    disp('由以上极点可知增广系统稳定')
else
    ree=r<=0;
    if ree==ones(size(ree))
        disp('由以上极点可知增广系统临界稳定')
    else
        disp('由以上极点可知增广系统不稳定')
    end
end   
%判断增广系统的可控性与可观性
Qc=cat(2,Bm,Am*Bm,Am^2*Bm,Am^3*Bm,Am^4*Bm,Am^5*Bm,Am^6*Bm);
Qo=cat(1,Cm,Cm*Am,Cm*Am^2,Cm*Am^3,Cm*Am^4,Cm*Am^5,Cm*Am^6);
rc=rank(Qc)
if rc==7
    disp('增广系统可控')
else
    disp('增广系统不可控')
end
ro=rank(Qo)
if ro==7
    disp('增广系统可观')
else
    disp('增广系统不可观')
end
%增广二次型最优控制
disp('增广二次型最优控制')
Qy=input('请设定输出加权对称常数矩阵Qy：');
Qy=diag(Qy)
Qu=input('请设定控制加权对称常数矩阵Qu：');
Qu=diag(Qu)
[K,P,E]=lqr(Am,Bm,Cm'*Qy*Cm,Qu); 
disp('黎卡提方程中矩阵P：')
P
disp('最优定常状态反馈矩阵K：')
K
disp('最优输出调节器的最优控制为：u*(t)=kx(t)')
disp('闭环系统的特征根：')
E
%状态观测器设计
length_of_Am=length(Am);
Om=obsv(Am,Cm);
rank_of_Om=rank(Om);
   if rank_of_Om==length(Am)
    %设置所需要的极点
   desired_poles=input('请设定期望极点：')
    %desired_poles=[-40 -41+0.32i -41-0.32i -42+1.25i -42-1.25i -43+13.48i -43-13.48i];
    %计算状态观测器增益矩阵L
    disp('状态观测器增益矩阵L:')
    L=place(Am',Cm',desired_poles)'
    disp('状态观测器系统状态空间模式表达式')
    Ao=Am-L*Cm
    Bo=[Bm L]
    Co=Cm
    Do=0
    else
      disp('不能设置状态观测器')
   end
disp('总的闭环系统的性质:')
   %带状态观测器的状态反馈系统
   R=zeros(2,7);
   %An=[Am,Bm*K;L*Cm,(Am-L*Cm+Bm*K)];
   An=[Am-Bm*K,Bm*K;zeros(7,7),(Am-L*Cm)];
   Bn=[Bm;Bm];
   Cn=[Cm R];
   Dn=0;
Qnc=cat(2,Bn,An*Bn,An^2*Bn,An^3*Bn,An^4*Bn,An^5*Bn,An^6*Bn,An^7*Bn,An^8*Bn,An^9*Bn,An^10*Bn,An^11*Bn,An^12*Bn,An^13*Bn);
Qno=cat(1,Cn,Cn*An,Cn*An^2,Cn*An^3,Cn*An^4,Cn*An^5,Cn*An^6,Cn*An^7,Cn*An^8,Cn*An^9,Cn*An^10,Cn*An^11,Cn*An^12,Cn*An^13);
rnc=rank(Qnc)
rno=rank(Qno)
if rnc==14
    disp('总的闭环系统可控')
else
    disp('总的闭环系统不可控')
end
if rno==14
    disp('总的闭环系统可观')
else
    disp('总的闭环系统不可观')
end
disp('  闭环传递函数： ')
sysk=ss(An,Bn,Cn,Dn);
Hk=tf(sysk)
%仿真结果
sim('f',60);%仿真模型名为f