Curso de Especialização em Automação Industrial
Grupo de Controle Automação e Robótica
GCAR/UFRGS
Controle Preditivo Baseado em Modelo
- Controle Preditivo Generalizado (GPC )
Prof. Dr. João Manoel Gomes da Silva Jr.
Vantagens
• Pode tratar sistemas instáveis
• Pode tratar sistemas de fase não mínima
• Pode utilizar pesos diferentes em cada instante na
função custo
maior flexibilidade de sintonia
Característica: utilização de um modelo baseado na
respresentação por função de transferência
Modelo
• Modelo CARMA (controller auto-regressive
moving average)
função de transferência
A( z 1 ) y (t )  z d B( z 1 )u (t  1)  C ( z 1 )e(t )
A( z 1 )  1  a1 z 1  a2 z 2    ana z na
B( z 1 )  b0  b1 z 1  b2 z 2    bnb z nb
C ( z 1 )  1  c1 z 1  a2 z 2    cnc z nc
e(t): ruído branco de média zero
u(t): sinal de controle
y(t): saída do processo
Modelo
• Modelo CARIMA
A( z 1 ) y (t )  B ( z 1 ) z d u (t  1)  C ( z 1 )
  1  z 1
Em geral C(z-1) é escolhido como 1
e(t )

Critério
J ( N1 , N 2 , N u ) 
N2
Nu
J  N1
j 1
2
2

(
j
)[
y
(
t

j
|
t
)

w
(
t

j
)]


(
j
)[

u
(
t

j

1
)]
ˆ


Predição
• Para definir a equação de predição, consideremos
primeiramente a seguinte equação diofantina
~
~
1  E j ( z 1 ) A( z 1 )  z  j Fj ( z 1 ) com A( z 1 )  A( z 1 )
Os polinômios Ej e Fj são determinados de forma
única e possuem grau (j-1) e na, respectivamente. Eles
~ 1
A
podem ser obtidos dividindo-se 1 por ( z ) até que o
resto possa ser fatorado como z  j Fj ( z 1 ) . O quociente
desta divisão será E j ( z 1 )
Predição
• Multiplicando-se a equação diofantina por z j E j ( z 1 ).
obtém-se o seguinte:
~
A( z 1 ) E j ( z 1 ) y (t  j )  E j ( z 1 ) B( z 1 )u (t  j  d  1) 
 E j ( z 1 )e(t  j )
(1  z  j Fj ( z 1 )) y(t  j )  E j ( z 1 ) B( z 1 )u(t  j  d  1)  E j ( z 1 )e(t  j )
~ 1
j
1
• Como 1  E j ( z ) A( z )  z Fj ( z )
1
y(t  j )  Fj ( z 1 ) y(t )  E j ( z 1 ) B( z 1 )u(t  j  d  1)  E j ( z 1 )e(t  j )
Predição
• Como o grau do polinômio E j ( z 1 ) é igual a j-1 os
termos de ruído são todos no futuro, logo a melhor
predição é dada por :
yˆ (t  j | t )  G j ( z 1 )u(t  j  d  1)  Fj ( z 1 ) y(t )
Onde:
G j ( z 1 )  E j ( z 1 ) B( z 1 )
Predição
• Considerando-se
F j ( z 1 )  f j , 0  f j ,1 z 1    f j ,na z na
E j ( z 1 )  e j , 0  e j ,1 z 1    e j , j 1 z ( j 1)
• Pode-se obter Gj iterativamente da seguinte maneira:
G j 1  E j 1 B  ( E j  f j , 0 z  j ) B
G j 1  G j  f j , 0 z  j B
N-Predições à frente
yˆ (t  d  1 | t )
yˆ (t  d  2 | t )

Gd 1u (t )  Fd 1 y (t )

Gd 2 u (t  1)  Fd 2 y (t )

yˆ (t  d  N | t )  Gd n u (t  N  1)  Fd  Ny (t )
y  Gu  F( z 1 ) y (t )  G ' ( z 1 )u (t  1)
Matrizes
 yˆ (t  d  1 | t ) 
 u (t ) 
 yˆ (t  d  2 | t ) 
 u (t  1) 
 u

y






 yˆ (t  d  N | t )
u (t  N  1)




 g0
 g
G 1
 
g
 N 1
Resposta ao
salto
0
 0
g0  0 




g N 2  g 0 
 Fd 1( z 1 ) 


(Gd 1 ( z 1 )  g 0 ) z
F



1
1
1
2
d 2( z ) 
1

(
G
(
z
)

g

g
z
)
z
1
d

2
0
1
F
(
z
)


G' ( z )  
  



F


