OKID (Observer Kalman Filter Identification)

OKID(Observer Kalman Filter Identification)는 시간 영역에서 비선형 시스템의 입력-출력 데이터를 이용하여 상태공간 이산시간(discrete-time) 선형 모델을 식별(identification)하는 알고리즘이다.

 

 

OKID는 1990년대 초 NASA의 Juang에 의해 처음 개발된 이래 다양한 항공기 모델을 식별하는 데 이용되어 왔으며, 완벽한 트림 조건이 아닌 경우나 센서 노이즈가 있는 경우에도 매우 효과적으로 모델을 식별할 수 있는 것으로 알려졌다.

OKID는 ERA(eigensystem realization algorithm)의 확장판으로서 ERA 알고리즘이 가진 두 가지 기본 제한 사항을 해결했다.

제한 사항이란 시스템의 초기값이 \(0\) 이어야 한다는 것과 시스템이 안정해야 한다는 것이다. 이러한 제한 사항은 알고리즘의 실제 적용에 있어서 어려움을 초래하는데, 예를 들면 실제 환경에서 항공기의 트림 상태는 일반적으로 돌풍 및 난기류와 같은 작은 대기 교란에도 불완전하여 시스템의 초기값이 \(0\) 이 안될 수가 있으며, 시스템의 안정성이 떨어질 경우 정상 상태로 진입하는데 상당한 시간이 소요되어 과도한 입출력 데이터를 저장하고 처리해야 해야 한다.

그렇다면 시스템의 안정성을 인위적으로 높일 방법은 없을까. 안정성이 높다면, 또는 감쇠가 빠르다면 시스템의 반응에서 초기값의 영향이 크게 감소하고 입출력 데이터의 양도 줄일 수 있으므로 ERA의 문제를 해결할 수 있을 것이다. 이것은 마치 피드백 제어기를 설계하는 문제와 같다. OKID에서는 다음에 설명하는 것처럼 대수적인 조작을 통해 피드백 효과를 줌으로써 안정성을 높이는 효과를 얻는 방법을 사용한다.

미지의 비선형 시스템을 특정 평형점에서 선형화한 모델이 식 (1)과 같이 이산시간 선형 시스템으로 표현된다고 하자.

 

\[ \begin{align} & \mathbf{x}_{k+1}=A \mathbf{x}_k+B \mathbf{u}_k \tag{1} \\ \\ & \mathbf{y}_k=C \mathbf{x}_k+D \mathbf{u}_k \end{align} \]

 

여기서 \(\mathbf{x}_k \in \mathbb{R}^n\), \(\mathbf{u}_k \in \mathbb{R}^p\), \(\mathbf{y}_k \in \mathbb{R}^q\), \(A \in \mathbb{R}^{n \times n}\), \(B \in \mathbb{R}^{n \times p}\), \(C \in \mathbb{R}^{q \times n}\), \(D \in \mathbb{R}^{q \times p}\) 이다.

식 (1)의 상태 업데이트 식에 \(G \mathbf{y}_k\) 를 더하고 빼도 식은 변함이 없지만 다음과 같이 출력 피드백 형식으로 바꿀 수 있다. 이것이 OKID의 핵심 아이디어이다.

 

\[ \begin{align} \mathbf{x}_{k+1} &= A \mathbf{x}_k+B \mathbf{u}_k+G \mathbf{y}_k-G \mathbf{y}_k \tag{2} \\ \\ &=(A+GC) \mathbf{x}_k+(B+GD) \mathbf{u}_k-G \mathbf{y}_k \\ \\ &= \tilde{A} \mathbf{x}_k+ \tilde{B} \mathbf{v}_k \end{align} \]

 

여기서 \(G \in \mathbb{R}^{n \times q}\) 를 관측기 게인(observer gain)이라고 하고,

 

\[ \begin{align} & \tilde{A}=A+GC \ \ \in \mathbb{R}^{n \times n} \tag{3} \\ \\ & \tilde{B}= [B+GD \ \ -G] \ \ \in \mathbb{R}^{n \times (p+q)} \\ \\ & \mathbf{v}_k= \begin{bmatrix} \mathbf{u}_k \\ \mathbf{y}_k \end{bmatrix} \ \ \in \mathbb{R}^{(p+q)} \end{align} \]

 

