% nominal.m % 1. to compare frequency response experiment result and physical model % 2. to find nominal plant P0(s) by curve fitting clear all; close all; load freqresp.mat; % inertia moment of driving side J_M = 1.11e-5 + 64.6e-6; % inertia moment of driven side J_L0 = 6.04e-6 + 64.6e-6 + 14.0e-6; J_L1 = J_L0 + 1.41e-4; % with small disc J_L2 = J_L0 + 5.38e-4; % with large disc % torsional spring constant of belt K_S = 250; %[Nm/rad] C_S = 0.02; C_L = 1e-5; % state space model of two inertia system Pp_a = [0, 1, -1; -K_S/J_M, -C_S/J_M, C_S/J_M; K_S/J_L1, C_S/J_L1, -(C_S+C_L)/J_L1]; Pp_b = [0, ; 1/J_M; 0]; Pp_c = [0, 1, 0]; Pp_d = 0; Pp1 = ss(Pp_a, Pp_b, Pp_c, Pp_d); % with small disc Pp_a = [0, 1, -1; -K_S/J_M, -C_S/J_M, C_S/J_M; K_S/J_L2, C_S/J_L2, -(C_S+C_L)/J_L2]; Pp2 = ss(Pp_a, Pp_b, Pp_c, Pp_d); % with large disc % bode plot of plant (with small disc) figure(1); bode(P1_g, 'k.-', P3_g, 'm.-', P5_g, 'c', Pp1, 'r'); legend('offset=0', 'offset=5', 'offset=10', 'model'); grid on; % bode plot of plant (with large disc) figure(2); bode(P2_g, 'k.-', P4_g, 'm.-', P6_g, 'c', Pp2, 'r'); legend('offset=0', 'offset=5', 'offset=10', 'model'); grid on; %axis([1, 2000, -1500, 180]); % phase h = get(1); %axis(h.Children(4),[1, 2000, -30, 70]); % gain Pave_g = (P3_g + P4_g + P5_g + P6_g)/4; % P1_g and P2_g are neglected P0 = fitfrd(Pave_g, 5, 1); % curve fitting P0 = balreal(P0); P0_g = frd(P0, w); figure(3); bode(Pave_g, 'b.-', P0, 'r'); legend('averaged freq. resp.', 'nominal plant', 'Location', 'SouthWest'); %axis([1, 2000, -1500, 180]); % phase h = get(3); %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 P0_g w f