Capítulo 3 - Métodos Determinísticos fileIdentiﬁcação de Sistemas Dinâmicos — Prof. Samir Martins Introdução I Métodos determinísticos !não dão tratamento especial

Identificação de Sistemas Dinâmicos — Prof. Samir Martins

Capítulo 3 - Métodos Determinísticos

Prof. Samir Martins

UFSJ-CEFET/MGPrograma de Pós Graduação em Engenharia Elétrica - PPGEL

São João del-Rei, 5 de setembro de 2018

1 / 61


Introdução

I Métodos determinísticos→ não dão tratamento especial aoruído presente nos dados, ainda que se aceite que os mesmosestejam contaminados;

I Métodos pouco imunes a ruídos;

I Requerem elevada relação sinal/ruído;

I Métodos estocásticos→ tratam o ruído presente nos dados(Capítulo 5);

2 / 61


Introdução





2 / 61


Introdução





2 / 61


Introdução





2 / 61


Alguns Casos Simples

I Um sistema de primeira ordem pode ser representado por:

H(s) =K

τs + 1, (1)

K e τ podem ser estimados por:

K =(y(∞)− y(0−)

A, (2)

y(τ) = 0,632(y(∞)− y(0−)) + y(0−), (3)

sendo A a amplitude do degrau aplicado.

3 / 61


Alguns Casos Simples

I Um sistema de primeira ordem pode ser representado por:

H(s) =K

τs + 1, (1)

K e τ podem ser estimados por:

K =(y(∞)− y(0−)

A, (2)

y(τ) = 0,632(y(∞)− y(0−)) + y(0−), (3)

sendo A a amplitude do degrau aplicado.

3 / 61


Casos Simples - Sistemas de Primeira Ordem

Figura 1: Resposta ao degrau (aplicado em t = 0s) de um sistema deprimeira ordem com constante de tempo τ .

4 / 61


Sistemas de Segunda Ordem Pouco Amortecidos

I Seja um possível modelo a ser identificado:

H(s) =Y (s)

U(s)=

Kω2n

s2 + 2ζωns + ω2n

(4)

I Para sistemas pouco amortecidos, o número de ciclos visíveispode ser aproximado por:

0.6ζ

(5)

I Exatidão do método depende da qualidade da aproximação√1− ζ2 ≈ 1→ ωn = ω (estimada diretamente do gráfico).

I K obtido como para sistemas de primeira ordem.

5 / 61




H(s) =Y (s)

U(s)=

Kω2n

s2 + 2ζωns + ω2n

(4)


0.6ζ

(5)



5 / 61




H(s) =Y (s)

U(s)=

Kω2n

s2 + 2ζωns + ω2n

(4)


0.6ζ

(5)



5 / 61




H(s) =Y (s)

U(s)=

Kω2n

s2 + 2ζωns + ω2n

(4)


0.6ζ

(5)



5 / 61


Estimação de ζ e ωn

Figura 2: Respostas ao degrau unitário. Sistemas com K = 1, ωn = 1,considerando ζ = 0,15 (–) e ζ = 0,4 (- -).

6 / 61



I Para a resposta menos amortecida:

ζ ≈ 0,64,5

= 0,13. (6)

I Temos aproximadamente quatro ciclos nos primeiros 25segundos→ "período médio":

T =254

ωn =2πT≈ 1,0. (7)

7 / 61



I Para a resposta menos amortecida:

ζ ≈ 0,64,5

= 0,13. (6)

I Temos aproximadamente quatro ciclos nos primeiros 25segundos→ "período médio":

T =254

ωn =2πT≈ 1,0. (7)

7 / 61



I Para a resposta mais amortecida:

ζ ≈ 0,6 (8)

I Apenas um ciclo é visível;

I Valor estimado 50% maior que o valor verdadeiro.

8 / 61




ζ ≈ 0,6 (8)



8 / 61




ζ ≈ 0,6 (8)



8 / 61


O Método de Sundaresan

I Pressupõe que o sistema possa ser aproximadosatisfatoriamente por:

H(s) =e−τd s

(τ1s + 1)(τ2s + 1)(9)

ou ainda:

H(s) =e−τd sω2

n

s2 + 2ζωns + ω2n

(10)

9 / 61



I Objetivo do método: determinação de parâmetros τd , τ1 e τ2 ouτd , ωn e ζ;

I Utiliza, para tal, a resposta do sistema ao degrau;

I Ganho posteriormente ajustado.

10 / 61






10 / 61






10 / 61


O método de Sundaresan

Figura 3: Resposta ao degrau típica de um sistema sobreamortecido comtempo morto.

11 / 61



I Para este método, tem-se:

m1 =

∫ ∞0

(1− y(t))dt . (11)

I H(s) e m1 estão relacionados por:

m1 = −dH(s)

ds|s=0 = τdτ1 + τ2. (12)

I A resposta ao degrau é dada por:

y(t) =

[1− τ1

τ1 − τ2e−

t−τdτ1 +

τ2

τ1 − τ2e−

t−τdτ2

]u(t − τd ), (13)

12 / 61




m1 =

∫ ∞0

(1− y(t))dt . (11)


m1 = −dH(s)

ds|s=0 = τdτ1 + τ2. (12)


y(t) =

[1− τ1

τ1 − τ2e−

t−τdτ1 +

τ2

τ1 − τ2e−

t−τdτ2

]u(t − τd ), (13)

12 / 61




m1 =

∫ ∞0

(1− y(t))dt . (11)


m1 = −dH(s)

ds|s=0 = τdτ1 + τ2. (12)


y(t) =

[1− τ1

τ1 − τ2e−

t−τdτ1 +

τ2

τ1 − τ2e−

t−τdτ2

]u(t − τd ), (13)

12 / 61


O Método de SundaresanI Igualando a segunda derivada da resposta ao degrau a zero

(ponto de inflexão):

ti = τd + αln η, (14)

I sendo η = τ1/τ2 e α = τ1τ2/(τ1 − τ2).I A inclinação da reta tangente ao ponto de inflexão pode ser

obtida por:

Mi =η1/(1−η)

α(η − 1). (15)

I tm pode ser expresso por:

tm = τd + α

[ln η +

η2 − 1η

]. (16)

13 / 61




ti = τd + αln η, (14)


obtida por:

Mi =η1/(1−η)

α(η − 1). (15)


tm = τd + α

[ln η +

η2 − 1η

]. (16)

13 / 61




ti = τd + αln η, (14)


obtida por:

Mi =η1/(1−η)

α(η − 1). (15)


tm = τd + α

[ln η +

η2 − 1η

]. (16)

13 / 61



I Combinando as equações anteriormente apresentadas,chega-se a:

Mi (tm −m1) =η1/(1−η)

η − 1ln η, (17)

I a qual pode ser escrita como:

λ = χe−χ, (18)

I sendo λ = (tm −m1)Mi e ainda

χ =ln ηη − 1

. (19)

14 / 61




Mi (tm −m1) =η1/(1−η)

η − 1ln η, (17)


λ = χe−χ, (18)


χ =ln ηη − 1

. (19)

14 / 61




Mi (tm −m1) =η1/(1−η)

η − 1ln η, (17)


λ = χe−χ, (18)


χ =ln ηη − 1

. (19)

14 / 61



Figura 4: Relação entre η e λ.

15 / 61


O Método de Sundaresan - Resumo

I determinar o ganho em regime permanente:

I determinar a área m1:

I determinar a inclinação da reta tangente ao ponto de inflexão dey(t). Este valor é Mi :

I determinar tm:

I determinar λ a partir de λ = (tm −m1)Mi :

16 / 61






I determinar tm:


16 / 61






I determinar tm:


16 / 61






I determinar tm:


16 / 61






I determinar tm:


16 / 61



I determinar η a partir da última Figura apresentada:

I determinar τ1, τ2 e τd por

τ1 =η

η1−η

Mi(20)

τ2 =η

11−η

Mi(21)

τd = m1 − τ1 − τ2. (22)

17 / 61



I determinar η a partir da última Figura apresentada:

I determinar τ1, τ2 e τd por

τ1 =η

η1−η

Mi(20)

τ2 =η

11−η

Mi(21)

τd = m1 − τ1 − τ2. (22)

17 / 61


O Método de Sundaresan - Caso Subamortecido

Figura 5: Resposta ao degrau típica de um sistema subamortecido comtempo morto.

18 / 61




y(t) = u(t − τd )

{1− ζωnet−τd

[ζ√

1− ζ2sen

(ωn(t − τd )

√1− ζ2

)+

+cos(ωn(t − τd )

√1− ζ2

)]}(23)

I Seguindo procedimentos análogos aos apresentados, obtém-se:

19 / 61



λ = (tm −m1)M1 =cos−1ζ√

1− ζ2e− ζcos−1(ζ)√

1−ζ2 , (24)

ωn =cos−1ζ√

1− ζ2

1tm −m1

, (25)

τd = m1 −2ζωn. (26)

20 / 61



λ = (tm −m1)M1 =cos−1ζ√

1− ζ2e− ζcos−1(ζ)√

1−ζ2 , (24)

ωn =cos−1ζ√

1− ζ2

1tm −m1

, (25)

τd = m1 −2ζωn. (26)

20 / 61



λ = (tm −m1)M1 =cos−1ζ√

1− ζ2e− ζcos−1(ζ)√

1−ζ2 , (24)

ωn =cos−1ζ√

1− ζ2

1tm −m1

, (25)

τd = m1 −2ζωn. (26)

20 / 61



I ζ Pode ser obtido a partir de λ por:

Figura 6: Relação entre ζ e λ

21 / 61




I determinar a área sombreada na resposta ao degrau, a áreaacima do valor de y(∞) deve ser subtraída. Este valor é m1:

I determinar a inclinação da tangente no ponto de inflexão dey(t). Este valor é Mi :

I determinar tm, que é a interseção da tangente com o valor emregime permanente de y(t):

I determinar λ

22 / 61







I determinar λ

22 / 61







I determinar λ

22 / 61







I determinar λ

22 / 61







I determinar λ

22 / 61



I determinar ζ graficamente:

I determinar ωn e τd .

23 / 61



I determinar ζ graficamente:

I determinar ωn e τd .

23 / 61


O Método de Sundaresan - Exemplo Experimental

I Identificação de Malha de controle de combustíveis;

I Sistema de controle que injeta até quatro combustíveis nafornalha de uma caldeira industrial.

I Devido à falta de medições e dados de fabricação, não se podeobter modelos pela física do processo;

24 / 61






24 / 61






24 / 61


A malha de combustíveis

Figura 7: Malha de controle de combustíveis de uma planta industrial real. osblocos "FIC"são controladores de vazão, ao passo que os blocos "subtrator","somador", "max/min"e √ representam circuitos com funções aritméticas.

25 / 61



Aspectos considerados na execução dos testes:

I Amplitude do sinal de excitação;I Ponto de operação;I Sentido do sinal de excitação.

Resultados obtidos pela aplicação do método de Sundaresan:

H1(s) =1,338e−1,9s

(3,406s + 1)(1,053s + 1)(27)

H2(s) =0,0182e−4,7s

s2 + 2(0,4)(0,228)s + 0,052(28)

H3(s) =0,189e−2s

s2 + 2(0,55)(0,346)s + 0,120(29)

26 / 61






H1(s) =1,338e−1,9s

(3,406s + 1)(1,053s + 1)(27)

H2(s) =0,0182e−4,7s

s2 + 2(0,4)(0,228)s + 0,052(28)

H3(s) =0,189e−2s

s2 + 2(0,55)(0,346)s + 0,120(29)

26 / 61






H1(s) =1,338e−1,9s

(3,406s + 1)(1,053s + 1)(27)

H2(s) =0,0182e−4,7s

s2 + 2(0,4)(0,228)s + 0,052(28)

H3(s) =0,189e−2s

s2 + 2(0,55)(0,346)s + 0,120(29)

26 / 61






H1(s) =1,338e−1,9s

(3,406s + 1)(1,053s + 1)(27)

H2(s) =0,0182e−4,7s

s2 + 2(0,4)(0,228)s + 0,052(28)

H3(s) =0,189e−2s

s2 + 2(0,55)(0,346)s + 0,120(29)

26 / 61


Figura 8: Entrada (acima) e respostas (abaixo) para a planta/modelo.Medição feita em indústria (–) e (- -) simulação utilizando H1(s).

27 / 61



28 / 61



29 / 61


Identificação em Malha Fechada

Método a seguir descrito apresentado originalmente por Yawana eSeborg (1982).

Figura 11: Sistema de controle em malha fechada.

Assume-se controlador puramente proporcional e:

H(s) =Ke−τd s

τs + 1(30)

30 / 61



Figura 12: Resposta ao degrau típica necessária para o método proposto porYawana e Seborg (1982).

31 / 61



I A função de transferência em malha fechada é dada por:

Y (s)

R(s)=

Kf e−τd s

τs + [Kf e−τd s + 1](31)

I Utilizando a seguinte aproximação de Padé:

e−τd s =1− 0,5τds1 + 0,5τds

(32)

I chega-se a:

32 / 61



I A função de transferência em malha fechada é dada por:

Y (s)

R(s)=

Kf e−τd s

τs + [Kf e−τd s + 1](31)

I Utilizando a seguinte aproximação de Padé:

e−τd s =1− 0,5τds1 + 0,5τds

(32)

I chega-se a:

32 / 61



Y (s)

R(s)=

K (1− 0,5τds)

τ2s2 + 2ζτs + 1(33)

sendo Kc o ganho do controlador, Kf = KcK e

K =Kf

Kf + 1(34)

τ =

[τdτ

2(Kf + 1)

]0,5

(35)

ζ =τ + 0,5τd (1− Kf )√

2τdτ(Kf + 1)(36)

33 / 61



Y (s)

R(s)=

K (1− 0,5τds)

τ2s2 + 2ζτs + 1(33)

sendo Kc o ganho do controlador, Kf = KcK e

K =Kf

Kf + 1(34)

τ =

[τdτ

2(Kf + 1)

]0,5

(35)

ζ =τ + 0,5τd (1− Kf )√

2τdτ(Kf + 1)(36)

33 / 61



I Os valores dos parâmetros do processo podem ser obtidos dasseguintes expressões:

K =y∞

Kc(A− y∞)(37)

τ =∆tπ

[ζ√

Kf + 1 +√ζ2(Kf + 1) + Kf

]√(1− ζ2)(Kf + 1), (38)

τd =∆t√

(1− ζ2)(Kf + 1)

π[ζ√

Kf + 1 +√ζ2(Kf + 1) + Kf

] , (39)

34 / 61



I Sendo que ζ pode ser avaliado de duas formas diferentes:

ζ =−ln

[y∞−ymyp1−y∞

]√π2 +

(ln[

y∞−ymyp1−y∞

])2(40)

ζ =−ln

[yp2−yinftyyp1−y∞

]√

4π2 +(

ln[

yp2−yinftyyp1−y∞

]) (41)

35 / 61



I O valor em regime permanente pode ainda ser estimado por:

y∞ ≈yp2yp1 − y2

m

yp2 + yp1 − 2ym. (42)

36 / 61


Identificação em Malha Fechada - Processo Real

I Objetivo→ obter modelos em malha aberta a partir de dadosreais coletados de uma planta industrial operando em malhafechada;

I Malha de corrente de um separador magnético.

I A partir dos dados medidos, obteve-se: yp1 = 33, yp2 = 28,ym = 16, ∆t = 7, y∞ = 22,8, com Kc = 2.

37 / 61






37 / 61






37 / 61



I O modelo em malha aberta obtido foi:

H1(s) =2,1923e−5,2853s

12,3681s + 1. (43)

I Ajustes realizados: aumento de K em 70% e aumento de ζpor 2,4:

H2(s) =3,7269e−3,4108s

23,2170s + 1(44)

38 / 61



Figura 13: Respostas em malha fechada em resposta a um degrau naentrada. (a) medido, (b) modelo obtido sem ajustes, (c) modelo obtido, apósajustes.

39 / 61


Identificação Usando Convolução

I Seja o somatório de convolução:

y(k) =∞∑j=0

h(j)u(k − j). (45)

Tomando-se as medições u(k) e y(k), pode-se escrever:

y(0) = h(0)u(0) + h(1)u(−1) + h(2)u(−2) + · · · (46)y(1) = h(0)u(1) + h(1)u(0) + h(2)u(−1) + · · · (47)y(2) = h(0)u(2) + h(1)u(1) + h(2)u(0) + · · · (48)

......

...

40 / 61



I Que também pode ser expresso em forma matricial:

y(0)y(1)y(2)

...

=

u(0) u(−1) u(−2) · · ·u(1) u(0) u(−1) · · ·u(2) u(1) u(0) · · ·

......

...

h(0)h(1)h(2)

...

, (49)

y = Uh (50)

I Se u(k) = 0 ∀k 6= 0, pode-se escrever:

y = u(0)I h, (51)h = y/u(0). (52)

41 / 61



I Se o sinal de saída estiver contaminado com ruído de medição(y(k) = y i (k) + e(k)):

h =yi + eu(0)

h =yi

u(0)+

eu(0)

(53)

I A segunda parcela da equação acima somente é eliminada se: i)u(0)→∞ ou; ii) e(k) = 0,∀ k .