이다. 관측기 게인 \(G\) 는 \(\tilde{A}\) 를 안정(stable)하게 되도록 선정한다. 식 (1)과 (2)는 수학적으로 동일한 식이다.

이제 시간스텝 \(k=i\) 부터 \(k=i+l_0\) 까지 식 (2)의 상태 업데이트 식을 전개해 보자.

 

\[ \begin{align} \mathbf{x}_{i+1} & = \tilde{A} \mathbf{x}_i+ \tilde{B} \mathbf{v}_i \tag{4} \\ \\ \mathbf{x}_{i+2} & = \tilde{A} \mathbf{x}_{i+1}+ \tilde{B} \mathbf{v}_{i+1} \\ &= \tilde{A}^2 \mathbf{x}_i+ \tilde{A} \tilde{B} \mathbf{v}_i+ \tilde{B} \mathbf{v}_{i+1} \\ \\ & \ \ \ \cdots \\ \\ \mathbf{x}_{i+l_0 } &= \tilde{A} \mathbf{x}_{i+l_0-1}+ \tilde{B} \mathbf{v}_{i+l_0-1} \\ &= \tilde{A}^{l_0 } \mathbf{x}_i+ \tilde{A}^{l_0-1} \tilde{B} \mathbf{v}_i+ \tilde{A}^{l_0-2} \tilde{B} \mathbf{v}_{i+1}+ \cdots + \tilde{B} \mathbf{v}_{i+l_0-1} \end{align} \]

 

 

 

식 (1)의 측정 방정식은 식 (4)에 의해서 시간스텝 \(k=i+l_0\) 에서 다음과 같이 된다.

 

\[ \begin{align} \mathbf{y}_{i+l_0 } &= C \mathbf{x}_{i+l_0 }+D \mathbf{u}_{i+l_0 } \tag{5} \\ \\ &=C \tilde{A}^{l_0 } \mathbf{x}_i+C \tilde{A}^{l_0-1} \tilde{B} \mathbf{v}_i+C \tilde{A}^{l_0-2} \tilde{B} \mathbf{v}_{i+1} \\ & \ \ \ \ + \cdots + C \tilde{B} \mathbf{v}_{i+l_0-1}+D \mathbf{u}_{i+l_0 } \end{align} \]

 

 

 

이제 \(i=0, 1, \cdots , l-1\) 에서 또는 시간스텝 \(k=l_0, l_0+1, \cdots , l_0+l-1\) 에서 식 (5)를 풀어 써보면 다음과 같다.

 

\[ \begin{align} \mathbf{y}_{l_0 } &= C \tilde{A}^{l_0 } \mathbf{x}_0+ C \tilde{A}^{l_0-1} \tilde{B} \mathbf{v}_0+C \tilde{A}^{l_0-2} \tilde{B} \mathbf{v}_1 \tag{6} \\ & \ \ \ \ \ + \cdots + C\tilde{B} \mathbf{v}_{l_0-1}+D \mathbf{u}_{l_0 } \\ \\ \mathbf{y}_{l_0+1} &=C \tilde{A}^{l_0 } \mathbf{x}_1+C \tilde{A}^{l_0-1} \tilde{B} \mathbf{v}_1+C \tilde{A}^{l_0-2} \tilde{B} \mathbf{v}_2 \\ & \ \ \ \ \ + \cdots + C\tilde{B} \mathbf{v}_{l_0 }+D \mathbf{u}_{l_0+1} \\ \\ \mathbf{y}_{l_0+2} &=C \tilde{A}^{l_0 } \mathbf{x}_2+C \tilde{A}^{l_0-1} \tilde{B} \mathbf{v}_2+C \tilde{A}^{l_0-2} \tilde{B} \mathbf{v}_3 \\ & \ \ \ \ \ + \cdots + C \tilde{B} \mathbf{v}_{l_0+1}+D \mathbf{u}_{l_0+2} \\ \\ & \ \ \ \cdots \\ \\ \mathbf{y}_{l_0+l-1} &=C \tilde{A}^{l_0 } \mathbf{x}_{l-1}+C \tilde{A}^{l_0-1} \tilde{B} \mathbf{v}_{l-1}+C \tilde{A}^{l_0-2} \tilde{B}\mathbf{v}_l \\ & \ \ \ \ \ + \cdots + C \tilde{B} \mathbf{v}_{l_0+l-2}+D \mathbf{u}_{l_0+l-1} \end{align} \]

 

 

 

