Linear Quadratic Regulator

Jeongsu Ahn·2025년 7월 25일

ControlBootcamp

목록 보기

6/7

Inverted Pendulum on Cart LQR Control

Reference: Prof Steve Brunton, Control Bootcamp
Youtube Link

위와 같은 카트에 달린 역진자 시스템을 LQR로 제어하는 예제를 실습해보자.

상태 피드백 제어

일반적으로 모든 상태를 피드백 받을 수 있다고 가정한다면, pole placement와 같은 제어기는 다음과 같이 정의된다.

$u = -Kx$

$\dot{x} = (A-BK)x$

그러나 이 제어 게인 $K$ 를 구하기 위해 극점을 설계하는 과정은 많은 노력이 필요하다. 이를 해결하기 위해 최적화의 개념을 도입하여 $K$ 를 도출하는 것이 LQR이다.

LQR 제어기 설계

최소화해야 하는 목적함수를 다음과 같이 정의한다:

$J = \int^\infty_0(x^TQx + u^TRu)dt$

무한대의 시간 동안 상태 오차를 최소화하며, 제어입력을 최소한으로 사용하는 것이 목적이다.

가중치 행렬 설정:

$Q = \text{diag}(1, 1, 10, 100)$

$R = 0.001$

여기서 $Q$ 는 대각행렬로, 각 상태변수에 대한 가중치를 나타낸다:

$Q_{11} = 1$ : 위치 $x$ 에 대한 가중치
$Q_{22} = 1$ : 속도 $v$ 에 대한 가중치
$Q_{33} = 10$ : 각도 $\theta$ 에 대한 가중치 (높은 우선순위)
$Q_{44} = 100$ : 각속도 $\omega$ 에 대한 가중치 (가장 높은 우선순위)

$R$ 은 제어입력에 대한 가중치로, 작은 값은 제어력 사용에 대한 제약이 적음을 의미한다.

시스템 모델링

역진자-카트 시스템의 선형화된 상태공간 모델은 다음과 같다:

$\dot{x} = Ax + Bu$

여기서 상태벡터는 $x = [x, \dot{x}, \theta, \dot{\theta}]^T$ 이며:

$x$ : 카트의 위치
$\dot{x}$ : 카트의 속도
$\theta$ : 진자의 각도 (수직으로부터의 편차)
$\dot{\theta}$ : 진자의 각속도

시스템 매개변수를 다음과 같이 설정한다:

$m = 1$ kg (진자의 질량)
$M = 5$ kg (카트의 질량)
$L = 2$ m (진자의 길이)
$g = 10$ m/s² (중력가속도)
$d = 1$ N·s/m (댐핑 계수)

이를 통해 시스템 행렬을 구성하면:

A = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & -\frac{d}{M} & \frac{mg}{M} & 0 \\ 0 & 0 & 0 & 1 \\ 0 & -\frac{d}{ML} & \frac{(m+M)g}{ML} & 0 \end{bmatrix}

B = \begin{bmatrix} 0 \\ \frac{1}{M} \\ 0 \\ \frac{1}{ML} \end{bmatrix}

Hamilton-Jacobi-Bellman (HJB) 방정식 이론

HJB 방정식의 개념

Hamilton-Jacobi-Bellman 방정식은 동적 계획법(Dynamic Programming)의 연속시간 버전으로, 최적 제어 문제를 해결하는 핵심 도구이다. 이는 Richard Bellman의 최적성 원리(Principle of Optimality)에 기반한다.

최적성 원리: "최적 정책의 어떤 부분도 최적이어야 한다"

즉, 시간 $t$ 에서 상태 $x(t)$ 에 있을 때, 앞으로의 최적 제어 전략은 현재 상태에서만 의존해야 하며, 과거 경로와는 무관해야 한다.

가치함수 (Value Function) 정의

현재 상태 $x$ 에서 시작하여 무한대까지의 최소 누적 비용을 나타내는 가치함수를 정의한다:

$V(x) = \min_{u(\cdot)} \int_0^{\infty} L(x(t), u(t)) dt$

주의: 여기서 $x(t)$ 는 초기조건 $x(0) = x$ 에서 시작하여 제어 정책 $u(\cdot)$ 에 의해 결정되는 시간에 따른 상태 궤적이다. 즉, $V(x)$ 는 초기 상태 $x$ 만의 함수이지만, 적분 내부의 $x(t)$ 는 시간에 따라 변하는 상태를 의미한다.

여기서 $L(x, u) = x^T Q x + u^T R u$ 는 순간 비용함수이다.

HJB 방정식의 물리적 의미

HJB 방정식은 다음과 같은 동적 일관성 조건을 수학적으로 표현한다:

"현재 시점에서의 최적 비용은 다음 순간의 최적 비용과 현재 순간의 비용의 합과 같아야 한다"

HJB 방정식 미분형태의 도출 과정

동적 일관성 조건 (Dynamic Consistency Condition)

HJB 방정식의 핵심 아이디어는 동적 일관성이다. 시간 $t$ 에서 상태 $x(t)$ 에 있을 때, 가치함수 $V(x(t))$ 는 다음 관계를 만족해야 한다:

V(x(t)) = \min_u \left[ \int_t^{t+dt} L(x(\tau), u(\tau)) d\tau + V(x(t+dt)) \right]

이는 "현재의 최적 비용 = 미소 시간 동안의 비용 + 다음 시점의 최적 비용"을 의미한다.

Taylor 전개를 통한 미분형태 도출

Step 1: 미소 시간 $dt$ 동안 상태 변화

시스템 동역학 $\dot{x} = f(x, u) = Ax + Bu$ 에 의해:

x(t+dt) = x(t) + f(x(t), u(t)) \cdot dt + O(dt^2)

여기서 $O$ 는 Big O 표기로 2차 이상의 High order term을 의미한다.

Step 2: 가치함수의 Taylor 전개

V(x(t+dt)) = V(x(t)) + \nabla V(x(t))^T [x(t+dt) - x(t)] + O(dt^2)

= V(x(t)) + \nabla V(x(t))^T f(x(t), u(t)) \cdot dt + O(dt^2)

Step 3: 적분 근사

미소 시간 동안의 비용:

\int_t^{t+dt} L(x(\tau), u(\tau)) d\tau \approx L(x(t), u(t)) \cdot dt

동적 일관성 조건 대입

위의 근사들을 동적 일관성 조건에 대입하면:

V(x(t)) = \min_u \left[ L(x(t), u(t)) \cdot dt + V(x(t)) + \nabla V(x(t))^T f(x(t), u(t)) \cdot dt \right]

양변에서 $V(x(t))$ 를 빼고 $dt$ 로 나누면:

0 = \min_u \left[ L(x(t), u(t)) + \nabla V(x(t))^T f(x(t), u(t)) \right]

HJB 방정식의 최종 미분형태

시간 의존성을 생략하여 일반적인 형태로 쓰면:

\boxed{0 = \min_u \left[ L(x, u) + \nabla V^T f(x, u) \right]}

여기서 $f(x, u) = Ax + Bu$ 는 시스템 동역학이다.

해밀토니안과 물리적 해석

HJB 방정식에서 대괄호 안의 식을 해밀토니안이라고 한다:

H(x, u, \nabla V) = L(x, u) + \nabla V^T f(x, u)

따라서 HJB 방정식은:

0 = \min_u H(x, u, \nabla V)

각 항의 의미:

$L(x, u) = x^T Q x + u^T R u$ :
- 현재 순간의 즉각적 비용
- 상태 오차와 제어 노력의 가중합
$\nabla V^T f(x, u)$ :
- 가치함수의 시간 변화율
- $\nabla V$ 는 "비용이 증가하는 방향"을 나타내는 그래디언트
- $f(x, u)$ 는 상태의 변화 방향
- 내적은 "비용 변화의 예상치"를 의미한다

HJB 방정식을 통한 LQR 제어 게인 도출

Step 1: HJB 방정식 구성

LQR 문제에서 HJB 방정식:

0 = \min_u \left[ x^T Q x + u^T R u + \nabla V^T (Ax + Bu) \right]

여기서 $\nabla V = \frac{\partial V}{\partial x}$ 는 가치함수의 그래디언트이다.

Step 2: 최적 제어 조건 (First-Order Optimality Condition)

제어 입력 $u$ 에 대해 최소화하기 위한 필요조건:

\frac{\partial}{\partial u} \left[ u^T R u + \nabla V^T B u \right] = 0

이를 계산하면:

2Ru + B^T \nabla V = 0

따라서 최적 제어:

u^*(x) = -\frac{1}{2}R^{-1} B^T \nabla V

Step 3: 가치함수의 이차형태 가정

선형-이차 문제의 특성상, 가치함수는 이차형태(Quadratic Form)를 가진다:

V(x) = x^T P x

여기서 $P$ 는 양정치 대칭행렬이다.

따라서 그래디언트:

\nabla V = \frac{\partial}{\partial x}(x^T P x) = 2Px

Step 4: LQR 제어 게인 도출

최적 제어를 다시 쓰면:

u^*(x) = -\frac{1}{2}R^{-1} B^T (2Px) = -R^{-1} B^T P x

상태 피드백 형태 $u = -Kx$ 와 비교하여:

\boxed{K = R^{-1} B^T P}

이것이 LQR 제어 게인이다.

Step 5: 대수 리카티 방정식 (Algebraic Riccati Equation) 도출

최적 제어 $u^*$ 를 HJB 방정식에 대입:

0 = x^T Q x + (Px)^T B R^{-1} B^T (Px) + (2Px)^T A x - (Px)^T B R^{-1} B^T (Px)

정리하면:

0 = x^T Q x + 2x^T P^T A x - x^T P^T B R^{-1} B^T P x

$P$ 가 대칭행렬이므로 $P^T = P$ 이고, $x$ 는 임의이므로:

\boxed{A^T P + P A - P B R^{-1} B^T P + Q = 0}

이것이 대수 리카티 방정식(ARE)이다.

리카티 방정식의 해석

물리적 의미

리카티 방정식의 각 항목은 다음과 같은 의미를 가진다:

