7.43 연립 상미분 방정식과 상태 공간 표현

1. 연립 상미분 방정식의 정의

연립 상미분 방정식(system of ordinary differential equations)은 여러 미지 함수와 그 도함수가 상호 결합된 미분 방정식의 집합이다. n개의 미지 함수 x_1(t), x_2(t), \ldots, x_n(t)에 대한 1차 연립 상미분 방정식의 일반적 형태는 다음과 같다.

\dot{x}_1 = f_1(t, x_1, x_2, \ldots, x_n)

\dot{x}_2 = f_2(t, x_1, x_2, \ldots, x_n)

\vdots

\dot{x}_n = f_n(t, x_1, x_2, \ldots, x_n)

이를 벡터 형태로 간결하게 표현하면 다음과 같다.

\dot{\mathbf{x}}(t) = \mathbf{f}(t, \mathbf{x}(t))

여기서 \mathbf{x}(t) = [x_1(t), x_2(t), \ldots, x_n(t)]^T \in \mathbb{R}^n은 상태 벡터(state vector)이고, \mathbf{f}: \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^n은 벡터장(vector field)이다.

2. 고차 미분 방정식의 연립 방정식 변환

임의의 n차 상미분 방정식은 n개의 1차 연립 상미분 방정식으로 변환할 수 있다. n차 방정식

y^{(n)} = F(t, y, \dot{y}, \ddot{y}, \ldots, y^{(n-1)})

에 대하여 다음과 같이 상태 변수를 정의한다.

x_1 = y, \quad x_2 = \dot{y}, \quad x_3 = \ddot{y}, \quad \ldots, \quad x_n = y^{(n-1)}

이렇게 하면 원래의 고차 방정식은 다음의 1차 연립 방정식으로 변환된다.

\dot{x}_1 = x_2

\dot{x}_2 = x_3

\vdots

\dot{x}_{n-1} = x_n

\dot{x}_n = F(t, x_1, x_2, \ldots, x_n)

이 변환은 로봇 동역학에서 2차 운동 방정식을 수치적으로 풀기 위한 표준적 전처리 단계이다.

3. 선형 연립 상미분 방정식

3.1 동차 선형 시스템

상수 계수를 갖는 동차 선형 연립 상미분 방정식은 다음과 같다.

\dot{\mathbf{x}}(t) = A\mathbf{x}(t)

여기서 A \in \mathbb{R}^{n \times n}은 시스템 행렬(system matrix)이다. 초기 조건 \mathbf{x}(0) = \mathbf{x}_0에 대한 해는 행렬 지수함수(matrix exponential)를 이용하여 다음과 같이 표현된다.

\mathbf{x}(t) = e^{At}\mathbf{x}_0

행렬 지수함수는 급수 전개로 정의된다.

e^{At} = I + At + \frac{(At)^2}{2!} + \frac{(At)^3}{3!} + \cdots = \sum_{k=0}^{\infty} \frac{(At)^k}{k!}

이 행렬 \Phi(t) = e^{At}를 상태 천이 행렬(state transition matrix)이라 하며, 다음 성질을 만족한다.

\Phi(0) = I, \quad \dot{\Phi}(t) = A\Phi(t), \quad \Phi(t_1 + t_2) = \Phi(t_1)\Phi(t_2), \quad \Phi^{-1}(t) = \Phi(-t)

3.2 고유값 분해를 이용한 해법

시스템 행렬 A의 고유값과 고유벡터를 이용하여 연립 방정식의 해를 구할 수 있다. A가 n개의 서로 다른 고유값 \lambda_1, \lambda_2, \ldots, \lambda_n을 갖고 대응하는 고유벡터가 \mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_n이면, 일반해는 다음과 같다.

\mathbf{x}(t) = c_1 e^{\lambda_1 t}\mathbf{v}_1 + c_2 e^{\lambda_2 t}\mathbf{v}_2 + \cdots + c_n e^{\lambda_n t}\mathbf{v}_n

여기서 상수 c_1, c_2, \ldots, c_n은 초기 조건에 의해 결정된다. A가 대각화 가능한 경우, A = PDP^{-1} (여기서 D = \text{diag}(\lambda_1, \ldots, \lambda_n), P = [\mathbf{v}_1 \; \mathbf{v}_2 \; \cdots \; \mathbf{v}_n])이면 행렬 지수함수는 다음과 같이 계산된다.

e^{At} = P e^{Dt} P^{-1} = P \, \text{diag}(e^{\lambda_1 t}, \ldots, e^{\lambda_n t}) \, P^{-1}

3.3 중복 고유값의 경우

고유값 \lambda가 대수적 중복도 m을 갖되 기하적 중복도가 m보다 작은 경우, 일반화된 고유벡터(generalized eigenvector)를 도입하여야 한다. 조르당 정규형(Jordan normal form) A = PJP^{-1}에서 \lambda에 대응하는 조르당 블록이 크기 k이면, 해당 성분의 해는 다음 형태를 갖는다.

e^{\lambda t}\left(\mathbf{v} + t\mathbf{w}_1 + \frac{t^2}{2!}\mathbf{w}_2 + \cdots + \frac{t^{k-1}}{(k-1)!}\mathbf{w}_{k-1}\right)

