[PSOC-11] 가우스 유사 스펙트럴 (GPM) 기반 최적제어

가우스 유사 스펙트럴 방법(GPM, Gauss pseudospectral method)에서는 \(N\) 개의 LG(Legendre-Gauss) 포인트를 콜로케이션 포인트로 사용하고, LG 포인트에 \(\tau_0=-1\) 을 포함한 점을 보간점으로 사용한다. 이산화 점은 보간점에 \(\tau_{N+1}=1\) 을 포함한 것이다. 따라서 가우스 유사 스펙트럴 방법은 \(N\) 개의 콜로케이션 포인트, \(N+1\) 개의 보간점와 \(N+2\) 개의 이산화 점을 사용한다.

 

 

LG 포인트는 \(N\) 차 르장드르 다항식 \(P_N (\tau)\) 의 해로 구성되어 있다. 가우스 유사 스펙트럴 방법에서는 상태변수 \(\mathbf{x}(\tau)\) 를 \( N\) 차 라그랑지 다항식으로 근사화한다.

 

\[ \mathbf{x}( \tau ) \approx \mathbf{X} (\tau)= \sum_{i=0}^N \mathbf{X}_i L_i^{(N)} (\tau) \tag{1} \]

 

여기서 \(\mathbf{X}_i= \mathbf{X}(\tau_i ) \in \mathbb{R}^n\) 이고 \(L_i^{(N)} (\tau)\) 는 \(N\) 차 라그랑지 다항식으로 다음과 같다.

 

\[ L_i^{(N)} (\tau)= \prod_{ \begin{matrix} k=0 \\ k \ne i \end{matrix} }^N \frac{\tau-\tau_k}{\tau_i-\tau_k } \tag{2} \]

 

 

 

\(\tau_0=-1\) 는 LG 포인트가 아님에 주의해야 한다. \(N\) 개의 LG 콜로케이션 포인트에서 상태변수 \(\mathbf{x}(\tau)\) 의 근사해는 다음과 같이 운동 미분방정식을 정확히 만족시켜야 한다.

 

\[ \begin{align} \frac{d \mathbf{x}( \tau)}{d \tau} & = \frac{(t_f-t_0 )}{2} \mathbf{f}(\mathbf{x}(\tau), \mathbf{u}(\tau), \tau, t_0, t_f ) \\ \\ \to \dot{\mathbf{X}}(\tau_j ) &= \sum_{i=0}^N \mathbf{X}_i \dot{L}_i^{(N)} (\tau_j ) \tag{3} \\ \\ & = \frac{(t_f-t_0)}{2} \mathbf{f}( \mathbf{X}_j, \mathbf{U}_j, \tau_j, t_0, t_f ), \ \ \ \ j=1, ..., N \end{align} \]

 

여기서 \(\mathbf{U}_i \in \mathbb{R}^p\) 이다. \(\dot{L}_i^{(N)} (\tau_j )=D_{ji}\) 로 치환하고 위 식을 행렬식으로 바꾸면 다음과 같다.

 

\[ \begin{bmatrix} D_{10} & D_{11} & \cdots & D_{1N} \\ \vdots & \vdots & \ddots & \vdots \\ D_{N0} & D_{N1} & \cdots & D_{NN} \end{bmatrix} \begin{bmatrix} \mathbf{X}_0^T \\ \mathbf{X}_1^T \\ \vdots \\ \mathbf{X}_N^T \end{bmatrix} = \frac{(t_f-t_0 )}{2} \begin{bmatrix} \mathbf{f}^T ( \mathbf{X}_1, \mathbf{U}_1, \tau_1, t_0, t_f ) \\ \vdots \\ \mathbf{f}^T (\mathbf{X}_N, \mathbf{U}_N, \tau_N, t_0, t_f) \end{bmatrix} \tag{4} \]

 

여기서 \([D_{ji} ] \in \mathbb{R}^{N \times (N+1) }\) 이다. 목적함수는 다음과 같이 \(N\) 개의 LG 포인트를 이용한 가우시안 쿼드래처로 근사화 된다.

 

\[ \begin{align} J & = \phi (\mathbf{x}(-1), \mathbf{x}(1), t_0, t_f )+ \frac{(t_f-t_0 )}{2} \int_{-1}^1 g( \mathbf{x}(\tau), \mathbf{u}( \tau), \tau, t_0, t_f) \ d \tau \\ \\ \to J &= \phi( \mathbf{X}_1, \mathbf{X}_{N+1}, t_0, t_f )+ \frac{(t_f-t_0 )}{2} \sum_{j=1}^N w_j \ g( \mathbf{X}_j, \mathbf{U}_j, \tau_j, t_0, t_f ) \tag{5} \end{align} \]

 

여기서 \(w_j\) 는 LG 쿼드래처 포인트의 가중치다. 최종 상태변수 값인 \(\mathbf{X}_{N+1}\) 은 보간점에 속한 변수가 아니기 때문에 따로 계산해야 한다. 우선 시스템 방정식을 \(\tau_{N+1}\) 까지 적분하면 다음과 같다.

 

