반응표면방법론

애니메이션 제목

2.1. 반응표면방법론(Response Surface Methodology, RSM)

반응표면방법론(Response Surface Methodology, RSM)은 실험설계 및 분석 기법 중 하나로, 여러 개의 독립변수가 종속변수에 미치는 영향을 조사하고 최적화하는 방법입니다. 반응표면방법론은 일반적으로 실험계획법(Design of Experiments, DOE)을 기반으로 하며, 종속변수의 값과 독립변수들의 값 사이에 있는 상호작용 관계를 모델링합니다. 이러한 모델링을 통해 독립변수들의 최적 조합을 찾아내고, 그 조합에서 최적의 종속변수 값을 예측할 수 있습니다. 반응표면방법론은 대부분의 산업분야에서 사용되며, 제품의 품질 향상, 제조과정의 최적화, 비용절감 등에 큰 도움을 줍니다. 또한, 실험횟수를 줄이면서도 효과적인 실험 결과를 얻을 수 있어, 실험 비용을 절감할 수 있는 장점이 있습니다.

반응표면분석

반응표면분석(response surface method)에서는 관심영역에 속한 임의의 계량인자(quantatitative factor)들의 값에서 반응변수 (예를 들면, 에너지 효율, 합격품 수율)의 예측이 반응표면 실험결과를 분석하는 목적입니다.

일단 반응표변분석에서는 두 개의 원인변수  $x_1$과 $x_2$, 그리고  반응변수 $y$의 모평균 $\eta = {\rm E}[y]$이 다음과 같은 함수관계를 가진다고 가정한다.

$$\eta = f(x_1, x_2)$$

위의 관계에서 함수 $f(x_1, x_2)$는 모르는 함수입니다.

반응표면분석의 목적은 반응변수의 모평균 $\eta$의 최대값$\eta^*$가 나타나는 범주형 원인변수의 수준 $x_1^*$와 $x_2^*$를 주어진 영역에서 찾는 것입니다.

$$\eta^*= \max {f(x_1, x_2)}=f(x_1^*, x_2^*)$$

함수 $f(x_1, x_2)$는 복잡한 형태를 가질 수 있지만 반응표변분석에서는 최대값의 주위에서 함수 $f(x_1, x_2)$를 이차함수로 근사하여 최적점을 찾습니다.

반응표면방법론을 통한 실험계획

1. 실험목적 : 각 요인에 따른 실험조건들 중에서 모평균의 값을 최적으로 하는 실험조건 찾기 및 재현성 검토

2. 실험설계

2.1. 요인배치법(분산분석을 통하여 유의한 요인을 선별하는 방법)으로 실험설계

2.2. 최대경사법(method of steepest ascent)으로 초기 실험원인변수의 범위를 최적점 근처로 이동시킬 수 있는 쉽고 경제적인 절차 설계

초기 실험에서 고려하는 원인변수(인자, 설명변수, 독립변수)들의 영역이 결과변수(반응변수, 종속변수)가 최적값을 가지는 영역에서 멀리 떨어져 있는 경우가 흔하므로 초기 실험에서는 최적의 영역으로 이동하기 위한 원인변수들의 방향을 알아내는 것이 중요

2.3. 초기 실험 설계

실험조건(원인변수, 요인의 교차표)에서 실험대상을 처리(treatment, 중재, intervention) 후, 결과변수(반응변수)의 모평균을 추정

2.4. 축차 실험 설계

선행 실험결과를 참조하여 다음 실험을 진행하는 실험, 1차모형과 최대경사법을 이용하여 최적점 근처로 실험점을 이동

2.5. 본 실험 설계 : 찾아낸 최적조건을 포함하는 원인변수의 영역에서 곡선을 나타내는 2차 다항 모형을 반응표면분석의 모형으로 가정하고 최적점을 찾기위한 실험 실행 –  2차다항식을 적합해야 하기 때문에 각 요인에 대하여 최소한 3개의 수준이 필요

2.6. 정준분석 설계

2차 다항 모형에서 계수의 추정값에 따른 반응표면의 모양은 4가지 모양으로 분류되는 데 그 중 어떤 모양에 해당하는지 판단하기 위해 정준분석(canonocal analysis)을 적용

3. 데이터 수집

4. 데이터 분석

4.1. 찾아낸 적합한 2차 다항식 모형식에서 최적점을 구함

반응표면방법론에 의한 실험절차

