EKF를 이용한 pH 공정의 공정변수예측

*이태철(학), 양대륙(정)

고려대학교 공과대학 화학공학과

Parameter Identification of pH Process Using EKF

Tae Chul Lee and Dae Ryook Yang

Department of Chemical Engineering, Korea University

서론

pH 공정은 비선형성이 심한 공정으로 모델링과 제어가 힘들다.  pH 공정의 동적 모사가 어려운 이유는 수용액상에서 산과 염기가 반응을 할 때 이 시스템은 분명히 하나의 화학 반응이 일어나는 반응기이다. 그러나 일반적으로 이러한 산염기 반응은 평형상태라 가정하기 때문에 반응상수를 무한대로 하고 반응기 모델을 얻을 수밖에 없다.  이러한 이유로 pH 공정의 동적 모델 개발은 어렵다. 이러한 경우 반응 속도식을 제외한 형태의 모델링이 시도되었으나 일반적 형태의 모델링이라고 볼 수는 없다.[1]

그림  pH neutralization process

또한 pH 공정의 동적 모델을 개발한 경우 이 모델의 비선형성이 매우 심하므로 이러한 모델을 이용한 공정 변수의 예측이 어려우므로 제어가 힘들다.

본 연구에서는 이러한 모델링의 어려움을 극복할 수 있는 reaction invariant를 이용한 모델을 통해 pH공정에서 공정변수를 RPEM(Recursive Prediction Error Method)를 이용해 추정해보았다.[1, 2]

본론

1. Reaction Invariant

Chemical reaction invariant란 시스템 내에서 화학반응에 의해 영향을 받지 않는 변수를 의미한다.  예를 들어 A→B 의 반응에서 cA+cB 변수는 화학반응에 의해 영향을 받지 않으며, 또한 cA+cB 변수를 이용해서 어떠한 화학반응이 일어났는지는 알 수 없다.  이러한 invariant는 일반적으로 시간이나 공간의 함수일 필요는 없다.[1]  Reaction invariant를 이용하여 화학반응식을 표현하는 것은 가능하며 이런 경우 invariant의 정보를 이용하여 화학반응식을 계산할 수 있으며, 화학반응식을 이용하여 invariant의 정보를 계산할 수도 있다.

이 reaction invariant를 이용하여 화학반응식을 구성하는 경우 새로운 변수를 이용하여 나타내는 것이 아니라 기존에 있는 반응식에 의한 제약조건을 나타내는 것이며, 이러한 reaction invariant를 이용하는 경우 수치계산을 매우 효과적으로 수행할 수 있다는 장점이 있다.

2. Dynamic Model for pH process[4]

본 연구에서 사용된 pH 중화장치는 그림 1과 같으며 각 stream에 대한 설명은 다음과 같다.  q1 : acid stream, q2 : buffer stream, q3 : base stream.  이 장치 내에서 일어나는 화학반응은 다음과 같다.