42 / 61


Estimação Determinística da Resposta ao Impulso

I Seja a função de transferência:

H(z) =Y (z)

U(z)=

z + 0,5z2 − 1,5z + 0,7

(54)

I H(z) simulada para uma entrada pseudo-aleatória;

I Resposta ao impulso obtida como anteriormente apresentado;

I Duas abordagens apresentadas: i) sem ruído e; ii) com ruído debaixíssima variância (σ2

e = 1× 10−5).

I h obtida de maneira determinística (h = U−1y);

43 / 61




H(z) =Y (z)

U(z)=

z + 0,5z2 − 1,5z + 0,7

(54)




e = 1× 10−5).


43 / 61




H(z) =Y (z)

U(z)=

z + 0,5z2 − 1,5z + 0,7

(54)




e = 1× 10−5).


43 / 61




H(z) =Y (z)

U(z)=

z + 0,5z2 − 1,5z + 0,7

(54)




e = 1× 10−5).


43 / 61




H(z) =Y (z)

U(z)=

z + 0,5z2 − 1,5z + 0,7

(54)




e = 1× 10−5).


43 / 61



Figura 14: Sinais de (a) entrada e (b) saída de H(z).

44 / 61



Figura 15: Resposta ao impulso de H(z) (o), juntamente com as respostasao impulso calculadas para o caso sem ruído (a) e com ruído (b).

