연속식 스티렌 중합 반응기의 비선형 모델예측제어

나상섭, 이현구

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

Nonlinear Model Predictive Control of a Continuous Styrene Polymerization Reactor

Sang-Seop Na and Hyun-Ku Rhee

School of Chemical Engineering & Institute of Chemical Processes,

Seoul National University

서 론

본 연구에서는 연속식 스티렌 중합 반응기의 전환율과 분자량을 제어하기 위하여 다항식 ARMA(auto-regressive moving average)모델에 기반한 비선형 모델예측제어(NLMPC) 알고리듬을 구현하고 제어기의 타당성을 검증하기 위하여 수치모사와 실험을 수행하였다. ARMA모델의 경우 공정의 비선형성이 다항식 형태로 표현되는 장점이 있으나 모델예측제어기에 적용하였을 때 목적함수를 최소화하는 최적화 단계에서 입력에 대하여 매우 높은 차수의 함수를 풀어야 하는 단점이 있다. 따라서 온라인 제어시 많은 계산량이 요구되는데 본 연구에서는 이를 개선하기 위하여 매 샘플링 시간마다 추정된 상태에 대해서 모델을 선형화 하는 기법을 적용하여 효과적인 제어가 이루어지도록 하였다. 또한 연속식 스티렌 중합 반응기에 대해서 온라인 밀도계와 점도계로부터 전환율과 분자량을 계산할 수 있는 상관식을 제안하고 제어에 이용하는 방안을 제시하고자 한다.

비선형 모델예측제어 알고리듬

MIMO(multiple input multiple output)계에 대한 다항식 ARMA모델의 정의는 다음과 같다.

[eqn: y_i `(k`)``=&``y_{i``0} ``+``SUM from {l=1} to {q_y }`SUM from {j = 1} to {n_y }`` { theta}`_{l`,`j } ^{1`,`i} ``y_l ``(k`-`j`) ``+``SUM from {l=1} to {q_u} `SUM from {j=1} to { n_u }``{ theta}`_{l`,`j } ^{2`,`i} ``u_l ``(k`-`j`)# &``+``SUM from {l=1} to {q_y} `SUM from {j=1} to { n_y }``SUM from {m=l} to {q_y} `SUM from {n=1} to { j }``{ theta}`_{l`,`j`,`m`,`n } ^{3`,`i} ``y_l ``(k`-`j`)``y_m``(k`-`n`)``]

(1)

[eqn: &``+``SUM from {l=1} to {q_u} `SUM from {j=1} to { n_u }``SUM from {m=l} to {q_u} `SUM from {n=1} to { j }``{ theta}`_{l`,`j`,`m`,`n } ^{4`,`i} ``u_l ``(k`-`j`)``u_m``(k`-`n`)# &``+``SUM from {l=1} to {q_y} `SUM from {j=1} to { n_y }``SUM from {m=1} to {q_u} `SUM from {n=1} to { n_u }``{ theta}`_{l`,`j`,`m`,`n } ^{5`,`i} ``y_l ``(k`-`j`)``u_m``(k`-`n`)``+`` CDOTS``]

식(1)에서 y , u는 각각 출력과 입력을 나타내며 qy는 출력의 개수, qu는 입력의 개수를 의미한다. 그리고 ny와 nu는 각각 출력과 입력의 시간지연 차수를 의미하는데 본 연구에서는 3차항 까지 고려하였다. 식(1)을 2×2 MIMO계에 대하여 전개한 뒤 Hernandez와 Arkun[1]이 제안한 delay coordinate method에 의하여 state-space realization을 하면 다음과 같이 유도된다.