여기서 \mathbf{w}_1, \mathbf{w}_2, \ldots, \mathbf{w}_{k-1}은 일반화된 고유벡터이다.

3.4 비동차 선형 시스템

비동차 선형 연립 방정식

\dot{\mathbf{x}}(t) = A\mathbf{x}(t) + \mathbf{b}(t)

의 해는 상태 천이 행렬을 이용한 변환 공식(variation of constants formula)으로 다음과 같이 표현된다.

\mathbf{x}(t) = e^{At}\mathbf{x}_0 + \int_0^t e^{A(t - \tau)}\mathbf{b}(\tau) \, d\tau

첫째 항 e^{At}\mathbf{x}_0은 초기 조건에 의한 자유 응답(free response)이고, 둘째 항은 외부 입력 \mathbf{b}(t)에 의한 강제 응답(forced response)이다. 강제 응답은 행렬 지수함수와 입력의 합성곱(convolution)으로 구성된다.

4. 상태 공간 표현

4.1 표준 형태

상태 공간 표현(state-space representation)은 입력-출력 관계를 포함하는 동적 시스템의 내부 표현 방식이다. 선형 시불변(LTI) 시스템의 상태 공간 표현은 다음 두 방정식으로 구성된다.

\dot{\mathbf{x}}(t) = A\mathbf{x}(t) + B\mathbf{u}(t) \quad \text{(상태 방정식)}

\mathbf{y}(t) = C\mathbf{x}(t) + D\mathbf{u}(t) \quad \text{(출력 방정식)}

여기서 각 변수와 행렬의 의미는 다음과 같다.

n을 시스템의 차수(order) 또는 상태 변수의 수라 하며, 이는 시스템의 동적 자유도를 나타낸다.

4.2 상태 변수의 물리적 의미

로봇 시스템에서 상태 변수는 시스템의 현재 상태를 완전히 기술하는 최소한의 변수 집합이다. 단일 관절 로봇의 경우, 관절 각도 \theta와 관절 각속도 \dot{\theta}를 상태 변수로 선택하는 것이 자연스럽다.

\mathbf{x} = \begin{bmatrix} \theta \\ \dot{\theta} \end{bmatrix}

n자유도 로봇 매니퓰레이터의 경우, 상태 벡터는 2n차원이 된다.

\mathbf{x} = \begin{bmatrix} \mathbf{q} \\ \dot{\mathbf{q}} \end{bmatrix} \in \mathbb{R}^{2n}

여기서 \mathbf{q} \in \mathbb{R}^n은 일반화 좌표 벡터이고, \dot{\mathbf{q}} \in \mathbb{R}^n은 일반화 속도 벡터이다.

4.3 상태 공간 표현의 비유일성

동일한 입출력 관계를 갖는 시스템에 대하여 상태 공간 표현은 유일하지 않다. 비특이 행렬 T \in \mathbb{R}^{n \times n}에 의한 상태 변환 \mathbf{z} = T\mathbf{x}를 수행하면, 변환된 상태 공간 표현은 다음과 같다.

\dot{\mathbf{z}} = TAT^{-1}\mathbf{z} + TB\mathbf{u}

\mathbf{y} = CT^{-1}\mathbf{z} + D\mathbf{u}

변환된 시스템 행렬 \bar{A} = TAT^{-1}, \bar{B} = TB, \bar{C} = CT^{-1}, \bar{D} = D는 원래 시스템과 동일한 입출력 특성을 갖는다. 이러한 변환을 유사 변환(similarity transformation)이라 하며, 시스템의 고유값은 유사 변환에 의해 불변이다.

5. 로봇 동역학의 상태 공간 변환

5.1 차 운동 방정식으로부터의 변환

로봇 매니퓰레이터의 운동 방정식은 일반적으로 다음의 2차 형태를 갖는다.

M(\mathbf{q})\ddot{\mathbf{q}} + C(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \boldsymbol{\tau}

여기서 M(\mathbf{q})는 관성 행렬, C(\mathbf{q}, \dot{\mathbf{q}})는 코리올리 및 원심력 행렬, \mathbf{g}(\mathbf{q})는 중력 벡터, \boldsymbol{\tau}는 관절 토크 벡터이다.

상태 변수를 \mathbf{x} = [\mathbf{q}^T, \dot{\mathbf{q}}^T]^T로 정의하면, 상태 방정식은 다음과 같다.

\dot{\mathbf{x}} = \begin{bmatrix} \dot{\mathbf{q}} \\ M(\mathbf{q})^{-1}\left[\boldsymbol{\tau} - C(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} - \mathbf{g}(\mathbf{q})\right] \end{bmatrix}

이는 비선형 상태 방정식 \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \mathbf{u})의 형태이며, 여기서 입력 \mathbf{u} = \boldsymbol{\tau}이다. 관성 행렬 M(\mathbf{q})가 양정치이므로 그 역행렬은 항상 존재하며, 상태 방정식이 명시적(explicit) 형태로 표현 가능함이 보장된다.

