自控系统数字仿真实验

1、实验一  面向微分方程的数值积分法仿真一、实验目的1掌握数值积分法的基本概念、原理及应用；2用龙格-库塔法解算微分方程，增加编写仿真程序的能力；3分析数值积分算法的计算步长与计算精度、速度、稳定性的关系；4. 对数值算法中的“病态问题”进行研究。二、实验内容1、已知系统微分方程及初值条件取步长h=0.1 ，试分别用欧拉方程法和RK4法求t=2h时的y值，并将求得的值与解析解比较(将三个解绘于同一坐标中，且用数值进行比较)，说明造成差异的原因。（编程完成；选用MATLAB ode函数完成。） 程序如下：t0=0;tf=2;h=0.1;y1=1;y2=1;y

2、3=1;t1=0;t2=0;t3=0;n=round(tf-t0)/h;for i=1:n %欧拉法 y1(i+1)=y1(i)+h*(2*h+y1(i); t1=t1,t1(i)+h;endfor i=1:n %四阶龙格-库塔法 k1=y2(i)+t2(i); k2=y2(i)+h*k1/2+t2(i)+h/2; k3=y2(i)+h*k2/2+t2(i)+h/2; k4=y2(i)+h*k3+t2(i)+h; y2(i+1)=y2(i)+h*(k1+2*k2+2*k3+k4)/6; t2=t2,t2(i)+h;endfor i=1:n y3(i+1)=2*exp(t3(i)-t3(i)-1

3、; t3=t3,t3(i)+h;endplot(t1,y1,'r',t2,y2,'g',t3,y3,'k')分析:红线为用欧拉法得到的结果，绿线为用四阶龙格库塔法得到的结果，蓝线为根据解析方程得到的结果。其差异原因主要有两个：（1）二者的方法不同，欧拉法是根据一阶微分方程计算得到的，龙格库塔法是根据四阶微分方程得到的；（2）由于步长取为0.1，所以得到的图像与解析解之间存在差异，若将步长取小，则得到的解将更靠近解析解。  2、已知系统的传递函数为在单位阶跃输入下，系统响应的解析解为 试分别用欧拉方程法和RK4法对系统进行仿真（

4、编程完成）：1）比较两种数值积分解与解析解的逼近程度；（绘图）2）改变步长，分析步长对数值解精度的影响；3）不断加大步长，分析计算稳定性的变化。程序如下： a=0 1 0; 0 0 1; -22.06 -27 -10;b=0;0;1; c=40.6 0 0;X1=0;0;0; t=0;Y1=0;X=0; u=1; Y2=0; Y3=0;X2=0;0;0; x=0; h=0.1;t0=0; tf=2; t1=0; t2=0; t3=0;N=(tf-t0)/h;for i=1:N k1=a*X1+b; k2=b+a*(h*k1/2+X1); k3=b+a*(h*k2/2+X1); k4=b+a*(

5、h*k3+X1); X1=X1+h*(k1+2*k2+2*k3+k4)/6; Y1=Y1,c*X1; t1=t1,t1(i)+h;endfor i=1:N x=X2(:,i)+h*(a*X2(:,i)+b*u); y=c*x; X2=X2,x; Y2=Y2,y; t2=t2,t2(i)+h;endfor i=1:N y=1.84-4.95*i*exp(-1.88*i)-1.5*exp(-1.88*i)-0.34*exp(-6.24*i); Y3=Y3,y; t3=t3,t3(i)+h;endplot(t1,Y1,'r',t2,Y2,'g',t3,Y3,'

6、b') 当h=0.01时，仿真图如下： 可见当步长很小时，数值积分解逼近解析解，其数值解得的精度也很高。 当h=0.1时，仿真图如下： 可见随着步长的增加，其数值积分解与解析解的逼近程度减小，精度降低。 当h=0.3时，仿真图如下： 可见其数值积分解已变得不稳定了，即步长不能取得太大3、求下图所示系统的阶跃响应y(t)的数值解。（v=1,k=1,开始时间t0=0,结束时间tf=10,h=0.25）分析k、v对系统响应的影响。（编程用RK4求解；Simulink）(修正：G(s)表达式的分母中，0.25s+1改为（0.25s+1）的二次方)    程序如下： k=1;

7、 a=conv(1 0 0,conv(0.25 1,0.25 1) b=2*k k; X0=0 0 0 0 ;v=1; n0=4; tf=10; h0=0.25; r=1; V=v; n=n0; T0=0; Tf=tf; h=h0; R=r; b=b/a(1); a=a/a(1); A=a(2: n+1); %首一化处理 A=rot90(rot90(eye(n-1,n); -fliplr(A); %形成能控标准型A阵 B=zeros(1,n-1),1' %形成输入阵B(n维列向量)m1=length(b); %分子系数向量维数m+1 C=fliplr(b),zeros(1,n-m1);

8、 %形成输出阵C(n维行向量) Ab=A-B*C*V; %形成闭环系数阵AbX=X0'y=0;t=T0; %赋初值，准备开始递推运算 N=round(Tf-T0)/h; %确定输出点数for i=1:N %四阶龙格-库塔法 k1=Ab*X+B*R; k2=Ab*(X+h*k1/2)+B*R; k3=Ab*(X+h*k2/2)+B*R; k4=Ab*(X+h*k3)+B*R; %求各阶次斜率 X=X+h*(k1+2*k2+2*k3+k4)/6; %求状态 y=y,C*X; %求输出并以向量形式保存 t=t,t(i)+h; %输出对应时刻以向量形式保存endt',y' %输

9、出数据形式结果plot(t,y) %输出曲线形式结果 k=0.5 K=1，v=1 K=1，V=0.5 K=1,v=5时， K=2，v=1时分析：当k取值增大，v值不变时，系统输出的波头增多，而且也变陡，稳态精度降低，当k增加到一定程度时系统便发散了（即不稳定了）。当v值增大，k值不变时，波头也是变多变陡，当v值增大到一定程度时系统便不稳定了。4、已知系统状态状态方程为（修正：x的系数矩阵中最后一个元素为-40，而非40）采用RK4法，步长分别取h=0.01、0.04、0.06，求系统的零输入响应，并绘图分析各状态变量的响应状态及产生的原因。（提示：病态系统）程序如下：a=-21 19 -20;

10、19 -21 20;40 -40 -40;x=1;0;-1;X=x;t=0;t0=0;tf=2;h=0.01;n=round(tf-t0)/h;for i=1:n x=X(:,i)+h*a*X(:,i); X=X,x; t=t,t(i)+h;endb=X(1,:);c=X(2,:);d=X(3,:);plot(t,b,'r',t,c,'g',t,d,'b') h=0.01时 H=0.02时 H=0.04时分析：如图，当h=0.01时，在t=0.2s以后系统输出便趋于平稳，当取h=0.02时，系统输出振荡剧烈，趋于稳定的时间也变长，当取h=0.04

11、后，系统输出呈发散振荡形式。当h=0.06后系统仍然是发散的，即当h的取值改变时，原先稳定的系统变得不稳定了，这便是病态系统。但是这个结果与书上的不同。三、思考题1、在进行仿真计算时，是否选用的数值积分法的阶次越高越好？ 不是，因为阶次越高，虽然精度高了，但其计算的复杂程度也越高，有时甚至无法完成计算。 2、选用数值积分法进行仿真的原则。(1) 精度 仿真结果的精度主要受三项误差影响： 截断误差：由算法本身的精度阶次所决定； 舍入误差：由计算机字长决定； 累积误差：由以上两项误差随计算时间长短累积情况决定。(2) 计算速度 取决于所用数值方法和计算步长。(3) 稳定性 数值稳定性主要

12、与计算步长h有关，不同的数值方法对h都有不同的稳定性限制范围，且与被仿真对象的时间常数有关。一般所选步长与系统最小时间常数有以下数量级关系： 而多步法、隐式算法有较好的数值稳定性，在对稳定性较注重时，应予以优先选择。 实验二第一部分  面向结构图的数值积分法仿真一、实验目的加深理解连续系统面向结构图仿真的原理及特点，进一步掌握数字积分法解算微分方程的方法，增加编写仿真程序的能力。二、实验内容1、用面向方框图的数字仿真方法对下列系统进行仿真。    2、求解下图所示系统在f=-1(t)阶跃扰动作用下第、第环节的动态过程。分别用面向框图的

13、数值积分法（RK4法）、MATLAB中有关系统建模的命令和Simulink三种方法求解。三实验程序及结果1.P=0,0.07,1,0.14;1,0.012,1,0;0,0.05,1,0.15;10,1,1,0;1,0.01,0,0.0008;WIJ=1,0,1;1,4,-1;2,1,1;3,2,1;3,5,-1;4,3,1;5,4,1n=5;Y0=1;Yt0=0 0 0 0 0;h=0.01;t=0;L1=5;T0=0;Tf=2;nout=4;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0

14、=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0) W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0'y=Y(nout);t=T0;N=round(Tf-T0)/(h*L1);for i=1:N; for j=1:L1; k1=Ab*Y+b1*1; k2=Ab*(Y+h*k1/2)+b1*1; k3=Ab*(Y+h*k2/2)+b1*1; k4=

15、Ab*(Y+h*k3)+b1*1; Y=Y+h*(k1+2*k2+2*k3+k4)/6; end; y=y,Y(nout); t=t,t(i)+L1*h; end; plot(t,y)2.P=1 0.0149 1 0; 1 0.00254 26.6667 0; 1 0.391 0 0.4199; 0 0.248 1 0; 1 2.88 1 0;WIJ=1 5 1; 2 1 1; 3 4 1; 4 2 1; 4 3 -1; 5 0 1; 5 4 -1;n=5;Y0=-1;Yt0=0 0 0 0 0;h=0.005;L1=10;T0=0;Tf=10;nout=5;A=diag(P(:,1);B=d

16、iag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0'y=Y(nout);t=T0;N=round(Tf-T0)/(h*L1);for i=1:N; for j=1:L1; K1=A

17、b*Y+b1*Y0; K2=Ab*(Y+h*K1/2)+b1*Y0; K3=Ab*(Y+h*K2/2)+b1*Y0; K4=Ab*(Y+h*K3)+b1*Y0; Y=Y+h*(K1+2*K2+2*K3+K4)/6; end;y=y,Y(nout);t=t,t(i)+h*L1;end;t',y'plot(t,y)grid on;xlabel('t/s');ylabel('Amplitude');axis(0 10 -0.08 0.02);结果：第二部分  面向结构图的离散相似法仿真一、实验目的   1掌握离散

18、相似法的基本概念和原理；2典型环节离散系数的求取及差分方程表示；3. 掌握非线性系统的数字仿真方法。二、实验内容1、已知控制系统结构图如图所示，设输入阶跃函数幅值Y0=10，滞环非线性参数s=1(滞环宽度），请用离散相似法编程和Simulink法对系统进行如下分析：1）不考虑非线性环节影响时，求解y(t)的阶跃响应；2）考虑非线性环节影响，其余参数不变，求解y(t)并与线性情况所得结果进行比较；3）改变的滞环非线性参数s，分析该非线性对系统的影响。2、系统结构图如图所示，先理论分析该系统是否会产生自激振荡，若会求出振荡的振幅和频率，并用离散相似法编程和Simulink法验证分析的结果

19、。三实验程序及结果1.(1)Simulink法结果：(2)MATLAB法P=1 10 5 10;1 0.5 1 0;1 0.1 1 0;0 1 1 0;WIJ=1 0 1;2 1 1;3 2 1;4 3 1;1 4 -1;X0=0 0 0 0;Z=0 0 0 6;S=0 0 0 1; h=0.01;L1=25;n=4;T0=0;Tf=10;nout=4;Y0=10;A=(P(:,1);B=(P(:,2);C=(P(:,3);D=(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0);W0(WI

20、J(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endendfor i=1:n if(A(i)=0); FI(i)=1; FIM(i)=h*C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0); FID(i)=D(i)/B(i); else end else FI(i)=exp(-h*A(i)/B(i); FIM(i)=(1-FI(i)*C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;FID(i)=

21、0; if(D(i)=0); FIM=(1-FI(i)*D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endend Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;N=length(t);for k=1:N-1 for l=1:L1 Ub=Uk; Uk=W*Y+W0*Y0; for i=1:n if(Z(i)=0) if(Z(i)=1) Uk(i)=satu(Uk(i),S(

22、i); end if(Z(i)=2) Uk(i)=dead(Uk(i),S(i); end if(Z(i)=3) Uk(i),Ubb(i)=backlash(Ubb(i),Uk(i),Ub(i),S(i); end end end Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; for i=1:n if(Z(i)=0) if(Z(i)=4) Y(i)=satu(Y(i),S(i); end if(Z(i)=5) Y(i)=dead(Y(

23、i),S(i); end if(Z(i)=6) Y(i),Ubb(i)=backlash(Ubb(i),Y(i),Yb(i),S(i); end end end end y=y,Y(nout);endplot(t,y)结果：S4=1S4=52. (1)Simulink法： (2)Matlab法P=0,1,1,0;1,1,1,0;2,1,1,0;WIJ=1 0 1;1 3 -1;2 1 1;3 2 1;n=3;Y0=10;X0=0 0 0 0;Yt0=0 0 0 0;h=0.01;L1=25;t0=0;tf=10;nout=3;Z=0 0 6;S=0 0 1;A=(P(:,1);B=(P(:,2

24、);C=(P(:,3);D=(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;sp3_3sp3_4plot(t,y)实验三 采样控制系统的数字仿真一、实验目的  1、掌握采样控制系统数字仿真的一般方法；2、掌握采样周期与计算步长的关系与确定方法。二、实验内容1、已知采样系统结构如图所示，分别用程序设计方法和Simulink法求系统的输出响应

25、。（取仿真时间:Tf=10；采样周期:T=0.2；计算步长：h=0.01）(注：第一个框图中，分子分母中z的指数均为-1；第二个框图中e的指数为-Ts)2、设某数字控制系统如图所示，采样周期为T=0.1s，初始状态。(注：第二个框图中e的指数为-Ts) （1）数字控制器为PI调节器 ，其中（2）数字控制器为PID调节器，其中要求： 1）用程序设计法、MATLAB控制工具箱时域响应分析函数、Simulink法求系统在单位阶跃输入信号作用下的输出响应；  2）比较两种控制器作用下系统的响应，由此得出微分调节的作用。三，程序及结果1.（1）matlab法：G=2.72

26、-1;F=0.717;P=0 1 1 0; 1 1 1 0;WIJ=1 0 1; 2 1 1;n=2;Y0=1;X0=0 0; h=0.01;T0=0;Tf=10;Ts=0.2;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n); for k=1:m if(WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endend for i=1:n if(A(i)=0) FI(i)=1; FIM

27、(i)=h*(C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0) FID(i)=D(i)/B(i); else end else FI(i)=exp(-h)*A(i)/B(i); FIM(i)=(1-FI(i)*(C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;FID(i)=0; if(D(i)=0) FIM(i)=(1-FI(i)*(D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)*A(i); FIC(i)=C(i)

28、/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(n,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=length(t0);t1=T0:Ts:Tf;M=length(t1);for k=1:M-1; U=uk; ek=Y0-x2; E=ek;E(1); uk=-F*U+G*E; for l=1:N-1 Ub=Uk; Uk=W*Y+W0*uk; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI

29、'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; y(N-1)*(k-1)+l),1)=Y(nout); end x2=Y(nout); endplot(1:1000,y);结果：2.（1）G=15.2 -15;F=-1;P=1 1 1 0; 1 4 1 0;WIJ=1 0 1; 2 1 1;n=2;Y0=1;X0=0 0; h=0.01;L1=25;T0=0;Tf=10;Ts=0.1;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,

30、1);W0=zeros(n,1);W=zeros(n,n); for k=1:m if(WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endend for i=1:n if(A(i)=0) FI(i)=1; FIM(i)=h*(C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0) FID(i)=D(i)/B(i); else end else FI(i)=exp(-h)*A(i)/B(i); FIM(i)=(1-FI(i)*(C

31、(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;FID(i)=0; if(D(i)=0) FIM(i)=(1-FI(i)*(D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)*A(i); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(2,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=l

32、ength(t0);t1=T0:Ts:Tf;M=length(t1);for k=1:M-1; U=uk; ek=Y0-x2; E=ek;E(1); uk=-F*U+G*E; for l=1:N-1 Ub=Uk; Uk=W*Y+W0*uk; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; y(N-1)*(k-1)+l),1)=Y(nout); end x2=Y(nout); endt=0:0.01:9.99;plot(t,y,'

33、-');grid on;xlabel('t/s');ylabel('Amplitude');hold on（2）matlab法G=55.2 -95 40;F=-1;P=1 1 1 0; 1 4 1 0;WIJ=1 0 1; 2 1 1;n=2;Y0=1;X0=0 0; h=0.01;L1=25;T0=0;Tf=10;Ts=0.1;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n); for k=1:m if(WIJ(k,2)=0);

34、W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endend for i=1:n if(A(i)=0) FI(i)=1; FIM(i)=h*(C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0) FID(i)=D(i)/B(i); else end else FI(i)=exp(-h)*A(i)/B(i); FIM(i)=(1-FI(i)*(C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;

35、FID(i)=0; if(D(i)=0) FIM(i)=(1-FI(i)*(D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)*A(i); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(3,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=length(t0);t1=T0:Ts:Tf;M=length(t1);for k=1:M-1; U=uk; ek

36、=Y0-x2; E=ek;E(1:2); uk=-F*U+G*E; for l=1:N-1 Ub=Uk; Uk=W*Y+W0*uk; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; y(N-1)*(k-1)+l),1)=Y(nout); end x2=Y(nout); endt=0:0.01:9.99;plot(t,y);grid on;xlabel('t/s');ylabel('Amplitude');

37、实验四  数字仿真技术的综合应用一、 实验目的使学生在系统地学习了计算机仿真的基本概念和方法后，对计算机仿真的全过程有一个清楚的概念和认识，能够将所学的知识应用到实际系统中去。二、实验内容    选用经典控制理论或现代控制理论中某一种你所熟知的控制器设计方法，为下列某一系统（三选一）进行控制器设计：  （三）直流电动机调速系统 1 系统物理模型、参数及要求1.1 系统物理模型本组设计的被控对象为前一课程设计直流拖动自动控制系统所采用的直流电机，其物理模型如下图1所示。图1 双闭环控制电流调速系统的特点是电机的转速和电流分别由两个独立的调节器分别

38、控制，且转速调节器的输出就是电流调节器的给定，因此电流环能够跟随转速的偏差调节电机电枢的电流。当转速低于给定转速时，转速调节器的积分作用使输出增加，即电流给定上升，并通过电流环调节使电机电枢电流增大，从而使电机获得加速转矩，电机转速上升。当实际转速高于给定转速时，转速调节器的输出减小，即给定电流减小，并通过电流环调节使电机电枢电流下降，电机将因为电磁转矩减小而减速。在当转速调节器饱和输出达到限幅值时，电流环即以最大电流限制实现电机的加速，使电机的启动时间最短。双闭环调速系统的原理框图如图2所示。图21.2 系统参数电机型号：DJ15minrad=n /1600 额定参数： ， ， ， ， 电枢

39、电阻：R=29，电枢电感：L=0.89mH电机飞轮惯量：GD=0.045N/m电枢回路电磁时间常数：T=0.03s, 系统的机电时间常数：T=0.034s电动机电势时间常数：C=0.105， 转矩常数：C=1.0Vmin/r电流反馈系数：=5V/A, 转速反馈系数：=0.003V/(rpm) 1.3 要求设计状态反馈控制器，使得系统的单位阶跃响应性能指标为：（1）调节时间小于2秒； （2）系统的超调量小于5%； （3）稳态误差小于1%。2 系统模型的建立2.1 模型抽象直流电机转矩和电枢电流的关系为： 电枢旋转产生反电动势e与旋转运动角速度的关系为： 根据牛顿第二定律列写运动平衡方程式为： 其

40、中b为电机摩擦系数，此处忽略不计。根据回路电压法列写电机电枢回路方程式为： 由于： ，可得： ， 其中，m为一个旋转体上的一个质点的质量，质量m为该质量的重量G和重力加速度g之比，R和D分别为旋转体的半径和直径，综合上两式可得： 从而可以得到电机电枢回路电压平衡和电机运动平衡的一组微分方程式 其中，摩擦系数 =b/9.55,此处忽略不计。设系统的状态变量为： ，以输入电压u为输入，转速n为输出。建立系统状态空间表达式为： ïïîïïíìúûùêëé=ú&#

41、251;ùêëé+úûùêëéúûùêëé=úûùêëé21212110012.10118.0-3.833358.32-xxyxxxx&&u由式带入数据进行计算化简可得：其中u代表给定电压值，代表系统内部电流变化，代表系统输出转速变化，y=x表示系统的转速输出。可得：A=-32.58 -0.118;8333.33 0; B=1.12;0; C=0 1; D=0。

42、2.2 所建模型的性能分析2.2.1 能空能观测性分析(1) 能空性分析MATLAB程序如下： A=-32.58 -0.118;8333.33 0;B=1.12;0;C=0 1;D=0;n=length(A) %求系统阶次qc=ctrb(A,B),nc=rank(qc) %求能控性判别矩阵及其秩if n=nc,disp('系统是可控的!'), %判断系统的能控性else disp('系统是不可控的!'),end运行结果为：>> 原系统的秩为：n = 2qc = 1.0e+003 * 0.0011 -0.0365 0 9.3333nc = 2系统是可控

43、的! (2) 能观测性分析MATLAB程序如下：A=-32.58 -0.118;8333.33 0;B=1.12;0;C=0 1;D=0;disp('原系统的秩为：');n=length(A) %求系统阶次qo=obsv(A,C),no=rank(qo) %求系统能观测性判别矩阵及其秩if n=no,disp('系统是可观测的!'), %判断系统的能观测性else disp('系统是不可观测的!'),end运行结果为：>>原系统的秩为：n = 2qo = 1.0e+003 * 0 0.0010 8.3333 0no = 2系统是可观测

44、的!2.2.2 系统稳定性分析(1) 利用MATLAB程序分析系统稳定性是系统正常工作的首要条件。只要系统的状态矩阵A的特征根全部具有负实部，系统就是状态稳定的。MATLAB程序如下:A=-32.58 -0.118;8333.33 0; B=1.12;0; C=0 1; D=0;z,p,k=ss2zp(A,B,C,D,1); Flags=0; n=length(A);for i=1:n; if real(p(i)>0; Flags=1; endenddisp('系统的零极点模型为:')z,p,kif Flags=1 disp('系统不稳定');else disp('系统是稳定的');end运行结果为：>> 系统的零极点模型为:z = Empty matrix: 0-by-1p = -16.2900+26.7949i -16.2900-26.7949ik = 9.3333e+00系统是稳定的(2) 利用零极点图分析MATLAB程序如下:A=-32.58 -0.118;8333.33 0;B=1.12;0;C