控制系统数字仿真与CAD-第四次实验-附完整代码

1. clear;
2. clear all;
3. clear close;
4. sa=2;sb=20;l=sb-sa;
5. e=1e-5;
6. x1=sa+(1-0.618)*l;
7. x2=sa+0.618*l;
8. q1=obj_asr(x1);
9. q2=obj_asr(x2);
10. count=0;
11. x_dot=[x1,x2];
12. while((sb-sa)/l>e)
13. if(q1>q2)
14.         sa=x1;
15.         x1=x2;
16.         q1=q2;
17.         x2=sa+0.618*(sb-sa);
18.         x_dot=[x_dot,x2];
19.         q2=obj_asr(x2);
20. else
21.         sb=x2;
22.         x2=x1;
23.         q2=q1;
24.         x1=sa+(1-0.618)*(sb- sa);
25.         x_dot=[x_dot,x1];
26.         q1=obj_asr(x1);
27. end
28. end
29. Kn_opt=(x1+x2)/2;
30. q_min=obj_asr(Kn_opt);
31. opt_asr_plot(Kn_opt);

1. function [result]=fun4(Kn)
2. tn=0.03;
3. Ki=0.269;
4. ti=0.067;
5. Ks=76;
6. Ts=0.00167;
7. R=6.58;
8. Tl=0.018;
9. Tm=0.25;
10. Ce=0.031;
11. Alpha=0.00337;
12. Beta=0.4;
13. A=[0 0 1 1 0];
14. B=[tn ti Ts Tl Tm*Ce];
15. C=[Kn Ki Ks 1/R R];
16. D=[Kn*tn Ki*ti 0 0 0];
17. c=8;
18. r=1;
19. W=[0 0 0 0 -Alpha;
20. 1 0 0 -Beta 0;
21. 0 1 0 0 0;
22. 0 0 1 0 -Ce;
23. 0 0 0 1 0];
24. W0=[1 0 0 0 0 ]';
25. block_num=5;
26. h=0.001; 
27. t_end=0.5;
28. t=0:h:t_end;
29. for k=1:block_num
30.     if (A(k)==0) %求积分或积分比例环节各系数
31.         FI(k)=1;
32.         FIM(k)=h*C(k)/B(k);
33.         FIJ(k)=h*h*C(k)/B(k)/2; 
34.         FIC(k)=1;
35.         FID(k)=0;
36.         if(D(k)~=0) %若为积分比例，修正fai d
37.             FID(k)=D(k)/B(k);
38.         end
39.     else
40.         FI(k)=exp(-h*A(k)/B(k)); % 求惯性或惯性比例环节各系数 
41.         FIM(k)=(1-FI(k))*C(k)/A(k);
42.         FIJ(k)=h*C(k)/A(k)-FIM(k)*B(k)/A(k); 
43.         FIC(k)=1;
44.         FID(k)=0;
45.         if(D(k)~=0)
46.             FIC(k)=C(k)/D(k)-A(k)/B(k);
47.             FID(k)=D(k)/B(k);
48.         end
49.     end
50. end
51. Y0=[0 0 0 0 0]'; 
52. n=length(t);
53. Y=Y0;
54. X=zeros(block_num,1);
55. result=Y;
56. Uk=zeros(block_num,1);
57. Ub=Uk; 
58. for m=1:(n-1)
59.     Ub=Uk;
60.     Uk=W*Y+W0*r;
61.     Uf=2*Uk-Ub;
62.     Udot=(Uk-Ub)/h;
63. %     X=FI'.*X+FIM'.*Uk;
64. %     Y=FIC'.*X+FID'.*Uf;
65.     X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;
66.     Y=FIC'.*X+FID'.*Uf;
67.     Y(1)=satur(Y(1),c);
68.     result=[result,Y];
69. end
70. 
71. end

1. function q=obj_asr(Kn)
2. result = fun4(Kn);
3. y_inf = 296.7350;
4. q=100*(max(result(5,:))-y_inf)/y_inf; 
5. end

1. function q=opt_asr_plot(Kn)
2. result = fun4(Kn);
3. h=0.001; 
4. t_end=0.5;
5. t=0:h:t_end;
6. figure
7. subplot(4,1,1),plot(t,result(4,:)),grid,title('Id')
8. subplot(4,1,2),plot(t,result(3,:)),grid,title('Ud')
9. subplot(4,1,3),plot(t,result(2,:)),grid,title('ACR')
10. subplot(4,1,4),plot(t,result(1,:)),grid,title('ASR')
11. figure;
12. plot(t,result(5,:),'r');grid,title('n') 
13. end

1. function [uo]=satur(ui,c)
2. if (abs(ui)<=c)
3.     uo=ui;
4. elseif (ui>c)
5.     uo=c;
6. else
7.     uo=-c;
8. end
9. end