| 查看: 1192 | 回复: 1 | ||
samemo木虫 (正式写手)
|
[求助]
求教鲁棒控制DK迭代算法
|
小弟科研苦于无人讨论 ，只能上贵版求教下面是我编写的程序，不过不知道为何在最后一步dksyn算法这里要跑很久， 并且还跑不出来；特此，求教对鲁棒控制中DK迭代算法有经验的同学, 望提出问题所在，不 吝感激！  %%%%%%%%%%%%%%%%%%%%%%%%%%% clear clc tic N=4; % number of teeth xF=0.744; Kt=462e6; % tangential cutting force coefficient(N/m) Kn=38.6e6; % normal cutting force coefficient(N/m) w0x=2350*2*pi; % angular natural frequencyx(rad/s) zetax=0.05; % relative dampingx(1) w0y=2350*2*pi; % angular natural frequencyy(rad/s) zetay=0.05; % relative dampingy(1) m_tx=0.015; % massx(kg) m_ty=0.015; % massy(kg) w_ax=1400*2*pi; % angular natural frequencyx(rad/s) zeta_ax=0.12; % relative dampingx(1) w_ay=1400*2*pi; % angular natural frequencyy(rad/s) zeta_ay=0.12; % relative dampingy(1) m_ax=0.14; % massx(kg) m_ay=0.14; % massy(kg) k_ambx=1.2566e-5; %parameter of the AMB model(Nm^2/A^2) k_amby=1.2566e-5; i_0=2.5; %parameter of the AMB model(A) v_0=0.5e-3; %parameter of the AMB model(m) aD=1; % radial depth of cut up_or_down=-1; % 1: up-milling,-1: down-milling if up_or_down ==1 % up-milling fist=0; % start angle fiex=acos(1-2*aD); % exit anlge elseif up_or_down== -1 % down-milling fist=acos(2*aD-1); % start angle fiex=pi; % exit angle end stx=600; % steps of spindle speed sty=60; % steps of depth of cut w_st=0e-3; % starting depth of cut (m) w_fi=3e-3; % final depth of cut (m) o_st=10e3; % starting spindle speed (rpm) o_fi=40e3; % final spindle speed (rpm) % computational parameters k=40; % number of discretization interval over one period intk=20; % number of numerical integration steps for Equation (37) m=k; % since time delay=time period wa=1/2; % since time delay =time period wb=1/2; % since time delay =time period % numerical integration of specific cutting force coefficient for i=1: k dtr=2*pi/N/k; % t,if=2/N hxx(i)=0; hxy(i)=0; hyx(i)=0; hyy(i)=0; for j=1:N % loop for toothj for h=1:intk % loop for numerical integration of hi fi(h)=(i-1)*dtr+(j-1)*2*pi/N+h*dtr/intk; if (fi(h)>fist)*(fi(h)<fiex) g(h)=1; % tooth is in the cut Fz(h)=xF*((0.2)*sin(fi(h)))^(xF-1);%线性化之后的一项 elseif (fi(h)>fist+2*pi)*(fi(h)<fiex+2*pi) g(h)=1; % tooth is in the cut %Fz(h)=xF*((0.2)*sin(fi(h)))^(xF-1);%线性化之后的一项 else g(h)=0; % tooth is out of cut Fz(h)=0; end end hxx(i)=hxx(i)+sum(Fz.*g.*(Kt.*cos(fi)+Kn.*sin(fi)).*sin(fi))/intk; hxy(i)=hxy(i)+sum(Fz.*g.*(Kt.*cos(fi)+Kn.*sin(fi)).*cos(fi))/intk; hyx(i)=hyx(i)+sum(Fz.*g.*(-Kt.*sin(fi)+Kn.*cos(fi)).*sin(fi))/intk; hyy(i)=hyy(i)+sum(Fz.*g.*(-Kt.*sin(fi)+Kn.*cos(fi)).*cos(fi))/intk; end end A=zeros(8,8); % matrix A A(1,3)=1; A(2,4)=1; A(3,1)=-w0x^2*(m_tx/m_ax)-w_ax^2; A(3,2)=w0x^2*(m_tx/m_ax); A(3,3)=-2*zetax*w0x*(m_tx/m_ax)-2*zeta_ax*w_ax; A(3,4)=2*zetax*w0x*(m_tx/m_ax); A(4,1)=w0x^2; A(4,2)=-w0x^2; A(4,3)=2*zetax*w0x; A(4,4)=-2*zetax*w0x; A(5,7)=1; A(6,8)=1; A(7,5)=-w0y^2*(m_ty/m_ay)-w_ay^2; A(7,6)=w0y^2*(m_ty/m_ay); A(7,7)=-2*zetay*w0y*(m_ty/m_ay)-2*zeta_ay*w_ay; A(7,8)=2*zetay*w0y*(m_ty/m_ay); A(8,5)=w0y^2; A(8,6)=-w0y^2; A(8,7)=2*zetay*w0y; A(8,8)=-2*zetay*w0y; Bt=zeros(8,2); % matrix Bt Bt(4,1)=1/m_tx; Bt(8,2)=1/m_ty; Ba=zeros(8,2); % matrix Ba Ba(3,1)=1/m_ax; Ba(7,2)=1/m_ay; Ct=zeros(2,8);%matrix Ct Ct(1,2)=1; Ct(2,6)=1; Ca=zeros(2,8);%matrix Ca Ca(1,1)=1; Ca(2,5)=1; k_ambx=1.2566e-5; %parameter of the AMB model(Nm^2/A^2) k_amby=1.2566e-5; i_0=2.5; %parameter of the AMB model(A) v_0=0.5e-3; %parameter of the AMB model(m) Ki=zeros(2,2);%matrix K_i Ki(1,1)=4*k_ambx*i_0/(v_0^2); Ki(2,2)=4*k_amby*i_0/(v_0^2); Ks=zeros(2,2);%matrix K_s Ks(1,1)=4*k_ambx*i_0^2/(v_0^3); Ks(2,2)=4*k_amby*i_0^2/(v_0^3); %the rational transfer functions G_d,where the size of A_d depends on the order of the pade approximation nmin=3e4;nmax=3.2e4; %controllers designed for range of spindle speed（rpm） tau_min=60/nmax/N;tau_max=60/nmin/N;%N for number of teeth; t_min for delay s=tf('s');%pade approximation sys=exp(-0.5*(tau_min+tau_max)*s); sysx=pade(sys,10);% the pade approximation for 10 order T=ss(sysx); A_d(1:10,1:10)=T.a;A_d(11:20,11:20)=T.a;B_d(1:10,1)=T.b;B_d(11:20,2)=T.b; C_d(1,1:10)=T.c;C_d(2,11:20)=T.c;D_d(1,1)=T.d;D_d(2,2)=T.d;%pade approximation %the rational transfer functions W_d num_wd=[(tau_max-tau_min)/2 0]; den_wd=[(tau_max-tau_min)/2/3.456 1]; W_wd=tf(num_wd,den_wd);%the weighting function W_wd(s) [A_w1,B_w1,C_w1,D_w1]=tf2ss(num_wd,den_wd); A_w(1,1)=A_w1;A_w(2,2)=A_w1;B_w(1,1)=B_w1;B_w(2,2)=B_w1; C_w(1,1)=C_w1;C_w(2,2)=C_w1;D_w(1,1)=D_w1;D_w(2,2)=D_w1; %the weighting function Kp=0.5e-6; %the gain of the weighting function(m/A) frh=7500; %the roll-off frequency(Hz) fph=1e5; %a pole frequency(Hz) num_KS=Kp*[1/(2*pi*frh) 1]; %num=Kp*[fph 2*pi*frh*fph]; den_KS=[1/(2*pi*fph) 1]; %den=[frh 2*pi*frh*fph]; W_KS=tf(num_KS,den_KS);%the weighting function W_KS(s) [A_KS1,B_KS1,C_KS1,D_KS1]=tf2ss(num_KS,den_KS); %T=ss(W_KS);T=ss(tf(num,den));A _KS=T.a;B_KS=T.b;C_KS=T.c;D_KS=T.d; A_KS(1,1)=A_KS1;A_KS(2,2)=A_KS1;B_KS(1,1)=B_KS1;B_KS(2,2)=B_KS1; C_KS(1,1)=C_KS1;C_KS(2,2)=C_KS1;D_KS(1,1)=D_KS1;D_KS(2,2)=D_KS1; %start to design controller K a_p=2e-3; %a_p is the maximal depth of cut for which stable cutting is desired( m) H=zeros(2,2); % average(zero-order) component of the Fourier sries expans ion according to Equations (20)–(21) H(1,1)=-1/k*sum(hxx);H(1,2)=-1/k*sum(hxy); H(2,1)=-1/k*sum(hyx);H(2,2)=-1/k*sum(hyy);%负号添加在H里面 %the state-space matrices of the generalized plant Ap(1:8,1:8)=A+Ba*Ks*Ca+0.5*a_p*Bt*H*(Ct-D_d*Ct); Ap(1:8,9:28)=-0.5*a_p*Bt*H*C_d;%a_p is the maximal depth of cut for which stable cutting is desired(m) Ap(9:28,1:8)=B_d*Ct; Ap(9:28,9:28)=A_d; Ap(29:30,1:8)=B_w*Ct; Ap(29:30,29:30)=A_w; Ap(31:32,31:32)=A_KS; Bp(1:8,1:2)=-0.5*a_p*Bt*H; Bp(1:8,3:4)=Bt; Bp(1:8,7:8)=Ba*Ki; Bp(31:32,7:8)=B_KS; Cp(1:2,1:8)=D_w*Ct; Cp(1:2,29:30)=C_w; Cp(3:4,1:8)=0.5*a_p*H*(Ct-D_d*Ct); Cp(3:4,9:28)=-0.5*a_p*H*C_d; Cp(5:6,31:32)=C_KS; Cp(7:8,1:8)=Ca; Dp(3:4,1:2)=-0.5*a_p*H; Dp(5:6,7:8)=D_KS; Dp(7:8,5:6)=[1 0;0 1]; for iu=1:8 [num_p((iu-1)*8+1:8*iu,  ,den_p(iu,]=ss2tf(Ap,Bp,Cp,Dp,iu);%将8输入8输出状态空间转换为传递函数，num_p为分子系数 end for g1=1:8 for g2=1:8 G(g2,g1)=tf(num_p(g2+(g1-1)*8,: ),den_p(g1, );%将传递函数表达出来end end delt_d=ucomplex('delt_d',0);delt_ap=ucomplex('delt_ap',0);%建立以0为圆心， 1为半径的不确定复数 delt_p(1,1)=delt_d;delt_p(2,2)=delt_d;delt_p(3,3)=delt_ap;delt_p(4,4)=delt _ap;%建立对角线上的模都小于1的不确定复数矩阵 I4=eye(4); A11=I4-[G(1,1),G(1,2),G(1,3),G(1,4);G(2,1),G(2,2),G(2,3),G(2,4);G(3,1),G(3 ,2),G(3,3),G(3,4);G(4,1),G(4,2),G(4,3),G(4,4)]*delt_p; A12=[G(1,5),G(1,6),G(1,7),G(1,8);G(2,5),G(2,6),G(2,7),G(2,8);G(3,5),G(3,6) ,G(3,7),G(3,8);G(4,5),G(4,6),G(4,7),G(4,8)]; P = [G(5,5),G(5,6),G(5,7),G(5,8);G(6,5),G(6,6),G(6,7),G(6,8);G(7,5),G(7,6 ),G(7,7),G(7,8);G(8,5),G(8,6),G(8,7),G(8,8)]+... [G(5,1),G(5,2),G(5,3),G(5,4);G(6,1),G(6,2),G(6,3),G(6,4);G(7,1),G(7,2) ,G(7,3),G(7,4);G(8,1),G(8,2),G(8,3),G(8,4)]*delt_p*inv(A11)*A12; %含有不确定性的[r i_c]到[z y]的传递函数；x,y两个方向，r和i_c都是2阶向量 [K,clp,bnd] = dksyn(P,2,2); %D-K iteration %%%%%%%%%%%%%%%%%%%%%%%%%%% 系统的对象图 http://dl.vmall.com/c01fse20t8 http://dl.vmall.com/c0ovxpsxdy [ Last edited by samemo on 2014-3-2 at 13:55 ] |
回复此楼» 猜你喜欢
081700，311，求调剂
已经有15人回复
-
一志愿北京化工085600 310分求调剂
已经有18人回复
-
085600材料与化工301分求调剂院校
已经有4人回复
-
材料与化工371求调剂
已经有14人回复
-
336材料与化工085600求调剂
已经有7人回复
-
材料334求调剂
已经有18人回复
-
331求调剂
已经有8人回复
-
332求调剂
已经有17人回复
-
一志愿南京航空航天大学 材料与化工329分求调剂
已经有4人回复
-
材料专硕322
已经有7人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- 关于迭代算法的收敛问题，求高人指点！
已经有3人回复
- matlab中利用LM算法解非线性方程组的问题
已经有5人回复
- 数学好的请进，请问这叫什么迭代算法？有具体描述。非线性迭代求解
已经有26人回复
samemo
木虫 (正式写手)
- 应助: 1 (幼儿园)
- 金币: 4089.5
- 散金: 20
- 红花: 4
- 帖子: 668
- 在线: 272.8小时
- 虫号: 1395679
- 注册: 2011-09-08
- 性别: GG
- 专业: 控制理论与方法
2楼2014-03-03 08:31:06
