Upload
others
View
1
Download
0
Embed Size (px)
Citation preview
INSTITUTO MILITAR DE ENGENHARIA
CAP LEONARDO OLIVEIRA DE ARAÚJO
IDENTIFICAÇÃO E CONTROLE DE ALGUMAS CLASSESDE SISTEMAS NÃO-ESTACIONÁRIOS
Dissertação de Mestrado apresentada ao Curso deMestrado em Engenharia Elétrica do Instituto Militarde Engenharia, como requisito parcial para obtenção dotítulo de Mestre em Ciências em Engenharia Elétrica.
Orientador: Paulo César Pellanda, Dr. ENSAECo-orientador: Roberto Ades, Dr PUC-Rio
Rio de Janeiro2006
c2006
INSTITUTO MILITAR DE ENGENHARIAPraça General Tibúrcio, 80-Praia VermelhaRio de Janeiro-RJ CEP 22290-270
Este exemplar é de propriedade do Instituto Militar de Engenharia, que poderá incluí-loem base de dados, armazenar em computador, microfilmar ou adotar qualquer forma dearquivamento.
É permitida a menção, reprodução parcial ou integral e a transmissão entre bibliote-cas deste trabalho, sem modificação de seu texto, em qualquer meio que esteja ou venhaa ser fixado, para pesquisa acadêmica, comentários e citações, desde que sem finalidadecomercial e que seja feita a referência bibliográfica completa.
Os conceitos expressos neste trabalho são de responsabilidade do(s) autor(es) e do(s)orientador(es).
A663 Araújo, Leonardo Oliveira deIdentificação e Controle de Algumas Classes de Sistemas
Não-estacionários / Leonardo Oliveira de Araújo. - Riode Janeiro : Instituto Militar de Engenharia, 2006.
146 p.: il, graf., tab.
Dissertação (mestrado) - Instituto Militar deEngenharia- Rio de Janeiro, 2006.
1. Sistemas, identificação e controle. 2. Sistemas não-estacionários. I. Título. II. Instituto Militar de Engen-haria.
CDD 621.319
2
INSTITUTO MILITAR DE ENGENHARIA
CAP LEONARDO OLIVEIRA DE ARAÚJO
IDENTIFICAÇÃO E CONTROLE DE ALGUMAS CLASSES DESISTEMAS NÃO-ESTACIONÁRIOS
Dissertação de Mestrado apresentada ao Curso de Mestrado em Engenharia Elétricado Instituto Militar de Engenharia, como requisito parcial para obtenção do título deMestre em Ciências em Engenharia Elétrica.
Orientador: Paulo César Pellanda, Dr. ENSAECo-orientador: Roberto Ades, Dr PUC-Rio
Aprovada em 31 de Janeiro de 2006 pela seguinte Banca Examinadora:
Paulo César Pellanda, Dr. ENSAE do IME - Presidente
Roberto Ades, Dr PUC-Rio do IME
Geraldo Magela Pinheiro Gomes, Dr. ENSAE do IME
Mario Cesar Mello Massa de Campos, Dr. ECP do CENPES
Rio de Janeiro2006
3
Àquele que é causa unívoca de todas as coisas, cha-mado por muitos nomes em tempos, locais e cul-turas diferentes: O Grande Arquiteto do Universo.
4
AGRADECIMENTOS
Aos amigos e orientadores, Maj Paulo César Pellanda e Maj Roberto Ades, exemplos
de dedicação, compromisso e notório saber, que sacrificaram inúmeras horas de lazer
familiar para possibilitar a existência deste trabalho.
Ao amigo Maj Juraci Ferreira Galdino pelas efetivas e claras orientações prestadas,
ultrapassando os limites de sua área de atuação científica específica, que possibilitou uma
abordagem de maior magnitude nesta dissertação.
Ao professor Marcelo Vieira Corrêa que contribuiu significativamente com elucidações
e arquivos de simulação para comparação de resultados sempre que solicitado, além de
permitir a exploração da tese de doutorado de sua autoria.
Ao Cel Ney Bruno por sua colaboração na modelagem do sistema do levitador mag-
nético, do qual é o idealizador.
À minha esposa, Rosane da Penha Andrade Araújo, pela compreensão e apoio du-
rante este período de extrema dedicação e aplicação ao estudos no qual fiquei envolvido.
Ao meu pai, José Alfredo de Araújo, e ao meu tio, Marlanfe Tavares Oliveira, pelos
sábios conselhos que me conduziram a este momento.
Ao Exército Brasileiro, em especial ao Instituto Militar de Engenharia, por esse
investimento tão significativo.
"Filho meu, não desprezes a instrução de tua mãe nem o ensinamento de teu pai,
ouça os teus maiores porque isto será uma coroa grandiosa para tua cabeça."
Provérbio, 1, 8
5
SUMÁRIO
LISTA DE ILUSTRAÇÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
LISTA DE TABELAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
LISTA DE SÍMBOLOS E ABREVIATURAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.1 Contexto e Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2 Objetivos da Dissertação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.3 Organização da Dissertação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2 FUNDAMENTOS TEÓRICOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 Projeção Estável e Instável de uma Realização . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.3 Redução de Ordem por Truncamento Balanceado . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4 Sistemas LPV e Quasi -LPV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.1 Modelagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.2 Análise de Estabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.5 Interpolação de Controladores Robustos de Estrutura Estimação-Controle . . 36
2.5.1 Formulação do Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.5.2 Estrutura Estimação-Controle Equivalente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5.3 Restrições Para a Seleção de Autovalores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.5.4 Parâmetro de Youla Dinâmico × Estático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.5.5 Continuação dos Subespaços Invariantes Selecionados . . . . . . . . . . . . . . . . . . . . 55
3 CONTRIBUIÇÕES PARA A IDENTIFICAÇÃO DE SISTEMAS 62
3.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.2 Representações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.3 Identificação de Sistemas Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3.1 Identificação usando Amostras Discretas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3.2 Identificação usando a Resposta em Freqüência . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.4 Identificação de Sistemas quasi -LPV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.4.1 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6
4 ESCALONAMENTO DE GANHOS EM ESTRUTURAS DO TIPO
ESTIMAÇÃO-CONTROLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.2 Escolha Sistemática dos Modelos de Síntese . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.2.1 Definição das Faixas de Linearização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.2.2 Transição entre Faixas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.3 Escalonamento de Ganhos por Imposição da Dinâmica de Malha Fechada . . . 81
4.3.1 Considerações sobre a Modelagem Não-Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.3.2 Controlador e Observador . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.3.3 Controlabilidade e Observabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.3.4 Ajuste do Rastreamento e do Ponto de Operação . . . . . . . . . . . . . . . . . . . . . . . . 85
4.4 Escalonamento de Ganhos por Cancelamento de Não-Linearidades . . . . . . . . . 86
4.4.1 Modelagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.4.2 Cancelamento das Não-Linearidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5 APLICAÇÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.1 Identificação de Modelos Lineares usando a Resposta Discreta no Tempo . . . 89
5.2 Identificação de Modelos Lineares usando a Resposta em Freqüência . . . . . . . 99
5.3 Identificação Quasi -LPV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.4 Controle de um Míssil Ar-Ar via Continuação de Subespaços Invariantes . . . . 112
5.5 Controle de um Sistema de Levitação Magnética por Técnicas de
Escalonamento de Ganhos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6 CONCLUSÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.1 Contribuições para a Identificação de Sistemas . . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.2 Contribuições para o Controle por Escalonamento de Ganhos . . . . . . . . . . . . . 127
6.3 Perspectivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
7 REFERÊNCIAS BIBLIOGRÁFICAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
8 APÊNDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
8.1 Apêndice 1: Modelo do Míssil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
8.2 Apêndice 2: Modelos Lineares Identificados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
8.3 Apêndice 3: Linearização do Modelo de um Levitador Magnético . . . . . . . . . . 145
7
LISTA DE ILUSTRAÇÕES
FIG.2.1 Sistemas em malha fechada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
FIG.2.2 Parametrização de Youla e a estrutura estimação-controle . . . . . . . . . . . . 40
FIG.4.1 Curva f(x) e sua linearização em x = 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
FIG.4.2 Curva f(x) e suas linearizações em x = 6: intervalo 4 ≤ x ≤ 8. . . . . . . . . 78
FIG.4.3 Erros nos processos de linearização. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
FIG.4.4 Transição entre faixas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
FIG.5.1 Sinal v aplicado as plantas físicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
FIG.5.2 Resposta em freqüência (módulo e fase) da dinâmica da planta . . . . . . . . 90
FIG.5.3 Curvas de respostas ao degrau para: planta e modelo. Curva de
erro entre as respostas do modelo e da planta: (A) Gm11(z), (B)
Gm12(z), (C) Gm
13(z) e (D) Gm14(z). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
FIG.5.4 Curvas de respostas ao degrau para: planta e modelo. Curva de
erro entre as respostas do modelo e da planta: (A) Gr15(z), (B)
Gr16(z), (C) Gm
17(z) e (D) Gm18(z). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
FIG.5.5 Resposta em freqüência (módulo e fase) da dinâmica da planta . . . . . . . . 95
FIG.5.6 Curvas de respostas ao degrau para: planta e modelo. Curva de
erro entre as respostas do modelo e da planta: (A) Gr21(z), (B)
Gr22(z) e (C) Gr
23(z). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
FIG.5.7 Curvas de respostas ao sinal v para: planta e modelo. Curva de
erro entre as respostas medida e estimada: Gr24(z) . . . . . . . . . . . . . . . . . . 98
FIG.5.8 Resposta em freqüência (módulo e fase): comparação entre a
planta e o modelo Gr32(s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
FIG.5.9 Resposta em freqüência (módulo e fase): comparação entre a
planta e o modelo Gm42(s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
FIG.5.10 Resposta em freqüência (módulo e fase): comparação entre a
planta e o modelo Gr52(s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
FIG.5.11 Resposta em freqüência (módulo e fase): comparação entre a
planta e o modelo Gm62(s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
FIG.5.12 Resposta em freqüência: comparação das amostras com os valores
fornecidos pelo modelo Gm72(s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8
FIG.5.13 Resposta em freqüência. Comparação das amostras com os valores
fornecidos pelos modelos (A) Gm81(s), (B) Gr
82(s) e (C) Gr83(s). . . . . . . . . 107
FIG.5.14 Entrada para validação do modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
FIG.5.15 Resposta da planta não-linear e do modelo identificado quasi -LPV
de um míssil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
FIG.5.16 Resposta da planta não-linear e do modelo identificado quasi -LPV
de um míssil entre 35s e 45s de vôo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
FIG.5.17 Estrutura de controle para o exemplo do míssil . . . . . . . . . . . . . . . . . . . . . . 113
FIG.5.18 Resposta ao degrau do sistema linearizado em malha fechada em
três pontos de operação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
FIG.5.19 Resposta do sistema não-linear em malha fechada para três tipos
de interpolação de controladores - aceleração . . . . . . . . . . . . . . . . . . . . . . . 119
FIG.5.20 Resposta do sistema não-linear em malha fechada para três tipos
de interpolação de controladores - aceleração: ampliação de parte
da FIG. 5.19. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
FIG.5.21 Resposta do sistema não-linear em malha fechada para três tipos
de interpolação de controladores - ângulo de ataque . . . . . . . . . . . . . . . . . 120
FIG.5.22 Resposta do sistema não-linear em malha fechada para três tipos
de interpolação de controladores - ângulo de ataque: ampliação
de parte da FIG. 5.22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
FIG.5.23 Resposta do sistema não-linear em malha fechada para três tipos
de interpolação de controladores - sinal de controle δc . . . . . . . . . . . . . . . 121
FIG.5.24 Resposta do sistema não-linear em malha fechada para três tipos
de interpolação de controladores - sinal de controle δc: ampliação
de parte da FIG. 5.23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
FIG.5.25 Comparação das saídas controladas com o sinal de referência. . . . . . . . . . 123
FIG.5.26 Comparação entre os sinais de controle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
FIG.8.1 Diagrama físico do míssil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
FIG.8.2 Sistema de levitação magnética de uma esfera de aço . . . . . . . . . . . . . . . . . 145
FIG.8.3 Valores medidos experimentalmente em laboratório: f(N) ×x(mm)× i(A) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
9
LISTA DE TABELAS
TAB.5.1 Dados referentes aos modelos identificados de GM1(s) . . . . . . . . . . . . . . . . 91
TAB.5.2 Valores dos coeficientes da FT contínua da planta utilizada . . . . . . . . . . . 94
TAB.5.3 Valores dos coeficientes da FT discreta da planta utilizada . . . . . . . . . . . . 94
TAB.5.4 Dados referentes aos modelos identificados . . . . . . . . . . . . . . . . . . . . . . . . . . 95
TAB.5.5 Custos referentes aos modelos identificados . . . . . . . . . . . . . . . . . . . . . . . . . 96
TAB.5.6 Valores da função custo da FT GM3(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
TAB.5.7 Valores da função custo da FT GM4(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
TAB.5.8 Valores da função custo da FT GM5(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
TAB.5.9 Valores da função custo da FT GM6(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
TAB.5.10 Valores da função custo da FT GM7(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
TAB.5.11 Valores dos coeficientes da FT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
TAB.5.12 Comparativo de custos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
TAB.5.13 Parâmetros do filtro W(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
TAB.5.14 Ganhos do controlador e do estimador aumentado . . . . . . . . . . . . . . . . . . 117
TAB.5.15 Pólos em malha fechada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
TAB.5.16 Distribuição dos pólos em malha fechada . . . . . . . . . . . . . . . . . . . . . . . . . . 118
TAB.8.1 Valores dos coeficientes da FT Gm21(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
TAB.8.2 Valores dos coeficientes da FT Gr21(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
TAB.8.3 Valores dos coeficientes da FT Gm22(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
TAB.8.4 Valores dos coeficientes da FT Gr22(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
TAB.8.5 Valores dos coeficientes da FT Gm23(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
TAB.8.6 Valores dos coeficientes da FT Gr23(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
TAB.8.7 Valores dos coeficientes da FT Gm24(z) = Gr
24(z) . . . . . . . . . . . . . . . . . . . . . 142
TAB.8.8 Valores dos coeficientes da FT Gm81(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
TAB.8.9 Valores dos coeficientes da FT Gm82(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
TAB.8.10 Valores dos coeficientes da FT Gr82(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
TAB.8.11 Valores dos coeficientes da FT Gm83(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
TAB.8.12 Valores dos coeficientes da FT Gr83(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
TAB.8.13 Diferenças entre as curvas do gráfico força × deslocamento (não-
linear) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
10
LISTA DE SÍMBOLOS E ABREVIATURAS
R - conjunto dos números reais
C - conjunto dos números complexos
Rn - conjunto dos vetores reais de dimensão n
Rm×n - conjunto das matrizes reais m× n
j - índice ou j=√−1
∈ e ∈ - pertence e não pertence
, - por definição
≈ e ∼= - aproximadamente igual a
× - produto cartesiano
∃ - existe
∀ - para todo
| - tal que
→ - tendendo a
=⇒ - implica em
co - envelope convexo: para Si ∈ Rm×n, i = 1, 2, ..., N,
co{S1, S2, ...SN} , {∑Ni=1 αiSi : αi ≥ 0,
∑Ni=1 αi = 1}
|α| - valor absoluto de α ∈ R ou de α ∈ C
spec(M) - espectro da matriz M ∈ Rn×n : {λi(M) : i = 1, ..., n}Re(α) - parte real de α ∈ C
In - matriz indentidade n× n
MT - transposta da matriz M
M−1 - inversa da matriz M
M † - pseudo-inversa, ou inversa de Penrose-Moore, da matriz M
σ(M) - valor singular máximo
σ(M) - valor singular mínimo de M
‖α‖2 - Norma 2 de α ∈ Cn
x(t) - vetor de sinais contínuos, onde para ∀t, x ∈ Rn
x(k) - sinal x ∈ Rn no instante (discreto) k
s - variável de Laplace
z - variável da transformada Z para os sistemas discretos
w - freqüência em rad/s
θ(t) - variável de interpolação
11
G(z) =
A B
C D
realização em espaço de estados da função de transferência discreta
G(z) = C(zI − A)−1B + D
G(s) =
A B
C D
realização em espaço de estados da função de transferência con-
tínua G(s) = C(sI − A)−1B + D
12
SIGLASLTI Linear Invariante no Tempo
LTV Linear Variante no Tempo
LPV Linear a Parâmetros Variáveis
LFT Transformação Linear Fracionária (do inglês Linear Fractional
Transformation)
LMI Desigualdade Matricial Linear
LQG Gaussiano Quadrático Linear
LQ Quadrático Linear
FT Função de Transferência
PRCBI Parameter Robust Control by Bayesian Identification
TF Transformada de Fourier
DFT Transformada discreta de Fourier (do inglês Discrete Fourier
Transform)
FFT Transformada rápida de Fourier (do inglês fast Fourier transform)
NARMAX Modelo não-linear auto-regressivo, de média móvel com entradas
exógenas (do inglês nonlinear autoregressive moving average model
with exogenous input)
13
RESUMO
Esta dissertação apresenta a identificação e o controle de uma classe de sistemasdinâmicos não-estacionários com representação quasi-LPV. Técnicas de identificação sãodefinidas e simuladas. Três metodologias do tipo escalonamento de ganhos clássico sãodesenvolvidas.
As duas primeiras técnicas de identificação exploradas são lineares. Uma fornececomo resultado uma função de transferência discreta e a outra resulta em uma função detransferência contínua. A terceira técnica permite a identificação de sistemas dinâmicosquasi -LPV e depende de uma identificação preliminar de pontos de operação do sistema.Esta última metodologia é aplicada com sucesso na identificação de um modelo não-linearrealista de um míssil ar-ar.
Em relação às técnicas de controle, a primeira contribuição do trabalho é possibi-litar uma escolha sistemática da família de modelos de síntese de controladores LTI aser interpolada, para uma planta não-linear. A segunda técnica apresentada utiliza umaabordagem quasi -LPV para modelar o sistema. Através de uma realimentação de esta-dos estimados, o objetivo é fixar os autovalores de matriz dinâmica de malha fechada. Aterceira metodologia busca o cancelamento do comportamento não-linear de uma mode-lagem restrita de sistemas não-estacionários.
Os resultados das simulações não-lineares indicam a eficiência das técnicas apresenta-das para a identificação e o controle de sistemas dinâmicos não-estacionários que possuamrepresentação quasi-LPV.
14
ABSTRACT
This dissertation deals with the identification and the control of nonlinear dynamicalsystems having a quasi-LPV representation. Techniques for identification are given andsimulated. Three methodologies of classic gain scheduling are developed.
The two first identification techniques are linear. One of them obtains a discretetransfer function as a solution and the other one results in a continuous transfer func-tion. The third method allows the identification of quasi-LPV dynamical systems andis dependent of a set of LTI plant known a priori. This technique is applied to theidentification of a realistic nonlinear model of an air-to-air missile.
As a with regard to control technique, the first contribution is a method systematicchoice of a set of LTI controller synthesis models to be scheduled for a nonlinear plant.The second technique use a quasi-LPV description of a system. The objetive is to fixthe eigenvalues of the closed-loop dynamic matrix by using a parameter dependent statefeedback law. The third control method aims to cancel nonlinearities for a restrictednonlinear plants state feedback.
The nonlinear simulation results indicate that the techniques presented in this workare efficient for the identification and the control of nonlinear dynamical systems havinga quasi-LPV representation.
15
1 INTRODUÇÃO
1.1 CONTEXTO E MOTIVAÇÃO
O desenvolvimento da engenharia de controle está intrinsecamente ligada à constru-
ção e análise de modelos que visam reproduzir, o mais fielmente possível, comportamentos
observados na natureza. Tais modelos possibilitam validar teorias e mostrar o funciona-
mento da dinâmica observada. Nos estudos de átomos, do sistema nervoso humano ou
de uma galáxia, por exemplo, modelos são fundamentais e indispensáveis.
Dentre os diversos modelos usados no decorrer da história da humanidade, o mo-
delo matemático destaca-se no atual contexto científico. Este modelo é definido por um
conjunto de equações diferenciais (tempo contínuo) ou equações de diferenças (tempo
discreto) que descrevem a dinâmica de uma planta física a ele relacionado. A complexi-
dade da modelagem, de forma geral, aumenta a medida que o nível de precisão requerida
na representação das dinâmicas envolvidas é maior.
Os modelos matemáticos permitem correlacionar causas e efeitos (AGUIRRE, 2000a),
possibilitando a análise dos dados observados. Porém, a modelagem exata e precisa de
toda a física de uma planta é de altíssima complexidade, senão impossível. No entanto,
dentro de tolerâncias pré-estabelecidas para cada caso, um modelo que tenha sucesso na
reprodução do comportamento de uma planta estudada é fundamental na engenharia,
em particular para as técnicas de controle.
Um sistema dinâmico tem sua análise feita geralmente nos domínios do tempo e da
freqüência. Por esse motivo, os métodos de identificação de plantas podem ser divididos
em (LJUNG, 1987):
• Métodos Paramétricos;
• Métodos Não-Paramétricos;
• Métodos no Domínio da Freqüência.
A técnica paramétrica assume que a estrutura do modelo do sistema é conhecida com
base nas leis físicas que regem a sua dinâmica, embora os parâmetros sejam considerados
desconhecidos. Por outro lado, a técnica não-paramétrica é necessária quando a estrutura
do modelo matemático não está disponível. O resultado dessa técnica não é um modelo
16
matemático, mas uma representação gráfica que caracteriza a dinâmica do sistema em
estudo (AGUIRRE, 2004). Por último, as técnicas no domínio da freqüência geram
modelos que mostram a evolução da dinâmica identificada nesse domínio, utilizando a
transformada de Fourier (HOUGEN, 1972).
Há ainda a possibilidade de se aplicar outras transformadas que possibilitem a iden-
tificação do sistema dinâmico analisado. Como exemplo, pode ser citada a identificação
utilizando a transformada Wavelets que é descrita em (COCA, 1995a,b, 1996a; BAKSHI,
1993), que são trabalhos precursores no desenvolvimento dessa técnica. Portanto, a mo-
delagem no domínio da freqüência é uma particularidade de uma modelagem no domínio
de uma determinada transformada.
Por serem mais simples, os modelos lineares foram fundamentais no desenvolvimento
de técnicas de identificação. Em contraste com a complexidade apresentada pelos mo-
delos não-lineares, o uso de modelos lineares é justificável (BILLINGS, 1980) pela sua
simplicidade, facilidade de obtenção e por possuírem um amplo ferramental matemático
na engenharia de controle. Esses modelos podem ser empregados em faixas operativas
restritas, ou seja, próximos dos pontos de operação onde o modelo foi validado.
A necessidade de obtenção de modelos mais precisos do comportamento de sis-
temas dinâmicos levou a uma modelagem que considera faixas mais amplas de operação
(BILLINGS, 1980). Em conseqüência, as equações obtidas nessa nova modelagem con-
sideram as não-linearidades dos modelos. Aponta-se ainda a dificuldade de identificar os
sistemas não-lineares e a impossibilidade da indicação de uma técnica que seja capaz de
apresentar soluções gerais aceitáveis. O desenvolvimento dos modernos computadores e
a disponibilização de um amplo ferramental matemático são aliados poderosos na mani-
pulação e análise dos dados medidos em plantas dessa natureza.
Em (CORRÊA, 2001), são ressaltados os estudos feitos na década de oitenta sobre a
identificação de sistemas dinâmicos não-lineares, dando particular atenção às técnicas de
representação NARMAX (do inglês Nonlinear Autoregressive Moving Average Model with
Exogenous Input) polinomial (LEONTARITIS, 1985), NARMAX racional (BILLINGS,
1989) e redes neurais (NARENDA, 1990). Na década seguinte, aquele mesmo autor
cita diversos trabalhos de identificação de plantas com modelos não-lineares1, bem como
apresenta as motivações e técnicas empregadas nos trabalhos desenvolvidos.
Neste contexto, o uso de informações a priori é enfatizada, pois provoca o surgimento
1Ver (BILLINGS, 1989; VALLVERDU, 1992; NOSHIRO, 1993; PRÖLL, 1994; JANG, 1994; ?; CU-
BILLOS, 1997; GENÇAY, 1997; TSOI, 1999) .
17
de uma nova classificação de modelagem. Segundo (CORRÊA, 2001), dependendo do
nível e/ou tipo de informação utilizada a priori, as técnicas de obtenção de modelos
passam a ser classificadas como (HERBERT, 1993; SJÖBERG, 1995; BOHLIN, 1995):
• modelagem caixa-branca;
• modelagem caixa-cinza;
• modelagem caixa-preta.
A modelagem caixa-branca, segundo (GARCIA, 1997), engloba os casos onde os
dados de entrada e saída do sistema são dispensáveis, pois o comportamento da dinâmica
do sistema é equacionado a partir do conhecimento da estrutura da planta. Nesse caso, os
parâmetros envolvidos possuem um significado físico. A modelagem caixa-preta, por sua
vez, é atribuída quando os modelos são gerados baseando-se exclusivamente na análise das
relações entre entradas e saídas do sistema dinâmico. Nesse contexto, não há informações
a priori, uma vez que a estrutura da planta é ignorada. Vantagens e desvantagens são
apontadas em ambos os tipos de modelagens, como pode ser visto em (POTTMANN,
1998) e (TIKHONOV, 1977) para caixas-brancas e caixas-pretas, respectivamente.
Entre esses dois casos extremos, a modelagem caixa-cinza agrega parte de ambos os
tipos: usa a relação entre entrada(s) e saída(s) medidas e informações a priori. Segundo
(JORGENSEN, 1995) essa modelagem é definida como: “a ciência de construção de
modelos que incorpora conhecimento a priori do sistema com um certo grau de incerteza
na seleção da estrutura da representação”.
Em (SJÖBERG, 1995), a modelagem em caixa-cinza é dividida em:
• Modelagem física: quando o modelo é determinado pelo conhecimento da estrutura
do sistema e sua dinâmica física/química, sendo estimados apenas os parâmetros,
ou parte destes, com o uso de medições;
• Modelagem semi-física: quando as informações medidas são usadas para sugerir
combinações não-lineares entre os sinais medidos, sendo que a escolha da estrutura
do modelo é fruto da análise destas informações.
Nota-se que na modelagem física, o modelo parametrizado e a modelagem caixa-
branca têm definições semelhantes. Sendo assim, as demais definições da classificação
podem ser entendidas como modelagem caixa-cinza.
As vantagens da modelagem caixa-cinza são relacionadas em (CORRÊA, 2001), com
base em trabalhos de outros autores, enfatizando:18
• os prováveis benefícios no projeto de controladores avançados que necessitam da
descrição adequada dos processos;
• o uso de informações a priori, que diminui efetivamente o número de parâmetros
a serem estimados, tornando o problema melhor condicionado com modelos menos
incertos, mesmo com relativa escassez de dados;
• a possibilidade do uso da informação a priori, pode vir a se tornar um meio para
escolha diante das diversas possibilidades de representações matemáticas na mode-
lagem de sistemas dinâmicos não-lineares.
Em (GOLDFREY, 1986; AGUIRRE, 2000b,a; CORRÊA, 2001), é evidenciado que
trabalhos em caixa-preta se mostram bem caracterizados localmente, mesmo que a planta
física identificada seja não-linear. A incorporação de informações a priori possibilita que
o modelo tenha uma faixa de correspondência, em relação à planta, mais abrangente,
com pouca perda de precisão e melhora a medida que toda a faixa de trabalho é excitada
pela entrada. Por fim, é dado ênfase ao papel das informações a priori na seleção de
representações matemáticas para o processo de modelagem.
Um outro desafio deste estudo é usar técnicas de controle já consagradas ou aqui
desenvolvidas para controlar as plantas cujo modelos apresentem não-linearidades ou
sejam lineares a parâmetros variáveis.
Um sistema linear pode ser considerado variante no tempo ou não-estacionário −Linear a Parâmetros Variantes (LPV) ou Linear Variante no Tempo (LTV, do inglês
Linear Time Variant) − quando um ou mais parâmetros do seu modelo variam ampla-
mente com o tempo e não podem ser considerados como parâmetros incertos. Muitas
vezes, os sistemas não-lineares comumente encontrados na prática podem ser tratados
como sistemas lineares não-estacionários.
Tradicionalmente, o controle de sistemas não-estacionários é realizado, na prática,
pelo uso de técnicas de escalonamento de ganhos (ou interpolação de controladores2). O
principal objetivo desses métodos é controlar um sistema que evolui num amplo domínio
de funcionamento, para o qual as técnicas de controle robusto Linear Invariante no Tempo
(LTI, do inglês Linear Time Invariant) se mostram ineficazes. Além de permitirem a
incorporação de propriedades de robustez em estabilidade face às incertezas do sistema,
2O termo em inglês gain-scheduling costuma ser traduzido para a língua portuguesa como “tabela-
mento de ganho” ou “escalonamento de ganho”. Neste trabalho é adotada a expressão “escalonamento de
ganho”. Para a tradução da palavra scheduling, é adotada a palavra “interpolação” ou “escalonamento”.
19
os controladores interpolados possuem outra característica favorável importante que é
a adaptação em tempo real do seu comportamento dinâmico, segundo a evolução dos
parâmetros (endógenos ou exógenos) que caracterizam as condições de funcionamento
do sistema. Essas técnicas ampliam o alcance dos métodos clássicos de controle robusto
LTI, que consideram somente as características lineares locais e condições particulares de
funcionamento do sistema. Este benefício da estratégia de controle por escalonamento de
ganhos é uma conseqüência da explícita utilização de informações adicionais importantes
oriundas da medida dos parâmetros variantes.
O escalonamento de ganhos é uma técnica bastante eficaz no controle de sistemas
não-estacionários. Desde as primeiras publicações sobre o assunto, há cerca de trinta
anos atrás, essa técnica tem sido aplicada em várias áreas. Apesar dessa prática ter se
difundido em larga escala, só recentemente o interesse teórico sobre ganhos escalonados
para sistemas LPV e sistemas não-lineares têm aumentado significativamente (SHAMMA,
1990; SHAHRUZ, 1992; KAMLER, 1995).
O desenvolvimento de novos métodos de interpolação de controladores tem desper-
tado grande interesse na comunidade científica pela ampla diversidade de suas aplicações
e pelo surgimento de técnicas modernas de controle robusto. Contudo, os métodos mais
utilizados no meio industrial são aqueles chamamos de “clássicos”, “convencionais” ou
“tradicionais”. Eles se baseiam em um conjunto de modelos LTI obtidos a partir da fi-
xação do parâmetro variante de um modelo originalmente LPV ou LTV por linearização
de um modelo não-linear em torno de uma família de pontos de operação3. Um con-
junto de técnicas de controle linear (LQG, PRLQG, PRCBI, H2, H∞, síntese µ , etc.,
e suas variantes) está, então, disponível para o projeto de uma família de controladores
LTI que ofereçam um compromisso razoável entre o desempenho e a robustez em esta-
bilidade em torno das condições de funcionamento dadas. Quanto as estratégias para a
interpolação, elas variam bastante segundo o método de síntese linear e a estrutura dos
controladores LTI escolhidos e são, na maioria das vezes, intuitivas e baseadas em uma
diretiva heurística principal: os controladores lineares são supostos suficientemente próxi-
mos para permitir transições suaves e para assimilar os comportamentos não-estacionários
provenientes das não-linearidades do sistema.
A interpolação é baseada na medida de um parâmetro que define os pontos de ope-
ração do sistema acarretando o ajuste da FT (Função de Transferência) que modela o
3Ver (SHAMMA, 1990; RUGH, 1991; KELLET, 1991; REICHERT, 1992; HYDER, 1993; NICHOLS,
1993; KAMLER, 1995; LAWRENCE, 1995; STILWELL, 1997) .
20
sistema ou dos elementos das matrizes A,B,C e D do modelo em espaço de estado. Entre
as estratégias de interpolação mais utilizadas, encontram-se aquelas que atuam sobre os
coeficientes de funções de transferência (NICHOLS, 1993), ou sobre os coeficientes de
realizações de estado (KELLET, 1991; HYDER, 1993), ou ainda, quando se tratar de
um conjunto de controladores lineares robustos H2 ou H∞, soluções de equações de Ric-
cati (REICHERT, 1992). O controlador interpolado convencional é, então, um sistema
não-estacionário (LPV ou não-linear), obtido por interpolação (linear, fuzzy, ou outra
qualquer) de controladores LTI em relação às variáveis de interpolação.
A interpolação clássica de controladores apresenta dificuldades teóricas no que se
refere à estabilidade e ao desempenho durante as transições entre os controladores locais,
em parte devido à forte dependência do comportamento do sistema global em relação à
estratégia de interpolação utilizada. Além disso, as etapas de síntese dos controladores
e a formulação da lei de interpolação são conduzidas separadamente, o que não garante
que o sistema não-estacionário em malha fechada seja estável e que o desempenho apre-
sentado em torno dos pontos de projeto prevaleçam de uma forma global. Simulações
exaustivas, inclusive simulações hardware-in-the-loop, são em geral efetuadas para avaliar
o comportamento do sistema controlado em regime não-estacionário.
Devido a estas dificuldades, poucos resultados teóricos em escalonamento de ganhos
convencionais têm sido publicados. Em (SHAHRUZ, 1992), um algoritmo para interpo-
lação linear de ganhos por realimentação de estados é proposto. Em (STILWELL, 2000),
os autores apresentam dois métodos para interpolação de controladores por escalona-
mento de ganhos: o primeiro apresenta fatores coprimos de funções de transferências e o
segundo é baseado na descrição em espaço de estado. Em (STILWELL, 1999), os ganhos
de realimentação e de observação de estados são interpolados no contexto do escalona-
mento de ganhos. Essas técnicas garantem, com certo grau de conservadorismo e sob
restrição de variação lenta do parâmetro, a estabilidade exponencial local do sistema em
malha fechada.
Atualmente, a interpolação de controladores é também tratada no contexto do con-
trole LPV (LU, 1992; BECKER, 1993; PACKARD, 1994; LU, 1995; BECKER, 1995;
APKARIAN, 1995a,b; WU, 1996; SCHERER, 1996; APKARIAN, 1998a,b; KAJIWARA,
1999; TUAN, 1999), onde funções de Lyapunov são utilizadas para definir a estabilidade
e o desempenho para uma larga faixa de variação do parâmetro. Na verdade os ter-
mos “controle LPV” ou “interpolação LPV” são utilizados para designar as técnicas de
interpolação de controladores onde as necessidades práticas da interpolação, ou seja, a
21
estabilidade e o desempenho, são satisfeitos de uma forma sistemática dentro de um con-
texto de otimização convexa sob restrições do tipo Desigualdade Matricial Linear (LMI
- do inglês Linear Matrix Inequality).
Os métodos de interpolação LPV distinguem-se das técnicas de interpolação, mais
pela sua maneira sistemática de tratar o problema, do que pelo seu objetivo. Do ponto
de vista conceitual, a interpolação do tipo LPV é, no entanto, muito diferente, já que as
questões relacionadas à estabilidade e ao desempenho em tempo variante são considera-
das diretamente na síntese dos controladores LPV. A tarefa mais exigente consiste em
resolver problemas de otimização do tipo LMI. Isso é relativamente fácil utilizando-se os
códigos de programação semi-definida disponíveis atualmente. Trata-se, na realidade, de
uma extensão dos métodos de controle robusto LTI, do tipo H2 e H∞, aos sistemas não-
estacionários LPV ou quasi -LPV. Em suma, as técnicas LPV são igualmente aplicáveis
a modelos por natureza LTV ou LPV, a modelos linearizados (parametrizados por va-
riáveis de interpolação) ou a modelos quasi -LPV. Estas vantagens explicam o repentino
interesse nos últimos anos por esse tipo de técnica (BECKER, 1993; PACKARD, 1994;
APKARIAN, 1995a; SCHERER, 1996; KAJIWARA, 1999).
Apesar dos esforços e desenvolvimentos teóricos recentes, alguns problemas delicados
persistem, não apenas na interpolação convencional, mas também na interpolação LPV.
Percebe-se que, entre as raras técnicas clássicas que garantem teoricamente a estabilidade
não-estacionária, uma boa parte é essencialmente baseada na interpolação de matrizes
de estado ou de ganhos de estruturas estimação-controle (observer-based structures) de
controladores LTI ótimos robustos ou de controladores ótimos LQ (Quadrático Linear)
ou LQG (Gaussiano Quadrático Linear). Uma dificuldade importante reside no fato
que o comportamento dinâmico dos controladores interpolados depende de uma forma
crítica das representações de estado adotadas para a família de controladores lineares
projetados sobre um conjunto de pontos de operação. Assim, as primeiras questões que
se apresentam são:
• Como escolher um conjunto de bases no espaço de estado que conduza a realizações
de controladores locais apropriados para a interpolação?
• Qual o melhor conjunto de pontos de operação sob a ótica da interpolação?
Em relação às técnicas LPV, algumas são potencialmente muito conservadoras, pois
utilizam funções de Lyapunov independentes dos parâmetros do sistema, toleram taxas
22
arbitrárias de variação desses parâmetros e exigem classes específicas e restritivas de repre-
sentações LPV (LU, 1992; BECKER, 1993; APKARIAN, 1995a,b; LU, 1995; SCHERER,
1996; APKARIAN, 2000; PELLANDA, 2002a). Outras, ao contrário, utilizam funções de
Lyapunov dependentes dos parâmetros, consideram limites realistas da taxa de variação
dos parâmetros e toleram uma dependência paramétrica geral do sistema (BECKER,
1995; WU, 1996; APKARIAN, 1998a,b; TUAN, 1999). Estas são muito pouco conser-
vadoras, mas introduzem uma grande complexidade na implementação do controlador
interpolado.
Esta dissertação aborda o estudo, o aprimoramento e o desenvolvimento de técnicas
de escalonamento de ganhos clássicas, com uma orientação particular para o caso da
realimentação de estados e da estrutura estimação-controle, aplicáveis em diversas áreas
no campo da engenharia de controle. Um enfoque especial é também dado ao desen-
volvimento de novos algoritmos de identificação do tipo caixa-preta aplicáveis a sistemas
lineares e do tipo caixa-cinza aplicáveis a sistemas não-lineares.
1.2 OBJETIVOS DA DISSERTAÇÃO
Os objetivos deste trabalho são:
a) Desenvolver novos algoritmos de identificação dos tipos caixa-preta e caixa-cinza,
para sistemas lineares e não-lineares (quasi -LPV), baseados no emprego de amostras
discretas dos sinais de entrada e saída da planta física nos domínios do tempo e da
freqüência;
b) Aprimorar a metodologia de escalonamento de ganhos clássico desenvolvido em
(PELLANDA, 2001) de modo a equipá-lo com um método de escolha sistemática
de um conjunto de pontos de síntese, a fim de aproximar o modelo não-linear da
planta em toda uma faixa de operação;
c) Desenvolver novas metodologias clássicas de escalonamento de ganhos baseadas na
realimentação de estados e na estrutura estimação-controle para algumas classes de
sistemas não-estacionários.
1.3 ORGANIZAÇÃO DA DISSERTAÇÃO
Além desta introdução, a dissertação está organizada em mais 6 capítulos, além de
um apêndice:
23
• O Capítulo 2 recorda os principais fundamentos que permitem entender o desen-
volvimento das contribuições desta dissertação.
• O Capítulo 3 apresenta as contribuições deste trabalho para a identificação de
sistemas. São apresentadas três novas metodologias baseadas na buscas de um
modelo que minimize um custo quadrático representativo do erro de estimação.
Duas delas são apropriadas para a identificação de modelos LTI e a outra é aplicável
a uma classe de sistemas não-lineares, fornecendo modelos não-estacionários do
tipo quasi -LPV. Todas as técnicas desenvolvidas apresentam solução analítica ao
problema de minimização do erro de estimação.
• O Capítulo 4 apresenta o processo de escalonamento de ganhos desenvolvido em
(PELLANDA, 2001) e propõe um método de escolha sistemática de um conjunto
de pontos de síntese que melhor aproxime o modelo não-linear da planta em toda
uma faixa de operação sendo, portanto, mais apropriado para o escalonamento de
ganhos. Apresenta também duas novas metodologias clássicas de escalonamento de
ganhos baseadas na realimentação de estados e na estrutura estimação-controle.
• No Capítulo 5 são ilustradas as técnicas apresentadas nos Capítulos 3 e 4 através
de diversas aplicações numéricas.
• O Capítulo 6 traz as conclusões tiradas a partir dos resultados obtidos ao longo
deste trabalho e aponta algumas perspectivas futuras.
• No Capítulo 7 são apresentadas as referências bibliográficas.
• No Apêndice encontram-se os modelos utilizados e identificados no Capítulo 5.
24
2 FUNDAMENTOS TEÓRICOS
2.1 INTRODUÇÃO
Neste capítulo são revistos alguns fundamentos teóricos e propriedades necessárias
para o entendimento dos trabalhos apresentados nos capítulos posteriores.
A aplicação das técnicas de identificação linear proposta nesta dissertação em plan-
tas físicas estáveis podem gerar modelos que apresentem pólos instáveis. Nestes casos,
também são gerados zeros localizados no plano s sobre esses pólos ou muito próximos
deles, acarretando um modelo de ordem não-mínima. Porém, como é sabido que as
funções identificadas são estáveis, torna-se interessante decompor o modelo identificado
na soma de suas partes estável e antiestável (ou instável), ignorando esta última e apli-
cando à primeira uma redução de ordem por truncamento balanceado. Assim, elimina-se
o problema da ordem não-mínima causado pela parte instável e se obtém um modelo
final estável de ordem reduzida. As Seções 2.2 e 2.3 apresentam o suporte teórico dessas
operações.
Na Seção 2.4 discute-se os conceitos sobre sistemas LPV e quasi -LPV, úteis para
o entendimento das técnicas propostas na Seção 3.4 e no Capítulo 4, que tratam, res-
pectivamente, da identificação e do controle por escalonamento de ganhos de sistemas
não-estacionários.
Por fim, a Seção 2.5 apresenta uma metodologia de interpolação clássica de contro-
ladores dinâmicos predeterminados sob a forma estimação-controle (PELLANDA, 2000,
2001, 2002b), onde os ganhos de realimentação e de estimação de estado são os parâme-
tros interpolados. Contudo, sabe-se que uma dificuldade da interpolação clássica é o fato
do comportamento dinâmico dos controladores interpolados depender fortemente das re-
presentações de estado adotadas para a família de controladores LTI projetados sobre o
conjunto de pontos de operação escolhido. A técnica apresentada uma seqüência de trans-
formações de similaridade a serem aplicadas à família original de controladores robustos
determinados a priori, de forma que o conjunto de controladores obtidos, equivalentes aos
primeiros do ponto de vista entrada-saída, tenha uma estrutura coerente, fisicamente in-
terpretável e de fácil interpolação. Depois da transformação, os controladores apresentam
uma estrutura do tipo estimação-controle, o que facilita a interpolação e a implementação
em tempo variante. O método pode ser aplicado em controladores discretos ou contínuos,
25
de ordem completa ou aumentada e, em particular, aos projetados através de técnicas de
controle robusto H2, H∞ e síntese µ, cuja interpolação é, em geral, de difícil obtenção
sem perdas consideráveis de desempenho.
2.2 PROJEÇÃO ESTÁVEL E INSTÁVEL DE UMA REALIZAÇÃO
Seja G(s) o modelo obtido, descrito pela seguinte realização de estado
G(s) =
A B
C D
(2.1)
onde A ∈ Rn×n, B ∈ Rn×p, C ∈ Rq×n e D ∈ Rq×p. Este modelo pode ser reescrito como
G(s) = G(s)+ + G(s)− (2.2)
onde G(s)+ e G(s)− correspondem, respectivamente, às projeções de G(s) nos semi-planos
da direita e esquerda do plano s.
Definindo as matrizes destas realizações, G(s)+ e G(s)−, como
G(s)− =
A11 B1
C1 D
(2.3)
e
G(s)+ =
A22 B2
C2 0
(2.4)
pode-se afirmar que:
G(s) = G(s)+ + G(s)− =
A11 0
0 A22
B1
B2
C1 C2 D
=
A B
C D
(2.5)
onde A11 ∈ Rk×k e as demais matrizes são de dimensões compatíveis sendo k o número
de autovalores estáveis.
Definindo P ∈ Rn×n como a matriz de transformação de similaridade tal que
P−1AP = A, P−1B = B, CP = C (2.6)
o problema de encontrar a decomposição desejada é resolvido ao se calcular P .
26
Aplicando uma decomposição de Schur na matriz A
[A11 A12
0 A22
]= UT AU = T (2.7)
onde UT U = I e T é a matriz triangular superior que possui os autovalores da matriz
A em sua diagonal com estes dispostos em ordem crescente, e sendo X uma matriz de
dimensões compatíveis tal que
A =
[I1 −X
0 I2
][A11 A12
0 A22
][I1 X
0 I2
]=
[A11 A12 −XA22
0 A22
][I1 X
0 I2
]
A =
[A11 A11X −XA22 + A12
0 A22
] (2.8)
chega-se a:
A11X −XA22 + A12 = 0 (2.9)
O problema se reduz, ao cálculo de uma solução X para a equação de Lyapunov em (2.9).
Adotando-se a seguinte partição: U =[
U1 U2
], de forma que U1 ∈ Rk×k, pode-se
reescrever o produto matricial P−1AP na forma[
I1 −X
0 I2
][UT
1
UT2
]A
[U1 U2
] [I1 X
0 I2
]=
[UT
1 −XUT2
UT2
]A
[U1 U1X + U2
]= P−1AP
(2.10)
Assim, tem-se também:[
B1
B2
]=
[UT
1 −XUT2
UT2
]B =
[(UT
1 −XUT2 )B
UT2 B
](2.11)
e [C1 C2
]= C
[U1 U1X + U2
]=
[CU1 C(U1X + U2)
](2.12)
o que define completamente a projeção estável G(s)− do modelo inicialmente identificado
G(s).
2.3 REDUÇÃO DE ORDEM POR TRUNCAMENTO BALANCEADO
Seja Gke(s), a parte estável, de ordem k, do modelo identificado G(s). O ajuste do
modelo Gke(s) aos dados fornecidos pelo sistema real pode acarretar uma ordem maior
27
que a necessária (k∗ = k−r). Uma função estimada não deve possuir uma ordem tão ele-
vada a ponto de acrescentar redundâncias, nem tão baixa a ponto de perder informações
relevantes. A resposta em freqüência pode ser utilizada para comparar o grau de proxi-
midade entre o sistema real e o modelo obtido, permitindo avaliar se a ordem do modelo
é compatível. Caso um modelo tenha sido encontrado com um valor de k elevado, existe
a possibilidade de se aplicar uma técnica para redução de sua ordem. Ou seja, pode se
encontrar um novo modelo de ordem reduzida, Gk−re (s), que preserve o ajuste anterior
de acordo com algum critério. O modelo final Gk−re (s) deve ser parcimonioso.
De acordo com (SKOGESTAD, 1997), uma realização balanceada é uma realização
assintoticamente estável, na qual os Gramianos de controlabilidade e observabilidade são
iguais e também são diagonais. Considerando a seguinte realização estável de G(s):
G(s) =
A B
C D
(2.13)
A realização em (2.13) é denominada balanceada se as soluções para as seguintes equações
de Lyapunov:
AP + PAT + BBT = 0
AT Q + QA + CT C = 0(2.14)
são tais que: P = Q = diag(σ1, σ2, . . . , σk) , Σ, onde σ1 ≥ σ2 ≥ σ3 ≥ . . . ≥ σk ≥ 0. Os
Gramianos de controlabilidade e observabilidade, respectivamente P e Q são definidos
por:
P ,∫∞0
eAtBBT eAT tdt
Q ,∫∞0
eAT tCT CeAtdt(2.15)
Neste caso∑
é dito Gramiano de G(s). Os σi são denominados de valores singulares de
Hankel de G(s), sendo definidos por (2.16):
σi , λ12i (PQ) (2.16)
onde {λi ; i = 1, . . . , k} são os autovalores de A. Cada valor σi está associado a um
estado xi da realização balanceada do sistema. O valor σi se constitui numa medida
da contribuição que o estado xi tem no comportamento de entrada e saída do sistema.
Em resumo, se σi >> σi+1, então o estado xi afeta o comportamento de entrada e saída
para o sistema muito mais do que xi+1, permitindo calcular uma estimativa de erro
em um modelo reduzido onde os r estados {xi+1, . . . , xk} tenham sido descartados. A
implementação para obtenção de modelos de ordem reduzida foi feita a partir da função
balmr do Matlab.28
2.4 SISTEMAS LPV E QUASI -LPV
2.4.1 MODELAGEM
Uma classe importante de sistemas dinâmicos não-estacionários pode ser represen-
tada por um conjunto de equações diferenciais não-lineares de ordem qualquer. Para uma
escolha apropriada dos vetores das variáveis de estado x(t) ∈ Rn, de entrada u(t) ∈ Rm
e de saída y(t) ∈ Rp, pode-se freqüentemente obter um modelo não-linear em relação
aos estados, mas linear em relação à entrada, que implique em uma equação matricial
diferencial de primeira ordem e uma equação matricial algébrica:
x(t) = A(θx, θp)x + B(θx, θp)u
y(t) = C(θx, θp)x + D(θx, θp)u(2.17)
As funções matriciais reais A(.), B(.), C(.) e D(.) são supostas contínuas e limi-
tadas, de dimensões compatíveis com as dimensões dos sinais e definem completamente a
dinâmica do sistema. Este tem uma característica não-linear e não-estacionária originada
pelas variáveis θx e θp:
• θx(x(t)) ∈ Rr1 é uma variável endógena, ou seja, que depende da dinâmica interna
do sistema e que o torna não-linear;
• θp(x(t)) ∈ Rr2 é um parâmetro exógeno, ou seja, que evolui no tempo de forma
independente da dinâmica interna do sistema.
Na síntese de controladores dentro da técnica de escalonamento clássico de ganhos,
a primeira etapa corresponde a obtenção de uma descrição linear aproximada do sistema
não-linear (2.17) que envolve um conjunto conveniente das variáveis de interpolação θ(t).
A maneira mais utilizada na prática consiste em:
• obter, via uma linearização Jacobiana clássica do modelo (2.17) em torno de um
conjunto de pontos de equilíbrio x(i)0 (u
(i)0 ), i = 1, 2, ..., um modelo linearizado:
x = A(θi)x + B(θi)u
y = C(θi)x + D(θi)u
parametrizado por
θi(t) =
[θ(x
(i)0 )
θp(t)
]∈ Rr, r = r1 + r2;
29
• definir uma trajetória nominal de x0(t) para o sistema e, supondo que θx(t) e
dθx(t)/dt são limitadas e independentes de x0(t) e dx0(t)/dt, derivar um modelo do
tipo LPVx = A(θ)x + B(θ)u
y = C(θ)x + D(θ)u(2.18)
onde o parâmetro e sua taxa de variação evoluem em domínios compactos, θ(t) ∈DΘ ⊂ Rr, θ(t) ∈ DΘd
⊂ Rr,∀t;
• eventualmente, escolher uma trajetória θ(t) ←− θ0(t) ou “congelar” o parâmetro
em um ponto dado θ(t) ←− θ0(t), para obter, respectivamente, um modelo LTV ou
LTI.
Nota-se que apesar das funções matriciais A(.), B(.), ... terem representações análo-
gas, em geral são diferentes daquelas que constam em (2.17).
Uma outra forma mais recente e mais direta de se chegar a um modelo similar àquele
de (2.18), a partir (2.17), consiste simplesmente em ignorar a etapa de linearização.
Escolhendo convenientemente a função θx(x(t)), reescreve-se o modelo numa forma onde
os termos não-lineares possam ser redefinido por um parâmetro variante unicamente em
função do tempo θx(t). De forma similar ao caso anterior, considera-se que as trajetórias
desse parâmetro são limitadas e independentes das trajetórias de x(t), o que desconecta
as funções matriciais A(.), B(.), ... do espaço de estados. Ele é então incluído na variável
de interpolação, juntamente com o parâmetro θp(t):
θ =
[θx(t)
θp(t)
](2.19)
Isso quer dizer que certos estados, ou funções dos estados, são classificados como
variáveis exógenas em certas partes do modelo, enquanto que em outras permanecem
como variáveis endógenas. Esta hipótese leva a um certo conservadorismo, mais ou menos
importante, na etapa de síntese dos controladores. Nesse caso particular, o modelo (2.18)
é denominado quasi -LPV.
Se uma lei de interpolação é satisfatória para todas as trajetórias no domínio DΘ ×DΘd
, ela é igualmente satisfatória para as trajetórias realistas dos estados que interferem
em θ. Contudo, o conservadorismo introduzido pela modelagem quasi -LPV é tão menos
desprezível quanto maior o número de estados implicados no parâmetro. Considerando,
por exemplo, o sistema não-linear
x1 = sen(x1) + x2, x2 = x1x2 + u
30
Uma representação quasi -LPV é
x = A(x)x + Bu =
[sen(x1)/x1 1
x2 0
]x +
[0
1
]u
com x1 6= 0 e θx(x) = x , [x1 x2]T . Esta representação é certamente mais conservadora
do que a seguinte:
x = A(x)x + Bu =
[sen(x1)/x1 1
0 x1
]x +
[0
1
]u
onde o parâmetro θx(x) é inteiramente definido por somente uma variável de estado
(θx(x) = x1) e, em conseqüência, a dimensão do espaço que pode incluir as trajetórias
não-realistas é menor. Enfim, um sistema pode ainda ser, pela sua própria natureza, LPV
ou LTV e nenhuma aproximação ou linearização suplementar é necessária para construir
o modelo (2.18).
O modelo LPV ou quasi -LPV (2.18) tolera uma dependência paramétrica bastante
geral, que engloba a maior parte das situações práticas. Esta propriedade requer a uti-
lização e o desenvolvimento de metodologias sofisticadas e complexas de análise e síntese
de leis de controle por escalonamento de ganho. No entanto, duas outras classes mais
restritivas de modelos LPV, obtidas da forma mais geral (2.18), são as vezes admissíveis
e mais adaptadas a certos métodos específicos de controle LPV. Embora estes modelos
não sejam diretamente utilizados neste trabalho, eles são apresentados a seguir.
Uma dessas classes refere-se aos modelos do tipo Transformação Linear Fracionária
(LFT, do inglês Linear Fractional Transformation). Uma dependência LFT do sistema
(2.18) em relação ao parâmetro θ(t) é definido como[
A(θ) B(θ)
C(θ) D(θ)
],
[A B2
C2 D22
]+
[Bθ
D2θ
]∆(θ)(I −Dθθ∆(θ))−1
[Cθ Dθ2
](2.20)
onde ∆(θ) é uma função matricial linear em θ. Este é suposto pertencente a um domínio
politópico
PΘ , co{Θ1, Θ2, ..., ΘL} (2.21)
onde Θi com i = 1, 2, ..., L designam os vértices do politopo PΘ. Tais modelos exigem,
em geral, hipóteses que simplifiquem sua obtenção, o que pode levar a um certo conser-
vadorismo.
A outra classe de representação LPV é a dos modelos politópicos. Se a função
matricial
31
S(θ) ,[
A(θ) B(θ)
C(θ) D(θ)
]
é afim em θ = [θ1, θ2, ..., θr]T , isto é,
S(θ) = S0 +∑r
l=1 θlSl
e as componentes escalares θl, l = 1, 2, ..., r, evoluem (por hipótese) independentemente
num domínio limitado DΘ que é um subconjunto do domínio politópico (2.21), DΘ ⊆ PΘ,
então o modelo admite um representação politópica
S(θ) = co{S1, ..., SL} = co{S(Θ1), ..., S(ΘL)} (2.22)
Este modelo não carrega nenhum conservadorismo se a condiçãoDΘ = PΘ é satisfeita.
Ao contrário, se PΘ representa um recobrimento politópico de DΘ o modelo é fortemente
conservador. Nenhuma atenção especial é dedicada à este tipo de representação neste
estudo.
2.4.2 ANÁLISE DE ESTABILIDADE
O critério de estabilidade por meio da localização dos autovalores da matriz de
dinâmica A para sistemas LTI não é suficiente para sistemas LPV. De fato, pode-se
mostrar, por exemplos simples, que mesmo quando todos os componentes da família de
sistemas LTI obtidos pelo congelamento do parâmetro são estáveis, ainda assim o sistema
LPV pode ser instável. Ou seja, para o sistema LPV
x(t) = A(θ(t))x(t) + B(θ(t))u(t)
y(t) = C(θ(t))x(t) + D(θ(t))u(t), ∀θ ∈ PΘ
(2.23)
onde A ∈ Rn×n, D ∈ Rp×m e PΘ é o domínio de variação paramétrica, a condição
Re(λi(A(θ))) < 0, i = 1, ...n, ∀θ ∈ PΘ (2.24)
não é suficiente para que se infira a estabilidade. Apesar disso, existe uma noção intuitiva
de que se a variação do parâmetro for suficientemente lenta, a estabilidade do conjunto de
sistemas LTI implica a estabilidade do sistema LPV. Nesse contexto se insere o teorema
a seguir.
Teorema 2.1. (DESOER, 1975) O sistema LPV 2.23 é globalmente estável se as
condições seguintes forem válidas:32
• Re(λi(A(θ))) < 0, i = 1, ...n, ∀θ ∈ PΘ,
• ∃α > 0 suficientemente pequeno tal que ‖ ddt
θ(t)‖ < α, ∀t ≥ 0, ∀θ ∈ PΘ.
Ainda que o Teorema 2.1 mostre que, sob condição de variação lenta do parâmetro,
a estabilidade do sistema LPV pode ser determinada a partir da estabilidade LTI, ele
não indica qual é o valor α que caracteriza a velocidade dessa variação. Assim, do ponto
de vista prático, não se costuma inferir sobre a estabilidade de um sistema LPV pela
estabilidade LTI. Faz-se, então, necessária a utilização de uma teoria mais abrangente,
como a teoria de estabilidade de Lyapunov (VIDYASAGAR, 1978; KHALIL, 1996). A
teoria de Lyapunov é bastante geral, de maneira que são abordados nesta seção apenas
alguns pontos básicos mais interessantes para o contexto deste estudo.
Seja o sistema autônomo geral
x = f(x, t), x ∈ Rn, t ≥ t0
x(t0) = x0
(2.25)
Supõe-se que sejam válidas as hipóteses de existência e unicidade para a solução do
sistema de equações diferenciais que regem o sistema (2.25). Além disso, supõe-se que
x0 = 0 seja um ponto de equilíbrio para o sistema4, ou seja, ∀t ≥ t0, f(0, t) = 0
Definição 2.1 (Estabilidade Assintótica Uniforme (Global)). A trajetória de equi-
líbrio x ≡ 0 é dita uniformemente assintoticamente estável se:
• ∀ε > 0, ∃δ(ε) tal que t ≥ t0 ≥ 0, ‖x(t0)‖ < δ ⇒ ‖x(t)‖ < ε,
• ∃c > 0 | ∀‖x(t0)‖ < c, x(t) → 0 quando t →∞, sendo que
∀ε > 0,∃T (c) | ∀t > t0 + T , ‖x(t)‖ < ε.
Se essas condições forem verificadas para c → ∞, então a trajetória de equilíbrio x ≡ 0
é dita globalmente uniformemente assintoticamente estável.
Definição 2.2 (Estabilidade Exponencial (Global)). A trajetória de equilíbrio x ≡ 0
é dita (globalmente) exponencialmente estável se:
• as condições de estabilidade assintótica uniforme (global) forem verificadas,
• ∃λ > 0 tal que ∀‖x(t0)‖ < c, ∃M(x(t0)) tal que t ≥ t0 ⇒ ‖x(t)‖ ≤ Me−λt
O teorema a seguir é central na teoria de Lyapunov.
4Essa hipótese não é restritiva, sendo satisfeita por simples translação (VIDYASAGAR, 1978, p.132).
33
Teorema 2.2. (VIDYASAGAR, 1978) O sistema (2.25) é globalmente uniformemente
assintoticamente estável se existe uma função continuamente diferenciável, V (t, x),
definida sobre R+ × Rn com valores em R, tal que, para todo x ∈ Rn e todo t positivo:
α1(‖x‖) ≤ V (t, x) ≤ α2(‖x‖) (2.26)
∂V
∂t+
∂V
∂xf(t, x) ≤ −α3(‖x‖) (2.27)
onde α1, α2 e α3 são funções definidas sobre R+ com valores em R+, contínuas, estrita-
mente crescentes, não limitadas e nulas na origem.
O primeiro membro de (2.27), que representa a derivada temporal de V , é chamado
de derivada de V ao longo das trajetórias do sistema (2.25) e é definido como
V (t, x) :=dV
dt=
∂V
∂t+
∂V
∂x
dx
dt=
∂V
∂t+
∂V
∂xf (2.28)
Pode-se notar que, pelo Teorema 2.2, V (t, 0) = 0 para todo t > 0.
O teorema a seguir é uma variação do Teorema 2.2. Ele faz uma extensão para funções
V que sejam continuamente diferenciáveis por partes. Por outro lado, ele particulariza
as funções αi para
αi(‖x‖) = λi‖x‖c, λi > 0, c > 0 (2.29)
Além disso, V passa a ser função apenas de x.
Teorema 2.3. (PETTERSSON, 1997) O sistema (2.25) é globalmente exponencialmente
estável se existe uma função continuamente diferenciável por partes, V (x), definida sobre
Rn com valores em R, tal que, para todo x ∈ Rn e todo t positivo:
λ1‖x‖c ≤ V (x) ≤ λ2‖x‖c (2.30)
V (x) ≤ −λ3‖x‖c (2.31)
com λ1, λ2, λ3 estritamente positivas e reais.
A função V que atende as condições dos Teoremas 2.2 e 2.3 é chamada de função
de Lyapunov. Essa função pode ser vista como uma medida generalizada da energia do
sistema. Existe uma idéia intuitiva por trás dos Teoremas 2.2 e 2.3 de que o sistema
será estável se a energia do sistema (ou a função V ) for decrescente ao longo de todas as
trajetórias desse sistema. Embora não esteja explicito, a estabilidade exponencial global
garante a estabilidade de um ponto de vista entrada/saída (KHALIL, 1996).
34
Vale destacar que, embora o Teorema 2.3 aponte para a estabilidade do sistema (2.25)
caso seja possível encontrar uma função de Lyapunov para esse sistema, ele não indica
como fazê-lo.
Seja, agora, a particularização do sistema autônomo (2.25) para o caso LPV
x = A(θ)x, x ∈ Rn, t ≥ t0
x(t0) = x0
(2.32)
Para analisar a estabilidade desse sistema, será utilizada uma função de Lyapunov V (θ, x)
quadrática
V (θ, x) = xT P (θ)x (2.33)
com P (θ) simétrica. Isso porque, para sistemas lineares, a estabilidade assintótica uni-
forme pode ser inferida, por meio do Teorema 2.3, utilizando-se uma função de Lya-
punov quadrática (VIDYASAGAR, 1978, p.183). Uma outra característica importante
dos sistemas lineares é que a estabilidade assintótica uniforme é equivalente à estabilidade
exponencial (VIDYASAGAR, 1978, p.170).
Para a função de Lyapunov dependente do parâmetro, V (θ, x), a derivada ao longo
das trajetórias é obtida por
V (θ, x) =∂V
∂θ
dθ
dt+
∂V
∂x
dx
dt(2.34)
Segundo o Teorema 2.3, para que V seja uma função de Lyapunov e, conseqüente-
mente, o sistema seja estável, ela deve satisfazer
V (θ, x) > 0, V (θ, x) < 0, x 6= 0
V (θ, 0) = 0 e V (θ, 0) = 0(2.35)
A primeira restrição é satisfeita quando a matriz P (θ) é positiva definida
P (θ) = P (θ)T > 0,∀θ ∈ PΘ (2.36)
A segunda implica, de acordo com (2.34), em
xT ∂P (θ)∂θ
dθdt
x + xT P (θ)(A(θ)x) + (xT AT (θ))P (θ)x =
xT [∂P (θ)∂θ
dθdt
+ P (θ)A(θ) + AT (θ)P (θ)]x < 0(2.37)
que é satisfeito quando
∂P (θ)
∂θ
dθ
dt+ P (θ)A(θ) + AT (θ)P (θ) < 0, ∀θ ∈ PΘ (2.38)
35
Assim, se for possível determinar uma matriz P dependente de θ e positiva definida
tal que a inequação matricial parametrizada (2.38) seja satisfeita, então V (θ, x) será uma
função de Lyapunov e o sistema (2.32) será estável.
A inequação matricial (2.38) é uma Desigualdade Linear Matricial Parametrizada,
ou PLMI (do inglês Parametrized Linear Matrix Inequality). Como o parâmetro θ é uma
função contínua do tempo, a PLMI tem que ser satisfeita para todos os infinitos valores
de θ.
Para sistemas LTI basta considerar a matriz P constante, de forma que o problema
de análise de estabilidade se reduz a determinar uma matriz P = P T > 0 tal que
PA + AT P < 0 (2.39)
A LMI (2.39) é condição necessária e suficiente para que a parte real dos pólos do
sistema LTI seja negativa (VIDYASAGAR, 1978; KAILATH, 1980).
2.5 INTERPOLAÇÃO DE CONTROLADORES ROBUSTOS DE ESTRUTURA
ESTIMAÇÃO-CONTROLE
No contexto do escalonamento de ganhos clássico, se os controladores lineares a serem
interpolados forem muito diferentes em pontos de operação vizinhos, uma variação rápida
é introduzida artificialmente na dinâmica de malha fechada, o que pode produzir um efeito
desestabilizante ou uma perda de desempenho. Então, para se obter um comportamento
desejável na interpolação de controladores lineares, estes devem obrigatoriamente ter
estruturas compatíveis.
Isto se torna crítico quando os dados são interpolados numa abordagem em espaço de
estado, pois o comportamento dinâmico do controlador escalonado pode depender forte-
mente das realizações adotadas para a família de controladores lineares projetada para
o conjunto escolhido de pontos de operação. Este fato é ilustrado em (LEITH, 1998a;
STILWELL, 1999, 2000). Em (STILWELL, 1999, 2000) os autores apresentam justifica-
tivas teóricas e condições suficientes para a alocação de controladores LTI de modo que
a estabilidade a tempo variante sempre exista5. Também apresentam condições (conser-
vadoras) para o cálculo de um limitante superior para a taxa de variação do parâmetro
de escalonamento, de forma que a estabilidade a tempo variante esteja assegurada. In-
felizmente, no contexto de interpolação em espaço de estado, esses métodos são restritos
a controladores de ordem completa, ou seja, de mesma ordem da planta. Além disso, as
5Ver também (?STILWELL, 1997) e os exemplos numéricos apresentados.
36
condições de estabilidade e desempenho são verificadas a posteriori, o que não elimina o
algo grau de empirismo na escolha dos pontos de operação.
Esses resultados mostram claramente que uma transição satisfatória depende não
somente da “distância” entre os pontos de operação, mas também da “proximidade” en-
tre os coeficientes dos respectivos controladores LTI. Deve-se notar que a interpolação
desses coeficientes representa também uma variação paramétrica para o sistema em malha
fechada e quanto maior a diferença entre eles para pontos vizinhos, maior será a sua taxa
de variação, o que prejudica a estabilidade a tempo variante.
Mesmo supondo que o conjunto de pontos de operação seja apropriadamente esco-
lhido, projetar um conjunto correspondente de controladores de modo a favorecer uma
interpolação suave, não é uma questão de fácil solução no contexto das técnicas clássi-
cas de escalonamento de ganhos, especialmente quando se trata de técnicas de controle
robusto dos tipos H2, H∞ e síntese µ. No caso particular da estrutura controlador-
observador, a interpolação dos ganhos de controle e de estimação é mais intuitiva, apesar
desses ganhos não serem as únicas variáveis a serem interpoladas. Os coeficientes do
controlador também dependem dos dados da realização do modelo do sistema a controlar
e devem evoluir de forma consistente com a sua dinâmica. Isto é, variações significativas
ou não-linearidades da planta devem ser adequadamente compensadas pelo ajuste do
controlador.
Nesta seção, são apresentadas as técnicas propostas em (PELLANDA, 2000, 2001,
2002b) para o cálculo de realizações equivalentes do tipo estimação-controle de um con-
junto de controladores estabilizantes arbitrários associados a um dado conjunto de pontos
de operação de um sistema. A família resultante de controladores de estrutura estimação-
controle pode ser interpolada mais facilmente. Apesar dessas técnicas serem gerais, de
modo a simplificar a apresentação, os casos dos controladores não estritamente próprios
e a tempo discreto são aqui omitidos, mas também a estes são aplicáveis tal metodologia.
Na literatura de controle pode-se encontrar várias expressões freqüentemente uti-
lizadas como termos equivalentes para “controlador de estrutura estimação-controle”:
observer-based structure, observer state (feedback) controller, LQG (form) controller,
estimator-controller structure (or form), observer-controller structure, state estimator-
state feedback structure, entre outras.
37
2.5.1 FORMULAÇÃO DO PROBLEMA
A definição a seguir é necessária para a compreensão do problema a ser proposto
nesta seção.
Definição 2.3 (Controlador de estrutura estimação-controle). (ZHOU, 1996)
Dada uma realização (A,B,C) de um sistema qualquer, com (A,B) estabilizável e (C,A)
detectável, então existem matrizes Kc e Kf tais que os autovalores de A−BKc e A−KfC
são estáveis. O controlador definido por u = K(s)y, onde
K(s) =
A−BKc −KfC Kf
−Kc 0
(2.40)
é chamado de “controlador de estrutura estimação-controle”.
Considere os sistemas em malha fechada descritos na FIG. 2.1, onde
G(s) =
[G11(s) G12(s)
G21(s) G22(s)
](2.41)
e
P (s) =
[P11(s) P12(s)
P21(s) P22(s)
](2.42)
são, respectivamente, o modelo de síntese nominal e aumentado. O sinal r é o sinal dereferência e os sinais w, e, z e y são, respectivamente, a entrada exógena, o controle, asaída controlada e a saída medida de P (s).
+ ( ) s P
( ) s G ( ) s W o
( ) s W i
( ) s K
w
r e y
z
u
+ ( ) s P
( ) s J
w r e y
z
u
( ) s Q
1 y 1 u
( ) s K e
(a) Controlador original (b) Controlador equivalent (a) - Controlador original (b) - Controlador equivalente
FIG. 2.1: Sistemas em malha fechada
O modelo nominal G22(s), suposto estritamente próprio sem perda de generalidade,
é definido pela representação estabilizável e detectável:
G22(s) =
A B
C 0
(2.43)
38
com A ∈ Rn×n, B ∈ Rn×m e C ∈ Rp×n. Supõem-se também que o sistema aumentado
P22(s) =
Ap Bp
Cp 0
(2.44)
com spec(Ap) ⊇ spec(A), é estabilizável e detectável e incorpora algumas dinâmicas
fictícias suplementares, notadamente as ponderações frequenciais e/ou os multiplicadores
dinâmicos estáveis Wi(s) e Wo(s). O sistema G22(s) corresponde à dinâmica física6 que se
deseja estabilizar através de um controlador de estrutura estimação-controle. O modelo
P22(s) se distingue de G22(s) pela possível presença de modos não-observáveis por y e/ou
não-controláveis por e. Mesmo que esses dois sistemas sejam equivalentes do ponto de
vista entrada-saída, neste capítulo, as notações (2.43) e (2.44) serão mantidas a fim de
possibilitar a distinção entre as dinâmicas física e aumentada.
Como já foi salientado anteriormente, não existe perda de generalidade ao se supor
que D = 0. De fato, se D 6= 0, é fácil formular um problema equivalente com D = 0
aplicando-se uma transformação linear fracionária apropriada ao controlador. Supondo
que K(s) seja um controlador para G22(s) com D = 0. Para um sistema em que D 6= 0,
o respectivo controlador pode ser reescrito como:
K(s) (I + DK(s))−1 (2.45)
O problema tratado nesta seção pode ser definido como:
Dados os sistemas G22(s) e P22(s) e um controlador estabilizante original (Figure 2.1-(a))
K(s) =
[AK BK
CK DK
](2.46)
onde AK ∈ RnK×nK , nK ≥ n, DK ∈ Rm×p e BK , CK são matrizes reais de dimensõescompatíveis, calcular uma transformação de similaridade, T , tal que o controlador
Ke(s) =
[T−1AKT T−1BK
CKT DK
](2.47)
equivalente a (2.46) do ponto de vista entrada-saída, apresente uma es-trutura explicitamente separada (FIG. 2.1-(b)) e que J11(s) tenha uma es-trutura estimação-controle em relação ao sistema físico G22(s) (FIG. 2.2).
O esquema mostrado na FIG. 2.2 explicita a estrutura procurada. Trata-se, portanto,
de buscar, através do cálculo de T , um controlador equivalente Ke(s) que apresente uma
6Eventualmente incluindo a dinâmica dos atuadores e dos sensores.
39
+ ( ) s P
( ) s G ( ) s W o
( ) s W i w
r e y
z
u
1 y 1 u
-K c
+ C + + B +
A
( ) s Q
_
K f
x ( ) s J ^ x ^
.
FIG. 2.2: Parametrização de Youla e a estrutura estimação-controle
estrutura do tipo LFT inferior em relação ao parâmetro de Youla Q(s) ∈ RH∞:
Ke(s) = J11(s) + J12(s)Q(s) [I + J22(s)Q(s)]−1 J21(s)
onde os coeficientes J11(s), J12(s), J21(s), J22(s) são os elementos da matriz transferência
entre [yT yT1 ]T e [uT uT
1 ]T
J(s) :=
[J11(s) J12(s)
J21(s) J22(s)
]=
A−BKc −KfC Kf B
−Kc 0 Im×m
−C Ip×p 0
(2.48)
Os ganhos Kc, Kf são escolhidos tais que A − BKc e A − KfC sejam estáveis7. A
representação de estado de J11(s) é a mesma vista em (2.40), com x sendo uma estimativa
assintótica dos reais estados x de G22(s). O erro x− x tende a 0 quando o tempo t tende
a infinito, para entradas exógenas nulas (w = 0).
Adotando a notação
Q(s) =
Aq Bq
Cq Dq
(2.49)
para uma realização mínima de Q, onde Aq ∈ Rnq , nq = nK − n e Bq, Cq, Dq têm
7É importante notar que essas restrições de estabilidade são sempre satisfeitas se existe uma tal matriz
T , isto é, se a equivalência entre (2.46) e (2.47) for assegurada.
40
dimensões compatíveis, obtém-se
Ke(s) =
A−BKc −KfC −BDqC BCq Kf + BDq
−BqC Aq Bq
−(Kc + DqC) Cq Dq
(2.50)
Assim, a função de transferência entre r e y pode ser escrita como
Try(s) = (I −G22Ke)−1G22 =
A + BDqC −B(Kc + DqC) BCq B
(Kf + BDq)C A−BKc −KfC −BDqC BCq 0
BqC −BqC Aq 0
C 0 0 0
(2.51)
O vetor de estados de (2.51) contém os estados do sistema (ou estados físicos), os estados
do observador e os estados do parâmetro de Youla: [xT xT xTq ]T .
Por abuso de linguagem, neste trabalho as representações de estado de J(s) e de
Ke(s), expressas por (2.48) e (2.50), são igualmente chamadas de controladores de estru-
tura estimação-controle.
Para facilitar a interpretação da dinâmica de malha fechada, aplica-se à (2.51) a
seguinte transformação de similaridade:
x
xq
x− x
=
I 0 0
0 0 I
−I I 0
x
x
xq
(2.52)
A nova representação de estado em malha fechada (2.53) inclui o erro de estimação do
estado físico, x− x, no vetor de estados:
Try(s) =
A−BKc BCq −B(Kc + DqC) B
0 Aq −BqC 0
0 0 A−KfC −B
C 0 0 0
(2.53)
Nesta representação, o principio da separação aparece claramente. Os autovalores
em malha fechada podem ser separados segundo os n autovalores do controlador, os n
autovalores do estimador e os nq autovalores do parâmetro de Youla: spec(A − BKc),
spec(A−KfC) e spec(Aq), respectivamente.
A função de transferência em malha fechada pode também ser representada a partir
41
do controlador original (2.46) :
Try(s) = (I −G22K)−1G22 =
A + BDKC BCK B
BKC AK 0
C 0 0
=
Acl Bcl
Ccl 0
(2.54)
cujo vetor de estado é[xT (T1x)T (T2xq)
T]T , com
[T1 T2
]= T (2.55)
Com a ajuda dessas definições e notações, pode-se deduzir as condições que garantem
a equivalência entrada-saída entre das representações (2.46) e (2.47).
2.5.2 ESTRUTURA ESTIMAÇÃO-CONTROLE EQUIVALENTE
A partir de (2.47) e (2.50), pode-se deduzir:
AKT − T
[A + BDKC 0
0 Aq
]− T
[BCK
0
]T +
[BKC 0
]= 0 (2.56)
T−1BK =
[Kf + BDq
Bq
](2.57)
CKT =[−(Kc + DqC) Cq
](2.58)
DK = Dq (2.59)
O problema se reduz, então, a resolver em T ∈ RnK×nK a EQ (2.56) e calcular Kc,
Kf , Bq, Cq e Dq usando (2.57), (2.58) e (2.59). Nota-se que Aq é desconhecido em (2.56)
e corresponde a uma variável adicional.
A solução teórica da EQ. (2.56) pode ser simplificada se ela for dissociada em dois
subproblemas. Adotando-se uma partição apropriada de T como em (2.55), obtém-se:
T1(A + BDKC)− AKT1 + T1BCKT1 −BKC = 0 (2.60)
e
(AK − T1BCK)T2 = T2Aq (2.61)
A equação de Riccati generalizada, não simétrica e retangular (2.60) pode ser ainda
reformulada como:
[−T1 I
]H︷ ︸︸ ︷[
A + BDKC BCK
BKC AK
] [I
T1
]= 0 (2.62)
42
Conseqüentemente, a matriz Hamiltoniana H associada a equação de Riccati (2.60)
é a própria matriz da dinâmica do sistema em malha fechada Acl, expressa em (2.54).
A equação de Riccati (2.60) pode ser resolvida em T1 ∈ RnK×n pela técnica clássica de
cálculo de subespaços invariantes que consiste em:
• Calcular um subespaço invariante associado a um conjunto de n autovalores,
spec(Λn), escolhidos entre 2n + nq autovalores de spec(Acl), isto é,[
A + BDKC BCK
BKC AK
] [U1
U2
]=
[U1
U2
]Λn (2.63)
onde U1 ∈ Rn×n e U2 ∈ RnK×n. Tal subespaço é facilmente calculado usando a
decomposição de Schur da matriz Acl.
• Calcular a solução
T1 = U2U−11 (2.64)
cuja existência é garantida quando todos os autovalores de malha fechada forem
distintos.
Usando o resultado anterior e de forma similar, o cálculo da solução T2 ∈ RnK×nq
da equação de Sylvester (2.61) se reduz a encontrar um subespaço invariante associado
com um conjunto de nq autovalores, spec(Aq), escolhidos entre os n + nq autovalores de
spec(AK − T1BCK).
Partindo de (2.47), (2.50) e (2.58), pode-se verificar que
AK − T1BCK = T
[A−KfC 0
−BqC Aq
]T−1 (2.65)
Pode-se, então, estabelecer a seguinte proposição:
Proposição 2.1. Os n autovalores escolhidos para o cálculo da solução T1 da EQ (2.60),
utilizando o método Hamiltoniano, são os n autovalores de realimentação de estados
associados à estrutura estimação-controle equivalente, isto é, spec(A−BKc). Além disso,
a dinâmica de malha fechada restante8 contém a dinâmica do estimador (A − KfC)
aumentada da dinâmica do parâmetro de Youla (Aq).
8Não incluída na seleção da solução de (2.60), isto é, a dinâmica de AK − T1BCK .
43
Esta proposição é uma conseqüência direta de (2.53) e (2.65), consistindo na base da
técnica aqui apresentada. A conclusão nela apresentada permitiu a definição de um algo-
ritmo para a seleção de autovalores necessária para o cálculo de controladores equivalentes
de estrutura estimação-controle (PELLANDA, 2001, 2002b).
Existe uma combinação de soluções segundo a escolha da partição dos autovalores
de malha fechada, primeiramente para o cálculo de T1 e depois para o cálculo de T2. Na
Seção 2.5.3 Algumas considerações adicionais a propósito das possíveis soluções de (2.60)
e (2.61)são apresentadas e discutidas.
Assim, dados um sistema de ordem n e um controlador de ordem nK , pode-se cal-
cular uma transformação linear T−1xK para os estados do controlador tal que eles sejam
separados em duas partes: uma partição correspondente à estimação dos estados físicos
e outra que corresponde aos estados do parâmetro de Youla.
2.5.3 RESTRIÇÕES PARA A SELEÇÃO DE AUTOVALORES
Controlabilidade e observabilidade
Existem várias soluções admissíveis para (2.60) e (2.61). Cada solução corresponde,
respectivamente, a um escolha particular dos n e dos nq autovalores do conjunto de
autovalores de malha fechada. No entanto, algumas combinações de seleção destes não
são admissíveis, pois podem não obedecer as seguintes propriedades:
Proposição 2.2. Considere-se a EQ (2.63). As duas propriedades duais seguintes são
verdadeiras:
(i) Se ∃ λ 6∈ spec(Λn) tal que λ é não-controlável pelo par (A,B), então U1 é singular.
(ii) Se ∃ λ ∈ spec(Λn) tal que λ é não-observável pelo par (A,C), então U2 apresenta
uma deficiência de posto de coluna.
Proposição 2.3. Considere as equações (2.60) e (2.61). Se ∃ λ ∈ spec(Aq) tal que λ é
não-observável por (A,C) ou não-controlável por (A,B), então[
T1 T2
]é singular.
As provas destas proposições são diretas utilizando as relações precedentes e as pro-
priedades de controlabilidade ou de observabilidade.
Pólos complexos conjugados
Outra possível restrição que pode reduzir o número de soluções admissíveis para
(2.60) e (2.61) está relacionado com a não separação dos pares de pólos complexos con-
jugados, ou seja, tais pares devem pertencer à mesma partição de autovalores de malha44
fechada, caso se pretenda obter somente coeficientes reais no cálculo do controlador equi-
valente. Note-se que esta escolha nem sempre é possível. Considerando, por exemplo,
o caso onde os pólos de malha fechada são todos complexos e n é ímpar, então os ga-
nhos terão coeficientes complexos. Para resolver este problema, é imperativo aumen-
tar judiciosamente a dinâmica do modelo do sistema com alguns modos reais (estáveis)
não-controláveis ou não-observáveis. Pode-se mostrar que uma condição necessária para
existência de uma parametrização real é que o número de autovalores de malha fechada
seja maior ou igual à 2(paridade(n)+paridade(nK − n).
Modos duplos
Segundo as Proposições 2.2 e 2.3, a seleção dos autovalores de spec(Acl), para designar
Λn em (2.63) e, em seguida, Aq em (2.61), é baseada numa análise de controlabilidade
e observabilidade modal. Como será discutido mais adiante, esta análise necessita do
cálculo dos autovetores associados aos modos da matriz dinâmica Acl e, por conseguinte,
da forma diagonal Λcl dessa matriz. Outro interesse de se utilizar a forma diagonal é a
diminuição do número de coeficientes de Aq a interpolar.
Entretanto, um problema difícil aparece caso existam modos duplos, ou "quase"
duplos, no espectro de malha fechada. Neste caso, Acl pode ser (quase) não-diagonalizável
e ter uma matriz de autovetores não-inversível ou mal condicionada. Conseqüentemente,
U = [UT1 UT
2 ]T em (2.63) e/ou T2 em (2.61) podem ter uma deficiência de posto de coluna,
dependendo das partições de spec(Acl) escolhidas para compor Λn e Aq.
A estrutura de bloco de Jordan oferece uma base de vetores bem condicionada e
seu cálculo para Acl seria uma solução para este problema. Contudo, tal estrutura é
numericamente difícil de ser determinada. Ao contrário, o cálculo de subespaços in-
variantes associados a formas bloco-diagonais é muito estável numericamente através da
manipulação de decomposições de Schur (GOLUB, 1996). O método consiste em trocar
sistematicamente as posições dos autovetores adjacentes através de rotações de Givens.
Assim, é possível deslocar para as primeiras linhas de Λcl9 cada autovalor isoladamente ou
cada par de modos duplos de uma forma de Schur triangular superior. Este procedimento
permite determinar, um por um, os seus autovetores associados e os elementos não nulos
fora da diagonal principal dos blocos não-diagonalizáveis. Executando-se, então, uma
seqüência de decomposições de Schur, cada uma delas tendo os autovalores ordenados
de uma forma apropriada, pode-se bloco-diagonalizar Acl e obter uma transformação de
9Para maiores detalhes, ver o Algoritmo 7.6.1 de (GOLUB, 1996).
45
similaridade bem condicionada:
AclUsc = UscΛcl (2.66)
com
Λcl = diag(λ1, · · · , λi, · · · , Λr, · · · , Λc, · · · , λnK+n−6) (2.67)
onde λi (i = 1, · · · , nK + n) correspondem aos autovalores distintos e/ou repetidos asso-
ciados a autovetores independentes, e
Λr =
[λr r
0 λ′r
](2.68)
e
Λc =
λc 0 c 0
0 λc 0 c
0 0 λ′c 0
0 0 0 λ′c
(2.69)
são blocos não-diagonalizáveis que correspondem, respectivamente, a dois autovalores
reais repetidos, λr ≈ λ′r, e a dois complexos, λc ≈ λ′c, com seus conjugados, λc ≈ λ′c.
Note-se que é possível ter autovalores repetidos, reais ou complexos, de multiplicidade
maior que 2.
As colunas de Usc em (2.66) que correspondem aos autovalores λi e aos autovalores
repetidos, λr, λc e λc, são seus autovetores associados. As demais colunas, isto é, aquelas
que correspondem aos autovalores repetidos, λ′r, λ′c e λ′c, são chamados de autovetores
generalizados. Um subconjunto qualquer de autovetores forma, sempre, uma base para
um subespaço invariante associado à Acl. Ao contrário, esta propriedade só é conservada
para um conjunto de colunas que contenha autovetores generalizados se os autovetores
respectivos estejam contidos neste conjunto. Ou seja, um conjunto que contenha autove-
tores generalizados associados a λ′r, λ′c e/ou λ′c, deveria conter também os autovetores
respectivos, associados a λr, λc e/ou λc. Entretanto, um bom condicionamento numérico
de T é geralmente obtido se nenhum autovetor generalizado for selecionado para compor
U em (2.63) e se nenhum modo duplo de malha fechada for selecionado para compor
spec(Aq) em (2.61). Então, é conveniente deixar os autovalores relativos aos autovetores
generalizados no conjunto de pólos do observador. Esta análise pode ser generalizada
para autovalores múltiplos.
Uma vez que a maioria das colunas de Usc são autovetores, o uso da forma bloco-
diagonal (2.67) facilita muito a análise de controlabilidade e observabilidade modal. Ao
46
contrário, uma outra tal análise seria muito mais complicada utilizando, por exemplo,
uma forma de Schur clássica do tipo triangular ou quase triangular superior.
Cálculo da forma real
O cálculo de subespaços invariantes exige a manipulação de números complexos
quando o espectro da matrix Hamiltoniana contém modos complexos. Neste caso, a
propagação de erros de arredondamento pode resultar em matrizes de transformação T
que não são rigorosamente reais. Entretanto, o controlador deverá ter exclusivamente
coeficientes reais. Por esta razão, é preferível partir de uma forma real de (2.66). Uma
alternativa à forma de Schur reordenada real é aplicar à (2.66) a seguinte transformação:
AclUscJ = UscJJ−1ΛclJ (2.70)
onde J é uma matriz bloco-diagonal cuja diagonal é composta por 1’s para os modos
reais e por blocos [1/√
2 −j1/√
2
1/√
2 j1/√
2
]
para modos complexos (repetidos ou distintos). Supondo-se que
Λcl = diag
λ1, λ2, λ2,
[λ3 r3
0 λ3
],
λ4 0 c4 0
0 λ4 0 c4
0 0 λ4 0
0 0 0 λ4
, · · ·
(2.71)
onde os autovalores λ1, λ2, λ3 e λ4 são, respectivamente, real distinto, complexo distinto,
real repetido e complexo repetido, obtém-se:
Wsc = UscJ = [u1 uR2 uI
2 u3 u3′ uR
4 uI4 uR
4
′uI
4
′ · · · ] (2.72)
onde o expoente R indica a parte real e o expoente I indica a parte imaginária, e
Λcl = J−1ΛclJ = diag
λ1,
[λR
2 λI2
−λI2 λR
2
],
[λ3 r3
0 λ3
],
λR4 λI
4 cR4 cI
4
−λI4 λR
4 −cI4 cR
4
0 0 λR4 λI
4
0 0 −λI4 λR
4
, · · ·
(2.73)
Assim as matrizes transformadas têm somente elementos reais.
47
Um problema de identificação de modo
Para resolver (2.56), é necessário particionar spec(Acl) em três subconjuntos, a saber:
o conjunto dos modos do controlador (ou de realimentação), o conjunto dos modos do
parâmetro de Youla e o conjunto dos modos do estimador. A etapa seguinte é o cálculo
de dois subespaços invariantes: um subespaço invariante de dimensão n da matriz Acl e
um subespaço invariante de dimensão nq da matriz AK − T1BCK . Estes subespaços são
associados, respectivamente, ao conjunto dos n pólos de realimentação de estados (2.63)
e ao conjunto dos nq pólos do parâmetro de Youla em (2.61).
Se esta abordagem facilita a compreensão conceitual do problema, por outro lado
uma dificuldade numérica pode ocorrer por causa da diferença entre as dimensões de
Acl e de AK . Esta dificuldade aparece notadamente por ocasião da identificação de
certos modos duplos de Acl, associados ou não a autovetores generalizados contidos no
spec(AK − T1BCK). A identificação desses modos demanda um estudo do paralelismo
dos respectivos autovetores. Entretanto, a diferença entre a dimensão dos autovetores em
(2.61) e a dimensão de autovetores em (2.63) torna este estudo complexo e se constitui em
um problema para distinção dos modos duplos. De fato, o cálculo de T em duas etapas,
T1 e T2, é um fator de ambigüidade quando dois modos são próximos um do outro, mesmo
que estes sejam diagonalizáveis. A propagação de erros de arredondamento no cálculo de
T1 pode também comprometer a boa escolha do spec(Aq).
Uma outra maneira de decompor (2.56) em dois problemas dissociados, superando
esta dificuldade, é reescrevê-la como:
[−T I
]
A + BDKC 0 BCK
0 Aq 0
BKC 0 AK
I
T
= 0 (2.74)
Nota-se que um subespaço invariante de dimensão nK é caracterizado por
A + BDKC 0 BCK
0 Aq 0
BKC 0 AK
V11
V12
V2
=
V11
V12
V2
ΛK (2.75)
onde ΛK é uma submatriz de Λcl em (2.71) ou de Λcl em (2.73). A matriz ΛK é associada
à um conjunto de n + nq modos, escolhidos entre aqueles da matriz Hamiltoniana, que
correspondem ao conjunto dos pólos de realimentação de estado e do parâmetro de Youla.
48
O conjunto dos autovalores de spec(ΛK) não deve conter autovalores de Λcl associados a
autovetores generalizados.
A solução de (2.75) implica na resolução de dois problemas:[
A + BDKC BCK
BKC AK
][V11
V2
]=
[V11
V2
]ΛK (2.76)
e
AqV12 = V12ΛK (2.77)
O primeiro problema, (2.76), é similar a (2.63) e consiste em encontrar um subespaço
invariante de Acl, de dimensão nK , associado ao spec(ΛK). O outro problema, (2.77),
consiste em encontrar um subespaço invariante de Aq de dimensão nK . O espaço invari-
ante de Aq de máxima dimensão é, naturalmente, o espaço inteiro que tem dimensão
nq < nK . Respeitando a partição escolhida dos espectros da Hamiltoniana, faz-se então
a escolha de uma base de dimensão nq para formar um tal espaço que é complementado
por vetores nulos a fim de satisfazer a restrição de dimensão imposta por ΛK . Isto é,
V12 = [0 · · · 0 e1 0 · · · 0 e2 0 · · · 0 enq 0 · · · 0] (2.78)
onde a submatriz [e1, · · · , enq ] pode ser escolhida como base canônica de Rnq de forma a
reduzir o número de coeficientes de Q(s) a interpolar. Este arranjo consiste em separar
nq pólos do parâmetro de Youla dentre os nK pólos escolhidos em (2.76). Uma vez que as
posições dos autovalores são completamente definidas por ΛK em (2.76) e (2.77), torna-se
fácil a sua identificação e, conseqüentemente, a alocação das colunas ei de V12 em (2.78).
Finalmente, a mudança de base T é obtida, sem ambigüidade, por
T = V2
[V11
V12
]−1
(2.79)
Considerando a nova formulação (2.76) à (2.79), as condições estabelecidas nas
Proposições 2.2 e 2.3 se tornam:
Proposição 2.4.
(i) Se ∃ λ 6∈ spec(ΛK) tal que λ é não-controlável por (A,B), então
[V11
V12
]é singular.
(ii) Se ∃ λ ∈ spec(ΛK) tal que λ é não-observável por (A,C), então V2 é singular.
49
(iii) Se ∃ λ ∈ spec(Aq) tal que λ é não-observável por (A,C) ou não-controlável por
(A,B), então
[V11
V12
]é singular.
As propriedades enunciadas nestas proposições mostram as restrições que devem ser
satisfeitas para a construção de uma solução inversível T .
Regras Adicionais de Seleção
Testes numéricos mostraram que melhores resultados podem ser obtidos no processo
de interpolação e estimação física, se as seguintes regras suplementares para a seleção da
partição dos autovalores de malha fechada forem respeitadas:
• Alocar em spec(A−BKc) os pólos que são menos controláveis, a fim de reduzir os
ganhos de realimentação de estado.
• Alocar em spec(A−KfC) os pólos mais rápidos e menos observáveis para ter um
estimador de estados eficiente e com ganhos reduzidos.
• Alocar em spec(Aq) pólos rápidos de tal maneira que o parâmetro de Youla se
comporte como uma transmissão direta no controlador.
Dois pontos particulares requerem atenção. O primeiro concerne a maneira de clas-
sificar os modos de acordo com as suas controlabilidade e observabilidade. Um método
útil para medir a controlabilidade de um modo λi em relação a uma entrada j de um
sistema (A,B,C) consiste em calcular
cij =| vT
i bj |‖vi‖ ‖bj‖
onde vi é o autovetor a esquerda associado a λi, bj é a j-ésima coluna de B e vTi bj
é chamado de fator de controlabilidade de λi. Este método foi aplicado em sistemas
físicos, tanto de ordens altas como de baixas (HAMDAN, 1989; MARTINS, 1990). A
controlabilidade global de λi pode ser medida por
ci =‖vT
i B‖‖vi‖
Se A tem todos os autovalores distintos e se supõe que Acl = A em (2.66), então cada
coluna ui de Usc e cada coluna vi de U−Tsc é, respectivamente, um autovetor normalizado
a direita e a esquerda associado ao autovalor λi. Então, ci é definido como:
ci = ‖B(i, :)‖ (2.80)50
onde B = U−1sc B.
Analogamente, têm-se os coeficientes de medida de observabilidade:
oi = ‖C(:, i)‖ (2.81)
onde C = CUsc.
A matriz B (C) é chamada de matriz controlabilidade (observabilidade) modal. Se
a i-ésima linha (coluna) de B (C) é nula, então λi é não-controlável (não-observável).
Examinando-se essas matrizes, pode-se classificar os modos em controláveis e observáveis,
controláveis e não-observáveis, não-controláveis e observáveis, e não-controláveis e não-
observáveis. Assim, a controlabilidade e a observabilidade absolutas dos modos são com-
pletamente definidas por B e C. Essas matrizes permitem classificar os modos segundo
a sua controlabilidade (observabilidade) relativa, uma vez que quanto menor for o coefi-
ciente ci (oi) associado a λi, menos controlável (observável) será o modo.
Supondo que A tenha um modo duplo real (ou complexo) λi = λ′i como em (2.68)
(ou (2.69)), pode-se dizer que:
• c′i (o′i) “grande” ⇐⇒ λi e λ′i são controláveis (observáveis);
• c′i (o′i) “pequeno” ⇐⇒ λ′i é não-controlável (não-observável);
• c′i (o′i) “pequeno” e ci (oi) “grande” =⇒ λi é controlável (observável);
• c′i (o′i) “pequeno” e ci (oi) “pequeno” ⇐⇒ λi é não-controlável (não-observável).
Este resultado é deduzido facilmente pela análise das relações de dependência linear
dos autovetores associados a λi e λ′i. Uma análise similar para os pólos múltiplos seria
uma tarefa mais difícil, porém tais extensões são possíveis. Entretanto, a existência de
pólos múltiplos em malha fechada é um fenômeno mais raro.
O segundo ponto consiste em fazer uma análise de controlabilidade e observabilidade
dos pólos de malha fechada respeitando as restrições impostas sobre o sistema em malha
aberta. Isto é, como definir B e C para calcular ci e oi utilizando (2.80) e (2.81). De fato,
de acordo com a Proposição 2.4, não há restrições impostas para a controlabilidade e a
observabilidade em malha fechada. E ainda, um problema numérico pode ocorrer para
identificar qual autovalor em malha fechada corresponde a um certo pólo não-controlável
ou não-observável, notadamente quando os espectros em malha fechada contêm modos
duplos.
Considerando as realizações em malha aberta e fechada, (2.43) e (2.54), pode-se
estabelecer a seguinte proposição cuja prova é trivial:51
Proposição 2.5. Supondo que CK e BK têm posto completo e considerando que
A = Acl, B =
[B 0
0 InK×nK
]e C =
[C 0
0 InK×nK
]
então, as seguintes propriedade duais são verdadeiras:
(i) λi é não-controlável por (A,B) ⇐⇒ λi é não-controlável por (A, B).
(ii) λi é não-observável por (A,C) ⇐⇒ λi é não-observável por (A, C).
Uma análise da observabilidade e da controlabilidade dos modos de malha fechada
do ponto de vista das entradas e saídas fictícias estabelecidas pelas matrizes B e C é
então justificada neste contexto.
Além dos conceitos de controlabilidade e observabilidade relativas, pode-se utilizar
também uma medida da contribuição de cada modo de malha fechada aos estados do
sistema para escolher aqueles que contribuem mais intensamente no comportamento
dinâmico da planta para compor o espectro do observador. O i-ésimo elemento do autove-
tor a esquerda vj mede o impacto do modo λj sobre a i-ésima variável de estado, enquanto
que o j-ésimo elemento do autovetor a direita ui pondera esta contribuição. O produto
destes dois elementos é adimensional, isto é, independente da escolha de unidades, e mede
a participação líquida. Assim, a matriz (2.82), chamada de matriz de participação, com-
bina os autovetores normalizados a direita e a esquerda e fornece uma maneira eficiente e
prática para medir a associação entre os modos e as variáveis de estado de malha fechada:
P = Usc · ∗U−Tsc = [pij] (2.82)
onde ·∗ designa a multiplicação elemento por elemento e pij é uma medida da participação
relativa do j-ésimo modo na i-ésima variável de estado. O coeficiente pij é chamado de
fator de participação (KUNDUR, 1994; ABED, 2000). A soma dos fatores de participação
que associados a uma variável de estado, ou a um modo qualquer, é igual a 1:∑n+nK
j=1 pij =
1, ∀i; e ∑n+nK
i=1 pij = 1, ∀j. É fácil mostrar que pij é igual à sensibilidade do autovalor
λj ao elemento aii da diagonal da matriz de estados em malha fechada, ∂λj
∂aii.
Se os estados de malha fechada são ordenados como em (2.54), a soma das amplitudes
dos n primeiros elementos da j-ésima coluna de P fornece a participação total de λj nos
n estados de G(s):
pj =n∑
i=1
|pij|, j = 1, · · · , n + nK (2.83)
52
Assim, o conceito de participação modal pode ser utilizado, no lugar dos conceitos de
controlabilidade e observabilidade modais relativas, através de uma escolha apropriada
da participação dos autovalores em malha fechada que respeita a dinâmica natural do
sistema físico.
A referência (PELLANDA, 2001) apresenta um algoritmo sistemático para a escolha
das partições de autovalores de malha fechada, necessária ao cálculo de controladores
equivalentes de estrutura estimação-controle, que leva em consideração as restrições ap-
resentadas nesta seção.
2.5.4 PARÂMETRO DE YOULA DINÂMICO × ESTÁTICO
É importante notar que quando nK é estritamente maior que n, o controlador equi-
valente engloba um parâmetro de Youla dinâmico e um estimador físico. Mas, quando
nK = n, Aq e T2 em (2.61) são vazios e o problema se reduz a solucionar (2.60) em T1 = T
e ao cálculo Kf , Kc e Dq utilizando as relações (2.57)-(2.59), ou seja,
Kf = T−1BK −BDK (2.84)
Kc = −CKT −DKC (2.85)
Q(s) = Dq = DK (2.86)
Em conseqüência, tem-se um parâmetro de Youla estático e uma estrutura de estimação-
controle que pode ser interpretada de duas formas:
• Estimação aumentada: o sistema aumentado P22 (não-mínimo) substitui o sistema
nominal G22; a dinâmica do sistema físico está incluída em P22; e os estados de G22
são estimados se o sentido físico é conservado na dinâmica aumentada.
• Estimação física: P22 substitui G22 com Wi(s) e Wo(s) estáticos (ou nulos) e Ap ∈Rn (ou P22=G22).
Enfim, é igualmente interessante notar, que o parâmetro de Youla dinâmico pode
ser entendido como um parâmetro de Youla estático com uma estimação aumentada.
Com efeito, calculando-se um controlador equivalente compreendendo um parâmetro Q
dinâmico, (Aq, Bq, Cq, Dq), e uma estrutura de estimação-controle, J(s) = f(A, B, C,
Kc, Kf ), pode-se interpretá-lo como um controlador que tem um parâmetro Q estático,
Dq, e uma estrutura de estimação-controle aumentada (ZHOU, 1996)
J(s) = f(A, B, C, Kc, Kf ) (2.87)53
onde A =
[A 0
0 Aq
], B =
[B
0
], C =
[C 0
], Kc =
[Kc −Cq
]e Kf =
[Kf
Bq
].
Em (2.87), J(s) é um controlador equivalente com uma estrutura de estimação-
controle para o sistema fictício aumentado
P (s) =
A 0 B
0 Aq 0
C 0 0
(2.88)
Considerando que um controlador equivalente para P22(s) é também um controlador
equivalente para P (s), que ambos os sistemas incluem as dinâmicas fictícias (Wi(s) e
Wo(s) no primeiro e Aq no segundo) e que geralmente Ap tem nK como dimensão (em
particular em problemas de síntese H∞ e µ), é interessante utilizar G22(s) em lugar de
P22(s) para o cálculo de controladores equivalentes. Obviamente, alguns inconvenientes,
como o cálculo de Aq, Bq e Cq, são evitados quando se utiliza P22(s). Além disso, o
cálculo de um controlador equivalente para P (s) seria impossível utilizando o método
Hamiltoniano10. No entanto, os controladores equivalentes que têm Q(s) dinâmico tam-
bém apresentam algumas vantagens que justificam a utilização do método mais geral
adotado:
• Enquanto que o parâmetro Aq é construído para se ter uma estimação física efi-
ciente, os filtros Wi(s) e Wo(s) são determinados independentemente dessas consi-
derações. Em outros termos, um controlador equivalente para P22(s) consideraria a
identificação da dinâmica fictícia tão importante quanto a identificação da dinâmica
física.
• Enquanto que G22(s) é, em geral, completamente controlável e observável, P22(s)
tem nK − n pólos não-observáveis ou não-controláveis, restringindo as possíveis
escolhas de modos.
• Devido ao modelo do sistema ser utilizado para construir o controlador equiva-
lente, tem-se a necessidade de interpolação da dinâmica fictícia. Pode-se mostra
que a matriz Aq, calculada por este método, tem propriedades que asseguram sua
estabilidade sobre intervalos de interpolação linear.
10Ver as Proposições 2.2, 2.3 e 2.4.
54
2.5.5 CONTINUAÇÃO DOS SUBESPAÇOS INVARIANTES SELECIONADOS
Os métodos clássicos de escalonamento de ganhos decompõem o projeto de um con-
trolador escalonado em vários projetos de controladores lineares (LEITH, 1998a,b). Esses
métodos oferecem estruturas de projeto abertas no sentido em que não existe nenhuma
restrição inerente a uma metodologia particular de síntese de controladores lineares e,
geralmente, o projetista é livre para utilizar os métodos de análise de estabilidade e
desempenho mais convenientes a uma aplicação particular.
Mesmo que a variável de interpolação seja uma função do tempo na implementação
dos controladores interpolados, ela é considerada como um parâmetro estacionário na
etapa de projeto. Considere-se aqui um conjunto de pontos de operação parametrizado
por uma variável de interpolação, θ ∈ R, que evolui num domínio compacto DΘ ⊂ R.Supõe-se que A(θ), B(θ) e C(θ) em (2.43) e Ap(θ), Bp(θ) e Cp(θ) em (2.44) sejam funções
contínuas em DΘ. Admite-se também que K(s, θi) em (2.46) são controladores LTI
estabilizantes projetados para θ = θi, i = 1, 2, · · · , r. O principal objetivo de um processo
de interpolação é definir uma lei de transição contínua e regular entre pontos de operação
adjacentes, θi e θi+1, ∀ i = 1, · · · , r− 1, de forma a preservar o desempenho obtido pelos
controladores LTI nas suas vizinhanças.
As leis de transição sempre introduzem distorções degradando a estabilidade ou o
desempenho nas zonas intermediárias. Entretanto, quando os coeficientes do controlador
evoluem de forma contínua e as amplitudes de variações são pequenas , as distorções
permanecem dentro de limites aceitáveis se o parâmetro do sistema varia lentamente.
Esta seção apresenta um método eficiente, baseado em representações do tipo estimação
e controle, para fazer frente a esse problema.Mais precisamente, o problema estudado nesta seção é:
Dados os sistemas G22(s, θ) (2.43) e P22(s, θ) (2.44), agora considerados dependentesdo parâmetro, e um conjunto inicial de controladores estabilizantes K(s, θi) (2.46),
calcular um conjunto de controladores equivalentes do tipo estimação-controle Ke(s, θi)(2.47),
tal que exista um caminho contínuo ligando a auto-estrutura (eigenstructure) do sistemade malha fechada entre os pontos de operação.
A lei de interpolação que atende os requisitos do problema é obtida por uma técnica
de continuação do tipo Euler-Newton. Esse processo permite o cálculo de um conjunto55
de controladores LTI equivalentes e compatíveis dinamicamente e assegura que existe uma
trajetória contínua que conecte as realizações do tipo estimação-controle.
Considere-se θ o parâmetro normalizado, θ , (θ − θi)/|θi+1 − θi|. Para θ ∈ [θi, θi+1]
tem-se θ ∈ [0, 1].
Em (LUI, 1997), os autores apresentam um método para o cálculo de autovalores de
uma matriz dada H(θ = 1) = H(1) e seus autovetores associados. A partir dos autova-
lores e autovetores conhecidos de uma matriz H(0) real, os autovalores e autovetores da
matriz parametrizada
H(θ) , (1− θ)H(0) + θH(1) (2.89)
são calculados separadamente de θ = 0 até θ = 1 e suas trajetórias são seguidas pelo uso
de técnicas de continuação. Em θ = 1 obtém-se então os dados para H(1). O cálculo
contínuo dos caminhos dos autovalores e autovetores (eigenpath) ao longo da trajetória
de H(θ) permite estabelecer uma correspondência entre as auto-estruturas de H(0) e
H(1), de onde vem o interesse desta técnica neste estudo.
Os autovalores de H(θ) são funções analíticas de θ, salvo para um número finito de
pontos onde alguns autovalores podem apresentar uma singularidade algébrica. Longe
destas singularidades, os autovetores podem ser escolhidos como funções analíticas de θ.
Estas singularidades ocorrem quando autovalores se cruzam ou são submetidos a bifur-
cações no caminho θ ∈ [0, 1]. Tais bifurcações são tipicamente encontradas quando um
modo faz uma transição entre real e complexo ou torna-se duplo, o que introduz dificul-
dades no cálculo, se eles não forem convenientemente manipulados. Algumas técnicas
para resolver esses problemas numéricos, dependendo da natureza da singularidade, são
também apresentados em (LUI, 1997). Mas, o mau condicionamento devido a autovetores
quase colineares aparece de uma forma freqüente durante o processo de continuação. Isto
vem do fato dos autovalores e autovetores serem calculados e seguidos independentemente
e constitui um inconveniente desse método de continuação.
Propõe-se a utilização de métodos similares para a obtenção dos conjuntos correspon-
dentes dos autovalores de matrizes Hamiltonianas adjacentes. Ao invés de seguir cada
modo de maneira independente, a idéia é seguir separadamente cada subespaço invariante
selecionado correspondentes a cada partição escolhida do espectro da Hamiltoniana. Essa
é uma maneira indireta de seguir um conjunto de modos simultaneamente. No contexto
aplicativo de controle deste estudo, esse método é mais confiável de um ponto de vista
computacional uma vez que os problemas de bifurcação e de mau condicionamento de-
vidos a autovetores quase colineares é consideravelmente reduzido. Um único problema
56
aparece quando os cruzamentos envolvem modos pertencentes a partições distintas, o que
é raro na prática.
Considere-se a matriz Hamiltoniana
H ,[
F R
S M
](2.90)
onde F ∈ Rn×n, M ∈ RnK×nK e R, S são reais de dimensões compatíveis. Considere-se
também os conjuntos de equações:[
F R
S M
] [U1
U2
]=
[U1
U2
]Λn (2.91)
[F R
S M
][V1
V2
]=
[V1
V2
]ΛK (2.92)
onde Λn ∈ Rn×n e ΛK ∈ RnK×nK são reais e bloco-diagonais como em (2.73), U1 ∈Rn×n e V2 ∈ RnK×nK são inversíveis, U2 ∈ RnK×n e V1 ∈ Rn×nK . Supõe-se que as
colunas de[
UT1 UT
2
]T
e de[
V T1 V T
2
]T
formam uma base para um subespaço de
H, respectivamente, de dimensão n e nK . Então, T1 = U2U−11 (∈ RnK×n) e T3 = V1V
−12
(∈ Rn×nK ) são, respectivamente, as soluções das equações de Riccati generalizadas, não-
simétricas e retangulares
T1F −MT1 + T1RT1 − S = 0 (2.93)
T3M − FT3 + T3ST3 −R = 0 (2.94)
É importante notar que as colunas de[
I T T1
]T
e de[
T T3 I
]T
geram também
subespaços invariantes de H, de dimensão n e nK , respectivamente. Então, dados H, T1
e T3 pode-se determinar Λn e ΛK . Conseqüentemente, tendo
H(θi) = Acl(θi)
Λn(θi)) = spec(A(θi)−B(θi)Kc(θi))
ΛK(θi)) = spec(A(θi)−B(θi)Kc(θi)) ∪ spec(Aq(θi))
H(θi+1) = Acl(θi+1)
(2.95)
basta executar uma continuação de T1(θi) e de T3(θi), calculadas no ponto de operação θi,
para determinar as dinâmicas correspondentes, Λn(θi+1) e ΛK(θi+1), no ponto de operação
adjacente θi+1.
Admita-se que (2.89) seja a matriz parametrizada associada a dois pontos de operação
adjacentes, H(θ = 0) = H(0) e H(θ = 1) = H(1). A equação de Riccati
F(T, θ) = T (θ)X(θ)− Y (θ)T (θ) + T (θ)Z(θ)T (θ)−W (θ) = 0 (2.96)57
onde θ ∈ [0, 1], corresponde a (2.93) se
T (θ) := T1(θ)
X(θ) := (1− θ)X(0) + θX(1) = F (θ) := A(θ) + B(θ)DK(θ)C(θ)
Y (θ) := (1− θ)Y (0) + θY (1) = M(θ) := AK(θ)
Z(θ) := (1− θ)Z(0) + θZ(1) = R(θ) := B(θ)CK(θ)
W (θ) := (1− θ)W (0) + θW (1) = S(θ) := BK(θ)C(θ)
(2.97)
e a (2.94) se
T (θ) := T3(θ)
X(θ) := (1− θ)X(0) + θX(1) = M(θ) := AK(θ)
Y (θ) := (1− θ)Y (0) + θY (1) = F (θ) := A(θ) + B(θ)DK(θ)C(θ)
Z(θ) := (1− θ)Z(0) + θZ(1) = S(θ) := BK(θ)C(θ)
W (θ) := (1− θ)W (0) + θW (1) = R(θ) := B(θ)CK(θ)
(2.98)
Executando a continuação de T em (2.96), efetua-se a continuação de T1 ou de T3,
de acordo com a escolha (2.97) ou (2.98). Para efetuar a continuação de T no intervalo
[θi, θi+1], é primeiramente necessário subdividi-lo em dois subintervalos da forma
0 = θ0 < θ1 < θ2 < · · · < θL = 1
Utiliza-se, em seguida, um método de continuação de Euler-Newton para calcular a
solução T (θl+1), (l = 0, · · · , L − 1), de F(T, θl+1) = 0, considerando que T (θl) é uma
solução conhecida de F(T, θl) = 0 em (2.96).
Aproximação de Euler
Para obter a solução de Riccati em θl+1, aplica-se o método de Newton na equação
F(T, θl+1) = 0 com a aproximação inicial T (θl+1)(0) = T (θl) + (θl+1 − θl)T (θl), onde
o ponto designa a derivada em relação ao parâmetro. Assim, diferenciando (2.96) em
θ = θl, obtém-se a equação de Sylvester:[T (θl)Z(θl)− Y (θl)
]T (θl) + T (θl)
[X(θl) + Z(θl)T (θl)
]+[
T (θl)X − Y T (θl) + T (θl)ZT (θl)− W]
= 0(2.99)
onde
X = X(1)−X(0)
Y = Y (1)− Y (0)
Z = Z(1)− Z(0)
W = W (1)−W (0)
(2.100)
58
A primeira iteração de Newton T (θl+1)(0) é então obtida resolvendo em T (θl) a
equação de Sylvester (2.99) e, em seguida, calculando o passo de Euler (θl+1 − θl)T (θl).
Algoritmo de Newton
O objetivo é solucionar de uma forma iterativa a equação algébrica
F(T, θl+1) = TX(θl+1)− Y (θl+1)T + TZ(θl+1)T −W (θl+1) = 0 (2.101)
a partir de uma condição inicial T = T (θl+1)(0). A expressão (2.101) pode ser desenvolvida
utilizando o teorema de Taylor. Desprezando os termos de ordem mais elevada, a forma
estendida de (2.101) para a k-ésima interação tem por aproximação de primeira ordem:
F(T (θl+1)
(k) + ∆T, θl+1
)≈
[T (θl+1)
(k)Z(θl+1)− Y (θl+1)]∆T +
∆T[X(θl+1) + Z(θl+1)T (θl+1)
(k)]
+ F(T (θl+1)(k), θl+1) = 0
(2.102)
Se T (θl+1)(k) fosse exata, então F e ∆T seriam nulas. Entretanto, como T (θl+1)
(k)
é somente uma aproximação de T (θl+1), o erro F é finito. Os valores atualizados são
calculados por
T (θl+1)(k+1) = T (θl+1)
(k) + ∆T
onde ∆T é a solução da equação de Sylvester (2.102). O processo é repetido até que o
erro ‖F‖ seja inferior a uma tolerância especificada.
A interação de Newton converge quadraticamente para T (θl+1) sob duas condições:
θl+1−θl suficientemente pequeno; e os caminhos dos autovalores e autovetores pertencente
a diferentes partições sejam suficientemente separados para todos os pontos da trajetória
de continuação.
Se a amostragem é grande pode acontecer que o algoritmo convirja para uma outra
solução. Essa situação é possível mesmo quando os caminhos das auto-estruturas são
aparentemente bem separados. Isso caracterizaria um salto de caminho e comprometeria
a compatibilidade entre as dinâmicas dos controladores LTI vizinhos.
Considere-se finalmente que o número de subintervalos tenha sido escolhido conve-
nientemente, mais que certos modos, correspondentes a partições distintas, se cruzam ou
se aproximam consideravelmente em um ponto dado da trajetória de continuação. En-
tão, o sistema linear implicado na solução da interação de Newton (2.102) pode ser mal
condicionado neste ponto e, conseqüentemente, a convergência do algoritmo de Newton
não é assegurada. A probabilidade que uma tal situação ocorra na prática é pequena.
Entretanto, se ela acontecer, será transparente ao método de continuação no sentido em
59
que o ponto de cruzamento é facilmente percebido. Uma vez que a partição dos modos
relativa à escolha inicial é recuperável em se calculando as matrizes Λn e ΛK no decorrer
de um processo de continuação, esse ponto será tratado da mesma forma que um modo
duplo no ponto de partida.
Algoritmo global
Um procedimento global para o cálculo de um controlador tendo uma estrutura de
estimação-controle no ponto de operação θi+1, a partir de um ponto vizinho θi, é obtido
pelo seguinte algoritmo:
Algoritmo 2.1. Continuação de subespaços invariantes
Etapa 1 : Escolher uma partição do spec(H(θi)), Λn(θi), Aq(θi) e ΛK(θi), e calcular
um controlador equivalente de uma estrutura de estimação-controle (Seções 2.5.2 a
2.5.4) para um ponto de operação θi.
Etapa 2 : Calcular as soluções T1(θi) e T3(θi) para as equações de Riccati (2.93) e
(2.94), de acordo com a partição escolhida de spec(H(θi)).
Etapa 3 : Executar uma continuação de Euler-Newton de T1(θi) e T3(θi) para obter os
correspondentes T1(θi+1) e T3(θi+1).
Etapa 4 : Calcular
Λn(θi+1) =
[I
T1(θi+1)
]H(θi+1)
[I
T1(θi+1)
]−1
(2.103)
e
ΛK(θi+1) =
[T3(θi+1)
I
]H(θi+1)
[T3(θi+1)
I
]−1
(2.104)
e bloco-diagonalizá-los como Λcl em (2.73).
Etapa 5 : Separar Aq(θi+1) de ΛK(θi+1) comparando os autovetores up e vj associados
a Λn(θi+1) e ΛK(θi+1), respectivamente. Isto é, calcular
cos(θpj) =| uT
p vj |‖up‖ ‖vj‖
para todo p, j (p = 1, 2, · · · , n e j = 1, 2, · · · , nK) para separar os vj que não
tenham um up paralelo correspondente.60
Etapa 6 : Utilizando Aq(θi+1) e ΛK(θi+1), calcular o controlador equivalente no ponto
de operação θi+1, ou seja, calcular T (2.79) em seguida calcular Kc, Kf , Bq, Cq e
Dq utilizando (2.57), (2.58) e (2.59).
Nota-se que aqui T1(θi+1) corresponde exatamente à primeira partição de T (θi+1)
em (2.55), enquanto que a segunda partição T2(θi+1) é determinada por Aq(θi+1). Por
outro lado, como a solução de uma equação de Riccati é independente da ordem dos
autovetores e autovalores, os ganhos Kc(θi+1) e Kf (θi+1) são independentes do arranjo
de Λn(θi+1) e ΛK(θi+1). Uma possível troca de ordem dos autovalores durante o processo
de diagonalização (Etapa 4) afetaria somente as posições das colunas de T2(θi+1). No
entanto, a ordem correta é facilmente recuperada analisando-se d proximidade entre os
autovalores de Aq(θi) e de Aq(θi+1).
Nota-se também que se nK = n, então T3(θ) = T−11 (θ). Em conseqüência, necessita-
se somente executar um prolongamento contínuo de T1(θ). Além disso, quando se tratar
de um conjunto de mais de dois controladores LTI, uma vez que a Etapa 1 é executada
no início do processo (por exemplo, i = 1), somente as Etapas 2 e 3 são necessárias para
determinar toda a familia de transformações lineares de estados (i = 2, · · · , r). Esse
procedimento permite calcular todo um conjunto de controladores equivalentes a partir
de uma escolha única da partição de autovalores de malha fechada e assegura que existe
uma trajetória contínua de transformações que conecta as suas realizações de estado do
tipo estimação-controle.
61
3 CONTRIBUIÇÕES PARA A IDENTIFICAÇÃO DE SISTEMAS
3.1 INTRODUÇÃO
Como foi tratado no Capítulo 1, a modelagem é classificada como: caixa-preta, caixa-
branca e caixa-cinza. Aborda-se neste capítulo a identificação dos três tipos de modela-
gens. Este enfoque visa utilizar as vantagens de combinar as modelagens de caixas-pretas
e caixas-brancas no que se refere ao uso de informações a priori. Segundo (CORRÊA,
2001), por ser um assunto relativamente novo, muitos dos problemas desta modelagem
estão em aberto. De acordo com (SOHLBERG, 1998) e (KARPLUS, 1976), classifica-se
a “tonalidade de cinza” em três níveis de intensidade, conforme o uso de informações
a priori no processo de modelagem. Esses tons se aplicam em sistemas econômicos,
hidráulicos e de poluição do ar, do cinza mais escuro para o mais claro, respectivamente.
Isso reflete a abrangência desta modelagem.
Nas modelagens onde as grandezas físicas/químicas que caracterizam a variação
paramétrica da dinâmica da planta podem ser medidas, como as saídas naturais ou a
partir de um sensoriamento aumentado do sistema, uma técnica é proposta neste capí-
tulo. Ela é paramétrica e permite a identificação de modelos quasi -LPV e uma classe de
modelos não-lineares. Um recurso desta técnica é a utilização de um parâmetro moni-
torado como argumento de uma função. Essa metodologia foi desenvolvida com o auxílio
da série de Taylor, determinando os coeficientes levantados mediante a aplicação da téc-
nica dos mínimos quadrados.
Vários aspectos importantes são discutidos no decorrer da descrição destas técnicas
de identificação. Porém, as linhas gerais citadas em (AGUIRRE, 2004) e (LJUNG, 1987)
são adotadas como as principais etapas a serem seguidas na identificação. São elas:
• Testes dinâmicos e coleta de dados: essa é a fase que envolve a geração de dados. Os
principais fatores estão ligados à escolha do sinal de excitação (quando possível), à
execução do teste e à escolha do período de amostragem. Em determinados casos,
a planta deve ser identificada quando em operação. Neste caso, se a entrada não
for capaz de excitar todos os modos, a identificação será restrita a uma faixa de
operação.
• Escolha da representação matemática: esta etapa está focada na representação da62
dinâmica, seja em funções de transferência ou em matrizes de espaço de estados.
Um abuso de linguagem é empregado na afirmativa anterior, pois uma das técnicas
de identificação deste estudo são aplicáveis a sistemas não-lineares, em geral.
• Determinação da estrutura do modelo: no caso de identificações lineares, é feita a
escolha do número de pólos e zeros. Como a técnica de modelos não-lineares ou
quasi -LPV dependem da identificação linear em pontos de operação da planta, a
estrutura dos modelos variantes no tempo fica condicionada, basicamente, à deter-
minação da estrutura dos modelos lineares.
• Estimação de parâmetros: trata-se da seleção do algoritmo e/ou dos métodos
numéricos para o cálculo dos parâmetros a serem utilizados.
• Validação do modelo: uma vez estabelecido um modelo ou uma família de modelos,
testes são usados para verificar se esses reproduzem as características da planta.
Além disso, no caso de uma família de modelos, uma comparação entre estes in-
dicará o que apresenta menor erro. Tal fato pode depender significativamente do
tipo de modelo empregado.
3.2 REPRESENTAÇÕES
Um sistema linear pode ser considerado variante no tempo ou não-estacionário - Li-
near a Parâmetros Variáveis (LPV) ou Linear Variante no Tempo (ou LTV, do inglês
Linear Time Variant) - quando um ou mais parâmetros do seu modelo variam signi-
ficativamente com o tempo e não podem ser considerados como parâmetros incertos. A
evolução destes sistemas pode ser matematicamente representada pela solução ϕ(t, u, x0)
de um sistema de equações diferenciais (COCA, 1996b):
x = f(x, u) (3.1)
onde x(t) é um vetor de dimensão n (n ∈ N) das variáveis de estado, u(t) é o vetor
de dimensão q (q ∈ N) das variáveis de entrada, x0 representa as condições iniciais e
f : Rn+q → Rn é um mapeamento não-linear e contínuo que representa a dinâmica do
sistema.
A saída do sistema é descrita por:
y = h(x, u) (3.2)
63
onde y é um vetor de dimensão p ∈ N. Quando f(.) e h(.) são funções lineares, (3.1) e
(3.2) formam um sistema de equações denominadas de espaço de estado. Se essas funções
são não-lineares, pode-se adotar a representação quasi -LPV, em um grande número de
casos.
3.3 IDENTIFICAÇÃO DE SISTEMAS LINEARES
As muitas técnicas de controle linear desenvolvidas, aliadas à simplicidade em se
trabalhar com modelos lineares, incentivam os métodos de identificação apresentados
nesta seção. São mostradas duas possibilidades de identificação de sistemas lineares
estáveis no presente trabalho. Em ambos os casos, não foi considerada a presença de
ruído nas medições amostradas.
A Seção 3.3.1 se baseia na transformada z e na função de transferência discreta.
Utilizam-se amostras discretas dos dados temporais de entrada e de saída da planta. O
resultado é a identificação dos coeficientes da função de transferência discreta e, conse-
qüentemente, os seus pólos e zeros.
A Seção 3.3.2 emprega a resposta em freqüência da planta para a determinação dos
coeficientes da função de transferência em regime contínuo. É possível se obter uma
seqüência discreta com valores em freqüência X(jw), a partir de valores amostrados de
um sinal vetorial x(t). Para melhor compreensão da obtenção da resposta em freqüência,
ver (VALLE, 2005) e (ADES, 2005).
O problema consiste em obter a função de transferência, discreta ou contínua, de
um sistema do tipo caixa-preta para entradas e saídas contínuas no tempo. As técnicas
empregadas foram desenvolvidas para sistemas SISO (uma entrada e uma saída). A
extensão à um sistema MIMO (múltiplas entradas e múltiplas saídas) pode ser feita
aplicando-se a metodologia a cada um dos elementos da matriz de transferência.
3.3.1 IDENTIFICAÇÃO USANDO AMOSTRAS DISCRETAS
Para um período de amostragem τ suficientemente pequeno, pode-se obter a
amostragem dos sinais de entrada e de saída da planta. Admite-se que o sinal apli-
cado à entrada do sistema excite todos os modos desta. Com este conjunto de dados é
possível equacionar sua dinâmica na forma de uma FT discreta:
G(z) =b0z
n + b1zn−1 + ... + bn−1z + bn
zn + a1zn−1 + ... + an−1z + an
=Y (z)
U(z)(3.3)
64
Sendo y(k) , y(kτ), manipulando (3.3) e aplicando a inversa da transformada z:
b0znU(z) + b1z
n−1U(z) + ... + bn−1zU(z) + bnU(z) =
= znY (z) + a1zn−1Y (z) + ... + an−1zY (z) + anY (z) ⇐⇒
b0u(k + n) + b1u(k + n− 1) + ... + bn−1u(k − 1) + bnu(k) =
= y(k + n) + a1y(k + n− 1) + ... + an−1y(k − 1) + any(k)
(3.4)
Fazendo k = k − n e isolando y(k) em (3.4), tem-se a seguinte equação diferença:
y(k) = −a1y(k − 1)− a2y(k − 2)− ...− any(k − n) + b0u(k) + ... + bnu(k − n) (3.5)
Escolhe-se então uma função custo, a ser minimizada, para uma identificação efi-
ciente. Pretende-se que as amostras discretas fornecidas pela saída da planta sejam
reproduzidas pelo modelo identificado. A função custo é calculada através da norma
dois da diferença entre os vetores de saída medido e aquele fornecido pelo modelo para a
mesma entrada aplicada na planta. Tem-se, então:
J = ‖yp(k)− y(k)‖2 (3.6)
onde y(k) é a saída do modelo linear estimada e yp(k) é a saída medida da planta.
Arbitrando-se "n" como a ordem da função de transferência discreta a ser identi-
ficada, sabendo que m À n e conhecendo-se as m amostras da entrada e da saída, é
possível equacionar a dinâmica desse sistema a partir de (3.5) como:
y(m) = −a1y(m− 1)− ... + b0u(m) + ... + bnu(m− n)
y(m− 1) = −a1y(m− 2)− ... + b0u(m− 1) + ... + bnu(m− n− 1)...
y(n + 1) = −a1y(n)− ...− any(1) + b0u(n + 1) + ... + bnu(1)
(3.7)
Definindo:
• Y =
y(m)
y(m− 1)...
y(n + 1)
• M =
y(m− 1) . . . y(m− n) u(m) . . . u(m− n)
y(m− 2) . . . y(m− n− 1) u(m− 1) . . . u(m− n− 1)... . . . ...
... . . . ...
y(n) . . . y(1) u(n + 1) . . . u(1)
65
• Θ =[−a1 . . . −an b0 . . . bn
]T
∈ R2n+1
O sistema (3.7) pode ser escrito de forma matricial como:
Y = MΘ (3.8)
O problema de identificação consiste em estimar os coeficientes Θ na equação ma-
tricial (3.8). Com m suficientemente grande, a técnica dos mínimos quadrados (pseudo-
inversa) pode ser aplicada. Assim, os coeficientes Θ podem ser determinados por:
MT MΘ = MT Y
Θ = (MT M)−1MT Y
Θ = M †Y
(3.9)
A matriz M e o vetor Y devem ser definidos com os valores medidos da planta (u(k)
e yp(k)). Com a obtenção de Θ, tem-se o modelo discreto linear identificado G(z).
Período de Amostragem e Sinal de Entrada
A escolha do período de amostragem τ de uma planta contínua no tempo requer uma
análise detalhada para uma representação fidedigna da dinâmica em tempo discreto, na
técnica de identificação proposta.
Quando a planta for contínua no tempo e o período de amostragem puder ser esco-
lhido, deve-se considerar o problema de condicionamento numérico da solução proposta.
Este fator está diretamente relacionado com o mapeamento entre a variável z, da FT
discreta, e a variável s, da FT contínua, dado por:
z = esτ (3.10)
De (3.10) é fácil ver que períodos de amostragem muito pequenos alocarão os pólos e
zeros em regiões muito próximas no plano z e com módulos tendendo a 1. Isso resulta em
alocação de pólos e zeros com valores decimais de significância relativa bastante reduzida,
ou seja, muito próximos. Como os coeficientes da FT são oriundos de operações aditivas
e multiplicativas com os pólos e zeros, para dois períodos diferentes, ambos próximos de
zero, os coeficientes das respectivas FT só se diferenciarão após várias em casas decimais.
Sabendo-se que os computadores tem um espaço limitado para guardar os dígitos de uma
variável, no caso os coeficientes da FT discreta, para τ → 0 pode haver arredondamentos
que comprometam uma identificação precisa. Para contornar este problema, períodos de
amostragem maiores devem ser privilegiados quando possíveis de serem escolhidos.66
Para um τ convenientemente escolhido, o sinal de entrada para identificação pode
ser um degrau unitário, sinal este que apresenta uma riqueza espectral suficiente para
excitar todos os modos da planta.
Quando τ for imposto e isso estiver causando os problemas aludidos, é possível tentar
minimizar estes efeitos com, pelo menos, dois recursos:
• multiplicar a entrada por uma constante;
• aplicar um somatório de pulsos retangulares finitos e defasados (trem de pulsos);
Uma combinação de ambas as técnicas pode ser implementa.
Aplicando a propriedade da homogeneidade de funções lineares, multiplicar o sinal
de entrada por uma constante tem igual efeito na resposta sem que haja alteração no
modelo. Entretanto, em termos operacionais, este recurso se mostra eficiente em alguns
casos. Quando o cálculo dos coeficientes for feito usando a forma recursiva da pseudo-
inversa, devido ao elevado número de amostras, o ajuste dos coeficientes dar-se-á por
uma diferença entre o valor medido e a estimação. Se a entrada for multiplicada por
uma constante que seja suficiente para ampliar esta diferença, o ajuste do coeficiente
convergirá com um menor número de amostras.
Um trem de pulsos irá gerar maior riqueza de amostras que permitirá um ajuste mais
preciso nas casas decimais de menor significância. Este ajuste fino apresenta melhores
resultados quando se usa a forma recursiva da pseudo-inversa.
Cabem duas últimas análises em relação ao modelo construído. A primeira é em re-
lação a estabilidade do modelo obtido. Por premissa a planta é estável e em conseqüência
não se admite modelos instáveis. Quando a identificação apresentar pólos que tornem o
modelo instável, deverá ser feita a projeção dentro do círculo unitário do plano z e de-
sprezar a parte instável, uma vez que uma planta estável não tem qualquer projeção fora
do círculo unitário do plano z (ver Seção 2.2). A segunda análise é relativa à alocação
dos pólos e zeros, pois ao se arbitrar a ordem da identificação, a função de transferência
discreta G(z) poderá ser não-mínima (ver Seção 2.3).
3.3.2 IDENTIFICAÇÃO USANDO A RESPOSTA EM FREQÜÊNCIA
Como abordado em (VALLE, 2005) e (ADES, 2005), uma vez obtida a resposta em
freqüência G(jω) de uma planta linear ou de um modelo linearizado em um ponto de
operação numa planta não-linear, é possível determinar a função de transferência (FT)
67
referente à dinâmica do sistema estudado. A abordagem da técnica aqui apresentada é
similar ao caso discreto anterior.
Considere a dinâmica de uma planta representada pela FT, G(s), na forma:
G(s) =b0s
n + b1sn−1 + ... + bn−1s + bn
sn + a1sn−1 + ... + an−1s + an
(3.11)
Deseja-se identificar um modelo que minimize a função custo:
J(Θ) = ‖G(jω)− G(Θ, jω)‖2 (3.12)
onde ω =[
ω1 ω2 . . . ωh
]T
∈ Rh.
O vetor ω contém as freqüências amostradas e selecionadas para a identificação do
sistema. Sugere-se que seja adotado um espaçamento logarítmico entre as freqüências
em ω. O vetor Θ representa os coeficientes a serem ajustados da função de transferência
identificada e G(s) o modelo identificado.
Dado que a resposta em freqüência G(jω) é conhecida, (3.11) pode ser particularizada
como:
G(jωi) = G(s)|s=jωi=
b0(jωi)n + b1(jωi)
n−1 + ... + bn−1(jωi) + bn
(jωi)n + a1(jωi)n−1 + ... + an−1(jωi) + an
(3.13)
Para evitar erros de interpretação nas equações desenvolvidas subseqüentes, usar-se-á
a notação Gjω , G(jω). De 3.13:
Gjωi
{(jωi)
n +∑n−1
p=0 [an−p(jωi)p]
}=
∑nq=0[bn−q(jωi)
q]
Gjωi(jωi)
n = −Gjωi
∑n−1p=0 [an−p(jωi)
p] +∑n
q=0[bn−q(jωi)q]
(3.14)
A EQ (3.14) pode ser escrita de forma matricial, para i ∈ {1, 2, ..., h}, como:
H ,
Gjω1 .(jω1)n
Gjω2 .(jω2)n
...
Gjωh.(jωh)
n
=
=
Gjωi(jω1)
n−1 . . . Gjωi(jω1)
0 (jω1)n . . . (jω1)
0
Gjωi(jω2)
n−1 . . . Gjωi(jω2)
0 (jω2)n . . . (jω2)
0
... . . . ...... . . . ...
Gjωi(jωh)
n−1 . . . Gjωi(jωh)
0 (jωh)n . . . (jωh)
0
−a1
...
−an
b0
...
bn
(3.15)
68
Reescrevendo (3.15) de forma a evitar coeficientes complexos na função de transfe-
rência:
H =
Re{Gjω1 .(jω1)n}
Im{Gjω1 .(jω1)n}
...
Re{Gjωh.(jωh)
n}Im{Gjωh
.(jωh)n}
= NΘ (3.16)
onde H ∈ R2h, Θ é o vetor dos coeficientes a ser encontrado e
N =
Re{(jω1)n−1} . . . Re{(jω1)
0} Re{(jω1)n} . . . Re{(jω1)
0}Im{(jω1)
n−1} . . . Im{(jω1)0} Im{(jω1)
n} . . . Im{(jω1)0}
Re{(jω2)n−1} . . . Re{(jω2)
0} Re{(jω2)n} . . . Re{(jω2)
0}Im{(jω2)
n−1} . . . Im{(jω2)0} Im{(jω2)
n} . . . Im{(jω2)0}
... . . . ...... . . . ...
Re{(jωh)n−1} . . . Re{(jωh)
0} Re{(jωh)n} . . . Re{(jωh)
0}Im{(jωh)
n−1} . . . Im{(jωh)0} Im{(jωh)
n} . . . Im{(jωh)0}
Novamente, a pseudo-inversa pode ser empregada em (3.16) para determinar o vetor
Θ:NT NΘ = NT H
Θ = (NT N)−1NT H
Θ = N †H
(3.17)
Desta forma, tem-se Θ =[−a1 . . . −an b0 . . . bn
]T
que corresponde à
seguinte FT
G(s) =b0s
n + b1sn−1 + ... + bn−1s + bn
sn + a1sn−1 + ... + an−1s + an
(3.18)
As observações feitas para o caso discreto, relativas a estabilidade e alocação de pólos
e zeros, são também válidas para este caso.
3.4 IDENTIFICAÇÃO DE SISTEMAS QUASI -LPV
Conforme explicado na Seção 2.4, um sistema linear pode ser considerado variante
no tempo ou não-estacionário − quasi -LPV ou Lineares Variantes no Tempo (LTV, do
inglês Linear Time Variant) − quando um ou mais parâmetros do seu modelo variam am-
plamente com o tempo e não podem ser considerados como parâmetros incertos. Muitas
69
vezes, os sistemas não-lineares, comumente encontrados na prática, podem ser tratados
como sistemas lineares não-estacionários.
Um sistema quasi -LPV é regido por equações diferenciais da forma:
Modelo M
{x(t) = A(θ(t))x(t) + B(θ(t))u(t)
y(t) = C(θ(t))x(t) + D(θ(t))u(t)
O parâmetro θ(t) varia continuamente no tempo dentro de um domínio paramétrico
PΘ. No caso escalar o parâmetro estará limitado por um valor máximo e mínimo, isto é,
θ ∈ [θm, θM ].
A técnica de identificação proposta pode ser estendida a todas as classe de plantas
cujas dinâmicas não-lineares possam ser reformuladas como sistemas quasi -LPV.
Para uma planta com as características descritas nos parágrafos anteriores e sendo
estável, pode-se construir um modelo quasi -LPV desse sistema. Para isso é necessário
identificar preliminarmente um conjunto de modelos lineares em tantos pontos de ope-
ração quanto necessários. Estes pontos corresponderão a valores pré-determinados do
parâmetro endógeno e/ou exógeno θ que determina aqueles pontos.
A técnica aqui descrita objetiva identificar um modelo quasi -LPV, que recebendo a
medição em tempo real do parâmetro θ e a entrada da planta, faz com que a saída do
modelo identificado convirja para a resposta da planta física.
A metodologia consiste em reescrever os diversos coeficientes do conjunto de mode-
los, obtidos na etapa de identificação dos pontos de operação da planta, como funções
polinomiais do parâmetro θ, considerado disponível para medição.
A identificação quasi -LPV tratada aqui pode ser classificada como caixa-cinza, de-
vido à necessidade do conhecimento prévio do parâmetro que torna a dinâmica do sis-
tema variante no tempo. Esta é a única informação prévia necessária para se levantar a
dinâmica física da planta.
3.4.1 METODOLOGIA
Considera-se que a dinâmica a pequenas pertubações de uma planta não-estacionária
possa ser representada por um conjunto de modelos lineares. Esses modelos podem
ser parametrizados por um sinal contínuo θ a ser medido. O problema consiste, numa
primeira etapa, em identificar os modelos lineares Gi(s), i ∈ {0, 1, ..., v}, v ∈ N, dos
pontos de operação de interesse da planta.
Numa segunda etapa, procura-se determinar, a partir dos modelos Gi(s), um novo
modelo quasi -LPV M, cujos coeficientes são funções polinomiais do parâmetro (sinal θ).70
Esse modelo M se torna linear com a determinação do valor de θ.
Adota-se para cada modelo Gi(s), uma FT de forma:
Gi(s) =bi0s
n + bi1sn−1 + ... + bi(n−1)s + bin
sn + ai1sn−1 + ... + ai(n−1)s + ain
(3.19)
para aij, bij ∈ R, j ∈ {0, 1, ..., n}.
Como o conjunto de modelos LIT Gi(s) foi obtido pela identificação do sistema, para
se comparar as entradas e as saídas do modelo com as da planta no ponto de operação i,
deve-se adicionar valores constantes a essas variáveis.
yi(t) = yi(t)− yNi
ui(t) = ui(t) + uNi
(3.20)
Os valores yNi e uNi ∈ R são as constantes que se referem ao ponto de operação i
da planta em relação ao modelo Gi(s). Os sinais y(t) e u(t) são a entrada e a saída da
planta, respectivamente. Os sinais y(t) e u(t) são a entrada e a saída, respectivamente,
correspondentes a aplicação no modelo linear.
No ponto de operação em regime estacionário, os estados do modelo são nulos. Isso
é incongruente com a pretensão da metodologia aqui aplicada, uma vez que os estados
não deverão ter nova referência (valores iniciais) a cada novo ponto de operação, mas
deverão ser contínuos, como são numa modelagem da planta, e guardar uma relação com
o modelo, ou seja, representar grandezas fixas definidas com o mesmo critério para todo o
conjunto de modelos lineares Gi(s), mesmo que estas não tenham significado físico. Isso
não ocorre com o conjunto Gi(s).
Assim sendo, deve-se compatibilizar os modelos lineares disponíveis, Gi(s), com as
necessidades de equacionamento proposto, Gi(s). Em outras palavras, deseja-se um con-
junto de modelos lineares cujos estados representem as mesmas grandezas em toda a faixa
de operação do parâmetro, sem deslocamento da entrada e da saída.
Para uma representação compatível com a metodologia proposta, uma manipulação
se faz necessária no conjunto de modelos LIT Gi(s). De (3.20):
Yi(s) = L[yi(t)]
Yi(s) = L[yi(t) + yNi]
Yi(s) = Y i(s) + yNi
s.
(3.21)
Analogamente, U(s) = U i(s) + uNi
s. Portanto:
G(s)i =Yi(s)
Ui(s)=
Y i(s) + yNi
s
U i(s)− uNi
s
=s.Y i(s) + yNi
s.U i(s)− uNi
(3.22)
71
obtém-se as FT relacionadas com a planta.
Para manter a coerência de notação, faz-se yNi = bi(n+1) e uNi = −ai(n+1). A nova
FT é dada por:
Gi(s) =bi0s
n+1 + bi1sn + ... + bins + bi(n+1)
sn+1 + ai1sn + ... + ains + ai(n+1)
(3.23)
Com j ∈ {0, 1, . . . , (n + 1)}, define-se:
aj =[
a1j a2j ... avj
]T
bj =[
b0j b1j ... bvj
]T
θj =[
θ1j θ2j ... θvj
]T
(3.24)
para θi ∈ R e i ∈ {1, 2, ..., v}.Com as definições acima, é possível estabelecer funções bases para se determinar a
dinâmica dependente do parâmetro para ai(θ) e bi(θ) que definem o modelo quasi -LPV:
G(s, θ) =b0(θ)s
n+1 + b1(θ)sn + ... + bn(θ)s + bn+1(θ)
sn+1 + a1(θ)sn + ... + an(θ)s + an+1(θ)(3.25)
O conceito de função de transferência dependente no parâmetro não tem sentido
quando este varia no tempo. A notação G(s, θ), onde o parâmetro θ pode assumir vários
valores em um dado domínio, é utilizada para indicar que um sistema G do tipo quasi -
LPV tendo como modelo freqüencial
G(s, θ) = C(θ)[sI − A(θ)]−1B(θ) + D(θ)
para valores estacionários de θ, ou seja, aqueles para os quais dθ/dt = 0, e possui como
modelo de estados
G(s, θ) =
A(θ) B(θ)
C(θ) D(θ)
e
[x(t)
y(t)
]=
[A(θ) B(θ)
C(θ) D(θ)
][x(t)
u(t)
]
para as trajetórias não-estacionárias de θ(t).
Particularmente, para definir as funções ai(θ) e bi(θ) em (3.25), será utilizado uma
base polinomial em θ do tipo f(θ) = α0.θp + α1.θ
p−1 + ... + αp, p ∈ N. Para chegar a
estas funções, é necessário arbitrar uma ordem de grandeza p para as funções a serem
aproximadas na variação dos coeficientes. Assim, pode-se equacionar:
aj =
θp1 θp−1
1 ... 1
θp2 θp−1
2 ... 1...
... . . . ...
θpk θp−1
k ... 1
α0
α1
...
αp
= Xαj (3.26)
72
Analogamente para bi. Mais uma vez, os coeficientes podem ser obtidos pela pseudo-
inversa.
αj = X†aj (3.27)
73
4 ESCALONAMENTO DE GANHOS EM ESTRUTURAS DO TIPO
ESTIMAÇÃO-CONTROLE
4.1 INTRODUÇÃO
Neste capítulo são apresentadas três metodologias para o controle de sistemas não-
estacionários que evoluem em grandes faixas de operação, utilizando interpolação clássica
de controladores dinâmicos predeterminados.
Na Seção 2.5 é apresentada uma metodologia (PELLANDA, 2000, 2001, 2002b) para
o cálculo de uma seqüência de transformações de similaridade a ser aplicada à família
original de controladores robustos determinados a priori, de forma que o conjunto de
controladores obtidos, equivalentes aos primeiros do ponto de vista entrada-saída, tenha
uma estrutura coerente, fisicamente interpretável e de fácil interpolação. Depois da trans-
formação, os controladores apresentam uma estrutura do tipo estimação-controle, o que
facilita a interpolação e a implementação em tempo variante. Contudo, a questão da
escolha dos pontos de projeto não é abordada nesta técnica, sendo assumido que eles são
convenientemente escolhidos. Uma abordagem para uma escolha sistemática do conjunto
de pontos de síntese, no contexto dessa técnica, é proposta na Seção 4.2.
A segunda técnica apresentada neste capítulo (Seção 4.3) constitui uma contribuição
deste trabalho e trata o problema do escalonamento de ganhos clássico de uma forma
distinta da primeira. O que se propõe é utilizar uma abordagem quasi -LPV para a
modelagem do sistema. A idéia principal é estabelecer uma realimentação de estados
estimados onde os ganhos são dependentes do parâmetro definido no modelo quasi -LPV.
Além disso, os ganhos são calculados como funções do parâmetro, de forma que o espectro
de autovalores da matriz dinâmica de malha fechada seja alocado em uma posição fixa
predefinida para todo o conjunto possível de trajetórias contínuas do parâmetro.
Por fim, é proposto um método que objetiva o cancelamento das não-linearidades
através da realimentação de estados (Seção 4.4). Assim como no caso anterior, é utilizada
uma técnica de alocação de autovalores em malha fechada. A técnica se aplica a uma
classe restrita de sistemas não-estacionários. Basicamente, a lei de controle proposta
busca anular os efeitos das dinâmicas indesejadas e impõe o controle sobre a parte linear
do modelo. Esse método também objetiva impedir que variações bruscas sejam inseridas
artificialmente na dinâmica de malha fechada.
74
4.2 ESCOLHA SISTEMÁTICA DOS MODELOS DE SÍNTESE
Na linearização clássica de um sistema não-linear x = f(x), onde f(x) : Rn → Rn é
um mapeamento não-linear e contínuo, é comumente utilizado o método de Taylor para
obtenção de um modelo linearizado em torno de um ponto de funcionamento escolhido
x = x. A resposta desse modelo linear é muito próxima da resposta do modelo não-
linear para pequenas perturbações em torno do ponto de operação considerado. Em
outras palavras, o método clássico de linearização atribui um coeficiente linear baseado
na derivada de primeira ordem e no valor da variável no ponto de operação previamente
escolhido (x = x): xδ = Axδ, onde xδ = x− x e
A =
∂f1(x,u)∂x1
∂f1(x,u)∂x2
... ∂f1(x,u)∂xn
∂f2(x,u)∂x1
∂f2(x,u)∂x2
... ∂f2(x,u)∂xn
...... . . . ...
∂fn(x,u)∂x1
∂fn(x,u)∂x2
... ∂fn(x,u)∂xn
x=x
(4.1)
Porém, no contexto clássico do escalonamento de controladores, utiliza-se um con-
junto de pontos de operação, uma vez que a técnica é aplicada em sistemas não-lineares
com ampla faixa de operação para os quais geralmente não existe controlador estabilizante
LTI. Ocorre que o comportamento dinâmico do sistema em malha fechada é fortemente
influenciado pela escolha desse conjunto.
Para exemplificar, supõe-se que o modelo tenha sido linearizado em um ponto onde
sua faixa de validade seja muito estreita para um dado critério de erro. Neste caso,
para grandes perturbações, onde o sistema evolua desse ponto de operação para um
ponto vizinho do conjunto escolhido, torna-se muito difícil manter a estabilidade ou
o desempenho do sistema em malha fechada através da interpolação dos respectivos
controladores LTI, a menos que o ponto vizinho esteja muito próximo. Ou seja, a não-
estacionariedade do sistema pode não ser devidamente compensada pela variação do
controlador.
Uma situação ainda mais desfavorável pode ocorrer se o sistema for conduzido, por
uma entrada exógena, a operar em um ponto intermediário entre dois pontos mal es-
colhidos, isto é, para aqueles onde os modelos linearizados são válidos para vizinhanças
muito restritas. Neste caso, o controlador poderá apresentar uma dinâmica completa-
mente incompatível com as especificações de projeto, podendo conduzir a uma perda de
estabilidade ou de desempenho.
A técnica de escolha do conjunto de modelos lineares de projeto aqui apresentada
75
se baseia na idéia de estabelecer faixas de operação onde se possa garantir um erro
máximo entre a dinâmica do modelo linearizado e a dinâmica do modelo não-linear na
faixa considerada. Ou seja, busca-se particionar o domínio global de operação em tantos
subintervalos de interpolação quantos sejam necessários para assegurar um erro máximo
estipulado. Com isto, contorna-se o problema da validade da linearização para vizin-
hanças estreitas de um ponto de operação, uma vez que toda a faixa para o qual o
modelo linear é calculado está dentro do erro previamente adotado.
Para ilustrar a técnica sugerida nesta seção e facilitar a compreensão da mesma,
considera-se um exemplo de uma função escalar, definida por f(x) = −x3− 5x2− 3x− 5,
para x ∈ R. A linearização de f(x) no ponto de operação x = 6 é dada por:
∂f(x)∂x
|x=6= [−3x2 − 10x− 3] |x=6=
= −3.62 − 10.6− 3 = −171 (4.2)
A tangente à curva f(x) para x = 6 tem o coeficiente angular igual a −171. A FIG.
4.1 mostra a curva f(x) e sua tangente, conforme a linearização (4.2).
0 2 4 6 8 10−2000
−1500
−1000
−500
0
500
1000
x
f(x)
Função não−linearLinearização Ponto
FIG. 4.1: Curva f(x) e sua linearização em x = 6.
Porém, existe a possibilidade de obtenção de uma função linear que aproxime não só
um ponto da curva não-linear f(x), mas toda uma faixa do domínio de x previamente
selecionada. Para tanto, é necessário calcular uma função f(x) = ax + b, f(x) : R→ R,para um dado intervalo no domíno, que minimize
∫(|f(x)− f(x)|2dx)
76
Uma solução simples é discretizar f(x) na faixa de interesse em m pontos, obtendo um
vetor de amostradas fm ∈ Rm.
Assim, deseja-se obter o coeficiente angular de f(x) que atenda:
fm = MΘ (4.3)
onde M ∈ Rm×2 e Θ = [a b]T ∈ R2. Assumindo que m > 2, o problema é convexo e tem
a solução dada pela teoria dos mínimos quadrados. Neste caso, o vetor fm contém os
pontos da equação não-linear f(x) para os valores discretos de x ∈ [x1, xm], onde x1 e xm
são, respectivamente, os valores inicial e final de x para o intervalo considerado. A matriz
M tem os seus elementos da primeira coluna M(1c) composto pelos valores discretos de
x no intervalo x ∈ [x1, xm], ou seja, M(1c) = [x1 x2 x3 ... xm]T . A ordenação seqüencial
de j ∈ {1, 2, ..., m} deve ser a mesma para a primeira coluna da matriz M e para o vetor
fm, isto é, o elemento xj ∈ M(1c) deverá estar na mesma linha que f(xj) ∈ fm em (4.3).
A segunda coluna de M é um vetor cujos elementos são iguais a 1.
Então, a solução Θ = Θ∗ que minimiza o erro médio quadrático para m > 2 é:
Θ∗ = M †fm (4.4)
Para o emprego da técnica apresentada na Seção (2.5), sugere-se aqui utilizar um
conjunto de modelos linearizados sobre uma família pré-definida de intervalos de opera-
ção, onde os coeficientes angulares são calculados conforme (4.4) em substituição àqueles
calculados pela forma Jacobiana tradicional.
No exemplo escalar apresentado, f(x) = −x3 − 5x2 − 3x − 5, a fim de comparar
as aproximações, arbitra-se x0 = 4 e xm = 8. Discretizando o valor de x no intervalo
de interesse, com taxa de amostragem de 0, 01, e considerando os valores de f(x) cor-
respondentes para esta discretização, obtém-se: a = −173, 4120 e b = 590, 6519. Na
linearização, o que interessa é o coeficiente angular, ou seja, a = −173, 4120, pois o valor
b é o deslocamento representativo da faixa de operação.
Novamente é apresentada a curva f(x) e sua tangente calculada para x = 6 em (4.2)
na FIG. 4.2 (linha contínua e tracejada, respectivamente). A linearização considerando o
coeficiente ponderado por todos os valores da faixa pré-selecionada, calculado conforme
(4.4), é representada pela linha com ponto e traço. Observa-se na FIG. 4.3 que o maior
erro em relação à curva não-linear, no intervalo representado, é menor para a linearização
proposta.
Esta técnica de linearização por faixa pode ser facilmente generalizada para o caso
onde o estado é multivariável. Contudo, resta ainda definir um processo sistemático de77
escolha dos intervalos de linearização.
4 4.5 5 5.5 6 6.5 7 7.5 8−900
−800
−700
−600
−500
−400
−300
−200
−100
x
f(x)
Função não−linearLinearização PontoLinearização Faixa
FIG. 4.2: Curva f(x) e suas linearizações em x = 6: intervalo 4 ≤ x ≤ 8.
4 4.5 5 5.5 6 6.5 7 7.5 80
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
x
Qua
drad
o do
Err
o
Linearização FaixaLinearização Ponto
FIG. 4.3: Erros nos processos de linearização.
4.2.1 DEFINIÇÃO DAS FAIXAS DE LINEARIZAÇÃO
Nesta seção é apresentado um algoritmo para a escolha sistemática dos intervalos de
linearização. O processo busca particionar o domínio de operação do sistema não-linear
em vários intervalos, de maneira que o maior erro em cada intervalo não ultrapasse um
valor máximo especificado. Como sugestão, pode-se adotar um valor máximo de uma
norma (por exemplo, H2 ou H∞) da diferença entre o modelo linearizado na faixa e os
modelos linearizados ponto a ponto. A linearização Jacobiana (ponto a ponto) do sistema
78
é executada sobre uma grade fina de pontos que abranja todo o domínio de operação.
A suposição inicial de que um único intervalo, igual a este domínio atende o critério
estabelecido (norma da diferença), é tomada como ponto de partida do algoritmo.
Caso o sistema linear calculado nesse intervalo não atenda ao critério, reduz-se a
faixa de linearização até obter-se a satisfação do valor máximo adotado para a norma
das diferenças. A não convergência do algoritmo implica em relaxamento da norma, isto
é, deve-se aumentar o valor arbitrado como tolerância.
Algoritmo 4.1. Faixas de Linearização de um Sistema Não-Linear
Etapa 1: Inicialização
• Escolher a norma que vai ser adotada como critério para verificação da pro-
ximidade da dinâmica entre os sistemas comparados.
• Arbitrar o valor máximo de referência Vref que será admitido para a norma
adotada.
• Fornecer os vetores M(1c) com pontos de linearização para cada estado, bem
como o conjunto de modelos lineares ou a equação geral que possa extrair estes.
• Assumir todo o domínio de operação como sendo um intervalo único inicial.
Etapa 2: Linearização na faixa
• Calcular Θ para cada estado que define a operação, de acordo com (4.4), onde
fm contém os valores calculados a partir de f(x), considerando a dada faixa
de linearização sob análise.
Etapa 3: Verificação dos erros
• Calcular a norma da diferença entre o modelo linearizado na faixa e os modelos
linearizados ponto a ponto definidos pelo vetor M(1c) e guardar a maior destas
diferenças Vmax.
• Caso Vref > Vmax, a largura da faixa sob análise é adequada e deve-se definir
os limites da faixa seguinte, retomando-se a Etapa 2. Caso contrário, a largura
deverá ser diminuída a partir da redução do seu limite superior e a Etapa 2 re-
tomada. Sugere-se aqui adotar um algoritmo do tipo bisseção até se encontrar
um intervalo que atenda o critério de erro. A nova faixa deve ser escolhida
79
como o complemento da união dos intervalos já escolhidos e que atendem esse
critério.
Etapa 4: Solução
• A solução é obtida quando todos os intervalos estiverem definidos e atenderem
o critério de erro.
Deve-se ressaltar que as faixas são escolhidas para cada estado que define a operação
do sistema não-linear.
4.2.2 TRANSIÇÃO ENTRE FAIXAS
Com a propriedade de maior aproximação na técnica de linearização proposta na
seção anterior, pode-se assumir que a norma do erro entre a dinâmica dos modelos não-
linear e linearizados por faixas é inferior à tolerância estabelecida. Esta garantia não
existe no processo de linearização clássica ou Jacobiana, a menos que um número muito
elevado de pontos seja utilizado, o que poderia inviabilizar a implementação prática do
controlador escalonado.
A proximidade das dinâmicas dos modelos linearizados por faixas e do modelo não-
linear não é, no entanto, suficiente para a assegurar um bom comportamento em malha
fechada. É necessário que haja também uma transição suave entre os controladores
projetados para cada faixa. A fim de minimizar a introdução de dinâmicas rápidas
indesejáveis ao sistema em malha fechada, propõe-se aqui uma interpolação particular à
metodologia apresentada na Seção 4.2.1.
Para cada modelo linear particular f(x) calculado especificamente para aproximar a
função não-linear f(x) em determinado intervalo do domínio de x, x0 ≤ x ≤ xm, devem
ser levantados os valores x = x desta variável que tornam verdadeira a expressão:
f(x)|[x0 , xm] = f(x) = a.x + b (4.5)
Como f(x) representa uma curva, o universo de soluções em x de (4.5) apresenta
pelo menos dois valores distintos (ver o exemplo FIG. 4.2). Caso seja obtido mais de dois
valores na solução dessa equação, devem ser considerados somente os valores extremos,
aqui denominados por xi0 e xim, ou seja, os valores que mais se aproximam de x0 e xm
para uma determinada faixa i. Considerando três intervalos adjacentes de linearização,
i ∈ {j − 1, j, j + 1}, j ∈ N, calcula-se, então, os respectivos controladores. Nos intervalos80
em que a planta operar dentro dos limites de [xi0, xim] se propõe adotar um controlador
fixo, ou seja, a interpolação não se faz necessária.
A interpolação deverá ocorrer quando a planta estiver operando dentro dos seguintes
intervalos de domínio de x: [x(j−1)m, xj0] e [xjm, x(j+1)0], conforme ilustra a FIG. 4.4. Um
método numérico pode ser incorporado ao Algoritmo 4.1 para que este forneça também
os valores de xi0 e xim de cada faixa calculada.
FIG. 4.4: Transição entre faixas.
4.3 ESCALONAMENTO DE GANHOS POR IMPOSIÇÃO DA DINÂMICA DE
MALHA FECHADA
Nesta seção, é desenvolvido um estudo voltado ao controle de sistemas não-
estacionários utilizando controladores de estrutura estimação-controle. As matrizes que
definem a dinâmica do sistema não-linear, assim como os ganhos de realimentação e de
estimação de estados, apresentam uma dependência quasi -LPV em relação aos parâme-
tros de interpolação. Estes são endógenos e podem também incluir uma parte exógena,
como mostra a Seção 2.4.1. O método proposto não é geral, sendo aplicável somente a
uma classe de sistemas não-estacionários. Contudo, uma grande quantidade de plantas
de natureza prática se enquadra nessa classe.
81
Teoricamente, não se pode definir pólos e zeros para sistemas não-estacionários. Por
esta razão, propõe-se impor a dinâmica de malha fechada através da alocação dos au-
tovalores de realimentação e de estimação dos estados para toda trajetória possível dos
parâmetros. Busca-se, assim, restringir as possíveis perdas de desempenho quando a
planta evolui em um largo domínio de operação. O objetivo desta metodologia é, então,
proporcionar uma resposta desejada em malha fechada independentemente do ponto ou
faixa de operação.
Uma característica relevante da metodologia de projeto de controladores escalonados
aqui apresentada é o fato de que a lei de interpolação (ou de evolução, neste caso) dos
ganhos é determinada pela própria dinâmica do modelo quasi -LPV da planta, através da
alocação do espectro de autovalores de malha fechada para o conjunto contínuo de mode-
los LTI gerados pela fixação do parâmetro. Desta forma, para variações suficientemente
lentas do parâmetro, a dinâmica não-estacionária do sistema em malha fechada é muito
próxima da desejada.
4.3.1 CONSIDERAÇÕES SOBRE A MODELAGEM NÃO-LINEAR
Este estudo trata das dinâmicas que podem ser modeladas, ainda que pelo uso de
uma aproximação por série de Taylor de ordem superior, como:
x = A(x, u)x + W (x, u) + B(x, u)u (4.6)
e a saída dada por
y = C(x, u)x (4.7)
onde x(t) ∈ Rn, u(t) ∈ Rm é o vetor de de entradas, y(t) ∈ Rp é o vetor de saídas, e
A(x, u), B(x, u) e C(x, u) são matrizes reais de dimensões compatíveis e ainda:
W (x, u) = Wx(x, u)x = Wu(x, u)u (4.8)
onde as matrizes W (x, u), Wx(x, u) e Wu(x, u) também são de dimensões compatíveis.
É importante notar que a matriz W (x, u) representa uma parte das equações da
dinâmica do sistema que pode ser reescrita como Wx(x, u)x ou Wu(x, u)u, de forma
que possa vir a ser somada à A(x, u), à B(x, u), ou ainda parcialmente a ambas as
matrizes sem que haja alteração da dinâmica do modelo. Para facilitar a exposição do
método, considera-se os casos extremos onde a matriz W (x, u) é incorporada a A(x, u)
ou a B(x, u), sem perda de generalidade, usando a notação:
F (x, u) = A(x, u) + Wx(x, u) e E(x, u) = B(x, u) + Wu(x, u)
82
Assim sendo, os modelos:
x = F (x, u)x + B(x, u)u (4.9)
e
x = A(x, u)x + E(x, u)u (4.10)
são equivalentes e representam a mesma dinâmica representada em (4.6).
4.3.2 CONTROLADOR E OBSERVADOR
Uma vez fixado o modelo no qual se pretende trabalhar, (4.9) ou (4.10), agora se
verifica a possibilidade de controlar estes modelos sem a necessidade de os linearizar. Em
outras palavras, buscam-se ganhos não-lineares de realimentação e estimação de estados,
cuja dinâmica varie de forma a manter os autovalores de malha fechada alocados segundo
um determinado critério de desempenho. Embora a dedução dos ganhos possa apresentar
uma maior complexidade matemática incluindo termos cruzados que não aparecem em
malha aberta, tanto nos estados quanto nas entradas, os procedimentos para o seu cálculo
seguem rigorosamente os passos do cálculo de alocação de pólos, bem conhecidos das
técnicas lineares.
Considerando uma realimentação de estados com um ganho dinâmico, a matriz
dinâmica de malha fechada é:
Ac(x, u) = A(x, u)− E(x, u)Kc1(x, u) (4.11)
ou
Fc(x, u) = F (x, u)−B(x, u)Kc2(x, u) (4.12)
onde Kc1 e Kc2 tem dimensões compatíveis com as do modelo quasi -LPV do sistema.
Neste contexto, é possível se empregar todos os cálculos demonstrados em
(KAILATH, 1980) para alocação dos pólos em malha fechada nos valores desejados,
embora a complexidade matemática seja maior. Com isto, a dependência da dinâmica
em malha fechada das variáveis de entrada e dos estados é tão mais atenuada quanto
mais lenta for a variação desses parâmetros. O comportamento do modelo passa a ser
similar ao de um sistema LTI cujos pólos estejam localizados nos valores desejados.
Um raciocínio análogo pode ser empregado para o caso do observador de estados:
Af (x, u) = A(x, u)−Kf1(x, u)C(x, u) (4.13)
83
ou
Ff (x, u) = F (x, u)−Kf2(x, u)C(x, u) (4.14)
onde Kf1 e Kf2 tem as dimensões compatíveis com sistema. Isso implica em um obser-
vador de estados para plantas não-estacionárias.
Combinando (4.11) com (4.13), ou (4.12) com (4.14), e considerando o modelo não-
linear como quasi -LPV, dispõe-se de todos os dados para o cálculo de um controlador do
tipo estimação-controle para a planta não-estacionária.
4.3.3 CONTROLABILIDADE E OBSERVABILIDADE
Embora os conceitos de controlabilidade e observabilidade de sistemas variantes no
tempo estejam relacionados com as propriedades dos seus respectivos Gramianos, a
hipótese de variação lenta dos parâmetros é considerada válida para aplicação da téc-
nica proposta. Então, pode-se avaliar a inversibilidade das matrizes de controlabilidade
e de observabilidade para cada ponto de uma dada trajetória dos parâmetros, para se
inferir sobre a controlabilidade e a observabilidade do sistema quasi -LPV. Observa-se
que, neste caso, o parâmetro é fortemente dependente dos estados e/ou das entradas, o
que em princípio não invalida a hipótese de variação lenta no contexto quasi -LPV.
Uma vez definido o domínio de operação da planta, pode-se definir as matrizes de
controlabilidade e observabilidade, C(x, u) e O(x, u), respectivamente. Logo, para que o
sistema seja instantaneamente controlável e observável em todo o domínio de variação
dos seus estados e possíveis entradas é necessário que existam as inversas dessas matrizes
para todo esse domínio contínuo. As matrizes C(x, u) e O(x, u) são calculadas segundo
as técnicas aplicáveis a sistemas LTI (KAILATH, 1980).
É mister citar que a escolha de como distribuir as equações que compõem a matriz
W (x, u) (4.8), segundo (4.9) ou (4.10), pode influir fortemente na condição de controla-
bilidade e/ou observabilidade instantâneas dos modelos adotados. Portanto, essa escolha
deve ser criteriosa no sentido de se evitar uma modelagem particular da planta que resulte
na inexistência das inversas de C(x, u), O(x, u) ou no seu mau condicionamento.
Além disso, esta escolha também repercute sobre a amplitude do sinal de controle
incidente na planta. De forma geral, uma vez fixados os autovalores desejados em malha
fechada, quanto maior for a distância dos autovalores do spec[A(x, u)] ou do spec[F (x, u)],
maior será a amplitude do sinal de controle para um dado sinal de entrada de referên-
cia. Em virtude desta observação, é interessante distribuir as equações que compõem
W (x, u) de forma que a matriz dinâmica em malha aberta tenha as trajetórias dos au-
84
tovalores o mais próxima possível dos autovalores desejados de malha fechada, caso isso
não comprometa a existência de C(x, u)−1 e O(x, u)−1.
4.3.4 AJUSTE DO RASTREAMENTO E DO PONTO DE OPERAÇÃO
Ajuste do Rastreamento
Para se obter o ajuste do rastreamento, deve-se proceder de forma similar à usada
para modelos LTI contínuos. Um cuidado especial deverá ser considerado quando existir
dependência das matrizes que compõem o cálculo do ganho de rastreamento (do bloco
em cascata da entrada), aqui dependentes do parâmetro, em relação à resposta: o valor
da saída deverá ser substituído pelo sinal de referência r a ser seguido, o que implica em
se calcular este ajuste como:
{C(x, r) [−A(x, r) + B(x, r)K(x, r)]−1 B(x, r)
}−1(4.15)
Esse procedimento é fundamental, pois a flutuação provocada pela mudança de opera-
ção da planta imposta por uma entrada impõe um ajuste na saída que introduz dinâmicas
indesejáveis ao sistema em malha fechada caso não se tome a referência r no cálculo desse
ganho variante.
Ajuste do Ponto de Operação
Em alguns casos, para o emprego desta metodologia, o modelo quasi -LPV definido
em (4.6) e (4.7) é uma aproximação polinomial. Esta aproximação poderá provocar um
desajuste em relação a um ponto de operação de referência, no caso do problema da
regulação. Como exemplo, para uma entrada nula, um dado modelo apresenta sua saída
em malha fechada diferente de zero, embora próxima deste valor. Neste caso, quando for
julgado necessário, um pequeno ajuste poderá ser feito.
a) Em regime estacionário no ponto de operação que se deseja o ajuste, verificar a
diferença dos valores de saída (ds) e comparar com o valor da referência (dr);
b) Somar à entrada a constante ds − dr. Não se espera que este ajuste resolva o
problema, pois está se trabalhando exatamente com valores nos quais a aproximação
feita é deficiente.
c) Simular e repetir o passo "b" até obter ds∼= dr.
85
4.4 ESCALONAMENTO DE GANHOS POR CANCELAMENTO DE NÃO-
LINEARIDADES
Nesta seção, é proposto um método que objetiva o cancelamento das não-linearidades
através da realimentação de estados. Assim como no caso anterior, é utilizada uma técnica
de alocação de autovalores em malha fechada. Basicamente, a lei de controle proposta
busca anular os efeitos das dinâmicas indesejadas através da reprodução da variação
paramétrica no controlador e impõe o controle sobre a parte linear do modelo. Esse
método também objetiva impedir que variações bruscas sejam inseridas artificialmente na
dinâmica de malha fechada, pois além de anular as não-linearidades através de variações
no controlador, o ganho permanece fixo.
Trabalhos como (ISIDORI, 1981, 1986, 1989; SCHOENWALD, 1992; BANASZUK,
1994; REILLY, 1996; BATTILOTTI, 2004) abordam a linearização por realimentação em
plantas do tipo:x = f(x) +
∑gi(x).ui
y1 = h1(x)...
ym = hm(x)
(4.16)
onde x = (x1, x2, . . . , xn) ∈ Rn são estados, ui ∈ R entradas e as saídas medidas são
yi ∈ R.Aplicando em (4.16) uma mudança de coordenadas do tipo:
z = Φ(x) (4.17)
para
Φ(x) =
φ1(x)
φ2(x)...
φn(x)
(4.18)
os trabalhos mencionados estudam e demonstram como tais sistemas podem ser linea-
rizados em torno de um ponto de operação. A técnica aqui apresentada difere desses
trabalhos mencionados. A abordagem enfatizada neste estudo é particular a uma classe
de plantas que possua um modelo no qual existe uma parte da dinâmica linear e outra
parte não-linear. Sob certas condições, é possível obter um sinal de controle que anule a
dinâmica não-linear, fazendo com que o sistema em malha fechada tenha um comporta-
mento linear.86
4.4.1 MODELAGEM
Este estudo foca modelos descritos que possam ser definidos na forma:
x = Ax + Ww(x, u) + Bu (4.19)
cuja saída é dada por
y = Cx + Du (4.20)
onde x(t) ∈ Rn, u(t) ∈ Rm é o vetor de de entradas, y(t) ∈ Rp é o vetor de saídas, w(.)
possui n linhas e 1 coluna, e A, W , B, C, D são matrizes reais de dimensões compatíveis.
Não são todos os modelos que possibilitam essa forma de equacionamento e que
podem ser controlados pela técnica proposta nesta seção, como é visto a seguir.
4.4.2 CANCELAMENTO DAS NÃO-LINEARIDADES
Neste estudo, há necessidade de que a tripla (A,B, C) seja controlável e observável.
Pode-se dizer que, para W = 0, o modelo LTI (A,B, C) é um caso particular do modelo
não-linear (4.19). Busca-se verificar se é possível anular as não-estacionariedades (w(.)),
somando ao sinal de controle da parte linear um sinal adicional, de forma que a equação
dinâmica de malha fechada seja da seguinte forma:
x = Ax + Ww(x, u) + B[r −K1x−K2w(x, u)] (4.21)
onde r é a entrada de referência, K1x + K2w(x, u) é o sinal de controle, com Ki ∈ Rm×n,
para i = 1, 2. Reagrupando os termos:
x = (A−BK1)x + (W −BK2)w(x, u) + Br (4.22)
Quando:
W = BK2 (4.23)
os efeitos das não-linearidades são anulados pelo sinal suplementar K2w(x, u). A solução
óbvia aponta para uma matriz B inversível.
Ainda que B não seja inversível, para solucionar (4.23), existe K2 que permite can-
celar as não-linearidades, quando as entradas forem desconectadas para cada efeito não-
linear ou parâmetro variável existente na dinâmica de um determinado estado xi. Em
outra palavras, se a dinâmica de xi é de natureza não-linear, quando existir uma entrada
ui associada ao estado xi exclusivamente, pode-se fazer com que a planta se comporte
como um sistema linear A, B, C e D via realimentação. Essa análise é conduzida caso87
a caso e é de fácil constatação, como ilustra o exemplo da Seção 5.5. Com isto a não-
linearidade é anulada, sendo possível aplicar todo o amplo ferramental disponibilizado
para as plantas lineares neste modelo.
88
5 APLICAÇÕES
5.1 IDENTIFICAÇÃO DE MODELOS LINEARES USANDO A RESPOSTA DIS-
CRETA NO TEMPO
Nos exemplos desta seção foi utilizado o sinal de entrada discreto v =
{8, 4, 1, 10, 12, 30,−1, 0, 10, 10,−20, 1, 5,−10,−10, 0, 0, 8,−3, 1, 2, 0} apresentado na FIG.
5.1. A largura dos pulsos de entrada para o conjunto v é de 1 segundo. O período de
amostragem usado para a identificação e o tempo total que a planta foi submetida ao
sinal de excitação são particulares a cada exemplo, estando nesses especificados.
Para a validação da identificação é aplicada à planta e ao modelo discreto obtido
uma entrada do tipo degrau unitário. Num tempo total de medição de 10 segundos, com
período de amostragem de 5 milisegundos, é avaliada a diferença entre os dois sinais de
saída com base na norma 2.
A taxa de amostragem de 0, 005 segundos, arbitrada, também é usada em todos os
exemplos desta seção para a medição dos dados de entrada e saída da planta identificada.
Os modelos equivalentes discretos deste capítulo foram obtidos utilizando a técnica
ZOH.
0 500 1000 1500 2000 2500 3000 3500−40
−30
−20
−10
0
10
20
30
40
Tempo (s)
FIG. 5.1: Sinal v aplicado as plantas físicas
89
Exemplo 1
Este exemplo é de uma função transferência de segunda ordem, estritamente própria.
A planta pode ser modelada como:
GM1(s) =10s + 10
s2 + 5s + 6(5.1)
A FT discreta equivalente à 5.1, para um período de amostragem de 0, 005s, é:
GM1d(z) =0, 0495z − 0, 0493
z2 − 1, 9752z + 0, 9753(5.2)
O correspondente diagrama de Bode da FT apresentada é:
10−2
10−1
100
101
102
0
0.5
1
1.5
2
2.5
|GM
1(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
−100
−80
−60
−40
−20
0
Fas
e (g
raus
)
Freqüência (rad/s)
FIG. 5.2: Resposta em freqüência (módulo e fase) da dinâmica da planta
Na TAB. 5.1, a primeira coluna traz a denominação do modelo calculado para um
período de amostragem de 0, 005 segundos: Gm(z) é o modelo identificado diretamente
pela metodologia proposta e Gr(z) é o equivalente modelo reduzido e/ou a projeção
estável (quando foi necessária o uso dessas técnicas). O tempo que o sinal de entrada é
aplicado está na segunda coluna. A terceira coluna mostra o ganho Kv ao qual a entrada
é multiplicada (Kv × v). A quarta coluna traz o custo respectivo. Os custos obtidos
são relativos ao sinal de validação, degrau unitário, imposto por 10 segundos ao modelo
tomado como planta e ao modelo identificado. Esta fato defini a quantidade de amostras
dos vetores utilizado no cálculo da função custo: 2001 amostras.
90
TAB. 5.1: Dados referentes aos modelos identificados de GM1(s)
Modelo Tempo (s) Kv Custo
Gm11(z) 10 0, 1 188, 2248
Gm12(z) 50 0, 1 128, 9441
Gm13(z) 50 1, 5 0, 9179
Gm14(z) 50 5 0, 0828
Gr15(z) 10 0, 1 83, 6623
Gr16(z) 50 0, 1 37, 0100
Gm17(z) 50 1, 5 0,1843
Gm18(z) 50 5 0, 01660
As FT discretas, relacionados na primeira coluna, são:
• Gm11(z) = 0,0495z−0,0486
z2−1,9610z+0,9615
• Gm12(z) = 0,0495z−0,0491
z2−1,9713z+0,9715
• Gm13(z) = 0,04950291940403z−0,04925515277825
z2−1,97514455776172z+0,97529312533261
• Gm14(z) = 0,04950290604650z−0,04925593005728
z2−1,97516022393298z+0,97530840121001
• Gr15(z) = 0,04950150028922z−0,04911644336826
z2−1,97229907240304z+0,97251800107981
• Gr16(z) = 0,04950354453251z−0,04921746938149
z2−1,97438668370996z+0,97455416477862
• Gm17(z) = 0,04950290757149z2−0,01626168765238z−0,03282947211210
z3−1,30864904330177z2−0,34115480929911z+0,65005087042431
• Gm18(z) = 0,04950290498158z2−0,01626181986433z−0,03282960139235
z3−1,30865165700456z2−0,34115487432492z+0,65005341878700
Os modelos identificados ou sua respectiva redução que não constam na TAB. 5.1
são apresentado no Apêndice 8.2 deste trabalho.
A saída medida da planta, a resposta do modelo e o erro entre estas, para um mesmo
sinal de entrada, dos modelos Gm1 1(z), Gm
12(z), Gm13(z) e Gm
14(z) estão na FIG. 5.3. A
entrada aplicada é o degrau unitário no primeiro segundo. Para os modelos Gr15(z),
Gr16(z), Gm
17(z) e Gm18(z), as curvas citadas estão na FIG. 5.4.
91
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
da
PlantaModelo
0 2 4 6 8 10−6
−4
−2
0
2
4
Tempo (s)
Err
o
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
da
PlantaModelo
0 2 4 6 8 10−4
−2
0
2
Tempo (s)
Err
o
(A) (B)
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
da
PlantaModelo
0 2 4 6 8 10−0.03
−0.02
−0.01
0
0.01
Tempo (s)
Err
o
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
daPlantaModelo
0 2 4 6 8 10−3
−2
−1
0
1x 10
−3
Tempo (s)
Err
o
(C) (D)
FIG. 5.3: Curvas de respostas ao degrau para: planta e modelo. Curva de erro entre asrespostas do modelo e da planta: (A) Gm
11(z), (B) Gm12(z), (C) Gm
13(z) e (D) Gm14(z).
92
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
da
PlantaModelo
0 2 4 6 8 10−3
−2
−1
0
1
2
Tempo (s)
Err
o
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
da
PlantaModelo
0 2 4 6 8 10−1.5
−1
−0.5
0
0.5
Tempo (s)
Err
o
(A) (B)
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
da
PlantaModelo
0 2 4 6 8 10−6
−4
−2
0
2x 10
−3
Tempo (s)
Err
o
0 2 4 6 8 100
20
40
60
Tempo (s)
Saí
daPlantaModelo
0 2 4 6 8 10−6
−4
−2
0
2x 10
−4
Tempo (s)
Err
o
(C) (D)
FIG. 5.4: Curvas de respostas ao degrau para: planta e modelo. Curva de erro entre asrespostas do modelo e da planta: (A) Gr
15(z), (B) Gr16(z), (C) Gm
17(z) e (D) Gm18(z).
93
Exemplo 2
Este exemplo é de uma função transferência de sétima ordem estritamente própria.
A planta contínua é descrita na TAB. 5.2 e seu equivalente discreto na TAB. 5.3.
TAB. 5.2: Valores dos coeficientes da FT contínua da planta utilizada
Graus em s Coeficientes do numerador Coeficientes do denominador
s7 0,1 1s6 4,753 16,82s5 90,01275 252,6691s4 863,2453825 2476,872788s3 4417,86362625 13659,19997939s2 11761,8216175 59044,409874979s1 14942,80380125 131354,017942688s0 6958,3200225 87580,350142506
TAB. 5.3: Valores dos coeficientes da FT discreta da planta utilizada
Graus em z Coeficientes do numerador Coeficientes do denominador
z7 0,1 1z6 -0,67580254321921 -6,91313222834198z5 1,95709308517528 20,48530559630960z4 -3,14831939251784 -33,72923096678250z3 3,03840030958530 33,32652110431409z2 -1,75917724122000 -19,76024986739179z1 0,56578095741999 6,51012567946567z0 -0,07797517522301 -0,91933931756652
94
O correspondente diagrama de Bode da FT apresentada é:
10−2
10−1
100
101
102
0
0.5
1
1.5
|G(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
−100
−50
0
50
100
Fas
e (g
raus
)
Freqüência (rad/s)
FIG. 5.5: Resposta em freqüência (módulo e fase) da dinâmica da planta
Na TAB. 5.4, a primeira coluna traz a enumeração das identificações realizadas para
a planta descrita. O tempo que o sinal de entrada é aplicado está na segunda coluna. A
terceira coluna mostra o sinal de identificação aplicado (Si): v para o conjunto de mesmo
nome já definido e d para o degrau unitário aplicado no primeiro segundo de medição.
A quarta coluna apresenta o sinal de validação do modelo identificado (Svd): v ou d.
Na quinta coluna encontra-se o ganho Kv ao qual o sinal v é multiplicado (Kv × v). A
sexta coluna traz a função de transferência discreta identificada, para um período de
amostragem τ . A função de transferência discreta de ordem reduzida, quando é possível,
está na sétima coluna. O número de amostras (Na) usados para o cálculo do valor da
função custo está apresentado na oitava coluna. Na nona coluna encontra-se o período
de amostragem τ em segundos.
TAB. 5.4: Dados referentes aos modelos identificados
Modelo Tempo (s) Si Svd Kv Gm(z) Gr(z) Na τ(s)
1.1 10 v d 0, 1 Gm21(z) Gr
21(z) 2001 0, 005
1.2 50 v d 5.109 Gm22(z) Gr
22(z) 10001 0, 005
1.3 50 v d 5 Gm23(z) Gr
23(z) 10001 0, 005
1.4 10 d v 0, 1 Gm24(z) Gr
24(z) 101 0, 1
95
Na TAB. 5.5, a primeira coluna traz o modelo identificado ou a redução do mesmo.
A ordem do modelo está na segunda coluna. A terceira coluna mostra o custo obtido
para o sinal de validação.
TAB. 5.5: Custos referentes aos modelos identificados
Gm(z) Ordem Custo
Gm21(z) 7 48, 5521
Gr21(z) 3 48, 4820
Gm22(z) 7 39, 9053
Gr22(z) 4 39.9048
Gm23(z) 15 1, 7758
Gr23(z) 4 1, 7740
Gm24(z) 4 0, 0082
Gr24(z) 4 0, 0082
Os únicos vetores para o cálculo da função custo que apresentam uma mesma quan-
tidade de amostras são relativos aos modelos 1.2 e 1.3 (10001 amostras). Logo, os valores
da terceira coluna não deve ser tomados como termo de comparação para avaliação do
melhor ajuste entre os modelos identificados. Contudo, nos gráficos apresentados, as
diferenças nos ajustes ficam evidentes e podem ser comparados.
As FT discretas identificadas, indicadas na TAB. 5.4, estão relacionadas no Apêndice
8.2 deste trabalho.
A saída medida da planta, as respostas dos modelos 1.1, 1.2, 1.3 e o erro entre a
medida e a estimação, para o degrau unitário aplicado no primeiro segundo às entradas,
estão na FIG. 5.6.
96
0 2 4 6 8 10−5
0
5
10
15
Tempo (s)
Saí
da
Planta
Modelo
0 2 4 6 8 10−5
0
5
10
Tempo (s)
Err
o
(A)
0 2 4 6 8 10−5
0
5
10
15
Tempo (s)
Saí
da
Planta
Modelo
0 2 4 6 8 10−4
−2
0
2
4
Tempo (s)
Err
o
(B)
0 2 4 6 8 10−5
0
5
10
15
Tempo (s)
Saí
da
Planta
Modelo
0 2 4 6 8 10−0.1
−0.05
0
0.05
0.1
0.15
Tempo (s)
Err
o
(C)
FIG. 5.6: Curvas de respostas ao degrau para: planta e modelo. Curva de erro entre asrespostas do modelo e da planta: (A) Gr
21(z), (B) Gr22(z) e (C) Gr
23(z).
97
0 5 10 15 20−2
−1
0
1
Tempo (s)S
aída
PlantaModelo
0 5 10 15 20−4
−2
0
2x 10
−3
Tempo (s)
Err
o
FIG. 5.7: Curvas de respostas ao sinal v para: planta e modelo. Curva de erro entre asrespostas medida e estimada: Gr
24(z)
Os exemplos aqui expostos evidenciam três fatores que contribuem para uma iden-
tificação com melhor ajuste. Dois destes já foram abordados na Seção 3.3: multiplicar
a entrada por uma constante e aplicar um somatório de degraus (trem de pulsos). O
terceiro fator fica claro no modelo 1.3. Neste caso, o maior número de pólos e zeros
contribuiu decisivamente para a queda do valor do custo11.
Na análise da Seção 3.3 foi mostrado como o período de amostragem τ influi no
condicionamento numérico da solução da técnica de identificação aqui proposta. Em
poucas palavras, um τ muito pequeno implica em poucos algarismos armazenados com-
putacionalmente que diferencie as posições de pólos e zeros, o que repercute diretamente
no cálculo dos coeficientes da FT. Outro fator, além do aumento do grau de liberdade
da solução, pode ser apontado na contribuição de um ajuste de custo mais adequado.
Para ajudar a clarificar o que aconteceu no modelo 1.3, pode-se considerar uma pequena
operação literal com 3 pólos (ou zeros).
Para |ni| ¿ 1, i ∈ N e ni ∈ C:
[z − (1− n1)].[z − (1− n2)].[z − (1− n3)] =
= [z2 − (2− n1 − n2)z + (1− n1 − n2 + n1n2)].[z − (1− n3)] =
= z3 − (3− n4)z2 + (3− n5)z − (1− n6)
(5.3)
11Comparar com o modelo 1.2, onde a constante que multiplica o sinal de entrada é significativamente
maior, τ e a quantidade de amostras consideradas permanecem inalterados.
98
onde n4 = n1 + n2 + n3, n5 = 2n1 + 2n2 + 2n3− n1n2, n6 = n1 + n2 + n3− n1n2− n1n3−n2n3 + n1n2n3.
Neste exemplo simples, para três valores de pólos (ou zeros) de módulo aproximada-
mente igual a 1, é explicito que os coeficientes irão gradualmente tendo uma variação
de ni devido as somas e multiplicações. No modelo 1.3, estas alterações provocaram um
gradual aumento das posições relativas das casas decimais dos coeficientes da FT dis-
creta levantada (ordem 15). Em conseqüência, obteve-se uma melhor estimativa da FT
identificada.
5.2 IDENTIFICAÇÃO DE MODELOS LINEARES USANDO A RESPOSTA EM
FREQÜÊNCIA
Em todos os exemplos acadêmicos desta seção, os dados amostrados foram obtidos
usando as funções de transferência que representam as plantas a serem identificadas. Para
a geração desses dados foi utilizado um vetor de freqüência, com 200 pontos, com o li-
mite superior de 6, 28308943338034.102 rad/seg e limite inferior de 0, 00009587379924.102
rad/seg, espaçados logaritmamente (base 10).
As tabelas com os valores dos custos das funções identificadas apresentam também
os resultados obtidos pela metodologia desenvolvida em (VALLE, 2005), cujos modelos
são indicados pela referência GRV .
Como nos exemplo precedentes, Gmi (s) é o modelo identificado e Gr
i (s) é o modelo
reduzido, quando possível, fazendo a projeção de Gmi (s) em uma soma de funções de
transferência estável e instável, porém desprezando-se em seguida a parte instável nesta
soma.
Exemplo 3
Este exemplo é de uma FT de segunda ordem estritamente própria. A planta tem a
seguinte modelagem:
GM3(s) =10s + 10
s2 + 5s + 6(5.4)
A TAB. 5.6 indica os valores dos custos calculados relativos a cada modelo identifi-
cado.
99
TAB. 5.6: Valores da função custo da FT GM3(s)
Modelo Ordem Custo
GRV (s) 2 8, 02914859270.10−3
GRV (s) 3 3, 11240181717.10−3
Gm31(s) 2 8, 6158694479621.10−4
Gm32(s) 3 6, 0310179151198.10−4
Gr32(s) 2 5, 9070667856717.10−4
Os modelos identificados nesta metodologia, indicados na TAB. 5.6, estão caracteri-
zados pelas FT apresentadas a seguir:
• Gm31(s) = 10,00000011166216s+9,99622103937097
s2+4,99963355793166s+5,99799440945547
• Gm32(s) = 10,00000000036061s2+4,84639011640609s−5,15181126472957
s3+4,48463901743252s2+3,42336739684553s−3,09121075307096
• Gr32(s) = 9,99999892482393s+9,99881332710704
s2+4,99988193165595s+5,99952113408498
A FIG. 5.8 corresponde à resposta em freqüência do modelo identificado Gr32(s), bem
como a resposta em freqüência da dinâmica de (5.4).
10−2
10−1
100
101
102
0
0.5
1
1.5
2
2.5
|GM
3(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
−100
−80
−60
−40
−20
0
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
FIG. 5.8: Resposta em freqüência (módulo e fase): comparação entre a planta e omodelo Gr
32(s).
100
Exemplo 4
Este exemplo é de uma FT de segunda ordem estritamente própria. A planta tem a
seguinte modelagem:
GM4(s) =5
s2 + 2s + 2(5.5)
A TAB. 5.7 indica os valores dos custos calculados relativos a cada modelo identifi-
cado.
TAB. 5.7: Valores da função custo da FT GM4(s)
Modelo Ordem Custo
GRV (s) 2 1, 665092433651023.10−4
GRV (s) 3 1, 542191623292908.10−4
Gm41(s) 2 8, 222305128202556.10−5
Gm42(s) 3 7, 0287026427989.10−5
Gr42(s) 2 8, 6557550461710.10−5
Os modelos identificados nesta metodologia, indicados na TAB. 5.7, estão caracteri-
zados pelas FT apresentadas a seguir:
• Gm41(s) = 0,00000000495838s+4,99994016380416
s2+1,99996475536416s+1,99997627463195
• Gm42(s) = 0,00000000001826s2+4,99999999963070s−0,88235020688053
s3+1,82352919843854s2+1,64705061487628s−0,35293969689422
• Gm42(s) = −0,00000112543695s+4,99999755011403
s2+1,99999977646143s+1,99999173147407
A figura 5.9 corresponde à resposta em freqüência do modelo Gm42(s), bem como a
resposta em freqüência da dinâmica de (5.5).
101
10−2
10−1
100
101
102
0
1
2
3
|GM
4(jw
)|Freqüência (rad/s)
10−2
10−1
100
101
102
−200
−150
−100
−50
0
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
FIG. 5.9: Resposta em freqüência (módulo e fase): comparação entre a planta e omodelo Gm
42(s).
Exemplo 5
Este exemplo é de uma função de transferência de primeira ordem biprópria. A
planta tem a seguinte modelagem:
GM5(s) =3s + 3
s + 4(5.6)
A TAB. 5.8 indica os valores do custos calculados relativos a cada modelo identificado.
TAB. 5.8: Valores da função custo da FT GM5(s)
Modelo Ordem Custo
GRV (s) 1 6, 39286642765.10−3
GRV (s) 2 4, 29458769714844.10−4
Gm51(s) 1 1, 343410171025567.10−4
Gm52(s) 2 1, 0734589783276.10−4
Gr52(s) 1 2, 902658582329.10−5
Os modelos identificados nesta metodologia, indicados na TAB. 5.8, estão caracteri-
zados pelas FT apresentadas a seguir:
• Gm51(s) = 2,99999999108632s+2,99991937462496
s+3,99996301493381
• Gm52(s) = 2,99999999993276s2+1,88570122262619s−1,11426791862194
s2+3,62856707328157s−1,48571674632358
102
• Gr52(s) = 2,99999999993276s+2,99998589668518
s+3,99999657765266
A FIG. 5.10 corresponde à resposta em freqüência do modelo Gr52(s), bem como à
resposta em freqüência da dinâmica de (5.6).
10−2
10−1
100
101
102
1
2
3|G
M5(
jw)|
Freqüência (rad/s)
10−2
10−1
100
101
102
0
10
20
30
40
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
FIG. 5.10: Resposta em freqüência (módulo e fase): comparação entre a planta e omodelo Gr
52(s).
Exemplo 6
Este exemplo é de uma função de transferência de primeira ordem biprópria e com
fase não mínima. A planta tem a seguinte modelagem:
GM6(s) =3s− 3
s + 4(5.7)
A TAB. 5.9 indica os valores dos custos calculados relativos a cada modelo identifi-
cado.
TAB. 5.9: Valores da função custo da FT GM6(s)
Modelo Ordem Custo
GRV (s) 1 1, 065477737936.10−2
GRV (s) 2 7, 160967685960582.10−4
Gm61(s) 1 3, 343933313214873.10−5
Gm62(s) 2 2, 144420689634599.10−5
Os modelos identificados nesta metodologia, indicados na TAB. 5.9, estão caracteri-
zados pelas FT apresentadas a seguir:103
• Gm61(s) = 2,99999999861816s−2,99998660896572
s+3,99999840291869
• Gm62(s) = 3,00000000001550s2−2,57142791023567s−0,42857004759117
s2+4,14285736307514s+0,57142905621459
A figura 5.11 corresponde a resposta em freqüência do modelo Gm62(s), bem como a
resposta em freqüência da dinâmica de (5.7).
10−2
10−1
100
101
102
0
1
2
3
|GM
6(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
0
50
100
150
200
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
FIG. 5.11: Resposta em freqüência (módulo e fase): comparação entre a planta e omodelo Gm
62(s).
Exemplo 7
Este exemplo é de uma FT de segunda ordem estritamente própria. A planta tem a
seguinte modelagem:
GM7(s) =10s− 20
s2 + 5s + 6(5.8)
A TAB. 5.10 indica os valores do custos calculados relativos a cada modelo.
TAB. 5.10: Valores da função custo da FT GM7(s)
Modelo Ordem Custo
GRV (s) 2 1, 369916888852.10−2
GRV (s) 3 7, 101343929997529.10−4
Gm71(s) 2 1, 872075101508796.10−4
Gm72(s) 3 8, 513035968792272.10−5
Os modelos identificados nesta metodologia, indicados na TAB. 5.10, estão caracte-
rizados pelas FT apresentadas a seguir:104
• Gm71(s) = 10,00000002291290s−19,99982589474207
s2+5,00001165296947s+5,99996457647605
• Gm72(s) = 10,00000000013230s2−17,06403865886867s−5,87187378726888
s3+5,29359613312761s2+7,46798279990970s+1,76156678510554
A FIG. 5.12 corresponde à resposta em freqüência do modelo Gm72(s), bem como à
resposta em freqüência da dinâmica de (5.8).
10−2
10−1
100
101
102
0
1
2
3
|GM
7(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
−100
0
100
200
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
FIG. 5.12: Resposta em freqüência: comparação das amostras com os valores fornecidospelo modelo Gm
72(s).
Exemplo 8
Este exemplo é de uma função de transferência de sétima ordem biprópria. A mode-
lagem da planta apresenta os seguintes coeficientes em sua função transferência:
TAB. 5.11: Valores dos coeficientes da FT
Graus em s Coeficientes do numerador Coeficientes do denominador
s7 0,1 1s6 4,753 16,82s5 90,01275 252,6691s4 863,2453825 2476,872788s3 4417,86362625 13659,19997939s2 11761,8216175 59044,409874979s1 14942,80380125 131354,017942688s0 6958,3200225 87580,350142506
105
Para uma melhor compreensão das identificações realizadas e dos custos obtidos, é
interessante verificar os Valores Singulares de Hankel associados à dinâmica da FT a ser
identificada. Estes valores foram obtidos usando o programa hksv do MATLAB 7.
Vh =
0, 79721577798123
0, 72170563888440
0, 26836386727286
0, 18257863242876
0, 00003739624354
0, 00003103517113
0, 00000680622801
Dentro da técnica sugerida, são identificados três modelos com ordens 4, 5 e 7. Os
coeficientes das FT dos modelos aqui levantados estão no Apêndice 8.2 deste trabalho.
TAB. 5.12: Comparativo de custos
Modelo Ordem Custo
GRV (s) 5 0,05609556083041GRV (s) 7 0,01753429577575Gm
81(s) 4 0, 00114170230673
Gm82(s) 5 0,03273684616308
Gr82(s) 4 0,00401288963516
Gm83(s) 7 0,42484749560738
Gr83(s) 6 0,37731525807160
Neste exemplo se pode verificar que o melhor ajuste ocorre quando a quantidade
de pólos do modelo identificado foi arbitrado em 4, ordem que o modelo usado como
planta pode obter boa aproximação. Para um maior número de pólos arbitrados, mesmo
quando este número ficou limitado à quantidade de pólos da FT original (TAB. 5.12), o
algoritmo se mostrou sensível ao aumento do grau de liberdade, deteriorando o ajuste da
FT identificada.
106
10−2
10−1
100
101
102
0
0.5
1
1.5
|G(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
−100
−50
0
50
100
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
(A)
10−2
10−1
100
101
102
0
0.5
1
1.5
|G(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
−100
−50
0
50
100
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
(B)
10−2
10−1
100
101
102
0
0.5
1
1.5
|G(jw
)|
Freqüência (rad/s)
10−2
10−1
100
101
102
−100
−50
0
50
100
Fas
e (g
raus
)
Freqüência (rad/s)
ModeloAmostras
(C)
FIG. 5.13: Resposta em freqüência. Comparação das amostras com os valoresfornecidos pelos modelos (A) Gm
81(s), (B) Gr82(s) e (C) Gr
83(s).
107
5.3 IDENTIFICAÇÃO QUASI -LPV
O presente exemplo baseia-se no modelo não-linear de um míssil ar-ar, onde a veloci-
dade foi considerada constante (M=3). O modelo está no Apêndice 8.1 deste trabalho.
Considerando que já foi levantado o conjunto de modelos LTI contínuos de 29 pontos
de operação, adotou-se a técnica descrita na Seção 3.4.
O conjunto de modelos lineares foi obtido variando o ângulo de ataque de 1o até
29o, de grau em grau. Uma vez estabelecidos os coeficientes que caracterizam o nume-
rador e o denominador de cada ponto de operação, foi levantada uma curva polinomial
parametrizada em α para cada um dos coeficientes das potências em s que caracterizam
a função transferência dos modelos lineares.
Assumindo como modelo genérico linear, no ponto de operação i, para i ∈{1, 2, ..., 29}, a FT assume a forma:
Gi(s) =bi0s
n + bi1sn−1 + ... + bi(n−1)s + bin
sn + ai1sn−1 + ... + ai(n−1)s + ain
(5.9)
Essa modelagem pode ser generalizada e reescrita em forma de equações matriciais,
partindo-se das formas canônicas. Neste exemplo, foi adotada a forma canônica obser-
vador.
xi(t) =
−ai1 1 . . . 0 0
−ai2 0 . . . 0 0...
... . . . ......
−ai(n−1) 0 . . . 0 1
−ain 0 . . . 0 0
x(t) +
bi1
bi2
...
bi(n−1)
bin
u(t)
y(t) =[
1 0 . . . 0 0]x(t)
(5.10)
onde x(t) ∈ Rn, aij e bij ∈ R, u(t) é a entrada e y(t) a saída do sistema SISO. As matrizes
têm dimensões compatíveis.
Generalizando o problema dentro da técnica proposta, a equação geral da planta
quasi -LPV, parametrizada pelo ângulo de ataque α, apresenta a seguinte dinâmica:
108
x(t) =
−a1(α) 1 . . . 0 0
−a2(α) 0 . . . 0 0...
... . . . ......
−a(n−1)(α) 0 . . . 0 1
−an(α) 0 . . . 0 0
x(t) +
b1(α)
b2(α)...
b(n−1)(α)
bn(α)
u(t)
y(t) =[
1 0 . . . 0 0]x(t)
(5.11)
ou ainda:
x(t) = A(α)x(t) + B(α)u(t)
y(t) = Cx(t)(5.12)
Neste exemplo, n = 4 e os polinômios identificados que definem os termos aj(α) e
bj(α) das matrizes A(α) e B(α) são:
a1(α) = 3, 3664.10−2α + 210, 5696
a2(α) = 1, 3113.10−5α4 − 0, 001889α3 − 0, 232776α2 + 19, 4452α + 2, 2594.104
a3(α) = 0, 001405α4 − 0, 1167696α3 − 38, 9809α2 + 3, 3576.103α + 6, 7673.103
a4(α) = −3, 0731.103α2 + 2, 7873.105α− 7, 2898.105
b1(α) = 0
b2(α) = −4, 5850.103
b3(α) = 0
b4(α) = −1, 1915.103α2 + 1, 0997.105α + 3, 1416.106
(5.13)
Para validar o modelo proposto, foi usada uma entrada, conforme a FIG. 5.14, e
comparada a primeira saída da planta, com a correspondente do modelo identificado.
109
0 10 20 30 40 50 60 70−40
−30
−20
−10
0
10
20
30
40
Tempo(s)
Âng
ulo
Com
anda
do (
º)
FIG. 5.14: Entrada para validação do modelo.
Como resultado, as figuras 5.15 e 5.16 apresentam as saídas medida e estimada.
0 10 20 30 40 50 60 70−150
−100
−50
0
50
100
150
Tempo(s)
Planta
Modelo
FIG. 5.15: Resposta da planta não-linear e do modelo identificado quasi -LPV de ummíssil.
110
35 35.5 36 36.5 37 37.5 38 38.5 39 39.5 40
−120
−100
−80
−60
−40
−20
0
20
40
Tempo(s)
Planta
Modelo
FIG. 5.16: Resposta da planta não-linear e do modelo identificado quasi -LPV de ummíssil entre 35s e 45s de vôo.
As pequenas diferenças de amplitude entre as respostas do modelo e da planta na FIG.
5.16 são notórias, bem como algumas diferenças de fase. A natureza desta metodologia é
calcada diretamente nos modelos lineares identificados nos pontos de operação da planta.
O que se constata nas curvas das figuras citadas está compatível com as aproximações
feitas de modelos variantes no tempo para LTI.
As diferenças apontadas também aparecem, mesmo em pequenos intervalos de afas-
tamento do ponto de operação, quando se utiliza modelos lineares, em um determinado
ponto de operação, e se compara a saída deste modelo com a resposta da planta. Con-
tudo, tais diferenças não influem significativamente no controle da planta, conforme será
demonstrado em exemplo deste mesmo capítulo na Seção 5.4.
Para reproduzir os resultados aqui apresentados, deverão ser tomados os cuidados de
evitar que a planta do míssil opere com ângulo de ataque com valores modulares inferior
a 2o. Com ângulo nesta faixa os modelos linearizados são instáveis, o que prejudica a
eficiência do método proposto.
111
É possível tentar obter uma função de ajuste para minimizar os erros apontados.
Porém não se pode garantir que essa dinâmica de ajuste seja facilmente incorporada
dentro de alguma técnica de controle já estabelecida.
5.4 CONTROLE DE UM MÍSSIL AR-AR VIA CONTINUAÇÃO DE SUBESPAÇOS
INVARIANTES
Nesta seção, é exemplificada a aplicação da técnica de escalonamento de ganhos
apresentada na Seção 2.5, bem como comparados os resultados da interpolação de con-
troladores quando é feita uma escolha sistemática dos modelos de síntese, como descrito
na Seção 4.2. Um problema realista de pilotage automática de um míssil ar-ar é uti-
lizado para esse fim. O seu modelo não-linear, assim como os dados numéricos, são
disponibilizados em (NICHOLS, 1993) e no Apêndice 8.1.
Mais especificamente, os resultados apresentados em (PELLANDA, 2001) através da
interpolação de realizações balanceadas dos controladores e dos ganhos da sua estrutura
estimação-controle, projetada sobre modelos linearizados em um conjunto de pontos de
operação, são comparados com os resultados desse projeto sobre um conjunto de modelos
linearizados em faixas de operação, conforme proposto na Seção 4.2.
O problema consiste em controlar a aceleração vertical ηc(t) usando um sinal de
controle δc(t) aplicado no profundor. A estrutura de controle em malha fechada é repre-
sentada na FIG. 5.17. A representação dos estados do sistema não-estacionário G(s, θ)
é mostrada em (8.1). Os estados do míssil são o ângulo de ataque α(t), a velocidade
angular em arfagem q(t), o ângulo do profundor δ(t) e sua derivada δ(t). A aceleração
vertical η(t) e a velocidade angular em arfagem são as saídas medidas. As condições
de equilíbrio do sistema são parametrizadas pelo ângulo de ataque (θ = |α|) caso seja
suposto que a velocidade do míssil é constante e igual à Mach 3. Como o modelo do
míssil é simétrico para α = 0, os controladores lineares K(s) são calculados para α ≥ 0
e interpolados em |α|.Em (STILWELL, 1997, 1999), os autores consideraram o método de modelagem
de valores singulares H∞ (H∞ loop shaping), introduzido em (MACFARLANE, 1990),
que incorpora um compromisso entre desempenho e estabilidade robusta para concepção
de controladores LTI. Tal exemplo permite ilustrar outros aspectos desse método, pois
a estrutura do piloto automático não é explicitamente mostrada no diagrama da FIG.
2.1-(a).
Os controladores globais C(s, α) , K(s, α)W (s, α) (FIG. 5.17) são construídos
112
FIG. 5.17: Estrutura de controle para o exemplo do míssil
combinando as funções de ponderação W (s, α) e os controladores H∞ iniciais K(s, α).
Considera-se a síntese LTI para α = [ 0o; 2, 64o ], α = [ 2, 64o; 10, 26o ] e α =
[ 10, 26o; 30o ], faixas selecionadas de operação do parâmetro utilizando o Algoritmo 4.1.
Os controladores K(s, α), calculados pela técnica de síntese H∞ em (MACFARLANE,
1990), dependem do ponto de funcionamento parametrizado por α e são da forma
K(s, α) =
Ap(α) + Bp(α)Kc1(α) + Kf1(α)Cp(α) −Kf1(α)
Kc1(α) 0
(5.14)
onde (Ap, Bp, Cp) são as matrizes de espaço de estado do sistema ponderado (shaped plant)
P (s, α) = W (s, α)G(s, α) e (Kc0, Kf0) são os ganhos de realimentação e do observador de
estado. A função de transferência do filtro é análoga àquela considerada em (STILWELL,
1997, 1999),
W (s, α) =
[g1
s−z1
s0
0 g2s−z2
s−p2
](5.15)
onde os parâmetros dependem de α e estão listados na TAB. 5.13.
TAB. 5.13: Parâmetros do filtro W(s)
α (degrees) g1 z1 g2 z2 p2
1a Faixa 0,15717 -174,23 0,16651 -7,6875 -0,689662a Faixa 0,14188 -183,80 0,20032 -5,6996 -0,876203a Faixa 0,23208 -125,83 0,34126 -3,276 -1,0574
Na aplicação desse método, considera-se primeiro o sistema ponderado P (s, α),
verifica-se quais são as faixas sobre as quais os controladores são projetados (α =
[ 0o; 2, 64o ]; α = [ 2, 64o; 10, 26o ]; α = [ 10, 26o; 30o ]) e calcula-se as transfor-
mações de similaridade T (α). Obtém-se os controladores aumentados Ke(s) que são
113
equivalentes aos controladores H∞ iniciais K(s) (n = nK = 6 e Q(s) = 0). Os con-
troladores globais C(s) = K(s)W (s) e Ce(s) = Ke(s)W (s) são de ordem 8. Os ganhos
dos controladores Kc e Kf , para esse caso, são apresentados na TAB. 5.14. Eles foram
calculados para as faixas de operação tendo sido obtidos através da seleção de partições
dos modos em malha fechada para a primeira faixa α = [ 0o; 2, 64o ] (Seção 2.5.3) e exe-
cutando um procedimento de continuação (Algoritmo 2.1 da Seção 2.5.5) para as faixas
restantes. Pode-se notar que, com algumas exceções para os ganhos do estimador, os
coeficientes são bastante próximos para diferentes faixas de operação.
Note que o filtro W (s) é conectado em cascata com o sistema fazendo parte da trans-
missão direta. Então, é conveniente conservar a estrutura para preservar as propriedades
de síntese LTI quando da interpolação. Este é o caso para Ce(s) uma vez que o filtro é
incorporado à posteriori ao controlador equivalente aumentado.
Ao contrário, nesse exemplo, os controladores LTI equivalentes que comportam pa-
râmetros de Youla dinâmicos e estimadores físicos só podem ser calculados para os con-
troladores globais C(s, α), com nK = 8, n = 4 e nq = 4. A estrutura do filtro seria, neste
caso particular, perdida no processo de interpolação de ganhos e de Q(s). Um contro-
lador interpolado baseado em um tal conjunto de controladores não seria apropriado para
inserção na malha. Entretanto, uma tal estrutura poderia ser explorada indiretamente
para melhorar a estimação física (PELLANDA, 2001).
Os valores dos pólos em malha fechada e sua distribuição para todos os controladores
e os três intervalos de operação são mostrados nas TAB. 5.15 et 5.16, com as seguintes
definiçõesCtr1 , spec(Ap −BpKc1), Obs1 , spec(Ap −Kf1Cp)
Ctr2 , spec(Ap −BpKc2), Obs2 , spec(Ap −Kf2Cp)
Ctr3 , spec(Ap −BpKc3), Obs3 , spec(Ap −Kf3Cp)
Alguns modos quase duplos interferiram na escolha da partição inicial dos pólos. Sua
distribuição é, no entanto, coerente para todos os pontos de operação.
Adota-se as estruturas de controle não-estacionária para Ce(s, α) = Ke(s, α)W (s, α).
Uma interpolação linear dos pólos, zeros, e ganhos de W (s, α) é utilizada de forma a
conservar as propriedades de modelagem dos valores singulares H∞ dentro das zonas de
transições. Três estratégias foram usadas para a interpolação de Ke(s, α) :
a) Os ganhos (Kc(α), Kf (α)) foram linearmente interpolados, enquanto que os out-
ros coeficientes matriciais de Ke(s, α), (Ap(α), Bp(α), Cp(α)), foram calculados via
conexão em série P (s, α) = W (s, α)G(s, α). O modelo linearizado em espaço de
114
estado parametrizado em α (A(α), B(α), C(α)) de G(s, α) foi suposto disponível
em tempo real e utilizado para atualizar o controlador.
b) Uma vez que os controladores LTI são estáveis para todos os pontos de opera-
ção, pôde-se utilizar uma interpolação linear das matrizes da realização balanceada
baseadas nos Gramianos de controlabilidade e observabilidade.
c) A interpolação da faixa 1 para a 2 e da faixa 2 para a 3 ocorre nos intervalos onde o
parâmetro α está entre [ 1, 8o; 4, 5o ] e [ 8, 82o; 14, 82o ], respectivamente. Estes
intervalos foram calculados de acordo com a técnica de linearização proposta. Fora
destes intervalos não há interpolação e os ganhos são mantidos constantes. As
matrizes A(α), B(α) e C(α), que fazem parte do controlador, são atualizadas em
todo o domínio, de forma análoga ao primeiro caso.
Note que o filtro W influencia duplamente Ce.
A FIG. 5.18 mostra a resposta LTI ao degrau para os três pontos de operação esco-
lhidos de forma arbitrária e considerados no controlador escalonado a ser comparado. O
objetivo do escalonamento de ganhos é manter o desempenho a tempo variante o mais
próximo possível do desempenho a tempo invariante apresentado nessa figura.
As FIG. 5.19 e 5.20 mostram as respostas do sistema não-linear em malha fechada,
com os três controladores interpolados, a uma seqüencia de entradas em degrau. A perda
de desempenho em relação ao comportamento estacionário da FIG. 5.18 é nítido, mas
ela é mais destoante para o piloto automático que usa a interpolação das matrizes. Para
este, as respostas no tempo apresentaram uma maior ultrapassagem e um menor fator
de amortecimento. Conforme já apresentado em (PELLANDA, 2001), a interpolação dos
ganhos tem um desempenho mais aproximado da resposta desejada, se comparado com a
interpolação das matrizes. Uma melhor aproximação é obtida com o terceiro controlador
que, além de interpolar os ganhos, usou o sistema de linearização de faixas e a inter-
polação ao fim de cada faixa com o início da faixa seguinte. Constata-se uma melhora
considerável na manutenção de desempenho em malha fechada, para este controlador
interpolado, em especial quando há grandes perturbações para valores maiores de acele-
ração, não sendo constatada piora no desempenho, em relação aos outros controladores,
em nenhum momento. De uma forma geral, a técnica proposta mostra uma menor ultra-
passagem e/ou um menor tempo de acomodação. A FIG. 5.20, que é uma ampliação da
anterior em um dado intervalo, mostra mais claramente esse resultado. Deve-se ressaltar
que a interpolação de matrizes ou de ganhos de controladores LTI projetados para mode-
115
los linearizados em pontos, via continuação de subespaços invariantes, são consideradas
técnicas satisfatórias.
Esse resultado é confirmado pelas FIG. 5.21 e 5.22, que mostram a evolução da ângulo
de ataque, que é também o parâmetro de interpolação. Verifica-se uma larga excursão
desse parâmetro.
Neste exemplo também é possível constatar que, mesmo com a melhoria de desem-
penho para o controlador linearizado em faixas com interpolação dos ganhos, o sinal de
controle não sofreu alteração significativa de amplitude. O fenômeno pode ser visto nas
FIG. 5.23 e 5.24. Observa-se uma discreta atenuação no sinal de controle com a técnica
proposta.
116
TAB. 5.14: Ganhos do controlador e do estimador aumentado
α [Kc]T Kf
0o
a2, 64o
27, 4335−0, 0315−5, 0115−0, 40460, 41230, 0016
1, 0072 −0, 9219−11, 0789 21, 8872−11, 1741 22, 2400−110, 4033 396, 7900−5, 1252 −7, 5098504, 5535 359, 0838
2, 64o
a10, 26o
34, 9753−0, 0337−5, 8022−0, 46900, 46790, 0018
0, 8618 −0, 3907−0, 4491 −0, 0152−0, 6287 0, 3866−191, 4988 433, 5060−5, 7902 −8, 4635677, 6276 385, 1673
10, 26o
a30o
54, 0926−0, 0313−8, 7852−0, 68480, 67510, 0024
0, 3372 0, 45670, 1308 −1, 4901
3, 436.10−3 −1, 1931−129, 5917 340, 2019−9, 1513 −10, 5533
1, 2129.103 −1, 142.102
TAB. 5.15: Pólos em malha fechada
Pólos α = [ 0o 2, 64o ] α = [ 2, 64o 10, 26o ] α = [ 10, 26o 30o ]
1, 2 −102, 23± j110, 07 −101, 59± j110, 40 −98, 588± j116, 213, 4 −104, 46± j107, 77 −104, 22± j108, 02 −102, 88± j109, 745 −97, 273 −88, 516 −107, 326 −20, 675 −20, 746 −36, 313
7, 8 −12, 931± j9, 7525 −14, 126± j11, 423 −16, 142± j11, 2989, 10 −10, 993± j8, 9727 −11, 818± j10, 304 −11, 710± j8, 832411 −0, 6897 −0, 8762 −1, 057212 −0, 6897 −0, 8762 −1, 0574
117
TAB. 5.16: Distribuição dos pólos em malha fechada
Pólos Ctr1 Obs1 Ctr2 Obs2 Ctr3 Obs3
1, 2 * * *3, 4 * * *5 * * *6 * * *
7, 8 * * *9, 10 * * *11 * * *12 * * *
Ctr1 , spec(Ap −BpKc1) Obs1 , spec(Ap −Kf1Cp)
Ctr2 , spec(Ap −BpKc2) Obs2 , spec(Ap −Kf2Cp)
Ctr3 , spec(Ap −BpKc3) Obs3 , spec(Ap −Kf3Cp)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.2
0
0.2
0.4
0.6
0.8
1
1.2
Ace
lera
ção
(g)
Tempo (s)
alfa=0 alfa=15alfa=30
FIG. 5.18: Resposta ao degrau do sistema linearizado em malha fechada em três pontosde operação
118
0 1 2 3 4 5 6 7 8 9 10−50
−40
−30
−20
−10
0
10
20
30
40
Tempo(s)
Interp Ganhos Fx
Interp Ganhos Pt
Interp Matrizes
Referência
FIG. 5.19: Resposta do sistema não-linear em malha fechada para três tipos deinterpolação de controladores - aceleração
4 4.2 4.4 4.6 4.8 5 5.2
0
5
10
15
20
25
30
35
40
Tempo(s)
Ace
lera
ção(
g)
Interp Ganhos FxInterp Ganhos PtInterp MatrizesReferência
FIG. 5.20: Resposta do sistema não-linear em malha fechada para três tipos deinterpolação de controladores - aceleração: ampliação de parte da FIG. 5.19.
119
0 1 2 3 4 5 6 7 8 9 10−25
−20
−15
−10
−5
0
5
10
15
20
25
Tempo(s)
Interp Ganhos (L Faixa)
Interp Ganhos (L Ponto)
Interp Matrizes
FIG. 5.21: Resposta do sistema não-linear em malha fechada para três tipos deinterpolação de controladores - ângulo de ataque
3.8 4 4.2 4.4 4.6 4.8 5 5.2 5.4
−20
−15
−10
−5
0
Tempo(s)
Interp Ganhos (L Faixa)
Interp Ganhos (L Ponto)
Interp Matrizes
FIG. 5.22: Resposta do sistema não-linear em malha fechada para três tipos deinterpolação de controladores - ângulo de ataque: ampliação de parte da FIG. 5.22
120
0 1 2 3 4 5 6 7 8 9 10−30
−20
−10
0
10
20
30
Tempo(s)
Interp Ganhos (L Faixa)Interp Ganhos (L Ponto)Interp Matrizes
FIG. 5.23: Resposta do sistema não-linear em malha fechada para três tipos deinterpolação de controladores - sinal de controle δc
4 4.2 4.4 4.6 4.8 5 5.2
0
5
10
15
20
25
Tempo(s)
Sin
al d
e C
ontr
ole
Interp Ganhos (L Faixa)Interp Ganhos (L Ponto)Interp Matrizes
FIG. 5.24: Resposta do sistema não-linear em malha fechada para três tipos deinterpolação de controladores - sinal de controle δc: ampliação de parte da FIG. 5.23
121
5.5 CONTROLE DE UM SISTEMA DE LEVITAÇÃO MAGNÉTICA POR TÉCNI-
CAS DE ESCALONAMENTO DE GANHOS
Nesta seção, explora-se o comportamento de um sistema do tipo levitador magnético,
em malha fechada, para três tipos de controladores:
• Controlador C1: controlador LTI Proporcional Derivativo (PD) projetado por alo-
cação de pólos de malha fechada.
• Controlador C2: controlador de ganho escalonado por imposição da dinâmica de
malha fechada, projetado conforme a técnica apresentada na Seção 4.3.
• Controlador C3: controlador de ganho escalonado por cancelamento de não-
linearidades, projetado conforme a técnica descrita na Seção 4.4.
O modelo em malha aberta do levitador magnético encontra-se descrito no Apêndice 8.3
desta dissertação.
O objetivo a ser atingido pelo sistema não-linear em malha fechada é reproduzir o
comportamento de malha fechada da planta linear com o controlador linear, que possui
um par de pólos alocados em −30, 51 ± j39, 97. O sinal de entrada, único para todas
simulações, assim como a resposta do sistema linear de referência (rref ) e do sistema
não-linear em malha fechada para os três controladores sob análise, são mostrados na
FIG. 5.25.
De acordo com o padrão de desempenho requerido, o controlador PD (C1) foi proje-
tado para operar com a planta linearizada no ponto x = 0 mm, com seus pólos de malha
fechada alocado em −30, 51±j39, 97. O vetor resposta deste sistema é denominado como
r1.
Os valores complexos −30, 51± j39, 97 adotados para a alocação de pólos no projeto
de C1 são também considerados para a alocação de autovalores nos projetos de C2 e
C3. Os vetores resposta destes sistemas em malha fechada são denominados por r2 e r3,
respectivamente.
A fim ter um padrão numérico de comparação de desempenho baseado na norma 2
da diferença ri−rref (i = 1, 2, 3) entre as respostas, utilizou-se um período de τ = 0, 02s.
Os dados numéricos para essa comparação são os seguintes:
‖r1 − rref‖2 = 0, 0110
‖r2 − rref‖2 = 0, 0035
‖r3 − rref‖2 = 0, 0057
122
0 1 2 3 4 5 6 7 8 9 10−8
−6
−4
−2
0
2
4x 10
−3
Tempo em s
Des
loca
men
to (
mm
)
Sist LinearC1C2C3Referência
FIG. 5.25: Comparação das saídas controladas com o sinal de referência.
Observa-se na FIG. 5.25 que quanto maior a perturbação em degrau na entrada, isto
é, quanto mais se afasta do ponto nominal x = 0, pior é o desempenho do controlador
LTI C1 em relação aos controladores escalonados C2 e C3. A variação no tempo destes
controladores tenta compensar os efeitos das não-linearidades. O resultado da figura vem
corroborar os dados numéricos das diferenças de normas das respostas.
Fazendo-se uma análise similar baseada na norma 2 dos sinais de controle (correntes
na entrada da planta - FIG. 5.26), verifica-se:
‖sC1‖2 = 18, 0265
‖sC2‖2 = 10, 9453
‖sC3‖2 = 11, 3550
onde sC1 , sC2 e sC3 são os vetores dos sinais de controle discretizados proporcionados
pelos respectivos controladores. Além disso, observa-se na FIG. 5.26, que os valores de
pico de cada um dos vetores são: sC1 = 8, 2367, sC2 = 3, 5973 e sC3 = 3, 1809. Isso
mostra que os controladores escalonados compensam de maneira satisfatória os efeitos
123
das não-linearidades sem, contudo, aumentar o esforço de controle e provocar riscos de
saturação.
0 1 2 3 4 5 6 7 8 9 10−10
−5
0
5
Tempo em s
C1C2C3
FIG. 5.26: Comparação entre os sinais de controle
Para o projeto dos controladores de ganho escalonado C2 e C3, foi considerada a
seguinte modelagem quasi -LPV:
f(x) = (3, 29553162.104x1 + 483, 237944595)/m
g(x) = (255, 11488730x1 + 2, 03988183)/m
A(x) =
[0 f(x)
g(x) 0
]
B(x) =
[0
g(x)
]
Aa(x) = A(x)2 + 61, 2A(x) + 2528, 5
[1 0
0 1
]
Be(x) = A(x)2 + 300A(x) + 22600I
124
C(x) =[
B(x) A(x)B(x)]
O(x) = I
Kc(x) =[
0 1]C(x)−1Aa(x)
Kf (x) = Be
[0 1
]T
KC31 = [665, 7681 10, 3501]
KC32 = [95246, 578 737, 326]
onde:
• A(x): matriz dinâmica de malha aberta;
• B(x): matriz de entrada;
• Aa(x): matriz do polinômio característico de realimentação de estados (utilizado
em C2);
• Be(x): matriz do polinômio característico do estimador (utilizado em C2);
• C(x): matriz controlabilidade dependente do parâmetro (utilizado em C2);
• O(x): matriz observabilidade dependente do parâmetro (utilizado em C2);
• Kc(x): ganho dinâmico de realimentação de estados (utilizado em C2);
• Kf (x): ganho dinâmico de estimação (utilizado em C2);
• KC31: ganho de realimentação de estados (utilizado em C3);
• KC32: ganho do sinal suplementar (utilizado em C3);
• m: massa da esfera.
Também, nos projetos de C2 e C3, verificou-se que os autovalores de malha fechada
do sistema quasi -LPV são constantes e iguais a −30, 60 ± j39, 90, para toda a faixa
de variação de x das simulações apresentadas nas figuras desta seção. Estes valores
diferem pouco dos valores desejados: −30, 51 ± j39, 97. Contudo, esta variação já era
esperada, devido à aproximação polinomial da função exponencial presente na dinâmica
do sistema.
125
6 CONCLUSÃO
Neste capítulo é apresentado um resumo, do ponto de vista metodológico, das con-
tribuições desta dissertação. É feita uma crítica dos resultados obtidos, bem como ap-
resentadas as perspectivas para futuras investigações, baseadas na experiência adquirida
durante o desenvolvimento deste trabalho.
6.1 CONTRIBUIÇÕES PARA A IDENTIFICAÇÃO DE SISTEMAS
O Capítulo 3 apresenta contribuições deste trabalho para a identificação de sistemas.
Inicialmente foi apresentada uma técnica que possibilita a identificação de sistemas
lineares do tipo caixa-preta usando vetores de amostras discretos da entrada e saída da
planta identificada. A manipulação dos dados permitiu a obtenção de uma FT discreta
que replica a saída medida. A técnica ajusta os coeficientes da FT pela minimização do
erro médio quadrático entre a saída estimada do modelo identificado e o sinal medido.
Problemas práticos relativos ao condicionamento numérico dos dados foram apontados e
sugestões para contorná-los foram apresentadas.
Uma segunda técnica de identificação, também aplicável a sistemas lineares do tipo
caixa-preta, foi proposta. Neste caso, além do vetor de amostras discretas da entrada e
da saída, foi necessário extrair destes dados a resposta em freqüência do sistema, como
apresentado em (ADES, 2005). O método identifica uma FT no domínio contínuo, cujos
coeficientes são ajustados pelo uso da técnica dos mínimos quadrados.
Em ambas metodologias sugeridas para identificação de modelos lineares constatou-
se a possibilidade de obtenção de FT não-mínimas e instáveis. Para solucionar esses
problemas, algumas considerações matemáticas na manipulação dos dados resultantes se
mostraram úteis e necessária, resultando na identificação de FT estáveis e mínimas.
A terceira técnica de identificação proposta se aplica a sistemas não-estacionários.
Parte-se da premissa que um conjunto suficientemente grande de modelos LTI é conhecido
para vários pontos de operação da planta. A metodologia é aplicada a plantas do tipo
caixa-cinza, onde é necessário conhecer o parâmetro que torna o sistema variante no
tempo e que seja possível a sua medição. Funções não-lineares dependentes do parâmetro
foram adotadas para os coeficientes da FT, que são também ajustados pela minimização
de um critério quadrático do erro. Ao fim do processo de identificação, obtém-se um
126
modelo quasi -LPV que se iguala aos modelos LTI para valores fixos do parâmetro nos
pontos de operação adotados.
Os métodos de identificação propostos possuem solução analítica e são de simples
aplicação, o que constitui uma grande vantagem. Além disso, os exemplos numéricos
indicaram a sua eficácia, tendo sido obtidos resultados análogos ou superiores aos apre-
sentados por outra técnica teoricamente mais sofisticada.
Finalmente, foi observado um aspecto numérico relevante, que merece destaque: o
ajuste do modelo identificado pelo primeiro método se mostrou conjuntamente sensível
ao tempo de amostragem e ao tipo de sinal de excitação.
6.2 CONTRIBUIÇÕES PARA O CONTROLE POR ESCALONAMENTO DE GA-
NHOS
No Capítulo 4 foi tratada a síntese de controladores para modelos não-lineares ou
LPV que evoluem em grandes faixas de operação. Três metodologias são apresentadas
nesse capítulo.
O primeiro método se baseia no trabalho desenvolvido em (PELLANDA, 2000, 2001,
2002b) e buscou dotá-lo de um procedimento sistemático de escolha da família de modelos
de síntese LTI. Foi apresentado um critério de escolha de subintervalos de interpolação e
de linearização para o modelo não-estacionário considerando toda o domínio de operação.
A técnica se baseia na garantia de um erro máximo entre a dinâmica dos modelos line-
arizados por faixa e a dinâmica do modelo não-linear. Para o conjunto de modelos LTI
calculados, é sintetizada uma família correspondente de controladores robustos. Uma se-
qüência de transformações de similaridade, calculada por um método de continuação de
subespaços invariantes, é aplicada aos controladores robustos originais de modo a obter
um conjunto equivalente com estrutura estimação-controle, que se prestam à interpo-
lação. A metodologia fornece uma transição suave de controladores, de forma que haja
pouca perda de desempenho.
A segunda técnica apresentada trata o problema do escalonamento de ganhos clássico
de uma forma distinta da primeira. O que se propõe é utilizar uma abordagem quasi -
LPV para a modelagem do sistema. A idéia principal é estabelecer uma realimentação
de estados estimados onde os ganhos são dependentes do parâmetro definido no modelo
quasi -LPV. Além disso, os ganhos são calculados como funções do parâmetro, de forma
que o espectro de autovalores da matriz dinâmica de malha fechada seja alocado em
uma posição fixa predefinida para todo o conjunto possível de trajetórias contínuas do
127
parâmetro.
Por fim, é proposto um método que objetiva o cancelamento das não-linearidades
através da realimentação de estados. Assim como no caso anterior, é utilizada uma
técnica de alocação de autovalores em malha fechada. A técnica se aplica a uma classe
restrita de sistemas não-estacionários. Basicamente, a lei de controle proposta busca
anular os efeitos das dinâmicas indesejadas e impõe o controle sobre a parte linear do
modelo. Esse método também objetiva impedir que variações bruscas sejam inseridas
artificialmente na dinâmica de malha fechada.
Todos os métodos de escalonamento de ganhos apresentados foram testados em sis-
temas não-lineares realistas, tendo sido constatadas as suas aplicabilidades e as vantagens
do seu uso. Contudo, as técnicas não são facilmente aplicáveis a sistemas mais complexos,
com um número maior de parâmetros ou de equações não-lineares nos seus modelos. Além
disso, duas delas tem aplicabilidade restrita a algumas classes de modelos não-lineares e
requerem adaptações judiciosas destes. Estas desvantagens são, no entanto, amenizadas
quando se leva em consideração que a maior parte das técnicas clássicas de escalonamento
de ganhos existentes não são facilmente aplicáveis a sistemas complexos.
6.3 PERSPECTIVAS
Com relação às técnicas apresentadas neste trabalho, existem alguns pontos que
merecem atenção e que representam possibilidades de trabalhos futuros.
• Os métodos de identificação desenvolvidos podem ainda ser testados em casos onde
o sinal medido esteja contaminado por ruído.
• Uma contribuição significativa à técnica desenvolvida na Seção 3.4 seria o estudo
da possibilidade de ser suprimida a monitoração do parâmetro que caracteriza o
modelo identificado, se este for endógeno. Devido a esta característica, em princípio,
a metodologia pode gerar o sinal internamente. Isso representaria uma identificação
caixa-preta para o modelo quasi -LPV.
• A identificação de modelos quasi -LPV é uma contribuição interessante deste tra-
balho. É possível desenvolver técnicas que identifiquem diretamente modelos LPV
com dependência LFT ou politópica em relação ao parâmetro. Esses modelos são
muito úteis para a síntese de controladores escalonados por técnicas LPV.
• No contexto da interpolação clássica de controladores preconcebidos, utiliza-se nor-
malmente uma transição linear. A busca de caminhos alternativos para a interpo-128
lação de ganhos da estrutura estimação-controle que conduzam a uma certa insen-
sibilidade às suas variações do ponto de vista entrada-saída é uma via de pesquisa
pouco explorada na literatura de controle.
• Uma comparação entre os três métodos de escalonamento de ganhos e com outras
técnicas não discutidas neste trabalho, é uma tarefa a ser executada.
• Estender as técnicas de escalonamento de ganhos para sistemas multivariáveis e
discretos pode ser uma complementação interessante.
129
7 REFERÊNCIAS BIBLIOGRÁFICAS
ABED, E. H., LINDSAY, D. e HASHLAMOUN, W. A. On participation factors for linearsystems. Automatica, 36:1489–1496, 2000.
ADES, R. e VALLE, R. C. Identificação de sistemas por conjuntos pré-selecionados viaimposição de pólos. VII SBAI/ II IEEE LARS, págs. 727–732, Set 2005.
AGUIRRE, L. A. A nonlinear dynamical approach to system identification. IEEE Circuitsand Systems Newsletter Society, 11(2):10–23,47, 2000a.
AGUIRRE, L. A. Introdução à Identificação de Sistemas: Técnicas Lineares e Não-Lineares Aplicada a Sistemas Reais. Editora UFMG, 2004.
AGUIRRE, L. A., DONOSO-GARCIA, P. F. e SANTOS-FILHO, R. Use of a prioriinformation in the identification of global nonlinear models - a case study using a buchconverter. IEEE Trans. Circuts Sust. I, 47(7):1081,1089, 2000b.
APKARIAN, P. e ADAMS, R. Advanced gain-scheduling techniques for uncertain sys-tems. IEEE Trans. Contr. Syst. Technology, 6:21–32, 1998a.
APKARIAN, P. e GAHINET, P. A convex charactezation of gain−scheduling H∞ con-trollers. IEEE Trans. Automat. Contr., 40(5):853–864, 1995a.
APKARIAN, P., GAHINET, P. e BECKER, G. Self-scheduled H∞ control of linearparameter-varying systems: A design example. Automatica, 31:1251–1261, Set 1995b.
APKARIAN, P., PELLANDA, P. e TUAN, H. Mixed H2/H∞ multi-channel linearparameter-varyng control in discrete time. System and Control Letters, 41:333–346,2000.
APKARIAN, P. e TUAN, H. D. Parameterized lmis in control theory. SIAM J. onControl and Optimization, 38:1241–1264, 1998b.
BAKSHI, B. R. e STEPHANOPOULOS, G. Wave-net: a multiresolution, hierarquicalneural network with localized learning. AIChE Journal, 39(1):57–91, 1993.
BANASZUK, A. e HAUSER, J. Approximate feedback linearization: Approximate inte-grating factors. Proceeding of the Americas Control Conferences, 6:2–10, Junho 1994.
BATTILOTTI, S. e CALIFANO, C. A constructive condition for dynamic feedback lin-earization. Systems e Control Lettlers, 52:329–338, 2004.
BECKER, G. Parameter-dependent control of an under-actuated mechanical system.Proc. IEEE Conf. Decision Contr., LA, 1995.
BECKER, G., PACKARD, A., PHILBRICK, D. e BALAS, G. Control of parametrically-dependent linear systems: A single quadratic Lyapunov approach. Proc. Amer. Contr.Conf., San Francisco, CA, págs. 2795–2799, 1993.
130
BILLINGS, S. A. Identification of nonlinear systems - a survey. Em In Proc IEE, art D,págs. 272–285, 1980.
BILLINGS, S. A., CHEN, S. e BACKHOUSE, M. J. The identification of linear andnon-linear models of a turbochaeged automotive diesel engine. Mechanical Systemsand Signal Processing, 3(2):123–142, 1989.
BOHLIN, T. e GRAEBE, S. T. Issues in nonlinear stochastic gray box identification. Int.J. of Adaptative Control and Signal Processing, 9(6):465, 1995.
COCA, D. e BILLINGS, S. A. Estimating derivatives from noisy data: A wavalet mul-tiresolution decomposition approach. Technical report, 1995a.
COCA, D. e BILLINGS, S. A. Smoothing chaotic signal for system identification: Awavalet multiresolution decomposition approach. Technical report, 1995b.
COCA, D. e BILLINGS, S. A. Continuous-time system identification for linear andnonlinear using wavelet decomposition. Technical report, 1996a.
COCA, D. e BILLINGS, S. A. A direct approach to identification of nonlinear differentialmodels from discrete data. Mechanical Systems and Sugnal Processing, 13(5):739–755,1996b.
CORRÊA, M. V. Identificação Caixa-Cinza de Sistemas Não-Lineares Utilizando Repre-sentações NARMAX Racionais e Polinomiais. Tese, Universidade Federal de MinasGerais, Belo Horizonte - MG, Brasil, Fevereiro 2001.
CUBILLOS, F. A. e LIMA, E. L. Identification and optimizing control of a roucherflotation circuit using an adaptable hybrid-neural model. Minerals Engineering, 10(7):707–721, 1997.
DESOER, C. e VIDYASAGAR, M. Feedback Systems: Input-Output Properties. AcademicPress, 1975.
GARCIA, C. Modelagem e simulação de processos industriais e de sistemas eletromecâni-cos. EDUSP, São Paulo, 1997.
GENÇAY, R. e LIU, T. Nonlinear modelling and prediction with feedforward and recur-rent networks. Physica D, (108):119–134, 1997.
GOLDFREY, K. R. Practical problems in identification. Goldfrey, K. and Jones, P., edi-tors, Signal Processing for Control, Lecture Notes in Control and Information Sciences,págs. 358–389, 1986.
GOLUB, G. H. e VAN LOAN, C. F. Matrix computtations. J. Hopkins University, 1996.
HAMDAN, A. M. A. e NAYFEH, A. H. Measures controllability and observability forfrst- and second-order linear systems. J. of Guid. Contr. and Dyn., Maio-Junho 1989.
HERBERT, J. A. e TULLEKEN, A. F. Grey-box modeling and identification usingphysical knowledge and bayesian techniques. Automatica, 29(2):285–308, 1993.
131
HOUGEN, J. O. Measurements and Control Applications for Practicing Engineers. Cah-mers Books, Massachusets, U.S.A., 1972.
HYDER, R. A. e GLOVER, K. The application of scheduled H∞ controllers to a vstolaircraft. IEEE Trans. Automat. Contr., 38(7):1021–39, 1993.
ISIDORI, A. Nonlinear Control Systems: An Introduction. Springer-Verlang, 1989.
ISIDORI, A., KRENER, A. J., GIORGI, C. G. e MONACO, S. Decoupling via feedback:a differential geometric approach. IEEE Trans. Automat. Control 26, págs. 331–345,1981.
ISIDORI, A., MOOG, C. H. e DE LUCA, A. A sufficient condition for full linearizationvia dynamic state feedback. Proceeding of the 25th IEEE CDC, págs. 203–208, 1986.
JANG, H. K. e KIM, K. J. Identification of loadspeaker nonlinearities using the NARMAXtechnique. J. Audio Eng. Soc., 42(1/2):50–59, 1994.
JORGENSEN, S. B. e HANGROS, K. M. Grey-box modeling for control: Qualitativemodels as a unifying framework. Int. J. of Adaptative Control and Signal Processing,9(6):547–562, 1995.
KAILATH, T. Linear Systems. Pretice-Hall, 1980.
KAJIWARA, H., APKARIAN, P. e GAHINET, P. LPV techniques for control of aninverted pendulum. IEEE Control Systems Magazine, 19:44–54, 1999.
KAMLER, I., PASCOAL, A. M., KHARGONEKAR, P. P. e COLEMAN, E. E. A velocityalgorithm for the implementation of gain-scheduled. Automatica, 31(8):1185–1191,1995.
KARPLUS, W. The Spectrum of Mathematical Modelling and System Simulation. Simu-lation os Systems. North-Holland Publishing Company. 1976.
KELLET, M. Continuos scheduling of H∞ controllers for a MS760 paris aircraft, in robustsystems design using H∞ and related methods. P. H. Hammond, págs. 197–219, 1991.
KHALIL, H. Nonlinear Systems. Prentice Hall, 1996.
KUNDUR, P. Power System Stability and Control. McGraw−Hill, Inc, 1994.LAWRENCE, D. A. e RUGH, W. J. Gain scheduling dynamic linear controllers for a
nonlinear plant. Automatica, 31(3):381–390, 1995.
LEITH, D. J. e LEITHEAD, W. E. Appropriate realization of MIMO gain scheduledcontrollers. Int. J. Contr., 70:13–50, 1998a.
LEITH, D. J. e LEITHEAD, W. E. Gain−scheduled and nonlinear systems: Dynamicanalysis by velocity based linearization families. Int. J. Contr, 70:289–317, 1998b.
LEONTARITIS, I. e BILLINGS, S. Input-output parametric models for non-linear sys-tems part I: Deterministic non-linear systems. Int. J. Control, 41(2):303–328, 1985.
132
LJUNG, L. System Identification, Theory for the Use. Prentice-Hall, New Jersey, 1987.
LU, W. M. e DOYLE, J. C. H∞ control of lft systems: An LMI approach. Proc. IEEEConf. Decision Contr., Tucson, AR, págs. 1997–2001, 1992.
LU, W. M. e DOYLE, J. C. H∞ control of nonlinear systems: A convex characterization.IEEE Trans. Aut. Control, 40:1668–1674, 1995.
LUI, S. H., KELLER, H. B. e KWOK, T. W. C. Homotopy method for the large, sparse,real nonsymmetric eigenvalue problem. SIAM, J. Matrix Anal. Appl., 18(2):312–333,1997.
MACFARLANE, D. C. e GLOVER, K. Robust controller design using normalized coprimefactor plant descriptions. Lecture Notes in Control and Information Sciences, 138, 1990.
MARTINS, N. e LIMA, L. T. G. Determination of suitable location dor power systemstabilizers and static VAR compensator for damping electromechanical oscillations inlarge scale power systems. IEEE Trans. Power Systems, 5(4):1455–1469, Setembro1990.
NARENDA, K. S. e PARTHASARATY, K. Identification and control of dynamical sys-tems using neural networks. IEEE Trans. Neural Net, 1:4–27, 1990.
NICHOLS, R. A., REICHERT, R. T. e RUGH, W. J. Gain sheduling for H∞ controllers:a flight control example. IEEE Trans. Contr. Systems Technology, 1(2):69–79, 1993.
NOSHIRO, M. F., LINKENS, D. e GOODE, K. Nonlinear identification of the Pco2
control system in man. Computer Method and Programs in Biomedicine, 40:189–202,1993.
PACKARD, A. Gain-scheduling via linear fractional transformations. System and ControlLetters, 22:79–92, 1994.
PELLANDA, P., APKARIAN, P. e TUAN, H. Missile autopilot design via a multichannelLFT/LPV control method. Int J. Robust and Nonlinear Control, 12(1):1–20, 2002a.
PELLANDA, P. C. Commande de Systèmes Instationnaires: Séquencement de Compen-sateurs et Commande LPV. Tese, École Nationale Supérieure de L’Aéronautique Etde L’Espace, Toulouse, França, 2001.
PELLANDA, P. C. e APKARIAN, P. Conception de Commandes Robustes, chapter UneMéthode d’Interpolation de Structures Estimation/Commande pour des Compensa-teurs H∞ e µ, págs. 269–315. Hermès Science Publications, Paris, 2002b.
PELLANDA, P. C., APKARIAN, P. e ALAZARD, D. Gain-scheduling through contin-uation of observer−based realizations - applications to H∞ e µ controllers. 39th IEEEConf. on Decision and Control, págs. 2787–2792, 2000.
PETTERSSON, S. e LENNARTSON, B. A converse theorem for exponential stabilityusing piecewise quadratic Lyapunov functions. Technical Report CTH/RT/I-97/008,Control Engineering Lab, Chalmers University of Technology, Gothenburg, Sweden,1997.
133
POTTMANN, M. e PEARSON, R. K. Block-oriented narmax models with output mul-tiplicities. AIChE Journal, 44(1):131–140, 1998.
PRÖLL, T. e KARIM, M. N. Model-predictive ph control using real-time NARX ap-proach. AIChE Journal, 40(2):269–282, 1994.
REICHERT, R. T. Dynamic scheduling of modern-robust-control autopilot designs formissiles. IEEE Contr. Syst. Magazine, 12(5):35–42, 1992.
REILLY, J. Development and Analysis of a Nonlinear Dynamic Inverse Control Strategy.Ph.d. thesis, Institute for Systems Research, 1996.
RUGH, W. J. Analytical framework for gain scheduling. IEEE Contr. Syst. Mag., 11(2):79–84, 1991.
SCHERER, C. Robust generalized H2 control for uncertain and LPV systems with generalscalings. Proc. IEEE Conf. on Decision and Control, Kobe, JP, págs. 3970–3975, 1996.
SCHOENWALD, D. A. e ÖZGÜNER, . Optimal control of feedback linearizable systems.Proceeding of the 31st Conference on Decision and Control, 10:15–30, Dezembro 1992.
SHAHRUZ, S. M. e BEHTASH, S. Design of controllers for linear parameter varyingsystems by the gain scheduling techinique. Journal Of Mathematical Analysis andApplications, 168(1):195–217, 1992.
SHAMMA, J. S. e ATHANS, M. Analysis of gain scheduled control for nonlinear plants.IEEE Trans. Automat. Contr., 35(8):898–907, 1990.
SJÖBERG, J., ZHANG, Q., LJUNG, L., BENVENISTE, A., DELYON, B., GLOREN-NEC, P., HJALMARSSON, H. e JUDITSKY, A. Nonlinear black-box modeling insystem identification: A unified overview. Automatica, 31(1691), 1995.
SKOGESTAD, S. e POSTLETHWAITE, I. Multivariable Feedback Control Analysis andDesign. John-Wiley and Sons Inc., 1997.
SOHLBERG, B. Supervision and Control for Industrial Processes - Using Grey BoxModels, Predictive Control and Fault Detection Methods. Springer-Verlag, 1998.
STILWELL, D. J. e RUGH, W. Stability preserving interpolation methods for the syn-thesis of gain scheduled controllers. Tech. Rep. JHU/ECE, págs. 97–18, 1997.
STILWELL, D. J. e RUGH, W. Interpolation of observer state feedback controllers forgain scheduling. IEEE Trans. Automat. Contr, 44(6):1225–1229, Junho 1999.
STILWELL, D. J. e RUGH, W. Stability preserving interpolation methods for the syn-thesis of gain scheduled controllers. Automatica, 36(5):665–671, Maio 2000.
TIKHONOV, A. N. e ARSENIN, V. Y. Solutions of ill-posed problems. Winston, Wash-ington, DC, 1977.
TSOI, A. C. e BACK, A. Discrete time recurrent neural network architectures: A unifyingreview. Neurocomputing, (15):183–223, 1999.
134
TUAN, H. D. e APKARIAN, P. Relaxations of parameterized LMIs with control appli-cations. Int J. Robust and Nonlinear Control, págs. 59–84, 1999.
VALLE, R. C. Identificação no domíno da freqüência por conjuntos pré-selecionadosvia imposição de pólos. Dissertação, Instituto Militar de Engenharia, Rio de Janeiro,Brasil, 2005.
VALLVERDU, M., KORENBERG, M. J. e CAMINAL, P. Model identification of theneural control of the cardiovascular system using NAXMAX models. IEEE Conference,1992.
VIDYASAGAR, M. Nonlinear Systems Analysis. Prentice Hall, 1978.
WU, F., YANG, X., PACKARD, A. e BECKER, G. Induced L2-norm control for LPVsystem with bounded parameter variations rates. Int J. Robust and Nonlinear Control,6:983–998, 1996.
ZHOU, K., DOULER, J. C. e GLOVER, K. Robust and Optimal Control. Prentice Hall,New Jersey, 1996.
135
8 APÊNDICES
136
8.1 APÊNDICE 1: MODELO DO MÍSSIL
O problema de guiagem do canal vertical de um míssil ar-ar, abordado no Capítulo
5, baseia-se no sistema não-linear apresentado em (REICHERT, 1992; NICHOLS, 1993)
e descrito na seqüência.
Os valores e as unidades das constantes presentes no modelo do míssil estão rela-
cionados abaixo (REICHERT, 1992; NICHOLS, 1993):
Kα = 0.7P0180Sπmvs
Kq = 0.7P0180SdπIy
Kz = 0.7P0S
mg
Ax = 0.7P0SCa
m
P0 = 973.3 lbs/ft2 - pressão estática a 20 000 ft
S = 0.44 ft2 - superfície de referência
m = 13.98 slugs - massa
vs = 1036.4 ft/s - velocidade do som a 20 000 ft
d = 0.75 ft - diâmetro
Iy = 182.5 slug.ft2 - momento de inércia em arfagem
Ca = −0.3 - coeficiente de arrasto
ζ = 0.7 - fator de amortecimento do atuador
ωa = 150 rad/s - freqüência natural não-amortecidada do atuador
g = 32.2 ft/s2 - constante de gravidade
an = 0.000103 graus−3 am = 0.000215 graus−3
bn = −0.00945 graus−2 bm = −0.0195 graus−2
cn = −0.1696 graus−1 cm = 0.051 graus−1
dn = −0.034 graus−1 dm = −0.206 graus−1
137
X cpδ X > 0 sistema instávelcp
vetor de velocidade
vetor de elevaçãoelevação da cauda
X < 0 sistema estávelcp
αângulode ataque
ângulo do profundor
centro de gravidade centro de pressão
FIG. 8.1: Diagrama físico do míssil
O modelo do eixo de elevação do míssil, ilustrado na FIG. 8.1, envolve o ângulo de
ataque α(t) (em ◦), a velocidade angular em arfagem q(t) (em ◦/s), o ângulo do profundor
δ(t) (em ◦) e sua derivada δ(t) (em ◦/s). A aceleração normal vertical η(t) (em g) e a
velocidade angular em arfagem são as saídas medidas. Uma descrição quasi-LPV do
modelo do míssil e do atuador é dada por:
α
q
δ
δ
=
Zα 1 Zδ 0
Mα 0 Mδ 0
0 0 0 1
0 0 −ω2a −2ζωa
α
q
δ
δ
+
0
0
0
ω2a
δc
[η
q
]=
[Nα 0 Nδ 0
0 1 0 0
]
α
q
δ
δ
(8.1)
onde
Zα = KαM cos α (anα2 + bn|α|+ cn(2−M/3))
Zδ = KαM cos α dn
Mα = KqM2 (amα2 + bm|α|+ cm(−7 + 8M/3))
Mδ = KqM2dm
Nα = KzM2 (anα2 + bn|α|+ cn(2−M/3))
Nδ = KzM2dn
(8.2)
A descrição não-linear acima representa um míssil voando a uma altitude de 20000
ft. É suposto verdadeiro o desacoplamento dos eixos de rumo e de rolagem.138
A dinâmica da planta pode ser parametrizada por α(t) e M(t), onde o número de
Mach M(t) é uma variável exógena. Devido à simetria do míssil em torno de α = 0,
controladores são projetados para α ≥ 0 e interpolados em |α(t)| e M(t), ou θ(t) =
[|α(t)|, M(t)]T .
Para simulações não-lineares e não-estacionárias, a trajetória temporal do número de
Mach é gerada por
M =1
vs
[−|η| sin(|α|) + AxM2 cos(α)
](8.3)
com M(0) = 4 como um perfil realista, como em (NICHOLS, 1993).
139
8.2 APÊNDICE 2: MODELOS LINEARES IDENTIFICADOS
• Gr11(z) = Gm
11(z)
• Gr12(z) = Gm
12(z)
• Gr13(z) = Gm
13(z)
• Gr14(z) = Gm
14(z)
• Gm15(z) = 0,04950150028922z2−0,01614040694779z−0,03271952599199
z3−1,30626734129036z2−0,34109576438472z+0,64772784779743
• Gm16(z) = 0,04950354453251z2−0,01622916422366z−0,03279766962922
z3−1,30800608617206z2−0,34113881328291z+0,64942398665823
• Gr17(z) = 0,04950290757149z−0.04925583555764
z2−1,97515832655574z+0,97530655124666
• Gr18(z) = 0,04950290498158z−0,04925599150817
z2−1,97516146313926z+0,97530960955677
TAB. 8.1: Valores dos coeficientes da FT Gm21(z)
Graus em z Coeficientes do numerador Coeficientes do denominador
z7 0,09994056772890 1z6 -0,08606556006252 -1,01583596455477z5 -0,04297535274090 -0,43006525404615z4 -0,00832012437565 -0,01622042321447z3 0,01493943961896 0,22056999286446z2 0,02370524493145 0,27562612447542z1 0,01484957953608 0,14494925398394z0 -0,01575333496272 -0,17486682191075
TAB. 8.2: Valores dos coeficientes da FT Gr21(z)
Graus em z Coeficientes do numerador Coeficientes do denominador
z3 0,09994056772890 1z2 -0,24285234034599 -2,57279809923296z1 0,19381783379975 2,15291103836761z0 -0,05084455783893 -0,57930627117792
140
TAB. 8.3: Valores dos coeficientes da FT Gm22(z)
Graus em z Coeficientes do numerador Coeficientes do denominador
z7 0,10000324302328 1z6 -0,21566066548438 -2,31166057975892z5 0,06083387070139 0,80896361264880z4 0,07706066967685 0,84999690396296z3 0,03450422028984 0,29375174201215z2 -0,02674762056180 -0,36566199002611z1 -0,06799132679234 -0,72341046003767z0 0,03799874353338 0,44804385533653
TAB. 8.4: Valores dos coeficientes da FT Gr22(z)
Graus em z Coeficientes do numerador Coeficientes do denominador
z4 0.10000324302328 1z3 -0.37294026756444 -3.88443512033357z2 0.52151902596859 5.65976920269300z1 -0.32411604060207 -3.66592768344410z0 0.07553430089659 0.89059892689331
141
TAB. 8.5: Valores dos coeficientes da FT Gm23(z)
Graus em z Coeficientes do numerador Coeficientes do denominador
z15 0,10000209457894 1z14 -0,10113143078328 -1,16639595695617z13 -0,04108313745854 -0,38781509208806z12 -0,00947105803537 -0,01059232042565z11 0,00847431750173 0,17247813392812z10 0,01623852890181 0,22506658278363z9 0,01900274846828 0,21876916117049z8 0,01691603774180 0,16413162943337z7 0,01301309603378 0,09908415635267z6 0,00654462155674 0,01800155941323z5 -0,00137797085309 -0,06518916979802z4 -0,00890146258836 -0,13137563427265z3 -0,01556653854061 -0,17838765081544z2 -0,01592297128421 -0,15452814065485z1 -0,00633343759552 -0,03429602746180z0 0,01960505691829 0,23115707077226
TAB. 8.6: Valores dos coeficientes da FT Gr23(z)
Graus em z Coeficientes do numerador Coeficientes do denominador
z4 0,10000209457894 1z3 -0,38158983368371 -3,97094580549617z2 0,54589140897210 5,91688491371750z1 -0,34699414016830 -3,92087879117162z0 0,08269064420559 0,97494190007441
TAB. 8.7: Valores dos coeficientes da FT Gm24(z) = Gr
24(z)
Graus em z Coeficientes do numerador Coeficientes do denominador
z4 0,09999996091961 1z3 0,07557471500680 -2,31788726012157z2 -0,34209668578165 2,76405692271227z 0,24004472255384 -1,81456386497634z0 -0,05412366295438 0,61261946201281
142
TAB. 8.8: Valores dos coeficientes da FT Gm81(s)
Graus em s Coeficientes do numerador Coeficientes do denominador
s4 0,100000000019 1s3 3,559261194897 4,88261197345s2 44,065817853652 159,79688533644s1 211,641443239558 375,74404367114s0 279,859916308965 3525,23239538712
TAB. 8.9: Valores dos coeficientes da FT Gm82(s)
Graus em s Coeficientes do numerador Coeficientes do denominador
s5 0,099999999997 1s4 3,545152400109 4,74152446811s3 43,565275974116 159,12434852911s2 205,483217599831 353,43948662078s1 250,871276216287 3475,00449293541s0 -40,040168369402 -476,00239481585
TAB. 8.10: Valores dos coeficientes da FT Gr82(s)
Graus em s Coeficientes do numerador Coeficientes do denominador
s4 0,099999999997 1s3 3,559300542403 4,87653626641s2 44,048978476956 159,78273845992s1 211,533722831513 375,01204147754s0 279,673443251354 3525,63554303942
143
TAB. 8.11: Valores dos coeficientes da FT Gm83(s)
Graus em s Coeficientes do numerador Coeficientes do denominador
s7 0,100000000075 1s6 3,660859658761 5,89859588099s5 47,663507019104 164,57268696623s4 255,748921460203 537,24106475674s3 486,827825951954 3879,50285078770s2 261,090865013284 3528,41880950224s1 -74,542118554929 -533,56060769679s0 -2,999610751776 -169,44791044336
TAB. 8.12: Valores dos coeficientes da FT Gr83(s)
Graus em s Coeficientes do numerador Coeficientes do denominador
s6 0,100000000075 1s5 3,684469738535 6,15803792308s4 48,605037385856 166,17034090025s3 267,971243650193 580,35263733440s2 554,996077026532 4030,07072414902s1 395,674141374724 4573,98858793540s0 17,437157086037 653,12433204538
144
8.3 APÊNDICE 3: LINEARIZAÇÃO DO MODELO DE UM LEVITADOR MAG-
NÉTICO
A FIG. 8.2 mostra fotos do sistema de levitação magnética de uma esfera de aço
implementado no Laboratório de Controle do Instituto Militar de Engenharia. A força
eletromagnética é usada para levitar uma esfera de metal. As variáveis utilizadas no
modelo matemático são mostradas na FIG. 8.2.
FIG. 8.2: Sistema de levitação magnética de uma esfera de aço
A equação de movimento da esfera, derivada da Segunda Lei de Newton, é encontrada
da seguinte maneira:∑
F = ma = mx =⇒ f(x, i)−mg = mx (8.4)
145
−4
−2
0
2
4
−1−0.5
00.5
10
2
4
6
8
10
x(mm)i(A)
For
ça(N
)
FIG. 8.3: Valores medidos experimentalmente em laboratório: f(N)× x(mm)× i(A)
onde f(x, i) é a força de sentido contrário à da gravidade ocasionada pelo campo eletro-
magnético do levitador, m é a massa e x é a posição da esfera, g é a aceleração da
gravidade e i é a corrente na bobina do levitador.
A não-linearidade do problema se encontra no fato de que, teoricamente, a força
eletromagnética decai numa relação inversamente proporcional ao quadrado da distância.
A relação exata é difícil de ser encontrada através de princípios físicos, tendo em vista a
grande complexidade da teoria eletromagnética.
Experimentalmente, foram levantados os dados constantes na FIG. 8.3 e na TAB.
8.13, referentes à força, deslocamento e corrente para o levitador magnético. Observa-se
a função f(x, i) é não-linear.
Dados do sistema:
a) distância nominal: x0 = 19 mm;
b) corrente nominal: i0 = 0, 750 A;
c) massa da esfera: m = 0, 346 kg;
d) f(x, i) não-linear: força em função do deslocamento x e da corrente i em relação
aos valores nominais (x0 e i0)
f(x, i) = k1eaxi + k2e
bx (8.5)
146
TAB. 8.13: Diferenças entre as curvas do gráfico força × deslocamento (não-linear)
x(mm) 4 3 2 1 0 -1 -2 -3 -4Intervalos (A) ∆f4 ∆f3 ∆f2 ∆f1 ∆f0 ∆f(−1) ∆f(−2) ∆f(−3) ∆f(−4)
-0,75 / -0,60 0,35 0,3 0,26 0,23 0,2 0,18 0,16 0,14 0,12-0,60 / -0,45 0,55 0,49 0,43 0,38 0,34 0,3 0,26 0,23 0,21-0,45 / -0,30 0,46 0,40 0,35 0,31 0,27 0,25 0,23 0,21 0,18-0,30 / -0,15 0,41 0,35 0,32 0,27 0,24 0,2 0,17 0,15 0,13-0,15 / 0,00 0,43 0,38 0,32 0,29 0,26 0,23 0,2 0,17 0,160,00 / 0,15 0,51 0,45 0,39 0,35 0,3 0,27 0,24 0,22 0,20,15 / 0,30 0,52 0,46 0,41 0,36 0,33 0,29 0,26 0,23 0,20,30 / 0,45 0,55 0,45 0,37 0,31 0,25 0,2 0,17 0,13 0,10,45 / 0,60 0,62 0,55 0,50 0,44 0,38 0,36 0,32 0,3 0,270,60 / 0,75 0,51 0,46 0,40 0,36 0,34 0,29 0,27 0,25 0,23
e) análise para o intervalo: −4 mm ≤ x ≤ 4 mm e −0, 750 A ≤ i ≤ 0, 750 A.
f) ponto de equilíbrio (peso = mg igual à força = f(x, i)) : f(x, i) = 3, 46 N , em
x = 0 mm e i = 0 A.
Tendo em vista este comportamento não-linear da força (em relação ao deslocamento
e à corrente), torna-se recomendável a linearização do sistema de levitação magnética em
torno de um ponto de operação.
Analiticamente, o comportamento de f(x, i) pode ser descrito por uma combinação
de exponenciais (uma vez que a força parece variar exponencialmente com o deslocamento
e quase linearmente com a corrente (ver gráficos de f(0, i) e f(x, 0, 75)).
Primeiramente, deve-se calcular os coeficientes k1, k2, a e b de (8.5):
a = 0, 1275 b = 0, 1332 k1 = 1, 9437 k2 = 3, 4843
Portanto, a força pode ser descrita pela expressão:
f(x, i) = 1, 9437e0,1275xi + 3, 4843e0,1332x (8.6)
Tendo em vista a teoria sobre linearização, pode-se aplicá-la, agora, à expressão de
f(x, i), de maneira a se obter o modelo linearizado para o levitador magnético.
- Primeiramente, determina-se o ponto de operação (x0, i0) como sendo (0, 0);
- Aplica-se a expansão em série de Taylor para f(x, i), mantendo-se apenas os termos
147
de primeira ordem:
f(x, i)− f(x0, i0) = Ka(x− x0) + Kb(i− i0)
f(x0, i0) = 1, 9437e0,1275x0 + 3, 4843e0,1332x = 3, 4843
Ka = ∂f∂x|x0,i0= k1aeax0i0 + k2be
bx0 = k1ai0 + k2b = 0, 4641N/mm = 464, 1N/m
Kb = ∂f∂i|x0,i0= k1e
ax0 = k1 = 1, 9437N/A
f(x, i) = 3, 4843 + 0, 4641(x− x0) + 1, 9437(i− i0)
f(x, i) = 3, 4843 + 464, 1x + 1, 9437i
(8.7)
onde i e x representam as variações de corrente e deslocamento, respectivamente, em
torno do ponto de operação (x0, i0).
148