Curso de Controle Discreto 2
Pe. Pedro M. Guimarães Ferreira S.J.
http://www.fplf.org.br/pedro_varios/
(Texto básico deste curso: Katsuiko Ogata, Discrete-time Control Systems. PrenticeHall, Second Edition, 1995)
CAPÍTULO 4: SÍNTESE DE SISTEMAS DE CONTROLE DE
TEMPO DISCRETO ATRAVÉS DE MÉTODOS CONVENCIONAIS
4-1 Introdução
Os métodos apresentados neste capítulo têm sido bem dominados desde os anos 1950’s.
O primeiro método é baseado no “root locus” já estudado nos sistemas de tempo
contínuo, uma ferramenta ainda muito útil, pois nos dá uma visualização do que
acontece com os pólos da função de transferência quando o ganho do sistema em malha
fechada varia.
O segundo método é baseado na resposta no domínio da frequência.
O terceiro método, ao contrário dos dois primeiros, que são gráficos, é um método
analítico, no qual tentamos obter um certo comportamento do sistema em malha
fechada, manipulando a função de transferência de pulso do controlador digital.
4-2 Mapeamento entre o plano s e o plano z
Como se sabe, tanto a estabilidade absoluta como a relativa de um sistema de tempo
contínuo depende da localização de seus pólos no plano s.
Tendo em vista que as variáveis s e z são relacionadas por z = eTs , é claro que os pólos
e zeros no plano z estão relacionados aos do plano s. Consequentemente, a estabilidade
dos sistemas lineares de tempo discreto também fica determinada pela localização dos
pólos no plano z.
Ora, definindo as partes real e imaginária:
s = σ + jω ,
(*)
T (σ + jω )
T σ jT ω
T σ j ( T ω + 2π k )
temos z = e
=e e =e e
.
(**)
Consequentemente, para cada s, há uma infinidade de valores de z, pois
k = 0; ±1; ±2; ±3;...
Se o sistema de tempo contínuo for estável, temos σ < 0 e se o sistema de tempo
discreto for estável, temos z = eTσ < 1 .
Observe-se que o eixo dos imaginários no plano s corresponde a z = 1 . Ou seja, o eixo
imaginário ( σ = 0 ) corresponde ao círculo unitário no plano z, e o semi-plano aberto da
esquerda para os sistemas de tempo contínuo corresponde ao interior do círculo unitário
para os sistemas de tempo discreto.
Faixas primárias e complementares
Considere os pontos no eixo imaginário do plano s. Seja ωs a frequência de
amostragem.
Consideremos (*) e (**) acima; quando um ponto se move no eixo imaginário do plano
1
1
s de − j ωs a j ωs , no plano z nós temos z = 1 e arg( z ) varia de −ωsT = −π a
2
2
ωsT = π no sentido trigonométrico, isto é, no sentido anti-horário, percorrendo todo o
1
1
1
círculo unitário, pois quando ω varia de − j ωs a j ωs , passa por ω = 0, o que
2
2
corresponde a ω T = 0.
Consequentemente, quando s percorre o eixo imaginário de −∞ a +∞ , o círculo
unitário no plano z é percorrido um número infinito de vezes. Isto implica que o semiplano da esquerda no plano s pode ser dividido em um número infinito de faixas de
largura ωs , conforme a figura 4-1abaixo. A faixa central é chamada primária, enquanto
que as outras são chamadas complementares.
Considere a faixa primária. Se traçarmos a sequência de pontos 1-2-3-4-5-1 no plano s,
conforme a figura 4-2(a) abaixo, esta trajetória é mapeada no plano z conforme a figura
4-2(b). Observe-se que os pontos 3 e 4 abaixo na figura (a) estão em - ∞ , sendo
mapeados na origem no plano z.
2
Na figura abaixo as imagens nos planos s e z para dois valores de σ . Observe-se que
quando a abacissa no plano s é menor que zero, o círculo correspondente tem raio
menor que 1; se estiver à direita, é maior, resultado da expressão z = eTs .
-Uma reta no plano s paralela ao eixo dos reais, com ordenada ω1 é mapeada no plano z
como uma linha radial que parte da origem e com ângulo ω1T no sentido
trigonométrico, como se pode ver na figura 4-5 abaixo.
-Uma região retangular no plano s é mapeada no plano z por uma região limitada por
duas linhas radiais e dois arcos circulares, conforme a figura 4-6 abaixo, observando-se
que o círculo no exterior da figura tem raio unitário.
3
-Uma linha radial no plano s, que corresponde, conforme sabido pela teoria de sistemas
de tempo contínuo, a uma taxa de amortecimento, é mapeada no plano z em uma espiral
que termina na origem. Para demonstrarmos este fato, seja
s = −ςωn + jωn 1 − ς 2 = −ςωn + jωd , onde ωd = ωn 1 − ς 2 : ver figura 4-7(a) abaixo.
No plano z esta linha se torna z = eTs = exp(−ςωnT + jωd T ) .
E tendo em vista a expressão de T em função da frequência de amostragem, vem
⎛
⎛
ω ⎞
ω
2πς ωd
2πς ωd ⎞
⎟ e arg( z ) = 2π d .
+ j 2π d ⎟ , donde z = exp ⎜ −
z = exp ⎜ −
⎜ 1 − ς 2 ωs
⎜ 1 − ς 2 ωs ⎟
ωs
ωs ⎟⎠
⎝
⎝
⎠
Por conseguinte, a magnitude de z diminui e a sua fase aumenta à medida que ωd
cresce, mantidos ς e T constantes, donde que z descreve uma espiral logarítmica,
como mostrado na figura 4-7(b) abaixo.
4
Note-se que para uma dada relação constante ωd / ωs , a magnitude de z torna-se uma
função somente de ς e a fase de z se torna constante.
Por exemplo, se a taxa de amortecimento for 0,3 e se ωd / ωs = 0,25, temos
⎛ 2π × 0,3
⎞
× 0, 25 ⎟ = 0, 61 e arg( z ) = 2π × 0, 25 = 900 .
z = exp ⎜ −
⎜ 1 − 0,32
⎟
⎝
⎠
⎛ 2π × 0,3
⎞
× 0,5 ⎟ = 0,3725 e
Com ωd / ωs = 0,5 , temos z = exp ⎜ −
⎜ 1 − 0,32
⎟
⎝
⎠
0
arg( z ) = 2π × 0,5 = 180 .
Assim, a espiral pode ser parametrizada em função de ωd / ωs (ver fig. 4-7(b)).
4-3 Análise da estabilidade de sistemas em malha fechada no plano z
Considere o seguinte sistema em malha fechada de tempo discreto dado pela sua
C ( z)
G( z)
transformada z:
=
.
(4-3)
R( z ) 1 + GH ( z )
É fácil verificar (faça!) que este sistema em malha fechada (SMF) tem G ( z ) no canal
direto e H ( z ) no canal de realimentação.
A estabilidade, como sabemos, pode ser analisada a partir da equação característica:
P ( z ) = 1 + GH ( z ) = 0 .
Analogamente aos sistemas de tempo contínuo, temos:
- Para que o sistema seja estável, as raízes da eq. característica devem estar no interior
do círculo unitário no plano z. (As raízes da eq. característica são os pólos do SMF).
Qualquer pólo fora do círculo unitário torna o SMF instável.
- Se existir um pólo simples em z = 1, o sistema é criticamente estável (também se diz
marginalmente estável). O sistema é também criticamente estável se existir um par de
pólos conjugados, com grau de multiplicidade um sobre o círculo unitário. Qualquer
pólo múltiplo no círculo unitário torna o sistema instável.
OBSERVACÃO: A estabilidade e a estabilidade crítica (marginal) são chamadas na
literatura especializada de BIBO estabilidade e BIBO estabilidade crítica,
respectivamente, a palavra BIBO significando “bounded input – bounded output”. Ou
seja, se a entrada do SMF for limitada (“bounded”), a saída também o será se o sistema
for BIBO estável, e analogamente para os outros dois casos.
Exemplo 4-2: Considere o SMF da figura 4-13 abaixo. Determine a estabilidade do
sistema quando K = 1. O canal direto, a partir do ganho K tem função de transferência
1 − e− s
1
.
G ( s) =
s s ( s + 1)
5
Solução: Vimos na eq. (3-58) que obtemos
⎡1 − e − s
1 ⎤ 0,3679 z + 0, 2642
(4-4)
G( z) = Z ⎢
⎥=
⎣ s s ( s + 1) ⎦ ( z − 0,3679)( z − 1)
Donde que a eq. característica é:
1 + G ( z ) = 0 ∴ ( z − 0,3679)( z − 1) + 0,3679 z + 0, 2642 = 0 , ou seja, z 2 − z + 0, 6321 = 0 .
Acham-se as raízes desta eq.: z1 = 0,5 + j 0, 6181 e z2 = 0,5 − j 0, 6181 . Donde,
z1 = z2 < 1 , e, portanto o SMF é estável.
Métodos para testar a estabilidade absoluta
Há três métodos para determinar a estabilidade de um sistema sem ter que resolver sua
eq. característica. Estudaremos dois deles. Observe-se que com o progresso dos
computadores e programas confiáveis como MATLAB, MAPPLE e MATHEMATICA,
as raízes de um polinômio podem ser calculadas com facilidade, o que torna menos
importante os métodos a seguir, os quais, entretanto ainda têm valor analítico,
permitindo, por exemplo, estabelecer relações entre coeficientes de um polinômio.
O método de Jury
Seja o polinômio característico
P ( z ) = a0 z n + a1 z n −1 + ... + an −1 z + an , a0 > 0
(Se a0 < 0 , multiplica-se a eq. por -1, o que não lhe altera as raízes).
Constrói-se a Tabela 4-1 abaixo
(4-5)
onde
bk =
an
an −1− k
a0
ak +1
, k = 0;1; 2;...; n − 1
6
ck =
bn −1 bn − 2− k
b0
bk +1
, k = 0;1; 2;...; n − 2
..........................................................
p p2 − k
, k = 0;1; 2
qk = 3
p0 pk +1
O critério de estabilidade de Jury nos diz que o polinômio característico dado em (4-5)
é estável se as seguintes condições forem satisfeitas:
1. an < a0
2. P ( z ) z =1 > 0
⎧ > 0 se n for par
3. P ( z ) z =−1 ⎨
⎩< 0 se n for impar
4. bn −1 > b0 , cn − 2 > c0 , ..., q2 > q0
(Observe-se que o teorema de Jury nos dá condição suficiente, mas não necessária).
Exemplo 4-3: Utilizando o critério de Jury, estabelecer as condições de estabilidade da
seguinte eq. característica: P ( z ) = a0 z 4 + a1 z 3 + a2 z 2 + a3 z + a4 = 0 , a0 > 0 .
Solução: As quatro condições acima nos dão imediatamente:
a4 < a0 , P(1) = a0 + a1 + a2 + a3 + a4 > 0,
P (−1) = a0 − a1 + a2 − a3 + a4 > 0, b3 > b0 , c2 > c0
Claro que b3 , b0 , c2 e c0 têm que ser calculados de acordo com as fórmulas acima.
(Faça-o!)
Exemplo 4-4: Examinar a estabilidade da seguinte equação característica:
P ( z ) = z 4 − 1, 2 z 3 + 0, 07 z 2 + 0,3z − 0, 08 = 0
Solução: Claro que a primeira condição -0,08 < 1 é satisfeita. Quanto à segunda
condição, P (1) = 1 − 1, 2 + 0, 07 + 0, 4 − 0, 08 = 0, 09 > 0 , também é satisfeita.
Terceira condição: P (−1) = 1 + 1, 2 + 0, 07 − 0,3 − 0, 08 = 1,89 > 0, também é satisfeita.
0,3
1
Quarta condição:
= −1, 024 ,
1 −0, 08
−0, 08 −0, 03
b0 =
= −0, 096 + 0, 03 = −0, 066 , a condição sendo satisfeita, pois temos
1
1, 2
b3 > b0 . Etc... Conclui-se que o sistema é estável.
Análise de estabilidade utilizando transformação bilinear e o critério de estabilidade
de Routh
s +1
Considere a transformação bilinear z =
. Esta transformação é dita bilinear, porque
s −1
ela nos dá zs − z = s + 1 , ou seja, uma equação linear em z quando s é constante e linear
7
em s quando z é constante. Desta última obtemos imediatamente s =
z +1
, uma relação
z −1
simétrica, como se vê.
Veremos agora que esta transformação mapeia o interior do círculo unitário no semiplano da esquerda e vice versa. Seja
s = σ + jω .
O interior do círculo unitário no plano z é dado, portanto por
s + 1 σ + jω + 1
(σ + 1) 2 + ω 2
z =
=
< 1 , ou seja,
< 1 ∴ (σ + 1) 2 + ω 2 < (σ − 1) 2 + ω 2
s − 1 σ + jω − 1
(σ − 1) 2 + ω 2
∴ 2σ < −2σ ∴σ < 0 , comprovando que o interior do círculo unitário no plano z é
efetivamente mapeado no semi-plano da esquerda no plano s. A demonstração da
recíproca pode ser feita simplesmente revertendo os passos da prova acima, ou então,
por processo mais trabalhoso, fazendo z = z e jθ com z < 1 , e substituindo na
expressão de s acima, para demonstrar que sua parte real é menor que zero.
Em vista deste resultado, para sabermos se um sistema de tempo discreto é estável,
podemos mapear z em s e daí verificar se as raízes da equação característica mapeada
tem suas raízes no semi-plano da esquerda.
Com efeito, seja
P ( z ) = a0 z n + a1 z −1 + ... + an −1 z + an = 0.
Utilizando a transformação bilinear, temos
n −1
⎛ s +1⎞
⎛ s +1⎞
⎛ s +1⎞
a0 ⎜
⎟ + a1 ⎜
⎟ + ... + an −1 ⎜
⎟ + an = 0.
⎝ s −1 ⎠
⎝ s −1 ⎠
⎝ s −1 ⎠
Multiplicando ambos os lados da eq. por ( s − 1) n , obtemos uma equação do tipo
Q( s ) = b0 s n + b1s n −1 + ... + bn −1s + bn = 0 .
Agora aplicamos o critério de Routh ou o de Hurwitz, estudados em sistemas de tempo
contínuo para concluir sobre a estabilidade, ou não, do sistema. Observe-se que o
critério de Routh, apesar de ser mais trabalhoso do que o de Hurwitz, permite concluir o
número de raízes no semi-plano da direita.
n
4-4 Análise da resposta nos regimes transitório e permanente
Considere-se a figura 4-14 abaixo
Suponha que a entrada seja um degrau unitário. Resposta típica de uma entrada em
degrau em sistemas de tempo contínuo é a que aparece na figura 4-16 abaixo, na qual
são definidos os seguintes parâmetros:
8
- Tempo de retardo (“delay time”) td : definido como o tempo para que a resposta atinja
a metade do valor em regime permanente.
- Tempo de subida (“rise time”) tr : definido como o tempo para que a resposta atinja de
10% a 90% do seu valor final, ou de 5% a 95%, ou de 0% a 100%, dependendo da
situação. Para sistemas sub-amortecidos, a faixa de 10% a 90% é usada comumente.
- Tempo de pico t p : definido como o tempo para que a resposta atinja o pico da
primeira ultrapassagem do valor final.
- Ultrapassagem (“overshoot”) máxima M p : medida acima do valor final,
correspondendo ao pico máximo atingido pela resposta.
- Tempo de assentamento (“settling time”) t s : é o tempo necessário para que a resposta
atinja e permaneça dentro de uma faixa, para mais ou para menos, em torno do valor
final, usualmente 2%. O tempo de assentamento é relacionado à maior constante de
tempo do sistema de controles.
Os parâmetros, ou especificações, acima são importantes, pois não basta apenas saber se
o valor final desejado é atingido. Assim, por exemplo, se M p é muito grande, o sistema
pode “queimar”.
Mas algumas das especificações acima não se aplicam a alguns sistemas: por exemplo,
se o sistema for super-amortecido, o tempo de pico e a ultrapassagem máxima não se
aplicam. E, por outro lado, outras especificações podem se aplicar, como por exemplo,
se houver um erro em regime permanente, é desejável que ele se mantenha dentro de
limites aceitáveis.
A resposta de um sistema a um delta de Kronecker, a um degrau, a uma rampa, etc.
podem ser facilmente obtidas utilizando o MATLAB, como veremos no exemplo
seguinte:
Exemplo 4-8: Seja o sistema de controle de tempo discreto definido por
C ( z)
0, 4673 z − 0,3393
= 2
.
R( z ) z − 1,5327 z + 0, 6607
Obter, utilizando o MATLAB, a resposta do sistema a um degrau unitário.
Abaixo os comandos do MATLAB Program 4-1.
(4-10)
9
O resultado obtido está na Figura 4-17 abaixo.
Análise do erro em regime permanente
No que se segue, investigaremos um tipo de erro em regime permanente, resultado da
incapacidade do sistema de seguir alguns tipos de entrada. Assim, por exemplo, um
sistema pode não apresentar erro em regime permanente quando a entrada é um degrau,
mas sim, se a entrada for uma rampa.
Seja um sistema de tempo contínuo cuja função de transferência em malha aberta é:
K (T s + 1)(Tb s + 1)...(Tm s + 1)
.
G(s) H (s) = N a
s (T1s + 1)(T2 s + 1)...(Tp s + 1)
(Este sistema é dito do tipo N porque tem um pólo na origem com grau de
multiplicidade N).
Como se viu no curso de sistemas de tempo contínuo, um sistema de tipo 0 (zero) dará
uma resposta com erro finito, diferente de zero em regime permanente a um degrau de
entrada, supondo que o sistema seja estável. Sistemas de tipo 1 darão resposta com erro
nulo em regime permanente a um degrau de entrada e erro finito, diferente de zero a
uma rampa de entrada. Aumentar o tipo do sistema pode criar dificuldade na
estabilização do mesmo.
Estas idéias podem ser estendidas para sistemas de tempo discreto da seguinte maneira:
10
Os sistemas de tempo discreto podem ser classificados igualmente conforme seu “tipo”
de acordo com o número de pólos em z = 1 na sua função de transferência em malha
aberta.
Suponha, com efeito, que a função de transferência em malha aberta de um sistema de
1
B( z )
B( z )
tempo discreto seja dada por
, onde
não contêm zero ou pólo em z
N
( z − 1) A( z )
A( z )
= 1. Então o sistema pode ser classificado de tipo 0;1;2;..., de acordo com o valor de N.
E o sentido é o mesmo do de tempo contínuo no que se refere à resposta em regime
permanente do SMF a entradas em degraus, rampas, parábolas, etc., todas em tempo
discreto.
Seja o sistema da figura 4-18 abaixo.
Suponhamos que o mesmo seja estável de modo que o teorema do valor final possa ser
aplicado para achar o valor da resposta em regime permanente.
Do diagrama de blocos, temos e(t ) = r (t ) − b(t ) .
Ora, do teorema do valor final (tabela 2-2, entrada 17), temos
(4-11)
lim e(kT ) = lim[(1 − z −1 ) E ( z )]
k →∞
z →1
⎡ G ( s) ⎤
Definamos G ( z ) = (1 − z −1 )Z ⎢ p ⎥ e
⎣ s ⎦
⎡ G p ( s) H ( s) ⎤
GH ( z ) = (1 − z −1 )Z ⎢
⎥.
s
⎣
⎦
C ( z)
GH ( z )
=
e
Então, do diagrama de blocos,
R( z ) 1 + GH ( z )
1
E( z) =
R( z ) .
(4-12)
1 + GH ( z )
Substituindo-se esta em (4-11), obtemos a expressão do erro em regime permanente
⎡
⎤
1
(4-13)
R( z ) ⎥ .
(steady-state): ess = lim ⎢ (1 − z −1 )
z →1
1 + GH ( z )
⎣
⎦
Tal como no caso de sistemas de tempo contínuo, é comum considerar-se especialmente
três tipos de entradas: degrau unitário, rampa unitária e aceleração unitária.
Degrau unitário:
r (t ) = 1(t ) (trata-se de degrau unitário aplicado em t = 0).
1
.
Da tabela 2-1 temos R( z ) =
1 − z −1
11
⎡
1
1 ⎤
1
Substituindo isto em (4-13), vem: ess = lim ⎢ (1 − z −1 )
= lim
.
−1 ⎥
z →1
z
→
1
1 + GH ( z ) 1 − z ⎦
1 + GH ( z )
⎣
(4-14)
Define-se a constante de erro de posição estática como K p = lim z →1 GH ( z ) .
Consequentemente, o erro em regime permanente devido a uma entrada em degrau
1
.
(4-15)
unitário é ess =
1+ K p
É claro que este erro em regime permanente se anula se só se K p = ∞ , o que implica
que GH ( z ) tenha um pólo em z = 1.
Rampa unitária:
r (t ) = t1(t ) , cuja transformada z é (ver Tabela 2-1):
Tz −1
R( z ) =
. Substituindo esta expressão em (4-13), temos:
(1 − z −1 ) 2
⎡
Tz −1
1
Tz −1 ⎤
−1
ess = lim ⎢(1 − z )
⎥ = lim
z →1 (1 − z −1 )(1 + GH ( z ))
z →1
1 + GH ( z ) (1 − z −1 ) 2 ⎦
⎣
T
1
(1 − z −1 )(1 + GH ( z ))
∴
=
= lim
lim
z →1 (1 − z −1 )(1 + GH ( z ))
ess z →1
T
Define-se a constante de erro de velocidade
(1 − z −1 )(1 + GH ( z ))
K v = lim
.
(4-16)
z →1
T
1
Por conseguinte, ess =
.
(4-17)
Kv
Se K v = ∞ , então o erro em regime permanente é nulo. Mas de (4-16) é claro que isto
implica que GH ( z ) tenha um pólo duplo em z = 1.
Parábola unitária:
1
1
é definida usualmente como r (t ) = t 21(t ) ; o fator é colocado de modo que a sua
2
2
derivada seja a rampa unitária. A transformada z (ver Tabela 2-1) é:
T 2 (1 + z −1 ) z −1
. Substituindo esta em (4-13), temos
R( z ) =
2(1 − z −1 )3
⎡
⎡
⎤
1
T 2 (1 + z −1 ) z −1 ⎤
T2
=
lim
ess = lim ⎢ (1 − z −1 )
⎥
⎢
⎥
1
3
−
−
1
2
z →1
1 + GH ( z ) 2(1 − z ) ⎦ z →1 ⎣ (1 − z ) (1 + GH ( z )) ⎦
⎣
⎡ (1 − z −1 ) 2 (1 + GH ( z )) ⎤
1
∴ = lim ⎢
⎥ . Definindo analogamente, como nos casos anteriores,
ess z →1 ⎣
T2
⎦
(1 − z −1 ) 2 (1 + GH ( z ))
,
z →1
T2
a constante de erro de aceleração K a = lim
(4-18)
1
.
(4-19)
Ka
O erro em regime permanente se torna nulo se só se K a = ∞ , o que implica, de acordo
com (4-18), que GH ( z ) tenha um pólo triplo em z = 1.
temos ess =
12
A Tabela 4-5 abaixo tem as três constantes de erros em regime permanente para
diversos tipos de configuração.
O primeiro e o segundo caso são diferentes, porque, como já vimos, em geral,
GH ( z ) ≠ G ( z ) H ( z ) ; e o 3º. e o 4º. caso são diferentes porque HG2 ( z ) ≠ H ( z )G2 ( z ) .
Resposta a distúrbios
Até agora estudamos a resposta da planta a sinais de referência que se deseja rastreados
em regime permanente com regime transitório satisfatório.
Passemos agora ao estudo da resposta do sistema a distúrbios. Considere a figura 4-19
abaixo. Na figura (a) o sinal de referência é suposto nulo, conforme indicado no próprio
diagrama de blocos. E, portanto o diagrama de blocos pode ser redesenhado como na
figura (b).
13
Então,
C ( z)
G( z)
=
.
N ( z ) 1 + GD ( z )G ( z )
(*)
Ora, se GD ( z )G ( z ) 1 para z < 1 , que é a região de estabilidade, então temos
C ( z)
1
≅
.
(**)
N ( z ) GD ( z )
Ora, remetendo-nos à figura (a), temos E ( z ) = R ( z ) − C ( z ) = −C ( z ) . Disto e da
1
expressão anterior, temos E ( z ) ≅ −
N ( z) .
GD ( z )
Portanto, quanto maior o ganho de GD ( z ) , tanto menor o erro E ( z ) .
Se GD ( z ) inclui um integrador (o que significa que GD ( z ) tem um pólo em z = 1),
então o erro em regime permanente devido a um distúrbio constante é nulo. Com efeito,
A
.
um distúrbio constante tem transformada z igual a N ( z ) =
1 − z −1
G ( z )
, onde G D ( z )
Ora, se GD ( z ) tem um pólo em z = 1, podemos escrever GD ( z ) = D
z −1
não tem zero em z = 1 (isto é, seu numerador não tem raiz em z = 1). Então, o erro em
⎡
− N ( z) ⎤
regime permanente é: ess = lim[(1 − z −1 ) E ( z )] = lim ⎢(1 − z −1 )
⎥
z →1
z →1
GD ( z ) ⎦
⎣
⎛
Az ⎞
⎜ z −1
⎟
z − 1 ⎟ = lim ⎛ z − 1 Az ⎞ = 0.
= lim ⎜
⎜
⎟
z →1
⎜ z G D ( z ) ⎟ z →1 ⎝ z G D ( z ) ⎠
⎜
⎟
⎝
z −1 ⎠
Na realidade, podemos chegar a esta conclusão sem precisar supor a aproximação (**)
acima. Com efeito, de (*), se GD ( z ) inclui um integrador e supondo que ele não seja
cancelado pelo numerador de G ( z ) , obtém-se a mesma conclusão. E note-se que não
cancelamento de raízes instáveis no produto GD ( z )G ( z ) é condição necessária para a
estabilidade do SMF, como é sabido do curso de sistemas de tempo contínuo.
14
Recordamos aqui, da teoria geral de sistemas lineares, que se um sistema é submetido
simultaneamente a um sinal de referência e a um distúrbio, a resposta vai ser a
superposição das duas respostas.
Considere o sistema da figura 4-20(a) abaixo, em que o distúrbio é aplicado na saída da
planta.
C( z)
E ( z)
1
=−
=
.
N ( z)
N ( z ) 1 + GD ( z )G ( z )
Consequentemente, para minimizar o efeito do distúrbio sobre o erro, o ganho de
GD ( z )G ( z ) deve ser tão grande quanto possível.
É fácil verificar que
GD ( z )G ( z )
C( z)
E( z)
=−
=−
,
N ( z)
N ( z)
1 + GD ( z )G ( z )
e neste caso, para minimizar o efeito do distúrbio, o ganho de GD ( z )G ( z ) deve ser tão
pequeno quanto possível.
Entretanto, para a figura (b), temos
É claro então que deve haver um compromisso, quando se trata de minimizar ambos os
efeitos. Usualmente, procura-se o compromisso principalmente na banda de frequências
que interessa.
Ou então, se as freqüências correspondentes ao sinal de referência e ao distúrbio forem
suficientemente disjuntas, as duas condições, isto é, ganhos de GD ( z )G ( z ) pequeno e
grande podem ser satisfeitas satisfatoriamente.
15
4-5 Projeto baseado no lugar das raízes (“root locus”)
O método do root locus para sistemas de tempo contínuo, já estudado no curso anterior,
pode ser aplicado sem alteração, a não ser o fato de que a fronteira entre a parte “boa”
(estável) no plano complexo é o interior do círculo unitário e não o semi plano aberto da
esquerda.
Seja o sistema da figura 4-21 abaixo.
A equação característica do sistema é 1 + G ( z ) H ( z ) = 0.
O MATLAB fornece o root locus, de modo que não precisaremos gastar tempo e
energia nos cálculos. Entretanto, é bom que dominemos o conceito, de resto simples, da
idéia do root locus, de modo a podermos interpretar os resultados obtidos no
computador. Definindo
F ( z ) = G ( z ) H ( z ) , a eq. característica é re-escrita
1 + F ( z) = 0 .
(*)
Observe-se que F ( z ) é a função de transferência em malha aberta. Da eq. acima, temos
F ( z ) = −1 , o que implica
F ( z) = 1 ,
arg[ F ( z )] = ±(2k + 1)π = ± (2k + 1)1800 , k = 1; 2;3;...
Recorda-se que os pólos e zeros de F ( z ) aparecem aos pares conjugados; como
conseqüência, o root locus é sempre simétrico com relação ao eixo real. Portanto, basta
construir a parte superior (acima do eixo dos reais) do root locus, a parte inferior sendo
simétrica.
Regras gerais para construir o lugar das raízes de um SMF
1ª. A partir da eq. (*) acima, escrevemos
K ( z + z1 )( z + z2 )...( z + zm )
1+
= 0.
(**)
( z + p1 )( z + p2 )...( z + pn )
Suporemos, como habitual no root locus, que K > 0.
A partir da forma fatorada acima do numerador e do denominador, plotamos no plano z
os zeros e pólos de F ( z ) .
2ª. Os pontos do root locus correspondentes a K = 0 são os pólos do sistema em malha
aberta, enquanto que os pontos do root locus correspondentes a K = ∞ são os zeros do
sistema em malha aberta. O root locus terá tantos ramos quantas são as raízes da eq.
característica. Cada ramo começa num pólo do sistema em malha aberta e termina num
zero do mesmo. A quase totalidade dos sistemas em controles são próprios, isto é, o
grau do numerador de F ( z ) é menor ou igual ao do denominador. E a maioria deles são
estritamente próprios, ou seja, o grau do numerador é estritamente menor do que o grau
16
do denominador de F ( z ) . Como consequência, neste caso, F ( z ) tem zero(s) no
infinito. E isto significa que o(s) ramo(s) correspondente(s) partindo do(s) respectivo(s)
pólo(s) vai (vão) para o infinito. Ou seja, olhando para a eq. (**) acima, vemos que n –
m ramos do root locus vão para o infinito.
3ª. Determinar os ramos do root locus no eixo dos reais; eles são constituídos por pólos
e zeros do sistema em malha aberta sobre o dito eixo. Para que um ponto do eixo real
pertença a um ramo do root locus, é necessário (e suficiente), que o número total de
zeros e pólos do sistema em malha aberta à sua direita seja ímpar.
4ª. Como dito acima, sendo n > m, há ramos que vão para o infinito; para determinar os
ângulos das assíntotas a estes ramos, usa-se a fórmula
±(2 N + 1) × 1800
Ângulo das assíntotas =
, N = 0;1; 2;...; n − m .
n−m
É claro que cada par de assíntotas é simétrico com relação ao eixo dos reais. O ponto em
que elas se interceptam no dito eixo é dado por
p + p + ... pn − ( z1 + z2 + ...zm )
.
(4-21)
σa = 1 2
n−m
5ª. A seguir, calculam-se os pontos “breakaway” e “break-in”, os quais ou estão no eixo
real ou ocorrem em pares de complexos conjugados. Se existe um ramo do root locus
entre dois pólos no eixo real, então existe um ponto breakaway no dito eixo entre os
dois pólos. E de modo semelhante, se existe um ramo do root locus entre dois zeros no
eixo real (um deles podendo estar no infinito), então existe pelo menos um ponto de
break-in entre os dois zeros.
KB( z )
A( z )
Seja F ( z ) =
. Da eq. (*) acima temos K = −
,
(4-22)
A( z )
B( z )
e os pontos breakaway e break-in são as raízes da eq.
dK
A '( z ) B( z ) − A( z ) B '( z )
=−
=0
(4-23)
dz
B2 ( z)
onde (') indica derivada com relação a z.
Todas as raízes de (4-23) correspondentes a valores positivos de K são pontos
breakaway ou break-in.
6ª. Em seguida, determinam-se os ângulos de partida e de chegada dos pontos
complexos do root locus. Conforme mostrado na figura 4-22, abaixo, o ângulo de
partida (ou de chegada) de um pólo (ou zero) do root locus é igual a 1800 menos a soma
dos ângulos dos vetores que ligam os outros pólos e zeros ao dito pólo (zero). Na
figura, temos um ponto de partida de um pólo indicado com um × no 2º. quadrante; no
caso, o sistema tem três pólos (indicados com × ) e um zero (indicado com uma
bolinha). Esta é a indicação universal para pólos e zeros no root locus.
7ª. A seguir, acham-se os pontos em que o root locus cruza o círculo unitário, o que é
resolvido facilmente fazendo z = 1 na eq. característica (*) ou (**), igualando-se as
partes real e imaginária, obtendo duas eqs. a partir das quais se podem obter z e o
respectivo valor de K.
17
8ª. Como diz o próprio nome, o “root locus” é o lugar das raízes da eq. característica do
SMF, ou seja, é o lugar dos pólos do SMF quando K varia de zero a infinito. Observe-se
que dado um certo K, para que ele satisfaça à eq. característica, terá que valer a relação,
( z + z1 )( z + z2 )...( z + zm )
1
de acordo com (**):
= . Para um dado K, é possível
( z + p1 )( z + p2 )...( z + pn ) K
localizar os pólos do SMF pelo método de tentativa-e-erro.
Cancelamento de pólos e zeros no produto G ( z ) H ( z )
É importante observar que os pólos que forem cancelados no produto G ( z ) H ( z ) não
aparecem no root locus. Assim, por exemplo, seja o SMF da figura 4-21 acima. Se
z+a
z+c
z+c
, e é a partir
G( z) =
e H ( z) =
, obtemos G ( z ) H ( z ) =
( z + a)( z + b)
( z + b)( z + d )
z+d
desta expressão que se constrói o root locus. Os pólos são obtidos da eq. característica
K ( z + c)
1 + G( z) H ( z) = 1 +
= 0, donde ( z + b)( z + d ) + K ( z + c) = 0 nos dão os
( z + b)( z + d )
pólos que aparecem no root locus quando K varia de zero a infinito.
Entretanto, note-se que a função de transferência do SMF é
K ( z + c)
G( z)
K ( z + c)( z + d )
( z + a )( z + b)
=
=
1 + G ( z ) H ( z ) 1 + K ( z + c) z + a ( z + a )( z + b)( z + d ) + K ( z + a )( z + c)
( z + a )( z + b) z + d
K ( z + c)( z + d )
, donde se vê que z = -a é um pólo do SMF, foi o
=
( z + a )[( z + b)( z + d ) + K ( z + c)]
pólo cancelado no produto G ( z ) H ( z ) .
Root locus de controle de sistemas digitais
No que se segue, investigaremos o efeito do ganho K e do período de amostragem T na
estabilidade de um SMF tal como o da figura 4-23 abaixo.
18
K
z
=K
.
−1
1− z
z −1
Vamos construir o root locus para três valores do período de amostragem T: 0,5; 1 e 2
segundos. Vamos determinar também o valor crítico de K para cada caso. E, finalmente,
vamos calcular os valores críticos dos pólos do SMF (isto é, a partir dos quais o SMF se
torna instável), para K = 2 em cada um dos três casos.
Calculemos primeiramente a transformada z de Gh ( s )G p ( s ) :
Suponha que o controlador digital é do tipo integral, ou seja, GD ( z ) =
⎡1 − e −Ts 1 ⎤
⎡ 1 ⎤
1 ⎤
⎡1
−1
−1
Z[Gh ( s )G p ( s )] = Z ⎢
⎥ = (1 − z )Z ⎢ s ( s + 1) ⎥ = (1 − z )Z ⎢ s − s + 1 ⎥ =
⎣
⎦
⎣
⎦
⎣ s s + 1⎦
−T
−T
z −1 ⎛ z
z ⎞
z −1
z − e − z +1 1− e
=
=
.
−
= 1−
⎜
−T ⎟
−T
z−e
z − e −T
z − e −T
z ⎝ z −1 z − e ⎠
Donde que a função de transferência de pulso do canal direto, que é a do sistema em
z 1 − e −T
malha aberta, é: G ( z ) = GD ( z )Z[Gh ( s )G p ( s )] = K
.
(4-25)
z − 1 z − e −T
z 1 − e −T
A equação característica é 1 + G ( z ) = 0 , ou seja, 1 + K
=0
(4-26)
z − 1 z − e −T
- Consideremos primeiramente T = 0,5 seg. A eq. (4-25) nos dá
0,3935 Kz
G( z) =
.
( z − 1)( z − 0, 6065)
Para construirmos o root locus, que está na figura 4-24 (a) abaixo, plotamos
primeiramente os pólos e zeros de G ( z ) .
19
Os pólos são indicados com × e o zero com um pequeno círculo, todos, neste caso, no
eixo dos reais. Usualmente se indica no root locus a evolução dos pólos, à medida que K
aumenta, setas indicando esta evolução, o que não é feito na figura.
Com K = 0, os pólos do SMF são os pólos de G ( z ) . O root locus tem dois ramos,
porque tem dois pólos. O ramo que começa em z = 1 caminha para a esquerda quando K
aumenta, enquanto que o ramo que começa em z = 0,6065 caminha para a direita
encontrando o outro no ponto “breakaway”, que é calculado, como vimos da eq. (4-22).
Os dois ramos se afastam simetricamente em relação ao eixo dos reais e se encontram
no ponto “break-in”, ambos estes pontos sendo obtidos pela mesma maneira:
( z − 1)( z − 0, 6065)
K =−
(4-27)
0,3935 z
e igualando a zero a derivada com relação a z:
dK
0,3935 z ( z − 0, 6065 + z − 1) − 0,3935( z − 1)( z − 0, 6065)
=−
=
dz
(0,3935 z ) 2
2 z 2 − 1, 6065 z − z 2 + 1, 6065 z − 0, 6065 z 2 − 0, 6065
=
= 0 , donde obtemos z = ±0, 7788 ;
0,3935 z 2
0,3935 z 2
o valor positivo é o ponto breakaway, enquanto que o negativo é o ponto “break-in”.
20
Substituindo estes valores em (4-27), achamos K = 0,1244 e K = 8.041,
respectivamente.
O root locus neste caso tem uma parte circular, como mostra a figura; na realidade, para
chegarmos a esta conclusão, precisamos levantar ponto por ponto, isto é, para cada valor
de K, determinar os pólos do SMF. Mas que esta curva seja um círculo, também pode
ser demonstrado analiticamente. (Tente!).
Como se vê na figura (a), o ramo da esquerda sobre o eixo real vai para o infinito,
saindo fora do círculo unitário, ou seja, o sistema torna-se instável a partir de um certo
valor de K. Para determinar este valor de K, utilizamos a eq. (4-26), que repetimos:
z 1 − e −T
( z − 1)( z − 0, 6065)
1+ K
= 0 ; substituindo T = 0,5 , obtemos K = −
e com z
−T
z −1 z − e
0,3935 z
(−2)(−1, 6065)
= -1, temos K = −
= 8,165.
−0,3935
- Vejamos agora o root locus quando T = 1. Reportamo-nos à figura (b) acima. Neste
0, 6321Kz
caso, da eq. (4-25), temos G ( z ) =
.
( z − 1)( z − 0,3679)
Calculando como no caso anterior a partir da eq., mas com os valores diferentes, achamse os pontos breakaway e break-in, a saber, z = 0,6065 e z = - 06065, respectivamente,
ou seja, simétricos com relação à origem. Os valores correspondentes de K são K =
0,2449 e 4,083. Como se poderia esperar, tendo em vista o caso anterior, ele é circular.
Utilizando o mesmo método que no caso anterior, calcula-se o valor crítico de K a partir
do qual o SMF se torna instável, obtendo-se K = 4,328.
- Passemos finalmente ao caso em que T = 2, reportando-nos à figura (c) acima. Neste
0,8647 Kz
caso, a eq. (4-25) nos dá: G ( z ) =
. Calculando os pontos de
( z − 1)( z − 0,1353)
breakaway e de break-in, temos (ver figura(c)) z = 0,3678 e z = -0,3678,
respectivamente, simétricos com relação à origem. Obtém-se o valor crítico de K a
partir do qual o SMF se torna instável: K = 2,626.
Observe-se que os root loci diminuem de diâmetro quando se aumenta o período T de
amostragem.
E verifica-se que, quanto maior o período de amostragem, menor se torna também o
valor de K a partir do qual o SMF se torna instável, ou seja, quanto menor a frequência
de amostragem, menor a margem de estabilidade do SMF.
Uma “receita de bolo” que se usa para o valor do período de amostragem é a seguinte: o
período de amostragem deve ser 8 ou 10 vezes menor do que o das oscilações na
resposta do sistema, se o sistema for sub-amortecido. Se o sistema for super-amortecido,
o período de amostragem deve ser 8 ou 10 vezes menor do que o do tempo de subida no
caso de resposta a um degrau.
Efetivamente, fazendo-se o tempo de amostragem cada vez menor, o que se torna
sempre mais possível com a “lei de Moore” (até quando?), o sistema de tempo discreto
tem um comportamento que se aproxima do de tempo contínuo.
Efeitos do período de amostragem T nas características do regime transitório:
Se o fator de amortecimento do SMF é ς , então o pólo correspondente no plano s acima
do eixo real é s = −ςωn + jωn 1 − ς 2 e sendo z = eTs , temos
(
)
z = exp T (−ςωn + jωn 1 − ς 2 ) ,
(4-29)
21
donde z = e −T ςωn
e arg[ z ] = T ωn 1 − ς 2 = T ωd = θ (em radianos).
(4-30)
Assim, se T = 0,5 e K = 2, é fácil verificar que o pólo do SMF com parte imaginária
positiva é: z = 0,4098 +j0,6623. Portanto, z = 0, 40982 + 0, 66232 = 0,7788, donde
e −T ςωn = 0,7788, ou seja,
T ςωn = 0,25.
0, 6623
= 58, 25D = 1, 0167 rad .
0, 4098
T ςωn
0, 25
=
, donde
Desta, de (4-30) e de (4-31), vem
2
1, 0167
T ωn 1 − ς
nos dá ς = 0,2388.
(4-31)
Por outro lado, arg[ z ] = arctg
(4-32)
ς
1− ς 2
= 0, 2459 , que
Voltemos à figura 4-23, cuja função de transferência de pulso do canal direto é G ( z ) e
G( z)
cuja função de transferência de pulso do SMF é
, onde G ( z ) é dado pela eq.
1 + G( z)
G( z)
Kz (1 − e −T )
(4-25), ou seja,
.
=
1 + G ( z ) ( z − 1)( z − e −T ) + Kz (1 − e −T )
Com T = 0,5 e K = 2, temos a transformada da resposta a um degrau unitário de
0,3935 × 2 z
entrada: C ( z ) =
R( z ) . Dividindo-se numerador e
( z − 1)( z − 0, 6065) + 0,3935 × 2 z
denominador da expressão por z 2 e tendo em vista a tabela 2-1 para a transformada do
0, 787 z −1
1
. A figura 4-26 (a)
degrau unitário, temos C ( z ) =
−1
−2
1 − 0,8195 z + 0, 6065 z 1 − z −1
abaixo dá a resposta do sistema.
22
Quando T = 1 e o mesmo K (= 2), temos a resposta ao degrau unitário:
1, 2642 z −1
1
, o resultado sendo mostrado na fig. (b) acima.
C ( z) =
−1
−2
1 − 0,1037 z + 0,3679 z 1 − z −1
Quando T = 2, com o mesmo K, temos: C ( z ) =
1, 7294 z −1
1
,o
−1
−2
1 + 0,5941z + 0,1353 z 1 − z −1
resultado estando na fig. (c).
Nas três figuras (a), (b) e (c), os pontos escuros correspondem a c(kT) e o tracejado é
uma interpolação. Comparando-se as figs. (a) e (c), fica evidente que a primeira é muito
mais precisa, ou seja, quanto menor o período de amostragem, melhor será a fidelidade
da resposta em tempo discreto. Se T → 0 , temos a resposta em tempo contínuo, da qual
nos aproximamos à medida que os processadores se tornam mais rápidos.
Nos três casos acima o teorema da amostragem (Shannon) é satisfeito, mas a terceira
resposta (c) é pouco satisfatória. Isto mostra, que não basta satisfazer ao dito teorema. A
“receita de bolo” aqui é a mesma já enunciada antes: 8 a 10 amostras por ciclo, se o
sistema é sub-amortecido, como é o caso acima, pois apresenta oscilação na resposta.
A seguir, investigaremos o efeito do período de amostragem no erro em regime
permanente. Consideraremos a rampa unitária como entrada nos três casos tratados. De
(4-25) com T = 0,5 e K =2, temos a função de transferência de pulso do canal direto:
0, 787 z
G( z) =
. Portanto a constante de erro de velocidade é dada por
( z − 1)( z − 0, 6065)
23
0, 787 z
⎡ z −1
⎤
= lim ⎢
⎥ = 4; portanto, o erro em regime
z
→
1
T
⎣ 0, 5 z ( z − 1)( z − 0, 6065) ⎦
1
= 0, 25 .
permanente é: ess =
Kv
1, 2642 z
Com T = 1 e K =2, temos G ( z ) =
, donde que a constante de erro de
( z − 1)( z − 0,3679)
K v = lim
(1 − z −1 )G ( z )
z →1
1, 2642 z
⎡ z −1
⎤
= lim ⎢
⎥=2
z →1
z →1
T
⎣ z ( z − 1)( z − 0, 3679) ⎦
e, portanto, o erro em regime permanente é igual a 0,5.
velocidade é K v = lim
(1 − z −1 )G ( z )
1, 7294 z
e pelo mesmo
( z − 1)( z − 0,1353)
procedimento nos casos anteriores obtemos o erro em regime permanente igual a 1.
A figura 4-27 abaixo mostra as três situações e, de novo, quanto menor o período de
amostragem, melhor a resposta.
Finalmente, com T = 2 e o mesmo K, temos G ( z ) =
Exemplo 4-9: Considere o SMF da figura 4-28 abaixo. Projetar um controlador digital
tal que os pólos dominantes no SMF tenham um fator de amortecimento ς = 0,5 e um
24
tempo de assentamento igual a 2 segundos. O período de amostragem é T = 0,2 seg.
Obter a resposta do sistema a um degrau de amplitude unitária, bem como a constante
de velocidade K v do sistema.
Solução: para sistemas padrão com um par de pólos dominantes, temos, conforme
estudado no curso de sistemas de tempo contínuo: tempo de assentamento = 2 =
4
4
=
.
ςωn 0,5ωn
(Aqui o tempo de assentamento está sendo definido como o tempo a partir do qual o
erro da resposta, em valor absoluto, é inferior a 2% . Se se quisesse ficar dentro de 5%,
o numerador da expressão acima seria igual a 3).
Donde que ωn = 4, que é a frequência natural não amortecida. Mas a frequência natural
real é a amortecida, isto é, ωd = ωn 1 − ς 2 = 4 1 − 0,52 = 3,464.
2π 2π
Já a frequência de amostragem é dada por ωs =
=
= 31,42.
0, 2
T
Observe-se que dividindo este número pelo anterior, isto é, 31,42/3,464 , obtemos
aproximadamente 9,07, ou seja, cerca de 9 amostragem por ciclo de oscilação
amortecida, um número que satisfaz à “regra de bolo” antes mencionada, o que nos leva
a concluir que o período de amostragem igual a 0,2 é satisfatório.
((
Repetimos aqui a eq. (4-29): z = exp T −ςωn + jωn 1 − ς 2
)) ; tendo em vista as
expressões de ωn e ωs , dadas logo acima, temos:
⎛ 2π
⎛
ωd ⎞
2π
2πς ωd ⎞
ωd .
⎟ = exp ⎜ −
⎟ e arg[ z ] = T ωd =
z = e −T ςωn = exp ⎜ −
ς
2 ⎟
⎜ ωs
⎜ 1 − ς 2 ωs ⎟
ω
ς
−
1
s
⎝
⎠
⎝
⎠
Ora, do enunciado do problema temos que ς = 0,5 e, como vimos, ωd = 3,464 e ωs =
31,42 . Substituindo-se estes valores acima, temos:
⎛ 2π × 0,5 3, 464 ⎞
z = exp ⎜ −
⎟ = e −0,4 = 0,6703 e
⎜ 1 − 0,52 31, 42 ⎟
⎝
⎠
2π
× 3, 464 = 0,6927 rad. = 39, 69D .
arg[ z ] =
31, 42
Portanto, z = 0, 6703arg[39, 69D ] = 0,5158 + j 0, 4281 .
A função de transferência de pulso da planta precedida do segurador é
⎡1 − e −0,2 s
⎡
⎤
1
1 ⎤
−1
G( z) = Z ⎢
⎥ = (1 − z )Z ⎢ s 2 ( s + 2) ⎥ . Decompondo em frações parciais
s ( s + 2) ⎦
⎣
⎦
⎣ s
25
o que está entre colchetes, usando a Tabela 2-1 e calculando, obtemos:
0, 01758( z + 0,876)
G( z) =
. Vamos, tentativamente, escolher um compensador GD ( z )
( z -1)( z - 0, 6703)
K ( z + b)
. Ora, se G ( z ) tiver um pólo estável, ele pode ser
para o sistema. Seja GD ( z ) =
z+a
K ( z − 0, 6703)
. Então a eq. característica do
cancelado. Então escolhamos GD ( z ) =
z+a
SMF será ( z − 1)( z + a) + 0, 01758 K ( z + 0,876) . Os parâmetros a e K serão escolhidos
de modo a obtermos uma resposta satisfatória. Com efeito, sejam -c e -d os pólos
desejados, de modo que o polinômio. característico desejado seja ( z + c)( z + d ) =
z 2 + (c + d ) z + cd .
(*)
Com a escolha de GD ( z ) acima, o polinômio característico do sistema é
( z − 1)( z + a) + 0, 01758K ( z + 0,876) = z 2 + (a − 1 + 0, 01758 K ) z − a + 0, 01758 × 0,876 K
= z 2 + (a − 1 + 0, 01758 K ) z − a + 0, 0154 K .
Igualando os coeficientes deste com os de (*) acima, temos a − 1 + 0, 01758K = c + d ,
−a + 0, 0154 K = cd . Somando membro a membro estas duas eqs., temos
c + d + cd + 1
0, 03298K − 1 = c + d + cd ∴ K =
e, por conseguinte, a = 0, 0154 K − cd . É
0, 03298
claro que c e d devem ser escolhidos de modo que as raízes de (*) fiquem no interior do
círculo unitário. Se escolhermos, por exemplo, um pólo duplo na origem, temos de (*),
c = d = 0 e, portanto, K = 1/ 0,03298 = 30,32 e a = 0,0154K = 0,4669. Com isto, temos
o controlador definido e o cálculo da resposta a um degrau unitário é rotina. Observe-se
que G(z) tem um pólo em 1 e, portanto, aplicando o teorema do valor final, obteremos
erro nulo em regime permanente.
4-6 Projeto com o método da resposta em frequência
Os conceitos da resposta em frequência têm a mesma importância em sistemas de tempo
discreto que em sistemas de tempo contínuo.
Supõe-se neste curso que o aluno já esteja familiarizado com os métodos de resposta em
frequência dos sistemas de tempo contínuo.
Sabemos que a resposta de um sistema linear a um sinal senoidal preserva a frequência
e modifica somente a amplitude e a fase do sinal de entrada.
Resposta de um sistema de tempo discreto a uma função senoidal de entrada.
A resposta em frequência de G ( z ) pode ser obtida substituindo-se z = e jωT em G ( z ) .
Seja o sistema de tempo discreto, que suporemos estável na figura 4-31abaixo:
O sinal de entrada é u (t ) = senωt , conforme indicado acima.
26
Então o sinal amostrado é u (kT ) = sen(kωT ) , cuja transformada z é:
zsenωT
U ( z ) = Z[ senkωT ] =
, conforme a entrada 14 da tabela 2-1. Então, a
jωT
( z − e )( z − e − jωT )
zsenωT
.
resposta do sistema é dada por X ( z ) = G ( z )U ( z ) = G ( z )
jωT
( z − e )( z − e − jωT )
az
az
+ termos devidos aos pólos
Obtemos uma expressão do tipo X ( z ) =
+
jωT
z −e
z − e − jωT
de G ( z ) , onde a e a são complexos conjugados.
(4-34)
jωT
Multiplicando ambos os lados desta eq. por ( z − e ) / z , obtemos
senωT
a ( z − e jωT ) z − e jωT
=
a
+
+
[termos devidos aos pólos de G ( z ) ].
z − e − jωT
z
z − e − jωT
Ora, quando z tende para e jωT , a segunda e a terceira parcelas do lado direito acima
tendem para zero, tendo em vista que G ( z ) é estável.
senωT
senωT
Donde, a = G ( z )
= G ( jωT ) jωT
− jωT
z −e
e − e − jωT
z = e jωT
G( z)
senωT
G (e jωT )
G (e − jωT )
e, portanto, a = −
.
=
cos ωT + jsenωT − cos ωT + jsenωT
2j
2j
Definamos G (e jωT ) = Me jθ ∴ G (e − jωT ) = Me − jθ .
Consequentemente, a eq. (4-34) pode ser escrita como
Me jθ
z
Me − jθ
z
−
+ termos devidos aos pólos de G ( z ) , ou ainda,
X ( z) =
jωT
2 j z −e
2 j z − e − jωT
= G ( jωT )
M ⎛ e jθ z
e − jθ z ⎞
−
⎜
⎟ + termos devidos aos pólos de G ( z ) .
2 j ⎝ z − e jωT z − e − jωT ⎠
A transformada inversa desta última é:
M jkωT jθ
(e e − e − jkωT e − jθ ) + Z−1[ termos devidos aos pólos de G ( z ) ] .
x(kT ) =
2j
O último termo do lado direito representa o regime transitório; sendo G ( z ) estável,
estes termos tenderão a zero. Portanto a resposta em regime permanente é:
M jkωT jθ
M j ( kωT +θ ) − j ( kωT +θ )
−e
(4-35)
xss (kT ) =
e e − e − jkωT e− jθ ) =
(
(e
)
2j
2j
M
= ( cos(kωT + θ ) + jsen(kωT + θ ) − cos(kωT + θ ) + jsen(kωT + θ ) )
2j
= Msen(kωT + θ ) ,
(4-36)
onde, como vimos, M é o ganho do sistema de tempo discreto quando submetido a uma
entrada senoidal e é dado por
M = M (ω ) = G (e jωT ) e
X ( z) =
(
)
θ = θ (ω ) = arg G ( e jωT ) .
Portanto, a eq. (4-36) pode ser escrita como
(
(
))
xss (kT ) = G (e jωT ) sen kωT + arg G ( e jωT ) .
Acabamos de provar que G (e jωT ) dá a magnitude e a fase da resposta em frequência de
G ( z ) . Portanto, para obter a resposta em frequência de G ( z ) , basta substituir z por
27
e jωT em G ( z ) . A função G (e jωT ) é comumente chamada de função de transferência
de pulso senoidal.
Observe-se que e j (ω + (2π / T ))T = e jωT e j 2π = e jωT , donde se conclui que a função G (e jωT )
é periódica, com período igual a T.
Exemplo 4-10: Considere o sistema dado por
x(kT ) = u (kT ) + ax((k − 1)T ) com 0 < a < 1 . Obter a resposta em regime permanente
quando a entrada u (kT ) é a senoide amostrada, ou seja, u (kT ) = Asen(kωT ) .
Solução: A transformada z da eq. que define o sistema é:
X ( z ) = U ( z ) + az −1 X ( z ) .
1
Definindo G ( z ) = X ( z ) / U ( z ) , temos G ( z ) =
. Substituindo z por e jωT em
−1
1 − az
1
1
G ( z ) , temos G (e jωT ) =
=
. Donde,
− jωT
1 − ae
1 − a cos ωT + jasenωT
1
asenωT
G (e jωT ) = M =
.
e arg[G (e jωT )] = θ = − arctg
1 − a cos ωT
1 + a 2 − 2a cos ωT
E, portanto, a resposta em regime permanente é:
A
asenωT ⎞
⎛
xss (kT ) = AMsen( kωT + θ ) =
sen ⎜ kωT − arctg
⎟
2
1 − a cos ωT ⎠
⎝
1 + a − 2a cos ωT
Transformação bilinear para o plano z
A fim de aplicar os métodos de frequência a sistemas de tempo discreto, uma
transformação de variável se faz necessária. Com efeito, tendo em vista que no plano z
temos z = e jωT , se tratarmos da resposta em frequência, a simplicidade dos diagrama de
Bode, que é logarítmico será perdida. Mais ainda, como o semi-plano da esquerda é
mapeado no interior do círculo unitário, os métodos frequenciais, que tratam do semiplano inteiro, não se aplicam ao plano z . Tendo em vista isto, uma mudança de variável
se impõe se quisermos usar o diagrama de Bode para sistemas de tempo discreto.
A transformação usada é:
1 + (T / 2) w
z=
,
(4-37)
1 − (T / 2) w
onde T é o período de amostragem. Da eq. acima, temos imediatamente
2 z −1
w=
.
(4-38)
T z +1
2 z −1
são mostrados na figura 4-32 abaixo.
Os mapeamentos z = eTs e w =
T z +1
28
Observe-se que na figura (a), temos uma faixa no semi-plano da esquerda, ao passo que
na figura (c), temos todo o semi-plano da esquerda. Na figura (b) o círculo é o unitário.
Através do mapeamento (4-38), G ( w) pode ser tratada como uma função de
transferência em w. Substituindo w por jν , técnicas convencionais de resposta em
frequência podem ser usadas para traçar o diagrama de Bode.
Faremos agora uma breve revisão do diagrama de Bode nesta seção, usando a
frequência ν como variável. A frequência fictícia ν e a frequência atual ω são
2 z −1
2 e jωT − 1 2 e jωT / 2 e jωT / 2 − 1
relacionadas por: w w= jν = jν =
=
=
; dividindo
T z + 1 z =e jωT T e jωT + 1 T e jωT / 2 e jωT / 2 + 1
numerador e denominador por e jωT / 2 , temos: jν =
2 e j (1/ 2)ωT − e − j (1/ 2)ωT 2
ωT
= jtg
, ou
− j (1/ 2)ωT
j (1/ 2)ωT
Te
+e
T
2
2 ωT
,
(4-39)
tg
2
T
que dá a relação entre a frequência atual ω e a frequência fictícia ν . Observe-se de (439) e também da figura 4-32 que, enquanto ω varia de -0,5 ωs a 0, a frequência fictícia
ν varia de −∞ a 0, e enquanto a primeira varia de 0 a 0,5 ωs , a segunda varia de 0 a ∞ .
A figura 4-33 abaixo mostra a relação entre a frequência fictícia ν multiplicada por T/2
e a frequência atual ω , quando esta varia de 0 a 0,5 ωs .
seja, ν =
29
Seja observado que quando ω T é pequeno, então tg
ωT
≅
ωT
e, de (4-39), temos
2
2
ν ≅ ω . Isto mostra que, para ω T pequeno, as funções de transferência G ( s) e G ( w) se
parecem. Note-se que esta é uma consequência da introdução do fator de escala 2/T em
(4-37). Isto significa que a função de transferência no plano w se aproxima da no plano
s quando T se aproxima de zero, algo sempre mais realista à medida que avança a
micro-eletrônica, de acordo com a lei de Moore.
Exemplo 4-11: Considere a função de transferência da figura 4-34 abaixo. O período de
amostragem é 0,1 seg. Calcular G ( w) .
⎡ 10 ⎤
⎡1 − e −Ts 10 ⎤
= (1 − z −1 )Z ⎢
Solução: A transformada z de G ( s ) é G ( z ) = Z ⎢
⎥
⎥
⎣ s s + 10 ⎦
⎣ s ( s + 10) ⎦
0, 6321
=
. Substituindo z dado por (4-37), obtemos, com T = 0,1
z − 0,3679
0, 6321
0, 6321(1 − 0, 05w)
(1 − 0, 05w)
G ( w) =
=
= 9, 241
.
1 + 0, 05w
0,
6321
+
0,
0684
w
w
+
9,
241
− 0,3679
1 − 0, 05w
Observe-se que a planta G ( s) tem um pólo em s = -10, enquanto que no plano w, o pólo
está em -9,241, não muito distante; quanto ao ganho, o da planta é 10 e o no plano w é
também 9,241, de novo não muito distante. Mas, em compensação, em G ( w) aparece
um zero em 20 = 2/T, o que não ocorre na planta. Mas, efetivamente, quando T tende a
zero, o zero em w = 2/T, se aproxima do infinito, que é um zero da planta G ( s) .
Diagramas de Bode
O diagrama de Bode, que nos dá o ganho e a fase, é largamente usado em sistemas
escalares (uma entrada e uma saída) de tempo contínuo, sendo de mais fácil aplicação
quando a função de transferência é dada na forma fatorada.
Tal como nos sistemas de tempo contínuo, o diagrama relativo nos dá o logaritmo de
G ( jν ) versus o logaritmo de ν . Como o logaritmo de um produto é a soma dos
logaritmos de fatores, a construção do diagrama de Bode fica simplificada se a função
de transferência for dada de forma fatorada.
Usando-se o diagrama de Bode, o controlador digital pode ser projetado usando os
métodos bem conhecidos usados em sistemas de tempo contínuo. Mas é importante
observar que pode haver diferença entre as magnitudes de G ( w) e a de G ( jν ) para
altas freqüências.
1 − 0, 05w
. O ganho de
Vejamos o exemplo anterior (4-11). Repetimos G ( w) = 9, 241
w + 9, 241
30
alta frequência é lim G ( jν ) = lim 9, 241
ν →∞
planta em alta frequência é
ν →∞
lim
10
ω →∞ jω +10
1 − 0, 05 jν
= 0,4621, enquanto que o ganho da
jν + 9, 241
= 0.
Há, portanto, uma diferença considerável. Mas é preciso não esquecer que ν = ∞ no
plano w corresponde a ω = 0,5ωs no plano s.
Vantagens do diagrama de Bode em projetos
São as seguintes as vantagens do diagrama de Bode para fins de projeto:
1. No diagrama de Bode a assíntota em baixas freqüências da curva de ganho é
indicativa de uma das constantes de erro estático K p , K v ou K a .
2. As especificações do regime transitório podem ser traduzidas na resposta em
frequência em termos de margem de ganho, margem de fase, banda passante, etc. Estas
especificações podem ser facilmente tratadas no diagrama de Bode. De modo particular,
as margens de ganho e de fase podem ser medidas diretamente no diagrama de Bode.
3. O projeto de um compensador (ou controlador) digital para satisfazer às
especificações de margem de ganho e margem de fase pode ser feito diretamente por
meio do diagrama de Bode de uma maneira simples e imediata.
Compensação por avanço de fase, atraso de fase e por avanço-atraso de fase.
- Recorda-se que a compensação por avanço de fase (“lead compensator”) aumenta a
banda passante do sistema, o sistema ficando com uma resposta mais rápida. Porém, tal
sistema pode ficar sujeito a ruídos de alta frequência, devido ao seu alto ganho para
freqüências maiores.
- Por sua vez, o sistema com compensação de atraso de fase (“lag compensator”) reduz
o ganho do sistema para altas freqüências, sem reduzir o ganho em baixas freqüências.
A banda passante do sistema é reduzida, e assim o sistema tem uma resposta mais lenta.
Em vista da redução do ganho em alta frequência, o ganho total do sistema pode ser
aumentado e assim o ganho em baixa frequência pode ser aumentado, aumentando
conseqüentemente a precisão do regime permanente. Além disso, qualquer ruído de alta
frequência é atenuado.
- Em algumas aplicações, um compensador de atraso de fase é colocado em cascata
(“em série”) com um compensador de avanço de fase, obtendo-se o compensador “lead
– lag”. Com este compensador, o ganho de baixa frequência pode ser aumentado (o que
significa um aumento na precisão da resposta em regime permanente) ao mesmo tempo
em que se consegue melhor banda passante e margens de estabilidade. Recorda-se
também que o controlador PID é um caso especial de compensador “lead-lag” (ou “laglead”, como dizem alguns autores). Com efeito, o controlador PD tem o mesmo
comportamento que o compensador “lead”, enquanto que o compensador PI se
comporta como um controlador “lag”.
Algumas observações a respeito do problema da quantização dos coeficientes
Do ponto de vista da implementação de micro-processadores, os compensadores “lead”
não apresentam problema especial de quantização, porque os pólos e zeros podem ser
bem separados uns dos outros; mas este não é o caso nos controladores “lead-lag”,
porque os zeros e pólos ficam perto uns dos outros, na realidade ficam próximos do
ponto z = 1.
Tendo em vista que os coeficientes do filtro devem ser realizados através de sistema
binário, se o número de bits empregado for insuficiente, as localizações dos zeros e
31
pólos do filtro podem não ser realizadas como desejado e assim o filtro pode não ser
realizado como desejado.
Ora, pequenos desvios na localização de pólos e zeros podem acarretar diferenças
significativas nas características da resposta em frequência do compensador e assim o
compensador pode não se comportar como desejado. Para minimizar este efeito da
quantização dos coeficientes, é necessário projetar o filtro de modo de modo que ele
seja pouco sujeito a este efeito indesejável.
Como a sensibilidade das raízes dos polinômios às variações dos parâmetros variam em
proporção direta com o grau do polinômio, não é recomendável a realização direta de
um filtro de ordem grande. (Lembra-se que a ordem de uma função de transferência é
igual ao grau do seu denominador). É preferível usar elementos de ordem baixa em
cascata ou em paralelo, como discutido na seção 3-6. Efetivamente, se escolhermos logo
de saída os pólos e zeros do compensador separados uns dos outros, o problema pode
ser evitado.
Nos compensadores analógicos seus pólos e zeros podem ser posicionados com
precisão. Ao convertermos o compensador analógico em digital, a versão digital do
compensador “lag” pode envolver bastante imprecisão na localização dos pólos e zeros.
(A coisa importante a ser lembrada é que os pólos e zeros do filtro no plano z devem
permanecer em pontos bem separados uns dos outros).
Procedimentos de projeto no plano w
Reportemo-nos à figura 4-35 abaixo.
Os procedimentos para um projeto no plano z são os seguintes:
1. Obter G ( z ) , que é a transformada z da planta precedida de um segurador (“hold”). A
1 + (T / 2) w
, isto é,
seguir, calcula-se G ( w) através da eq. (4-37): z =
1 − (T / 2) w
G ( w) = G ( z ) z =[1+ (T / 2) w]/[1−(T / 2) w] . É fundamental que o período de amostragem T seja
escolhido apropriadamente. Uma “regra de bolo” é amostrar com uma frequência 10
vezes maior do que a banda passante do sistema em malha fechada. Há que notar que
embora o controle digital e o processamento de sinais usem procedimentos semelhantes
na amostragem de sinais de tempo contínuo, as freqüências de amostragem são muito
diferentes. No processamento de sinais, estas freqüências são muito altas, enquanto que
em controle digital elas são, em geral, baixas. Esta diferença entre as freqüências de
amostragem nas duas áreas se deve à diferença dos processos dinâmicos envolvidos e
aos diferentes compromissos (“trade-offs”) nas duas áreas).
2. Fazer w = jν e obter o diagrama de Bode de G ( jν ) .
3. Medir, a partir do diagrama de Bode, as constantes de erro estático, a margem de fase
e a margem de ganho.
4. Supondo que o ganho de baixa frequência do controlador de tempo discreto (ou
32
digital) seja igual a um, determinar o ganho do sistema, satisfazendo às exigências para
se obter um erro estático dado. A seguir, usando técnicas convencionais de projeto para
sistemas de tempo contínuo, determinar os pólos e zeros da função de transferência do
controlador GD ( w) . A função de transferência em malha aberta (canal direto do SMF) é
GD ( w)G ( w) .
5. Transformar a função de transferência GD ( w) em GD ( z ) através de (4-38), a saber,
2 z −1
w=
. GD ( z ) é a função de transferência de pulso do controlador digital.
T z +1
6. Realizar a função de transferência GD ( z ) através de um algoritmo computacional.
Observações:
- G ( w) é uma função de transferência de fase não mínima (isto é, tem pelo menos um
zero na parte “ruim” do plano complexo).
– O eixo horizontal no diagrama de Bode, o das freqüências, é distorcido no plano w.
Exemplo 4-12: Considere o sistema de controle digital dado na figura 4-36 abaixo.
Projete um controlador digital no plano w com margens de ganho de pelo menos 10 dB
e margem de fase de 50 0 e uma constante de erro estático de velocidade K v = 2 / seg. O
período de amostragem é T = 0,2.
Solução:
⎡1 − e −0,2 s K ⎤
⎡ K ⎤
K ( z + 0,9356)
−1
= 0, 01873
G( z) = Z ⎢
⎥ = (1 − z )Z ⎢ 2
⎥
( z − 1)( z − 0,8187)
s ( s + 1) ⎦
⎣ s ( s + 1) ⎦
⎣ s
K (0, 01873z + 0, 01752)
= 2
.
z − 1,8187 z + 0,8187
1 + (T / 2) w 1 + 0,1w
, obtendo
=
A seguir, substituímos na expressão acima z =
1 − (T / 2) w 1 − 0,1w
⎡
⎤
⎛ 1 + 0,1w ⎞
+ 0, 01752 ⎥
K ⎢ 0, 01873 ⎜
⎟
⎝ 1 − 0,1w ⎠
⎣
⎦
G ( w) =
2
⎛ 1 + 0,1w ⎞
⎛ 1 + 0,1w ⎞
⎜ 1 − 0,1w ⎟ − 1,8187 ⎜ 1 − 0,1w ⎟ + 0,8187
⎝
⎠
⎝
⎠
K (−0, 000333w2 − 0, 09633w + 0,9966)
=
w2 + 0,9969 w
(*)
33
w ⎞⎛
w⎞
⎛
K ⎜1 +
⎟ ⎜1 − ⎟
300 ⎠ ⎝ 10 ⎠
. Tentaremos um compensador “lead”; se ele não resolver o
≅ ⎝
w( w + 1)
problema, tentaremos um outro de ordem maior. Suponhamos, para simplificar, que o
controlador GD ( w) tenha ganho unitário para a faixa de freqüências baixas e tenha a
1+τ w
, 0 < α < 1.
seguinte forma: GD ( w) =
(**)
1 + ατ w
A função de transferência de malha aberta é então
1 + τ w K (−0, 000333w2 − 0, 09633w + 0,9966)
.
GD ( w)G ( w) =
1 + ατ w
w2 + 0,9969 w
Ora, quando z → 1 , w → 0 . Portanto, K v = lim wGD ( w)G ( w) , e da expressão acima,
w→0
temos K v ≅ K. Mas do enunciado do problema, temos K = 2. Substituindo em (*) e na
expressão subsequente, temos portanto
w ⎞⎛
w⎞
⎛
2 ⎜1 +
1− ⎟
⎟⎜
2
2(−0, 000333w − 0, 09633w + 0,9966)
300 ⎠⎝ 10 ⎠
G ( w) =
≅ ⎝
. A figura 4-37
2
w( w + 1)
w + 0,9969 w
abaixo apresenta o diagrama de Bode para o sistema.
Os ganhos, em decibéis, são marcados do lado esquerdo e as fases, em graus, do lado
direito, como de costume. Os ganhos e fases de G ( jν ) são apresentados em linhas
tracejadas.
Recorda-se que a margem de ganho é o acréscimo de ganho, em decibéis, que se pode
dar à malha aberta tal que o SMF se torne instável.
E por outro lado, a margem de fase é o ângulo que pode ser acrescentado à fase do
34
sistema em malha aberta tal que o SMF se torne instável.
A margem de ganho é obtida medindo a diferença de ganho, com relação a 0 dB, na
frequência em que a assíntota do diagrama de fase atinge - 1800 e é, como se pode ver
acima, igual a 14,5 dB.
A margem de fase é obtida medindo a diferença de fase com relação a - 1800 na
frequência em que a assíntota do diagrama de ganhos atinge 0 dB, sendo neste caso,
igual a 30 0 .
As duas frequências mencionadas acima são chamadas de frequências de “crossover”.
O compensador deve ser tal, portanto, que aumente a margem de fase em 20 0 . Mas este
aumento modifica a curva de ganho do diagrama de Bode, a respectiva frequência de
“crossover” (que são as duas frequências mencionadas acima) sendo deslocada para a
direita. Segue-se um processo de tentativa e erro. O que se segue, apresentado pelo livro
texto deste curso, é um método um tanto sofisticado, que não precisa ser conhecido
pelos alunos. O importante é que se chegue, por tentativa e erro, a um valor satisfatório.
O método apresentado consiste no seguinte: vamos dar um aumento de fase de 28 0 ,
tentativamente. Seja então φm = = 280 . Sabemos da teoria estudada no curso anterior
1−α
, o que dá α = 0, 361, onde este α é o da expressão (**). Da mesma
que senφm =
1+ α
expressão (**) temos as “frequências de esquina” (“corner frequencies”) ν = 1/ τ e
ν = 1/(ατ ) . Sabemos também que φm ocorre na média geométrica das duas freqüências
de esquina, ou seja, em ν = 1/( ατ ) . A modificação (aumento) na curva de ganho
devida à inclusão do termo (1 + τ jν ) /(1 + ατ jν ) é
1 + τ jν
1 + ατ jν
=
ν =1/( ατ )
1
α
.
A seguir, achamos a frequência onde o ganho do sistema sem compensador é igual a
⎛ 1 ⎞
⎛ 1 ⎞
-20log ⎜
⎟ = - 4,425dB.
⎟ = -20log ⎜
0,361
⎝ α⎠
⎝
⎠
Para achar a frequência onde o ganho é - 4,425dB, substitui-se w = jν na expressão de
G ( w) , obtendo G ( jν ) =
⎛ ν ⎞
2 1+ ⎜
⎟
⎝ 300 ⎠
2
⎛ν ⎞
1+ ⎜ ⎟
⎝ 10 ⎠
ν 1 +ν 2
2
. Por tentativa e erro, acha-se que
com ν = 1,7 , o ganho é aproximadamente – 4,4 dB. Escolhemos, portanto esta
frequência para a nova frequência de “crossover”, ν c .
1
1
, e, portanto, τ =
= 0,979, donde ατ = 0,3534 .
Ora, ν c = 1, 7 =
ατ
1, 7 α
1+τ w
Portanto o compensador “lead” é GD ( w) =
.
(4-40)
1 + ατ w
As curvas de ganho e de fase tanto de GD ( jν ) como as da função de transferência de
malha aberta GD ( jν )G ( jν ) aparecem no diagrama de Bode em linhas cheias. E do
diagrama se conclui que a margem de fase é 500 e a margem de ganho é 14 dB,
satisfazendo ao desejado no enunciado.
2 z −1
O controlador dado em (4-40) tem que ser re-escrito no plano z através de w =
=
T z +1
35
⎛ z −1 ⎞
1 + 0,979 ⎜ 10
⎟
z −1
⎝ z + 1 ⎠ = 2,3798 z − 1,9387 .
10
. Então, temos GD ( z ) =
z +1
z − 0,5589
⎛ z −1 ⎞
1 + 0,3534 ⎜ 10
⎟
⎝ z +1⎠
Donde que a função de transferência em malha aberta é:
2,3798 z − 1,9387 0, 03746( z + 0,9356)
GD ( z )G ( z ) =
z − 0,5589
( z − 1)( z − 0,8187)
0, 0891z 2 + 0, 0108 z − 0, 0679
= 3
. E a função de transferência do SMF é:
z − 2,3776 z 2 + 1,8352 z − 0, 4576
C ( z)
0, 0891z 2 + 0, 0108 z − 0, 0679
= 3
R ( z ) z − 2, 2885 z 2 + 1,846 z − 0,5255
0, 089( z + 0,9357)( z − 0,8145)
=
.
( z − 0,8126)( z − 0, 7379 − j 0,3196)( z − 0, 7379 + j 0,3196)
Observe-se que o zero em 0,8145 “quase” cancela o pólo em 0,8126 (um quase
cancelamento que não preocupa, porque o modo é estável) e assim os pólos complexos
conjugados atuam como pólos dominantes, o sistema comportando-se aproximadamente
como um sistema de segunda ordem.
Para obter a resposta no domínio do tempo, podemos usar o MATLAB. A resposta ao
degrau unitário pode ser obtida com os comandos do MATLAB Program 4-2 abaixo.
O resultado aparece na Figura 4-38 abaixo. Da figura, obtemos uma ultrapassagem
(“overshoot”) de aproximadamente 20% e um tempo de assentamento (“settling time”)
de cerca de 4 seg.
36
4-6 MÉTODO DE PROJETO ANALÍTICO
A razão principal da desvantagem dos controladores analógicos é que todos eles, sejam
pneumáticos, hidráulicos, etc., têm limitações físicas. Tais limitações podem ser quase
completamente ignoradas quando se projetam controladores digitais. Assim, muitos
esquemas de controle que eram impossíveis com controladores analógicos, tornaram-se
possíveis com controladores digitais.
Nesta seção estudaremos os controladores “deadbeat”, que fazem com que a resposta do
SMF apresente erro nulo em regime permanente, tempo mínimo de assentamento e
ausência de oscilações entre os instantes de amostragem depois que se atinge o regime
permanente.
Considere o sistema da Figura 4-39(a) abaixo.
Deseja-se projetar um controlador digital GD ( z ) de modo que o SMF responda com o
menor tempo de assentamento possível, com erro nulo em regime permanente, seja a
entrada um degrau, uma rampa ou uma parábola. Mais ainda, deseja-se que não existam
oscilações entre os instantes de amostragem depois que a resposta do sistema atingir o
regime permanente.
Definamos a transformada z da planta precedida do segurador de ordem zero:
⎡1 − e −Ts
⎤
G( z) = Z ⎢
G p ( s ) ⎥ . A função de transferência de malha aberta se torna
⎣ s
⎦
GD ( z )G ( z ) , conforme mostrado na figura 4-39(b).
Suponhamos que se queira que a função de transferência em malha fechada seja F ( z ) ,
GD ( z )G ( z )
C ( z)
dada. Então, temos: F ( z ) =
=
.
(4-41)
R( z ) 1 + GD ( z )G ( z )
a0 z N + a1 z N −1 + ... + aN
, ou
zN
F ( z ) = a0 + a1 z −1 + ... + aN z − N ,
Seja F ( z ) =
com N a determinar. De (4-41), temos: GD =
(4-42)
F ( z)
.
G ( z )[1 − F ( z )]
(4-43)
37
O controlador deve ser fisicamente realizável, o que impõe certas condições:
1. O grau do numerador deve ser menor ou igual ao grau do denominador, pois do
contrário, seria necessário conhecer as entradas futuras para produzir a resposta
presente.
2. Se a planta G p ( s ) tiver algum retardo e − Ls , então o SMF deve ter, pelo menos, o
mesmo retardo, pois do contrário, a saída do SMF teria que responder antes de uma
entrada ser dada, o que é impossível para um sistema fisicamente realizável.
3. Se G ( z ) for expandido em uma série em z −1 , o termo de potência mais baixa da
expansão de F ( z ) em z −1 deve ser pelo menos tão grande quanto a de G ( z ) : ver (441). Assim, por exemplo, se a expansão da série de G ( z ) começar com o termo em
z −1 , então o primeiro termo da expansão de F ( z ) em (4-42) deve ser zero, ou seja, a0 =
0, isto é, a expansão de F ( z ) será a1 z −1 + a2 z −2 + ... + aN z − N , onde N ≥ n , sendo n a
ordem da planta.
Além das condições de realizabilidade física, é preciso prestar atenção à estabilidade do
SMF. Como sabemos, estabilidade implica (condição necessária, mas não suficiente!)
que não haja cancelamento de pólos e zeros instáveis no produto das funções de
transferência da planta, segurador e controlador.
Investiguemos agora o que acontece com F ( z ) quando G ( z ) tiver um pólo instável ou
criticamente estável, isto é, um pólo fora ou sobre o círculo unitário. Seja α este pólo e
G ( z)
, sendo que G1 ( z ) não inclui termo que cancele z − α .
definamos G ( z ) = 1
z −α
Então a função de transferência do SMF é dada, tendo em vista (4-41), por
G ( z)
GD ( z ) 1
z −α .
F ( z) =
(4-44)
G1 ( z )
1 + GD ( z )
z −α
1
z −α
Além disso, 1 − F ( z ) =
=
, ou seja, α é um zero da
G1 ( z ) z − α + GD ( z )G1 ( z )
1 + GD ( z )
z −α
função de transferência do erro (diferença entre o sinal de entrada no SMF e a resposta
da planta), fato bem conhecido do curso anterior.
Observe-se ainda de (4-44) que se G ( z ) tiver um zero fora ou sobre o círculo unitário,
ele não pode ser cancelado (pois isto desestabilizaria o SMF), o que significa que todo
zero de G ( z ) que esteja fora ou sobre o círculo unitário, aparece em F ( z ) .
Procedamos com o projeto. Da figura 4-39(b), temos:
E ( z ) = R ( z ) − C ( z ) = R ( z )[1 − F ( z )] .
(4-45)
1
Para uma entrada em degrau unitário, temos R( z ) =
;
1 − z −1
Tz −1
para uma entrada em rampa unitária, temos R( z ) =
;
(1 − z −1 ) 2
e para uma entrada em parábola unitária, isto é, r (t ) = 0,5t 21(t ) , temos:
T 2 z −1 (1 + z −1 )
.
2(1 − z −1 )3
E assim, em geral, a transformada z de entradas de forma polinomial é:
R( z ) =
38
P( z )
,
(4-46)
(1 − z −1 ) q +1
onde P( z ) é um polinômio em z −1 .
Note-se que para um degrau unitário, temos P( z ) = 1 e q = 0; para uma rampa unitária,
temos P( z ) = T z −1 e q = 1 e para uma parábola unitária, temos
P( z ) = 0,5T 2 z −1 (1 + z −1 ) e q = 2.
P( z )[1 − F ( z )]
Substituindo esta em (4-45), temos: E ( z ) =
.
(4-47)
(1 − z −1 ) q +1
Para assegurar que a resposta do sistema atinja o regime permanente em número finito
de passos, mantendo então o erro igual a zero, E ( z ) deve ser um polinômio em z −1 com
um número finito de termos. Então, em vista de (4-47), escolhemos F ( z ) tal que
1 − F ( z ) = (1 − z −1 ) q +1 N ( z ) ,
(4-48)
−1
onde N ( z ) é um polinômio em z com um número finito de termos.
Consequentemente, E ( z ) = P ( z ) N ( z ) .
(4-49)
Finalmente, substituindo isto em (4-43), temos:
F ( z)
GD ( z ) =
.
(4-50)
G ( z )(1 − z −1 ) q +1 N ( z )
Com este compensador, obtém-se erro nulo nos instantes de amostragem após um
número finito de instantes.
Quanto aos tempos entre os instantes de amostragem, para garantir que o erro também
seja nulo, isto é, para que não haja “ondulações” entre os instantes de amostragem, é
preciso que
c(t ≥ nT ) = constante, para entrada em degrau;
c(t ≥ nT ) = constante, para entrada em rampa;
(4-50*)
c(t ≥ nT ) = constante para entrada em parábola.
Para garantir estas condições, observe-se que a planta é de tempo contínuo, portanto o
controle que aciona a planta, isto é, u (t ) , deve ser ou constante ou monotonicamente
crescente, ou decrescente, em regime permanente.
Alguns comentários:
1. Tendo em vista que F ( z ) é um polinômio em z −1 , fica claro que todos os pólos do
SMF estão em z = 0. Estes múltiplos pólos na origem são muito sensíveis à variação de
parâmetros. Efetivamente, o sistema de controles “dead-beat” não é “robusto”, ou seja,
a perturbação de parâmetros destrói a solução “dead-beat”, o erro tendendo a zero, mas
assintoticamente, e não em número finito de passos.
2. Observe-se que um sistema que apresente ótimo tempo de assentamento e erro nulo
para um certo tipo de entradas, pode revelar um desempenho muito ruim quando se
muda o tipo de entrada. É neste contexto que se torna importante o estudo do chamado
“controle ótimo”. Mas seja dito também que em número muito grande de processos
industriais, a solução do controle ótimo é dispensável, precisamente porque o tipo de
entrada é conhecido.
3. O período de amostragem, T, não entrou nos cálculos deste problema. Mas é claro
que quanto menor T, mais rapidamente o erro atingirá o valor nulo.
R( z ) =
39
Exemplo 4-13: Considere o problema da mesma figura 4-39(a) acima, onde
1
. Achar um compensador GD ( z ) tal que o SMF exiba uma resposta
G p ( s) =
s ( s + 1)
“dead-beat” para um degrau de entrada. O período de amostragem é T = 1 seg. A
seguir, usando este mesmo compensador, calcule a resposta do sistema a uma rampa de
entrada.
Solução: Seguindo os passos do estudo teórico acima, calculemos
⎡1 − e −Ts
⎡ 1 ⎤
1 ⎤
−1
G( z) = Z ⎢
⎥ = (1 − z )Z ⎢ 2
⎥
⎣ s ( s + 1) ⎦
⎣ s s ( s + 1) ⎦
⎡ z −1
⎤ 0,3679(1 + 0, 7181z −1 ) z −1
1
1
=
.
(4.51)
= (1 − z ) ⎢
−
+
−1 2
−1
1 − 0,3679 z −1 ⎥⎦
(1 − z −1 )(1 − 0,3679 z −1 )
⎣ (1 − z ) 1 − z
Passamos ao diagrama de blocos da figura (b), definindo a função de transferência do
C ( z)
GD ( z )G ( z )
=
. Agora observe-se de (4-51) que se G ( z ) for
SMF: F ( z ) =
R( z ) 1 + GD ( z )G ( z )
−1
expandido em uma série em z −1 , o primeiro termo da série será 0,3679 z −1 . Portanto
F ( z ) deve começar com um termo em z −1 .
Reportemo-nos à eq. (4-42), lembrando que N ≥ n . Como a planta é de 2ª. ordem,
tentemos:
F ( z ) = a1 z −1 + a2 z −2 .
(4-52)
Como a entrada é um degrau, temos da eq. (4-48):
1 − F ( z ) = (1 − z −1 ) N ( z ) .
(4-53)
Como vimos, tendo G ( z ) um pólo em z = 1, a estabilidade do SMF exige que 1 − F ( z )
tenha um zero em z = 1. Mas de (4-53) vemos que esta condição já é satisfeita.
Agora, observe-se que a saída de um segurador de ordem zero é constante entre os
instantes de amostragem; em consequência, a 1ª. condição de (4-50*) é satisfeita. Donde
que: U ( z ) = b0 + b1 z −1 .
(4-53*)
Ora, da figura (b) é claro que: U ( z ) =
C ( z)
G( z)
=
E tendo em vista (4-51), temos: U ( z ) = F ( z )
C ( z ) R( z )
R( z ) G ( z )
= F ( z)
R( z )
G( z)
.
1
(1 − z −1 )(1 − 0,3679 z −1 )
1 − z −1 0,3679(1 + 0, 7181z −1 ) z −1
(1 − 0,3679 z −1 )
.
(4-53**)
0,3679(1 + 0, 7181z −1 ) z −1
Ora, para que U ( z ) seja uma série com apenas dois termos, conforme (4-53*), é
(4-54)
preciso que F ( z ) = (1 + 0, 7181z −1 ) z −1 F1 ,
= F ( z)
onde F1 é uma constante. Então, de (4-53**),
U ( z ) = 2,1781(1 − 0,3679 z −1 ) F1 .
(4-55)
Agora vamos determinar N ( z ) , F ( z ) e F1 . Substituindo (4-52) em (4-53), temos
1 − a1 z −1 − a2 z −2 = (1 − z −1 ) N ( z ) .
Ora, se dividirmos esta eq. por 1 − z −1 , o quociente do lado esquerdo é
1 + (1 − a1 ) z −1 = N ( z )
(4-56)
e o resto é (1 − a1 − a2 ) z −2 . Ora, o resto deve ser nulo, o que implica
1 − a1 − a2 = 0.
(4-57)
40
Mas das eqs. (4-52) e (4-54), temos
F ( z ) = a1 z −1 + a2 z −2 = (1 + 0, 7181z −1 ) z −1 F1 .
Portanto, a1 + a2 z −1 = (1 + 0, 7181z −1 ) F1 .
Dividindo-se esta eq. por 1 + 0, 7181z −1 , obtém-se no lado esquerdo o quociente a1 e o
resto (a2 − 0, 7181a1 ) z −1 .
Igualando o quociente com F1 e zerando o resto, temos
F1 = a1 e a2 − 0, 7181a1 = 0.
Desta última e de (4-57), temos a1 = 0,582 e a2 = 0,418. Portanto,
F ( z ) = 0,582 z −1 + 0, 418 z −2
e F1 = 0,582 ,
(4-58)
(4-59)
sendo que a eq. (4-56) nos dá N ( z ) = 1 + 0, 418 z −1 .
(4-60)
A função de transferência de pulso do controlador digital é determinada a partir de
(4-50), com (4-51), (4-54) e (4-60):
(1 + 0, 7181z −1 ) z −1 (0,582)
F ( z)
GD ( z ) =
=
G ( z )(1 − z −1 ) N ( z ) 0,3679(1 + 0, 7181z −1 )
(1 − z −1 )(1 + 0, 418 z −1 )
(1 − z −1 )(1 − 0,3679 z −1 )
1,582 − 0,582 z −1
.
1 + 0, 418 z −1
Com este controlador digital, obtém-se a seguinte função de transferência de pulso do
C ( z)
0,582( z + 0, 7181)
SMF:
= F ( z ) = 0,582 z −1 + 0, 418 z −2 =
.
R( z )
z2
E a resposta a um degrau unitário é obtida por
1
C ( z ) = F ( z ) R( z ) = (0,582 z −1 + 0, 418 z −2 )
= 0,582 z −1 + z −2 + z −3 + z −4 + .... . E,
−1
1− z
portanto, c(0) = 0; c(1) = 0,582; c(k ) = 1 para k = 2;3; 4;...
=
Por outro lado, de (4-55) com F1 = 0,582 , temos U ( z ) = 2,1781× 0,582(1 − 0,3679 z −1 )
= 1,582 − 0,582z −1 . Portanto, o controle u (k ) torna-se nulo para k ≥ 2 , como desejável
e, consequentemente, não há ondulações entre os intervalos de amostragem na resposta,
depois de atingido o regime permanente.
A figura 4-40(a) abaixo mostra c(k ) e u (k ) versus k e c(t ) versus t no caso de a entrada
ser um degrau unitário.
41
Vamos agora investigar, como requerido no enunciado, o comportamento da resposta do
sistema a uma entrada em rampa unitária:
z −1
C ( z ) = F ( z ) R ( z ) = (0,582 z −1 + 0, 418 z −2 )
(1 − z −1 ) 2
= 0,582 z −2 + 1,582 z −3 + 2,582 z −4 + 3,582 z −5 + ...
Por sua vez, o controle é dado, em vista das eqs. (4-51) e (4-59) por:
C ( z ) F ( z ) R( z ) F ( z ) z −1
0,582 z −1 + 0, 418 z −2
z −1
=
=
=
U ( z) =
G( z )
G( z )
G( z ) (1 − z −1 ) 2 0,3679(1 + 0, 7181z −1 ) z −1 (1 − z −1 ) 2
(1 − z −1 )(1 − 0,3679 z −1 )
z −1
= (1,582 − 0,582 z )
= 1,582 z −1 + z −2 + z −3 + z −4 + ...
−1 2
(1 − z )
Ou seja, o sinal u (k ) torna-se constante a partir de k = 2, o que implica que a resposta
em regime permanente não terá oscilações, depois de atingido o regime permanente.
A figura 4-40(b) acima mostra c( k ) e u ( k ) versus k e c(t ) versus t para o caso de a
entrada ser uma rampa unitária. Observe-se da primeira figura que há um erro na
1
resposta, mas que se mantém constante. E este erro é, como sabemos, ess =
, onde
Kv
−1
⎡1 − z −1
⎤
⎡
⎤
F ( z)
=
K v = lim ⎢
GD ( z )G( z ) ⎥ = lim ⎢(1 − z −1 )
−
1
z →1
(1 − z ) N ( z ) ⎥⎦
⎣ T
⎦ z →1 ⎣
0,582 z −1 + 0, 418 z −2
lim
= 0, 7052 , donde ess = 1, 418 , como indicado na 1ª. figura de 4z →1
1 + 0, 418 z −1
40(b).
Exemplo 4-14: Considere o mesmo problema que o anterior, agora a constante de
velocidade K v é especificada, igual a 4 seg −1 . (Por conta disto, como veremos, o tempo
de assentamento será maior). Mas, além disso, desejamos que o tempo de assentamento
seja o mínimo possível que satisfaça a estas especificações. O período de amostragem é
T = 1. Projetar um controlador digital que satisfaça a estas especificações.
42
Solução: A transformada z da planta precedida pelo segurador de ordem zero já foi
0,3679(1 + 0, 7181z −1 ) z −1
obtida no exemplo anterior em (4-51): G( z ) =
. E tal como no
(1 − z −1 )(1 − 0,3679 z −1 )
GD ( z )G ( z )
C ( z)
=
. E como o primeiro termo da
exemplo anterior, temos F ( z ) =
R( z ) 1 + GD ( z )G ( z )
expansão de GD ( z ) é 0,3679z −1 , F ( z ) deve também começar com um termo em z −1 :
F ( z ) = a1 z −1 + a2 z −2 + ... + aN z − N , onde N ≥ n = 2 . Em vista das restrições adicionais
deste problema, vamos tentar N = 3. (No exercício anterior resolvemos o problema com
N = 2). Então, seja F ( z ) = a1 z −1 + a2 z −2 + a3 z −3 .
(4-61)
(Se não obtivermos resultado satisfatório, tentaremos N > 3). Da equação (4-48), temos,
em vista do fato que a entrada é um degrau: 1 − F ( z ) = (1 − z −1 ) N ( z ) .
(4-62)
Por outro lado, do final do problema anterior e do enunciado deste, temos
⎡1 − z −1
⎤
⎡
⎤ F (1)
F ( z)
=
K v = lim ⎢
=4.
GD ( z )G( z ) ⎥ = lim ⎢(1 − z −1 )
−1
→
z →1
1
z
(1 − z ) N ( z ) ⎥⎦ N (1)
⎣
⎣ T
⎦
Mas de (4-62) temos F (1) =1.
1
Disto e da expressão acima, temos K v =
= 4.
(4-63)
N (1)
Tendo em vista que se requer que a resposta do sistema não apresente oscilações entre
os intervalos de amostragem depois de atingido o regime permanente, devemos ter
U ( z ) = b0 + b1 z −1 + b2 z −2 + b( z −3 + z −4 + z −5 + ...) . Mas a planta tem um integrador,
portanto b = 0, ou seja, U ( z ) = b0 + b1 z −1 + b2 z −2 .
C ( z ) C ( z ) R( z )
R( z )
=
= F ( z)
Mas da fig. 4-39 (b), temos U ( z ) =
G( z ) R( z ) G( z )
G( z )
1 − 0,3679 z −1
= F ( z)
.
0,3679(1 + 0, 7181z −1 ) z −1
Para que U ( z ) seja uma série em z −1 com três termos, F ( z ) deve ter a seguinte forma:
F ( z ) = (1 + 0, 7181z −1 ) z −1 F1 ( z ) ,
(4-64)
onde F1 ( z ) deve ser um polinômio do 1º. grau em z −1 . E, portanto de U ( z ) acima, temos
U ( z ) = 2, 7181(1 − 0,3679 z −1 ) F1 ( z ) .
(4-65)
Ora, de (4-61) e (4-62), temos 1 − F ( z ) = 1 − a1 z −1 − a2 z −2 − a3 z −3 = (1 − z −1 ) N ( z ) .
Se dividirmos esta eq. por 1 − z −1 , o quociente de 1 − a1 z −1 − a2 z −2 − a3 z −3 é
1 + (1 − a1 ) z −1 + (1 − a1 − a2 ) z −2 e o resto é (1 − a1 − a2 − a3 ) z −3 . Então,
N ( z ) = 1 + (1 − a1 ) z −1 + (1 − a1 − a2 ) z −2
e o resto deve ser nulo, ou seja, 1 − a1 − a2 − a3 = 0.
(4-66)
(4-67)
observando-se que de (4-63) devemos ter N (1) = 1/ 4 . Portanto, substituindo z −1 =1 em
(4-68)
(4-66), temos 2a1 + a2 = 2,75.
Além disso, a eq. (4-64) pode ser re-escrita como
F ( z ) = a1 z −1 + a2 z −2 + a3 z −3 = (1 + 0, 7181z −1 ) z −1 F1 ( z ) , ou seja,
(4-68*)
a1 + a2 z −1 + a3 z −2 = (1 + 0, 7181z −1 ) F1 ( z ) .
A divisão desta última eq. por 1 + 0, 7181z −1 dá o quociente a1 + (a2 − 0, 7181a1 ) z −1 e o
43
resto igual a [a3 − 0, 7181(a2 − 0, 7181a1 )]z −2 . Ora, então,
F1 ( z ) = a1 + (a2 − 0, 7181a1 ) z −1 e
(4-68**)
a3 − 0, 7181(a2 − 0, 7181a1 ) = 0.
(4-69)
As eqs. (4-67), (4-68) e (4-69) nos dão: a1 = 1, 26184 , a2 = 0, 22633 e a3 = −0, 48816 .
E, portanto, de (4-68*), F ( z ) = 1, 26184 z −1 + 0, 22633z −2 − 0, 48816 z −3 ,
enquanto que de (4-68**), vem F1 ( z ) = 1, 26184 − 0, 67979 z −1 .
E por outro lado, a eq. (4-66) nos dá N ( z ) = 1 − 0, 26184 z −1 − 0, 48817 z −2 .
Com estes resultados, obtemos de (4-50):
(1 − 0,5387 z −1 )(1 − 0,3679 z −1 )
GD ( z ) = 3, 4298
(1 − 08418 z −1 )(1 + 0,5799 z −1 )
e a resposta é: C ( z ) = 1, 2618 z −1 + 1, 4882 z −2 + z −3 + z −4 + ... , ou seja,
c(1) = 1, 2628 , c(2) = 1, 4882 , c(k ) = 1 para k = 3; 4;5;... ,
donde se vê que a resposta ao degrau unitário tem uma ultrapassagem máxima de
aproximadamente 50% e o tempo de assentamento é 3 seg.
Observe-se que da eq. (4-65) nós temos:
U ( z ) = 3, 4298 − 3,1096 z −1 + 0, 6798 z −2 ,
ou seja, o controle u (k ) torna-se nulo para k ≥ 3 , conforme a hipótese.
Consequentemente, não há ondulações entre amostragens, depois de alcançado o regime
permanente.
A Figura 4-41 abaixo mostra c(k ) versus k , u (k ) versus k e u (t ) versus t para uma
entrada em degrau unitário. Observe-se que a tentativa N = 3 na eq. (4-61) revelou-se
satisfatória.
44
Vamos agora investigar a resposta do sistema a uma rampa unitária. Obtemos, depois de
muitos cálculos: C ( z ) = 1, 2618 z −2 + 2, 75 z −3 + 3, 75 z −4 + ... ,
enquanto que: U ( z ) = 3, 4298 z −1 − 3, 202 z −2 + z −3 + z −4 + z −5 + ... ,
ou seja, o sinal de controle torna-se constante para k > 2, o que implica que a resposta
do sistema não exibirá oscilações entre amostragens depois que o sistema atinge o
regime permanente.
A Figura 4-42 abaixo mostra a resposta do sistema e o controle, sendo que este é dado
tanto nos instantes de amostragem, como em tempo contínuo. Seja notado ainda que o
erro em regime permanente é ess = 1/ K v = 1/ 4 , como indicado na 1ª. figura.
Comparando este exemplo, 4-14, com o anterior, 4-13, verificamos que o 2º. melhora a
resposta à rampa unitária a expensas do tempo de assentamento, isto é, precisa-se de
uma amostragem a mais para se atingir o regime permanente. Mas quanto à resposta ao
degrau, o primeiro é melhor, pois tem um tempo de assentamento menor e não tem
ultrapassagem.
45
Problema com soluções:
Problema A-4-1: Mostre que a configuração dos pólos perto de z = 1 no plano z é
similar à dos pólos perto de s = 0 no plano s.
Solução: z = eTs . Perto da origem no plano s temos z = eTs = 1 + Ts + 0,5T 2 s 2 + .... ,
donde z − 1 ≈ Ts , provando o resultado.
Problema A-4-2: Considere o sistema descrito por:
y (k ) − 0, 6 y (k − 1) − 0,81y (k − 2) + 0, 67 y (k − 3) − 0,12 y (k − 4) = x(k ) , onde y é a
resposta e x é a entrada. Estudar a estabilidade do sistema pelo método de Jury.
Solução: A função de transferência de pulso do sistema é:
z4
Y ( z)
1
=
=
.
X ( z ) 1 − 0, 6 z −1 − 0,81z −2 + 0, 67 z −3 − 0,12 z −4 z 4 − 0, 6 z 3 − 0,81z 2 + 0, 67 z − 0,12
Portanto o polinômio característico é: P ( z ) = z 4 − 0, 6 z 3 − 0,81z 2 + 0, 67 z1 − 0,12 .
Aplicando o critério de Jury para estabilidade, temos:
1. −0,12 < 1 . OK.
2. P (1) = 1 − 0, 6 − 0,81 + 0, 67 − 0,12 = 0,14 > 0 . OK.
3. P (−1) = 1 + 0, 6 − 0,81 − 0, 67 − 0,12 = 0 , a condição não sendo satisfeita, pois deveria
ser > 0. Mas vemos que o sistema tem um pólo em z = -1.
−0,12
1
−0,12 0, 67
4. b3 =
= −0,9856 e b2 =
= −0,598 . Donde b3 > b2 . OK.
1
−0,12
1
−0, 6
5. c2 =
b3
b0
b0
b3
=
−0,9856
−0,598
−0,598
−0,9856
= 0, 6138 ,
46
b2 −0,9856 −0,5196
=
= −0,5834 . Ora c2 > c0 . OK.
b1
−0,598 −0,9072
Portanto concluímos que todas as condições de Jury, à exceção da 3ª, são satisfeitas.
Pode-se verificar que o pólo em -1 tem grau de multiplicidade 1 e portanto o sistema é
criticamente estável.
c0 =
b3
b0
Problema A-4-3: Considere a seguinte eq. característica
P ( z ) = z 3 − 1,3z 2 − 0, 08 z + 0, 24 = 0 . Determine se alguma(s) raiz(es) da eq. está fora do
círculo unitário. Usar a transformação bilinear e o critério de Routh.
⎛ w +1⎞
⎛ w +1⎞
⎛ w +1⎞
Solução: ⎜
⎟ − 1,3 ⎜
⎟ − 0, 08 ⎜
⎟ + 0, 24 = 0 , ou seja, após as operações
⎝ w −1 ⎠
⎝ w −1 ⎠
⎝ w −1 ⎠
óbvias, temos: −0,14w3 + 1, 06w2 + 5,1w + 1,98 = 0 . Como os coeficientes do polinômio
não têm o mesmo sinal, já podemos conclui que este polinômio não é Hurwitz, isto é,
não tem todas as raízes na parte boa do plano complexo, isto é, não tem todas as raízes
no semi-plano aberto da esquerda. Mas o critério de Routh permite determinar quantas
são as raízes no semi-plano da direita. Dividindo toda a eq. acima por -0,14 , temos:
w3 − 7,571w2 − 36, 43w − 14,14 = 0 ,
(4-71)
3
cujo tablóide de Routh é:
2
w3
w2
w
w0
1
−36, 43
−7,571 −14,14
.
0
−38,3
0
−14,14
Recorda-se que o termo relativo a w é
−7, 571 × ( −36, 43) + 14,14
−7, 571
= −38, 29 ≅ −38, 3 , ao
−38,3 × (−14,14)
= −14,14
−38,3
Como se vê, há uma troca de sinal na segunda coluna, donde se conclui que há uma raiz
no semi-plano da direita do plano w, e, portanto uma raiz fora do círculo unitário no
plano z.
passo que o termo correspondente a w0 é
Y ( z)
0, 787 z −1
Problema A-4-4: Considere o sistema dado por
=
=
U ( z ) 1 − 0,8195 z −1 + 0, 6065 z −2
0, 787 z
. O tempo de amostragem é T = 0,5 seg. Usando o MATLAB,
2
z − 0,8195 z + 0, 6065
plotar a resposta do sistema a uma rampa unitária até k = 20.
Solução: ver o MATLAB Program 4-3 abaixo, seguido da resposta na Figura 4-43.
47
Problema A-4-5: Demonstrar que se a eq. característica de um SMF for
KB( z )
1+
= 0 , onde B( z ) e A( z ) são coprimos e não contêm K como fator, então os
A( z )
pontos de “breakaway” e “break-in” podem ser determinados pelas raízes de
dK
A '( z ) B( z ) − A( z ) B '( z )
=−
= 0 , onde ' indica derivada com relação a z.
dz
B2 ( z)
(4-72)
Prova: Re-escrevendo a eq. característica: f ( z ) = A( z ) + KB( z ) = 0 .
Suponha que a eq. acima tenha uma raiz múltipla de ordem r. Seja z1 esta raiz. Então,
temos: f ( z ) = ( z − z1 ) r ( z − z2 )...( z − z p ) . E é claro que
df ( z )
= 0 . Isto significa que
dz z = z1
48
as raízes múltiplas de f ( z ) satisfazem à eq.
df ( z )
= 0 , ou seja,
dz
df ( z )
= A '( z ) + KB '( z ) = 0 ,
(4-73)
dz
A '( z )
. Este valor de K conduz a múltiplas raízes da eq. característica.
donde K = −
B '( z )
A '( z )
Substituindo este valor de K em (4-72), temos f ( z ) = A( z ) −
B ( z ) = 0 , ou seja,
B '( z )
A '( z ) B( z ) − A( z ) B '( z ) = 0.
(4-74)
Quando esta eq. é resolvida, obtêm-se os valores de z onde ocorrem múltiplas raízes.
A( z )
, e, portanto,
Por outro lado, da eq. (4-72), nós temos K = −
B( z )
dK
A '( z ) B( z ) − A( z ) B '( z )
=−
. Mas quando igualamos isto a zero, obtemos a eq.
dz
B2 ( z)
(4-74). Portanto, os pontos breakaway e break-in podem são determinados a partir de
dK
= 0 , concluindo aprova. Mas deve ser notado que nem todas as soluções desta eq.
dz
são pontos breakaway ou break-in. Será se o valor correspondente de K for real.
Problema A-4-7: Obtenha o lugar das raízes (root locus) no plano z do diagrama de
blocos na Figura 4-45 abaixo para os seguintes períodos de amostragem: T = 1; 2 e 4
seg.
Solução: Reportando-nos ao exemplo 3-5, em que não havia K no numerador da planta,
K [(T − 1 − e −T ) z −1 + (1 − e −T − Te −T ) z −2
temos G ( z ) =
.
(4-75)
(1 − z −1 )(1 − e −T z −1 )
1. Para T = 1, da eq. de cima, temos G ( z ) =
K [(1 − 1 − e −1 ) z −1 + (1 − e −1 − e −1 ) z −2
(1 − z −1 )(1 − e −1 z −1 )
0,3679 K ( z + 0, 7181)
.
(4-75*)
( z − 1)( z − 0,3679)
Para a construção do root locus, utilizam-se as regras, (re)vistas depois de (4-20). Em
primeiro lugar, plotam-se os pólos e os zeros da função de transferência, conforme a
=
49
Figura 4-46(a) abaixo, observando-se que a função tem um zero no infinito.
Consequentemente, o root locus tem dois ramos. Para K = 0, os pólos do SMF
coincidem com os do sistema em malha aberta, que são os da função de transferência
(4-75*). À medida que K aumenta, os pólos caminham um para a direita e o outro para
a esquerda no eixo real, de acordo com a 3ª. regra, isto é, cada ponto sobre o eixo real,
para pertencer a um ramo do root locus, deve ter à sua direita um número impar de
pólos mais zeros, até se encontrarem no ponto breakaway, calculado de acordo com a
5ª. regra, deixando então o eixo real, numa curva que é simétrica com relação ao eixo
real (pois os pólos ocorrem aos pares conjugados), até se encontrarem de novo no eixo
real, no ponto break-in, também calculado de acordo com a 5ª. regra, caminhando um
ramo para a direita até encontrar o zero finito e o outro para a esquerda, para o infinito,
observando-se de novo a 3ª. regra, mencionada acima. Obtemos 0,6479 para o ponto de
breakaway e – 2,0841 para o break-in. De acordo com a 6ª. regra, calculam-se os
ângulos de partida (do breakaway) e de chegada (no break-in) no eixo real. Obtém-se
900 em ambos os casos. Observe-se ainda que, de acordo com a 8ª. regra, todo ponto
50
( z − 1)( z − 0,3679)
, a qual é obtida a
0,3679( z + 0, 7181)
partir da eq. característica 1 + KG ( z ) = 0 , que também pode ser usada na construção do
root locus. A partir disto pode-se verificar, ponto por ponto, que os pontos do root locus
fora do eixo real formam um círculo. Seja z = x + jy . O círculo unitário tem por eq.
x 2 + y 2 = 1 , enquanto o root locus, com centro em -0,7181 e raio igual a 0,7181 +
0,3679 = 1,086 tem por eq. ( x + 0, 7181) 2 + y 2 = 1,086 2 . Destas duas eqs. achamos os
valores de x e y em que os círculos se interceptam e, substituindo no valor de K acima,
encontramos K = 2,3925 , que é o ganho crítico de estabilidade conforme indicado na
figura (a). Ou seja, para K > 2,3925 , o SMF é instável.
1,1353K ( z + 0,5232)
2. Para T = 2, temos de (4-75): G ( z ) =
. O root locus é mostrado
( z − 1)( z − 0,1353)
na fig. (b) acima.
3, 0183K ( z + 0,301)
3. Finalmente para T = 4, temos da eq. (4-75): G ( z ) =
. O root
( z − 1)( z − 0, 0183)
locus está na figura (c). Observe-se que, neste caso, o SMF é estável para todo K.
do root locus satisfaz à eq. (4-24), , isto é, K =
Problema A-4-8: Seja o sistema de controles da Figura 4-47 abaixo. O período de
amostragem é T = 1 seg. Projetar um controlador PI tal que os pólos dominantes do
SMF tenham um fator de amortecimento ς igual a 0,5 e o número de amostras por ciclo
de oscilações senoidais amortecidas igual a 10. Obter a resposta do sistema a um degrau
unitário de entrada.
Solução:
⎡1 − e −Ts e −2 s ⎤
1 ⎤
(1 − e −1 ) z −1
−1
−2 ⎡
−1
−2
=
−
=
−
G( z) = Z ⎢
(1
z
)
z
Z
(1
z
)
z
⎥
⎢ s ( s + 1) ⎥
(1 − z −1 )(1 − e −1 z −1 )
⎣
⎦
⎣ s s + 1⎦
0, 6321z −3
0, 6321
=
= 2
. O controlador digital PI tem, como sabemos, a
−1
1 − 0,3679 z
z ( z − 0,3679)
K z − K p + KI z
1
seguinte função de transferência GD ( z ) = K p + K I
= p
−1
z −1
1− z
Kp
z−
(K + KI )z − K p
K p + KI
= p
= (K p + KI )
.
(4-75**)
z −1
z −1
51
E, portanto, a função de transferência em malha aberta é
⎛
Kp ⎞
(K p + KI ) ⎜ z −
⎟
⎜
K p + K I ⎟⎠
0, 6321
⎝
.
GD ( z )G ( z ) =
2
z −1
z ( z − 0,3679)
(4-75***)
A seguir plotamos no plano complexo os pólos do sistema de malha aberta, conforme
(4-75***), na Figura 4-48 abaixo, o zero finito sendo desconhecido por enquanto, pois
depende de K p e K I .
Tendo em vista que queremos 10 amostras por ciclo de oscilações senoidais
amortecidas, o pólo dominante em malha fechada no semi-plano superior deve ficar
num ângulo de 360 (= 3600 ÷ 10 ) = 2π /10 , conforme indicado na figura.
Re-escrevemos as eqs. (4-1) e (4-2):
⎛ −2πς ω ⎞
ω
d
⎟ e arg[ z ] = 2π d .
(4-1bis) e (4-2bis)
z = exp ⎜
2
⎜ 1 − ς ωs ⎟
ω
s
⎝
⎠
Então desta última, temos, para o pólo dominante: 2π
ωd
ω
= 2π /10 , ou seja, d = 0,1.
ωs
ωs
⎛ −2π (0,5) 1 ⎞
Visto que ς = 0,5, temos da eq. (4-1), repetida acima: z = exp ⎜
⎟
⎜ 1 − 0,52 10 ⎟
⎝
⎠
−0,3628
= 0,6958 , o que no permite plotar o pólo dominante P na figura, ou seja,
=e
z = 0, 6958arg(360 ) = 0,5629 + j 0, 409 .
Ora, de acordo com a 6ª. regra, que se aplica a qualquer ponto do root locus, e não
somente aos pontos complexos de partida (breakaway) e de chegada (break-in), a soma
dos ângulos do ponto P aos pólos do sistema em malha aberta menos a soma dos
ângulos aos zeros finitos deve somar 1800 . Então, tendo em vista que o pólo na origem
é duplo, o ângulo ao zero finito do sistema em malha aberta é:
360 + 360 + 136,90 + 64, 620 − 1800 = 93,520 , conforme indicado na figura.
E isto determina a posição do zero em 0,5881 sobre o eixo real. Por conseguinte, de (4-
52
75**),
Kp
K p + KI
= 0,5881.
(4-76)
z − 0,5881
, onde K = K p + K I . O ganho K é
z −1
z − 0,5881
0, 6321
= 1,
determinado pela condição de magnitude: K
2
z − 1 z ( z − 0,3679) z =0,5629+ j 0,409
Então o controlador PI é: GD ( z ) = K
o que nos dá K = 0,507, ou seja, K p + K I = 0,507 .
(4-77)
Desta eq. e de (4-76), chega-se a K p = 0, 2982 e K I = 0, 2088 . E, portanto o
controlador PI é GD ( z ) = 0,507
1 − 0,5881z −1
e a função de transferência da malha
1 − z −1
aberta é:
⎛ 1 − 0,5881z −1 0, 6321z −3 ⎞
⎛ 1 − 0,5881z −1
⎞
z −3
=
GD ( z )G ( z ) = 0,507 ⎜
0,3205
⎜
−1
−1 ⎟
−1
−1 ⎟
1 − 0,3679 z ⎠
1 − 0,3679 z ⎠
⎝ 1− z
⎝ 1− z
, a função de transferência do SMF sendo
C ( z)
0,3205 z −3 − 0,1885 z −4
=
.
R( z ) 1 − 1,3679 z −1 + 0,3679 z −2 + 0,3205 z −3 − 0,1885 z −4
A resposta c(kT ) pode ser obtida facilmente conforme o MATLAB Program 4-4, que
apresenta os comandos usados, o resultado aparecendo na Figura 4-49, logo a seguir.
53
Problema A-4-9: Considere o sistema da Figura 4-50 abaixo. Deseja-se projetar um
controlador digital tal que os pólos dominantes do SMF tenham um fator de
amortecimento ς = 0,5 . Queremos que o número de amostras por ciclo da oscilação
senoidal amortecida seja igual a 8. O período de amostragem é T = 0,2 . Usando o
método do lugar das raízes, determinar a função de transferência de pulso do
controlador digital. Obter a resposta do sistema projetado a uma entrada em degrau
unitário bem como a constante de erro de velocidade, K v .
Solução: Repetimos as eqs. (4-1) e (4-2), que já foram repetidas no problema anterior:
⎛ −2πς ω ⎞
d
⎟,
(4-78)
z = e −T ωn = exp ⎜
⎜ 1 − ς 2 ωs ⎟
⎝
⎠
arg[ z ] = T ωd = 2π
ωd
=θ .
ωs
(4-78*)
Tendo em vista que queremos 8 amostras por ciclo na oscilação senoidal amortecida, o
pólo dominante do sistema em malha fechada deve estar localizado sobre uma linha que
passe pela origem e forme ângulo de 450 (= 3600 / 8) com o eixo dos reais. Então, temos
ω
ω 1
π
arg[ z ] = 450 = = 2π d , o que dá d = .
(4-79)
ωs
ωs 8
4
54
Tendo em vista que T = 0,2 , temos ωs =
2π 2π
=
= 10π , e portanto,
T
0, 2
1
10π
= 3,927 . Substituindo (4-79) em (4-78) e usando o valor dado do
8
8
amortecimento, obtemos z = e −0,4535 = 0, 6354 . Portanto o ponto P, que é o pólo
ωd = ω s =
dominante na parte superior do plano tem magnitude 0,6534 e ângulo de 450 , ou seja, o
ponto P é dado por z = 0, 4493 + j 0, 4493 , tal como indicado na Figura 4-51 abaixo.
A seguir, calculamos a função de transferência de pulso da planta precedida de um
⎡1 − e −Ts
⎡ 1 ⎤
1 ⎤
amostrador de ordem zero: G ( z ) = Z ⎢
= (1 − z −1 )Z ⎢ 2
⎥
⎥
⎣ s ( s + 1) ⎦
⎣ s s ( s + 1) ⎦
⎡ 0, 2 z −1
⎤
0, 01873(1 + 0,9356 z −1 ) z −1
1
1
= (1 − z ) ⎢
−
+
⎥ = (1 − z −1 )(1 − 0,8187 z −1 )
−1 2
−1
1 − e −0,2 z −1 ⎦
⎣ (1 − z ) 1 − z
0, 01873( z + 0,9356)
=
.
( z − 1)( z − 0,8187)
−1
(4-79*)
Plotamos os pólos e zero da figura 4-51, conforme (4-79*). Tendo em vista que o ponto
P é o pólo dominante desejado do SMF, o ângulo com que o controlador deve
contribuir é 140, 790 + 129, 430 − 17,970 − 1800 = 72, 250 , onde o primeiro ângulo do lado
esquerdo da igualdade é o formado a partir do pólo z = 1, o segundo é o formado a partir
do pólo z = 0,8187 e o terceiro a partir do zero z = -0,9356 , a soma das contribuições
dos ângulos da planta e do compensador devendo ser igual a 1800 .
s+a
Seja o controlador da forma: GD ( z ) = K
.
s+b
Escolhemos o zero do controlador de modo a cancelar o pólo z = 0,8187. A partir disto e
da condição do ângulo vista acima, calculamos o pólo do compensador: z = 0,1595. E
55
assim, temos o compensador procurado, a menos de K ainda a ser determinado:
1 − 0,8187 z −1
.
(4-79**)
GD ( z ) = K
1 − 0,1595 z −1
Consequentemente, a função de transferência de pulso de malha aberta é:
1 − 0,8187 z −1 0, 01873(1 + 0,9356 z −1 ) z −1
GD ( z )G ( z ) = K
.
1 − 0,1595 z −1 (1 − z −1 )(1 − 0,8187 z −1 )
K pode ser determinado a partir da condição de magnitude, a saber, tendo havido
0, 01873( z + 0,9356)
cancelamento de pólo e zero: K
= 1 , ou seja,
( z − 1)( z − 0,1595) z =0,4493+ j 0,4493
K = 13,934 . De (4-79**), temos o compensador procurado:
1 − 0,8187 z −1
, a função de transferência de pulso do SMF sendo
GD ( z ) = 13,934
1 − 0,1595 z −1
GD ( z )G ( z )
0, 261z −1 + 0, 2442 z −2
C ( z)
=
=
facilmente calculável:
.
R( z ) 1 + GD ( z )G ( z ) 1 − 0,8985 z −1 + 0, 4073z −2
A Figura 4-52 abaixo mostra c(kT ) , a resposta do sistema no domínio do tempo a um
degrau unitário. Finalmente, a constante de velocidade de erro K v é determinada por
⎡1 − z −1
⎤
⎡1 − z −1 0, 261(1 + 0, 9356 z −1 ) z −1 ⎤
K v = lim ⎢
GD ( z )G ( z ) ⎥ = lim ⎢
⎥ = 3, 005 .
−1
−1
z →1
⎣ T
⎦ z →1 ⎣ 0, 2 (1 − 0,1595 z )(1 − z ) ⎦
Problema A-4-10: Considere o sistema da Figura 4-53 abaixo. Suponha que as
especificações de desempenho são dadas em termos de margem de fase, margem de
ganho, erro estático de velocidade, etc. Estabelecer os procedimentos para projetar
compensadores “lead” e “lag” pelo método frequencial.
56
Solução:
1+τ w
, com
1 + ατ w
0 < α < 1 . Consequentemente, a função de transferência em malha aberta é:
1+τ w
1+τ w
GD ( w)G ( w) = K D
G ( w) =
G1 ( w) , onde G1 ( w) = K D G ( w) . Os passos
1 + ατ w
1 + ατ w
seguintes não estão necessariamente na ordem que se deve seguir sempre. Esta depende
dos dados do problema. Os passos seguintes nem sempre precisam ser usados, de novo
dependendo dos dados do problema. Mas é bom tê-los presentes sempre que se faz um
projeto.
1. Determinar o ganho K D para satisfazer a especificação a respeito da velocidade de
erro, observando que z → 1 se só se w → 0 .
2. Plotar o diagrama de Bode de G1 ( w) = K D G ( w) com o ganho K D determinado
acima. Medir a margem de fase.
3. Determinar o avanço necessário de fase φ a ser adicionado ao sistema.
1. Compensador “lead”. Como sabemos, ele tem a forma GD ( w) = K D
4. Acrescentar 50 a 120 a φ , tendo em vista o deslocamento da frequência de
“crossover”. Definir este novo avanço de fase como φm . Determinar o efeito de
1−α
.
atenuação α a partir de senφm =
1+ α
5. Determinar a freqência onde a magnitude do sistema descompensado G1 ( jν ) é igual
(
)
a 20 log 1/ α . Escolher esta frequência como a nova frequência de “crossover”. Esta
(
)
frequência corresponde a ν m = 1/ ατ e o valor máximo do deslocamento (“shift”) da
fase φm corresponde a esta frequência.
6. Determinar as frequências de canto (“corner frequencies”) do compensador “lead” do
seguinte modo: zero do compensador “lead” igual a
1
τ
; pólo do compensador “lead”
1
.
ατ
7. Conferir a margem de ganho para verificar se é satisfatória. Se não o for, repetir o
processo, modificando os pólos e zeros do compensador até que um resultado
satisfatório seja obtido.
1+τ w
, onde
2. Compensador “lag”. O compensador tem agora a forma GD ( w) = K D
1 + βτ w
β > 1 . A função de transferência do sistema em malha aberta é
1+τ w
1+τ w
GD ( w)G ( w) = K D
G ( jw) =
G1 ( jw) , onde G1 ( w) = K D G ( w) . Valem aqui
1 + βτ w
1 + βτ w
igual a
57
também as observações acima (compensador “lead”) sobre ordem e importância dos
passos.
1. Determinar o ganho K D para satisfazer a especificação a respeito da velocidade de
erro.
2. Se o sistema descompensado G1 ( w) não satisfizer às especificações de ganhos de
margem e de fase, ache a frequência para a qual o ângulo de fase da função de
transferência de malha aberta seja igual a - 1800 mais a margem de fase desejada. A
esta deve-se somar 50 a 120 para compensar pela redução de fase do compensador
“lag”. Escolha esta frequência como a nova frequência de “crossover”.
3. Para evitar prejuízos devidos à margem de fase em vista do compensador “lag”, o
pólo e zero do compensador devem ser colocados bem abaixo da nova frequência de
“crossover”. Portanto, escolha-se a frequência de “corner” ν = 1/ τ (correspondendo ao
zero do compensador “lag”) uma década abaixo da nova frequência de “crossover”.
4. Determinar a atenuação necessária para trazer a curva de magnitude a 0 dB no novo
ganho da frequência de “crossover”. Notando-se que esta atenuação é igual a -20 log β ,
determine o valor de β . E então a outra frequência de “corner”, correspondendo ao
pólo do compensador “lag” é determinado por ν = 1/( βτ ) .
5. Uma vez que o compensador é projetado no plano w, ele deve ser transformado no
compensador no plano z, GD ( z ) . Note-se que no plano z, o pólo e o zero do
compensador estão próximos um do outro, e na realidade, próximos de z = 1. Note-se
ainda que, tendo em vista que os coeficientes do filtro (compensador) devem ser
realizados por palavras digitais, se o número de bits usado for insuficiente, a localização
do zero e do pólo do filtro podem não ficar onde se deseja e o compensador terá um
desempenho diferente do esperado.
Problema A-4-11: Projetar um controlador digital para o sistema da Figura 4-54 abaixo.
Usar o método do diagrama de Bode no plano w. As especificações do projeto são as
seguintes: margem de fase de 550 , margem de ganho de pelo menos 10 dB e erro de
velocidade de 5/seg., período de amostragem T = 0,1. Depois de projetar o
compensador, construa o root locus, localize os pólos do SMF no root locus e determine
o número de amostragens por ciclo de oscilações senoidais amortecidas.
Solução: A transformada z da planta precedida pelo segurador de ordem zero é
⎡1 − e −Ts
⎡
⎤
1
1 ⎤
= (1 − z −1 )Z ⎢ 2
G( z) = Z ⎢
⎥
⎥=
⎣ s ( s + 2) ⎦
⎣ s s ( s + 2) ⎦
1 + 0,9355 z −1
z + 0,9355
0, 004683 z
= 0, 004683
. Passamos para o
−1
−1
(1 − z )(1 − 0,8187 z )
( z − 1)( z − 0,8187)
−1
58
1 + (Tw / 2) 1 + 0, 05w
, obtendo, após
=
1 − (Tw / 2) 1 − 0, 05w
0,5(1 + 0, 001666 w)(1 − 0, 05w)
simplificações, G ( w) =
.
w(1 + 0,501w)
O diagrama de Bode de G(w) está na Figura 4-55 abaixo com w = jν .
w
1+
1+τ w
a ,
= KD
Escolhemos um controlador da forma GD ( w) = K D
w
1 + ατ w
1+
b
onde a = 1/ τ e b = 1/(ατ ) . Então a função de transferência de malha aberta é
1 + ( w / a) 0,5(1 + 0, 001666 w)(1 − 0, 05w)
GD ( w)G ( w) = K D
. É dado do problema que
w(1 + 0,501w)
1 + ( w / b)
K v = 5/seg. Então, K v = lim[ wGD ( w)G ( w)] = 0,5 K D , donde K D = 10.
plano w através da transformação z =
w→0
Usando agora técnica convencional de projeto (passos 2 a 4 do compensador “lag”),
achamos que a função de transferência do controlador é dada por
w ⎞
⎛
⎜ 1 + 1,994 ⎟
(4-79**)
GD ( w) = 10 ⎜
⎟.
⎜ 1+ w ⎟
⎜
12,5 ⎟⎠
⎝
Portanto, a função de transferência do sistema em malha aberta é:
w ⎞
⎛
⎜ 1 + 1,994 ⎟ 0,5(1 + 0, 001666 w)(1 − 0, 05w)
.
GD ( w)G ( w) = 10 ⎜
⎟
w(1 + 0,501w)
⎜ 1+ w ⎟
⎜
12,5 ⎟⎠
⎝
Traçando o diagrama de Bode (abaixo), verificamos que esta função de transferência
nos dá uma margem de fase de aproximadamente 550 e uma margem de ganho de
aproximadamente 12,4 dB. Portanto, todas as especificações foram obtidas.
59
Só nos resta obter GD ( z ) através da transformação bilinear w =
2 z −1
. Substituindo
T z +1
este valor em (4-79**) acima, obtemos, após alguns cálculos:
⎛ 1 − 08187 z −1 ⎞
. E a função de transferência em malha aberta é:
GD ( z ) = 42, 423 ⎜
−1 ⎟
⎝ 1 − 0, 2308 z ⎠
0,1987( z + 0,9355)
.
( z − 0, 2308)( z − 1)
A Figura 4-56 abaixo mostra o root locus do sistema.
GD ( z )G ( z ) =
Foram traçados no root locus os lugares dos pólos quando ς = 0,5 e quando ς = 0, 6 ,
(fazendo variar a freqência ν ?). Dá para estimar que o fator de amortecimento ς é
aproximadamente igual a 0,55. A linha que liga a origem ao pólo no plano superior
pode ser medida no root locus, obtendo-se um ângulo de 370 . Portanto o número de
amostragem por ciclo de oscilações senoidais amortecidas é 3600 / 370 = 9,73.
Problema A-4-12: Seja o sistema de controle digital da Figura 4-57 abaixo. Projetar um
controlador digital no plano w de modo que a margem de fase seja de 500 e a margem
de ganho seja de pelo menos 10 dB. O período de amostragem é T = 0,1. Após projetar
o controlador, obter a constante de erro de velocidade, K v . Obter também a resposta do
sistema a uma entrada em degrau unitário.
2
−1
−1
⎡1 − e −Ts 1 ⎤
⎡1⎤
−1
−1 T (1 + z ) z
= (1 − z )Z ⎢ 3 ⎥ = (1 − z )
, esta última
Solução: G ( z ) = Z ⎢
2⎥
2(1 − z −1 )3
⎣s ⎦
⎣ s s ⎦
60
0, 005( z + 1)
.
( z − 1) 2
1 + 0, 05w
A seguir, utilizando a transformação bilinear, temos z =
, donde que
1 − 0, 05w
1 − 0, 05w
1 − 0, 05 jν
, ou seja, G ( jν ) =
. A partir disto, podemos obter o
G ( w) =
2
( jν ) 2
w
diagrama de Bode para G ( jν ) , que se encontra na Figura 4-58 abaixo.
expressão tendo sido obtida da Tabela 2-1. Donde G ( z ) =
Observe-se que a margem de fase do sistema em malha aberta pode ser medida igual a
aproximadamente - 20 . Então é necessário adicionar um controlador “lead” de modo a
obter as margens de fase e ganho requeridas.
Utilizando métodos convencionais de projeto (ver Problema A-4-10), pode-se verificar
que o seguinte compensador lead satisfaz às especificações:
1+ w
w +1
= 64
GD ( w) = 4
. Com a adição deste compensador, a frequência do
w + 16
1 + w /16
ganho de crossover passa a ser ν = 4. Observe-se que o maior avanço de fase que este
1−α
compensador pode produzir é dado por arcsenφm , sendo senφm =
, onde α é a
1+ α
relação dos termos independentes de z no numerador e denominador da função de
transferência do compensador e neste caso é:
1 − 1/16
= arcsen(0,8824) = 61,930 . A margem de fase obtida do diagrama
φm = arcsen
1 + 1/16
de Bode é 50,620 e a margem de ganho, como se pode ver, é de aproximadamente
igual a 13 dB. Portanto, as especificações são satisfeitas.
z − 0,9048
Agora, através da transformação bilinear, obtemos: GD ( z ) = 37,333
. (4-79)
z − 0,1111
E a função de transferência de malha aberta é:
0,1867(1 − 0,9048 z −1 )(1 + z −1 ) z −1
GD ( z )G ( z ) =
.
(1 − 0,1111z −1 )(1 − z −1 ) 2
A constante de erro de velocidade é obtida através da fórmula:
61
⎡1 − z −1
K v = lim ⎢
⎤
GD ( z )G ( z ) ⎥ ; efetuando, obtém-se K v = ∞ , ou seja, o erro em regime
⎣ T
⎦
permanente quando a entrada é uma rampa é nulo, o que já seria de se esperar, tendo em
vista que a planta tem um duplo integrador.
A função de transferência de pulso do SMF é:
C ( z ) 0,1867 z −1 + 0, 0178 z −2 − 0,1689 z −3
=
. A Figura 4-59 abaixo mostra a resposta do
1 − 1,9244 z −1 + 1, 24 z −2 − 0, 28 z −3
R( z )
sistema a uma entrada em degrau unitário.
z →1
Observe que de (4-79), o controlador tem um zero em 0,9048, portanto próximo ao
duplo pólo em z = 1. Quando isto acontece, isto é, um zero próximo a um pólo, a
resposta apresenta uma longa “cauda”, como se vê da fig. 4-59.
Problema A-4-13: Seja o SMF da Figura 4-60 abaixo.
Observe-se que a planta tem um retardo de 5 seg. Supõe-se que o tempo de amostragem
T seja igual ao de retardo da planta. A resposta desejada do sistema a um degrau
unitário de entrada é mostrada abaixo na Figura 4-61(a).
62
A resposta vai de 0 (no quinto segundo) ao valor final em 10 segundos, não havendo
ultrapassagem nem erro em regime permanente. O tempo de assentamento é de 15
segundos (medido de t = 0 a t = 15 segundos). Deseja-se que não haja ondulações
depois que o tempo de assentamento é alcançado. Projetar um controlador digital
GD ( z ) .
Solução: A transformada z da planta precedida do segurador de ordem zero é:
⎡1- eTs e −5 s ⎤
⎤
1
0,3935 z −2
−1
−1 ⎡
(1
)
Z
−
=
z
z
=
G( z) = Z ⎢
⎥
⎢ s (10 s + 1) ⎥ 1 − 0, 6065 z −1 . O sistema é
⎣
⎦
⎣ s 10 s + 1 ⎦
estável, com os pólos simples em z = 0 e z = 0,6065. A função de transferência de pulso
C ( z)
GD ( z )G ( z )
=
= F ( z) .
é:
(4-80)
R( z ) 1 + GD ( z )G ( z )
Da figura (a), vemos que h(1 − e −0,1(t −5) ) ; em t = 15 temos h(1 − e −1 ) = 1,582. Da figura
(a), que mostra que o controle é “dead-beat”, temos nos instantes discretos: c(0) = 0 ,
c(1) = 0 , c(2) = h(1 − e−0,5 ) = 0,6225 , c(k ) = 1 para k = 3;4;5;... Disto, temos:
1
C ( z ) = 0, 6225 z −2 + z −3 + z −4 + z −5 + ... = 0, 6225 z −2 + z −3
=
1 − z −1
0, 6225 z −2 + 0,3775 z −3
.
(4-80*)
1 − z −1
1
e em visa de (4-80*), temos:
Mas de (4-80) temos C ( z ) = F ( z ) R( z ) = F ( z )
1 − z −1
F ( z ) = 0, 6225 z −2 + 0,3775 z −3 .
63
F ( z)
.
(4-80**)
G ( z )[1 − F ( z )]
Mas de (4-48) lembra-se que 1 − F ( z ) = (1 − z −1 ) N ( z ) , onde N ( z ) é um polinômio em
z −1 , ou seja, 1 − 0, 6225 z −2 − 0,3775 z −3 = (1 − z −1 ) N ( z ) , obtendo-se
N ( z ) = 1 + z −1 + 0,3775 z −2 , donde que 1 − F ( z ) = (1 − z −1 )(1 + z −1 + 0,3775 z −2 ) .
Utilizando (4-80**), obtemos, após alguns cálculos triviais:
1,582(1 − 0,3678 z −2 )
.
GD ( z ) =
(1 − z −1 )(1 + z −1 + 0,3775 z −2 )
C ( z)
Vejamos o controle que aciona a planta U ( z ) =
. Substituindo os resultados
G( z)
Ora, GD ( z ) =
1 − 0,3678 z −2
anteriores, obtemos, após alguns cálculos: U ( z ) = 1,582
1 − z −1
= 1,582 + 1,582 z −1 + z −2 + z −3 + z −4 + ... Donde se vê que o controle u (k ) é constante
para k > 1, o que faz que não haja oscilacões entre as amostragens depois de atingido o
regime permanente. u (k ) versus k é mostrado na figura 4-61(b) acima.
Problema A-4-14: Seja o sistema na figura 4-62 abaixo. Projetar um controlador digital
GD ( z ) tal que o SMF tenha um tempo de assentamento mínimo com erro em regime
permanente nulo para uma entrada em rampa unitária. Além disso, não deve haver
oscilacões entre os instantes de amostragem, depois de atingido o regime permanente.
Dado T = 1 seg. Depois de projetado o controlador, verificar a resposta do sistema a
uma entrada em delta de Kronecker e a um degrau unitário.
−1
−1
⎡1 − e −Ts 1 ⎤
⎡ 1 ⎤ (1 + z ) z
−1
(1
)Z
−
Solução: G ( z ) = Z ⎢
=
z
=
2⎥
⎢⎣ s 3 ⎥⎦ 2(1 − z −1 ) 2 ,
⎣ s s ⎦
esta última igualdade sendo obtida usando a Tabela 2-1. Definamos a função de
C ( z)
GD ( z )G ( z )
=
= F ( z) .
transferência de pulso:
R( z ) 1 + GD ( z )G ( z )
Observe-se que se G ( z ) for expandido como uma série em z −1 , o primeiro termo da
série será 0,5 z −1 . Consequentemente, F ( z ) é uma série que começa com z −1 :
F ( z ) = a1 z −1 + a2 z −2 + ....aN z − N , ou seja, queremos um controle “dead-beat”, onde
N ≥ n , n sendo a ordem do sistema. No caso presente, n = 2.
Tendo em vista que a entrada é uma rampa, temos da eq. (4-48):
1 − F ( z ) = (1 − z −1 ) 2 N ( z ) .
(4-81)
Observe-se que G ( z ) tem um pólo duplo crítico em z = 1. Em virtude da estabilidade
do SMF, então 1 − F ( z ) deve ter um duplo zero em z = 1 para cancelar este pólo da
64
planta e prover para o controle dead beat. Mas acontece que, de (4-81), esta condição já
é satisfeita.
Tendo em vista que não devem existir oscilacões na resposta em regime permanente,
então U ( z ) deve ser uma série em z −1 com a seguinte característica:
U ( z ) = b0 + b1 z −1 + b2 z −2 + ....bN −1 z − N +1 + b( z − N + z − N −1 + z − N − 2 + ....)
Mas como N ≥ n e n = 2, e tendo a planta um integrador duplo, então b acima deve ser
nulo, pois do contrário, a resposta cresceria parabolicamente, e não linearmente.
Consquentemente, temos: U ( z ) = b0 + b1 z −1 + b2 z −2 + ....bN −1 z − N +1 .
C ( z ) C ( z ) R( z )
R( z )
.
Ora, do diagrama de blocos, temos U ( z ) =
=
= F ( z)
G ( z ) R( z ) G ( z )
G( z)
Mas da tabela 2-1 temos R( z ) =
Tz −1
, com T = 1. Donde,
(1 − z −1 ) 2
z −1
2(1 − z −1 ) 2
2
= F ( z)
.
−1 2
−1
−1
(1 − z ) (1 + z ) z
1 + z −1
Mas para que U ( z ) seja uma série em z −1 com um número finito de termos, F ( z ) tem
que ser divisível por 1+ z −1 , ou seja,
F ( z ) = (1 + z −1 ) F1 ( z ) ,
(4-82)
U ( z) = F ( z)
onde F1 ( z ) é um polinômio em z −1 com um número finito de termos. Donde,
U ( z ) = 2 F1 ( z ) .
(4-83)
De (4-81) e (4-82) vemos que F ( z ) deve ter termos pelo menos de z −3 . Então,
(4-83*)
suporemos que F ( z ) = a1 z −1 + a2 z −2 + a3 z −3 .
Esta expressão de F ( z ) envolve o número mínimo de termos, o transiente terminará em
3 amostragens. Para determinar os coeficientes de F ( z ) acima, temos da eq. (4-81):
1 − a1 z −1 − a2 z −2 − a3 z −3 = (1 − z −1 ) 2 N ( z ) .
Dividindo o lado esquerdo da eq. acima por (1 − z −1 ) 2 , o resto deve ser nulo. Obtemos:
N ( z ) = 1 + (2 − a1 ) z −1
(4-83**)
e o resto nulo, a saber,
[2(2 − a1 ) − (1 + a2 )]z −2 − (2 − a1 + a3 ) z −3 = 0 para todo z, o que implica:
2(2 − a1 ) − (1 + a2 ) = 0
e
(4-84)
2 − a1 + a3 = 0.
(4-85)
Por outro lado, de (4-82), nós temos a1 z −1 + a2 z −2 + a3 z −3 = (1 + z −1 ) F1 ( z ) .
Dividindo o lado esquerdo desta eq. por 1 + z −1 , o resto deve ser nulo. Obtemos então
F1 ( z ) = a1 z −1 + (a2 − a1 ) z −2 ,
(4-85*)
a1 − a2 + a3 = 0 .
(4-86)
Resolvendo as eqs. (4-84), (4-85) e (4-86), temos:
a1 = 1, 25 ; a2 = 0,5 e a3 = −0, 75 .
(4-86*)
Portanto, obtemos de (4-83**): N ( z ) = 1 + 0, 75 z −1
e de (4-85*): F1 ( z ) = 1, 25 z −1 (1 − 0, 6 z −1 ) .
−1
−2
(4-86**)
−3
−1
−1
De (4-83*) e de (4-86*): F ( z ) = 1, 25 z + 0,5 z − 0, 75 z = 1, 25 z (1 + z )(1 − 0, 6 z −1 ) .
65
De (4-50), através de cálculo rotineiro: GD ( z ) =
2,5(1 − 0, 6 z −1 )
.
(1 + 0, 75 z −1 )
Então, por cálculos simples:
C ( z ) = F ( z ) R( z ) = 1, 25 z −2 + 3z −3 + 4 z −4 + 5 z −5 + ... , ou seja,
c(0) = 0 ; c(1) = 0 ; c(2) = 1, 25 ; c(k ) = k para k = 3;4;5;...
Note-se que de (4-83) e de (4-86**), temos U ( z ) = 2,5 z −1 − 1,5 z −2 .
Ou seja, o controle torna-se nulo para k > 2. Consequentemente, não há oscilacões entre
as amostragens no regime permanente. A Figura 4-63 abaixo mostra a resposta e o
controle plotados contra o tempo, sendo que quanto ao controle, ele é mostrado também
depois de passar pelo segurador de ordem zero.
Investiguemos agora as respostas deste sistema a uma entrada que seja um delta de
Kronecker e a um degrau unitário.
Ora, quando a entrada é um delta de Kronecker, temos:
C ( z ) = F ( z ) = 1, 25 z −1 + 0,5 z −2 − 0, 75 z −3 .
F ( z)
E o controle é, efetuando cálculos simples: C ( z ) =
= 2,5 − 6,5 z −1 + 5,5 z −2 − 1,5 z −3 .
G( z)
Note-se que o controle se torna nulo para k > 3. Portanto, não há oscilacões entre as
amostragens a partir de t = 4T = 4.
Quando a entrada é um degrau, temos:
1
C ( z ) = F ( z ) R( z ) = (1, 25 z −1 + 0,5 z −2 − 0, 75 z −3 )
= 1, 25 z −1 + 1, 75 z −2 + z −3 + z −4 + ...
−1
1− z
Observa-se que a ultrapassagem máxima é de 75%. Calculamos também
R( z )
U ( z) = F ( z)
= 2,5 − 4 z −1 + 1,5 z −2 , donde vemos que o controle se anula a partir de
G( z)
k = 3, instante a partir do qual não há oscilacões entre amostragens.
66
A Figura 4-64(a) abaixo mostra a resposta e o controle quando a entrada é um delta de
Kronecker, enquanto que (b) mostra resposta e controle quando a entrada é um degrau
unitário.
67
CAPÍTULO 5: ANÁLISE NO ESPACO DE ESTADO
5-1 Introdução
Nos capítulos anteriores tratamos de sistemas invariantes no tempo com uma entrada e
uma saída. Para tais sistemas, os métodos usados nos dois capítulos anteriores são
bastante úteis.
Os métodos usados não se aplicam, entretanto, a sistemas não lineares, exceto quando
estes são muito simples (aos quais se aplicam as “funções descreventes”) e não se
aplicam ao problema do controle ótimo e do controle adaptativo que são, geralmente,
não lineares.
Ora, os modernos sistemas de controle podem ter muitas entradas e respostas, as quais
podem ser inter-relacionadas de um modo complicado. Os métodos do espaço de
estado, que estudaremos neste capítulo, são apropriados, mas não exclusivos, para
estudar tais sistemas. E têm também a vantagem de se aplicarem, com as devidas
modificações, a sistemas não lineares.
O método do espaço de estado
Ele é baseado na descrição do sistema em n eqs. de diferença (ou diferencial, se se trata
de sistemas de tempo contínuo) de 1ª. ordem, as quais podem ser agrupadas, dando
lugar a uma eq. vetorial-matricial de 1ª. ordem. O uso da notação vetorial-matricial
simplifica bastante a representação matemática do sistema de equações.
Estado
O estado de um sistema é o menor conjunto de variáveis (chamadas variáveis de estado)
tais que, o conhecimento destas variáveis em t = t0 , juntamente com o conhecimento da
entrada para t ≥ t0 , determina completamente o comportamento do sistema para t ≥ t0 .
Variáveis de Estado
As variáveis de estado de um sistema dinâmico são as variáveis constituindo o menor
número de variáveis que determinam o estado do sistema, tal como referido acima. Se
pelo menos n variáveis x1 ; x2 ; ....; xn são necessárias para descrever completamente o
comportamento do sistema, então estas variáveis são um conjunto de variáveis de
estado. É importante observar que as variáveis de estado não precisam ser quantidades
fisicamente mensuráveis ou observáveis. Na prática, porém, é muitas vezes conveniente
escolher quantidades facilmente mensuráveis, se isto for possível, quando se faz
necessária a realimentação de variáveis mensuráveis.
Vetor de estado
Se n variáveis de estado são necessárias para descrever completamente o
comportamento de um sistema, então estas n variáveis podem ser consideradas como as
n componentes de um vetor x . Tal vetor é chamado vetor de estado. (Observe-se que
doravante os vetores e matrizes são representados em negrito).
Espaço de estado
O espaço de dimensão n cujos eixos coordenados são constituídos por x1 ; x2 ; ....; xn é
chamado um espaço de estado. Qualquer estado do sistema pode ser representado por
um ponto no espaço de estado.
Equações no espaço de estado
Na análise do espaço de estado, nós trabalharemos com três tipos de variáveis que
68
integram a modelagem de um sistema dinâmico: variáveis de entrada (“inputs”),
variáveis de saída ou resposta (“outputs”) e variáveis de estado.
Para sistemas variantes no tempo (lineares ou não-lineares), de tempo discreto, as eqs.
de estado podem ser escritas da seguinte forma:
x(k + 1) = f [x(k ), u(k ), k ] ,
enquanto que a eq. da resposta é dada por
y (k ) = g[x(k ), u(k ), k ] .
Para sistemas lineares variantes no tempo, temos
x(k + 1) = G (k )x(k ) + H (k )u(k )
y (k ) = C(k )x(k ) + D(k )u(k ) .
Suporemos sempre que o estado x(k ) tem dimensão n, a resposta y (k ) tem dimensão m
e a entrada u(k ) tem dimensão r e as matrizes têm as dimensões apropriadas. O fato de
as matrizes G , H , C e D dependerem explicitamente de k indicam que o sistema é
variante no tempo. Caso k não aparecesse explicitamente como argumento nestas
matrizes, o sistema seria invariante no tempo.
Ou seja, as eqs. de sistemas lineares invariantes no tempo, que nos ocuparão neste
capítulo, são dadas por:
x(k + 1) = Gx(k ) + Hu(k ) ,
(5-1)
y (k ) = Cx(k ) + Du(k ) .
(5-2)
A Figura 5-1(a) acima apresenta o diagrama de blocos correspondentes às eqs. (5-1) e
(5-2). A figura (b) se refere a sistemas de tempo contínuo, que deveria ter sido estudado
no curso anterior. Note-se que as configurações são basicamente iguais, o retardo z −1 I
no tempo discreto sendo substituído pela integral no caso do tempo contínuo.
69
5-2 Representação no espaço de estado de sistemas de tempo discreto
Formas canônicas para equações de tempo discreto no espaço de estado
Há muitas técnicas para obter representação em espaço de estado de sistemas de tempo
discreto.
Seja o sistema de tempo discreto descrito por
y (k ) + a1 y (k − 1) + a2 y (k − 2) + ... + an y ( k − n) = b0u ( k ) + b1u (k − 1) + ... + bnu (k − n) . (5-5)
Esta eq. pode ser escrita, equivalentemente, na forma de uma função de transferência de
Y ( z ) b0 + b1 z −1 + ... + bn z − n
=
(5-6)
pulso, a saber,
U ( z ) 1 + a1 z −1 + ... + an z − n
Y ( z ) b0 z n + b1 z n −1 + ... + bn
= n
.
(5-7)
U ( z)
z + a1 z n −1 + ... + an
Há muitas maneiras de realizar a representação em espaço de estado para o sistema de
tempo discreto descrito pelas eqs. (5-5), (5-6) ou (5-7). Serão apresentadas aqui as
seguintes representações:
- Forma canônica controlável
- Forma canônica observável
- Forma canônica diagonal
- Forma canônica de Jordan
ou:
Forma canônica controlável
A representação de espaço de estado do sistema de tempo discreto dado pelas eqs.
acima pode ser apresentada pelas seguintes eqs. vetoriais:
1
0
.... 0 ⎤ ⎡ x1 (k ) ⎤ ⎡ 0 ⎤
⎡ x1 (k + 1) ⎤ ⎡ 0
⎢
⎥ ⎢
0
1
.... 0 ⎥⎥ ⎢⎢ x2 (k ) ⎥⎥ ⎢ 0 ⎥
⎢ x2 (k + 1) ⎥ ⎢ 0
⎢ ⎥
⎢ ..... ⎥ = ⎢ ....
(5-8)
....
.... .... .... ⎥ ⎢ .... ⎥ + ⎢....⎥ u (k )
⎢
⎥ ⎢
⎥⎢
⎥ ⎢ ⎥
0
0
.... 1 ⎥ ⎢ xn −1 (k ) ⎥ ⎢ 0 ⎥
⎢ xn −1 (k + 1) ⎥ ⎢ 0
⎢⎣ xn (k + 1) ⎥⎦ ⎢⎣ − an − an −1 − an − 2 .... −a1 ⎥⎦ ⎢⎣ xn (k ) ⎥⎦ ⎢⎣ 1 ⎥⎦
⎡ x1 (k ) ⎤
⎢ x (k ) ⎥
⎢ 2
⎥
(5-9)
y (k ) = [bn − anb0 bn −1 − an −1b0 .... b2 − a2b0 b1 − a1b0 ] ⎢ .... ⎥ + b0u (k ) .
⎢
⎥
⎢ xn −1 (k ) ⎥
⎢⎣ xn (k ) ⎥⎦
A prova deste resultado básico deveria ter sido vista no curso anterior. Se você está
duvidando, faça um exemplo simples, de sistema de ordem 2, para se convencer.
Agora, observe-se que se revertermos a ordem das variáveis, isto é, se definirmos:
⎡ xˆ1 (k ) ⎤ ⎡ 0 0 .... 0 1 ⎤ ⎡ x1 (k ) ⎤
⎢ xˆ (k ) ⎥ ⎢ 0 0 .... 1 0 ⎥ ⎢ x (k ) ⎥
⎢ 2 ⎥=⎢
⎥⎢ 2 ⎥,
(5-9*)
⎢ .... ⎥ ⎢... ... ... ... ...⎥ ⎢ .... ⎥
⎢
⎥ ⎢
⎥
⎥⎢
⎣ xˆn (k ) ⎦ ⎣ 1 0 ... 0 0 ⎦ ⎣ xn (k ) ⎦
então as eqs. acima podem ser escritas da seguinte forma:
70
⎡ xˆ1 (k + 1) ⎤ ⎡ −a1
⎢ˆ
⎥ ⎢
⎢ x2 (k + 1) ⎥ ⎢ 1
⎢ xˆ3 (k + 1) ⎥ = ⎢ 0
⎢
⎥ ⎢
⎢ .... ⎥ ⎢ ...
⎢⎣ xˆn (k + 1) ⎥⎦ ⎢⎣ 0
y (k ) = [b1 − a1b0
− a2 ... − an −1
0
...
0
1
...
0
...
...
...
0
...
1
b2 − a2b0
− an ⎤ ⎡ xˆ1 (k ) ⎤ ⎡ 1 ⎤
0 ⎥⎥ ⎢⎢ xˆ2 (k ) ⎥⎥ ⎢⎢ 0 ⎥⎥
0 ⎥ ⎢ xˆ3 (k ) ⎥ + ⎢ 0 ⎥ u (k )
⎥ ⎢ ⎥
⎥⎢
... ⎥ ⎢ ... ⎥ ⎢...⎥
0 ⎥⎦ ⎢⎣ xˆn (k ) ⎥⎦ ⎢⎣ 0 ⎥⎦
⎡ xˆ1 (k ) ⎤
⎢ xˆ (k ) ⎥
... bn − an b0 ] ⎢ 2 ⎥ + b0u (k ) .
⎢ .... ⎥
⎢
⎥
⎣ xˆn (k ) ⎦
(5-10)
(5-11)
Forma canônica observável
O sistema dado pelas eqs. (5-5), (5-6) ou (5-7) pode ser posto na seguinte forma:
⎡ x1 (k + 1) ⎤ ⎡ 0 0 ... 0 −an ⎤ ⎡ x1 (k ) ⎤ ⎡ bn − anb0 ⎤
⎢ x (k + 1) ⎥ ⎢ 1 0 ... 0 −a ⎥ ⎢ x (k ) ⎥ ⎢b − a b ⎥
n −1 ⎥ ⎢ 2
⎢ 2
⎥ ⎢
⎥ ⎢ n −1 n −1 0 ⎥
⎢ ..... ⎥ = ⎢.... .... ... .... .... ⎥ ⎢ .... ⎥ + ⎢
⎥ u (k )
(5-12)
....
⎢
⎥ ⎢
⎥⎢
⎥ ⎢
⎥
⎢ xn −1 (k + 1) ⎥ ⎢ 0 0 ... 0 − a2 ⎥ ⎢ xn −1 (k ) ⎥ ⎢ b2 − a2b0 ⎥
⎢⎣ xn (k + 1) ⎥⎦ ⎢⎣ 0 0 ... 1 − a1 ⎥⎦ ⎢⎣ xn (k ) ⎥⎦ ⎢⎣ b1 − a1b0 ⎥⎦
⎡ x1 (k ) ⎤
⎢ x (k ) ⎥
⎢ 2
⎥
y (k ) = [ 0 0 .... 0 1] ⎢ .... ⎥ + b0u (k ) ,
⎢
⎥
⎢ xn −1 (k ) ⎥
⎢⎣ xn (k ) ⎥⎦
a prova deste resultado também devendo ter sido vista no curso anterior.
(5-13)
Observe-se que a matriz da eq. de estado (5-12) é a transposta da matriz de (5-8). Mais
ainda, os elementos do vetor que multiplica a resposta (entrada) em (5-8)- (5-9) são os
mesmos que multiplica a e entrada (resposta) em (5-12) e (5-13). Dizemos então que o
par ((5-8), (5-9)) é dual do par ((5-12), (5-13)), ou que os pares são duais um do outro.
Note-se também que se a ordem das variáveis de estado for invertida tal como em
(5-9*), então as eqs. (5-12) e (5-13) são transformadas em:
1 ... 0 0 ⎤ ⎡ xˆ1 (k ) ⎤ ⎡ bn − anb0 ⎤
⎡ xˆ1 (k + 1) ⎤ ⎡ −a1
⎢ xˆ (k + 1) ⎥ ⎢ − a
0 ... 0 0 ⎥⎥ ⎢⎢ xˆ2 (k ) ⎥⎥ ⎢⎢bn −1 − an −1b0 ⎥⎥
⎢ 2
⎥ ⎢ 2
⎥ u (k )
⎢ ..... ⎥ = ⎢ .... .... ... .... ....⎥ ⎢ .... ⎥ + ⎢
....
(5-14)
⎥
⎢
⎥ ⎢
⎥⎢
⎥ ⎢
⎢ xˆn −1 (k + 1) ⎥ ⎢ − an −1 0 ... 0 1 ⎥ ⎢ xˆn −1 (k ) ⎥ ⎢ b2 − a2b0 ⎥
⎢⎣ xˆn (k + 1) ⎥⎦ ⎢⎣ − an
0 ... 0 0 ⎥⎦ ⎢⎣ xˆn (k ) ⎥⎦ ⎢⎣ b1 − a1b0 ⎥⎦
71
⎡ xˆ1 (k ) ⎤
⎢ xˆ (k ) ⎥
⎢ 2
⎥
y (k ) = [1 0 .... 0 0] ⎢ .... ⎥ + b0u (k ) .
⎢
⎥
⎢ xˆn −1 (k ) ⎥
⎢⎣ xˆn (k ) ⎥⎦
(5-15)
Forma canônica diagonal
Se os pólos da função de transferência de pulso dada pelas eqs. (5-6) e (5-7) são todos
distintos, então a representação em estado de estado pode ser posta na forma diagonal:
0 ⎤ ⎡ x1 (k ) ⎤ ⎡ 1 ⎤
⎡ x1 (k + 1) ⎤ ⎡ p1 0 ... 0
⎢ x (k + 1) ⎥ ⎢ 0 p ... 0
0 ⎥⎥ ⎢⎢ x2 (k ) ⎥⎥ ⎢ 1 ⎥
2
⎢ 2
⎥ ⎢
⎢ ⎥
⎢ ..... ⎥ = ⎢.... .... ... .... .... ⎥ ⎢ .... ⎥ + ⎢....⎥ u (k ) ,
(5-16)
⎢
⎥ ⎢
⎥⎢
⎥ ⎢ ⎥
⎢ xn −1 (k + 1) ⎥ ⎢ 0 0 ... pn −1 0 ⎥ ⎢ xn −1 (k ) ⎥ ⎢ 1 ⎥
⎢⎣ xn (k + 1) ⎥⎦ ⎢⎣ 0 0 ... 0
pn ⎥⎦ ⎢⎣ xn (k ) ⎥⎦ ⎢⎣ 1 ⎥⎦
⎡ x1 (k ) ⎤
⎢ x (k ) ⎥
⎢ 2
⎥
(5-17)
y (k ) = [ c1 c2 ... cn −1 cn ] ⎢ .... ⎥ + b0u (k ) ,
⎢
⎥
⎢ xn −1 (k ) ⎥
⎢⎣ xn (k ) ⎥⎦
onde os pi ' s são os polos da função de transferência e os ci ' s são os resíduos da
expansão da função de transferência em frações parciais, isto é,
⎛ Y ( z)
⎞
( z − pi ) ⎟ . Este resultado também deveria ter sido provado no curso
ci = lim ⎜
z → pi U ( z )
⎝
⎠
anterior.
Forma canônica de Jordan
Se a função de transferência de pulso tiver um pólo de ordem m em z = p1 e todos os
outros pólos forem distintos, então obtemos as seguintes eqs. de estado e de resposta:
⎡ x1 (k + 1) ⎤ ⎡ p1
⎢ x (k + 1) ⎥ ⎢ 0
⎢ 2
⎥ ⎢
⎢
⎥ ⎢ ...
....
⎢
⎥ ⎢
(
1)
x
k
+
m
⎢
⎥=⎢0
⎢ xm+1 (k + 1) ⎥ ⎢ 0
⎢
⎥ ⎢
....
⎢
⎥ ⎢ ...
⎢ x (k + 1) ⎥ ⎢ 0
⎣ n
⎦ ⎣
1 0 ... 0
p1 1 ... 0
... ... ... ...
0
0
0 ...
0 ...
p1
0
... ... ... ...
0 0 0 0
0
0
...
...
...
...
0 ...
pm+1 ...
...
0
...
0
0 ⎤ ⎡ x1 ( k ) ⎤ ⎡ 0 ⎤
0 ⎥⎥ ⎢⎢ x2 ( k ) ⎥⎥ ⎢ 0 ⎥
⎢ ⎥
... ⎥ ⎢ .... ⎥ ⎢...⎥
⎥⎢
⎥ ⎢ ⎥
0 ⎥ ⎢ xm (k ) ⎥ + ⎢ 0 ⎥ u (k ) , (5-18)
0 ⎥ ⎢ xm+1 (k ) ⎥ ⎢ 1 ⎥
⎥⎢
⎥ ⎢ ⎥
... ⎥ ⎢ .... ⎥ ⎢ 1 ⎥
pn ⎦⎥ ⎢⎣ xn (k ) ⎥⎦ ⎢⎣ 1 ⎥⎦
onde a matriz acima tem quatro blocos: o primeiro superior à esquerda tem p1 ' s na
diagonal e 1’s na primeira sobre-diagonal, os outros elementos nulos; o segundo bloco
superior à direita tem todos os elementos nulos; o primeiro boco inferior à esquerda
também tem todos os elementos nulos e finalmente o segundo bloco inferior à direita
tem os outros pólos da função de transferência na diagonal e os outros elementos nulos.
A eq. da resposta é, tal como no caso anterior:
72
⎡ x1 (k ) ⎤
⎢ x (k ) ⎥
⎢ 2
⎥
(5-19)
y (k ) = [ c1 c2 ... cn −1 cn ] ⎢ .... ⎥ + b0u (k ) .
⎢
⎥
⎢ xn −1 (k ) ⎥
⎢⎣ xn (k ) ⎥⎦
A matriz de (5-18) é uma forma canônica de Jordan. O bloco com 1’s na sobre-diagonal
é chamado um bloco de Jordan. Se a função de transferência de pulso tiver q pólos com
grau de multiplicidade maior que 1, a matriz terá q blocos de Jordan.
Y ( z)
z +1
= 2
. Achar as formas canônicas
U ( z ) z + 1,3z + 0, 4
controlável, observável e diagonal (ou de Jordan, se for o caso) para este sistema.
Solução: A forma canônica controlável é escrita por simples inspeção da função de
transferência:
1 ⎤ ⎡ x1 (k ) ⎤ ⎡0 ⎤
⎡ x1 ( k + 1) ⎤ ⎡ 0
⎡ x1 (k ) ⎤
=
+
u
k
=
(
)
y
(
k
)
1
1
,
[
]
⎢ x (k + 1) ⎥ ⎢ −0, 4 −1,3⎥ ⎢ x (k ) ⎥ ⎢1 ⎥
⎢ x (k ) ⎥ .
⎦⎣ 2 ⎦ ⎣ ⎦
⎣ 2
⎦ ⎣
⎣ 2 ⎦
A forma canônica observável também é escrita por simples inspeção e é transposta da
anterior:
⎡ x1 (k + 1) ⎤ ⎡0 −0, 4 ⎤ ⎡ x1 (k ) ⎤ ⎡1⎤
⎡ x1 (k ) ⎤
⎢ x (k + 1) ⎥ = ⎢1 −1,3 ⎥ ⎢ x (k ) ⎥ + ⎢1⎥ u (k ) , y (k ) = [ 0 1] ⎢ x (k ) ⎥ .
⎦⎣ 2 ⎦ ⎣ ⎦
⎣ 2
⎦ ⎣
⎣ 2 ⎦
Exemplo 5-1: Seja o sistema
Para escrever a forma canônica diagonal, precisamos achar os pólos e os respectivos
resíduos da função de transferência. Neste caso, obtemos
Y ( z)
5/3
−2 / 3
. Donde temos imediatamente:
=
+
U ( z ) z + 0,5 z + 0,8
0 ⎤ ⎡ x1 (k ) ⎤ ⎡1⎤
⎡ x1 (k + 1) ⎤ ⎡ −0,5
⎡ x1 (k ) ⎤
,
=
+
u
k
=
−
(
)
y
(
k
)
5
/
3
2
/
3
[
]
⎢ x (k + 1) ⎥ ⎢ 0
⎢ x (k ) ⎥ .
−0,8⎥⎦ ⎢⎣ x2 (k ) ⎥⎦ ⎢⎣1⎥⎦
⎣ 2
⎦ ⎣
⎣ 2 ⎦
As representações no espaço de estado não são únicas
Seja o sistema representado por
x(k + 1) = Gx(k ) + Hu(k ) ,
(5-1bis)
y (k ) = Cx(k ) + Du(k ) .
(5-2 bis)
Definamos um novo vetor xˆ (k ) por meio de
x(k ) = Pxˆ (k ) ,
(5-22)
onde P é uma matriz não singular, a saber, inversível e, consequentemente, com o
determinante diferente de zero .
Substituindo esta em (5-1bis), vem:
Pxˆ (k + 1) = GPxˆ (k ) + Hu(k ) .
(5-23)
−1
Multiplicando ambos os lados desta última por P , vem
xˆ (k + 1) = P −1GPxˆ (k ) + P −1Hu(k )
(5-24)
ˆ , P −1H = H
ˆ .
Definamos P −1GPxˆ = G
Então a eq. (5-24) pode ser escrita na forma
ˆ ˆ (k ) + Hu
ˆ (k ) .
xˆ (k + 1) = Gx
(5-25)
Por outro lado, substituindo (5-22) em (5-2bis), temos
y (k ) = CPxˆ (k ) + Du(k ) .
(5-26)
73
ˆ e D=D
ˆ , temos
Definindo CP = C
ˆ ˆ (k ) + Du
ˆ (k ) ,
y (k ) = Cx
(5-27)
Então acabamos de mostrar que a realização (5-1bis) e (5-2bis) é equivalente à
realização dada por (5-25) e (5-27), pois ambas correspondem ao mesmo par de entrada
e saída u(k ) e y (k ) .
Tendo em vista que a matriz P é qualquer (mas não singular), conclui-se que existe
uma infinidade de realizações (ou representações) no espaço de estado de um mesmo
sistema.
5-3 Resolvendo equações discretas no espaço de estado
Nesta seção resolveremos a eq. (5-1bis) por dois métodos, recurrência e transformada z.
Em geral, as eqs. de diferença são mais fáceis de serem resolvidas do que as eqs.
diferenciais.
Vejamos primeiramente o método da recurrência:
Da eq. (5-1bis) temos imediatamente:
x(1) = Gx(0) + Hu(0) ,
x(2) = Gx(1) + Hu(1) = G 2 x(0) + GHu(0) + Hu(1) ,
x(3) = Gx(2) + Hu(2) = G 3 x(0) + G 2 Hu(0) + GHu(1) + Hu(2) ,
....................................................................................................
A generalização é imediata:
k −1
x(k ) = G k x(0) + ∑ G k − j −1Hu( j ) , k = 1;2;3;....
(5-30)
j =0
E a resposta do sistema é dada por (5-2bis), ou seja,
k −1
y (k ) = CG k x(0) + C∑ G k − j −1Hu( j ) + Du(k )
(5-31)
j =0
Matriz de transição de estado
Considere a eq. homogênea de (5-1bis):
x(k + 1) = Gx(k )
(5-32)
De (5-30) com u(k ) = 0, temos
x(k ) = G k x(0) .
A matriz G k é chamada matriz de transição de estado e também matriz fundamental e
costuma ser representada por Ψ (k ) . (Há uma diferença entre matriz de transição de
estado e matriz fundamental em sistemas variantes no tempo). Ou seja, temos, da eq.
logo acima:
x(k ) = Ψ(k )x(0) ,
(5-33)
onde Ψ (k ) satisfaz evidentemente à condição
Ψ (k + 1) = GΨ (k ) , com Ψ (0) = I
(5-34)
Em termos de matriz de transição de estado a eq. (5-30) pode ser escrita na seguinte
forma:
k −1
x(k ) = Ψ ( k )x(0) + ∑ Ψ (k − j − 1)Hu( j ) .
(5-36)
j =0
O somatório acima é uma convolução discreta, que também pode ser escrita na forma:
k −1
x(k ) = Ψ ( k )x(0) + ∑ Ψ ( j )Hu(k − j − 1) .
(5-37)
j =0
74
Substituindo estas eqs. em (5-31), temos:
k −1
y (k ) = CΨ (k )x(0) + C∑ Ψ (k − j − 1)Hu( j ) + Du(k )
(5-38)
j =0
k −1
= CΨ (k )x(0) + C∑ Ψ ( j ) Hu( k − j − 1) + Du( k ) .
(5-39)
j =0
Método da transformada z para a solução de equações de estado discretas
Consideremos de novo a eq.
x(k + 1) = Gx(k ) + Hu(k ) .
Tomando a transformada z de ambos os lados desta eq., temos
zX( z ) − zx(0) = GX( z ) + HU( z ) , ou seja,
( zI − G ) X( z ) = zx(0) + HU( z ) . Donde,
X( z ) = ( zI − G ) −1 zx(0) + ( zI − G ) −1 HU( z ) .
Tomando a transformada inversa desta, temos
x(k ) = Z-1[( zI − G ) −1 z ]x(0) + Z-1[( zI − G ) −1 HU( z )] .
Comparando esta com (5-30), temos
G k = Z-1[( zI − G ) −1 z ]
k −1
e
∑G
k − j −1
Hu( j ) = Z-1[( zI − G ) −1 HU( z )] ,
(5-1bis)
(5-41)
(5-42)
(5-43)
(5-44)
j =0
onde k = 1; 2; 3;...
Exemplo 5-2: Obter a matriz de transição de estado do sistema (5-1bis) e (5-2bis), com
1⎤
⎡ 0
⎡1⎤
G=⎢
, H = ⎢ ⎥ , C = [1 0] e D = 0 . A seguir, calcule o estado x(k ) e a
⎥
⎣ −0,16 −1⎦
⎣1⎦
⎡ x (0) ⎤ ⎡ 1 ⎤
resposta y (k ) com u (k ) = 1 para k = 0;1;2;... e o estado inicial x(0) = ⎢ 1 ⎥ = ⎢ ⎥ .
⎣ x2 (0) ⎦ ⎣ −1⎦
Solução: De (5-43), temos Ψ (k ) = G k = Z-1[( zI − G ) −1 z ] .
z +1
1
⎡
⎤
−1
⎢
−1 ⎤
( z + 0, 2)( z + 0,8) ( z + 0, 2)( z + 0,8) ⎥
⎡ z
⎢
⎥
=
Ora, ( zI − G ) −1 = ⎢
⎥
z
−0,16
⎢
⎥
⎣ 0,16 z + 1⎦
⎢ ( z + 0, 2)( z + 0,8) ( z + 0, 2)( z + 0,8) ⎥
⎣
⎦
−1/ 3
⎡ 4/3
+
⎢ z + 0, 2 z + 0,8
=⎢
⎢ −0,8 / 3 0,8 / 3
⎢ z + 0, 2 + z + 0,8
⎣
5/3
+
z + 0, 2
−1/ 3
+
z + 0, 2
−5 / 3 ⎤
z + 0,8 ⎥
⎥ . Donde,
4/3 ⎥
z + 0,8 ⎥⎦
1 z
5 z
5 z
⎡ 4 z
⎤
−
⎢ 3 z + 0, 2 − 3 z + 0,8
3 z + 0, 2 3 z + 0,8 ⎥
⎥
Ψ (k ) = G k = Z-1[( zI-G ) −1 z ] = Z-1 ⎢
0,8 z
1 z
4 z ⎥
⎢ 0,8 z
⎢ − 3 z + 0, 2 + 3 z + 0,8 − 3 z + 0, 2 + 3 z + 0,8 ⎥
⎣
⎦
75
1
5
5
⎡ 4
⎤
k
k
k
(
0,
2)
(
0,8)
(
0,
2)
(−0,8) k ⎥
−
−
−
−
−
⎢ 3
3
3
3
=⎢
(5-45)
⎥
⎢ − 0,8 (−0, 2) k + 0,8 (−0,8) k − 1 (−0, 2) k + 4 (−0,8) k ⎥
⎢⎣ 3
⎥⎦
3
3
3
Agora vamos calcular o estado: X( z ) = ( zI − G ) −1 zx(0) + ( zI − G ) −1 HU ( z )
= ( zI − G ) −1[ zx(0) + HU ( z )] . É dado do enunciado do problema que
1
z
. Ora, X( z ) = ( zI − G ) −1[ zx(0) + HU ( z )] ; mas
=
U ( z) =
−1
1− z
z −1
2
⎤
⎡ z ⎤ ⎡ z
⎢
⎥
⎢
⎥
⎡ z ⎤ z −1
z −1 ⎥ .
zx(0) + HU ( z ) = ⎢ ⎥ + ⎢
⎥=⎢ 2
⎣− z ⎦ ⎢ z ⎥ ⎢ − z + 2 z ⎥
⎣⎢ z − 1 ⎦⎥ ⎢⎣ z − 1 ⎦⎥
22
25 ⎤
⎡ 17
z
z⎥
⎢ − 6 z
⎡
⎤
( z 2 + 2) z
9
18
+
+
⎢
⎥
⎢ ( z + 0, 2)( z + 0,8)( z − 1) ⎥
z + 0, 2 z + 0,8 z − 1 ⎥
⎢
⎢
⎥
=
E, portanto, X( z ) =
; donde,
⎢ 3, 4
⎢
⎥
17, 6
7 ⎥
(− z 2 + 1,84 z ) z
⎢ 6 z − 9 z 18 z ⎥
⎢
⎥
z
+
z
+
z
−
(
0,
2)(
0,8)(
1)
⎣
⎦
+
+
⎢
⎥
⎣ z + 0, 2 z + 0,8 z − 1 ⎦
22
25 ⎤
⎡ 17
− (−0, 2) k + (−0,8) k +
⎢
6
9
18 ⎥
x(k ) = Z-1[ X( z )] = ⎢
⎥,
3,
4
17,
6
7
k
k
⎢
(−0, 2) −
(−0,8) + ⎥
⎢⎣ 6
9
18 ⎥⎦
22
25 ⎤
⎡ 17
k
k
⎢ − 6 (−0, 2) + 9 (−0,8) + 18 ⎥
y (k ) = Cx(k ) = [1 0] ⎢
⎥
⎢ 3, 4 (−0, 2) k − 17, 6 (−0,8) k + 7 ⎥
⎢⎣ 6
9
18 ⎥⎦
17
22
25
= − (−0, 2) k + (−0,8) k + .
6
9
18
Cálculo de ( zI − G ) −1
Este cálculo consome usualmente, se for feito manualmente, bastante tempo. Há
métodos analíticos e computacionais para calcular ( zI − G ) −1 .
Como é bem sabido, ( zI − G ) −1 pode ser expresso como função da adjunta de zI − G e
do determinante da mesma, a saber,
adj ( zI − G )
(5-46)
( zI − G ) −1 =
zI − G
Note-se que o determinante zI − G pode ser expresso como
zI − G = z n + a1 z n −1 + a2 z n − 2 + ...an
Pode-se demonstrar que
adj ( zI − G ) = Iz n −1 + H1 z n − 2 + H 2 z n −3 + ... + H n −1 , onde
H1 = G + a1I
H 2 = GH1 + a2 I
(5-47)
(5-48)
76
.......................
(5-49)
H n −1 = GH n − 2 + an −1I
H n = GH n −1 + an I = 0 .
Os coeficientes ai ' s podem ser obtidos ou pelo cálculo do determinante (5-47) ou pelas
fórmulas (lembrando que o traço de uma matriz é a soma dos elementos da sua diagonal
principal):
1
1
1
a1 = −tr (G ) , a2 = − tr (GH1 ) , a3 = − tr (GH 2 ) ,..., an = − tr (GH n −1 ) ,
(5-50)
n
2
3
Procedemos sequencialmente da seguinte forma:
1
a1 = −tr (G ) , H1 = G + a1I , a2 = − tr (GH1 ) , H 2 = GH1 + a2 I , etc.
2
Conhecidos os ai ' s e os H i ' s , calcula-se ( zI − G ) −1 , usando (5-46) - (5-48).
0 ⎤
⎡ 0,1 0,1
⎢
Exemplo 5-3: Determinar ( zI − G ) , onde G = ⎢0,3 −0,1 −0, 2 ⎥⎥ . Obter também G k .
⎢⎣ 0
−0,3 ⎥⎦
0
Solução: Neste caso o determinante é de ordem 3 e pode ser obtido facilmente. Mas
usaremos o método alternativo proposto, para treinar.
Da primeira eq. de (5-50), temos a1 = −tr (G ) = -(0,1-0,1-0,3) = 0,3.
−1
0 ⎤
⎡ 0,1 0,1
⎢
Da primeira eq. de (5-49), temos H1 = G + a1I = ⎢ 0,3 −0,1 −0, 2 ⎥⎥ +
⎢⎣ 0
−0,3 ⎥⎦
0
0 ⎤
⎡ 0,3 0
⎢ 0 0,3 0 ⎥
⎢
⎥
⎢⎣ 0
0 0,3⎥⎦
0 ⎤
⎡ 0, 4 0,1
= ⎢⎢ 0,3 0, 2 −0, 2 ⎥⎥ .
⎢⎣ 0
0
0 ⎥⎦
1
Da 2ª. eq. de (5-50), temos a2 = − tr (GH1 ) =
2
⎧ ⎡ 0,1 0,1
0 ⎤ ⎡0, 4 0,1
0 ⎤⎫
⎡ 0, 07 0, 03 −0, 02 ⎤
1 ⎪⎢
1 ⎢
⎪
⎥
⎢
⎥
− tr ⎨ ⎢0,3 −0,1 −0, 2 ⎥ ⎢ 0,3 0, 2 −0, 2 ⎥ ⎬ = − tr ⎢ 0, 09 0, 01 0, 02 ⎥⎥ = -0,04.
2 ⎪
2
⎪
⎢
⎥
⎢
⎥
⎢⎣ 0
0
0 ⎥⎦
0
0
−
0,3
0
0
0
⎣
⎦
⎣
⎦
⎩
⎭
Mas da 2ª. eq. (5-49), temos H 2 = GH1 + a2 I ; substituindo os valores obtidos, temos:
⎡ 0, 03 0, 03 −0, 02 ⎤
H 2 = ⎢⎢ 0, 09 −0, 03 0, 02 ⎥⎥ . E então, de novo de (5-50), temos a3 = -0,012.
⎢⎣ 0
0
−0, 04 ⎥⎦
Com estes elementos, já podemos calcular ( zI − G ) −1 . Mas é bom conferir a última eq.
(5-49), a saber, GH 2 − 0, 012I = 0 .
De (5-46)-(5-48), vem:
77
z + 0,1
0,1
−0, 02
⎡
⎤
⎢ ( z + 0, 2)( z − 0, 2) ( z + 0, 2)( z − 0, 2) ( z + 0, 2)( z − 0, 2) ⎥
⎢
⎥
z
z
⎢
0,3
−
0,1
−
0,
2(
−
0,1)
⎥
( zI − G ) −1 = ⎢
. Expandindo
( z + 0, 2)( z − 0, 2) ( z + 0, 2)( z − 0, 2)
( z + 2)( z − 2) ⎥
⎢
⎥
⎢
⎥
1
0
0
⎢
⎥
z + 0,3
⎣
⎦
cada uma das funções racionais e usando a Tabela 2-1, se for necessário, achamos
G k = Z-1[( zI − G ) −1 z ]
⎡ 0, 25(−0, 2) k + 0, 75(0, 2) k
⎢
= ⎢ −0, 75(−0, 2) k + 0, 75(0, 2) k
⎢⎣
0
−0, 25(−0, 2) k + 0, 25(0, 2) k
0, 5(−0, 2) k − 0,1(0, 2) k − 0, 4( −0, 3) k ⎤
0, 75(−0, 2) k + 0, 25(0, 2) k
−1, 5( −0, 2) k − 0,1(0, 2) k + 1, 6(−0, 3) k ⎥
0
⎥
⎥⎦
( −0, 3) k
5-4 MATRIZ DE TRANFERÊNCIA DE PULSO
Consideremos novamente as eqs. de estado
x(k + 1) = Gx(k ) + Hu(k ) ,
y (k ) = Cx(k ) + Du(k ) .
De (5-42), com x(0) = 0 , obtemos
(5-1bis)
(5-2bis)
X( z ) = ( zI − G ) −1 HU( z ) ,
Y( z ) = [C( zI − G ) −1 H + D]U( z ) = F ( z ) , onde
F( z ) = C( zI − G ) −1 H + D
(5-60)
é a chamada matriz de função de transferência de pulso. Trata-se evidentemente de uma
matriz m × r, lembrando as dimensões das matrizes em questão.
adj ( zI − G )
. Donde que,
Mas sabemos que ( zI − G ) −1 =
zI − G
F( z ) =
Cadj ( zI − G )H
+D.
zI − G
É claro que os pólos de F( z ) são as raízes de zI − G = 0, esta sendo então a eq.
característica do sistema de tempo discreto, ou seja,
z n + a1 z n −1 + a2 z n − 2 + ...an −1 z + an = 0.
Transformações de similaridade
Na seção 5.2 vimos que há um número infinito de representações em estado de estado
de um mesmo sistema, ou seja, se definirmos
x(k ) = Pxˆ (k ) ,
sendo P uma matriz não singular, obtemos o sistema
ˆ ˆ (k ) + Hu
ˆ (k )
xˆ (k + 1) = Gx
(5-61)
ˆ ˆ (k ) + Du
ˆ (k ) ,
y (k ) = Cx
(5-62)
onde, como vimos,
ˆ ; P −1H = H
ˆ e D = D̂ .
ˆ ; CP = C
P −1GPx = G
Ora, a função de transferência de pulso de (5-61) e (5-62) é
ˆ ( zI − G
ˆ ) −1 H
ˆ +D
ˆ.
Fˆ ( z ) = C
78
Mas substituindo os valores acima, temos
Fˆ ( z ) = CP( zI − P −1GP) −1 P −1H + D = C( zI − G ) −1 H + D = F( z ) .
Ou seja, acabamos de demonstrar que a matriz de transferência não se altera quando
fazemos uma transformação de similaridade, ou seja, transformação de estado da forma
x(k ) = Pxˆ (k ) , com P não singular, um resultado já esperado, uma vez que a entrada e a
resposta do sistema não se alteram quando fazemos transformação de estado.
E é fácil verificar que o polinômio característico não se altera, pois
ˆ = P −1 zI − G P = zI − G , ou seja, os auto-valores de G são invariantes sob
zI − G
transformação de similaridade.
5-5 Discretização de equações de espaço de estado de tempo contínuo
Recordamos que
∞
1
1
Ak t k
.
e At = I + At + A 2t 2 + .... + A k t k + .... = ∑
2!
k!
k =0 k !
A série convergindo, pode ser diferenciada termo a termo, dando:
d At
A 3t 2
A k t k −1
+ .... +
+ ...
e = A + A 2t +
2!
(k − 1)!
dt
(*)
⎛
⎞
1
A k t k −1
= A ⎜ I + At + A 2t 2 + .... +
+ ... ⎟ = Ae At
2!
(k − 1)!
⎝
⎠
⎛
⎞
1
A k t k −1
= ⎜ I + At + A 2t 2 + .... +
+ ... ⎟ A = e At A .
2!
(k − 1)!
⎝
⎠
A seguinte propriedade, bem conhecida das funções exponenciais, também vale para
exponenciais de matrizes:
e At e As = e A (t + s ) .
Com efeito,
At As
e e
⎛ ∞ A k t k ⎞⎛ ∞ A k s k
= ⎜∑
⎟⎜ ∑
⎝ k =0 k ! ⎠⎝ k =0 k !
⎞ ∞ k ⎛ k t i s k −i ⎞ ∞ k (t + s ) k
= e A (t + s ) .
⎟ = ∑A ⎜∑
⎟ = ∑A
k!
⎠ k =0 ⎝ i =0 i !(k − i )! ⎠ k =0
(Verificar a segunda e terceira igualdades, usando (*) acima).
É claro que se tivermos s = -t, vem imediatamente
e At e − At = e− At e At = e A ( t −t ) = e A 0 = e0 = I .
É fácil provar também que e( A + B )t = e At eBt se AB = BA e, em geral,
e( A + B ) t ≠ e At eBt se AB ≠ BA .
Passemos ao problema desta seção, isto é, a discretização do sistema de tempo contínuo.
O sinal de tempo contínuo é amostrado a cada T unidades de tempo e mantido o seu
valor até a próxima amostragem, ou seja, trata-se de um amostrador de ordem zero,
seguido de um segurador (holder).
Neste ponto recorda-se que no caso de sistemas de tempo contínuo, temos:
x (t ) = Ax(t ) + Bu(t ) ,
(5-63)
a solução sendo dada por:
t
x(t ) = e At x(0) + ∫ e A (t −τ ) Bu(τ )dτ ,
(5-64)
0
onde x(0) é obviamente o estado inicial.
Temos então
79
u(t ) = u(kT ) para kT ≤ t < kT + T .
Ora, de (5-64) com t = (k + 1)T e t = kT vem, respectivamente:
( k +1)T
∫
x((k + 1)T ) = e A ( k +1)T x(0) + e A ( k +1)T
e− Aτ Bu(τ )dτ ,
(5-69)
(5-70)
0
kT
e x(kT ) = e AkT x(0) + e AkT ∫ e − Aτ Bu(τ )dτ .
(5-71)
0
Multiplicando esta última eq. por e AT e subtraindo de (5-70), temos, tendo em vista os
limites de integração das integrais:
( k +1)T
x((k + 1)T ) = e AT x(kT ) + e A ( k +1)T
∫
e − Aτ Bu(τ )dτ .
kT
Mas em vista de (5-69), u(τ ) = u(kT ) = constante no integrando desta última equação.
Donde, fazendo t = τ − kT , temos
T
x((k + 1)T ) = e AT x(kT ) + e AT ∫ e − At Bu(kT )dt .
0
Definindo agora λ = T − t , vem:
0
T
T
0
x((k + 1)T ) = e AT x(kT ) − ∫ e Aλ Bu(kT )d λ = e AT x(kT ) + ∫ e Aλ Bu(kT )d λ .
Definamos G (T ) = e AT
(5-72)
(5-73)
⎛T
⎞
(5-74)
H (T ) = ⎜ ∫ e Aλ d λ ⎟ B ,
⎝0
⎠
com o que a eq. (5-72) se torna:
x((k + 1)T ) = G (T )x(kT ) + H (T )u(kT ) .
(5-75)
E, claro,
y (kT ) = Cx(kT ) + Du(kT ) .
(5-76)
Finalmente, observe-se que se a matriz A for não singular, então de (5-74) temos
também:
⎛T
⎞
H (T ) = ⎜ ∫ e Aλ d λ ⎟ B = A −1 (e AT − I )B = (e AT − I ) A −1B .
⎝0
⎠
A eq. de estado (5-75) é dita equivalente com amostrador de ordem zero à eq. de estado
original (5-63). Note-se ainda de (5-73) que se T 1 , então G (T ) ≈ G (0) = e A 0 = I .
e
Y ( s)
1
=
. Obter a representação de
U ( s) s + a
estado deste sistema e em seguida discretize-o, isto é, obtenha a representação de estado
discreta do sistema, bem como a função de transferência.
Solução: Obtemos imediatamente
x = −ax + u e y = x . (Neste caso as formas canônica controlável, observável e
diagonal coincidem).
Ora, das eqs. (5-73) e (5-74), temos
Exemplo 5-4: Seja o sistema dado por G ( s ) =
G (T ) = e
− aT
T
e H (T ) = ∫ e
0
− aλ
1 − e − aT
.
dλ =
a
Portanto a versão discreta das eqs. de estado é:
80
x(k + 1) = e − aT x(k ) +
1 − e − aT
u (k ) ,
a
y (k ) = x(k ) .
E em vista de (5-60), a função de transferência de pulso para este sistema é:
1 − e − aT
(1 − e − aT ) z −1
.
F ( z ) = C ( zI − G ) −1 H = ( z − e − aT ) −1
=
a
a (1 − e − aT z −1 )
Podemos conferir este resultado calculando a função de transferência pelo método
usado nos capítulos anteriores, a saber, calculando a função de transferência da planta
precedida de um segurador de ordem zero:
⎛ 1 − e −Ts 1 ⎞
⎛ 1 ⎞ (1 − e − aT ) z −1
−1
G( z) = Z ⎜
= F ( z ) , conferindo.
⎟ = (1 − z )Z ⎜
⎟=
− aT −1
⎝ s ( s + a ) ⎠ a (1 − e z )
⎝ s s+a⎠
Exemplo 5-5: Obter as eqs. de tempo discreto, a resposta do sistema e a função de
transferência de pulso, com período de amostragem T = 1 do seguinte sistema de tempo
Y ( s)
1
=
.
contínuo: G ( s ) =
U ( s ) s ( s + 2)
Solução: O sistema acima pode ser representado pelas seguintes eqs., como se pode
⎡ x ⎤ ⎡0 1 ⎤ ⎡ x1 ⎤ ⎡ 0 ⎤
⎡x ⎤
verificar facilmente: ⎢ 1 ⎥ = ⎢
+ ⎢ ⎥ u , y = [1 0] ⎢ 1 ⎥ .
⎢
⎥
⎥
⎣ x2 ⎦ ⎣0 −2 ⎦ ⎣ x2 ⎦ ⎣1 ⎦
⎣ x2 ⎦
A eq. de estado a ser obtida tem a forma: x((k + 1)T ) = G (T )x(kT ) + H (T )u (kT ) ,
onde G (T ) e H (T ) são dadas pelas eqs. (5-73) e (5-74). Obtém-se sem muita
dificuldade: G (T ) = e
H (T ) =
(∫
T
0
)
AT
⎡1 0,5(1 − e −2T ) ⎤
=⎢
⎥,
e −2T
⎣0
⎦
⎛ T ⎡1 0,5(1 − e −2T ) ⎤ ⎞ ⎡0 ⎤
AT
e dt B = ⎜ ∫ ⎢
⎥dt ⎟⎟ ⎢ ⎥ =
⎜ 0 0
e −2T
⎦ ⎠ ⎣1 ⎦
⎝ ⎣
⎡ ⎛
e −2T − 1 ⎞ ⎤
+
0,5
T
⎢ ⎜
⎟⎥
2 ⎠ ⎥ . Então,
⎢ ⎝
⎢⎣ 0,5(1 − e −2T ) ⎥⎦
⎡ ⎛
e −2T − 1 ⎞ ⎤
⎡ x1 ((k + 1)T ) ⎤ ⎡1 0,5(1 − e −2T ) ⎤ ⎡ x1 (kT ) ⎤ ⎢0,5 ⎜ T +
⎟⎥
2 ⎠ ⎥ u (kT ) ,
⎥⎢
⎢ x ((k + 1)T ) ⎥ = ⎢
⎥+⎢ ⎝
−2T
x
(
kT
)
e
⎣ 2
⎦ ⎣0
⎦ ⎢
⎦⎣ 2
−2T
⎣ 0,5(1 − e ) ⎥⎦
⎡ x (kT ) ⎤
y (kT ) = [1 0] ⎢ 1
⎥.
⎣ x2 (kT ) ⎦
Quando fazemos T = 1, estas eqs. se tornam
⎡ x1 (k + 1) ⎤ ⎡1 0, 4323⎤ ⎡ x1 (k ) ⎤ ⎡0, 2838⎤
⎡ x1 (k ) ⎤
⎢ x (k + 1) ⎥ = ⎢0 0,1353 ⎥ ⎢ x (k ) ⎥ + ⎢ 0, 4323⎥ u (k ) e y (k ) = [1 0] ⎢ x (k ) ⎥ .
⎦⎣ 2 ⎦ ⎣
⎦
⎣ 2
⎦ ⎣
⎣ 2 ⎦
A função de transferência de pulso é obtida da eq. (5-60):
F ( z ) = C( zI − G ) −1 H + D
81
0, 4323
⎡ 1
⎤
⎢
z − 1 ( z − 1)( z − 0,1353) ⎥ ⎡ 0, 2838⎤
⎡ z − 1 −0, 4323 ⎤ ⎡ 0, 2838⎤
⎢
⎥⎢
=
= [1 0] ⎢
1
0
[
]
⎥ ⎢ 0, 4323⎥
⎥
z
−
0
0,1353
1
⎢
⎥ ⎣ 0, 4323⎦
⎣
⎦ ⎣
⎦
0
⎢
⎥
z − 0,1353
⎣
⎦
−1
−2
0, 2838 z + 0,1485 0, 2838 z + 0,1485 z
.
=
=
( z − 1)( z − 0,1353) (1 − z −1 )(1 − 0,1353 z −1 )
−1
Método de discretização de sistemas de tempo contínuo usando o MATLAB
Dado o sistema de tempo contínuo x (t ) = Ax(t ) + Bu(t ) , trata-se de obter o sistema de
tempo discreto correspondente x(k + 1) = Gx(k ) + Hu(k ) .
O comando do MATLAB correspondente é:
[G,H]=c2d(A,B,T) ,
onde T é o período de amostragem, o qual deve ser especificado em segundos. Se se
quer quatro casas decimais, use-se o format short; caso se queira maior precisão, usar o
format long. Se não se especifica o formato, o MATLAB usa o format short.
Vejamos o exemplo em que o sistema de tempo contínuo é dado por
1 ⎤ ⎡ x1 ⎤ ⎡0 ⎤
⎡ x1 ⎤ ⎡ 0
(5-77)
⎢ x ⎥ = ⎢ −25 −4 ⎥ ⎢ x ⎥ + ⎢1 ⎥ u
⎦⎣ 2⎦ ⎣ ⎦
⎣ 2⎦ ⎣
Os comandos para obter G e H estão no retângulo abaixo:
Como se vê no quadro acima, o tempo de amostragem usado foi T = 0,05 seg. Nos
retângulos seguintes os períodos de amostragem são 0,2 e 1 seg., respectivamente.
Como se verifica, as matrizes G e H mudam bastante quando se muda muito o tempo
de amostragem.
1 0 0⎤
⎡ 0
⎡ 0 ⎤
⎢ 20, 601 0 0 0 ⎥
⎢ ⎥
⎥ , B = ⎢ −1 ⎥ . Com o período
Vejamos outro exemplo, em que A = ⎢
⎢ 0
⎢ 0 ⎥
0 0 1⎥
⎢
⎥
⎢ ⎥
⎣ −0, 4905 0 0 0 ⎦
⎣0,5⎦
82
de amostragem T = 0,05, obtemos o seguinte resultado de acordo com o retângulo
abaixo:
Ao discretizarmos um sistema de tempo contínuo, obtemos o estado e a resposta nos
instantes de amostragem. Ora, muitas vezes estamos também interessados no valor do
estado e da resposta entre os instantes de amostragem. É disso que trataremos no final
deste capítulo.
Seja o sistema de tempo contínuo:
x = Ax + Bu, y = Cx + Du .
Suponhamos que a entrada u passe por um amostrador e um segurador de ordem zero.
Vimos que u(τ ) = u(kT ) para kT ≤ τ < kT + T .
Sabemos que a solução do sistema de tempo contínuo acima é:
t
x(t ) = e A (t −t0 ) x(t0 ) + ∫ e A (t −τ ) Bu(τ )dτ , onde t0 é o instante inicial.
t0
Para obtermos o valor do estado entre os instantes de amostragem, façamos na
expressão acima
t = kT + ∆T , onde 0 < ∆T < T e t0 =kT. Obtemos:
x(kT + ∆T ) = e A∆T x(kT ) +
kT +∆T
∫
e A ( kT +∆T −τ ) Bu(kT )dτ .
kT
Definamos λ = kT + ∆T − τ .
Donde que quando τ = kT , λ = ∆T e quando τ = kT + ∆T , λ = 0 ; por outro lado,
d λ = − dτ . Portanto, temos
x(kT + ∆T ) = e
A∆T
0
x(kT ) −
∫e
Aλ
Bu(kT )d λ = e
A∆T
∆T
x(kT ) +
∆T
Definimos:
G (∆T ) = e A∆T ,
⎛ ∆T
⎞
H (∆T ) = ⎜ ∫ e Aλ u(kT )d λ ⎟ B .
⎝0
⎠
Obtemos, portanto
x(kT + ∆T ) = G (∆T )x(kT ) + H (∆T )u(kT ) .
Por sua vez, a resposta é:
y (kT + ∆T ) = Cx(kT + ∆T ) + Du(kT ) , ou seja,
y (kT + ∆T ) = CG (∆T )x(kT ) + [CH (∆T ) + D]u(kT ) .
∫e
Aλ
Bu(kT )d λ .
0
(5.78)
(5-79)
(5-80)
(5-81)
83
Exemplo 5-6: Seja o sistema discutido no Exemplo (5-5). Obter a eq. de estado
discretizada. Obter as expressões para T = 1 e ∆ T = 0,5 seg.
Solução: Vimos no exemplo 5-5 que
⎡ ⎛
e −2T − 1 ⎞ ⎤
⎡1 0,5(1 − e −2T ) ⎤
⎢ 0,5 ⎜ T +
⎟⎥
G (T ) = ⎢
2 ⎠ ⎥ . Então de (5-80) e (5-81).
⎥ e H (T ) = ⎢ ⎝
−2T
e
⎣0
⎦
⎢⎣ 0,5(1 − e −2T ) ⎥⎦
temos
⎡ x1 (kT + ∆T ) ⎤ ⎡1 0,5(1 − e
⎢ x (kT + ∆T ) ⎥ = ⎢
e −2 ∆T
⎣ 2
⎦ ⎣0
−2 ∆T
⎡ ⎛
e −2 ∆T − 1 ⎞ ⎤
) ⎤ ⎡ x1 (kT ) ⎤ ⎢0,5 ⎜ ∆T +
⎟⎥
+⎢ ⎝
2
⎥⎢
⎠ ⎥ u (kT ) ,
⎥
⎦ ⎣ x2 (kT ) ⎦ ⎢
−2 ∆T
) ⎥⎦
⎣ 0,5(1 − e
⎡ ⎛
e −2 ∆T − 1 ⎞ ⎤
∆
+
0,5
T
⎡1 0,5(1 − e −2 ∆T ) ⎤ ⎡ x1 (kT ) ⎤
⎢ ⎜
⎟⎥
+ [1 0] ⎢ ⎝
y (kT + ∆T ) = [1 0] ⎢
2
⎥⎢
⎠
⎥
−2 ∆T
⎥ u (kT ) .
x
(
kT
)
0
e
⎦
⎣
⎦⎣ 2
⎢⎣ 0,5(1 − e −2 ∆T ) ⎥⎦
Para T = 1 e ∆ T = 0,5 seg., temos
⎡ x1 (k + 0,5) ⎤ ⎡1 0,3161⎤ ⎡ x1 (k ) ⎤ ⎡ 0, 092 ⎤
⎢ x (k + 0,5) ⎥ = ⎢0 0,3679 ⎥ ⎢ x (k ) ⎥ + ⎢ 0,3161⎥ u (k ) ,
⎦⎣ 2 ⎦ ⎣
⎦
⎣ 2
⎦ ⎣
⎡ x (kT ) ⎤
y (k + 0,5) = [1 0,3161] ⎢ 1
⎥ + 0, 092u (k ) .
⎣ x2 (kT ) ⎦
Problemas e soluções
Problema 1: Achar a forma canônica controlável do seguinte sistema:
Y ( z)
z −1 + 5 z −2
.
=
U ( z ) 1 + 4 z −1 + 3z −2
Y ( z)
z +5
= 2
.
Solução: Multiplicando numerador e denominador por z 2 , temos
U ( z) z + 4z + 3
⎡ x (k + 1) ⎤ ⎡ 0 1 ⎤ ⎡ x1 (k ) ⎤ ⎡0 ⎤
Portanto, de (5-7)-(5-9), temos ⎢ 1
⎥=⎢
⎥ + ⎢ ⎥ u (k ) ,
⎥⎢
⎣ x2 (k + 1) ⎦ ⎣ −3 −4 ⎦ ⎣ x2 (k ) ⎦ ⎣1 ⎦
⎡ x (k ) ⎤
y (k ) = [5 1] ⎢ 1 ⎥ .
⎣ x2 (k ) ⎦
Problema A-5-6: Obter uma representação de estado do sistema da Figura 5-4. O
período de amostragem é T = 1 seg.
84
Solução: Primeiramente obtemos a transformada z da função de transferência em malha
aberta (ou de “feedforward”, ou do canal direto)
⎛ 1 − e− s
⎛
⎞ 0,3679( z + 0, 7181)
1 ⎞
1
−1
.
G( z) = Z ⎜
⎟ = (1 − z )Z ⎜ 2
⎟=
⎝ s ( s + 1) ⎠ ( z − 1)( z − 0,3679)
⎝ s s ( s + 1) ⎠
Expandindo esta em frações parciais, temos:
1
0, 6321
z −1
0, 6321z −1
−
=
−
G( z) =
. A Figura 5-5 abaixo mostra o
z − 1 z − 0,3679 1 − z −1 1 − 0,3679 z −1
diagrama de blocos do sistema. Escolhemos a saída de cada unidade de retardo como
variável de estado, conforme mostrado na figura.
Do diagrama de blocos acima temos:
zX 1 ( z ) = X 1 ( z ) − [ X 1 ( z ) − 0, 6321X 2 ( z )] + U ( z ) ,
zX 2 ( z ) = 0,3679 X 2 ( z ) − [ X 1 ( z ) − 0, 6321X 2 ( z )] + U ( z ) ,
Y ( z ) = X 1 ( z ) − 0, 6321X 2 ( z ) .
E destas eqs. temos imediatamente:
x1 (k + 1) = 0, 6321x2 (k ) + u (k ) ,
x2 (k + 1) = − x1 (k ) + x2 (k ) + u (k ) ,
y (k ) = x1 (k ) − 0, 6321x2 (k ) .
Ou, em forma vetorial,
⎡ x1 (k + 1) ⎤ ⎡ 0 0, 6321⎤ ⎡ x1 (k ) ⎤ ⎡1⎤
⎡ x (k ) ⎤
+ ⎢ ⎥ u (k ) , y (k ) = [1 −0, 6321] ⎢ 1 ⎥ .
⎢ x (k + 1) ⎥ = ⎢ −1
⎢
⎥
⎥
1 ⎦ ⎣ x2 (k ) ⎦ ⎣1⎦
⎣ 2
⎦ ⎣
⎣ x2 (k ) ⎦
Teorema A-5-9: Seja a eq. característica de uma matriz quadrada de dimensão n dada
por λ I − A = λ n + a1λ n −1 + ... + an −1λ + an = 0 . Demonstre que a matriz A satisfaz á eq.
característica, isto é, A n + a1A n −1 + ... + an −1A + an I = 0 .
(Este é o teorema de Cayley-Hamilton, que deve ter sido dado no curso anterior ou no
curso de Álgebra Linear).
Definição: (Acontece que a eq. característica não é necessariamente a equação de menor
grau satisfeita pela matriz A ). O polinômio de menor grau tal que A seja sua raiz é
denominado polinômio mínimo. Seja φ (t ) este polinômio:
85
φ (λ ) = λ m + a1λ m−1 + .... + am−1λ + am , com m ≤ n .
Então, por definição, φ ( A) = A m + a1A m −1 + .... + am −1A + am I = 0 .
Propriedade A-5-10: Seja d (λ ) o maior divisor comum de todos os elementos de
adj (λ I − A) e suponha que ele seja mônico, isto é, o coeficiente do seu termo de mais
adj (λ I − A)
.
alto grau é 1. Então o polinômio mínimo é dado por: φ (λ ) =
d (λ )
Propriedade A-5-11: Se os auto-valores de uma matriz são todos distintos, então o
polinômio mínimo coincide com o polinômio característico da matriz. Se algum autovalor tiver grau de multiplicidade maior que 1 e houver um bloco de Jordan associado a
ele, então o polinômio mínimo também coincide com o polinômio característico da
matriz. Se, porém não houver bloco de Jordan associado ao auto-valor com grau de
multiplicidade maior que 1, então o polinômio tem grau menor do que o polinômio
característico da matriz.
⎡1 0 ⎤
Assim, por exemplo, no primeiro caso, se A1 = ⎢
⎥ , então
⎣0 2⎦
λ −1 0
= (λ − 1)(λ − 2) = λ 2 − 3λ + 2 . É fácil verificar que A12 − 3A1 + 2I = 0 e é
λ −2
0
claro que não existem números a e b tais que aA1 + bI = 0 , portanto o polinômio
mínimo coincide com o polinômio característico da matriz.
⎡1 1⎤
2
Seja agora o caso de um bloco de Jordan: A 2 = ⎢
⎥ ∴ λ I − A 2 = (λ − 1) e é fácil
0
1
⎣
⎦
verificar que A 22 − 2 A 2 + I = 0 ; é fácil ver também que não existem números a e b tais
que aA1 + bI = 0 .
Finalmente seja o caso em que o auto-valor tem multiplicidade 2, mas não há bloco de
⎡1 0 ⎤
Jordan: A 3 = ⎢
⎥ ; o polinômio característico é igual ao do caso anterior, ou seja,
⎣0 1 ⎦
(λ − 1) 2 ; mas neste caso, com a = 1 e b = -1, temos aA 3 + bI = 0 , ou seja, o polinômio
mínimo é λ − 1 , com grau menor que o do polinômio característico.
Y ( s)
ω2
= 2
. Obtenha
U ( s) s + ω 2
representação de estado deste sistema de tempo contínuo, discretize-o e obtenha a
representação de estado de estado do sistema discretizado. Finalmente, obtenha também
a função de transferência de pulso.
Solução: Da definição do sistema temos imediatamente a eq. diferencial:
y + ω 2 y = ω 2u .
y
Definamos as variáveis de estado x1 = y, x2 = . Ou seja, em forma vetorial:
Problema A-5-14: Seja o seguinte sistema, um oscilador:
ω
⎡ x1 ⎤ ⎡ 0 ω ⎤ ⎡ x1 ⎤ ⎡ 0 ⎤
⎡ x1 ⎤
⎢ x ⎥ = ⎢ −ω 0 ⎥ ⎢ x ⎥ + ⎢ω ⎥ u , y = [1 0] ⎢ x ⎥ .
⎦⎣ 2⎦ ⎣ ⎦
⎣ 2⎦ ⎣
⎣ 2⎦
Em vista de A e B dados acima, de (5-73) e (5-74) temos:
G = e AT = L-1[( sI − A) −1 ]
86
ω ⎤
⎡ s
−1
2
2
2
⎛
⎞
⎢
sen cos ωT ⎤
⎡ s −ω ⎤
s + ω 2 ⎥ ⎡ cos ωT
-1 s + ω
⎟
= L-1 ⎜ ⎢
L
=
,
⎢
⎥=⎢
⎥
⎜ ⎣ω s ⎦ ⎟
−ω
cos ωT ⎦⎥
s ⎥ ⎣ − sen cos ωT
⎢
⎝
⎠
⎢⎣ s 2 + ω 2 s 2 + ω 2 ⎥⎦
T
⎛ T ⎡ cos ωT senωT ⎤ ⎞ ⎡ 0 ⎤ ⎡1 − cos ωT ⎤
H = ∫ e Aλ B = ⎜ ∫ ⎢
⎥ dλ ⎟ ⎢ ⎥ = ⎢
⎥.
0
0
⎝ ⎣ − senωT cos ωT ⎦ ⎠ ⎣ω ⎦ ⎣ senωT ⎦
Temos então
⎡ x1 ((k + 1)T ) ⎤ ⎡ cos ωT senωT ⎤ ⎡ x1 (kT ) ⎤ ⎡1 − cos ωT ⎤
⎢ x ((k + 1)T ) ⎥ = ⎢ − senωT cos ωT ⎥ ⎢ x (kT ) ⎥ + ⎢ senωT ⎥ u (kT ) ,
⎦⎣ 2
⎦
⎣ 2
⎦ ⎣
⎦ ⎣
(
)
⎡ x (kT ) ⎤
y (kT ) = [1 0] ⎢ 1
⎥.
⎣ x2 (kT ) ⎦
A função de transferência de pulso é:
−1
− senωT ⎤ ⎡1 − cos ωT ⎤
⎡ z − cos ωT
−1
F ( z ) = C( zI − G ) H + D = [1 0] ⎢
z − cos ωT ⎥⎦ ⎢⎣ senωT ⎥⎦
⎣ senωT
senωT ⎤ ⎡1 − cos ωT ⎤
⎡ z − cos ωT
1
= 2
[1 0] ⎢
z − cos ωT ⎥⎦ ⎢⎣ senωT ⎥⎦
z − 2 z cos ωT + 1
⎣ − senωT
(1 − cos ωT )( z + 1) (1 − cos ωT )(1 + z −1 ) z −1
=
.
= 2
z − 2 z cos ωT + 1
1 − 2 z −1 cos ωT + z −2
Seja observado que se obtém o mesmo resultado achando a transformada z da função de
transferência de tempo contínuo,
ω2
precedida de um segurador de ordem zero.
s2 + ω 2
Problema A-5-15: Seja o sistema da figura 5-8(a) abaixo, o qual tem polos complexos.
Este sistema é marginalmente estável, como sabemos do primeiro curso. A figura 5-8(b)
é uma versão discretizada do sistema de tempo contínuo. Este sistema também é
marginalmente estável. Suponha que o sistema seja submetido a um degrau unitário na
entrada. Mostre que o sistema discretizado pode ter oscilações ocultas quando o período
de amostragem T assume certo valor.
s2 1
s
= 2
, ou seja,
2
s +4 s s +4
y(t)=cos2t . Esta resposta é mostrada na figura 5-9(a) abaixo.
Solução: A resposta do sistema de tempo contínuo é Y ( s ) =
87
A função de transferência de pulso do sistema discretizado é:
⎛ 1 − e −Ts s 2 ⎞
Y ( z)
1 − z −1 cos 2T
⎛ s ⎞
−1
−1
. Portanto a
= Z⎜
=
−
=
−
(1
z
)Z
(1
z
)
⎟
⎜ 2
⎟
2
U ( z)
1 − 2 z −1 cos 2T + z −2
⎝s +4⎠
⎝ s s +4⎠
resposta do sistema ao degrau unitário será: (1 − z −1 )
1 − z −1 cos 2T
1
=
−1
−2
1 − 2 z cos 2T + z 1 − z −1
1 − z −1 cos 2T
.
1 − 2 z −1 cos 2T + z −2
É claro da tabela 2-1, entrada 15, que a resposta é y (kT ) = cos 2kT , com ω = 2.
1
⎛π ⎞
⎛π ⎞
Se, por exemplo, T = π seg., temos y (kT ) = y ⎜ k ⎟ = cos ⎜ k ⎟ , ou seja,
4
⎝2 ⎠
⎝2 ⎠
y (0) = 1; y (1) = 0; y (2) = −1; y (3) = 0; y (4) = 1;... Um gráfico destes resultados é
apresentado na figura 5-9(b) acima, no qual se vê que a resposta é oscilatória.
(1 − z −1 )(1 − z −1 ) 1
1
Se porém T = π , então temos Y ( z ) =
=
, e portanto a
−1
−2
−1
1− 2z − z
1− z
1 − z −1
resposta no domínio do tempo é y (k ) = 1 para k = 1; 2;3;... , mostrado na figura (c).
Esta resposta, ao contrário da anterior, pode induzir a idéia de que y (k ) se mantenha
constante, ou seja, com T = π temos oscilações ocultas.
88