최적설계 | 제약조건 문제 수치해법(3): 순차 이차계획법 SQP

1. 순차 이차계획법

   앞서 다룬 순차 선형계획법은 일반적인 제약조건 최적화 문제를 풀이를 위한 단순하고 직관적인 알고리즘이지만,  상황에 따라 정확한 최적해에 수렴하지 않는 등 강건성이 부족하다는 단점이 있다. 반면에 순차 이차계획법(Sequantial Quadratic Programming, SQP)은 강건성을 갖춘 결과를 보여 비교적 널리 사용되고 있다. 해당 알고리즘은 ▲탐색방향을 결정하는 문제와 ▲이동거리를 결정하는 문제를 포함함다. 

2. 탐색방향 결정

   순차 이차계획법에서 탐색방향을 결정하는 문제는 비선형 목적함수와 제약조건 함수에 대해서 선형화된 근사를 사용한다. 아래와 같이 선형화된 목적함수에 탐색방향 이계항을 더하여 목적함수와 탐색방향의 길이를 최소화하는 이차계획문제로 만든다. 해당 이차계획문제는 복록하므로 최소값이 존재한다면 전역적으로 유일하다. 아래 링크로 접속하면 비선형·제약조건 최적설계문제를 선형화하는 방법에 대해 확인할 수 있다. 
 

$$ \begin{align} \mathrm{minimize}~~~~&\overline{f} = \mathbf{c}^T \mathbf{d} + \frac{1}{2}\mathbf{d}^T \mathbf{d} \\\\ \mathrm{subject~to}~~~~&\mathbf{N}^T\mathbf{d} = \mathbf{e} \\\\ &\mathbf{A}^T\mathbf{d} \leq \mathbf{b} \end{align}$$

 

최적설계 | 제약조건 문제 수치해법(1): 설계문제 선형화

 

