CONTROLE PREDITIVO NÃO-LINEAR PRÁTICO APLICADO A UM CSTR

ADRIANO S. M. BRANDÃO

Programa de Pós-graduação em Engenharia de Automação e Sistemas, Universidade Federal de Santa

Catarina
Florianópolis – Santa Catarina - Brasil
E-mail: asmbrandao@gmail.com

Resumo Figura 1. Representação simplificada do reator.

Este artigo apresenta um estudo da aplicação de um con-

trolador preditivo não-linear em duas diferentes varia-
ções no controle da concentração e temperatura em um
reator CSTR São apresentadas as formulações matemá-
ticas dos controladores utilizados e alguns detalhes rele-
vantes na implementação dos controladores. Por fim, são
apresentados resultados e análises do comportamento de
ambos controladores em diferentes cenários.

1. Introdução

O controle preditivo baseado em modelo (MPC) é ambiente desprezíveis. A Figura 1 apresenta um es-
uma técnica de controle avançado de sucesso na in- quema simplificado do reator, evidenciando as variá-
dústria, e aplicável em sistemas multivariáveis e com veis de interesse, enquanto o Quadro 1 apresenta um
restrições. Neste tipo de controle, é necessário um mo- descritivo de cada variável. O objetivo deste trabalho
delo que pode ser linear ou não-linear. Do ponto de é controlar a concentração de B (𝐶𝐵 ) e a temperatura
vista prático, existem várias dificuldades para a apli- (𝑇) na corrente de descarga, manipulando as vazões de
cação do MPC com modelos não-lineares (NMPC), carga do reator (Fin ) e da camisa (Fj,in ).
como a necessidade de obter a solução de um pro-
blema de otimização não convexo. O controle NMPC Quadro 1. Descrição das principais variáveis do processo.
é atualmente um campo aberto para pesquisas Variável Descrição
(CAMACHO; BORDONS, 2007). O estudo de pro-
𝐅𝐢𝐧 Vazão volumétrica na carga do reator
postas como o PNMPC de (PLUCENIO, 2010), como
o realizado em (YANG, 2016) é interessante para via- 𝐂𝐀,𝐢𝐧 Concentração de A na carga do reator
bilizar a aplicação destas técnicas em processos com- 𝐓𝐢𝐧 Temperatura na carga do reator
plexos como o CSTR. 𝐅𝐣,𝐢𝐧 Vazão volumétrica na carga da camisa
No meio industrial, o reator tipo tanque agitado 𝐓𝐣,𝐢𝐧 Temperatura na carga da camisa
(CSTR) é amplamente utilizado em processos que
𝐅𝒐𝒖𝒕 Vazão na descarga do reator
ocorrem de forma não isotérmica com reações. O
comportamento dinâmico deste tipo de reator é alta- 𝐂𝐀 Concentração de A no reator
mente não-linear, devido ao alto acoplamento entre as 𝑪𝑩 Concentração de B no reator
variáveis do sistema, tornando interessante a aplicação 𝑻 Temperatura no reator
de controladores preditivos não-lineares. 𝐓𝒋 Temperatura na camisa
Este trabalho consiste na simulação de duas variantes
do PNMPC aplicadas a um CSTR, bem como uma
comparação do desempenho dos controladores pro- 2.1 Modelo matemático
postos.
A partir das leis de conservação de massa e ener-
gia, foram obtidas as seguintes equações diferenciais
2. Descrição do processo para representar o CSTR:

O reator perfeitamente agitado (continuous stirred- Fin (𝑡) = F𝑜𝑢𝑡 (𝑡) (1)
tank reactor model, CSTR) é um tipo de reator ideal,
no qual assume-se que a mistura dos componentes 𝑑(𝐶𝐴 (𝑡)) CA,in (𝑡). Fin (𝑡) CA (𝑡). Fout (𝑡)
dentro do reator é perfeita, não existindo gradientes de = −
𝑑𝑡 Aa . ℎ Aa . ℎ (2)
temperatura ou composição. Neste processo, consi- −
𝐸
dera-se uma reação irreversível exotérmica simples e −𝐶𝐴 (𝑡). 𝑘0 . 𝑒 𝑅.𝑇(𝑡)

em fase líquida A→B, com cinética de primeira ordem 𝑑(𝐶𝐵 (𝑡)) C𝐵 (𝑡). Fout (𝑡)
e a presença de uma camisa térmica. O nível do reator =−
𝑑𝑡 Aa . ℎ (3)
é considerado constante e as perdas térmicas para o −
𝐸
+𝐶𝐴 (𝑡). 𝑘0 . 𝑒 𝑅.𝑇(𝑡)
 𝑑(𝑇(𝑡)) Fin (t). Tin (t) F𝑜𝑢𝑡 (t). T(t) Tabela 4. Restrições operacionais do reator.
