Identificação de sistemas

Prévia do material em texto

Universidade Estadual do Oeste do Paraná
UNIOESTE/Campus de Foz do Iguaçu
Centro de Engenharias e Ciências Exatas - CECE
Identificação de Sistemas
Notas de Aulas
Prof. Dr. Samuel da Silva
Foz do Iguaçu, 2009.
Prefácio
Este texto apresenta as notas de aulas da disciplina de Identificação de Sis-
temas oferecida por mim no Centro de Engenharias e Ciências Exatas da
Universidade Estadual do Oeste do Paraná. Uma vez que existe uma escas-
sez de livros sobre este tema na biblioteca da UNIOESTE, eu espero que
este texto sirva como um guia simples, conciso e com boa qualidade gráfica1
para auxiliar no estudo por parte dos alunos dos cursos de graduação em
engenharias elétrica e mecânica da UNIOESTE. Deve ficar claro que esta
apostila não tem como objetivo substituir os excelentes livros na área, como
por exemplo [6] , [1] ou [4]. Os estudantes são encorajados a se aventurar em
busca de mais informações nestes tópicos. Espero contar com o apoio dos
alunos e demais colaboradores para melhorar este texto gradualmente, sendo
assim, sugestões, correções e comentários são muito bem vindos. Boa leitura
e estudo!
Prof. Dr. Samuel da Silva
fevereiro de 2009.
1O texto foi redigido com o LATEX2ε.
2
Sumário
Lista de Figuras 5
1 Introdução 8
1.1 Exemplos de aplicação . . . . . . . . . . . . . . . . . . . . . . 9
1.2 Histórico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Etapas de um problema de identificação . . . . . . . . . . . . 10
2 Introdução à Probabilidade 12
2.1 Espaço amostral . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.1 Álgebra de conjuntos . . . . . . . . . . . . . . . . . . . 13
2.2 Probabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.1 Probabilidade condicional . . . . . . . . . . . . . . . . 14
2.2.2 Independência . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.3 Regra de Bayes . . . . . . . . . . . . . . . . . . . . . . 15
2.2.4 Probabilidade total . . . . . . . . . . . . . . . . . . . . 15
2.3 Variável aleatória . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Função distribuição de probabilidade . . . . . . . . . . . . . . 16
2.5 Função densidade de probabilidade . . . . . . . . . . . . . . . 16
2.6 Variáveis aleatórias Gaussianas . . . . . . . . . . . . . . . . . 17
2.7 Momentos estatísticos . . . . . . . . . . . . . . . . . . . . . . 17
2.8 Correlação e covariância . . . . . . . . . . . . . . . . . . . . . 18
2.9 Processo estocástico . . . . . . . . . . . . . . . . . . . . . . . . 19
2.10 Processo estacionário . . . . . . . . . . . . . . . . . . . . . . . 21
2.11 Processo ergódico . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.12 Ruídos brancos e coloridos . . . . . . . . . . . . . . . . . . . . 22
2.12.1 Ruído branco . . . . . . . . . . . . . . . . . . . . . . . 22
2.12.2 Ruído colorido . . . . . . . . . . . . . . . . . . . . . . . 23
2.13 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3
3 Representações Matemáticas Lineares 24
3.1 Função de transferência . . . . . . . . . . . . . . . . . . . . . . 24
3.1.1 Caso contínuo - domínio s . . . . . . . . . . . . . . . . 25
3.1.2 Caso discreto - domínio z . . . . . . . . . . . . . . . . 27
3.2 Função de resposta ao impulso . . . . . . . . . . . . . . . . . . 28
3.2.1 Caso contínuo - Integral de convolução . . . . . . . . . 28
3.2.2 Caso discreto - Somatório de convolução . . . . . . . . 29
3.3 Realização no espaço de estados . . . . . . . . . . . . . . . . . 29
3.4 Representações matemáticas lineares de tempo discreto . . . . 32
3.4.1 Modelo FIR . . . . . . . . . . . . . . . . . . . . . . . . 34
3.4.2 Modelo auto-regressivo . . . . . . . . . . . . . . . . . . 35
3.4.3 Modelo auto-regressivo com entradas externas . . . . . 35
3.4.4 Modelo auto-regressivo com média móvel e entradas
externas . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.4.5 Modelo de erro na saída . . . . . . . . . . . . . . . . . 37
3.4.6 Modelo de Box-Jenkins . . . . . . . . . . . . . . . . . . 37
3.5 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4 Métodos Não-Paramétricos 40
4.1 Métodos determinísticos . . . . . . . . . . . . . . . . . . . . . 40
4.1.1 Sistema de 1.o ordem . . . . . . . . . . . . . . . . . . . 41
4.1.2 Sistema de 2.o ordem . . . . . . . . . . . . . . . . . . . 41
4.1.3 Método de Sundaresan . . . . . . . . . . . . . . . . . . 47
4.1.4 Identificação usando convolução . . . . . . . . . . . . . 49
4.1.5 Identificação da função de resposta em frequência . . . 52
4.2 Métodos estocásticos . . . . . . . . . . . . . . . . . . . . . . . 54
4.2.1 Método das correlações . . . . . . . . . . . . . . . . . . 54
4.2.2 Densidade espectral de potência . . . . . . . . . . . . . 57
4.2.3 Janelas . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2.4 Periodograma ponderado . . . . . . . . . . . . . . . . . 64
4.2.5 Periodograma de Welch . . . . . . . . . . . . . . . . . 66
4.2.6 Estimadores H1, H2 e Hv . . . . . . . . . . . . . . . . . 67
4.2.7 Função de coerência . . . . . . . . . . . . . . . . . . . 68
4.3 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5 Métodos Paramétricos 76
5.1 Estimador de mínimos quadrados (MQ) . . . . . . . . . . . . . 78
5.2 Propriedades do estimador MQ . . . . . . . . . . . . . . . . . 87
5.2.1 Ortogonalidade . . . . . . . . . . . . . . . . . . . . . . 87
5.2.2 Polarização . . . . . . . . . . . . . . . . . . . . . . . . 88
5.2.3 Variância . . . . . . . . . . . . . . . . . . . . . . . . . 90
4
5.3 Estimador de Yule-Walker . . . . . . . . . . . . . . . . . . . . 91
5.4 Estimação de parâmetros de modelos ARX . . . . . . . . . . . 93
5.5 Estimador estendido MQ . . . . . . . . . . . . . . . . . . . . . 94
5.6 Método da variável instrumental (VI) . . . . . . . . . . . . . . 96
5.7 Estimador MQ generalizado . . . . . . . . . . . . . . . . . . . 98
5.8 Forma recursiva do estimador MQ . . . . . . . . . . . . . . . . 100
5.9 Métodos de subespaço . . . . . . . . . . . . . . . . . . . . . . 102
5.9.1 Parâmetros de Markov de um sistema . . . . . . . . . . 102
5.9.2 Algoritmo de realização de autosistemas (ERA) . . . . 104
5.9.3 Método da exponencial complexa . . . . . . . . . . . . 107
5.10 Comentários finais . . . . . . . . . . . . . . . . . . . . . . . . 110
5.11 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6 Determinação de Estruturas, Verificação e Validação de Mo-
delos 114
6.1 Determinação de ordem de modelos . . . . . . . . . . . . . . . 115
6.2 Validação de modelos . . . . . . . . . . . . . . . . . . . . . . . 116
Referências Bibliográficas 118
5
Lista de Figuras
2.1 Realizações de sinais medidos em um processo estocástico. . . 20
2.2 Exemplo de um processo estacionário. . . . . . . . . . . . . . . 21
2.3 Distribuição de partes de um processo estacionário. . . . . . . 22
3.1 Representação de um sistema dinâmico linear, invariante com
o tempo e estacionário com ruído colorido adicionado a saída. 32
4.1 Resposta ao degrau unitário (A = 1) de um sistema de pri-
meira ordem. . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.2 Exemplo de resposta ao degrau unitário para um sistema com
um grau de liberdade. . . . . . . . . . . . . . . . . . . . . . . 43
4.3 Resposta de sistema subamortecido evidenciando amplitudes
sucessivas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4 Resposta livre do sistema. . . . . . . . . . . . . . . . . . . . . 46
4.5 Espectro de um sinal senoidal continuo. . . . . . . . . . . . . . 61
4.6 Espectro de um sinal senoidal amostrado. . . . . . . . . . . . . 61
4.7 Aplicação de janelas para reduzir o efeito de leakage. . . . . . 62
4.8 Característica espectrais de janelas. . . . . . . . . . . . . . . . 62
4.9 Efeito do janelamento no espectro de amplitude do sinal. . . . 64
4.10 Efeitodo janelamento na energia do sinal. . . . . . . . . . . . 65
4.11 Rotina no matlab para cálculo da densidade espectral de po-
tência usando o periodograma ponderado por uma janela Han-
ning. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.12 Densidade espectral de potência do sinal. . . . . . . . . . . . . 68
4.13 Sinais de entrada e saída do sistema . . . . . . . . . . . . . . . 69
4.14 Estimativa da FRF e da função de coerência. . . . . . . . . . . 70
4.15 Estimativa da FRF e da função de coerência. . . . . . . . . . . 70
4.16 Resposta ao degrau unitário. . . . . . . . . . . . . . . . . . . . 71
4.17 Resposta ao impulso unitário. . . . . . . . . . . . . . . . . . . 71
4.18 Resposta ao impulso unitário. . . . . . . . . . . . . . . . . . . 72
4.19 Resposta ao degrau unitário. . . . . . . . . . . . . . . . . . . . 72
4.20 Resposta ao degrau unitário. . . . . . . . . . . . . . . . . . . . 73
6
4.21 Resposta ao degrau y1. . . . . . . . . . . . . . . . . . . . . . . 73
4.22 Resposta ao degrau y2. . . . . . . . . . . . . . . . . . . . . . . 74
4.23 Resposta ao degrau do sinal de pressão de um tubo. . . . . . . 74
5.1 Propriedade de ortogonalidade do estimador MQ. . . . . . . . 88
7
Capítulo 1
Introdução
Nos dias atuais a competitividade faz com que as indústrias necessitem
cada vez mais de análise e controle de seus processos, sendo assim, existe
uma busca crescente pela obtenção de modelos matemáticos que representem
processos físicos. Três abordagens podem ser usadas para se identificar um
modelo de um processo:
Caixa Branca: consiste da análise físico matemática de um processo uti-
lizando o conhecimento de leis físicas conhecidas, no geral envolvendo
equações de balanço e relações constitutivas de propriedades. Exem-
plos de métodos caixa branca de modelagem incluem o uso de métodos
dos elementos finitos, elementos de contorno, volumes finitos, etc. Esta
abordagem pode se demonstrar extremamente complexa e cara em al-
gumas situações e é mais usada em projetos quando se desconhece
totalmente qual é a planta que se tem em mãos.
Caixa Preta: consiste em obter um modelo matemático de um sistema já
existente e em operação utilizando para isto medidas de entradas, saí-
das e/ou estados do sistema. Pode-se assumir pouco conhecimento do
sistema.
Caixa Cinza: método mais recente onde se realiza uma união entre técnicas
puramente analíticas de modelagem (caixa branca) com ferramentas
baseadas em dados experimentais das plantas (caixa preta). Depen-
dendo do tipo de aplicação, o método pode ser "cinza escuro"ou "cinza
claro". Esta categoria corresponde ao clássico problema de ajuste de
modelos, ou solução do problema inverso, usando dados experimentais
para correção de parâmetros físicos de modelos.
O foco deste curso será em métodos de identificação do tipo caixa preta.
Um problema de identificação nesta classe pode compreender dois tipos dife-
8
rentes de abordagem. A primeira é a estimação de parâmetros, que consiste
na determinação de grandezas físicas não-observáveis a partir de grandezas
mensuráveis de entrada e/ou saídas. O resultado da estimação é um modelo
linear ou não-linear, determínistico ou estocástico, contínuo ou discreto, es-
tático ou dinâmico, etc. Outra abordagem se refere a filtragem de dados que
é a determinação de estados a partir de dados de entrada/saída. Em um
problema de filtragem linear, se conhece o modelo que pode ser representado
por uma função de transferência ou realização no espaço de estados. Um
problema de filtragem não-linear já envolve o desconhecimento deste mo-
delo. Exemplos são o filtro de Kalman e outros observadores de estado. Este
curso será concentrado em estimação de parâmetros em especial de sistemas
lineares e monovariáveis via dados de entrada e saída.
A estimação pode ser realizada de duas maneiras: offline ou em bate-
lada, quando se adquire todos os dados da planta e depois se processa para
identificação do modelo, ou online, usando algum método recursivo onde o
modelo é atualizado de tempos em tempos.
1.1 Exemplos de aplicação
O grande interesse em aplicações de identificação de sistemas (ID) se deve
a enorme facilidade de se adquirir dados de processos. O avanço tecnológico
dos últimos anos permitiu que fossem acessíveis placas de aquisição de dados
com inúmeros canais, altas taxas de amostragem, etc. tudo a um preço baixo.
Este monte de informação precisa ser processado de alguma forma e neste
ponto é que a utilização de ID tem se mostrado uma ferramenta valiosa para
engenheiros.
As aplicações de ID são as mais variadas possíveis se destacando:
• Controle de processos, especialmente envolvendo ferramentas sofisticadas
de controle preditivo para ajuste de plantas complexas.
• Análise dinâmica de plantas existentes.
• Detecção, diagnóstico e prognóstico de falhas.
• Processamento de sinais.
• Sistemas biológicos.
• Economia.
9
entre inúmeras outras aplicações em todas as áreas do conhecimento.
Uma correta utilização de métodos de ID pode se transformar em enorme lu-
cro, economia, segurança, etc. sendo estes fatores motivadores para o estudo
e aprendizado de ID visando aplicações em ambiente profissional.
1.2 Histórico
O desenvolvimento de ID pode ser dividido de forma didática da seguinte
forma, visando facilitar o seu estudo:
até 1960: desenvolvimento de toda a base estatística envolvida em métodos
ID, em especial com as contribuições de Gauss envolvendo o método
dos mínimos quadrados (LSM).
1960-1970: neste período foram desenvolvidas praticamente todas as prin-
cipais ferramentas de ID que serão estudadas durante este curso.
1970-1985: nesta fase foi realizada uma consolidação das principais fer-
ramentas deste curso com intensa aplicação em indústrias diversas.
Destaca-se o aparecimento de softwares dedicados a ID e desenvolvi-
mento de placas de aquisição sofisticadas.
1985 até hoje: utilização de métodos frequenciais e construção de modelos
relevantes para aplicações de controle robusto. Também é destaque o
desenvolvimento de métodos ID para obtenção de modelos em malha
fechada.
1990 até hoje: observa-se uma popularização de técnicas de ID conside-
rando que as plantas são não-lineares com aplicações bem sucedidas de
ferramentas como redes neurais artificiais, lógica fuzzy, sistemas espe-
cialistas, algoritmos genéticos, séries de Volterra, etc.
Este curso visa dar um enfoque no ferramental básico que foi desenvolvido
na década de 60. Com esta base os estudantes interessados são capazes de
"decolar"rumo à outros horizontes de métodos de ID, basta um pouco de
"transpiração".
1.3 Etapas de um problema de identificação
Um problema de identificação de um modelo matemático de um sistema
envolve no geral as seguintes etapas:
10
1. A primeira etapa é aquisição dos dados. Para este procedimento
deve-se escolher quais saídas e entradas são necessárias para cons-
trução de um modelo matemático, quantos N pares de entrada e/ou
saída são necessários, quais as taxas de amostragem devem ser usa-
das, número de canais, etc. O resultado deste processo é um conjunto
ZN = {u(t), y(t)}Nt=1. A partir destas medições é obtida uma função
representando o sistema gˆn
(
t, ZN
)
. Com está função estima-se a saída
yˆ como sendo:
yˆ = gˆ
(
t, ZN
)
(1.1)
A função gˆn(t, ZN) é escrita em termos de um número finito de parâ-
metros θˆ. Este θˆ representa a estrutura (forma) do modelo. Deve-se
notar que o desejo é que θˆ tenda ao valor e a estrutura real do sistema
θ, o que nem sempre pode ser facilmente obtido1.
2. A segunda etapa é decidir qual estrutura θˆ deve ser empregada.
3. Em seguida deve-se escolher os valores dos parâmetros θˆ de tal forma
que yˆ tenda ao valor real medido y, para isto deve-se ter que θˆ deve
tender a estrutura e os parâmetros reais θ. Esta etapa consiste de um
problema de otimização escritoda seguinte maneira:
Minimize
N∑
i=1
∥∥y(t)− gˆn(t, ZN)∥∥ (1.2)
4. Por fim a última etapa consiste em realizar uma validação e verificação
(V&V) do modelo. Este processo consiste em analisar eventuais erros
de modelagem, polarização e variância. Diversas técnicas são empre-
gadas em V&V, como por exemplo, comparação direta da saída do
modelo com dados experimentais, análise de correlação de resíduos,
etc. Deve ficar claro que cada uma desta técnicas devem ser realizadas
se empregando um conjunto de dados de ensaios diferentes dos usados
no procedimento de identificação do modelo.
Diversas técnicas podem ser empregadas nas etapas explicadas anterior-
mente. Durante este curso cada um destes passos serão estudados em detalhes
se apresentando os diferentes métodos utilizados.
1Em algumas aplicações pode até ser impóssivel.
11
Capítulo 2
Introdução à Probabilidade
Antes de começar a estudar métodos de ID é necessário revisar alguns
conceitos básicos de probabilidades que são utilizados nestas formulações. O
claro entendimento dos conceitos apresentados neste capítulo facilita muita
a ampla compreensão de ID, uma vez que praticamente todos os métodos a
serem estudados neste curso são estocásticos por natureza e devem ser descri-
tos por modelos envolvendo ruído. A saída de modelos estocástico envolvem
variáveis aleatórias. A seguir alguns teoremas e definições são apresentados
sobre o tema, muitos já estudados em cursos básicos de estatística e outros
mais avançados e tópicos de cursos introdutórios em processos estocásticos.
Processo estocástico se refere a uma grandeza que evolui no tempo sob ação
de algum mecanismo que é totalmente aleatório. Este capítulo é um resumo
conciso de notas de aulas do docente e que estão contidos em excelentes livros
introdutórios sobre o tema como [2].
2.1 Espaço amostral
Um espaço amostral Ω é definido como um conjunto de todos os possíveis
resultados.
Exemplo 2.1 Ω = {cara, coroa}; Ω = {1, 2, 3, 4, 5, 6}; Ω = {x ∈ <; 0 ≤
x ≤ 1}
Por outro lado, um evento é qualquer subconjunto de Ω pertencente a
uma classe de subconjuntos. Nestes casos ω é um ponto genérico do espaço
amostral Ω e φ é um evento impossível.
12
2.1.1 Álgebra de conjuntos
Para trabalhar com espaços amostrais são necessárias certas operações.
Imagine dois conjuntos E = {1, 2, 3} e F = {1, 2, 3, 5, 6}, assim:
União: E ∪ F = {1, 2, 3, 5, 6}, lê-se E + F ou E união com F .
Subconjunto: Note que E ⊂ F , o que significa que E é um subconjunto de
F .
Intersecção: E ∩ F = {1, 2, 3}, lê-se E e F .
Complemento: Ec é o conjunto de tudo que não está em E. Por exemplo,
E ∪Ec = Ω, o que significa que todo o espaço amostral é complemento
de E, ou ainda E ∩ Ec = φ, caso o complemento não tenha nenhum
elemento em comum.
Disjuntos: Dois conjuntos são disjuntos ou mutuamente exclusivos se estes
não tem nenhum elemento em comum, ou seja se E ∩ F for igual a φ.
Diferenças: E − F = E ∩ F c ou F − E = F ∩ Ec.
Com base nestes resultados associamos sempre: resultados de experimen-
tos com eventos e estes com conjuntos.
Exemplo 2.2 Seja A um conjunto de resultados x (face cara de uma moeda)
e B um conjunto de resultados y (face coroa de uma moeda). Neste caso nosso
espaço amostral é Ω = {cara, coroa}. temos que:
1. Ocorre x ou y ⇒ A ou B ⇒ A ∪B.
2. Ocorre x e y ⇒ A e B ⇒ A ∩B.
3. Não ocorre x ⇒ não A ⇒ A¯ ou Ac.
4. x e y mutuamente exclusivos ⇒ ou A ou B ⇒ A ∩B = φ
Uma coleção de subconjuntos de Ω formam um campo (ou anel) de Borel:
F. Um campo de Borel F é um campo não vazio que satisfaz:
• Se A ∈ F então Ac ∈ F.
• Se A,B ∈ F então A ∪B ∈ F.
• Ω ∈ F.
Um campo de Borel tem as seguintes propriedades:
• φ ∈ F.
• Se A e B ∈ F então A ∩B ∈ F e A−B = A ∩Bc ∈ F.
13
2.2 Probabilidade
O conceito de probabilidade é usado para dar peso a conjuntos diferen-
tes. Matematicamente, probabilidade é qualquer função real definida sobre
a classe F de Borel P : F→ < tal que para A,B ∈ F tem-se os axiomas:
• P (A) ≥ 0.
• P (Ω) = 1.
• Se A∩B = φ, ou seja A e B são mutuamente exclusivos, então P (A∪
B) = P (A) + P (B).
As propriedades de probabilidade são dadas por:
• P (φ) = 0.
• P (Ac) = 1− P (A).
• P (A ∪B) = P (A) + P (B)− P (A ∩B).
• Se A ⊂ B então P (A) ≤ P (B).
A partir destes conceitos iniciais é possível definir um espaço de proba-
bilidades como sendo formado por três variáveis (Ω,F, P ), sendo Ω o espaço
amostral, F uma classe de evento (anel de Borel) e P uma função de proba-
bilidade.
2.2.1 Probabilidade condicional
Seja B um evento tal que P (B) > 0 então para A:
P (A|B) = P (A ∩B)
P (B)
⇒ P (A ∩B) = P (A|B)P (B) (2.1)
sendo que a probabilidade condicional é lida como A probabilidade de A
dado que ocorreu B é a divisão da probabilidade de intersecção entre A e B
pela probabilidade de ocorrer B.
Exemplo 2.3 Qual a probabilidade de selecionar 3 cartas sem reposição e
obter três azes em um baralho com 52 cartas. Resposta: 4
52
3
51
2
50
= 24
132600
.
14
2.2.2 Independência
Independência é um conceito extremamente importante no estudo de pro-
cessos estocásticos. Dois evento A ∈ F e B ∈ F com P (A) > 0 e P (B) > 0
são ditos serem idependentes se:
P (A ∩B) = P (A)P (B) (2.2)
Se A e B são idenpendentes, então:
• P (A|B) = P (A).
• Se Ac e Bc, A e Bc, Ac e B são independentes.
2.2.3 Regra de Bayes
A regra de Bayes dá origem a um excelente estimador conhecido como
estimador de Bayes. Sua definição é:
P (Aj|B) = P (B|Aj)P (Aj)∑n
i=1 P (B|Ai)P (Ai)
(2.3)
2.2.4 Probabilidade total
Sejam A1, A2, · · · , An eventos mutuamente exclusivos (disjuntos), ou seja,
A1 ∩A2 ∩ · · ·An = φ com P (Ai) > 0 então a probabilidade total é dada por:
P (B) =
n∑
i=1
P (B|Ai)P (Ai) (2.4)
2.3 Variável aleatória
A primeira coisa a esclarecer é que uma variável aleatória X não é uma
variável e sim uma função com domínio no espaço amostral Ω, assim:
X : Ω→ < tal que: {ω ∈ Ω, X(ω) ∈ (−∞, x] ∈ F} (2.5)
Se Ω é enumerável1 a variável aleatória é discreta, caso contrário é contí-
nua. A variável aleatória é um evento de importância única e que associa-se
uma probabilidade a este X(ω).
1Um espaço amostral ser enumerável significa que ele é um conjunto infinito com cor-
respondência biunívoca com números inteiros positivos.
15
2.4 Função distribuição de probabilidade
Uma função distribuição de probabilidade (FDP2) é uma função cumula-
tiva monótona crescente definida como:
FX(x) , P ({X(Ω) ≤ x}) (2.6)
sendo {X(Ω) ≤ x} = {ω ∈ Ω, X(Ω) ∈ (−∞, x]} ∈ F, ou seja, é um
evento.
Exemplo 2.4 Variável aleatória tempo de espera t de um ônibus no PTI. A
função cumulativa neste caso pode ser dada por:
P (T (ω) ≤ t) = FT (t) (2.7)
Esta função é cumulativa e monótona crescente, pois se você espera um
tempo t → ∞ a probabilidade de não ocorrer o evento (ônibus chegar) é
pequena.
Uma FDP tem as seguintes propriedades:
• FX(−∞) = 0.
• FX(+∞) = 1.
• Se x2 > x1, então FX(x2) ≤ FX(x1) (função monótona não decres-
cente).
2.5 Função densidade de probabilidade
Caso uma FDP seja diferenciável é mais comum o emprego da função
densidade de probabilidade (PDF) pX(x):
pX(x) ,
d
dx
FX(x) (2.8)
Existem diferentes tipos de PDF, como exemplo a distribuição uniforme,
normal, etc. Se integrarmos uma PDF pX(x) obtemos a FDP. Observa-se
também que:
FX(x) =
∫ x
−∞
pX(y)dy = 1 (2.9)
2Não é o que vocês estão pensando!
16
Uma vez que FX(−∞) = 0. Como FX(x) é uma função monótona cres-
cente tem-se que pX(x) ≥ 0.
Um erro comum que os estudantes cometem é associar PDF com variáveis
aleatórias. Uma PDF não tem nada a ver com aleatoriedade. Na verdade
PDF é uma ferramenta matemática poderosa que é utilizada para caracterizar
variáveis aleatórias.
2.6 Variáveis aleatórias Gaussianas
Todos os testes de ID envolvem certo grau de aleatoriedade, nestes casosé impossível prever ruídos futuros, porém o comportamento estatístico de um
ruído podem ser modelado. Variáveis aleatórias Gaussianas são excelentes
modelos de ruídos. As variáveis aleatórias podem ser contínuas ou discretas,
exemplo:
Variáveis aleatórias contínuas: exponencial, gama, Cauchy.
Variáveis aleatórias discretas: binomial, geométrica, Poisson.
O primeiro conceito importante a reforçar é que uma variável aleatória
Gaussiana não necessariamente tem distribuição normal. Uma variável ale-
atória também pode ser bidimensional. Z ∼ N(0, 1) → é variável aleatória
normal 0− 1:
pZ(z) =
1√
2pi
exp
(
−z
2
2
)
(2.10)
com P{|Z| < 1} = 0.68.
2.7 Momentos estatísticos
Momentos estatísticos são métricas envolvendo variáveis aleatórias, onde
se pode fazer uma analogia com teorema dos momentos estudados em mecâ-
nica geral.
O momento estatístico de 1.o ordem, também conhecido como esperança
matemática é dado por3:
E(X) , x¯ =
∫
Ω
X(ω)dP (ω) (2.11)
3Considerando que a variável aleatória X é necessariamente integrável.
17
=
∫ +∞
−∞
xdFX(x)
=
∫ +∞
−∞
xpX(x)dx < +∞
Estimar a esperança matemática na prática pode ser difícil, uma vez que
se pode desconhecer o valor da PDF do sistema em análise.
Generalizando o momento estatístico para ordem k obtém-se uma função
geradora de momentos:
mk , E[Xk] =
∫ +∞
−∞
xkpX(x)dx < +∞ (2.12)
m1 é a média (x¯, µ), m2 é o segundo momento e assim por diante. Cada
um destes parâmetros servem para descrever características estatísticas de
um sinal. Conhecendo o valor da média µ pode-se escrever o momento cen-
trado de ordem k como:
ck , E[(X − µ)k] =
∫ +∞
−∞
(x− µ)kpX(x)dx (2.13)
c1 ≡ 0, já o momento centrado de ordem 2 (c2) é a variância (σ2), sendo
σ o desvio padrão. O momento centrado de 3.o ordem é conhecido como
kurtosis e descreve de forma simplista se uma distribuição é mais ou menos
achatada. Já o momento centrado de 4.o ordem é chamado de skewness é
descreve o grau de simetria de uma distribuição.
2.8 Correlação e covariância
Para um vetor de variáveis aleatórias x ∈ <n, a matriz de covariância é
uma grandeza de dimensão n× n definida por:
cov[x] = E[(x− E[x])(x− E[x])T ] (2.14)
que pode ser escrita como:
cov[x] = E[xxT ]− E[x]E[xT ] (2.15)
Já a correlação é definida como:
Rxy = E[xy] (2.16)
Para quantificar o grau de correlação (semelhança) entre dois sinais é
utilizado o coeficiente de correlação:
18
ρxy =
cov(x,y)
σx, σy
(2.17)
Se este valor for 0 não há tendência, agora se for 1 ou próxima de 1 os
dois sinais são altamente correlacionados. A partir destas ideias se define a
função de autocorrelação (FAC) de x(t) é:
rxx(τ, t) = E[x(t)x
∗(t+ τ)] (2.18)
se x(t) for real x∗(t+ τ) = x(t+ τ). A função de autocovariância de x(t)
é:
cxx(τ, t) = rxx(τ, t)− E[x(t)]E[x∗(t+ τ)] (2.19)
Já a função de correlação cruzada (FCC) entre dois sinais x(t) e y(t) é:
rxy(τ, t) = E[x(t)y
∗(t+ τ)] (2.20)
A função de covariância cruzada de x(t) e y(t) é:
cxy(τ, t) = rxy(τ, t)− E[x(t)]E[y∗(t+ τ)] (2.21)
As funções de autocorrelação e correlação cruzada são métricas que
quantificam o grau de semelhança entre dois sinais ou o mesmo sinal (no
caso de autocorrelação) em função da defasagem entre eles. Estes conceitos
são muito usados em métodos de ID em várias etapas:
• Podem ser usados para auxiliar se existe correlação entre dados can-
didatos de entrada e saída, uma vez que pode ser extremamente com-
plicado se definir na prática quem são os sinais de entrada e saída em
plantas complexas.
• Empregados em processos de verificação e validação de modelos.
• Usados para caracterizar sinais aleatórios e ruídos brancos ou coloridos.
2.9 Processo estocástico
Um processo estocástico, pode ser expresso graficamente por um con-
junto de testes com amostras aleatórias discretas xk[n] com k = 1, 2, .., K
realizações e n = 1, 2, ..., N pontos cada. Neste caso só é possível analisar
as características médias deste processo. A fig. (2.1) mostra um exemplo
gráfico de processo estocástico.
19
Fig. 2.1: Realizações de sinais medidos em um processo estocástico.
As métricas utilizadas para descrever as características de processos esto-
cásticos são os momentos estatísticos. Entre os momentos estatísticos mais
importantes se destacam as funções de autocorrelação (FAC) rxx(n,m)4 que
na forma geral podem ser estimadas por:
Rxx(n,m) = lim
K→∞
1
K
K∑
k=1
xk[n]xk[n+m] (2.22)
e as funções de correlações cruzadas (FCC)
RFx(n,m) = lim
K→∞
1
K
K∑
k=1
Fk[n]xk[n+m] (2.23)
sendo m o número de atrasos temporais e n o número de amostras, sendo
que no caso contínuo são descritos por τ . É interessante notar que a FAC é
a média do produto entre xk[n] e xk[n+m]5 e a FCC é a média do produto
entre duas seqüência diferentes Fk[n] e xk[n+m].
4No caso apresentado aqui elas são definidas para variáveis discretas (sequencias).
5Estas equações são idênticas ao caso contínuo apresentados nas eqs. (2.18) e (2.20).
20
2.10 Processo estacionário
Um processo estocástico é dito estacionário se suas propriedades estatís-
ticas não variam com o tempo (se mantém constante). A fig. (2.2) apresenta
a resposta de um sistema estacionário. Caso se divida este sinal em várias
partes e se calcule a distribuição de probabilidade em cada uma destas par-
tes irá se constatar que a distribuição estatística é a mesma, conforme a fig.
(2.3). A estacionariedade é relacionada a invariância do sistema e é uma
propriedade do processo e não do sinal em si.
Fig. 2.2: Exemplo de um processo estacionário.
2.11 Processo ergódico
Um processo estocástico é dito ergódico quando as propriedades médias
calculadas no tempo para qualquer realização são iguais às propriedades cal-
culadas a partir das médias do conjunto. Assim as FAC e FCC de processos
estacionários e ergódicos se tornam dependentes apenas dos atrasos m6, por-
tanto rxx(m,n) = rxx(m) e rxy(m,n) = rxy(m). Existem vários métodos
temporais para se estimar as correlações (pois dificilmente elas são conheci-
das por serem baseadas na definição de um limite). Um dos métodos mais
conhecidos é o método de Levinson-Durbin. A rigor deveríamos utilizar os
termos função de autocovariância e função de covariância cruzada, que são
iguais as FAC e FCC, mas retirando o efeito da média.
6No caso contínuo o atraso é τ .
21
Fig. 2.3: Distribuição de partes de um processo estacionário.
2.12 Ruídos brancos e coloridos
Ruído em ID corresponde as parcelas dos sinais que não é de interesse
para determinar as características dinâmicas do sistema. Em geral o ruído
apresenta natureza aleatória, como por exemplo dois casos especiais: o ruído
branco e o ruído colorido.
2.12.1 Ruído branco
Um ruído branco e(k) é um sinal puramente aleatório caracterizado por
uma FAC que satisfaz ree(m) = 0, ∀m 6= 0. O espectro de um ruído branco
tem potência em todas as frequências. É importante observar que ruído
branco apresenta uma característica estocástica. Isto significa na prática
que em um ruído branco toda a banda de frequencias tem importância. O
ruído ser branco não significa conhecimento sobre sua distribuição. Um ruído
branco e(k) pode ter diferentes distribuições de probabilidade. O interessante
em ID é que os resíduos7 tenham características de ruído branco, ou seja, que
7Lembrando que resíduo é a diferença entre a saída de um modelo identificado com os
dados experimentais.
22
possua sua FAC ree(m) = 0, ∀m 6= 0, sendo que no atraso zero a FAC possui
valores não nulos. É impossível determinar o valor em uma amostra k + 1.
2.12.2 Ruído colorido
Um ruído colorido já é um sinal que contém uma parte deterministica e
uma parte aleatória (branca). Um ruído colorido pode ser modelado como
um ruído branco sendo processado por um filtro, por exemplo, um processo
auto-regressivo commédia móvel. Ao contrário do ruído branco, um ruído
colorido não contém potência em todas as frequências. A maior parte da
energia de um ruído colorido está contida em bandas estreitas de frequencia.
2.13 Exercícios
Ex. 2.1 Suponha que a ocorrência ou não de chuva dependa das condições
de tempo do dia exatamente anterior. Admita-se que se chove hoje, choverá
amanhã com probabilidade 0.7 e que se não chover hoje, choverá amanhã
com probabilidade 0.4. Sabendo-se que choveu hoje, calcule a probabilidade
de que choverá depois de amanhã.
Ex. 2.2 Utilizando o Matlab ou Scilab gere um sinal aleatório u(k) com
distribuição normal ou uniforme. Use os comando randn e rand. Verifique
se sua FAC é realmente ruu(m) = 0, ∀m 6= 0, estimando a FAC usando
m = 20 atrasos.
Ex. 2.3 O arquivo dados1.mat fornecido pelo docente contém dois sinais no-
meados de y1 e y2. Qual destes sinais é ruído branco e qual é ruído colorido?
Justifique sua resposta usando o Matlab8. Não se esqueça de retirar o efeito
da média no sinal no momento de calcular as FACs.
Ex. 2.4 O arquivo dados2.mat contém três sinais x1, x2 e x3. Defina qual
destes sinais são correlacionados entre si. Use o mesmo programa imple-
mentado anteriormente. Aqui deverá ser estimada a função de correlação
cruzada entre estes sinais.
8Dica: Utilize os comandos no Matlab iddata (para escrever no formato do toolbox
ident), detrend (para remover tendências) e covf (para calcular a matriz de covariâncias)
para auxiliar nesta escolha. Use o help. A FAC do sinal pode ser obtida por R(4,:) saída
do comando covf
23
Capítulo 3
Representações Matemáticas
Lineares
A ideia dos métodos de ID é representar uma planta ou processo
físico qualquer por um modelo matemático, seja ele paramétrico ou não-
paramétrico. Por modelo não-paramétrico se entende uma representação
onde as informações da planta são completamente descritas por um gráfico,
por exemplo: resposta ao impulso, resposta em frequência, resposta ao de-
grau unitário, etc. Já um modelo paramétrico é um modelo que descreve a
partir de equações o comportamento dinâmico e/ou estático de um sistema.
Muitos dos conceitos deste capítulo não são novos para os estudantes, muitos
vistos em cursos de Sistemas Dinâmicos e Controle. Este é o caso da repre-
sentação linear mais básica usada em ID, a função de transferência1. Porém
outras representações também são possíveis, como a resposta em frequência
e a resposta ao impulso, assim como modelos em tempo discreto2 represen-
tados por equações a diferenças, como AR, ARX, ARMA, etc. Grande parte
das técnicas de ID a serem vistas neste curso são baseadas em modelos line-
ares discretos. Sendo assim, a meta desta capítulo é revisar representações
matemáticas clássicas de modelos e introduzir o uso de modelos AR, ARX,
ARMA, ARMAX, etc. em análise de sistemas dinâmicos lineares.
3.1 Função de transferência
A função de transferência representa uma relação entre dados de entrada
e saída em um sistema. Por definição, função de transferência H(s) é a
1Nos seus casos contínuos e discretos.
2Este tópico, provavelmente, deve ser o mais novo para todos os alunos, portanto, aqui
deverá se exigir maior esforço no momento de estudar.
24
transformada de Laplace da resposta ao impulso h(t) de um sistema:
H(s) =
∫ +∞
−∞
h(t)e−stdt (3.1)
Na eq. (3.1) o limitante inferior da integral é −∞ se assumindo que a
resposta ao impulso do sistema possa ser não-causal. Uma vez conhecida
H(s) se conhece totalmente o sistema de interesse sendo possível calcular
rapidamente a resposta deste frente a outras excitações. Para se obter a
função de transferência podemos adotar várias abordagens, e. g.:
• Aplicando a definição da transformada de Laplace na equações diferen-
cial que descreve o fenômeno físico de interesse.
• Aplicando técnicas de análise espectral ou outros métodos nos dados
de entrada/saída de uma planta e estimando sua respectiva função de
transferência3.
• Aplicando a definição da transformada de Laplace, caso se conheça a
resposta ao impulso h(t) do sistema4.
O problema típico de modelagem de sistemas dinâmicos é encontrar qual
a função de transferência que pode representar bem o comportamento de
uma planta visando análise ou controle. Uma função de transferência, nor-
malmente é escrita em função de uma razão de polinômios, seja em s ou
z.
3.1.1 Caso contínuo - domínio s
Uma função de transferência de sistema contínuo em uma forma polinô-
mial pode ser escrita na forma:
H(s) =
N(s)
D(s)
=
∑q
i=0 bis
i∑n
j=0 ajs
j
=
Y (s)
U(s)
(3.2)
sendo Y (s) e U(s) as transformadas de Laplace dos sinais de entrada
e saída do sistema. Deve ficar claro para o estudante que apesar de ser
teoricamento possível descrever uma função de transferência envolvendo U(s)
e Y (s) está abordagem não é na prática usada, uma vez que transformada
de Laplace não se aplica a sinais e sim para sistemas. Para sinais o melhor é
usar a transformada de Fourier.
3Esta técnica será vista mais adiante neste curso.
4Neste caso deve-se determinar a priori quem é h(t).
25
As raízes do polinômio N(s) (numerador de H(s)) são chamadas de zeros
e são dados por s = z1, z2, . . . , zq. De forma semelhante as raízes do polinômio
D(s) (denominador de H(s)) são conhecidos como pólos, s = p1, p2, . . . , pn.
Estes pólos e zeros podem ter multiplicidade.
A ordem de uma função de transferência é determinada pela ordem do
número de pólos s, ou seja, ordem n. Para sistemas físicos reais a ordem n
é sempre igual ou maior do que a ordem q (n ≥ q) para sistemas próprios
ou n > q para sistemas estritamente próprios. Um sistema ter mais zeros
do que pólos significa que a sua resposta ao impulso começa a responde em
instante menores do que zeros, o que leva ao conceito de sistema não-causal,
ou seja, um sistema passa altas para frequencias ω →∞.
Os pólos e zeros podem ser reais ou complexos. Se forem complexos eles
aparecem em pares complexos conjugados, ou seja, α ± jβ. Em sistemas
dinâmicos reais sempre aparecem pólos complexos conjugados, uma vez que
equações diferenciais são descritas por parâmetros reais (massa, resistência,
capacitância, amortecimento, etc.). Os pólos de uma função de transferência
trazem informações importantes sobre a dinâmica de um sistema, como por
exemplo, amortecimento ξ e frequencia natural ωn, além de informações sobre
a estabilidade do sistema. Se pelo menos um pólo tiver parte real Re[] maior
do que zero:
Re[pi] > 0, i− 1, · · · , n (3.3)
tal função H(s) é instável. No caso de zeros Re[zi] > 0, i − 1, · · · , n
a função de transferência é chamada de fase miníma. Uma função de trans-
ferência pode ser decomposta na forma:
H(s) =
∏q
i=1 (s− zi)∏n
j=1 (s− pj)
(3.4)
Sendo o resíduo de H(s) escrito como:
Ji = H(s)(s− pi)|s=pi (3.5)
O conceito de resíduo é bastante usado em decomposição modal, que
por sua vez pode ser usada para obter a resposta temporal de um sistema.
A decomposição modal, também conhecida como decomposição em frações
parciais, se refere a descrever um sistema como um somatório de módulos
básicos, conhecidos como modos. Matematicamente, a decomposição modal
é equivalente a forma canônica de Jordan. É possível que em um sistema
real existam alguns modos que sejam mais importantes do que outros no que
se refere a características dinâmicas que estes representam. Sendo assim, é
fundamental se definir de início qual a faixa de frequências pretende-se extrair
26
alguma informação dos sinais para então propor um modelo contendo os
modos que são capazes de descrever o comportamento que se deseja modelar.
Por exemplo, imagine que você deseja amortecer oscilações indesejadas em
um sistema elétrico qualquer. Um modelo completo para este sistema pode
compreender diferentes modos. Se a sua meta é se concentrar em uma faixa
de frequências estreita, obviamente que seu modelo deve conteros modos que
são dominantes e importantes naquela faixa.
O processo de decomposição modal compreende basicamente um desaco-
plamento dos modos dominantes do sistema. Existem três casos:
1. pólos reais repetidos.
2. pólos complexos não repetidos.
3. pólos com multiplicidade maior do que um.
Apenas para ilustração o caso para sistemas com pólos complexos não
repetidos é descrito por:
H(s) =
n∑
i=1
(
Ji
(s− pi) +
J∗i
(s− p∗i )
)
(3.6)
sendo o símbolo ∗ usado para representar o complexo conjugado5.
3.1.2 Caso discreto - domínio z
O caso discreto é equivalente a utilizar a definição da tranformada z em
cima da resposta ao impulso discreta h(k):
H(z) =
+∞∑
k=−∞
h(k)z−k (3.7)
Mais uma vez se define um caso geral assumindo que h(k) representa
a resposta ao impulso em tempo discreto t = kT , sendo T o período de
amostragem, de um sistema não-causal. Análogo ao domínio contínuo s, este
sistema também pode ser descrito em termos dos parâmetros de uma equação
a diferenças, relações de entrada U(z) e saída Y (z), ou entre relação de pólos
e zeros no domínio z.
O caso discreto da função de transferênciaH(z) é diretamente relacionado
aos conceitos de modelos em tempo discreto. A rigor H(z) é uma represen-
tação relacionada no domínio da frequencia , uma vez que a transformada
5a função matlab usada para calcular os resíduo é residue.
27
de Fourier de sequencias é um caso particular da transformada envolvendo
a variável z, que é tomada em todo plano complexo. Esta relação é especial
z = ejωn. A transformada z toma todos os valores do plano complexo e não
apenas na circunferência de raio unitário. No domínio de tempo discreto é
possível propor uma representação muito útil e poderosa em ID: H(q), sendo
q−1 um operador de atraso.
3.2 Função de resposta ao impulso
Aplicando a transformada inversa de Laplace na resposta na função de
transferência H(s) obtém-se a função de resposta ao impulso (IRF) h(t):
h(t) =
1
2pii
∫ +∞
−∞
H(s)estds (3.8)
Uma série de métodos pode ser empregado para estimar a IRF h(t) a
partir dos sinais de entrada e saída da planta. Alguns destes métodos serão
revistos no capítulo 4 deste texto.
3.2.1 Caso contínuo - Integral de convolução
Uma das grandes vantagens do uso de transformadas em aplicações de
engenharia é a seguinte relação:
Y (s) = H(s)U(s) (3.9)
No domínio temporal isto equivale a uma operação de convolução envol-
vendo a IRF h(t) e o sinal de excitação u(t). No domínio contínuo a operação
de convolução é escrita na forma integral:
y(t) =
∫ +∞
−∞
h(τ)u(t− τ)dτ =
∫ +∞
−∞
u(τ)h(t− τ)dτ = h ∗ u (3.10)
sendo τ um atraso temporal variante entre −∞ até +∞. Assim caso se
conheça a IRF e o sinal de excitação u(t) é possível descrever a resposta do
sistema sem necessariamente se resolver a equação diferencial governante do
processo físico em questão.
28
3.2.2 Caso discreto - Somatório de convolução
No caso discreto a resposta y(n) de um sistema com uma função de res-
posta ao impulso h(k) pode ser calculada pelo somatório de convolução con-
siderando o caso causal:
y(k) =
∞∑
j=0
h(j)u(k − j) (3.11)
sendo u(k) o sinal de excitação. Caso a resposta ao impulso tenha du-
ração infinita, h(k) é chamada função de resposta ao impulso infinita (IIR).
Normalmente, todo o sistema dinâmico (com realimentação de pólos) é es-
sencialmente do tipo IIR. Entretanto, muitas vezes são bastante amortecidos
e h(k) pode ter uma curta duração, podendo ser aproximada por uma res-
posta ao impulso finita (FIR). No caso FIR o modelo da eq. (3.12) pode ser
truncado em N amostras:
y(k) =
N∑
j=0
h(j)u(k − j) (3.12)
O somatório de convolução envolvendo um modelo FIR é uma represen-
tação discreta muito importante em ID e será revisto mais a frente. A partir
da representação da eq. (3.12) pode-se propor métodos para estimação de
h(k), assumindo conhecidos a entrada u(k) e a saída y(k), como por exemplo,
método dos mínimos quadrados ou método das covariâncias.
3.3 Realização no espaço de estados
A utilização de funções de transferência para representação de sistemas é
mais indicada para sistemas do tipo SISO6. O caso de um sistema multivariá-
vel (MIMO) deve empregar o conceito de matriz de função de transferência
H(s). Mesmo assim, a função de transferência, seja ela SISO ou MIMO, não
fornece informação sobre o que acontece dentro do sistema. Em outras pala-
vras, não fornece informação sobre as variáveis de estado x de um sistema.
Um bom modelo que fornece informações internas do sistema é a realização
no espaço de estados. Esta representação é temporal e muito útil em várias
situações:
• Representar sistemas MIMO.
6Single Input-Single Output.
29
• Representar sistemas não-lineares.
• Diversas técnicas para projeto de controladores robustos para plantas
descritas no espaço de estados.
• Utilização de técnicas de filtragem para estimação de estados.
Um caso típico de um modelo no espaço de estados é descrito pela ex-
pressão a seguir:
x˙ = Ax + Bu (3.13)
y = Cx + Du
sendo x ∈ <n o vetor de estados com dimensão n, x˙ a derivada do vetor
de estados, u ∈ <r o vetor de entradas formado por r funções temporais, y
o vetor de saídas medidas e A, B, C, D são as matrizes que representam a
realização do sistema sendo chamadas de matriz dinâmica, de entradas, de
medidas e de transmissão direta, respectivamente. A informação interna do
sistema está contida no vetor de estados.
A matriz dinâmica A contém toda a informação contida nos pólos, como
frequencias naturais e fatores de amortecimento, já as matrizes B e C contém
informações sobre os zeros e os modos de um sistema. A matriz de transmis-
são direta D é usada nos casos onde a resposta medida depende diretamente
da matriz de entrada, por exemplo quando se mede aceleração. As matrizes
A, B, C, D são normalmente assumidas serem constantes, mas podem variar
no tempo dependendo do tipo de aplicação.
O problema de identificação de uma realização no espaço de estados com-
preende encontrar quais matrizes A, B, C e D são necessárias para descrição
de um modelo do comportamento dinâmico uma vez conhecido o valores de
u e y. Deve ficar claro para o estudante que a o vetor de estados x não ne-
cessariamente precisa ser uma grandeza mensurável. Normalmente, o vetor
de estados não contém um significado físico e em muitas aplicações deve ser
estimado usando algum filtro, conhecido como observador de estados. Ob-
servadores de estado clássicos são o de Luenberg (determinístico) e o filtro
de Kalman (estocático). Outra questão interessante é que existem infinitas
realizações no espaço de estados para o conjunto de dados de entrada e saída
de uma planta, ou seja, inúmeras matrizes A, B, C, D podem representar
o mesmo sistema. Isto significa que estas podem ser decompostas de formas
diferentes, como por exemplo: realização modal, realização balanceada, reali-
zação na forma canônica, etc. Estas transformações lineares da realização no
espaço de estados são úteis para resolver e minimizar questões relacionadas a
30
convergência e condicionamento numérico que são encontrados em problemas
de análise e projeto de controladores.
Inúmeros métodos de ID podem ser empregados para obtenção de uma
realização no espaço de estados de uma planta. Métodos clássicos envolvem
o uso de realização de autosistemas. Um método poderoso nesta classe é o
método ERA, que nada mais é do que uma formulação baseada na construção
de uma matriz de Toeplitz utilizando a resposta ao impulso IRF estimada de
dados de entrada e saída. Com a matriz de Toeplitz montada realiza-se uma
decomposição em valores singulares e aplicando algumas definições chega-
se a uma realização no espaço de estados de um modelo para um processo
físico real. Este método foi proposto por Juang [3] visando aplicações em
identificação de modeloas para estruturasaeroespaciais, como satélites.
Também é importante observar que é possível realizar transformações
entre realizações diferentes, por exemplo, a partir do conhecimento da reali-
zação no espaço de estados se obter a função de transferência correspondente:
H(s) = C (Is−A)−1 B + D (3.14)
sendo I a matriz de identidade com ordem n.
Uma variante discreta para a realização no espaço e estados é descrita
por:
x(k + 1) = Φx(k) + Γu(k) (3.15)
y(k) = Cx(k)
sendo x(k) ∈ <n o vetor de estados em tempo discreto k, u(k) ∈ <r é o
vetor de entradas com dimensão r no instante k, y(k) ∈ <p é o vetor de saída
com dimensão p e as matrizes Φ, Γ e C a realização no espaço de estados
em tempo discreto. Também é possível ter uma matriz de transmissão direta
para o caso discreto, que foi omitida aqui por simplicidade. Similarmente ao
caso contínuo a matriz de função de transferência H(z) no domínio discreto
pode ser calculada a partir do conhecimento da realização discreta:
H(z) = C (zI−Φ)−1 Γ (3.16)
Conhecida uma realização discreta e possível transforma-la para contí-
nuo ou vice-versa conhecendo-se qual o período de amostragem do sistema
T . Inúmeros métodos podem ser empregados para isto, como método da in-
variância da resposta impulsiva, método de Tustin (transformação bilinear),
etc.7
7Consulte os comandos c2d e d2c no Matlab para informação adicional.
31
3.4 Representações matemáticas lineares de
tempo discreto
Este é o ponto mais importante deste capítulo: representações de siste-
mas lineares em tempo discreto. Duas razões são usadas para justificar esta
afirmação:
• Essencialmente todos os processos físicos são contínuos, porém os dados
que são mensuráveis destes processos são normalmente amostrados e
armazenados de forma discreta. Nada mais natural do que tratá-los na
forma em que são processados.
• Como poderá ser visto no decorrer desta apostila, a totalidade dos mé-
todos de ID a serem apresentados neste curso irá se basear em modela-
gem para sistemas dinâmicos estocásticos. Neste caso será fundamental
que o modelo contenha parcelas que sejam capazes de modelar todos
as variáveis aleatórias que possam aparecer nestes processos. Modelos
discretos para ruídos são temas clássicos e com muita coisa já proposta
e testada em processos reais.
Esta seção tem como meta apresentar as representações matemáticas li-
neares em tempo discreto mais importantes e usadas em ID.
Considere um sistema linear, invariante com o tempo e estacionário que
possa ser representados como na fig. (3.1)
Fig. 3.1: Representação de um sistema dinâmico linear, invariante com o
tempo e estacionário com ruído colorido adicionado a saída.
A transformada z da saída do sistema pode ser descrita por:
Y (z) = G(z)U(z) +X(z) = G(z)U(z) +H(z)W (z) (3.17)
32
sendo G(z) a função de transferência da parte determinística do sistema,
X(z) um processo estocástico estacionário com densidade espectral racional
(ruído colorido) X(z) = H(z)W (z), ou seja, um ruído branco W (z) proces-
sado por uma função de transferência H(z) representado a parte estocástica
do sistema. Todas as funções de transferência são racionais e estáveis. Em
particular H(z) tem as seguintes propriedades:
• H(z)−1 é estável, ou seja, o sistema e de fase mínima.
• limz→∞H(z) = 1, ou seja, o sistema é causal.
Supondo o caso da parte determinística e estocástica conterem pólos em
comum, pode-se propor um modelo geral no domínio z do tipo:
Y (z) =
z−dB(z)
F (z)A(z)
+
C(z)
D(z)A(z)
w(z) (3.18)
sendo A(z), B(z), C(z), D(z), F (z) polinômios em z, cuja raízes são pó-
los e zeros das partes estocásticas H(z) e determinísticas G(z) do sistema e d
o eventual atraso de transporte puro do sistema. Note que um sistema deste
tipo tem memória, ou seja, para descrever uma entrada y(k) são necessárias
informações ocorridas em instantes passados. Obviamente trabalhar no do-
mínio z não é interessante, uma vez que seria mais proveitoso trabalhar no
tempo discreto diretamente com os sinais coletados experimentalmente y(k)
e u(k). Sendo assim, aplicando a transformada z inversa usando o operador
de atraso temporal8q−1 :
A
(
q−1
)
y(k) = q−d
B (q−1)
F (q−1)
u(k) +
C (q−1)
D (q−1)
w(k) (3.19)
sendo u(k) o sinal de entrada, y(k) o sinal de saída e w(k) um ruído
branco. Os polinômios que definem o modelo sãod dados por:
A
(
q−1
)
= 1 + a1q
−1 + · · ·+ anaq−na
B
(
q−1
)
= b0 + b1q
−1 + · · ·+ bnbq−nb
C
(
q−1
)
= 1 + c1q
−1 + · · ·+ cncq−nc
D
(
q−1
)
= 1 + d1q
−1 + · · ·+ dndq−nd
F
(
q−1
)
= 1 + f1q
−1 + · · ·+ fnf q−nf
A partir do modelo geral da eq. (3.19) diferentes modelos são obtidos na
literatura. A seguir uma descrição sucinta destes modelos.
8Note que q−1y(k) = y(k − 1)
33
3.4.1 Modelo FIR
Uma importante representação de saída envolve a soma de convolução,
que já foi vista anteriormente. A convolução entre a IRF h(k) o sinal de
entrada u(k) pode ser escrito como9:
y(k) =
∞∑
j=0
h(j)u(k − j) + w(k) (3.20)
O modelo acima é conhecido como modelo de resposta infinita ao impulso
(IIR). Caso o sistema seja estável, o que significa que existe um valor de M
tal que h(k) possa ser truncado de forma a ter |h(k)| < �, ∀k > M , sendo �
um valor pequeno tendendo a zero, é possível truncar e obter um modelo de
resposta finita ao impulso:
y(k) =
M∑
j=0
h(j)u(k − j) + w(k) (3.21)
M é o número de termos da resposta ao impulso h(k). O modelo FIR
pode ser obtido com a eq. (3.19) considerando que A(q−1) = C(q−1) =
D(q−1) = F (q−1) = 1 e B(q−1) é um polinômio de ordem M , assim nb = M .
Com isto um modelo FIR é descrito como:
y(k) = B
(
q−1
)
u(k − d) + w(k) (3.22)
Importante observar que este modelo contém pólos somente na origem.
Sendo assim, um sistema dinâmico só pode ser representado exatamente por
um modelo FIR, se este for estável e sua resposta ao impulso decair rapi-
damente, uma vez que a ordem do modelo nb está associada a quantidade
de amostras temporais da resposta ao impulso. A ideia do modelo FIR é
descrever um sistema dinâmico por um modelo envolvendo apenas zeros em
posições arbitrárias e pólos somente na origem. Uma vez que o modelo não
contém pólos, não há problemas de estabilidade, o que faz com que este tipo
de representação linear seja muito usada em controle adaptativo de processos.
Uma forma poderosa de expandir modelos FIR é empregando uma base
ortogonal, ou seja, descrever um modelo FIR considerando a imposição de
uma função ortogonal com pólos correspondentes a dinâmica dominante do
sistema. Com isto, a ordem do polinômio B(q−1) pode ser drasticamente
reduzida e o problema de estimação dos parâmetros deste polinômio pode
9Assumindo que a saída y(k) seja contaminada como um ruído branco w(k).
34
ser conduzido de forma numérica mais simples10. O artigo [8] mostra uma
aplicação de modelos FIR com uma base ortogonal conhecida como filtro
de Kautz aplicando os resultados utilizando o método das covariâncias que
deverá ser estudado no próximo capítulo deste texto.
3.4.2 Modelo auto-regressivo
O modelo mais simples que pode ser descrito pela eq. (3.19) é o modelo
auto-regressivo (AR). Um modelo AR é função apenas do polinômio A (q−1)
e descrito por:
A
(
q−1
)
y(k) = w(k) (3.23)
Um modelo AR pode ser usado para descrever um modelo com entradas
desconhecidas. Uma vez que o ruído w(k) aparece explicitamente no modelo
AR, este é classificado na classe de modelos com erro na equação. Note que:
y(k) =
1
A (q−1)
w(k) (3.24)
Este modelo contém apenas pólos em posições arbitrárias. Diversos mé-
todos podem ser usados para determinar o polinômio A(q−1), entre eles o
método dos mínimos quadrados e o estimador de Yule-Walker. Estes méto-
dos serão estudados em detalhes no capítulo 5 desta apostila. Um comando
Matlab que pode ser empregado na modelagem via AR é a função ar. É
muito comumusar a representação modelo AR(na), ou seja um modelo AR
com ordem na.
3.4.3 Modelo auto-regressivo com entradas externas
A estrutura de ummodelo auto-regressivo com entradas externas (ARX11)
corresponde a seguinte equação:
A
(
q−1
)
y(k) = B
(
q−1
)
u(k − d) + w(k) (3.25)
que também representa um modelo com erro na saída. Interessante ob-
servar que:
y(k) =
B (q−1)
A (q−1)
u(k − d) + 1
A (q−1)
w(k) =
B (q−1)
A (q−1)
u(k − d) + e(k) (3.26)
10Um dos temas propostos para o projeto do curso será a implementação de um modelo
FIR em base ortogonal.
11AutoRegressive with eXogenous inputs.
35
assim o ruído e(k) que é adicionado a saída do modelo AR não é branco
e sim um ruído branco w(k) filtrado por um modelo AR. Em um modelo
ARX(na, nb) há somente relações algébricas entre a previsão da saída e as
entradas passadas. Um modelo AR já é descrito por pólos e zeros, ressal-
tando que o ruído pode ser processado por um filtro AR e não faz parte do
regressor. O capítulo 5 contém métodos que serão estudados para obtenção
dos polinômios A(q−1) e B(q−1) de um modelo ARX. Deve ser observado que
o modelo AR é um caso particular de modelo ARX.
3.4.4 Modelo auto-regressivo com média móvel e entra-
das externas
O modelo ARMAX12 é descrito por:
A
(
q−1
)
y(k) = B
(
q−1
)
u(k − d) + C (q−1)w(k) (3.27)
ou ainda como:
y(k) =
B (q−1)
A (q−1)
u(k − d) + C (q
−1)
A (q−1)
w(k) (3.28)
O modelo ARMAX(na, nb, nc) também é um modelo de erro na equação.
A perturbação não é branca e modelada como sendo um ruído branco w(k)
processado por um filtro ARMA:
e(k) =
C (q−1)
A (q−1)
w(k) (3.29)
O modelo de ruído é mais completo do que o caso ARX. Por sua vez,
a estimação do modelo exige estimativa dos três polinômios A(q), B(q) e
C(q). Um problema de ordem prática é a montagem do regressor para es-
timar C(q). Este problemas pode ser resolvido via mínimos quadrados com
matriz estendida ou via variáveis instrumentais e será estudado em detalhes
no capítulo 5 desta apostila.
Um caso particular do modelo ARMAX é o modelo ARMA, quando há
desconhecimento do sinal de excitação13 de uma planta u(k):
A
(
q−1
)
y(k) = C
(
q−1
)
w(k) (3.30)
Um modelo ARMA(na, nc) descreve como um ruído branco w(k) afeta
a saída y(k) com um filtro formado por A(q) e C(q). Modelos ARMA são
12AutoRegressive Moving Average with eXogenous inputs.
13O que pode ser muito comum em aplicações industriais.
36
muito usados como modelos de séries temporais, ou seja, situações em que
somente se tem um sinal medido, y(k). Exemplo: série temporal de preços
de um produto, temperatura de um lago ao longo de uma década, etc.
3.4.5 Modelo de erro na saída
Um modelo de erro na saída (OE14), também conhecido como erro com
estrutura paralela é usado quando o único ruído que afeta um sistema é um
ruído do tipo branco. Assim:
y(k) =
B (q−1)
F (q−1)
u(k − d) + w(k) (3.31)
Para este modelo ser estável as raízes de F (q) devem estar dentro do
CRU. Um modelo OE pode ser interpretado como a superposição de um
ruído branco com um processo ARX sem ruído.
3.4.6 Modelo de Box-Jenkins
O modelo de Box-Jenkins é dado por:
y(k) =
B (q−1)
F (q−1)
u(k − d) + C (q
−1)
D (q−1)
w(k) (3.32)
sendo que este modelo também é do tipo OE, porém o ruído adicionado é
colorido. Além disto a parte determinística e a parte estocástico do modelo
não têm pólos em comum, ao contrário do modelo ARMAX.
3.5 Exercícios
Ex. 3.1 Escreva todos os modelos AR, ARX, ARMA, ARMAX, OE e Box-
Jenkins usando representações por diagramas de blocos, tanto no domínio z
quando usando o operador de atraso q−1.
Ex. 3.2 Seja uma função de transferência H(z) assintoticamente estável.
Usando H(z) e um sinal de entrada qualquer u(k) obtém-se por simulação
y(k). A este sinal é acrescentado ruído branco w(k) a fim de representar
possíveis erros de medição. Escolha uma das representações discretas que
melhor descreveria este sistema. Que modificação seria necessária fazer no
ruído para poder utilizar uma representação ARX?
14Output error
37
Ex. 3.3 Considere um sistema discreto descrito pela seguinte equação a di-
ferenças:
y(n) = x(n) + ax(n− 1) + bx(n)x(n− 1) (3.33)
sendo x(n) um sinal de entrada, y(n) um sinal de saída e a uma constante
real não-nula.
a) Determine o sistema quanto as propriedades e em função do valor da cons-
tante real b, justificando a sua resposta: memória, linearidade, invâriância
com o tempo discreto, estabilidade e causalidade.
b) Se o sistema for linear e invariante com o tempo discreto para algum valor
de b, calcule a resposta ao impulso deste sistema.
c) Que tipo de representação discreta tem este sistema?
Ex. 3.4 Considere um sistema discreto linear e invariante com o tempo dis-
creto descrito pela equação a diferenças?
y(n) + 0.5y(n− 1) = x(n) + x(n− 2). (3.34)
a) Descreve este sistema com uma representação discreta adequada, ressal-
tando os coeficientes de seus polinômios.
b) Calcule a resposta ao impulso causal.
c) Calcule a resposta ao impulso não-causal, supondo h(n) = 0 para n ≥ 3.
d) Calcule a função de transferência H(z) deste sistema.
e) Calcule a resposta do sistema para a entrada x(n) = ejnpi/2.
f) Calcule a resposta do sistema causal para a entrada:
x(n) = 2δ(n) + δ(n− 1) (3.35)
sendo δ(n) um impulso unitário em um instante n.
Ex. 3.5 Considere um sistema discreto descrito por:
y(n) =
2∑
k=0
bkx(n− k) + c (3.36)
sendo x(n) uma sequencia de entrada, y(n) uma sequencia de saída, c
uma constante real finita e bk (k = 0, 1, 2) constantes reais finitas e não
nulas.
a) Descreve este sistema com uma representação discreta adequada, ressal-
tando os coeficientes de seus polinômios.
38
b) Determine o sistema quanto as propriedades e em função do valor da
constante real c, justificando a sua resposta: memória, linearidade, invâri-
ância com o tempo discreto, estabilidade e causalidade.
c) Se o sistema for linear e invariante com o tempo discreto para algum valor
de b, calcule a resposta ao impulso deste sistema.
d) Calcule a função de transferência para este sistema H(z).
Ex. 3.6 Considere um sistema linear e invariante com o tempo discreto des-
crito pela seguinte equação a diferenças:
y(n) + ay(n− 1) = x(n) + bx(n− 1) (3.37)
com condições iniciais nulas. Suponha que as constantes a e b são reais
e diferentes de zero.
a) Calcule a resposta ao impulso causal.
b) Calcule a resposta ao impulso não-causal.
c) Calcule a função de transferência H(z).
d) Qual a condição para o sistema causal ser estável? Jusitfique.
Ex. 3.7 Considere um sistema linear e invariante com o tempo com uma
IRF do tipo h(n) = anu(n), sendo u(n) a resposta ao degrau unitário e com
x(n) = u(n) − u(n − N). Com base nestas informações pede-se a obten-
ção da resposta y(n) calculando como convolução discreta entre h(n) e x(n).
Obtenha o mesmo resultado usando a transoformada z.
39
Capítulo 4
Métodos Não-Paramétricos
Como já ressaltado anteriormente, os métodos não-paramétricos tem
como meta obter uma representação gráfica do sistema a partir de dados
de entrada e/ou saída de uma planta. Esta representação pode ser uma
resposta ao degrau, impulso, resposta em frequência, etc. O termo não-
paramétrico se refere ao tipo de saída obtida com esta classe de métodos:
uma curva gráfica. Esta resposta, por sua vez, pode ser identificada atra-
vés de duas classes de métodos1: determínisticos e estocásticos. Este capí-
tulo tem como objetivo apresentar alguns métodos comuns de obtenção de
modelos não-paramétricos. Inicialmente, são revistos alguns métodos clássi-
cos de identificação de características de sistemas de primeira e segunda or-
dem usando técnicas determínisticas, além do procedimento de identificação
usando a soma de convolução. Na sequencia, os métodos estocásticos, en-
volvendoa equação de Wiener-Hopf, são apresentados. Em cada uma destas
classes, tanto os procedimentos temporais quanto espectrais são explorados.
4.1 Métodos determinísticos
Os métodos determinísticos assumem que a quantidade de ruído presente
nos dados são desprezíveis e tratam as entradas/saídas do sistema como si-
nais de processos puramente determínisticos. Vários métodos podem ser
usados na obtenção das características de sistemas físicos usando estes mé-
todos, sendo que os mais comuns são revistos a seguir. Estas formulações se
mostram úteis em muitas aplicações industriais, desde que atenda a certos
requisitos, como ter relação sinal/ruído relativamente alta.
1Estes métodos se referem ao tipo de algoritmo usado na estimativa do modelo não-
paramétrico.
40
4.1.1 Sistema de 1.o ordem
Talvez o caso mais simples em identificação seja considerar a dinâmica de
sistemas de 1.o ordem. Um sistema típico destes tem a seguinte função de
transferência2:
H(s) =
Y (s)
U(s)
=
K
τs+ 1
(4.1)
sendoK o ganho do sistema e τ uma constante de tempo. O conhecimento
destes dois parâmetros identifica o modelo como um todo. Uma forma de
obter este modelo é aplicar como sinal de entrada um degrau de amplitude
A. Caso o sistema tenha um nível de ruído baixo estes os parâmetros K e τ
podem ser estimados facilmente. O ganho K pode ser expressado como:
K =
y(∞)− y(0)
A
(4.2)
sendo y(∞) o valor em regime permanente do sinal y(t). A constante
de tempo τ pode ser determinada a partir do tempo em que a resposta ao
degrau atinge 63.2 % da amplitude em regime permanente:
y(τ) = 0.632 (y(∞)− y(0)) + y(0) (4.3)
O eventual atraso de transporte no sistema pode ser estimado com este
método assumindo que ele é igual ao período entre a aplicação do degrau e o
início da resposta, mas suas considerações são muito restritivas. A fig. (4.1)
apresente uma exemplo de estimação gráfica do ganho K e da constante de
tempo τ a partir do conhecimento da resposta ao degrau.
Com a fig. (4.1) observa que o sistema possui um τ = 0.1 segundos e um
ganho K = 10, aplicando as eqs. (4.3) e (4.2).
4.1.2 Sistema de 2.o ordem
Um sistema típico de segunda ordem de um sistema subamortecido3 tem
a seguinte função de transferência:
H(s) =
Y (s)
U(s)
=
Kω2n
s2 + 2ξωns+ ω2n
(4.4)
sendo que ωn é a frequência natural do sistema e ξ o fator de amorteci-
mento. Este dois parâmetros juntamente com o ganho K definem totalmente
2Que provavelmente já deve ter sido exaustivamente estudada.
3Que apresenta resposta oscilante e contém pólos complexos conjugados.
41
Step Response
Time (sec)
A
m
pl
itu
de
0 0.1 0.2 0.3 0.4 0.5 0.6
0
2
4
6
8
10
12
System: H
Time (sec): 0.0994
Amplitude: 6.3
System: H
Time (sec): 0.577
Amplitude: 9.97
Fig. 4.1: Resposta ao degrau unitário (A = 1) de um sistema de primeira
ordem.
o sistema. O caso de interesse aqui se refere a sistemas pouco amortecidos, o
que na prática significa que ξ << 1. Também se assume por simplicidade que
K = 1. A resposta temporal ao degrau unitário de um sistema de segunda
ordem pode ser descrita por:
y(t) = 1− 1√
1− ξ2 e
−ξωntsen
(
ωnt
√
1− ξ2 + tan−1
(√
1− ξ2
ξ
))
(4.5)
Note que na fig. (4.2) são descritos alguns parâmetros que descrevem o
comportamento dinâmico de um sistema e podem ser usados para analisar
qualitativamente se um sistema tem comportamento adequado ou não, de
acordo com especificações de projeto. Uma destas medidas é o sobre-sinal,
mais conhecido pelo termo em inglês overshoot OS. Este valor é dado pelo
máximo valor da resposta menos o valor desta quando o sistema entra em
regime permanente
OS = ymax (t)− 1 = exp
(
−ξpi√
1− ξ2
)
, (4.6)
o overshoot ocorre exatamente em um tempo de pico tp descrito como
tp =
pi
ωn
√
1− ξ2 . (4.7)
42
0 0.5 1 1.5 2 2.5
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
x 10−3
X: 0.1013
Y: 0.001778
Tempo [s]
x
(t)
X: 1.759
Y: 0.0009952
OS
t
s
tp
Fig. 4.2: Exemplo de resposta ao degrau unitário para um sistema com um
grau de liberdade.
Outra característica importante é o período de oscilações Td dado por
Td =
2pi
ωn
√
1− ξ2 = 2tp. (4.8)
Por fim o tempo de ajuste, ts, define o tempo em que a resposta do
sistema atinge o regime permanente dentro de um intervalo de ±5%4. Uma
aproximação para ts pode ser escrita como
ts =
3
ωnξ
. (4.9)
É importante observar que a partir das equações anteriores é possível
identificar um sistema com um determinado fator de amortecimento ξ e
freqüência natural ωn de acordo com os parâmetros de tempo de ajuste,
overshoot, período de oscilações e tempo de pico para conduzirem a uma
resposta com características e forma desejada.
Outro método determinístico de identificação envolve usar a resposta li-
vre. A resposta livre de um sistema de 2.o ordem levemente amortecido pode
4Há definições para ts quando este intervalo é ±3%.
43
ser descrita no domínio temporal por:
y(t) = Y e−ξωntsen (ωdt+ φ) (4.10)
sendo Y a amplitude da saída, ωd = ωn
√
1− ξ2 a frequencia natural
amortecida e φ = tan−1
(√
1−ξ2
ξ
)
.
Utilizando a resposta ao degrau é possível identificar o tempo de esta-
belecimento, tempo de pico, overshoot, ξ, ωn, etc. O mesmo pode ser feito
usando a resposta livre (resposta a condições iniciais ou resposta ao impulso).
A utilização da resposta ao degrau neste caso já foi estudada no curso de
controle linear I, portanto, a enfase básica aqui será apresentar o decremento
logarítmico.
O decremento logarítmico δ é definido como o logarítmo natural da razão
de duas amplitudes sucessivas. Considere a resposta y(t) do caso subamor-
tecido (0 < ξ < 1) visto na fig. (4.3). O decremento logarítmico δ é escrito
como:
δ = ln
(
y(t)
y(t+ td)
)
, (4.11)
sendo td = 2pi
ωd
o período entre duas oscilações sucessivas, onde ωd é a
freqüência angular natural amortecida.
Para um caso geral tem-se:
δ = ln
(
y0
y1
)
= ln
(
y1
y2
)
= ln
(
yn−2
yn−1
)
, (4.12)
sendo n o número de oscilações realizadas. A Eq. (4.12) pode ser rescrita
da forma:
eδ =
y0
y1
=
y1
y2
=
yn−2
yn−1
=
yn−1
yn
. (4.13)
Notando que y0
yn
= y0
y1
y1
y2
y2
y3
· · · yn−1
yn
pode-se escrever a relação:
enδ =
y0
yn
. (4.14)
Com isto obtém-se uma nova expressão para o decremento logarítmico δ
em função do número de ciclos n realizados no movimento oscilatório
δ =
1
n
ln
(
y0
yn
)
. (4.15)
Lembrando que a resposta de um sistema subamortecido é do tipo
44
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
−8
−6
−4
−2
0
2
4
6
8
10
x 10−3
Tempo [s]
x
(t)
 [m
]
x1
x2
x3
x0
x4
Fig. 4.3: Resposta de sistema subamortecido evidenciando amplitudes suces-
sivas.
y(t) = Y e−ξωntsen (ωdt+ φ) (4.16)
Substituindo a Eq. (4.16) na Eq. (4.12) obtém-se a seguinte equação
δ = ln
(
y0
y1
)
= ln
(
Y e−ξωnt0sen (ωdt0 + φ)
Y e−ξωnt1sen (ωdt1 + φ)
)
, (4.17)
sendo t1 = t0 + td, onde td = 2piωd . Após algumas manipulações algébricas
na Eq. (4.17) chega-se a expressão do decremento logarítmico δ em função
do fator de amortecimento ξ
δ =
2piξ√
1− ξ2 , (4.18)
Ou ainda da forma
ξ =
δ√
4pi2 + δ2
(4.19)
Assim se conheço duas amplitudes sucessivas y0 e y1, ou se uma amplitude
y0 e uma amplitude yn após n ciclos, posso calcular o decremento logarítmico
δ entre elas e estimar com a Eq. (4.19) o fator de amortecimento ξ do sistema.
45
Exemplo 4.1 Considere um sistema de segunda-ordem representado por um
oscilador do tipo massa-mola-amortecedor com massa m = 20kg e desloca-
mento inicial y0 = 0.01 m. A fig. (4.4) mostra a resposta livre deste sistema.
Estime os coeficientes equivalentes de rigidez e amortecimento viscosodeste
sistema.
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
−8
−6
−4
−2
0
2
4
6
8
10
x 10−3
X: 0.06
Y: 0.004993
Tempo [s]
x
(t)
 [m
]
X: 0
Y: 0.01
Fig. 4.4: Resposta livre do sistema.
Solução: Considerando duas amplitudes sucessivas y0 = 0.01 m e y1 = 0.005
m, mostradas na fig. (4.4), o decremento logarítmico é calculado a seguir:
δ = ln
(
y0
y1
)
= ln
(
0.01
0.005
)
= 0.693.
Com o δ calculado emprega-se a Eq. (4.19) para se estimar o fator de amor-
tecimento ξ
ξ = δ√
4pi2+δ2
= 0.693√
4pi2+(0.693)2
= 0.11.
Como o fator de amortecimento ξ está entre 0 e 1, este sistema é subamor-
tecido. Sabendo que o período entre as duas oscilações sucessivas é td = 0.06
s, também visto na fig. (4.4), pode-se calcular a freqüência angular natural
amortecida
ωd =
2pi
td
= 104.7 rad/s.
46
O valor da freqüência angular natural é dado por
ωn =
ωd√
1−ξ2
= 104.7√
1−(0.1)2 = 105.3 rad/s.
A rigidez do sistema pode ser escrita lembrando que ωn =
√
k
m
, o que leva a
k=mω2n = (20) (105.3)
2 = 2.22× 105 N/m.
Já o coeficiente de amortecimento viscoso é estimado por:
c=2mωnξ = 2(20)(105.3)(0.11) = 4.63× 102 N.s/m
4.1.3 Método de Sundaresan
O método de Sundaresan é uma formulação clássica baseada em estima-
tivas de H(s) a partir da curva de resposta ao degrau. O sistema de interesse
deve ser escrito das seguintes formas:
H(s) =
e−τds
(τ1s+ 1) (τ2s+ 1)
(4.20)
H(s) =
e−τdsω2n
s2 + 2ξωns+ ω2n
(4.21)
A eq. (4.20) é mais indicada para representar sistemas dinâmicos lineares
sobreamortecidos, ou seja, sistemas que não oscilam quando se aplicada a
resposta ao degrau ou impulso. Isto significa que a função de transferência
H(s) é descrita apenas por pólos reais. Já a eq. (4.21) representa o caso de
um sistema dinâmico subamortecido, ou seja, que oscila devido a presença
de pólos complexos conjugados presentes em H(s). Em ambos os casos se
assume que possa existir um atraso de transporte τd no sistema.
O caso sobreamortecido tem os seguintes passos:
1. Determine o ganho de regime permanente do sinal dividindo a variação
do sinal de saída pela amplitude do degrau aplicado à entrada.
2. Determine a área dada pela integral:
m1 =
∫ ∞
0
(1− y(t)) dt (4.22)
o valor de m1 pode ser obtido com o comando area no matlab, uma
vez que y(t) é obtido experimentalmente.
3. Determine graficamente a inclinação Mi da tangente no ponto de
inflexão de y(t).
47
4. Determine tm que é a interseção da tangente com o valor em regime
permanente de y(t).
5. Determine λ a partir de λ = (tm −m1)Mi.
6. Determine η a partir da equação5:
λ = χe−χ (4.23)
χ =
lnη
η − 1 (4.24)
7. Por fim τ1, τ2 e τd são determinados a partir de:
τ1 =
η
η
1−η
Mi
(4.25)
τ2 =
η
1
1−η
Mi
(4.26)
τd = m1 − τ1 + τ2 (4.27)
Para o caso de sistema subamortecido o método de Sundaresan envolve
as etapas:
1. Determinação do ganho em regime permanente, dividindo a variação
do sinal de saída pela amplitude do degrau aplicado à entrada.
2. Deteminar a área do sinal m1 usando algum procedimento gráfico:
m1 =
∫ ∞
0
(1− y(t)) dt (4.28)
3. Estimar a inclinação da tangente no ponto de inflexão de y(t). Este
valor é Mi.
4. Determinar tm que a interseção da tangente como valor em regime
permanente do sinal.
5. Determine λ = (tm −m1)Mi.
5Que pode ser resolvida graficamente a partir da relação entre η e λ.
48
6. Determine ξ a partir da solução gráfica da equação:
λ =
cos−1ξ√
1− ξ2 e
−ξcos−1ξ√
1−ξ2 (4.29)
7. Determine ωn e τd a partir das expressões:
ωn =
cos−1ξ√
1− ξ2
1
tm −m1 (4.30)
τd = m1 − 2ξ
ωn
(4.31)
4.1.4 Identificação usando convolução
Até agora, os métodos não-paramétricos determinísticos estudados eram
baseados em excitar o sistema físico com entradas particulares, como entrada
degrau ou impulso. Na identificação de um processo real esta consideração
pode ser uma barreira prática, ou por ser impossível aplicar tais tipos de
sinais em uma planta ou pela necessidade de obtenção de um modelo mate-
mático do sistema a partir de dados de operação.
Uma forma muito útil de contornar esta limitação é empregar a soma de
convolução com o objetivo de identificar a resposta ao impulso discreto de
um sistema, h(k), a partir de sinais de entrada u(k) e saída da planta y(k).
Como visto anteriormente, a soma de convolução assumindo que o sistema é
do tipo IIR e causal é calculada por:
y(k) =
∞∑
j=0
h(j)u(k − j) (4.32)
Com a eq. (4.32) pode-se escrever o seguinte conjunto de equações:
y(0) = h(0)u(0) + h(1)u(−1) + h(2)u(−2) + . . .
y(1) = h(0)u(1) + h(1)u(0) + h(2)u(−1) + . . .
y(2) = h(0)u(2) + h(1)u(1) + h(2)u(0) + . . .
...
...
...
Ou em uma forma mais elegante usando matrizes:
49

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)
...
 (4.33)