45 / 61


Identificação no Domínio da Frequência

I Resposta em regime permanente de um sistema linear excitadopor uma senóide, também será uma senóide de mesmafrequência.

I Amplitude e fase da resposta dependem do sistema, e podemser alteradas

I Amplitude→ multiplicada por |H(jω)|

I Fase→ atrasada de φ radianos, sendo H(jω) = |H(jω)|ejφ.

46 / 61







46 / 61







46 / 61







46 / 61



I Procedimento possível para estimar H(jω)→ excitar o sistemacom senoides de diferentes frequências.

I A partir daí, obtém-se o ganho e fase.

I Problema→ nem sempre se pode aplicar sinais senoidais aoprocesso a ser identificado

I Problema→ teste de longa duração.

47 / 61







47 / 61







47 / 61







47 / 61



I Aplicando-se a transformada de Fourier à integral deconvolução:

H(jω) =Y (jω)

U(jω). (55)

I Sinais "abertos"(resposta a um degrau - sistema estável deprimeira ordem)→ não são absolutamente integráveis(condições de Dirichlet);

I Diferenciação destes sinais: forma alternativa para se obtersinais "fechados", não alterando a resposta em frequência:

48 / 61




H(jω) =Y (jω)

U(jω). (55)



48 / 61




H(jω) =Y (jω)