= −
𝑑𝑡 Aa . ℎ Aa . ℎ Variável Mínimo Máximo Unidade
−
𝐸 (4) 𝑪𝒃 150 300 mol.m-3
𝑈. 𝐴. (𝑇(𝑡) − 𝑇𝑗 (𝑡)) 𝐶𝐴 (𝑡). Δ𝐻. 𝑘0 . 𝑒 𝑅.𝑇(𝑡)
− − 𝑻 314 320 K
𝐶𝑝 . ρ. Aa . ℎ ρ. 𝐶𝑝
𝑭𝒊𝒏 0 6e-4 m3.s-1
𝑭𝒋 0 10e-4 m3.s-1
𝑑 (𝑇𝑗 (𝑡)) Fj,in (t). (Tj,in (t) − T𝑗 (t)) ∆𝑭𝒊𝒏 -0.5e-4 0.5e-4 m3.s-1
= ∆𝑭𝒋 -0.5e-4 0.5e-4 m3.s-1
𝑑𝑡 Vj
(5)
𝑈. 𝐴. (𝑇(𝑡) − 𝑇𝑗 (𝑡))
+ 3. Controle proposto
𝐶𝑝,𝑗 . ρj . Vj

O MPC é uma família de controladores nos quais a

2.2 Parâmetros do modelo e de operação ação de controle é calculada a partir da minimização
Na Tabela 1 são apresentados os parâmetros utilizados de uma função objetivo que considera o erro no segui-
nas equações 1,2,3,4 e 5, enquanto que na Tabela 2 são mento de referência em uma projeção do comporta-
apresentadas as condições iniciais dos estados. Na Ta- mento do processo em um horizonte deslizante e os
bela 3, é mostrado o ponto de operação inicial a ser movimentos das variáveis manipuladas, com suas res-
considerado para a condição inicial aproximadamente pectivas ponderações. Uma das grandes vantagens do
no estado estacionário. Todos os parâmetros do mo- MPC, quando comparado a outros controladores, é
delo foram adaptados de (CORRIOU, 2004). Na Ta- que a abordagem do caso multivariado com restrições
bela 4, são apresentadas as restrições que devem ser é apenas uma extensão do caso SISO (CAMACHO;
atendidas pelo sistema de controle a ser implemen- BORDONS, 2007). Na função objetivo, a predição é
tado. calculada para um horizonte de predição, sendo que o
controlador projeta incrementos nas variáveis manipu-
Tabela 1. Parâmetros do modelo. ladas em um horizonte de controle, que deve ser me-
Parâmetro Valor Unidade nor que o horizonte de predição. A equação 6 apre-
𝑨𝒂 1,1310 m² senta uma das formas de se definir a função objetivo
𝒌𝟎 2x107 s-1 do controlador para o problema multivariado.
𝑬 7x104 J.mol-1 2
𝑹 8,314462 J.K-1.mol-1 ̂ ‖ + ‖𝚫𝑼‖2𝑅
𝐽 = ‖𝑾 − 𝒀 (6)
𝑄
𝚫𝑯 -7x104 J.mol-1
𝛒 1000 kg.m-3 Sendo:
𝑪𝒑 3000 J.kg-1.K-1 𝑊1 𝑤𝑖 (𝑡 + 𝑁1 )
𝑽𝒋 0,1 m-3 𝑊2 𝑤𝑖 (𝑡 + 𝑁1 + 1)
𝑾 = [ ⋮ ] , 𝑊𝑖 = [ ] , 𝑖 = 1. . 𝑁𝑦 ,
𝑼 900 W.m-2.K-1 ⋮
𝑨 20 m² 𝑊𝑁𝑦 𝑤𝑖 (𝑡 + 𝑁2 )
𝛒𝐣 1000 kg.m-3 𝑌̂1 𝑦̂𝑖 (𝑡 + 𝑁1 |𝑡)
J.kg-1.K-1 ̂
𝑪𝒑,𝒋 4200 ̂ = 𝑌2 , 𝑌̂𝑖 = [𝑦̂𝑖 (𝑡 + 𝑁1 + 1|𝑡)] , 𝑖 = 1. . 𝑁𝑦 ,
𝒀
𝒉 0,5 m ⋮ ⋮
̂
𝑌
[ 𝑁𝑦 ] 𝑦̂𝑖 (𝑡 + 𝑁2 |𝑡)
Δ𝑈1 Δ𝑢𝑖 (𝑡)
Tabela 2. Condições iniciais. Δ𝑈2 Δ𝑢𝑖 (𝑡 + 1)
𝚫𝑼 = [ ⋮ ] , 𝚫𝑼𝒊 = [ ] , 𝑖 = 1. . 𝑁𝑢 ,
Saída Valor Unidade ⋮
𝑪𝒂 (𝟎) 3586 mol.m-3 Δ𝑈𝑁𝑢 Δ𝑢𝑖 (𝑡 + 𝑁3 − 1)
𝑪𝒃 (𝟎) 156,979 mol.m-3 𝛿1 . 𝑰𝑵𝒚 𝟎𝑵𝒚 … 𝟎𝑵𝒚 𝟎𝑵𝒚
𝑻(𝟎) 314,279 K 𝟎𝑵𝒚 𝛿2 . 𝑰𝑵𝒚 … 𝟎𝑵𝒚 𝟎𝑵𝒚
𝑻𝒋 (𝟎) 314,877 K 𝑸= ⋮ ⋮ ⋱ ⋮ ⋮ ,
𝟎𝑵𝒚 𝟎𝑵𝒚 … 𝛿𝑁𝑦 −1 . 𝑰𝑵𝒚 𝟎𝑵𝒚
[ 𝟎𝑵 𝟎𝑵 … 𝟎𝑵 𝛿𝑁 . 𝑰𝑵 ]
𝒚 𝒚 𝒚 𝑦 𝒚
Tabela 3. Ponto de operação inicial.