\[ \mathbf{x}(\tau_{N+1} )=\mathbf{x}(\tau_0 )+\frac{(t_f-t_0 )}{2} \int_{-1}^1 \mathbf{f}( \mathbf{x}(\tau), \mathbf{u}(\tau), \tau, t_0, t_f ) \ d\tau \tag{6} \]

 

여기서 적분항을 \(N\) 개의 LG 포인트를 이용하여 가우시안 쿼드래처로 근사화하면 다음과 같이 최종 상태변수 값 \(\mathbf{X}_{N+1}\) 을 구할 수 있다.

 

\[ \mathbf{X}_{N+1} = \mathbf{X}_0 + \frac{(t_f-t_0 )}{2} \sum_{j=1}^N w_j \mathbf{f}( \mathbf{X}_j, \mathbf{U}_j, \tau_j, t_0, t_f ) \tag{7} \]

 

경계조건은 다음과 같다.

 

\[ \begin{align} & \psi( \mathbf{x}(-1), \mathbf{x}(1), t_0, t_f )=0 \\ \\ \to \ & \psi (\mathbf{X}_1, \mathbf{X}_{N+1}, t_0, t_f )=0 \tag{8} \end{align} \]

 

경로제약 조건은 \(N\) 개의 LG 포인트에서 다음 식을 만족해야 한다.

 

\[ \begin{align} & \frac{(t_f-t_0 )}{2} \mathbf{C}( \mathbf{x} (\tau), \mathbf{u}(\tau), \tau, t_0, t_f ) \le 0 \\ \\ \to \ & \frac{(t_f-t_0)}{2} \mathbf{C}( \mathbf{X}_j, \mathbf{U}_j, \tau_j, t_0, t_f ) \le 0 , \ \ \ j=1, ..., N \tag{9} \end{align} \]

 

가우스 유사 스펙트럴에서 유도된 비선형 프로그래밍 문제(NLP)는 식 (4), (7), (8), (9) 의 대수적 제약 조건에 따라 식 (5)의 목적함수를 최소화하는 것이다. 이 때 최적화 대상 변수는 다음과 같다.

 

\[ (\mathbf{X}_0, \ \mathbf{X}_1, \ ... , \ \mathbf{X}_N ), \ ( \mathbf{U}_1, \ \mathbf{U}_2, \ ... , \ \mathbf{U}_N ), \ t_0, \ t_f \tag{8} \]

 

여기서 \(t_0, \ t_f\) 는 최적제어 문제에 따라서 변수가 될 수도 아닐 수도 있다. 또한 초기 제어값 \(\mathbf{U}_0\) 와 최종 제어값 \(\mathbf{U}_{N+1}\) 은 최적화 대상변수가 아니기 때문에 계산되지 않는다.

콜로케이션 포인트에서 제어입력 값 \(\mathbf{U}_i\) 를 계산한 후 다음과 같이 \((N-1)\) 차 라그랑지 보간법을 이용하여 제어입력의 시간 함수를 구할 수 있다.

 

\[ \mathbf{u}(\tau) \approx \mathbf{U}(\tau)= \sum_{i=1}^N \mathbf{U}_i L_i^{(N-1)} (\tau) \tag{11} \]

 

여기서

 

\[ L_i^{(N-1)} (\tau)= \prod_{ \begin{matrix} k=1 \\ k \ne i \end{matrix} }^N \frac{\tau-\tau_k}{\tau_i-\tau_k } \tag{12} \]

 

이다.

 

 

