ERA (Eigensystem Realization Algorithm)

Ho-Kalman 식별 알고리즘에서는 임펄스 반응(impulse response)을 이용하여 마코프 파라미터(Markov parameters)를 측정하였다. 그렇다면 일반적인 입출력 데이터를 이용하여 마코프 파라미터를 획득하는 방법은 없을까.

 

 

다음과 같이 미지의 이산시간(discrete-time) 선형 시스템이 있다고 하자.

 

$$\begin{array}{r}\text{(1)}& & {\mathbf{x}}_{k+1}=A{\mathbf{x}}_k+B{\mathbf{u}}_k\\ & {\mathbf{y}}_k=C{\mathbf{x}}_k+D{\mathbf{u}}_k\end{array}$$

 

여기서 ${\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}$ 이다.

시간스텝 $k$ 에서 출력은 다음과 같이 계산할 수 있다.

 

$$\begin{array}{r}\text{(2)}& {\mathbf{y}}_0& =C{\mathbf{x}}_0+D{\mathbf{u}}_0\\ {\mathbf{y}}_1& =C{\mathbf{x}}_1+D{\mathbf{u}}_1\\ & =C(A{\mathbf{x}}_0+B{\mathbf{u}}_0)+D{\mathbf{u}}_1\\ & =CA{\mathbf{x}}_0+CB{\mathbf{u}}_0+D{\mathbf{u}}_1\\ \\ {\mathbf{y}}_2& =C{\mathbf{x}}_2+D{\mathbf{u}}_2\\ & =C(A{\mathbf{x}}_1+B{\mathbf{u}}_1)+D{\mathbf{u}}_2\\ & =C{A}^2{\mathbf{x}}_0+CAB{\mathbf{u}}_0+CB{\mathbf{u}}_1+D{\mathbf{u}}_2\\ \\ & \cdots \\ \\ {\mathbf{y}}_k& =C{A}^k{\mathbf{x}}_0+\sum _{i=0}^{k-1}C{A}^{k-1-i}B{\mathbf{u}}_i+D{\mathbf{u}}_k\end{array}$$

 

초기값이 $0$ 이라고 가정하면, 즉 ${\mathbf{x}}_0=0$ 이면 시간스텝 $k$ 에서 출력은 다음과 같이 된다.

 

$$\begin{array}{}\text{(3)}& {\mathbf{y}}_k=\sum _{i=0}^{k-1}C{A}^{k-1-i}B{\mathbf{u}}_i+D{\mathbf{u}}_k\end{array}$$

 

식 (3)을 $k=0,1,...,l-1$ 에서 행렬 형식으로 쓰면 다음과 같다.

 

$$\begin{array}{}\text{(4)}& {\mathbf{y}}_{0:l-1}=YU\end{array}$$

여기서

$$\begin{array}{r}\text{(5)}& & {\mathbf{y}}_{0:l-1}=[{\mathbf{y}}_0{\mathbf{y}}_1{\mathbf{y}}_2\cdots {\mathbf{y}}_{l-1}]\in {\mathbb{R}}^{q\times l}\\ & Y=[DCBCAB\cdots C{A}^{l-2}B]\in {\mathbb{R}}^{q\times pl}\\ \\ & U=\left[\begin{array}{ccccc}{\mathbf{u}}_0& {\mathbf{u}}_1& {\mathbf{u}}_2& \cdots & {\mathbf{u}}_{l-1}\\ & {\mathbf{u}}_0& {\mathbf{u}}_1& \cdots & {\mathbf{u}}_{l-2}\\ & & {\mathbf{u}}_0& \cdots & {\mathbf{u}}_{l-3}\\ & & & \ddots & ⋮\\ & & & & {\mathbf{u}}_0\end{array}\right]\in {\mathbb{R}}^{pl\times l}\end{array}$$

 

이다. $U$ 와 ${\mathbf{y}}_{0:l-1}$ 는 입출력 데이터로만 이루어진 행렬이고 마코프 파라미터로 구성된 행렬 $Y$ 는 추정할 대상이다.

식 (5)에 의하면 총 식의 개수는 $q\times l$ 이고 추정해야 할 미지수 개수는 $q\times pl$ 이다. 따라서 미지수의 개수가 식의 개수보다 많기 때문에 $p=1$ 즉 입력이 $1$ 개일 때를 제외하고는 $Y$ 는 무수히 많은 해가 존재하며 유일하게 계산해 낼 수 없다. 입력이 1개 이더라도 $Y={\mathbf{y}}_{0:l-1}{U}^{-1}$ 이어야 하므로 $U^{-1}$ 가 존재해야 하기 때문에 ${\mathbf{u}}_0\ne 0$ 이어야 하고 입력이 충분히 복잡해야 한다.

