리아푸노프 방정식

수치해석 시뮬레이션 최적설계 마지막 수정: 2026-07-29 04:29:11

1. 개요[편집]

리아푸노프 방정식(Lyapunov equation)은 미지 대칭행렬 $P$ 에 대한 선형 행렬방정식으로, 선형 시스템의 안정성 판정과 에너지 계산에 쓰인다. 연속시간판과 이산시간판이 각각 이렇게 생겼다.

A^{\top}P + PA + Q = 0 \qquad\text{(연속)}, \qquad A^{\top}PA - P + Q = 0 \qquad\text{(이산)}

$A$ 와 $Q$ 가 주어지고 $P$ 를 구한다. 여기서 결정적으로 중요한 사실은 $P$ 에 대해 선형이라는 것이다. 리카티 방정식이 $PBR^{-1}B^\top P$ 라는 이차항 때문에 고유값 문제나 반복법을 동원해야 하는 것과 달리, 리아푸노프 방정식은 원리적으로 연립 1차방정식일 뿐이다. 실제로 리카티를 뉴턴법으로 풀 때 매 반복마다 푸는 부분 문제가 바로 이 방정식이다.

이름 때문에 랴푸노프 지수와 헷갈리기 쉬운데, 둘은 다른 대상이다. 지수 쪽은 비선형 동역학계에서 이웃한 궤적이 벌어지는 속도(카오스 판정)를 재는 스칼라 지표이고, 여기서 다루는 방정식은 선형계의 이차형식 에너지함수를 구성하는 행렬방정식이다. 같은 사람(알렉산드르 리아푸노프, 1892년 학위논문) 이름을 딴 것뿐이다.¹

2. 안정성 정리와의 연결[편집]

이 방정식이 왜 이 모양인지는 리아푸노프의 직접법에서 나온다. $\dot{x} = Ax$ 에 대해 에너지 후보로 이차형식

V(x) = x^{\top}Px, \qquad P \succ 0

를 잡고 궤적을 따라 미분하면

\dot{V} = \dot{x}^{\top}Px + x^{\top}P\dot{x} = x^{\top}\left(A^{\top}P + PA\right)x

가 된다. 여기서 $A^\top P + PA = -Q$ 로 두면 $\dot V = -x^\top Q x$ 이므로, $Q \succ 0$ 이면 $\dot V < 0$ 이 되어 원점이 점근 안정이다. 즉 방정식의 해 $P$ 가 곧 리아푸노프 함수의 계수행렬이다.

정리로 정리하면 이렇다.

$Q \succ 0$ 일 때, 리아푸노프 방정식이 유일한 해 $P \succ 0$ 을 갖는 것과 $A$ 가 후르비츠(모든 고유값의 실수부가 음수)인 것은 동치다. 이산시간에서는 $A$ 가 슈어(모든 고유값의 절댓값이 1 미만)인 것과 동치다.

이 동치성이 실용적으로 대단한 이유는, 안정성 판정을 고유값 계산 없이 할 수 있게 해 주기 때문이다. 아무 $Q \succ 0$ (대개 $Q = I$ )를 잡고 방정식을 푼 뒤 $P$ 가 정부호인지만 촐레스키 분해로 확인하면 끝이다. 게다가 이 논리는 비선형계로 확장된다 — 평형점 근처에서 선형화한 $A$ 로 $P$ 를 구하면 그 $V = x^\top P x$ 가 비선형계의 국소 리아푸노프 함수가 되고, $\dot V < 0$ 이 성립하는 등위집합이 곧 끌림 영역의 보수적 추정치가 된다.

3. 유일성 조건과 적분 표현[편집]

$Q$ 가 정부호가 아니어도 해 자체의 유일성은 따로 논할 수 있다. 연속시간 리아푸노프 방정식은 $A$ 의 고유값 $\lambda_i$ 에 대해

\lambda_i + \lambda_j \neq 0 \quad (\forall i, j)

일 때, 그리고 그때만 유일한 해를 갖는다. 이산시간에서는 조건이 $\lambda_i \lambda_j \neq 1$ 로 바뀐다. 즉 고유값이 허수축에 대칭으로(혹은 단위원에 대해 역수 짝으로) 놓이면 해가 없거나 무한히 많다. 순허수 고유값 한 쌍만 있어도 $\lambda + \bar\lambda = 0$ 이라 바로 걸린다.

