기계항공시스템해석 MATLAB 실습 - 시스템구성및시간 / 주파수응답그래프 - 박사과정서종상 azuresky@snuackr Tel:0-880-194 301-113 Vehicle Dynamics & Control Laboratory
Set Transfer Function and State Space 1) tf Transfer Function From = tf(numerator, Denominator) Numerator ( 분자 ), Denominator( 분모 ) 를통하여 Transfer Function 을 Define 한다 ) ss State Space From = ss(a, B, C, D) State Space A,B,C,D 로부터 State Space Form 을유도 Ex) Transfer Function 을 Define TF = 1 ss ( + 05) system_num=[1]; system_den=[1 05 0]; system=tf(system_num,system_den);
Transfer Function State Space 3) tfss [A,B,C,D] = tfss( 분자, 분모 ) Numerator ( 분자 ), Denominator( 분모 ) 를통하여 State Space Matrix A,B,C,D 유도 4) sstf (Numerator, Denominator)= sstf(a, B, C, D) Numerator ( 분자 ), Denominator( 분모 ) 를통하여 Transfer Function (Num,Den) 유도
Transfer Function 활용법 System1=tf(num1,den1); System=tf(num,den); System 간곱셈, 덧셈등의계산이가능 system_num1=[1]; system_den1=[1 05 0]; system1=tf(system_num1,system_den1); system_num= [1 35]; system_den= [1 58]; system=tf(system_num,system_den); OLS=system1*system;
선형시불변모델의시간응답 함수 내용 step 계단입력에대한선형시간불변모델의응답을구한다 Impulse 충격량입력에대한선형시간불변모델의응답을구한다 Initial 상태방정식으로표현된모델의초기조건에의한과도응답을구한다 lsim 임의의입력에대한선형시간불변모델의응답을구한다 time response란선형시간불변모델에시간 t의함수로표현되는입력이가해질경우에응답을시간 t에대해구하는것을말한다 보통 time response를구하는목적은 delay time, rise time, peak time, max overshoot, settling time 등과같은과도상태에서제어시스템의성능을나타내는값들을계산하기위함이다
Step Response step : 단위계단입력에대한응답을얻고자할경우사용된다 step(num,den) % num, den으로정의된 transfer function에대한 step response를그린다 step(num,den,t) % num, den으로정의된 transfer function에대한 step response를시간 t에대해그린다 step(a,b,c,d) % 상태공간으로정의된시스템에대한 step response를그린다 step(a,b,c,d,ui) % 상태공간으로정의된시스템의 input ui에대한 step response를그린다 step(a,b,c,d,ui,t) % 상태공간으로정의된시스템에 input ui에대한 step response를시간 t에대해그린다
Step Response Ex 1 1) m 파일 editor window ( 명령어입력 ) clear; clc; % t=0:001:10; num=[0 0 5]; den=[1 4 5]; step(num,den); % 단위계단응답 % step(num,den,t); 시간 t 에대한단위계단응답 grid on; title('unit step response of G(s)=5/(s^+4s+5)')
Step Response Ex 1 ) 결과 14 unit step response of G(s)=5/(s +4s+5) 1 1 Amplitude 08 06 04 0 0 0 05 1 15 5 3 Time (sec)
Step Response Ex 1) m 파일 editor window ( 명령어입력 ) clear; clc; A=[-1-1; 65 0]; B=[1 1; 1 0]; C=[1 0; 0 1]; D=[0 0; 0 0]; figure(1) step(a,b,c,d,1) % input u1에대한step response grid on; title('step response plot : input =u1(u=0)'); figure() step(a,b,c,d,) % input u에대한 step response grid on; title('step response plot : input =u(u1=0)'); figure(3) step(a,b,c,d) grid on; % input u1,u 에대한 step response
Step Response Ex ) 결과 04 step response plot : input =u1(u=0) 03 step response plot : input =u(u1=0) 0 To: Out(1) 0 0 To: Out(1) 01 0-0 -01 Amplitude -04 Amplitude -0 15 15 To: Out() 1 05 To: Out() 1 05 0 0 4 6 8 10 1 Time (sec) Step Response From: In(1) 04 From: In() 0 0 4 6 8 10 1 Time (sec) 0 To: Out(1) 0-0 Amplitude -04 15 To: Out() 1 05 0 0 4 6 8 10 1 0 4 6 8 10 1 Time (sec)
Impulse Response impulse : 단위충격량입력에대한응답을얻고자할경우사용된다 impulse(num,den) % num, den으로정의된 transfer function에대한 impulse response를그린다 impulse(num,den,t) % num, den으로정의된 transfer function에대한 impulse response를시간 t에대해그린다 impulse( A,B,C,D) % 상태공간으로정의된시스템에대한 impulse response를그린다 impulse(a,b,c,d,iu) % 상태공간으로정의된시스템의 input iu에대한 impulse response를그린다 impulse(a,b,c,d,iu,t) % 상태공간으로정의된시스템에 input iu에대한 impulse response를시간 t에대해그린다
Initial Response initial : 상태방정식으로표현된모델의초기조건에의한과도응답을얻고자할경우 사용된다 step 이나 impulse는외부입력에대한응답을구한다 이때초기조건은모두 0으로되는데, 초기조건에대한응답을구하고싶을때는 initial 을사용한다 즉, initial은외부입력이없다고하고단지초기조건에대한응답만을구한다 이때에는모델이상태방정식으로표현된모델이어야만하는데그이유는전달함수의경우에는이미전달함수를유도하기위해초기조건을모두 0으로한다는가정을사용했기때문이다
Initial Response Ex 1 Initial( A,B,C,D,x0) % 상태공간으로정의된시스템의초기값 x0 에대한 response 를그린다 Initial(A,B,C,D,x0, t) % 상태공간으로정의된시스템에 input iu 에대한 impulse response 를시간 t 에대해그린다 1) m 파일 editor window ( 명령어입력 ) clear; clc; A=[-1-1; 65 0]; B=[1 1; 1 0]; C=[1 0; 0 1]; D=[0 0; 0 0]; x0=[1;0]; % x1 의초기값 1, x 의초기값 0 figure(1) initial(a,b,c,d,x0); grid on; title('initial response plot');
Initial Response Ex 1 ) 결과 1 initial response plot 05 To: Out(1) 0-05 Amplitude -1 1 To: Out() 0-1 - 0 4 6 8 10 1 Time (sec)
lsim Response lsim : 임의의입력에대한선형시간불변모델의응답을구한다 step 이나 impulse 는미리정해진외부입력에대한응답을구한다 하지만 lsim 은임의의입력에대한응답을구할수있다 lsim(num,den,u,t) % transfer function으로정의된시스템의 input u에대한 response를시간 t에대해그린다 lsim(num,den,u,t,x0) % transfer function으로정의된시스템의초기값 x0, input u에대한 response를시간 t에대해그린다 lsim(a,b,c,d,u,t) % 상태공간으로정의된시스템의 input u에대한 response를시간 t에대해그린다 lsim(a,b,c,d,u,t,x0) % 상태공간으로정의된시스템의초기값 x0, input u에대한 response를시간 t에대해그린다
lsim Response Ex 1 1) m 파일 editor window ( 명령어입력 ) clear; clc; t=0:001:10; A=[-1-1; 65 0]; B=[1 1; 1 0]; C=[1 0; 0 1]; D=[0 0; 0 0]; u=[0*t+; 0*t+5]; % input u1=*1(t), u=5*1(t) x0=[5;0]; % 초기값 x1=5, x=0 figure(1) lsim (A,B,C,D,u,t) % 시간 t에동안상태방정식으로정의된시스템의 input u에대한 response grid on; title('linear simulation response plot'); figure() lsim (A,B,C,D,u,t,x0) % 시간 t에동안상태방정식으로정의된시스템의 input u와초기값 x0에대한 response grid on; grid on; title('linear simulation response plot');
lsim Response Ex 1 ) 결과 6 linear simulation response plot 6 linear simulation response plot 4 4 To: Out(1) To: Out(1) 0 0 - Amplitude To: Out() - 1 10 8 6 4 Amplitude To: Out() -4 0 15 10 5 0 0 1 3 4 5 6 7 8 9 10 0 0 1 3 4 5 6 7 8 9 10 Time (sec) Time (sec)
Plotting bode diagrams with MATLAB 1) From Transfer Function [Magnitude, Phase, Frequency values]=bode(numerator, Denominator, W); Numerator ( 분자 ), Denominator( 분모 ) 를통하여 Bode 의 Magnitude 와 Phase 를구한다 ) From State Space [Magnitude, Phase, Frequency values]=bode(a, B, C, D, W); State Space 의각행렬 A, B, C, D 를통하여 Bode 의 Magnitude 와 Phase 를구한다
Bode Plot Example TF = 1 ss ( + 05) %% System Define system_num=[1]; system_den=[1 05 0]; %% Generates logarithmic points w=logspace(0,3,100); %% Bode [mag,phase,w]=bode(system_num,system_den); % bode(system_num,system_den); %% Magnitude magdb = 0*log10(mag); %% Plot figure(1) subplot(11) semilogx(w/(*pi),magdb,'b','linewidth',3);hold on; grid on; xlabel('frequency [rad/s]') ylabel('magnitude [db]') subplot(1) semilogx(w/(*pi),phase,'b','linewidth',3);hold on; grid on; xlabel('frequency [rad/s]') ylabel('phase [deg]') Magnitude [db] Phase [deg] 50 0-50 10-3 10-10 -1 10 0 10 1 Frequency [rad/s] -50-100 -150-00 10-3 10-10 -1 10 0 10 1 Frequency [rad/s]
Suspension Example m y Design Considerations 1 Ride Quality Sprung mass acceleration : y k b Rattle space Suspension Deflection : y x m1 x 3 Tire Force Vibration Tire Deflection : x u k 1 P u Spring Stiffness : k Damping Ratio : b Tire Stiffness : k 1 Suspension Design Parameters
Dynamic Equations Free Body Diagram k ( y x) by ( x ) y m m1 x k ( y x) by ( x ) k ( u x) 1 Dynamic Equations m x= k ( y x) + b( y x ) + k ( u x) 1 1 m y = k ( y x) b( y x )
Laplace Transform Laplace Transform [ m s + bs + ( k + k )] X () s = ( bs + k ) Y () s + k U () s 1 1 1 [ m s + bs + k ] Y () s = ( bs + k ) X () s Displacement of Mass Y() s k ( bs + k ) = U ( s) m m s ( m m ) bs [( k m ( k k ) m ] s k bs k k 1 4 3 1 + 1+ + 1+ 1+ + 1 + 1 X() s k ( m s + bs + k ) = U ( s) m m s ( m m ) bs [( k m ( k k ) m ] s k bs k k 1 4 3 1 + 1+ + 1+ 1+ + 1 + 1 Design Considerations s Y() s s k ( bs + k ) 1 1() = = 4 3 U( s) mm 1 s + ( m1+ m) bs + [( km1+ ( k1+ k) m] s + kb 1s+ k1k G s Sprung mass acceleration : y G Y() s X() s kms 1 () s = = 4 3 U ( s) m1m s + ( m1+ m) bs + [( km1+ ( k1+ k) m] s + k1bs + k1k Suspension Deflection : y x X() s U() s mm s ( m + m ) bs k ( m + m ) s 4 3 1 1 1 3() = = 4 3 U ( s) m1m s + ( m1+ m) bs + [( km1+ ( k1+ k) m] s + k1bs + k1k G s Tire Deflection : x u
State Equation General Form of State Equation x = Ax + Bu y = Cx + Du Dynamic Equations m z = k ( z z ) + b( z z ) + k ( u z ) 1 u s u s u 1 u The State variables ( x= z, y = z ) u s m ( ) ( ) zs = k zs zu b z s z u x = z z : Suspension Deflection 1 x = z : absolute velocity of sprung mass x = z u : Tire Deflection 3 x = z : absolute velocity of unsprung mass 4 s s u u u 1 st order State equations x = z z = x x 1 3 s u 4 k b k b b x = ( z z ) ( z z ) = x x + x m m m m m u s u s u 1 x = z u = x u 4 k b k k b k b x 4 = ( z z ) + ( z z ) ( z u) = x + x x m m m m m m m 1 1 s u s u u 1 3 4 1 1 1 1 1 1 1 4 x
System Matrix Matrix Form of State equations (system matrix) 0 1 0 1 x 1 k b b x1 0 0 x m m m x 0 = + u x 3 0 0 0 1 x 3 1 x 4 k b k1 b x4 0 m1 m1 m1 m1 Matrix Form of State equations (output matrix) k b b y = x = x x + x : Sprung mass acceleration 1 1 4 m m m y = z z = x : Suspension Deflection s u 1 y = z u = x : Tire Deflection 3 u 3 y k b b 0 m m m x 1 1 x y 1 0 0 0 = x 3 y 3 x4 0 0 1 0
MATLAB Simulation using Laplace Transform Suspension Parameters m1=55; m=400; b=1000; k1=180000; k=18000; % unsprung mass % sprung mass % damping ratio % stiffness of Tire % stiffness of sping Displacement of Mass (Transfer function) % Transfer Function of sprung mass displacement num_s=[k1*b k1*k]; den=[m1*m (m1+m)*b [k*(m1+m)+k1*m] k1*b k1*k]; % Transfer Function of sprung mass displacement num_u=[k1*m k1*b k1*k]; Design Considerations (Transfer function) % Transfer Function of sprung mass acceleration num_1=[k1*b k1*k 0 0]; % Transfer Function of suspension deflection num_=[-k1*m 0 0]; % Transfer Function of tire deflection num_3=[-m1*m -(m1+m)*b -k*(m1+m)]; printsys(num_1,den) % print system transfer function
MATLAB Simulation using Laplace Transform Making Input functions t=0:001:0; % 시간을정의 %sine 함수 u1=01*sin(0*t); % sine 함수를이용한자갈길 u=00*sin(4*t)+00*abs(sin(4*t)); % abs() : 절대값함수 % 과속방지턱 u3=005*sin(*pi/0*(t-5))+abs(005*sin(*pi/0*(t-5))); Simulation & Plot result % Linear simulation y=lsim(num_s,den,u1,t); x=lsim(num_u,den,u1,t); plot(t,u1,t,x+03,t,y+1); grid on; title(' 시스템의응답 '); xlabel('time [sec]'); ylabel('displacement'); legend('u(t)','x(t)+03','y(t)+1');
Simulation Results Simulation result(1) Simulation result() : 자갈길 m m m 1 m 1 Simulation result(3) : 과속방지턱 Deflection of Elements m m 1
MATLAB Simulation using State Equation Suspension Parameters m1=55; m=400; b=1000; k1=180000; k=18000; % unsprung mass % sprung mass % damping ratio % stiffness of Tire % stiffness of sping State Equation % Define State eqations A=[ 0 1 0-1; -k/m -b/m 0 b/m; 0 0 0 1 k/m1 b/m1 -k1/m1 -b/m1]; B=[0; 0; -1; 0]; C=[-k/m -b/m 0 b/m; 1 0 0 0; 0 0 1 0]; D=[0; 0; 0]; Making Input functions t=0:001:0; %sine 함수 u1=01*sin(0*t); % 시간을정의 % sine 함수를이용한자갈길 u=00*sin(4*t)+00*abs(sin(4*t)); % abs() : 절대값함수 % 과속방지턱 u3=005*sin(*pi/0*(t-5))+abs(005*sin(*pi/0*(t-5))); % 입력의미분함수구하기 : State equation 의입력 u3_d=diff(u3)/001; u3_d=[u3_d u3_d(length(u3_d))];
MATLAB Simulation using State Equation Simulation & plot result % Linear simulation y=lsim(a,b,c,d,u3_d,t); % n-row, 3-column figure subplot(311) plot(t,y(:,1),'linewidth',); grid on; ylabel('acceleration[m/s^]'); title('design Considerations'); subplot(31) plot(t,y(:,),'linewidth',); grid on; ylabel('suspension deflection[m]'); subplot(313) plot(t,y(:,3),'linewidth',); grid on; ylabel('tire deflection[m]'); xlabel('time[sec]'); y 가 n 행 3 열의데이터
Simulation Results 자갈길 과속방지턱
Parametric Study using for loop Parameters & Input function m1=55; m=400; b=1000; k1=180000; t=0:001:0; % unsprung mass % sprung mass % damping ratio % stiffness of Tire % 시간을정의 % sine 함수를이용한자갈길 u1=00*sin(4*t)+00*abs(sin(4*t)); % abs() : 절대값함수 Parametric study k=[1000 18000 4000 30000]; % 변화시킬 parameter 를배열로정의 for i=1:1:4 % index의설정, 시작 : 간격 : 끝 num1=[k1*b k1*k(i)]; % 전달함수를정의 den1=[m1*m (m1+m)*b [m1*k(i)+(k1+k(i))*m] k1*b k1*k(i)]; num=[k1*m k1*b k1*k(i)]; den=[m1*m (m1+m)*b [m1*k(i)+(k1+k(i))*m] k1*b k1*k(i)]; y(:,i)=lsim(num1,den1,u1,t); x(:,i)=lsim(num,den,u1,t); end % index 에부여하여 matrix 로저장 figure plot(t,y); % matrix의출력 grid on; title(' 차체의응답 (y(t)'); xlabel('time [sec]'); ylabel('displacement'); legend('k_=1000','k_=18000','k_=4000','k_=30000');
Displacement of Sprungmass
Bode Plot Transfer Function Transfer Function % Generates logarithmic points w=logspace(0,3,100); % Bode % [mag,phase,w]=bode(num_s,den,w); [mag,phase,w]=bode(num_u,den,w); % Magnitude magdb = 0*log10(mag); % Plot figure(1) subplot(11) semilogx(w/(*pi),magdb,'b','linewidth',3);hold on; grid on; xlabel('frequency [rad/s]') ylabel('magnitude [db]') subplot(1) semilogx(w/(*pi),phase,'b','linewidth',3);hold on; grid on; xlabel('frequency [rad/s]') ylabel('phase [deg]')
Bode Plot State Space State Space % Generates logarithmic points w=logspace(0,3,100); % Bode % [mag,phase,w]=bode(num_s,den); [mag,phase,w]=bode(num_u,den); % bode(system_num,system_den); % Magnitude magdb = 0*log10(mag); % Plot figure(1) subplot(11) semilogx(w/(*pi),magdb,'b','linewidth',3);hold on; grid on; xlabel('frequency [rad/s]') ylabel('magnitude [db]') subplot(1) semilogx(w/(*pi),phase,'b','linewidth',3);hold on; grid on; xlabel('frequency [rad/s]') ylabel('phase [deg]')
Save and Load the Data 1) Save save 파일명 (Workspace 에있는변수명 ) -ascii Workspace 에있는변수를 ascii code 로된 파일명 으로저장한다 ) Load load 파일명 -ascii ascii code 로된 파일명 으로된파일을같은이름의 Workspace 변수로저장한다 * 단이때문자열이포함되어서는안된다 save resulttxt y -ascii % y 변수를 resulttxt 란파일을만들어서저장 clear all; clc; load resulttxt % resulttxt 를불러옴
Help function 3) Help help 명령어 Workspace 에있는변수를 ascii code 로된 파일명 으로저장한다 save resulttxt y -ascii % y 변수를 resulttxt 란파일을만들어서저장 clear all; clc; load resulttxt % resulttxt 를불러옴