[MPC] MPC를 위한 두가지 QP 모델 - 2

MPC(model predictive control) 문제를 최적화 문제인 QP(quadratic program)문제로 변환할 때 널리 사용되는 일반적인 방법( https://pasus.tistory.com/229 )은 모델이 조밀해져서 문제의 구조가 손실되는 단점이 있다.

 

고속 MPC(fast MPC)에서는 QP 문제로 변환할 시 적절한 변수 재정렬을 사용하여 희소(sparse)구조를 최대로 이용할 수 있도록 특별한 방법으로 변환하며, 최적화 기법을 적용할 시 warm start, fixed iteration, early termination등의 휴리스틱 기법을 이용하여 MPC 계산량을 대폭 줄이는 방법을 사용한다. fast MPC알고리즘은 스탠퍼드의 Boyd 교수와 그 제자의 논문인 'Fast Model Predictive Control Using Online Optimization'에 자세히 설명되어 있다. 여기서는 fast MPC에서 사용된 QP 변환 모델을 유도해 보고자 한다.

MPC는 다음과 같은 제약조건을 갖는 선형 시스템에서

 

\[ \begin{align} & \mathbf{x}_{t+1}=A \mathbf{x}_t+B \mathbf{u}_t \tag{1} \\ \\ & \mathbf{y}_t=C \mathbf{x}_t \\ \\ & \ \ \ \ \ \mathbf{u}_{min} \le \mathbf{u}_{t+i} \le \mathbf{u}_{max}, \ \ \ i=0, ... , N-1 \tag{2} \\ \\ & \ \ \ \ \ \mathbf{y}_{min} \le \mathbf{y}_{t+i} \le \mathbf{y}_{max}, \ \ \ i=1, ... , N \end{align} \]

 

매 시간 스텝마다 다음 목적함수가 일정 성능 예측구간 \([t, \ t+N ]\) 에서 최적화되도록 제어입력의 시퀀스 \( \mathbf{u}_{t:t+N-1}^*=( \mathbf{u}_t^*, \mathbf{u}_{t+1}^*, ... , \mathbf{u}_{t+N-1}^* )\) 를 계산하는 문제다. 여기서 \(\mathbf{x}_t \in \mathbb{R}^n\), \( \mathbf{u}_t \in \mathbb{R}^p\) 이다.

 

\[ \min_{\mathbf{u}_{t:t+N-1}} J = \mathbf{x}_{t+N}^T Q_f \mathbf{x}_{t+N} + \sum_{i=0}^{N-1} \mathbf{x}_{t+i}^T Q \mathbf{x}_{t+i} + \mathbf{u}_{t+i}^T R \mathbf{u}_{t+i} \tag{3} \]

 

식 (1)에 의하면 시간 스텝(time step) \(t\) 부터 \(t+N\) 까지 미래의 상태를 다음과 같이 계산할 수 있다.

 

\[ \begin{align} & -B \mathbf{u}_t+ \mathbf{x}_{t+1}=A \mathbf{x}_t \tag{4} \\ \\ & -A \mathbf{x}_{t+1}-B \mathbf{u}_{t+1}+ \mathbf{x}_{t+2}=0 \\ \\ & -A \mathbf{x}_{t+2}-B \mathbf{u}_{t+2}+ \mathbf{x}_{t+3}=0 \\ \\ & \ \ \ \ \ \ \ \ \ \ \cdots \\ \\ & -A \mathbf{x}_{t+N-1}-B \mathbf{u}_{t+N-1}+ \mathbf{x}_{t+N}=0 \end{align} \]

 

식 (2)의 제약조건 수식은 다음과 같이 재구성할 수 있다.

 

\[ \begin{align} &\begin{bmatrix} I \\ -I \end{bmatrix} \mathbf{u}_{t+i} \le \begin{bmatrix} \mathbf{u}_{max} \\ -\mathbf{u}_{min} \end{bmatrix} , \ \ \ i=0, ... , N-1 \tag{5} \\ \\ & \begin{bmatrix} C \\ -C \end{bmatrix} \mathbf{x}_{t+i} \le \begin{bmatrix} \mathbf{y}_{max} \\ -\mathbf{y}_{min} \end{bmatrix} , \ \ \ i=1, ... , N \end{align} \]

 

식 (3)의 목적함수도 다음과 같이 행렬 형식으로 표현할 수 있다.

 

\[ \begin{align} J &= \mathbf{x}_t^T Q \mathbf{x}_t + \begin{bmatrix} \mathbf{x}_{t+1} & \mathbf{x}_{t+2} & \cdots & \mathbf{x}_{t+N-1} & \mathbf{x}_{t+N} \end{bmatrix} \begin{bmatrix} Q & & & & \\ & Q & & & \\ & & ⋱ & & \\ & & & Q & \\ & & & & Q_f \end{bmatrix} \begin{bmatrix} \mathbf{x}_{t+1} \\ \mathbf{x}_{t+2} \\ \mathbf{x}_{t+3} \\ \vdots \\ \mathbf{x}_{t+N} \end{bmatrix} \tag{6} \\ \\ & \ \ \ \ \ \ \ \ + \begin{bmatrix} \mathbf{u}_t & \mathbf{u}_{t+1} & \cdots & \mathbf{u}_{t+N-2} & \mathbf{u}_{t+N-1} \end{bmatrix} \begin{bmatrix} R & & & & \\ & R & & & \\ & & ⋱ & & \\ & & & R & \\ & & & & R \end{bmatrix} \begin{bmatrix} \mathbf{u}_t \\ \mathbf{u}_{t+1} \\ \mathbf{u}_{t+2} \\ \vdots \\ \mathbf{u}_{t+N-1} \end{bmatrix} \end{align} \]

 

이제 최적화 변수 \(\mathbf{z}\) 를 다음과 같이 정의한다.

 

\[ \mathbf{z}= \begin{bmatrix} \mathbf{u}_t \\ \mathbf{x}_{t+1} \\ \mathbf{u}_{t+1} \\ \vdots \\\mathbf{u}_{t+N-1} \\ \mathbf{x}_{t+N} \end{bmatrix} \ \ \in \mathbb{R}^{N(n+p)} \tag{7} \]

 

그러면 식 (6)의 목적함수를 다음과 같이 최적화 변수의 2차함수로 바꿀 수 있다.

 

\[ \tilde{J}= \mathbf{z}^T H \mathbf{z} \tag{8} \]

 

여기서 \(\mathbf{x}_t^T Q \mathbf{x}_t\) 는 최적화 대상이 아니기 때문에 제외하였으며,

 

\[ H= \begin{bmatrix} R & & & & & & \\ & Q & & & & & \\ & & R & & & & \\ & & & ⋱ & & & \\ & & & & Q & & \\ & & & & & R & \\ & & & & & & Q_f \end{bmatrix} \ \ \in \mathbb{R}^{N(n+p) \times N(n+p) } \]

 

이다. 식 (5)의 제약조건도 다음과 같이 최적화 변수의 부등식으로 바꿀 수 있다.

 

\[ P \mathbf{z} \le \mathbf{h} \tag{9} \]

 

여기서

 

\[ P= \begin{bmatrix} I & & & \\ & C & & \\ & & ⋱ & \\ & & & C \\-I & & & \\ & -C & & \\ & & ⋱ & \\ & & & -C \end{bmatrix} \ \ \in \mathbb{R}^{2N(n+p) \times N(n+p)} , \ \ \ \mathbf{h}= \begin{bmatrix} \mathbf{u}_{max} \\ \mathbf{x}_{max} \\ \vdots \\ \mathbf{x}_{max} \\ -\mathbf{u}_{min} \\ -\mathbf{x}_{min} \\ \vdots \\ -\mathbf{x}_{min} \end{bmatrix} \]

 

이다.

 

 

시스템 방정식 (4)는 다음과 같이 최적화 변수의 등식으로 바꿀 수 있다.

 

\[ D \mathbf{z} = \mathbf{b} \tag{10} \]

 

여기서

 

\[ \begin{align} & D= \begin{bmatrix} -B & I & 0 & 0 & 0 & 0 & & & & 0 \\ 0 & -A & -B & I & 0 & 0 & & & & \\ & & & -A & -B & I & & & & & \\ & & & & & & ⋱ & & & \\ & & & & & & & I & 0 & 0 \\ & & & & & & & -A & -B & I \end{bmatrix} \ \ \in \mathbb{R}^{Nn \times N(n+p) } , \\ \\ & \mathbf{b}= \begin{bmatrix} A\mathbf{x}_t \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0\end{bmatrix} \ \ \in \mathbb{R}^{Nn} \end{align} \]

 

이다. 식 (8)-(10)에 의하면 식 (1)-(3)의 MPC 최적화 문제는 다음과 같이 QP 최적화 문제로 변환된다.

 

\[ \begin{align} & \min_{\mathbf{z}} \tilde{J} = \mathbf{z}^T H \mathbf{z} \tag{11} \\ \\ & \mbox{subject to } \ \ P \mathbf{z} \le \mathbf{h}, \ \ \ \ D\mathbf{z}= \mathbf{b} \end{align} \]

 

위 최적화 문제의 해는 \(\mathbf{x}_t\) 값에 대한 \(\mathbf{z}^*=(\mathbf{u}_t^*, \mathbf{x}_{t+1}^*, \mathbf{u}_{t+1}^* , ... , \mathbf{x}_{t+N}^*)\) 으로서 제어입력과 상태변수의 최적 시퀀스이다.

 

 

식 (11)의 QP 문제의 행렬과 벡터의 구조를 보면 널리 사용되는 일반적인 방법과는 달리 MPC 모델의 구조가 훼손되지 않았고 희소(sparse)행렬임을 알 수 있다. 따라서 범용 QP 솔버 대신 이를 최대로 이용할 수 있는 특화된 QP 솔버를 만들 수 있다.