U(jω). (55)



48 / 61



H(jω) =F{y(t)}F{u(t)}

=jωY (jω)

jωU(jω). (56)

I Para sinais de tempo discreto:

y∗(k) = y(k − 1)− y(k) (57)u∗(k) = u(k − 1)− u(k) (58)

49 / 61



H(jω) =F{y(t)}F{u(t)}

=jωY (jω)

jωU(jω). (56)

I Para sinais de tempo discreto:

y∗(k) = y(k − 1)− y(k) (57)u∗(k) = u(k − 1)− u(k) (58)

49 / 61



I E a estimativa da resposta em frequência é:

H(ejω) =Y ∗(ejω)

U∗(ejω), (59)

I sendo que o argumento ejω indica a transformada de Fourier desinais discretos.

I A transformada Discreta de Fourier para um sinal y(k), com Namostras, é definida como:

Y (ejω) =1√N

N∑k=1

y(k)e−jωk , ω =2πkN

, k = 1, . . . ,N. (60)

50 / 61


Identificaçao no Domínio da Frequência

I Se y(t) é composto de uma parcela ideal (sem ruído - y i (t)) e deruído de medição e(t), tem-se:

H(jω) =Y i (jω)

U(jω)+

E(jω)

U(jω). (61)

I Para que e(t) não afete a estimativa da resposta emfrequência→ e(t) = 0,∀ t .

