27.3.3.1. 자동 생성된 24x24 상태 전이 행렬의 희소(Sparse) 연산 코드 구조
EKF의 공분산 예측 단계(predictCovariance)는 센서 노이즈와 시간의 흐름에 따라 기체 상태의 불확실성이 어떻게 팽창하는가를 연산하는 핵심 과정이다. 수학적으로 이 팽창은 24개의 기체 상태 변수들이 서로 얽혀 있는 미분 방정식의 야코비안 행렬(Jacobian Matrix, \mathbf{F})을 통해 계산된다.
본 절에서는 576(24 \times 24)개의 요소를 가진 거대한 상태 전이 행렬 연산이, 연산량이 제한된 임베디드 코어(예: STM32)에서 어떻게 실시간으로 처리되는지, 파이썬(SymPy) 기반의 스크립트로 자동 생성(Autogenerated)된 C++ 희소(Sparse) 연산 코드의 구조적 특징을 분석한다.
1. 컴퓨팅 부하의 근원: 꽉 찬 전배열(Dense Matrix) 연산의 파국
리카티(Riccati) 미분 방정식에 기반한 공분산 팽창 수식은 다음과 같이 근사된다.
\mathbf{P}_{k+1} = \mathbf{F}_k \mathbf{P}_k \mathbf{F}_k^T + \mathbf{Q}
- \mathbf{P}_k: 현재 상태의 24 \times 24 공분산 행렬
- \mathbf{F}_k: 24 \times 24 상태 전이 야코비안 행렬 (State Transition Jacobian Matrix)
- \mathbf{Q}: 24 \times 24 프로세스 노이즈 행렬
만약 \mathbf{F}_k 행렬을 일반적인 2차원 부동소수점 배열(float F[24][24])로 메모리에 할당하고 3중 for 문을 돌려 행렬 곱셈 연산(\mathbf{O}(N^3))을 수행한다면, 루프 한 번에 대략 24 \times 24 \times 24 \approx 13,824 번의 부동소수점 곱셈-덧셈(FMA: Fused Multiply-Add) 연산이 필요하다.
더 큰 문제는 \mathbf{F}_k 안의 요소들을 채워 넣는 과정 자체다. 예를 들어 위치(Position)를 쿼터니언(Quaternion)의 각 항(q_0, q_1, q_2, q_3)으로 편미분한 결과식 안에는 삼각함수 파생형과 제곱 연산이 수십 개 얽혀 있다. 이를 무식하게 실시간 연산으로 돌리면 Pixhawk의 CPU 점유율은 100%를 찍고 제어 루프가 붕괴(Watchdog Timeout)될 것이다.
2. 희소성(Sparsity) 규칙과 0 행렬요소(Zero-Entry) 가지치기
다행스럽게도 EKF의 상태 변수 간에는 ‘전혀 물리적 연관성이 없는’ 변수 쌍이 다수 존재한다.
예를 들어, “지구 자기장의 세기(인덱스 16~18)“가 변한다고 해서 “가속도계의 센서 바이어스(인덱스 13~15)“가 변할 리가 없다. 따라서 이 두 변수를 서로 편미분한 야코비안 행렬의 해당 셀(Cell) 값은 무조건 수학적으로 0 이다.
PX4 개발팀은 24 \times 24 행렬(576개 요소를 가짐) 내부를 뜯어본 결과, 물리적 연관성이 없어 영구히 0 값을 갖는 요소가 전체의 80% 이상을 차지하는 극단적인 희소 행렬(Sparse Matrix) 형태임을 파악했다.
0이 곱해지면 결과는 무조건 0이며, 이 연산을 CPU 클럭을 써가며 수행해 봤자 아무런 의미가 없다. 따라서 PX4 알고리즘 개발의 철학은 “0이 되는 위치는 아예 코드에 쓰지도 않고 건너뛴다” 로 수렴한다.
3. SymPy 심볼릭 연산 기반 C++ 코드 오토제너레이션(Auto-Generation)
PX4 ECL(Estimation and Control Library) 폴더 안 깊숙한 곳(src/lib/ecl/EKF/python/)에는 파이썬(Python)으로 작성된 기괴한 스크립트들이 존재한다. 이 스크립트들은 컴파일러가 아니라, 수학자 역할을 하는 봇(Bot)이다.
3.1 파이썬 스크립트의 동작 원리
이 스크립트는 파이썬의 대수학 라이브러리인 SymPy (Symbolic Python) 를 호출한다.
- 기체의 24상태 비선형 동역학 수식(쿼터니언 회전, 중력, 가속도 변환 등)을 수학 기호(Symbol) 문자열로 입력한다.
- SymPy 라이브러리에게 24 \times 24 편미분(Jacobian) 연산을 시킨다.
- \mathbf{F}_k \mathbf{P}_k \mathbf{F}_k^T 수식을 기호 상으로 전부 전개(Expand)한 뒤, 인수가 0으로 날아가는 항을 대수학적으로 모조리 소거한다.
- 살아남은 수식들을 분석하여 공통으로 반복되는 부분식(Common Subexpression)을 묶어 임시 변수(
t시리즈)로 치환한다. - 최종적으로 이 최적화된 수학 공식을
C/C++문법의 문자열 텍스트로 인쇄(Generate)하여.cpp소스 코드 텍스트 파일로 저장시킨다.
3.2 자동 생성된 타원형(Elliptical) 하드코딩 코드 분석
실제 predictCovariance 함수 내의 코드를 보면 개발자가 짰다고는 상상할 수 없는 기계적인 하드코딩(Hard-coding) 블록이 펼쳐져 있다. for문 같은 반복 루프는 전혀 보이지 않고, 오직 평면적인 행렬 직접 접근 연산만 576줄(혹은 그 이상) 늘어서 있다.
// 자동 생성된 희소 연산 기반의 공분산 예측 코드 패턴
// 1. 공통 부분식(Common Subexpressions, CSE)을 위한 임시 스택 변수 선언
// 쿼터니언과 델타 앵글 값들이 교차 결합되는 공통 수식들을 미리 한 번만 계산하여 캐싱함
const float t0 = q0 * q0;
const float t1 = q1 * q1;
const float t2 = q2 * q2;
const float t3 = q3 * q3;
const float t10 = q0 * q1 * 2.0f;
// ... (t0 부터 t400 이상까지 무수한 서브 수식 변수들 나열) ...
// 2. 루프문(For-loop)이 완전히 제거된 직접 대입 연산 (Loop Unrolling)
// 예측 후의 공분산 행렬 nextP의 개별 원소(0행 0열부터)에 직접 계산 값 꽂아넣기
nextP(0, 0) = P(0,0) + dt*( P(0,0)*t25 + P(0,1)*t26 /*...중략...*/ );
// 야코비안이 0인 항은 애초에 코드 생성 시 누락됨 (Sparsity 이용)
// 예를 들어, 자기장 편미분 항이 0인 인덱스의 경우 수식 자체가 짧거나 아예 P 값을 복사만 함
nextP(16, 16) = P(16, 16) + dt * q_mag_noise; // 다른 상태변수와 연관이 없다면 이렇게 단순해짐!
// 상삼각(Upper-triangular) 행렬만 연산 (대칭 행렬의 절반만 연산하여 속도 2배 향상)
nextP(1, 0) = nextP(0, 1); // 미리 계산한 상단 값을 덮어쓰기만 수행
4. 구조적 장단점 및 최적화 철학 소결
이런 오토제너레이션 기반의 하드코딩 기법은 C++ 코드로서 일반적인 가독성(Readability)은 최악이다. 개발자가 에러를 발견해도 C++ 코드를 직접 고칠 수 없고, 원본 Python SymPy 스크립트 수식을 수정한 후 다시 코드를 제너레이트(Generate)해야 하는 구조적 단절성을 초래한다.
하지만 그 대가로 얻은 극강의 런타임 성능(Runtime Performance) 은 모든 단점을 압도한다.
반복문을 타면서 발생하는 분기 예측 실패(Branch Misprediction)나 점프 오버헤드가 제로(0)가 되며, \mathbf{O}(N^3) 의 행렬 곱셈 부하가 비단순 덧셈/곱셈 수백 번으로 축소되어 시간 복잡도가 사실상 임베디드 O(1) 상수에 수렴한다.
이는 “비행의 안전과 실시간성(Real-time Operation)은 그 어떤 소프트웨어 공학적 우아함보다도 최우선시된다“라는 무인항공기 생태계의 철학을 가감 없이 보여주는 설계다. 이 무식하지만 가장 확실하고 빠른 코드의 사슬 덕분에, EKF 필터는 \mathbf{P} 행렬을 부풀리는 무거운 짐을 매번 지치지 않고 들어 올릴 수 있다.
다음 절에서는 \mathbf{O}(N^3) 의 부담을 털고 일어난 이 거대한 공분산 행렬이, 유한한 부동소수점(float) 늪 속에서 대칭성(Symmetry)을 수호하기 위해 펼치는 소프트웨어적 구명 활동 로직에 대해 해부한다.