1. 반응값을 나타내는 변수를 반응변수로, 계량인자를 설명변수로 간주하여 실험 자료를 이용하여 회귀분석을 실시.

2. 관심영역의 최적조건 근처에서 곡선효과 존재

3. 곡선효과를 반영하는 가장 간결한 모형인 이차 다항모형 가정

4. 적절한 모형 찾기.

5. 예측치를 최적으로 하는 최적조건을 관심영역에서 찾고, 최적조건에서의 재현성 검토하기.

6. 반응표면분석의 절차는 다음과 같이 크게 3단계로 구성된다.

7. 2수준 일부실시법에 의해서 핵심인자들을 선별하기.

8. 선별된 핵심인자들에 대한 축차적인 실험 설계(중심점을 갖는 2수준 요인 배치법)와 분석에 의해서 최적조건 근처의 설명변수들의 영역 으로 이동하기. (최대경사법 적용)

$$Y=\beta_0 +\sum_{i=1}^{p}\beta_i X_i+\varepsilon$$

9. 최적조건 근처에서 2차 모형을 가정, CCD(중심합성설계)에 의한 실험설계 및 실험자료의 회귀분석을 통한 적절한 모형 찾기 및 최적조건 찾기와 재현성 검토.

$$Y=\beta_0 +\sum_{i=1}^{k}\beta_i X_i+\sum_{i=1}^{k}\beta_{ii} X_i^2+\sum\sum_{i\lt j}^{k}\beta_{ij} X_i X_j+\varepsilon$$

2.2. 반응표면분석 모형

생산공정 등에서 많이 사용되고 있는 회귀분석의 응용 중에 반응표면분석(response surface analysis)이 있습니다. 이 방법에서는 분석에 포함된 설명변수들의 값에 따라 반응변수의 변화가 만드는 반응표면에 대한 통계적 분석를 실시하게 되는데, 보통 반응변수와 설명변수들의 관계를 일차 또는 이차 다항회귀모형으로 적합합니다. 그리고, 일반회귀분석에서와는 다르게 추정량의 유의성등에 대한 분석보다는 반응표면의 형태에 대한 분석과 반응변수의 최적조건을 찾는 것이 주요 목적입니다.

2차 다항모형을 반응표면분석에 사용하는 이유는 반응변수의 값이 최대가 되는 독립변수들의 값을 구하기 위함입니다. 이렇게 반응변수의 값이 최대가 되는 독립변수들의 값을 최적점(optimum point)이라고 부릅니다.

설명변수가 1개인 이차모형의 최적점과 극점(정상점)

먼저 간단한 설명변수가 1개인 이차모형을 살펴보자. 설명변수가 1개인 이차모형을 적합하여 다음과 같은 추정식을 얻었다고 한다면

$$\hat{Y}=\hat{\beta}_0 +\hat{\beta}_1 X+\hat{\beta}_2 X^2$$

위 식은 $X$의 이차식이므로 $Y$의 최적조건은 이 식을 $X$에 대해 미분하여 구해집니다. 즉,

$$\dfrac{d\hat{Y}}{dX}=\hat{\beta}_1 +2\hat{\beta}_2 X=0$$

에서 극점(정상점)

$$X_m =-\dfrac{\hat{\beta}_1}{2\hat{\beta}_2}$$

을 얻을 수 있으며, 이 점에서 $\hat{Y}$의 최대(최소)값은

$$\hat{Y}_m =\hat{\beta}_0 -\dfrac{\hat{\beta}_1^2}{4\hat{\beta}_2}$$

로 주어짐을 알수 있습니다. 그리고 이차 도함수가 $\dfrac{d^2\hat{Y}}{dX^2}=2\hat{\beta}_2$이므로, $\hat{\beta}_2$이 양수이면 $\hat{Y}_m$은 최소값이 되고, 음수이면 $\hat{Y}_m$은 최대값이 됩니다.

설명변수가 2개인 일차모형의 회귀평면

설명변수가 2개인 일차모형의 다음과 같은 추정식이 있다면

$$\hat{Y}=\hat{\beta}_0 +\hat{\beta}_1 X_1 +\hat{\beta}_2 X_2$$

위의 식은 $X_1$와 $X_2$의 일차함수이므로 $Y$의 최적점은 구할 수 없고 단지 증가 또는 감소하는 방향만을 구할 수 있습니다.위의 식에서 $\hat{Y}$의 값을 고정시키면 $X_1$과 $X_2$는 다음과 같은 직선의 식을 형성하는데