이제 앞선 포스팅(https://pasus.tistory.com/281)에 있는 최적제어 예제 (a)를 GPM을 이용하여 풀어보자. 최적제어 문제를 정적 파라미터 최적화 문제로 변환되는 과정을 자세히 살펴보기 위해서 콜로케이션 포인트의 숫자를 작게 잡도록 한다. 여기서는 3개의 LG 콜로케이션 포인트를 사용한다. LG 콜로케이션 포인트는 다음과 같다.

 

N = 3; % number of collocation points
points = GaussPoints(N,'LG'); % LG collocation points (N)
col_point = points' 

col_point = 
   -0.7746    0    0.7746

 

보간점은 콜로케이션 포인트에 \(-1\) 점을 추가한 것이다.

 

int_point = [-1 col_point] % interpolating points (N+1)

int_point = 
   -1.0000   -0.7746    0    0.7746

 

식 (1)에 의하면 상태변수는 보간점을 이용하여 다음과 같이 이산화된다.

 

\[ \begin{align} & x_1 (\tau) \approx X_1 (\tau)= \sum_{i=0}^3 X_{i1} L_i^{(3)} (\tau) \\ \\ & x_2 (\tau) \approx X_2 (\tau)= \sum_{i=0}^3 X_{i2} L_i^{(3)} (\tau) \end{align} \]

 

3차 라그랑지 다항식은 다음과 같다.

 

\[ L_i^{(3)} (\tau)= \prod_{ \begin{matrix} k=0 \\ k \ne i \end{matrix} }^3 \frac{\tau-\tau_k}{\tau_i-\tau_k } \]

 

콜로케이션 포인트에서 미분방정식을 만족해야 하므로 식 (3)에 의하면 다음 식을 유도할 수 있다.

 

\[ \begin{align} & \dot{x}_1 (\tau)=x_2 (\tau) \\ \\ & \dot{x}_2(\tau)=u(\tau) \\ \\ \to \ & \dot{X}_1 (\tau_j )= \sum_{i=0}^3 X_{i1} \dot{L}_i^{(3)} (\tau_j )=X_{j2} \\ \\ & \dot{X}_2 (\tau_j )= \sum_{i=0}^3 X_{i2} \dot{L}_i^{(3)} (\tau_j ) =U_j, \ \ \ \ \ j=1, 2, 3 \end{align} \]

 

식 (4)에 의하면 위 식의 행렬 형식은 다음과 같다.

 

\[ \begin{bmatrix} D_{10} & D_{11} & D_{12} & D_{13} \\ D_{20} & D_{21} & D_{22} & D_{23} \\ D_{30} & D_{31} & D_{32} & D_{33} \end{bmatrix} \begin{bmatrix} X_{01} & X_{02} \\ X_{11} & X_{12} \\ X_{21} & X_{22} \\ X_{31} & X_{32} \end{bmatrix} = \begin{bmatrix} X_{12} & U_1 \\ X_{22} & U_2 \\ X_{32} & U_3 \end{bmatrix} \]

 

 

Dji = LagrangePolyDerivative(col_point,int_point) % Differentiation matrix

Dji = 
   -3.0000    2.5000    0.5820   -0.0820
    1.5000   -2.8637    1.0000    0.3637
   -3.0000    5.0820   -4.5820    2.5000
    

 

목적함수는 콜로케이션 포인트를 이용한 가우시안 쿼드래처로 근사화 시킨다.

 

\[ \begin{align} & J_1= \frac{1}{2} \int_{-1}^1 u^2 (\tau) \ d \tau \\ \\ \to \ & J_1= \frac{1}{2} \sum_{j=1}^3 w_j U_j^2 \end{align} \]

 

 

gweights = GaussianQuadrature(N,'LG') % LG quadrature weights

gweights = 0.5556  0.8889  0.5556
    

 

초기조건과 최종조건은 다음과 같다.

 

\[ \begin{align} & x_1 (-1)=2, \ x_2 (-1)=1, \ x_1 (1)=0, \ x_2 (1)=0 \\ \\ \to \ & X_{01}=2, \ X_{02}=1, \ X_{41}=0,\ X_{42}=0 \end{align} \]

 

여기서 \(X_{41}\) 과 \(X_{42}\) 는 식 (7)로 계산해야 한다.

 

\[ \begin{align} X_{41} &= X_{01}+ \sum_{j=1}^3 w_j X_{j2} \\ \\ X_{42} &= X_{02}+ \sum_{j=1}^3 w_j U_j \end{align} \]

 

위 식을 종합하면 다음과 같은 정적 파라미터 최적화 문제를 구성할 수 있다.

 

\[ \begin{align} & \min J_1 = \frac{1}{2} \sum_{j=1}^3 w_j U_j^2 \\ \\ & \mbox{subject to }: \\ \\ & \ \ \ \ \ \begin{bmatrix} D_{10} & D_{11} & D_{12} & D_{13} \\ D_{20} & D_{21} & D_{22} & D_{23} \\ D_{30} & D_{31} & D_{32} & D_{33} \end{bmatrix} \begin{bmatrix} X_{01} & X_{02} \\ X_{11} & X_{12} \\ X_{21} & X_{22} \\ X_{31} & X_{32} \end{bmatrix} = \begin{bmatrix} X_{12} & U_1 \\ X_{22} & U_2 \\ X_{32} & U_3 \end{bmatrix} \\ \\ & \ \ \ \ \ X_{11}=2, \ X_{12}=1, \\ \\ & \ \ \ \ \ 0 = X_{01}+ \sum_{j=1}^3 w_j X_{j2}, \ \ \ \ \ 0= X_{02}+ \sum_{j=1}^3 w_j U_j \end{align} \]

 

매트랩의 비선형 최적화 문제 솔버인 fmincon 함수를 이용하여 위 문제를 풀면 다음과 같다.

 

\[ \begin{align} & \begin{bmatrix} X_{01} & X_{02} \\ X_{11} & X_{12} \\ X_{21} & X_{22} \\ X_{31} & X_{32} \end{bmatrix} = \begin{bmatrix} 2 & 1 \\ 2.1070 & -0.0127 \\ 1.2500 & -1.7500 \\ 0.0930 & -0.7873 \end{bmatrix} \\ \\ & \begin{bmatrix} U_1 \\ U_2 \\ U_3 \end{bmatrix} = \begin{bmatrix} -3.9857 \\ -0.5000 \\ 2.9857 \end{bmatrix} \end{align} \]

 

아래 그림은 라그랑지 보간 다항식을 이용하여 상태변수와 제어입력의 반응을 그린 것이다. 해석적으로 구한 상태변수와 제어입력의 반응과 일치함을 알 수 있다.