3. 강하함수 계산

   순차 이차계획법에서는 이동거리를 결정하기 위해 목적함수 대신 강하함수를 정의하여 이를 최소화한다. 탐색방향과 설계점이 결정된 상태에서의 강하함수는 이동거리에 대한 일변수 함수이므로 구간감소법을 이용해 이동거리를 결정한다. 비제약조건 문제에서는 수렴 상태를 확인하기 위해 목적함수가 강하함수로서 사용되었다. 하지만 제약조건 문제에서는 목적함수에 제약조건 위배에 대한 벌칙을 반영한 강하함수를 정의하여 이를 기반으로 수렴 상태를 확인한다. 쉐니크니 강하함수(Pshenichny's descent function)는 최적화 공학 문제에 일반적으로 사용되는 함수이며, 임의의 점에서 쉐니크니 강하함수는 아래와 같이 정의된다. f(x_k)는 설계점 x_k에서의 목적함수 값이고, 벌칙인자 R은 순차 이차계획법 시작 시 사용자가 지정하는 양수이며, V(x_k)는 이차계획문제를 이루는 모든 제약조건 중 최대 제약조건 위배량이다.

 

$$ \begin{align} \Phi \left( \mathbf{x}^{\left(k\right)} \right) = f \left( \mathbf{x}^{\left(k\right)} \right) + RV \left( \mathbf{x}^{\left(k\right)} \right) \end{align}$$

 
   벌칙인자는 최적화 반복 과정 중에 새롭게 갱신될 수 있는데, 매 단계에서 벌칙인자가 모든 라그랑지 승수의 합과 크거나 같은지 확인해야 한다. 이는 알고리즘이 수렴하기 위한 필요조건이며, 다음과 같이 나타낼 수 있다.
 

$$ \begin{align} R \geq r_k \end{align}$$

 

   위 식에서 r_k는 k번째 단계에서 이차계획문제의 모든 라그랑지 승수의 합이다. 등호제약조건의 라그랑지 승수는 부호에 제약이 없으므로 절대값을 사용한다. 반면에 부등호제약조건의 라그랑지 승수는 음수가 아니어야 한다는 것을 기억하자.

 

$$ \begin{align} r_k = \sum_{i=1}^{p} \left | v_i^{\left(k\right)}\right | + \sum_{i=1}^{m} u_i^{\left(k\right)} \end{align}$$

 
   R_k가 k번째 단계에서의 벌칙인자 값이라면, 순차 이차계획법 알고리즘이 수렴하기 위한 필요조건은 벌칙인자를 다음과 같이 선택하면 만족한다. 벌칙함수는 탐색방향 결정 단계에서 계산되며, 탐색방향이 결정된 상태인 이동거리 결정 단계에서는 그 값이 고정된다.
 

$$ \begin{align} R=\mathrm{max} \left( R_k,~r_k\right) \end{align}$$

 
   k번째 단계에서 최대 제약조건 위배량 인자 V_k는 벌칙인자와 마찬가지로 음수가 아니어야 하며, 해당 단계의 설계점에서 제약조건 함수값을 이용해 아래와 같이 결정된다. 등호제약조건은 0이 아니면 위배된 것이기 때문에 절대값을 사용하였다. 만약 모든 제약조건이 해당 설계점에서 만족되면 최대 제약조건 위배량 인자는 0이 된다.
 

$$ \begin{align} V\left( \mathbf{x}^{\left(k\right)} \right) = \mathrm{max} \left\{ 0;~~\left|h_1\left( \mathbf{x}^{\left(k\right)} \right)\right|,\left|h_2\left( \mathbf{x}^{\left(k\right)} \right)\right|,\cdots,\left|h_p\left( \mathbf{x}^{\left(k\right)} \right)\right|;~~g_1\left( \mathbf{x}^{\left(k\right)} \right),g_2\left( \mathbf{x}^{\left(k\right)} \right),\cdots,g_m\left( \mathbf{x}^{\left(k\right)} \right) \right\} \end{align}$$

 

4. 이동거리 결정

   탐색방향과 현재 설계점이 결정되면 갱신되는 설계는 다음과 같이 이동거리에 대한 일변수 함수가 된다.
 

$$ \begin{align} \mathbf{x}^{\left(k+1\right)}=\mathbf{x}^{\left(k\right)}+\alpha\mathbf{d}^{\left(k\right)} \end{align}$$

 
   강하함수에 갱신된 설계를 대입하면 강하함수 또한 이동거리에 대한 일변수 함수가 된다. 이동거리는 아래와 같은 강하함수를 최소화하는 문제를 풀이함으로써 결정된다. 이동거리 결정을 위한 풀이 방법으로는 구간감소법이 있으나, 대부분의 알고리즘에서는 부정확 이동거리 탐색법을 채택하고 있다. 아래 링크로 접속하면 이동거리를 결정하기 위한 구간감소법과 부정확 이동거리 탐색법을 확인할 수 있다.

 

$$ \begin{align} \Phi \left(\alpha\right) = \Phi \left(\mathbf{x}^{\left(k\right)}+\alpha\mathbf{d}^{\left(k\right)}\right) \end{align}$$

 

최적설계 | 비제약조건 문제 수치해법(3): 구간감소법

최적설계 | 제약조건 문제 수치해법(4): 부정확 이동거리 탐색법

 

5. 예제

   다음 비선형·제약조건 문제에 대해서 순차 이차계획문제를 적용하기 위해 초기 설계점 (40, 0.5)에서 이차계획문제를 생성해보자. 그 다음 해당 설계점에서의 탐색방향과, 이동거리에 대한 일변수 함수인 강하함수를 정의해보자. 이때 벌칙인자의 초기값은 1로 선정한다.

 

$$ \begin{align} \mathrm{minimize}~~~~&f = x_1^2 + 320x_1x_2 \\\\ \mathrm{subject~to}~~~~ &g_1 = \frac{1}{100}\left( x_1 - 60x_2\right) \leq 0 \\\\ &g_2 = 1 - \frac{x_1\left(x_1-x_2\right)}{3600} \leq 0 \\\\ &g_3 = -x_1 \leq 0 \\\\ &g_4 = -x_2 \leq 0 \end{align}$$

 

5.1. 탐색방향 결정

   탐색방향을 결정하기 위해 위 문제를 초기 설계점에서 선형화한 뒤 탐색방향의 길이를 최소화하기 위한 이차항을 더하면 다음과 같은 이차계획문제를 얻을 수 있다. 목적함수에서 상수는 최적해에 영향을 미치지 않으므로 생략하였다.
 

$$ \begin{align} \mathrm{minimize}~~~~&\overline{f} = 240d_1 + 12800d_2 + 0.5\left(d_1^2+d_2^2\right) \\\\ \mathrm{subject~to}~~~~ &\overline{g}_1 = 0.01d_1 -0.6d_2+0.1 \leq 0 \\\\ &\overline{g}_2 = -0.022d_1 +0.011d_2 +0.561 \leq 0 \\\\ &\overline{g}_3 = -d_1-40 \leq 0 \\\\ &\overline{g}_4 = -d_2-0.5 \leq 0 \end{align}$$

 

   위 이차계획문제를 행렬식으로 나타내면 다음과 같다.

 

$$ \begin{align} \mathrm{minimize}~~~~&\overline{f} = \frac{1}{2} \mathbf{d}^T \mathbf{P} \mathbf{d} + \mathbf{q}^T \mathbf{d} \\\\ \mathrm{subject~to}~~~~ &\mathbf{G} \mathbf{d} \leq \mathbf{h} \end{align}$$

 

$$ \begin{align} &\mathbf{P} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ \end{bmatrix},~~~~ \mathbf{q} = \begin{bmatrix} 240 \\ 12800 \\ \end{bmatrix}\\\\ &\mathbf{G} = \begin{bmatrix} 0.01 & -0.6 \\ -0.022 & 0.011 \\ -1 & 0 \\ 0 & -1 \end{bmatrix},~~~~ \mathbf{h} = \begin{bmatrix} -0.1 \\ -0.561 \\ 40 \\ 0.5 \end{bmatrix} \end{align}$$

 

   위 이차계획문제에 qpsolvers.solve_qp 함수를 이용하면 아래와 같은 프로그램을 작성하여 풀이할 수 있으며, 그 결과는 다음과 같다.

 

$$ \begin{align} \therefore~\mathbf{d}^{\left(0\right)} = \begin{bmatrix} 25.798 \\ 0.597 \\ \end{bmatrix} \end{align}$$

# 라이브러리 추가
import numpy as np
from qpsolvers import solve_qp

# 이차계획문제 가격계수 행렬 및 벡터 정의
P = np.array([[1.0, 0.0], [0.0, 1.0]])
q = np.array([240.0, 12800.0])

# 선형계획문제 부등호제약조건 행렬 정의
G = np.array([[0.01, -0.6], [-0.022, 0.011], [-1.0, 0.0], [0.0, -1.0]])
h = np.array([-0.1, -0.561, 40.0, 0.5])

# 이차계획문제 최적화 알고리즘
x = solve_qp(P, q, G, h, solver="cvxopt")

# 이차계획문제 최적 결과 출력
print(x)

최적설계 | 파이썬 기반 이차계획문제 알고리즘 qpsolvers.solve_qp

 

5.2. 강하함수 계산

   강하함수를 계산하기 위해 벌칙인자와 최대 제약조건 위배량 인자를 계산한다. 이차계획문제에 대하여 라그랑주 함수를 정의하고 KKT 최적성 조건을 적용하여 나타내면 다음과 같다.

 

$$ \begin{align} L~=~&\overline{f} + u_1\overline{g}_1 + u_2\overline{g}_2 + u_3\overline{g}_3 + u_4\overline{g}_4 \\\\ =~&240d_1 + 12800d_2 + 0.5\left(d_1^2 + d_2^2\right) \\\\ &+ u_1\left(0.01d_1-0.6d_2+0.1\right) \\\\ &+ u_2\left(-0.022d_1+0.011d_2+0.561\right) \\\\ &+ u_3\left(-d_1-40\right) \\\\ &+ u_4\left(-d_2-0.5\right) \end{align}$$

 

$$ \begin{align} \frac{\partial L}{\partial d_1} &= 240 + d_1 + 0.01u_1 -0.022u_2-u_3 = 0 \\\\ \frac{\partial L}{\partial d_2} &= 12800 + d_2 -0.6u_1 + 0.011u_2-u_4 = 0 \end{align}$$

 

   부등호제약조건에 탐색방향 벡터를 대입하였을 때 첫 번째와 두 번째 부등호제약조건은 활성, 세 번째와 네 번째 부등호제약조건은 만족이므로, 세 번째와 네 번째 부등호제약조건의 라그랑주 승수는 0의 값을 갖는다. KKT 최적성 조건에 이를 반영하고 탐색방향 벡터를 대입하면 다음과 같은 이차방정식을 얻을 수 있다.

 

$$ \begin{align} &0.01u_1-0.022u_2+265.798=0 \\\\ -&0.6u_1+0.011u_2+12800.597=0 \\\\ \\\\ \therefore~~&u_1=21736.968,~~u_2=21962.167 \end{align}$$

 

   이로부터 벌칙인자를 갱신하면 다음과 같다. 해당 값은 이동거리 결정 단계에서 고정된다는 것에 유의하자.

 

$$ \begin{align} R=\mathrm{max}\left( R_0,~r_0\right)=\mathrm{max}\left( 1,~43699.135\right)=43699.135 \end{align}$$

 

   마지막으로 제약조건에 초기 설계점을 대입하여 최대 제약조건 위배량 인자를 계산하면 다음과 같다.

 

$$ \begin{align} V\left(\mathbf{x}^{\left(0\right)}\right) &= \mathrm{max}\left\{ 0;~~g_1\left(\mathbf{x}^{\left(0\right)}\right),g_2\left(\mathbf{x}^{\left(0\right)}\right),g_3\left(\mathbf{x}^{\left(0\right)}\right),g_4\left(\mathbf{x}^{\left(0\right)}\right) \right\} \\\\ &= \mathrm{max}\left\{ 0;~~0.1,0.561,-40,-0.5\right\} \\\\ &= 0.561 \end{align}$$

 

   위 연산을 기반으로 해당 단계에서의 강하함수 값을 계산하면 다음과 같다.

 

$$ \begin{align} \Phi \left( \mathbf{x}^{\left(0\right)} \right) &= f \left( \mathbf{x}^{\left(0\right)} \right) + RV \left( \mathbf{x}^{\left(0\right)} \right) \\\\ &= 8000 + 43699.135 \times 0.561 \\\\ &= 32515.215 \end{align}$$

 

5.3. 이동거리 결정

   구간감소법이나 부정확 이동거리 탐색법을 활용해 해당 단계에서 강하함수를 최소화하는 이동거리를 탐색한다. 만약 탐색을 위한 최소 간격을 0.1로 선정하였다면 아래와 같이 설계점을 갱신하여 강하함수의 값을 계산한다. 이때 강하함수의 계산은 앞서 보인 과정과 동일하다.

 

$$ \begin{align} \mathbf{x}^{\left(0,1\right)} &= \mathbf{x}^{\left(0\right)} +\alpha_{0,1}\mathbf{d}^{\left(0\right)} \\\\ &= \begin{bmatrix} 40 \\ 0.5 \end{bmatrix} + 0.1 \times \begin{bmatrix} 25.798 \\ 0.597 \end{bmatrix} \\\\ &= \begin{bmatrix} 42.580 \\ 0.560 \end{bmatrix} \end{align}$$

 

$$ \begin{align} \Phi \left( \mathbf{x}^{\left(0,1\right)} \right) = f \left( \mathbf{x}^{\left(0,1\right)} \right) + RV \left( \mathbf{x}^{\left(0,1\right)} \right) \end{align}$$

 

 

 

[함께 읽으면 좋은 페이지]

 

최적설계 | 이차계획법

최적설계 | 제약조건 문제 수치해법(5): 제약최속강하법 CSD

 

 
 
 
 

참고문헌
- Arora, J. S. (2016). Introduction to optimum design. Elsevier.