식 (6)을 행렬 형식으로 정리하면 다음과 같다.

 

\[ \mathbf{y}_{l_0:l_0+l-1}=C \tilde{A}^{l_0 } \mathbf{x}_{0:l-1}+ \bar{Y}V \tag{7} \]

여기서

\[ \begin{align} & \mathbf{y}_{l_0:l_0+l-1}=[ \mathbf{y}_{l_0} \ \ \ \mathbf{y}_{l_0+1} \ \ \ \mathbf{y}_{l_0+2} \ \ \ \cdots \ \ \ \mathbf{y}_{l_0+l-1} ] \ \ \in \mathbb{R}^{q \times l} \tag{8} \\ \\ & \mathbf{x}_{0:l-1} =[ \mathbf{x}_0 \ \ \ \mathbf{x}_1 \ \ \ \mathbf{x}_2 \ \ \ \cdots \ \ \ \mathbf{x}_{l-1} ] \ \ \in \mathbb{R}^{n \times l} \\ \\ & \bar{Y} = [D \ \ \ C \tilde{B} \ \ \ C\tilde{A}\tilde{B} \ \ \ \cdots \ \ \ C\tilde{A}^{l_0-2} \tilde{B} \ \ \ C\tilde{A}^{l_0-1} \tilde{B} ] \ \ \in \mathbb{R}^{q \times [(p+q) l_0+p]} \\ \\ & V= \begin{bmatrix} \mathbf{u}_{l_0 } & \mathbf{u}_{l_0+1} & \mathbf{u}_{l_0+2} & \cdots & \mathbf{u}_{l_0+l-1} \\ \mathbf{v}_{l_0-1} & \mathbf{v}_{l_0 } & \mathbf{v}_{l_0+1} & \cdots & \mathbf{v}_{l_0+l-2} \\ \vdots & \cdots & \vdots & \ddots & \vdots \\ \mathbf{v}_1 & \mathbf{v}_2 & \mathbf{v}_3 & \cdots & \mathbf{v}_l \\ \mathbf{v}_0 & \mathbf{v}_1 & \mathbf{v}_2 & \cdots & \mathbf{v}_{l-1} \end{bmatrix} \ \ \mathbb{R}^{[(p+q) l_0+p] \times l} \end{align} \]

 

이다. \(\bar{Y}\) 는 관측기 마코프(Markov) 파라미터로 구성된 행렬이다.

여기서 \(l_0\) 를 \(\tilde{A}^{l_0 } \approx 0\) 이 되도록 충분히 크게 잡으면 식 (7)은 다음과 같이 된다.

 

\[ \mathbf{y}_{l_0:l_0+l-1} \approx \bar{Y}V \tag{9} \]

 

\( \mathbf{y}_{l_0:l_0+l-1}\) 와 \(V\) 는 입출력 데이터로만 이루어진 행렬이고 마코프 파라미터 행렬 \(\bar{Y}\) 는 시스템을 구성하는 행렬로서 추정할 대상이다.

식 (9)에 의하면 초기값에 의한 영향이 사라진다. 또한 식 (9)에 의하면 총 식의 개수는 \(q \times l\) 이고 추정해야 할 미지수의 개수는 \(q \times [(p+q) l_0+p]\) 이다. 따라서 데이터의 길이 \(l\) 을 \((p+q) l_0+p\) 이상으로 잡으면, 즉 \(l \ge (p+q) l_0+p\) 이면 식의 개수가 미지수의 개수보다 커지므로 행렬 \(\bar{Y}\) 를 다음과 같이 근사적으로 계산할 수 있다.

 