5.2 선형화와 선형 상태 공간 모델

비선형 상태 방정식을 평형점(equilibrium point) (\mathbf{x}_e, \mathbf{u}_e) 주위에서 1차 테일러 전개하면 선형 근사 모델을 얻는다. 편차 변수 \delta\mathbf{x} = \mathbf{x} - \mathbf{x}_e, \delta\mathbf{u} = \mathbf{u} - \mathbf{u}_e에 대하여 다음이 성립한다.

\delta\dot{\mathbf{x}} \approx A\delta\mathbf{x} + B\delta\mathbf{u}

여기서 야코비안 행렬 A와 B는 다음과 같이 정의된다.

A = \frac{\partial \mathbf{f}}{\partial \mathbf{x}}\bigg\vert_{(\mathbf{x}_e, \mathbf{u}_e)}, \quad B = \frac{\partial \mathbf{f}}{\partial \mathbf{u}}\bigg\vert_{(\mathbf{x}_e, \mathbf{u}_e)}

이 선형화된 상태 공간 모델은 평형점 근방에서의 로봇 시스템 거동을 분석하고, 선형 제어 기법(예: LQR, 극배치법)을 적용하기 위한 기초가 된다.

6. 안정성 해석과 고유값의 역할

6.1 고유값과 시스템 안정성

선형 시불변 시스템 \dot{\mathbf{x}} = A\mathbf{x}의 안정성은 시스템 행렬 A의 고유값에 의해 완전히 결정된다.

• 모든 고유값의 실수부가 음이면(\text{Re}(\lambda_i) < 0, \forall i), 시스템은 점근 안정(asymptotically stable)하다.
• 모든 고유값의 실수부가 영 이하이고, 실수부가 영인 고유값의 대수적 중복도와 기하적 중복도가 같으면, 시스템은 안정(stable)하다.
• 실수부가 양인 고유값이 하나라도 존재하면, 시스템은 불안정(unstable)하다.

6.2 고유값의 물리적 해석

고유값이 \lambda = \sigma + j\omega일 때, 실수부 \sigma는 응답의 감쇠율을 결정하고, 허수부 \omega는 진동 주파수를 결정한다.

• \sigma < 0: 감쇠 응답 (시간 상수 \tau = -1/\sigma)
• \sigma > 0: 발산 응답
• \omega \neq 0: 주파수 \omega/(2\pi)의 진동 성분 존재

로봇 관절 시스템에서 고유값 분석은 과도 응답의 감쇠 특성과 진동 주파수를 예측하는 데 직접적으로 활용된다.

7. 이산 시간 상태 공간 표현

디지털 제어 시스템에서는 연속 시간 모델을 이산화(discretization)하여 사용한다. 샘플링 주기 T_s에 대한 이산 시간 상태 공간 모델은 다음과 같다.

\mathbf{x}[k+1] = A_d \mathbf{x}[k] + B_d \mathbf{u}[k]

\mathbf{y}[k] = C_d \mathbf{x}[k] + D_d \mathbf{u}[k]

영차 유지(zero-order hold, ZOH) 방식으로 이산화하면, 이산 시스템 행렬은 다음과 같이 결정된다.

A_d = e^{AT_s}, \quad B_d = \left(\int_0^{T_s} e^{A\tau} d\tau\right) B = A^{-1}(A_d - I)B

여기서 두 번째 등식은 A가 비특이일 때 성립한다. 이산 시간 시스템의 안정성은 이산 시스템 행렬 A_d의 고유값 크기가 모두 1보다 작은지(\vert\lambda_i\vert < 1, \forall i)에 의해 판정된다.

8. 요약

연립 상미분 방정식은 다자유도 로봇 시스템의 동적 거동을 기술하는 기본 틀이다. 고차 미분 방정식을 1차 연립 방정식으로 변환함으로써 통일된 수학적 체계를 갖추게 되며, 이것이 상태 공간 표현의 토대가 된다. 상태 공간 표현은 시스템 행렬, 입력 행렬, 출력 행렬, 직접 전달 행렬의 네 가지 행렬로 시스템을 완전히 기술하며, 다입력 다출력(MIMO) 시스템으로의 자연스러운 확장이 가능하다. 행렬 지수함수와 고유값 분석은 시스템의 시간 응답과 안정성을 체계적으로 파악하는 핵심 도구이며, 로봇 제어기 설계의 이론적 근간을 형성한다.

참고 문헌

• Khalil, H. K. (2002). Nonlinear Systems (3rd ed.). Prentice Hall.
• Chen, C. T. (1999). Linear System Theory and Design (3rd ed.). Oxford University Press.
• Ogata, K. (2010). Modern Control Engineering (5th ed.). Prentice Hall.
• Craig, J. J. (2005). Introduction to Robotics: Mechanics and Control (3rd ed.). Pearson.
• Siciliano, B., Sciavicco, L., Villani, L., & Oriolo, G. (2009). Robotics: Modelling, Planning and Control. Springer.
• Strang, G. (2016). Introduction to Linear Algebra (5th ed.). Wellesley-Cambridge Press.

version: 0.1.0