51 / 61


Estimação da Resposta em Frequência - Um CasoSimulado

I Considere que:

H(z) =Y (z)

U(z)=

z + 0,5z2 − 1,5z + 0,7

(62)

foi excitada por um sinal pseudo-aleatório.

I O número complexo H(ejωi ) = α + jβ foi obtido dividindo-seY (ejωi ) por U(ejωi ):

|H(ejωi )| =√α2 + β2 (63)

fase[H(ejωi )

]= tan−1

[β

α

]. (64)

52 / 61


Estimação da Resposta em Frequência - Um CasoSimulado

Figura 16: (a) ganho e (b) fase das respostas em frequência do sistemaoriginal (traço contínuo) e estimada a partir de y(k) e u(k) (tracejado).Procedimento sensível a problemas numéricos, especialmente em elevadasfrequências.

53 / 61


Estimação da Resposta em Frequência - Um CasoReal

I Dados coletados da planta piloto de bombeamento de água;

I A variação do sinal de saída é relativamente pequenacomparada à faixa total de operação (menos de 0,1V, em umfundo de escala de 5V);

I O teste dinamico não começa com a planta "desligada";

I O tempo de amostragem utilizado foi de Ts = 1,044s.

54 / 61







54 / 61







54 / 61







54 / 61



Figura 17: (a) pulso de entrada, u(k), e (b) resposta de vazão da planta,y(k).

