Question
SH 선배님께 받은 Easy_Model 분석과 관련해 질문드립니다.
MATLAB, RecurDyn 파일과 함께 설명 자료도 함께 첨부해주셔서 참고해 분석하며 MATLAB을 이용한 기구학 분석을 익히는 Exercise를 진행 중입니다.
코드를 보며 이해하는 중 질문이 몇 가지 생겨, 승현 선배 혹은 파일을 확인한 후 답변 해주실 수 있는 다른 분들께 질문드립니다.
질문 1: a_Easy_Kinematic_Model을 통해 어떠한 시스템인지 이해했습니다. 하지만 b_Easy_Energy로 넘어오며 분석의 목적을 파악하는데 어려움이 생겼습니다. PPT 속 과정을 통해 최종적으로 추출하려는 값들이 무엇인지 알고싶습니다. 특히 My_RK4_easy의 역할을 이해하는데 도움이 필요합니다 (MATLAB에 정말 미숙하기 때문에 전체적인 흐름을 더 자세히 설명해주시면 감사하겠습니다).
질문 2: b_Easy_Energy의 Lagrange Equation에 대입하는 과정에서 Link1의 A(원점)에 위치한 spring의 potential energy를 계산하게 되는데, 이 스프링은 a_Easy_Kinematic_Model에서는 불필요해 언급이 없는 것으로 이해하면 되는걸까요?
% 질문 2
PE = (1/2*K*(alpha - deg2rad(60))^2); % PE stored in the torsional spring at A
KE1 = (1/2*I1*alpha_dot^2); % KE of link1(rot)
KE2 = (1/2*I2*beta_dot^2 + 1/2*m2*(v_M(1)^2 + v_M(2)^2)); % KE of link2(rot+trans)
KE = KE1 + KE2;
MATLAB
복사
Answer
Q 1. [b_Easy_Energy]를 통해 최종적으로 추출하려는 값과 [My_RK4_easy]에 대한 설명
최종적으로 추출하려는 값
[a_Easy_Kinematic_Model]에 나타나는 모델 (Two-Link Model)은 alpha라는 각도 변수의 값이 결정되면, 두 링크의 상태가 결정되는 모델임을 알 수 있음 (beta 또한 alpha에 대한 함수이기 때문에).
즉, 해당 코드는 Two-Link Model의 링크 위치를 alpha에 대해 기술하고, alpha의 값을 변화해가면서 기구가 어떻게 거동하는지를 기구학 해석 (Kinematics Analysis)을 통해 확인한 것임.
기구학 해석을 통해 Two-Link Model이 적절한 거동을 취함을 확인하였고,
이후 동역학 해석 (Dynamics Analysis)을 수행하여 힘을 입력으로 넣어주었을 때 해당 시스템이 어떻게 거동하는지 분석하는 과정이 필요함.
동역학 해석을 수행하기 위해서는 운동방정식 (Equation of Motion)을 기술해야 함.
Two-Link Model 처럼 둘 이상의 물체가 움직이는 경우, Lagrangian Mechanics에 기반하여 운동방정식을 기술하는 것이 일반적임.
이 방법으로 기술된 운동방정식을 Euler-Lagrange Equation이라고 말하며 다음과 같이 정의됨:
여기서 은 라그랑지안 (Lagrangian), 는 일반화 좌표 (Generalized Coordinates), 는 일반화 외력 (Generalized External Forces)이다.
일반화 좌표는 해당 시스템의 상태를 결정하는 변수라고 볼 수 있으며, Two-Link Model에서는 alpha가 이에 해당함.
Euler-Lagrange Equation을 이용해서 운동방정식을 기술하기 위해서는 해당 시스템의 라그랑지안을 구해야 하며, 라그랑지안은 다음과 같이 정의됨:
여기서 는 운동 에너지 (Kinematic Energy)이고 는 포텐셜 에너지 (Potential Energy)이다.
[b_Easy_Energy]는 Two-Link Model의 운동 에너지와 포텐셜 에너지를 MATLAB Symbolic을 이용하여 수식으로 정의하고 (코드 맨 윗줄에 syms를 이용하여 Symbolic 변수가 정의 됨),
PE = (1/2*K*(alpha - deg2rad(60))^2);% Total PE: stored in the torsional spring at A
KE1 = (1/2*I1*alpha_dot^2); % KE of link1 (rot)
KE2 = (1/2*I2*beta_dot^2 + 1/2*m2*(v_M(1)^2 + v_M(2)^2)); % KE of link2 (rot + trans)
KE = KE1 + KE2; % Total KE
MATLAB
복사
Symbolic 미분과 관련된 함수 (diff)를 이용하여 Euler-Lagrange Equation을 구하는 스크립트임.
L_eq = KE - PE; % Lagrangian by definition
L_eq1 = diff(L_eq, alpha_dot); % Generalized momentum conjugate to alpha [ d L / d q_dot ]
L_eq2 = diff(L_eq1, alpha)*alpha_dot + diff(L_eq1, alpha_dot)*alpha_2dot; % [ d / d t ] with chain rule
L_eq3 = diff(L_eq, alpha); % [ d L / d q ]
Eq = collect(L_eq2 - L_eq3, [alpha_2dot, alpha_dot]); % Sort L_eq2 by derivs of alpha
[coeffs_A_B, terms] = coeffs(Eq, alpha_2dot); % Coeffcients of alpha_2dot
Q = diff(r_C(1), alpha)*(-F); % Generalized external force
MATLAB
복사
즉, 해당 스크립트를 실행시킴으로써 최종적으로 Two-Link Model에 대한 운동방정식을 Euler-Lagrange Equation 형태로 얻게 됨.
변수 Eq가 Euler-Lagrange Equation의 좌변이고 Q가 우변임.
[My_RK4_easy]에 대한 설명
우선 [My_RK4_easy] 함수 스크립트는 기본적으로 가독성이 떨어져서 이해하기 어려운게 맞음.
좀비 코드를 만들어낸 승현이를 탓하자
아마 지금은 좀 더 깔끔하게 코드 짤 줄 알테지만, 과거에 대충 짜둔거 수정하기 귀찮아서 그대로 둔 듯…
우선 [My_RK4_easy] 함수 스크립트는 크게 두 개의 함수로 이루어졌으며, 첫 번째 함수 [My_RK4_easy]는 내가 MATLAB 강의 때도 설명 했고 일반적으로 수치 해석 수업에서 배우는 4차 Runge-Kutta (RK4)임.
function x_new = My_RK4_easy(x, h, F)
% x: current state variables, h: time step size, F: external force input
% x_new: next state variables
k1 = Dyn_Eqn(x, F);
k2 = Dyn_Eqn(x + 0.5*k1*h, F);
k3 = Dyn_Eqn(x + 0.5*k2*h, F);
k4 = Dyn_Eqn(x + k3*h, F);
x_new = x + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h;
end
MATLAB
복사
아마 이해하기 어려운 부분은 운동방정식을 기술해둔 [Dyn_Eqn] 일 것인데,
여기서 링크 길이나 MOI (Moment of Inertia) 같은 변수의 특정 값을 정의하는 부분 싹 다 날리고, 핵심적인 부분만 보면 다음과 같음:
function f = Dyn_Eqn(x, F)
% State Variable: alpha & alpha_dot
alpha = x(:,1);
alpha_dot = x(:,2);
% State Matrix
A = I1 + 0.5000*m2*(4.5000*L1^2*sin(alpha)^2 - (0.5000*L1^4*cos(alpha)^2*sin(alpha)^2)/(L2^2*((L1^2*cos(alpha)^2)/L2^2 - 1))) - (I2*L1^2*sin(alpha)^2)/(L2^2*((L1^2*cos(alpha)^2)/L2^2 - 1));
B = (0.5000*m2*(9*L1^2*cos(alpha)*sin(alpha) + (L1^4*cos(alpha)*sin(alpha)^3)/(L2^2*((L1^2*cos(alpha)^2)/L2^2 - 1)) - (L1^4*cos(alpha)^3*sin(alpha))/(L2^2*((L1^2*cos(alpha)^2)/L2^2 - 1)) - (L1^6*cos(alpha)^3*sin(alpha)^3)/(L2^4*((L1^2*cos(alpha)^2)/L2^2 - 1)^2)) - 0.5000*m2*(4.5000*L1^2*cos(alpha)*sin(alpha) + (0.5000*L1^4*cos(alpha)*sin(alpha)^3)/(L2^2*((L1^2*cos(alpha)^2)/L2^2 - 1)) - (0.5000*L1^4*cos(alpha)^3*sin(alpha))/(L2^2*((L1^2*cos(alpha)^2)/L2^2 - 1)) - (0.5000*L1^6*cos(alpha)^3*sin(alpha)^3)/(L2^4*((L1^2*cos(alpha)^2)/L2^2 - 1)^2)) - (I2*L1^2*cos(alpha)*sin(alpha))/(L2^2*((L1^2*cos(alpha)^2)/L2^2 - 1)) - (I2*L1^4*cos(alpha)*sin(alpha)^3)/(L2^4*((L1^2*cos(alpha)^2)/L2^2 - 1)^2))*alpha_dot^2 + 0.5000*K*(2*alpha - 2.0944);
Q = 2*F*L1*sin(alpha);
% State Equation
f(:,1) = alpha_dot;
f(:,2) = -B/A + Q/A;
end
MATLAB
복사
즉, MATLAB에서 [My_RK4_easy]를 이용하여 수치적으로 운동방정식을 풀기 위해,
앞서 수립한 운동방정식을 상태 방정식 형태 (State-Space Form)로 다시 써둔 것에 불과함을 알 수 있음.
상태 변수에 해당하는 x는 alpha와 alpha_dot으로 이루어져 있음.
A는 운동방정식에서 alpha_2dot의 계수에 해당하며 [b_Easy_Energy] 에서 이미 구했으니 복사하여 사용함.
B는 운동방정식에서 alpha_2dot의 계수 이외의 항에 해당하며 이 또한 [b_Easy_Energy] 에서 이미 구했으니 복사하여 사용함.
Q는 운동방정식에서 일반화 외력에 해당하며 [b_Easy_Energy] 에서 이미 구했으니 복사하여 사용함 (다만 힘의 방향을 고려하여 F 부호 결정).
f는 운동방정식의 상태 방정식 형태에 해당 하며 이에 대한 정의를 잘 모를 경우 별도의 공부가 필요함 (원본 코드에서 if 문 위 아래 내용이 같아 생략해도 무방).
[My_RK4_easy]와 같은 과정이 운동방정식을 수치 적분을 통해 풀어나가는 가장 기본적인 형태임.
Q 2. [a_Easy_Kinematic_Model]에서 스프링 존재에 따른 포텐셜 에너지 관련 내용이 없는 이유
[a_Easy_Kinematic_Model]의 목적은 단순히 기구학 해석 (Kinematics Analysis)에 있음!
앞서 Q 1 답변 초반에 이야기 했던 것처럼, [a_Easy_Kinematic_Model]는 기구가 어떻게 거동하는지 분석하는 기구학 해석이 주된 목적임.
승현이 답변
Kinematic에서는 단순히 구동하는 데에 있어서 크기나 길이를 확인하거나 각도가 접힐때 간섭이 안생기는지 등등을 파악했던거라서 spring이 들어가지 않음.
근데, 이제 스프링을 넣고 리커다인과 동일한 거동을 보고싶으면 저 lagrange를 runge kutta 해서 시간에 따른 각도를 구하고 이 각도를 kinematic model.m에 다시 넣으면 됨.
즉, 동역학 해석을 수행하기에 앞서, 구상하고 있는 시스템이 정의한 변수 (alpha)가 변화함에 따라 어떻게 거동하는 지 시각적으로 확인하기 위한 스크립트라고 볼 수 있음.
힘을 고려하지 않고 alpha 값에 따른 두 링크의 위치와 자세가 어떻게 변하는지 확인하기 위한 코드이기에 스프링 존재 여부는 중요하지 않음.
동역학 해석 결과로 얻어지는 시간에 따른 alpha 프로파일을 이용하여, 동역학 해석 결과를 시각화 하기 위해 사용할 수 있긴함.