diff --git a/matlab/identifyFxFyStage.m b/matlab/identifyFxFyStage.m index e920648..32ad10c 100644 --- a/matlab/identifyFxFyStage.m +++ b/matlab/identifyFxFyStage.m @@ -14,23 +14,21 @@ function [mot1,mot2]=identifyFxFyStage() obj.currstep=iddata([zeros(10,1); obj.currstep.data(:,2)],[zeros(10,1); obj.currstep.data(:,3)],50E-6); f=load(strcat(path,sprintf('full_bode_mot%d.mat',motid))); - obj.frq=f.frq*2*pi; %convert from Hz to rad/s + obj.w=f.frq*2*pi; %convert from Hz to rad/s + if motid==2 + f.db_mag(1:224)=f.db_mag(225); % reset bad values at low frequencies + end obj.mag=10.^(f.db_mag/20); %mag not in dB obj.phase=f.deg_phase*pi/180; %phase in rad response = obj.mag.*exp(1j*obj.phase); - obj.meas= idfrd(response,obj.frq,0); + obj.meas= idfrd(response,obj.w,0); end function tfc=currstep(obj) opt=tfestOptions; opt.Display='off'; tfc = tfest(obj.currstep, 2, 0,opt); - den=tfc.Denominator; - num=tfc.Numerator; - k=num(1); - w0=sqrt(num(1)); - damp=den(2)/2/w0; - s=sprintf('k:%g w0:%g damp:%g\n',k,w0,damp); + s=str2ndOrd(tfc); %disp(s); %h = stepplot(tf1); %l=obj.currstep.OutputData @@ -45,6 +43,36 @@ function [mot1,mot2]=identifyFxFyStage() setoptions(h,'FreqUnits','Hz','Grid','on'); end + function s=str2ndOrd(tf) + den=tf.Denominator; + num=tf.Numerator; + k=num(1)/den(3); + w0=sqrt(den(3)); + damp=den(2)/2/w0; + s=sprintf('k:%g w0:%g damp:%g\n',k,w0,damp); + end + + function tf=fyModel() + num=[ 2.36527033e+11, 1.17108082e+13, 3.62387303e+17]; + den=[ 1.00000000e+00, 6.64495206e+03, 2.12777376e+07, ... + 1.23728427e+10, 3.07054470e+13, 1.72592127e+15, ... + 5.39888656e+17]; + + tf=idtf(num,den); + end + + function tf=fxModel() + num=[ 1.23284092e+11, 4.14791803e+13, 1.18702926e+18, ... + 2.96296718e+20, 2.67179357e+24, 4.04662786e+26, ... + 1.59131515e+30, 1.02778572e+32, 1.64551888e+35]; + den=[ 1.00000000e+00, 6.93892369e+03, 3.17041055e+07, ... + 7.66104262e+10, 2.36504992e+14, 2.23054854e+17, ... + 5.12578678e+20, 2.04416512e+23, 3.27771400e+26, ... + 4.77145416e+28, 3.85452959e+31, 1.28911178e+33, ... + 9.52157664e+34]; + tf=idtf(num,den); + end + function mot=fyStage() mot=loadData('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/MXTuning/18_10_02/',1); mot.tfc=currstep(mot); @@ -52,12 +80,13 @@ function [mot1,mot2]=identifyFxFyStage() opt=tfestOptions; opt.Display='off'; opt.initializeMethod='iv'; - opt.WeightingFilter=[1,5;20,670]*(2*pi); % Hz->rad/s conversion + opt.WeightingFilter=[1,5;30,670]*(2*pi); % Hz->rad/s conversion figure(); - mot.tf4_2 = tfest(mot.meas, 4, 2, opt); - mot.tf6_4 = tfest(mot.meas, 6, 4, opt); - h=bodeplot(mot.meas,'r',mot.tf4_2,'b',mot.tf6_4,'g'); + mot.tf2_0 = tfest(mot.meas, 2, 0, opt);disp(str2ndOrd(mot.tf2_0)); + mot.tf_mdl= fyModel(); + %h=bodeplot(mot.meas,'r',mot.tf4_2,'b',mot.tf6_4,'g'); + h=bodeplot(mot.meas,'r',mot.tf2_0,'b',mot.tf_mdl,'g',mot.w); setoptions(h,'FreqUnits','Hz','Grid','on'); end @@ -69,16 +98,14 @@ function [mot1,mot2]=identifyFxFyStage() opt=tfestOptions; opt.Display='off'; opt.initializeMethod='iv'; - opt.WeightingFilter=[1,5;20,670]*(2*pi); % Hz->rad/s conversion + opt.WeightingFilter=[1,4;10,670]*(2*pi); % Hz->rad/s conversion figure(); - mot.tf4_2 = tfest(mot.meas, 4, 2, opt); - mot.tf6_4 = tfest(mot.meas, 6, 4, opt); + mot.tf2_0 = tfest(mot.meas, 2, 0, opt);disp(str2ndOrd(mot.tf2_0)); mot.tf13_9 = tfest(mot.meas, 13, 9, opt); - num=[5.42637491e-24 1.45926254e-21 5.20861422e-17 9.92527094e-15 1.16707977e-10 1.31240975e-08 7.03191689e-05 3.08626613e-03 7.32824533e+00]; - den=[2.01035570e-35,2.33560078e-31,9.14767135e-28,2.52369732e-24,7.42150891e-21,6.89695386e-18,1.65017156e-14,5.77522779e-12,1.08386286e-08,1.13336206e-06,1.27552247e-03,2.19776479e-02,1.00000000e+00]; - mot.tf_py = idtf(num,den); - h=bodeplot(mot.meas,'r',mot.tf4_2,'b',mot.tf6_4,'g',mot.tf13_9,'m',mot.tf_py,'b'); + mot.tf_mdl = fxModel(); + %h=bodeplot(mot.meas,'r',mot.tf4_2,'b',mot.tf6_4,'g',mot.tf13_9,'m',mot.tf_py,'b'); + h=bodeplot(mot.meas,'r',mot.tf2_0,'b',mot.tf_mdl,'g',mot.w); setoptions(h,'FreqUnits','Hz','Grid','on'); %controlSystemDesigner('bode',1,mot.tf_py); % <<<<<<<<< This opens a transferfiûnction that can be edited @@ -91,10 +118,16 @@ end function f=SCRATCH() + [m1,m2]=identifyFxFyStage(); + controlSystemDesigner(1,m2.tf_py); % <<<<<<<<< This opens a transferfiûnction that can be edited %identification toolbox systemIdentification + + + + %opt=tfestOptions('Display','off'); %opt=tfestOptions('Display','on','initializeMethod','svf'); %opt=tfestOptions('Display','on','initializeMethod','iv','WeightingFilter',[]); @@ -183,7 +216,5 @@ function f=SCRATCH() controlSystemDesigner('bode',mot2); - [m1,m2]=identifyFxFyStage(); - controlSystemDesigner(1,m2.tf_py); % <<<<<<<<< This opens a transferfiûnction that can be edited end \ No newline at end of file diff --git a/python/MXTuning.py b/python/MXTuning.py index 6537354..ed7c1d5 100755 --- a/python/MXTuning.py +++ b/python/MXTuning.py @@ -61,20 +61,25 @@ class MXTuning(Tuning): def bode_model_plot(self, mot,base): self.bode_full_plot(mot,base) fig=plt.gcf() + _N=1.#E-3 # normalization factor: -> moves 3 decades to right but has factors around 1 + # s -> ms + # Hz -> kHz + # rad/s -> rad/ms if mot==1: - db_mag1=17.3 #dB - mag1=10**(db_mag1/20) - f1=6.5 #Hz - w1=f1*2*np.pi #rad/sec + #identify matlab: k:0.671226 w0:134.705 damp:0.191004 + mag1=0.671226 #10**(db_mag1/20) + db_mag1=20*np.log10(mag1)#dB + w1=134.705*_N #rad/sec + f1=w1/2/np.pi; # ca. 6.5Hz T1=1/w1 - d1=.7 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 + d1=0.19 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 num1=np.poly1d([mag1]) den1 = np.poly1d([T1**2,2*T1*d1,1]) #first resonance frequency f2=np.array([197,199]) d2=np.array([.02,.02])#daempfung - w2=f2*2*np.pi #rad/sec + w2=f2*2*np.pi*_N #rad/sec T2=1/w2 num2 = np.poly1d([T2[0]**2,2*T2[0]*d2[0],1]) den2 = np.poly1d([T2[1]**2,2*T2[1]*d2[1],1]) @@ -83,34 +88,36 @@ class MXTuning(Tuning): #current loop 2nd order approx - f4=900. - d4=1 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 - w4=f4*2*np.pi #rad/sec - T4=1/w4 - num4 = np.poly1d([1.]) - den4 = np.poly1d([T4**2,2*T4*d4,1]) - #mdl= signal.lti(num4, den4) #num denum - #bode(mdl) + #identification with matlab: k=1, w0=8725, d=0.75 + dc=.75 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 + wc=8725.*_N # rad/sec + #...but phase lag seems to have earlier effect -> reduce wc + wc*=.5 # rad/sec + fc=wc/2/np.pi # ca 1388Hz + Tc=1/wc + numc = np.poly1d([1.]) + denc = np.poly1d([Tc**2,2*Tc*dc,1]) - num=num1*num2*num4#*num3 - den=den1*den2*den4#*den3 + num=num1*num2*numc#*num3 + den=den1*den2*denc#*den3 mdl= signal.lti(num, den) #num denum - print num,den + print num + print den print mdl elif mot==2: - # basic 1/s^2 system with damping an d resonance - db_mag1=17.3 #dB - mag1=10**(db_mag1/20) - f1=4.5 #Hz - w1=f1*2*np.pi #rad/sec - d1=.3 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 + #identify matlab: k:1.7282 w0:51.069 damp:0.327613 + mag1=1.7282 #10**(db_mag1/20) + db_mag1=20*np.log10(mag1)#dB + w1=51.069*_N #rad/sec + f1=w1/2/np.pi; # ca. 6.5Hz T1=1/w1 - num1 = np.poly1d([mag1]) + d1=0.32 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 + num1=np.poly1d([mag1]) den1 = np.poly1d([T1**2,2*T1*d1,1]) - #first resonance frequency + #resonance frequency f2=np.array([57.8,61.8]) - d2=np.array([.05,.055])#daempfung + d2=np.array([.08,.095])#daempfung w2=f2*2*np.pi #rad/sec T2=1/w2 num2 = np.poly1d([T2[0]**2,2*T2[0]*d2[0],1]) @@ -118,9 +125,9 @@ class MXTuning(Tuning): mdl= signal.lti(num2, den2) #num denum #bode(mdl) - #second resonance frequency - f3=np.array([138,151]) - d3=np.array([.04,.03])#daempfung + #resonance frequency + f3=np.array([136,148]) + d3=np.array([.05,.05])#daempfung w3=f3*2*np.pi #rad/sec T3=1/w3 num3 = np.poly1d([T3[0]**2,2*T3[0]*d3[0],1]) @@ -128,9 +135,9 @@ class MXTuning(Tuning): #mdl= signal.lti(num3, den3) #num denum #bode(mdl) - #second resonance frequency + #resonance frequency f4=np.array([410,417]) - d4=np.array([.015,.02])#daempfung + d4=np.array([.015,.015])#daempfung w4=f4*2*np.pi #rad/sec T4=1/w4 num4 = np.poly1d([T4[0]**2,2*T4[0]*d4[0],1]) @@ -138,8 +145,9 @@ class MXTuning(Tuning): #mdl= signal.lti(num3, den3) #num denum #bode(mdl) - f5=np.array([228,230]) - d5=np.array([.03,.03])#daempfung + #resonance frequency + f5=np.array([230,233]) + d5=np.array([.04,.04])#daempfung w5=f5*2*np.pi #rad/sec T5=1/w5 num5 = np.poly1d([T5[0]**2,2*T5[0]*d5[0],1]) @@ -147,18 +155,22 @@ class MXTuning(Tuning): #current loop 2nd order approx - fc=900. - dc=1 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 - wc=fc*2*np.pi #rad/sec + #identification with matlab: k=1, w0=8725, d=0.75 + dc=.75 # daempfung =1 -> keine resonanz -> den1= np.poly1d([T1,1])**2 + wc=8725.*_N # rad/sec + #...but phase lag seems to have earlier effect -> reduce wc + wc*=.5 # rad/sec + fc=wc/2/np.pi # ca 1388Hz Tc=1/wc numc = np.poly1d([1.]) denc = np.poly1d([Tc**2,2*Tc*dc,1]) - #mdl= signal.lti(num4, den4) #num denum - #bode(mdl) num=num1*num2*num3*num4*num5*numc den=den1*den2*den3*den4*den5*denc mdl= signal.lti(num, den) #num denum + print num + print den + print mdl bode(mdl) w=np.logspace(0,np.log10(2000),1000)*2*np.pi w,mag,phase = signal.bode(mdl,w)