[eqn: x_1 (k)``=``y_1 `(k) ~,~x_2 (k)``=``y_2 `(k)~,~ &x_3 (k)``=``y_1 `(k-1) ~,~x_4 (k)``=``y_2 `(k-1)# CDOTS~,~&x_13 (k)``=``u_1 `(k-3) ~,~x_14 (k)``=``u_2 `(k-3)``]

(2)

[eqn: x `(k+1)``=``A CDOT x `(k) + B``]

(3)

[eqn: y (k) ``=`` BMATRIX { {1}& {0}& {CDOTS}& {CDOTS}& {CDOTS}& {0}# {0}& {1}& {CDOTS}& {CDOTS}& {CDOTS}& {0} } CDOT x (k)``]

(4)

식(3)에서 행렬 A와 B는 다음과 같이 정의되는데 이를 간략히 표현하면 식(5)와 같다.

[eqn: &x `(k+1`) ``=``F`[x `(k`)~,~u `(k`)`]# &y `(k+1`) ``=``h`[x `(k`) `] ``]           [eqn: A``=`` BMATRIX { {0}& { }& { }& {CDOTS}& {0}& { }& { }& { }# {0}& { }& { }& {CDOTS}& {0}& { }& { }& { }# {1}& { }& { }& {CDOTS}& {0}& { }& { }& { }# {0}& {1}& { }& {CDOTS}& {0}& { }& { }& { }# { }& { }& { }& { }& {DOTSVERT}& { }& { }& { }# { }& { }& { }& { }& {DOTSVERT}& { }& { }& { }# { }& { }& { }& { }& { }& {0}& {CDOTS}& {0}# { }& { }& { }& { }& { }& {0}& {CDOTS}& {0}# { }& { }& { }& { }& { }& {1}& {CDOTS}& {0}# { }& { }& { }& { }& { }& { }& { }& {DOTSVERT}# }] , [eqn: B``=``BMATRIX { {p_1 [ x (k) `,`u (k)]}# {p_2 [ x (k) `,`u (k)]}# {0}# {DOTSVERT}# {0}# u_1 (k)# u_2 (k)# {0}# {DOTSVERT}# {0} }]

(5)

행렬 B에서 p1과 p2는 각각의 출력에 대해서 식별을 통하여 얻어지는  다항식을 의미한다.

외란의 효과를 고려하고 steady-state offset을 제거하기 위하여 처음 두 개의 state(x1(k),x2(k))에만 process noise가 가해진다고 가정하고 disturbance state variable[eqn: {x}^w]을 도입하면 식(6)과 같은 augmented state model이 얻어진다.

[eqn: BMATRIX { {x `(k) }## {{x}^w (k)}``} ``&=``BMATRIX { { F ``( x `(k-1)``,``u `(k-1)`) }## {{x}^w (k-1)}``} ``=`` BMATRIX { {G}& {0}## { I}& {0 }``} CDOT w `(k -1) ## y `(k) ``&=``h ``(x `(k-1)``) ``+``{x}^w `(k-1) ``+``nu `(k)]식(6)에서 w는 process white noise이며 □는 measurement white noise이다.

(6)

식(6)의 state model에 대해서 EKF기법을 도입하면 다음과 같은 implementation equation과 linear approximation equation이 유도된다.

nonlinear implementation equation

[eqn: BMATRIX { { hat x `(k/k-1`)}## {{hat x}^w `(k/k-1`)}``} ``&=`` BMATRIX { F `(`{hat x `(k-1/k-1`)``,`` u`(k-1`)``)}## {{ hat x}^w `(k-1/k-1`)}``}## hat y `(k/k`) ``&=``h `( hat x `(k/k-1`)``)``+``{hat x}^w `(k/k-1`)``]	linear approximation equation

(7)

[eqn: BMATRIX { {hat x `(k/k`)}## {{hat x}^w `(k/k`)}``} ``=`` BMATRIX { {hat x `(k/k-1`)}## {{hat x}^w `(k/k-1`)}``}``+``L`(k`) ``(y_m `(k) - hat y `(k/k`)``)``]식(8)에서 L(k)는 estimation gain이며 ym은 측정값을 의미한다. 본 연구에서는 목적함수를 다음과 같이 정의하였다.

(8)

식(9)에서 p와 m은 각각 prediction horizon[eqn: J ``=`` &SUM from { {i}=1} to P LEFT { LAMBDA _y ``LEFT [ ``y (k+i`) ``-``r`(k+i`) RIGHT ] RIGHT }^2 ``+``SUM from {{j}=1} to m LEFT { LAMBDA_u ``DELTA u `(k+j-1`)``RIGHT}^2``]과 control horizon을 나타내며 □y, □u는 각각 출력과 입력에 대한 가중치를 의미한다. 그리고 r은 목표치를 나타낸다.

(9)

output prediction equation은 state model과 estimation equation을 현재의 state와 입력에 대해서 선형화 하여 구하였으며 이를 식(9)에 대입하면 quadratic형태로 유도된다. 본 연구에서는 구현된 모델예측제어 알고리듬을 연속식 스티렌 중합 반응기의 제어에 적용하기 위하여 입력을 냉각수의 유입온도와 공급물의 유량으로 하고 출력을 단량체의 전환율과 중량 평균 분자량으로 정의하였다. 수치모사에서는 다음과 같은 공정의 제약조건이 있다고 가정하였다.

● jacket inlet temperature(u1) range : [eqn: 55℃`` <= `` T_jin`` <= ``85℃]

● feed flow rate(u2) range : [eqn: 5□/min`` <= `` q_f`` <= ``30□/min``]

● conversion(y1) range : [eqn: 0`` <= `` X`` <= ``0.5``]

● molecular weight(y2) range : [eqn: 10000`` <= `` Mw`` <= ``50000``]

● [eqn: LEFT | DELTA T_jin ``RIGHT |`` <= ``5℃``] , [eqn: LEFT | DELTA q_f ``RIGHT |`` <= ``2□/min``]

시스템 식별

Table 1과 같은 반응조건에 대해서 연속식 스티렌 중합 반응기를 식별하는 실험을 수행하였다. 실험에서는 확률이 0.1인 pseudo random multi-level input signal을 사용하였으며 중합 반응기의 시간지연이 매우 크기 때문에 입력신호가 한 번 결정되면 최소 40분을 유지하도록 하였다. 온라인 밀도계와 점도계를 설치하고 이들 자료로부터 전환율과 분자량을 계산하기 위하여 상관식을 이용하였다. 전환율의 경우 물질 수지를 고려하여 식(10)을 유도할 수 있다.

[eqn: X``&=``LEFT{ LEFT ( 1 over rho_d `-` 1 over rho_md RIGHT ) `-`W_sf ` LEFT ( 1 over rho_sd `-`1 over rho_md RIGHT ) RIGHT}/ { LEFT ( 1 over rho_pd `-` 1 over rho_md RIGHT ) }]

(10)

식(10)에서 □d가 온라인 밀도계에서 측정되는 반응물의 밀도이며 하첨자 m, s, p는 각각 단량체, 용매, 고분자를 의미한다. 그리고 Wsf는 용매의 무게분율을 나타낸다. 분자량의 경우 intrinsic viscosity와 specific viscosity간의 관계는 Lyon's equation을 적용하였으며 이 식으로부터 intrinsic viscosity를 계산하고 Mark-Houwink equation을 적용하여 분자량을 계산하였다. Lyon's equation은 식(11)과 같다.

[eqn: eta_sp over { C ``[eta] } ``=``exp`(`{k_H `[eta] `C} over {1`-`b`C }`)``]

(11)

Table 1의 조건하에서 실험을 실시하여 얻은 식별의 결과를 Figure 1에 나타내었다. 전환율과 분자량은 20분마다 시료를 채취하여 각각 무게측정법과 GPC(Gel Permeation Chromatography)를 통하여 분석한 결과와 식(10)및 식(11)로부터 온라인으로 계산된 결과들을 비교하여 나타내었는데 전 영역에서 매우 잘 일치함을 알 수 있다. 또한 수치모사 결과와 식별된 ARMA모델로부터 추정된 결과도 실제 플랜트의 출력과 거의 유사함을 알 수 있다. ARMA모델의 식별은 온라인으로 계산된 출력값들을 사용하였으며 아래와 같이 결정되었다.

Figure 1. The experimental result of the identification for the continuous Styrene polymerization reactor.

y1(k) = 0.01 +0.49y1(k-1) -0.0015u12(k-2) +0.32y1(k-2)

-0.14y1(k-1)u2(k-2) +0.07y1(k-3) +0.02u1(k-3)u1(k-1)

-0.008y2(k-3)u1(k-1) -0.01y2(k-3)

y2(k) = -0.002 +0.58y2(k-1) +0.37y2(k-2) -0.15y2(k-3)u1(k-2)

-0.03u2(k-2) +2.19y1(k-1)y2(k-1) -1.17y1(k-3)y2(k-1)

+0.06y2(k-1)u1(k-2) -0.02u1(k-3)u2(k-2)

-0.07y2(k-3)u1(k-1)

Table 1 Reference Conditions for the continuous Styrene Polymerization used in  experiment.

Items	Values
Initial Charge Monomer(Styrene) Solvent(Toluene) Initiator(AIBN)	400□ 400□ 8.0g
Feed Conditions Monomer : Solvent (volume ratio) Initiator	1:1 10g/□
Operating Conditions Reactor Temperature Feed flow rate	50 ∼ 85℃ 5 ∼ 30□/min

수치 모사

식별된 다항식 ARMA모델을 이용하여 전환율과 분자량을 제어하는 수치모사를 수행하였다. Sampling time은 2분이며 prediction horizon과 control horizon은 각각 15와 3으로 하였다. Figure 2는 전환율은 0.1로 일정하게 유지시키고 분자량의 목표치를 25000→20000→25000으로 변경했을 때의 제어결과이며 Figure 3은 전환율과 분자량의 목표치를 각각 0.1→0.1→0.15와 25000→20000→20000으로 변경했을 때의 제어 결과를 나타낸다. 제어기 출력이 공정의 제약조건을 만족하면서 빠르게 수렴되는 것을 알 수 있고 두 가지 경우 모두 목표치를 잘 추적하고 있음을 알 수 있다.

Figure 2. The NLMPC simulation result for the case in which only the setpoint of Mw is changed.

Figure 3. The NLMPC simulation result of the case in which setpoints of X and Mw are changed.

결 론

본 연구에서는 다항식 ARMA모델을 이용하여 연속식 스티렌 중합 반응기를 실험을 통하여 식별하였는데 강한 비선형성에도 불구하고 정밀하게 식별할 수  있었다. 그리고 매 샘플링 시간마다 모델을 선형화 하는 기법을 도입하여 온라인 제어에 효과적인 비선형 모델예측제어 알고리듬을 개발하여 이를 구현하고 수치모사를 통하여 성능을 검토하였는데 만족스러운 결과를 얻을 수 있었다. 이러한 제어방법은 실제 생산공정에 적용하기에 매우 효과적일 것으로 판단된다.

감 사

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

참고 문헌

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

2. J. H. Lee and N. L. Ricker, "Extended Kalman Filter Based Nonlinear Model Predictive Control", Ind. Eng. Chem. Res., 33, pp. 1530(1994).

3. B. R. Maner and F. J. Doyle Ⅲ, "Polymerization Reactor Control Using Autoregressive-Plus Volterra-Based MPC", AIChE J., 43(7), PP. 1763-1784(1997).

4. M.-J. Park, S.-S. Na, S.-M. Ahn and H.-K. Rhee, "Measurement Technique for the Physical Properties in Polymerization Reactors using on-line Densitometer and Viscometer", Proc. 13th KACC, 2, pp. 1967-1970, Pusan, Korea(Oct. 1998).