%% ABOUT THIS FILE % ---------------------------------------------------------------------- % Type-A with Dummy Mass for KAGRA % Coded by K. Okutomi on 2016.06.28 % ---------------------------------------------------------------------- %% PRELIMINARY clear all; % Clear workspace close all; % Close plot windows %% cd('/Users/kokioktm/Documents/suspension/MATLAB/script/TypeA'); % Change directory addpath('../../utility'); % Add path to utilities addpath('../../oklab'); addpath('controlmodel'); addpath('susmodel'); g = 9.81; pp = 2.0*pi; %% IMPORT SUSPENSION MODEL matfile='TypeA161114mdl'; load([matfile,'.mat']); %% TUNING DAMPER % This part compensates the failure in converting the structural damping to % viscous damping. %% %%%%% NO-CONTROL MODEL %%%%% %% IMPORT SERVO FILTERS addpath('parameter'); % Add path to servo TypeA_paramNoCtrl170105; % NO CONTROL MODE rmpath ('parameter'); % Remove path to servo %% IMPORT SIMULINK MODEL mdlfile='TypeA_ctrl161209'; st =linmod(mdlfile); invl =strrep(st.InputName, [mdlfile,'/'],''); outvl =strrep(st.OutputName,[mdlfile,'/'],''); sysc0 =ss(st.a,st.b,st.c,st.d,'inputname',invl,'outputname',outvl); %% FREQUENCY freq =logspace(-2,2,1001); freq2=logspace(-3,1,1001); %% PLOT TRANSFER FUNCTION OF THE SUSPENSION MECHANICAL RESPONSE %{ TF_F0 = calcbodesus(sysc0,... {'injLF0','LVDT_LF0';'injTF0','LVDT_TF0';'injYF0','LVDT_YF0'}, freq); TF_BF = calcbodesus(sysc0,... {'injLBF','LVDT_LBF';'injTBF','LVDT_TBF';'injVBF','LVDT_VBF';... 'injRBF','LVDT_RBF';'injPBF','LVDT_PBF';'injYBF','LVDT_YBF'}, freq); TFtest = calcbodesus(sysc0,... {'injYBF','LVDT_YF0'}, freq2); TF_GAS = calcbodesus(sysc0,... {'injGASF0','LVDT_GASF0';'injGASF1','LVDT_GASF1';'injGASF2','LVDT_GASF2';... 'injGASF3','LVDT_GASF3';'injGASBF','LVDT_GASBF'},freq2); %} %% LF0 actuation close all bodesusplotcmpopt3(sysc0,... {'injLF0','LVDT_LF0';'injLF0','LVDT_LBF';'injLF0','dispLDP'},... freq,... 'ylim',[1e-6,1e2],... 'unit1','N','unit2','m'); saveas(gca, 'fig/tfTypeA_DP_actLF0_160728.png') %% YF0 actuation bodesusplotcmpopt3(sysc0,... {'injYF0','LVDT_YF0';'injYF0','LVDT_YBF';'injYF0','dispYDP'},... freq2,... 'ylim',[1e-4,1e2],... 'unit1','(NEm)','unit2','rad'); saveas(gca, 'fig/tfTypeA_DP_actYF0_160728.png') %% LBF actuation close all bodesusplotcmpopt3(sysc0,... {'injLBF','LVDT_LF0';'injLBF','LVDT_LBF';'injLBF','dispLDP'},... freq,... 'ylim',[1e-6,1e2],... 'unit1','N','unit2','m'); saveas(gca, 'fig/tfTypeA_DP_actLBF_160728.png') %% YBF actuation close all bodesusplotcmpopt3(sysc0,... {'injYBF','LVDT_YF0';'injYBF','LVDT_YBF';'injYBF','dispYDP'},... freq2,... 'ylim',[1e-6,1e2],... 'unit1','(NEm)','unit2','rad'); % saveas(gca, 'fig/tfTypeA_DP_actYBF_160728.png') %% GASF0 actuation close all bodesusplotcmpopt3(sysc0,... {'injGASF0','LVDT_GASF0';'injGASF0','LVDT_GASBF';'injGASF0','dispVDP'},... freq,... 'ylim',[1e-3,1e3],... 'unit1','N','unit2','m'); saveas(gca, 'fig/tfTypeA_DP_actGASF0_160728.png') %} %% PLOT OPEN-LOOP TRANSFER FUNCTION OF FEEDBACK CONTROLS %% LF0 control Loop ipservo = myzpk([0], [3., 3.], 300*pp); OLTF_LF0 = bodesusplotOLTF(sysc0,... 'injLF0','LVDT_LF0',ipservo,freq,... 'ylim',[1e-4,1e2],... 'unit1','1'); saveas(gca, 'fig/OLtf_TypeA_DP_LF0_160728.png') %% YF0 control loop ipservoY = myzpk([0], [1.e0, 1.e0], 10*pp); OLTF_YF0 = bodesusplotOLTF(sysc0,... 'injYF0','LVDT_YF0',ipservoY,freq2,... 'ylim',[1e-2,1e3],... 'unit1','1'); saveas(gca, 'fig/OLtf_TypeA_DP_YF0_160728.png') %% YBF control loop close all bfservoY = myzpk([0], [0.3, 0.3], 5*pp); OLTF_YF0 = bodesusplotOLTF(sysc0,... 'injYBF','LVDT_YBF',bfservoY,freq,... 'ylim',[1e-4,1e2],... 'unit1','1'); saveas(gca, 'fig/OLtf_TypeA_DP_YBF_160824.png') %% GASF0 control loop close all gasf0servo = myzpk([0], [5., 5.], 100*pp); OLTF_YF0 = bodesusplotOLTF(sysc0,... 'injGASF0','LVDT_GASF0',gasf0servo,freq,... 'ylim',[1e-4,1e2],... 'unit1','1'); saveas(gca, 'fig/OLtf_TypeA_DP_GASF0_160728.png') %% %%%%% CONTROL MODEL %%%%% %% IMPORT SERVO FILTERS addpath('parameter'); % Add path to servo TypeA_DP_paramAllCtrl160728; % ALL CONTROL MODE rmpath ('parameter'); % Remove path to servo %% IMPORT SIMULINK MODEL mdlfile='TypeA_DP_ctrl160728'; % typeA_DP ver.160713 st =linmod(mdlfile); invl =strrep(st.InputName, [mdlfile,'/'],''); outvl =strrep(st.OutputName,[mdlfile,'/'],''); sysc =ss(st.a,st.b,st.c,st.d,'inputname',invl,'outputname',outvl); %% PLOT NYQUIST DIAGRAM % set(0,'defaultFigurePosition',[50, 50, 850, 650]); % set(0,'defaultAxesLineWidth',1.5); % set(0,'defaultLineLineWidth',1.5); nyq = nyquistoptions; nyq.Grid = 'on'; nyq.ShowFullContour = 'on'; fig = figure; % nyqLF0 = nyquistplot(sysc('LVDT_LF0','injLF0'),nyq); % nyqYF0 = nyquistplot(sysc('LVDT_YF0','injYF0'),nyq); % nyqYBF = nyquistplot(sysc('LVDT_YBF','injYBF'),nyq); nyqGASF0 = nyquistplot(sysc('LVDT_GASF0','injGASF0'),nyq); a=get(gca,'Children'); data=get(a(1),'Children'); set(data(1), 'LineWidth',3); set(data(2), 'LineWidth',3); %% PLOT CLOSED LOOP TRANSFER FUNCTION %% LF0 actuation bodesusplotcmpopt3(sysc,... {'injLF0','LVDT_LF0';'injLF0','LVDT_LBF';'injLF0','dispLDP'},... freq,... 'ylim',[1e-6,1e2],... 'unit1','N','unit2','m'); saveas(gca, 'fig/CLTF_TypeA_DP_actLF0_160728.png') %% YF0 actuation bodesusplotcmpopt3(sysc,... {'injYF0','LVDT_YF0';'injYF0','LVDT_YBF';'injYF0','dispYDP'},... freq2,... 'ylim',[1e-4,1e2],... 'unit1','(NEm)','unit2','rad'); saveas(gca, 'fig/CLTF_TypeA_DP_actYF0_160728.png') %% YBF actuation close all bodesusplotcmpopt3(sysc,... {'injYBF','LVDT_YF0';'injYBF','LVDT_YBF';'injYBF','dispYDP'},... freq2,... 'ylim',[1e-6,1e2],... 'unit1','(NEm)','unit2','rad'); saveas(gca, 'fig/CFTF_TypeA_DP_actYBF_160728.png') %% GASF0 actuation close all bodesusplotcmpopt3(sysc,... {'injGASF0','LVDT_GASF0';'injGASF0','LVDT_GASBF';'injGASF0','dispVDP'},... freq,... 'ylim',[1e-3,1e3],... 'unit1','N','unit2','m'); saveas(gca, 'fig/CLTF_TypeA_DP_actGASF0_160728.png') %% PLOT VIBRATION ISOLARION RATIO bodesusplotcmpopt4({sysc,sysc0},... 'accLGnd','dispLDP',... freq,'legend',{'Control ON','Control OFF'},... 'ylim',[1e-3,1e3],... 'unit1','m','unit2','m'); saveas(gca, 'fig/VIR_TypeA_DP_toLDP_160728.png') %% bodesusplotcmpopt3(sysc0,... {'accLGnd','LVDT_LF0';'accLGnd','LVDT_LBF';'accLGnd','dispLDP'},... freq,... 'title','Displacement from Ground Acceleration, Control OFF',... 'ylim',[1e-8,1e4],... 'unit1','(m/s^2)','unit2','m'); bodesusplotcmpopt3(sysc,... {'accLGnd','LVDT_LF0';'accLGnd','LVDT_LBF';'accLGnd','dispLDP'},... freq,... 'title','Displacement from Ground Acceleration, Control ON',... 'ylim',[1e-8,1e4],... 'unit1','(m/s^2)','unit2','m'); % saveas(gca, 'fig/VIR_TypeA_DP_NoCtrl_toLDP_160728.png') %% NOISE % NOISE data_LVDT = importdata('../../noise/LVDTnoiseADC_disp.dat'); % um/rtHz noise_LVDT = interp1(data_LVDT(:,1),data_LVDT(:,2)*1e-6,freq')'; % m/rtHz noiseLVDT_LF0 = noise_LVDT/sqrt(3/2); % [m/rtHz] data_yawBF = importdata('../../noise/psd_yawBF.txt'); noiseYawBF = interp1(data_yawBF(:,1),data_yawBF(:,2),freq')'; % data_GEO = importdata('../../noise/GEOnoiseproto_vel.dat'); % um/sec/rtHz % noise_GEO = interp1(data_GEO(:,1),data_GEO(:,2)*1e-6,freq')';% m/sec/rtHz % n_GEO_LF0 = noise_GEO/sqrt(3/2); % [m/rtHz] data_seismic = importdata('../../noise/KamiokaSeismicHighNoise.dat'); % m/rtHz noiseSeismic = interp1(data_seismic(:,1),data_seismic(:,2),freq')';% m/rtHz % % %% MODEL % addpath('servofilter'); % Add path to servo % typeA_RT_servo_lock_160427; % LOCK ACQUISITION MODE % rmpath ('servofilter'); % Remove path to servo % st =linmod(mdlfile); % invl =strrep(st.InputName, [mdlfile,'/'],''); % outvl =strrep(st.OutputName,[mdlfile,'/'],''); % sysc =ss(st.a,st.b,st.c,st.d,'inputname',invl,'outputname',outvl); %% Noise propagation close all bodesusplotcmpopt3(sysc,... {'noiseLVDT_LBF','LVDT_LBF';'noiseLVDT_YBF','LVDT_YBF'},... freq2,... 'ylim',[1e-3,1e3]);%,... % 'unit1','N','unit2','m'); % saveas(gca, 'fig/CLTF_TypeA_DP_actGASF0_160728.png') %% bodesusplotcmpopt4({sysc,sysc0},... 'noiseLVDT_LBF','LVDT_LBF',... freq2,'legend',{'Control ON','Control OFF'},... 'ylim',[1e-3,1e3],... 'unit1','m','unit2','m'); bodesusplotcmpopt4({sysc,sysc0},... 'noiseLVDT_YBF','LVDT_YBF',... freq2,'legend',{'Control ON','Control OFF'},... 'ylim',[1e-3,1e3],... 'unit1','rad','unit2','rad'); %% CALC [TFnLVDT_BFL,~]=bodesus(sysc0,'noiseLVDT_LF0','LVDT_LBF',freq); [TFnLVDT_BFY,~]=bodesus(sysc0,'noiseLVDT_LF0','LVDT_LBF',freq); % [magsG,~]=bodesus(sysc,'n_ACC_LF0','disp_LDM',freq); [magsS,~]=bodesus(sysc,'accLGnd','dispLDP',freq); [mag0S,~]=bodesus(sysc0,'accLGnd','dispLDP',freq); % magsL=magsL.*n_LVDT_LF0; % magsG=magsG.*n_GEO_LF0; magsS=magsS.*freq.*freq*pp*pp.*noiseSeismic; mag0S=mag0S.*freq.*freq*pp*pp.*noiseSeismic; % magsT=sumpsd({magsL,magsG,magsS}); % rmsS=makerms(freq,magsS); % rmsSv=makerms(freq,magsS.*freq*pp); % rms0S=makerms(freq,mag0S); % rms0Sv=makerms(freq,mag0S.*freq*pp); % rmsSy=makerms(freq,magsSy); % rmsSyv=makerms(freq,magsSy.*freq*pp); % rms0Sy=makerms(freq,mag0Sy); % rms0Syv=makerms(freq,mag0Sy.*freq*pp); % rmsT=makerms(freq,magsT); % rmsTv=makerms(freq,magsT.*freq*pp); %% K. Okutomi editted below [magcSyawBF,~]=bodesus(sysc,'noiseLVDT_YBF','LVDT_YBF',freq); [mag0SyawBF,~]=bodesus(sysc0,'noiseLVDT_YBF','LVDT_YBF',freq); [magcSyawDP,~]=bodesus(sysc,'noiseLVDT_YBF','dispYDP',freq); [mag0SyawDP,~]=bodesus(sysc0,'noiseLVDT_YBF','dispYDP',freq); %% specYawBF=magcSyawBF.*freq.*freq*pp*pp.*noiseYawBF; spe0YawBF=mag0SyawBF.*freq.*freq*pp*pp.*noiseYawBF; specYawDP=magcSyawDP.*freq.*freq*pp*pp.*noiseYawBF; spe0YawDP=mag0SyawDP.*freq.*freq*pp*pp.*noiseYawBF; %% rmscYawBF=makerms(freq,specYawBF); rms0YawBF=makerms(freq,spe0YawBF); rmscYawDP=makerms(freq,specYawDP); rms0YawDP=makerms(freq,spe0YawDP); %% YBF control loop close all bfservoY = myzpk([0], [0.3, 0.3], 5*pp); OLTF_YF0 = bodesusplotOLTF(sysc0,... 'injYBF','LVDT_YBF',bfservoY,freq,... 'ylim',[1e-4,1e2],... 'unit1','1'); % TRANSFER FUNCTION [magYBF,phsYBF] = mybode(sysc0('LVDT_YBF','injYBF')*bfservoY,freq); TFsusYBF = magYBF.*exp(1j*phsYBF*pi/180); [magFltYBF,phsFltYBF] = mybode(bfservoY,freq); TFfltYBF = magFltYBF.*exp(1j*phsFltYBF*pi/180); OLTFybf = -1*TFsusYBF.*TFfltYBF; magOLTFybf = abs(OLTFybf); phsOLTFybf = angle(OLTFybf)*180/pi; % feedback = OLTFybf./TFYBF./(1-OLTFybf).* %% Requirement % reqBS=freq; % for i=1:length(freq) % if freq(i)<8 % reqBS(i)=0; % else % reqBS(i)=1e-15*(freq(i)/10)^(-1.5); % end % end %% PLOT mypsdplotopt({noiseYawBF,spe0YawBF,specYawBF,rms0YawBF,rmscYawBF},freq,... 'legend',{'Seismic','Control OFF','Control ON',... 'Control OFF RMS','Control ON RMS'},... 'color',{'k-','b-','r-','c--','m--'},... 'ylim',[1e-18,1e-6],'ylabel','Magnitude [rad/rtHz] or [rad]') % saveas(gca,'fig/psdTypeA_DP_AllCtrl_LDP160728.png') %% mypsdplotopt({noiseYawBF,spe0YawDP,specYawDP,rms0YawDP,rmscYawDP},freq,... 'legend',{'Seismic','Control OFF','Control ON',... 'Control OFF RMS','Control ON RMS'},... 'color',{'k-','b-','r-','c--','m--'},... 'ylim',[1e-18,1e-6],'ylabel','Magnitude [rad/rtHz] or [rad]') % saveas(gca,'fig/psdTypeA_DP_AllCtrl_1mmAssymYDP160728.png') %% % syspole = poledamping({sysc0, sysc}); % saveas(gca, 'fig/decaytimeTypeA_DP_AllCtrl160728.png') % mypsdplotopt({magsG,magsL,magsS,magsT,rmsT,rmsTv},freq,... % 'title','Blending at 0.05 Hz',... % 'legend',{'GEO','LVDT','Seismic','total','RMS',... % 'RMS velocity'},... % 'color',{'b-','g-','k-','r-','r--','m--'},... % 'ylim',[1e-12,1e-5],'ylabel','Magnitude [m/rtHz] or [m] or [m/sec]') % saveas(gca,'figure/typeA_RT_IPdamp_DCservoON_160427','epsc') % % %% COMPARE % mypsdplotopt({magsG,magsL,magsS,magsT,mag0S,reqBS},freq,... % 'title','Noise coupling',... % 'legend',{'geophone','LVDT','seismic','Controlled total'... % 'No controlled','BS requirement'},... % 'color',{'g-','b-','k-','r-','m--','k--'},... % 'ylim',[1e-20,1e-5],'ylabel','Magnitude [m/rtHz]') % saveas(gca,'figure/typeA_RT_IPdamp_highfreq_noise_160427','epsc') %