$$X_2 =-\dfrac{\hat{\beta}_1}{\hat{\beta}_2}X_1 +c$$

($c$는 $c=\dfrac{(\hat{Y}-\hat{\beta}_0)}{\hat{\beta}_2}$으로 정의되는 상수), 이 직선 위에 있는 모든 $X_1$과 $X_2$값들의 조합에서는 $\hat{Y}$이 같은 값을 가지며, 이 직선의 기울기와 다른 방향으로 $X_1$과 $X_2$의 값을 움직이면 $\hat{Y}$의 값은 증가 또는 감소하게 됩니다. 이 중에서 위 직선의 기울기와 직각이 되는 방향을 최대상승(steepest ascent) 또는 최대하강(steepest descent) 방향이라 하며 이 방향의 기울기는 $\dfrac{\hat{\beta}_2}{\hat{\beta}_1}$으로 주어집니다.

설명변수가 $k$개 있는 이차모형의 최적점과 극점(정상점)

일반적으로 반응표면분석에서 많이 사용되는 모형은 설명변수가 $k$개 있는 다음과 같은 이차모형입니다.

$$Y=\beta_0 +\sum_{i=1}^{k}\beta_i X_i+\sum_{i=1}^{k}\beta_{ii} X_i^2+\sum\sum_{i\lt j}^{k}\beta_{ij} X_i X_j+\varepsilon$$

그리고, 벡터와 행렬을 이용하면 위의 모형은

$$Y=\beta_0 +\boldsymbol{x}^{\prime}\boldsymbol{\beta} +\boldsymbol{x}^{\prime}\boldsymbol{B}\boldsymbol{x}+\varepsilon \tag{7.9}$$

과 같이 나타낼 수 있습니다. 위의 식에서 $\boldsymbol{x}=(X_1,X_2,\cdots,X_k)^{\prime},\ \boldsymbol{\beta}=(\beta_1,\beta_2,\cdots,\beta_k)^{\prime}$이며 $k \times k$ 행렬 $\boldsymbol{B}$는 다음과 같이 정의됩니다.

$$\boldsymbol{B}=\begin{bmatrix}\beta_{11} & \frac{1}{2}\beta_{12} & \cdots & \frac{1}{2}\beta_{1k} \\\frac{1}{2}\beta_{12} & \beta_{22} & \cdots & \frac{1}{2}\beta_{2k} \\\cdots & \cdots & \cdots & \cdots \\\frac{1}{2}\beta_{1k} & \frac{1}{2}\beta_{2k} & \cdots & \beta_{kk}\end{bmatrix}$$

모형 (7,9)는 이차모형이므로 설명변수들의 상관관계가 그렇게 크지 않아 각 설명변수를 식 (7,8)과 같이 변환한 후 중회귀모형을 위한 계산공식을 사용하여 최소제곱추정량을 구합니다. 이와 같이 적합된 회귀식을

$$\hat{Y}=\hat{\beta}_0 +\boldsymbol{x}^{\prime}\hat{\boldsymbol{\beta}} +\boldsymbol{x}^{\prime}\hat{\boldsymbol{B}}\hat{\boldsymbol{x}} \tag{7.10}$$

이라 합니다. ($\hat{\boldsymbol{B}}$은 행렬 $\boldsymbol{B}$의 각 원소를 추정량으로 대체한 행렬이다.) 반응변수 $Y$의 최적값은 역시 식(7.10)을 설명변수에 대하여 미분하여 구해지는데, 설명변수가 2개 이상인 경우에는 도함수 벡터로 주어집니다. 즉, 일차 도함수벡터를 구하여 $\boldsymbol{0}$으로 두면

$$\dfrac{\partial\hat{\boldsymbol{Y}}}{\partial \boldsymbol{x}}=\hat{\boldsymbol{\beta}} +2\hat{\boldsymbol{B}}\boldsymbol{x}=\boldsymbol{0}$$

이 되고, 여기서 극점(정상점)

$$\boldsymbol{x}_s = -\dfrac{1}{2}\hat{\boldsymbol{B}}^{-1}\hat{\boldsymbol{\beta}} \tag{7.11}$$

을 얻을 수 있으며 이 점을 정상점(stationary point)이라 합니다.

2.3. 최대경사법

데이터변환

원인변수 X의 범위를 (-1, +1)로 변환합니다.

선형회귀식