Entrada Valor inicial Unidade 𝜆1 . 𝑰𝑵𝒖 𝟎𝑵𝒖 … 𝟎𝑵𝒖 𝟎𝑵𝒖

𝑪𝒂,𝒊𝒏 (𝟎) 3900 mol.m-3 𝟎𝑵𝒖 𝜆2 . 𝑰𝑵𝒖 … 𝟎𝑵𝒖 𝟎𝑵𝒖
𝑭𝒊𝒏 (𝟎) 3x10-4 m3.s-1 𝑹= ⋮ ⋮ ⋱ ⋮ ⋮ .
𝑻𝒊𝒏 (𝟎) 295 K 𝟎𝑵𝒖 𝟎𝑵𝒖 … 𝜆𝑁𝑢 −1 . 𝑰𝑵𝒖 𝟎𝑵𝒖
𝑭𝒋 (𝟎) 5x10-4 m3.s-1 [ 𝟎𝑵𝒖 𝟎𝑵𝒖 … 𝟎𝑵𝒖 𝜆𝑁𝑢 . 𝑰𝑵𝒖 ]
𝑻𝒋,𝒊𝒏 (𝟎) 320 K
 Quadro 2. Descrição das variáveis do problema de controle. ⃖
𝑌 Valores de 𝒀 anteriores ao instante
Variável Descrição atual.
Predição para o instante 𝑡 + 𝑘 do va- ⃐
𝑈 Valores de 𝑼 anteriores ao instante
̂𝒊 (𝒕 + 𝑘|𝒕)
𝒚
lor da i-ésima variável controlada, re- atual.
alizada no instante 𝑡.
Existem algumas abordagens diferentes para resolver
𝜹𝒊 Ponderação da i-ésima variável con-
trolada. o problema do MPC não-linear. Uma abordagem in-
tuitiva é a utilização de técnicas de otimização não-
𝚫𝒖𝒊 (𝒕) Valor do incremento da i-ésima vari-
ável manipulada no instante 𝑡. linear para se obter o mínimo da função objetivo, mas
este tipo de solução é complexa. Outra abordagem é a
𝝀𝒊 Ponderação da i-ésima variável ma-
nipulada. utilização de aproximações lineares para a função da
predição, possibilitando utilizar a QP para obter o mí-
𝒘𝒊 (𝒕) Referência a ser seguida pela i-ésima
nimo da função objetivo. Esta é a abordagem utilizada
variável controlada no instante 𝑡.
𝑵𝒚 Número de variáveis controladas no PNMPC, que será a estratégia de controle aplicada
neste trabalho.
𝑵𝒖 Número de variáveis manipuladas
𝑵𝟏 Início do horizonte de predição As variáveis controladas a serem consideradas neste
𝑵𝟐 Fim do horizonte de predição trabalho são a concentração de B (𝐶𝐵 ) e a temperatura
𝑵𝟑 Horizonte de controle (𝑇) na corrente de descarga, as variáveis manipuladas
são as vazões de carga do reator (Fin ) e da camisa
No MPC linear, o modelo utilizado para a predição do (Fj,in ). Também são considerados dois distúrbios me-
comportamento do processo é um modelo linear, por- didos: Temperatura na carga do reator (Tin ) e Tempe-
tanto as predições feitas a partir deste modelo podem ratura na carga da camisa (𝑇j,in ). Também foram in-
ser expressas na forma: cluídas no controlador as restrições operacionais apre-
sentadas na Tabela 4.
̂ 𝑙𝑖𝑛 = 𝑭 + 𝑮. 𝚫𝑼
𝒀 (7)
3.1 O PNMPC
Esta forma de escrever as predições é interessante
pois, após a substituição da equação 7 na equação 6, a O Practical NonlinearModel Predictive Control
função objetivo pode ser escrita na forma padrão da (PNMPC), proposto por (PLUCENIO, 2010), utiliza
programação quadrática (QP): uma aproximação linear do modelo não-linear de pre-
dição obtida a partir da série de Talylor ou Maclaurin.
1 Nesta aproximação, é possível escrever a predição na
𝐽= . 𝚫𝑼𝑇 . 𝑯. 𝚫𝑼 + 𝒃𝑇 . 𝚫𝑼 + 𝑥0 (8) mesma forma da equação 7:
2

