해밀턴-자코비-벨만 방정식 1 (Deterministic case) (feat. LQR)

1. Optimal Control problem

최적 제어(Optimal Control)는 아래의 그림에서 제어 입력(control) $\mathbf{u}(t)$을 디자인하여 cost $J$를 최소화시키면서 상태(state) $\mathbf{x}(t_0)$에서 $\mathbf{x}(t_f)$까지 도달하는 것을 목표로 한다.

여기서 위 시스템의 dynamics와 cost function을 다음과 같이 상정해보자.

$$
\begin{align}
    \frac{d}{dt}\mathbf{x}&=f(\mathbf{x}(t), \mathbf{u}(t)) \label{eq1} \tag{1} \\
    J(\mathbf{x}(t), \mathbf{u}(t), t_0)&=Q(\mathbf{x}(t_f), t_f)+\int_{t_0}^{t_f} \mathcal{L}(\mathbf{x}(\tau), \mathbf{u}(\tau)) d\tau \label{eq2} \tag{2} \\
\end{align}
$$

(2)의 $Q$는 말기 상태에 도달했을 때에 Terminal cost이고 $\mathcal{L}$은 도달하면서 받을 일시적 cost이다. 
참고로 (1)에서 $f$가 linear이고 (2)의 $J$가 'Sum of Square'처럼 quadratic form인 케이스가 바로 LQR(Linear Quadratic Regulator) 문제이다.

그리고 아래와 같이 각 시점에서의 최적 제어를 찾아 목표에 달성한 optimal cost function의 존재에 대해 생각해볼 수 있다. 

$$V(\mathbf{x}(t), t) = \min_{\mathbf{u}(t)} J(\mathbf{x}(t), \mathbf{u}(t), t) \label{eq3} \tag{3}$$


제어 이론에서는 이 optimal cost function를 value function 혹은 cost-to-go function 이라고 칭한다. 시스템이 위의 optimal cost를 달성할 수 있으면 최적화 프로세스를 끝내고 실제의 제어 입력을 투입할 수 있기 때문에 "cost-to-go"라고 명명했다고 한다.
개념의 혼동을 막기위해 덧붙이자면 강화학습에서는 $J$를 value function, $V$를 optimal value function이라 칭한다. 여기서는 제어 이론의 관점에서 서술하고 있으므로 제어 이론의 용어를 따르기로 한다. 즉 여기서 칭하는 value function은 강화학습의 optimal value function과 동치이다.

2. Hamilton-Jacobi-Bellman(HJB) Equation

해밀턴 자코비 벨만 방정식(HJB equation)은 제어 입력의 최적성(optimality)에 대해 필요 충분 조건을 제시함으로써 최적 제어 이론에서 현 시점 가장 중요한 방정식으로 자리매김하고 있다. 일반적으로 해밀턴 자코비 벨만 방정식은 식 (3)의 value function에 대한 비선형 편미분 방정식 형태를 띄고 있다. 즉 이 방정식의 해는 value function 자신이 된다는 의미이다. 일단 이 해가 알려지면 해밀턴 자코비 방정식 내부의 해밀토니안(Hamiltonian)에 대해 최대 or 최솟값을 구해서 최적 제어를 얻을 수 있게 된다. 

---

* Note : "Hamiltonian" in Pontryagin's Minimum(Maximum) Principle

해밀턴 자코비 벨만 방정식을 논의하기 전에 먼저 최적 제어 문제에서 사용되는 해밀토니안(Hamiltonian)의 존재에 대해 살펴볼 필요가 있다. 여기서 해밀토니안은 식 (1)과 같은 dynamical system을 위한 최적 제어 문제를 해결하기 위해 사용되는 함수로서 Lev Pontryagin의 Minimum(Maximum) principle에서 처음 등장한다. 식 (2)의 $J$를 cost function으로 정의하면 *Minimum* principle이고 reward function으로 정의하면 *Maximum* Principle이다. 위에서 cost로 정의했기 때문에 여기선 Minimum principle로 논의를 전개하도록 한다. Pontryagin Minimum principle은 여기서 주 논의 대상이 아니기 때문에 간단히만 언급하자면, 최적 제어 문제는 Hamiltonian을 최적화하도록 제어를 선택하면 자연스레 해결된다는 것을 증명한 이론이다. 