1
1
1
( N 1)
N
 d N ( z ) 
)z 
Gd  N ( z )  g 0  g1 z    g N 1 z
Resposta Livre
y  Gu  F( z 1 ) y (t )  G ' ( z 1 )u (t  1)
• Como os dois últimos termos desta equação dependem
apenas do passado podemos agrupá-los em f:
y  Gu  f
f =resposta livre
~ 1
f j1  z (1  A( z ))f j  B( z 1 )u(t  d  j )
Com: f 0  y (t )
e u (t  j )  0
para
j0
Critério
J  (Gu  f  w )T (Gu  f  w )  u T u
w  [ w(t  d  1) w(t  d  2)  w(t  d  N )]T
1 T
J  u Hu  b T u  f 0
2
Com:
H  2(G T G  λI)
b T  2(f  w) T G
f 0  (f  w) T (f  w)
Solução Ótima
u  H 1b  (G T G  λI) 1 G T (w  f)
• Estratégia de horizonte deslizante:
u(t )  K(w  f)
K=primeira linha da matrix (G T G  λI) 1 G T
Esquema Geral
Observações
• Note que a ação de controle é calculada
com relação a erros (preditos) futuros e não
passados
• Matriz a ser invertida de ordem N
(=horizonte de predição)
• EPSAC: Nu=1  limitado com relação a
otimalidade.
Exemplo
Sistema de 1a ordem:
e(t )
(1  az ) y (t )  (b0  b1 z )u (t  1) 

com a=0.8, bo=0.4, b1=0.6, N1=1, N=Nu=3
1
1
~
A( z 1 )  A( z 1 )(1  z 1 )  1  1.8 z 1  0.8 z 2
E1 ( z 1 )  1 F1 ( z 1 )  1.8  0.8 z 1
E2  1  1.8 z 1
F2  2.44  1.44 z 1
E 3 1  1.8 z 1  2.44 z 2
F3  2.952  1.952z 1
Exemplo
G1  0.4  0.6 z 1
G2  0.4  1.32 z 1  1.08z 2
G3  0.4  1.32 z 1  2.056z 2  1.464z 3
0
0   u (t )  
0 .6  u ( t
 yˆ (t  1 | t )   0.4
 yˆ (t 1. 20 8| t)u 
( t  11 .)32
 2 . 4 04 .4y ( t ) 0 1 . 4 
4 uy ((tt  11))  1 . 0 8  u ( t 

 

 
 yˆ (t  3 | t )  2.056 1.32 0.4 u (t  2)  1 . 4 6 4  u ( t 
    
 0.6u (t  1)  1.8 y (t )  0.8 y (t  1)

  1.08u (t  1)  2.44 y (t )  1.44 y (t  1) 


1.464u (t  1)  2.952 y (t )  1.952 y (t  1)

f
Exemplo
• Considerando:   0.8
0.286 0.147 
 0.133
(G T G  λI) 1 G T   0.154  0.165 0.286 


 0.029  0.154 0.1334
Exemplo
u (t )   0.604u (t  1)  1.371y (t )  0.805 y (t  1) 
 0.133w(t  1)  0.286w(t  2)  0.147w(t  3)
Considerando uma referência constante:
u (t )  0.396u (t  1)  0.604u (t  2)  1.371y (t )  0.805 y (t  1) 
0.133w(t  1)  0.286w(t  2)  0.147w(t  3)
ou seja, o sinal de controle é função de saídas e controles
passados e referências futuras
Exemplo
Exemplo
Cálculo sem passar pela solução da eq. Diofantina  usar
resposta ao salto
j
j 1
g j   ai g j 1   bi
i 1
com
i 0
g k  0, k  0
g 0  b0  0.4
g1  a1 g 0  b0  b1  1.32
g 2  a1 g1  a2 g 2  a3 g 0  b0  b1  2.056
0
0
 0.4
G   1.32 0.4 0 


2.056 1.32 0.4
Exemplo
A resposta livre pode ser calculada através do modelo
considerando que todos os controles futuros são iguais a u(t-1):
y (t )
 0.8 y (t  1)  0.4u (t  1)  0.6u (t  2)
y (t  1) 
0.8 y (t )  0.4u (t )  0.6u (t  1)
y(t  1)  1.8 y(t )  0.8 y(t 1)  0.4u(t )  0.6u(t 1)
f (t  1)

1.8 y (t )  0.8 y (t  1)  0.6u (t  1)
f (t  2)  1.8 f (t  1)  0.8 y (t )  2.44 y (t )  1.44 y (t  1)  1.08u (t  1)
f (t  3) 
1.8 f (t  2)  0.8 f (t  1) 

2.952 y (t )  1.952 y (t  1)  1.464u (t  1)