Sendo:
̂ = 𝑓𝑛𝑙 (𝒀
𝒀 ⃖ ,𝑼
⃐ , 𝚫𝑼) ≅ 𝑭 + 𝑮𝑷𝑵𝑴𝑷𝑪 . 𝚫𝑼 (10)
𝑯 = 2. (𝑮𝑇 . 𝑸. 𝑮 + 𝑹)
𝒃 = 2. 𝑮𝑇 . 𝑸𝑻 . (𝑭 − 𝑾)
Sendo
É interessante a formulação apresentada na equação 8,
pois a solução de problemas de programação quadrá- ⃖ ,𝑼
𝑭 = 𝑓𝑛𝑙 (𝒀 ⃐ , 𝚫𝑼 = 𝟎𝑁 .𝑁 𝑥1 )
𝑢 3
tica é bem consolidada na literatura. Existem técnicas
eficientes para a resolução de QP’s com restrições, o 𝜕𝑌̂1 𝜕𝑌̂1 𝜕𝑌̂1
que é interessante para um bom desempenho do con- ⋯
𝜕Δ𝑈1 𝜕Δ𝑈2 𝜕Δ𝑈𝑁𝑢
trolador e incentiva o uso deste tipo de formulação do 𝜕𝑌̂2 𝜕𝑌̂2 𝜕𝑌̂2
problema de controle preditivo. ̂
𝜕𝒀 …
𝑮𝑷𝑵𝑴𝑷𝑪 = | = 𝜕Δ𝑈1 𝜕Δ𝑈2 𝜕Δ𝑈𝑁𝑢
𝜕𝚫𝑼 𝑼=𝑼
0 ⋮ ⋮ ⋱ ⋮
A estratégia de controle adotada neste trabalho foi o 𝜕𝑌̂𝑁 𝜕𝑌̂𝑁𝑦 𝜕𝑌̂𝑁
𝑦 𝑦
controle preditivo baseado em modelo não-linear, que ⋯
é caracterizado pela utilização de um modelo não-li- [𝜕Δ𝑈1 𝜕Δ𝑈2 𝜕Δ𝑈𝑁𝑢 ]
near para as predições, que podem ser escritas na 𝜕𝑌̂𝑖
=
forma: 𝜕Δ𝑈𝑗
𝜕𝑦̂𝑖 (𝑡 + 𝑁1 |𝑡)
0 ⋯ 0
̂ = 𝑓𝑛𝑙 (𝒀
⃖ ,𝑼
⃐ , 𝚫𝑼) 𝜕Δ𝑢1 (𝑡)
𝒀 (9)
𝜕𝑦̂𝑖 (𝑡 + 𝑁1 + 1|𝑡) 𝜕𝑦̂𝑖 (𝑡 + 𝑁1 + 1|𝑡)
… 0
𝜕Δ𝑢1 (𝑡) 𝜕Δ𝑢1 (𝑡 + 1) ,
Sendo: ⋮
⋮ ⋮ ⋱
𝒇𝒏𝒍 Função não-linear que calcula as pre- 𝜕𝑦̂𝑖 (𝑡 + 𝑁2 |𝑡) 𝜕𝑦̂𝑖 (𝑡 + 𝑁2 |𝑡) 𝜕𝑦̂𝑖 (𝑡 + 𝑁2|𝑡)
⋯
dições. [ 𝜕Δ𝑢1 (𝑡) 𝜕Δ𝑢1 (𝑡 + 1) 𝜕Δ𝑢1 (𝑡 + 𝑁3 − 1)]
𝑖 = 1. . 𝑁𝑦 , 𝑗 = 1. . 𝑁𝑢
Ou seja, a partir do Jacobiano de 𝒀 ̂ (𝑮𝑷𝑵𝑴𝑷𝑪 ) com re- de uma matriz 𝑮𝑷𝑵𝑴𝑷𝑪 ponderada a partir de duas ou-
lação a 𝚫𝑼, e da avaliação do modelo não-linear para tras matrizes 𝑮𝟏 e 𝑮𝟐 . O Quadro 4 apresenta um algo-
⃖ , ⃐𝑼, 𝚫𝑼 = 𝟎𝑁 .𝑁 𝑥1 )), é
ações de controle nulas (𝑓𝑛𝑙 (𝒀 ritmo para o cálculo da matriz dinâmica em quatro
𝑢 3
obtido um modelo linearizado que permite escrever a passos.
função objetivo do MPC como expresso na equação 8.
Desta forma, a ação de controle pode ser obtida pela Quadro 4. Algoritmo para obter GPNMPC ponderada.
resolução da QP, de forma similar ao MPC linear. No
PNMPC, é obtida uma nova matriz 𝑮𝑷𝑵𝑴𝑷𝑪 a cada pe-
ríodo de amostragem do controlador. A cada período de amostragem do controlador:

3.2 Obtendo 𝑭 e 𝑮𝑷𝑵𝑴𝑷𝑪 - Calcular a matriz 𝐺1 , utilizando o procedimento

mostrado no Quadro 3, utilizando 𝚫𝑼 = 𝟎𝑁𝑢 .𝑁3𝑥1 ;
As derivadas necessárias para o Jacobiano podem ser - Utilizando a QP, obter o 𝚫𝑼∗ que minimiza a fun-
obtidas analiticamente, mas dada a diversidade e com- ção objetivo na forma expressa na equação 8;
plexidade dos modelos não lineares, (PLUCENIO, - Calcular a matriz 𝐺2 , utilizando o procedimento
2010) sugere a utilização de derivadas numéricas e mostrado no Quadro 3, utilizando 𝚫𝑼 = 𝚫𝑼∗ ;
propõe um algoritmo para o cálculo de 𝑭 e 𝑮𝑷𝑵𝑴𝑷𝑪 ,
- 𝐺𝑃𝑁𝑀𝑃𝐶 = 0,5. 𝐺1 + 0,5. 𝐺2 .
que é apresentado a seguir adaptado para o caso mul-
tivariado.

O 𝑮𝑷𝑵𝑴𝑷𝑪 obtido pelo algoritmo do Quadro 4 pode

Quadro 3. Algoritmo para obter F e GPNMPC
ser emtão utilizado para obter as ações de controle a
serem aplicadas no processo. Uma interpretação grá-
Sendo: fica do procedimento e mostrada na Figura 2, para o
caso monovariável (𝑦 = 𝑓(𝑥)).
𝑮𝑖𝑷𝑵𝑴𝑷𝑪 a i-ésima coluna da matriz
𝑮𝑷𝑵𝑴𝑷𝑪;
̂ 𝑖 a i-ésima iteração do vetor 𝑌̂;
𝒀 Figura 2. Interpretação gráfica do gradiente ponderado
𝜀𝑖 um incremento pequeno na i-ésima va- para o caso SISO.
riável manipulada. Uma sugestão para o
valor de 𝜀𝑖 é 𝒖𝒊 (𝒕 − 1)/100;
𝚫𝑼𝑗𝑖 vetor coluna de dimensão 𝑁𝑢 . 𝑁3 x1
com zeros em todas as posições, exceto na
posição i, onde vale 𝜀𝑗 ;

A cada período de amostragem do controlador:

̂ 0 =𝑓𝑛𝑙 (𝒀
- Obter o vetor 𝒀 ⃖ , ⃐𝑼, 𝚫𝑼 = 𝟎𝑁 .𝑁 𝑥1 );
𝑢 3
- Fazer 𝑭 = 𝒀̂ ;
0

Para i de 1 a 𝑁𝑢
Para j de 1 a 𝑁3
(𝑖−1).𝑁 +𝑗
3
- Obter 𝚫𝑼𝑖
-𝒀̂ (𝑖−1).𝑁3 +𝑗 Fonte: (PLUCENIO, 2010)

⃖ , ⃐𝑼, 𝚫𝑼(𝑖−1).𝑁3+𝑗 )
= 𝑓𝑛𝑙 (𝒀
- 𝑮𝑷𝑵𝑴𝑷𝑪3
(𝑖−1).𝑁 +𝑗 3.4 Correção da predição
= (𝒀̂ (𝑖−1).𝑁3+𝑗 − 𝒀
̂ 0 )⁄𝜀𝑖 ;
Os algoritmos de MPC linear, como o DMC e o GPC
Fim Para utilizam mecanismos para fazer correções nas predi-
Fim Para ções dos modelos do processo a cada período de amos-
tragem, com o objetivo de obter erro de predição nulo
no regime permanente. Para o PNMPC, é utilizada
3.3 𝑮𝑷𝑵𝑴𝑷𝑪 ponderada uma versão explícita da correção feita no GPC, so-
mando-se à resposta livre em cada ponto do horizonte
A depender do grau de não-linearidade do processo a de predição a integral do erro de predição filtrado:
ser controlado, a matriz 𝑮𝑷𝑵𝑴𝑷𝑪 obtida pelo procedi-
𝜂(𝑘) = 𝜂(𝑘 − 1). (1 + 𝑓𝑑 )
mento apresentado na seção anterior pode não repre- (11)
sentar adequadamente o comportamento das saídas. −𝑓𝑑 . 𝜂(𝑘 − 2) + 𝑒(𝑘). 𝑘𝑖
Desta forma, (PLUCENIO, 2010) propõe a utilização
Sendo: Tabela 5. Faixas utilizadas na normalização.
𝑒(𝑘) = 𝑦(𝑘) − 𝑦̂(𝑘|𝑘 − 1); Variável Referência Referência Unidade
𝑓𝑑 e 𝑘𝑖 parâmetros do filtro. Inferior Superior
𝒙 𝒙𝒊𝒏𝒇 𝒙𝒔𝒖𝒑
Em (PLUCENIO, 2010), é apresentada uma forma de 𝑪𝒂 3000 4500 mol.m-3
selecionar os coeficientes do filtro pela seleção da di- 𝑪𝒃 0 500 mol.m-3
nâmica desejada para o filtro bem como um exemplo
𝑻 280 320 K
do impacto deste filtro na reposta do sistema de con-
𝑻𝒋 280 390 K
trole na presença de ruídos. Neste trabalho, foram uti-
lizados os coeficientes 𝑓𝑑 = 0 𝑘𝑖 = 1. 𝑭𝒊𝒏 0 6e-4 m3.s-1
𝑻𝒊𝒏 280 310 K
𝑭𝒋 0 10e-4 m3.s-1
3.5 Normalização
𝑻𝒋,𝒊𝒏 300 340 K
Uma operação interessante a ser incluída nos MPC’s
é a normalização das entradas e saídas. Este procedi- 3.6 Sintonia
mento torna as ordens de grandeza das matrizes e ve-
tores utilizados na função objetivo similares, facili- A sintonia dos parâmetros do controlador preditivo foi
tando a definição das ponderações das variáveis mani- realizada empiricamente, após a análise do comporta-
puladas e controladas e possivelmente melhorando a mento do sistema em malha aberta, sendo apresenta-
performance do algoritmo de resolução da QP. A nor- dos na Tabela 6 os parâmetros selecionados.
malização pode ser realizada pela seguinte equação:
Tabela 6. Parâmetros de sintonia do controlador PNMPC
𝑥 − 𝑥𝑖𝑛𝑓
𝑥𝑛 = (12) Parâmetro Valor
𝑥𝑠𝑢𝑝 − 𝑥𝑖𝑛𝑓
𝑵𝟏 1
𝑵𝟐 400
Sendo: 𝑵𝟑 3
𝑥 a variável não normalizada; 𝑻𝒔 35 segundos
𝑥𝑖𝑛𝑓 o valor inferior da variável não norma- 𝝀𝑭𝒊𝒏 1
lizada na faixa considerada; 𝝀𝑭𝒋 1
𝑥𝑠𝑢𝑝 o valor inferior da variável não norma-
𝜹𝑪𝒃 1
lizada na faixa considerada;
𝜹𝑻 1
𝑥𝑛 a variável normalizada.
3.7 Avaliação do modelo não-linear
O processo inverso à normalização é realizado pela
equação:
Para o cálculo das predições do modelo de EDO’s
não-lineares, foi utilizado o método de Euler, descrito
𝑥 = 𝑥𝑛 . (𝑥𝑠𝑢𝑝 − 𝑥𝑖𝑛𝑓 ) + 𝑥𝑖𝑛𝑓 (13) pela equação 16.

No caso especial da variável a ser normalizada ser um 𝑿(𝑘 + 1) = 𝑿(𝑘) + ℎ. 𝑓 (16)

valor incremental, como as ações de controle (Δ𝑈), a
normalização é realizada segundo a equação:
Sendo:
𝑿 os estados do modelo não-linear;
Δ𝑥
Δ𝑥𝑛 = (14) 𝑓 a função não-linear, considerando o mo-
𝑥𝑠𝑢𝑝 − 𝑥𝑖𝑛𝑓 delo escrito na forma 𝑿̇ = 𝑓(𝑡, 𝑿);
ℎ o passo de integração.
Similarmente, o processo inverso é dado pela equa-
ção: É intuitivo selecionar o passo de integração igual ao
tempo de amostragem do MPC, de forma que a cada
Δ𝑥 = Δ𝑥𝑛 . (𝑥𝑠𝑢𝑝 − 𝑥𝑖𝑛𝑓 ) (15) iteração o método de Euler calcula as predições a se-
rem utilizadas na função objetivo. Como a parametri-
zação do modelo do CSTR apresenta ordens de gran-
A Tabela 5 apresenta as faixas utilizadas na normali- deza muito discrepantes, a utilização do método de
zação das variáveis do modelo do CSTR. Euler só é possível quando se considera um passo de
integração pequeno. Ao utilizar passos de integração
igual ao tempo de amostragem do controlador, o mé-
todo não consegue convergir e produz uma solução in-
correta. Nesta situação, o modelo é classificado como
stiff. Após testes variand0-se o passo de integração, foi
selecionado um passo de 1 segundo para a avaliação Quadro 7. Ações tomadas na Situação 2.
do modelo.
Tempo /s Descrição
1000 𝑤𝑇 passa a ser 315,5 K
4. Simulações
10000 ℎ elevado em 0.1 m
Tanto o modelo do processo real quanto o controlador 15000 𝑤𝐶𝐵 passa a ser 220 mol.m-3
PNMPC foram implementados como S-functions no 20000 𝑤𝑇 elevado em 1 K
Matlab® Simulink®, sendo utilizado o solver ode15s.
23000 ℎ reduzido em 0.2 m
Inicialmente foi realizada uma simulação em malha
aberta, para avaliação do comportamento dinâmico do
processo (Situação 0). Para comparar o desempenho  Situação 3
do PNMPC e PNMPC com 𝑮𝑷𝑵𝑴𝑷𝑪 ponderada, fo- Com o sistema em malha fechada, foram aplicadas as
ram realizadas simulações em quatro situações dife- mesmas mudanças de referência da situação 2, mas
rentes (1 a 4) com ambos controladores, sendo estas também foram introduzidas variações na temperatura
descritas a seguir. da carga da camisa e na concentração de A na carga
do reator. Também foi introduzido um ruído branco na
 Situação 0 temperatura na carga do reator. O Quadro 8 apresenta
as alterações realizadas. O tempo de simulação foi de
Com o sistema em malha aberta, foram aplicadas as
3.104 segundos.
ações descritas no Quadro 5. O tempo de simulação
foi de 9.104 segundos.
Quadro 8. Ações tomadas na Situação 3.
Quadro 5. Ações tomadas na Situação 0.

Tempo /s Descrição Tempo /s Descrição

1000 𝐹𝑖𝑛 elevado em 0,5e-4 m3.s-1 1000 𝑤𝑇 passa a ser 315,5 K
9000 𝐹𝑖𝑛 elevado em 0,5e-4 m3.s-1 10000 𝐶𝐴,𝑖𝑛 elevado em 5 mol.m-3
17000 𝐹𝑖𝑛 elevado em 0,5e-4 m3.s-1 15000 𝑤𝐶𝐵 passa a ser 220 mol.m-3
24000 𝐹𝑖𝑛 elevado em 0,5e-4 m3.s-1 20000 𝑤𝑇 elevado em 1 K
32000 𝐹𝑖𝑛 reduzido em 2,0e-4 m3.s-1 25000 𝑇𝑗,𝑖𝑛 elevado em 2 K
48000 𝐹𝑗,𝑖𝑛 elevado em 0,5e-4 m3.s-1
58000 𝐹𝑗,𝑖𝑛 elevado em 0,5e-4 m3.s-1  Situação 4
Na situação final, foram aplicadas as mesmas ações
68000 𝐹𝑗,𝑖𝑛 reduzido em 1,5e-4 m3.s-1 da situação 3 e foram realizadas simulações para ho-
rizonte de controle igual a 2 e 1.
 Situação 1
Com o sistema em malha fechada, foram realizadas Também foi calculada nas situações de 1 a 4 a integral
mudanças na referência para as duas variáveis contro- do erro quadrático (ISE) para auxiliar na comparação
ladas, bem como introduzidos distúrbios na concen- do desempenho entre o PNMPC e o PNMPC com pon-
tração de A na carga, apresentados no Quadro 6. O deração. O ISE é calculado pela seguinte expressão:
tempo de simulação foi de 6.104 segundos.
𝑛𝑖𝑡
Quadro 6. Ações tomadas na Situação 1. 2
ISEi = ∑(wi (t + j) − yi (t + j)) , 𝑖 = 1. . 𝑁𝑦
Tempo /s Descrição 𝑗=1
1000 𝑤𝐶𝐵 passa a ser 200,0 mol.m-3
15000 𝐶𝐴,𝑖𝑛 elevado em 10% de 𝐶𝐴,𝑖𝑛 (0) Sendo 𝑛𝑖𝑡 o número de iterações da simulação.
26000 𝑤𝑇 passa a ser 316 K
30000 𝑤𝐶𝐵 reduzido em 30,0 mol.m-3 5. Resultados
35000 𝐶𝐴,𝑖𝑛 reduzido em 10% de 𝐶𝐴,𝑖𝑛 (0)
A seguir são apresentados, na situação 0, resultados e
40000 𝑤𝑇 reduzido em 1 K
comentários sobre o comportamento dinâmico do pro-
50000 𝐶𝐴,𝑖𝑛 reduzido em 10% de 𝐶𝐴,𝑖𝑛 (0) cesso na simulação, e os resultados e análises do de-
sempenho das duas propostas de controle nas situa-
 Situação 2 ções de 1 a 4.
Com o sistema em malha fechada, foram realizadas
mudanças nas referências para as duas variáveis con-
troladas, além de se variar o nível do reator, que é con-
siderado constante no modelo do controlador. O Qua-
dro 7 apresenta as alterações realizadas. O tempo de
simulação nesta situação foi de 3.104 segundos.
 Figura 4. Gráficos da situação 1.
Figura 3. Gráficos da situação 0

Tabela 7. ISE calculado para a situação 1.

PNMPC com
PNMPC
ponderação
Cb 1.0357e+05 9.6479e+04
T 122.0356 102.4610
 Figura 6. Gráficos da situação 3.
Figura 5. Gráficos da situação 2.

Tabela 8. ISE calculado para a situação 2. Tabela 9. ISE calculado para a situação 3.

PNMPC com PNMPC com

PNMPC PNMPC
ponderação ponderação
Cb 1.6763e+05 1.5262e+05 Cb 1.5240e+05 1.4144e+05
T 96.2622 77.8919 T 104.3836 89.6115
 Figura 7. Gráficos da situação 4.
5.1 Situação 0
Observando-se a variação das variáveis controladas ao
serem aplicados os degraus, é possível observar o ga-
nho variante e os tempos de resposta diferentes entre
as variáveis controladas. Esta simulação é útil para a
escolha do horizonte de controle pois permite observar
o tempo de resposta do sistema.

5.2 Situação 1
Analisando o gráfico da resposta das variáveis contro-
ladas, observa-se que os controladores conseguem se-
guir as referências e rejeitar perturbações, sendo que a
concentração de B parece ser mais sensível aos distúr-
bios aplicados. O PNMPC ponderado apresentou
comportamento muito similar ao PNMPC, mas obser-
vando o ISE, ele é ligeiramente melhor.

5.3 Situação 2
A partir da análise dos gráficos e do ISE, é possível
afirmar que os controladores também apresentam de-
sempenho similar, com o PNMPC ponderado um
pouco melhor. Também é notável o impacto que a va-
riação do nível do processo teve na resposta, sendo
que a vazão da carga da camisa saturou ao compensar
o impacto da variação do nível em certo momento. A
vazão da carga do reator esteve muito próxima de sa-
turar. Isto evidencia o esforço necessário para contro-
lar o sistema quando existem diferenças entre o pro-
cesso e o modelo do controlador. Nesta situação tam-
bém é possível notar que o controlador movimenta a
vazão da carga da camisa com muito mais velocidade
que a vazão de carga do reator.

5.4 Situação 3
Neste exemplo, é notável que o controlador utiliza a
vazão da carga da camisa para compensar o ruído de
alta frequência, já que esta variável manipulada inte-
rage com o sistema com uma dinâmica mais rápida.
Os dois controladores novamente apresentam desem-
penho similar, sendo o PNMPC levemente superior
quando se utiliza o ISE como critério de avaliação.

5.5 Situação 4
A análise tanto do ISE quando do gráfico da resposta
nesta situação permite inferir que a melhor resposta
foi obtida quando o horizonte de controle é unitário, o
que é curioso pois esta solução possui menor custo
Tabela 10. ISE calculado para a situação 4. computacional de cálculo. Apesar disto, quando com-
PNMPC com parado com a resposta da Situação 3, fica claro que o
PNMPC horizonte de controle igual a 3 apresenta uma melhor
ponderação
resposta. Comparando o PNMPC ponderado e o
Cb 1.9611e+05 2.1739e+05 PNMPC, mais uma vez não houveram diferenças sig-
N3 = 2
T 126.3923 160.2671 nificativas, mas desta vez, quando considerado o ISE,
Cb 1.7281e+05 2.0497e+05 o PNMPC apresentou desempenho levemente supe-
N3 = 1 rior. Também é notável que os controladores propaga-
T 71.4218 69.2881
ram o ruído nas variáveis manipulada, o que não é de-
sejado por provocar desgaste de atuadores.
 6. Conclusão

O PNMPC se mostrou capaz de controlar o

CSTR, mas o controle apresentou sensibilidade a
ruídos nos distúrbios não medidos. Os parâmetros
do filtro da correção da predição não foram devi-
damente ajustados e poderiam melhorar a res-
posta do controlador. Não foi notada uma me-
lhora significativa quando se compara o PNMPC
ponderado com o PNMPC. A estratégia de con-
trole estudada neste trabalho se mostrou uma al-
ternativa interessante quando é necessário aplicar
controle preditivo não-linear, por conseguir traba-
lhar em ambiente não-linear e utilizar as ferra-
mentas para resolução de problemas lineares, que
apresentam menor custo computacional.

Referências

PLUCENIO, Agustinho. DESENVOLVIMENTO

DE TÉCNICAS DE CONTROLE NÃO LINEAR
PARA ELEVAÇÃO DE FLUIDOS
MULTIFÁSICOS. 2010. 275 f. Tese (Doutorado) -
Curso de Engenharia de Automação e Sistemas, Uni-
versidade Federal de Santa Catarina, Florianópolis,
2010. Disponível em: <https://reposito-
rio.ufsc.br/xmlui/handle/123456789/93796>. Acesso
em: 11 dez. 2016.

YANG, Bryant Bruce Picon. CONTRIBUIÇÕES

AOS ALGORITMOS PNMPC E DMC COM
PASSO DE ESTABILIDADE. 2016. 132 f. Disser-
tação (Mestrado) - Curso de Engenharia de Automa-
ção e Sistemas, Universidade Federal de Santa Cata-
rina, Florianópolis, 2016.

CAMACHO, E. F.; BORDONS, C.. Model Predictive

control. Advanced Textbooks In Control And Sig-
nal Processing, 2007. Springer London.

CORRIOU, JEAN-PIERRE. Process Control: Theory

and Applications. London: SpringerVerlag,2004. p
684.

MOTA, M. V. de C; VIEIRA, F. T. ANÁLISE DO

REATOR DO TIPO TANQUE AGITADO (CSTR):
COMPORTAMENTO DINÂMICO NA PARTIDA E
ESTABILIDADE DA OPERAÇÃO
ESTACIONÁRIA. COBEQIC. Campinas, jul. 2015.