$$
H(\mathbf{x}(t), \mathbf{u}(t), \mathbf{\lambda}(t),t) = \mathbf{\lambda}^{T}(t)f(\mathbf{x}(t), \mathbf{u}(t),t) + \mathcal{L}(\mathbf{x}(t), \mathbf{u}(t)) \label{eq4} \tag{4}
$$

$$
H(\mathbf{x}^*(t), \mathbf{u}^*(t), \lambda^*(t), t) \leq H(\mathbf{x}^*(t), \mathbf{u}(t), \lambda^*(t), t) \label{eq5} \tag{5}
$$

여기서 $\mathbf{x}^*$는 최적 상태 trajectory고 $\lambda^{*}$는 최적 라그랑지안(lagrangian) 승수 벡터이다. 식 (5)에서의 함의는 제어 입력 $\mathbf{u}$를 제외한 다른 input을 넣은 작거나 같은 해밀토니안이 있다면 그 해밀토니안에 들어간 제어 입력 $\mathbf{u}^*$가 $\mathbf{u}$에 비해 더 최적화된 제어 입력이라는 뜻이다.

식 (4)가 바로 해밀토니안으로서 라그랑지안 형식으로 표현된 최적화 대상 함수(objective function)이다.  이 노트에선 이러한 해밀토니안의 개념만 킵해두고 가자.

---

2.1 Hamilton-Jacobi-Bellman(HJB) Equation formula

Richard E. Bellman은 이 최적화된 해밀토니안이 value function의 시간에 대한 미분값에 음수를 취한 것과 같아야함을 제시하였는데, 이것이 바로 해밀턴 자코비 벨만 방정식(HJB Equation)이다.

$$
\begin{align}
-\frac{\partial V}{\partial t} 
&= H_{opt} \\
&= \min_{\mathbf{u}(t)} \left ( \lambda^T(t) f(\mathbf{x}, \mathbf{u})+\mathcal{L}(\mathbf{x}, \mathbf{u}) \right ) \\
&= \min_{\mathbf{u}(t)} \left ( \left( \frac{\partial V}{\partial x} \right )^T f(\mathbf{x}, \mathbf{u})+\mathcal{L}(\mathbf{x}, \mathbf{u}) \right ) \label{eq6} \tag{6}
\end{align}
$$

아래에서 이 해밀턴 자코비 벨만 방정식이 어떻게 유도되었는지 살펴보도록 한다.

2.2 Derivation of the HJB Equation

기술할 해밀턴 자코비 벨만 방정식 유도를 위해선 Dynamic programming과 테일러 전개에 대한 사전 지식이 필요하다. 이것들을 간단하게 리뷰해보자.

---

* Note : Dynamic Programming과 Bellman Optimality

Dynamic programming은 Bellman Optimality로부터 비롯되는 알고리즘이라고 할 수 있다.
Bellman Optimality는 위에서 봤던 최적 제어 문제를 아래의 그림처럼 임의의 중간 지점 $t$를 두어도 $t_0$ ~ $t$까지의 value function과 $t$ ~ $t_f$까지의 value function의 합으로 $t_0$~$t_f$의 value function을 만들 수 있다는 성질을 뜻한다.

식 (3)의 value function 표현에서 terminal point까지 변수로 만들어 수식으로 표현하면 아래와 같다.

$$
V(\mathbf{x}(t_0), t_0, t_f)=V(\mathbf{x}(t_0), t_0, t)+V(\mathbf{x}(t), t, t_f) \label{eq7} \tag{7}
$$

이 Bellman Optimality를 일반화시키면 어떠한 길고 복잡한 trajectory를 갖는 system을 상대적으로 더 간단한 작은 time step으로 쪼개어 풀어도 최적성이 유지됨을 시사한다. 그리고 이것이 정확히 dynamic programming의 프로세스다.

---

---

* Note: 테일러 전개(Taylor Expansion)

