비선형 ARMA 모델을 이용한 연속식 PS 중합 반응기의 제어

오창근, 나상섭, 이현구

서울대학교 응용화학부·화학공정신기술연구소

Control of a Continuous Styrene Polymerization Reactor

Using Nonlinear ARMA Model

Chang-Keun Oh, Sang-Seop Na and Hyun-Ku Rhee

School of Chemical Engineering & Institute of Chemical Processes,

Seoul National University

서 론

중합 반응계는 복잡한 반응기구와 심한 비선형성 때문에 선형 모델에 기초를 둔 제어기로는 만족스러운 제어결과를 얻기 어렵다. 신경회로망 같은 비선형 모델을 사용할 경우, 추정해야 할 파라메터가 많고 모델 연구의 어려움으로 구현이 쉽지 않은데 비하여 비선형 ARMA 모델은 입출력 정보로부터 제어 대상계를 용이하게 식별할 수 있는 장점이 있다[3]. 또한 steepest descent method는 초기 수렴속도는 느리지만, 가속이 가능하며 안정적인 최적화 알고리즘이다[5].

본 연구에서는 연속식 PS (Polystyrene) 중합반응계의 수학적 모델식들을 유도하여 이들을 플랜트로 가정한 뒤, 이로부터  자켓 온도와 전환율이 각각 입력과 출력이 되는 비선형 ARMA 모델을 구성하였다. 그리고 이 모델을 바탕으로 steepest descent method를 이용해서 최적화를 수행하는 모델 예측 제어기를 설계하고, 정상상태에서 전환율을 제어하는 수치모사를 병행하여 비선형 ARMA 모델의 타당성을 검증하였다.

반응기구와 모델식