$A$ 가 후르비츠이면 해를 적분으로 명시할 수 있다.

P = \int_{0}^{\infty} e^{A^{\top}t}\, Q\, e^{At}\, dt

$A$ 가 안정이므로 지수행렬이 지수적으로 죽어 적분이 수렴한다. 미분해서 대입하면 방정식을 만족함이 바로 확인된다. 이산시간판은 급수

P = \sum_{k=0}^{\infty} \left(A^{\top}\right)^{k} Q\, A^{k}

이다. 이 표현이 물리적 의미를 준다. $Q = C^\top C$ 로 두면 $x_0^\top P x_0 = \int_0^\infty \|y(t)\|^2 dt$ , 즉 초기 상태에서 출발한 자유 응답 출력의 총 에너지다. $P$ 는 에너지 저장량을 재는 행렬이었던 것이다.

4. 수치 해법[편집]

4.1. 벡터화 — 하면 안 되는 방법[편집]

크로네커 곱을 쓰면 이 행렬방정식을 평범한 선형계로 펼 수 있다. $\mathrm{vec}(AXB) = (B^\top \otimes A)\mathrm{vec}(X)$ 항등식을 적용하면

\left(I \otimes A^{\top} + A^{\top} \otimes I\right)\mathrm{vec}(P) = -\,\mathrm{vec}(Q)

가 된다. 계수행렬은 $n^2 \times n^2$ 이고, 그 고유값이 정확히 $\lambda_i + \lambda_j$ 라서 앞서 본 유일성 조건이 “이 행렬이 정칙이다”와 같은 말임을 바로 보여 준다. 이론적 설명으로는 훌륭하다.

문제는 비용이다. $n^2$ 미지수를 가우스 소거법으로 풀면 $O((n^2)^3) = O(n^6)$ 연산에 $O(n^4)$ 메모리다. $n = 100$ 짜리 장난감 문제도 $10^{12}$ 연산에 계수행렬만 $10^8$ 개 성분이다. 대칭성을 이용해 미지수를 $n(n+1)/2$ 로 줄여도 지수는 그대로다. 벡터화는 유도용이고, 계산에는 쓰지 않는다.

4.2. 바텔스-스튜어트 알고리즘[편집]

실전 표준은 1972년 바텔스(R. H. Bartels)와 스튜어트(G. W. Stewart)의 알고리즘이며, $O(n^3)$ 에 끝난다. 전략은 “ $A$ 를 삼각형으로 만들어 놓고 대입해 나간다”이다.

실 슈어 분해 $A = U T U^\top$ 을 계산한다( $U$ 직교, $T$ 준상삼각 — 복소 고유값 쌍이 $2\times2$ 블록으로 남는다). QR 분해 기반 QR 알고리즘으로 $O(n^3)$ .
방정식을 변환한다. $\tilde P = U^\top P U$ , $\tilde Q = U^\top Q U$ 로 두면 $T^\top \tilde P + \tilde P T + \tilde Q = 0$ .
$T$ 가 (준)삼각이므로 $\tilde P$ 의 열(또는 블록)을 후진 대입으로 하나씩 푼다. 각 단계는 $1\times1$ 또는 $2\times2$ 스칼라 방정식이고, 앞서 구한 열들이 우변에 누적된다.
$P = U \tilde P U^\top$ 로 되돌린다.

슈어 분해가 지배적 비용이라 전체가 $O(n^3)$ , 메모리는 $O(n^2)$ 다. LAPACK의 dtrsyl(변환된 삼각계 풀이)과 SLICOT의 SB03MD가 이 경로를 구현하며, MATLAB lyap/dlyap, SciPy solve_lyapunov도 전부 여기에 얹혀 있다. 해머링(Hammarling) 변형은 $Q = C^\top C$ 일 때 $P$ 를 거치지 않고 촐레스키 인수 $R$ ( $P = R^\top R$ )을 직접 계산해, 반올림 때문에 $P$ 가 정부호성을 잃는 사고를 원천 차단한다. 그라미안처럼 조건수가 나쁜 경우에 특히 중요하다.²

4.3. 대규모 문제[편집]