여기선 테일러 전개에 대해 간단하게 알아보고자 한다.  
테일러 전개는 특정 지점에서 잘 모르는 함수에 대해 근사하여 알아내는 수학적 도구이다. 만약 1) target이 되는 지점과 다른 source에서의 함수 값을 알고, 2) 그 source에서의 n계 도함수 값을 알 수 있다면 테일러 전개를 활용하여 target 지점의 함수값을 근사할 수 있다. 여기서 target과 source의 거리가 멀어질수록, 근사의 정확도를 위해 더 높은 n계 도함수가 필요하다.
해당 함수를 $f$, target 지점을 $a$, source지점을 $x$라 할 때 수식 표현은 다음과 같다.

$$
\begin{align}
    f(a) 
    &= \sum_{n=0}^{\infty}\frac{f^{(n)}(x)}{n!}(a-x)^n \\
    &= f(x)+ \frac{f^{'}(x)}{1!}(a-x) + \frac{f^{''}(x)}{2!}(a-x)^2 + \cdots \label{eq8} \tag{8}\\
\end{align}
$$

만약 $f$가 2변수 함수일 경우 각 변수의 target 지점이 $a$, $b$, source지점이 $x$, $y$가 되면 아래와 같다.

$$
\begin{align}
    f(a,b) 
    &= \sum_{n_1=0}^{\infty}\sum_{n_2=0}^{\infty}\frac{(a-x)^{n_1}(b-y)^{n_2}}{n_1! n_2!} \frac{\partial^{n_1+n_2} f}{\partial x^{n_1} \partial y^{n_2}} (x,y) \\
    &= f(x,y)+ \frac{1}{1!} \left ( (a-x) f_{x}(x,y) +(b-y) f_{y}(x,y) \right ) \\ &+ \frac{1}{2!} \left ( (a-x)^2 f_{xx}(x,y) + (b-y)^2 f_{yy}(x,y) + 2(a-x)(b-y) f_{xy}(x,y) \right ) \\ & + H.O.T \label{eq9} \tag{9}\\
\end{align}
$$

이런식으로 다변수 함수에 대해서도 일반화할 수 있다.

테일러 전개는 머신러닝 논문에서의 모델링 과정, 수학적인 증명 과정, Theorem 도출 과정 등 곳곳에서 다양하게 활용된다. 이때 함수에 관한 특별한 이슈가 있지 않는 이상 보통 계산의 복잡성 대비 근사치에 대한 낮은 효율로 인해 통상 1계 도함수까지만 사용해서 근사식으로 사용한다. 그러나 시스템이 확률미분방정식(Stochastic Differential Equation)인 경우에서는 2계 도함수까지 근사해야 정석이 되는데, 이는 Stochastic case의 HJB 방정식을 다루는 다음 포스팅에서 자세히 다룬다.

1변수 및 2변수 함수의 target 지점이 각각 source 에서 작은 단위로 떨어져 있는 곳이라 할 때 주로 표현되는 일반적인 함수 근사식은 그러므로 아래와 같다.

$$
f(x+\Delta x) = f(x)+ f^{'}(x) \Delta x + H.O.T \label{eq10} \tag{10}
$$

$$
f(x+\Delta x, y+\Delta y) = f(x,y) + f_x(x,y)\Delta x + f_y(x,y)\Delta y + H.O.T \label{eq11} \tag{11} 
$$

---

이제 필요한 툴들은 모두 갖추었으니 해밀턴 자코비 벨만 방정식을 유도해보자. 여기서의 사용되는 시간은 $t_0<t<t_f$이다. 식 (3)에 식 (2)를 대입하여 시작하면 다음과 같다.

$$
\begin{align}
V(\mathbf{x}(t),t)
&= \min_{\mathbf{u}[t, t_f]} \left ( \int_{t}^{t_f} \mathcal{L}(\mathbf{x}(\tau), \mathbf{u}(\tau)) d\tau + Q(\mathbf{x}(t_f), t_f)\right ) \label{eq12} \tag{12} \\
&= \min_{\mathbf{u}[t, t_f]} \left ( \int_{t}^{t+dt} \mathcal{L}(\mathbf{x}(\tau), \mathbf{u}(\tau)) d\tau + \min_{\mathbf{u}[t+dt, t_f]} \left ( \int_{t+dt}^{t_f} \mathcal{L}(\mathbf{x}(\tau), \mathbf{u}(\tau)) d\tau \right ) + Q(\mathbf{x}(t_f), t_f) \right ) \label{eq13} \tag{13}\\
&= \min_{\mathbf{u}[t, t+dt]} \left ( \int_{t}^{t+dt} \mathcal{L}(\mathbf{x}(\tau), \mathbf{u}(\tau)) d\tau + V(\mathbf{x}(t+dt),t+dt) \right ) \label{eq14} \tag{14} \\
&= \min_{\mathbf{u}[t, t+dt]} \left (\mathcal{L}(\mathbf{x}(t), \mathbf{u}(t))dt + V(\mathbf{x}(t), t)+dV(\mathbf{x}(t), t) \right ) \label{eq15} \tag{15} \\
&= \min_{\mathbf{u}[t, t+dt]} (\mathcal{L}(\mathbf{x}(t), \mathbf{u}(t)) dt + V(\mathbf{x}(t),t)+ V_x (\mathbf{x}(t), t)^T d\mathbf{x} + V_t (\mathbf{x}(t),t) dt) \label{eq16} \tag{16} \\
\therefore -V_t (\mathbf{x}(t),t)&= \min_{\mathbf{u}[t, t+dt]} (\mathcal{L}(\mathbf{x}(t), \mathbf{u}(t)) + V_x (\mathbf{x}(t), t)^T f(\mathbf{x},\mathbf{u})) \label{eq17} \tag{17}
\end{align}
$$

식 (13)은 앞서 정의한 dynamic programming에 의해 전개되었다. 식 (15)의 $\mathcal{L}$ 파트는 일반적인 리만 적분의 형식으로 직사각형 왼쪽 모서리를 사용하였고 $dV$파트는 식 (11)의 2변수 테일러 전개식의 결과로 식 (16)이 된다. 식 (17)에서는 $V(\mathbf{x},t)$를 양변에 상계시켜주고($V$는 $\mathbf{u}$를 input으로 받지 않기 때문에 min 오퍼레이터로부터 제약 받지 않는다.) $d\mathbf{x}$를 식 (1)에 따라 $f(\mathbf{x}, \mathbf{u}) dt$로 치환시킨 다음 양변을 $dt$로 나눠준 결과이다. 

이렇게 해밀턴 자코비 벨만 방정식이 유도됨을 살펴볼 수 있다.

2.3 How to use HJB Equation

해밀턴 자코비 벨만 방정식을 이용해 최적 제어를 구하는 문제는 다음과 같은 과정을 갖는다.

1. $V$에 대한 ansatz를 추론한다.
2. ansatz로 추론한 $V$로 $\frac{\partial V}{\partial t}$, $\frac{\partial V}{\partial x}$를 구하고, 주어진  system과 cost function을 이용해 HJB 방정식을 정의한다.
3. First order condition으로 $\frac{\partial H}{\partial u}=0$를 풀어 $u^{*}$ 형태를 구한다
4. 최적 제어 형태 $u^{*}$를 HJB 방정식에 대입한 뒤 편미분 방정식을 풀어 최적 제어 $u^{*}$에 붙은 파라미터를 구한다.

해밀턴 자코비 벨만 방정식을 이용해 위와 같은 프로세스를 거치면 선형은 물론 비선형 시스템에 대한 최적 제어를 구할 수 있다. 해밀턴 자코비 벨만 방정식의 특장점은 상태 방정식이 ODE 뿐 아니라 SDE(Stochastic Differential Equation)형식으로 전개되어 확률론적 최적 제어(Stochastic Optimal Control)문제로 확장되어도 그대로 적용해도 된다는 점이다. 본 컨텐츠를 상세히 포스팅한 것도 그 동기이지만 그래서 dynamics가 SDE로 모델링되는 자산 가격이 직,간접적으로 상태(state)로서 포함되는 다양한 금융 문제(ex. 포트폴리오 최적화, 옵션의 공정 가격, 최적 매수 매도 호가 지정 등등)를 다루는 논문들에서 해밀턴 자코비 벨만 방정식을 이용하여 해를 찾기도 한다. 

한편 HJB 방정식은 상태 $\mathbf x$와 제어 입력 $\mathbf u$의 차원이 커질수록 Computation cost가 기하급수적으로 늘어나고, 또한 $H$가 $\mathbf u$에 대한 convexity가 보장이 되어야만 위 프로세스의 3번처럼 First order condition을 사용하여 analytic solution을 구할 수 있다. 금융 논문에서는 그래서 보통 sum of squared 형식의 quadratic 함수나 exponential 형태 기반의 효용 함수(Utility function)를 cost function으로 사용한다. 이외에 n차 미분으로 global minimum에 해당하는 analytic solution을 구하기 어려운 문제 형태에서는 *신경망*을 이용해 numerical solution을 내기도 한다.

앞서 식 (1)의 $f$가 linear이고 식 (2)의 $J$가 quadratic form인 케이스가 LQR(Linear Quadratic Regulator) 문제가 된다고 언급했었다. 이 경우엔 First Order Condition을 이용해 Analytic solution을 구할 수 있는데, 본 포스팅의 마지막 섹션으로 아래에서 LQR 예제를 통해 해밀턴 자코비 벨만 방정식이 최적 제어 도출에 어떻게 활용되는지 살펴보도록 한다.

3. LQR example

problem)

$$
\frac{d}{dt}\mathbf{x}
=A \mathbf{x} + B \mathbf{u}
$$

$$
J(\mathbf{x}(t_0), \mathbf{u}, t_0) = \int_{t_0}^{t_f} (\mathbf{u}^T R \mathbf{u} + \mathbf{x}^T Q \mathbf{x})dt + \mathbf{x}^T(t_f) Q \mathbf{x}(t_f)
$$

$$
\mathbf{x}(t_0) = \mathbf{x}_0
$$

계산의 편의를 위해 $R$, $Q$는 symmetric matrix로 가정한다.

 

solution)

1. $V$에 대한 ansatz를 추론한다.

위의 cost function $J$의 형태에서 어떤 최적 제어가 정해지면 value function은 
$$V = \mathbf{x}^T P \mathbf{x}$$
와 같이 어떤 parameter matrix $P$를 둔 위의 형태로 ansatz를 갖게 될 것을 추론할 수 있다. 여기서 계산의 편의를 위해 $P$도 symmetric matrix로 가정한다.

2. ansatz로 추론한 $V$로 $\frac{\partial V}{\partial t}$, $\frac{\partial V}{\partial x}$를 구하고, 주어진  system과 cost function을 이용해 HJB 방정식을 정의한다.

$$\mathcal{L(\mathbf{x},\mathbf{u})}= \mathbf{u}^T R \mathbf{u} + \mathbf{x}^T Q \mathbf{x}$$
$$
f(\mathbf{x},\mathbf{u}) = A \mathbf{x}+B\mathbf{u}
$$
이므로
$$
\begin{align}
-\frac{\partial V}{\partial t}  = \mathbf{0}
&= \min_{\mathbf{u}(t)} H \\
&= \min_{\mathbf{u}(t)} \left (\frac{\partial V}{\partial x}^T f(\mathbf{x}, \mathbf{u})+\mathcal{L}(\mathbf{x}, \mathbf{u})\right ) \\
&= \min_{\mathbf{u}(t)} \left ( \underbrace{2x^TP}_{\frac{\partial V}{\partial x}} \underbrace{( A \mathbf{x} + B \mathbf{u})}_{f} + \underbrace{\mathbf{u}^T R \mathbf{u} + \mathbf{x}^T Q \mathbf{x}}_{\mathcal{L}} \right )\\
&= \min_{\mathbf{u}(t)} \left ( \mathbf{u}^T R \mathbf{u} + 2 \mathbf{x}^T P B \mathbf{u} + \mathbf{x}^T Q \mathbf{x} + 2 \mathbf{x}^T P A \mathbf{x} \right )
\end{align}
$$

3. First order condition으로 $\frac{\partial H}{\partial u}=0$를 풀어 $u^{*}$ 형태를 구한다.

$$
\begin{align}
\frac{\partial H}{\partial u}&= 0 \\
&= 2 \mathbf{u}^T R + 2 \mathbf{x}^T P B
\end{align}
$$
$$
\\
=> \mathbf{u}^T R = - \mathbf{x}^T P B \\
$$
$$
\therefore \mathbf{u}^{*T} = -\mathbf{x}^T P B R^{-1}, \qquad \mathbf{u}^{*} = -R^{-1} B^T P^T \mathbf{x}
$$

4. 최적 제어 형태 $\mathbf{u}^{*}$를 HJB 방정식에 대입한 뒤 편미분 방정식을 풀어 최적 제어에 붙은 파라미터(P)를 구한다.

3번에서 구한 최적 제어는 H의 최소값을 보장하는 닫힌 해(closed-form solution)이다. 그러므로 위 2번에서 구한 HJB 방정식의 $\mathbf{u}$에다 최적 제어 $\mathbf{u}^{*}$를 대입하면 아래와 같다.

$$
\begin{align}
\mathbf{0} 
&= - \left [ \left( -\mathbf{x}^T PBR^{-1} \right ) R \left ( -R^{-1}B^T P ^T \mathbf{x} \right ) + 2 \mathbf{x}^T PB \left( -R^{-1}B^T P^T \mathbf{x}  \right )  + \mathbf{x}^T Q \mathbf{x} + 2 \mathbf{x}^T PA \mathbf{x} \right ] \\
&= - \left[ \mathbf{x}^T PBR^{-1}B^T P ^T \mathbf{x} - 2 \mathbf{x}^T PBR^{-1}B^T P^T \mathbf{x} + \mathbf{x}^T Q \mathbf{x} + 2 \mathbf{x}^T PA \mathbf{x}\right ] \\
&= - \left[ \mathbf{x}^T \left ( PBR^{-1}B^T P ^T - 2PBR^{-1}B^T P^T+ Q + 2PA \right) \mathbf{x}\right ] \\
&= Q + PA + A^TP + PBR^{-1}B^T P ^T - 2PBR^{-1}B^T P^T
\end{align}
$$

참고로 여기서 4번의 마지막 식 $ 0 = Q + PA + A^TP + PBR^{-1}B^T P ^T - 2PBR^{-1}B^T P^T$을 algebraic riccati equation이라 한다.

algebraic riccati equation에서 $P$를 구하는 방법에 대해 여기서 더 진도를 나가는 건 LQR 문제에 한정된 지엽적인 풀이이므로 간단하게만 언급하자면, 풀이 알고리즘이 여러 개가 있지만 대표적으로 *Hamiltonian Matrix* $\mathcal{H} =\begin{bmatrix} A & -BR^{-1}B^T \\ -Q & -A^T \end{bmatrix}$ 로 놓고 Hamiltonian matrix의 Eigenvalue matrix를 만든 뒤 2개의 sub-matrix로 갈라 $U_1, U_2$를 만들고 $P=U_2 U_1^{-1}$로 구하는 방법이 있다. 그러면 결국 이 $P$를 이용하여 최적 제어에 대해 수치적으로 알 수 있게 된다.

이런 식으로 다른 문제에 적용할 때에 시스템의 상황 별로 4번까지 진행하여 나온 방정식을 풀어 최적 제어를 구해야한다.

4. Conclusion

본 포스팅에선 HJB 방정식과 그 유도, 사용법에 대해 살펴보았다. HJB 방정식은 dynamic problem에서 최적제어를 구할 때 사용되며, system이 linear일 때, nonlinear 일때, deterministic dynamics일 때, stochastic dynamics일 때 모두 사용되는 유용한 툴이다. 본 포스팅에서는 deterministic dynamics를 갖는 시스템에서 HJB 방정식을 유도했는데 금융과 같은 분야는 예를 들면 자산 가격과 같은 것을 상태 방정식으로 하여 Stochastic dynamics를 가지므로 다른 형태의 HJB 방정식이 유도된다. 이 때 Ito calculus가 사용되는데, 그것까지 설명하면 너무 포스팅이 길어져 부득이 2부까지 포스팅하게 되었다. 다음 포스팅에서는 Stochastic case의 HJB 방정식에 대해 살펴보겠다.

 

Reference

1. Bechhoefer, J. (2021). Control Theory for Physicists. Cambridge University Press.

2. https://www.youtube.com/watch?v=JCMOq0taHW8&t=832s 

3. https://www.youtube.com/watch?v=dP8yzkznnK8&t=265s