clear all; close all; % format of the data file: % 1st column ... frequency [Hz] % 2nd column ... gain from T_M to omega_M % 3rd column ... phase [deg] from T_M to omega_M % 4th column ... gain from T_M to omega_L % 5th column ... phase [deg] from T_M to omega_L % dimension of plant % input ... torque T_M [Nm] % output ... speed omega_M [rpm] load frdata_amp2_offset10_fixed.dat; load frdata_amp3_offset10_fixed.dat; load frdata_amp4_offset10_fixed.dat; load frdata_amp3_offset5_fixed.dat; load frdata_amp3_offset0_fixed.dat; load frdata_amp10_offset20_fixed.dat; P1_jw = frdata_amp2_offset10_fixed(:,2).*exp(i*pi/180.0*frdata_amp2_offset10_fixed(:,3)); P2_jw = frdata_amp3_offset10_fixed(:,2).*exp(i*pi/180.0*frdata_amp3_offset10_fixed(:,3)); P3_jw = frdata_amp4_offset10_fixed(:,2).*exp(i*pi/180.0*frdata_amp4_offset10_fixed(:,3)); P4_jw = frdata_amp3_offset5_fixed(:,2) .*exp(i*pi/180.0*frdata_amp3_offset5_fixed(:,3)); P5_jw = frdata_amp3_offset0_fixed(:,2) .*exp(i*pi/180.0*frdata_amp3_offset0_fixed(:,3)); P6_jw = frdata_amp10_offset20_fixed(:,2) .*exp(i*pi/180.0*frdata_amp10_offset20_fixed(:,3)); f = frdata_amp2_offset10_fixed(:,1); % frequency grid is commonly used for frequency response exp. w = 2*pi*f; P1_g = frd(P1_jw, w); P2_g = frd(P2_jw, w); P3_g = frd(P3_jw, w); P4_g = frd(P4_jw, w); P5_g = frd(P5_jw, w); P6_g = frd(P6_jw, w); % phisycal model J_M = 1.73E-05; J_L = 2.04E-04; K_S = 523; %[Nm/rad] C_S = 0.01; C_L = 1e-5; Pp_a = [0, 1, -1; -K_S/J_M, -C_S/J_M, C_S/J_M; K_S/J_L, C_S/J_L, -(C_S+C_L)/J_L]; Pp_b = [0, ; 1/J_M; 0]; Pp_c = [0, 1, 0]; Pp_d = 0; Pp = ss(Pp_a, Pp_b, Pp_c, Pp_d); Pp_g = Pp; figure(1); bode(P1_g, 'b.-', P2_g, 'k.-', P3_g, 'g.-', P4_g, 'm.-', P5_g, 'y.-', P6_g, 'c', Pp_g, 'r--'); legend('amp=2, offset=10', 'amp=3, offset=10', 'amp=4, offset=10', 'amp=3, offset=5', 'amp=3, offset=0', 'amp=10, offset=20', 'physical model', 'Location', 'SouthWest'); grid on; axis([1, 2000, -1500, 180]); % phase h = get(1); axis(h.Children(4),[1, 2000, -30, 70]); % gain Pave_g = (P1_g + P2_g + P3_g + P4_g + P6_g)/5; P0 = fitfrd(Pave_g, 8, 1); P0 = balreal(P0); P0_g = frd(P0, w); figure(2); bode(Pave_g, 'b.-', Pp_g, 'k', P0_g, 'r') legend('averaged freq. resp.', 'physical model', 'nominal plant', 'Location', 'SouthWest'); axis([1, 2000, -1500, 180]); % phase h = get(2); axis(h.Children(4),[1, 2000, -30, 70]); % gain print -depsc nominal.eps; % save as EPS format print -djpeg -r0 nominal.jpg; % save as JPEG format print -dpng -r0 nominal.png; % save as PNG format save nominal.mat P0 P1_g P2_g P3_g P4_g P5_g P6_g Pp_g P0_g w f