다음은 반응변수 $y$를 에측하는 모형으로 두 개의 원인변수, $x_1$과 $x_2$를 가지는 다음과 같은 일차선형 모형은 다음과 같습니다.

$$Y=\beta_0 +\beta_1 X_1 + \beta_2 X_2 + \varepsilon$$

일차선형 모형의 회귀계수들을 추정하면 반응변수의 평균. ${\rm E}[y \mid x]$에 대한 다음과 같은 예측식을 얻을 수 있습니다.

$$\widehat{{\rm E}[y \mid x]}=\hat{y}=\hat{\beta}_0+\hat{\beta}_1 x_1+\hat{\beta}_2 x_2$$

최대경사법

최대경사법(method of steepest ascent)은 공간에서 다변량 함수의 변화가 가장 크게 변하는 방향을 찾는 방법입니다. 이 방법은 경사하강법(기울기하강볍, gradient descent method)의 반대 방법입니다. 반응변수가 가장 크게 증가하는 방향을 원인변수로 만들어지는 2차원 공간 $(x_1,x_2)$에서 찾기 위해 이변량 함수 $f(x_1,x_2)$에 대한 기울기 벡터(gradient),  $\nabla f$를 정의합니다. 기울기 벡터는 다음과 같이 각 축에 대한 편미분(부분미분, partial derivative)로 이루어진 벡터입니다. $\nabla f$는 주어진 점에서 함수 $f$가 가장 빨리 증가하는 방향을 의미합니다. 따라서 기울기 벡터 $\nabla f$를 최대경사(steepest ascent) 벡터라고 부릅니다.

$$\nabla f =\nabla f(x_1,x_2) =\begin{bmatrix}
\dfrac{\partial f(x_1,x_2)}{\partial x_1} \\
\dfrac{\partial f(x_1,x_2)}{\partial x_2}
\end{bmatrix}$$

기울기 벡터 $\nabla f$에 대한 기하학적 의미는 다음과 같은 그림으로 나타낼 수 있습니다.

선형회귀의 2차원 공간 2차원 공간 $(x_1,x_2)$에서 반응변수가 가장 크게 증가하는 방향, 즉 최대경사 방향 $\nabla f$은 다음과 같이 구합니다.

$$\nabla f =\nabla f(x_1,x_2) =\begin{bmatrix}
\dfrac{\partial f(x_1,x_2)}{\partial x_1} \\
\dfrac{\partial f(x_1,x_2)}{\partial x_2}
\end{bmatrix}=\begin{bmatrix}
\dfrac{\partial (\hat{\beta}_0+\hat{\beta}_1 x_1+\hat{\beta}_2 x_2)}{\partial x_1} \\
\dfrac{\partial (\hat{\beta}_0+\hat{\beta}_1 x_1+\hat{\beta}_2 x_2)}{\partial x_2}
\end{bmatrix}=\begin{bmatrix}
\hat{\beta}_1 \\
\hat{\beta}_2
\end{bmatrix}$$

2.4. 정준분석 (stationary analysis)

반응변수 $Y$가 이 점에서 최대값 또는 최소값을 갖는지에 대한 판정을 역시 이차 도함수를 이용하는데, 이 경우에는

$$\dfrac{\partial^2\hat{\boldsymbol{Y}}}{\partial\boldsymbol{x}\partial\boldsymbol{x}^{\prime}}=2\hat{\boldsymbol{B}}$$

과 같이 행렬로 주어집니다. 행렬 $\hat{\boldsymbol{B}}$를 이용한 판정은 아래에서 설명되고 있습니다. 정상점 주위에서 반응변수가 생성하는 반응표면의 형태에 대한 분석을 정준분석(canonical analisys)이라 하며, 행렬 $\hat{\boldsymbol{B}}$에 대한 고유값분석(eigenvalue analisys)과 반응표면의 등고선그림(contour plot)을 이용합니다.

반응변수 $Y$가 정상점에서 갖는 값은

$$\hat{Y}_s =\hat{\beta}_0 +\boldsymbol{x}_s^{\prime}\hat{\boldsymbol{\beta}}+\boldsymbol{x}_s^{\prime}\hat{\boldsymbol{B}}\boldsymbol{x}_s$$

로 주어지며, 정준분석을 하기 위하여 설명변수들의 원점을 정상점 $\boldsymbol{x}_s$로 옮가는 변환을 합니다. 즉, 새로운 벡터 $\tilde{\boldsymbol{x}}$를

$$\tilde{\boldsymbol{x}}=\boldsymbol{x}-\boldsymbol{x}_s$$

