DAS 5945 - Técnicas de
Controle Avançado aplicadas à
Indústria do Petróleo e Gás.
aula 12
2008.2
1. Estruturas de Controle
- Controle de Razão (Ratio Control)
- Controle de Tanques
- Ganho Escalonado
- Controle adaptativo
- Controle em Cascata
- Preditor de Smith
2. Introdução ao Controle Preditivo Baseado em
Modelo – CPBM ou MBPC(Model Based Predictive
Control)
Controle de Razão (Ratio Control)
Uma necessidade comum em controle de processos é
manipular a razão entre a vazão f de uma substância em
relação a vazão total q.
f
R=
q
q
f
Processo
c=
q
e
Controle
Estrutura A
f
q
Processo
r
+
f
e
Controle
Estrutura B
+
f f = rq
Controle de Tanques
tanque
LT
qsaída
qentrada
O controle do nível é efetuado de forma a absorver as variações da
vazão de entrada mantendo a vazão de saída com variações suaves.
Uma técnica é utilizar controle proporcional com
ec=e*abs(e)
Ganho Escalonado (Gain Schedulling)
Escalonamento
Look-up table
Yr
e
+
-
Controlador
u
Processo
y
Controle adaptativo (Adaptive Control)
Regulador Auto-ajustável
Especificação
Projeto do
Controlador
Yr
+
e
-
Controlador
Estimador
u
Processo
Controle em Cascata (Cascade Control)
Exemplo:
Temperatura
desejada
Yr
+
e
Vazão de vapor
desejada
Cp
-
u1
+
-
Vazão
Temperatura
Abertura de
Válvula
Cs
u2
Temperatura
Processo
Y
Preditor de Smith
Controlador
yr
e
+
-
C
Processo
y
p
u
~
p
e
− Ls
-+
+
+
O Controlador C é projetado como se não houvesse
atraso de transporte.
CONTROLE PREDITIVO BASEADO EM MODELO
• Uso explícito de um modelo para predizer
as saídas do processo em instantes de
tempo futuro. (horizonte de predição)
CPBM
• Cálculo de uma sequência de controle
(horizonte de controle) que minimiza uma
função objetivo.
• A cada instante o horizonte é deslocado
em direção ao futuro. Aplica-se somente o
primeiro sinal da sequência calculada em
cada iteração.
Estrutura básica dos controladores Preditivos
Trajetória da
referência
Entradas e
saídas passadas
Modelo
para
predição
Saídas
previstas
Erros
futuros
Entradas
Futuras
Otimizador
( Entradas Futuras )
Função Custo
Restrições
Modelos para predição:
Modelos não paramétricos: ( Plantas estáveis )
1-Resposta ao Impulso Unitário (MAC, EPSAC)
2-Resposta ao Degrau Unitário (DMC)
Modelo paramétricos:
1- Função de Transferência
(GPC,UPC,EPSAC,EHAC,MUSMAR ou MURHAC)
2- Equações de Estado (PFC)
3- Modelo não linear (NN, Fuzzy, PNMPC)
Modelo baseado na resposta ao impulso unitário
y (k ) = ∑i =1 hi u (k − i )
∞
Os termos hi são os valores nos instantes k da resposta
ao impulso .
Coeficientes hi
Impulse Response
0.35
0.3
Amplitude
0.25
0.2
0.15
0.1
0.05
0
0
5
10
15
Time (sec)
20
25
Resp. ao degrau
•
0
0
•
0
0
•
0.3000
1.0000
•
0.3000
1.0000
•
0.2100
2.0000
•
0.5100
2.0000
•
0.1470
3.0000
•
0.6570
3.0000
•
0.1029
4.0000
•
0.7599
4.0000
•
0.0720
5.0000
•
0.8319
5.0000
•
0.0504
6.0000
•
0.8824
6.0000
•
0.0353
7.0000
•
0.9176
7.0000
•
0.0247
8.0000
•
0.9424
8.0000
•
0.0173
9.0000
•
0.9596
9.0000
•
0.0121 10.0000
•
0.9718 10.0000
•
0.0085 11.0000
•
0.9802 11.0000
u(k)
1 2 3 ...
y (k ) =
∑
∞
h u (k − i)
i =1 i
y (4) = h1u (4 − 1) + h2u (4 − 2) + h3u (4 − 3) + h4u (4 − 4)
y (4) = h1 *1 + h2 *1 + h3 *1 + h4 * 0
y (4) = 0 *1 + 0.3 *1 + .21 *1 + .147 * 0
Y (4) = 0.51
K-N
Para sistemas estáveis existe um N tal que hN ≅ 0
y (k ) = ∑i =1 hi u (k − i ) → y (k ) = ∑i =1 hi u (k − i )
∞
N
Modelo FIR (Finite Impulse Response)
Modelo baseado na resposta ao degrau
y (k ) = ∑i =1 g i ∆u (k − i )
∞
y (0) = 0
y (4) = g1∆u (4 − 1) + g 2 ∆u (4 − 2) + g 3 ∆u (4 − 3)
y (4) = g1∆u (3) + g 2 ∆u (2) + g 3 ∆u (1)
y (4) = g1 (u (3) − u (2)) + g 2 (u (2) − u (1)) + g 3 (u (1) − u (0))
y (4) = g1 0 + g 2 0 + g 31
y (4) = 0.51
TAREFA
Obter o gráfico de y(k) utilizando o modelo baseado em
reposta ao degrau unitário.
∆u (k )
1.0
u (k )
1
1 − z −1
.2 z −1
1 − 0.8 z −1
u (k )
0.5
1.0
∆u (k )
9
17
k
0.5
17
9
k
TAREFA
solução
clear all
close all
z=tf('z',1);
G=.2*z^-1/(1-.8*z^-1);
Ng=70;
[rd]=step(G,Ng+1);
g=rd(2:Ng+1);
du=[0.5;zeros(7,1);0.5;zeros(7,1);-.5;zeros(Ng-17,1)];
N=length(du)
for k=1:N+1
som=0;
for i=1:Ng
j=k-i;
if j>0
som=som+g(i)*du(j)
else
break
end
end
y(k)=som;
End
plot(y)
u=[0.5*ones(8,1);ones(8,1);.5*ones(8,1)];
yft(1)=0;
for s=2:25
yft(s)=.8*yft(s-1)+.2*u(s-1);
End
Hold
plot(yft,'kd')
Em geral os algoritmos MBPC utilizam a seguinte
estratégia:
1- Faz-se uma predição da saída para os próximos p
instantes de amostragem considerando que não
haverá alteração na ação de controle (∆u(k+j)=0) para
j=0 até m-1.
2- Calcula-se os próximos ∆u(k+j), j=0 até m-1, de
forma a minimizar uma função objetivo.
A função objetivo típica é uma medida da diferença
entre as referências futuras e as predições futuras em
um horizonte (horizonte de predição) mais o esforço
de controle em outro horizonte (horizonte de controle).
Aplica somente
∆u (k )
1
1 − z −1
∆u
Calcula
u (k )
Processo
é o vetor com incrementos de ação
futura (m elementos) que minimiza uma
função custo.
⎡ ∆u (k ) ⎤
⎢ ∆u (k + 1) ⎥
⎥
∆u = ⎢
M
⎥
⎢
⎥
⎢
⎣∆u (k + m − 1)⎦
MPC calcula o vetor ∆u mas aplica somente
o primeiro elemento do vetor ∆u, i.e., ∆u(k).
Dynamic Matrix Control - DMC
Cutler e Ramaker, Shell Oil Co. (1980)
Tipo de Controle Preditivo bem aceito na Indústria
Petroquímica.
Principais vantagens:
• O modelo é a resposta ao degrau (de fácil obtenção)
• Trata processos multivariáveis com facilidade.
• Fácil compreensão e ajuste.
Principais desvantagens:
• Trabalha somente com plantas estáveis (Existe
versão para operar com processos integrativos)
DMC – Montagem da Função Objetivo
1º. Passo = Escrever as predições
y (k ) = ∑i =1 g i ∆u (k − i ) + n(k )
∞
∞
~
y (k + j / k ) = ∑i =1 g i ∆u (k + j − i ) + n~ (k + j / k )
Significa o valor da saída y j tempos de
amostragem na frente do instante atual k.
Escrevendo as predições...
∞
~
y (k + j / k ) = ∑i =1 g i ∆u (k + j − i ) + n~ (k + j / k )
∞
j
~
y (k + j / k ) = ∑i =1 g i ∆u (k + j − i ) + ∑i = j +1 g i ∆u (k + j − i ) + n~ (k + j / k )
Contribuição de
Contribuição de
∆u futuro k, k+1, k+2...
∆u passado k-1, k-2, k-3...
Perturbação
~
y (k + j / k ) = ∆y (k + j / k )(∆u f ) + y∆pu f =0 (k + j / k ) + n~ (k + j / k )
A perturbação para todos os instantes futuros é considerada igual a
perturbação no instante k.
n~ (k + j / k ) = n~ (k / k ) = y m (k ) − ~
y (k / k ) para todo j ∈ [1, p ]
Escrevendo as predições...
~
Futuro ∆y (k + j / k )
~
passado y (k + j / k ) ∆u =0
j
∞
~
y (k + j / k ) = ∑i =1 g i ∆u (k + j − i ) + ∑i = j +1 g i ∆u (k + j − i )
+ y (k ) − ~
y (k / k )
m
perturbação
futuro
passado
j
∞
~
y (k + j / k ) = ∑i =1 g i ∆u (k + j − i ) + ∑i = j +1 g i ∆u (k + j − i )
+ y m (k ) − ∑i =1 g i ∆u (k − i )
∞
perturbação
Escrevendo as predições ...
∑ g ∆u(k + j − i) + f (k + j)
f (k + j ) = ∑
g ∆u (k + j − i ) − ∑ g ∆u (k − i ) + y (k )
f (k + j ) = ∑ g ∆u (k − i ) − ∑ g ∆u (k − i ) + y (k )
f (k + j ) = ∑ ( g − g )∆u (k − i ) + y (k )
~
y (k + j / k ) =
j
i =1 i
∞
∞
i = j +1 i
i =1 i
∞
i =1
∞
i =1
m
∞
j +i
i =1 i
j +i
i
m
m
g j +i − g i = 0 para i = N
f (k + j ) = ∑i =1 ( g j +i − g i )∆u (k − i ) + y m (k )
N
DMC – Escrevendo as predições ...
⎡ g1
⎢g
⎢ 2
⎢ g3
⎢
j
∑i =1 g i ∆u (k + j − i) = ⎢ g 4
⎢ g5
j = 1: p
⎢
⎢ M
⎢g
⎣ p
0
g1
0
0
0
0
g2
g3
g1
g2
0
g1
g4
M
g p −1
g3
M
g p−2
g2
M
g p −3
j
~
y (k + j / k ) = ∑i =1 g i ∆u (k + j − i ) + f (k + j )
~
Y px1 = G pxm ∆u mx1 + f px1
⎤ ⎡ ∆u (k ) ⎤
⎥⎢
⎥
u
(
k
1
)
∆
+
⎥⎢
⎥
0 ⎥ ⎢ ∆u (k + 2) ⎥
⎥⎢
⎥
0 ⎥ ⎢ ∆u (k + 3) ⎥
g1 ⎥ ⎢ ∆u (k + 4) ⎥
⎥⎢
⎥
M ⎥⎢
M
⎥
g p − m +1 ⎥⎦ ⎢⎣∆u (k + m − 1)⎥⎦
0
0
Representação AP para as predições (válida para
sistemas lineares e não lineares)
~
Y = F + G∆u
onde
~
F = Y∆u =0
~
∂Y
G=
∂∆u
⎡ ∂y (k + 1)
⎢ ∂∆u (k )
⎢
⎢ ∂y (k + 2)
G = ⎢ ∂∆u (k )
⎢
M
⎢ ∂y (k + p )
⎢
⎣ ∂∆u (k )
⎤
⎥
⎥
∂y (k + 2)
⎥
0
L
⎥
∂∆u (k + 1)
⎥
0
M
L
∂y (k + 2)
∂y (k + p ) ⎥
L
⎥
∂∆u (k + 1)
∂∆u (k + m − 1) ⎦
0
L
0
Representação AP para as predições (válida
para sistemas lineares e não lineares)
⎡ ∂y (k + 1)
⎢ ∂∆u (k )
⎡ y (k + 1) ⎤ ⎡ y (k + 1) ⎤
⎢
⎢ y (k + 2) ⎥ ⎢ y (k + 2) ⎥
⎢ ∂y (k + 2)
⎢
⎥=⎢
⎥
+ ⎢ ∂∆u (k )
⎢
⎥
⎢
⎥
M
M
⎢
M
⎢
⎥ ⎢
⎥
⎢
∂y (k + p )
⎣ y (k + p)⎦ ⎣ y (k + p )⎦ ∆u =0
⎢
⎣ ∂∆u (k )
Jacobiano de
⎤
⎥
∆u (k ) ⎤
⎥⎡
⎢ ∆u (k + 1) ⎥
∂y (k + 2)
⎥
0
L
⎥
⎥⎢
∂∆u (k + 1)
⎢
⎥
M
⎥
0
M
L
⎢
⎥
∂y (k + 2)
∂y (k + p ) ⎥ ⎣∆u (k + m − 1)⎦
L
⎥
∂∆u (k + 1)
∂∆u (k + m − 1) ⎦
0
L
~
Y(∆u)
0
Exercício:
Para os sistema representado pela FT
y( z)
bz −1
=
u ( z ) 1 − az −1
1- Encontre a representação das predições para p=3 e
m=3 a partir da definição AP.
Lembre:
u (k ) = u (k − 1) + ∆u (k )
u (k + 1) = u (k ) + ∆u (k + 1) = u (k − 1) + ∆u (k ) + ∆u (k + 1)
u (k + 2) = u (k + 1) + ∆u (k + 2) = u (k − 1) + ∆u (k ) + ∆u (k + 1) + ∆u (k + 2)
2- Mostre que a matriz G obtida é igual aquela obtida
com os coeficientes da resposta ao degrau (sistema
linear).
Solução:
y (k + 1 / k ) ?
y( z)
bz −1
=
u ( z ) 1 − az −1
y ( k + 1 / k ) = ay (k ) + bu (k )
y (k ) = ay (k − 1) + bu (k − 1)
y ( k + 1 / k ) = a[ay (k − 1) + bu (k − 1)] + b[u (k − 1) + ∆u (k )]
y ( k + 1 / k ) = f1 (u (k − 1), y (k − 1)) + b∆u ( k )
y (k + 2 / k ) ?
y ( k + 2 / k ) = ay (k + 1) + bu ( k + 1)
y ( k + 2 / k ) = a[ f1 (u (k − 1), y (k − 1)) + b∆u ( k )] + b[u (k − 1) + ∆u ( k ) + ∆u (k + 1)]
y (k + 2 / k ) = f 2 ( y (k − 1), u (k − 1)) + ab∆u (k ) + b∆u (k + 1)
De forma análoga ...
y (k + 3 / k ) ?
y ( k + 3 / k ) = ay ( k + 2) + bu (k + 2)....
y (k + 3 / k ) = f 3 ( y (k − 1), u (k − 1)) + a 2b∆u (k ) + ab∆u (k + 1) + b∆u (k + 2)
Solução:
⎡ ∂y ( k + 1)
⎢ ∂∆u (k )
⎡ y ( k + 1) ⎤ ⎡ y ( k + 1) ⎤
⎢
⎢ y ( k + 2) ⎥ ⎢ y ( k + 2) ⎥
⎢ ∂y (k + 2)
⎢
⎥=⎢
⎥
+ ⎢ ∂∆u (k )
⎢
⎥ ⎢
⎥
M
M
⎢
M
⎢
⎥ ⎢
⎥
⎣ y (k + p )⎦ ⎣ y (k + p )⎦ ∆u =0 ⎢ ∂y (k + p )
⎢
⎣ ∂∆u (k )
⎤
⎥
∆u (k ) ⎤
⎥⎡
∂y (k + 2)
⎥ ⎢ ∆u (k + 1) ⎥
L
0
⎥
⎥⎢
∂∆u (k + 1)
⎢
⎥
M
⎥
M
L
0
⎢
⎥
∂y (k + 2)
∂y (k + p ) ⎥ ⎣∆u (k + m − 1)⎦
L
⎥
∂∆u (k + 1)
∂∆u (k + m − 1) ⎦
0
L
0
∂y (k + 2)
∂y (k + 1)
∂y (k + 3)
= b,
= ab + b,
= a 2b + ab + b
∂∆u (k )
∂∆u ( k )
∂∆u (k )
∂y (k + 2)
∂y (k + 2)
=b
= ab + b,
∂∆u (k + 1)
∂∆u (k + 1)
∂y ( k + 3)
=b
∂∆u (k + 2)
1)
b
0
0⎤ ⎡ ∆u (k ) ⎤
⎡ y (k + 1) ⎤ ⎡ f1 ( y (k − 1), u (k − 1)) ⎤
⎡
⎢ y (k + 2)⎥ = ⎢ f ( y (k − 1), u (k − 1))⎥
⎢ (ab + b)
⎥ ⎢ ∆u (k + 1) ⎥
b
0
+
⎥
⎢
⎥ ⎢ 2
⎥
⎢
⎥⎢
2
⎢⎣ y (k + 3) ⎥⎦ ⎢⎣ f 3 ( y (k − 1), u (k − 1)) ⎥⎦ ∆u =0 ⎢⎣(a b + ab + b) (ab + b) b ⎥⎦ ⎢⎣∆u (k + 2) ⎥⎦
Exercício:
2- Mostre que a matriz G obtida é igual aquela obtida com os coeficientes
da resposta ao degrau (sistema linear).
y (k ) = ay (k − 1) + bu ( k − 1), y (1) = 0, u (k ) = 1 for k ≥ 1
Obtenção dos coeficientes da resposta ao degrau unitário...
g1 = y ( 2) = ay (1) + bu (1) = b
g 2 = y (3) = ay (2) + bu (2) = ab + b
g 3 = y (4) = ay (3) + bu (3) = a(ab + b) + b = a 2b + ab + b
2)
⎡ g1
G = ⎢⎢ g 2
⎢⎣ g 3
0
g1
g2
0⎤
b
⎡
0 ⎥⎥, G = ⎢⎢ (ab + b )
⎢⎣ a 2b + ab + b
g1 ⎥⎦
(
)
0
b
(ab + b )
0⎤
0⎥⎥
b ⎥⎦
DMC (continuação)
1º. Passo– Forma compacta de escrever as predições
~
Y = G∆u + f
2º. Passo = Escrever a função objetivo
Para o caso SISO ...
Função objetivo:
p
min J
∆u
m
J = ∑ (~
y (k + j / k ) − w(k + j / k ) ) + λ ∑ ∆u (k + j − 1) 2
j =1
2
j =1
Função objetivo:
min J
∆u
p
m
J = ∑ (~
y (k + j / k ) − w(k + j / k ) ) + λ ∑ ∆u (k + j − 1) 2
2
j =1
j =1
Escrevendo na forma matricial ...
w(k + j / k ) for j = 1 : p → W px1
~
~
y (k + j / k ) → Y = G∆u + f
J = (G ∆ u + f − W ) T (G ∆ u + f − W ) + λ∆ u T ∆ u
J = (G ∆ u + f − W ) T (G ∆ u + f − W ) + λ ∆ u T ∆ u
J = (∆ u T G T + f T − W T )(G ∆ u + f − W ) + λ ∆ u T ∆ u
∂J
= 2G T G ∆ u + 2G T f − 2G T W + 2 λ ∆ u = 0
∂∆u
∂J
=0
∂∆u
∂2 J
T
=
2
G
G + 2λI é positiva definida
2
∂∆u
(G T G + λ I ) ∆ u = G T ( W − f )
−1
∆ u = (G G + λ I ) G (W − f)
T
T
Procedimento do algoritmo:
1- Inicialização:
Definir g (Nx1)= vetor com coeficientes da resposta ao degrau do
modelo da planta.
p = horizonte de predição
m = horizonte de controle
λ = ponderação da ação de controle
α = filtro de referência
G(pxm) = Matriz com os coeficientes da resposta ao degrau g.
Calcular
K 1xp = (G T G + λ I ) −1 G T (1 , :)
Montar os vetores
Ua = Vetor com as entradas passadas (iniciar com Ua(:)=urp
Ya = Vetor com saídas passadas (iniciar com Ya(:)=
Procedimento:
2- No instante k:
A-Calcular: y(1) com os valores anteriores:u(k-1), u(k-2)... y(k-1), y(k2)...etc.
B-Atualizar o vetor referências W(p,1)
C-Atualizar Vetor com saídas passadas Y(na,1) colocando y(k) no
topo e deslocando os demais valores para baixo.
D-Calcular: o vetor f
∆u = K(W − f)
u (k ) = u (k − 1) + du
du(k)
E-Atualizar
o vetor Ua deslocando os elementos para
baixo e colocando u(k) no topo.
F-Voltar a 2-A com k=k+1
Parâmetros importantes do DMC
Métodos para escolha de p.
O valor de p (número de iterações do horizonte de
predição) pode ser selecionado como segue:
Como leva cerca de 4τ para um processo de 1a.
ordem atingir regime permanente,
4τ
pTs = 4τ → p =
Ts
Onde τ é a constante de tempo dominante do processo.
t acom.
t acom.
< Ts <
1)
15
6
Se
t acom. = pTs
Tacom.=tempo de acomodação ( 95%)
Isermann
pTs
pTs
< Ts <
→
15
6
−1− d
6 < p < 15
bz
2) G ( z ) =
, a=e
−1
1 − az
Ts
Ts
1
ln(a) = − , ln( ) =
a
τ
τ
−
Ts
τ
4τ 4
=
ln(1 / a)
Ts
p= 4
ln(1 / a)
3) 20 < p < 30
Coleman Brosilow and Babu Joseph
O valor de p (horizonte de predição) determina o
tamanho da matriz G.
pTs = 4τ
Aumentar Ts para diminuir p pode
significar um intervalo de tempo muito
grande entre aplicações de controle, o
que não é recomendável.
p muito grande pode significar ênfase no regime
permanente e menos ênfase no transitório !
p muito pequeno pode significar ênfase no transitório
sem olhar para onde vai o regime permanente.!
Existe um compromisso ÓTIMO !
O horizonte de controle m define o grau de liberdade
sobre u(t) futuro de forma a obter o ótimo da função
objetivo.
Um valor grande para m deixa o problema para ser
resolvido no futuro.
Um valor pequeno força a solução imediata.
Processos complexos necessitam de m maior pois a
solução ótima pode necessitar de mais manipulações
de controle.
m pequeno reduz o tamanho da matriz (GtG +λI) a ser
invertida.
1) 1 ≤ m ≤ 5
Coleman Brosilow and Babu Joseph
TAREFA
Água Fria
Gás
Água Quente
y( z)
0.2713 z −3
H ( z) =
=
u ( z ) 1 − 0.8351z −1
Fazer o Controle da Temperatura (y(k))
manipulando a abertura (u(k) da válvula de controle
utilizando o algoritmo DMC.
Exemplo do Livro MPC de E. F. Camacho and C. Bordons
DMC
y
2
Saída
referência
1.5
1
0.5
0
0
20
40
60
80
100
120
140
160
180
tempo
du
1.5
u(k)
du(k)
1
0.5
0
-0.5
-1
0
20
40
60
80
100
tempo
120
140
160
180