55 / 61



Figura 18: Módulos das transformadas de Fourier dos sinais de entrada (a) esaída (b).

56 / 61



Foram utilizados apenas os dados correspondentes até a frequênciade ω = 0,015rad/s.

Figura 19: Resposta em frequência estimada. Representação gráfica deH(ejω). (a) ganho e (b) fase.

57 / 61


Comentários Finais

Figura 20: Resposta ao degrau de um sistema de segunda ordemsubamortecido, altamente contaminada com ruído.

58 / 61


Comentários Finais

I A partir da Figura anterior, torna-se impossível aplicar métodosdeterminísticos para obtenção de modelos.

I Considerando o exemplo do ajuste de reta, com ruído,tomando-se dois pontos, chega-se a:

[θ1θ2

]=

[2,872× 10−1

−1,091× 10−3

], (65)

sendo y = θ1 + θ2x .

59 / 61


Comentários Finais



[θ1θ2

]=

[2,872× 10−1

−1,091× 10−3

], (65)


59 / 61


Comentários Finais



[θ1θ2

]=

[2,872× 10−1

−1,091× 10−3

], (65)


59 / 61


Comentários Finais

Figura 21: Dados medidos com ruído estão marcados com cruzes, e o traçoé a reta estimada - y = 2,872× 10−1 − 1,091× 10−3x

60 / 61


Agradecimentos

MUITO OBRIGADO!

Prof. Samir MartinsDEPEL/UFSJ

[email protected]

61 / 61

Documents

Capítulo 3 - Métodos Determinísticos fileIdentiﬁcação de Sistemas Dinâmicos — Prof. Samir Martins Introdução I Métodos determinísticos !não dão tratamento especial