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

1. 제약최속강하법

   앞서 다루었던 최속강하법에서 탐색방향은 목적함수가 가장 빠르게 감소하는 방향, 즉 목적함수의 경사도 벡터와 반대 방향이었다. 그러나 이전과 다르게 최적 설계문제에 제약조건이 추가되면 탐색방향은 목적함수가 감소하는 방향이면서도 선형화된 제약조건 또한 만족해야 한다. 이는 최속강하법에서 제약조건을 고려하여 탐색방향을 수정하기 때문에 이를 제약최속강하법(Constrained Steepest Descent, CSD)이라고 부른다. 제약최속강하법의 풀이와 알고리즘은 순차 이차계획법과 동일하다.
 

2. 알고리즘

제약최속강하법 알고리즘

   제약최속강하법을 시작하기에 앞서 종료 판정기준 두 가지를 정의한다. 첫 번째 기준은 (1) 최대 제약조건 위배량에 대한 허용 오차로, 해당 설계점에서의 등호제약조건과 부등호제약조건 위배량 중 가장 큰 값이 해당 기준보다 작거나 같아야 한다. 두 번째 기준은 (2) 탐색방향 벡터의 길이를 0으로 만들기 위한 것으로, 탐색방향 벡터의 길이가 해당 기준보다 작거나 같아야 한다. 반복 과정을 종료하기 위해서는 두 종료 판정기준 모두를 만족해야 한다. 이제 벌칙인자 초기값과 초기 설계점을 시작으로 종료 판정기준을 만족할 때까지 탐색을 반복한다. 현재 설계점에서 목적함수와 제약조건의 함수값과 경사도를 계산하고, 선형 테일러 급수 전개를 통해 문제를 국소적으로 선형화한다. 탐색방향 벡터에 대해 이차계획문제를 구성하여 풀이하여 탐색방향 벡터를 구하고, 제약조건의 함수값으로 최대 제약조건 위배량을 계산한다. 종료 판정기준을 만족하지 않는다면 라그랑지 승수의 합과 벌칙인자를 비교하여 벌칙인자의 값을 갱신하고, 부정확 이동거리 탐색법을 적용해 최대한 멀리 이동한다. 새로운 설계점을 갱신하고 종료 판정기준을 만족할 때가지 탐색을 반복한다.
 

3. 예제

   초기 설계점 (1, 1)에서 다음 설계문제를 최적화해보자.
 

$$ \begin{align} \mathrm{minimize}~~~~&f = x_1^2 + x_2^2 - 3x_1x_2 \\\\ \mathrm{subject~to}~~~~ &g_1 = \frac{1}{6}x_1^2 + \frac{1}{6}x_2^2 - 1 \leq 0 \\\\ &g_2 = -x_1 \leq 0 \\\\ &g_3 = -x_2 \leq 0 \end{align}$$

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

$$ \begin{align} \mathrm{minimize}~~~~&\overline{f} = -d_1 -d_2 + 0.5\left(d_1^2+d_2^2\right) \\\\ \mathrm{subject~to}~~~~ &\overline{g}_1 = \frac{1}{3}d_1 + \frac{1}{3}d_2 - \frac{2}{3} \leq 0 \\\\ &\overline{g}_2 = -d_1-1 \leq 0 \\\\ &\overline{g}_3 = -d_2-1 \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} -1 \\ -1 \\ \end{bmatrix}\\\\ &\mathbf{G} = \begin{bmatrix} \frac{1}{3} & \frac{1}{3} \\ -1 & 0 \\ 0 & -1 \\ \end{bmatrix},~~~~ \mathbf{h} = \begin{bmatrix} \frac{2}{3} \\ 1 \\ 1 \\ \end{bmatrix} \end{align}$$

 
   위 이차계획문제에 qpsolvers.solve_qp 함수를 이용하면 아래와 같은 프로그램을 작성하여 풀이할 수 있으며, 그 결과는 다음과 같다. 해당 탐색방향 벡터는 알고리즘의 종료 판정조건을 만족하지 않는다.
 

$$ \begin{align} \therefore~\mathbf{d}^{\left(0\right)} = \begin{bmatrix} 1.0 \\ 1.0 \\ \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([-1.0, -1.0])

# 선형계획문제 부등호제약조건 행렬 정의
G = np.array([[1/3, 1/3], [-1.0, 0.0], [0.0, -1.0]])
h = np.array([2/3, 1.0, 1.0])

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

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

 
   현재 설계점에서 최대 제약조건 위배량 인자를 계산하면 다음과 같다. 해당 최대 제약조건 위배량 인자는 알고리즘의 종료 판정조건을 만족한다.

 