라 정의하며, 식(7.10)에 대입히면

$$\eqalign{\hat{Y}&=\hat{\beta}_0 +(\boldsymbol{x}_s +\tilde{\boldsymbol{x}})^{\prime}\hat{\boldsymbol{\beta}}+(\boldsymbol{x}_s +\tilde{\boldsymbol{x}})^{\prime}\hat{\boldsymbol{B}}(\boldsymbol{x}_s +\tilde{\boldsymbol{x}})
\cr&=\hat{\beta}_0+\boldsymbol{x}_s^{\prime}\hat{\boldsymbol{\beta}}+\boldsymbol{x}_s^{\prime}\hat{\boldsymbol{B}}\boldsymbol{x}_s+\tilde{\boldsymbol{x}}^{\prime}\hat{\boldsymbol{\beta}}+2\tilde{\boldsymbol{x}}^{\prime}\hat{\boldsymbol{B}}\boldsymbol{x}_s +\tilde{\boldsymbol{x}}^{\prime}\hat{\boldsymbol{B}}\tilde{\boldsymbol{x}}\cr&=\hat{Y}_s +\tilde{\boldsymbol{x}}^{\prime}\hat{\boldsymbol{B}}\tilde{\boldsymbol{x}}+\tilde{\boldsymbol{x}}^{\prime}(\hat{\boldsymbol{\beta}}+2\hat{\boldsymbol{B}}\boldsymbol{x}_s)\cr&=\hat{Y}_s +\tilde{\boldsymbol{x}}^{\prime}\hat{\boldsymbol{B}}\tilde{\boldsymbol{x}}} \tag{7.12}$$

의 식을 얻는다. (마지막 등식은 식(7.11)에 의하여 $\hat{\boldsymbol{\beta}}+2\hat{\boldsymbol{B}}\boldsymbol{x}_s=\boldsymbol{0}$이므로 성립) 그러므로, 이 식에서 행렬 $\hat{\boldsymbol{B}}$이 양정치이면(즉, 모든 $\boldsymbol{0}$ 아닌 벡터 $\boldsymbol{w}$에 대해 $\boldsymbol{w}^{\prime}\hat{\boldsymbol{B}}\boldsymbol{w} \gt 0$) 정상점에서의 값 $\hat{Y}_s$은 최소값이 되고, 음정치이면 $\hat{Y}_s$은 최대값이 됨을 알 수 있다.행렬의 양정치 여부를 판별하는 간단한 방법은 행렬의 고유값을 구해보는 것으로, 고유값이 모두 양수이면 양정치이고, 모두 음수이면 음정치가 됩니다. 그러므로 $\lambda_1,\lambda_2,\cdots,\lambda_k$를 행렬 $\hat{\boldsymbol{B}}$의 $k$개 고유값이라 하면 다음과 같이 판별할 수 있습니다.

1) 모든 $\lambda_j$가 양수이면 정상점 $\boldsymbol{x}_s$는 최소점이 되며, $\hat{Y}$은 최소값 $\hat{Y}_s$을 갖습니다.
2) 모든 $\lambda_j$가 음수이면 정상점 $\boldsymbol{x}_s$는 최대점이 되며, $\hat{Y}$은 최대값 $\hat{Y}_s$을 갖습니다.
3) $\lambda_j$가 일부는 양수이고, 일부는 음수이면, $\hat{Y}$은 $\boldsymbol{x}_s$에서 벗어나는 방향에 따라 중가 또는 감소하게 되는데, 이 때 $\boldsymbol{x}_s$를 안장점(안부점, saddle point)이라 합니다.

정준분석을 이용한 정상점의 판단

$$\eqalign{\hat{Y}&=\hat{Y}_s +\tilde{\boldsymbol{x}}^{\prime}\hat{\boldsymbol{B}}\tilde{\boldsymbol{x}}\cr&=\hat{Y}_s+\lambda_1(\boldsymbol{\upsilon}_1^{\prime}\tilde{\boldsymbol{x}})^2+\lambda_2(\boldsymbol{\upsilon}_2^{\prime}\tilde{\boldsymbol{x}})^2 +\cdots +\lambda_k(\boldsymbol{\upsilon}_k^{\prime}\tilde{\boldsymbol{x}})^2\cr&=\hat{Y}_s +\lambda_1 w_1^2 +\lambda_2 w_2^2 +\cdots +\lambda_k w_k^2} \tag{7.13}$$

