회귀적 예측오차방법의 매개변수의 영향에 관한 연구

임순혁,김홍찬,양대륙

고려대학교 화학공학과

Research on the effect of the parameter of RPE Method

Soon Hyuk Im, Hong Chan Kim, Dae Ryook Yang

Department of Chemical Engineering, Korea University

1.서론

증류공정의 혼합물의 비점차를 이용하여 물질을 분리하는 방법은 그 자체가 열역학적인 측면에서 비선형적이고, 탑상부와 하부제품의 농도조절은 상호간섭이 심하여 각각의 농도제어에는 어려움이 많다. 또한 원료의 조성과 외부환경의 변화 등을 통해 끊임없이 외란이 유입된다. 이러한 공정자체의 비선형성과 농도제어환 사이의 상호간섭이 존재하는 단일 증류공정에 대하여, 공정특성을 적절히 표현할 수 있는 물리적 의미를 가진 공정모델을 개발하고 이를 이용한 적응예측제어를 수행하고자 한다. 증류공정에 대한 엄밀모델의 경우 일반적으로 복잡한 열역학적 관계식과 함께 방대한 수의 미분방정식으로 표현된다. 본 연구에서는 엄밀모델과 유사한 정확도를 가지며 계산시간을 단축시키기 위한 수단인 차수감소모델, 특히 배열법(Collocation method)을 사용한다. 배열법은 불연속적인 변수들을 연속적인 것으로 가정하여 편미분 방정식으로 표현하고, 내삽함수를 도입하여 공정수식을 배열점에서 재구성하는 것으로, 공정의 물리적 수지식을 그대로 보존하며 여기에 수학적 근사만이 도입되었다는 점에서 주목할 만하다. 본 연구에서는 모델예측기법으로 사용하는 회귀적 예측방법의 변수변화가 적응예측제어에 미치는 영향을 살펴보고자 한다.

2.본론

이성분 증류공정의 엄밀모델을 차수감소기법에 적용하여 단순모델로 만들기 위한 기본 가정들은 다음과 같다.

□Constant molar overflow

□Liquid holdup in each tray is constant

□Vapor holdup is negligible

□Vapor leaves a tray at equilibrium with liquid

□The holdups in the accumulators are constant

이성분계 증류공정의 동특성방정식은 다음과 같다.

응축기

[eqn: scale 80 rm H_c {dx_1}over {dt} =V_r (y_1 -x_1 )]	(1)

각 단에서,

[eqn: scale 80 rm H_m {dx_j } over {dt} = V_m (y_j -y_j-1 ) +L_m (x_j-1 - x_j ) +q F_1 x_f +(1-q)F_2 x_f]	(2)

여기서 j=NF인경우 F1=F, j=NF-1인 경우 F2=F, 다른 모든경우 0

재비기에서

[eqn: scale 80 rm H_b {dx_NT} over {dt} = L_s x_NT-1 - V_s y_NT-1 -Bx_NT]	(3)

다른 관계식