$Q$ : 상태 비용 - 시스템이 목표 상태에서 벗어나는 정도
$A^T P + P A$ : 개루프 시스템의 안정성 기여도
$-P B R^{-1} B^T P$ : 피드백 제어의 안정화 효과

해의 존재성과 유일성

ARE의 해 $P$ 가 존재하고 유일하기 위한 조건:

제어가능성: $(A, B)$ 쌍이 제어 가능해야 함
관측가능성: $(A, Q^{1/2})$ 쌍이 관측 가능해야 함
가중치 조건: $Q \geq 0$ , $R > 0$

이 조건들이 만족되면, ARE는 유일한 양정치 해 $P > 0$ 을 가진다.

LQR에서의 특수성

이차 가치함수: LQR 문제에서는 $V(x) = x^T P x$ (이차형태)이므로:

\nabla V = \frac{\partial}{\partial x}(x^T P x) = 2Px

이는 선형 그래디언트를 가지므로 HJB 방정식이 해석적으로 풀린다.

명시적 해: 이차 가치함수 덕분에 HJB 방정식이 대수 방정식 (리카티 방정식)으로 변환되어, 수치적 또는 해석적으로 쉽게 풀 수 있다.

이것이 LQR이 최적 제어 이론의 입문으로 자주 사용되는 이유이다 - 일반적인 HJB 방정식은 편미분방정식이지만, LQR에서는 대수방정식이 된다.

MATLAB 구현

시스템 매개변수 및 행렬 정의

% System parameters
m = 1;    % pendulum mass (kg)
M = 5;    % cart mass (kg)  
L = 2;    % pendulum length (m)
g = 10;   % gravity (m/s^2)
d = 1;    % damping coefficient (N*s/m)

% State-space matrices
A = [0      1              0           0;
     0    -d/M           m*g/M         0;
     0      0              0           1;
     0   -d/(M*L)  (m+M)*g/(M*L)      0];

B = [0; 1/M; 0; 1/(M*L)];

% Weight matrices
Q = diag([1, 1, 10, 100]);  % State weights
R = 0.001;                  % Control weight

LQR 제어 게인 계산

% Solve Algebraic Riccati Equation
P = care(A, B, Q, R);

% Calculate LQR gain
K = R\(B'*P);  % Equivalent to inv(R)*B'*P

% Display results
fprintf('LQR Control Gain K:\n');
disp(K);

% Closed-loop eigenvalues
A_cl = A - B*K;
eig_cl = eig(A_cl);
fprintf('Closed-loop eigenvalues:\n');
disp(eig_cl);

시뮬레이션 및 시각화

% Initial conditions
x0 = [0; 0; 0.2; 0];  % Initial angle: 0.2 rad (≈ 11.5°)

% Time vector
t = 0:0.01:10;

% Closed-loop system simulation
sys_cl = ss(A_cl, B, eye(4), 0);
[y, t] = initial(sys_cl, x0, t);

% Extract states
x_pos = y(:,1);     % Cart position
x_vel = y(:,2);     % Cart velocity  
theta = y(:,3);     % Pendulum angle
theta_dot = y(:,4); % Pendulum angular velocity

% Calculate control input
u = -K * y';

% Plotting
figure('Position', [100, 100, 1200, 800]);

% State responses
subplot(2,3,1);
plot(t, x_pos, 'LineWidth', 2);
title('Cart Position', 'FontSize', 12);
xlabel('Time (s)', 'FontSize', 10);
ylabel('Position (m)', 'FontSize', 10);
grid on;

subplot(2,3,2);
plot(t, x_vel, 'LineWidth', 2);
title('Cart Velocity', 'FontSize', 12);
xlabel('Time (s)', 'FontSize', 10);
ylabel('Velocity (m/s)', 'FontSize', 10);
grid on;

subplot(2,3,3);
plot(t, theta*180/pi, 'LineWidth', 2);
title('Pendulum Angle', 'FontSize', 12);
xlabel('Time (s)', 'FontSize', 10);
ylabel('Angle (deg)', 'FontSize', 10);
grid on;

subplot(2,3,4);
plot(t, theta_dot, 'LineWidth', 2);
title('Pendulum Angular Velocity', 'FontSize', 12);
xlabel('Time (s)', 'FontSize', 10);
ylabel('Angular Velocity (rad/s)', 'FontSize', 10);
grid on;

subplot(2,3,5);
plot(t, u, 'LineWidth', 2);
title('Control Input', 'FontSize', 12);
xlabel('Time (s)', 'FontSize', 10);
ylabel('Force (N)', 'FontSize', 10);
grid on;

% Phase portrait (angle vs angular velocity)
subplot(2,3,6);
plot(theta*180/pi, theta_dot, 'LineWidth', 2);
title('Phase Portrait (θ vs θ̇)', 'FontSize', 12);
xlabel('Angle (deg)', 'FontSize', 10);
ylabel('Angular Velocity (rad/s)', 'FontSize', 10);
grid on;

sgtitle('LQR Control of Inverted Pendulum', 'FontSize', 14, 'FontWeight', 'bold');

Jeongsu Ahn

될 때까지 하는 사람

이전 포스트

Inverted Pendulum on a Cart

다음 포스트