의 정준형태(canonical form)로 주어짐을 보일 수 있습니다. 위의 식에서 $w_j$는 $w_j = \boldsymbol{\upsilon}_j^{\prime}\tilde{\boldsymbol{x}}$로 정의되며 정준변수(canonical variable)라 부릅니다. 즉, 원래 변수를 고유벡터들로 회전하여 얻어지는 새로운 변수라 할 수 있습니다. 따라서, 고유값의 부호로 새로운 변수의 축 방향으로의 증가와 감소 여부를 알 수 있고, 고유값의 크기로는 경사의 정도를 파악 할 수 있습니다. 유의할 점은 설명변수들의 단위가 다르면 식 (7.13)에서 고유값들의 크기를 상대적으로 비교할 수 없으므로, 정준분석을 할 때는 설명변수들을 식 (7.8)으로 표준화하여 사용할 필요가 있습니다.

설명변수의 개수가 작은 경우에는 등고선 그림과 3차원 그림을 이용하면 편리합니다. 등고선 그림(contour plot)은 반응표면에서 $\hat{Y}$의 값이 같은 점들을 연결하여 선으로 나타낸 것으로, $\hat{Y}$이 증가 또는 감소하는 방향과 속도(경사)를 쉽게 파악 할 수 있게 해줍니다. 그러나, 설명변수 개수가 많은 경우에는 2차원 또는 3차원 그림에 축이 설정 되지 않는 다른 변수들의 값을 고정시켜야 하므로 단면만을 봐야하는 제약이 따릅니다. 아래그림은 등고선 그림으로 나타낸 반응표면의 여러가지 형태를 보여주고 있습니다.

2.5. 중심합성설계(CCD)

반응표면분석에서 원인변수(독립변수)의 개수가 $k$개이면 한 원인변수(요인)에 대하여 최소한 3개의 변수값(수준, 카테고리, 집단, 범주)이 필요합니다. 최소한 3개의 수준이 필요한 이유는 2차 다항식을 특정하기 위해서는 최소한 3개의 데이터가 필요하기 떄문입니다. 따라서, 요인의 개수가 $k$이면 실험점(범주)의 개수는 최소한 $3^k$개가 필요합니다. 즉, 실험점의 개수는  요인의 수에 따라 기하급수적으로 늘어납니다. 따라서 비용이 수반하는 실험의 횟수를 결정하는 실험점은 경제적으로 설계되어야 합니다. 

반응표면분석에서는 실험을 단계적으로(sequentially) 실행하면서 최적점 근처에서 더 많은 실험을 수행할 수 있는 실험설계를 합니다. 그에 따른 실험절차는 다음과 같습니다.

1차 모형 적합을 위한 1단계 실험

– 먼저 중심점(center points)에서 $n_0$개의 실험을 수행하고 각 요인의 상자점들(factional points)에서 $2^k$번의 실험을 실시

– 이러한 기초 실험은 최대경사법을 이용하는 중간 과정에서 언제나 수행해야 하는 실험이다.

– 기초실허은 1차 모형을 적합하기 위한 실험설계이다.

– 기초 실험에서 실험 구간에 최적점이 포함되어 있다고 판단되면 2차 모형 적합을 위한 2단계 실험을 실시한다.

– 1단계 실험의 개수는 상자점 $2^k$개

2차 모형 적합을 위한 2단계 실험

–  최적점이 가까워진 경우 2차 모형 적합을 위해 1단계 실험의 실험점들에 추가적인 실험점을 더하여 2단계 실험을 진행한다.

– 2차 모형을 적합하기 위한 실험으로 각 요인에 대하여 2개의 축점(axial points, star points)을 추가한다. 기초실험에서 수행한 실험점들과 합해져서 각 요인마다 3개의 수준을 가지는 실험점을 만든다.

– 2단계 실험에서 추가한 축점은 2차 모형의 효율적인 추정을 고려하여 선택한다.

– 2단계 실험의 개수는 축점 $2k$개

반응표면분석에서 위와 같이 최적점 탐색을 하는 경우 사용되는 대표적인 실험설계는 중심합성설계(CCE : Central Composite design)이다.

그림 5.1 구형 계획법을 이용한 중심합성설계

중심합성설계에서 축차적으로 실험을 실시한 경우 $k$개의 요인을 고려하면 실험점들은 그 특성에 따라서 다음과 같이 나타납니다.

2.6. 균등정밀 중심합성설계