[eqn: scale 80 rm &L_s =L_r +qF#rm&V_s =V_r - (1-q) F#rm&y_j =K_j x_j+1]	(4)

차수감소모델

배열법은 정류부와 탈거부에 대하여 연속 공간좌표 z를 설정하고, 모든변수를 z에 대하여 연속으로 가정한다. 다항식의 내삽함수로써 Cubic spline을 사용하고, 단순모델식을 임의의 배열점에서 재구성하면 다항식 근사에 의한 이성분 증류공정의 차수감소모델을 얻을 수 있다. 차수감소모델의 구성에 필요한 내삽함수의 일반적인 형태는 다음과 같이 쓸수 있다.

[eqn: scale 80 rm x(z) = SUM from { { p}=1} to n S_p (z) x(z_p )]	(5)

수식표현의 복잡성을 덜기 위해, 다음의 관계를 정의한다.

[eqn: scale 80 rm &a_p (z) =S_p (z- DELTA z) scale 80 rm&b_p (z) =S_p (z+DELTAz)#scale 80 rm&c_p (z) =S_p (z) - S_p (z-DELTAz)scale 80 rm&d_p (z) =S_p (z) - S_p (z+DELTA z)]	(6)

여기서 □z는 단간격이고, 정류부에서는 1/(NF-2)이고, 탈거부에서는 1/(NT-NF +1) 이다.

응축기

[eqn: scale 80 rm {dx_r (z_1 )} over {dt}= V_r over H_c LEFT [ y_r (z_1 ) -x_r (z_1 ) RIGHT ]]	(7)

정류단

[eqn: scale 80 rm {dx_r (z_p ) } over dt = 1 over H_r SUM from { { j}=1} to nr} rmc_j (z_p ) [V_r y_r (z_j ) -L_r x_r (z_j ) ]], p=2,…,nr	(8)

탈거단와 재비기

[eqn: scale 80 rm {dx_s (z_p ) } over dt = 1 over H_s SUM from { { j}=1} to ns d_j (z_p ) [ V_s y_s (z_j ) - L_s x_s (z_j ) ]], p=1,…,ns-1	(9)

여기서  p=1일때, Hs =Hb

경계조건으로,

[eqn: scale 80 rm&y_s (z_1 ) = x_s (z_1 )#rm&L_s x_s (z_ns ) = L_r x_r (z_nr ) +q Fx_f#rm&V_s y_s (z_ns ) = V_r y_r (z_nr ) - (1-q)F x_f]	(10)

회귀적 예측오차방법(recursive prediction error method, RPE)

차수감소모델이 정상상태와 동적상태에서 단순모델을 적절히 근사하였더라도, 여러 가지 가정과 단순화를 통하여 이루어진 것이므로 실제 공정의 운전상황을 완전하게 예측하는 것을 기대하기는 힘들다. 따라서 차수감소모델이 포함하고 있는 미측정변수를 조정변수로 취하고, 공정조업자료로부터 이 값들을 추정하는 실시간대 모델인식을 수행하는 방법을 택한다. 차수감소모델이 액상 몰분율에 대하여 선형이고 공정변수는 비선형적으로 포함되어 있으므로, 이와 같은 경우에 확장 칼만휠터를 사용하여 모델인식을 전개한다.

공정변수 □를 가진 일반적인 선형상태모델은 다음과 같이 쓸 수 있다.

[eqn: scale 80 rm &x(t+1) = A(theta)x(t) + B(theta)u(t) + v(t) #scale 80 rm &y(t) = C(theta)x(t) + e(t) #scale 80 rm&~~~E[v(t) v(s) ^T ] = R_1 delta_ts,~~ E[e(t) e(s) ^T ] = R_2 delta_ts, ~~E[v(t) e(s) ^T ] = R_12 delta_ts]	(11)

여기서 공정의 잡음에 관한 정보인 R1, R2, R12를 정확하게 예측하지 못하면 추정변수가 발산하거나 bias되는 경향이 있는데 이를 보완하기 위해 정상상태 Kalman이득을 변수화하여 모델에 첨부한 innovation predictor model을 사용한다.

[eqn: scale 80 rmx(t+1) &= A(theta) x(t) + B(theta)u(t) + {bar G } (theta) epsilon(t) #scale 80 rm y(t) &= C(theta)x(t) + epsilon(t)]	(12)

일반적인 시불변계의 선형모델은 입출력자료의 선형휠터를 통한 휠터링에 의해 다음과 같이 쓸 수 있다.

[eqn: scale 80 rm&phi(t+1,theta) = {tilde F} (theta) phi (t,theta) + {tilde G (theta)} z(t)# scale 80 rm&~~~hat y (t|theta)```` = {tilde H} (theta) phi(t, theta)]	(13)

공정변수 □가 예측치[eqn: scale 80 rmhat y (t | theta)]에 주는 영향을 수식화하기 위해 모델이 공정변수 □로 미분이 가능하다고 하고, 다음과 같이 정의한다.

[eqn: scale 80 rm&left [d over {d theta} hat y (t|theta) right] ^T = PSI (t, theta)]	(14)

복잡성을 덜기 위해 다음과 같은 행렬을 도입한다.

[eqn: scale 80 rm& PARTIAL over {partial theta} left[ tilde F (theta) phi + tilde G (theta) z right ] = M(theta, phi, z)#scale 80 rm& ~~~~~~~~partial over {partial theta} left [ tilde H (theta) phi right] = D (theta, phi)]	(15)

Innovation 예측모델에 회귀적 예측오차방법을 적용한 알고리듬은 다음과 같다.

[eqn: scale 80 rm&epsilon (t) ``= y(t) - hat y (t) #scale 80 rm&hat LAMBDA (t) = hat LAMBDA (t-1) + gamma(t) [ epsilon(t) epsilon supT (t) - hat LAMBDA (T-1)] #scale 80 rm&`R(t) = `R(t-1) + gamma(t) [ PSI(t){hat LAMBDA}^-1 (t) PSI ^T (t) - R(t-1)] #scale 80 rm&hat theta (t) `= hat theta (t-1) + gamma(t) R ^-1 (t) PSI(t) {hat LAMBDA} ^-1 (t) epsilon(t) #scale 80 rm& hat x (t+1) = A sub t hat x (t) + B _ t u(t) + G _t epsilon(t) #scale 80 rm& hat y (t+1) = C _t hat x (t+1) #scale 80 rm&W(t+1) = [A_t - G_t C_t ]W(t) + {bar M} _t - G_t D_t #scale 80 rm&PSI(t+1) = W ^T (t+1) C_t ^T + D ^T (hat theta (t) , hat x (t+1)) #scale 80 rm&~~~~~where, A _t = A(hat theta (t)), B_t = B( hat theta (t)) ,C_t = C (hat theta (t)), G_t = G (hat theta (t)) #scale 80 rm&~~~~~~~~~~```{bar M} _t = {bar M} ( hat theta (t), hat x (t), u(t), epsilon (t)), D _t = D ( hat theta (t), hat x (t))]	(16)

3.결과 및 토의

공정모사 data는 다음과 같다.

표  물-메탄올 이성분계의 운전조건

No. of trays Feed tray Normal holups	13 7th from the condenser Hc=15 mol, Hr=2.5 mol, Hs=3.5 mol, HB=23 mol
Feed flow rate Feed condition Top product flow rate Reflux flow rae	216 mol/hr xF=0.5, q=1.001 108 mol/hr 60 mol/hr
Distibution constants from top tray to the reboiler	1.1247 1.2087 1.2940 1.3752 1.4470 1.5063 1.5523 1.5537 1.5582 1.5730 1.6227 2.1007 3.0235 4.9596
xF change	T<20,xF=0.5: T<200,xF=0.45: T<400,xF=0.5:

위 RPE 알고리듬에서 공정변수 □, 다시 말하면 배열점에서의 분배계수 4개와탈거단에서의 liquid holdup을 예측함에 있어서 대각 Kalman 이득행렬(G), 공분산행렬(□), 예측오차의 자승으로 된 기준함수의 이차미분(R) 등의 초기값에 대한 민감도와 시간에 따른 feed 농도의 변화에 대한 수렴여부 등을 분석하였다. 첫번째는 feed의 농도 변화가 없는 경우이다. 우선 3개의 초기값(G(0),□(0),R(0))은 작을수록 실제값에 수렴하는 속도가 빠르지만 전반부에서 진동을 나타내는 반면 값이 증가할수록 진동은 약해져서 없어지고 수렴속도가 현저히 떨어지는 경향을 보여준다.  그림1은 예측값 중의 하나인 증류탑 하단에서의 분배계수에 대한 여러 초기값의 대각 Kalman 이득행렬의 영향을 보여주고 그림2는 공분산행렬의 영향을 보여준다.

그림1. Kalman 이득행렬의 변화에 따른 탑하부 분배계수의 수렴 경향

그림2. 공분산행렬의 변화에 따른 탑하부 분배계수의 수렴 경향

시간에 따른 feed 농도의 변화가 있을 때의 RPE 알고리듬에 의한 이 feed의 농도의 예측을 살펴보면 위의 농도변화가 없는 경우와 각 초기값의 경향은 그 크기만 다를 뿐 흡사함을 알 수 있다.  그림3은 기준함수의 이차미분(R)의 초기값에 따른 feed의 농도 예측의 영향을 보여주고 있으며 그림4는 condition number의 reset 값에 대한 feed의 농도 예측의 영향을 보여주고 있다.

그림3. 기준함수의 이차미분(R)의 초기값의 변화에 대한 xF 추정 경향

그림4. covariance resetting의 변화에 대한 xF 추정 경향

여기서 중요하게 살펴보아야 할 변수는 condition number이다. condition number의 reset값 역시 위에 언급했듯이 다른 변수와 같은 경향을 보이나 중간정도의 값에서는 다른 경향을 보인다. 즉 아주 낮은 값에서 시작했을 경우에 처음에는 진동이 많다가 서서히 줄어들면서 xF 의 추적이 느려지나 어느 한계값에 가까워지면 추정이 다시 빨라지게 된다. 그러나 이 한계값을 넘어서게 되면 추정이 오른쪽으로 이동된 모양을 보이면서 xF를 벗어나는 정도가 심해진다.

4.참고문헌

1) 김흥식, 1990, Ph.D. Thesis, 서강대학교

2) L.Ljung, T. S□derstr□m, Theory and Practice of Recirsive Identification, The MIT Press, 1983

3) Mohinder S. Grewal, Angus P. Andrews, Kalman Filtering, Prentice-Hall, 1993