A princípio a eq. (4.33) pode ser empregada para estimação da resposta ao
impulso h(k) do sistema. Em termos práticos se admite que a IRF h(k) possa
ser truncada em J amostras. Se for possível coletar N amostras temporais
de y(k) e u(k) tal que N >> J , a eq. (4.33) pode ser rescrita:

y(0)
y(1)
y(2)
...
y(N)

=

u(0) u(−1) u(−2) · · · u(−J)
u(1) u(0) u(−1) · · · u(−J + 1)
...
...
... . . .
...
...
...
... . . .
...
u(N) u(N − 1) · · · · · · u(−J +N)


h(0)
h(1)
h(2)
...
h(J)

(4.34)
Escrevendo de forma compacta:
y = Uh (4.35)
sendo que y é um vetor N × 1, U é uma matriz N × J e h é um vetor
J × 1. A matriz U é conhecida por matriz de Toeplitz.
Na eq. (4.35) a incógnita é o vetor h. Assim, a IRF h pode ser solucionada
via o método dos mínimos quadrados padrão, assim h = U−1y. Entretanto,
no geral, a relação da eq. (4.35) define uma matriz U de regressores em u(k)
que pode ser muito grande e mal condicionada numericamente, dificultando o
cálculo da inversa visando a estimativa de h(k). Para melhorar o condiciona-
mento numérico de U o sinal de excitação deve ser suficientemente ativo, ou
seja, capaz de excitar bem toda a dinâmica do sistema. Um sinal de entrada
deste tipo tem duas características interessante em ID:
Numérica: A inversão da matriz de Toeplitz U pode ser realizada de ma-
neira mais fácil, sem problemas de origem numérica.
Física: O sinal de excitação foi tal que conseguiu excitar de forma eficiente
toda a dinâmica de interesse no processo de identificação da planta.
Apesar da eq. (4.35) não ser apenas restrita a sinais de entrada u(k)
particulares, imagine o caso de se aplicar uma resposta ao impulso unitário
50
u(k) = 1 para k = 0 e u(k) = 0 para todo k 6= 0, então a matriz de Toeplitz
será diagonal, o que facilita bastante a identificação do sistema:
y = u(0)Ih (4.36)
sendo I a matriz identidade N × J . Assim:
h = y/u(0) (4.37)
o que obviamente equivale a dizer que h = y. Portanto a resposta ao
impulso, assim como a resposta ao degrau, tem característica de um sinal
suficientemente ativo, o que ajuda a explicar por que estes sinais são bastante
usados em análise de sistemas dinâmicos.
Infelizmente, o método apresentado aqui ainda é determinístico uma vez
que se assume que a relação sinal/ruído é alta. Caso se considere que o sinal
de saída y(k) = ye(k) + e(k), sendo que ye(k) é o sinal determinístico real e
e(k) um ruído de medição que pode ser do tipo branco ou colorido. Com isto
a eq. (4.37) torna-se:
h =
ye + e
u(0)
=
ye
u(0)
+
e
u(0)
(4.38)
A equação anterior mostra então que a IRF h(k) para um sistema é
composta por duas parcelas, a primeira real e a segunda causada pelo ruído,
conhecida como polarização. Assim, caso