중심합성설계에서 2차 다항식을 추정하기 위한 2단계 실험에서 추가하는 축점(axial points)을 어떻게 선택하느냐가 중요합니다.

축점을 배치할 때 중요한 고려사항은 2차 다항식에서 얻은 추정값들의 분산이 각 실험접들에서 동일하게 나타나게 하는 것입니다. 이러한 성질을 회전 가능성(Rotatibility) 라고 부른다. 이러한 회전 가능성이 만족하면 고려한 모든 실험점들에서 구한 예측값들의 정도(precision)가 같다는 의미입니다.

회전가능성이 중요한 이유는 반응표면분석이 반응의 최적점을 찾는 실험이고 최적점의 위치는 알 수 없으므로 모든 방향에 대한 예측값의 정도를 동일하게 설정하는 것이 합리적이기 때문입니다.

2차 다항식도 선형모형에 속하므로 다음과 같은 선형모형을 고려할 때

$$y=X\beta + e$$

최소제곱법으로 얻은 회귀계수의 추정량을 $\hat \beta$이라고 하자. 만일 관심이 있는 실험점 $x$에서 반응값의 예측치는 $\hat y=x^t\hat \beta$이다. 또한 예측치의 분산은 다음과 같이 주어집니다.

$${\rm Var}[\hat y]= {\rm Var}[\hat y | x] = \sigma^2x^t(X^tX)^{-1}x$$

“중심합성설계가 회전 가능하다”는 의미는 실험에서 사용한 모든 실험점들에서 반응변수 예측값의 분산이 동일하다는 것을 의미합니다. 

$$ {\rm Var}[\hat y | x_i] = {\rm Var}[\hat y | x_j]$$

여기서,  모든 $i$, $j$에 대해 등식이 성립

이렇게 실험에서 고려한 모든 실험점에서 예측값의 분산이 같은 실험설계를 ‘회전가능하다’(rotatable) 라고 말하며 일반적으로 균등정밀 중심합성설계(uniform precision CCD)라고 부릅니다.

중심합성설계가 회전가능하게 되는 조건은 일반적으로 2차 실험에서 추가되는 축점의 길이, $\alpha$가 결정합니다. 상자점의 수가 F개인 경우 회전가능한 중심합성설계를 얻기 위해서는 축점길이($\alpha$)를 $F^{\frac{1}{4}}$로 설정하면 됩니다.

다음은 요인의 개수 $k$에 대한 균등 정밀 중심합성설계의 각 실험점의 개수와 축점의 길이 $\alpha$를 나타낸 표입니다.

균등 정밀 중심합성설계의 각 실험점의 개수와 축점의 거리

참고로 중심점 (0, 0)에서의 실험의 수 $n_0$는 일반적으로 3개에서 5개입니다.

2.7. 구형 실험설계

반응표면분석에서 회전 가능한 계획법은 실험에서 고려한 모든 실험점들의 분산을 같게 하는 실험계획입니다. 하지만 실제 현장에서는 분산이 모두 동일한 경우 보다는 모든 실험점이 원점에서 같은 거리에 있는 경우를 선호하는 경우도 있습니다.

고려한 모든 실험점이 중심점으로부터 거리가 동일한 계획을 구형 계획법(spherical design)이라고 합니다. 요인의 개수가 $k$개인 경우, 구형 계획법에서 축점의 거리는 $\alpha=\sqrt{k}$입니다.

2.8. Box-Benken 실험설계

Box-Benken 설계는 2차 모형 적합을 위한 효율적인 3수준 실험설계로 3수준 요인배치법의 일부 실험조건에서만 실험을 실시하는 것입니다.

예를 들어 3개의 요인가 있다면 관심영역 상자점에서 2개의 요인를 먼저 선택하고 $2^2$개의 (±1, ±1)  수준에서 실험을 하고 나머지 인자는 수준의 중앙값인 0으로 고정합니다.

이렇게 고려하는 인자의 개수가 $k$개인 경우, 2개의 인자를 선택하는 조합의 수는 $\dfrac{k(k-1)}{2}$개이다. 또한 각 조합마다 4개의 실험점(±1, ±1)  이 추가되므로 총 실험점의 개수는 다음과 같습니다.

$$\dfrac{k(k-1)}{2}+n_0$$

Box-Benken 설계는 선택된 2개의 조합이 블럭으로 나타나는 블럭 일부실시법입니다.

요인의 개수가 3개인 경우 Box-Benken 설계

