24小时热门版块排行榜

返回列表

samemo

木虫 (正式写手)

应助: 1 (幼儿园)
金币: 4094.5
散金: 20
红花: 4
帖子: 667
在线: 272.6小时
虫号: 1395679
注册: 2011-09-08
性别: GG
专业: 控制理论与方法

[求助] 求教鲁棒控制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 ]

回复此楼