[eqn: 개시반응 ~:~ I ~buildrel rarrow { k_d }~ 2 varphi CDOT ~~~~~~~~# varphi CDOT ~ + ~ M ~ buildrel rarrow {k_i} ~R_1 CDOT]

[eqn: 성장반응~ :~ R_1 CDOT ~+~ M ~ buildrel rarrow { k_p} ~ R_2 CDOT ~~~# R_j CDOT ~+~ M ~ buildrel rarrow { k_p } ~ R_j+1 CDOT ~]

[eqn: 종결반응 ~:~ R_i CDOT ~+~ R_j CDOT ~ buildrel rarrow {k_tc } ~ P_i+j ~ (결합)]

[eqn: 연쇄이동반응 ~:~ R_i CDOT ~+~ M ~ buildrel rarrow {k_trm }~ P_i ~+~ R_1 CDOT # R_i CDOT ~+~ S ~ buildrel rarrow {k_trs } ~ P_i ~+~ S CDOT ~~]

위와 같이 고분자 중합 반응 기구는 개시반응, 성장반응, 종결반응 및 연쇄이동반응으로 구성되며, styrene 중합 반응의 경우에는 결합에 의한 종결 반응만이 존재한다. 이를 바탕으로 연속식 중합 반응기내의 개시제, 단량체 및 용매에 대한 물질 수지식, 활성고분자 농도와 비활성 고분자 농도 각각의 0차, 1차 및 2차 모멘트에 대한 수지식, 그리고 에너지 수지식을 다음과 같이 얻을 수 있다.

[eqn: d(IV) over dt ~=~ I_f q_f - I q + (-k_d ) IV~ ,~~~~~~~~` d(SV) over dt ~=~ S_f q_f - S q +(-k_trs S G_0 ) V` # d(MV) over dt ~=~ M_f q_f - M q + left { -2fk_d I -(k_p + k_trm ) MG_0 right } V`~~~~~~~~~~~~~~~~~~~~~~~~~~]

[eqn: { d ( G_0 V ) } over dt &=~ -qG_0 + ( 2`f``k_d `I - k_tc `G_0 ^2 )V # { d ( G_1 V ) } over dt &=~ -qG_1 + left { 2`f`k_d`I + k_p`MG_0` -k_tc`G_0 G_1 + (k_trm M + k_trs S)(G_0 - G_1 ) right } V# { d (G_2 V) } over dt &=~ -qG_2 + left { 2`f``k_d `I + k_p M ( G sub 0 + 2 G sub 1 ) - k_tc G sub 0 G sub 2 + ( k_trm M + k_trs S)(G sub 0 - G sub 2 ) right } V`]    (1)    [eqn: { d (F_0 V ) } over dt &=~ -qF_0 + left { 0.5 k_tc G sub 0 ^2 + (k_trm M + k_trs S ) G sub 0 RIGHT }V # { d (F_1 V ) } over dt &=~ -qF_1 + left { k_tc G_0 G_1 + ( k_trm M + k_trs S ) G_1 right } V # { d (F_2 V ) } over dt &=~ -qF_2 + left { k_tc ( G_0 G_2 + G_1 ^2 ) + (k_trm M + k_trs S ) G_2 right } V`]

[eqn: {d ( ( rho C_p ) _m V T )} over dt ~=~ ( rho C_p )_m q_f T_f - (rho C_p )_m q T + ( - Delta H_p ) k_p M G_0 V - UA ( T - T_j )]          (2)여기서 G0, G1, 및 G2는 각각 활성 고분자 농도의 0차, 1차, 및 2차 모멘트, F0, F1, 및 F2는 각각 비활성 고분자 농도의 0차, 1차, 및 2차 모멘트이고, □Hp는 성장반응의 반응열, U는 총괄 열전달 계수, A는 열전달 면적을 나타내며, qf는 공급물의 유량, q는 배출물의 유량을 의미한다.

본 연구에서는 젤 효과를 고려하기 위해 Hamer등(1981)이 사용한 다음과 같은 종결 반응의 속도상수식을 사용하였다[2].

[eqn: & g(X,T ) ≡ k_t over k_{t`` 0} = exp [ -2(Ax+Bx^2 + Cx^3 ) ] # & A = 2.57 - 5.05×10^-3 T , ~ B = 9.56 - 1.76×10^-2 T, ~C = -3.03 + 7.85×10^-3 T # & x=X(1-f_s )]         (3)

여기서 X는 단량체의 전환율, x는 부피 분율이 fs인 용매가 존재할 때의 단량체 전환율을 나타낸다.

본 연구에서 사용한 속도상수들은 문헌에 보고된 값들을 인용하여 사용하였으며 그 값들을 Table 1에 수록하였다.

Table 1. Parameter values for rate constants

rate constant	freqency factor	activation energy(cal/mol)
kd(l/s)	5.95×[eqn: 10^13`]	29,600
kp(l/mols)	1.06×[eqn: 10^7`]	7,068
kto(l/mols)	1.25×[eqn: 10^9`]	1,675
ktrs(l/mols)	3.148×[eqn: 10^9`]	16,264
ktrm(l/mols)	1.186×[eqn: 10^7`]	11,767

시스템 식별

비선형 ARMA 모델의 형태는 다음과 같이 정의된다.

[eqn: y(k) = & theta_0 + SUM from { i=1} to { n_y } theta_1,i `y(k-i ) + SUM from { i=rho } to {n_u } theta_2,i `u(k-i) + SUM from { i=1} to {n_y } SUM from { j=1} to i theta_3,(i,j) y(k-i) y(k-j) # & + SUM from { i=rho} to {n_y } SUM from { j=rho} to i theta_4,(i,j) u(k-i) u(k-j) + SUM from { i=1} to {n_y } SUM from { j=rho } to { n_u } theta_5,(i,j) y(k-i) u(k-j) + ...]

여기서 ny는 출력의 차수, nu는 입력의 차수, □는 입·출력간의 시간 지연이고, □는 추정해야 할 파라메터들이다.

본 연구에서는 Table 2에 주어진 반응조건하에서 (1), (2)의 방정식을 실제 플랜트로 가정하였다. 이 때 샘플링 시간은 1분이고, 입력 u는 무차원화된 자켓의 입구온도( u = (T-333.15)/(353.15-333.15) )이며, 출력 y는 전환율이다.

PS 중합반응계의 식별을 위하여 사용한 입·출력 데이터는 Fig. 1과 같다. 입력은 최소 switching time이 20인 변형된 Random Binary Signal (RBS)이며, 이에

Table 2. Reaction Conditions

efficiency ofinitiator	weight of initiator	fraction of solvent	feed temp.	feed flow rate
0.6	10g/l	0.5	293.15K	10ml

대한 플랜트와 ARMA 모델의 출력을 비교하였는데 Percentage Variance Accounted For (VAF)값이 98.3으로 잘 일치하는 결과를 얻었다.

사용된 ARMA 모델의 구조는 ny=3, nu=3으로 하여 전개하였으며, stepwise algorithm을 이용하여 출력과 상관관계가 큰 8개의 변수를 선택하였다[1].

[eqn: y(k) =& theta_0 + theta_1 y(k-1) + theta_2 y(k-2) + theta_3 y(k-3) + theta_4 u(k-2) + theta_5 u(k-3) # & + theta_6 u(k-1)u(k-2) + theta_7 y(k-2) u(k-1) + theta_8 y(k-3)u(k-1)]             (4)

파라메터들은 RLS방법을 이용해서 추정하였으며, forgetting factor □=0.983으로 하여 계산하였다. 계산된 결과는 다음과 같다.

□0 = 0.0008	□1 = 0.3286	□2 = 0.3286	□3 = 0.3286	□4 = -0.0003
□5 = -0.0003	□6 = 0.0013	□7 = 0.0023	□8 = 0.0023

제어 알고리즘

본 연구에서는 목적 함수를 다음과 같이 정의하여 모델 예측 제어 알고리즘을 적용하였다.

[eqn: E = 1 over 2 LEFT ( Gamma_y SUM from { {l}=1} to {p} [ y_pred (k+l|k) - y_sp (k+l)]^2 + Gamma_u SUM from { {l} = 0 } to {m-1} [ u(k+l) - u(k+l-1) ]^2 RIGHT )]            (5)

여기서 □y와 □u는 각각 출력과 입력의 가중치이고, p는 예측범위이며, m은 제어범위이다.

이 때 u(k), ..., u(k+m-1)이 서로 독립적이고, u(k+m-1)=u(k+m)=... 라고 가정하면, steepest descent method를 이용하여 다음과 같이 u(k), ..., u(k+m-1)을 반복계산을 통하여 얻을 수 있다.

[eqn: & {bar u}(k) _ i+1 = {bar u}(k)_i - eta {partial`E} over {partial {bar u}(k)_i } + alpha( {bar u } (k) _ i - {bar u} (k) _i-1 ) ~~~~~~ until `````E <= epsilon ## & {partial E} over {partial {bar u}(k)_i } = Gamma_y ``J_y left{ {bar y } (k+1)_i - {bar y_sp} (k+1) right} + Gamma_u left{ 2{bar u}(k) _i - {bar u}(k-1) _i - {bar u}(k+1) _i right} ## & { bar u} (k) _i = { [ u(k)`` u(k+1)`` `...` ``u(k+p-1) ]^T} _i ## & { bar y } (k+1) _ i = { [ y(k+1)`` y(k+2)`` `...` ``y(k+p) ]^T} _i]         (6)

여기서 □는 step size, □는 moment constant이고, Jy는 [eqn: it {bar y} (k+1) _i]의 Jacobian을 의미한다. E가 적절한 오차인 □에 도달하거나 반복계산횟수가 미리 정해진 한계값에 이르면, 계산을 중단하고 얻어진 u(k)를 플랜트에 적용하였다. 단, u(k)가 하한이나 상한에 도달하면 상한값이나 하한값을 u(k)로 바꾸어 적용하였다[4].

수치 모사

본 수치모사에서 사용한 제어파라메터 값들은  다음과 같다.

□(초기값)=0.01,		□의 증가속도=1.05,		□의 감소속도=0.7
□=0.9,	□=0.00001 (또는 최대 반복 계산횟수 = 100)
umin = 0,			umax = 1,		□umax = 0.2

목표값을 0.1로 유지하다가 0.12로 계단변화를 시켰을 때의 제어결과를 Fig. 2에 나타내었다. 입력에는 편차가 0.02, 출력에는 편차가 0.0008인 white noise를 주고, m=20, p=7로 하여 수치모사를 수행하였다. Fig. 2에서 관찰되는 바와 같이   제어기가 전환율을 비교적 정확히 제어하고 있음을 알 수 있다.

Fig. 1. I/O data for the identification of continuous PS polymerization system.	Fig. 2. Simulation result of the reactor control for a setpoint change (m=20, p=7, □y=0.9, □u=0.1).

결 론

본 연구에서는 연속식 PS 중합반응계에 대하여 기존에 알려져 있는 수학적 모델식을 실제 플랜트로 가정하고, 변형된 RBS로 얻은 입출력 자료를 이용해 SISO 비선형 ARMA모델을 구성하였다. 그리고, 이 비선형 ARMA 모델을 근거로 steepest descent method를 이용해서 최적화를 하는 모델 예측 제어 알고리즘을 구성하고, 임의의 목표값에 대하여 수치모사를 수행하여 모델의 타당성을 검증하였다. 이러한 제어방법은 구현이 용이한  비선형 ARMA 모델을 도입하여 간단한 steepest descent method만으로 최적화를 하여 제어를 할 수 있다는 장점을 지니고 있다. 앞으로는 위의 제어방법을 MIMO 공정에 대해 확장하는 방법을 연구하고, 실제 실험을 통하여 검증할 예정이다.

감 사

본 연구는 학술진흥재단과 포항공과대학교 공정산업의 지능자동화 연구센터(ARC)를 통한 한국과학재단 우수연구센터의 연구비지원에 의하여 수행되었으며, 이에 대하여 심심한 감사의 뜻을 표합니다.

참고 문헌

1. Draper, N. R. and H. Smith, Applied Regression Analysis, Wiley, New York, pp. 163-195(1966).

2. Hamer, J. W., T. A. Akarmov, and W. H. Ray, "The dynamic behavior of continuous polymerization reactors - II. Nonisothermal solution homopolymerization  and copolymerization in a CSTR", Chem. Eng. Sci., 36, pp. 1897-1914(1981).

3. Hern□ndez, E. and Y. Arkun, "Control of Nonlinear Systems Using Polynomial ARMA Models", AIChE J., 39(3), pp. 446-460(1993).

4. Zhan, J. and M. Ishida, "The Multi-step Predictive Control of Nonlinear SISO Processes with a Neural Model Predictive Control Method", Computers & Chem. Engng, 21(2), pp. 201-210(1997).

5. Na, S.-S., C.-K. Oh and H.-K. Rhee, □Identification of Continuous Styrene Polymerization Reactor System using Nonlinear ARMA Model□, Theories and Applications of Chemical Engineering, 4(1), pp. 349-352(1998).