$n$ 이 $10^5$ 를 넘어가는 반이산화 PDE 계에서는 $P$ 를 통째로 저장하는 것 자체가 불가능하다. 이때는 $Q = BB^\top$ 처럼 낮은 계수인 성질을 이용해, 해도 낮은 계수 인수 $P \approx ZZ^\top$ ( $Z$ 는 $n \times r$ , $r \ll n$ )로만 근사한다. 대표적으로 ADI 반복(교대 방향 음함수법)과 그 저계수판인 LR-ADI, 그리고 크리로프 부분공간법에 기반한 사영법이 있다. 그라미안의 고유값이 대개 급격히 감소한다는 사실이 이 근사를 정당화한다.

일반화판인 실베스터 방정식 $AX + XB = C$ 도 같은 슈어-후진대입 구조로 풀리며( $A$ , $B$ 둘 다 슈어화), 리아푸노프 방정식은 $B = A^\top$ 인 특수한 경우다.

5. 그라미안과 응용[편집]

가장 널리 쓰이는 응용은 그라미안(Gramian) 계산이다. $\dot x = Ax + Bu$ , $y = Cx$ 에 대해

A W_c + W_c A^{\top} + BB^{\top} = 0, \qquad A^{\top}W_o + W_o A + C^{\top}C = 0

의 해가 각각 가제어 그라미안 $W_c$ 와 가관측 그라미안 $W_o$ 다. $W_c \succ 0$ 이면 가제어, $W_o \succ 0$ 이면 가관측이라는 판정이 되고, 정량적으로는 $W_c$ 가 “그 상태로 데려가는 데 드는 최소 입력 에너지”의 역행렬을, $W_o$ 가 “그 상태가 출력에 남기는 에너지”를 잰다.

여기서 축소차수모델의 평형 절단(balanced truncation)이 나온다. $W_c$ 와 $W_o$ 가 같은 대각행렬이 되도록 좌표를 바꾸면 그 대각 성분이 한켈 특이값 $\sigma_i = \sqrt{\lambda_i(W_c W_o)}$ 이고, 작은 $\sigma_i$ 에 대응하는 상태는 “밀기도 어렵고 보이지도 않는” 방향이라 잘라낸다. 이 방법이 좋은 이유는 오차에 보장된 상계 $\|G - G_r\|_\infty \le 2\sum_{i>r}\sigma_i$ 가 존재하기 때문이다. 주성분 분석 계열 POD가 데이터에서 에너지를 재는 것과 달리, 여기서는 입출력 관계에서 직접 잰다.

$H_2$ 노름도 리아푸노프 방정식 한 번으로 나온다.

\|G\|_2^2 = \mathrm{tr}\left(C W_c C^{\top}\right) = \mathrm{tr}\left(B^{\top} W_o B\right)

백색 잡음이 들어올 때 출력의 정상상태 분산이 곧 이 값이라, 확률적 성능 지표로 쓰인다. 강건 제어의 $H_2$ 설계와 칼만 필터의 정상상태 공분산 해석이 모두 이 식 위에 서 있다. 참고로 필터 공분산 전파식에서 측정 갱신 항을 빼면 정확히 리아푸노프 방정식이 되고, 측정 갱신을 넣는 순간 이차항이 생겨 리카티가 된다 — 선형과 이차를 가르는 것이 피드백의 유무다.³

6. 관련 문서[편집]

리아푸노프의 1892년 논문 운동 안정성의 일반 문제는 러시아어로 나와 서방에 알려지는 데 반세기가 걸렸다. 1949년 프랑스어 번역, 1992년에야 영어 완역본이 나왔다. 그 사이 제어공학자들은 “안정하려면 고유값을 계산하라”는 말만 반복하고 있었다. ↩
$P$ 를 다 구한 뒤 촐레스키를 돌리면 이미 늦는 경우가 있다. 조건수가 $10^{16}$ 근처인 그라미안은 부동소수점에서 음의 고유값을 갖고 나오며, 이걸 그대로 평형 절단에 넣으면 한켈 특이값이 NaN이 된다. 인수를 직접 전파하는 설계는 칼만 필터의 제곱근 필터와 완전히 같은 발상이다. ↩
그래서 리카티를 뉴턴법으로 풀면 리아푸노프 방정식을 반복마다 하나씩 푸는 형태가 된다. “비선형 문제를 선형 문제의 수열로 바꾼다”는 수치해석의 상투 수단이 여기서도 그대로 나온다. 대신 반복마다 $O(n^3)$ 이 들어서, 슈어 방법으로 한 방에 푸는 쪽이 이기는 경우가 대부분이다. ↩