하지만 만약 행렬 $A$ 가 안정(stable)하여서 어떤 시간스텝 $k\ge l_0$ 에서는 $A^k\approx 0$ 이 된다면 $k\ge l_0$ 구간에서는 식 (3)을 다음과 같이 근사화할 수 있다.

$$\begin{array}{r}\text{(6)}& {\mathbf{y}}_k& =\sum _{i=0}^{k-1}C{A}^{k-1-i}B{\mathbf{u}}_i+D{\mathbf{u}}_k\\ & =\sum _{i=0}^{k-l_0-1}C{A}^{k-1-i}B{\mathbf{u}}_i+\sum _{i=k-l_0}^{k-1}C{A}^{k-1-i}B{\mathbf{u}}_i+D{\mathbf{u}}_k\\ \\ & \approx \sum _{i=k-l_0}^{k-1}C{A}^{k-1-i}B{\mathbf{u}}_i+D{\mathbf{u}}_k,k\ge l_0\end{array}$$

 

식 (6)에 의하면 식 (4)와 (5)는 다음과 같이 근사화 된다.

 

$$\begin{array}{}\text{(7)}& {\mathbf{y}}_{0:l-1}\approx YU\end{array}$$

여기서

$$\begin{array}{r}\text{(8)}& & {\mathbf{y}}_{0:l-1}=[{\mathbf{y}}_0{\mathbf{y}}_1{\mathbf{y}}_2\cdots {\mathbf{y}}_{l_0}\cdots {\mathbf{y}}_{l-1}]\in {\mathbb{R}}^{q\times l}\\ & Y=[DCBCAB\cdots C{A}^{l_0-1}B]\in {\mathbb{R}}^{q\times p(l_0+1)}\\ \\ & U=\left[\begin{array}{ccccccc}{\mathbf{u}}_0& {\mathbf{u}}_1& {\mathbf{u}}_2& \cdots & {\mathbf{u}}_{l_0}& \cdots & {\mathbf{u}}_{l-1}\\ & {\mathbf{u}}_0& {\mathbf{u}}_1& \cdots & {\mathbf{u}}_{l_0-1}& \cdots & {\mathbf{u}}_{l-2}\\ & & {\mathbf{u}}_0& \cdots & {\mathbf{u}}_{l_0-2}& \cdots & {\mathbf{u}}_{l-3}\\ & & & \ddots & ⋮& \cdots & ⋮\\ & & & & {\mathbf{u}}_0& \cdots & {\mathbf{u}}_{l-l_0-1}\end{array}\right]\in {\mathbb{R}}^{p(l_0+1)\times l}\end{array}$$

 

이다. 식 (8)에 의하면 총 식의 개수는 $q\times l$ 이고 추정해야 할 미지수 개수는 $q\times p(l_0+1)$ 이다. 따라서 식 (8)에서 데이터의 길이 $l$ 을 $p(l_0+1)$ 이상으로 잡으면, 즉 $l\ge p(l_0+1)$ 이면 식의 개수가 미지수의 개수보다 커지므로 행렬 $Y$ 를 다음과 같이 근사적으로 계산할 수 있다.

$$\begin{array}{}\text{(9)}& Y={\mathbf{y}}_{0:l-1}{U}^{+}\end{array}$$

 

여기서 $U^{+}$ 는 $U$ 의 유사 역행렬(pseudo inverse matrix)이다. 행렬 $U$ 는 입력 데이터 ${\mathbf{u}}_i$ 로 구성되어 있으므로 ${\mathbf{u}}_i$ 는 $U$ 의 유사 역행렬이 존재할 수 있도록 충분히 복잡하게 설계되어야 한다.

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

 

 

식 (6)에 의하면 근사화에 의한 오차는 $l_0$ 가 커질 수록 작아질 것이다. 하지만 $U$ 의 사이즈도 같이 커지기 때문에 유사 역행렬 $U^{+}$ 를 계산하는데 어려움이 따른다. 만약 시스템이 덜 안정적일 경우, 즉 감쇄가 느리게 진행될 경우, 정수 $l_0$ 를 큰 값으로 정해야 하기 때문에 ERA 는 실용성이 떨어진다.

그렇다면 인위적으로 시스템의 안정성을 높일 방법은 없을까. 그리고 초기값이 $0$ 이어야 한다는 제약을 푸는 방법은 없을까.