요인과 블록에 따른 위치값

2.9. 실험설계에 따른 실험점의 개수

요인의 개수가 $k$개인 경우 각 실험계획법에 대한 실험점의 개수는 다음의 표에 있는 공식으로 계산할 수 있습니다. $n_c$는 중심점 (0,0)에서 실험의 개수입니다

2차 반응표면분석 모형 적합을 위한 실험설계의 실험점 개수 비교

2.10. 반응표면방법론(Response Surface Methodology)을 수행하는 데 사용되는 파이썬 패키지

“rsm”은 Response Surface Methodology를 수행하는 데 사용되는 Python package입니다. 이 패키지는 실험계획부터 모델링, 최적화 및 시각화까지 반응표면 방법론 분석을 위한 다양한 기능을 제공합니다. 주요 기능은 다음과 같습니다.

(1) 실험계획 생성: 중심합성법(Central Composite Design) 및 박스-벤킨법(Box-Behnken Design) 등의 실험계획법 생성 기능을 제공합니다.

(2) 모델링: 실험 데이터에 대한 반응표면 모델을 적합하는 기능을 제공합니다. 이를 위해 다중선형회귀분석, 이차방정식 모델, 다항식 모델, 신경망 모델 등을 사용할 수 있습니다.

(3) 분산분석(ANOVA): 반응표면 모델의 항들이 유의한지 검정하는 기능을 제공합니다.

(4) 최적화: 반응표면 모델을 사용하여 반응변수의 최적 조건을 찾는 기능을 제공합니다. 이를 위해 그리드 탐색, 단순 최적화 알고리즘 등을 사용할 수 있습니다.

(5) 시각화: 실험 데이터와 반응표면 모델을 시각화하여 분석결과를 확인하는 기능을 제공합니다. 이를 위해 2D 및 3D 플롯, 등고선 플롯 등을 사용할 수 있습니다.

“rsm” 패키지는 NumPy, SciPy, Matplotlib, Pandas 등의 패키지를 기반으로 하며, PyPI(Python Package Index)에서 제공됩니다. 패키지의 문서에는 설치 방법 및 사용 예제가 포함되어 있습니다.

https://pypi.org/project/rsm/

3.2. 함수

=ROWS(F2:F2) : 지정된 배열 또는 범위에 있는 행의 개수.

4.1 용어

반응표면방법론(Response surface methodology, RSM)

통계학에서 반응표면방법론(RSM)은 여러 설명변수(원인변수)와 하나 이상의 반응변수(결과변수)간의 관계를 탐색합니다.  이 방법론은 1951년 George E. P. Box와 K. B. Wilson에 의해 도입되었습니다. RSM의 주요 아이디어는 일련의 설계된 실험을 사용하여 최적의 응답을 얻는 것입니다. Box와 Wilson은 이를 위해 2차 다항식 모델을 사용할 것을 제안합니다. 그들은 이 모델이 근사치일 뿐이라는 것을 인정하지만 이 모델을 추정하고 추정값을 적용하기가 쉽습니다. 그래서 반응표면방법론은 프로세스에 대해 거의 알려지지 않은 경우에도 사용가능합니다. RSM과 같은 통계적 접근 방식을 사용하여 운영 요소를 최적화하여 특수 물질의 생산을 극대화할 수 있습니다. 최근에는 제형 최적화를 위해 적절한 실험설계(DoE)를 사용하는 RSM이 광범위하게 사용되었습니다. 기존의 방법과 달리 프로세스 변수 간의 상호 작용은 통계적 기법으로 확인할 수 있습니다.

Reference

 WIkipedia : Response surface methodology

중심합성설계(Central Composite Design)

통계학에서 중심합성설계(Central Composite Design)는 반응표면방법론(Response Surface Mehthodology)에 포함된 유용한 실험설계입니다. 중심합성설계는  완전한  3수준 요인에 따른 실험을 하지 않고 반응변수에 대한 2차(quadratic) 모형을 구축하기 위한 실험설계입니다.

설계된 실험이 수행된 후 선형회귀가 적용되며 때로는 결과를 얻기 위해 반복적으로 실험이 수행됩니다. 코드화된 변수는 중심합성설계를 할 때 자주 사용됩니다.

Reference

 WIkipedia : Central Composite Design

4.2. 참조

https://ilovedata.github.io/teaching/experiment/response/steepdescent.html

https://angeloyeo.github.io/2020/08/16/gradient_descent.html