[eqn: rm H_2 CO_3 &<-> HCO_3^- + H^+ # rm HCO_3^- &<-> CO_3^2- + H^+ # rm H_2 O &<-> OH^- + H^+]

반응평형상수는 다음과 같이 주어진다.

[eqn: K_a1 = rm {[HCO_3^- ] [H^+ ]} over {[H_2 CO_3 ]}, it K_a2 = rm {[CO_3^2- ] [H^+ ]} over {[HCO_3^- ]}, it K_w = rm [H^+ ][OH^- ]]

Reaction invariant를 이용한 pH 공정의 화학적 평형을 구성한다.  본 시스템의 경우 두 개의 invariant를 각각의 입력 흐름에 대해 정의한다.

[eqn: W_ai &= rm [H^+ ]_i - [OH^- ]_i - [HCO_3^- ]_i - 2[CO_3^2- ]_i # W_bi &= rm [H_2 CO_3 ]_i + [HCO_3^- ]_i + [CO_3^2- ]_i]

Invariant Wa는 전하와 관련된 항이며, Wb는 carbonate ion과 관련된 항이다.

이러한 invariant를 이용한 pH 공정의 동특성 방정식은 다음과 같다.

[eqn: V {dW_a4} over { dt} &= q_1 (W_a1 - W_a4 ) + q_2 (W_a2 - W_a4 ) + q_3 (W_a3 - W_a4 ) # V {dW_b4} over { dt} &= q_1 (W_b1 - W_b4 ) + q_2 (W_b2 - W_b4 ) + q_3 (W_b3 - W_b4 ) #]

위의 식에서 다음과 같은 상태변수와 입력 출력을 정의하여 비선형 상태방정식을 구성할 수 있다.

[eqn: {bold x} = [W_a4 W_b4 ]^T , u=q_3 , y=pH]

[eqn: {bold {dot x}} = - 1 over V bmatrix{ {q_1 + q_2 } & 0 # 0 &{q_1 + q_2} } {bold x} - 1 over V bmatrix{ W_a3 - x_1 # W_b3 - x_2 } u + bmatrix{q_1 W_a1 + q_2 W_a2 # q_1 W_b1 + q_2 W_b2 }]

[eqn: c(x,y) = x_1 + 10^{y-14} - 10^-y + x_2 {1+2 times 10^{y-pK2}} over { 1 + 10^{pK1-y} + 10^{y-pK2}}]

3. RPEM

실제 공정을 모델을 이용해 추정하는데는 한계가 있고 또한 측정할 수 없는 공정 파라미터 있는 경우 파라미터의 추정을 위해  인식이라는 과정을 거쳐야 한다.  이러한 모델인식을 도입하여 제어하는 경우 모델인식을 하지 않는 경우보다 더 나은 성능을 보여준다.  이 경우 측정되는 출력변수의 값을 계속적으로 이용하여 주어진 모델의 변수를 향상시켜 가는 방법이 대부분이다.  가장 대표적인 방법으로는 Kalman filtering이 있고 모델이 비선형인 경우에는 Extended Kalman filtering을 이용한다.  이러한 EKF도 또한 공정잡음에 대한 정확한 정보가 없는 경우 추정변수가 발산하거나 bias되는 경우가 생기게 된다.  이러한 경우를 방지하기 위해 RPEM(Recursive Prediction Error Method)을 사용한다.[2]

RPEM은 innovation model을 이용하여 다음의 회귀적 예측오차 방법의 알고리즘을 구성할 수 있다.

[eqn: epsilon (t)= y(t) - hat y (t)]

[eqn: hat LAMBDA (t) = hat LAMBDA (t-1) + gamma(t) [ epsilon(t) epsilon^T (t) - hat LAMBDA (t-1)]]

[eqn: R(t) = `R(t-1) + gamma(t) [ PSI(t){hat LAMBDA}^-1 (t) PSI^T (t) - R(t-1)]]

[eqn: hat theta (t)= hat theta (t-1) + gamma(t) R^-1 (t) PSI(t) {hat LAMBDA} ^-1 (t) epsilon(t)]

[eqn: W(t+1) = [A_t - G_t C_t ]W(t) + {bar M} _t - G_t D_t]

[eqn: PSI(t+1) = W ^T (t+1) C_t ^T + D ^T (hat theta (t) , hat x (t+1))]

[eqn: hat x (t+1) &= A sub t hat x (t) + B _ t u(t) + G _t epsilon(t) # hat y (t+1) &= C _t hat x (t+1)]

[eqn: bar M (theta,x,u,epsilon) == partial over {partial theta} LEFT [ A(theta)x+B(theta)u+G(theta)epsilon RIGHT ] vert_{theta= hat theta }]

[eqn: D(theta, x)= partial over {partial theta} [C(theta)x] vert_{theta=hat theta}]

본 연구에서는 상태방정식과 출력식의 비선형성이 매우 강하므로 위의 방법을 그대로 사용할 수 없다.  예측출력은 equation solver를 이용하여 계산하며 비선형 모델을 선형화하여 사용하지 않고 비선형의 상태로 이용하여 알고리즘을 구성하여 사용하였다.

그림  공정 변수의 예측

결론

Table 1에 주어진 데이터를 이용하고 Wa2와 Wb2를 공정파라미터라 선정한 후 RPEM을 수행하였다.  공정을 PRBS를 이용하여 excitation시킨 후 RPEM을 적용하여 35분 동안 5초의 샘플링 시간을 이용하여 420개의 데이터를 이용하였다.  RPEM을 수행하여 모델의 파라미터 및 상태변수, 출력변수에 대한 예측치에 대한 결과는 그림 2와 3에 주어져 있다.  결과에서 보듯이 RPEM을 이용하여 파라미터 인식을 뛰어나게 잘하고 있음을 보여주고 있다.

현재 연구는 time invariant 공정 파라미터에 대한 파라미터 인식을 적용하였으나 추후 연구과제는 time varying 공정파라미터에 대한 인식이라 할 수 있다.

그림  상태변수 및 출력변수의 예측값과 실제값의 비교

Table  Nominal plant parameters

Parameters	Nominal values	Parameters	Nominal values
[eqn: V]	[eqn: 2900 ml]	[eqn: q_3]	[eqn: 15.8 ~rm [ml/s]]
[eqn: K_a_1]	[eqn: 4.47 times 10^-7]	[eqn: pH_4]	[eqn: 7.00]
[eqn: K_a_2]	[eqn: 5.62 times 10^-11]	[eqn: W_a_1]	[eqn: 0.003 ~rm [M]]
[eqn: [q_1 ]]	[eqn: 0.003 ~rm [M] ~HNO_3]	[eqn: W_b_1]	[eqn: 5.00 times 10^-5 rm ~[M]]
	[eqn: 5.0 times 10^-5 rm ~[M] ~H_2 CO_3]	[eqn: W_a_2]	[eqn: -0.03 ~rm [M]]
[eqn: [q_2 ]]	[eqn: 0.03 ~rm [M] ~NaHCO_3]	[eqn: W_b_2]	[eqn: 0.03 rm ~[M]]
[eqn: [q_3 ]]	[eqn: 0.003 ~rm [M] ~NaOH]	[eqn: W_a_3]	[eqn: -3.05 times 10^-3 ~rm [M]]
	[eqn: 5.0 times 10^-5 rm ~[M] ~NaHCO_3]	[eqn: W_b_3]	[eqn: 5.00 times 10^-5 rm ~[M]]
[eqn: q_1]	[eqn: 16.6 ~rm [ml/s]]	[eqn: W_a_4]	[eqn: -4.50 times 10^-4 ~rm [M]]
[eqn: q_2]	[eqn: 0.55 ~rm [ml/s]]	[eqn: W_b_4]	[eqn: 5.50 times 10^-4 ~rm [M]]

참고문헌

1.	Waller, K. V. and Makila, P. M. : Ind. Eng. Chem. Process Des. Dev., 20, 1 (1981)

2.	L. Ljung : "Theory and Practice of Recursive Identification", the MIT Press, 1983

3.	Gaustafsson, T. F. and Waller, K. V. :  Chem Eng. Sci., 38, 389 (1992)

4.	Heson, M. A.  and Seborg, D. E.  : Int. J Adapative Control and Signal Processing, 11, 171 (1997)