$$ \begin{align} V\left(\mathbf{x}^{\left(k\right)}\right) &= \mathrm{max}\left\{ 0;~~g_1\left(\mathbf{x}^{\left(k\right)}\right),g_2\left(\mathbf{x}^{\left(k\right)}\right),g_3\left(\mathbf{x}^{\left(k\right)}\right) \right\} \\\\ &= \mathrm{max}\left\{ 0;~~-\frac{2}{3}, -1, -1\right\} \\\\ &= 0 \end{align}$$

 

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

$$ \begin{align} L~=~&\overline{f} + u_1\overline{g}_1 + u_2\overline{g}_2 + u_3\overline{g}_3 \\\\ =~&-d_1 -d_2 + 0.5\left(d_1^2 + d_2^2\right) \\\\ &+ u_1\left(\frac{1}{3}d_1 + \frac{1}{3}d_2 -\frac{2}{3}\right) \\\\ &+ u_2\left(-d_1 -1\right) \\\\ &+ u_3\left(-d_2 -1\right) \end{align}$$

 

$$ \begin{align} \frac{\partial L}{\partial d_1} &= -1 + d_1 + \frac{1}{3}u_1 -u_2 = 0 \\\\ \frac{\partial L}{\partial d_2} &= -1 + d_2 + \frac{1}{3}u_1 -u_3 = 0 \end{align}$$

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

$$ \begin{align} u_1=6 \end{align}$$

 
   이로부터 벌칙인자를 갱신하면 다음과 같다. 
 

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

   부정확 이동거리 탐색법을 활용해 해당 단계에서 강하함수를 최소화하는 이동거리를 탐색한다. 이를 위해 현재 설계대안에서 강하함수를 계산하면 다음과 같다.

 

$$ \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) \\\\ &= -1 + 6 \times 0 \\\\ &= -1 \end{align}$$

 

   기본 시험 이동거리를 시작으로 이동거리 탐색을 실시한다. 현재 단계에서 주어진 탐색방향을 따라 시험하는 설계점은 다음과 같이 계산된다.

 

$$ \begin{align} \mathbf{x}^{\left(k+1,0\right)} &= \mathbf{x}^{\left(k\right)} + t_0 \mathbf{d}^{\left(k\right)} \\\\ &= \begin{bmatrix} 1 \\ 1 \end{bmatrix} + 1 \times \begin{bmatrix} 1 \\ 1 \end{bmatrix} \\\\ &= \begin{bmatrix} 2 \\ 2 \end{bmatrix} \end{align}$$

   강하조건을 만족하는지 확인하기 위해 시험 설계점에서 강하함수를 계산한다. 먼저 시험 설계점에서 최대 제약조건 위배량 인자를 계산하면 다음과 같다.
 

$$ \begin{align} V\left(\mathbf{x}^{\left(k+1,0\right)}\right) &= \mathrm{max}\left\{ 0;~~g_1\left(\mathbf{x}^{\left(k+1,0\right)}\right),g_2\left(\mathbf{x}^{\left(k+1,0\right)}\right),g_3\left(\mathbf{x}^{\left(k+1,0\right)}\right) \right\} \\\\ &= \mathrm{max}\left\{ 0;~~\frac{1}{3},-2,-2\right\} \\\\ &= \frac{1}{3} \end{align}$$

   위 연산을 기반으로 시험 설계점에 대해서 강하함수를 계산하면 다음과 같다.

 

$$ \begin{align} \Phi \left( \mathbf{x}^{\left(k+1,0\right)} \right) &= f \left( \mathbf{x}^{\left(k+1,0\right)} \right) + RV \left( \mathbf{x}^{\left(k+1,0\right)} \right) \\\\ &= -4 + 6 \times \frac{1}{3} \\\\ &= -2 \end{align}$$

   시험 설계점에 대해서 강하조건을 검토한다. 시험 설계점에서의 강하함수를 초기 설계점에서의 강하함수와 비교하였을 때 그 값이 감소하였으므로, 즉 강하조건을 만족하므로 해당 시험 이동거리는 수용 가능하다. 
 

$$ \begin{align} \mathbf{\Phi} \left(\mathbf{x}^{\left(k+1,0\right)}\right) + t_0 \gamma \left\| \mathbf{d}^{\left(k\right)} \right\|^2 &= -2 + 1 \times 0.5 \times \left( 1^2 + 1^2 \right) \\\\ &= -1 \\\\ &\leq \mathbf{\Phi} \left(\mathbf{x}^{\left(k\right)}\right) \\\\ &= -1 \end{align}$$

   따라서 해당 단계에서의 이동거리는 다음과 같이 결정된다. 새로운 설계점 (2, 2)에서 탐색방향과 이동거리를 구하여 최적점을 구하기 위한 반복을 수행한다.
 

$$ \begin{align} \alpha_0 = t_0 = 1 \end{align}$$

 
 

 

 

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

 

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

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

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

 

 

 

 

 

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