MATLAB을 사용하여 바퀴 형 도립 진자의 운동 방정식 도출

본 기사에서는 MATLAB/Symbolic Math Toolbox를 사용하여 바퀴형 도립진자의 운동방정식을 도출해 나갑니다.
이 기사를 통해 MATLAB을 이용하여 복잡한 미분 방정식을 풀 수 있다는 것을 실감해 주시면 좋겠습니다.

전치

트랜지스터 기술 7월호 에 차륜형 도립 진자의 수학 모델의 상세한 도출이 쓰여져 있으므로, 이 기사에서는 도중식의 해설은 하지 않습니다. 7월호에는 수학이나 역학에 관한 정중한 해설이 많이 막혀 있으므로, 1권 가지고 두어 손해는 없다고 생각합니다. 사자. 웃음

바퀴 형 도립 진자의 각 파라미터

MATLAB에 도출해 주는 운동 방정식에 포함되는 「위치 에너지」와 「운동 에너지」에 대해서, 「진자」와 「바퀴」의 2개의 관점으로부터 도출해 갑니다.

바퀴 형 도립 진자의 물리량

〇바퀴(wheel)의 물리량

바퀴의 무게 중심 좌표를 $(x_w, y_w)$로 한다.

무게 중심 좌표

\left\{
\begin{array}{ll}
x_w = r_w(\theta_p + \theta_w) \\
y_w = 0
\end{array}
\right.

무게 중심 속도

\left\{
\begin{array}{ll}
\dot{x}_w = r_w(\dot{\theta}_p + \dot{\theta}_w) \\
\dot{y}_w = 0
\end{array}
\right.

병진 운동 에너지

K_{w1} = \frac{1}{2}m_w\dot{x}_w^2 +\frac{1}{2}m_w\dot{y}_w^2 \\
= \frac{1}{2}m_wr_w^2(\dot{\theta}_p + \dot{\theta}_w)^2

회전 운동 에너지

K_{w2} = \frac{1}{2}I_w(\dot{\theta}_p + \dot{\theta}_w)^2

위치 에너지

U_w = m_wgy_w = 0

〇진자(pendulum)의 물리량

진자의 무게 중심 좌표를 $(x_p, y_p)$로 한다.

무게 중심 좌표

\left\{
\begin{array}{ll}
x_p = x_w + r_psin(\theta_p) = r_w(\theta_p + \theta_w) + r_psin(\theta_p) \\
y_p = r_pcos(\theta_p)
\end{array}
\right.

무게 중심 속도

\left\{
\begin{array}{ll}
\dot{x}_p = r_w(\dot{\theta}_p + \dot{\theta}_w) + r_p\dot{\theta}_pcos(\theta_p) \\
\dot{y}_p = -r_p\dot{\theta}_psin(\theta_p)
\end{array}
\right.

병진 운동 에너지

K_{p1} = \frac{1}{2}m_p\dot{x}_p^2 + \frac{1}{2}m_p\dot{y}_p^2 \\
= \frac{1}{2}m_p\bigl(r_w(\dot{\theta}_p + \dot{\theta}_w) + r_p\dot{\theta}_pcos(\theta_p)\bigr)^2 + \frac{1}{2}m_p\bigl(-r_p\dot{\theta}_psin(\theta_p)\bigr)^2

회전 운동 에너지

K_{p2} = \frac{1}{2}I_p\dot{\theta}_p^2

위치 에너지

U_p = m_pgy_p = m_pgr_pcos(\theta_p)

라그랑주 운동 방정식의 도출

모델의 운동 에너지와 위치 에너지의 도출이 가능하면,
MATLAB/Symbolic Math Toolbox를 사용하여, 바퀴형 도립 진자의 라그랑주의 운동 방정식을 도출해 갑니다.

lagrange.m

clear; format compact   %initialize

syms m_w r_w theta_w dtheta_w ddtheta_w I_w real % wheel
syms m_p r_p theta_p dtheta_p ddtheta_p I_p real % pendulum
syms theta_w dtheta_w ddtheta_w phi_p dphi_p ddphi_p real % radian
syms f_w g real % torque, gravity

q   = [theta_w theta_p]';
dq  = [dtheta_w dtheta_p]';
ddq = [ddtheta_w ddtheta_p]';
f   = [f_w 0]';

%--------------------------------------------------------

x_w = r_w*(q(1) + q(2));
y_w = 0;

dx_w = diff(x_w, q(1))*dq(1) + diff(x_w, q(2))*dq(2);
dy_w = diff(x_w, q(1))*dq(1) + diff(y_w, q(2))*dq(2);

x_p = x_w + r_p*sin(q(2));
y_p = r_p*cos(q(2));

dx_p = diff(x_p, q(1))*dq(1) + diff(x_p, q(2))*dq(2);
dy_p = diff(x_p, q(1))*dq(1) + diff(y_p, q(2))*dq(2);

%--------------------------------------------------------

K_w1 = 0.5*m_w*dx_w^2 + 0.5*m_w*dy_w^2;
K_w2 = 0.5*I_w*(dq(1) + dq(2))^2;
U_w = m_w*g*y_w;

K_p1 = 0.5*m_p*dx_p^2 + 0.5*m_p*dy_p^2;
K_p2 = 0.5*I_p*dq(2)^2;
U_p = m_p*g*y_p;

L = (K_w1 + K_w2 + K_p1 + K_p2) - (U_w + U_p);

%--------------------------------------------------------

N = length(q);
for i = 1:N
    dLq(i) = diff(L, dq(i));

    temp = 0;
    for j = 1:N
        temp = temp + diff(dLq(i),dq(j))*ddq(j) + diff(dLq(i),q(j))*dq(j);
    end
    ddLq(i) = temp;

    eq(i) = ddLq(i) - diff(L,q(i)) - f(i);
end

eq = simplify(eq')

m 파일을 실행하면 바퀴와 진자에 대한 라그랑주 운동 방정식을 도출할 수 있습니다.

출력 결과

>> lagrange
eq =
                                                                                                            ddtheta_p*(I_w + m_w*r_w^2 + m_p*r_w*(r_w + r_p*cos(theta_p)) - m_p*r_p*r_w*sin(theta_p)) - dtheta_p*(dtheta_p*m_p*r_p*r_w*cos(theta_p) + dtheta_p*m_p*r_p*r_w*sin(theta_p)) - f_w + ddtheta_w*(I_w + 2*m_p*r_w^2 + 2*m_w*r_w^2)
 I_p*ddtheta_p + I_w*ddtheta_p + I_w*ddtheta_w + ddtheta_p*m_p*r_p^2 + ddtheta_p*m_p*r_w^2 + ddtheta_w*m_p*r_w^2 + ddtheta_p*m_w*r_w^2 + ddtheta_w*m_w*r_w^2 - g*m_p*r_p*sin(theta_p) - dtheta_p^2*m_p*r_p*r_w*sin(theta_p) + 2*ddtheta_p*m_p*r_p*r_w*cos(theta_p) + ddtheta_w*m_p*r_p*r_w*cos(theta_p) - ddtheta_w*m_p*r_p*r_w*sin(theta_p)

MATLAB의 출력 결과는 수식으로 볼 수 없기 때문에 markdown의 수식 모드를 사용하여 성형했습니다.

위의 식이 바퀴, 아래의 식이 진자의 운동 방정식을 나타냅니다.

\ddot{\theta}_p(I_w + m_wr_w^2 + m_pr_w(r_w + r_pcos(\theta_p)) - m_pr_pr_wsin(\theta_p)) - \dot{\theta}_p(\dot{\theta}_pm_pr_pr_wcos(\theta_p) + \dot{\theta}_pm_pr_pr_wsin(\theta_p)) - f_w + \ddot{\theta}_w(I_w + 2m_pr_w^2 + 2m_wr_w^2) = 0\\

I_p\ddot{\theta}_p + I_w\ddot{\theta}_p + I_w\ddot{\theta}_w + \ddot{\theta}_pm_pr_p^2 + \ddot{\theta}_pm_pr_w^2 + \ddot{\theta}_wm_pr_w^2 + \ddot{\theta}_pm_wr_w^2 + \ddot{\theta}_wm_wr_w^2 - m_pgr_psin(\theta_p) - \dot{\theta}_p^2m_pr_pr_wsin(\theta_p) + 2\ddot{\theta}_pm_pr_pr_wcos(\theta_p) + \ddot{\theta}_wm_pr_pr_wcos(\theta_p) - \ddot{\theta}_wm_pr_pr_wsin(\theta_p) = 0

참고문헌

트랜지스터 기술 7월호

도립 진자로 배우는 제어 공학