\[\begin{align} \bar{Y} &= \mathbf{y}_{l_0:l_0+l-1} V^+ \tag{10} \\ \\ &= \mathbf{y}_{l_0:l_0+l-1} V^T (VV^T )^{-1} \end{align} \]

 

여기서 \(V^+\) 는 \(V\) 의 유사 역행렬(pseudo inverse matrix)이다. 행렬 \(V\) 의 유사 역행렬이 존재하려면 시스템의 입력이 모든 운동 모드를 여기(excitation)할 정도로 충분히 복잡해야 한다.

참고로 항공기 식별을 위한 비행 테스트에서는 스텝(step), 펄스(pulse) 또는 싱글렛(singlet), 더블렛(doublet) 등 세 가지 유형의 제어 입력이 많이 사용된다.

 

 

이제 식 (8)을 참고하여 마코프 파라미터 행렬 \(\bar{Y}\) 를 다음과 같이 분할한다.

 

\[\begin{align} \bar{Y} &= [D \ \ \ C\tilde{B} \ \ \ C\tilde{A}\tilde{B} \ \ \ \cdots \ \ \ C\tilde{A}^{l_0-2} \tilde{B} \ \ \ C\tilde{A}^{l_0-1} \tilde{B} ] \tag{11} \\ \\ &= [ \bar{Y}_0 \ \ \ \bar{Y}_1 \ \ \ \bar{Y}_2 \ \ \ \cdots \ \ \ \bar{Y}_{l_0-1} \ \ \ \bar{Y}_{l_0 } ] \end{align} \]

 

여기서 식 (3)을 이용하면 \(\bar{Y}\) 의 구성 성분을 다음과 같이 쓸 수 있다.

 

\[\begin{align} \bar{Y}_0 &= D \tag{12} \\ \\ \bar{Y}_i &=C \tilde{A}^{i-1} \tilde{B} \\ \\ &=[ C(A+GC)^{i-1} (B+GD) \ \ \ \ -C(A+GC)^{i-1} G] \\ \\ &=[\tilde{\bar{Y}}_i^{(1) } \ \ \ \ -\tilde{\bar{Y}}_i^{(2) } ] \end{align} \]

 

그러면 식 (12)에서 다음과 같이 원래 시스템의 마코프 파라미터를 계산해 낼 수 있다.

 

\[\begin{align} Y_1 &= CB=C(B+GD)-CGD \tag{13} \\ &= \tilde{\bar{Y}}_1^{(1) }- \tilde{\bar{Y}}_1^{(2) } D \\ \\ Y_2 &= CAB=C(A+GC)(B+GD)-CGCB-C(A+GC)GD \\ &= \tilde{\bar{Y}}_2^{(1) }- \tilde{\bar{Y}}_1^{(2) } Y_1- \tilde{\bar{Y}}_2^{(2) } D \end{align} \]

 

식 (13)을 이용하면 모든 원래 시스템의 마코프 파라미터를 계산할 수 있다.

 

\[\begin{align} Y_0 &= D= \bar{Y}_0 \tag{14} \\ \\ Y_i & =CA^{i-1} B \\ \\ &=\tilde{\bar{Y}}_i^{(1) }- \sum_{k=1}^{l_0} \tilde{\bar{Y}}_k^{(2)} Y_{i-k} \end{align} \]

 

식 (14)에서 계산한 마코프 파라미터를 Ho-Kalman 식별 알고리즘에 대입하면 미지의 시스템 (1)의 행렬 \(A, B, C, D\) 와 시스템 차수를 식별할 수 있다. 이와 같은 식별 방법을 OKID라고 한다.

OKID 의 장점은 적절한 관측기 게인 \(G\) 를 선정하여 시스템의 안정성을 증대 시킴으로써 감쇠가 작은 시스템에서도 마코프 파라미터의 수를 줄일 수 있고 이를 통해 더 작은 항켈(Hankel) 행렬을 사용할 수 있으므로 식별 알고리즘의 계산이 용이해 진다는 것이다.